AlgSubShowerSR Class Reference

#include <AlgSubShowerSR.h>

Inheritance diagram for AlgSubShowerSR:
AlgBase

List of all members.

Public Member Functions

 AlgSubShowerSR ()
virtual ~AlgSubShowerSR ()
virtual void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void Trace (const char *c) const

Private Member Functions

void CalculateEnergyVertexAngle (CandSubShowerSRHandle &cch, Int_t det)
void SubShowerID (CandSubShowerSRHandle &cch)

Private Attributes

Double_t fPECut
Double_t fECut
Double_t fMaxDiffCut
Double_t fVtxDiffCut
Double_t fOut2RmPHCut
Double_t fTrkFraCut
Double_t fTrkLikeStpCut
Double_t fXtkFraCut
Double_t fMaxXTalkStpPH
Bool_t fDoEMFit
Double_t fMIPperGeV
Double_t maxStpPHinSS
Double_t maxdiff
Double_t vtxdiff
Double_t trkfra
Double_t avgstpperpln
Double_t xtkfra
Double_t avgph
Double_t out2Rmph
Double_t bothfitprob
Double_t stpdismean
Double_t stpdisrms
std::list< Double_t > LongPl
std::list< Double_t > LongPh
std::list< Double_t > LongFlu
std::list< Double_t > LongPhpe
std::list< Double_t > TranStp
std::list< Double_t > TranPh
std::list< Double_t > TranFlu
Double_t loweststp
Double_t uppeststp
Double_t begz
Double_t endz
Double_t Theta_rad
Double_t totw
Double_t Eguess
Int_t nn0l
Int_t nn0t
TH1F * flspl
TH1F * sflspl
TF1 * clusterlong
TF1 * clustertran
TF1 * chi2func
TH1F * LongProf
TH1F * LongHist
TH1F * TranProf
TH1F * TranHist
TH1F * LongProfb
TH1F * LongHistb
TH1F * TranProfb
TH1F * TranHistb

Detailed Description

Definition at line 20 of file AlgSubShowerSR.h.


Constructor & Destructor Documentation

AlgSubShowerSR::AlgSubShowerSR (  ) 

Definition at line 54 of file AlgSubShowerSR.cxx.

References Chi2(), LongShwProfSamAng(), and TranShwProfSamAng().

00054                                :
00055   fPECut(1.5),fECut(0.1),fMaxDiffCut(1),fVtxDiffCut(2),
00056   fOut2RmPHCut(0.3),fTrkFraCut(0.8),fTrkLikeStpCut(1.3),
00057   fXtkFraCut(0.8),fMaxXTalkStpPH(2.0),fDoEMFit(1),fMIPperGeV(18.5)
00058 {
00059   flspl = new TH1F("flspl","flspl",500,-0.5,499.5); 
00060   sflspl = new TH1F("sflspl","sflspl",500,-0.5,499.5);
00061   clusterlong = new TF1("clusterlong",LongShwProfSamAng,0,1,2);
00062   clustertran = new TF1("clustertran",TranShwProfSamAng,0,1,2);
00063   chi2func = new TF1("chi2func",Chi2,0,5000,1);
00064   
00065   LongProf = new TH1F("longprof","Long Prof",1,0,1);
00066   LongHist = new TH1F("longhist","Long Hist",1,0,1);
00067   TranProf = new TH1F("tranprof","Tran Prof",1,0,1);
00068   TranHist = new TH1F("tranhist","Tran Hist",1,0,1);
00069   LongProfb = new TH1F("longprofb","Long Prof Both",1,0,1);
00070   LongHistb = new TH1F("longhistb","Long Hist Both",1,0,1);
00071   TranProfb = new TH1F("tranprofb","Tran Prof Both",1,0,1);
00072   TranHistb = new TH1F("tranhistb","Tran Hist Both",1,0,1);
00073 }

AlgSubShowerSR::~AlgSubShowerSR (  )  [virtual]

Definition at line 76 of file AlgSubShowerSR.cxx.

References chi2func, clusterlong, clustertran, flspl, LongFlu, LongHist, LongHistb, LongPh, LongPhpe, LongPl, LongProf, LongProfb, sflspl, TranFlu, TranHist, TranHistb, TranPh, TranProf, TranProfb, and TranStp.

00077 {
00078 
00079   delete flspl;
00080   delete sflspl;
00081   delete clusterlong;
00082   delete clustertran;
00083   delete chi2func;
00084 
00085   delete LongProf;
00086   delete LongHist;
00087   delete TranProf;
00088   delete TranHist;
00089   delete LongProfb;
00090   delete LongHistb;
00091   delete TranProfb;
00092   delete TranHistb;
00093 
00094   LongPl.clear();
00095   LongPh.clear();
00096   LongFlu.clear();
00097   LongPhpe.clear();
00098   TranStp.clear();
00099   TranPh.clear();
00100   TranFlu.clear();
00101 
00102 }


Member Function Documentation

void AlgSubShowerSR::CalculateEnergyVertexAngle ( CandSubShowerSRHandle cch,
Int_t  det 
) [private]

Definition at line 150 of file AlgSubShowerSR.cxx.

References avgph, avgstpperpln, begz, EMFluctuation::CalcLongFluc(), EMFluctuation::CalcTranFluc(), D_EFF, EC_EFF, Eguess, endz, flspl, fMIPperGeV, fPECut, CandHandle::GetDaughterIterator(), Calibrator::GetMIP(), CandRecoHandle::GetNStrip(), CandSubShowerSRHandle::GetPlaneView(), Calibrator::Instance(), Msg::kDebug, Msg::kError, CalDigitType::kPE, CalDigitType::kSigCorr, PlaneView::kU, PlaneView::kV, Msg::kVerbose, PlaneView::kX, PlaneView::kY, LongFlu, LongPh, LongPhpe, LongPl, loweststp, Munits::m, maxdiff, maxStpPHinSS, Munits::ms, MSG, n, nn0l, nn0t, out2Rmph, RM_EFF, CandSubShowerSRHandle::SetAvgDev(), CandRecoHandle::SetEndPlane(), CandRecoHandle::SetEndZ(), CandSubShowerSRHandle::SetEnergy(), CandSubShowerSRHandle::SetSlope(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), sflspl, ST_WID, T_EFF, Theta_rad, totw, TranFlu, TranPh, TranStp, trkfra, uppeststp, vtxdiff, X0_EFF, and xtkfra.

Referenced by RunAlg().

