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	pam-fi	1.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