AlgShowerEM Class Reference

#include <AlgShowerEM.h>

Inheritance diagram for AlgShowerEM:
AlgBase AlgReco

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

Detailed Description

Definition at line 18 of file AlgShowerEM.h.


Constructor & Destructor Documentation

AlgShowerEM::AlgShowerEM (  ) 

Definition at line 44 of file AlgShowerEM.cxx.

00045 {
00046 }

AlgShowerEM::~AlgShowerEM (  )  [virtual]

Definition at line 49 of file AlgShowerEM.cxx.

00050 {
00051 }


Member Function Documentation

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

Read in info:

Implements AlgBase.

Definition at line 54 of file AlgShowerEM.cxx.

References AlgReco::Calibrate(), CandContext::GetDataIn(), Registry::GetDouble(), Registry::GetInt(), Calibrator::GetMIP(), CandHandle::GetVldContext(), Calibrator::Instance(), Msg::kDebug, Msg::kError, CalStripType::kGeV, CalStripType::kMIP, CalDigitType::kPE, CalDigitType::kSigCorr, Msg::kVerbose, Munits::m, MSG, and n.

00055 {
00056   
00057   MSG("ShowerEM", Msg::kDebug) << "Starting AlgShowerEM::RunAlg()" << endl;
00058   Calibrator& calibrator=Calibrator::Instance();
00059 
00060   assert(ch.InheritsFrom("CandShowerEMHandle"));
00061   CandShowerEMHandle &csh = dynamic_cast<CandShowerEMHandle &>(ch);
00062   assert(cx.GetDataIn());
00063   assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00064   const TObjArray *tary = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00065   
00066   Int_t nstp = 0;
00067   for (Int_t i=0; i<=tary->GetLast(); i++) {
00068     TObject *tobj = tary->At(i);
00069     assert(tobj->InheritsFrom("CandClusterHandle"));
00070     CandClusterHandle *clusterhandle = dynamic_cast<CandClusterHandle*>(tobj);
00071     nstp+=clusterhandle->GetNStrip();
00072   }
00073      
00074   if(nstp==0){
00075     MSG("ShowerEM", Msg::kError) 
00076       << "Found empty daughter list in CandShowerEM" << endl;
00077     csh.SetDirCosU(0.);
00078     csh.SetDirCosV(0.);
00079     csh.SetDirCosZ(1.);
00080     csh.SetVtxPlane(-1);
00081     csh.SetVtxZ(-1.);
00082     csh.SetVtxU(0.);
00083     csh.SetVtxV(0.);
00084     csh.SetVtxT(0.);
00085     csh.SetEnergy(0);
00086     csh.SetEigenVectors(0);
00087     csh.SetEigenValues(0);
00088     csh.SetOutPH(0);
00089     csh.SetAvgDev(0);
00090     csh.SetShwStatus(0);
00091     return;
00092   }
00093 
00094   Int_t ShwStatus = 1;
00095   Double_t fVertex[3] = {0};
00096   Double_t fVtxPlane = -1;
00097   Double_t fAngle[4] = {0};
00098 
00099   Double_t *fEigenVector = new Double_t[8];
00100   fEigenVector[0]=0; fEigenVector[1]=0; fEigenVector[2]=0; fEigenVector[3]=0;
00101   fEigenVector[4]=0; fEigenVector[5]=0; fEigenVector[6]=0; fEigenVector[7]=0;
00102   Double_t *fEigenValue = new Double_t[4];
00103   fEigenValue[0]=0; fEigenValue[1]=0; fEigenValue[2]=0; fEigenValue[3]=0;
00104   Double_t *fOutPH = new Double_t[5]; //0-Rad, 1-Lon, 2-Tmax, 3-MaxPl, 4-Redo
00105   fOutPH[0]=0; fOutPH[1]=0; fOutPH[2]=0; fOutPH[3]=0; fOutPH[4]=0;
00106   Double_t *fAvgDev = new Double_t[4];
00107   fAvgDev[0]=-1.; fAvgDev[1]=-1.; fAvgDev[2]=-1.; fAvgDev[3]=-1.;
00108 
00109   Double_t fMaxPlane = -1;
00110   Double_t fMaxStrip = -1;
00111   Double_t fMinStrip = 999;
00112   Double_t fMinTime = 999999999;
00113 
00114   Int_t fIter = ac.GetInt("NIter");   
00115   Double_t fCutOff = ac.GetDouble("CutOff");
00116   Double_t fShwFrac = ac.GetDouble("ShwFrac");
00117   Double_t fPECut = ac.GetDouble("PECut");
00118   Double_t fRadCut = ac.GetDouble("RadCut");
00119 
00120   Int_t fMaxBigPlaneSep = ac.GetInt("MaxBigPlaneSep");
00121   Int_t fIsolatedPlaneCut = ac.GetInt("IsolatedPlaneCut");
00122   Int_t fNGapPlanes = ac.GetInt("NGapPlanes");
00123   Int_t fNGapStrips = ac.GetInt("NGapStrips");
00124   Double_t fMaxAvgDev = ac.GetDouble("MaxAvgDev");
00125   Double_t fMaxOffFrac = ac.GetDouble("MaxOffFrac");
00126   Double_t fAvgDevCutForLargeOutliers = ac.GetDouble("AvgDevCutForLargeOutliers");
00127   Double_t fVtxBegPlaneDiffCut = ac.GetDouble("VtxBegPlaneDiffCut");
00128   Double_t fVtxPlaneDiffCut = ac.GetDouble("VtxPlaneDiffCut");
00129   Double_t fTrkFraCut = ac.GetDouble("TrkFraCut");
00130   Double_t fUVDiffCut = ac.GetDouble("UVDiffCut");
00131   Double_t mip2gev = ac.GetDouble("Mip2GeV");
00132 
00133   MSG("ShowerEM", Msg::kVerbose) 
00134     << "Configured Quantities: " 
00135     << "NIter = " << fIter << "\n"
00136     << "CutOff = " << fCutOff << "\n"
00137     << "ShwFrac = " << fShwFrac << "\n"
00138     << "PECut = " << fPECut << "\n"
00139     << "RadCut = " << fRadCut << "\n"    
00140     << "MaxBigPlaneSep = " << fMaxBigPlaneSep << "\n"
00141     << "IsolatedPlaneCut = " << fIsolatedPlaneCut << "\n"
00142     << "NGapPlanes = " << fNGapPlanes << "\n"
00143     << "NGapStrips = " << fNGapStrips << "\n"
00144     << "MaxAvgDev = " << fMaxAvgDev << "\n"
00145     << "MaxOffFrac = " << fMaxOffFrac << "\n"
00146     << "AvgDevCutForLargeOutliers = " << fAvgDevCutForLargeOutliers << "\n"
00147     << "VtxBegPlaneDiffCut = " << fVtxBegPlaneDiffCut << "\n"
00148     << "VtxPlaneDiffCut = " << fVtxPlaneDiffCut << "\n"
00149     << "TrkFraCut = " << fTrkFraCut << "\n"
00150     << "UVDiffCut = " << fUVDiffCut << "\n"
00151     << "Mip2GeV = " << mip2gev << "\n"
00152     << endl;
00153 
00154   //important hard coded values - don't change!
00155   Double_t Ec_eff = 21.6297;
00156   Double_t X0_eff = 0.0401542;
00157   Double_t d_eff = 0.06;
00158   Double_t d = .0595;
00159   Double_t t = .0411;
00160   
00162   Int_t *plane = new Int_t[nstp];
00163   Int_t *strip = new Int_t[nstp];
00164   Double_t *ph = new Double_t[nstp];
00165   Double_t *z = new Double_t[nstp];
00166   Double_t *tpos = new Double_t[nstp];
00167   Int_t *out = new Int_t[nstp];
00168   for(Int_t i=0;i<nstp;i++) out[i] = 0; 
00169   
00170   Float_t su[2][2] = {{0}};
00171   Float_t sv[2][2] = {{0}};
00172   Int_t woutstpu = -1;
00173   Int_t woutstpv = -1;
00174   Double_t totph = 0.;
00175 
00176   CandStripHandle *usefulStrip=0;
00177   Int_t icnt = 0;
00178   for (Int_t i=0; i<=tary->GetLast(); i++) {
00179     TObject *tobj = tary->At(i);
00180     assert(tobj->InheritsFrom("CandClusterHandle"));
00181     CandClusterHandle *clusterhandle = dynamic_cast<CandClusterHandle*>(tobj);
00182     CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator());
00183     while (CandStripHandle *stp = stripItr()) {
00184       plane[icnt] = stp->GetPlane();
00185       strip[icnt] = stp->GetStrip();
00186       ph[icnt] = calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
00187       z[icnt] = stp->GetZPos();
00188       tpos[icnt] = stp->GetTPos();
00189       if(stp->GetCharge(CalDigitType::kPE)<=fPECut) out[icnt] = 1;
00190       else totph += ph[icnt];
00191       usefulStrip = stp;
00192       icnt++;
00193     }
00194   }
00195   UgliGeomHandle ugh(*usefulStrip->GetVldContext());
00196   
00198   
00199   Int_t maxpl = -1;
00200   Double_t maxz = 0;
00201   Int_t secmaxpl = -1;
00202   Int_t maxstp = -1;
00203   Double_t maxtpos[2] = {0,0};
00204   Int_t maxplu = -1;
00205   Double_t maxzu = 0;
00206   Int_t maxstpu = 0;
00207   Double_t maxtposu = 0;
00208   Int_t begstpu = 0;
00209   //Double_t begtposu = 0;
00210   Int_t begplu = -1;
00211   Double_t begzu = 0;
00212 
00213   Int_t maxplv = -1;
00214   Double_t maxzv = 0;  
00215   Int_t maxstpv = 0;
00216   Double_t maxtposv = 0;
00217   Int_t begstpv = 0;
00218   //Double_t begtposv = 0;
00219   Int_t begplv = -1;
00220   Double_t begzv = 0;
00221   
00222   Float_t totwu = 0.;
00223   Float_t totwv = 0.;
00224   Float_t normposu = 0;
00225   Float_t normposv = 0;
00226   
00227   Double_t out2Rmphu = 1;
00228   Double_t out2Rmphv = 1;
00229   
00230   Double_t Thetau = 0;
00231   Double_t Theta_radu = 0;
00232   Double_t Thetav = 0;
00233   Double_t Theta_radv = 0;
00234   Double_t Thetau1 = 0;
00235   Double_t Theta_radu1 = 0;
00236   Double_t Thetav1 = 0;
00237   Double_t Theta_radv1 = 0;
00238 
00239   Double_t avgdevu = 0.;
00240   Double_t wavgdevu = 0.;
00241   Double_t avgdevv = 0.;
00242   Double_t wavgdevv = 0.;
00243   Double_t *evalu = new Double_t[2];
00244   Double_t *evalv = new Double_t[2];
00245   TVector2 *evu0 = new TVector2(0.,0.);
00246   TVector2 *evu1 = new TVector2(0.,0.);
00247   TVector2 *evv0 = new TVector2(0.,0.);
00248   TVector2 *evv1 = new TVector2(0.,0.);
00249   for (Int_t i = 0; i<2; i++) {
00250     evalu[i] = 0.;
00251     evalv[i] = 0.;
00252   }
00253   bool redou = false;
00254   bool redov = false;
00255   bool oncemoreu = true;
00256   bool oncemorev = true;
00257   bool remax = true;
00258   Int_t countu = 0;
00259   Int_t countv = 0;
00260   int Nu = 0;
00261   int Nv = 0;
00262   int Nu1 = 0;
00263   int Nv1 = 0;
00264   double trkfra = 0.;
00265   double uvdif = 0.;
00266   while(((out2Rmphu>fCutOff||redou||oncemoreu)&&countu<fIter)
00267         ||((out2Rmphv>fCutOff||redov||oncemorev)&&countv<fIter)){
00268     while (remax){
00269       remax = false;
00270       Double_t Eguess = totph*mip2gev;
00271       if(Eguess==0) {ShwStatus = 0; return;}
00272       Double_t tmax = log(Eguess*1000./Ec_eff)-0.858; //t_max formula
00273       Double_t tmaxpl = tmax*X0_eff/d_eff;
00274 
00275       //find max plane
00276       TH1F *flspl = new TH1F("flspl","flspl",500,-0.5,499.5);
00277       TH1F *sflspl = new TH1F("sflspl","sflspl",500,-0.5,499.5);
00278       for(Int_t i=0;i<nstp;i++){
00279         if (out[i]<=0) flspl->Fill(plane[i],ph[i]);
00280       }
00281       maxpl = Int_t(flspl->GetBinCenter(flspl->GetMaximumBin()));
00282       Int_t getz = 0;
00283       for(Int_t i=0;i<nstp;i++){
00284         if (out[i]<=0&&plane[i]!=maxpl) sflspl->Fill(plane[i],ph[i]);
00285         else if (getz<1&&out[i]<=0&&plane[i]==maxpl) {
00286           maxz = z[i];
00287           getz++;
00288         }
00289       }
00290       secmaxpl = Int_t(sflspl->GetBinCenter(sflspl->GetMaximumBin()));      
00291 
00292       if (abs(maxpl-secmaxpl)>fMaxBigPlaneSep*tmaxpl){
00293         Bool_t pos = true;
00294         Bool_t neg = true;
00295         for(Int_t i=1;i<=fIsolatedPlaneCut;i++){
00296           if(neg&&flspl->GetBinContent(flspl->FindBin(maxpl-i))<=0) neg=true;
00297           else neg=false;
00298           if(pos&&flspl->GetBinContent(flspl->FindBin(maxpl+i))<=0) pos=true;
00299           else pos=false;
00300         }
00301         if(pos||neg){
00302           for(Int_t i=0;i<nstp;i++){
00303             if (out[i]<=0&&plane[i]==maxpl) {
00304               out[i] = 4;
00305             }
00306           }
00307           remax = true; 
00308         }
00309         MSG("ShowerEM", Msg::kDebug)  << "Using second maximum: "<<maxpl
00310                                       <<" "<<secmaxpl<<" "<<tmaxpl<<endl;
00311       }
00312 
00313       delete flspl;
00314       delete sflspl;
00315 
00316       if (maxpl*secmaxpl<=0) {ShwStatus = 0; return;}
00317     }  
00318 
00319 
00320     //find boundary of primary cluster:
00321     //plane:
00322     TH1F *flspl = new TH1F("flspl","flspl",500,-0.5,499.5);
00323     for(Int_t i=0;i<nstp;i++){
00324       if (out[i]<=0) flspl->Fill(plane[i],ph[i]);
00325     }
00326     maxpl = Int_t(flspl->GetBinCenter(flspl->GetMaximumBin()));
00327     
00328     Int_t k = 0;
00329     Bool_t FoundGap = false;
00330     while(!FoundGap){
00331       Bool_t allEmpty = true;
00332       for(Int_t i=0;i<fNGapPlanes;i++){
00333         if(allEmpty&&flspl->GetBinContent(flspl->FindBin(maxpl-k-i))<=0) {
00334           allEmpty=true;
00335         }
00336         else allEmpty = false;
00337       }
00338       if(allEmpty) FoundGap=true;
00339       k++;
00340     }
00341     
00342     Int_t l = 0;
00343     FoundGap = false;
00344     while(!FoundGap){
00345       Bool_t allEmpty = true;
00346       for(Int_t i=0;i<fNGapPlanes;i++){
00347         if(allEmpty&&flspl->GetBinContent(flspl->FindBin(maxpl+l+i))<=0) {
00348           allEmpty=true;
00349         }
00350         else allEmpty = false;
00351       }
00352       if(allEmpty) FoundGap=true;
00353       l++;
00354     }
00355     
00356     for(Int_t i=0;i<nstp;i++){
00357       if (out[i]<=0&&((plane[i]<=(maxpl-k))||(plane[i]>=(maxpl+l)))) {
00358         out[i] = 2;
00359       }
00360     }
00361     delete flspl;
00362     //Int_t npl = l+k-1;
00363     
00364     //strip:
00365     Int_t stplb[2] = {200,200};
00366     Int_t stpub[2] = {0,0};
00367     Int_t maxstpa[2] = {0,0};
00368     for (Int_t j = 0; j<2; j++){
00369       double maxph[2] = {0,0};
00370       for (int pl = maxpl-k+1;pl<maxpl+l;pl++){
00371         TH1F *flsstp = new TH1F("flsstp","flsstp",201,-0.5,200.5);
00372         if(j==0&&((pl<249&&pl%2==0)||(pl>=249&&pl%2==1))){
00373           for(Int_t i=0;i<nstp;i++){
00374             if(plane[i]==pl&&out[i]<=0) flsstp->Fill(strip[i],ph[i]);
00375           }
00376         }
00377         else if (j==1&&((pl<249&&pl%2==1)||(pl>=249&&pl%2==0))){
00378           for(Int_t i=0;i<nstp;i++){
00379             if(plane[i]==pl&&out[i]<=0) flsstp->Fill(strip[i],ph[i]);
00380           }
00381         }
00382         maxstp = Int_t(flsstp->GetBinCenter(flsstp->GetMaximumBin()));
00383         if((j==0&&((pl<249&&pl%2==0)||(pl>=249&&pl%2==1)))||
00384            (j==1&&((pl<249&&pl%2==1)||(pl>=249&&pl%2==0)))){
00385           if (maxph[j] < flsstp->GetBinContent(flsstp->GetMaximumBin())){
00386           maxph[j] = flsstp->GetBinContent(flsstp->GetMaximumBin());
00387           maxstpa[j] = maxstp;
00388           }
00389         }
00390           
00391         Int_t m[2] = {0};
00392         FoundGap = false;
00393         while(!FoundGap){
00394           Bool_t allEmpty = true;
00395           for(Int_t i=0;i<fNGapStrips;i++){
00396             if(allEmpty&&flsstp->GetBinContent(flsstp->FindBin(maxstp-m[j]-i))<=0){
00397               allEmpty=true;
00398             }
00399             else allEmpty = false;
00400           }
00401           if(allEmpty) FoundGap=true;
00402           m[j]++;
00403         }
00404         
00405         Int_t n[2] = {0};
00406         FoundGap = false;
00407         while(!FoundGap){
00408           Bool_t allEmpty = true;
00409           for(Int_t i=0;i<fNGapStrips;i++){
00410             if(allEmpty&&flsstp->GetBinContent(flsstp->FindBin(maxstp+n[j]+i))<=0){
00411               allEmpty=true;
00412             }
00413             else allEmpty = false;
00414           }
00415           if(allEmpty) FoundGap=true;
00416           n[j]++;
00417         }
00418         for(Int_t i=0;i<nstp;i++){
00419           if((j==0&&((pl<249&&pl%2==0)||(pl>=249&&pl%2==1)))||
00420              (j==1&&((pl<249&&pl%2==1)||(pl>=249&&pl%2==0)))){
00421             if (plane[i]==pl&&out[i]<=0
00422                 &&((strip[i]<=(maxstp-m[j]))||(strip[i]>=(maxstp+n[j])))) {
00423               out[i] = 1;         
00424             }
00425           }
00426         }
00427         if (maxstp>0&&(maxstp-m[j])<stplb[j]) stplb[j] = maxstp-m[j];
00428         if (maxstp>0&&(maxstp+n[j])>stpub[j]) stpub[j] = maxstp+n[j];
00429         delete flsstp;
00430       }
00431     }
00432     int ctpl = 0;
00433     int ctstp[2] = {0,0};
00434     for(Int_t i=0;i<nstp;i++){
00435       if(out[i]<=0){
00436         if(plane[i]==maxpl&&ctpl<1){
00437           maxz = z[i];
00438           ctpl++;
00439         }
00440         for (int j=0; j<2; j++){
00441           if((j==0&&((plane[i]<249&&plane[i]%2==0)
00442                      ||(plane[i]>=249&&plane[i]%2==1)))||
00443              (j==1&&((plane[i]<249&&plane[i]%2==1)
00444                      ||(plane[i]>=249&&plane[i]%2==0)))){
00445             if(strip[i]==maxstpa[j]&&ctstp[j]<1){
00446               maxtpos[j] = tpos[i];
00447               ctstp[j]++;
00448             }
00449           }
00450         }
00451       }
00452     }
00453     //      MSG("ShowerEM", Msg::kDebug)  <<stplb[0]<<" "<<stplb[1]<<" "<<stpub[0]<<" "<<stpub[1]<<" "<<maxz<<" "<<maxtpos[0]<<" "<<maxtpos[1]<<endl;    
00454 
00455     //iterate over U strips
00456     while((out2Rmphu>fCutOff||redou||oncemoreu)&&countu<fIter){
00457       redou = false;
00458       oncemoreu = false;
00459       totwu = 0;
00460       countu++;
00461       if(woutstpu!=-1&&out[woutstpu]==0) out[woutstpu] = 1;
00462       else if (out[woutstpu]<0) out[woutstpu] = 5;
00463       TH1F *flsplu = new TH1F("flsplu","flsplu",500,-0.5,499.5);      
00464       TH2F *u2d = new TH2F("u2d","u2d",l+k+1,maxz-(k+0.5)*d,maxz+(l+0.5)*d,
00465                            stpub[0]-stplb[0]+1,maxtpos[0]-(maxstpa[0]-stplb[0]+0.5)*t,
00466                            maxtpos[0]+(stpub[0]-maxstpa[0]+0.5)*t);
00467       for(Int_t i=0;i<nstp;i++){
00468 
00469         if(((plane[i]<249&&plane[i]%2==0)
00470             ||(plane[i]>=249&&plane[i]%2==1))&&out[i]<=0){
00471           flsplu->Fill(plane[i],ph[i]);
00472           u2d->Fill(z[i],tpos[i],ph[i]);
00473         }
00474       }
00475       maxplu = Int_t(flsplu->GetBinCenter(flsplu->GetMaximumBin()));
00476       Int_t nplu = 0;
00477       for (Int_t pl = maxpl-k+1;pl<maxpl+l;pl++){
00478         if (flsplu->GetBinContent(flsplu->FindBin(pl))>0) nplu++;
00479       }
00480       TProfile *up = u2d->ProfileX();
00481       TF1 *p1 = new TF1("p1","pol1");
00482       up->Fit("p1","NQ");
00483       Double_t slopeu = p1->GetParameter(1);
00484       //Double_t eslopeu = p1->GetParError(1);
00485       Double_t chi2u = p1->GetChisquare();
00486       Int_t j = 0;
00487       Double_t sumb = 0.;
00488       while(sumb==0&&j<=500) {
00489         j++;
00490         if(flsplu->GetBinContent(j)>0) {
00491           begplu = Int_t(flsplu->GetBinCenter(j));
00492           sumb += flsplu->GetBinContent(j);
00493         }
00494       }
00495       
00496       delete flsplu;
00497       delete u2d;
00498       delete up;
00499       delete p1;
00500       Float_t maxstpphu = 0.;
00501       Float_t begstpphu = 0.;
00502       Int_t nstpmaxplu = 0;
00503       for(Int_t i=0;i<nstp;i++){
00504         if(out[i]<=0){
00505           if(plane[i]==maxplu&&ph[i]>maxstpphu){
00506             maxstpu = strip[i];
00507             maxstpphu = ph[i];
00508             maxtposu = tpos[i];
00509           }
00510           if(plane[i]==maxplu){
00511             maxzu = z[i];
00512             nstpmaxplu++;
00513           }
00514           if(plane[i]==begplu&&ph[i]>begstpphu){
00515             begzu = z[i];
00516             begstpu = strip[i];
00517             begstpphu = ph[i];
00518           }
00519         }
00520       }
00521       
00522       Double_t cntu = 0;
00523       for(Int_t i=0;i<nstp;i++){
00524         if(out[i]<=0){  
00525           if((plane[i]<249&&plane[i]%2==0)||(plane[i]>=249&&plane[i]%2==1)){
00526             
00527             Double_t dtpos = tpos[i] - maxtposu;
00528             Double_t dz = z[i] - maxzu;
00529             su[0][0]+=dz*dz*ph[i];
00530             su[0][1]+=dtpos*dz*ph[i];
00531             su[1][1]+=dtpos*dtpos*ph[i];
00532             totwu += ph[i];
00533             cntu+=1;
00534             normposu += dtpos*dtpos+dz*dz;
00535             
00536           }
00537         }
00538       }
00539       su[1][0]=su[0][1];
00540       
00541       if(totwu*normposu!=0){
00542         
00543         TMatrixDSym mu(2);
00544         Double_t *angu = new Double_t[2];
00545         
00546         for(Int_t i=0;i<2;i++) angu[i] = -999;
00547         for (Int_t i=0;i<2;i++){
00548           for (Int_t j=0;j<2;j++){
00549             mu(i,j)=su[i][j]/totwu/normposu;
00550           }
00551         }
00552         
00553         if(mu.IsA()){
00554           TMatrixDSymEigen eigen(mu);
00555           TMatrixD EVect1 = eigen.GetEigenVectors();      
00556           TVectorD EVal1 = eigen.GetEigenValues();
00557           evu0 = new TVector2(EVect1(0,0),EVect1(1,0));
00558           evu1 = new TVector2(EVect1(0,1),EVect1(1,1));
00559           angu[0] = evu0->Phi()*180/TMath::Pi();
00560           angu[1] = evu1->Phi()*180/TMath::Pi();
00561           evalu[0] = EVal1(0);
00562           evalu[1] = EVal1(1);
00563         }
00564         else {
00565           evalu[0] = 0.;
00566           evalu[1] = 0.;
00567         }
00568         for (Int_t i=0;i<2;i++){
00569           if (angu[i]>90&&angu[i]<270) angu[i]-=180;
00570           if (angu[i]>270&&angu[i]<360) angu[i]-=360;
00571           if (angu[i]==90||angu[i]==270) angu[i]-=0.01;
00572         }
00573 
00574         Thetau = angu[0];
00575         Thetau1 = angu[1];
00576         Theta_radu = (Thetau/180)*TMath::Pi();
00577         Theta_radu1 = (Thetau1/180)*TMath::Pi();
00578         Double_t *shwstpc = new Double_t[nstp];
00579         Double_t *shwstpc1 = new Double_t[nstp];
00580         Double_t *shwstpc2 = new Double_t[nstp];
00581         Double_t *wshwstpc = new Double_t[nstp];
00582         for(Int_t i=0;i<nstp;i++) shwstpc[i] = -1.;
00583         for(Int_t i=0;i<nstp;i++) shwstpc1[i] = -1.;
00584         for(Int_t i=0;i<nstp;i++) shwstpc2[i] = -1.;
00585         for(Int_t i=0;i<nstp;i++) wshwstpc[i] = -1.;
00586         
00587         Double_t wdevmaxu = 0.;
00588         woutstpu = -1;
00589         avgdevu = 0.;
00590         wavgdevu = 0.;
00591         Double_t avgdevu1 = 0.;
00592         Double_t avgdevu2 = 0.;
00593         Double_t wavgdevu1 = 0.;
00594         Double_t wavgdevu2 = 0.;
00595         Int_t instpu = 0;
00596         //Int_t instpu1 = 0;
00597         Double_t winstpu = 0;
00598         Float_t offstpu = 0;
00599         Float_t offstpu2 = 0;
00600         //Float_t offstpu3 = 0;
00601         out2Rmphu = 0;
00602         for(Int_t i=0;i<nstp;i++){
00603           if(out[i]<=0){
00604             if((plane[i]<249&&plane[i]%2==0)||(plane[i]>=249&&plane[i]%2==1)){
00605               shwstpc[i]=fabs(tpos[i]
00606                               -(maxtposu+tan(Theta_radu)
00607                                 *(z[i]-maxzu)))*TMath::Cos(Theta_radu);
00608               shwstpc2[i]=fabs(tpos[i]
00609                                -(maxtposu+tan(Theta_radu1)
00610                                  *(z[i]-maxzu)))*TMath::Cos(Theta_radu1);
00611               avgdevu1 += shwstpc[i];
00612               avgdevu2 += shwstpc2[i];
00613               wavgdevu1 += shwstpc[i]*ph[i];
00614               wavgdevu2 += shwstpc2[i]*ph[i];
00615               instpu++;
00616               winstpu+=ph[i];
00617               if (shwstpc[i]>1.*t) offstpu++;
00618               if (shwstpc2[i]>1.*t) offstpu2++;
00619               //              if (shwstpc2[i]>6.*t) offstpu3++;
00620             }
00621           }
00622         }
00623         avgdevu1/=instpu;
00624         avgdevu2/=instpu;
00625         wavgdevu1/=winstpu;
00626         wavgdevu2/=winstpu;
00627         avgdevu = avgdevu1;
00628         wavgdevu = wavgdevu1;
00629         offstpu/=instpu;
00630         offstpu2/=instpu;
00631         for(Int_t i=0;i<nstp;i++){
00632           if((fabs(Thetau)>fabs(Thetau1))
00633              &&((avgdevu1>fMaxAvgDev*t&&avgdevu2>fMaxAvgDev*t)
00634                 ||(offstpu>fMaxOffFrac&&offstpu2>fMaxOffFrac))
00635              &&strip[i]==maxstpu&&nstpmaxplu==1) {
00636             out[i] = 4;
00637             oncemoreu = true;
00638             MSG("ShowerEM", Msg::kDebug) <<"u maxstp removed"<<endl;
00639           }
00640         }
00641         if(fabs(tan(Theta_radu)-slopeu)>fabs(tan(Theta_radu1)-slopeu)
00642            &&(wavgdevu1>wavgdevu2/1000.||nplu==1)){
00643           Thetau = angu[1];      
00644           Thetau1 = angu[0];
00645           avgdevu = avgdevu2;
00646           wavgdevu = wavgdevu2;
00647           for(Int_t i=0;i<nstp;i++){
00648             shwstpc[i] = shwstpc2[i];
00649           }
00650         }
00651         
00652         Theta_radu = (Thetau/180)*TMath::Pi();
00653         Theta_radu1 = (Thetau1/180)*TMath::Pi();
00654         for(Int_t i=0;i<nstp;i++){
00655           if(out[i]<=0){
00656             if((plane[i]<249&&plane[i]%2==0)||(plane[i]>=249&&plane[i]%2==1)){
00657               wshwstpc[i] = shwstpc[i]*ph[i];
00658               if(shwstpc[i]>fRadCut*t) {
00659                 out2Rmphu += ph[i];
00660                 if(wshwstpc[i]>wdevmaxu){ 
00661                   wdevmaxu = wshwstpc[i];
00662                   woutstpu = i;         }
00663               }
00664             }
00665           }
00666           else if(out[i]==1){
00667             if((plane[i]<249&&plane[i]%2==0)||(plane[i]>=249&&plane[i]%2==1)){
00668               shwstpc[i] = fabs(tpos[i]-(maxtposu+tan(Theta_radu)*(z[i]-maxzu)))*TMath::Cos(Theta_radu);;
00669               wshwstpc[i] = shwstpc[i]*ph[i];
00670               if(shwstpc[i]<=fRadCut*t) {
00671                 out[i] = -1;
00672                 redou = true;
00673               }
00674             }
00675           }
00676         }
00677         MSG("ShowerEM", Msg::kDebug)  <<"avgu "<< avgdevu <<" "<< avgdevu1 
00678                                       <<" "<<avgdevu2<<" "<<wavgdevu1 <<" "
00679                                       <<wavgdevu2<<" "<<offstpu<<" "
00680                                       <<offstpu2<<" "<<tan(Theta_radu)<<" "
00681                                       <<slopeu<<" "<<chi2u<<endl;
00682         if (avgdevu>fAvgDevCutForLargeOutliers*t){
00683           for(Int_t i=0;i<nstp;i++){
00684             if(((plane[i]<249&&plane[i]%2==0)
00685                 ||(plane[i]>=249&&plane[i]%2==1))&&out[i]<=0){
00686               if((shwstpc[i]-avgdevu)>fRadCut*t) {
00687                 if (out[i]==0) out[i] = 1;  
00688                 else if (out[i]<0) out[i] = 5;
00689                 oncemoreu = true;
00690               }
00691             }
00692           }
00693         }
00694         out2Rmphu /= totwu;
00695         Nu = 0;
00696         Nu1 = 0;
00697         for (int pl = maxpl-k+1;pl<maxpl+l;pl++){
00698           int nstppl = 0;
00699           for(int i=0;i<nstp;i++){
00700             if(plane[i]==pl&&((plane[i]<249&&plane[i]%2==0)
00701                               ||(plane[i]>=249&&plane[i]%2==1))
00702                &&out[i]<=0){
00703               nstppl++;
00704             }
00705           }
00706           if (nstppl>0) Nu++; 
00707           if (nstppl==1) Nu1++; 
00708         }
00709       }
00710     }
00711     
00712     //iterate over V strips    
00713     while((out2Rmphv>fCutOff||redov||oncemorev)&&countv<fIter){
00714       redov = false;
00715       oncemorev = false;
00716       totwv = 0;
00717       countv++;
00718       if(woutstpv!=-1&&out[woutstpv]==0) out[woutstpv] = 1;
00719       else if (out[woutstpv]<0) out[woutstpv] = 5;
00720       TH1F *flsplv = new TH1F("flsplv","flsplv",500,-0.5,499.5);
00721       TH2F *v2d = new TH2F("v2d","v2d",l+k+1,maxz-(k+0.5)*d,maxz+(l+0.5)*d,
00722                            stpub[1]-stplb[1]+1,maxtpos[1]-(maxstpa[1]-stplb[1]+0.5)*t,
00723                            maxtpos[1]+(stpub[1]-maxstpa[1]+0.5)*t);
00724       for(Int_t i=0;i<nstp;i++){
00725         if(((plane[i]>=249&&plane[i]%2==0)
00726             ||(plane[i]<249&&plane[i]%2==1))&&out[i]<=0){
00727           flsplv->Fill(plane[i],ph[i]);
00728           v2d->Fill(z[i],tpos[i],ph[i]);
00729         }
00730       }
00731       maxplv = Int_t(flsplv->GetBinCenter(flsplv->GetMaximumBin()));
00732       Int_t nplv = 0;
00733       for (Int_t pl = maxpl-k+1;pl<maxpl+l;pl++){
00734         if (flsplv->GetBinContent(flsplv->FindBin(pl))>0) nplv++;
00735       }
00736       TProfile *vp = v2d->ProfileX();
00737       TF1 *p1 = new TF1("p1","pol1");      
00738       vp->Fit("p1","NQ");
00739       Double_t slopev = p1->GetParameter(1);      
00740       //Double_t eslopev = p1->GetParError(1);      
00741       Double_t chi2v = p1->GetChisquare();
00742       Int_t j = 0;
00743       Double_t sumb = 0.; 
00744       while (sumb==0&&j<=500) {
00745         j++;
00746         if(flsplv->GetBinContent(j)>0) {
00747           begplv = Int_t(flsplv->GetBinCenter(j));
00748           sumb += flsplv->GetBinContent(j);
00749         }
00750       }
00751       delete flsplv;
00752       delete v2d;
00753       delete vp;
00754       delete p1;
00755       Float_t maxstpphv = 0.;
00756       Float_t begstpphv = 0.;
00757       Int_t nstpmaxplv = 0;
00758       for(Int_t i=0;i<nstp;i++){
00759         if(out[i]<=0){
00760           if(plane[i]==maxplv&&ph[i]>maxstpphv){
00761             maxstpv = strip[i];
00762             maxstpphv = ph[i];
00763             maxtposv = tpos[i];
00764           }
00765           if(plane[i]==maxplv){
00766             maxzv = z[i];
00767             nstpmaxplv++;
00768           }
00769           if (plane[i]==begplv&&ph[i]>begstpphv){
00770             begzv = z[i];
00771             begstpv = strip[i];
00772             begstpphv = ph[i];
00773           }
00774         }
00775       }
00776       
00777       Double_t cntv = 0;
00778       for(Int_t i=0;i<nstp;i++){
00779         if(out[i]<=0){
00780           if ((plane[i]>=249&&plane[i]%2==0)||(plane[i]<249&&plane[i]%2==1)){
00781             Double_t dtpos = tpos[i] - maxtposv;
00782             Double_t dz = z[i] - maxzv;
00783             sv[0][0]+=dz*dz*ph[i];
00784             sv[0][1]+=dtpos*dz*ph[i];
00785             sv[1][1]+=dtpos*dtpos*ph[i];
00786             totwv += ph[i];
00787             cntv+=1;
00788             normposv += dtpos*dtpos+dz*dz;
00789           }
00790         }
00791       }
00792       sv[1][0]=sv[0][1];
00793 
00794       if(totwv*normposv!=0){
00795         TMatrixDSym mv(2);
00796         Double_t *angv = new Double_t[2];
00797         
00798         for(Int_t i=0;i<2;i++) angv[i] = -999;
00799         for (Int_t i=0;i<2;i++){
00800           for (Int_t j=0;j<2;j++){
00801             mv(i,j)=sv[i][j]/totwv/normposv;
00802           }
00803         }
00804         
00805         if(mv.IsA()){
00806           TMatrixDSymEigen eigen(mv);
00807           TMatrixD EVect2 = eigen.GetEigenVectors();
00808           TVectorD EVal2 = eigen.GetEigenValues();      
00809           evv0 = new TVector2(EVect2(0,0),EVect2(1,0));
00810           evv1 = new TVector2(EVect2(0,1),EVect2(1,1));
00811           angv[0] = evv0->Phi()*180/TMath::Pi();
00812           angv[1] = evv1->Phi()*180/TMath::Pi();
00813           evalv[0] = EVal2(0);
00814           evalv[1] = EVal2(1);
00815         }
00816         else {
00817           evalv[0] = 0.;
00818           evalv[1] = 0.;
00819         }   
00820         for (Int_t i=0;i<2;i++){
00821           if (angv[i]>90&&angv[i]<270) angv[i]-=180;
00822           if (angv[i]>270&&angv[i]<360) angv[i]-=360;
00823           if (angv[i]==90||angv[i]==270) angv[i]-=0.01;
00824         }
00825         
00826         Thetav = angv[0];
00827         Thetav1 = angv[1];
00828         Theta_radv = (Thetav/180)*TMath::Pi();
00829         Theta_radv1 = (Thetav1/180)*TMath::Pi();
00830         Double_t *shwstpc = new Double_t[nstp];
00831         Double_t *shwstpc1 = new Double_t[nstp];
00832         Double_t *shwstpc2 = new Double_t[nstp];
00833         Double_t *wshwstpc = new Double_t[nstp];
00834         for(Int_t i=0;i<nstp;i++) shwstpc[i] = -1.;
00835         for(Int_t i=0;i<nstp;i++) shwstpc1[i] = -1.;
00836         for(Int_t i=0;i<nstp;i++) shwstpc2[i] = -1.;
00837         for(Int_t i=0;i<nstp;i++) wshwstpc[i] = -1.;
00838         
00839         Double_t wdevmaxv = 0.;
00840         woutstpv = -1;
00841         avgdevv = 0.;
00842         wavgdevv = 0.;
00843         Double_t avgdevv1 = 0.;
00844         Double_t avgdevv2 = 0.;
00845         Double_t wavgdevv1 = 0.;
00846         Double_t wavgdevv2 = 0.;
00847         Int_t instpv = 0;
00848         //Int_t instpv1 = 0;
00849         Double_t winstpv = 0;
00850         Float_t offstpv = 0;
00851         Float_t offstpv2 = 0;
00852         //      Float_t offstpv3 = 0;
00853         out2Rmphv = 0;
00854         for(Int_t i=0;i<nstp;i++){
00855           if(out[i]<=0){
00856             if((plane[i]<249&&plane[i]%2==1)||(plane[i]>=249&&plane[i]%2==0)){
00857               shwstpc[i]=fabs(tpos[i]
00858                               -(maxtposv+tan(Theta_radv)
00859                                 *(z[i]-maxzv)))*TMath::Cos(Theta_radv);
00860               shwstpc2[i]=fabs(tpos[i]
00861                                -(maxtposv+tan(Theta_radv1)
00862                                  *(z[i]-maxzv)))*TMath::Cos(Theta_radv1);
00863               avgdevv1 += shwstpc[i];
00864               avgdevv2 += shwstpc2[i];
00865               wavgdevv1 += shwstpc[i]*ph[i];
00866               wavgdevv2 += shwstpc2[i]*ph[i];
00867               instpv++;
00868               winstpv+=ph[i];
00869               if (shwstpc[i]>1) offstpv++;
00870               if (shwstpc2[i]>1) offstpv2++;
00871               //              if (shwstpc2[i]>6) offstpv3++;
00872             }
00873           }
00874         }
00875         avgdevv1/=instpv;
00876         avgdevv2/=instpv;
00877         wavgdevv1/=winstpv;
00878         wavgdevv2/=winstpv;
00879         avgdevv = avgdevv1;
00880         wavgdevv = wavgdevv1;
00881         offstpv/=instpv;
00882         offstpv2/=instpv;
00883         for(Int_t i=0;i<nstp;i++){
00884           if ((fabs(Thetav)>fabs(Thetav1))
00885               &&((avgdevv1>fMaxAvgDev*t&&avgdevv2>fMaxAvgDev*t)
00886                  ||(offstpv>fMaxOffFrac&&offstpv2>fMaxOffFrac))
00887               &&strip[i]==maxstpv&&nstpmaxplv==1) {
00888             out[i] = 4;
00889             oncemorev = true;
00890             MSG("ShowerEM", Msg::kDebug) <<"v maxstp removed"<<endl;
00891           }
00892         }
00893         if(fabs(tan(Theta_radv)-slopev)>fabs(tan(Theta_radv1)-slopev)
00894            &&(wavgdevv1>wavgdevv2/1000.||nplv==1)){
00895           Thetav = angv[1];      
00896           Thetav1 = angv[0];
00897           avgdevv = avgdevv2;
00898           wavgdevv = wavgdevv2;
00899           for(Int_t i=0;i<nstp;i++){
00900             shwstpc[i] = shwstpc2[i];
00901           }
00902         }
00903         
00904         Theta_radv = (Thetav/180)*TMath::Pi();
00905         Theta_radv1 = (Thetav1/180)*TMath::Pi();
00906 
00907         for(Int_t i=0;i<nstp;i++){
00908           if(out[i]<=0){
00909             if((plane[i]<249&&plane[i]%2==1)||(plane[i]>=249&&plane[i]%2==0)){
00910               wshwstpc[i] = shwstpc[i]*ph[i];
00911               if(shwstpc[i]>fRadCut*t) {
00912                 out2Rmphv += ph[i];
00913                 if(wshwstpc[i]>wdevmaxv){ 
00914                   wdevmaxv = wshwstpc[i];
00915                   woutstpv = i;         }
00916               }
00917             }
00918           }
00919           else if(out[i]==1){
00920             if((plane[i]<249&&plane[i]%2==1)||(plane[i]>=249&&plane[i]%2==0)){
00921               shwstpc[i] = fabs(tpos[i]-(maxtposv+tan(Theta_radv)*(z[i]-maxzv)))*TMath::Cos(Theta_radv);;
00922               wshwstpc[i] = shwstpc[i]*ph[i];
00923               if(shwstpc[i]<=fRadCut*t) {
00924                 out[i] = -1;
00925                 redov = true;
00926               }
00927             }
00928           }
00929         }
00930         MSG("ShowerEM", Msg::kDebug)  <<"avgv "<< avgdevv <<" "<<avgdevv1 
00931                                       <<" "<<avgdevv2<<" "<<wavgdevv1 <<" "
00932                                       <<wavgdevv2<<" "<<offstpv<<" "
00933                                       <<offstpv2<<" "<<tan(Theta_radv)<<" "
00934                                       <<slopev<<" "<<chi2v<<endl;
00935         if (avgdevv>fAvgDevCutForLargeOutliers*t){
00936           for(Int_t i=0;i<nstp;i++){
00937             if(((plane[i]<249&&plane[i]%2==1)
00938                 ||(plane[i]>=249&&plane[i]%2==0))&&out[i]<=0){
00939               if((shwstpc[i]-avgdevv)>fRadCut*t) {
00940                 if (out[i]==0) out[i] = 1;  
00941                 else if (out[i]<0) out[i] = 5;
00942                 oncemorev = true;
00943               }
00944             }
00945           }
00946         }
00947         out2Rmphv /= totwv;
00948         Nv = 0;
00949         Nv1 = 0;
00950         for (int pl = maxpl-k+1;pl<maxpl+l;pl++){
00951           int nstppl = 0;
00952           for(int i=0;i<nstp;i++){
00953             if(plane[i]==pl&&((plane[i]<249&&plane[i]%2==1)
00954                               ||(plane[i]>=249&&plane[i]%2==0))
00955                &&out[i]<=0){
00956               nstppl++;
00957             }
00958           }
00959           if (nstppl>0) Nv++; 
00960           if (nstppl==1) Nv1++; 
00961         }       
00962       }
00963     }
00964     if ((Nu+Nv)>0){
00965       trkfra = double(Nu1+Nv1)/double(Nu+Nv);
00966       uvdif = double(abs(Nu-Nv))/double(Nu+Nv);
00967     }
00968     else {
00969       trkfra = -1;
00970       uvdif = -1;
00971     }
00972     Double_t shwfra = (totwu+totwv)/totph;
00973     
00974     if(totwu*totwv==0||shwfra<fShwFrac||(countu==fIter&&out2Rmphu>fCutOff)
00975        ||(countv==fIter&&out2Rmphv>fCutOff)) ShwStatus = 0;
00976     if(ShwStatus==0) return;
00977 
00978     if (trkfra>fTrkFraCut) ShwStatus = 2;
00979     if (uvdif>fUVDiffCut) ShwStatus = 3;  
00980   
00981     Double_t Eguess = (totwu+totwv)*mip2gev;
00982     if(Eguess==0) {ShwStatus = 0; return;}
00983     Double_t tmax = log(Eguess*1000./Ec_eff)-0.858;
00984     Double_t tmaxz = tmax*X0_eff;
00985     
00986     Double_t begz = begzv;
00987     Double_t begpl = begplv;
00988     if(begzu<begzv) {
00989       begz = begzu;
00990       begpl = begplu;
00991     }
00992     Double_t vtxz = maxz-tmaxz/sqrt(TMath::Power(tan(Theta_radu),2)+TMath::Power(tan(Theta_radv),2)+1.);
00993 
00994     MSG("ShowerEM", Msg::kDebug) <<shwfra<<" "<<Eguess<<" "<<tmaxz<<" "
00995                                  <<vtxz<<" "<<begz<<" "<<maxz<<endl;
00996     
00997     if ((vtxz-begz)<=fVtxBegPlaneDiffCut*d) {
00998       fVertex[2] = begz;
00999       fVtxPlane = begpl;
01000     }
01001     else {
01002       fVertex[2] = vtxz;
01003       fVtxPlane = begpl+(vtxz-begz)/d;
01004       for(Int_t i=0;i<nstp;i++){
01005         if(((z[i]-vtxz)<-fVtxPlaneDiffCut*d)
01006            &&((plane[i]<249&&plane[i]%2==0)
01007               ||(plane[i]>=249&&plane[i]%2==1))
01008            &&out[i]<=0) {
01009           out[i] = 3;
01010           oncemoreu = true;
01011         }
01012         else if(((z[i]-vtxz)<-fVtxPlaneDiffCut*d)
01013                 &&((plane[i]<249&&plane[i]%2==1)
01014                    ||(plane[i]>=249&&plane[i]%2==0))
01015                 &&out[i]<=0) {
01016           out[i] = 3;
01017           oncemorev = true;
01018         }
01019       }
01020     }
01021   }
01022 
01023   fVertex[0] = maxtposu+tan(Theta_radu)*(fVertex[2]-maxzu);
01024   fVertex[1] = maxtposv+tan(Theta_radv)*(fVertex[2]-maxzv);
01025   fAngle[0] = Theta_radu;
01026   fAngle[1] = Theta_radv;
01027   fAngle[2] = Theta_radu1;
01028   fAngle[3] = Theta_radv1;
01029   fEigenValue[0] = evalu[0];
01030   fEigenValue[1] = evalv[0];
01031   fEigenValue[2] = evalu[1];
01032   fEigenValue[3] = evalv[1];
01033   fAvgDev[0] = avgdevu;
01034   fAvgDev[1] = avgdevv;
01035   fAvgDev[2] = wavgdevu;
01036   fAvgDev[3] = wavgdevv;
01037   fEigenVector[0] = evu0->X();
01038   fEigenVector[1] = evu0->Y();
01039   fEigenVector[2] = evv0->X();
01040   fEigenVector[3] = evv0->Y();
01041   fEigenVector[4] = evu1->X();
01042   fEigenVector[5] = evu1->Y();
01043   fEigenVector[6] = evv1->X();
01044   fEigenVector[7] = evv1->Y();
01045 
01046   icnt = 0;
01047   Int_t jcnt = 0;
01048 
01049   for (Int_t i=0; i<=tary->GetLast(); i++) {
01050     TObject *tobj = tary->At(i);
01051     assert(tobj->InheritsFrom("CandClusterHandle"));
01052     CandClusterHandle *clusterhandle = dynamic_cast<CandClusterHandle*>(tobj);
01053     CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator());
01054     while (CandStripHandle *stp = stripItr()) {
01055       if(out[icnt]<=0){
01056         //get min/max number of strips
01057         if(stp->GetStrip()>fMaxStrip) fMaxStrip = stp->GetStrip();
01058         if(stp->GetStrip()<fMinStrip) fMinStrip = stp->GetStrip();
01059         if(stp->GetPlane()>fMaxPlane) fMaxPlane = stp->GetPlane();
01060         //find minimum strip time at vertex plane
01061         PlexStripEndId stripid = stp->GetStripEndId();
01062         if(stp->GetPlane()==int(fVtxPlane+0.5)){
01063           if(stp->GetTime()<fMinTime) fMinTime = stp->GetTime();
01064         }
01065         csh.AddDaughterLink(*stp);
01066         jcnt++;
01067       }
01068       else{
01069         fOutPH[out[icnt]-1] += calibrator.
01070           GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
01071       }
01072       icnt++;
01073     }
01074     csh.AddCluster(clusterhandle);
01075   }
01076   
01077   csh.SetVtxPlane(int(fVtxPlane+0.5));
01078   csh.SetVtxZ(fVertex[2]);  
01079   csh.SetVtxU(fVertex[0]);
01080   csh.SetVtxV(fVertex[1]);
01081   csh.SetVtxT(fMinTime);
01082 
01083   double dudz = TMath::Tan(fAngle[0]);
01084   double dvdz = TMath::Tan(fAngle[1]);
01085   double dsdz = sqrt(dudz*dudz + dvdz*dvdz + 1.);
01086   csh.SetDirCosU(dudz/dsdz);
01087   csh.SetDirCosV(dvdz/dsdz);
01088   csh.SetDirCosZ(1./dsdz);
01089 
01090   csh.SetEigenVectors(fEigenVector);
01091   csh.SetEigenValues(fEigenValue);
01092   csh.SetAvgDev(fAvgDev);  
01093   csh.SetOutPH(fOutPH);
01094   csh.SetShwStatus(ShwStatus);
01095 
01096   //SetUV(&csh);
01097   Calibrate(&csh);
01098 
01099   Double_t TotMip = csh.GetCharge(CalStripType::kMIP);
01100   Double_t TotEnergy = csh.GetCharge(CalStripType::kGeV);
01101 
01102   MSG("ShowerEM", Msg::kVerbose) << endl
01103                                  << " TotMIP = " << TotMip
01104                                  << endl << "Setting shower energy to " 
01105                                  << TotEnergy << endl;
01106   
01107   csh.SetEnergy(TotEnergy);
01108   
01109 }

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

Reimplemented from AlgBase.

Definition at line 1112 of file AlgShowerEM.cxx.

01113 {
01114 }


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1