00152 {
00153 
00154   Int_t nstp      = cch.GetNStrip();
00155   Int_t begpl     = 500;
00156   Int_t endpl     = -1;
00157   CandStripHandleItr stripItr(cch.GetDaughterIterator());  
00158   while (CandStripHandle *stp = stripItr()){
00159     Int_t pln = stp->GetPlane();
00160     if(pln<begpl) begpl = pln;
00161     if(pln>endpl) endpl = pln;
00162   }
00163 
00164   std::vector<Int_t> plane(nstp,0);
00165   std::vector<Int_t> strip(nstp,0);
00166   std::vector<Double_t> ph(nstp,0);
00167   std::vector<Double_t> ph_pe(nstp,0);
00168   std::vector<Double_t> z(nstp,0);
00169   std::vector<Double_t> tpos(nstp,0);
00170   std::vector<Double_t> time(nstp,0);
00171 
00172   totw            = 0.;
00173   Double_t totph  = 0.;
00174   Int_t Ns0       = 0;
00175   Int_t icnt      = 0;
00176 
00177   maxStpPHinSS    = 0;
00178   loweststp       = 200;
00179   uppeststp       = -200;
00180   Double_t lowestpl = begpl;
00181   Double_t uppestpl = endpl;
00182 
00183   Calibrator& calibrator=Calibrator::Instance();
00184 
00185   stripItr.Reset();
00186   while (CandStripHandle *stp = stripItr()) {
00187     
00188     plane[icnt] = stp->GetPlane();
00189     strip[icnt] = stp->GetStrip();
00190     z[icnt] = stp->GetZPos();
00191     tpos[icnt] = stp->GetTPos();
00192     time[icnt] = stp->GetTime();
00193 
00194     ph[icnt] = calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
00195     totw += ph[icnt];    
00196     
00197     ph_pe[icnt] = stp->GetCharge(CalDigitType::kPE);
00198     totph += ph_pe[icnt];
00199 
00200     if(tpos[icnt]/T_EFF*100.>uppeststp) {
00201       uppeststp = tpos[icnt]/T_EFF*100.;
00202       uppestpl = plane[icnt];
00203     }
00204     if(tpos[icnt]/T_EFF*100.<loweststp) {
00205       loweststp = tpos[icnt]/T_EFF*100.;
00206       lowestpl = plane[icnt];
00207     }
00208     if(ph[icnt]>maxStpPHinSS) maxStpPHinSS = ph[icnt];
00209     if (ph_pe[icnt]<=fPECut) Ns0++;
00210     icnt++;
00211   }
00212 
00213   maxdiff      = 0.;
00214   vtxdiff      = 0.;
00215   trkfra       = 0;
00216   avgstpperpln = 0;
00217   xtkfra       = 0;
00218   avgph        = 0;  
00219   out2Rmph     = 0;
00220 
00221   Int_t maxpl      = -1;
00222   Double_t maxz    = 0.;
00223   begz             = 0.;
00224   endz             = 30.;
00225   Int_t maxstp     = -1;
00226   Double_t maxtpos = 0.;
00227   Int_t begstp     = -1;
00228   Double_t mint    = 999999999;
00229   Float_t normpos  = 0;
00230   
00231   Double_t Theta      = 0;
00232   Theta_rad           = 0;
00233   Double_t Theta1     = 0;
00234   Double_t Theta_rad1 = 0;
00235   Double_t avgdev     = 0.;
00236   Double_t wavgdev    = 0.;
00237   Double_t avgdev1    = 0.;
00238   Double_t wavgdev1   = 0.;
00239   Double_t avgdev2    = 0.;
00240   Double_t wavgdev2   = 0.;
00241   Double_t eval[2]    = {0};
00242 
00243   Eguess = totw*2./fMIPperGeV;
00244   if(Eguess==0) {
00245     MSG("SubShowerSR",Msg::kError) << "Eguess = 0" << endl; 
00246     return;
00247   }
00248 
00249   Double_t tmax = TMath::Log(Eguess*1000./EC_EFF)-0.858;
00250   Double_t tmaxpl = tmax*X0_EFF/D_EFF;
00251   flspl->Reset();
00252   sflspl->Reset();
00253 
00254   for(Int_t i=0;i<nstp;i++){
00255     flspl->Fill(plane[i],ph[i]);
00256   }
00257   maxpl = Int_t(flspl->GetBinCenter(flspl->GetMaximumBin()));
00258   Int_t getz = 0;
00259   for(Int_t i=0;i<nstp;i++){
00260     if (plane[i]!=maxpl) sflspl->Fill(plane[i],ph[i]);
00261     else if (getz<1&&plane[i]==maxpl) {
00262       maxz = z[i];
00263       getz++;
00264     }
00265   }
00266   double secmaxpl = Int_t(sflspl->GetBinCenter(sflspl->GetMaximumBin()));
00267   maxdiff = TMath::Abs(maxpl-secmaxpl)/tmaxpl;
00268   vtxdiff = TMath::Abs(maxpl-begpl)/tmaxpl;
00269   MSG("SubShowerSR",Msg::kVerbose)
00270     << "vtxdiff " <<vtxdiff << " tmaxpl " << tmaxpl <<endl;
00271   
00272   Int_t Ns         = 0;
00273   Int_t Ns1        = 0;
00274   Double_t stpden  = 0.;
00275   Double_t stpsp   = 0.;
00276   double phplbeg   = 0.;
00277   double phplend   = 0.;
00278   Double_t avgbeg  = 0.;
00279   Double_t posdif  = 0.;
00280   double preavg    = 0.;
00281   double prephpl   = 0.;
00282   double prenstppl = 0.;
00283 
00284   for (int pl = begpl;pl<=endpl;pl++){
00285     int nstppl  = 0;
00286     double low  = 200.;
00287     double high = -200.;
00288     double phpl = 0.;
00289     double avg  = 0.;
00290     for(int i=0;i<nstp;i++){
00291       if(plane[i]==pl) {
00292         nstppl++;
00293         phpl += ph[i];
00294         avg +=tpos[i]*ph[i]; 
00295         if(tpos[i]/T_EFF*100.>high) high = tpos[i]/T_EFF*100.;
00296         if(tpos[i]/T_EFF*100.<low) low = tpos[i]/T_EFF*100.;
00297       }
00298     }
00299     if (nstppl>0) {
00300       Ns++;
00301       avgstpperpln += Double_t(nstppl);
00302       stpden += double(nstppl)*phpl/(high-low+1.);
00303       avg = avg/T_EFF*100./phpl;    
00304       if (nstppl==1) Ns1++;
00305       if(pl>begpl && phplbeg!=0. && preavg!=0.) {
00306         double avgdif = ( (phpl + prephpl)*(avg - preavg) / 
00307                           (2.*double(nstppl + prenstppl) ) );
00308         stpsp += TMath::Abs(avgdif);
00309         if (avgdif>=0) posdif += 1;
00310         MSG("SubShowerSR",Msg::kVerbose) 
00311           << "avg " << avg << " preavg " << preavg << " nstppl " << nstppl
00312           << " prenstppl " << prenstppl << " phpl " << phpl << " prephpl "
00313           << " avgdif " << avgdif << " posdif " << posdif 
00314           << " stpsp " << stpsp << endl;
00315       }
00316       if(pl>=begpl&&phplbeg==0.) {
00317         phplbeg = phpl;
00318         avgbeg = avg;
00319       }
00320       if(pl<=endpl) phplend = phpl;
00321       preavg = avg;
00322       prephpl = phpl;
00323       prenstppl = nstppl;
00324     }
00325   }
00326   stpden = stpden/totw;
00327   
00328   if (Ns>1) {
00329     stpsp = stpsp/(totw-(phplbeg+phplend)/2.);
00330     posdif = posdif/double(Ns-1);
00331   }
00332   else {
00333     stpsp = 100.;
00334     posdif = 0.5;
00335   }  
00336   MSG("SubShowerSR",Msg::kVerbose) << "stpden " << stpden << " stpsp " 
00337                                    << stpsp << " posdif " << posdif << endl;
00338 
00339   if(endpl==begpl) {
00340     avgstpperpln = 999;
00341     trkfra = 0;
00342   }
00343   else {
00344     avgstpperpln /= (Double_t(endpl - begpl)/2. + 1.);
00345     trkfra = double(Ns1)/double(Ns);
00346   }  
00347   xtkfra = double(Ns0)/double(nstp);
00348   avgph = totph/double(nstp);
00349   Float_t maxstpph = 0.;
00350   Float_t begstpph = 0.;
00351   Int_t nstpmaxpl  = 0;
00352 
00353   for(Int_t i=0;i<nstp;i++){
00354     if(plane[i]==maxpl&&ph[i]>maxstpph){
00355       maxstp = strip[i];
00356       maxstpph = ph[i];
00357       maxtpos = tpos[i];
00358     }
00359     if(plane[i]==maxpl){
00360       maxz = z[i];
00361       nstpmaxpl++;
00362     }
00363     if(plane[i]==begpl){
00364       if(ph[i]>begstpph){
00365         begz = z[i];
00366         begstp = strip[i];
00367         begstpph = ph[i];
00368       }
00369       if(time[i]<mint) mint = time[i];
00370     }
00371     Bool_t gotendz = false;
00372     if(!gotendz&&plane[i]==endpl){
00373       endz = z[i];
00374       gotendz = true;
00375     }
00376   }
00377 
00378   Float_t ss[2][2] = {{0}};
00379   for (int m = 0;m<2;m++){
00380     for (int n = 0;n<2;n++){
00381       ss[m][n] = 0.;
00382     }
00383   }
00384 
00385   Double_t cnt = 0;
00386   for(Int_t i=0;i<nstp;i++){
00387     Double_t dtpos = tpos[i] - maxtpos;
00388     Double_t dz = z[i] - maxz;
00389     ss[0][0]+=dz*dz*ph[i];
00390     ss[0][1]+=dtpos*dz*ph[i];
00391     ss[1][1]+=dtpos*dtpos*ph[i];
00392     cnt+=1;
00393     normpos += dtpos*dtpos+dz*dz;    
00394   }
00395   ss[1][0]=ss[0][1];
00396     
00397   if(totw*normpos!=0){
00398     
00399     TMatrixDSym ms(2);
00400     for (int m = 0;m<2;m++){
00401       for (int n = 0;n<2;n++){
00402         ms[m][n] = 0.;
00403       }
00404     }
00405 
00406     Double_t ang[2];    
00407     for(Int_t i=0;i<2;i++) ang[i] = -999;
00408     for (Int_t i=0;i<2;i++){
00409       for (Int_t j=0;j<2;j++){
00410         ms(i,j)=ss[i][j]/totw/normpos;
00411       }
00412     }
00413 
00414     if(ms.IsA()){
00415       TMatrixDSymEigen eigen(ms);
00416       TMatrixD EVect1 = eigen.GetEigenVectors();          
00417       TVectorD EVal1 = eigen.GetEigenValues();
00418       TVector2 ev0(EVect1(0,0),EVect1(1,0));
00419       TVector2 ev1(EVect1(0,1),EVect1(1,1));
00420       ang[0] = ev0.Phi()*180/TMath::Pi();
00421       ang[1] = ev1.Phi()*180/TMath::Pi();
00422       eval[0] = EVal1(0);
00423       eval[1] = EVal1(1);
00424     }
00425     else {
00426       eval[0] = 0.;
00427       eval[1] = 0.;
00428     }
00429     for (Int_t i=0;i<2;i++){
00430       if (ang[i]>90&&ang[i]<270) ang[i]-=180;
00431       if (ang[i]>270&&ang[i]<=360) ang[i]-=360;
00432       if (ang[i]==90||ang[i]==270) ang[i]-=0.01;
00433     }
00434     
00435     Theta = ang[0];
00436     Theta1 = ang[1];
00437     Theta_rad = (Theta/180)*TMath::Pi();
00438     Theta_rad1 = (Theta1/180)*TMath::Pi();
00439     double vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00440     double vs1 = maxtpos+tan(Theta_rad1)*(begz-maxz);
00441     MSG("SubShowerSR",Msg::kVerbose) << "vs " << vs << " vs1 " << vs1 
00442                                      << " avgbeg " << avgbeg*T_EFF/100. 
00443                                      << " Theta " << Theta << endl;
00444     if(ss[0][0]==0. || (TMath::Abs(Theta)>TMath::Abs(Theta1) && 
00445                         TMath::Abs(avgdev*T_EFF/100.-vs) > 2 * 
00446                         TMath::Abs(avgbeg*T_EFF/100.-vs1) && 
00447                         (stpden<0.5 || (stpsp>0.5 && 
00448                                         (posdif>0. && posdif<1.)
00449                                         || Ns<=2)))){
00450       double swap = Theta;
00451       Theta = Theta1;
00452       Theta1 = swap;
00453     }
00454     Theta_rad = (Theta/180)*TMath::Pi();
00455     Theta_rad1 = (Theta1/180)*TMath::Pi();
00456     vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00457     MSG("SubShowerSR",Msg::kVerbose) 
00458       << "vs " << vs << " avgbeg " << avgbeg*T_EFF/100. 
00459       << " Theta " << Theta << endl;
00460     if(vs>(uppeststp+(uppeststp-loweststp)*1.5/2.)*T_EFF/100.||
00461        vs<(loweststp-(uppeststp-loweststp)*1.5/2.)*T_EFF/100.
00462        ||vs<-4.||vs>4.) {
00463       double bigang = Theta;
00464       if(TMath::Abs(tan(Theta_rad))<TMath::Abs(tan(Theta_rad1))) bigang = Theta1;
00465       Theta = 0.;
00466       Theta1 = TMath::Abs(tan(bigang))/tan(bigang)*89.99;
00467       Theta_rad = (Theta/180)*TMath::Pi();
00468       Theta_rad1 = (Theta1/180)*TMath::Pi();
00469     }
00470     MSG("SubShowerSR",Msg::kVerbose) 
00471       << "vs " << vs << " avgbeg " << avgbeg*T_EFF/100. 
00472       << " Theta " << Theta << endl;
00473     
00474     std::vector<Double_t> shwstpc(nstp,0);
00475     std::vector<Double_t> shwstpc1(nstp,0);
00476     std::vector<Double_t> shwstpc2(nstp,0);
00477     std::vector<Double_t> wshwstpc(nstp,0);
00478     for(Int_t i=0;i<nstp;i++) shwstpc[i] = -1.;
00479     for(Int_t i=0;i<nstp;i++) shwstpc1[i] = -1.;
00480     for(Int_t i=0;i<nstp;i++) shwstpc2[i] = -1.;
00481     for(Int_t i=0;i<nstp;i++) wshwstpc[i] = -1.;
00482         
00483     avgdev = 0.;
00484     wavgdev = 0.;
00485 
00486     Float_t offstp = 0;
00487     Float_t offstp2 = 0;
00488 
00489     out2Rmph = 0;
00490     for(Int_t i=0;i<nstp;i++){
00491       shwstpc[i]=TMath::Abs(tpos[i]
00492                       -(maxtpos+tan(Theta_rad)
00493                         *(z[i]-maxz)))*TMath::Cos(Theta_rad);
00494       shwstpc2[i]=TMath::Abs(tpos[i]
00495                        -(maxtpos+tan(Theta_rad1)
00496                          *(z[i]-maxz)))*TMath::Cos(Theta_rad1);
00497       if (shwstpc[i]>2*RM_EFF/100.) out2Rmph+=ph[i];
00498       if (shwstpc[i]>1) offstp++;
00499       if (shwstpc2[i]>1) offstp2++;
00500       avgdev1 += shwstpc[i];
00501       avgdev2 += shwstpc2[i];
00502       wavgdev1 += shwstpc[i]*ph[i];
00503       wavgdev2 += shwstpc2[i]*ph[i];
00504     }
00505 
00506     out2Rmph/=totw;
00507     avgdev1/=nstp;
00508     avgdev2/=nstp;
00509     wavgdev1/=totw;
00510     wavgdev2/=totw;
00511     avgdev = avgdev1;
00512     wavgdev = wavgdev1;
00513     offstp/=nstp;
00514     offstp2/=nstp;
00515   } 
00516   double vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00517 
00518   //check that vertex from slope is within 25cm of 
00519   //subshower transverse extent to prevent unphysical values
00520   if(vs<loweststp*T_EFF/100.-0.25 || 
00521      vs>uppeststp*T_EFF/100.+0.25) {
00522     vs = avgbeg*T_EFF/100.;
00523     Theta_rad = 0.;
00524   }
00525   if(cch.GetPlaneView()==PlaneView::kX || 
00526      cch.GetPlaneView()==PlaneView::kU) {
00527     cch.SetVtxU(vs); 
00528   }
00529   if(cch.GetPlaneView()==PlaneView::kY || 
00530      cch.GetPlaneView()==PlaneView::kV) {
00531     cch.SetVtxV(vs);
00532   }
00533   cch.SetVtxZ(begz);
00534   cch.SetEndZ(endz);
00535   cch.SetVtxT(mint);
00536   cch.SetVtxPlane(begpl);
00537   cch.SetEndPlane(endpl);
00538   cch.SetEnergy(Eguess/2.);
00539   cch.SetSlope(tan(Theta_rad));
00540   cch.SetAvgDev(wavgdev); 
00541 
00542   LongPl.clear();
00543   LongPh.clear();
00544   LongFlu.clear();
00545   LongPhpe.clear();
00546   TranStp.clear();  
00547   TranPh.clear();
00548   TranFlu.clear();
00549 
00550   EMFluctuation EMFlu(Eguess);
00551   nn0l = 0;
00552   nn0t = 0;
00553 
00554   if(cos(Theta_rad)!=0.){
00555     Bool_t xs = false;
00556     Int_t dpl = 2;
00557     if(det==1&&(begpl-249)*(endpl-249)<0) xs = true;
00558     for(Int_t j=begpl;j<=endpl;j+=dpl){
00559       dpl = 2;
00560       if(xs&&(begpl-249)*(j-249)<0) {
00561         dpl = 3;
00562         xs = false;
00563       }
00564       Double_t totphPl = 0.;
00565       Double_t totphPl_pe = 0.;
00566       Double_t plz = -1.;
00567       Double_t nstpPl = 0;
00568       for(Int_t i=0;i<nstp;i++){
00569         if(plane[i]==j && ph[i]>0) {
00570           totphPl+=ph[i];
00571           totphPl_pe+=ph_pe[i];
00572           plz = z[i];
00573           Double_t zt_vtx_corrected[2];
00574           zt_vtx_corrected[0] = z[i]-begz;
00575           zt_vtx_corrected[1] = tpos[i]-vs-(z[i]-begz)*tan(Theta_rad);
00576           nstpPl++;
00577         }
00578       }
00579       LongPl.push_back(plz-begz);
00580       LongPh.push_back(totphPl);
00581       LongPhpe.push_back(totphPl_pe);
00582       if(totphPl>0.) nn0l++;
00583     }
00584 
00585     double z_adj = 0.;
00586     std::list<Double_t>::iterator PlIt = LongPl.begin();
00587     std::list<Double_t>::iterator PhIt = LongPh.begin();
00588     std::list<Double_t>::iterator PhpeIt = LongPhpe.begin();
00589     double t0 = (*PhIt);
00590     PhIt++;
00591     double t1 = (*PhIt);
00592     PhIt = LongPh.begin();
00593     if(vtxdiff>1.25) z_adj = Int_t((1-vtxdiff)*tmaxpl);
00594     if(t0>t1/2.&&t0<t1) z_adj = 1;
00595     else if(t0>=t1) {
00596       z_adj = 2;
00597       LongPl.push_front(0);
00598       LongPh.push_front(0.);
00599       LongPhpe.push_front(0.);
00600     }
00601 
00602     Int_t inEM = 0;
00603     PlIt = LongPl.begin();
00604     while(PlIt!=LongPl.end()){
00605       if(z_adj!=0) {
00606         (*PlIt) += z_adj*D_EFF/100.;
00607         begz    += z_adj*D_EFF/100.;
00608         endz    += z_adj*D_EFF/100.;
00609       }
00610       double fluem = 0.;
00611       if ((*PlIt)>=0.) {
00612         fluem = EMFlu.CalcLongFluc((*PlIt)/cos(Theta_rad));
00613         inEM++;
00614       }
00615       double error = 0.;
00616       if ((*PhpeIt)<=0.001) error = fluem; 
00617       else {
00618         MSG("SubShowerSR",Msg::kVerbose)
00619           << "fluem " << fluem << " *PhpeIt " << (*PhpeIt) << endl;
00620         error = TMath::Sqrt(TMath::Power(fluem,2) + 
00621                             TMath::Power(1./TMath::Sqrt(*PhpeIt),2));
00622       }
00623       if(inEM==1) error = TMath::Sqrt(TMath::Power(error,2)+1.);
00624       MSG("SubShowerSR",Msg::kVerbose) 
00625         << (*PlIt) << " long ph pe " << (*PhpeIt) 
00626         << " error " << error << " " << fluem << endl;
00627       LongFlu.push_back(error);
00628       PlIt++;
00629       PhpeIt++;      
00630     }
00631     MSG("SubShowerSR",Msg::kVerbose)<<"z_adj "<<z_adj<<endl; 
00632 
00633     PlIt = LongPl.end();
00634     PlIt--;
00635     Double_t miss_z = (*PlIt);
00636     PlIt = LongPl.begin();
00637     t0 = *(LongPl.begin());
00638     while(miss_z-t0<2*tmax*X0_EFF/100.){
00639       miss_z+=2*D_EFF/100.;
00640       LongPl.push_back(miss_z);
00641       LongPh.push_back(0.);
00642       Double_t miss_flu = EMFlu.CalcLongFluc(miss_z/cos(Theta_rad));
00643       MSG("SubShowerSR",Msg::kVerbose)<<"miss_flu "<<miss_flu<<endl;
00644       LongFlu.push_back(miss_flu);
00645     }
00646     MSG("SubShowerSR",Msg::kVerbose)<< "0 loweststp " << loweststp
00647                                     << " uppeststp " << uppeststp << endl;
00648     Double_t Thetabd = atan(ST_WID/(Double_t(LongPl.size()) * 
00649                                     D_EFF*2.))*180./TMath::Pi();
00650     Double_t dThetai = 1./180.*TMath::Pi();
00651     if(TMath::Abs(Thetabd)<5.) dThetai = TMath::Abs(Thetabd)/180.*TMath::Pi()/5.;
00652     Double_t maxph0 = 0.;
00653     Double_t maxvs0 = vs;
00654     Double_t maxtheta0 = Theta_rad;
00655     for(Int_t k=-6;k<=6;k++){
00656       for(Int_t j=-5;j<=5;j++){
00657         Double_t totph0 = 0.;
00658         for(Int_t i=0;i<nstp;i++){
00659           Double_t zt_vtx_corrected[2];
00660           zt_vtx_corrected[0] = z[i]-begz;
00661           zt_vtx_corrected[1] = tpos[i] - (vs+ST_WID/100.*k/4.) - 
00662             (z[i]-begz) * tan(Theta_rad+j*dThetai);
00663           if(TMath::Abs(zt_vtx_corrected[1]-0.)<=T_EFF/2./100.&&ph[i]>0) {
00664             totph0+=ph[i];
00665           }
00666         }
00667         if(totph0>maxph0){
00668           maxph0 = totph0;
00669           maxvs0 = vs+ST_WID/100.*k/4.;
00670           maxtheta0 = Theta_rad+j*dThetai;
00671           MSG("SubShowerSR",Msg::kVerbose)
00672             <<"Theta_rad "<<Theta_rad<<" vs "<<vs
00673             <<" maxtheta0 "<<maxtheta0<<" maxvs0 "<<maxvs0
00674             <<" maxph0 "<<maxph0<<" j "<<j<<" k "<<k
00675             <<" dThetai "<<dThetai<<endl;
00676         }
00677       }
00678     }
00679     MSG("SubShowerSR",Msg::kDebug) <<"Theta_rad "<<Theta_rad
00680                                    <<" vs "<<vs
00681                                    <<" maxtheta0 "<<maxtheta0
00682                                    <<" maxvs0 "<<maxvs0<<endl;
00683     Theta_rad = maxtheta0;
00684     vs = maxvs0;
00685     for(Int_t j=int((loweststp-uppeststp)/2.)-1;
00686         j<=int((uppeststp-loweststp)/2.)+1;j++){
00687 
00688       Double_t totphStp = 0.;
00689       Double_t totphStp_pe = 0.;
00690       Double_t stptpos = j*T_EFF/100.;
00691 
00692       Double_t nstpTr = 0.;
00693       for(Int_t i=0;i<nstp;i++){
00694         Double_t zt_vtx_corrected[2];
00695         zt_vtx_corrected[0] = z[i]-begz;
00696         zt_vtx_corrected[1] = tpos[i]-vs-(z[i]-begz)*tan(Theta_rad);
00697 
00698         if(TMath::Abs(zt_vtx_corrected[1]-stptpos)<=T_EFF/2./100.&&ph[i]>0) {
00699           totphStp+=ph[i];
00700           totphStp_pe+=ph_pe[i];
00701           nstpTr++;
00702         }
00703       }
00704       TranStp.push_back(stptpos);
00705       TranPh.push_back(totphStp);
00706       if(totphStp>0.) nn0t++;
00707       double fluem = EMFlu.CalcTranFluc(stptpos*cos(Theta_rad));
00708       double error = 0.;
00709       if (totphStp_pe<=0.001) error = fluem;
00710       else {
00711         MSG("SubShowerSR",Msg::kVerbose)<<"fluem "<<fluem
00712                                         <<" totphStp_pe "<<totphStp_pe<<endl;
00713         error = TMath::Sqrt(TMath::Power(fluem,2) + 
00714                             TMath::Power(1./TMath::Sqrt(totphStp_pe),2)); 
00715         MSG("SubShowerSR",Msg::kVerbose)<<"tran ph "<<totphStp_pe
00716                                         <<" error "<<error<<" "<<fluem<<endl;
00717       }
00718       TranFlu.push_back(error);     
00719     }
00720 
00721     std::list<Double_t>::iterator StpIt = TranStp.begin();
00722     Double_t miss_tn = *(TranStp.begin());
00723     StpIt = TranStp.end();
00724     StpIt--;
00725     Double_t miss_tp = *StpIt;
00726 
00727     while(miss_tp-miss_tn<4*RM_EFF/100.){
00728       miss_tn-=T_EFF/100.;
00729       miss_tp+=T_EFF/100.;
00730       TranStp.push_front(miss_tn);
00731       TranStp.push_back(miss_tp);
00732       TranPh.push_front(0.);
00733       TranPh.push_back(0.);
00734       Double_t miss_flun = EMFlu.CalcLongFluc(miss_tn*cos(Theta_rad));
00735       Double_t miss_flup = EMFlu.CalcLongFluc(miss_tp*cos(Theta_rad));
00736       MSG("SubShowerSR",Msg::kVerbose)<<"miss_flun "<<miss_flun
00737                                       <<" miss_flup "<<miss_flup<<endl;
00738       TranFlu.push_front(miss_flun);
00739       TranFlu.push_back(miss_flup);
00740       loweststp--;
00741       uppeststp++;
00742     }
00743   }
00744 }

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

