#include //////////////////////////////////////////////////////////// /// 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()"<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 << "===============" <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 << "===============" <