AlgShowerEM Class Reference

#include <AlgShowerEM.h>

Inheritance diagram for AlgShowerEM:

List of all members.

Public Member Functions
	AlgShowerEM ()
virtual	~AlgShowerEM ()
virtual void	RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void	Trace (const char *c) const

---

Constructor & Destructor Documentation

AlgShowerEM::AlgShowerEM (   )

Definition at line 43 of file AlgShowerEM.cxx. { }

AlgShowerEM::~AlgShowerEM (   )  [virtual]

Definition at line 48 of file AlgShowerEM.cxx. { }

---

Member Function Documentation

void AlgShowerEM::RunAlg (  AlgConfig &  ac, CandHandle &  ch, CandContext &  cx )  [virtual]

Read in info: Implements AlgBase. Definition at line 53 of file AlgShowerEM.cxx. References abs(), CandShowerHandle::AddCluster(), CandHandle::AddDaughterLink(), AlgReco::Calibrate(), CandRecoHandle::GetCharge(), CandStripHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), Registry::GetDouble(), Registry::GetInt(), Calibrator::GetMIP(), CandClusterHandle::GetNStrip(), CandStripHandle::GetPlane(), CandStripHandle::GetStrip(), CandStripHandle::GetStripEndId(), CandStripHandle::GetTime(), CandStripHandle::GetTPos(), CandHandle::GetVldContext(), CandStripHandle::GetZPos(), Calibrator::Instance(), MSG, CandShowerEMHandle::SetAvgDev(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandShowerEMHandle::SetEigenValues(), CandShowerEMHandle::SetEigenVectors(), CandShowerHandle::SetEnergy(), CandShowerEMHandle::SetOutPH(), CandShowerEMHandle::SetShwStatus(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), and CandRecoHandle::SetVtxZ(). { MSG("ShowerEM", Msg::kDebug) << "Starting AlgShowerEM::RunAlg()" << endl; Calibrator& calibrator=Calibrator::Instance(); assert(ch.InheritsFrom("CandShowerEMHandle")); CandShowerEMHandle &csh = dynamic_cast<CandShowerEMHandle &>(ch); assert(cx.GetDataIn()); assert(cx.GetDataIn()->InheritsFrom("TObjArray")); const TObjArray *tary = dynamic_cast<const TObjArray*>(cx.GetDataIn()); Int_t nstp = 0; for (Int_t i=0; i<=tary->GetLast(); i++) { TObject *tobj = tary->At(i); assert(tobj->InheritsFrom("CandClusterHandle")); CandClusterHandle *clusterhandle = dynamic_cast<CandClusterHandle*>(tobj); nstp+=clusterhandle->GetNStrip(); } if(nstp==0){ MSG("ShowerEM", Msg::kError) << "Found empty daughter list in CandShowerEM" << endl; csh.SetDirCosU(0.); csh.SetDirCosV(0.); csh.SetDirCosZ(1.); csh.SetVtxPlane(-1); csh.SetVtxZ(-1.); csh.SetVtxU(0.); csh.SetVtxV(0.); csh.SetVtxT(0.); csh.SetEnergy(0); csh.SetEigenVectors(0); csh.SetEigenValues(0); csh.SetOutPH(0); csh.SetAvgDev(0); csh.SetShwStatus(0); return; } Int_t ShwStatus = 1; Double_t fVertex[3] = {0}; Double_t fVtxPlane = -1; Double_t fAngle[4] = {0}; Double_t *fEigenVector = new Double_t[8]; fEigenVector[0]=0; fEigenVector[1]=0; fEigenVector[2]=0; fEigenVector[3]=0; fEigenVector[4]=0; fEigenVector[5]=0; fEigenVector[6]=0; fEigenVector[7]=0; Double_t *fEigenValue = new Double_t[4]; fEigenValue[0]=0; fEigenValue[1]=0; fEigenValue[2]=0; fEigenValue[3]=0; Double_t *fOutPH = new Double_t[5]; //0-Rad, 1-Lon, 2-Tmax, 3-MaxPl, 4-Redo fOutPH[0]=0; fOutPH[1]=0; fOutPH[2]=0; fOutPH[3]=0; fOutPH[4]=0; Double_t *fAvgDev = new Double_t[4]; fAvgDev[0]=-1.; fAvgDev[1]=-1.; fAvgDev[2]=-1.; fAvgDev[3]=-1.; Double_t fMaxPlane = -1; Double_t fMaxStrip = -1; Double_t fMinStrip = 999; Double_t fMinTime = 999999999; Int_t fIter = ac.GetInt("NIter"); Double_t fCutOff = ac.GetDouble("CutOff"); Double_t fShwFrac = ac.GetDouble("ShwFrac"); Double_t fPECut = ac.GetDouble("PECut"); Double_t fRadCut = ac.GetDouble("RadCut"); Int_t fMaxBigPlaneSep = ac.GetInt("MaxBigPlaneSep"); Int_t fIsolatedPlaneCut = ac.GetInt("IsolatedPlaneCut"); Int_t fNGapPlanes = ac.GetInt("NGapPlanes"); Int_t fNGapStrips = ac.GetInt("NGapStrips"); Double_t fMaxAvgDev = ac.GetDouble("MaxAvgDev"); Double_t fMaxOffFrac = ac.GetDouble("MaxOffFrac"); Double_t fAvgDevCutForLargeOutliers = ac.GetDouble("AvgDevCutForLargeOutliers"); Double_t fVtxBegPlaneDiffCut = ac.GetDouble("VtxBegPlaneDiffCut"); Double_t fVtxPlaneDiffCut = ac.GetDouble("VtxPlaneDiffCut"); Double_t fTrkFraCut = ac.GetDouble("TrkFraCut"); Double_t fUVDiffCut = ac.GetDouble("UVDiffCut"); Double_t mip2gev = ac.GetDouble("Mip2GeV"); MSG("ShowerEM", Msg::kVerbose) << "Configured Quantities: " << "NIter = " << fIter << "\n" << "CutOff = " << fCutOff << "\n" << "ShwFrac = " << fShwFrac << "\n" << "PECut = " << fPECut << "\n" << "RadCut = " << fRadCut << "\n" << "MaxBigPlaneSep = " << fMaxBigPlaneSep << "\n" << "IsolatedPlaneCut = " << fIsolatedPlaneCut << "\n" << "NGapPlanes = " << fNGapPlanes << "\n" << "NGapStrips = " << fNGapStrips << "\n" << "MaxAvgDev = " << fMaxAvgDev << "\n" << "MaxOffFrac = " << fMaxOffFrac << "\n" << "AvgDevCutForLargeOutliers = " << fAvgDevCutForLargeOutliers << "\n" << "VtxBegPlaneDiffCut = " << fVtxBegPlaneDiffCut << "\n" << "VtxPlaneDiffCut = " << fVtxPlaneDiffCut << "\n" << "TrkFraCut = " << fTrkFraCut << "\n" << "UVDiffCut = " << fUVDiffCut << "\n" << "Mip2GeV = " << mip2gev << "\n" << endl; //important hard coded values - don't change! Double_t Ec_eff = 21.6297; Double_t X0_eff = 0.0401542; Double_t d_eff = 0.06; Double_t d = .0595; Double_t t = .0411; Int_t *plane = new Int_t[nstp]; Int_t *strip = new Int_t[nstp]; Double_t *ph = new Double_t[nstp]; Double_t *z = new Double_t[nstp]; Double_t *tpos = new Double_t[nstp]; Int_t *out = new Int_t[nstp]; for(Int_t i=0;i<nstp;i++) out[i] = 0; Float_t su[2][2] = {{0}}; Float_t sv[2][2] = {{0}}; Int_t woutstpu = -1; Int_t woutstpv = -1; Double_t totph = 0.; CandStripHandle *usefulStrip; Int_t icnt = 0; for (Int_t i=0; i<=tary->GetLast(); i++) { TObject *tobj = tary->At(i); assert(tobj->InheritsFrom("CandClusterHandle")); CandClusterHandle *clusterhandle = dynamic_cast<CandClusterHandle*>(tobj); CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator()); while (CandStripHandle *stp = stripItr()) { plane[icnt] = stp->GetPlane(); strip[icnt] = stp->GetStrip(); ph[icnt] = calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr)); z[icnt] = stp->GetZPos(); tpos[icnt] = stp->GetTPos(); if(stp->GetCharge(CalDigitType::kPE)<=fPECut) out[icnt] = 1; else totph += ph[icnt]; usefulStrip = stp; icnt++; } } UgliGeomHandle ugh(*usefulStrip->GetVldContext()); Int_t maxpl = -1; Double_t maxz = 0; Int_t secmaxpl = -1; Int_t maxstp = -1; Double_t maxtpos[2] = {0,0}; Int_t maxplu = -1; Double_t maxzu = 0; Int_t maxstpu = 0; Double_t maxtposu = 0; Int_t begstpu = 0; //Double_t begtposu = 0; Int_t begplu = -1; Double_t begzu = 0; Int_t maxplv = -1; Double_t maxzv = 0; Int_t maxstpv = 0; Double_t maxtposv = 0; Int_t begstpv = 0; //Double_t begtposv = 0; Int_t begplv = -1; Double_t begzv = 0; Float_t totwu = 0.; Float_t totwv = 0.; Float_t normposu = 0; Float_t normposv = 0; Double_t out2Rmphu = 1; Double_t out2Rmphv = 1; Double_t Thetau = 0; Double_t Theta_radu = 0; Double_t Thetav = 0; Double_t Theta_radv = 0; Double_t Thetau1 = 0; Double_t Theta_radu1 = 0; Double_t Thetav1 = 0; Double_t Theta_radv1 = 0; Double_t avgdevu = 0.; Double_t wavgdevu = 0.; Double_t avgdevv = 0.; Double_t wavgdevv = 0.; Double_t *evalu = new Double_t[2]; Double_t *evalv = new Double_t[2]; TVector2 *evu0 = new TVector2(0.,0.); TVector2 *evu1 = new TVector2(0.,0.); TVector2 *evv0 = new TVector2(0.,0.); TVector2 *evv1 = new TVector2(0.,0.); for (Int_t i = 0; i<2; i++) { evalu[i] = 0.; evalv[i] = 0.; } bool redou = false; bool redov = false; bool oncemoreu = true; bool oncemorev = true; bool remax = true; Int_t countu = 0; Int_t countv = 0; int Nu = 0; int Nv = 0; int Nu1 = 0; int Nv1 = 0; double trkfra = 0.; double uvdif = 0.; while(((out2Rmphu>fCutOff||redou||oncemoreu)&&countu<fIter) ||((out2Rmphv>fCutOff||redov||oncemorev)&&countv<fIter)){ while (remax){ remax = false; Double_t Eguess = totph*mip2gev; if(Eguess==0) {ShwStatus = 0; return;} Double_t tmax = log(Eguess*1000./Ec_eff)-0.858; //t_max formula Double_t tmaxpl = tmax*X0_eff/d_eff; //find max plane TH1F *flspl = new TH1F("flspl","flspl",500,-0.5,499.5); TH1F *sflspl = new TH1F("sflspl","sflspl",500,-0.5,499.5); for(Int_t i=0;i<nstp;i++){ if (out[i]<=0) flspl->Fill(plane[i],ph[i]); } maxpl = Int_t(flspl->GetBinCenter(flspl->GetMaximumBin())); Int_t getz = 0; for(Int_t i=0;i<nstp;i++){ if (out[i]<=0&&plane[i]!=maxpl) sflspl->Fill(plane[i],ph[i]); else if (getz<1&&out[i]<=0&&plane[i]==maxpl) { maxz = z[i]; getz++; } } secmaxpl = Int_t(sflspl->GetBinCenter(sflspl->GetMaximumBin())); if (abs(maxpl-secmaxpl)>fMaxBigPlaneSep*tmaxpl){ Bool_t pos = true; Bool_t neg = true; for(Int_t i=1;i<=fIsolatedPlaneCut;i++){ if(neg&&flspl->GetBinContent(flspl->FindBin(maxpl-i))<=0) neg=true; else neg=false; if(pos&&flspl->GetBinContent(flspl->FindBin(maxpl+i))<=0) pos=true; else pos=false; } if(pos||neg){ for(Int_t i=0;i<nstp;i++){ if (out[i]<=0&&plane[i]==maxpl) { out[i] = 4; } } remax = true; } MSG("ShowerEM", Msg::kDebug) << "Using second maximum: "<<maxpl <<" "<<secmaxpl<<" "<<tmaxpl<<endl; } delete flspl; delete sflspl; if (maxpl*secmaxpl<=0) {ShwStatus = 0; return;} } //find boundary of primary cluster: //plane: TH1F *flspl = new TH1F("flspl","flspl",500,-0.5,499.5); for(Int_t i=0;i<nstp;i++){ if (out[i]<=0) flspl->Fill(plane[i],ph[i]); } maxpl = Int_t(flspl->GetBinCenter(flspl->GetMaximumBin())); Int_t k = 0; Bool_t FoundGap = false; while(!FoundGap){ Bool_t allEmpty = true; for(Int_t i=0;i<fNGapPlanes;i++){ if(allEmpty&&flspl->GetBinContent(flspl->FindBin(maxpl-k-i))<=0) { allEmpty=true; } else allEmpty = false; } if(allEmpty) FoundGap=true; k++; } Int_t l = 0; FoundGap = false; while(!FoundGap){ Bool_t allEmpty = true; for(Int_t i=0;i<fNGapPlanes;i++){ if(allEmpty&&flspl->GetBinContent(flspl->FindBin(maxpl+l+i))<=0) { allEmpty=true; } else allEmpty = false; } if(allEmpty) FoundGap=true; l++; } for(Int_t i=0;i<nstp;i++){ if (out[i]<=0&&((plane[i]<=(maxpl-k))||(plane[i]>=(maxpl+l)))) { out[i] = 2; } } delete flspl; //Int_t npl = l+k-1; //strip: Int_t stplb[2] = {200,200}; Int_t stpub[2] = {0,0}; Int_t maxstpa[2] = {0,0}; for (Int_t j = 0; j<2; j++){ double maxph[2] = {0,0}; for (int pl = maxpl-k+1;pl<maxpl+l;pl++){ TH1F *flsstp = new TH1F("flsstp","flsstp",201,-0.5,200.5); if(j==0&&((pl<249&&pl%2==0)||(pl>=249&&pl%2==1))){ for(Int_t i=0;i<nstp;i++){ if(plane[i]==pl&&out[i]<=0) flsstp->Fill(strip[i],ph[i]); } } else if (j==1&&((pl<249&&pl%2==1)||(pl>=249&&pl%2==0))){ for(Int_t i=0;i<nstp;i++){ if(plane[i]==pl&&out[i]<=0) flsstp->Fill(strip[i],ph[i]); } } maxstp = Int_t(flsstp->GetBinCenter(flsstp->GetMaximumBin())); if((j==0&&((pl<249&&pl%2==0)||(pl>=249&&pl%2==1)))|| (j==1&&((pl<249&&pl%2==1)||(pl>=249&&pl%2==0)))){ if (maxph[j] < flsstp->GetBinContent(flsstp->GetMaximumBin())){ maxph[j] = flsstp->GetBinContent(flsstp->GetMaximumBin()); maxstpa[j] = maxstp; } } Int_t m[2] = {0}; FoundGap = false; while(!FoundGap){ Bool_t allEmpty = true; for(Int_t i=0;i<fNGapStrips;i++){ if(allEmpty&&flsstp->GetBinContent(flsstp->FindBin(maxstp-m[j]-i))<=0){ allEmpty=true; } else allEmpty = false; } if(allEmpty) FoundGap=true; m[j]++; } Int_t n[2] = {0}; FoundGap = false; while(!FoundGap){ Bool_t allEmpty = true; for(Int_t i=0;i<fNGapStrips;i++){ if(allEmpty&&flsstp->GetBinContent(flsstp->FindBin(maxstp+n[j]+i))<=0){ allEmpty=true; } else allEmpty = false; } if(allEmpty) FoundGap=true; n[j]++; } for(Int_t i=0;i<nstp;i++){ if((j==0&&((pl<249&&pl%2==0)||(pl>=249&&pl%2==1)))|| (j==1&&((pl<249&&pl%2==1)||(pl>=249&&pl%2==0)))){ if (plane[i]==pl&&out[i]<=0 &&((strip[i]<=(maxstp-m[j]))||(strip[i]>=(maxstp+n[j])))) { out[i] = 1; } } } if (maxstp>0&&(maxstp-m[j])<stplb[j]) stplb[j] = maxstp-m[j]; if (maxstp>0&&(maxstp+n[j])>stpub[j]) stpub[j] = maxstp+n[j]; delete flsstp; } } int ctpl = 0; int ctstp[2] = {0,0}; for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if(plane[i]==maxpl&&ctpl<1){ maxz = z[i]; ctpl++; } for (int j=0; j<2; j++){ if((j==0&&((plane[i]<249&&plane[i]%2==0) ||(plane[i]>=249&&plane[i]%2==1)))|| (j==1&&((plane[i]<249&&plane[i]%2==1) ||(plane[i]>=249&&plane[i]%2==0)))){ if(strip[i]==maxstpa[j]&&ctstp[j]<1){ maxtpos[j] = tpos[i]; ctstp[j]++; } } } } } // MSG("ShowerEM", Msg::kDebug) <<stplb[0]<<" "<<stplb[1]<<" "<<stpub[0]<<" "<<stpub[1]<<" "<<maxz<<" "<<maxtpos[0]<<" "<<maxtpos[1]<<endl; //iterate over U strips while((out2Rmphu>fCutOff||redou||oncemoreu)&&countu<fIter){ redou = false; oncemoreu = false; totwu = 0; countu++; if(woutstpu!=-1&&out[woutstpu]==0) out[woutstpu] = 1; else if (out[woutstpu]<0) out[woutstpu] = 5; TH1F *flsplu = new TH1F("flsplu","flsplu",500,-0.5,499.5); TH2F *u2d = new TH2F("u2d","u2d",l+k+1,maxz-(k+0.5)*d,maxz+(l+0.5)*d, stpub[0]-stplb[0]+1,maxtpos[0]-(maxstpa[0]-stplb[0]+0.5)*t, maxtpos[0]+(stpub[0]-maxstpa[0]+0.5)*t); for(Int_t i=0;i<nstp;i++){ if(((plane[i]<249&&plane[i]%2==0) ||(plane[i]>=249&&plane[i]%2==1))&&out[i]<=0){ flsplu->Fill(plane[i],ph[i]); u2d->Fill(z[i],tpos[i],ph[i]); } } maxplu = Int_t(flsplu->GetBinCenter(flsplu->GetMaximumBin())); Int_t nplu = 0; for (Int_t pl = maxpl-k+1;pl<maxpl+l;pl++){ if (flsplu->GetBinContent(flsplu->FindBin(pl))>0) nplu++; } TProfile *up = u2d->ProfileX(); TF1 *p1 = new TF1("p1","pol1"); up->Fit("p1","NQ"); Double_t slopeu = p1->GetParameter(1); //Double_t eslopeu = p1->GetParError(1); Double_t chi2u = p1->GetChisquare(); Int_t j = 0; Double_t sumb = 0.; while(sumb==0&&j<=500) { j++; if(flsplu->GetBinContent(j)>0) { begplu = Int_t(flsplu->GetBinCenter(j)); sumb += flsplu->GetBinContent(j); } } delete flsplu; delete u2d; delete up; delete p1; Float_t maxstpphu = 0.; Float_t begstpphu = 0.; Int_t nstpmaxplu = 0; for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if(plane[i]==maxplu&&ph[i]>maxstpphu){ maxstpu = strip[i]; maxstpphu = ph[i]; maxtposu = tpos[i]; } if(plane[i]==maxplu){ maxzu = z[i]; nstpmaxplu++; } if(plane[i]==begplu&&ph[i]>begstpphu){ begzu = z[i]; begstpu = strip[i]; begstpphu = ph[i]; } } } Double_t cntu = 0; for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if((plane[i]<249&&plane[i]%2==0)||(plane[i]>=249&&plane[i]%2==1)){ Double_t dtpos = tpos[i] - maxtposu; Double_t dz = z[i] - maxzu; su[0][0]+=dz*dz*ph[i]; su[0][1]+=dtpos*dz*ph[i]; su[1][1]+=dtpos*dtpos*ph[i]; totwu += ph[i]; cntu+=1; normposu += dtpos*dtpos+dz*dz; } } } su[1][0]=su[0][1]; if(totwu*normposu!=0){ TMatrixDSym mu(2); Double_t *angu = new Double_t[2]; for(Int_t i=0;i<2;i++) angu[i] = -999; for (Int_t i=0;i<2;i++){ for (Int_t j=0;j<2;j++){ mu(i,j)=su[i][j]/totwu/normposu; } } if(mu.IsA()){ TMatrixDSymEigen eigen(mu); TMatrixD EVect1 = eigen.GetEigenVectors(); TVectorD EVal1 = eigen.GetEigenValues(); evu0 = new TVector2(EVect1(0,0),EVect1(1,0)); evu1 = new TVector2(EVect1(0,1),EVect1(1,1)); angu[0] = evu0->Phi()*180/TMath::Pi(); angu[1] = evu1->Phi()*180/TMath::Pi(); evalu[0] = EVal1(0); evalu[1] = EVal1(1); } else { evalu[0] = 0.; evalu[1] = 0.; } for (Int_t i=0;i<2;i++){ if (angu[i]>90&&angu[i]<270) angu[i]-=180; if (angu[i]>270&&angu[i]<360) angu[i]-=360; if (angu[i]==90||angu[i]==270) angu[i]-=0.01; } Thetau = angu[0]; Thetau1 = angu[1]; Theta_radu = (Thetau/180)*TMath::Pi(); Theta_radu1 = (Thetau1/180)*TMath::Pi(); Double_t *shwstpc = new Double_t[nstp]; Double_t *shwstpc1 = new Double_t[nstp]; Double_t *shwstpc2 = new Double_t[nstp]; Double_t *wshwstpc = new Double_t[nstp]; for(Int_t i=0;i<nstp;i++) shwstpc[i] = -1.; for(Int_t i=0;i<nstp;i++) shwstpc1[i] = -1.; for(Int_t i=0;i<nstp;i++) shwstpc2[i] = -1.; for(Int_t i=0;i<nstp;i++) wshwstpc[i] = -1.; Double_t wdevmaxu = 0.; woutstpu = -1; avgdevu = 0.; wavgdevu = 0.; Double_t avgdevu1 = 0.; Double_t avgdevu2 = 0.; Double_t wavgdevu1 = 0.; Double_t wavgdevu2 = 0.; Int_t instpu = 0; //Int_t instpu1 = 0; Double_t winstpu = 0; Float_t offstpu = 0; Float_t offstpu2 = 0; //Float_t offstpu3 = 0; out2Rmphu = 0; for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if((plane[i]<249&&plane[i]%2==0)||(plane[i]>=249&&plane[i]%2==1)){ shwstpc[i]=fabs(tpos[i] -(maxtposu+tan(Theta_radu) *(z[i]-maxzu)))*TMath::Cos(Theta_radu); shwstpc2[i]=fabs(tpos[i] -(maxtposu+tan(Theta_radu1) *(z[i]-maxzu)))*TMath::Cos(Theta_radu1); avgdevu1 += shwstpc[i]; avgdevu2 += shwstpc2[i]; wavgdevu1 += shwstpc[i]*ph[i]; wavgdevu2 += shwstpc2[i]*ph[i]; instpu++; winstpu+=ph[i]; if (shwstpc[i]>1.*t) offstpu++; if (shwstpc2[i]>1.*t) offstpu2++; // if (shwstpc2[i]>6.*t) offstpu3++; } } } avgdevu1/=instpu; avgdevu2/=instpu; wavgdevu1/=winstpu; wavgdevu2/=winstpu; avgdevu = avgdevu1; wavgdevu = wavgdevu1; offstpu/=instpu; offstpu2/=instpu; for(Int_t i=0;i<nstp;i++){ if((fabs(Thetau)>fabs(Thetau1)) &&((avgdevu1>fMaxAvgDev*t&&avgdevu2>fMaxAvgDev*t) ||(offstpu>fMaxOffFrac&&offstpu2>fMaxOffFrac)) &&strip[i]==maxstpu&&nstpmaxplu==1) { out[i] = 4; oncemoreu = true; MSG("ShowerEM", Msg::kDebug) <<"u maxstp removed"<<endl; } } if(fabs(tan(Theta_radu)-slopeu)>fabs(tan(Theta_radu1)-slopeu) &&(wavgdevu1>wavgdevu2/1000.||nplu==1)){ Thetau = angu[1]; Thetau1 = angu[0]; avgdevu = avgdevu2; wavgdevu = wavgdevu2; for(Int_t i=0;i<nstp;i++){ shwstpc[i] = shwstpc2[i]; } } Theta_radu = (Thetau/180)*TMath::Pi(); Theta_radu1 = (Thetau1/180)*TMath::Pi(); for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if((plane[i]<249&&plane[i]%2==0)||(plane[i]>=249&&plane[i]%2==1)){ wshwstpc[i] = shwstpc[i]*ph[i]; if(shwstpc[i]>fRadCut*t) { out2Rmphu += ph[i]; if(wshwstpc[i]>wdevmaxu){ wdevmaxu = wshwstpc[i]; woutstpu = i; } } } } else if(out[i]==1){ if((plane[i]<249&&plane[i]%2==0)||(plane[i]>=249&&plane[i]%2==1)){ shwstpc[i] = fabs(tpos[i]-(maxtposu+tan(Theta_radu)*(z[i]-maxzu)))*TMath::Cos(Theta_radu);; wshwstpc[i] = shwstpc[i]*ph[i]; if(shwstpc[i]<=fRadCut*t) { out[i] = -1; redou = true; } } } } MSG("ShowerEM", Msg::kDebug) <<"avgu "<< avgdevu <<" "<< avgdevu1 <<" "<<avgdevu2<<" "<<wavgdevu1 <<" " <<wavgdevu2<<" "<<offstpu<<" " <<offstpu2<<" "<<tan(Theta_radu)<<" " <<slopeu<<" "<<chi2u<<endl; if (avgdevu>fAvgDevCutForLargeOutliers*t){ for(Int_t i=0;i<nstp;i++){ if(((plane[i]<249&&plane[i]%2==0) ||(plane[i]>=249&&plane[i]%2==1))&&out[i]<=0){ if((shwstpc[i]-avgdevu)>fRadCut*t) { if (out[i]==0) out[i] = 1; else if (out[i]<0) out[i] = 5; oncemoreu = true; } } } } out2Rmphu /= totwu; Nu = 0; Nu1 = 0; for (int pl = maxpl-k+1;pl<maxpl+l;pl++){ int nstppl = 0; for(int i=0;i<nstp;i++){ if(plane[i]==pl&&((plane[i]<249&&plane[i]%2==0) ||(plane[i]>=249&&plane[i]%2==1)) &&out[i]<=0){ nstppl++; } } if (nstppl>0) Nu++; if (nstppl==1) Nu1++; } } } //iterate over V strips while((out2Rmphv>fCutOff||redov||oncemorev)&&countv<fIter){ redov = false; oncemorev = false; totwv = 0; countv++; if(woutstpv!=-1&&out[woutstpv]==0) out[woutstpv] = 1; else if (out[woutstpv]<0) out[woutstpv] = 5; TH1F *flsplv = new TH1F("flsplv","flsplv",500,-0.5,499.5); TH2F *v2d = new TH2F("v2d","v2d",l+k+1,maxz-(k+0.5)*d,maxz+(l+0.5)*d, stpub[1]-stplb[1]+1,maxtpos[1]-(maxstpa[1]-stplb[1]+0.5)*t, maxtpos[1]+(stpub[1]-maxstpa[1]+0.5)*t); for(Int_t i=0;i<nstp;i++){ if(((plane[i]>=249&&plane[i]%2==0) ||(plane[i]<249&&plane[i]%2==1))&&out[i]<=0){ flsplv->Fill(plane[i],ph[i]); v2d->Fill(z[i],tpos[i],ph[i]); } } maxplv = Int_t(flsplv->GetBinCenter(flsplv->GetMaximumBin())); Int_t nplv = 0; for (Int_t pl = maxpl-k+1;pl<maxpl+l;pl++){ if (flsplv->GetBinContent(flsplv->FindBin(pl))>0) nplv++; } TProfile *vp = v2d->ProfileX(); TF1 *p1 = new TF1("p1","pol1"); vp->Fit("p1","NQ"); Double_t slopev = p1->GetParameter(1); //Double_t eslopev = p1->GetParError(1); Double_t chi2v = p1->GetChisquare(); Int_t j = 0; Double_t sumb = 0.; while (sumb==0&&j<=500) { j++; if(flsplv->GetBinContent(j)>0) { begplv = Int_t(flsplv->GetBinCenter(j)); sumb += flsplv->GetBinContent(j); } } delete flsplv; delete v2d; delete vp; delete p1; Float_t maxstpphv = 0.; Float_t begstpphv = 0.; Int_t nstpmaxplv = 0; for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if(plane[i]==maxplv&&ph[i]>maxstpphv){ maxstpv = strip[i]; maxstpphv = ph[i]; maxtposv = tpos[i]; } if(plane[i]==maxplv){ maxzv = z[i]; nstpmaxplv++; } if (plane[i]==begplv&&ph[i]>begstpphv){ begzv = z[i]; begstpv = strip[i]; begstpphv = ph[i]; } } } Double_t cntv = 0; for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if ((plane[i]>=249&&plane[i]%2==0)||(plane[i]<249&&plane[i]%2==1)){ Double_t dtpos = tpos[i] - maxtposv; Double_t dz = z[i] - maxzv; sv[0][0]+=dz*dz*ph[i]; sv[0][1]+=dtpos*dz*ph[i]; sv[1][1]+=dtpos*dtpos*ph[i]; totwv += ph[i]; cntv+=1; normposv += dtpos*dtpos+dz*dz; } } } sv[1][0]=sv[0][1]; if(totwv*normposv!=0){ TMatrixDSym mv(2); Double_t *angv = new Double_t[2]; for(Int_t i=0;i<2;i++) angv[i] = -999; for (Int_t i=0;i<2;i++){ for (Int_t j=0;j<2;j++){ mv(i,j)=sv[i][j]/totwv/normposv; } } if(mv.IsA()){ TMatrixDSymEigen eigen(mv); TMatrixD EVect2 = eigen.GetEigenVectors(); TVectorD EVal2 = eigen.GetEigenValues(); evv0 = new TVector2(EVect2(0,0),EVect2(1,0)); evv1 = new TVector2(EVect2(0,1),EVect2(1,1)); angv[0] = evv0->Phi()*180/TMath::Pi(); angv[1] = evv1->Phi()*180/TMath::Pi(); evalv[0] = EVal2(0); evalv[1] = EVal2(1); } else { evalv[0] = 0.; evalv[1] = 0.; } for (Int_t i=0;i<2;i++){ if (angv[i]>90&&angv[i]<270) angv[i]-=180; if (angv[i]>270&&angv[i]<360) angv[i]-=360; if (angv[i]==90||angv[i]==270) angv[i]-=0.01; } Thetav = angv[0]; Thetav1 = angv[1]; Theta_radv = (Thetav/180)*TMath::Pi(); Theta_radv1 = (Thetav1/180)*TMath::Pi(); Double_t *shwstpc = new Double_t[nstp]; Double_t *shwstpc1 = new Double_t[nstp]; Double_t *shwstpc2 = new Double_t[nstp]; Double_t *wshwstpc = new Double_t[nstp]; for(Int_t i=0;i<nstp;i++) shwstpc[i] = -1.; for(Int_t i=0;i<nstp;i++) shwstpc1[i] = -1.; for(Int_t i=0;i<nstp;i++) shwstpc2[i] = -1.; for(Int_t i=0;i<nstp;i++) wshwstpc[i] = -1.; Double_t wdevmaxv = 0.; woutstpv = -1; avgdevv = 0.; wavgdevv = 0.; Double_t avgdevv1 = 0.; Double_t avgdevv2 = 0.; Double_t wavgdevv1 = 0.; Double_t wavgdevv2 = 0.; Int_t instpv = 0; //Int_t instpv1 = 0; Double_t winstpv = 0; Float_t offstpv = 0; Float_t offstpv2 = 0; // Float_t offstpv3 = 0; out2Rmphv = 0; for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if((plane[i]<249&&plane[i]%2==1)||(plane[i]>=249&&plane[i]%2==0)){ shwstpc[i]=fabs(tpos[i] -(maxtposv+tan(Theta_radv) *(z[i]-maxzv)))*TMath::Cos(Theta_radv); shwstpc2[i]=fabs(tpos[i] -(maxtposv+tan(Theta_radv1) *(z[i]-maxzv)))*TMath::Cos(Theta_radv1); avgdevv1 += shwstpc[i]; avgdevv2 += shwstpc2[i]; wavgdevv1 += shwstpc[i]*ph[i]; wavgdevv2 += shwstpc2[i]*ph[i]; instpv++; winstpv+=ph[i]; if (shwstpc[i]>1) offstpv++; if (shwstpc2[i]>1) offstpv2++; // if (shwstpc2[i]>6) offstpv3++; } } } avgdevv1/=instpv; avgdevv2/=instpv; wavgdevv1/=winstpv; wavgdevv2/=winstpv; avgdevv = avgdevv1; wavgdevv = wavgdevv1; offstpv/=instpv; offstpv2/=instpv; for(Int_t i=0;i<nstp;i++){ if ((fabs(Thetav)>fabs(Thetav1)) &&((avgdevv1>fMaxAvgDev*t&&avgdevv2>fMaxAvgDev*t) ||(offstpv>fMaxOffFrac&&offstpv2>fMaxOffFrac)) &&strip[i]==maxstpv&&nstpmaxplv==1) { out[i] = 4; oncemorev = true; MSG("ShowerEM", Msg::kDebug) <<"v maxstp removed"<<endl; } } if(fabs(tan(Theta_radv)-slopev)>fabs(tan(Theta_radv1)-slopev) &&(wavgdevv1>wavgdevv2/1000.||nplv==1)){ Thetav = angv[1]; Thetav1 = angv[0]; avgdevv = avgdevv2; wavgdevv = wavgdevv2; for(Int_t i=0;i<nstp;i++){ shwstpc[i] = shwstpc2[i]; } } Theta_radv = (Thetav/180)*TMath::Pi(); Theta_radv1 = (Thetav1/180)*TMath::Pi(); for(Int_t i=0;i<nstp;i++){ if(out[i]<=0){ if((plane[i]<249&&plane[i]%2==1)||(plane[i]>=249&&plane[i]%2==0)){ wshwstpc[i] = shwstpc[i]*ph[i]; if(shwstpc[i]>fRadCut*t) { out2Rmphv += ph[i]; if(wshwstpc[i]>wdevmaxv){ wdevmaxv = wshwstpc[i]; woutstpv = i; } } } } else if(out[i]==1){ if((plane[i]<249&&plane[i]%2==1)||(plane[i]>=249&&plane[i]%2==0)){ shwstpc[i] = fabs(tpos[i]-(maxtposv+tan(Theta_radv)*(z[i]-maxzv)))*TMath::Cos(Theta_radv);; wshwstpc[i] = shwstpc[i]*ph[i]; if(shwstpc[i]<=fRadCut*t) { out[i] = -1; redov = true; } } } } MSG("ShowerEM", Msg::kDebug) <<"avgv "<< avgdevv <<" "<<avgdevv1 <<" "<<avgdevv2<<" "<<wavgdevv1 <<" " <<wavgdevv2<<" "<<offstpv<<" " <<offstpv2<<" "<<tan(Theta_radv)<<" " <<slopev<<" "<<chi2v<<endl; if (avgdevv>fAvgDevCutForLargeOutliers*t){ for(Int_t i=0;i<nstp;i++){ if(((plane[i]<249&&plane[i]%2==1) ||(plane[i]>=249&&plane[i]%2==0))&&out[i]<=0){ if((shwstpc[i]-avgdevv)>fRadCut*t) { if (out[i]==0) out[i] = 1; else if (out[i]<0) out[i] = 5; oncemorev = true; } } } } out2Rmphv /= totwv; Nv = 0; Nv1 = 0; for (int pl = maxpl-k+1;pl<maxpl+l;pl++){ int nstppl = 0; for(int i=0;i<nstp;i++){ if(plane[i]==pl&&((plane[i]<249&&plane[i]%2==1) ||(plane[i]>=249&&plane[i]%2==0)) &&out[i]<=0){ nstppl++; } } if (nstppl>0) Nv++; if (nstppl==1) Nv1++; } } } if ((Nu+Nv)>0){ trkfra = double(Nu1+Nv1)/double(Nu+Nv); uvdif = double(abs(Nu-Nv))/double(Nu+Nv); } else { trkfra = -1; uvdif = -1; } Double_t shwfra = (totwu+totwv)/totph; if(totwu*totwv==0||shwfra<fShwFrac||(countu==fIter&&out2Rmphu>fCutOff) ||(countv==fIter&&out2Rmphv>fCutOff)) ShwStatus = 0; if(ShwStatus==0) return; if (trkfra>fTrkFraCut) ShwStatus = 2; if (uvdif>fUVDiffCut) ShwStatus = 3; Double_t Eguess = (totwu+totwv)*mip2gev; if(Eguess==0) {ShwStatus = 0; return;} Double_t tmax = log(Eguess*1000./Ec_eff)-0.858; Double_t tmaxz = tmax*X0_eff; Double_t begz = begzv; Double_t begpl = begplv; if(begzu<begzv) { begz = begzu; begpl = begplu; } Double_t vtxz = maxz-tmaxz/sqrt(TMath::Power(tan(Theta_radu),2)+TMath::Power(tan(Theta_radv),2)+1.); MSG("ShowerEM", Msg::kDebug) <<shwfra<<" "<<Eguess<<" "<<tmaxz<<" " <<vtxz<<" "<<begz<<" "<<maxz<<endl; if ((vtxz-begz)<=fVtxBegPlaneDiffCut*d) { fVertex[2] = begz; fVtxPlane = begpl; } else { fVertex[2] = vtxz; fVtxPlane = begpl+(vtxz-begz)/d; for(Int_t i=0;i<nstp;i++){ if(((z[i]-vtxz)<-fVtxPlaneDiffCut*d) &&((plane[i]<249&&plane[i]%2==0) ||(plane[i]>=249&&plane[i]%2==1)) &&out[i]<=0) { out[i] = 3; oncemoreu = true; } else if(((z[i]-vtxz)<-fVtxPlaneDiffCut*d) &&((plane[i]<249&&plane[i]%2==1) ||(plane[i]>=249&&plane[i]%2==0)) &&out[i]<=0) { out[i] = 3; oncemorev = true; } } } } fVertex[0] = maxtposu+tan(Theta_radu)*(fVertex[2]-maxzu); fVertex[1] = maxtposv+tan(Theta_radv)*(fVertex[2]-maxzv); fAngle[0] = Theta_radu; fAngle[1] = Theta_radv; fAngle[2] = Theta_radu1; fAngle[3] = Theta_radv1; fEigenValue[0] = evalu[0]; fEigenValue[1] = evalv[0]; fEigenValue[2] = evalu[1]; fEigenValue[3] = evalv[1]; fAvgDev[0] = avgdevu; fAvgDev[1] = avgdevv; fAvgDev[2] = wavgdevu; fAvgDev[3] = wavgdevv; fEigenVector[0] = evu0->X(); fEigenVector[1] = evu0->Y(); fEigenVector[2] = evv0->X(); fEigenVector[3] = evv0->Y(); fEigenVector[4] = evu1->X(); fEigenVector[5] = evu1->Y(); fEigenVector[6] = evv1->X(); fEigenVector[7] = evv1->Y(); icnt = 0; Int_t jcnt = 0; for (Int_t i=0; i<=tary->GetLast(); i++) { TObject *tobj = tary->At(i); assert(tobj->InheritsFrom("CandClusterHandle")); CandClusterHandle *clusterhandle = dynamic_cast<CandClusterHandle*>(tobj); CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator()); while (CandStripHandle *stp = stripItr()) { if(out[icnt]<=0){ //get min/max number of strips if(stp->GetStrip()>fMaxStrip) fMaxStrip = stp->GetStrip(); if(stp->GetStrip()<fMinStrip) fMinStrip = stp->GetStrip(); if(stp->GetPlane()>fMaxPlane) fMaxPlane = stp->GetPlane(); //find minimum strip time at vertex plane PlexStripEndId stripid = stp->GetStripEndId(); if(stp->GetPlane()==int(fVtxPlane+0.5)){ if(stp->GetTime()<fMinTime) fMinTime = stp->GetTime(); } csh.AddDaughterLink(*stp); jcnt++; } else{ fOutPH[out[icnt]-1] += calibrator. GetMIP(stp->GetCharge(CalDigitType::kSigCorr)); } icnt++; } csh.AddCluster(clusterhandle); } csh.SetVtxPlane(int(fVtxPlane+0.5)); csh.SetVtxZ(fVertex[2]); csh.SetVtxU(fVertex[0]); csh.SetVtxV(fVertex[1]); csh.SetVtxT(fMinTime); double dudz = TMath::Tan(fAngle[0]); double dvdz = TMath::Tan(fAngle[1]); double dsdz = sqrt(dudz*dudz + dvdz*dvdz + 1.); csh.SetDirCosU(dudz/dsdz); csh.SetDirCosV(dvdz/dsdz); csh.SetDirCosZ(1./dsdz); csh.SetEigenVectors(fEigenVector); csh.SetEigenValues(fEigenValue); csh.SetAvgDev(fAvgDev); csh.SetOutPH(fOutPH); csh.SetShwStatus(ShwStatus); //SetUV(&csh); Calibrate(&csh); Double_t TotMip = csh.GetCharge(CalStripType::kMIP); Double_t TotEnergy = csh.GetCharge(CalStripType::kGeV); MSG("ShowerEM", Msg::kVerbose) << endl << " TotMIP = " << TotMip << endl << "Setting shower energy to " << TotEnergy << endl; csh.SetEnergy(TotEnergy); }

void AlgShowerEM::Trace (  const char *  c  )  const [virtual]

Reimplemented from AlgBase. Definition at line 1111 of file AlgShowerEM.cxx. { }

---

The documentation for this class was generated from the following files:

• AlgShowerEM.h
• AlgShowerEM.cxx

---