Implements AlgBase.

Definition at line 106 of file AlgSubShowerSR.cxx.

References CalculateEnergyVertexAngle(), fDoEMFit, fECut, fMaxDiffCut, fMaxXTalkStpPH, fMIPperGeV, fOut2RmPHCut, fPECut, fTrkFraCut, fTrkLikeStpCut, fVtxDiffCut, fXtkFraCut, CandContext::GetDataIn(), VldContext::GetDetector(), Registry::GetDouble(), Registry::GetInt(), Msg::kDebug, Detector::kFar, PlaneView::kUnknown, MSG, and SubShowerID().

00107 {
00108 
00109   MSG("SubShowerSR",Msg::kDebug) << "In AlgSubShowerSR::RunAlg" << endl;
00110   assert(ch.InheritsFrom("CandSubShowerSRHandle"));
00111   CandSubShowerSRHandle &cch = dynamic_cast<CandSubShowerSRHandle &>(ch);
00112   assert(cx.GetDataIn());
00113   assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00114   const TObjArray *tary = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00115   if(tary->GetLast()<=-1) return;
00116   Int_t detector = 0;
00117   for (Int_t i=0; i<=tary->GetLast(); i++) {
00118     TObject *tobj = tary->At(i);
00119     assert(tobj->InheritsFrom("CandStripHandle"));
00120     CandStripHandle *striphandle = dynamic_cast<CandStripHandle*>(tobj);
00121     cch.AddDaughterLink(*striphandle);
00122     if(cch.GetPlaneView()==PlaneView::kUnknown) 
00123       cch.SetPlaneView(striphandle->GetPlaneView());
00124     if (i==0) {
00125       const VldContext *vldcptr = striphandle->GetVldContext();
00126       assert(vldcptr);
00127       VldContext vldc = *vldcptr;                                         
00128       if(vldc.GetDetector()==Detector::kFar) detector = 1;
00129     }
00130   }
00131 
00132   fPECut         = ac.GetDouble("PECut");  
00133   fECut          = ac.GetDouble("ECut");
00134   fMaxDiffCut    = ac.GetDouble("MaxDiffCut");
00135   fVtxDiffCut    = ac.GetDouble("VtxDiffCut");
00136   fOut2RmPHCut   = ac.GetDouble("Out2RmPHCut");
00137   fTrkFraCut     = ac.GetDouble("TrkFraCut");
00138   fTrkLikeStpCut = ac.GetDouble("TrkLikeStpCut");
00139   fXtkFraCut     = ac.GetDouble("XtkFraCut");
00140   fMaxXTalkStpPH = ac.GetDouble("MaxXTalkStpPH");
00141   fDoEMFit       = ac.GetInt("DoEMFit");
00142   fMIPperGeV     = ac.GetDouble("MIPperGeV");
00143 
00144   CalculateEnergyVertexAngle(cch,detector);
00145   SubShowerID(cch);
00146 
00147 }

