CaloElectron/src/CaloShape.cpp

#include <CaloShape.h>
////////////////////////////////////////////////////////////
/// INITIALIZATIONS & GLOBAL PARAMETERS
////////////////////////////////////////////////////////////

CaloShape_parameters * CaloShape_parameters::_parameters = 0; 

/**
 * Set external parameters to default values
 */
void CaloShape_parameters::SetDefault(){
    return;
};
////////////////////////////////////////////////////////////
/// CLASS IMPLEMENTATION
////////////////////////////////////////////////////////////


void CaloShape::Init(){
    for(int ip=0; ip<23; ip++){
        qtot[ip]=0.;
        delta[ip]=0.;
        rms[ip]=0.;
        skew[ip]=0.;
        kurt[ip]=0.;
        qtrk[ip]=0.;
        mean[ip]=0.;
        nstrip[ip]=0;
    }
    qmax    = 0.;
    qtotmax = 0.;
    pmax    = -1;
}


bool CaloShape::Set(PamLevel2* l2, int ntr, int view, Bool_t usemechanicalalignement){
    if(!l2)return false;
    CaloTrkVar *c = l2->GetTrack(ntr)->GetCaloTrack();
    if(!c)return false;
    CaloLevel1* l1 = l2->GetCaloLevel1();
    if(!l1)return false;    
    Set(l1,c,view,usemechanicalalignement);
    return true;
};

