AlgTrackSR Class Reference

#include <AlgTrackSR.h>

Inheritance diagram for AlgTrackSR:
AlgBase AlgReco AlgTrack

List of all members.

Public Member Functions

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

Private Member Functions

Int_t FindNeighborIndex (Track2DSR *track, Int_t currentIndex, Int_t direction)
Int_t FindTimingDirection (CandTrackSRHandle &cth, CandStripHandleItr &stripItr, Track2DSR *utrack, Track2DSR *vtrack, Double_t *fitparm, Double_t &maxPathLength, Double_t &timefitchi2)
void FindStripsInTrack (Track2DSR *uTrack, Track2DSR *vTrack, CandTrackSRHandle &cth)

Private Attributes

Int_t fParmNExtraStrip
Double_t fParmMinPlanePE
Int_t fParmIsCosmic
Detector::Detector_t fDetector
Double_t fParmMinStripPE
Double_t fParmMaxClusterSizeTimeFit
Double_t fParmMaxTimeFitRes
Int_t fParmMinTimeFitSize
Double_t fParmMaxTimeFitOutlier

Detailed Description

Definition at line 26 of file AlgTrackSR.h.


Constructor & Destructor Documentation

AlgTrackSR::AlgTrackSR (  ) 

Definition at line 64 of file AlgTrackSR.cxx.

00065 {
00066 }

AlgTrackSR::~AlgTrackSR (  )  [virtual]

Definition at line 69 of file AlgTrackSR.cxx.

00070 {
00071 }


Member Function Documentation

Int_t AlgTrackSR::FindNeighborIndex ( Track2DSR track,
Int_t  currentIndex,
Int_t  direction 
) [private]

Definition at line 723 of file AlgTrackSR.cxx.

References Track2DSR::GetCluster(), Track2DSR::GetSlope(), TrackClusterSR::GetStripList(), Msg::kVerbose, MSG, and tc.

Referenced by FindStripsInTrack().

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 }

void AlgTrackSR::FindStripsInTrack ( Track2DSR uTrack,
Track2DSR vTrack,
CandTrackSRHandle cth 
) [private]

Definition at line 534 of file AlgTrackSR.cxx.

References CandHandle::AddDaughterLink(), FindNeighborIndex(), fParmIsCosmic, fParmMinStripPE, fParmNExtraStrip, CandStripHandle::GetCharge(), Track2DSR::GetCluster(), Track2DSR::GetLast(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), Track2DSR::GetSlope(), TrackClusterSR::GetStripList(), CandStripHandle::GetTPos(), TrackClusterSR::GetTPos(), TrackClusterSR::GetTPosError(), TrackClusterSR::GetZPos(), Msg::kDebug, PlaneView::kU, PlaneView::kUnknown, PlaneView::kV, Msg::kVerbose, MSG, CandTrackHandle::SetInShower(), CandTrackHandle::SetTrackPointError(), TMath::Sort(), and tc.

Referenced by RunAlg().

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 }

Int_t AlgTrackSR::FindTimingDirection ( CandTrackSRHandle cth,
CandStripHandleItr &  stripItr,
Track2DSR utrack,
Track2DSR vtrack,
Double_t *  fitparm,
Double_t &  maxPathLength,
Double_t &  timefitchi2 
) [private]

Definition at line 765 of file AlgTrackSR.cxx.

References bfld::AsString(), UgliStripHandle::ClearFiber(), digit(), fDetector, fParmMaxClusterSizeTimeFit, fParmMaxTimeFitOutlier, fParmMaxTimeFitRes, fParmMinTimeFitSize, Track2DSR::GetBegPlane(), CandHandle::GetDaughterIterator(), CandTrackHandle::GetdS(), PlexSEIdAltL::GetEnd(), Track2DSR::GetEndPlane(), UgliStripHandle::GetHalfLength(), CandHandle::GetNDaughters(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandTrackHandle::GetT(), CandStripHandle::GetTime(), CandTrackSRHandle::GetTimeWeight(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), CandHandle::GetVldContext(), CandTrackHandle::IsInShower(), Msg::kDebug, Detector::kNear, StripEnd::kNegative, PlaneView::kU, StripEnd::kUnknown, PlaneView::kV, Msg::kVerbose, MSG, PropagationVelocity::Velocity(), LinearFit::Weighted(), and UgliStripHandle::WlsPigtail().

Referenced by RunAlg().

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 }

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

Implements AlgBase.

Definition at line 74 of file AlgTrackSR.cxx.

References bfld::AsString(), AlgReco::Calibrate(), Track2DSR::Dup(), fDetector, FindStripsInTrack(), FindTimingDirection(), fParmIsCosmic, fParmMaxClusterSizeTimeFit, fParmMaxTimeFitOutlier, fParmMaxTimeFitRes, fParmMinPlanePE, fParmMinStripPE, fParmMinTimeFitSize, fParmNExtraStrip, CandRecoHandle::GetBegPlane(), CandStripHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Track2DSR::GetDirection(), Registry::GetDouble(), CandRecoHandle::GetEndPlane(), Registry::GetInt(), GetMomFromRange(), CandHandle::GetNDaughters(), CandStripHandle::GetPlane(), Track2DSR::GetPlaneView(), UgliGeomHandle::GetScintPlnHandle(), CandStripHandle::GetTime(), CandStripHandle::GetTPos(), UgliGeomHandle::GetTransverseExtent(), UgliPlnHandle::GetZ0(), UgliGeomHandle::GetZExtent(), CandStripHandle::GetZPos(), PlexPlaneId::IsValid(), UgliScintPlnHandle::IsValid(), Msg::kDebug, CalDigitType::kPE, PlaneView::kU, PlaneView::kV, Msg::kVerbose, MSG, AlgTrack::SetdS(), and AlgTrack::SetUVZT().

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 }

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

Reimplemented from AlgBase.

Definition at line 527 of file AlgTrackSR.cxx.

00528 {
00529 }


Member Data Documentation

Definition at line 40 of file AlgTrackSR.h.

Referenced by FindTimingDirection(), and RunAlg().

Int_t AlgTrackSR::fParmIsCosmic [private]

Definition at line 39 of file AlgTrackSR.h.

Referenced by FindStripsInTrack(), and RunAlg().

Definition at line 42 of file AlgTrackSR.h.

Referenced by FindTimingDirection(), and RunAlg().

Definition at line 45 of file AlgTrackSR.h.

Referenced by FindTimingDirection(), and RunAlg().

Double_t AlgTrackSR::fParmMaxTimeFitRes [private]

Definition at line 43 of file AlgTrackSR.h.

Referenced by FindTimingDirection(), and RunAlg().

Double_t AlgTrackSR::fParmMinPlanePE [private]

Definition at line 38 of file AlgTrackSR.h.

Referenced by RunAlg().

Double_t AlgTrackSR::fParmMinStripPE [private]

Definition at line 41 of file AlgTrackSR.h.

Referenced by FindStripsInTrack(), and RunAlg().

Definition at line 44 of file AlgTrackSR.h.

Referenced by FindTimingDirection(), and RunAlg().

Definition at line 37 of file AlgTrackSR.h.

Referenced by FindStripsInTrack(), and RunAlg().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1