loon: AlgTrackSR.cxx Source File

00001 
00002 // $Id: AlgTrackSR.cxx,v 1.71 2007/02/04 06:03:24 rhatcher Exp $
00003 //
00004 // AlgTrackSR.cxx
00005 //
00006 // Begin_Html<img src="../../pedestrians.gif" align=center>
00007 // <a href="../source_warning.html">Warning for beginners</a>.<br> 
00008 //
00009 // AlgTrackSR is a concrete TrackSR Algorithm class.
00010 //
00011 // Author:  R. Lee 2001.02.26
00012 // Revised: B. Rebel 2004.03.31
00013 //
00014 // Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and 
00015 // <a href="../CandTrackSR.html"> CandTrackSR Classes</a> (part of
00016 // <a href="../index.html">The MINOS Class User Guide</a>)End_Html
00018 
00019 #include <cassert>
00020 
00021 #include "Algorithm/AlgFactory.h"
00022 #include "Algorithm/AlgHandle.h"
00023 #include "Algorithm/AlgConfig.h"
00024 #include "Candidate/CandContext.h"
00025 //#include "./AlgTrackSR.h"
00026 #include "CandTrackSR/AlgTrackSR.h"
00027 #include "CandTrackSR/CandTrackSRHandle.h"
00028 #include "CandTrackSR/TrackClusterSR.h"
00029 #include "CandTrackSR/Track2DSR.h"
00030 #include "Conventions/CalTimeType.h"
00031 #include "Conventions/Mphysical.h"
00032 #include "Conventions/Munits.h"
00033 #include "Conventions/PlaneView.h"
00034 #include "DataUtil/GetMomFromRange.h"
00035 #include "MessageService/MsgService.h"
00036 #include "MinosObjectMap/MomNavigator.h"
00037 #include "Navigation/NavKey.h"
00038 #include "Navigation/NavSet.h"
00039 #include "Navigation/XxxItr.h"
00040 #include "Plex/PlexStripEndId.h"
00041 #include "Plex/PlexPlaneId.h"
00042 #include "RecoBase/ArrivalTime.h"
00043 #include "RecoBase/CandSliceHandle.h"
00044 #include "RecoBase/CandStripHandle.h"
00045 #include "RecoBase/CandStripListHandle.h"
00046 #include "RecoBase/LinearFit.h"
00047 #include "RecoBase/PropagationVelocity.h"
00048 #include "UgliGeometry/UgliGeomHandle.h"
00049 #include "UgliGeometry/UgliScintPlnHandle.h"
00050 #include "Validity/VldContext.h"
00051 #include "TMath.h"
00052 
00053 NavKey KeyOnStripTPos( const CandStripHandle *cth )
00054 {
00055   return cth->GetTPos();
00056 }
00057 
00058 ClassImp(AlgTrackSR)
00059 
00060 //______________________________________________________________________
00061   CVSID("$Id: AlgTrackSR.cxx,v 1.71 2007/02/04 06:03:24 rhatcher Exp $");
00062 
00063 //______________________________________________________________________
00064 AlgTrackSR::AlgTrackSR()
00065 {
00066 }
00067 
00068 //______________________________________________________________________
00069 AlgTrackSR::~AlgTrackSR()
00070 {
00071 }
00072 
00073 //______________________________________________________________________
00074 void AlgTrackSR::RunAlg(AlgConfig &ac, CandHandle &ch, CandContext &cx)
00075 {
00076 
00077 
00078   MSG("Alg", Msg::kDebug) << "Starting AlgTrackSR::RunAlg()" << endl;
00079 
00080   // get parameters
00081   fParmNExtraStrip = ac.GetInt("NExtraStrip");
00082   fParmMinPlanePE = ac.GetDouble("MinPlanePE");
00083   fParmIsCosmic = ac.GetInt("IsCosmic");
00084   fParmMinStripPE = ac.GetDouble("MinStripPE");
00085   fParmMaxClusterSizeTimeFit = ac.GetDouble("MaxClusterSizeTimeFit");
00086   fParmMaxTimeFitRes = ac.GetDouble("MaxTimeFitRes");
00087   fParmMinTimeFitSize = ac.GetInt("MinTimeFitSize");
00088   fParmMaxTimeFitOutlier = ac.GetDouble("MaxTimeFitOutlier");
00089 
00090   assert(ch.InheritsFrom("CandTrackSRHandle"));
00091   CandTrackSRHandle &cth = dynamic_cast<CandTrackSRHandle &>(ch);
00092 
00093   assert(cx.GetDataIn());
00094 
00095   assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00096 
00097  
00098   const TObjArray *tary =
00099     dynamic_cast<const TObjArray*>(cx.GetDataIn());
00100 
00101   Track2DSR *utrack=0;
00102   Track2DSR *vtrack=0;
00103   CandStripHandle *strip = 0;
00104 
00105   Int_t ntrackdigit=0;
00106 
00107   const CandSliceHandle *slice = 0;
00108 
00109   for (Int_t i=0; i<=tary->GetLast(); i++) {
00110     TObject *tobj = tary->At(i);
00111     if (tobj->InheritsFrom("CandSliceHandle")) {
00112       slice = dynamic_cast<CandSliceHandle*>(tobj);
00113       continue;
00114     }
00115     Track2DSR *track = dynamic_cast<Track2DSR*>(tobj);
00116     assert(track);
00117     if(!utrack && track->GetPlaneView()==PlaneView::kU) utrack = track;
00118     if(!vtrack && track->GetPlaneView()==PlaneView::kV) vtrack = track;
00119 
00120   }//end loop to find 2D tracks
00121 
00122   FindStripsInTrack(utrack, vtrack, cth);
00123 
00124   CandStripHandleItr trackStripItr(cth.GetDaughterIterator());
00125   while( (strip = trackStripItr()) ){
00126     MSG("AlgTrackSR", Msg::kDebug) << strip->GetPlane() << " " << strip->GetTPos() << endl;
00127   }
00128   MSG("AlgTrackSR", Msg::kDebug) << endl;
00129 
00130   assert(slice);
00131   CandStripHandleItr sliceItr(slice->GetDaughterIterator());
00132   Int_t maxplane=0;
00133   while( (strip = sliceItr()) ){
00134     if(strip->GetPlane()>maxplane) maxplane = strip->GetPlane();
00135   }
00136   assert(maxplane);
00137 
00138   //there can be at most 485 planes in either detector - round up
00139   //to 500
00140   Float_t planepe[500];
00141   Float_t zpe[500];
00142   for (int i=0; i<=maxplane; i++) {
00143     planepe[i] = 0.;
00144     zpe[i] = 0.;
00145   }
00146 
00147   sliceItr.Reset();
00148 
00149   while( (strip = sliceItr()) ){
00150     if(strip->GetPlane()>=0){
00151       planepe[strip->GetPlane()] += strip->GetCharge(CalDigitType::kPE);
00152       zpe[strip->GetPlane()] = strip->GetZPos();
00153     }
00154   }
00155 
00156   CandStripHandleItr stripItr(cth.GetDaughterIterator());
00157   CandStripHandle *firststrip = stripItr.Ptr();
00158 
00159   assert(firststrip);
00160 
00161   const VldContext *vldcontext = firststrip->GetVldContext();
00162   assert(vldcontext);
00163 
00164   //get the detector type from the VldContext
00165   fDetector = vldcontext->GetDetector();
00166 
00167   UgliGeomHandle ugh(*vldcontext);
00168 
00169   // calculate transverse positions
00170 
00171   MSG("AlgTrackSR", Msg::kDebug) << "in detector type " 
00172                                  << Detector::AsString(fDetector) 
00173                                  << endl
00174                                  << " u track direction = " 
00175                                  << utrack->GetDirection()
00176                                  << " v track direction = " 
00177                                  << vtrack->GetDirection()
00178                                  << endl;
00179 
00180   if(utrack->GetDirection()>0.) cth.SetTimeSlope(1.);
00181   else cth.SetTimeSlope(-1.);
00182 
00183   CandRecoHandle *candreco = dynamic_cast<CandRecoHandle *>(&cth);
00184   assert(candreco);
00185   
00186   Int_t begplane = candreco->GetBegPlane();
00187   Int_t endplane = candreco->GetEndPlane();
00188   
00189   cth.SetVtxPlane(begplane);
00190   cth.SetVtxZ(cth.GetZ(begplane));
00191   cth.SetVtxU(cth.GetU(begplane));
00192   cth.SetVtxV(cth.GetV(begplane));
00193   cth.SetDirCosU(cth.GetDirCosU(begplane));
00194   cth.SetDirCosV(cth.GetDirCosV(begplane));
00195   cth.SetDirCosZ(cth.GetDirCosZ(begplane));
00196   
00197   cth.SetEndPlane(endplane);
00198   cth.SetEndZ(cth.GetZ(endplane));
00199   cth.SetEndU(cth.GetU(endplane));
00200   cth.SetEndV(cth.GetV(endplane));
00201   cth.SetEndDirCosU(cth.GetDirCosU(endplane));
00202   cth.SetEndDirCosV(cth.GetDirCosV(endplane));
00203   cth.SetEndDirCosZ(cth.GetDirCosZ(endplane));
00204 
00205   SetUVZT(&cth);
00206 
00207   Int_t nshower = 0;
00208   Int_t ntrackstrip = 0;
00209   Double_t aveTime=0;
00210   while( (strip = stripItr()) ){
00211     nshower = cth.IsInShower(strip);
00212     aveTime += strip->GetTime(); 
00213     if(nshower<=1){
00214       ++ntrackstrip;
00215       ntrackdigit += strip->GetNDaughters();
00216     }
00217   }//end loop to see how many track strips there are
00218 
00219   if(cth.GetNDaughters()>0) aveTime /=cth.GetNDaughters();
00220   cth.SetVtxT(aveTime);
00221   cth.SetEndT(aveTime);
00222 
00223   stripItr.Reset();
00224   
00225   cth.SetNTrackStrip(ntrackstrip);
00226   
00227   //timewalk correction constants
00228   //  Double_t c1 = ac.GetDouble("TimeWalk1"); //-14.45;
00229   //  Double_t c2 = ac.GetDouble("TimeWalk2"); //7.498
00230   //  Double_t c3 = ac.GetDouble("TimeWalk3"); //-1.566;
00231   Double_t maxPathLength = 0.;
00232   Double_t timefitchi2 = 0.;
00233   Double_t fitparm[] = {0.,0.};
00234  
00235   
00236   Int_t ntimedigit = FindTimingDirection(cth, stripItr, utrack, vtrack, 
00237                                          fitparm, maxPathLength, timefitchi2);
00238  
00239   if(ntimedigit<=fParmMinTimeFitSize){
00240     fitparm[0]=aveTime;
00241     fitparm[1]=0;
00242     timefitchi2=0;
00243   }
00244 
00245   stripItr.Reset();
00246 
00247   if(!fParmIsCosmic 
00248      || fitparm[1]*utrack->GetDirection()>0.)
00249     cth.SetTimeSlope(TMath::Abs(fitparm[1]));
00250   else
00251     cth.SetTimeSlope(-1.*TMath::Abs(fitparm[1]));
00252   
00253   if(fitparm[1]<0. || !fParmIsCosmic){ 
00254     // in beam mode, reconstruct tracks backward, but now want track going forward
00255     cth.SetTimeOffset(fitparm[0]+maxPathLength*fitparm[1]);
00256     
00257     // direction changed from initial guess
00258     MSG("AlgTrackSR",Msg::kDebug) << "DIRECTION CHANGED from " 
00259                                   << utrack->GetDirection() << endl;
00260     Int_t vtxplane = cth.GetVtxPlane();
00261     Int_t endplane = cth.GetEndPlane();
00262     cth.SetVtxPlane(endplane);
00263     cth.SetEndPlane(vtxplane);
00264     
00265     MSG("AlgTrackSR",Msg::kDebug) << "new vertex and end planes = " 
00266                                   << cth.GetVtxPlane() << " " 
00267                                   << cth.GetEndPlane() << endl;
00268     SetUVZT(&cth);
00269   } 
00270   else
00271     cth.SetTimeOffset(fitparm[0]);
00272   
00273   cth.SetNTrackDigit(ntrackdigit);
00274   cth.SetNTimeFitDigit(ntimedigit);
00275   cth.SetTimeFitChi2(timefitchi2);
00276   
00277   MSG("AlgTrackSR",Msg::kVerbose) << "timeslope = " << cth.GetTimeSlope() << endl;
00278   MSG("AlgTrackSR",Msg::kVerbose) << "ntrackdigit = " << cth.GetNTrackDigit() 
00279                                   << endl;
00280   MSG("AlgTrackSR",Msg::kVerbose) << "ntimefitdigit = " << cth.GetNTimeFitDigit() 
00281                                   << endl;
00282   MSG("TrackSR",Msg::kVerbose) << "timefitchi2 = " << cth.GetTimeFitChi2() << endl;
00283 
00284   begplane = cth.GetBegPlane();
00285   endplane = cth.GetEndPlane();
00286 
00287   Track2DSR *utrack0 = utrack->Dup();
00288   Track2DSR *vtrack0 = vtrack->Dup();
00289 
00290   cth.SetUTrack(utrack0);
00291   cth.SetVTrack(vtrack0);
00292 
00293   cth.SetVtxPlane(begplane);
00294   cth.SetVtxZ(cth.GetZ(begplane));
00295   cth.SetVtxU(cth.GetU(begplane));
00296   cth.SetVtxV(cth.GetV(begplane));
00297   cth.SetVtxT(cth.GetTimeOffset());
00298   cth.SetDirCosU(cth.GetDirCosU(begplane));
00299   cth.SetDirCosV(cth.GetDirCosV(begplane));
00300   cth.SetDirCosZ(cth.GetDirCosZ(begplane));
00301 
00302   cth.SetEndPlane(endplane);
00303   cth.SetEndZ(cth.GetZ(endplane));
00304   cth.SetEndU(cth.GetU(endplane));
00305   cth.SetEndV(cth.GetV(endplane));
00306   cth.SetEndT(cth.GetTimeOffset()+cth.GetTimeSlope()*cth.GetdS());
00307   cth.SetEndDirCosU(cth.GetDirCosU(endplane));
00308   cth.SetEndDirCosV(cth.GetDirCosV(endplane));
00309   cth.SetEndDirCosZ(cth.GetDirCosZ(endplane));
00310 
00311   Int_t plane0 = -999;
00312   Int_t plane1 = -999;
00313   Int_t idir = (cth.GetDirCosZ()>0. ? 1 : -1);
00314   Double_t minTPos = 0.;
00315   Double_t maxTPos = 0.;
00316   Double_t minZ    = 0.;
00317   Double_t maxZ    = 0.;
00318   ugh.GetTransverseExtent(PlaneView::kU, minTPos, maxTPos);
00319   ugh.GetZExtent(minZ, maxZ);
00320 
00321   for (int i=0; i<maxplane; i++) {
00322     if (planepe[i]>=fParmMinPlanePE && planepe[i+1]>=fParmMinPlanePE) {
00323       if (idir>0) {
00324         if (plane0<0) {
00325           plane0 = i;
00326         }
00327         if (plane1<0 || i+1>plane1) {
00328           plane1 = i+1;
00329         }
00330       } else {
00331         if (plane1<0) {
00332           plane1 = i;
00333         }
00334         if (plane0<0 || i+1>plane0) {
00335           plane0 = i+1;
00336         }
00337       }
00338     }
00339   }
00340   MSG("TrackSR",Msg::kDebug) << "Vertex and End Planes >= " << fParmMinPlanePE 
00341                              << " pe: " << plane0 << " " << plane1 << endl;
00342 
00343   if(fParmIsCosmic){
00344     Int_t thisplane = cth.GetVtxPlane()-idir;
00345     Double_t thisu = cth.GetVtxU();
00346     Double_t thisv = cth.GetVtxV();
00347     Double_t thisz = cth.GetVtxZ();
00348     Double_t thisx = 0.707*(thisu-thisv);
00349     Double_t thisy = 0.707*(thisu+thisv);
00350     Double_t thisdu = cth.GetDirCosU()/cth.GetDirCosZ();
00351     Double_t thisdv = cth.GetDirCosV()/cth.GetDirCosZ();
00352     Double_t thisdx = 0.707*(thisdu-thisdv);
00353     Double_t thisdy = 0.707*(thisdu+thisdv);
00354     Int_t nplaneskip = 0;
00355     while(thisplane>=0 && thisplane<=500 && thisplane*idir>=plane0*idir 
00356           && thisx<=maxTPos && thisx>=minTPos 
00357           && thisy<=maxTPos && thisy>=minTPos
00358           && thisu<=maxTPos && thisu>=minTPos
00359           && thisv<=maxTPos && thisv>=minTPos 
00360           && nplaneskip<=3
00361           ){
00362 
00363       PlexPlaneId plexplane(fDetector,thisplane);
00364       
00365       if( plexplane.IsValid() ){
00366         UgliScintPlnHandle planeid = ugh.GetScintPlnHandle(plexplane);
00367         if( planeid.IsValid() ){
00368           Double_t dz = (Double_t)planeid.GetZ0()-thisz;
00369           if (thisx<=maxTPos && thisx>=minTPos 
00370               && thisy<=maxTPos && thisy>=minTPos
00371               && thisu<=maxTPos && thisu>=minTPos
00372               && thisv<=maxTPos && thisv>=minTPos){ 
00373             thisz = (Double_t)planeid.GetZ0();
00374             thisx += dz*thisdx;
00375             thisy += dz*thisdy;
00376             thisu += dz*thisdu;
00377             thisv += dz*thisdv;
00378             cth.SetVtxU(thisu);
00379             cth.SetVtxV(thisv);
00380             cth.SetVtxZ(thisz);
00381             cth.SetVtxPlane(thisplane);
00382             cth.SetU(thisplane,thisu);
00383             cth.SetV(thisplane,thisv);
00384             MSG("TrackSR",Msg::kDebug) << "SETTING VERTEX, (u,v,z) = (" << thisu 
00385                                        << "," << thisv << "," << thisz << ")" << endl;
00386           }
00387           if (planepe[thisplane]>0.) {
00388             nplaneskip = 0;
00389           } else {
00390             nplaneskip++;
00391           }
00392         }
00393       }
00394       thisplane -= idir;
00395     }
00396     thisplane = cth.GetEndPlane()+idir;
00397     thisu = cth.GetEndU();
00398     thisv = cth.GetEndV();
00399     thisz = cth.GetEndZ();
00400     thisx = 0.707*(thisu-thisv);
00401     thisy = 0.707*(thisu+thisv);
00402     thisdu = cth.GetEndDirCosU()/cth.GetEndDirCosZ();
00403     thisdv = cth.GetEndDirCosV()/cth.GetEndDirCosZ();
00404     thisdx = 0.707*(thisdu-thisdv);
00405     thisdy = 0.707*(thisdu+thisdv);
00406     nplaneskip = 0;
00407     while(thisplane>=0 && thisplane<=500 && thisplane*idir<=plane1*idir 
00408           && thisx<=maxTPos && thisx>=minTPos 
00409           && thisy<=maxTPos && thisy>=minTPos
00410           && thisu<=maxTPos && thisu>=minTPos
00411           && thisv<=maxTPos && thisv>=minTPos
00412           && nplaneskip<=3
00413           ){
00414 
00415       PlexPlaneId plexplane(fDetector,thisplane);
00416       if( plexplane.IsValid() ){
00417         UgliScintPlnHandle planeid = ugh.GetScintPlnHandle(plexplane);
00418         if (planeid.IsValid()) {
00419           Double_t dz = (Double_t)planeid.GetZ0()-thisz;
00420           if (thisx<=maxTPos && thisx>=minTPos 
00421               && thisy<=maxTPos && thisy>=minTPos
00422               && thisu<=maxTPos && thisu>=minTPos
00423               && thisv<=maxTPos && thisv>=minTPos){
00424             thisz = (Double_t)planeid.GetZ0();
00425             thisx += dz*thisdx;
00426             thisy += dz*thisdy;
00427             thisu += dz*thisdu;
00428             thisv += dz*thisdv;
00429             cth.SetEndU(thisu);
00430             cth.SetEndV(thisv);
00431             cth.SetEndZ(thisz);
00432             cth.SetEndPlane(thisplane);
00433             cth.SetU(thisplane,thisu);
00434             cth.SetV(thisplane,thisv);
00435             MSG("TrackSR",Msg::kDebug) << "SETTING END, (u,v,z) = (" << thisu 
00436                                        << "," << thisv << "," << thisz << ")" << endl;
00437           }
00438           if (planepe[thisplane]>0.) {
00439             nplaneskip = 0;
00440           } else {
00441             nplaneskip++;
00442           }
00443         }
00444       }
00445       thisplane += idir;
00446     }
00447   }
00448 
00449 
00450   // traces calculated assuming far detector
00451 
00452   Double_t thistrace[2][4] = {{0.,0.,0.,0.},{0.,0.,0.,0.}}; // u,v,x,y
00453   Double_t thistracez[2][4] = {{0.,0.,0.,0.},{0.,0.,0.,0.}};
00454 
00455   Double_t thisx[2][4] = {{cth.GetVtxU(),cth.GetVtxV(),0.,0.},
00456                            {cth.GetEndU(),cth.GetEndV(),0.,0.}};
00457   Double_t thisz[2] = {cth.GetVtxZ(),cth.GetEndZ()};
00458 
00459   Double_t thisdcos[2][4] = {{-cth.GetDirCosU(),-cth.GetDirCosV(),
00460                               0.,0.},
00461                              {cth.GetEndDirCosU(),cth.GetEndDirCosV(),
00462                               0.,0.}};
00463 
00464   Double_t thisdcosz[2] = {-cth.GetDirCosZ(),cth.GetEndDirCosZ()};
00465 
00466   Double_t mintrace[2] = {999.,999.};
00467   Double_t mintracez[2] = {999.,999.};
00468 
00469   for (Int_t i=0; i<2; i++) {
00470     for (Int_t j=0; j<4; j++) {
00471       if (j==2) {
00472         thisx[i][j] = .70710678*(thisx[i][0]-thisx[i][1]);
00473         thisdcos[i][j] = .70710678*(thisdcos[i][0]-thisdcos[i][1]);
00474       }     
00475       if (j==3) {
00476         thisx[i][j] = .70710678*(thisx[i][0]+thisx[i][1]);
00477         thisdcos[i][j] = .70710678*(thisdcos[i][0]+thisdcos[i][1]);
00478       }
00479       if (thisx[i][j]<=maxTPos && thisx[i][j]>=minTPos
00480           && thisz[i]>=minZ && thisz[i]<=maxZ) {
00481         if (thisdcos[i][j]>0.) {
00482           thistrace[i][j] = (maxTPos-thisx[i][j])/TMath::Abs(thisdcos[i][j]);
00483         } 
00484         else if( thisdcos[i][j] < 0.) {
00485           thistrace[i][j] = (maxTPos+thisx[i][j])/TMath::Abs(thisdcos[i][j]);
00486         }
00487         thistracez[i][j] = thistrace[i][j]*TMath::Abs(thisdcosz[i]);
00488 // check in z direction
00489         Double_t delz = 0.;
00490         if (thisdcosz[i]>0) {
00491           delz = maxZ-thisz[i];
00492         } else {
00493           delz = thisz[i]-minZ;
00494         }
00495         if (delz<thistracez[i][j] && thisdcosz[i] != 0.) {
00496           thistrace[i][j] = delz/TMath::Abs(thisdcosz[i]);
00497           thistracez[i][j] = delz;
00498         }
00499       }
00500       if (thistrace[i][j]<mintrace[i]) {
00501         mintrace[i] = thistrace[i][j];
00502         mintracez[i] = thistracez[i][j];
00503       }
00504     }
00505   }
00506 
00507   cth.SetVtxTrace(mintrace[0]);
00508   cth.SetVtxTraceZ(mintracez[0]);
00509   cth.SetEndTrace(mintrace[1]);
00510   cth.SetEndTraceZ(mintracez[1]);
00511 
00512   SetdS(&cth);
00513   Double_t range = cth.GetRange();
00514   if (range>0.) {
00515     //    cth.SetMomentum(0.105658*exp(1.0363*log(range)-4.383));
00516   // taken from PDG R/M vs p/M plot for Fe
00517   // until June 2006 used to be  cth.SetMomentum(0.048 + 1.660e-3*range + 3.057e-8*range*range);
00518      cth.SetMomentum(GetMomFromRange(range));   
00519   } else {
00520     cth.SetMomentum(0.); // range not valid, set momentum to 0 GeV/c
00521   }
00522   Calibrate(&cth);
00523 
00524 }
00525 
00526 //______________________________________________________________________
00527 void AlgTrackSR::Trace(const char * /* c */) const
00528 {
00529 }
00530 
00531 //----------------------------------------------------------------------
00532 //figure out which strips belong in the track
00533 
00534 void AlgTrackSR::FindStripsInTrack(Track2DSR *uTrack, Track2DSR *vTrack,
00535                                    CandTrackSRHandle &cth)
00536 {
00537   MSG("AlgTrackSR", Msg::kDebug) << "in find strips in track" << endl;
00538   
00539   Track2DSR *track = 0;
00540   Int_t trkCtr = 0;
00541 
00542   TrackClusterSR *tc = 0;
00543   CandStripHandle *strip = 0;
00544   Int_t nstripexp = 0;
00545   PlaneView::PlaneView_t stripView = PlaneView::kUnknown;
00546   
00547   Int_t itcbef=0,itcaft=0;
00548   TrackClusterSR *tcbef = 0;
00549   TrackClusterSR *tcaft = 0;
00550   Double_t localslope = 0.;
00551   Double_t localtpos = 0.;
00552 
00553   Int_t numStrips = 0;
00554   
00555   Float_t stripTPos[192];
00556   Float_t stripCharge[192];
00557   Int_t sortedStrips[192];
00558   
00559   Int_t indxbest = 0;
00560   Double_t distbest = 0;
00561   Double_t dist = 0.;
00562   Double_t distsum = 0.;
00563  
00564   Int_t nshower=0;
00565   
00566   while( trkCtr<2 ){
00567     
00568     track = uTrack;
00569     if(trkCtr==1) track = vTrack;
00570 
00571     ++trkCtr;
00572 
00573     //loop over the clusters in the track
00574     for(Int_t itc=0; itc<track->GetLast()+1; ++itc){
00575       tc = track->GetCluster(itc);
00576       numStrips = tc->GetStripList()->GetLast()+1;
00577       cth.SetTrackPointError(tc->GetPlane(),tc->GetTPosError());
00578       // assume 4:1 aspect ratio for scintillator strips
00579       nstripexp = TMath::Max(1,(Int_t)(TMath::Abs(track->GetSlope(itc))/4.)+1);
00580       nstripexp += fParmNExtraStrip;
00581       
00582       //if this cluster has less than the maximum number of strips allowed
00583       //take all the strips
00584       if(numStrips<=nstripexp){
00585         
00586         MSG("AlgTrackSR", Msg::kDebug) << "cluster on plane " << tc->GetPlane()
00587                                        << " has " 
00588                                        << numStrips
00589                                        << "/" << nstripexp << " strips" << endl;
00590 
00591         //loop over the strips in the cluster
00592         for(Int_t istrip=0; istrip<numStrips; ++istrip){
00593           
00594           strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(istrip));
00595           stripView = strip->GetPlaneView();
00596           
00597           if(stripView == PlaneView::kU || stripView == PlaneView::kV){
00598             cth.AddDaughterLink(*strip);
00599             cth.SetInShower(strip,0);
00600           }
00601         }
00602       }
00603       else{ 
00604         
00605         //more than the expected number of strips,
00606         //only add strips which match fit
00607         
00608         MSG("AlgTrackSR", Msg::kDebug) << "cluster on plane " << tc->GetPlane()
00609                                          << " has " 
00610                                          << tc->GetStripList()->GetLast()+1
00611                                          << "/" << nstripexp << " strips" << endl;
00612 
00613         itcbef=0;
00614         itcaft=0;
00615         
00616         //get the bounding clusters.  if on an endpoint cluster get the two
00617         //closest to it
00618         if(itc == 0){
00619           itcbef = FindNeighborIndex(track, itc, 1);
00620           itcaft = FindNeighborIndex(track, itcbef, 1);
00621         }
00622         else if(itc == track->GetLast() ){
00623           itcbef = FindNeighborIndex(track, itc, -1);
00624           itcaft = FindNeighborIndex(track, itcbef, -1);
00625         }
00626         else{
00627           itcbef = FindNeighborIndex(track, itc,-1);
00628           itcaft = FindNeighborIndex(track, itc, 1);
00629         }
00630         
00631         tcbef = track->GetCluster(itcbef);
00632         tcaft = track->GetCluster(itcaft);
00633         
00634         MSG("AlgTrackSR", Msg::kVerbose) << itc
00635                                          << " looking at cluster on plane "
00636                                          << tc->GetPlane() 
00637                                          << " prev cluster on plane "
00638                                          << tcbef->GetPlane()
00639                                          << " netx cluster on plane "
00640                                          << tcaft->GetPlane() << endl;
00641 
00642         //get the slope from the clusters bounding the current one
00643         localslope = tcaft->GetTPos()-tcbef->GetTPos();
00644         if(tcaft->GetZPos()!=tcbef->GetZPos())
00645           localslope /= tcaft->GetZPos()-tcbef->GetZPos();
00646 //      else MSG("AlgTrackSR", Msg::kWarning) << "cannot find local slope "
00647 //                                            << "clusters have same zpos  "
00648 //                                            << localslope << endl;
00649         
00650         //get the expected tpos at the current cluster
00651         localtpos = tcbef->GetTPos();
00652         localtpos += localslope*(tc->GetZPos()-tcbef->GetZPos());
00653 
00654         MSG("AlgTrackSR", Msg::kVerbose) << "local slope = " << localslope
00655                                          << " local tpos = " << localtpos
00656                                          << endl;
00657         
00658         
00659         //find the set of strips in the cluster of size nstripexp which is closest
00660         //to the expect tpos
00661         
00662         //first order the strips by tpos
00663         for(Int_t i = 0; i<numStrips; ++i){
00664           strip = dynamic_cast<CandStripHandle *>(tc->GetStripList()->At(i));
00665           stripTPos[i] = strip->GetTPos();
00666           stripCharge[i] = strip->GetCharge();
00667           MSG("AlgTrackSR", Msg::kVerbose) << stripTPos[i] << endl;
00668         }
00669 
00670         //sort the array of TPos values in ascending order
00671         TMath::Sort(numStrips,stripTPos,sortedStrips,false);
00672 
00673         indxbest = 0;
00674         distbest = 0;
00675         dist = 0.;
00676         distsum = 0.;
00677         for(Int_t indx=0; indx<=numStrips-nstripexp; ++indx){
00678 
00679           distsum = 0.;
00680 
00681           for(Int_t istrip=indx; istrip<indx+nstripexp; ++istrip){
00682 
00683             MSG("AlgTrackSR", Msg::kVerbose) << indx << " " << istrip 
00684                                              << " " << sortedStrips[istrip]
00685                                              << " " << stripTPos[sortedStrips[istrip]]
00686                                              << endl;
00687 
00688             dist = TMath::Abs(stripTPos[sortedStrips[istrip]]-localtpos);
00689             if(fParmIsCosmic==1)
00690               dist *= (0.3+exp(-0.2*stripCharge[sortedStrips[istrip]]));
00691             
00692             distsum += dist;
00693           }
00694           if (!indx || distsum<distbest) {
00695             distbest = distsum;
00696             indxbest = indx;
00697           }
00698         }
00699         
00700         // count how many hit strips above 1 pe
00701         nshower=0;
00702         for(int i= 0; i<=numStrips-nstripexp; ++i){
00703           strip = dynamic_cast<CandStripHandle *>(tc->GetStripList()->At(i));
00704           if(strip->GetCharge()>fParmMinStripPE) ++nshower;
00705         }
00706 
00707         for(Int_t istrip=indxbest; istrip<indxbest+nstripexp; istrip++) {
00708           strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(sortedStrips[istrip]));
00709           MSG("AlgTrackSR", Msg::kVerbose) << strip->GetPlane() << " " 
00710                                            << strip->GetTPos() << endl;
00711           cth.AddDaughterLink(*strip);
00712           cth.SetInShower(strip,nshower-nstripexp);
00713         }
00714       }//end if more than expected number of strips in cluster
00715     }//end loop over track clusters
00716   }//end loop over tracks
00717 
00718 
00719   return;
00720 }
00721 
00722 //----------------------------------------------------------------------
00723 Int_t AlgTrackSR::FindNeighborIndex(Track2DSR *track, 
00724                                     Int_t currentIndex,
00725                                     Int_t direction)
00726 {
00727   MSG("AlgTrackSR", Msg::kVerbose) << "in find neighbor index" << endl;
00728   
00729   Bool_t found = false;
00730   Int_t neighborIndex = currentIndex+direction;
00731   Int_t nstrip = 0;
00732   Int_t dstrip = 0;
00733   TrackClusterSR *tc = 0;
00734 
00735   //loop over clusters in the track starting from the one
00736   //right next to the current cluster - take the next cluster in the list
00737   //which has the expected number of strips in it based on the slope at 
00738   //that cluster
00739   while(!found 
00740         && neighborIndex>=0 && neighborIndex<=track->GetLast()
00741         ){
00742 
00743     MSG("AlgTrackSR", Msg::kVerbose) << "current index = " << currentIndex
00744                                      <<" neighbor index = " << neighborIndex
00745                                      << endl;
00746     
00747     tc = track->GetCluster(neighborIndex);
00748     dstrip = tc->GetStripList()->GetLast()+1;
00749     dstrip -= TMath::Max(1,(Int_t)(TMath::Abs(track->GetSlope(neighborIndex))/4.));
00750     if(dstrip<=nstrip) found = true;
00751           
00752     ++nstrip;
00753     neighborIndex += direction;
00754   }
00755 
00756   neighborIndex -= direction;
00757   assert(neighborIndex>=0 && neighborIndex<= track->GetLast());
00758   
00759   return neighborIndex;
00760 }
00761 
00762 //---------------------------------------------------------------------
00763 //figure out if the track is going forward or backwards with respect to
00764 //z.  return the number of points used to determine the direction
00765 Int_t AlgTrackSR::FindTimingDirection(CandTrackSRHandle &cth,
00766                                       CandStripHandleItr &stripItr, 
00767                                       Track2DSR *utrack,
00768                                       Track2DSR *vtrack,
00769                                       Double_t *fitparm,
00770                                       Double_t &maxPathLength,
00771                                       Double_t &timefitchi2)
00772 {
00773 
00774   MSG("AlgTrackSR", Msg::kVerbose) << "in find timing direction" << endl;
00775 
00776   fitparm[0]=0;
00777   fitparm[1]=0;
00778   stripItr.Reset();
00779   Int_t plane = 0;
00780   maxPathLength = 0.;
00781   //  Double_t digitpe = 0.;
00782   // Double_t logadc = 0.;
00783   const Double_t min_wgt = 0.4;
00784 
00785   CandStripHandle *strip = 0;
00786   CandDigitHandle *digit = 0;
00787 
00788   StripEnd::StripEnd_t stripEnd = StripEnd::kUnknown;
00789 
00790   //get VldContext and UgliGeomHandle
00791   CandStripHandle *firststrip = stripItr.Ptr();
00792 
00793   assert(firststrip);
00794 
00795   const VldContext *vldcontext = firststrip->GetVldContext();
00796   assert(vldcontext);
00797 
00798   UgliGeomHandle ugh(*vldcontext);
00799 
00800   //arrays for determining the timing direction - limit to
00801   //1000 entries - really shouldnt need more than that to 
00802   //be able to figure out if the event is going north or 
00803   //south
00804   Double_t time[1000];
00805   Double_t pathLength[1000];
00806   Double_t weight[1000];
00807   Int_t digitCtr = 0;
00808 
00809   for(Int_t i = 0; i<1000; ++i){
00810     time[i] = 0.;
00811     pathLength[i] = 0.;
00812     weight[i] = 0.;
00813   }
00814   
00815   //need to put in a check to make sure you dont have more planes in one
00816   //view than the other
00817 
00818   //get the endpoints for each track
00819   Int_t uBeg = utrack->GetBegPlane();
00820   Int_t uEnd = utrack->GetEndPlane();
00821   Int_t vBeg = vtrack->GetBegPlane();
00822   Int_t vEnd = vtrack->GetEndPlane();
00823 
00824   Int_t trackBeg = TMath::Min(uBeg,vBeg);
00825   Int_t trackEnd = TMath::Max(uEnd,vEnd);
00826   Int_t direction = 1;
00827 
00828   //if the track was loaded up going north to south
00829   if(uBeg>uEnd){
00830     direction = -1;
00831     trackBeg = TMath::Max(uBeg,vBeg);
00832     trackEnd = TMath::Min(uEnd,vEnd);
00833   }
00834   if(TMath::Abs(uBeg-vBeg)>3.){
00835     if(direction>0) trackBeg = TMath::Max(uBeg,vBeg) - 1;
00836     else if(direction<0) trackBeg = TMath::Min(uBeg,vBeg) + 1;
00837   } 
00838   if(TMath::Abs(uEnd-vEnd)>3.){
00839     if(direction>0) trackEnd = TMath::Min(uEnd,vEnd) + 1;
00840     if(direction<0) trackEnd = TMath::Max(uEnd,vEnd) - 1;
00841   }
00842   MSG("AlgTrackSR", Msg::kDebug) << "u end points = " << uBeg
00843                                  << " " << uEnd << " v end points = "
00844                                  << vBeg << " " << vEnd << endl
00845                                  << " track end points = " << trackBeg
00846                                  << " " << trackEnd << endl;
00847 
00848   Double_t halfLength = 0.;
00849   Double_t apos = 0.;
00850   Double_t distFromCenter = 0.;
00851   Double_t fiberDist = 0.;
00852   
00853   //loop over all strips in the track
00854   while( (strip = stripItr()) && digitCtr<1000){
00855     
00856     plane = strip->GetPlane();
00857 
00858     MSG("AlgTrackSR", Msg::kDebug) << "strip from plane " << plane 
00859                                    << " " << Detector::AsString(fDetector) << endl;
00860 
00861     //only look at double ended strips if in the far detector
00862     //and strips from planes that have orthogonal measurements 
00863     if( (strip->GetNDaughters() == 2 || fDetector == Detector::kNear)  
00864        && plane*direction>=trackBeg*direction 
00865        && plane*direction<=trackEnd*direction){
00866       
00867       MSG("AlgTrackSR", Msg::kVerbose) << "strip good for timing" << endl;
00868 
00869       //get the iterator over the digits in the strip
00870       CandDigitHandleItr digitItr(strip->GetDaughterIterator());
00871 
00872       while( (digit = digitItr()) && digitCtr<1000){
00873         
00874         stripEnd = digit->GetPlexSEIdAltL().GetEnd();
00875               
00876         pathLength[digitCtr] = cth.GetdS()-cth.GetdS(plane);
00877         maxPathLength = TMath::Max(maxPathLength, 
00878                                    (double)(cth.GetdS()-cth.GetdS(plane)));
00879         
00880         
00881         //time is recorded in seconds, but i prefer working in ns.
00882         //get the time from the track rather than the strips
00883         time[digitCtr] = cth.GetT(plane,stripEnd)*1.e9;
00884 
00885         //cant use the GetT method for Near detector right now because the near 
00886         //detector has multiple digits for each strip end and GetT only returns one 
00887         //time - gotta fix that.
00888         if(fDetector == Detector::kNear){
00889           UgliStripHandle stripHandle = ugh.GetStripHandle(strip->GetStripEndId());
00890           halfLength = stripHandle.GetHalfLength();
00891           apos = 0.;
00892           fiberDist = 0.;
00893           
00894           if(strip->GetPlaneView() == PlaneView::kV) apos = cth.GetU(plane);
00895           else if(strip->GetPlaneView() == PlaneView::kU) apos = cth.GetV(plane);
00896           distFromCenter = 0.;
00897           if(strip->GetPlaneView() == PlaneView::kV) 
00898             distFromCenter = (stripEnd == StripEnd::kNegative) ? apos : -apos;
00899           else if(strip->GetPlaneView() == PlaneView::kU) 
00900             distFromCenter = (stripEnd == StripEnd::kNegative) ? -apos : apos;
00901           
00902           fiberDist = (halfLength + distFromCenter + stripHandle.ClearFiber(stripEnd) 
00903                        + stripHandle.WlsPigtail(stripEnd));
00904           
00905           time[digitCtr] = strip->GetTime(digit->GetPlexSEIdAltL().GetEnd()) - fiberDist/PropagationVelocity::Velocity();
00906           time[digitCtr] *= 1.e9;
00907           
00908           //digitpe = digit->GetCharge(CalDigitType::kPE);
00909           //    if (digitpe>0.0) logadc =  TMath::Log(digitpe)/2.3;
00910           //  else {
00911           //  logadc = 0.0;
00912           // MSG("AlgTrackSR", Msg::kWarning)
00913           //    << "FindTimingDirection() digitpe was "
00914           //   << digitpe << ", avoid TMath::Log()." << endl;
00915           //  }
00916           //      time[digitCtr] -= (c1*logadc
00917           //         +c2*logadc*logadc
00918           //         +c3*logadc*logadc*logadc);
00919         }//end if near detector
00920 
00921         //dont use hits in showers for calculating timing
00922         if(cth.IsInShower(strip)<=fParmMaxClusterSizeTimeFit)
00923           weight[digitCtr] = TMath::Min(cth.GetTimeWeight(digit),min_wgt);
00924         
00925 
00926         MSG("AlgTrackSR", Msg::kDebug) << digitCtr 
00927                                        << " " << pathLength[digitCtr]
00928                                        << " " << time[digitCtr]
00929                                        << " " << strip->GetTime(stripEnd)*1.e9
00930                                        << " " << weight[digitCtr]
00931                                        << endl;
00932         
00933         ++digitCtr;
00934         
00935       }//end loop over digits in strip
00936     }//end if double ended strip
00937   }//end loop over strips in track
00938 
00939   Double_t efitparm[2];
00940   Double_t maxresid = fParmMaxTimeFitRes*1.e9+1.;
00941   Double_t resid = 0.;
00942   Int_t ctr = digitCtr;
00943   Int_t imaxresid = 0;
00944   Int_t nremove = 0;
00945   Bool_t dofit=false;
00946   while(maxresid>fParmMaxTimeFitRes
00947         && digitCtr-nremove>fParmMinTimeFitSize 
00948         && (1.*nremove)<fParmMaxTimeFitOutlier*digitCtr
00949         ){
00950     dofit=true;
00951     LinearFit::Weighted(ctr,pathLength,time,weight,fitparm,efitparm);
00952     maxresid = 0.;
00953     imaxresid = 0;
00954     for(int i=0; i<ctr; ++i){
00955       resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
00956       if(weight[i]>0. && TMath::Abs(resid)>maxresid){
00957         maxresid = TMath::Abs(resid);
00958         imaxresid = i;
00959       }
00960     }
00961     
00962     MSG("AlgTrackSR",Msg::kDebug) << "constrained fit, dt/ds slope, offset = " 
00963                                   << fitparm[1] << " " << fitparm[0] << endl;
00964     MSG("AlgTrackSR",Msg::kDebug) << "maximum residual " << maxresid << " " 
00965                                   << pathLength[imaxresid] << " " 
00966                                   << time[imaxresid] << " " 
00967                                   << weight[imaxresid] << endl;
00968     
00969     if(maxresid>fParmMaxTimeFitRes){
00970       MSG("AlgTrackSR",Msg::kVerbose) << "  removing" << endl;
00971       weight[imaxresid] = 0.;
00972       nremove++;
00973     }
00974     
00975     timefitchi2 = 0.;
00976     for(int i=0; i<ctr; ++i){
00977       resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
00978       timefitchi2 += weight[i]*resid*resid;
00979     }
00980     MSG("AlgTrackSR",Msg::kDebug) << "chi2 = " << timefitchi2 
00981                                   << "   n = " << digitCtr-nremove << endl;
00982   }//end loop to find timing fit and remove outliers
00983   
00984   timefitchi2 = 0.;
00985   if(dofit){
00986     for(int i=0; i<ctr; ++i){
00987       if(weight[i]>0.){
00988         MSG("AlgTrackSR",Msg::kDebug) << "TIMEFIT " << i << " " 
00989                                       <<pathLength[i] << " " 
00990                                       << time[i] << " " << weight[i] << endl;
00991         resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
00992       timefitchi2 += weight[i]*resid*resid;
00993       }
00994     }
00995   }
00996     
00997   MSG("AlgTrackSR",Msg::kDebug) << " TIMEFIT " << fitparm[1] 
00998                                 << " chi^2/ndf = " 
00999                                 << timefitchi2/(1.*(digitCtr-nremove)) 
01000                                 << " " << TMath::Abs(uBeg-uEnd)
01001                                 << " " << TMath::Abs(vBeg-vEnd) << endl;
01002   
01003   //check the chi^2 for the fit - if it is really bad dont change the
01004   //initial guess at the direction for the event, ie just make sure
01005   //that the slope in time vs pathlength is positive;
01006   if(fitparm[1]<0. && timefitchi2>=10.*(digitCtr-nremove))
01007     fitparm[1] *= -1.;
01008 
01009   //return to units of seconds
01010   fitparm[0] *= 1.e-9;
01011   fitparm[1] *= 1.e-9;
01012 
01013   return digitCtr-nremove;
01014 }