bool CaloShape::Set(CaloLevel1* ca1, CaloTrkVar* ca2, int view, Bool_t usemechanicalalignement){

    Reset();
//      cout << " CaloEvent::Set()"<<endl; 

    if(!ca1)return false;
    CaloStrip cstrip;
    cstrip = CaloStrip(ca1,usemechanicalalignement);  

    for(int ih=0;ih<ca1->istrip;ih++){
            
        int iv=-1;
        int ip=-1;
        int is=-1;
        ca1->DecodeEstrip(ih,iv,ip,is);

        // ----------------------
        // track intersection
        // ----------------------
        int is0 = -100;
        if(ca2) is0 = ca2->tibar[ip][iv] - 1;

        if(iv==view){
            cstrip.Set(view,ip,is);
            float fxy = 0.;
            if( !view ) fxy = (float) (cstrip.GetX());
            else        fxy = (float) (cstrip.GetY());
            float fq = (float) (cstrip.GetE());         
            nstrip[ip]++;
            qtot[ip]  += fq;
            delta[ip] += fq*fxy;
            rms[ip]   += fq*fxy*fxy;
            skew[ip]  += fq*fxy*fxy*fxy;
            kurt[ip]  += fq*fxy*fxy*fxy*fxy;
            if( is==is0 || is==is0-1 || is==is0+1 )qtrk[ip]+=fq;
            if(fq>qmax)qmax = fq;
                
        }       
    }// end loop over hit strips

/*      cout << "===============" <<endl; */
/*      cout << qtot[4] << endl;  */
/*      cout << delta[4] << endl;  */
/*      cout << rms[4] << endl;  */
/*      cout << skew[4] << endl;  */
/*      cout << kurt[4] << endl;  */

    for(int ip=0; ip<22; ip++){

        float mu1t = 0.;
        float mu2t = 0.;
        float mu3t = 0.;
        float mu4t = 0.;
        if( qtot[ip]>0. ){

            mu1t = delta[ip]/qtot[ip];
            mu2t = rms[ip]/qtot[ip];
            mu3t = skew[ip]/qtot[ip];
            mu4t = kurt[ip]/qtot[ip];           

/*              delta[ip] = mu1t; */
/*              rms[ip] = mu2t - mu1t*mu1t; */
/*              skew[ip] = mu3t - 3*mu1t*mu2t + 2*mu1t*mu1t*mu1t; */
/*              kurt[ip] = mu4t - 4*mu3t*mu1t + 6*mu2t*mu1t*mu1t - 3*mu1t*mu1t*mu1t*mu1t; */

            mean[ip] = mu1t;
            // ----------------------
            // track intersection
            // ----------------------
            float fxy0 = mu1t;
            if(ca2)fxy0 = (float) (ca2->tbar[ip][view]);
            //
            delta[ip] = mu1t - fxy0;
            rms[ip]   = mu2t - 2*mu1t*fxy0 + fxy0*fxy0;
            skew[ip]  = mu3t - 3*mu2t*fxy0 + 3*mu1t*fxy0*fxy0 - fxy0*fxy0*fxy0;
            kurt[ip]  = mu4t - 4*mu3t*fxy0 + 6*mu2t*fxy0*fxy0 - 4*mu1t*fxy0*fxy0*fxy0 + fxy0*fxy0*fxy0*fxy0;

        }

//          if( nstrip[ip]>2 ){
        nstrip[22] += nstrip[ip];
        qtot[22]   += qtot[ip];
        delta[22]  += qtot[ip] * delta[ip];
        rms[22]    += qtot[ip] * rms[ip];
        skew[22]   += qtot[ip] * skew[ip];
        kurt[22]   += qtot[ip] * kurt[ip];
        mean[22]   += qtot[ip] * mean[ip];
        qtrk[22]   += qtrk[ip];
//          }

        if(rms[ip]>0.)rms[ip] = sqrt( rms[ip] );
        else          rms[ip] = 0.;
        if(rms[ip]>0.)skew[ip] = skew[ip]/pow(rms[ip],3.);
//      if(rms[ip]>0.)skew[ip] = skew[ip];
        else          skew[ip] = 0.;
        if(rms[ip]>0.)kurt[ip] = kurt[ip]/pow(rms[ip],4.) - 3.;
//      if(rms[ip]>0.)kurt[ip] = kurt[ip]; 
        else          kurt[ip] = 0.;

        if( qtot[ip]>qtotmax ){
            qtotmax = qtot[ip];
            pmax = ip;
        }
           
    }   

    if(qtot[22]>0.)mean[22] = mean[22]/qtot[22];
    else           mean[22] = 0.;
    if(qtot[22]>0.)delta[22] = delta[22]/qtot[22];
    else           delta[22] = 0.;
    if(rms[22]>0.)rms[22] = sqrt( rms[22]/qtot[22] );
    else          rms[22] = 0.;
    if(rms[22]>0.)skew[22] = skew[22]/( qtot[22]*pow(rms[22],3.) );
//    if(rms[22]>0.)skew[22] = skew[22]/( qtot[22] ) ;
    else          skew[22] = 0.;
    if(rms[22]>0.)kurt[22] = kurt[22]/( qtot[22]*pow(rms[22],4.) ) - 3.;
//    if(rms[22]>0.)kurt[22] = kurt[22]/( qtot[22] );
    else          kurt[22] = 0.;
        

/*      cout << "===============" <<endl; */
/*      cout << qtot[4]  << endl;  */
/*      cout << delta[4] << endl;  */
/*      cout << rms[4]   << endl;  */
/*      cout << skew[4]  << endl;  */
/*      cout << kurt[4]  << endl;  */

    return true;

};
///
void CaloShape::Print(){

    cout << endl;
//    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
    cout << setw(3)<< "IP";
    cout << setw(5)<< "N";
    cout << setw(13)<< "Q";
    cout << setw(13)<< "DX";
    cout << setw(13)<< "RMS";
    cout << setw(13)<< "skew";
    cout << setw(13)<< "kurt";
    cout << endl;
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
    for(int ip=0; ip<=22; ip++){
        cout << setw(3)<< ip;
        cout << setw(5)<< nstrip[ip];
        cout << setw(13)<< qtot[ip];
        cout << setw(13)<< delta[ip];
        cout << setw(13)<< rms[ip];
        cout << setw(13)<< skew[ip];
        cout << setw(13)<< kurt[ip];
        cout << endl;
    }
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
    cout << endl;
}


ClassImp(CaloShape);
ClassImp(CaloShape_parameters);