void AlgSubShowerSR::SubShowerID ( CandSubShowerSRHandle cch  )  [private]

Definition at line 748 of file AlgSubShowerSR.cxx.

References avgph, avgstpperpln, begz, bothfitprob, chi2func, clusterlong, clustertran, D_EFF, MuELoss::e, Eguess, endz, err(), fDoEMFit, fECut, fMaxDiffCut, fMaxXTalkStpPH, fOut2RmPHCut, fPECut, fTrkFraCut, fTrkLikeStpCut, fVtxDiffCut, fXtkFraCut, CandSubShowerSRHandle::GetEnergy(), CandRecoHandle::GetNStrip(), Msg::kDebug, ClusterType::kEMLike, ClusterType::kHadLike, ClusterType::kTrkLike, ClusterType::kUnknown, Msg::kVerbose, ClusterType::kXTalk, LongFlu, LongHist, LongHistb, LongPh, LongPl, LongProf, LongProfb, loweststp, maxdiff, maxStpPHinSS, MSG, nn0l, nn0t, out2Rmph, CandSubShowerSRHandle::SetClusterID(), CandSubShowerSRHandle::SetProbEM(), T_EFF, Theta_rad, totw, TranFlu, TranHist, TranHistb, TranPh, TranProf, TranProfb, TranStp, trkfra, uppeststp, vtxdiff, and xtkfra.