1	#include <CaloShape.h>
2	////////////////////////////////////////////////////////////
3	/// INITIALIZATIONS & GLOBAL PARAMETERS
4	////////////////////////////////////////////////////////////
5
6	CaloShape_parameters * CaloShape_parameters::_parameters = 0;
7
8	/**
9	* Set external parameters to default values
10	*/
11	void CaloShape_parameters::SetDefault(){
12	return;
13	};
14	////////////////////////////////////////////////////////////
15	/// CLASS IMPLEMENTATION
16	////////////////////////////////////////////////////////////
17
18
19	void CaloShape::Init(){
20	for(int ip=0; ip<23; ip++){
21	qtot[ip]=0.;
22	delta[ip]=0.;
23	rms[ip]=0.;
24	skew[ip]=0.;
25	kurt[ip]=0.;
26	qtrk[ip]=0.;
27	mean[ip]=0.;
28	nstrip[ip]=0;
29	}
30	qmax = 0.;
31	qtotmax = 0.;
32	pmax = -1;
33	}
34
35
36	bool CaloShape::Set(PamLevel2* l2, int ntr, int view, Bool_t usemechanicalalignement){
37	if(!l2)return false;
38	CaloTrkVar *c = l2->GetTrack(ntr)->GetCaloTrack();
39	if(!c)return false;
40	CaloLevel1* l1 = l2->GetCaloLevel1();
41	if(!l1)return false;
42	Set(l1,c,view,usemechanicalalignement);
43	return true;
44	};
45
46	bool CaloShape::Set(CaloLevel1* ca1, CaloTrkVar* ca2, int view, Bool_t usemechanicalalignement){
47
48	Reset();
49	// cout << " CaloEvent::Set()"<<endl;
50
51	if(!ca1)return false;
52	CaloStrip cstrip;
53	cstrip = CaloStrip(ca1,usemechanicalalignement);
54
55	for(int ih=0;ih<ca1->istrip;ih++){
56
57	int iv=-1;
58	int ip=-1;
59	int is=-1;
60	ca1->DecodeEstrip(ih,iv,ip,is);
61
62	// ----------------------
63	// track intersection
64	// ----------------------
65	int is0 = -100;
66	if(ca2) is0 = ca2->tibar[ip][iv] - 1;
67
68	if(iv==view){
69	cstrip.Set(view,ip,is);
70	float fxy = 0.;
71	if( !view ) fxy = (float) (cstrip.GetX());
72	else fxy = (float) (cstrip.GetY());
73	float fq = (float) (cstrip.GetE());
74	nstrip[ip]++;
75	qtot[ip] += fq;
76	delta[ip] += fq*fxy;
77	rms[ip] += fqfxyfxy;
78	skew[ip] += fqfxyfxy*fxy;
79	kurt[ip] += fqfxyfxyfxyfxy;
80	if( is==is0 \|\| is==is0-1 \|\| is==is0+1 )qtrk[ip]+=fq;
81	if(fq>qmax)qmax = fq;
82
83	}
84	}// end loop over hit strips
85
86	/* cout << "===============" <<endl; */
87	/* cout << qtot[4] << endl; */
88	/* cout << delta[4] << endl; */
89	/* cout << rms[4] << endl; */
90	/* cout << skew[4] << endl; */
91	/* cout << kurt[4] << endl; */
92
93	for(int ip=0; ip<22; ip++){
94
95	float mu1t = 0.;
96	float mu2t = 0.;
97	float mu3t = 0.;
98	float mu4t = 0.;
99	if( qtot[ip]>0. ){
100
101	mu1t = delta[ip]/qtot[ip];
102	mu2t = rms[ip]/qtot[ip];
103	mu3t = skew[ip]/qtot[ip];
104	mu4t = kurt[ip]/qtot[ip];
105
106	/* delta[ip] = mu1t; */
107	/* rms[ip] = mu2t - mu1tmu1t; /
108	/* skew[ip] = mu3t - 3mu1tmu2t + 2mu1tmu1tmu1t; /
109	/* kurt[ip] = mu4t - 4mu3tmu1t + 6mu2tmu1tmu1t - 3mu1tmu1tmu1tmu1t; /
110
111	mean[ip] = mu1t;
112	// ----------------------
113	// track intersection
114	// ----------------------
115	float fxy0 = mu1t;
116	if(ca2)fxy0 = (float) (ca2->tbar[ip][view]);
117	//
118	delta[ip] = mu1t - fxy0;
119	rms[ip] = mu2t - 2mu1tfxy0 + fxy0*fxy0;
120	skew[ip] = mu3t - 3mu2tfxy0 + 3mu1tfxy0fxy0 - fxy0fxy0*fxy0;
121	kurt[ip] = mu4t - 4mu3tfxy0 + 6mu2tfxy0fxy0 - 4mu1tfxy0fxy0fxy0 + fxy0fxy0fxy0fxy0;
122
123	}
124
125	// if( nstrip[ip]>2 ){
126	nstrip[22] += nstrip[ip];
127	qtot[22] += qtot[ip];
128	delta[22] += qtot[ip] * delta[ip];
129	rms[22] += qtot[ip] * rms[ip];
130	skew[22] += qtot[ip] * skew[ip];
131	kurt[22] += qtot[ip] * kurt[ip];
132	mean[22] += qtot[ip] * mean[ip];
133	qtrk[22] += qtrk[ip];
134	// }
135
136	if(rms[ip]>0.)rms[ip] = sqrt( rms[ip] );
137	else rms[ip] = 0.;
138	if(rms[ip]>0.)skew[ip] = skew[ip]/pow(rms[ip],3.);
139	// if(rms[ip]>0.)skew[ip] = skew[ip];
140	else skew[ip] = 0.;
141	if(rms[ip]>0.)kurt[ip] = kurt[ip]/pow(rms[ip],4.) - 3.;
142	// if(rms[ip]>0.)kurt[ip] = kurt[ip];
143	else kurt[ip] = 0.;
144
145	if( qtot[ip]>qtotmax ){
146	qtotmax = qtot[ip];
147	pmax = ip;
148	}
149
150	}
151
152	if(qtot[22]>0.)mean[22] = mean[22]/qtot[22];
153	else mean[22] = 0.;
154	if(qtot[22]>0.)delta[22] = delta[22]/qtot[22];
155	else delta[22] = 0.;
156	if(rms[22]>0.)rms[22] = sqrt( rms[22]/qtot[22] );
157	else rms[22] = 0.;
158	if(rms[22]>0.)skew[22] = skew[22]/( qtot[22]*pow(rms[22],3.) );
159	// if(rms[22]>0.)skew[22] = skew[22]/( qtot[22] ) ;
160	else skew[22] = 0.;
161	if(rms[22]>0.)kurt[22] = kurt[22]/( qtot[22]*pow(rms[22],4.) ) - 3.;
162	// if(rms[22]>0.)kurt[22] = kurt[22]/( qtot[22] );
163	else kurt[22] = 0.;
164
165
166	/* cout << "===============" <<endl; */
167	/* cout << qtot[4] << endl; */
168	/* cout << delta[4] << endl; */
169	/* cout << rms[4] << endl; */
170	/* cout << skew[4] << endl; */
171	/* cout << kurt[4] << endl; */
172
173	return true;
174
175	};
176	///
177	void CaloShape::Print(){
178
179	cout << endl;
180	// cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
181	cout << setw(3)<< "IP";
182	cout << setw(5)<< "N";
183	cout << setw(13)<< "Q";
184	cout << setw(13)<< "DX";
185	cout << setw(13)<< "RMS";
186	cout << setw(13)<< "skew";
187	cout << setw(13)<< "kurt";
188	cout << endl;
189	cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
190	for(int ip=0; ip<=22; ip++){
191	cout << setw(3)<< ip;
192	cout << setw(5)<< nstrip[ip];
193	cout << setw(13)<< qtot[ip];
194	cout << setw(13)<< delta[ip];
195	cout << setw(13)<< rms[ip];
196	cout << setw(13)<< skew[ip];
197	cout << setw(13)<< kurt[ip];
198	cout << endl;
199	}
200	cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
201	cout << endl;
202	}
203
204
205	ClassImp(CaloShape);
206	ClassImp(CaloShape_parameters);
207