Referenced by RunAlg().

00748                                                           {
00749 
00750   ClusterType::ClusterType_t id = ClusterType::kUnknown;
00751 
00752   if((xtkfra>fXtkFraCut || avgph<2*fPECut) && maxStpPHinSS<fMaxXTalkStpPH ) 
00753     id = ClusterType::kXTalk;
00754   else if(trkfra>fTrkFraCut || avgstpperpln<fTrkLikeStpCut)
00755     id = ClusterType::kTrkLike;
00756   else if(cch.GetEnergy()<fECut) 
00757     id = ClusterType::kHadLike;
00758   else if(maxdiff>fMaxDiffCut) 
00759     id = ClusterType::kHadLike;
00760   else if(vtxdiff>fVtxDiffCut) 
00761     id = ClusterType::kHadLike;
00762   else if(out2Rmph>fOut2RmPHCut) 
00763     id = ClusterType::kHadLike;
00764   else id = ClusterType::kEMLike;
00765 
00766   cch.SetClusterID(id);
00767   MSG("SubShowerSR",Msg::kVerbose)<<"id "<<id<<" E "<<cch.GetEnergy()<<endl;
00768 
00769   bothfitprob = 0.;
00770   Bool_t makeplots = false;
00771   MSG("SubShowerSR",Msg::kDebug) 
00772     << "nn0l " << nn0l << " LongPl.size() " << LongPl.size() 
00773     << " nn0t " << nn0t << " TranStp.size() " << TranStp.size() << endl;
00774 
00775   //checks before trying EM fit:
00776   if(fDoEMFit && 
00777      (id==ClusterType::kEMLike || 
00778       (id==ClusterType::kHadLike && cch.GetEnergy()>=fECut)) &&
00779      LongPl.size()-nn0l <= 3 && TranStp.size()-nn0t <= 3 ){
00780 
00781     MSG("SubShowerSR",Msg::kDebug)<<"in fitting "<<endl;
00782 
00783     std::list<Double_t>::iterator PlIt = LongPl.begin();
00784     std::list<Double_t>::iterator PhIt = LongPh.begin();
00785     std::list<Double_t>::iterator PluIt = LongFlu.begin();
00786     std::list<Double_t>::iterator StpIt = TranStp.begin();
00787     std::list<Double_t>::iterator StphIt = TranPh.begin();
00788     std::list<Double_t>::iterator StpluIt = TranFlu.begin();
00789 
00790     MSG("SubShowerSR",Msg::kVerbose) 
00791       << "funct boundary " << ((loweststp-uppeststp)/2.-1)*T_EFF-T_EFF/2. 
00792       << " " << ((uppeststp-loweststp)/2.+1)*T_EFF+T_EFF/2. << endl;
00793     
00794     Double_t thetarange = Theta_rad*0.05;
00795     if(thetarange<10./180.*TMath::Pi()) thetarange = 10./180.*TMath::Pi();
00796     Double_t erange = Eguess*0.05;
00797     if(erange<0.5) erange = 0.5;
00798     
00799     clusterlong->SetRange(0.-D_EFF/2.,(endz-begz+(*PlIt))*100+D_EFF);
00800     clusterlong->SetParameter(0,Eguess);
00801     clusterlong->SetParLimits(0,Eguess-erange,Eguess+erange);
00802     clusterlong->SetParameter(1,Theta_rad);
00803     clusterlong->SetParLimits(1,Theta_rad-thetarange,Theta_rad+thetarange);
00804     
00805     clustertran->SetRange(((loweststp-uppeststp)/2.-1)*T_EFF,
00806                           ((uppeststp-loweststp)/2.+1)*T_EFF);
00807     clustertran->SetParameter(0,Eguess);
00808     clustertran->SetParLimits(0,Eguess-erange,Eguess+erange);
00809     clustertran->SetParameter(1,Theta_rad);
00810     clustertran->SetParLimits(1,Theta_rad-thetarange,Theta_rad+thetarange);
00811     
00812     LongProf->Reset();
00813     LongProf->SetBins(int((endz-begz+(*PlIt))*100+2*D_EFF),
00814                       0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);    
00815     LongHist->Reset();
00816     LongHist->SetBins(int((endz-begz+(*PlIt))*100+2*D_EFF),
00817                       0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00818     TranProf->Reset();
00819     TranProf->SetBins(int((uppeststp-loweststp+2)*T_EFF),
00820                       ((loweststp-uppeststp)/2.-1)*T_EFF,
00821                       ((uppeststp-loweststp)/2.+1)*T_EFF);
00822     TranHist->Reset();
00823     TranHist->SetBins(int((uppeststp-loweststp+2)*T_EFF),
00824                       ((loweststp-uppeststp)/2.-1)*T_EFF,
00825                       ((uppeststp-loweststp)/2.+1)*T_EFF);
00826     LongProfb->Reset();
00827     LongProfb->SetBins(int((endz-begz+(*PlIt))*100+2*D_EFF),
00828                        0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);    
00829     LongHistb->Reset();
00830     LongHistb->SetBins(int((endz-begz+(*PlIt))*100+2*D_EFF),
00831                        0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00832     TranProfb->Reset();
00833     TranProfb->SetBins(int((uppeststp-loweststp+2)*T_EFF),
00834                        ((loweststp-uppeststp)/2.-1)*T_EFF,
00835                        ((uppeststp-loweststp)/2.+1)*T_EFF);
00836     TranHistb->Reset();
00837     TranHistb->SetBins(int((uppeststp-loweststp+2)*T_EFF),
00838                        ((loweststp-uppeststp)/2.-1)*T_EFF,
00839                        ((uppeststp-loweststp)/2.+1)*T_EFF);
00840 
00841     MSG("SubShowerSR",Msg::kVerbose) << "loweststp " << loweststp 
00842                                      << " uppeststp " << uppeststp << endl;
00843 
00844     Double_t Minchi2Both = 100000;
00845     Double_t BestEBoth = Eguess;
00846     Double_t BestThetaBoth = Theta_rad;
00847     double bestnormlb = 1;
00848     double bestnormtb = 1;
00849     Bool_t fittedboth = false;
00850     Double_t dE = 0.1;
00851     Double_t Eini = Eguess-erange;
00852     Double_t Ein = Eini;
00853 
00854     MSG("SubShowerSR",Msg::kDebug) 
00855       << "done initial stuff " <<Eguess << " " << Ein 
00856       << " " << Theta_rad << " " << thetarange << endl;
00857 
00858     while(Ein>=Eini&&Ein<=Eguess+erange){
00859       MSG("SubShowerSR",Msg::kVerbose)<<"Ein "<<Ein<<endl;
00860       clusterlong->SetParameter(0,Ein);
00861       Double_t dTheta = 2.5/180.*TMath::Pi();
00862       Double_t Thetaini = Theta_rad-thetarange;
00863       Double_t Thetain = Thetaini;
00864       while(Thetain>=Thetaini&&Thetain<=Theta_rad+thetarange){
00865         MSG("SubShowerSR",Msg::kVerbose)<<"Thetain "<<Thetain<<endl;
00866         if(cos(Thetain)!=0.){
00867           clusterlong->SetParameter(1,Thetain);
00868           Double_t norml = 0.0;
00869           std::vector<Double_t> funcl(LongPl.size(),0);
00870           std::vector<Double_t> errl(LongPl.size(),0);
00871           Int_t nfpl = 0;
00872           UInt_t j = 0;
00873           Int_t nlong = 0;
00874           PhIt = LongPh.begin();
00875           while(PhIt!=LongPh.end()){
00876             if((*PhIt)>0.) nlong++;
00877             PhIt++;
00878           }
00879           PhIt = LongPh.begin();
00880           PlIt = LongPl.begin();
00881           while(PlIt!=LongPl.end()){
00882             funcl[j] = 0.;
00883             if((*PlIt)>=0.) funcl[j] = clusterlong->Eval((*PlIt)*100.);
00884             MSG("SubShowerSR",Msg::kVerbose)<<"funcl[j] "<<funcl[j]<<endl;
00885             norml += funcl[j];
00886             if (funcl[j]>0.||(*PlIt)<0.) nfpl++;
00887             PlIt++;
00888             j++;
00889           }
00890           MSG("SubShowerSR",Msg::kVerbose)<<"nfpl "<<nfpl<<endl;
00891           Double_t chi2Long = 0.;
00892           PlIt = LongPl.begin();
00893           PhIt = LongPh.begin();
00894           PluIt = LongFlu.begin();
00895           j = 0;
00896           MSG("SubShowerSR",Msg::kVerbose)
00897             << "nfpl " << nfpl << " norml " << norml << " nlong " << nlong
00898             << " Ein " << Ein << " Thetain " << Thetain << endl;
00899 
00900           if(double(nfpl)/double(LongPl.size())>0.5 && 
00901              norml>1e-3 && nlong>1){
00902             
00903             while(PlIt!=LongPl.end()){
00904               Double_t data = (*PhIt)/totw*norml;
00905               Double_t err = funcl[j]*(*PluIt);
00906               if ((*PluIt)==0.0) err = 0.0001;
00907               UInt_t sub = j;
00908               while(funcl[sub]<0.05*norml&&j==0) {
00909                 sub++;
00910                 err = (funcl[sub]+funcl[j])/2.*(*PluIt);
00911                 MSG("SubShowerSR",Msg::kVerbose)<<"err "<<err<<" sub "<<endl;
00912                 if(err>0.0001) break;
00913                 else if(sub==LongPl.size()-1||(sub-j)>=3){
00914                   err = 0.0001;
00915                   break;
00916                 }
00917               }
00918               sub = j;
00919               while(funcl[sub]<0.05*norml&&j>0) {
00920                 sub--;
00921                 err = (funcl[sub]+funcl[j])/2.*(*PluIt);
00922                 MSG("SubShowerSR",Msg::kVerbose)<<"err "<<err<<" sub "<<endl;
00923                 if(err>0.0001) break;
00924                 else if(sub==0||(j-sub)>=3){
00925                   err = 0.0001;
00926                   break;
00927                 }
00928               }
00929               
00930               MSG("SubShowerSR",Msg::kDebug) 
00931                 << " LongPl.size() " << LongPl.size() 
00932                 << " Eguess " <<Eguess << " norml " << norml 
00933                 << " Theta " << Theta_rad << " j " << j 
00934                 << " LongPl.at(j) " << (*PlIt) 
00935                 << " err " << funcl[j] << " " << funcl[j-1] 
00936                 << " " << funcl[j+1] << " " << (*PluIt) 
00937                 << " err " << err << " data - funcl " << data-funcl[j] << endl;
00938                       
00939               if(err<0.0001) err = 0.0001;
00940               MSG("SubShowerSR",Msg::kVerbose)
00941                 <<"err "<<err<<" funcl[j] "<<funcl[j]<<endl;
00942               
00943               if(j==0) 
00944                 errl[j] = TMath::Sqrt(TMath::Power(err,2) + 
00945                                       TMath::Power((funcl[j+1] - 
00946                                                     funcl[j])/2.,2));
00947               else if(j==LongPl.size()-1) 
00948                 errl[j] = TMath::Sqrt(TMath::Power(err,2) + 
00949                                       TMath::Power((funcl[j] - 
00950                                                     funcl[j-1])/2.,2));
00951               else 
00952                 errl[j] = TMath::Sqrt(TMath::Power(err,2) + 
00953                                       TMath::Power((funcl[j+1] - 
00954                                                     funcl[j-1])/4.,2));
00955               
00956               MSG("SubShowerSR",Msg::kVerbose) 
00957                 << "data-funcl[j] " << data-funcl[j] 
00958                 << " errl[j] " << errl[j] << endl;
00959               
00960               chi2Long += ( TMath::Power(data-funcl[j],2) / 
00961                             TMath::Power(errl[j],2) );
00962               
00963               MSG("SubShowerSR",Msg::kVerbose)
00964                 << "chi2long " << (*PlIt) << " " << data 
00965                 << " " << funcl[j] << " " << (*PluIt) << " " << err 
00966                 << " " << TMath::Power(data-funcl[j],2)/TMath::Power(err,2) 
00967                 << " " << chi2Long << endl;
00968               
00969               PlIt++;
00970               PhIt++;
00971               PluIt++;
00972               j++;
00973             }
00974 
00975             MSG("SubShowerSR",Msg::kVerbose) 
00976               << "long fit " << chi2Long/(LongPl.size()-1) 
00977               << "norml " << LongPl.size() << " " << norml << endl;
00978           }
00979           else chi2Long = 100000;
00980           
00981           Double_t normt = 0.0;
00982           std::vector<Double_t> funct(TranStp.size(),0);
00983           std::vector<Double_t> errt(TranStp.size(),0);
00984           Int_t nfpt = 0;
00985           j = 0;
00986           Int_t ntran = 0;
00987           StphIt = TranPh.begin();
00988           while(StphIt!=TranPh.end()){
00989             if(*StphIt>0.) ntran++;
00990             StphIt++;
00991           }
00992           StphIt = TranPh.begin();
00993           StpIt = TranStp.begin();
00994           while(StpIt!=TranStp.end()){
00995             funct[j] = clustertran->Eval(*StpIt*100.);
00996             normt += funct[j];
00997             if(funct[j]>0.) nfpt++;
00998             StpIt++;
00999             j++;
01000           }
01001           
01002           Double_t chi2Tran = 0.;
01003           MSG("SubShowerSR",Msg::kVerbose)<<"nfpt "<<nfpt<<" "<<normt<<endl;
01004           if(double(nfpt)/double(TranStp.size())>0.5 && 
01005              normt>1e-3 && ntran>1){
01006             
01007             StpIt = TranStp.begin();
01008             StphIt = TranPh.begin();
01009             StpluIt = TranFlu.begin();
01010             j = 0;
01011             
01012             while(StpIt!=TranStp.end()){
01013               Double_t data = *StphIt/totw*normt;
01014               Double_t err = funct[j]*(*StpluIt);
01015               if (*StpluIt==0.0) err = 0.0001;
01016               UInt_t sub = j;
01017               while(funct[sub]<0.05*normt&&*StpIt<=0) {
01018                 sub++;
01019                 err = (funct[sub]+funct[j])/2.*(*StpluIt);
01020                 if(err>0.0001) break;
01021                 else if(sub==TranStp.size()-1||(sub-j)>=3){
01022                   err = 0.0001;
01023                   break;
01024                 }               
01025               }
01026               sub = j;
01027               while(funct[sub]<0.05*normt && *StpIt>0) {
01028                 sub--;
01029                 err = (funct[sub]+funct[j])/2.*(*StpluIt);
01030                 if(err>0.0001) break;
01031                 else if(sub==0||(j-sub)>=3){
01032                   err = 0.0001;
01033                   break;
01034                 }
01035               }
01036               MSG("SubShowerSR",Msg::kDebug) 
01037                 << "funct " << j << " " << *StpIt << " " << *StpluIt 
01038                 << " err " << " " << funct[j] << " " << funct[j-1]
01039                 << " " << funct[j+1] << " err " << err 
01040                 << " data - funcl " << data-funcl[j] << endl;
01041 
01042               if(err<0.0001) err = 0.0001;
01043               if(j==0) 
01044                 errt[j] = TMath::Sqrt(TMath::Power(err,2) + 
01045                                       TMath::Power((funct[j+1] - 
01046                                                     funct[j])/2.,2));
01047               else if(j==TranStp.size()-1) 
01048                 errt[j] = TMath::Sqrt(TMath::Power(err,2) + 
01049                                       TMath::Power((funct[j] - 
01050                                                     funct[j-1])/2.,2));
01051               else 
01052                 errt[j] = TMath::Sqrt(TMath::Power(err,2) + 
01053                                       TMath::Power((funct[j+1] - 
01054                                                     funct[j-1])/4.,2));
01055               
01056               MSG("SubShowerSR",Msg::kVerbose) << "data-funcl[j] " << 
01057                 data-funcl[j] << " errt[j] " << errt[j] << endl;
01058 
01059               chi2Tran += ( TMath::Power(data-funct[j],2) / 
01060                             (TMath::Power(errt[j],2)) );
01061               
01062               MSG("SubShowerSR",Msg::kVerbose) << "chi2tran " <<
01063                 TMath::Power(data-funcl[j],2)/TMath::Power(err,2) 
01064                                                << " " << chi2Tran << endl;
01065               
01066               StpIt++;
01067               StphIt++;
01068               StpluIt++;
01069               j++;
01070             }
01071             
01072             MSG("SubShowerSR",Msg::kVerbose) 
01073               << "tran fit " << chi2Tran/(TranStp.size()-1) 
01074               << "normt " << TranStp.size() << " " << normt << endl;
01075           }
01076           else chi2Tran = 100000;
01077 
01078           bothfitprob = 0.;
01079           Double_t chi2Both = 0.;
01080           
01081           MSG("SubShowerSR",Msg::kVerbose) << "chi2Tran " <<chi2Tran 
01082                                            << " chi2Long " << chi2Long << endl;
01083           
01084           chi2Both = chi2Tran+chi2Long;
01085           if((double(nfpl)/double(LongPl.size())>0.5&&norml>1e-3&&nlong>1)&&
01086              (double(nfpt)/double(TranStp.size())>0.5&&normt>1e-3&&ntran>1)&&
01087              chi2Both<=Minchi2Both){
01088             Minchi2Both = chi2Both;
01089             BestEBoth = Ein;
01090             BestThetaBoth = Thetain;
01091             bestnormtb = normt;
01092             bestnormlb = norml;
01093             fittedboth = true;
01094             if(makeplots){
01095               PlIt = LongPl.begin();
01096               PhIt = LongPh.begin();
01097               PluIt = LongFlu.begin();
01098               j = 0;
01099               while(PlIt!=LongPl.end()){
01100                 double data = (*PhIt)/totw*norml;
01101                 Int_t bin = LongProfb->FindBin((*PlIt)*100.);
01102                 LongProfb->SetBinContent(bin,funcl[j]);           
01103                 LongProfb->SetBinError(bin,errl[j]);
01104                 LongHistb->SetBinContent(bin,data);
01105                 LongHistb->SetBinError(bin,0.0);
01106                 PlIt++;
01107                 PhIt++;
01108                 PluIt++;
01109                 j++;
01110               }
01111               StpIt = TranStp.begin();
01112               StphIt = TranPh.begin();
01113               StpluIt = TranFlu.begin();
01114               j = 0;
01115 
01116               while(StpIt!=TranStp.end()){
01117                 double data = *StphIt/totw*normt;
01118                 Int_t bin = TranProfb->FindBin(*StpIt*100.);
01119                 TranProfb->SetBinContent(bin,funct[j]);
01120                 TranProfb->SetBinError(bin,errt[j]);
01121                 TranHistb->SetBinContent(bin,data);
01122                 TranHistb->SetBinError(bin,0.);
01123                 StpIt++;
01124                 StphIt++;
01125                 StpluIt++;
01126                 j++;
01127               }
01128             }
01129           }
01130         }
01131         Thetain+=dTheta;
01132       }
01133       Ein+=dE;
01134     }
01135 
01136     if(fittedboth){
01137       if(LongPl.size()>1&&TranStp.size()>1){
01138         Int_t ndfLong = LongPl.size();
01139         Int_t ndfTran = TranStp.size();
01140         Int_t ndfBoth = ndfLong + ndfTran;
01141         chi2func->SetParameter(0,ndfBoth);
01142         if(fittedboth && Minchi2Both/ndfBoth<100.)      
01143           bothfitprob = 1.0 - chi2func->Integral(0.0,Minchi2Both);
01144         else bothfitprob = 0.0;
01145         
01146         MSG("SubShowerSR",Msg::kDebug) 
01147           << "Best Fit:  Both- " << BestEBoth << " GeV " 
01148           << BestThetaBoth << " rad with chi2/ndof " 
01149           << Minchi2Both/ndfBoth << " ndf " << ndfBoth 
01150           << " of probability " << bothfitprob 
01151           << " comparing to " << Eguess << " GeV " 
01152           << Theta_rad << " rad" << endl;
01153 
01154       }
01155 
01156       if(makeplots){
01157         gStyle->SetOptStat(0);
01158         char name[256];
01159         TCanvas *can = new TCanvas();
01160         can->cd();
01161         can->Divide(2,1);
01162         can->cd(1);
01163         LongHistb->SetMaximum(bestnormlb);
01164         LongHistb->Draw("EP");
01165         LongProfb->SetMarkerColor(4);
01166         LongProfb->Draw("EPsame");
01167         can->cd(2);
01168         TranHistb->SetMaximum(bestnormtb);
01169         TranHistb->Draw("EP");
01170         TranProfb->SetMarkerColor(4);
01171         TranProfb->Draw("EPsame");
01172         if(bestnormlb+bestnormtb>0){
01173           sprintf(name,"Fit%f_%fb.eps",Eguess,Theta_rad);
01174           can->Print(name);
01175         }
01176         delete can;
01177       }
01178     }
01179   }
01180   else {
01181     bothfitprob = 0.;
01182   }
01183   MSG("SubShowerSR",Msg::kVerbose) << "long " << LongPl.size() 
01184                                    << " tran " <<TranStp.size()
01185                                    << " id " << id 
01186                                    << " eng " << cch.GetEnergy()
01187                                    << " nstrip " << cch.GetNStrip() 
01188                                    << " Prob " << bothfitprob << endl;
01189   cch.SetProbEM(Float_t(bothfitprob));
01190 }

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

Reimplemented from AlgBase.

Definition at line 1193 of file AlgSubShowerSR.cxx.

01194 {
01195 }


Member Data Documentation

Double_t AlgSubShowerSR::avgph [private]

Definition at line 52 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::avgstpperpln [private]

Definition at line 50 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::begz [private]

Definition at line 68 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::bothfitprob [private]

Definition at line 54 of file AlgSubShowerSR.h.

Referenced by SubShowerID().

TF1* AlgSubShowerSR::chi2func [private]

Definition at line 80 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

Definition at line 78 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

Definition at line 79 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

Double_t AlgSubShowerSR::Eguess [private]

Definition at line 72 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::endz [private]

Definition at line 69 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Bool_t AlgSubShowerSR::fDoEMFit [private]

Definition at line 43 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

Double_t AlgSubShowerSR::fECut [private]

Definition at line 35 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

TH1F* AlgSubShowerSR::flspl [private]

Definition at line 76 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and ~AlgSubShowerSR().

Double_t AlgSubShowerSR::fMaxDiffCut [private]

Definition at line 36 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

Double_t AlgSubShowerSR::fMaxXTalkStpPH [private]

Definition at line 42 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

Double_t AlgSubShowerSR::fMIPperGeV [private]

Definition at line 44 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and RunAlg().

Double_t AlgSubShowerSR::fOut2RmPHCut [private]

Definition at line 38 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

Double_t AlgSubShowerSR::fPECut [private]

Definition at line 34 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), RunAlg(), and SubShowerID().

Double_t AlgSubShowerSR::fTrkFraCut [private]

Definition at line 39 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

Double_t AlgSubShowerSR::fTrkLikeStpCut [private]

Definition at line 40 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

Double_t AlgSubShowerSR::fVtxDiffCut [private]

Definition at line 37 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

Double_t AlgSubShowerSR::fXtkFraCut [private]

Definition at line 41 of file AlgSubShowerSR.h.

Referenced by RunAlg(), and SubShowerID().

std::list<Double_t> AlgSubShowerSR::LongFlu [private]

Definition at line 60 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR().

TH1F* AlgSubShowerSR::LongHist [private]

Definition at line 83 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

TH1F* AlgSubShowerSR::LongHistb [private]

Definition at line 87 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

std::list<Double_t> AlgSubShowerSR::LongPh [private]

Definition at line 59 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR().

std::list<Double_t> AlgSubShowerSR::LongPhpe [private]

Definition at line 61 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and ~AlgSubShowerSR().

std::list<Double_t> AlgSubShowerSR::LongPl [private]

Definition at line 58 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR().

TH1F* AlgSubShowerSR::LongProf [private]

Definition at line 82 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

TH1F* AlgSubShowerSR::LongProfb [private]

Definition at line 86 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

Double_t AlgSubShowerSR::loweststp [private]

Definition at line 66 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::maxdiff [private]

Definition at line 47 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::maxStpPHinSS [private]

Definition at line 46 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Int_t AlgSubShowerSR::nn0l [private]

Definition at line 73 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Int_t AlgSubShowerSR::nn0t [private]

Definition at line 74 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::out2Rmph [private]

Definition at line 53 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

TH1F* AlgSubShowerSR::sflspl [private]

Definition at line 77 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and ~AlgSubShowerSR().

Double_t AlgSubShowerSR::stpdismean [private]

Definition at line 55 of file AlgSubShowerSR.h.

Double_t AlgSubShowerSR::stpdisrms [private]

Definition at line 56 of file AlgSubShowerSR.h.

Double_t AlgSubShowerSR::Theta_rad [private]

Definition at line 70 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::totw [private]

Definition at line 71 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

std::list<Double_t> AlgSubShowerSR::TranFlu [private]

Definition at line 64 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR().

TH1F* AlgSubShowerSR::TranHist [private]

Definition at line 85 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

TH1F* AlgSubShowerSR::TranHistb [private]

Definition at line 89 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

std::list<Double_t> AlgSubShowerSR::TranPh [private]

Definition at line 63 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR().

TH1F* AlgSubShowerSR::TranProf [private]

Definition at line 84 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

TH1F* AlgSubShowerSR::TranProfb [private]

Definition at line 88 of file AlgSubShowerSR.h.

Referenced by SubShowerID(), and ~AlgSubShowerSR().

std::list<Double_t> AlgSubShowerSR::TranStp [private]

Definition at line 62 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR().

Double_t AlgSubShowerSR::trkfra [private]

Definition at line 49 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::uppeststp [private]

Definition at line 67 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::vtxdiff [private]

Definition at line 48 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().

Double_t AlgSubShowerSR::xtkfra [private]

Definition at line 51 of file AlgSubShowerSR.h.

Referenced by CalculateEnergyVertexAngle(), and SubShowerID().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1