AlgTrackCamList Class Reference

#include <AlgTrackCamList.h>

Inheritance diagram for AlgTrackCamList:
AlgBase

List of all members.

Public Member Functions

 AlgTrackCamList ()
virtual ~AlgTrackCamList ()
virtual void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void Trace (const char *c) const
void RunTheFinder (CandSliceHandle *slice)
void FormTheHits (CandSliceHandle *slice)
void FormTheClusters ()
void IDTrkAndShwClusters ()
void FormTriplets ()
void FindAllAssociations ()
void FindPreferredJoins ()
void FindMatchedJoins ()
void FirstUVComparison ()
void Form2DTracks ()
void Join2DTracks ()
void SecondUVComparison ()
void Form3DTracks (CandSliceHandle *slice)
void LookForHitsAcrossGap (TrackCam *Trk)
void ExtendTrack (TrackCam *Trk)
void FillGapsInTrack (TrackCam *Trk)
void NearDetectorTriplets ()
void ClearUp ()
void MatchUV (TrackCam *trku, TrackCam *trkv)

Private Attributes

vector< HitCam * > AllHitBank [500]
vector< HitCam * > HitBank [500]
vector< HitCam * > HitsInTracks [2]
vector< ClusterCam * > ClusterBank [500]
vector< ClusterCam * > ClusterList [500]
vector< TrackSegmentCam * > SegmentBank [500]
vector< TrackSegmentCam * > NDSegmentBank [500]
vector< TrackSegmentCam * > ViewSegBank [2]
vector< TrackSegmentCam * > TempTrack [2]
vector< TrackSegmentCam * > PossibleJoins [2]
vector< TrackCam * > FinalTrackBank [2]
VldContextvldc
int NumModules
int PlanesInModule
int ModuleType
const double StripWidth
double PECut
double PECut2
bool CambridgeAnalysis

Detailed Description

Definition at line 23 of file AlgTrackCamList.h.


Constructor & Destructor Documentation

AlgTrackCamList::AlgTrackCamList (  ) 

Definition at line 52 of file AlgTrackCamList.cxx.

00052                                  :
00053   StripWidth(4.108e-2)
00054 {
00055 }

AlgTrackCamList::~AlgTrackCamList (  )  [virtual]

Definition at line 60 of file AlgTrackCamList.cxx.

00061 {
00062 }


Member Function Documentation

void AlgTrackCamList::ClearUp (  ) 

Definition at line 5353 of file AlgTrackCamList.cxx.

References AllHitBank, ClusterBank, ClusterList, FinalTrackBank, HitBank, HitsInTracks, NDSegmentBank, SegmentBank, TempTrack, and ViewSegBank.

Referenced by RunAlg().

05354 {
05355   // Reset containers and tidy memory usage
05356   unsigned int i, j;
05357 
05358   TempTrack[0].clear(); TempTrack[1].clear();
05359   ViewSegBank[0].clear(); ViewSegBank[1].clear();
05360   HitsInTracks[0].clear(); HitsInTracks[1].clear(); 
05361 
05362   for (i=0; i<2; ++i) {
05363     for(j=0; j<FinalTrackBank[i].size(); ++j) {delete FinalTrackBank[i][j];}
05364     FinalTrackBank[i].clear();
05365   }
05366 
05367 
05368   for (i=0; i<500; ++i) {
05369     for(j=0; j<AllHitBank[i].size(); ++j) {if(AllHitBank[i][j]) delete AllHitBank[i][j];}
05370     AllHitBank[i].clear();
05371     HitBank[i].clear();
05372 
05373     for(j=0; j<ClusterBank[i].size(); ++j) {if(ClusterBank[i][j]) delete ClusterBank[i][j];}
05374     ClusterBank[i].clear();
05375 
05376     for(j=0; j<ClusterList[i].size(); ++j) {if(ClusterList[i][j]) delete ClusterList[i][j];}
05377     ClusterList[i].clear();
05378 
05379     for(j=0; j<SegmentBank[i].size(); ++j) {if(SegmentBank[i][j]) delete SegmentBank[i][j];}
05380     SegmentBank[i].clear();
05381 
05382     for(j=0; j<NDSegmentBank[i].size(); ++j) {if(NDSegmentBank[i][j]) delete NDSegmentBank[i][j];}
05383     NDSegmentBank[i].clear();
05384   }
05385 
05386 }

void AlgTrackCamList::ExtendTrack ( TrackCam Trk  ) 

Definition at line 4913 of file AlgTrackCamList.cxx.

References TrackCam::AddHit(), AllHitBank, TrackCam::GetBegDir(), TrackCam::GetBegPlane(), TrackCam::GetBegTPos(), TrackCam::GetBegZPos(), TrackCam::GetEndDir(), TrackCam::GetEndPlane(), TrackCam::GetEndTPos(), TrackCam::GetEndZPos(), TrackCam::GetPlaneView(), HitCam::GetTPos(), HitCam::GetZPos(), HitsInTracks, Msg::kVerbose, ModuleType, MSG, PlanesInModule, and StripWidth.

Referenced by Form3DTracks().

04914 {
04915   MSG("AlgTrackCamList", Msg::kVerbose) << " Attempting to extend track at beginning and end " << endl;
04916   
04917   int PlanesSinceAdded, CurrentPlane(0), Increment(0);
04918   unsigned int nhits;
04919   double TrkDir(0.), TrkTPos(0.), TrkZPos(0.), PredictedTPos, extrem1, extrem2;
04920   double StripCentre, CurrentZPos, MinDistanceToHit, StripCharge;
04921 
04922   double Tolerance; int MaxPlanes;
04923   int PlaneView=Trk->GetPlaneView();
04924 
04925   bool ConsiderHit;
04926 
04927   if(ModuleType==1 || ModuleType==3) {MaxPlanes=8; Tolerance=2.*StripWidth;}
04928   else {MaxPlanes=10; Tolerance=4.*StripWidth;} // Leigh: Changed from 40 
04929 
04930   // Extend beginning first, then the end.
04932   for(int direction=0; direction<2; ++direction) {
04933   
04934     PlanesSinceAdded=1; 
04935     if(direction==0) {CurrentPlane=Trk->GetBegPlane()-1; Increment=-1;}
04936     else if(direction==1) {CurrentPlane=Trk->GetEndPlane()+1; Increment=1;}
04937 
04938     
04939     while(PlanesSinceAdded<=MaxPlanes && CurrentPlane>=0 && CurrentPlane<=PlanesInModule) { 
04940       // Don't do anything if there are no hits, or we aren't in the same view.
04941       nhits=AllHitBank[CurrentPlane].size();
04942       if(nhits>0 && AllHitBank[CurrentPlane][0]->GetPlaneView()==PlaneView) {
04943 
04944         if(direction==0) {
04945           TrkDir=Trk->GetBegDir(); 
04946           TrkTPos=Trk->GetBegTPos(); 
04947           TrkZPos=Trk->GetBegZPos(); 
04948         }
04949         else if(direction==1) {
04950           TrkDir=Trk->GetEndDir(); 
04951           TrkTPos=Trk->GetEndTPos(); 
04952           TrkZPos=Trk->GetEndZPos(); 
04953         }
04954 
04955         HitCam* CurrentHit = 0;
04956         CurrentZPos = AllHitBank[CurrentPlane][0]->GetZPos();
04957 
04958         MinDistanceToHit = Tolerance;
04959         
04960         // Loop over the hits on the plane
04961         for(unsigned int i=0; i<nhits; ++i) {
04962 
04963           // Don't add a hit that is already in another track in the slice
04964           ConsiderHit=true;
04965           for(unsigned int p=0; p<HitsInTracks[PlaneView].size(); ++p) {
04966             if(AllHitBank[CurrentPlane][i]==HitsInTracks[PlaneView][p]) {ConsiderHit=false; break;}
04967           }
04968           if(ConsiderHit==false) {continue;}
04969 
04970 
04971           // Find the most likely hit to add
04972           PredictedTPos = TrkTPos + TrkDir*(CurrentZPos-TrkZPos);
04973           extrem1 = PredictedTPos + (0.055*TrkDir);
04974           extrem2 = PredictedTPos - (0.055*TrkDir);
04975           
04976           StripCentre=AllHitBank[CurrentPlane][i]->GetTPos();
04977           StripCharge=AllHitBank[CurrentPlane][i]->GetCharge();
04978           
04979           MSG("AlgTrackCamList", Msg::kVerbose) << "Prediction " << direction << " " << PredictedTPos 
04980                                                 << " Actual " << StripCentre   << " CurrentZPos " << CurrentZPos << " trkdir " << TrkDir << endl;
04981           
04982           if(fabs(StripCentre-PredictedTPos)<MinDistanceToHit 
04983              && (ModuleType==1 || ModuleType==3 || StripCharge>2 || (fabs(StripCentre-PredictedTPos)<StripWidth && PlanesSinceAdded<20)) ) {
04984             MinDistanceToHit = fabs(StripCentre-PredictedTPos);
04985             CurrentHit = AllHitBank[CurrentPlane][i];
04986           }
04987           else if(fabs(StripCentre-extrem1)<MinDistanceToHit 
04988                   && (ModuleType==1 || ModuleType==3 || StripCharge>2 || (fabs(StripCentre-extrem1)<StripWidth && PlanesSinceAdded<20)) ) {
04989             MinDistanceToHit = fabs(StripCentre-extrem1);
04990             CurrentHit = AllHitBank[CurrentPlane][i];
04991           }
04992           else if(fabs(StripCentre-extrem2)<MinDistanceToHit 
04993                   && (ModuleType==1 || ModuleType==3 || StripCharge>2 || (fabs(StripCentre-extrem2)<StripWidth && PlanesSinceAdded<20)) ) {
04994             MinDistanceToHit = fabs(StripCentre-extrem2);
04995             CurrentHit = AllHitBank[CurrentPlane][i];
04996           }
04997         }
04998         
04999         // Add best hit we have found
05000         if(CurrentHit) {
05001           Trk->AddHit(CurrentHit); 
05002           HitsInTracks[PlaneView].push_back(CurrentHit);
05003 
05004           PlanesSinceAdded=0;
05005 
05006           MSG("AlgTrackCamList", Msg::kVerbose) << "Track extended " << direction << " TPos " << CurrentHit->GetTPos() 
05007                                                 << " ZPos " << CurrentHit->GetZPos() << endl;
05008         }
05009       }
05010       else {PlanesSinceAdded++;}
05011       
05012       CurrentPlane+=Increment;
05013     }
05014   }
05016 
05017 
05018   return;
05019 }

void AlgTrackCamList::FillGapsInTrack ( TrackCam Trk  ) 

Definition at line 5025 of file AlgTrackCamList.cxx.

References TrackCam::AddHit(), AllHitBank, TrackCam::GetBegPlane(), HitCam::GetCharge(), TrackCam::GetEndPlane(), TrackCam::GetEntries(), TrackCam::GetHit(), HitCam::GetPlane(), TrackCam::GetPlaneView(), HitCam::GetTPos(), HitCam::GetZPos(), HitsInTracks, Msg::kVerbose, ModuleType, MSG, size, and StripWidth.

Referenced by Form3DTracks().

05026 {
05027   MSG("AlgTrackCamList", Msg::kVerbose) << " Attempting to fill gaps in track " << endl;
05028 
05029 
05030   // First store track hits in plane order
05031   vector<HitCam*> TrackHits[490];
05032   int PlaneView=Trk->GetPlaneView();
05033 
05034   for(unsigned int i=0; i<Trk->GetEntries(); ++i) {
05035     HitCam* Hit = Trk->GetHit(i);
05036     TrackHits[Hit->GetPlane()].push_back(Hit);
05037   }
05038 
05039   
05040   // Find gaps
05041   int PrevPlaneWithHit=Trk->GetBegPlane();
05042   int NextPlaneWithHit=Trk->GetEndPlane();
05043 
05044   double PrevTPos, NextTPos, PrevZPos, NextZPos, CurrentZPos;
05045   double PredictedGradient, PredictedTPos, MinDistanceToHit, StripCharge;
05046   double Tolerance; bool ConsiderHit;
05047 
05048   if(ModuleType==1 || ModuleType==3) {Tolerance=1.5*StripWidth;}
05049   else {Tolerance=3.*StripWidth;}
05050 
05051   for(int i=Trk->GetBegPlane()+1; i<Trk->GetEndPlane(); i++) {
05052     // Don't do anything if there are no hits, or we aren't in the same view.
05053     if(AllHitBank[i].size()>0 && AllHitBank[i][0]->GetPlaneView()==PlaneView) {
05054       
05055       // Find the hit planes on either side of a gap.
05056       if(TrackHits[i].size()>0) {PrevPlaneWithHit=i; continue;}
05057 
05058       for(int j=i+2; j<=Trk->GetEndPlane(); j++) {
05059         if(TrackHits[j].size()>0) {NextPlaneWithHit=j; break;}
05060       }
05061  
05062       // Now try and fill the gap. Limit to maximum gap size?
05063       if(PrevPlaneWithHit!=NextPlaneWithHit) {
05064         CurrentZPos=AllHitBank[i][0]->GetZPos();
05065 
05066         PrevZPos=TrackHits[PrevPlaneWithHit][0]->GetZPos();
05067         NextZPos=TrackHits[NextPlaneWithHit][0]->GetZPos();
05068 
05069         PrevTPos=TrackHits[PrevPlaneWithHit][0]->GetTPos();
05070         NextTPos=TrackHits[NextPlaneWithHit][0]->GetTPos();
05071 
05072         PredictedGradient=(NextTPos-PrevTPos)/(NextZPos-PrevZPos);
05073         PredictedTPos=PrevTPos+((CurrentZPos-PrevZPos)*PredictedGradient);
05074 
05075         HitCam* CurrentHit = 0;
05076 
05077         MinDistanceToHit=Tolerance;
05078 
05079         for(unsigned int k=0; k<AllHitBank[i].size(); ++k) {
05080           HitCam* hit = AllHitBank[i][k];
05081 
05082 
05083           // Don't add a hit that is already in another track in the slice
05084           ConsiderHit=true;
05085           for(unsigned int p=0; p<HitsInTracks[PlaneView].size(); ++p) {
05086             if(hit==HitsInTracks[PlaneView][p]) {ConsiderHit=false; break;}
05087           }
05088           if(ConsiderHit==false) {continue;}
05089           
05090 
05091           // Find the most likely hit to add
05092           StripCharge=hit->GetCharge();
05093 
05094           if(fabs(PredictedTPos-hit->GetTPos())<MinDistanceToHit 
05095              && (ModuleType==1 || ModuleType==3 || StripCharge>2 || fabs(PredictedTPos-hit->GetTPos())<StripWidth)) {
05096 
05097             MinDistanceToHit=fabs(PredictedTPos-hit->GetTPos());
05098             CurrentHit=hit;
05099           }
05100         }
05101       
05102         if(CurrentHit) {
05103           Trk->AddHit(CurrentHit); 
05104           TrackHits[i].push_back(CurrentHit);
05105           HitsInTracks[PlaneView].push_back(CurrentHit);
05106 
05107           MSG("AlgTrackCamList", Msg::kVerbose) << "Gap filled at TPos " << CurrentHit->GetTPos()
05108                                                 << " ZPos " << CurrentHit->GetZPos() << endl;
05109         }
05110       }
05111     }
05112   }
05113   
05114   // Clean up
05115   for(int i=0; i<490; ++i) {TrackHits[i].clear();}
05116   
05117   return;
05118 }

void AlgTrackCamList::FindAllAssociations (  ) 

Definition at line 984 of file AlgTrackCamList.cxx.

References TrackSegmentCam::GetAssocSegBeg(), TrackSegmentCam::GetAssocSegEnd(), TrackSegmentCam::GetBegPlane(), TrackSegmentCam::GetBegTPos(), TrackSegmentCam::GetEndPlane(), TrackSegmentCam::GetEndTPos(), TrackSegmentCam::GetNAssocSegBeg(), TrackSegmentCam::GetNAssocSegEnd(), Msg::kVerbose, MSG, NumModules, PlanesInModule, and SegmentBank.

Referenced by RunTheFinder().

00985 {
00986   // For each triplet formed, we consider nearby triplets and see whether
00987   // the triplets are associated. 
00988 
00989   // For each triplet, we simply find the nearby other triplets that have 
00990   // a compatible beginning/end position/direction.
00991 
00992   // Later we look more that one triplet away, to find strings of 
00993   // "preferred associations" that might indicate a track.
00994 
00995 
00996   int i=0;
00997 
00998 
00999   // Loop over the planes
01001   for(int Module=0; Module<NumModules; ++Module) {
01002 
01003     for(int Plane=3; Plane<PlanesInModule-1; ++Plane) {
01004       i=Plane + Module*(PlanesInModule+1);
01005         
01006       if((i+2)<500 && SegmentBank[i].size()>0 && SegmentBank[i+2].size()>0) {
01007         
01008         // Look at the segments on plane i
01009         for(unsigned int k0=0; k0<SegmentBank[i].size(); ++k0) {
01010           
01011           // For each of these, look at the segments on plane i+2
01012           for(unsigned int kp=0; kp<SegmentBank[i+2].size(); ++kp) {
01013             
01014             // Check associations. 
01015             if(SegmentBank[i][k0]->IsAssoc(SegmentBank[i+2][kp])) {
01016               
01017               // If segments are associated, we add to list of ALL possible joins
01018               SegmentBank[i][k0]->AddAssocSegToEnd(SegmentBank[i+2][kp]);
01019               SegmentBank[i+2][kp]->AddAssocSegToBeg(SegmentBank[i][k0]);
01020             }
01021           }
01022         }
01023       }
01024     } // End loop over planes in module
01025   } // End loop over modules
01027 
01028 
01029 
01030 
01031   // Print out list of all associations
01033   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF ALL SEG ASSOCIATIONS *** " << endl;
01034 
01035   // Print out list of triplets
01036   for(int i=0; i<500; ++i) {
01037     for(unsigned int j=0;j<SegmentBank[i].size(); ++j) {
01038       TrackSegmentCam* segment = SegmentBank[i][j];
01039 
01040       MSG("AlgTrackCamList", Msg::kVerbose) << "--- Plane " << i << ", Seg Number " << j 
01041                                             << ", BegPlane " << segment->GetBegPlane() 
01042                                             << ", EndPlane " << segment->GetEndPlane() 
01043                                             << ", BegTPos " << segment->GetBegTPos() 
01044                                             << ", EndTPos " << segment->GetEndTPos() 
01045                                             <<endl;
01046  
01047       MSG("AlgTrackCamList", Msg::kVerbose) << " Beg: " <<endl;
01048       for(unsigned int k=0; k<segment->GetNAssocSegBeg(); ++k) {
01049         TrackSegmentCam* segtmp = segment->GetAssocSegBeg(k);
01050         MSG("AlgTrackCamList", Msg::kVerbose) << "  Assoc number=" << k
01051                                               << "  begpln: " << segtmp->GetBegPlane()
01052                                               << "  begtpos: " << segtmp->GetBegTPos()
01053                                               << "  endpln: " << segtmp->GetEndPlane()
01054                                               << "  endtpos: " << segtmp->GetEndTPos() << endl;
01055       }
01056       
01057       MSG("AlgTrackCamList", Msg::kVerbose) << " End: " <<endl;
01058       for(unsigned int k=0; k<segment->GetNAssocSegEnd(); ++k) {
01059         TrackSegmentCam* segtmp = segment->GetAssocSegEnd(k);
01060         MSG("AlgTrackCamList", Msg::kVerbose) << "  Assoc number=" << k
01061                                               << "  begpln: " << segtmp->GetBegPlane()
01062                                               << "  begtpos: " << segtmp->GetBegTPos()
01063                                               << "  endpln: " << segtmp->GetEndPlane()
01064                                               << "  endtpos: " << segtmp->GetEndTPos() << endl;           
01065       }
01066     }
01067   }
01069 
01070   return;
01071 }

void AlgTrackCamList::FindMatchedJoins (  ) 

Definition at line 1762 of file AlgTrackCamList.cxx.

References TrackSegmentCam::AddMatchSegToBeg(), TrackSegmentCam::AddMatchSegToEnd(), TrackSegmentCam::AddSegment(), TrackSegmentCam::GetBegDir(), TrackSegmentCam::GetBegPlane(), TrackSegmentCam::GetBegTPos(), TrackSegmentCam::GetBegZPos(), TrackSegmentCam::GetEndDir(), TrackSegmentCam::GetEndPlane(), TrackSegmentCam::GetEndTPos(), TrackSegmentCam::GetEndZPos(), TrackSegmentCam::GetEntries(), TrackSegmentCam::GetMatchSegBeg(), TrackSegmentCam::GetMatchSegEnd(), TrackSegmentCam::GetNMatchSegBeg(), TrackSegmentCam::GetNMatchSegEnd(), TrackSegmentCam::GetNPrefSegEnd(), TrackSegmentCam::GetSeedSegment(), TrackSegmentCam::IsAssoc(), Msg::kVerbose, ModuleType, MSG, NumModules, PlanesInModule, SegmentBank, StripWidth, and ViewSegBank.

Referenced by RunTheFinder().

01763 {
01764   // Having made all the preferred assocations, we actually match triplets 
01765   // together to form longer segments of the track.
01766 
01767   // We use the idea of seed segments to identify the last segment in a chain
01768   // of segments.
01769 
01770 
01771 
01772   int i;
01773   bool JoinFlag;
01774 
01775   // Loop over the planes
01776   for(int Module=0; Module<NumModules; ++Module) {
01777     for(int Plane=3; Plane<PlanesInModule-1; ++Plane) {
01778   
01779       i=Plane + Module*(PlanesInModule+1);
01780       
01781       // Loop over the triplets centered on plane i
01782       for(unsigned int j=0; j<SegmentBank[i].size(); ++j) {
01783         
01784         // Firstly, consider isolated segments. These must be their own seed segments.
01785         if(SegmentBank[i][j]->GetNAssocSegBeg()==0 && SegmentBank[i][j]->GetNAssocSegEnd()==0) {
01786           ViewSegBank[SegmentBank[i][j]->GetPlaneView()].push_back(SegmentBank[i][j]);
01787           SegmentBank[i][j]->SetSeedSegment(SegmentBank[i][j]);
01788         }
01789 
01790 
01791         // Then, consider those segments which have some preferred associations
01792         if(SegmentBank[i][j]->GetNPrefSegBeg()>0 || SegmentBank[i][j]->GetNPrefSegEnd()>0) {
01793           JoinFlag=false;
01794 
01795           // If there is one pref assoc at beginning (always false for first segment in module)
01796           if(SegmentBank[i][j]->GetNPrefSegBeg()==1) {
01797             TrackSegmentCam* Segm = SegmentBank[i][j]->GetPrefSegBeg(0);
01798 
01799             // Check that we are considering a simple chain of segments
01800             if(Segm->GetNPrefSegEnd()==1) {
01801               // Get the last segment in the chain of segments
01802               TrackSegmentCam* SeedSeg = Segm->GetSeedSegment();
01803 
01804               SeedSeg->AddSegment(SegmentBank[i][j]); 
01805               SegmentBank[i][j]->SetSeedSegment(SeedSeg);
01806               SegmentBank[i][j]->SetUID(-1);
01807               JoinFlag=true;
01808             }
01809           }
01810 
01811           // Now consider those triplets which aren't just one pref assoc at beg and one at end
01812           if(JoinFlag==false) {
01813             ViewSegBank[SegmentBank[i][j]->GetPlaneView()].push_back(SegmentBank[i][j]);
01814             SegmentBank[i][j]->SetSeedSegment(SegmentBank[i][j]);
01815             
01816             // Loop over the preferred associations at beginning
01817             for(unsigned int k=0; k<SegmentBank[i][j]->GetNPrefSegBeg(); ++k) {
01818               TrackSegmentCam* SeedSeg = SegmentBank[i][j]->GetPrefSegBeg(k)->GetSeedSegment();
01819               TrackSegmentCam* MatchSeg = SegmentBank[i][j];
01820 
01821               if(SeedSeg->GetBegPlane()<=MatchSeg->GetBegPlane() && SeedSeg->GetEndPlane()<=MatchSeg->GetEndPlane()) {
01822 
01823                 if( (SeedSeg->GetEntries()>3 && MatchSeg->GetEntries()>3)
01824                     || 4.*fabs( (SeedSeg->GetBegTPos()-SeedSeg->GetEndTPos())/(SeedSeg->GetEndPlane()-SeedSeg->GetBegPlane()) -
01825                                 (MatchSeg->GetBegTPos()-MatchSeg->GetEndTPos())/(MatchSeg->GetEndPlane()-MatchSeg->GetBegPlane()) ) < 10*StripWidth ) {
01826                   MatchSeg->AddMatchSegToBeg(SeedSeg);
01827                   SeedSeg->AddMatchSegToEnd(MatchSeg);
01828                 }
01829               }
01830 
01831             }
01832           }
01833         }
01834       } // End loop over triplets centered on plane i
01835       
01836     } // End loop over planes in module
01837 
01838   } // End loop over modules
01840 
01841 
01842 
01843   // Form matched associations between non-adjacent segments
01845   vector<TrackSegmentCam*> TempSeg1;
01846   vector<TrackSegmentCam*> TempSeg2;
01847 
01848   const int npts0=4; int npts; double inv_npts; int Module; int Plane;
01849   bool CheckCoilHole;
01850   
01851   
01852   // Loop over the planeviews U and V
01853   for(int view=0; view<2; ++view) {
01854     TempSeg1.clear(); TempSeg2.clear();
01855 
01856  
01857     // Loop over entries stored in ViewSegBank, above
01858     for(unsigned int i=0; i<ViewSegBank[view].size(); ++i) {
01859       TrackSegmentCam* Seg1 = ViewSegBank[view][i];
01860 
01861       // If no matched entries at end of segment
01862       if(Seg1->GetNMatchSegEnd()==0) {
01863         JoinFlag=false;
01864 
01865         
01866         // For FD, attempt to match segments across the coil hole
01867         // For ND, need to be careful making associations, so we don't join two separate tracks together
01869 
01870         // Check which module we are in
01871         Module=int(Seg1->GetEndPlane()/(PlanesInModule+1));
01872         npts=npts0;
01873 
01874         // These are the tpos values required to locate the FD coil hole
01875         if( (ModuleType==1 && Seg1->GetEndTPos()>-0.62 && Seg1->GetEndTPos()<0.62) 
01876             || ( 1+(Seg1->GetEndPlane()-Seg1->GetBegPlane())/2>5) ) {
01877           
01878           // Calculate range of possible association and configure
01879           // - worst case scenario - track crosses entire width of coil hole, tcoil = ~0.6m
01880           // - it will take a distance in z, zcoil = tcoil * dz/dt to do this
01881           // - zcoil is equivalent to nplanes = zcoil/0.0594 = tcoil / ( 0.0594 * dt/dz )
01882           // - npts = nplanes/2.0 = 1.0 / ( 0.198 * dt/dz)
01883           // - where dt/dz = Seg1->GetEndDir();
01884           inv_npts=0.198*Seg1->GetEndDir(); //was inv_npts=0.5*Seg1->GetEndDir();
01885           
01886           if(inv_npts<0) {inv_npts=-inv_npts;}
01887           if(inv_npts<0.01) {inv_npts=0.01;}
01888           npts=int(1./inv_npts);
01889 
01890           if(npts<13) {npts=13;}    //replaced 8 with 13
01891           if(npts>100) {npts=100;}  //replaced 18 with 30 - possibly this should be a function of the gradient of the track?
01892         }
01893 
01894         
01895         // Loop over this range of possible association
01896         for(int j=1; j<npts; ++j) {
01897           Plane=Seg1->GetEndPlane()+(2*j);
01898           
01899           if( JoinFlag==false && Plane<((PlanesInModule+1)*(Module+1)) ) {
01900             
01901             // Loop over segments on each possible plane
01902             for(unsigned int k=0; k<SegmentBank[Plane].size(); ++k) {
01903               TrackSegmentCam* Seg2 = SegmentBank[Plane][k];
01904 
01905               // Check association, first across coil-hole 
01906               if(ModuleType==1 && (( Seg2->GetBegTPos()>-0.62 && Seg2->GetBegTPos()<0.21
01907                                      && Seg1->GetEndTPos()>-0.21 && Seg1->GetEndTPos()<0.62)
01908                                    || (Seg2->GetBegTPos()>-0.21 && Seg2->GetBegTPos()<0.62
01909                                        && Seg1->GetEndTPos()>-0.62 && Seg1->GetEndTPos()<0.21)) ) {CheckCoilHole=true;}
01910               else {CheckCoilHole=false;}
01911               
01912 
01913               // Then main association checks
01914               if( (Seg2->GetSeedSegment()==Seg2 && Seg2->GetNMatchSegBeg()==0) 
01915                   && (Seg2->GetBegPlane()-Seg1->GetEndPlane()>0 || (Seg1->GetEntries()>5 && Seg2->GetEntries()>5) )
01916                   && (j<npts0 || CheckCoilHole==true || (1+(Seg1->GetEndPlane()-Seg1->GetBegPlane())/2>5) ) ) {
01917                 
01918                 if( Seg1->IsAssoc(Seg2) ) {
01919 
01920                   // ND nearby 2D track protection
01922                   if( ModuleType!=2 || (fabs(Seg1->GetEndDir()-Seg2->GetBegDir())<0.25 && 
01923                                         ((Seg1->GetEntries()<11 && Seg2->GetEntries()<11 
01924                                           && (ModuleType!=2 || Seg2->GetBegPlane()-Seg1->GetEndPlane()<50)) 
01925                                          || Seg2->GetBegPlane()-Seg1->GetEndPlane()<11) ) ) {
01926 
01927                     // Leigh: Now for some additional direction protection for the ND.
01928                     if(ModuleType==2){
01929                       // Perform a linear extrapolation from the end of Seg 1 to the start of Seg 2 to check if the
01930                       // transverse positions agree. If not, then we don't have a match.
01931                       double distZ = Seg2->GetBegZPos() - Seg1->GetEndZPos();
01932                       double expectedTPos = Seg1->GetEndTPos() + Seg1->GetEndDir()*distZ;
01933                       if(fabs(expectedTPos - Seg2->GetBegTPos()) < 5*StripWidth){
01934 
01935                         TempSeg1.push_back(Seg1);
01936                         TempSeg2.push_back(Seg2);
01937                     
01938                         JoinFlag=true;
01939                     
01940                         MSG("AlgTrackCamList", Msg::kVerbose) << "   NON-ADJ ASSOC: " 
01941                                                           << Seg1->GetBegPlane() << " " << Seg1->GetEndPlane() << " (" 
01942                                                           << Seg1->GetBegTPos() << ", " << Seg1->GetEndTPos() << ") "
01943                                                           << " enddir " << Seg1->GetEndDir()
01944                                                           << " Entries " << Seg1->GetEntries() << " "
01945                                                           << Seg2->GetBegPlane() << " " << Seg2->GetEndPlane() << " (" 
01946                                                           << Seg2->GetBegTPos() << ", " << Seg2->GetEndTPos() << ") " 
01947                                                           << " begdir " << Seg2->GetBegDir() << " Entries " << Seg2->GetEntries() << endl; 
01948                                         }
01949                     }
01950                                 else{
01951                       TempSeg1.push_back(Seg1);
01952                       TempSeg2.push_back(Seg2);
01953                     
01954                       JoinFlag=true;
01955                     
01956                       MSG("AlgTrackCamList", Msg::kVerbose) << "   NON-ADJ ASSOC: " 
01957                                                           << Seg1->GetBegPlane() << " " << Seg1->GetEndPlane() << " (" 
01958                                                           << Seg1->GetBegTPos() << ", " << Seg1->GetEndTPos() << ") "
01959                                                           << " enddir " << Seg1->GetEndDir()
01960                                                           << " Entries " << Seg1->GetEntries() << " "
01961                                                           << Seg2->GetBegPlane() << " " << Seg2->GetEndPlane() << " (" 
01962                                                           << Seg2->GetBegTPos() << ", " << Seg2->GetEndTPos() << ") " 
01963                                                           << " begdir " << Seg2->GetBegDir() << " Entries " << Seg2->GetEntries() << endl; 
01964                     }
01965                   }
01967                 }
01968               }
01969             }
01970           }
01971         } 
01972       }
01973     } // End loop over entries in ViewSegBank
01974     
01975 
01976     // Make the associations
01977     for(unsigned int i=0; i<TempSeg1.size(); ++i) {
01978       TrackSegmentCam* Seg1 = TempSeg1[i];
01979       TrackSegmentCam* Seg2 = TempSeg2[i];
01980 
01981       if(Seg1->GetBegPlane()<=Seg2->GetBegPlane() && Seg1->GetEndPlane()<=Seg2->GetEndPlane()) {
01982         
01983         if( (Seg1->GetEntries()>3 && Seg2->GetEntries()>3)
01984             || 4.*fabs( (Seg1->GetBegTPos()-Seg1->GetEndTPos())/(Seg1->GetBegPlane()-Seg1->GetEndPlane()) -
01985                         (Seg2->GetBegTPos()-Seg2->GetEndTPos())/(Seg2->GetBegPlane()-Seg2->GetEndPlane()) ) < 10*StripWidth ) {
01986           Seg2->AddMatchSegToBeg(Seg1);
01987           Seg1->AddMatchSegToEnd(Seg2);
01988         }
01989       }
01990 
01991     }
01992 
01993   } // End loop over planeviews
01995 
01996 
01997 
01998   // Print out list of segments
02000   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF MATCHED SEG ASSOCIATIONS *** " << endl;
02001 
02002   for(int View=0; View<2; ++View) {
02003 
02004     MSG("AlgTrackCamList", Msg::kVerbose) << "View: " << View << endl;
02005 
02006     for(unsigned int i=0; i<ViewSegBank[View].size(); ++i) {
02007       TrackSegmentCam* segment = (TrackSegmentCam*)(ViewSegBank[View][i]);
02008       MSG("AlgTrackCamList", Msg::kVerbose) << "---  Seg0 number =" << i
02009                                             << "  begpln: "  << segment->GetBegPlane()
02010                                             << "  begtpos: " << segment->GetBegTPos()
02011                                             << "  endpln: "  << segment->GetEndPlane()
02012                                             << "  endtpos: " << segment->GetEndTPos() << endl;
02013 
02014       MSG("AlgTrackCamList", Msg::kVerbose) << "     Beg: " << endl;
02015       for(unsigned int j=0; j<segment->GetNMatchSegBeg(); ++j) {
02016         TrackSegmentCam* segtmp = segment->GetMatchSegBeg(j);
02017         MSG("AlgTrackCamList", Msg::kVerbose) << "  Match number=" << j
02018                                               << "  begpln: "  << segtmp->GetBegPlane()
02019                                               << "  begtpos: " << segtmp->GetBegTPos()
02020                                               << "  endpln: "  << segtmp->GetEndPlane()
02021                                               << "  endtpos: "  << segtmp->GetEndTPos() << endl;
02022       }
02023       
02024       MSG("AlgTrackCamList", Msg::kVerbose) << "     End: " << endl;
02025       for(unsigned int j=0; j<segment->GetNMatchSegEnd(); ++j) {
02026         TrackSegmentCam* segtmp = segment->GetMatchSegEnd(j);
02027         MSG("AlgTrackCamList", Msg::kVerbose) << "  Match number=" << j
02028                                               << "  begpln: "  << segtmp->GetBegPlane()
02029                                               << "  begtpos: " << segtmp->GetBegTPos()
02030                                               << "  endpln: "  << segtmp->GetEndPlane()
02031                                               << "  endtpos: " << segtmp->GetEndTPos() << endl;           
02032       }
02033     }
02034   }
02036 
02037 
02038 }

void AlgTrackCamList::FindPreferredJoins (  ) 

Definition at line 1078 of file AlgTrackCamList.cxx.

References TrackSegmentCam::AddPrefSegToBeg(), TrackSegmentCam::AddPrefSegToEnd(), TrackSegmentCam::GetAssocSegBeg(), TrackSegmentCam::GetAssocSegEnd(), TrackSegmentCam::GetBegPlane(), TrackSegmentCam::GetBegTPos(), TrackSegmentCam::GetEndPlane(), TrackSegmentCam::GetEndTPos(), TrackSegmentCam::GetNAssocSegBeg(), TrackSegmentCam::GetNAssocSegEnd(), TrackSegmentCam::GetNPrefSegBeg(), TrackSegmentCam::GetNPrefSegEnd(), TrackSegmentCam::GetPrefSegBeg(), TrackSegmentCam::GetPrefSegEnd(), TrackSegmentCam::GetTmpTrkFlag(), TrackSegmentCam::IsAssoc(), Msg::kVerbose, Munits::m, ModuleType, MSG, NumModules, PlanesInModule, SegmentBank, TrackSegmentCam::SetTmpTrkFlag(), and StripWidth.

Referenced by RunTheFinder().

01079 {
01080   // Having made all possible associations between triplets above, we 
01081   // select those associations that are most track-like and so are preferred.
01082   
01083   // For a given triplet, we know which triplets are associated with the
01084   // beginning and which are associated at the end. If there are associations
01085   // between these beginning and end triplets too, then we are quite likely 
01086   // to be considering track-like segments.
01087 
01088 
01089 
01090   int i;
01091 
01092   vector<TrackSegmentCam*> mTempSeg;
01093   vector<TrackSegmentCam*> pTempSeg;
01094   
01095   bool JoinFlag; bool mJnFlag; bool pJnFlag; 
01096 
01097   double GradientTolerance=8.;
01098   if(ModuleType==2) {GradientTolerance=5.;}
01099 
01100 
01101   // Loop over the planes
01103   for(int Module=0; Module<NumModules; ++Module) {
01104 
01105     for(int Plane=3; Plane<PlanesInModule-1; ++Plane) {
01106       i=Plane + Module*(PlanesInModule+1);
01107       
01108 
01109       // SET THE TEMPORARY TRACK FLAGS
01110      
01111       // Consider previous adjacent triplets   ( i-2 <-> ) i <-> i+2
01112       //                                                  Seg
01113       // If there are associations at both ends of triplet, set triplet's 
01114       // temporary track flag to 1. 
01115 
01116       // If one of beginning associated triplets also starts at lower plane 
01117       // than current triplet, set the triplet's temporary track flag to 2
01119       mTempSeg.clear();
01120 
01121       // Loop over the triplets centered on plane i
01122       for(unsigned int j=0; j<SegmentBank[i].size(); ++j) {
01123         TrackSegmentCam* Seg = SegmentBank[i][j];
01124 
01125         // Check to see if there are associations at end
01126         if(Seg->GetNAssocSegEnd()>0) {
01127           mTempSeg.push_back(Seg);
01128           
01129           // Check to see if there are associations at beginning
01130           if(Seg->GetNAssocSegBeg()>0) {
01131             
01132             // We now know that triplet has possible associations at both ends
01133             Seg->SetTmpTrkFlag(1);
01134 
01135             // Loop over the associations at beginning and set temp track flags
01136             // for segments on plane i
01137             for(unsigned int k=0; k<Seg->GetNAssocSegBeg(); ++k) {
01138               if(Seg->GetTmpTrkFlag()<2
01139                  && (Seg->GetAssocSegBeg(k)->GetBegPlane() < Seg->GetBegPlane()) ) 
01140                 {Seg->SetTmpTrkFlag(2); break;}
01141             }
01142           }
01143         }
01144       }
01146 
01147 
01148 
01149       // Subsequent adjacent triplets   i <-> i+2 ( <-> i+4 )
01150       //                                      Segp
01151       // If there are associations at both ends of triplet centered on 
01152       // plane i+2, set triplet's temporary track flag to 1.
01153 
01154       // If one of the end associated triplets also ends at higher plane
01155       // than i+2 triplet, set i+2 triplet's temporary track flag to 2
01157       pTempSeg.clear();
01158 
01159       // Loop over the triplets centered on plane i+2
01160       for(unsigned int j=0; j<SegmentBank[i+2].size(); ++j) {
01161         TrackSegmentCam* Segp = SegmentBank[i+2][j];
01162 
01163         // Check to see if there are associations at beginning
01164         if(Segp->GetNAssocSegBeg()>0) {
01165           pTempSeg.push_back(Segp);
01166           
01167           // Check to see if there are associations at end
01168           if(Segp->GetNAssocSegEnd()>0) {
01169             
01170             // We now know that triplet has possible associations at both ends
01171             Segp->SetTmpTrkFlag(1);
01172             
01173             // Loop over associations at end and set temp track flags
01174             // for segments on plane i+2
01175             for(unsigned int k=0; k<Segp->GetNAssocSegEnd(); ++k) {
01176               if(Segp->GetTmpTrkFlag()<2
01177                  && (Segp->GetAssocSegEnd(k)->GetEndPlane() > Segp->GetEndPlane()) ) 
01178                 {Segp->SetTmpTrkFlag(2); break;}
01179             }
01180           }
01181         }
01182       }
01184       
01185   
01186 
01187       // Look for preferred joins   ( i-2 <-> ) i <-> i+2
01189       // Loop the over segments on plane i+2, which we know have possible beginning associations
01190       for(unsigned int j=0; j<pTempSeg.size(); ++j) {
01191         JoinFlag=false;
01192 
01193         // Loop over these beginning associations. If one beginning associated
01194         // triplet already has temp track flag 2, set JoinFlag to true.
01195         for(unsigned int k=0; k<pTempSeg[j]->GetNAssocSegBeg(); ++k) {
01196           if(pTempSeg[j]->GetAssocSegBeg(k)->GetTmpTrkFlag()>1) {JoinFlag=true; break;}   
01197         }
01198 
01199         // If JoinFlag is true, and a beginning assoc segment has temp track flag 1 (i.e. possible 
01200         // beginning and end associations), set the temp track flag to 2 for this beg assoc segment
01201         if(JoinFlag==true) {
01202           // Loop over beginning associations again and set temp track flags
01203           // for segments on plane i
01204           for(unsigned int k=0; k<pTempSeg[j]->GetNAssocSegBeg(); ++k) {
01205             if(pTempSeg[j]->GetAssocSegBeg(k)->GetTmpTrkFlag()==1) {
01206               pTempSeg[j]->GetAssocSegBeg(k)->SetTmpTrkFlag(2); 
01207             }
01208           }
01209         }
01210       }
01212 
01213 
01214 
01215       // Look for preferred joins    i <-> i+2 ( <-> i+4 )
01217       // Loop over segments on plane i, which we know have possible end associations
01218       for(unsigned int j=0; j<mTempSeg.size(); ++j) {
01219         JoinFlag=false;
01220 
01221         // Loop over these end associations. If one end associated triplet already
01222         // has temp track flag 2, set JoinFlag to true.
01223         for(unsigned int k=0; k<mTempSeg[j]->GetNAssocSegEnd(); ++k) {
01224           if(mTempSeg[j]->GetAssocSegEnd(k)->GetTmpTrkFlag()>1) {JoinFlag=true; break;}   
01225         }
01226         
01227         // If JoinFlag is true, and an end assoc segment has temp track flag 1 (i.e. possible 
01228         // beginning and end associations), set the temp track flag to 2 for this end assoc segment
01229         if(JoinFlag==true) {
01230           // Loop over end associations again and set temp track flags
01231           // for segments on plane i+2
01232           for(unsigned int k=0; k<mTempSeg[j]->GetNAssocSegEnd(); ++k) {
01233             if(mTempSeg[j]->GetAssocSegEnd(k)->GetTmpTrkFlag()==1) {
01234               mTempSeg[j]->GetAssocSegEnd(k)->SetTmpTrkFlag(2); 
01235             }
01236           }
01237         }
01238       }
01240 
01241 
01242       // Now have temporary track flags set to 2 for all segments on planes i and i+2
01243       // which look like very promising track segments      
01245 
01246 
01247 
01248 
01249       // MAKE THE PREFERRED JOINS BETWEEN TRIPLETS ON PLANE i AND PLANE i+2
01251       if(mTempSeg.size()>0 && pTempSeg.size()>0) {
01252       
01253         // Look for doubly preferred joins    ( i-2 <-> ) i <-> i+2 ( <-> i+4 )
01255 
01256         // Loop over segments on plane i, which we know have possible end associations
01257         for(unsigned int j=0; j<mTempSeg.size(); ++j) {
01258           TrackSegmentCam* Seg = mTempSeg[j];
01259 
01260           // For these, loop over possible end associations
01261           for(unsigned int k=0; k<Seg->GetNAssocSegEnd(); ++k) {
01262             TrackSegmentCam* Segp = Seg->GetAssocSegEnd(k);
01263  
01264             // If both segments given temp track flag 2 above check assocations 
01265             // between triplets on either side
01266             if(Seg->GetTmpTrkFlag()>1 && Segp->GetTmpTrkFlag()>1) {
01267 
01268               // For each Segp, see if there is a Segm association
01269               mJnFlag=false;
01270               for(unsigned int l=0; l<Seg->GetNAssocSegBeg(); ++l) {
01271                 TrackSegmentCam* Segm = Seg->GetAssocSegBeg(l);
01272                 if(Segm->IsAssoc(Segp) ) {mJnFlag=true; break;}
01273               }
01274 
01275               // For each end assoc of Segp, see if it is assoc with Seg
01276               pJnFlag=false;
01277               for(unsigned int l=0; l<Segp->GetNAssocSegEnd(); ++l) {
01278                 TrackSegmentCam* Segp2 = Segp->GetAssocSegEnd(l);
01279                 if(Seg->IsAssoc(Segp2) ) {pJnFlag=true; break;}
01280               }
01281 
01282               // Make the preferred assocations
01283               if(mJnFlag==true && pJnFlag==true){
01284                 if(fabs((Seg->GetBegTPos()-Seg->GetEndTPos())-
01285                         (Segp->GetBegTPos()-Segp->GetEndTPos()))<GradientTolerance*StripWidth
01286                    && Seg->GetBegPlane()<=Segp->GetBegPlane()
01287                    && Seg->GetEndPlane()<=Segp->GetEndPlane() ) {
01288                   Seg->AddPrefSegToEnd(Segp); Segp->AddPrefSegToBeg(Seg);}
01289 
01290                 else {MSG("AlgTrackCamList", Msg::kVerbose) << "Seg " << Seg->GetBegPlane() << " " << Seg->GetEndPlane() << " "
01291                                                             << Seg->GetBegTPos() << " " << Seg->GetEndTPos() << " Segp " 
01292                                                             << Segp->GetBegPlane() << " " << Segp->GetEndPlane() << " "
01293                                                             << Segp->GetBegTPos() << " " << Segp->GetEndTPos() << endl;}
01294               }
01295             }
01296           }
01297         }
01299      
01300 
01301 
01302         // Look for singly preferred joins (1)    ( i-2 <-> ) i <-> i+2 
01304 
01305         // Loop over segments on plane i, which we know have possible end associations
01306         for(unsigned int j=0; j<mTempSeg.size(); ++j) {
01307           TrackSegmentCam* Seg = mTempSeg[j];
01308           
01309           // If segment given temp track flag 2 above 
01310           if(Seg->GetTmpTrkFlag()>1) {
01311             JoinFlag=false;
01312 
01313             // Check to see if we will already have made the preferred association above
01314             for(unsigned int k=0; k<Seg->GetNAssocSegEnd(); ++k) {
01315               TrackSegmentCam* Segp = Seg->GetAssocSegEnd(k);
01316               
01317               if(Segp->GetTmpTrkFlag()>1) {JoinFlag=true; break;}
01318             }
01319 
01320             // If have not made the association, can check assocations 
01321             // between triplets on either side
01322             if(JoinFlag==false) {
01323               for(unsigned int k=0; k<Seg->GetNAssocSegEnd(); ++k) {
01324                 TrackSegmentCam* Segp = Seg->GetAssocSegEnd(k);
01325 
01326                 // For each Segp, see if there is a Segm associated
01327                 mJnFlag=false;
01328                 for(unsigned int l=0; l<Seg->GetNAssocSegBeg(); ++l) {
01329                   TrackSegmentCam* Segm = Seg->GetAssocSegBeg(l);
01330                   if(mJnFlag==false && Segm->IsAssoc(Segp) ) {mJnFlag=true; break;}
01331                 }
01332                
01333                 // Make the preferred assocations
01334                 if(mJnFlag==true){
01335                   if(fabs((Seg->GetBegTPos()-Seg->GetEndTPos())-
01336                           (Segp->GetBegTPos()-Segp->GetEndTPos()))<GradientTolerance*StripWidth
01337                      && Seg->GetBegPlane()<=Segp->GetBegPlane()
01338                      && Seg->GetEndPlane()<=Segp->GetEndPlane() ) {
01339                     Seg->AddPrefSegToEnd(Segp); Segp->AddPrefSegToBeg(Seg);}
01340 
01341                   else {MSG("AlgTrackCamList", Msg::kVerbose) << "Seg " << Seg->GetBegPlane() << " " << Seg->GetEndPlane() << " "
01342                                                               << Seg->GetBegTPos() << " " << Seg->GetEndTPos() << " Segp " 
01343                                                               << Segp->GetBegPlane() << " " << Segp->GetEndPlane() << " "
01344                                                               << Segp->GetBegTPos() << " " << Segp->GetEndTPos() << endl;}
01345                 }
01346               }
01347             }
01348           }
01349         }
01351 
01352 
01353 
01354         // Look for singly preferred joins (2)    i <-> i+2 ( <-> i+4 )
01356 
01357         // Loop the over segments on plane i+2, which we know have possible beginning associations
01358         for(unsigned int j=0; j<pTempSeg.size(); ++j) {
01359           TrackSegmentCam* Segp = pTempSeg[j];
01360 
01361           // If segment given temp track flag 2 above
01362           if(Segp->GetTmpTrkFlag()>1) {
01363             JoinFlag=false;
01364 
01365             // Check to see if we will already have made the preferred association above
01366             for(unsigned int k=0; k<Segp->GetNAssocSegBeg(); ++k) {
01367               TrackSegmentCam* Seg = Segp->GetAssocSegBeg(k);
01368               
01369               if(Seg->GetTmpTrkFlag()>1) {JoinFlag=true; break;}
01370             }
01371 
01372             // If have not made the association, can check assocations 
01373             // between triplets on either side
01374             if(JoinFlag==false) {
01375               for(unsigned int k=0; k<Segp->GetNAssocSegBeg(); ++k) {
01376                 TrackSegmentCam* Seg = Segp->GetAssocSegBeg(k);
01377 
01378                 // For each Seg, see if there is a Segp2 associated
01379                 pJnFlag=false;
01380                 for(unsigned int l=0; l<Segp->GetNAssocSegEnd(); ++l) {
01381                   TrackSegmentCam* Segp2 = Segp->GetAssocSegEnd(l);
01382                   if(pJnFlag==false && Seg->IsAssoc(Segp2) ) {pJnFlag=true; break;}
01383                 }
01384                
01385                 // Make the preferred assocations
01386                 if(pJnFlag==true){
01387                   if(fabs((Seg->GetBegTPos()-Seg->GetEndTPos())-
01388                           (Segp->GetBegTPos()-Segp->GetEndTPos()))<GradientTolerance*StripWidth
01389                      && Seg->GetBegPlane()<=Segp->GetBegPlane()
01390                      && Seg->GetEndPlane()<=Segp->GetEndPlane() ) {
01391                   Seg->AddPrefSegToEnd(Segp); Segp->AddPrefSegToBeg(Seg);}
01392 
01393                   else {MSG("AlgTrackCamList", Msg::kVerbose) << "Seg " << Seg->GetBegPlane() << " " << Seg->GetEndPlane() << " "
01394                                                               << Seg->GetBegTPos() << " " << Seg->GetEndTPos() << " Segp " 
01395                                                               << Segp->GetBegPlane() << " " << Segp->GetEndPlane() << " "
01396                                                               << Segp->GetBegTPos() << " " << Segp->GetEndTPos() << endl;}
01397                 }
01398               }
01399             }
01400           }
01401         }
01403 
01404 
01405 
01406         // Look for other joins we have missed    i <-> i+2 
01408 
01409         // Loop over segments on plane i, which we know have possible end associations
01410         for(unsigned int j=0; j<mTempSeg.size(); ++j) {
01411           TrackSegmentCam* Seg = mTempSeg[j];
01412           
01413           // For these, loop over the end assocations
01414           for(unsigned int k=0; k<Seg->GetNAssocSegEnd(); ++k) {
01415             TrackSegmentCam* Segp = Seg->GetAssocSegEnd(k);
01416             
01417             // If segments and its end assoc segment were given temp track flag 2 above
01418             if(Seg->GetTmpTrkFlag()<2 && Segp->GetTmpTrkFlag()<2) {
01419               
01420               // Look at each possible end assoc for Seg to see if we've already
01421               // made the assocation
01422               mJnFlag=false;
01423               for(unsigned int l=0; l<Seg->GetNAssocSegEnd(); ++l) {
01424                 if(Seg->GetAssocSegEnd(l)->GetTmpTrkFlag()>1) {mJnFlag=true; break;}
01425               }
01426 
01427               // Look at each possible beginning assoc for Segp to see if we've already
01428               // made the assocation
01429               pJnFlag=false;
01430               for(unsigned int l=0; l<Segp->GetNAssocSegBeg(); ++l) {
01431                 if(Segp->GetAssocSegBeg(l)->GetTmpTrkFlag()>1) {pJnFlag=true; break;}
01432               }
01433                       
01434               // Make the preferred associations if it hasn't been done by one of the loops above
01435               if(mJnFlag==false && pJnFlag==false) {
01436                 if(fabs((Seg->GetBegTPos()-Seg->GetEndTPos())-
01437                         (Segp->GetBegTPos()-Segp->GetEndTPos()))<GradientTolerance*StripWidth
01438                    && Seg->GetBegPlane()<=Segp->GetBegPlane()
01439                    && Seg->GetEndPlane()<=Segp->GetEndPlane() ) {
01440                 Seg->AddPrefSegToEnd(Segp); Segp->AddPrefSegToBeg(Seg);}
01441 
01442                 else {MSG("AlgTrackCamList", Msg::kVerbose) << "Seg " << Seg->GetBegPlane() << " " << Seg->GetEndPlane() << " "
01443                                                             << Seg->GetBegTPos() << " " << Seg->GetEndTPos() << " Segp " 
01444                                                             << Segp->GetBegPlane() << " " << Segp->GetEndPlane() << " "
01445                                                             << Segp->GetBegTPos() << " " << Segp->GetEndTPos() << endl;}
01446               }
01447             }       
01448           }
01449         }
01450         
01452 
01453 
01454       }
01456 
01457 
01458 
01459       // Clear flags ready to consider next planes
01461       for(unsigned int j=0; j<mTempSeg.size(); ++j) {
01462         mTempSeg[j]->SetTmpTrkFlag(0);
01463       }
01464       
01465       for(unsigned int j=0; j<pTempSeg.size(); ++j) {
01466         pTempSeg[j]->SetTmpTrkFlag(0);
01467       }
01469 
01470 
01471     } // End loop over planes in module
01472 
01473   } // End loop over modules
01475 
01476   
01477   
01478   
01479   // Look for any other preferred associations we may have missed.
01480   // Find segments with no preferred end association and look a bit
01481   // further to find preferred associations.
01482 
01483   // Not for ND
01484   if(ModuleType==1 || ModuleType==3) {
01485 
01486     int NewPlane; int Increment;
01487     double SegTPos, NearbySegTPos; 
01488     bool AssocsStopHere;
01489     bool ConsiderSegment;
01490 
01491     for(int Module=0; Module<NumModules; ++Module) {
01492       for(int Plane=3; Plane<PlanesInModule-1; ++Plane) {
01493         i=Plane + Module*(PlanesInModule+1);
01494           
01495         // Loop over all segments
01496         for(unsigned int j=0; j<SegmentBank[i].size(); ++j) {
01497           TrackSegmentCam* Seg0 = SegmentBank[i][j];
01498         
01499           // Loop over the segments with no preferred end association
01501           if(Seg0->GetNPrefSegEnd()==0 && Seg0->GetNPrefSegBeg()>0) {
01502           
01503             // Look at nearby segments to see if these already continue the associations.
01505             AssocsStopHere=true;
01506             SegTPos=0.5*(Seg0->GetBegTPos()+Seg0->GetEndTPos());
01507 
01508             for(unsigned int k=0; k<SegmentBank[i].size(); ++k) {
01509               TrackSegmentCam* NearbySeg = SegmentBank[i][k];
01510             
01511               if(NearbySeg==Seg0) {continue;}
01512             
01513               NearbySegTPos=0.5*(NearbySeg->GetBegTPos()+NearbySeg->GetEndTPos());
01514 
01515               if(fabs(NearbySegTPos-SegTPos)<1 && 
01516                  NearbySeg->GetNPrefSegEnd()>0 && NearbySeg->GetNPrefSegBeg()>0) {AssocsStopHere=false; break;}
01517             }
01518 
01519             if(AssocsStopHere==false) {break;}
01521 
01522 
01523             JoinFlag=false;
01524 
01525             // Consider the segment, plus all the segments, SegBeg, with which 
01526             // it has a preferred beginning association.
01528             for(int k=-1; k<int(Seg0->GetNPrefSegBeg()); ++k) {
01529               // First time, consider Seg0, rather than one of its preferred beginning associations.
01530               TrackSegmentCam* SegBeg = 0;
01531               if(k<0) {SegBeg=Seg0;}
01532               else{SegBeg = Seg0->GetPrefSegBeg(k);}
01533               Increment=4;
01534 
01535               // Look up to 12 planes ahead within the same module.
01536               while(Increment<=14) {
01537                 NewPlane=SegBeg->GetEndPlane()+Increment;
01538               
01539                 if(NewPlane>=(Module+1)*(PlanesInModule+1)) {break;}
01540 
01541                 // See if SegBeg is associated with a segment on this plane, NextSeg.
01542                 // If so, also see if SegBeg is associated with one of NextSeg's 
01543                 // preferred end associated segments.
01544 
01545                 // If so, make the preferred associations between SegBeg and NextSeg.
01546                 for(unsigned int l=0; l<SegmentBank[NewPlane].size(); ++l) {
01547                   TrackSegmentCam* NextSeg = SegmentBank[NewPlane][l];
01548 
01549 
01550                   // Check segments are not now linked in a chain of preferred associations.
01552                   ConsiderSegment=true;
01553                   for(unsigned int m=0; m<NextSeg->GetNPrefSegBeg(); ++m) {
01554                     TrackSegmentCam* PrefSeg = NextSeg->GetPrefSegBeg(m);
01555 
01556                     if(PrefSeg->GetTmpTrkFlag()==1) {NextSeg->SetTmpTrkFlag(1); ConsiderSegment=false; break;}
01557                   }
01558                   if(ConsiderSegment==false) {continue;}
01560 
01561 
01563                   if( SegBeg->IsAssoc(NextSeg) ) {
01564                   
01565                     for(unsigned int m=0; m<NextSeg->GetNPrefSegEnd(); ++m) {
01566                       TrackSegmentCam* SegEnd = NextSeg->GetPrefSegEnd(m);
01567                     
01568                       if( SegBeg->IsAssoc(SegEnd) ) {
01569                         JoinFlag=true;
01570                       
01571                         SegBeg->AddPrefSegToEnd(NextSeg);
01572                         NextSeg->AddPrefSegToBeg(SegBeg);
01573 
01574                         SegBeg->SetTmpTrkFlag(1); NextSeg->SetTmpTrkFlag(1);
01575 
01576                         MSG("AlgTrackCamList", Msg::kVerbose) << "Made missing Pref assoc, SegBeg "
01577                                                               << SegBeg->GetBegPlane() << ", " << SegBeg->GetEndPlane() << ", "
01578                                                               << SegBeg->GetBegTPos() << ", " << SegBeg->GetEndTPos()
01579                                                               << " NextSeg " << NextSeg->GetBegPlane() << ", " 
01580                                                               << NextSeg->GetEndPlane() << ", " << NextSeg->GetBegTPos() 
01581                                                               << ", " << NextSeg->GetEndTPos() << endl;
01582                         break;
01583                       }
01584                     }
01585                   }
01587                 }
01588                 Increment+=2;
01589               }
01590 
01591               // Reset TmpTrkFlags
01592               for(int p=i; p<(1+Module)*PlanesInModule; ++p) {
01593                 for(unsigned int q=0; q<SegmentBank[p].size(); ++q) {SegmentBank[p][q]->SetTmpTrkFlag(0);}
01594               }
01595 
01596               if(JoinFlag==true) {break;}
01597             }
01598           }
01600         
01601 
01602         
01603           // Now loop over the segments with no preferred beginning association
01605           if(Seg0->GetNPrefSegBeg()==0 && Seg0->GetNPrefSegEnd()>0) {
01606           
01607             // Look at nearby segments to see if these already continue the associations.
01609             AssocsStopHere=true;
01610             SegTPos=0.5*(Seg0->GetBegTPos()+Seg0->GetEndTPos());
01611 
01612             for(unsigned int k=0; k<SegmentBank[i].size(); ++k) {
01613               TrackSegmentCam* NearbySeg = SegmentBank[i][k];
01614             
01615               if(NearbySeg==Seg0) {continue;}
01616             
01617               NearbySegTPos=0.5*(NearbySeg->GetBegTPos()+NearbySeg->GetEndTPos());
01618 
01619               if(fabs(NearbySegTPos-SegTPos)<1 && 
01620                  NearbySeg->GetNPrefSegEnd()>0 && NearbySeg->GetNPrefSegBeg()>0) {AssocsStopHere=false; break;}
01621             }
01622 
01623             if(AssocsStopHere==false) {break;}
01625 
01626 
01627             JoinFlag=false;
01628 
01629             // Consider the segment, plus all the segments, SegEnd, with which 
01630             // it has a preferred end association.
01632             for(int k=-1; k<int(Seg0->GetNPrefSegEnd()); ++k) {
01633               // First time, consider Seg0, rather than one of its preferred end associations.
01634               TrackSegmentCam* SegEnd = 0;
01635               if(k<0) {SegEnd=Seg0;}
01636               else{SegEnd = Seg0->GetPrefSegEnd(k);}
01637               Increment=4;
01638 
01639               // Look up to 12 planes behind within the same module.
01640               while(Increment<=14) {
01641                 NewPlane=SegEnd->GetBegPlane()-Increment;
01642               
01643                 if(NewPlane<=(Module)*(PlanesInModule+1)) {break;}
01644 
01645                 // See if SegEnd is associated with a segment on this plane, PrevSeg.
01646                 // If so, also see if SegEnd is associated with one of PrevSeg's 
01647                 // preferred beginning associated segments.
01648 
01649                 // If so, make the preferred associations between SegEnd and PrevSeg.
01650                 for(unsigned int l=0; l<SegmentBank[NewPlane].size(); ++l) {
01651                   TrackSegmentCam* PrevSeg = SegmentBank[NewPlane][l];
01652 
01653 
01654                   // Check segments are not now linked in a chain of preferred associations.
01656                   ConsiderSegment=true;
01657                   for(unsigned int m=0; m<PrevSeg->GetNPrefSegEnd(); ++m) {
01658                     TrackSegmentCam* PrefSeg = PrevSeg->GetPrefSegEnd(m);
01659 
01660                     if(PrefSeg->GetTmpTrkFlag()==1) {PrevSeg->SetTmpTrkFlag(1); ConsiderSegment=false; break;}
01661                   }
01662                   if(ConsiderSegment==false) {continue;}
01664         
01665 
01667                   if( PrevSeg->IsAssoc(SegEnd) ) {
01668                   
01669                     for(unsigned int m=0; m<PrevSeg->GetNPrefSegBeg(); ++m) {
01670                       TrackSegmentCam* SegBeg = PrevSeg->GetPrefSegBeg(m);
01671                             
01672                       if( SegBeg->IsAssoc(SegEnd) ) {
01673                         JoinFlag=true;
01674                       
01675                         SegEnd->AddPrefSegToBeg(PrevSeg);
01676                         PrevSeg->AddPrefSegToEnd(SegEnd);
01677 
01678                         SegEnd->SetTmpTrkFlag(1); PrevSeg->SetTmpTrkFlag(1);
01679 
01680                         MSG("AlgTrackCamList", Msg::kVerbose) << "Made missing Pref assoc, SegEnd "
01681                                                               << SegEnd->GetBegPlane() << ", " << SegEnd->GetEndPlane() << ", "
01682                                                               << SegEnd->GetBegTPos() << ", " << SegEnd->GetEndTPos()
01683                                                               << " PrevSeg " << PrevSeg->GetBegPlane() << ", " 
01684                                                               << PrevSeg->GetEndPlane() << ", " << PrevSeg->GetBegTPos() 
01685                                                               << ", " << PrevSeg->GetEndTPos() << endl;
01686                         break;
01687                       }
01688                     }
01689                   }
01691                 }
01692                 Increment+=2;
01693               }
01694 
01695               // Reset TmpTrkFlags
01696               for(int p=i; p>Module*(PlanesInModule); --p) {
01697                 for(unsigned int q=0; q<SegmentBank[p].size(); ++q) {SegmentBank[p][q]->SetTmpTrkFlag(0);}
01698               }
01699 
01700               if(JoinFlag==true) {break;}
01701             }
01702           }
01704 
01705         }
01706       }
01707     }
01708 
01709   }
01711  
01712   
01713 
01714 
01715   // Print out list of preferred associations
01717   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF PREF SEG ASSOCIATIONS *** " << endl;
01718 
01719   // Print out list of triplets
01720   for(int i=0; i<500; ++i) {
01721     for(unsigned int j=0;j<SegmentBank[i].size(); ++j) {
01722       TrackSegmentCam* segment = SegmentBank[i][j];
01723 
01724       MSG("AlgTrackCamList", Msg::kVerbose) << " ----Plane " << i << ", Seg Number " << j 
01725                                             << ", BegPlane " << segment->GetBegPlane() 
01726                                             << ", EndPlane " << segment->GetEndPlane() 
01727                                             << ", BegTPos " << segment->GetBegTPos() 
01728                                             << ", EndTPos " << segment->GetEndTPos() 
01729                                             <<endl;
01730  
01731       MSG("AlgTrackCamList", Msg::kVerbose) << " Beg: " << endl;
01732       for(unsigned int k=0; k<segment->GetNPrefSegBeg(); ++k) {
01733         TrackSegmentCam* segtmp = segment->GetPrefSegBeg(k);
01734         MSG("AlgTrackCamList", Msg::kVerbose) << "  Pref number=" << k
01735                                               << "  begpln: " << segtmp->GetBegPlane()
01736                                               << "  begtpos: " << segtmp->GetBegTPos()
01737                                               << "  endpln: " << segtmp->GetEndPlane()
01738                                               << "  endtpos: " << segtmp->GetEndTPos() << endl;
01739       }
01740       
01741       MSG("AlgTrackCamList", Msg::kVerbose) << " End: " << endl;
01742       for(unsigned int k=0; k<segment->GetNPrefSegEnd(); ++k) {
01743         TrackSegmentCam* segtmp = segment->GetPrefSegEnd(k);
01744         MSG("AlgTrackCamList", Msg::kVerbose) << "  Pref number=" << k
01745                                               << "  begpln: " << segtmp->GetBegPlane()
01746                                               << "  begtpos: " << segtmp->GetBegTPos()
01747                                               << "  endpln: " << segtmp->GetEndPlane()
01748                                               << "  endtpos: " << segtmp->GetEndTPos() << endl;           
01749       }
01750     }
01751   }
01753   
01754   return;
01755 }

void AlgTrackCamList::FirstUVComparison (  ) 

Definition at line 2046 of file AlgTrackCamList.cxx.

References CambridgeAnalysis, ClusterBank, TrackSegmentCam::GetBegPlane(), ClusterCam::GetBegTPos(), TrackSegmentCam::GetBegTPos(), TrackSegmentCam::GetCluster(), VldContext::GetDetector(), TrackSegmentCam::GetEndPlane(), ClusterCam::GetEndTPos(), TrackSegmentCam::GetEndTPos(), TrackSegmentCam::GetEntries(), ClusterCam::GetNDFlag(), ClusterCam::GetShwFlag(), ClusterCam::GetTrkFlag(), ClusterCam::GetTrkPlnFlag(), TrackSegmentCam::GetUID(), HitBank, Detector::kNear, Msg::kVerbose, MSG, PlanesInModule, ClusterCam::SetNDFlag(), TrackSegmentCam::SetUID(), StripWidth, ViewSegBank, and vldc.

Referenced by RunTheFinder().

02047 {
02048   // Having made sections of the track, from 'matched' segments, we compare
02049   // the sections in the U view with those in the V view, looking for compatibity.
02050 
02051   // First all clusters are investigated to mark the segments as track-like or
02052   // shower-like. Then track-like segments which have 'overlapping' partners in the 
02053   // opposite view are marked by setting their UID values to 2
02054 
02055   MSG("AlgTrackCamList", Msg::kVerbose) << " *** Comparing views (1) *** " << endl;
02056 
02057 
02058   int shwctr, plnctr;
02059   unsigned int nsegments=0;
02060   unsigned int nsegmentsV=0;
02061   vector<TrackSegmentCam*> mysegbnk[2];
02062 
02063 
02064   // Determine track-like segments
02066   for(int View=0; View<2; ++View) {
02067 
02068     // Loop over all segments in the current view
02069     nsegments = ViewSegBank[View].size();
02070     for(unsigned int i=0; i<nsegments; ++i) {
02071 
02072       TrackSegmentCam* Seg = ViewSegBank[View][i];
02073       shwctr=0; plnctr=0;
02074       
02075       // Loop over all the clusters in the current segment
02076       const unsigned int ncluster = Seg->GetEntries();
02077       for(unsigned int j=0;j<ncluster;++j) {
02078         ClusterCam* cluster = Seg->GetCluster(j);
02079         
02080         // Find how track-like or shower-like the cluster is
02081         if(cluster->GetTrkFlag()>1) { 
02082           if(cluster->GetShwFlag()<1) {shwctr++;}
02083         }
02084         if(cluster->GetTrkPlnFlag()>0) {plnctr++;}
02085       }
02086       
02087       // Set the UID accordingly. This depends greatly on the initial values 
02088       // of the shower flags used when identifying track and shower clusters.
02089       if(CambridgeAnalysis==false) {Seg->SetUID(1);}
02090 
02091       if( plnctr>1 && Seg->GetEntries()>4 ) { // Just added
02092         if(CambridgeAnalysis==true && shwctr>0) {Seg->SetUID(1);} 
02093         if(shwctr>2) {Seg->SetUID(2);} 
02094       }
02095     }
02096   }
02098 
02099 
02100   // Loop over segments in U view
02102   nsegments = ViewSegBank[0].size();
02103   for(unsigned int i=0; i<nsegments; ++i) {
02104 
02105     TrackSegmentCam* Segu = ViewSegBank[0][i];
02106     
02107     // For these, loop over segments in V view
02108     nsegmentsV = ViewSegBank[1].size();
02109     for(unsigned int j=0; j<nsegmentsV; ++j) {
02110 
02111       TrackSegmentCam* Segv = ViewSegBank[1][j];
02112 
02113       // Select pairs of segments, one from each view, which overlap by more than 1 plane
02114       // and which both have UID>0 (i.e. were marked as track-like above)
02115       if( Segu->GetBegPlane()+1<Segv->GetEndPlane() 
02116           && Segu->GetEndPlane()>Segv->GetBegPlane()+1 
02117           && Segu->GetUID()>0 && Segv->GetUID()>0 ) {
02118         if(Segv->GetUID()>1) {mysegbnk[0].push_back(Segu); }
02119         if(Segu->GetUID()>1) {mysegbnk[1].push_back(Segv); }
02120       }
02121     }
02122   }
02124 
02125   
02126   // Loop over the segments we just selected and set their UIDs to 2
02128   for(int View=0; View<2; ++View) {
02129     nsegments = mysegbnk[View].size();
02130     for(unsigned int i=0; i<nsegments; ++i) {
02131       mysegbnk[View][i]->SetUID(2);
02132     }
02133   }
02135 
02136 
02137 
02138   // Mark clusters that aren't going to be available for ND +/-10 plane track finding
02140   if(vldc->GetDetector()==Detector::kNear) {
02141     
02142     // First mark all those for which we will probably find a +/-2 plane track
02144     for(int View=0; View<2; ++View) {
02145       for(unsigned int i=0; i<mysegbnk[View].size(); ++i) {
02146         TrackSegmentCam* Seg = mysegbnk[View][i];
02147 
02148         for(unsigned int j=0; j<Seg->GetEntries(); ++j) {
02149           Seg->GetCluster(j)->SetNDFlag(0);  
02150         }
02151       }
02152     }
02154 
02155 
02156     // Then look specifically for clusters that are in +/-2 plane sections of detector
02158     double BegTPos, EndTPos, Centre, HitTPos;
02159     double Tolerance=2.*StripWidth;
02160 
02161     // Loop over all clusters
02162     for(int i=0; i<PlanesInModule; ++i) {
02163       for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
02164         ClusterCam* Clust = ClusterBank[i][j];
02165 
02166         // Get cluster properties       
02167         BegTPos=Clust->GetBegTPos();
02168         EndTPos=Clust->GetEndTPos();
02169         Centre=0.5*(BegTPos+EndTPos);
02170 
02171 
02172         // Now compare with the hits +/-2 planes away
02173         if(Clust->GetNDFlag()>0 && i>1) {
02174           for(unsigned int k=0; k<HitBank[i-2].size(); ++k) {
02175             HitTPos=HitBank[i-2][k]->GetTPos();
02176 
02177             if(fabs(HitTPos-Centre)<Tolerance || fabs(HitTPos-BegTPos)<Tolerance 
02178                || fabs(HitTPos-EndTPos)<Tolerance) {Clust->SetNDFlag(0); break;}
02179           }
02180         }
02181 
02182         if(Clust->GetNDFlag()>0) {
02183           for(unsigned int k=0; k<HitBank[i+2].size(); ++k) {
02184             HitTPos=HitBank[i+2][k]->GetTPos();
02185             
02186             if(fabs(HitTPos-Centre)<Tolerance || fabs(HitTPos-BegTPos)<Tolerance 
02187                || fabs(HitTPos-EndTPos)<Tolerance) {Clust->SetNDFlag(0); break;}
02188           }
02189         }
02190       }
02191     }
02193 
02194   } // End ND only section
02196 
02197 
02198   // Print out the results
02200   for(int View=0; View<2; ++View) {
02201     MSG("AlgTrackCamList", Msg::kVerbose) << "  View =" << View << endl;
02202 
02203     nsegments = ViewSegBank[View].size();
02204     for(unsigned int i=0; i<nsegments; ++i) {
02205 
02206       TrackSegmentCam* Seg = ViewSegBank[View][i];
02207 
02208       MSG("AlgTrackCamList", Msg::kVerbose) << "  Seg number=" << i 
02209                                             << "  begpln: "  << Seg->GetBegPlane()
02210                                             << "  begtpos: " << Seg->GetBegTPos()
02211                                             << "  endpln: "  << Seg->GetEndPlane()
02212                                             << "  endtpos: " << Seg->GetEndTPos() 
02213                                             << "  UID: "     << Seg->GetUID() 
02214                                             << "  clusters: "<< Seg->GetEntries()
02215                                             << endl;
02216     }
02217   }
02219   
02220   return;
02221 }

void AlgTrackCamList::Form2DTracks (  ) 

Definition at line 2746 of file AlgTrackCamList.cxx.

References TrackSegmentCam::AddSegment(), TrackSegmentCam::GetBegPlane(), TrackSegmentCam::GetBegTPos(), TrackSegmentCam::GetCluster(), TrackSegmentCam::GetEndPlane(), TrackSegmentCam::GetEndTPos(), TrackSegmentCam::GetEntries(), TrackSegmentCam::GetMatchSegBeg(), TrackSegmentCam::GetMatchSegEnd(), TrackSegmentCam::GetNMatchSegBeg(), TrackSegmentCam::GetNMatchSegEnd(), TrackSegmentCam::GetScore(), TrackSegmentCam::GetTmpTrkFlag(), ClusterCam::GetTrkFlag(), ClusterCam::GetTrkPlnFlag(), TrackSegmentCam::GetUID(), id, Msg::kVerbose, Msg::kWarning, MSG, PossibleJoins, TrackSegmentCam::SetTmpTrkFlag(), ClusterCam::SetTrkFlag(), TrackSegmentCam::SetUID(), size, and ViewSegBank.

Referenced by RunTheFinder().

02747 {
02748   // Of the segments identified as good above, we identify the segments which
02749   // are good seeds for the track i.e. from which we can propagate back and forth
02750   // along matched segments to find a long track.
02751  
02752   // Then, for the first seed segment, we try to propagate backwards and forwards,
02753   // marking the segments we use with different TmpTrkFlag settings.
02754 
02755   // For a seed segment, Seg0, we firstly move in the direction of increasing plane
02756   // number. We mark all the segments we use with TmpTrkFlag 1. The seed segment is
02757   // labelled with TmpTrkFlag 3.
02758 
02759   // Some paths will lead further than others. We are interested in the longest paths, 
02760   // and move back towards Seg0 along the longest paths, marking the segments used 
02761   // with TmpTrkFlag 2.
02762 
02763   // Having done this, we can carry out a similar procedure, but initially moving 
02764   // backwards from Seg0 to lower plane numbers, before returning to Seg0 again.
02765 
02766   // So, this arrangement of segments...
02767   //
02768   //       Seg  Seg                                  Seg
02769   //                  Seg                      Seg
02770   //                        Seg   Seg0   Seg 
02771   //                  Seg                      Seg   Seg   Seg
02772   //                                                 
02773   //               Seg                               Seg
02774   //                                                       Seg   Seg
02775 
02776   // Is labelled as follows...
02777   //
02778   //       2     2                                    1 
02779   //                   2                        1 
02780   //                         2     3      2  
02781   //                   1                        2     2     2 
02782   //                                                 
02783   //                1                                 2 
02784   //                                                        2     2 
02785   //
02786 
02787   // It is possible that there are multiple paths along segments with TmpTrkFlag 2, 
02788   // as above. The best path is selected by obtaining a score for each, based on
02789   // straightness and length.
02790 
02791 
02792 
02793   // Define containers
02794   vector<TrackSegmentCam*> BestSeedSegments;
02795   vector<TrackSegmentCam*> SeedSegments;
02796 
02797   vector<TrackSegmentCam*> Temp;
02798   vector<TrackSegmentCam*> mTemp;
02799   vector<TrackSegmentCam*> pTemp;
02800   vector<TrackSegmentCam*> BegTemp;
02801   vector<TrackSegmentCam*> EndTemp;
02802   vector<TrackSegmentCam*> TempBank1;
02803   vector<TrackSegmentCam*> TempBank2;
02804 
02805   vector <vector<TrackSegmentCam*> > BegBank;
02806   vector <vector<TrackSegmentCam*> > EndBank;
02807 
02808   vector <vector<TrackSegmentCam*> > TempBegBank;
02809   vector <vector<TrackSegmentCam*> > TempEndBank;
02810 
02811 
02812   bool Cont; bool tmpFlag; bool TrackFlag;
02813   int ntrks; int Counter; int nplane; int npts;
02814   int ObjCounter; int ObjVectorCounter; const int fObjVectorMax=100000;
02815   int id; double Score; double TopScore;
02816 
02817   unsigned int MostClusters;
02818   bool AlreadyAdded;
02819 
02820 
02821   // Main loop is over the views, first for U, then V
02824   if(ViewSegBank[0].size()>0 && ViewSegBank[1].size()>0) {
02825     for(int View=0; View<2; ++View) {
02826     
02827       Cont=true; ntrks=0;
02828 
02829       while(Cont==true) {
02830         Cont=false;
02831         BegBank.clear(); EndBank.clear();
02832 
02833 
02834 
02835         // First selection of SEED SEGMENTS - segment from which we can move
02836         // along matched associations to find majority of the track.
02837         //
02838         // First guesses at seed segments are stored in SeedSegments Container
02840         BestSeedSegments.clear(); SeedSegments.clear();
02841 
02842         // Loop over the entries in ViewSegBank
02843         for(unsigned int i=0; i<ViewSegBank[View].size(); ++i) {
02844           TrackSegmentCam* Seg = ViewSegBank[View][i];
02845 
02846           if(Seg->GetUID()>0 && Seg->GetUID()<3) {
02847             Counter=0;
02848             
02849             // Loop over the constituent clusters and count the number that are track-like
02850             for(unsigned int j=0; j<Seg->GetEntries(); ++j) {
02851               ClusterCam* Clust = Seg->GetCluster(j);
02852               if(Clust->GetTrkPlnFlag()>0 && Clust->GetTrkFlag()==2) {Counter++;}
02853             }
02854             
02855             // If we've already made a track and this segment doesn't contain 3 track-like
02856             // clusters, don't consider it any further
02857             if(ntrks>0 && Counter<3) {Seg->SetUID(0);}
02858           }
02859 
02860           // Temporarily store all the likely seed segments. Those with UID 2 are most track-like.
02861           if(Seg->GetUID()==2) {BestSeedSegments.push_back(Seg);}
02862           if(Seg->GetUID()==1) {SeedSegments.push_back(Seg);}
02863         }
02864 
02865 
02866         // Hopefully we have segments with UID 2 (most track-like), else take those with UID 1
02867         // (just track-like) and move them from SeedSegments to BestSeedSegments. 
02868         // Can then empty SeedSegments.
02869         if(BestSeedSegments.size()==0) {BestSeedSegments=SeedSegments;}
02870         
02871         SeedSegments.clear();
02872 
02873         SeedSegments=BestSeedSegments;
02875 
02876 
02877 
02878 
02879         // Now, using our initial seed segments, try propagating back and forth
02880         // from each segment to refine our choice.
02881 
02882         // Final SEED SEGMENTS are stored BestSeedSegments Container
02884         BestSeedSegments.clear();
02885 
02886 
02887         // If we have some seed segments to work with, loop over them
02888         if(SeedSegments.size()>0) { 
02889           nplane=-1;
02890 
02891           for(unsigned int i=0; i<SeedSegments.size(); ++i) {
02892 
02893             // See how far we can propagate backwards to lower plane numbers
02895             mTemp.clear(); BegTemp.clear(); 
02896 
02897             // First put seed segment into BegTemp then begin while loop
02898             BegTemp.push_back(SeedSegments[i]); TempBank1.clear(); 
02899 
02900             while(BegTemp.size()>0) {
02901               Temp.clear();
02902               
02903               // Copy any segments into a temporary holder and loop over them,
02904               // gradually propagating back to lower plane numbers
02905               Temp=BegTemp; BegTemp.clear();
02906             
02907               for(unsigned int j=0; j<Temp.size(); ++j) {
02908                 tmpFlag=false;
02909                 
02910                 // For each segment, loop over their beginning matches
02911                 for(unsigned int k=0; k<Temp[j]->GetNMatchSegBeg(); ++k) {
02912                   TrackSegmentCam* Segb = Temp[j]->GetMatchSegBeg(k);
02913                   
02914                   // If matched segment is track-like and has zero TmpTrkFlag, 
02915                   // push it back into BegTemp, so the while loop continues.
02916                   if(Segb->GetUID()>-1 && Segb->GetUID()<3) {
02917                     tmpFlag=true;
02918                     
02919                     if(Segb->GetTmpTrkFlag()<1) {
02920                       BegTemp.push_back(Segb);
02921                       Segb->SetTmpTrkFlag(1); // Temporarily mark segments used
02922                       TempBank1.push_back(Segb);
02923                     }
02924                   }
02925                 }
02926                 
02927                 // Must have gone as far as we can go with this seed segment. Store segment at
02928                 // lowest plane.
02929                 if(tmpFlag==false) {mTemp.push_back(Temp[j]);}
02930               }
02931             } // End while BegTemp.size()>0
02932 
02933             // Set the flags back to zero, for when we carry out actual propagation, below.
02934             for(unsigned int j=0; j<TempBank1.size(); ++j) {TempBank1[j]->SetTmpTrkFlag(0);}
02936 
02937 
02938             
02939             // Now see how far we can propagate forwards to higher plane numbers
02941             pTemp.clear(); EndTemp.clear();
02942 
02943             // Put the seed segment into EndTemp then start while loop
02944             EndTemp.push_back(SeedSegments[i]); TempBank1.clear();
02945             
02946             while(EndTemp.size()>0) {
02947               Temp.clear();
02948               
02949               // Copy any segments into a temporary holder and loop over them,
02950               // gradually propagating to higher plane numbers
02951               Temp=EndTemp; EndTemp.clear();
02952              
02953               for(unsigned int j=0; j<Temp.size(); ++j) {
02954                 tmpFlag=false;
02955                 
02956                 // For each segment, loop over their end matches
02957                 for(unsigned int k=0; k<Temp[j]->GetNMatchSegEnd(); ++k) {
02958                   TrackSegmentCam* Sege = Temp[j]->GetMatchSegEnd(k);
02959 
02960                   // If the matched segment is track-like and has zero TmpTrkFlag, 
02961                   // push back into EndTemp, so the while loop continues.
02962                   if(Sege->GetUID()>-1 && Sege->GetUID()<3) {
02963                     tmpFlag=true;
02964                     
02965                     if(Sege->GetTmpTrkFlag()<1) {
02966                       EndTemp.push_back(Sege);
02967                       Sege->SetTmpTrkFlag(1); // Temporarily mark segments used
02968                       TempBank1.push_back(Sege);
02969                     }
02970                   }
02971                 }
02972                 
02973                 // Must have gone as far as we can go with this seed segment. Store segment at
02974                 // highest plane.
02975                 if(tmpFlag==false) {pTemp.push_back(Temp[j]);}
02976               }
02977             } // End while EndTemp.size()>0
02978             
02979             // Set the flags back to zero, for when we carry out actual propagation, below.
02980             for(unsigned int j=0; j<TempBank1.size(); ++j) {TempBank1[j]->SetTmpTrkFlag(0);}
02982             
02983             
02984             
02985             // Find the maximum number of planes we can span by propagating the seed segment
02986             // back and forth.
02988             npts=-1;
02989             for(unsigned int j=0; j<mTemp.size(); ++j) {
02990               for(unsigned int k=0; k<pTemp.size(); ++k) {
02991                 if((1+pTemp[k]->GetEndPlane()-mTemp[j]->GetBegPlane())>npts) {
02992                   npts=1+pTemp[k]->GetEndPlane()-mTemp[j]->GetBegPlane();
02993                 }
02994               }
02995             }
02996             
02997             SeedSegments[i]->SetNPlanes(npts);
02998             if(npts>nplane) {nplane=npts;}
03000           
03001           } // End loop over seed segments stored above
03002 
03003 
03004           // Store the seed segments that lead to the biggest propagation in BestSeedSegments
03006           if(nplane>0) {
03007             for(unsigned int i=0; i<SeedSegments.size(); ++i) {
03008               if(SeedSegments[i]->GetNPlanes()>nplane-5) {BestSeedSegments.push_back(SeedSegments[i]);}       
03009             }
03010           }
03012           
03013           
03014         } // End finding final seed segments
03016         
03017 
03018         // Order contents of BestSeedSegments
03020         SeedSegments=BestSeedSegments; BestSeedSegments.clear();
03021         
03022         for(unsigned int i=0; i<SeedSegments.size(); ++i) {
03023           MostClusters=0;
03024           TrackSegmentCam* LargestSeg = 0;        
03025           
03026           for(unsigned int j=0; j<SeedSegments.size(); ++j) {
03027             AlreadyAdded=false;
03028             
03029             for(unsigned int k=0; k<BestSeedSegments.size(); ++k) {
03030               if(BestSeedSegments[k]==SeedSegments[j]) {AlreadyAdded=true;}
03031             }
03032             if(AlreadyAdded) {continue;}
03033 
03034             if(SeedSegments[j]->GetEntries()>MostClusters) {
03035               MostClusters=SeedSegments[j]->GetEntries();
03036               LargestSeg=SeedSegments[j];
03037             }
03038           }
03039 
03040           if(LargestSeg) {BestSeedSegments.push_back(LargestSeg);}
03041         }
03043 
03044         
03045 
03046         // Now use the final seed segments to actually PROPAGATE back and forth, 
03047         // setting the TmpTrkFlags for the segments used in the propagation.
03049         TempBank2.clear();
03050 
03051         // Loop over the seed segments, which are stored in BestSeedSegments
03052         for(unsigned int i=0; i<BestSeedSegments.size(); ++i) {
03053           TrackSegmentCam* Seg0 = BestSeedSegments[i];
03054 
03055           MSG("AlgTrackCamList", Msg::kVerbose) << " Seed BegPlane" << Seg0->GetBegPlane() 
03056                                                 << " Endplane " << Seg0->GetEndPlane()
03057                                                 << " Begtpos " << Seg0->GetBegTPos() 
03058                                                 << " Endtpos " << Seg0->GetEndTPos()
03059                                                 << " UID " << Seg0->GetUID()
03060                                                 << " entries " << Seg0->GetEntries() << endl;
03061 
03062 
03063           // If we make a track, we will definitely include the seed segment, so set
03064           // it's TmpTrkFlag to 3.
03065           Seg0->SetTmpTrkFlag(3); 
03066           TrackFlag=false;
03067 
03068           TempBank1.clear(); TempBank1.push_back(Seg0); 
03069 
03070 
03071           // Track backwards (1)
03073           mTemp.clear(); BegTemp.clear();
03074 
03075           // Push the seed segment into BegTemp and begin while loop
03076           BegTemp.push_back(Seg0);
03077 
03078           while(BegTemp.size()>0) {
03079             Temp.clear();
03080 
03081             // Copy any segments into a temporary holder and loop over them,
03082             // gradually propagating back to lower plane numbers
03083             Temp=BegTemp; BegTemp.clear();
03084 
03085             for(unsigned int j=0; j<Temp.size(); ++j) {
03086               tmpFlag=false;
03087               TrackSegmentCam* Segtmp = Temp[j];
03088               
03089               // For each segment, loop over their beginning matches
03090               for(unsigned int k=0; k<Segtmp->GetNMatchSegBeg(); ++k){
03091                 TrackSegmentCam* Segb = Segtmp->GetMatchSegBeg(k);
03092                 
03093                 // If the matched segment is track-like and has zero TmpTrkFlag, 
03094                 // push back into BegTemp, so the while loop continues.
03095                 if(Segb->GetUID()>-1 && Segb->GetUID()<3) {
03096                   tmpFlag=true;
03097                   
03098                   if(Segb->GetTmpTrkFlag()<1) {
03099                     BegTemp.push_back(Segb);
03100                     Segb->SetTmpTrkFlag(1); // Mark the segments with TmpTrkFlag 1
03101                     TempBank1.push_back(Segb);
03102                   }
03103                 }
03104               }
03105 
03106               // Must have gone as far as we can go with this seed segment. Store segment at
03107               // lowest plane.
03108               if(tmpFlag==false) {mTemp.push_back(Segtmp);}
03109             }
03110           }
03111           
03112           // Note the first plane of the potential 2D track
03113           nplane=999;
03114           for(unsigned int j=0; j<mTemp.size(); ++j) {
03115             if(mTemp[j]->GetBegPlane()<nplane) {nplane=mTemp[j]->GetBegPlane();}
03116           }
03117 
03118 
03119           // For the longest possible paths, move back towards the seed segment,
03120           // setting TmpTrkFlags to 2 for the segments used.
03122 
03123           // Find the segments reached by the longest propagations
03124           BegTemp.clear();
03125           for(unsigned int j=0; j<mTemp.size(); ++j) {
03126             if(mTemp[j]->GetBegPlane()<nplane+5) {BegTemp.push_back(mTemp[j]);}
03127           }
03128 
03129 
03130           // From these segments, propagate forwards and set the flags.
03131           while(BegTemp.size()>0) {
03132             Temp.clear();
03133             
03134             for(unsigned int j=0; j<BegTemp.size(); ++j) {
03135               if(BegTemp[j]->GetTmpTrkFlag()<2) {
03136                 BegTemp[j]->SetTmpTrkFlag(2); // Mark the segments with TmpTrkFlag 2
03137                 Temp.push_back(BegTemp[j]);
03138               }
03139               
03140               if(BegTemp[j]->GetTrkFlag()<1) {
03141                 BegTemp[j]->SetTrkFlag(1);    // Also mark the segments with TrkFlag 1
03142                 TempBank2.push_back(BegTemp[j]);
03143                 TrackFlag=true;
03144               }
03145             }
03146             BegTemp.clear();
03147             
03148             
03149             for(unsigned int j=0; j<Temp.size(); ++j) {
03150               for(unsigned int k=0; k<Temp[j]->GetNMatchSegEnd(); ++k) {
03151                 if(Temp[j]->GetMatchSegEnd(k)->GetTmpTrkFlag()>0) {
03152                   BegTemp.push_back(Temp[j]->GetMatchSegEnd(k));
03153                 }
03154               }
03155             }
03156           }
03158           
03159 
03160 
03161           // Track forwards (1)
03163           pTemp.clear(); EndTemp.clear();
03164           
03165           
03166           // Push seed segment into EndTemp and begin while loop
03168           EndTemp.push_back(Seg0);
03169 
03170           while(EndTemp.size()>0) {
03171             Temp.clear();
03172          
03173             // Copy any segments into a temporary holder and loop over them,
03174             // gradually propagating to higher plane numbers
03175             Temp=EndTemp; EndTemp.clear();
03176 
03177             for(unsigned int j=0; j<Temp.size(); ++j) {
03178               tmpFlag=false;
03179               TrackSegmentCam* Segtmp = Temp[j];
03180               
03181               // For each segment, loop over their end matches
03182               for(unsigned int k=0; k<Segtmp->GetNMatchSegEnd(); ++k){
03183                 TrackSegmentCam* Sege = Segtmp->GetMatchSegEnd(k);
03184                 
03185                 // If the matched segment is track-like and has zero TmpTrkFlag, 
03186                 // push back into BegTemp, so the while loop continues.
03187                 if(Sege->GetUID()>-1 && Sege->GetUID()<3) {
03188                   tmpFlag=true;
03189                   
03190                   if(Sege->GetTmpTrkFlag()<1) {
03191                     EndTemp.push_back(Sege);
03192                     Sege->SetTmpTrkFlag(1); // Mark the segments with TmpTrkFlag 1
03193                     TempBank1.push_back(Sege);
03194                   }
03195                 }
03196               }
03197 
03198               // Must have gone as far as we can go with this seed segment. Store segment at
03199               // highest plane.
03200               if(tmpFlag==false) {pTemp.push_back(Segtmp);}
03201             }
03202           }
03203 
03204           // Note the end plane of potential 2D track     
03205           nplane=-999;
03206           for(unsigned int j=0; j<pTemp.size(); ++j) {
03207             if(pTemp[j]->GetEndPlane()>nplane) {nplane=pTemp[j]->GetEndPlane();}
03208           }
03209 
03210 
03211           // For the longest possible paths, move back towards the seed segment,
03212           // setting TmpTrkFlags to 2 for the segments used.
03214 
03215           // Find the segments reached by the longest propagations
03216           EndTemp.clear();
03217           for(unsigned int j=0; j<pTemp.size(); ++j) {
03218             if(pTemp[j]->GetEndPlane()>nplane-5) {EndTemp.push_back(pTemp[j]);}
03219           }
03220 
03221           
03222           // From these segments, propagate back and set the flags.
03223           while(EndTemp.size()>0) {
03224             Temp.clear();
03225             
03226             for(unsigned int j=0; j<EndTemp.size(); ++j) {
03227               if(EndTemp[j]->GetTmpTrkFlag()<2) {
03228                 EndTemp[j]->SetTmpTrkFlag(2); // Mark the segments with TmpTrkFlag 2
03229                 Temp.push_back(EndTemp[j]);
03230               }
03231               
03232               if(EndTemp[j]->GetTrkFlag()<1) {
03233                 EndTemp[j]->SetTrkFlag(1);    // Also mark the segments with TrkFlag 1
03234                 TempBank2.push_back(EndTemp[j]);
03235                 TrackFlag=true;
03236               }
03237             }
03238             EndTemp.clear();
03239             
03240             for(unsigned int j=0; j<Temp.size(); ++j) {
03241               for(unsigned int k=0; k<Temp[j]->GetNMatchSegBeg(); ++k) {
03242                 if(Temp[j]->GetMatchSegBeg(k)->GetTmpTrkFlag()>0) {
03243                   EndTemp.push_back(Temp[j]->GetMatchSegBeg(k));
03244                 }
03245               }
03246             }
03247           }
03248 
03249           // End setting TmpTrkFlags
03251 
03252 
03253 
03254 
03255           // From paths of segments with TmpTrkFlag==2, find the best set of
03256           // segments at the beginning and the best set of segments at the end
03258 
03259           // If propagation from the seed segment has given us a good potential 2D track
03260           MSG("AlgTrackCamList", Msg::kVerbose) << "TrackFlag " << TrackFlag << endl;
03261 
03262           if(TrackFlag==true) { 
03263 
03264             // Track backwards (2) and find best set of beginning segments
03266             TempBegBank.clear(); ObjVectorCounter=0;
03267             vector<TrackSegmentCam*> TmpBegSegBank;
03268 
03269             TmpBegSegBank.push_back(Seg0);
03270             TempBegBank.push_back(TmpBegSegBank);
03271             ObjCounter=1;
03272             
03273 
03274             // Store the series of beginning segments that were flagged above during the propagation.
03275             // Each entry in TempBegbank will be a vectors of segments (a series of beginning segments 
03276             // flagged above).
03277 
03278             // TempBegBank will contain vectors of segments representing every possible path back through
03279             // the segments marked with TmpTrkFlag 2 (possible beginning sections for the 2D track)
03281             while(ObjCounter>0 && ObjVectorCounter<fObjVectorMax) {
03282               ObjCounter=0;
03283               npts=TempBegBank.size();
03284 
03285               for(int j=0; j<npts; ++j) {
03286                 TrackSegmentCam* Segtmp = TempBegBank[j].back();
03287                 mTemp.clear();
03288                 
03289                 for(unsigned int k=0; k<Segtmp->GetNMatchSegBeg(); ++k) {
03290                   TrackSegmentCam* Segbeg = Segtmp->GetMatchSegBeg(k);
03291                   if(Segbeg->GetTmpTrkFlag()>1) {mTemp.push_back(Segbeg);};
03292                 }
03293 
03294                 if(mTemp.size()>0) {
03295                   for(unsigned int k=1; k<mTemp.size(); ++k) {
03296                     if(ObjVectorCounter<fObjVectorMax) { 
03297                       vector<TrackSegmentCam*> NewObj = TempBegBank[j];
03298                       
03299                       NewObj.push_back(mTemp[k]);
03300                       TempBegBank.push_back(NewObj);
03301 
03302                       ObjVectorCounter++;
03303                     }
03304                     else {MSG("AlgTrackCamList", Msg::kWarning) << " *** CRAZY MEMORY GUZZLING *** " << endl; break;}
03305                   }
03306 
03307                   if(ObjVectorCounter>=fObjVectorMax) {break;}
03308 
03309                   TempBegBank[j].push_back(mTemp[0]); ObjCounter++;
03310                 }
03311               }
03312             } // End while ObjCounter>0
03314 
03315             // Find beginning section with best score and store it
03316             MSG("AlgTrackCamList", Msg::kVerbose) << "beginning score " << endl;
03317 
03318             id=-1; TopScore=-1.;
03319             for(unsigned int j=0; j<TempBegBank.size(); ++j) {
03320               Score = Seg0->GetScore(&TempBegBank[j],0);
03321               if(Score>TopScore) {TopScore=Score; id=j;}
03322             }
03323 
03324             if(id!=-1) {BegBank.push_back(TempBegBank[id]);}
03325             TempBegBank.clear();
03327 
03328 
03329 
03330 
03331             // Track forwards (2) and find best set of end segments
03333             TempEndBank.clear(); ObjVectorCounter=0;
03334             vector<TrackSegmentCam*> TmpEndSegBank;
03335 
03336             TmpEndSegBank.push_back(Seg0);
03337             TempEndBank.push_back(TmpEndSegBank);
03338             ObjCounter=1;
03339             
03340 
03341             // Store the series of end segments that were flagged above during the propagation.
03342             // Each entry in TempEndbank will be a vectors of segments (a series of end segments 
03343             // flagged above).
03344 
03345             // TempBegBank will contain vectors of segments representing every possible path forward through
03346             // the segments marked with TmpTrkFlag 2 (possible end sections for the 2D track)
03348             while(ObjCounter>0 && ObjVectorCounter<fObjVectorMax) {
03349               ObjCounter=0;
03350               npts=TempEndBank.size();
03351               
03352               for(int j=0; j<npts; ++j) {
03353                 TrackSegmentCam* Segtmp = TempEndBank[j].back();
03354 
03355                 pTemp.clear();
03356                 
03357                 for(unsigned int k=0; k<Segtmp->GetNMatchSegEnd(); ++k) {
03358                   TrackSegmentCam* Segend = Segtmp->GetMatchSegEnd(k);
03359                   if(Segend->GetTmpTrkFlag()>1) {pTemp.push_back(Segend);};
03360                 }
03361 
03362                 if(pTemp.size()>0) {
03363                   for(unsigned int k=1; k<pTemp.size(); ++k) {
03364                     if(ObjVectorCounter<fObjVectorMax) {
03365                       vector<TrackSegmentCam*> NewObj = TempEndBank[j];
03366                       
03367                       NewObj.push_back(pTemp[k]);
03368                       TempEndBank.push_back(NewObj);
03369 
03370                       ObjVectorCounter++;
03371                     }
03372                     else {MSG("AlgTrackCamList", Msg::kWarning) << " *** CRAZY MEMORY GUZZLING *** " << endl; break;}
03373                   }
03374 
03375                   if(ObjVectorCounter>=fObjVectorMax) {break;}
03376 
03377                   TempEndBank[j].push_back(pTemp[0]); ObjCounter++;
03378                 }
03379               }
03380             } // End while ObjCounter>0
03382 
03383             // Find end section with best score and store it
03384             MSG("AlgTrackCamList", Msg::kVerbose) << "end score " << endl;
03385 
03386             id=-1; TopScore=-1.;
03387             for(unsigned int j=0; j<TempEndBank.size(); ++j) {
03388               Score = Seg0->GetScore(0,&TempEndBank[j]);
03389               if(Score>TopScore) {TopScore=Score; id=j;}
03390             }
03391 
03392             if(id!=-1) {EndBank.push_back(TempEndBank[id]);}
03393             TempEndBank.clear(); // maybe clear individual vectors of segments
03395           }
03397 
03398           
03399           // Clear up
03400           for(unsigned int j=0; j<TempBank1.size(); ++j) {TempBank1[j]->SetTmpTrkFlag(0);}
03401 
03402         } // End loop over BestSeedSegments
03404         for(unsigned int j=0; j<TempBank2.size(); ++j) {TempBank2[j]->SetTrkFlag(0);}
03405 
03406 
03407 
03408 
03409 
03410         // Using the best beginning sections and the best end sections from all 
03411         // seed segments, find the BEST COMPLETE 2D track
03413         if(BegBank.size()>0 && EndBank.size()>0) {
03414 
03415           // Calculate a score, including the seed segment and best beginning/end paths,
03416           // to find the 'best' complete track
03417           MSG("AlgTrackCamList", Msg::kVerbose) << "overall score " << endl;
03418 
03419           id=-1; TopScore=-1.;
03420           for(unsigned int i=0; i<BegBank.size(); ++i) {
03421             Score = BegBank[i][0]->GetScore(&BegBank[i],&EndBank[i]);
03422             if(Score>TopScore) {TopScore=Score; id=i;}
03423           }
03424 
03425           // If we have found a best track, add together the segments to form
03426           // a segment that is the 2D track, giving this UID 3. Mark segments  
03427           // used by setting their UID to -1.
03428           if(1+id>0) {
03429             // Get the seed segment
03430             TrackSegmentCam* Seg0 = BegBank[id][0];
03431 
03432             // Set segment UIDs to -1 if they are part of the track, so they
03433             // can't be used in another track.
03434             for(unsigned int i=1; i<BegBank[id].size(); ++i) {
03435               TrackSegmentCam* Segbeg = BegBank[id][i];
03436               Seg0->AddSegment(Segbeg);
03437               Segbeg->SetUID(-1);
03438             }
03439 
03440             for(unsigned int i=1; i<EndBank[id].size(); ++i) {
03441               TrackSegmentCam* Segend = EndBank[id][i];
03442               Seg0->AddSegment(Segend);
03443               Segend->SetUID(-1);
03444             }
03445             
03446             // Mark long segment as a 2D track by setting its UID to 3
03447             Seg0->SetUID(3);
03448 
03449             // Mark all the clusters in 2D track with TrkFlag 3
03450             for(unsigned int i=0; i<Seg0->GetEntries(); ++i) {
03451               Seg0->GetCluster(i)->SetTrkFlag(3);
03452             }
03453 
03454             // Store for use in Join2D tracks method
03455             PossibleJoins[View].push_back(Seg0);
03456 
03457             // Can now look for more 2D tracks
03458             Cont=true; ntrks++;
03459           }
03460         }
03462 
03463         // Clean up
03464         BegBank.clear(); EndBank.clear();
03465 
03466         
03467       } // End while Cont==true loop
03468     } // End loop over planeviews
03469   } // End check to see if there is anything in ViewSegBank
03472   
03473 
03474   
03475   // Print out the list of 2D tracks
03477   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF 2D TRACKS *** " << endl;
03478   for(int View=0; View<2; ++View) {
03479 
03480     MSG("AlgTrackCamList", Msg::kVerbose) << "View: " << View << endl;
03481     for(unsigned i=0; i<ViewSegBank[View].size(); ++i) {
03482       TrackSegmentCam* Seg = ViewSegBank[View][i];
03483       if(Seg->GetUID()>2) {
03484         MSG("AlgTrackCamList", Msg::kVerbose) << "  i=" << i
03485                                               << "  begpln: "  << Seg->GetBegPlane()
03486                                               << "  begtpos: " << Seg->GetBegTPos()
03487                                               << "  endpln: "  << Seg->GetEndPlane()
03488                                               << "  endtpos: " << Seg->GetEndTPos() << endl;
03489       }
03490     }
03491   }
03493 
03494   return;
03495 }

void AlgTrackCamList::Form3DTracks ( CandSliceHandle slice  ) 

Definition at line 4191 of file AlgTrackCamList.cxx.

References ClusterCam::AddHit(), TrackCam::AddHit(), CambridgeAnalysis, ClusterList, MuELoss::e, ExtendTrack(), FillGapsInTrack(), FinalTrackBank, TrackCam::GetBegPlane(), TrackSegmentCam::GetBegPlane(), ClusterCam::GetBegTPos(), ClusterCam::GetCharge(), HitCam::GetCharge(), TrackSegmentCam::GetCluster(), TrackSegmentCam::GetEndPlane(), TrackCam::GetEndPlane(), ClusterCam::GetEndTPos(), TrackSegmentCam::GetEntries(), ClusterCam::GetEntries(), TrackCam::GetEntries(), ClusterCam::GetHit(), TrackCam::GetHit(), ClusterCam::GetNDFlag(), ClusterCam::GetPlane(), HitCam::GetPlane(), ClusterCam::GetPlaneView(), ClusterCam::GetShwPlnFlag(), HitCam::GetTPos(), ClusterCam::GetTrkPlnFlag(), HitCam::GetZPos(), ClusterCam::GetZPos(), HitsInTracks, id, ClusterCam::IsTrkAssoc(), Msg::kVerbose, LookForHitsAcrossGap(), Munits::m, MatchUV(), ModuleType, MSG, TrackCam::SetPartner(), HitCam::SetTrkFlag(), size, StripWidth, and TempTrack.

Referenced by RunTheFinder().

04192 {
04193   // The 3D track segments are now stored in the TempTracks container, with
04194   // U and V segments stored separately. In this method we look at the consituent
04195   // clusters and hits in these segments, performing fits to choose the final
04196   // strips for the track.
04197 
04198 
04199   MSG("AlgTrackCamList", Msg::kVerbose) << " *** Formation of 3D tracks *** " << endl;
04200 
04201   int ktemp, View;
04202   double swx, swy, swx2, swxy, sw, ProjectedTPos, TPosWindow(0.);
04203   double w, x, y, m, c;
04204   int UStrips, VStrips, TrackStrips, TrackPlanes;
04205 
04206   int begPlaneView[2], endPlaneView[2];
04207   int begplane1, begplane2, endplane1, endplane2, maxplane, Plane;
04208   double dScore, Score; int id;
04209 
04210   int ShwOrTrkList[500];
04211   int ViewList[500];
04212   vector<HitCam*> HitList[500];
04213   vector<TrackCam*> FinalTrackTempHolder[2];
04214   for(int k=0; k<500; ++k) {ShwOrTrkList[k]=-1; ViewList[k]=-1;}
04215 
04216 
04217 
04218   // Main loop is over the pairs of track segments, stored above
04220   for(unsigned int i=0; i<TempTrack[0].size(); ++i) {
04221     MSG("AlgTrackCamList", Msg::kVerbose) << " making track " << i << endl;
04222     
04223     
04224     // Find the beginning and end planes for each of this pair of segments.
04225     // Order them so that we have:
04226     // begplane2, begplane 1, endplane 1, endplane 2.
04228     begplane2=-1; endplane2=-1; begplane1=-1; endplane1=-1;
04229     
04230     for(View=0; View<2; ++View) {
04231       TrackSegmentCam* seg = TempTrack[View][i];
04232       if(begplane2<0 || seg->GetBegPlane()<begplane2) { begplane1=begplane2; begplane2=seg->GetBegPlane(); }
04233       else { begplane1=seg->GetBegPlane(); }
04234       
04235       if(endplane2<0 || seg->GetEndPlane()>endplane2) { endplane1=endplane2; endplane2=seg->GetEndPlane(); }
04236       else { endplane1=seg->GetEndPlane(); }
04237 
04238       begPlaneView[View]=-1; endPlaneView[View]=-1;
04239     }
04240 
04241     maxplane=1+endplane2;
04242     
04243     MSG("AlgTrackCamList", Msg::kVerbose) << begplane2 << " (" << begplane1 << ") -> (" << endplane1 << ") " << endplane2 << endl;
04245  
04246 
04247     // Unpack the clusters from the segments.
04248     // For each track, should only have one cluster on a plane.
04250     for(View=0; View<2; ++View) {
04251       TrackSegmentCam* Seg = TempTrack[View][i];
04252       
04253       for(unsigned int j=0; j<Seg->GetEntries(); ++j) {
04254         ClusterCam* Clust = Seg->GetCluster(j);
04255 
04256         Plane=Clust->GetPlane();
04257         if(Plane>=begplane2 && Plane<maxplane) {ClusterList[Plane].push_back(Clust);}
04258       }
04259     }
04261     
04262     
04263     
04264     // Trim the segments. Check we haven't tracked much further
04265     // in one view than the other, and adjust as necessary.
04266 
04267     // Maximum discrepancy allowed at each end is 3 planes. As only
04268     // clusters on planes greater than begplane1 are considered.
04270     int PrevPln, NextPln;
04271     bool LeavesAtEdge=false;
04272 
04273     // Check hasn't just left edge of detector in one view.
04274     if(ModuleType==1 && (ClusterList[begplane1][0]->GetEndTPos()>3.9 
04275                          || ClusterList[begplane1][0]->GetBegTPos()<-3.9)) {begplane1=begplane2+1; LeavesAtEdge=true;}
04276     if(ModuleType==1 && (ClusterList[endplane1][0]->GetEndTPos()>3.9 
04277                          || ClusterList[endplane1][0]->GetBegTPos()<-3.9)) {endplane1=endplane2-1; LeavesAtEdge=true;}
04278 
04279     if(!LeavesAtEdge) {
04280       // For the beginning
04282       // Jump back 4 planes before begplane1 (unless this takes us past begplane2)
04283       begplane1=begplane1-4; PrevPln=-1;
04284       Plane=begplane1; if(Plane<begplane2) {Plane=begplane2;}
04285     
04286       // Working forwards from this plane, find the first plane that has
04287       // a cluster.
04288       while(Plane<maxplane) {
04289         if(ClusterList[Plane].size()>0) {PrevPln=Plane; break;}
04290         else {Plane++;}
04291       }
04292     
04293       // Look at the next 2 planes in the same view as this cluster
04294       // and see if they have clusters too. Try and veto geometries such as
04295       //
04296       //                bp1
04297       // u       x   x   o   o
04298       // v         o   x   *   o
04299       //          bp2      * == shw-like, or gap
04300       //
04301       // In this case, begplane1 is set to begplane 2, so the lone hit at 
04302       // begplane 2 is later vetoed.
04303       if(PrevPln>-1 && PrevPln+4<maxplane && ClusterList[PrevPln+2].size()==0) {
04304         if(ClusterList[PrevPln+4].size()>0) {
04305           ClusterCam* Clust = ClusterList[PrevPln+4][0];
04306         
04307           if(Clust->GetTrkPlnFlag()>0 && Clust->GetShwPlnFlag()==0) {
04308             PrevPln=begplane1;
04309           }
04310         }
04311         begplane1=PrevPln;
04312       }
04313     
04314     
04315       // Similarly, for the end 
04317       // Jump forward 4 planes after endplane1 (unless this takes us past endplane2)
04318       endplane1=endplane1+4; NextPln=-1;
04319       Plane=endplane1; if(Plane>maxplane-1) {Plane=maxplane-1;}
04320     
04321       // Working backwards from this plane, find the first plane that has
04322       // a cluster.
04323       while(Plane>begplane2-1) {
04324         if(ClusterList[Plane].size()>0) {NextPln=Plane; break;}
04325         else {Plane--;}
04326       }
04327     
04328       // Look at the next 2 planes in the same view as this cluster
04329       // and see if they have clusters too, as above.
04330       if(NextPln>-1 && NextPln-4>(begplane2-1) && ClusterList[NextPln-2].size()==0) {
04331         if(ClusterList[NextPln-4].size()>0) {
04332           ClusterCam* Clust = ClusterList[NextPln-4][0];
04333         
04334           if(Clust->GetTrkPlnFlag()>0 && Clust->GetShwPlnFlag()==0) {
04335             NextPln=endplane1;
04336           }
04337         }
04338         endplane1=NextPln;
04339       }
04341     }
04343 
04344     
04345 
04346     // Count the strips and record planeviews, number of occupied planes, etc.
04348     UStrips=0; VStrips=0; TrackStrips=0;
04349 
04350     // Look at the region over which the track is now defined
04351     for(int k=begplane2; k<maxplane; ++k) {  
04352 
04353       if( //  (k>=begplane1 && k<=endplane1) &&
04354           ClusterList[k].size()>0 ) {
04355 
04356         ClusterCam* Clust = ClusterList[k][0];
04357         
04358         View=Clust->GetPlaneView();
04359         ViewList[k]=View;
04360         ShwOrTrkList[k]=1;
04361         
04362         // Record first and last planes in each view.
04363         if(begPlaneView[View]<0 || k<begPlaneView[View]) {begPlaneView[View]=k;}
04364         if(endPlaneView[View]<0 || k>endPlaneView[View]) {endPlaneView[View]=k;}
04365         
04366         if(View==0) {UStrips++;}
04367         if(View==1) {VStrips++;}
04368         
04369         if(Clust->GetTrkPlnFlag()>0) {TrackStrips++;}
04370       }
04371     }
04372     
04373     
04374     MSG("AlgTrackCamList", Msg::kVerbose) 
04375       << " U: (" << begPlaneView[0] << "," << endPlaneView[0] << ") " 
04376       << " V: (" << begPlaneView[1] << "," << endPlaneView[1] << ") " 
04377       << " UStrips " << UStrips << " VStrips " << VStrips << " TrackStrips " << TrackStrips << endl;
04379     
04380     
04381 
04382     
04383     // If we have at least 3 strips in each view and at least 2 clusters
04384     // on track-like planes, proceed with the fits.
04386     if(UStrips>2 && VStrips>2 && (CambridgeAnalysis==false || TrackStrips>2) ) {
04387 
04388 
04389       // FLAG PLANES which we think are track-like
04391       for(View=0; View<2; ++View) {
04392         for(Plane=begplane2; Plane<maxplane; ++Plane) {
04393 
04394           if(ViewList[Plane]==View && ShwOrTrkList[Plane]>0) {
04395             ClusterCam* Clust0 = ClusterList[Plane][0];
04396 
04397             // Always set flag to 3 for the +/- 5 plane ND track segments
04398             if(Clust0->GetNDFlag()==2) {ShwOrTrkList[Plane]=3;}
04399 
04400 
04401             // Find next planes, in the same view, where we have a cluster, both forwards
04402             // and backwards.
04404             ktemp=0; PrevPln=-1;
04405             while(Plane-ktemp>begplane2) {
04406               ktemp++;
04407               if(ViewList[Plane-ktemp]==View) {PrevPln=ktemp; break;}
04408             }
04409 
04410             ktemp=0; NextPln=-1;
04411             while(Plane+ktemp<maxplane-1) {
04412               ktemp++;
04413               if(ViewList[Plane+ktemp]==View) {NextPln=ktemp; break;}
04414             }
04416 
04417 
04418             // If there are planes forwards and backwards in the same view...
04420             if(PrevPln>0 && NextPln>0) {
04421               ClusterCam* NextClust = ClusterList[Plane+NextPln][0];
04422               ClusterCam* PrevClust = ClusterList[Plane-PrevPln][0];
04423               
04424               // If either both clusters are within range for this view, or they are the
04425               // next plane along and the current plane is a track plane.
04426               if( (PrevClust->GetPlane()>begPlaneView[View]
04427                    || (Clust0->GetPlane()-PrevClust->GetPlane()<3 && Clust0->GetTrkPlnFlag()>0))
04428                   && (NextClust->GetPlane()<endPlaneView[View]
04429                       || (NextClust->GetPlane()-Clust0->GetPlane()<3 && Clust0->GetTrkPlnFlag()>0) )) {
04430 
04431 
04432                 // If the clusters have track-like association, identify plane as
04433                 // being track like, by setting entry in ShwOrTrkList to 3.
04434                 if(Clust0->IsTrkAssoc(PrevClust,NextClust)>1) {
04435                   ShwOrTrkList[Plane]=3;
04436 
04437                   // If there is a good overlap of strip numbers between the three clusters
04438                   if( (NextClust->GetBegTPos()-Clust0->GetEndTPos()>-1.1*StripWidth 
04439                        && Clust0->GetBegTPos()-PrevClust->GetEndTPos()>-1.1*StripWidth)
04440                       || (NextClust->GetEndTPos()-Clust0->GetBegTPos()<1.1*StripWidth
04441                           && Clust0->GetEndTPos()-PrevClust->GetBegTPos()<1.1*StripWidth) ) {
04442 
04443                     // If forwards and backwards planes are adjacent in view,
04444                     // identify these planes as track-like too
04445                     if( Clust0->GetPlane()-PrevClust->GetPlane()<3 || Clust0->GetNDFlag()==2) {ShwOrTrkList[Plane-PrevPln]=3;}
04446                     if( NextClust->GetPlane()-Clust0->GetPlane()<3 || Clust0->GetNDFlag()==2) {ShwOrTrkList[Plane+NextPln]=3;}
04447                   }
04448                 }
04449 
04450                 // Otherwise, if cluster is small, can identify its plane as track-like.
04451                 // Don't set anything for forwards and backwards planes.
04452                 else {
04453                   if(Clust0->GetEndTPos()-Clust0->GetBegTPos()<2.1*StripWidth 
04454                      && Clust0->GetTrkPlnFlag()>0) {ShwOrTrkList[Plane]=3;}
04455                 }
04456               }
04457             }
04459 
04460           }
04461         } // End loop over planes
04462       } // End loop over views
04464 
04465 
04467       MSG("AlgTrackCamList", Msg::kVerbose) << " *** SORT TRACK HITS *** " << endl;
04468       
04469       for(Plane=begplane2;Plane<maxplane; ++Plane) {
04470         if(ViewList[Plane]>-1) {
04471           ClusterCam* Clust = ClusterList[Plane][0];
04472           MSG("AlgTrackCamList", Msg::kVerbose) << Clust->GetPlane() 
04473                                                 << " " << Clust->GetBegTPos() << " " << Clust->GetEndTPos()
04474                                                 << " " << ShwOrTrkList[Plane] << endl;
04475         }
04476       }
04478 
04479 
04480 
04481       // Perform GLOBAL LINEAR FIT
04483       int ShowerStart, ShowerEnd, FitStart, FitEnd, PrevTrkPln, NextTrkPln;
04484 
04485       for(View=0; View<2; ++View) {
04486         for(Plane=begplane2; Plane<maxplane; ++Plane) {
04487 
04488 
04489           // Currently ShwOrTrkList flag==1 means the plane was part of the 3D track.
04490           // ShwOrTrkList flag==3 means that we also flagged the plane as being track-like, above.
04491           //
04492           // Here, we look at the planes we haven't flagged as track-like, find track-like planes
04493           // on either side and carry out a linear fit through this 'shower-like' region.
04494           // 
04495           // Using the fit, clusters in this shower-like region are identified, then hits.
04497           if(ViewList[Plane]==View && ShwOrTrkList[Plane]<2) {
04498             MSG("AlgTrackCamList", Msg::kVerbose) << " *** GLOBAL FIT ***   Plane=" << Plane << endl;
04499 
04500             // Look for next planes in the view that we have flagged as track-like, 
04501             // forwards and backwards.
04503             ktemp=0; NextTrkPln=-1;
04504             while(Plane+ktemp<endplane1 && Plane+ktemp<maxplane-1) {
04505               ktemp++;
04506               if(ViewList[Plane+ktemp]==View && ShwOrTrkList[Plane+ktemp]>2) {NextTrkPln=ktemp; break;}
04507             }
04508             
04509             ktemp=0; PrevTrkPln=-1;
04510             while(Plane-ktemp>begplane1 && Plane-ktemp>begplane2) {
04511               ktemp++;
04512               if(ViewList[Plane-ktemp]==View && ShwOrTrkList[Plane-ktemp]>2) {PrevTrkPln=ktemp; break;}
04513             }
04515 
04516             
04517             // Find the track-like planes which are just outside of the shower-like region ('ShowerStart' 
04518             // and 'ShowerEnd' planes away), as well as the planes which are up to 7 planes either side 
04519             // of the shower-like region ('FitStart' and 'FitEnd' planes away). These are the 
04520             // track used in the linear fit through the shower.
04522             if(PrevTrkPln>0) {FitStart=-PrevTrkPln-7; ShowerStart=-PrevTrkPln;} 
04523             else {FitStart=begplane1-Plane+1; ShowerStart=FitStart;}
04524 
04525             if(Plane+FitStart<begplane2) {FitStart=-Plane+begplane2;}
04526             if(Plane+ShowerStart<begplane2) {ShowerStart=-Plane+begplane2;}
04527 
04528             if(NextTrkPln>0) {FitEnd=1+NextTrkPln+7; ShowerEnd=1+NextTrkPln;}
04529             else {FitEnd=endplane1-Plane; ShowerEnd=FitEnd;}
04530 
04531             if(Plane+FitEnd>maxplane) {FitEnd=maxplane-Plane;}
04532             if(Plane+ShowerEnd>maxplane) {ShowerEnd=maxplane-Plane;}
04534 
04535 
04536             ClusterCam* ClustSeed = ClusterList[Plane][0];
04537 
04538             m=0.; c=0.5*(ClustSeed->GetBegTPos()+ClustSeed->GetEndTPos());
04539             swx=0.; swy=0.; swx2=0.; swxy=0.; sw=0.;
04540 
04541 
04542             // Carry out the linear fit through the shower, obtaining gradient and intercept.
04544             for(int k=FitStart; k<FitEnd; ++k) {
04545               if(ViewList[Plane+k]==View) {
04546                 ClusterCam* ClustTemp = ClusterList[Plane+k][0];
04547                 
04548                 for(unsigned int l=0; l<ClustTemp->GetEntries(); ++l) {
04549                   HitCam* HitTemp = ClustTemp->GetHit(l);
04550 
04551                   x=HitTemp->GetZPos();
04552                   y=HitTemp->GetTPos();
04553                   w=HitTemp->GetCharge()/ClustTemp->GetCharge();
04554                   
04555                   if(ShwOrTrkList[Plane+k]>2) {w=5.*w;}
04556                   else {w=0.5*w;}
04557                   
04558                   swx+=w*x; swy+=w*y; swx2+=w*x*x; swxy+=w*x*y; sw+=w;
04559                 }
04560               }
04561             }
04562             double denom=sw*swx2-swx*swx;
04563             if(sw>0. && fabs(denom)>1.e-10) {
04564               m=(sw*swxy-swx*swy)/denom;
04565               c=(swy*swx2-swx*swxy)/denom;
04566             }
04568 
04569 
04570             // Now loop over the planes in the shower-like region
04572             for(int k=ShowerStart; k<ShowerEnd; ++k) {
04573               if(ViewList[Plane+k]==View && ShwOrTrkList[Plane+k]<2) {
04574                 // ***First cluster in ClusterList is that from segment in 3D track***
04575                 ClusterCam* ClustTemp = ClusterList[Plane+k][0];
04576                 
04577                 // Projected strip number in shower
04578                 ProjectedTPos = m*ClustTemp->GetZPos()+c;
04579 
04580                 // Calculate window around projected strip number
04581                 if(m>=0.) {TPosWindow = (0.6*StripWidth)+0.02*m;}
04582                 if(m<=0.) {TPosWindow = (0.6*StripWidth)-0.02*m;}
04583 
04584                 MSG("AlgTrackCamList", Msg::kVerbose) << " Plane=" << Plane+k << " ProjectedTPos=" << ProjectedTPos 
04585                                                       << " TPosWindow=" << TPosWindow << " (" << ClustTemp->GetBegTPos() 
04586                                                       << " " << ClustTemp->GetEndTPos() << ") " << endl;
04587 
04588 
04589                 // If  1. we have found track-like planes on either side of the shower-like region,
04590                 // or  2. we haven't found any of these track-like planes (entire track in this view must be shower),
04591                 // or  3. current cluster overlaps quite nicely with projected strip number from the fit,
04592                 //
04593                 // then proceed with trying to find the best hits from within the cluster.
04595                 if( (PrevTrkPln>0 && NextTrkPln>0) || (PrevTrkPln<0 && NextTrkPln<0)
04596                     || (ClustTemp->GetEndTPos()>ProjectedTPos-(1.6*StripWidth) && ClustTemp->GetBegTPos()<ProjectedTPos+(1.6*StripWidth)) ) {
04597 
04598                   id=-1; dScore=0.; Score=999.9;
04599 
04600                   // Loop over hits in cluster and get a 'score' for each. Score is based
04601                   // on distance from projected tpos, and we want to minimise this.
04602                   for(unsigned int l=0; l<ClustTemp->GetEntries(); ++l) {
04603                     HitCam* HitTemp = ClustTemp->GetHit(l);
04604                     
04605                     dScore=HitTemp->GetTPos()-ProjectedTPos;
04606                     if(dScore<0) {dScore=-dScore;}
04607 
04608                     if(dScore<Score) {id=l; Score=dScore;}
04609                   }
04611 
04612         
04613                   // If we have found a 'best hit'
04615                   if((1+id)>0) {
04616                     HitCam* Hit = ClustTemp->GetHit(id);
04617                     
04618                     // Create new cluster and store it. 
04619                     // ***Last cluster in ClusterList will be that found by the fit.***
04620                     ClusterCam* ClustNew = new ClusterCam(Hit);
04621                     ClusterList[Plane+k].push_back(ClustNew);
04622                     
04623                     ShwOrTrkList[Plane+k]=2;
04624                     
04625                     // If any other hits in the cluster found by the fit are also good, 
04626                     // add these to the new cluster too.
04627                     for(unsigned int l=0; l<ClustTemp->GetEntries(); ++l) {
04628                       HitCam* HitTemp = ClustTemp->GetHit(l);
04629                         
04630                       if( (HitTemp->GetTPos()-Hit->GetTPos()>0
04631                            && HitTemp->GetTPos()-Hit->GetTPos()<TPosWindow)
04632                           || (HitTemp->GetTPos()-Hit->GetTPos()<0
04633                               && HitTemp->GetTPos()-Hit->GetTPos()>-TPosWindow) ) {
04634                           
04635                         ClustNew->AddHit(HitTemp);
04636                       }
04637                     }
04638 
04639                     for(unsigned int l=0; l<ClustNew->GetEntries(); ++l) {
04640                       MSG("AlgTrackCamList", Msg::kVerbose) << "   ...global add hit=" << ClustNew->GetHit(l)->GetTPos() << " (" << Plane+k << ") " << endl;
04641                     }
04642                   }
04644 
04645                 }
04646               }
04647             }
04649           }
04650         }  // End loop over planes
04651       } // End loop over views
04653 
04654 
04655 
04656       // Now, Perform LOCAL LINEAR FITS
04658       for(View=0; View<2; ++View) {
04659         for(Plane=begplane2; Plane<maxplane; ++Plane) {
04660 
04661           // If plane has been declared track-like
04662           if(ViewList[Plane]==View && ShwOrTrkList[Plane]>1) {
04663             MSG("AlgTrackCamList", Msg::kVerbose) << " *** FIT (2) ***   Plane=" << Plane << endl;
04664 
04665             // Get last cluster from ClusterList. This will have been flagged immediately as track-like,
04666             // or have been stored as a result of a linear fit through a shower-like region.
04667             // Using the cluster, calculate a projected strip number and window for use in the local fit.
04668             ClusterCam* clr = ClusterList[Plane].back();
04669             ProjectedTPos=0.5*(clr->GetBegTPos()+clr->GetEndTPos()); 
04670             TPosWindow=0.6*StripWidth;
04671 
04672             // Loop over the hits in this cluster
04674             const unsigned int nhits = clr->GetEntries();
04675             if(nhits>1){
04676               for(unsigned int k=0; k<nhits; ++k) {clr->GetHit(k)->SetTrkFlag(1);}
04677               
04678               swx=0.0; swy=0.0; swx2=0.0; swxy=0.0; sw=0.0;
04679 
04680               // For these, loop over the hits in nearby clusters (within 4 planes of the current plane)
04681               // and carry out the linear fit
04683               for(int k=-4; k<5; ++k){
04684                 if( (Plane+k)>=begplane2 && (Plane+k)<(maxplane) && ViewList[Plane+k]==View && ShwOrTrkList[Plane+k]>1 ) {
04685                   ClusterCam* clrtmp = ClusterList[Plane+k].back(); // Again get last entry 
04686 
04687                   for(unsigned int l=0; l<clrtmp->GetEntries(); ++l){
04688                     HitCam* hittmp = clrtmp->GetHit(l);
04689                     x=hittmp->GetZPos(); 
04690                     y=hittmp->GetTPos(); 
04691                     w=hittmp->GetCharge()/clrtmp->GetCharge();
04692 
04693                     swx+=w*x; swy+=w*y; swx2+=w*x*x; swxy+=w*x*y; sw+=w;
04694                   }
04695                 }
04696               }
04697               double denom=sw*swx2-swx*swx;
04698               if(sw>1.0 && fabs(denom)>1.e-10){
04699                 m=(sw*swxy-swx*swy)/denom;
04700                 c=(swy*swx2-swx*swxy)/denom;
04701                 ProjectedTPos=m*clr->GetZPos()+c; 
04702                 if(m>=0.0) TPosWindow = (0.6*StripWidth)+0.02*m; 
04703                 if(m<=0) TPosWindow = (0.6*StripWidth)-0.02*m; 
04704               }
04706             }
04708             
04709             MSG("AlgTrackCamList", Msg::kVerbose) << "  local plane=" << Plane << " ProjectedTPos=" << ProjectedTPos 
04710                                                   << " TPosWindow=" << TPosWindow << endl;
04711 
04712 
04713             // Using the results of the fit, obtain a score for each hit. Again, this is
04714             // based on distance from the projected strip number. We want to pick the hit which
04715             // minimises this distance.
04717             id=-1; dScore=0.; Score=999.;
04718             
04719             for(unsigned int k=0; k<nhits; ++k){
04720               HitCam* hit = clr->GetHit(k);  
04721               dScore=hit->GetTPos()-ProjectedTPos; if(dScore<0) dScore=-dScore;
04722               
04723               if(dScore<Score){id=k; Score=dScore;}
04724             }
04725             
04726             // If we have found a good hit, store it and flag the plane as part of the 
04727             // final track
04728             if(1+id>0){  
04729               HitCam* hit = clr->GetHit(id); 
04730               HitList[Plane].push_back(hit); ShwOrTrkList[Plane]=4;
04731 
04732               // Include any other hits which also match projection within the window
04733               for(unsigned int k=0; k<nhits; ++k){
04734                 HitCam* hittmp = clr->GetHit(k);
04735                 if( ( hittmp->GetTPos()-hit->GetTPos()>0 && hittmp->GetTPos()-hit->GetTPos()<TPosWindow)
04736                     || ( hittmp->GetTPos()-hit->GetTPos()<0 && hittmp->GetTPos()-hit->GetTPos()>-TPosWindow) ){
04737                   HitList[Plane].push_back(hittmp);
04738                 }
04739               }
04740             }
04741 
04742             for(unsigned int l=0; l<HitList[Plane].size(); ++l){
04743               HitCam* hittmp = HitList[Plane][l];
04744               MSG("AlgTrackCamList", Msg::kVerbose) << "   ...local add hit=" << hittmp->GetTPos() << " (" << Plane << ") " << endl;
04745             }
04747 
04748           } // End demand that plane has been flagged as track-like
04749         } // End loop over planes
04750       } // End loop over views
04752       
04753 
04754 
04755       // *** FINALLY, FORM THE TRACKS ***
04757       // Count the number of planes flagged as being part of final track
04758       TrackPlanes=0;
04759       
04760       for(Plane=begplane2; Plane<maxplane; ++Plane) {
04761         if(ShwOrTrkList[Plane]>3) {TrackPlanes++;}
04762       }
04763 
04764       MSG("AlgTrackCamList", Msg::kVerbose) << " TrackPlanes " << TrackPlanes << endl;
04765 
04766       if(CambridgeAnalysis==false || TrackPlanes>5) {
04767         for(View=0; View<2; ++View) {
04768           
04769           TrackCam* Trk = new TrackCam(slice);
04770           FinalTrackTempHolder[View].push_back(Trk);
04771           
04772           for(Plane=begplane2; Plane<maxplane; ++Plane) {
04773             if(ViewList[Plane]==View && ShwOrTrkList[Plane]>3) {
04774               
04775               for(unsigned int k=0; k<HitList[Plane].size(); ++k) {
04776                 HitList[Plane][k]->SetTrkFlag(2);
04777                 Trk->AddHit(HitList[Plane][k]);
04778               }
04779             }
04780           }
04781         }
04782         
04783         TrackCam* Trku = FinalTrackTempHolder[0].back();
04784         TrackCam* Trkv = FinalTrackTempHolder[1].back();
04785         
04786         Trku->SetPartner(Trkv); 
04787         Trkv->SetPartner(Trku);
04788         
04789         for(unsigned int p=0; p<Trku->GetEntries(); ++p) {HitsInTracks[0].push_back(Trku->GetHit(p));}
04790         for(unsigned int p=0; p<Trkv->GetEntries(); ++p) {HitsInTracks[1].push_back(Trkv->GetHit(p));}
04791       }
04792       
04793     } // End U/V/TrackStrips check
04795 
04796     for(int k=0; k<500; ++k) {
04797       ShwOrTrkList[k]=-1; 
04798       ViewList[k]=-1;
04799       ClusterList[k].clear();
04800       HitList[k].clear();
04801     }
04802 
04803   } // End loop over temporary tracks   
04805  
04806  
04807 
04808   // Check tracks don't contain same hits as a previous track
04810   bool GoodTrack; bool ThreePlanes[2];
04811   int DuplicateStrips[2]; int Planes[2]; int ThisPlane;
04812   vector<HitCam*> HitsInPrevTracks[2];
04813 
04814   for(unsigned int p=0; p<FinalTrackTempHolder[0].size(); ++p) {
04815 
04816     TrackCam* Trku = FinalTrackTempHolder[0][p];
04817     TrackCam* Trkv = FinalTrackTempHolder[1][p];
04818 
04820     if(ModuleType==1) {LookForHitsAcrossGap(Trku); LookForHitsAcrossGap(Trkv);}
04821 
04822     ExtendTrack(Trku); ExtendTrack(Trkv);
04823     FillGapsInTrack(Trku); FillGapsInTrack(Trkv);
04825 
04826 
04828     GoodTrack=false; 
04829 
04830     if(Trku->GetEntries()>2 && Trkv->GetEntries()>2) {
04831       DuplicateStrips[0]=0; DuplicateStrips[1]=0;
04832 
04833       // For u
04834       Planes[0]=-999; Planes[1]=-999; ThreePlanes[0]=false;
04835 
04836       for(unsigned int p=0; p<Trku->GetEntries(); ++p) {
04837         HitCam* TrackHit = Trku->GetHit(p);
04838 
04839         for(unsigned int q=0; q<HitsInPrevTracks[0].size(); ++q) {
04840           if(TrackHit==HitsInPrevTracks[0][q]) {DuplicateStrips[0]++; break;}
04841         }
04842 
04843         ThisPlane=TrackHit->GetPlane();
04844         if(Planes[0]<0) {Planes[0]=ThisPlane;}
04845         else if(Planes[1]<0 && Planes[0]!=ThisPlane) {Planes[1]=ThisPlane;}
04846         else if(Planes[0]!=ThisPlane && Planes[1]!=ThisPlane) {ThreePlanes[0]=true;}
04847 
04848         if(DuplicateStrips[0]>2) {break;}
04849       }
04850 
04851 
04852       // For v
04853       Planes[0]=-999; Planes[1]=-999; ThreePlanes[1]=false;
04854 
04855       for(unsigned int p=0; p<Trkv->GetEntries(); ++p) {
04856         HitCam* TrackHit = Trkv->GetHit(p);
04857 
04858         for(unsigned int q=0; q<HitsInPrevTracks[1].size(); ++q) {
04859           if(TrackHit==HitsInPrevTracks[1][q]) {DuplicateStrips[1]++; break;}
04860         }
04861 
04862         ThisPlane=TrackHit->GetPlane();
04863         if(Planes[0]<0) {Planes[0]=ThisPlane;}
04864         else if(Planes[1]<0 && Planes[0]!=ThisPlane) {Planes[1]=ThisPlane;}
04865         else if(Planes[0]!=ThisPlane && Planes[1]!=ThisPlane) {ThreePlanes[1]=true;}
04866 
04867         if(DuplicateStrips[1]>2) {break;}
04868       }
04869 
04870       if(ThreePlanes[0] && ThreePlanes[1] && DuplicateStrips[0]<3 && DuplicateStrips[1]<3) {GoodTrack=true;}    
04871     }
04873 
04874 
04875     if(GoodTrack) {
04876       FinalTrackBank[0].push_back(Trku);
04877       FinalTrackBank[1].push_back(Trkv);
04878 
04879       for(unsigned int q=0; q<Trku->GetEntries(); ++q) {HitsInPrevTracks[0].push_back(Trku->GetHit(q));}
04880       for(unsigned int q=0; q<Trkv->GetEntries(); ++q) {HitsInPrevTracks[1].push_back(Trkv->GetHit(q));}
04881 
04882       if(ModuleType==2) {MatchUV(Trku,Trkv);}
04883     }
04884   }
04885 
04886   FinalTrackTempHolder[0].clear(); FinalTrackTempHolder[1].clear();
04887   HitsInPrevTracks[0].clear(); HitsInPrevTracks[1].clear();
04889 
04890 
04891 
04892   // Print out list of final tracks           
04894   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF TRACKS *** " << endl;
04895   for(View=0; View<2; ++View) {
04896     for(unsigned int i=0; i<FinalTrackBank[View].size(); ++i) {
04897       TrackCam* TrkTemp = FinalTrackBank[View][i];
04898       
04899       MSG("AlgTrackCamList", Msg::kVerbose) 
04900         << "  " << View << " " << i
04901         << "  " << TrkTemp->GetBegPlane() << "  " << TrkTemp->GetEndPlane() << endl;
04902     }
04903   }
04905 
04906   return;
04907 }

void AlgTrackCamList::FormTheClusters (  ) 

Definition at line 283 of file AlgTrackCamList.cxx.

References ClusterCam::AddHit(), ClusterBank, HitBank, ClusterCam::IsHitAssoc(), NumModules, and PlanesInModule.

Referenced by RunTheFinder().

00284 {
00285   // For each plane where we stored hits, we look to form 1D clusters 
00286   // of these hits. This is controlled by the IsHitAssoc method in ClusterCam.
00287   // Pointers to the ClusterCam objects are stored in plane order in a 
00288   // clusterbank
00289   
00290 
00291   bool AddingHits;
00292   int i;
00293 
00294   for(int Module=0; Module<NumModules; ++Module) {
00295     for(int Plane=0; Plane<PlanesInModule+1; ++Plane) {
00296       i=Plane + Module*(PlanesInModule+1);
00297       
00298       // Loop over the hits on plane i
00299       for(unsigned int j=0; j<HitBank[i].size(); ++j) {
00300         
00301         // Make a cluster for the first hit on the plane
00302         if(HitBank[i][j]->GetUID()==0) {
00303           ClusterCam* Clust = new ClusterCam(HitBank[i][j]);
00304           ClusterBank[HitBank[i][j]->GetPlane()].push_back(Clust);
00305 
00306           // Change UID from 0 to 1 to show that hit has been added       
00307           HitBank[i][j]->SetUID(1); 
00308           AddingHits=true;
00309           
00310           // Loop over other hits on plane to form clusters
00311           while(AddingHits==true) {
00312             AddingHits=false;
00313             
00314             for(unsigned int k=0; k<HitBank[i].size(); ++k) {
00315               if(HitBank[i][k]->GetUID()==0 && Clust->IsHitAssoc(HitBank[i][k])) {
00316                 Clust->AddHit(HitBank[i][k]); 
00317                 HitBank[i][k]->SetUID(1);
00318                 AddingHits=true;
00319               }
00320             }
00321           }
00322 
00323         }
00324       } // End loop over hits on plane
00325     } // End loop over planes in module
00326   } // End loop over modules
00327 
00328   return;
00329 }

void AlgTrackCamList::FormTheHits ( CandSliceHandle slice  ) 

Definition at line 219 of file AlgTrackCamList.cxx.

References AllHitBank, CandHandle::GetDaughterIterator(), HitBank, Msg::kVerbose, CandDeMuxDigit::kXTalk, ModuleType, MSG, PECut, and PECut2.

Referenced by RunTheFinder().

00220 {
00221   // Create and (store pointers to) HitCam objects for all strips in slice
00222   // meeting PH and XTalk requirements. HitCam class is essentially a wrapper
00223   // for CandStripHandles
00224 
00225   bool IsXTalk[2];
00226   int Counter; double StripCharge;
00227   
00228   // Loop over all strips in the slice
00229   TIter SliceItr(slice->GetDaughterIterator());
00230   while(CandStripHandle* SlcStrip = dynamic_cast<CandStripHandle*>(SliceItr())) {
00231     
00232     // Default Xtalk flag to false, as cannot identify ND XTalk
00233     IsXTalk[0]=false; IsXTalk[1]=false;
00234   
00235     // Identify XTalk strips for FD
00236     if(ModuleType==1) {
00237       IsXTalk[0]=true; IsXTalk[1]=true; Counter=0;
00238       
00239       TIter DigItr(SlcStrip->GetDaughterIterator());
00240       while(CandDeMuxDigitHandle* Digit = (CandDeMuxDigitHandle*)(DigItr())) {
00241         if( !((Digit->GetDeMuxDigitFlagWord()<8)
00242               && (Digit->GetDeMuxDigitFlagWord() & CandDeMuxDigit::kXTalk)==(CandDeMuxDigit::kXTalk)) ) {
00243           IsXTalk[Counter]=false;
00244         }
00245         if(Counter==0) {Counter=1;}
00246       }
00247     }
00248     
00249     // Add the hits to the bank, where they are stored in plane order
00250     if(IsXTalk[0]==false || IsXTalk[1]==false) {
00251       StripCharge=SlcStrip->GetCharge();
00252 
00253       if(StripCharge>PECut2 || StripCharge>PECut) {
00254         HitCam* hit = new HitCam(SlcStrip);
00255 
00256         if(StripCharge>PECut) {HitBank[SlcStrip->GetPlane()].push_back(hit);}
00257         if(StripCharge>PECut2) {AllHitBank[SlcStrip->GetPlane()].push_back(hit);}
00258       }
00259     }
00260   }
00261   
00262 
00263   // Print out the list of hits
00264   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF HITS *** " << endl;
00265   for(int i=0; i<500; ++i){
00266     for(unsigned int j=0; j<HitBank[i].size(); ++j) {
00267       MSG("AlgTrackCamList", Msg::kVerbose)
00268         << " pln="  << HitBank[i][j]->GetPlane() 
00269         << " strp=" << HitBank[i][j]->GetStrip() 
00270         << " chg="  << HitBank[i][j]->GetCharge()
00271         << " time=" << HitBank[i][j]->GetTime() << endl;
00272     }
00273   }
00274 
00275   return;
00276 }

void AlgTrackCamList::FormTriplets (  ) 

Definition at line 492 of file AlgTrackCamList.cxx.

References ClusterBank, ClusterCam::GetBegTPos(), ClusterCam::GetEndTPos(), ClusterCam::GetPlane(), Munits::km, Msg::kVerbose, MSG, NumModules, PlanesInModule, SegmentBank, and size.

Referenced by RunTheFinder().

00493 {
00494   // Treating U and V views separately, we consider associations between clusters
00495   // on nearby planes. We look for groups of three clusters, each on different 
00496   // planes, that can be joined together in a small track segment, or triplet. 
00497 
00498   // All possible triplets of the following form are created:
00499   
00500   // If we sit on plane 0, and there are no possible hits on a plane marked X, m 
00501   // means a lower plane number than our origin plane, and p means a higher plane 
00502   // number.
00503    
00504   // m 0 p         The simplest triplet
00505   
00506   // m X 0 p       One gap
00507   // m 0 X p
00508   
00509   // m X X 0 p     Two gaps
00510   // m X 0 X p
00511   // m 0 X X p
00512 
00513 
00514   int TripAssocNum;
00515   bool AlternateTriplets;
00516   bool JoinFlag;
00517   int i;
00518   
00519   // Begin loop to make the triplets
00520   for(int Module=0; Module<NumModules; ++Module) {
00521 
00522     for(int Plane=2; Plane<PlanesInModule-1; ++Plane) {
00523       i=Plane + Module*(PlanesInModule+1);
00524 
00525       // Loop over clusters on plane. This is our plane 0, the origin.
00526       for(unsigned int k0=0; k0<ClusterBank[i].size(); ++k0) {
00527 
00528         if(ClusterBank[i][k0]->GetTrkFlag()>0) {
00529           JoinFlag=false;
00530 
00531 
00532           // Look for triplets of form m 0 p
00534           if((i-2)>=0 && ClusterBank[i-2].size()>0 && (i+2)<500 && ClusterBank[i+2].size()>0) {
00535           
00536             // Look at the clusters on plane i-2...
00537             for(unsigned int km=0; km<ClusterBank[i-2].size(); ++km) {
00538             
00539               // ...and the clusters on plane i+2
00540               for(unsigned int kp=0; kp<ClusterBank[i+2].size(); ++kp) {
00541               
00542                 if(ClusterBank[i-2][km]->GetTrkFlag()>0 && ClusterBank[i+2][kp]->GetTrkFlag()>0) {
00543 
00544                   // Check the track-like association of the three clusters
00545                   TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+2][kp]);
00546 
00547                   // If the association is good, make the triplet
00548                   if(TripAssocNum>0) {
00549                     TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-2][km], ClusterBank[i][k0], ClusterBank[i+2][kp]);
00550 
00551                     SegmentBank[i].push_back(seg0);
00552                     JoinFlag=true;
00553 
00554                     // Indicate that triplet is exceptionally track-like, by setting flag to 2.
00555                     if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
00556                   }
00557                 }
00558               }
00559             }
00560           }
00562         
00563 
00564         
00565           // Look for triplets of form m <-> 0 p
00567           if((i+2)<500 && ClusterBank[i+2].size()>0) {
00568 
00569         
00570             // Look for triplets of form m X 0 p  
00572             if( Plane>3 && (i-4)>=0 && ClusterBank[i-4].size()>0) {
00573 
00574               // Look at the clusters on plane i+2...
00575               for(unsigned int kp=0; kp<ClusterBank[i+2].size(); ++kp) {
00576               
00577                 // ...and the clusters on plane i-4
00578                 for(unsigned int km=0; km<ClusterBank[i-4].size(); ++km) {
00579 
00580                   if(ClusterBank[i-4][km]->GetTrkFlag()>0 && ClusterBank[i+2][kp]->GetTrkFlag()>0) {
00581 
00582                     // Check the track-like association of the three clusters
00583                     TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i+2][kp]);
00584                   
00585                     // If the association is good, check whether we can make the triplet
00586                     if(TripAssocNum>0) {
00587                       AlternateTriplets=false;
00588                       
00589                       // Check to see if there is also an alternate triplet possibility with clusters on plane i-2
00590                       // i.e. Want to make sure that the X in "m X 0 p" is correct.
00591 
00592                       if(AlternateTriplets==false && ClusterBank[i-2].size()>0) {
00593                         // Look at the clusters on plane i-2
00594                         for(unsigned int ktmp=0; ktmp<ClusterBank[i-2].size(); ++ktmp) {
00595                           if(AlternateTriplets==false && ClusterBank[i-2][ktmp]->GetTrkFlag()>0
00596                              && ClusterBank[i-2][ktmp]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i][k0])>0 
00597                              && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp],ClusterBank[i+2][kp])>0) 
00598                             {AlternateTriplets=true;}
00599                         }
00600                       }
00601                       
00602                       
00603                       // If everything is good, make the triplet
00604                       if(AlternateTriplets==false) {
00605                         TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-4][km],ClusterBank[i][k0],ClusterBank[i+2][kp]);
00606                         
00607                         SegmentBank[i].push_back(seg0);
00608                         JoinFlag=true;
00609         
00610                         // Indicate that triplet is exceptionally track-like, by setting flag to 2.
00611                         if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
00612                       }
00613                     }
00614                   }
00615                 }
00616               }
00617             }
00619 
00620 
00621             // Look for triplets of form m X X 0 p  
00623             if( Plane>5 && (i-6)>=0 && ClusterBank[i-6].size()>0) {
00624 
00625               // Look at the clusters on plane i+2...
00626               for(unsigned int kp=0; kp<ClusterBank[i+2].size(); ++kp) {
00627               
00628                 // ...and look at the clusters on plane i-6
00629                 for(unsigned int km=0; km<ClusterBank[i-6].size(); ++km) {
00630                   
00631                   if(ClusterBank[i-6][km]->GetTrkFlag()>0 && ClusterBank[i+2][kp]->GetTrkFlag()>0) {
00632 
00633                     // Check the track-like association of the three clusters
00634                     TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-6][km],ClusterBank[i+2][kp]);
00635         
00636                     // If the association is good, check whether we can make the triplet
00637                     if(TripAssocNum>0) {
00638                       AlternateTriplets=false;
00639                       
00640                       // Check to see if there are alternate triplet possibilities with clusters on plane i-2 or i-4
00641                       // i.e. Want to make sure that the Xs in "m X X 0 p" are correct.
00642 
00643                       // Check first X
00644                       if(AlternateTriplets==false && ClusterBank[i-2].size()>0) {
00645                         // Look at any clusters on plane i-2
00646                         for(unsigned int ktmp=0; ktmp<ClusterBank[i-2].size(); ++ktmp) {
00647                           if(AlternateTriplets==false && ClusterBank[i-2][ktmp]->GetTrkFlag()>0
00648                              && ClusterBank[i-2][ktmp]->IsTrkAssoc(ClusterBank[i-6][km],ClusterBank[i][k0])>0 
00649                              && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp],ClusterBank[i+2][kp])>0) 
00650                             {AlternateTriplets=true;}
00651                         }
00652                       }
00653                       
00654                       // Check second X, plane i-4
00655                       if(AlternateTriplets==false && ClusterBank[i-4].size()>0) {
00656                         // Look at any clusters on plane i-4
00657                         for(unsigned int ktmp=0; ktmp<ClusterBank[i-4].size(); ++ktmp) {
00658                           if(AlternateTriplets==false && ClusterBank[i-4][ktmp]->GetTrkFlag()>0
00659                              && ClusterBank[i-4][ktmp]->IsTrkAssoc(ClusterBank[i-6][km],ClusterBank[i][k0])>0 
00660                              && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-4][ktmp],ClusterBank[i+2][kp])>0) 
00661                             {AlternateTriplets=true;}
00662                         }
00663                       }
00664                       
00665                       // Check both Xs together, planes i-2 and i-4
00666                       if(AlternateTriplets==false && ClusterBank[i-2].size()>0 && ClusterBank[i-4].size()>0) {
00667                         // Look again at any clusters on plane i-4...
00668                         for(unsigned int ktmp=0; ktmp<ClusterBank[i-4].size(); ++ktmp) {
00669                           
00670                           // ...and look again any clusters on plane i-2
00671                           for(unsigned int ktmp1=0; ktmp1<ClusterBank[i-2].size(); ++ktmp1) {
00672                             if(AlternateTriplets==false 
00673                                && ClusterBank[i-4][ktmp]->GetTrkFlag()>0 && ClusterBank[i-2][ktmp1]->GetTrkFlag()>0
00674                                && ClusterBank[i-4][ktmp]->IsTrkAssoc(ClusterBank[i-6][km],ClusterBank[i-2][ktmp1])>0
00675                                && ClusterBank[i-2][ktmp1]->IsTrkAssoc(ClusterBank[i-4][ktmp],ClusterBank[i][k0])>0
00676                                && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp1],ClusterBank[i+2][kp])>0)
00677                               {AlternateTriplets=true;}
00678                           }
00679                           
00680                         }
00681                       }
00682                       
00683                       // If everything is good, make the triplet
00684                       if(AlternateTriplets==false) {
00685                         TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-6][km],ClusterBank[i][k0],ClusterBank[i+2][kp]);
00686                         
00687                         SegmentBank[i].push_back(seg0);
00688                         JoinFlag=true;
00689 
00690                         // Indicate that triplet is exceptionally track-like, by setting flag to 2.
00691                         if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
00692                       }
00693                     }
00694                   }
00695                 }
00696               }
00697             }
00699 
00700           }
00702 
00703 
00704 
00705           // Look for triplets of form m 0 <-> p
00707           if((i-2)>0 && ClusterBank[i-2].size()>0) {
00708 
00709         
00710             // Look for triplets of form m 0 X p  
00712             if( Plane<PlanesInModule-3 && (i+4)<500 && ClusterBank[i+4].size()>0) {
00713 
00714               // Look at the clusters on plane i-2...
00715               for(unsigned int km=0; km<ClusterBank[i-2].size(); ++km) {
00716               
00717                 // ... and the clusters on plane i+4
00718                 for(unsigned int kp=0; kp<ClusterBank[i+4].size(); ++kp) {
00719 
00720                   if(ClusterBank[i-2][km]->GetTrkFlag()>0 && ClusterBank[i+4][kp]->GetTrkFlag()>0) {
00721 
00722                     // Check the track-like association of the three clusters
00723                     TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+4][kp]);
00724                    
00725                     // If the association is good, check whether we can make the triplet
00726                     if(TripAssocNum>0) {
00727                       AlternateTriplets=false;
00728                       
00729                       // Check to see if there is also an alternate triplet possibility with clusters on plane i+2
00730                       // i.e. Want to make sure that the X in "m 0 X p" is correct.
00731 
00732                       if(AlternateTriplets==false && ClusterBank[i+2].size()>0) {
00733                         // Look at the clusters on plane i+2
00734                         for(unsigned int ktmp=0; ktmp<ClusterBank[i+2].size(); ++ktmp) {
00735                           if(AlternateTriplets==false && ClusterBank[i+2][ktmp]->GetTrkFlag()>0
00736                              && ClusterBank[i+2][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+4][kp])>0 
00737                              && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+2][ktmp])>0) 
00738                             {AlternateTriplets=true;}
00739                         }
00740                       }
00741                       
00742                       
00743                       // If everything is good, make the triplet
00744                       if(AlternateTriplets==false) {
00745                         // Store segment on empty plane too
00746                         for(int j=0; j<2; ++j) {
00747                           TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-2][km],ClusterBank[i][k0],ClusterBank[i+4][kp]);
00748                           
00749                           SegmentBank[i+2*j].push_back(seg0);
00750                           JoinFlag=true;
00751                         }
00752                         
00753                         // Indicate that triplet is exceptionally track-like, by setting flag to 2.
00754                         if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
00755                       }
00756                     }
00757                   }
00758                 }
00759               }
00760             }
00762 
00763 
00764             // Look for triplets of form m 0 X X p  
00766             if(Plane<PlanesInModule-5 && (i+6)<500 && ClusterBank[i+6].size()>0) {
00767 
00768               // Look at the clusters on plane i-2...
00769               for(unsigned int km=0; km<ClusterBank[i-2].size(); ++km) {
00770               
00771                 // ...and look at the clusters on plane i+6
00772                 for(unsigned int kp=0; kp<ClusterBank[i+6].size(); ++kp) {
00773 
00774                   if(ClusterBank[i-2][km]->GetTrkFlag()>0 && ClusterBank[i+6][kp]->GetTrkFlag()>0) {
00775                     // Check the track-like association of the three clusters
00776                     TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+6][kp]);
00777                     
00778                     // If the association is good, check whether we can make the triplet
00779                     if(TripAssocNum>0) {
00780                       AlternateTriplets=false;
00781         
00782                       // Check to see if there are alternate triplet possibilities with clusters on plane i+2 or i+4
00783                       // i.e. Want to make sure that the Xs in "m 0 X X p" are correct.
00784               
00785                       // Check first X, plane i+2
00786                       if(AlternateTriplets==false && ClusterBank[i+2].size()>0) {
00787                         // Look at any clusters on plane i+2
00788                         for(unsigned int ktmp=0; ktmp<ClusterBank[i+2].size(); ++ktmp) {
00789                           if(AlternateTriplets==false && ClusterBank[i+2][ktmp]->GetTrkFlag()>0
00790                              && ClusterBank[i+2][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+6][kp])>0 
00791                              && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+2][ktmp])>0) 
00792                             {AlternateTriplets=true;}
00793                         }
00794                       }
00795                       
00796                       // Check second X, plane i+4
00797                       if(AlternateTriplets==false && ClusterBank[i+4].size()>0) {
00798                         // Look at any clusters on plane i+4
00799                         for(unsigned int ktmp=0; ktmp<ClusterBank[i+4].size(); ++ktmp) {
00800                           if(AlternateTriplets==false && ClusterBank[i+4][ktmp]->GetTrkFlag()>0
00801                              && ClusterBank[i+4][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+6][kp])>0 
00802                              && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+4][ktmp])>0) 
00803                             {AlternateTriplets=true;}
00804                         }
00805                       }
00806                       
00807                       // Check both Xs together, planes i+2 and i+4
00808                       if(AlternateTriplets==false && ClusterBank[i+2].size()>0 && ClusterBank[i+4].size()>0) {
00809                         // Look again at any clusters on plane i+2...
00810                         for(unsigned int ktmp=0; ktmp<ClusterBank[i+2].size(); ++ktmp) {
00811                           
00812                           // ...and look again at any clusters on plane i+4
00813                           for(unsigned int ktmp1=0; ktmp1<ClusterBank[i+4].size(); ++ktmp1) {
00814                             if(AlternateTriplets==false 
00815                                && ClusterBank[i+2][ktmp]->GetTrkFlag()>0 &&ClusterBank[i+4][ktmp1]->GetTrkFlag()>0
00816                                && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+2][ktmp])>0
00817                                && ClusterBank[i+2][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+4][ktmp1])>0
00818                                && ClusterBank[i+4][ktmp1]->IsTrkAssoc(ClusterBank[i+2][ktmp],ClusterBank[i+6][kp])>0)
00819                               {AlternateTriplets=true;}
00820                           }
00821                           
00822                         }
00823                       }
00824                       
00825                       // If everything is good, make the triplet
00826                       if(AlternateTriplets==false) {
00827                         // Store segment on empty planes too
00828                         for(int j=0; j<3; ++j) {
00829                           TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-2][km],ClusterBank[i][k0],ClusterBank[i+6][kp]);
00830                           
00831                           SegmentBank[i+2*j].push_back(seg0);
00832                           JoinFlag=true;
00833                         }
00834                         // Indicate that triplet is exceptionally track-like, by setting flag to 2.
00835                         if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
00836                       }
00837                     }
00838                   }
00839                 }
00840               }
00841             }
00843 
00844           }
00846 
00847 
00848 
00849           // Look for triplets of form m X 0 X p
00851           if(Plane>3 && Plane<PlanesInModule-3 && (i-4)>=0 && ClusterBank[i-4].size()>0 && (i+4)<500 && ClusterBank[i+4].size()>0) {
00852 
00853             // Look at clusters on planes i-4 and i+4
00854             for(unsigned int km=0; km<ClusterBank[i-4].size(); ++km) {
00855               for(unsigned int kp=0; kp<ClusterBank[i+4].size(); ++kp) {
00856 
00857                 if(ClusterBank[i-4][km]->GetTrkFlag()>0 && ClusterBank[i+4][kp]->GetTrkFlag()>0) {
00858                   // Check the track-like association of the three clusters
00859                   TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i+4][kp]);
00860                   
00861                   // If the association is good, check whether we can make the triplet
00862                   if(TripAssocNum>0) {
00863                     AlternateTriplets=false;
00864                     
00865                     // Check to see if there are alternate triplet possibilities with clusters on plane i-2 or i+2
00866                     // i.e. Want to make sure that the Xs in "m X 0 X p" are correct.
00867                     
00868                     // Check first X, plane i-2
00869                     if(AlternateTriplets==false && ClusterBank[i-2].size()>0) {
00870                       // Look at any clusters on plane i-2
00871                       for(unsigned int ktmp=0; ktmp<ClusterBank[i-2].size(); ++ktmp) {
00872                         if(AlternateTriplets==false && ClusterBank[i-2][ktmp]->GetTrkFlag()>0
00873                            && ClusterBank[i-2][ktmp]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i][k0])>0 
00874                            && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp],ClusterBank[i+4][kp])>0) 
00875                           {AlternateTriplets=true;}
00876                       }
00877                     }
00878                     
00879                     // Check second X, plane i+2
00880                     if(AlternateTriplets==false && ClusterBank[i+2].size()>0) {
00881                       // Look at any clusters on plane i+2
00882                       for(unsigned int ktmp=0; ktmp<ClusterBank[i+2].size(); ++ktmp) {
00883                         if(AlternateTriplets==false  && ClusterBank[i+2][ktmp]->GetTrkFlag()>0
00884                            && ClusterBank[i+2][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+4][kp])>0 
00885                            && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i+2][ktmp])>0) 
00886                           {AlternateTriplets=true;}
00887                       }
00888                     }
00889                     
00890                     // Check both Xs together, planes i-2 and i+2
00891                     if(AlternateTriplets==false && ClusterBank[i-2].size()>0 && ClusterBank[i+2].size()>0) {
00892                       
00893                       // Loop again at any clusters on plane i-2
00894                       for(unsigned int ktmp=0; ktmp<ClusterBank[i-2].size(); ++ktmp) {
00895                         
00896                         // Look again at any clusters on plane i+2
00897                         for(unsigned int ktmp1=0; ktmp1<ClusterBank[i+2].size(); ++ktmp1) {
00898                           if(AlternateTriplets==false 
00899                              && ClusterBank[i-2][ktmp]->GetTrkFlag()>0 && ClusterBank[i+2][ktmp1]->GetTrkFlag()>0
00900                              && ClusterBank[i-2][ktmp]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i][k0])>0 
00901                              && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp],ClusterBank[i+2][ktmp1])>0 
00902                              && ClusterBank[i+2][ktmp1]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+4][kp])>0)
00903                             {AlternateTriplets=true;}
00904                         }
00905                       }
00906                     }
00907 
00908                     // If everything is good, make the triplet
00909                     if(AlternateTriplets==false) {
00910                       // Store segment on empty planes too
00911                       for(int j=0; j<2; ++j) {
00912                         TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-4][km],ClusterBank[i][k0],ClusterBank[i+4][kp]);
00913                         
00914                         SegmentBank[i+2*j].push_back(seg0);
00915                         JoinFlag=true;
00916                       }
00917                       // Indicate that triplet is exceptionally track-like, by setting flag to 2.
00918                       if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
00919                     }               
00920                   }
00921                 }
00922               }
00923             }
00924           }
00926 
00927 
00928           // Set the exceptionally track-like flag for lone hits that form part of a triplet
00929           if(JoinFlag==true && (ClusterBank[i][k0]->GetEndTPos()==ClusterBank[i][k0]->GetBegTPos())) {
00930             ClusterBank[i][k0]->SetTrkFlag(2);
00931           }
00932 
00933           
00934         }
00935 
00936       } // End loop over clusters on plane
00937     
00938     } // End loop over planes in module
00939 
00940   } // End loop over modules
00941 
00942 
00943   // Print out list of hits and 1D clusters
00944   MSG("AlgTrackCamList", Msg::kVerbose) << " *** 2ND LIST OF CLUSTERS *** " << endl;
00945   for(int i=0; i<500; ++i){
00946     for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
00947       MSG("AlgTrackCamList", Msg::kVerbose)
00948         << " plane="    << ClusterBank[i][j]->GetPlane() 
00949         << " begtpos="  << ClusterBank[i][j]->GetBegTPos() 
00950         << " endtpos="  << ClusterBank[i][j]->GetEndTPos() 
00951         << " trkflag="  << ClusterBank[i][j]->GetTrkFlag()
00952         << " shwpln="   << ClusterBank[i][j]->GetShwPlnFlag() 
00953         << " trkpln="   << ClusterBank[i][j]->GetTrkPlnFlag()
00954         << " shwflag="  << ClusterBank[i][j]->GetShwFlag() << endl;
00955     }
00956   }
00957 
00958 
00959   // Print out list of triplets
00960   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF TRIPLETS *** " << endl;
00961   for(int i=0; i<500; ++i) {
00962     for(unsigned int j=0;j<SegmentBank[i].size(); ++j) {
00963       MSG("AlgTrackCamList", Msg::kVerbose) << " Segment Num " << j << ", Plane " << i <<endl;
00964  
00965       ClusterCam* clr0 = SegmentBank[i][j]->GetCluster(0);
00966       ClusterCam* clr1 = SegmentBank[i][j]->GetCluster(1);
00967       ClusterCam* clr2 = SegmentBank[i][j]->GetCluster(2);
00968 
00969       MSG("AlgTrackCamList", Msg::kVerbose) 
00970         << " (" << clr0->GetPlane() << "," << clr0->GetBegTPos() << "->" << clr0->GetEndTPos() << ") "
00971         << " (" << clr1->GetPlane() << "," << clr1->GetBegTPos() << "->" << clr1->GetEndTPos() << ") "
00972         << " (" << clr2->GetPlane() << "," << clr2->GetBegTPos() << "->" << clr2->GetEndTPos() << ") " << endl;
00973     }
00974   }
00975   
00976   return;
00977 }

void AlgTrackCamList::IDTrkAndShwClusters (  ) 

Definition at line 336 of file AlgTrackCamList.cxx.

References ClusterBank, OscFit::GetCharge(), ClusterCam::GetCharge(), ClusterCam::GetEntries(), ClusterCam::IsShwAssoc(), Msg::kVerbose, MSG, NumModules, PECut, PlanesInModule, ClusterCam::SetShwPlnFlag(), and ClusterCam::SetTrkPlnFlag().

Referenced by RunTheFinder().

00337 {
00338   // Look at the 1D clusters formed and use simple techniques to get an idea
00339   // of how track-like or shower-like each cluster is.
00340   // Then look at all clusters on the plane and see how track-like or shower-like
00341   // the plane is.
00342 
00343   //  Detector::Detector_t Detector = vldc->GetDetector();
00344 
00345   int i, k0;
00346   int nclust0, nclust1, nhits0, nhits1, ShwAssocNum;
00347   vector<ClusterCam*> TempClust0;
00348   vector<ClusterCam*> TempClust1;
00349 
00350 
00351   // Simple checks to set initial track flags
00353   for(int Module=0; Module<NumModules; ++Module) {
00354     for(int Plane=0; Plane<PlanesInModule+1; ++Plane) {
00355       i=Plane + Module*(PlanesInModule+1);
00356       
00357       for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
00358         if(ClusterBank[i][j]->GetCharge()>PECut || ClusterBank[i][j]->GetDigits()>1 ) {
00359           ClusterBank[i][j]->SetTrkFlag(1);
00360         }
00361       }
00362     }
00363   }
00365 
00366 
00367   
00368   // Identify the shower-like clusters
00370   for(int Module=0; Module<NumModules; ++Module) {
00371     for(int Plane=1; Plane<PlanesInModule+1; ++Plane) {
00372       i=Plane + Module*(PlanesInModule+1);
00373       
00374       for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
00375         ClusterCam* Clust0 = ClusterBank[i][j];
00376 
00377         nhits0=0; nhits1=0; nclust0=0; nclust1=0;
00378         TempClust0.clear(); TempClust1.clear(); 
00379         
00380         
00381         // Look for shower associations in nearby planes
00383         // Loop over nearby planes
00384         for(int k=-1; k<2; ++k) {
00385           k0=i+(2*k);
00386           
00387           // Check k0 is a valid nearby plane
00388           if(k0>Module*(PlanesInModule+1) && k0<(Module+1)*(PlanesInModule+1)) {
00389 
00390             // Loop over clusters on plane k0
00391             for(unsigned int k1=0; k1<ClusterBank[k0].size(); ++k1) {
00392               ClusterCam* Clust1 = ClusterBank[k0][k1];
00393 
00394               // Check shower-like associations
00395               if(Clust0->GetCharge()>1. && Clust1->GetCharge()>1.) {
00396                 ShwAssocNum=Clust0->IsShwAssoc(Clust1);
00397                 
00398                 // Store the association data and the clusters
00399                 if(ShwAssocNum>0) {
00400                   nhits0+=Clust1->GetEntries(); 
00401                   nclust0+=1;
00402                   TempClust0.push_back(Clust1);
00403                 }
00404                 
00405                 if(ShwAssocNum>1) {
00406                   nhits1+=Clust1->GetEntries(); 
00407                   nclust1+=1;
00408                   TempClust1.push_back(Clust1);
00409                 }
00410               }
00411             } // End loop over clusters on nearby plane
00412           }
00413         } // End loop over nearby planes
00415 
00416 
00417         // Use the above results to set the shower flags
00418         // Has important implications for finding tracks embedded in showers
00420         if(nclust1>1 && nhits1>4) {
00421           for(unsigned int k=0; k<TempClust1.size(); ++k) {
00422             TempClust1[k]->SetShwFlag(1);
00423           }
00424         }
00425 
00426         if(nclust0>4 && nhits0>5) {
00427           for(unsigned int k=0; k<TempClust0.size(); ++k) {
00428             TempClust0[k]->SetShwFlag(1);
00429           }
00430         }
00432 
00433       }
00434     }
00435   }
00437 
00438 
00439 
00440   // Identify track-like and shower-like planes
00442   double Charge;
00443 
00444   for(int Module=0; Module<NumModules; ++Module) {
00445     for(int Plane=0; Plane<PlanesInModule+1; ++Plane) {
00446       i=Plane + Module*(PlanesInModule+1);
00447       Charge=0.;
00448 
00449       // Get total charge in clusters on plane
00450       for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
00451         Charge+=ClusterBank[i][j]->GetCharge();
00452       }
00453 
00454       if(Charge>0.) {
00455         for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
00456           ClusterCam* Clust = ClusterBank[i][j];
00457 
00458           // Set the plane flags for the clusters // should be <3
00459           if(Clust->GetEntries()<7 && (Clust->GetCharge()/Charge)>0.5) {Clust->SetTrkPlnFlag(1);}
00460           if(Clust->GetEntries()>1 && (Clust->GetCharge()/Charge)>0.5) {Clust->SetShwPlnFlag(1);}
00461         }
00462       }
00463     }
00464   }
00466   
00467 
00468 
00469   // Print out list of hits and 1D clusters
00470   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF CLUSTERS *** " << endl;
00471   for(int i=0; i<500; ++i){
00472     for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
00473       MSG("AlgTrackCamList", Msg::kVerbose)
00474         << " plane="    << ClusterBank[i][j]->GetPlane() 
00475         << " begtpos="  << ClusterBank[i][j]->GetBegTPos() 
00476         << " endtpos="  << ClusterBank[i][j]->GetEndTPos() 
00477         << " trkflag="  << ClusterBank[i][j]->GetTrkFlag()
00478         << " shwpln="   << ClusterBank[i][j]->GetShwPlnFlag() 
00479         << " trkpln="   << ClusterBank[i][j]->GetTrkPlnFlag()
00480         << " shwflag="  << ClusterBank[i][j]->GetShwFlag() << endl;
00481     }
00482   }
00483 
00484   return;
00485 }

void AlgTrackCamList::Join2DTracks (  ) 

Definition at line 3502 of file AlgTrackCamList.cxx.

References TrackSegmentCam::AddSegment(), TrackSegmentCam::GetBegPlane(), TrackSegmentCam::GetBegTPos(), TrackSegmentCam::GetBegZPos(), SntpHelpers::GetCluster(), TrackSegmentCam::GetCluster(), TrackSegmentCam::GetEndPlane(), TrackSegmentCam::GetEndTPos(), TrackSegmentCam::GetEndZPos(), TrackSegmentCam::GetEntries(), ClusterCam::GetPlane(), TrackSegmentCam::GetUID(), TrackSegmentCam::IsAssoc(), Msg::kVerbose, ModuleType, MSG, Particle::Overlap(), PossibleJoins, TrackSegmentCam::SetUID(), size, and ViewSegBank.

Referenced by RunTheFinder().

03503 {
03504   MSG("AlgTrackCamList", Msg::kVerbose) << " *** JOIN 2D TRACKS *** " << endl;
03505 
03506   double ShortestDist; double LongestAllowedDist=5.;
03507   double DeltaTPos, DeltaZPos, Dist;
03508   int DeltaPlane, Length, UID, NumAgree;
03509 
03510   TrackSegmentCam* Seg1ToAdd=0;
03511   TrackSegmentCam* Seg2ToAdd=0;
03512   bool PossibleAssocs, Cont, Overlap, DiffGrad;
03513 
03514   if(ModuleType==2) {LongestAllowedDist=1.;} //Leigh: Changed from 2.0m
03515 
03516 
03517   // Deal with two nearby tracks in ND  
03519   if(ModuleType==2 && PossibleJoins[0].size()==2 && PossibleJoins[1].size()==2) {
03520     Overlap=false;
03521 
03522     // Check first to see if the 2D tracks overlap and share clusters
03524     for(int View=0; View<2; ++View) {
03525       NumAgree=0; 
03526 
03527       for(unsigned int i=0; i<PossibleJoins[View][0]->GetEntries(); ++i) {
03528         ClusterCam* Clust = PossibleJoins[View][0]->GetCluster(i);
03529         
03530         for(unsigned int j=0; j<PossibleJoins[View][1]->GetEntries(); ++j) {
03531           if(PossibleJoins[View][1]->GetCluster(j)==Clust) {Overlap=true; break;}
03532         }
03533         if(Overlap) {break;}
03534 
03535       }
03536       if(Overlap) {break;}
03537     }
03539 
03540     
03541 
03542     // Check for very different gradients
03544     DiffGrad=false;
03545 
03546     for(int View=0; View<2; ++View) {
03547       if(PossibleJoins[View][0]->GetBegPlane()<PossibleJoins[View][1]->GetBegPlane()) {
03548         if(fabs(PossibleJoins[View][0]->GetEndDir()
03549                 -PossibleJoins[View][1]->GetBegDir())>0.7) {DiffGrad=true;}
03550       }
03551       else {
03552         if(fabs(PossibleJoins[View][1]->GetEndDir()
03553                 -PossibleJoins[View][0]->GetBegDir())>0.7) {DiffGrad=true;}  
03554       }
03555       if(DiffGrad) {break;}
03556     }
03557 
03558     if(DiffGrad) {
03559       PossibleJoins[0].clear(); 
03560       PossibleJoins[1].clear();
03561       return;
03562     }
03564 
03565     
03566 
03568     if(Overlap==false) {
03569       // Check associations apply in both views
03570       if( (PossibleJoins[0][0]->IsAssoc(PossibleJoins[0][1])
03571            || PossibleJoins[0][1]->IsAssoc(PossibleJoins[0][0]))
03572           && (PossibleJoins[1][0]->IsAssoc(PossibleJoins[1][1])
03573               || PossibleJoins[1][1]->IsAssoc(PossibleJoins[1][0])) ) {
03574       
03575         PossibleJoins[0][0]->AddSegment(PossibleJoins[0][1]); 
03576         PossibleJoins[0][1]->SetUID(-1);
03577 
03578         PossibleJoins[1][0]->AddSegment(PossibleJoins[1][1]); 
03579         PossibleJoins[1][1]->SetUID(-1);        
03580       }
03581 
03582       // Clear up and leave Join2DTracks method
03583       PossibleJoins[0].clear(); PossibleJoins[1].clear();
03584     
03585       return;
03586     }
03588   }
03590 
03591 
03592 
03594   else if(ModuleType==2 && PossibleJoins[0].size()>2 && PossibleJoins[1].size()>2 
03595           && PossibleJoins[0].size()==PossibleJoins[1].size()) {
03596     
03597     const unsigned int Size = PossibleJoins[0].size();
03598 
03599     bool UView[100], VView[100]; bool AllMatch=true;
03600     int BegPlane, EndPlane;
03601 
03602     if(Size<100) {
03603       for(unsigned int i=0; i<Size; ++i) {UView[i]=false; VView[i]=false;}
03604 
03605       for(unsigned int i=0; i<Size; ++i) {
03606         BegPlane=PossibleJoins[0][i]->GetBegPlane();
03607         EndPlane=PossibleJoins[0][i]->GetEndPlane();
03608         
03609         for(unsigned int j=0; j<Size; ++j) {
03610           
03611           if(VView[j]) {continue;}
03612 
03613           if(abs(BegPlane-PossibleJoins[1][j]->GetBegPlane())<8 
03614              && abs(EndPlane-PossibleJoins[1][j]->GetEndPlane())<8) {
03615             
03616             UView[i]=true; VView[j]=true; break;
03617           }
03618         }
03619       }
03620       
03621       for(unsigned int i=0; i<Size; ++i) {if(UView[i]==false || VView[i]==false) {AllMatch=false;}}
03622       
03623       if(AllMatch) {
03624         PossibleJoins[0].clear(); 
03625         PossibleJoins[1].clear(); 
03626         return;
03627       }
03628     }
03629   }
03631 
03632 
03633 
03634   // First, make the best possible joins between the tracks
03636   for(int View=0; View<2; ++View) {
03637     PossibleAssocs=true;
03638 
03639     while(PossibleAssocs) {
03640       PossibleAssocs=false; 
03641       ShortestDist=LongestAllowedDist; Seg1ToAdd=0; Seg2ToAdd=0; 
03642 
03643       // Consider all pairs of 2D tracks and find the single best association.
03644       for(unsigned int i=0; i<PossibleJoins[View].size(); ++i) {
03645         TrackSegmentCam* Seg1 = PossibleJoins[View][i];
03646 
03647         if(Seg1->GetUID()>2) {
03648 
03649           for(unsigned int j=0; j<PossibleJoins[View].size(); ++j) {
03650             TrackSegmentCam* Seg2 = PossibleJoins[View][j];
03651           
03652             if(Seg2->GetUID()>2 && Seg1!=Seg2) {         
03653 
03654               // Calculate distance between segments
03655               DeltaTPos=fabs(Seg2->GetBegTPos()-Seg1->GetEndTPos());
03656               DeltaZPos=fabs(Seg2->GetBegZPos()-Seg1->GetEndZPos());
03657               DeltaPlane=Seg2->GetBegPlane()-Seg1->GetEndPlane();
03658 
03659               Dist=pow(DeltaTPos*DeltaTPos + DeltaZPos*DeltaZPos,0.5);
03660 
03661               if(DeltaPlane>=0 && Dist<ShortestDist)  {
03662                 if(Seg1->IsAssoc(Seg2)) {Seg1ToAdd=Seg1; Seg2ToAdd=Seg2; ShortestDist=Dist;}
03663               }
03664               else if(DeltaPlane<=0 && Dist<ShortestDist) {
03665                 if(Seg2->IsAssoc(Seg1)) {Seg1ToAdd=Seg1; Seg2ToAdd=Seg2; ShortestDist=Dist;}
03666               }
03667             }
03668           }
03669         }
03670       }
03671 
03672 
03673       // If a good association is found, add the segments and look for the next association.
03674       if(ShortestDist<LongestAllowedDist && Seg2ToAdd!=0 && Seg1ToAdd!=0) { 
03675         MSG("AlgTrackCamList", Msg::kVerbose) << "Missing 2D track assocation: Seg1 " 
03676                                               << Seg1ToAdd->GetBegPlane() << " " << Seg1ToAdd->GetEndPlane() 
03677                                               << " (" << Seg1ToAdd->GetBegTPos() << "," << Seg1ToAdd->GetEndTPos() 
03678                                               << "), Seg2 " << Seg2ToAdd->GetBegPlane() 
03679                                               << " " << Seg2ToAdd->GetEndPlane() << " (" << Seg2ToAdd->GetBegTPos() 
03680                                               << "," << Seg2ToAdd->GetEndTPos() << endl;
03681 
03682         Seg1ToAdd->AddSegment(Seg2ToAdd); Seg2ToAdd->SetUID(-1); PossibleAssocs=true;
03683       }
03684     }
03686 
03687 
03688 
03689     // Now find the longest possible 2D track and remove 2d tracks that overlap.
03691     Cont=true; 
03692 
03693     while(Cont) {
03694       Cont=false;
03695 
03696       TrackSegmentCam* BestTrk = 0;
03697       Length=0;
03698 
03699       for(unsigned int i=0; i<PossibleJoins[View].size(); ++i) {
03700         TrackSegmentCam* Seg = PossibleJoins[View][i];
03701         if(Seg->GetUID()>0 && Seg->GetUID()<4 && (Seg->GetEndPlane()-Seg->GetBegPlane())>Length) {
03702           BestTrk=Seg; Length=Seg->GetEndPlane()-Seg->GetBegPlane();
03703         }
03704       }
03705 
03706 
03707       if(BestTrk) {
03708         BestTrk->SetUID(4);
03709 
03710         for(unsigned int i=0; i<BestTrk->GetEntries(); ++i) {
03711           ClusterCam* BestTrkClust = BestTrk->GetCluster(i);
03712         
03713           if(BestTrkClust->GetPlane()==BestTrk->GetBegPlane() || BestTrkClust->GetPlane()==BestTrk->GetEndPlane()) {continue;}
03714 
03715           for(unsigned int j=0; j<PossibleJoins[View].size(); ++j) {
03716             TrackSegmentCam* Seg = PossibleJoins[View][j];
03717             NumAgree=0; 
03718 
03719             if(Seg->GetUID()<0 || Seg->GetUID()>3) {continue;}
03720   
03721             for(unsigned int k=0; k<Seg->GetEntries(); ++k) {
03722               ClusterCam* SegClust = Seg->GetCluster(k);
03723 
03724               if(SegClust==BestTrkClust) {NumAgree++;}
03725 
03726               if(NumAgree>2) {Seg->SetUID(-1); break;}
03727             }
03728           }
03729         }
03730       }
03731 
03732       // Conditions for continuing here
03733       for(unsigned int i=0; i<PossibleJoins[View].size(); ++i) {
03734         UID=PossibleJoins[View][i]->GetUID();
03735 
03736         if(UID>0 && UID<4) {Cont=true;}
03737       }
03738       
03739     }
03741 
03742    
03743     // Reset UIDs
03744     for(unsigned int i=0; i<PossibleJoins[View].size(); ++i) {
03745       if(PossibleJoins[View][i]->GetUID()==4) {PossibleJoins[View][i]->SetUID(3);}
03746     }
03747 
03748     PossibleJoins[View].clear();
03749   }
03751 
03752 
03753 
03754   // Print out the list of 2D tracks
03756   for(int View=0; View<2; ++View) {
03757     MSG("AlgTrackCamList", Msg::kVerbose) << "View: " << View << endl;
03758     for(unsigned i=0; i<ViewSegBank[View].size(); ++i) {
03759       TrackSegmentCam* Seg = ViewSegBank[View][i];
03760       if(Seg->GetUID()>2) {
03761         MSG("AlgTrackCamList", Msg::kVerbose)  << "  i=" << i
03762                                                << "  begpln: "  << Seg->GetBegPlane()
03763                                                << "  begtpos: " << Seg->GetBegTPos()
03764                                                << "  endpln: "  << Seg->GetEndPlane()
03765                                                << "  endtpos: " << Seg->GetEndTPos() << endl;
03766       }
03767     }
03768   }
03770 
03771   return;
03772 }

void AlgTrackCamList::LookForHitsAcrossGap ( TrackCam Trk  ) 

Definition at line 5125 of file AlgTrackCamList.cxx.

References TrackCam::AddHit(), AllHitBank, MuELoss::e, TrackCam::GetBegPlane(), TrackCam::GetBegTPos(), TrackCam::GetBegZPos(), HitCam::GetCharge(), TrackCam::GetEndPlane(), TrackCam::GetEndTPos(), TrackCam::GetEndZPos(), TrackCam::GetEntries(), TrackCam::GetHit(), HitCam::GetPlane(), TrackCam::GetPlaneView(), HitCam::GetTPos(), HitCam::GetZPos(), HitsInTracks, min, and StripWidth.

Referenced by Form3DTracks().

05126 {
05127   // Try and add missing hits that are just across the SM gap
05128   
05129   int PlanesSinceAdded, CurrentPlane(0), HitsAdded(0);
05130   double CurrentZPos, TargetTPos;
05131   double TrkBegZPos, TrkBegTPos, TrkBegDir(0.);
05132   double TrkEndZPos, TrkEndTPos, TrkEndDir(0.);
05133   double Tolerance(0.);
05134 
05135   double MinDistanceToHit, DistanceToCurrentHit;
05136   double centre, extrem1, extrem2;
05137   unsigned int nhits; 
05138 
05139   double z, t, w, sw, swx, swx2, swy, swyx;
05140   bool ConsiderHit;
05141   
05142   int PlaneView=Trk->GetPlaneView();
05143 
05144   // If we want to look back into SM 1
05146   if(Trk->GetBegPlane()>248 && Trk->GetBegPlane()<=256) {
05147               
05148     if(PlaneView==0) {CurrentPlane=248;} else if(PlaneView==1) {CurrentPlane=247;}
05149     PlanesSinceAdded=0;
05150     TrkBegZPos=Trk->GetBegZPos();
05151     TrkBegTPos=Trk->GetBegTPos();
05152  
05153 
05154     // Get beginning direction
05156     z=0.; t=0.; sw=0.; swx=0.; swx2=0.; swy=0.; swyx=0.;
05157      
05158     // Include first few hits in track
05160     for(unsigned int i=0; i<Trk->GetEntries(); ++i) {
05161       HitCam* hit = Trk->GetHit(i);
05162 
05163       if(hit->GetPlane()<Trk->GetBegPlane()+8) {
05164         z=hit->GetZPos(); t=hit->GetTPos(); w=hit->GetCharge();
05165         sw+=w; swx+=w*z; swx2+=w*z*z; swy+=w*t; swyx+=w*t*z;  
05166       }
05167     }
05168     // Use result so far to define a tolerance
05169     double denom=swx*swx-sw*swx2;
05170     if(fabs(denom)>1.e-10) {TrkBegDir = (swx*swy-sw*swyx)/denom; Tolerance = 0.2 + fabs(TrkBegDir);}
05172 
05173               
05174     // Look at some of the relevant hits across the gap
05176     for(int i=CurrentPlane; i>CurrentPlane-8; i-=2) {
05177       for(unsigned int j=0; j<AllHitBank[i].size(); ++j) {
05178         HitCam* hit = AllHitBank[i][j];
05179           
05180         z=hit->GetZPos(); t=hit->GetTPos(); w=hit->GetCharge();
05181 
05182         // If hit is within tolerance region, include with weight 0.5
05183         if(fabs(t-(TrkBegTPos+TrkBegDir*(z-TrkBegZPos)))<Tolerance) {
05184           sw+=w; swx+=w*z; swx2+=w*z*z; swy+=w*t; swyx+=w*t*z;  
05185         }
05186       }
05187     }
05189     // Get results of total fit
05190     denom=swx*swx-sw*swx2;
05191     if(fabs(denom)>1.e-10) {TrkBegDir = (swx*swy-sw*swyx)/denom;}
05193 
05194                       
05195 
05196     // Look to add the hits
05198     while(CurrentPlane>0 && PlanesSinceAdded<10 && HitsAdded<4) {
05199                 
05200       nhits=AllHitBank[CurrentPlane].size();
05201                 
05202       if(nhits>0) {
05203         CurrentZPos=AllHitBank[CurrentPlane][0]->GetZPos();
05204 
05205         HitCam* BestHit = 0;
05206         MinDistanceToHit=3.*StripWidth;
05207 
05208         for(unsigned int j=0; j<nhits; ++j) {
05209           // Don't add a hit that is already in another track in the slice
05210           ConsiderHit=true;
05211           for(unsigned int p=0; p<HitsInTracks[PlaneView].size(); ++p) {
05212             if(AllHitBank[CurrentPlane][j]==HitsInTracks[PlaneView][p]) {ConsiderHit=false; break;}
05213           }
05214           if(ConsiderHit==false) {continue;}
05215 
05216 
05217           // Compare transverse positions
05218           TargetTPos=TrkBegTPos+TrkBegDir*(CurrentZPos-TrkBegZPos);
05219 
05220           centre=AllHitBank[CurrentPlane][j]->GetTPos();
05221           extrem1=centre+(0.055*TrkBegDir);
05222           extrem2=centre-(0.055*TrkBegDir);
05223                     
05224           DistanceToCurrentHit=min(fabs(centre-TargetTPos),
05225                                    min(fabs(extrem1-TargetTPos),fabs(extrem2-TargetTPos)));
05226 
05227           if(DistanceToCurrentHit<MinDistanceToHit) {
05228             MinDistanceToHit=DistanceToCurrentHit; 
05229             BestHit=AllHitBank[CurrentPlane][j]; 
05230             PlanesSinceAdded=0;
05231           }
05232         }
05233 
05234         if(BestHit) {Trk->AddHit(BestHit);}
05235       }
05236       else {PlanesSinceAdded+=2;}
05237       CurrentPlane-=2;
05238     }
05240   }
05242 
05243 
05244 
05245 
05246   // If we want to look forward into SM 2
05248   if(Trk->GetEndPlane()>=242 && Trk->GetEndPlane()<250) {
05249         
05250     if(PlaneView==0) {CurrentPlane=251;} else if(PlaneView==1) {CurrentPlane=250;}
05251     PlanesSinceAdded=0;
05252     TrkEndZPos=Trk->GetEndZPos();
05253     TrkEndTPos=Trk->GetEndTPos();
05254         
05255     
05256     // Get end direction
05258     z=0.; t=0.; sw=0.; swx=0.; swx2=0.; swy=0.; swyx=0.;
05259      
05260     // Include last few hits in track
05262     for(unsigned int i=0; i<Trk->GetEntries(); ++i) {
05263       HitCam* hit = Trk->GetHit(i);
05264 
05265       if(hit->GetPlane()>Trk->GetEndPlane()-8) {
05266         z=hit->GetZPos(); t=hit->GetTPos(); w=hit->GetCharge();
05267         sw+=w; swx+=w*z; swx2+=w*z*z; swy+=w*t; swyx+=w*t*z; 
05268       }
05269     }
05270     // Use result so far to define a tolerance
05271     double denom=swx*swx-sw*swx2;
05272     if(fabs(denom)>1.e-10) {TrkEndDir = (swx*swy-sw*swyx)/denom; Tolerance = 0.2 + fabs(TrkEndDir);}
05274               
05275 
05276     // Look at some of the relevant hits across the gap
05278     for(int i=CurrentPlane; i<CurrentPlane+8; i+=2) {
05279       for(unsigned int j=0; j<AllHitBank[i].size(); ++j) {
05280         HitCam* hit = AllHitBank[i][j];
05281           
05282         z=hit->GetZPos(); t=hit->GetTPos(); w=hit->GetCharge();
05283 
05284         // If hit is within tolerance region, include with weight 0.5
05285         if(fabs(t-(TrkEndTPos+TrkEndDir*(z-TrkEndZPos)))<Tolerance) {
05286           sw+=w; swx+=w*z; swx2+=w*z*z; swy+=w*t; swyx+=w*t*z;  
05287         }
05288       }
05289     }
05291 
05292     // Get results of total fit
05293     denom=swx*swx-sw*swx2;
05294     if(fabs(denom)>1.e-10) {TrkEndDir = (swx*swy-sw*swyx)/denom;}
05296   
05297     
05298 
05299     // Look at hits in region
05301     while(CurrentPlane<500 && PlanesSinceAdded<10 && HitsAdded<4) {
05302 
05303       nhits=AllHitBank[CurrentPlane].size();
05304         
05305       if(nhits>0) {
05306         CurrentZPos=AllHitBank[CurrentPlane][0]->GetZPos();
05307 
05308         HitCam* BestHit = 0;
05309         MinDistanceToHit=3.*StripWidth;
05310                   
05311         for(unsigned int j=0; j<nhits; ++j) {
05312           // Don't add a hit that is already in another track in the slice
05313           ConsiderHit=true;
05314           for(unsigned int p=0; p<HitsInTracks[PlaneView].size(); ++p) {
05315             if(AllHitBank[CurrentPlane][j]==HitsInTracks[PlaneView][p]) {ConsiderHit=false; break;}
05316           }
05317           if(ConsiderHit==false) {continue;}
05318 
05319 
05320           // Compare transverse positions
05321           TargetTPos=TrkEndTPos+TrkEndDir*(CurrentZPos-TrkEndZPos);
05322 
05323           centre=AllHitBank[CurrentPlane][j]->GetTPos();
05324           extrem1=centre+(0.055*TrkEndDir);
05325           extrem2=centre-(0.055*TrkEndDir);
05326                     
05327           DistanceToCurrentHit=min(fabs(centre-TargetTPos),
05328                                    min(fabs(extrem1-TargetTPos),fabs(extrem2-TargetTPos)));
05329 
05330           if(DistanceToCurrentHit<MinDistanceToHit) {
05331             MinDistanceToHit=DistanceToCurrentHit; 
05332             BestHit=AllHitBank[CurrentPlane][j]; 
05333             PlanesSinceAdded=0;
05334           }
05335         }
05336         
05337         if(BestHit) {Trk->AddHit(BestHit);}
05338       }
05339       else {PlanesSinceAdded+=2;}
05340       CurrentPlane+=2;
05341     }
05343   }
05345 
05346   return;
05347 }

void AlgTrackCamList::MatchUV ( TrackCam trku,
TrackCam trkv 
)

Definition at line 5392 of file AlgTrackCamList.cxx.

References TrackCam::AddHit(), AllHitBank, TrackCam::GetBegDir(), TrackCam::GetBegTPos(), TrackCam::GetBegZPos(), OscFit::GetCharge(), TrackCam::GetEndDir(), TrackCam::GetEndTPos(), TrackCam::GetEndZPos(), TrackCam::GetEntries(), TrackCam::GetHit(), HitCam::GetPlane(), HitCam::GetTPos(), HitBank, HitsInTracks, PECut, HitCam::SetUID(), and size.

Referenced by Form3DTracks().

05393 {
05394   // Initial check and declarations
05395   if(PECut<1.) {return;}
05396 
05397   vector<HitCam*> UTrack[500];
05398   vector<HitCam*> VTrack[500];
05399 
05400 
05401   // Store track hits  
05403   for(unsigned int i=0; i<trku->GetEntries(); ++i) {
05404     HitCam* hit = trku->GetHit(i);
05405     UTrack[hit->GetPlane()].push_back(hit);
05406   }
05407 
05408   for(unsigned int i=0; i<trkv->GetEntries(); ++i) {
05409     HitCam* hit = trkv->GetHit(i);
05410     VTrack[hit->GetPlane()].push_back(hit);
05411   }
05413 
05414 
05415   // Find levels of overlap
05417   int MinUPlane=500, MaxUPlane=-20;
05418   int MinUPlaneGTPECut=500, MaxUPlaneGTPECut=-20;
05419 
05420   int MinVPlane=500, MaxVPlane=-20;
05421   int MinVPlaneGTPECut=500, MaxVPlaneGTPECut=-20;
05422 
05423   bool GTPECut=false, ConsiderHit=true, LowPH=false;
05424   double PredictedTPos=0; int LastPlaneAdded;
05425 
05426   for(int i=0; i<500; ++i) {
05427     // For U
05428     if(UTrack[i].size()>0) {
05429       if(i<MinUPlane) {MinUPlane=i;}
05430       if(i>MaxUPlane) {MaxUPlane=i;}
05431 
05432       GTPECut=false;
05433       
05434       for(unsigned int j=0; j<UTrack[i].size(); ++j) {
05435         if(UTrack[i][j]->GetCharge()>PECut) {GTPECut=true; break;}
05436       }
05437 
05438       if(GTPECut) {
05439         if(i<MinUPlaneGTPECut) {MinUPlaneGTPECut=i;}
05440         if(i>MaxUPlaneGTPECut) {MaxUPlaneGTPECut=i;}    
05441       }
05442     }
05443 
05444     // For V
05445     if(VTrack[i].size()>0) {
05446       if(i<MinVPlane) {MinVPlane=i;}
05447       if(i>MaxVPlane) {MaxVPlane=i;}
05448 
05449       GTPECut=false;
05450 
05451       for(unsigned int j=0; j<VTrack[i].size(); ++j) {
05452         if(VTrack[i][j]->GetCharge()>PECut) {GTPECut=true; break;}
05453       }
05454 
05455       if(GTPECut) {
05456         if(i<MinVPlaneGTPECut) {MinVPlaneGTPECut=i;}
05457         if(i>MaxVPlaneGTPECut) {MaxVPlaneGTPECut=i;}
05458       }
05459     }
05460   }
05462 
05463 
05464 
05465 
05466   // Carry out the PH dependent UV match
05467   //
05468   // Two cases:
05469   // 1. Have extended too far in one view, including low PH strips
05470   //
05471   // 2. Have missed hits in one view and need to add them
05472 
05473   // Beginning
05475   if(abs(MinUPlane-MinVPlane)>8) {
05476 
05477     if(MinUPlane<MinVPlane) {
05478       // Removing hits from beginning of u track
05480       if(MinVPlane-MinUPlaneGTPECut<8) {
05481         
05482         for(int i=MinUPlane; i<MinUPlaneGTPECut; ++i) {
05483           for(unsigned int j=0; j<UTrack[i].size(); ++j) {
05484             UTrack[i][j]->SetUID(-1);
05485           }
05486         }
05487       }
05489 
05490 
05491       // Adding hits to beginning of v track
05493       else {
05494         LastPlaneAdded=MinVPlane;
05495 
05496         for(int i=MinVPlane-1; i>=MinUPlane-6; --i) {
05497           if(i<0 || i>120 || LastPlaneAdded-i>30) {break;}
05498 
05499           if((AllHitBank[i].size()==1 || HitBank[i].size()==1) && AllHitBank[i][0]->GetPlaneView()==1) {
05500 
05501             PredictedTPos=trkv->GetBegTPos()+trkv->GetBegDir()*(AllHitBank[i][0]->GetZPos()-trkv->GetBegZPos());
05502 
05503             HitCam* hit = 0;
05504             if(HitBank[i].size()==1) {hit=HitBank[i][0]; LowPH=false;}
05505             else {hit=AllHitBank[i][0]; LowPH=true;}
05506 
05507             // Don't add a hit that is already in another track in the slice
05508             ConsiderHit=true;
05509             for(unsigned int p=0; p<HitsInTracks[1].size(); ++p) {
05510               if(hit==HitsInTracks[1][p]) {ConsiderHit=false; break;}
05511             }
05512 
05513             // Check to see if it is genuinely in a clean +/-5 plane instrumented part of detector
05514             if(ConsiderHit) {
05515               for(int q=-2; q<=2; ++q) {
05516                 if(!LowPH) {for(int q=-1; q<=1; ++q) {if(q!=0 && HitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
05517                 else {for(int q=-1; q<=1; ++q) {if(q!=0 && AllHitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
05518               }
05519             }
05520 
05521             if(ConsiderHit==false) {continue;}
05522 
05523             if(fabs(hit->GetTPos()-PredictedTPos)<0.4) {trkv->AddHit(hit); LastPlaneAdded=i;}
05524           }      
05525         }
05526       }
05528     }
05529 
05530     else {
05531       // Removing hits from beginning of V track
05533       if(MinUPlane-MinVPlaneGTPECut<8) {
05534         
05535         for(int i=MinVPlane; i<MinVPlaneGTPECut; ++i) {
05536           for(unsigned int j=0; j<VTrack[i].size(); ++j) {
05537             VTrack[i][j]->SetUID(-1);
05538           }
05539         }
05540       }
05542 
05543 
05544       // Adding hits to beginning of u track
05546       else {
05547         LastPlaneAdded=MinUPlane;
05548 
05549         for(int i=MinUPlane-1; i>=MinVPlane-6; --i) {
05550           if(i<0 || i>120 || LastPlaneAdded-i>30) {break;}
05551 
05552           if((AllHitBank[i].size()==1 || HitBank[i].size()==1) && AllHitBank[i][0]->GetPlaneView()==0) {
05553 
05554             PredictedTPos=trku->GetBegTPos()+trku->GetBegDir()*(AllHitBank[i][0]->GetZPos()-trku->GetBegZPos());
05555 
05556             HitCam* hit = 0;
05557             if(HitBank[i].size()==1) {hit=HitBank[i][0]; LowPH=false;}
05558             else {hit=AllHitBank[i][0]; LowPH=true;}
05559 
05560             // Don't add a hit that is already in another track in the slice
05561             ConsiderHit=true;
05562             for(unsigned int p=0; p<HitsInTracks[0].size(); ++p) {
05563               if(hit==HitsInTracks[0][p]) {ConsiderHit=false; break;}
05564             }
05565 
05566             // Check to see if it is genuinely in a clean +/-5 plane instrumented part of detector
05567             if(ConsiderHit) {
05568               if(!LowPH) {for(int q=-1; q<=1; ++q) {if(q!=0 && HitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
05569               else {for(int q=-1; q<=1; ++q) {if(q!=0 && AllHitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
05570             }
05571 
05572             if(ConsiderHit==false) {continue;}
05573 
05574             if(fabs(hit->GetTPos()-PredictedTPos)<0.4) {trku->AddHit(hit); LastPlaneAdded=i;}
05575           }
05576         }
05577       }
05579     }
05580   }
05582 
05583 
05584   // End
05586   if(abs(MaxUPlane-MaxVPlane)>8) {
05587 
05588     if(MaxUPlane>MaxVPlane) {
05589       // Removing hits from end of u track
05591       if(MaxUPlaneGTPECut-MaxVPlane<8) {
05592         
05593         for(int i=MaxUPlaneGTPECut+1; i<=MaxUPlane; ++i) {
05594           for(unsigned int j=0; j<UTrack[i].size(); ++j) {
05595             UTrack[i][j]->SetUID(-1);
05596           }
05597         }
05598       }
05600       
05601 
05602       // Adding hits to end of v track
05604       else {
05605         LastPlaneAdded=MaxVPlane;
05606 
05607         for(int i=MaxVPlane; i<MaxUPlane+6; ++i) {
05608           if(i<0 || i>120 || i-LastPlaneAdded>30) {break;}
05609 
05610           if((AllHitBank[i].size()==1 || HitBank[i].size()==1) && AllHitBank[i][0]->GetPlaneView()==1) {
05611 
05612             PredictedTPos=trkv->GetEndTPos()+trkv->GetEndDir()*(AllHitBank[i][0]->GetZPos()-trkv->GetEndZPos());
05613 
05614             HitCam* hit = 0;
05615             if(HitBank[i].size()==1) {hit=HitBank[i][0]; LowPH=false;}
05616             else {hit=AllHitBank[i][0]; LowPH=true;}
05617 
05618             // Don't add a hit that is already in another track in the slice
05619             ConsiderHit=true;
05620             for(unsigned int p=0; p<HitsInTracks[1].size(); ++p) {
05621               if(hit==HitsInTracks[1][p]) {ConsiderHit=false; break;}
05622             }
05623 
05624             // Check to see if it is genuinely in a clean +/-5 plane instrumented part of detector
05625             if(ConsiderHit) {
05626               if(!LowPH) {for(int q=-1; q<=1; ++q) {if(q!=0 && HitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
05627               else {for(int q=-1; q<=1; ++q) {if(q!=0 && AllHitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
05628             }
05629 
05630             if(ConsiderHit==false) {continue;}
05631 
05632             if(fabs(hit->GetTPos()-PredictedTPos)<0.4) {trkv->AddHit(hit); LastPlaneAdded=i;}
05633           }      
05634         }
05635       }
05637     }
05638     
05639     else {
05640       // Removing hits from end of v track
05642       if(MaxVPlaneGTPECut-MaxUPlane<8) {
05643         
05644         for(int i=MaxVPlaneGTPECut+1; i<=MaxVPlane; ++i) {
05645           for(unsigned int j=0; j<VTrack[i].size(); ++j) {
05646             VTrack[i][j]->SetUID(-1);
05647           }
05648         }
05649       }
05651       
05652 
05653       // Adding hits to end of u track
05655       else {
05656         LastPlaneAdded=MaxUPlane;
05657 
05658         for(int i=MaxUPlane; i<MaxVPlane+6; ++i) {
05659           if(i<0 || i>120 || i-LastPlaneAdded>30) {break;}
05660 
05661           if((AllHitBank[i].size()==1 || HitBank[i].size()==1) && AllHitBank[i][0]->GetPlaneView()==0) {
05662 
05663             PredictedTPos=trku->GetEndTPos()+trku->GetEndDir()*(AllHitBank[i][0]->GetZPos()-trku->GetEndZPos());
05664 
05665             HitCam* hit = 0;
05666             if(HitBank[i].size()==1) {hit=HitBank[i][0]; LowPH=false;}
05667             else {hit=AllHitBank[i][0]; LowPH=true;}
05668 
05669             // Don't add a hit that is already in another track in the slice
05670             ConsiderHit=true;
05671             for(unsigned int p=0; p<HitsInTracks[0].size(); ++p) {
05672               if(hit==HitsInTracks[0][p]) {ConsiderHit=false; break;}
05673             }
05674 
05675             // Check to see if it is genuinely in a clean +/-5 plane instrumented part of detector
05676             if(ConsiderHit) {
05677               if(!LowPH) {for(int q=-1; q<=1; ++q) {if(q!=0 && HitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
05678               else {for(int q=-1; q<=1; ++q) {if(q!=0 && AllHitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
05679             }
05680 
05681             if(ConsiderHit==false) {continue;}
05682 
05683             if(fabs(hit->GetTPos()-PredictedTPos)<0.4) {trku->AddHit(hit); LastPlaneAdded=i;}
05684           }      
05685         }
05686       }
05688     }
05689   }
05691 
05692 
05693   // Clear up and return
05694   for(int i=0; i<500; ++i) {UTrack[i].clear(); VTrack[i].clear();}
05695   
05696   return;
05697 }

void AlgTrackCamList::NearDetectorTriplets (  ) 

Definition at line 2228 of file AlgTrackCamList.cxx.

References TrackSegmentCam::AddAssocSegToBeg(), TrackSegmentCam::AddAssocSegToEnd(), TrackSegmentCam::AddMatchSegToBeg(), TrackSegmentCam::AddMatchSegToEnd(), TrackSegmentCam::AddSegment(), ClusterBank, TrackSegmentCam::GetAssocSegBeg(), TrackSegmentCam::GetAssocSegEnd(), TrackSegmentCam::GetBegPlane(), ClusterCam::GetBegTPos(), TrackSegmentCam::GetBegTPos(), TrackSegmentCam::GetCluster(), TrackSegmentCam::GetEndPlane(), ClusterCam::GetEndTPos(), TrackSegmentCam::GetEndTPos(), TrackSegmentCam::GetEntries(), TrackSegmentCam::GetNAssocSegBeg(), TrackSegmentCam::GetNAssocSegEnd(), ClusterCam::GetPlane(), TrackSegmentCam::GetPlaneView(), TrackSegmentCam::GetSeedSegment(), TrackSegmentCam::GetUID(), TrackSegmentCam::IsAssoc(), Munits::km, Msg::kVerbose, MSG, NDSegmentBank, PlanesInModule, ClusterCam::SetNDFlag(), TrackSegmentCam::SetSeedSegment(), TrackSegmentCam::SetUID(), size, and ViewSegBank.

Referenced by RunTheFinder().

02229 {
02230   MSG("AlgTrackCamList", Msg::kVerbose) << " MAKE ND TRIPLETS " << endl;
02231 
02232   int TripAssocNum; bool AlternateTriplets;
02233   vector<TrackSegmentCam*> NDViewSegBank[2];
02234   vector<TrackSegmentCam*> TempNDViewSegBank[2];
02235 
02236 
02237   // MAKE THE TRIPLETS
02239   for(int i=10; i<PlanesInModule-9; ++i) {
02240 
02241     // Only consider Fully Instrumented planes
02242     if((i-1)%10!=0 && (i-6)%10!=0) {continue;}
02243 
02244     // Loop over clusters on plane. This is our plane 0, the origin.
02245     for(unsigned int k0=0; k0<ClusterBank[i].size(); ++k0) {
02246     
02247       if(ClusterBank[i][k0]->GetTrkFlag()>0 && ClusterBank[i][k0]->GetNDFlag()>0) {
02248         
02249         // Look for triplets of form m 0 p
02251         if( (i-10)>=0 && ClusterBank[i-10].size()>0 && (i+10)<=PlanesInModule && ClusterBank[i+10].size()>0 ) {
02252           
02253           // Look at the clusters available on plane i-10...
02254           for(unsigned int km=0; km<ClusterBank[i-10].size(); ++km) {
02255             
02256             // ...and the clusters available on plane i+10
02257             for(unsigned int kp=0; kp<ClusterBank[i+10].size(); ++kp) {
02258               
02259               if( ClusterBank[i-10][km]->GetTrkFlag()>0 && ClusterBank[i-10][km]->GetNDFlag()>0
02260                   && ClusterBank[i+10][kp]->GetTrkFlag()>0 && ClusterBank[i+10][kp]->GetNDFlag()>0 ) {
02261                 
02262                 // Check the track-like association of the three clusters
02263                 TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][km],ClusterBank[i+10][kp],10);
02264         
02265                 // If the association is good, make the triplet
02266                 if(TripAssocNum>0) {
02267                   TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-10][km], ClusterBank[i][k0], ClusterBank[i+10][kp]);
02268                   for(unsigned int j=0; j<seg0->GetEntries(); ++j) {seg0->GetCluster(j)->SetNDFlag(2);}
02269 
02270                   NDSegmentBank[i].push_back(seg0);
02271                 }
02272               }
02273             }
02274           }
02275         }
02277 
02278 
02279 
02280         
02281         // Look for triplets of form m X 0 p  
02283         if( (i-20)>=0 && ClusterBank[i-20].size()>0 && (i+10)<=PlanesInModule && ClusterBank[i+10].size()>0 ) {
02284 
02285           // Look at the clusters on plane i+10...
02286           for(unsigned int kp=0; kp<ClusterBank[i+10].size(); ++kp) {
02287               
02288             // ...and the clusters on plane i-20
02289             for(unsigned int km=0; km<ClusterBank[i-20].size(); ++km) {
02290 
02291               if( ClusterBank[i-20][km]->GetTrkFlag()>0 && ClusterBank[i+10][kp]->GetTrkFlag()>0
02292                   && ClusterBank[i-20][km]->GetNDFlag()>0 && ClusterBank[i+10][kp]->GetNDFlag()>0 ) {
02293 
02294                 // Check the track-like association of the three clusters
02295                 TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i+10][kp],10);
02296                   
02297                 // If the association is good, check whether we can make the triplet
02298                 if(TripAssocNum>0) {
02299                   AlternateTriplets=false;
02300                       
02301                   // Check to see if there is also an alternate triplet possibility with clusters on plane i-2
02302                   // i.e. Want to make sure that the X in "m X 0 p" is correct.
02303 
02304                   if(AlternateTriplets==false && ClusterBank[i-10].size()>0) {
02305                     // Look at the clusters on plane i-10
02306                     for(unsigned int ktmp=0; ktmp<ClusterBank[i-10].size(); ++ktmp) {
02307                       if(AlternateTriplets==false && ClusterBank[i-10][ktmp]->GetTrkFlag()>0
02308                          && ClusterBank[i-10][ktmp]->GetNDFlag()>0
02309                          && ClusterBank[i-10][ktmp]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i][k0],10)>0 
02310                          && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][ktmp],ClusterBank[i+10][kp],10)>0) 
02311                         {AlternateTriplets=true;}
02312                     }
02313                   }
02314                                       
02315                   // If everything is good, make the triplet
02316                   if(AlternateTriplets==false) {
02317                     TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-20][km],ClusterBank[i][k0],ClusterBank[i+10][kp]);
02318                     for(unsigned int j=0; j<seg0->GetEntries(); ++j) {seg0->GetCluster(j)->SetNDFlag(2);}
02319 
02320                     NDSegmentBank[i].push_back(seg0);
02321                   }
02322                 }
02323               }
02324             }
02325           }
02326         }
02328 
02329 
02330 
02331 
02332         // Look for triplets of form m 0 X p  
02334         if( (i-10)>=0 && ClusterBank[i-10].size()>0 && (i+20)<=PlanesInModule && ClusterBank[i+20].size()>0 ) {
02335 
02336           // Look at the clusters on plane i-10...
02337           for(unsigned int km=0; km<ClusterBank[i-10].size(); ++km) {
02338               
02339             // ... and the clusters on plane i+20
02340             for(unsigned int kp=0; kp<ClusterBank[i+20].size(); ++kp) {
02341 
02342               if( ClusterBank[i-10][km]->GetTrkFlag()>0 && ClusterBank[i+20][kp]->GetTrkFlag()>0
02343                   && ClusterBank[i-10][km]->GetNDFlag()>0 && ClusterBank[i+20][kp]->GetNDFlag()>0 ) {
02344 
02345                 // Check the track-like association of the three clusters
02346                 TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][km],ClusterBank[i+20][kp],10);
02347                    
02348                 // If the association is good, check whether we can make the triplet
02349                 if(TripAssocNum>0) {
02350                   AlternateTriplets=false;
02351                       
02352                   // Check to see if there is also an alternate triplet possibility with clusters on plane i+2
02353                   // i.e. Want to make sure that the X in "m 0 X p" is correct.
02354 
02355                   if(AlternateTriplets==false && ClusterBank[i+10].size()>0) {
02356                     // Look at the clusters on plane i+10
02357                     for(unsigned int ktmp=0; ktmp<ClusterBank[i+10].size(); ++ktmp) {
02358                       if(AlternateTriplets==false && ClusterBank[i+10][ktmp]->GetTrkFlag()>0 
02359                          && ClusterBank[i+10][ktmp]->GetNDFlag()>0
02360                          && ClusterBank[i+10][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+20][kp],10)>0 
02361                          && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][km],ClusterBank[i+10][ktmp],10)>0) 
02362                         {AlternateTriplets=true;}
02363                     }
02364                   }
02365                                       
02366                   // If everything is good, make the triplet
02367                   if(AlternateTriplets==false) {
02368                     // Store segment on empty plane too
02369                     for(int j=0; j<2; ++j) {
02370                       TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-10][km],ClusterBank[i][k0],ClusterBank[i+20][kp]);
02371                       for(unsigned int k=0; k<seg0->GetEntries(); ++k) {seg0->GetCluster(k)->SetNDFlag(2);}
02372   
02373                       NDSegmentBank[i+10*j].push_back(seg0);
02374                     }
02375                   }
02376                 }
02377               }
02378             }
02379           }
02380         }
02382 
02383 
02384 
02385 
02386         // Look for triplets of form m X 0 X p
02388         if( (i-20)>=0 && ClusterBank[i-20].size()>0 && (i+20)<=PlanesInModule && ClusterBank[i+20].size()>0 ) {
02389 
02390           // Look at clusters on planes i-20 and i+20
02391           for(unsigned int km=0; km<ClusterBank[i-20].size(); ++km) {
02392             for(unsigned int kp=0; kp<ClusterBank[i+20].size(); ++kp) {
02393 
02394               if( ClusterBank[i-20][km]->GetTrkFlag()>0 && ClusterBank[i+20][kp]->GetTrkFlag()>0
02395                   && ClusterBank[i-20][km]->GetNDFlag()>0 && ClusterBank[i+20][kp]->GetNDFlag()>0 ) {
02396 
02397                 // Check the track-like association of the three clusters
02398                 TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i+20][kp],10);
02399                   
02400                 // If the association is good, check whether we can make the triplet
02401                 if(TripAssocNum>0) {
02402                   AlternateTriplets=false;
02403                     
02404                   // Check to see if there are alternate triplet possibilities with clusters on plane i-10 or i+10
02405                   // i.e. Want to make sure that the Xs in "m X 0 X p" are correct.
02406                     
02407                   // Check first X, plane i-10
02408                   if(AlternateTriplets==false && ClusterBank[i-10].size()>0) {
02409                     // Look at any clusters on plane i-10
02410                     for(unsigned int ktmp=0; ktmp<ClusterBank[i-10].size(); ++ktmp) {
02411                       if(AlternateTriplets==false && ClusterBank[i-10][ktmp]->GetTrkFlag()>0 
02412                          && ClusterBank[i-10][ktmp]->GetNDFlag()>0
02413                          && ClusterBank[i-10][ktmp]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i][k0],10)>0 
02414                          && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][ktmp],ClusterBank[i+20][kp],10)>0) 
02415                         {AlternateTriplets=true;}
02416                     }
02417                   }
02418                     
02419                   // Check second X, plane i+10
02420                   if(AlternateTriplets==false && ClusterBank[i+10].size()>0) {
02421                     // Look at any clusters on plane i+10
02422                     for(unsigned int ktmp=0; ktmp<ClusterBank[i+10].size(); ++ktmp) {
02423                       if(AlternateTriplets==false  && ClusterBank[i+10][ktmp]->GetTrkFlag()>0 
02424                          && ClusterBank[i+10][ktmp]->GetNDFlag()>0
02425                          && ClusterBank[i+10][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+20][kp],10)>0 
02426                          && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i+10][ktmp],10)>0) 
02427                         {AlternateTriplets=true;}
02428                     }
02429                   }
02430                     
02431                   // Check both Xs together, planes i-10 and i+10
02432                   if(AlternateTriplets==false && ClusterBank[i-10].size()>0 
02433                      && ClusterBank[i+10].size()>0) {
02434                       
02435                     // Loop again at any clusters on plane i-10
02436                     for(unsigned int ktmp=0; ktmp<ClusterBank[i-10].size(); ++ktmp) {
02437                         
02438                       // Look again at any clusters on plane i+10
02439                       for(unsigned int ktmp1=0; ktmp1<ClusterBank[i+10].size(); ++ktmp1) {
02440                         if(AlternateTriplets==false 
02441                            && ClusterBank[i-10][ktmp]->GetTrkFlag()>0 && ClusterBank[i+10][ktmp1]->GetTrkFlag()>0
02442                            && ClusterBank[i-10][ktmp]->GetNDFlag()>0  && ClusterBank[i+10][ktmp1]->GetNDFlag()>0 
02443                            && ClusterBank[i-10][ktmp]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i][k0],10)>0 
02444                            && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][ktmp],ClusterBank[i+10][ktmp1],10)>0 
02445                            && ClusterBank[i+10][ktmp1]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+20][kp],10)>0)
02446                           {AlternateTriplets=true;}
02447                       }
02448                     }
02449                   }
02450 
02451                   // If everything is good, make the triplet
02452                   if(AlternateTriplets==false) {
02453                     // Store segment on empty planes too
02454                     for(int j=0; j<2; ++j) {
02455                       TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-20][km],ClusterBank[i][k0],ClusterBank[i+20][kp]);
02456                       for(unsigned int k=0; k<seg0->GetEntries(); ++k) {seg0->GetCluster(k)->SetNDFlag(2);}
02457         
02458                       NDSegmentBank[i+10*j].push_back(seg0);
02459                     }
02460                   }
02461                 }
02462               }
02463             }
02464           }
02465         }
02467 
02468       }
02469     }
02470   }
02472 
02473 
02474 
02475   // Print out list of triplets
02477   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF ND TRIPLETS *** " << endl;
02478   for(int i=0; i<500; ++i) {
02479     for(unsigned int j=0;j<NDSegmentBank[i].size(); ++j) {
02480       MSG("AlgTrackCamList", Msg::kVerbose) << " Segment Num " << j << ", Plane " << i <<endl;
02481       
02482       ClusterCam* clr0 = NDSegmentBank[i][j]->GetCluster(0);
02483       ClusterCam* clr1 = NDSegmentBank[i][j]->GetCluster(1);
02484       ClusterCam* clr2 = NDSegmentBank[i][j]->GetCluster(2);
02485       
02486       MSG("AlgTrackCamList", Msg::kVerbose) << " (" << clr0->GetPlane() << "," 
02487                                             << clr0->GetBegTPos() << "->" << clr0->GetEndTPos() << ") "
02488                                             << " (" << clr1->GetPlane() << "," << clr1->GetBegTPos() 
02489                                             << "->" << clr1->GetEndTPos() << ") "
02490                                             << " (" << clr2->GetPlane() << "," << clr2->GetBegTPos() 
02491                                             << "->" << clr2->GetEndTPos() << ") " << endl;
02492     }
02493   }
02495 
02496 
02497 
02498 
02499   // Now loop over these segments and check associations
02501   for(int i=10; i<PlanesInModule-9; ++i) {
02502     
02503     if((i+10)<=PlanesInModule && NDSegmentBank[i].size()>0 && NDSegmentBank[i+10].size()>0) {
02504         
02505       // Look at the segments on plane i
02506       for(unsigned int k0=0; k0<NDSegmentBank[i].size(); ++k0) {
02507         TrackSegmentCam* Seg0 = NDSegmentBank[i][k0];
02508 
02509         // For each of these, look at the segments on plane i+10
02510         for(unsigned int kp=0; kp<NDSegmentBank[i+10].size(); ++kp) {
02511           TrackSegmentCam* NextSeg = NDSegmentBank[i+10][kp];
02512 
02513           // Check associations. 
02514           if(Seg0->IsAssoc(NextSeg)) {
02515               
02516             // If segments are associated, add to list of all associations
02517             Seg0->AddAssocSegToEnd(NextSeg);
02518             NextSeg->AddAssocSegToBeg(Seg0);
02519           }
02520         }
02521       }
02522     }
02523   }
02525   
02526 
02527 
02528   // Print out list of all associations
02530   MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF ALL ND TRIPLET ASSOCIATIONS *** " << endl;
02531 
02532   // Print out list of triplets
02533   for(int i=0; i<500; ++i) {
02534     for(unsigned int j=0;j<NDSegmentBank[i].size(); ++j) {
02535       TrackSegmentCam* segment = NDSegmentBank[i][j];
02536 
02537       MSG("AlgTrackCamList", Msg::kVerbose) << "--- Plane " << i << ", Seg Number " << j 
02538                                             << ", BegPlane " << segment->GetBegPlane() 
02539                                             << ", EndPlane " << segment->GetEndPlane() 
02540                                             << ", BegTPos " << segment->GetBegTPos() 
02541                                             << ", EndTPos " << segment->GetEndTPos() 
02542                                             <<endl;
02543       
02544       MSG("AlgTrackCamList", Msg::kVerbose) << " Beg: " <<endl;
02545       for(unsigned int k=0; k<segment->GetNAssocSegBeg(); ++k) {
02546         TrackSegmentCam* segtmp = segment->GetAssocSegBeg(k);
02547         MSG("AlgTrackCamList", Msg::kVerbose) << "  Assoc number=" << k
02548                                               << "  begpln: " << segtmp->GetBegPlane()
02549                                               << "  begtpos: " << segtmp->GetBegTPos()
02550                                               << "  endpln: " << segtmp->GetEndPlane()
02551                                               << "  endtpos: " << segtmp->GetEndTPos() << endl;
02552       }
02553       
02554       MSG("AlgTrackCamList", Msg::kVerbose) << " End: " <<endl;
02555       for(unsigned int k=0; k<segment->GetNAssocSegEnd(); ++k) {
02556         TrackSegmentCam* segtmp = segment->GetAssocSegEnd(k);
02557         MSG("AlgTrackCamList", Msg::kVerbose) << "  Assoc number=" << k
02558                                               << "  begpln: " << segtmp->GetBegPlane()
02559                                               << "  begtpos: " << segtmp->GetBegTPos()
02560                                               << "  endpln: " << segtmp->GetEndPlane()
02561                                               << "  endtpos: " << segtmp->GetEndTPos() << endl;           
02562       }
02563     }
02564   }
02566 
02567 
02568 
02569   // Look for chains of associations and make ND 2D tracks
02571   bool JoinFlag;
02572   
02573   for(int i=10; i<PlanesInModule-9; ++i) {
02574     // Loop over the triplets centered on plane i
02575     for(unsigned int j=0; j<NDSegmentBank[i].size(); ++j) {
02576 
02577       TrackSegmentCam* Seg = NDSegmentBank[i][j];
02578 
02579 
02580       // Firstly, consider isolated segments. These must be their own seed segments.
02581       if(Seg->GetNAssocSegBeg()==0 && Seg->GetNAssocSegEnd()==0) {
02582         NDViewSegBank[Seg->GetPlaneView()].push_back(Seg);
02583         Seg->SetSeedSegment(Seg);
02584       }
02585 
02586 
02587       // Then, consider those segments which have some associations
02588       if(Seg->GetNAssocSegBeg()>0 || Seg->GetNAssocSegEnd()>0) {
02589         JoinFlag=false;
02590 
02591         // If there is one assoc at beginning (always false for first segment in module)
02592         if(Seg->GetNAssocSegBeg()==1) {
02593           TrackSegmentCam* Segm = Seg->GetAssocSegBeg(0);
02594 
02595           // Check that we are considering a simple chain of segments
02596           if(Segm->GetNAssocSegEnd()==1) {
02597             // Get the last segment in the chain of segments
02598             TrackSegmentCam* SeedSeg = Segm->GetSeedSegment();
02599 
02600             SeedSeg->AddSegment(Seg); 
02601             Seg->SetSeedSegment(SeedSeg);
02602             Seg->SetUID(-1);
02603             JoinFlag=true;
02604           }
02605         }
02606 
02607         // Now consider those triplets which aren't just one assoc at beg and one at end
02608         if(JoinFlag==false) {
02609           Seg->SetSeedSegment(Seg);
02610           NDViewSegBank[Seg->GetPlaneView()].push_back(Seg);
02611                   
02612           // Loop over the associations at beginning
02613           for(unsigned int k=0; k<Seg->GetNAssocSegBeg(); ++k) {
02614             TrackSegmentCam* SeedSeg = Seg->GetAssocSegBeg(k)->GetSeedSegment();
02615             
02616             Seg->AddMatchSegToBeg(SeedSeg);
02617             SeedSeg->AddMatchSegToEnd(Seg);
02618           }
02619         }
02620       }
02621 
02622     } 
02623   }
02625 
02626 
02627 
02628 
02629   // Print out the list of 2D tracks
02631   MSG("AlgTrackCamList", Msg::kVerbose)  << " *** LIST OF ND 2D TRACKS *** " << endl;
02632   for(int View=0; View<2; ++View) {
02633 
02634     MSG("AlgTrackCamList", Msg::kVerbose) << "View: " << View << endl;
02635     for(unsigned i=0; i<NDViewSegBank[View].size(); ++i) {
02636       TrackSegmentCam* Seg = NDViewSegBank[View][i];
02637       MSG("AlgTrackCamList", Msg::kVerbose) << "  i=" << i
02638                                             << "  begpln: "  << Seg->GetBegPlane()
02639                                             << "  begtpos: " << Seg->GetBegTPos()
02640                                             << "  endpln: "  << Seg->GetEndPlane()
02641                                             << "  endtpos: " << Seg->GetEndTPos() << endl;
02642     }
02643   }
02645 
02646 
02647 
02648 
02649   // Now match to existing 2D tracks
02650   // If no match, store separately in ViewSegBank
02652   bool Match;
02653 
02654   // For ND U segments
02656   for(unsigned int i=0; i<NDViewSegBank[0].size(); ++i) {
02657     TrackSegmentCam* NDSegu = NDViewSegBank[0][i];
02658     Match=false;
02659       
02660     // Consider the other U segments
02661     for(unsigned int j=0; j<ViewSegBank[0].size(); ++j) {
02662       TrackSegmentCam* Segu = ViewSegBank[0][j];
02663 
02664       if(Segu->GetUID()>1) {
02665 
02666         if(NDSegu->GetEndPlane()<Segu->GetEndPlane()) {Match=NDSegu->IsAssoc(Segu);}
02667         else {Match=Segu->IsAssoc(NDSegu);}
02668 
02669         if(Match) {
02670           if(NDSegu->GetBegPlane()<Segu->GetBegPlane()) {Segu->AddMatchSegToBeg(NDSegu);}
02671           else {Segu->AddMatchSegToEnd(NDSegu);}
02672 
02673           MSG("AlgTrackCamList", Msg::kVerbose) << "U Match, NDbegplane " << NDSegu->GetBegPlane() 
02674                                                 << " Endplane " << NDSegu->GetEndPlane()
02675                                                 << " Begtpos " << NDSegu->GetBegTPos() 
02676                                                 << " Endtpos " << NDSegu->GetEndTPos()
02677                                                 << " 2d track begplane " << Segu->GetBegPlane() 
02678                                                 << " Endplane " << Segu->GetEndPlane()
02679                                                 << " Begtpos " << Segu->GetBegTPos() 
02680                                                 << " Endtpos " << Segu->GetEndTPos() << endl;
02681         }
02682       }
02683     }
02684 
02685     if(Match) {NDSegu->SetUID(2);}
02686     else {NDSegu->SetUID(2); TempNDViewSegBank[0].push_back(NDSegu);}
02687   }
02689 
02690 
02691   // For ND V segments
02693   for(unsigned int i=0; i<NDViewSegBank[1].size(); ++i) {
02694     TrackSegmentCam* NDSegv = NDViewSegBank[1][i];
02695     Match=false;
02696         
02697     // Consider the other V segments
02698     for(unsigned int j=0; j<ViewSegBank[1].size(); ++j) {
02699       TrackSegmentCam* Segv = ViewSegBank[1][j];
02700 
02701       if(Segv->GetUID()>1) {
02702 
02703         if(NDSegv->GetEndPlane()<Segv->GetEndPlane()) {Match=NDSegv->IsAssoc(Segv);}
02704         else {Match=Segv->IsAssoc(NDSegv);}
02705 
02706         if(Match) {
02707           if(NDSegv->GetBegPlane()<Segv->GetBegPlane()) {Segv->AddMatchSegToBeg(NDSegv);}
02708           else {Segv->AddMatchSegToEnd(NDSegv);}
02709 
02710           MSG("AlgTrackCamList", Msg::kVerbose) << "V Match, NDbegplane " << NDSegv->GetBegPlane() 
02711                                                 << " Endplane " << NDSegv->GetEndPlane()
02712                                                 << " Begtpos " << NDSegv->GetBegTPos() 
02713                                                 << " Endtpos " << NDSegv->GetEndTPos()
02714                                                 << " 2d track begplane " << Segv->GetBegPlane() 
02715                                                 << " Endplane " << Segv->GetEndPlane()
02716                                                 << " Begtpos " << Segv->GetBegTPos() 
02717                                                 << " Endtpos " << Segv->GetEndTPos() << endl;
02718         }
02719       }
02720     }
02721 
02722     if(Match) {NDSegv->SetUID(2);}
02723     else {NDSegv->SetUID(2); TempNDViewSegBank[1].push_back(NDSegv);}
02724   }
02725 
02726 
02727   // Now store those segments marked above in ViewSegBank
02728   for(int View=0; View<2; ++View) {
02729     for(unsigned int i=0; i<TempNDViewSegBank[View].size(); ++i) {
02730       ViewSegBank[View].push_back(TempNDViewSegBank[View][i]);
02731     }
02732   }
02734 
02735 
02736   // Clean up
02737   for(int i=0; i<2; ++i) {NDViewSegBank[i].clear(); TempNDViewSegBank[i].clear();}
02738 
02739   return;
02740 }

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

Implements AlgBase.

Definition at line 67 of file AlgTrackCamList.cxx.

References CandHandle::AddDaughterLink(), CambridgeAnalysis, ClearUp(), FinalTrackBank, Registry::Get(), AlgFactory::GetAlgHandle(), CandContext::GetCandRecord(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), AlgFactory::GetInstance(), Registry::GetInt(), CandContext::GetMom(), CandHandle::GetNDaughters(), RecMinos::GetVldContext(), Detector::kCalDet, Msg::kError, Registry::KeyExists(), Detector::kFar, Detector::kNear, Msg::kWarning, CandTrackCam::MakeCandidate(), ModuleType, MSG, NumModules, PECut, PECut2, PlanesInModule, RunTheFinder(), CandContext::SetCandRecord(), CandContext::SetDataIn(), CandHandle::SetName(), CandHandle::SetTitle(), size, and vldc.

00068 {
00069   assert(cx.GetDataIn());
00070 
00071   // Check for CandSliceListHandle input
00072   if (cx.GetDataIn()->InheritsFrom("CandSliceListHandle")) {
00073     const CandSliceListHandle *cslh = dynamic_cast<const CandSliceListHandle*>(cx.GetDataIn());
00074     const MomNavigator *mom = cx.GetMom();
00075     
00076     CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00077     assert(candrec);
00078     
00079     vldc = (VldContext*)(candrec->GetVldContext());
00080     
00081     
00082     // Create the new tracklist
00083     CandTrackCamListHandle& tracklist = dynamic_cast<CandTrackCamListHandle&>(ch);
00084     
00085     if( !cslh || cslh->GetNDaughters()<1 ) {
00086       // Require number of CandSlices to be non-zero to do anything more
00087       MSG("AlgTrackCamList", Msg::kWarning) << " !cslh || cslh->GetNDaughters()<1 " << endl;
00088       return;
00089     }
00090  
00091     // Configure for correct detector
00092     switch(vldc->GetDetector()){
00093     case Detector::kFar:
00094       ModuleType=1; NumModules=2; PlanesInModule=248; break;
00095     case Detector::kNear:
00096       ModuleType=2; NumModules=1; PlanesInModule=120; break;
00097     case Detector::kCalDet:
00098       ModuleType=3; NumModules=1; PlanesInModule=60; break;
00099     default:
00100       ModuleType=-1; NumModules=0; PlanesInModule=0; break;
00101     }
00102      
00103  
00104     // Unpack AlgConfig
00105     CambridgeAnalysis=false;
00106     if( ac.KeyExists("CambridgeAnalysis") ) {CambridgeAnalysis = ac.GetInt("CambridgeAnalysis");}
00107   
00108  
00109     // Iterate over the slices
00110     TIter sliceItr(cslh->GetDaughterIterator());
00111      
00112     while (CandSliceHandle *slice = dynamic_cast<CandSliceHandle*>(sliceItr())) {
00113 
00114       // Find tracks in slice
00116       if(ModuleType==1 || ModuleType==3) {PECut=1.8; PECut2=1.8;}
00117       else if(ModuleType==2) {PECut=2.; PECut2=0.5;}
00118       RunTheFinder(slice);
00119 
00120         
00121       // Setup AlgTrackCam
00123       AlgFactory &af = AlgFactory::GetInstance();
00124       TString algtrackconfig = "default";
00125       if( ac.KeyExists("AlgTrackCamConfig") ) {
00126         const char* tmps;
00127         if( ac.Get("AlgTrackCamConfig",tmps)) { algtrackconfig = tmps; }
00128       }
00129       AlgHandle ah_trk = af.GetAlgHandle("AlgTrackCam", algtrackconfig.Data());
00130 
00131         
00132       // Create the CandContext
00134       CandContext cx0(this, mom);
00135       cx0.SetCandRecord(candrec);
00136         
00137       unsigned int NTracksFound = FinalTrackBank[0].size();
00138 
00139       // Run the finder again with different PECut settings
00140       if(CambridgeAnalysis==false && NTracksFound==0) {
00141         ClearUp(); PECut=0.5; PECut2=0.5;
00142         RunTheFinder(slice); 
00143         NTracksFound = FinalTrackBank[0].size();
00144       }
00145   
00146       // Sanity Check
00147       if(NTracksFound!= FinalTrackBank[1].size()) {
00148         MSG("AlgTrackCamList", Msg::kError) <<"Different Number of U and V Tracks Found!"<<endl;
00149         if( NTracksFound>=FinalTrackBank[1].size()) {NTracksFound = FinalTrackBank[1].size();}
00150       }
00151 
00152       //Loop over the U and V track found
00154       for(unsigned int i=0; i<NTracksFound; ++i) {
00155         TrackCam* tracku = FinalTrackBank[0][i];
00156         TrackCam* trackv = FinalTrackBank[1][i];
00157         
00158         //Fill the CandContext for this track
00160         TObjArray trackin;
00161         trackin.Add(tracku);
00162         trackin.Add(trackv);
00163         trackin.Add(slice);
00164         cx0.SetDataIn(&trackin);
00165           
00166         //Create the CandTrack
00167         CandTrackCamHandle cth = CandTrackCam::MakeCandidate(ah_trk, cx0);
00168         cth.SetName(TString("CandTrackCamHandle"));
00169         cth.SetTitle(TString("Created by AlgTrackCamList"));
00170           
00171         // Add CandTrack to CandTrackList
00173         if(cth.GetNDaughters()>0) {tracklist.AddDaughterLink(cth);}
00174       }
00175       
00176       ClearUp();
00177 
00178     } // End loop over slices
00179       
00180   } // End if inherits from CandSliceListHandle
00181   
00182 }

void AlgTrackCamList::RunTheFinder ( CandSliceHandle slice  ) 

Definition at line 189 of file AlgTrackCamList.cxx.

References FindAllAssociations(), FindMatchedJoins(), FindPreferredJoins(), FirstUVComparison(), Form2DTracks(), Form3DTracks(), FormTheClusters(), FormTheHits(), FormTriplets(), IDTrkAndShwClusters(), Join2DTracks(), Msg::kVerbose, ModuleType, MSG, NearDetectorTriplets(), and SecondUVComparison().

Referenced by RunAlg().

00190 {
00191   // Configure algorithm for the relevant detector
00192   // For ND, initial algorithm is applied only to calorimeter,
00193   // spectrometer is treated separately
00194   MSG("AlgTrackCamList", Msg::kVerbose) << "FINDING TRACKS IN SLICE" << endl;
00195 
00196   // Run the methods
00197   FormTheHits(slice);
00198   FormTheClusters(); 
00199   IDTrkAndShwClusters();
00200   FormTriplets();
00201   FindAllAssociations();
00202   FindPreferredJoins();
00203   FindMatchedJoins();
00204   FirstUVComparison();
00205   if(ModuleType==2) {NearDetectorTriplets();}
00206   Form2DTracks();
00207   Join2DTracks();
00208   SecondUVComparison();
00209   Form3DTracks(slice);
00210 
00211   return;
00212 }

void AlgTrackCamList::SecondUVComparison (  ) 

Definition at line 3779 of file AlgTrackCamList.cxx.

References TrackSegmentCam::AddSegment(), TrackSegmentCam::GetBegPlane(), TrackSegmentCam::GetBegTime(), TrackSegmentCam::GetBegTPos(), TrackSegmentCam::GetCluster(), TrackSegmentCam::GetEndPlane(), TrackSegmentCam::GetEndTime(), TrackSegmentCam::GetEndTPos(), TrackSegmentCam::GetEntries(), TrackSegmentCam::GetPartner(), TrackSegmentCam::GetUID(), id, TrackSegmentCam::IsAssoc(), Msg::kVerbose, ModuleType, MSG, NumModules, PlanesInModule, TrackSegmentCam::SetPartner(), ClusterCam::SetTrkFlag(), TrackSegmentCam::SetUID(), size, TempTrack, and ViewSegBank.

Referenced by RunTheFinder().

03780 {
03781   // Look at the 2D tracks formed and try to match 2D tracks in the U view
03782   // with those in the V view, based on a rough overlap in plane numbers and time.
03783   // Store those that match.
03784 
03785   // 2D tracks are only ever formed within a single SuperModule. At this point
03786   // we try to match FD 2D tracks across the SM gap. Using the 2D tracks we have
03787   // just stored, we check association across the gap. Any final ambiguity is
03788   // resolved by getting a score for each possible match. 
03789 
03790   // If things look good, U/V partners are set and 2D tracks are joined across the gap.
03791 
03792   // Any other possible 2D track matches between views are then considered and the
03793   // best matches chosen by obtaining a score characterising the degree of overlap.
03794 
03795 
03796 
03797   vector<TrackSegmentCam*> TempTrackSeg[2];
03798   int Module; bool Cont;
03799 
03800   int id, idm, idp; 
03801   double TopScore, Score, NClusters, TopClusters;
03802   int begplane1(0), begplane2(0), endplane1(0), endplane2(0), dplane1, dplane2;
03803   int NumUTrks(0), NumVTrks(0);
03804   
03805   // Store the 2D track segments which roughly overlap in terms of planes
03806   // and time.
03808  
03809   // Loop over all u segments
03810   for(unsigned int i=0; i<ViewSegBank[0].size(); ++i) {
03811     TrackSegmentCam* Segu = ViewSegBank[0][i];
03812  
03813     // If segment is a 2D track   
03814     if(Segu->GetUID()>2) {
03815 
03816       // For each of these, loop over 2D tracks in V view
03817       for(unsigned int j=0; j<ViewSegBank[1].size(); ++j) {
03818         TrackSegmentCam* Segv = ViewSegBank[1][j];
03819         
03820         if(Segv->GetUID()>2) {
03821         
03822           // Check overlap between u and v segments
03823           if( (abs(Segv->GetBegPlane()-Segu->GetBegPlane())<10 || abs(Segv->GetEndPlane()-Segu->GetEndPlane())<10 )
03824               && ( (Segv->GetEndPlane()-Segu->GetBegPlane())>2 && (Segu->GetEndPlane()-Segv->GetBegPlane())>2 )
03825               && ( (Segv->GetEndTime()-Segu->GetBegTime())>-100 && (Segu->GetEndTime()-Segv->GetBegTime())>-100 ) ) {
03826             
03827             TempTrackSeg[0].push_back(Segu);
03828             TempTrackSeg[1].push_back(Segv);
03829           }
03830         }
03831       }
03832     }
03833   }
03834 
03835   // Track combinations we may have missed above:
03837   // Simple combination
03838   TrackSegmentCam* USeg=0; TrackSegmentCam* VSeg=0;
03839 
03840   for(unsigned int i=0; i<ViewSegBank[0].size(); ++i) {if(ViewSegBank[0][i]->GetUID()>2) {NumUTrks++; USeg=ViewSegBank[0][i];}}
03841   for(unsigned int i=0; i<ViewSegBank[1].size(); ++i) {if(ViewSegBank[1][i]->GetUID()>2) {NumVTrks++; VSeg=ViewSegBank[1][i];}}
03842 
03843   if(TempTrackSeg[0].size()==0 && TempTrackSeg[1].size()==0 
03844      && NumUTrks==1 && NumVTrks==1 && USeg && VSeg) {
03845     TempTrackSeg[0].push_back(USeg); 
03846     TempTrackSeg[1].push_back(VSeg);
03847   }
03849 
03850 
03852   // Essentially same as above, but if we have also found a junky second track in one view
03853   if((ModuleType==1 || ModuleType==3) && TempTrackSeg[0].size()==0 && TempTrackSeg[1].size()==0) {
03854     // First combination
03855     if(NumUTrks==1 && NumVTrks==2) {
03856       TrackSegmentCam* VSeg2=0;
03857 
03858       for(unsigned int i=0; i<ViewSegBank[1].size(); ++i) {
03859         if(ViewSegBank[1][i]->GetUID()>2) {VSeg2=ViewSegBank[1][i]; break;}
03860       }
03861       
03862       if(USeg && VSeg && VSeg2) {
03863         if(abs(int(USeg->GetEntries())-int(VSeg->GetEntries()))<30 && abs(int(USeg->GetEntries())-int(VSeg2->GetEntries()))>30) {
03864           TempTrackSeg[0].push_back(USeg); 
03865           TempTrackSeg[1].push_back(VSeg);
03866         }
03867         else if(abs(int(USeg->GetEntries())-int(VSeg->GetEntries()))>30 && abs(int(USeg->GetEntries())-int(VSeg2->GetEntries()))<30) {
03868           TempTrackSeg[0].push_back(USeg); 
03869           TempTrackSeg[1].push_back(VSeg2);       
03870         }
03871       }      
03872     }
03873     // Equivalent combination
03874     else if(NumUTrks==2 && NumVTrks==1) {
03875       TrackSegmentCam* USeg2=0;
03876 
03877       for(unsigned int i=0; i<ViewSegBank[0].size(); ++i) {
03878         if(ViewSegBank[0][i]->GetUID()>2) {USeg2=ViewSegBank[0][i]; break;}
03879       }
03880       
03881       if(VSeg && USeg && USeg2) {
03882         if(abs(int(VSeg->GetEntries())-int(USeg->GetEntries()))<30 && abs(int(VSeg->GetEntries())-int(USeg2->GetEntries()))>30) {
03883           TempTrackSeg[0].push_back(USeg); 
03884           TempTrackSeg[1].push_back(VSeg);
03885         }
03886         else if(abs(int(VSeg->GetEntries())-int(USeg->GetEntries()))>30 && abs(int(VSeg->GetEntries())-int(USeg2->GetEntries()))<30) {
03887           TempTrackSeg[0].push_back(USeg2); 
03888           TempTrackSeg[1].push_back(VSeg);        
03889         }
03890       }      
03891     }
03892   }
03894 
03895 
03897 
03898 
03899 
03900 
03901   // Print out this initial overlap information
03903   MSG("AlgTrackCamList", Msg::kVerbose) << " ** Guesses at tracks (only showing for UID>2)"<<endl;
03904   MSG("AlgTrackCamList", Msg::kVerbose) << "================================" << endl;
03905   for(unsigned int i=0; i<TempTrackSeg[0].size(); ++i) {
03906     TrackSegmentCam* Seg = TempTrackSeg[0][i];
03907     TrackSegmentCam* SegPartner = TempTrackSeg[1][i];
03908     if(Seg->GetUID()>2 && SegPartner->GetUID()>2)
03909       {
03910         MSG("AlgTrackCamList", Msg::kVerbose) << "  Seg    " 
03911                                               << "  begpln: "  << Seg->GetBegPlane()
03912                                               << "  begtpos: " << Seg->GetBegTPos()
03913                                               << "  endpln: "  << Seg->GetEndPlane()
03914                                               << "  endtpos: " << Seg->GetEndTPos() << endl;
03915         MSG("AlgTrackCamList", Msg::kVerbose) << "  Partner" 
03916                                               << "  begpln: "  << SegPartner->GetBegPlane()
03917                                               << "  begtpos: " << SegPartner->GetBegTPos()
03918                                               << "  endpln: "  << SegPartner->GetEndPlane()
03919                                               << "  endtpos: " << SegPartner->GetEndTPos() << endl;
03920         MSG("AlgTrackCamList", Msg::kVerbose) << "================================" << endl;
03921       }
03922   }
03924 
03925 
03926 
03927 
03928   // Special FD specific section to consider matching 2D tracks across the SM gap
03930   if(NumModules>1) {
03931 
03932     // Initialise container: Create one vector of empty segment vectors for each SM
03933     vector< vector<TrackSegmentCam*> > TempMod[2]; 
03934     for(unsigned int i=0; i<2; ++i) {
03935       for(unsigned int j=0; j<2; ++j) {vector<TrackSegmentCam*> TempVector; TempMod[i].push_back(TempVector);}
03936     }
03937 
03938     // Loop over the roughly overlapping 2D tracks we stored above
03940     for(unsigned int i=0; i<TempTrackSeg[0].size(); ++i) {
03941       TrackSegmentCam* Segu = TempTrackSeg[0][i];
03942       TrackSegmentCam* Segv = TempTrackSeg[1][i];
03943       
03944       // Check which module the U 2D track ends in
03945       Module = int((Segu->GetEndPlane())/(PlanesInModule+1));
03946       
03947       // Store U and V 2D tracks if they are suitably near the SM gap
03948       if( Segu->GetBegPlane()<Module*(PlanesInModule+1)+20
03949           && Segv->GetBegPlane()<Module*(PlanesInModule+1)+20 && Module>0) {
03950    
03951         TempMod[1][0].push_back(Segu);
03952         TempMod[1][1].push_back(Segv);
03953       }
03954 
03955       if( Segu->GetEndPlane()>(Module+1)*(PlanesInModule+1)-20
03956           && Segv->GetEndPlane()>(Module+1)*(PlanesInModule+1)-20 && Module<1) {
03957         
03958         TempMod[0][0].push_back(Segu);
03959         TempMod[0][1].push_back(Segv);
03960       }
03961     }
03963 
03964 
03965 
03966     // If there are U and V 2D tracks on either side of the SM gap, look for associations.
03967     // Also calculate a score based on the degree of overlap we would have if a join was made.
03969     if(TempMod[0][0].size()>0 && TempMod[1][0].size()>0) {
03970       Cont=true;
03971 
03972       while(Cont==true) {
03973         Cont=false;
03974         idm=-1; idp=-1; TopScore=0.; NClusters=0; TopClusters=0;
03975 
03976         // Loop over the 2D tracks in SM 1
03977         for(unsigned int i=0; i<TempMod[0][0].size(); ++i) {  
03978           TrackSegmentCam* Segum = TempMod[0][0][i];
03979           TrackSegmentCam* Segvm = TempMod[0][1][i];
03980         
03981           // For these, loop over the 2D tracks in SM 2
03982           if( (Segum->GetPartner()==Segvm && Segvm->GetPartner()==Segum)
03983               || (Segum->GetPartner()==0 && Segvm->GetPartner()==0
03984                   && Segum->GetUID()>2 && Segvm->GetUID()>2) ) {
03985             
03986             for(unsigned int j=0; j<TempMod[1][0].size(); ++j) {
03987               TrackSegmentCam* Segup = TempMod[1][0][j];
03988               TrackSegmentCam* Segvp = TempMod[1][1][j];
03989 
03990               // Check the association of the 2D tracks across the gap
03992               if(Segup->GetUID()>2 && Segvp->GetUID()>2) {
03993                 if(Segum->IsAssoc(Segup) && Segvm->IsAssoc(Segvp)) {
03994                   
03995                   if(Segum->GetBegPlane()<=Segvm->GetBegPlane()) {
03996                     begplane2=Segum->GetBegPlane();
03997                     begplane1=Segvm->GetBegPlane();
03998                   }
03999 
04000                   if(Segvm->GetBegPlane()<=Segum->GetBegPlane()) {
04001                     begplane2=Segvm->GetBegPlane();
04002                     begplane1=Segum->GetBegPlane();
04003                   }
04004 
04005                   if(Segup->GetEndPlane()>=Segvp->GetEndPlane()) {
04006                     endplane2=Segup->GetEndPlane();
04007                     endplane1=Segvp->GetEndPlane();
04008                   }
04009 
04010                   if(Segvp->GetEndPlane()>=Segup->GetEndPlane()) {
04011                     endplane2=Segvp->GetEndPlane();
04012                     endplane1=Segup->GetEndPlane();
04013                   }               
04014 
04015                   dplane1=endplane2-begplane2;
04016                   dplane2=endplane1-begplane1;
04017                   
04018                   Score=double(dplane2*dplane2)/double(dplane1);
04019 
04020                   NClusters=int(Segum->GetEntries()+Segvm->GetEntries()+Segup->GetEntries()+Segvp->GetEntries());
04021 
04022                   if(Score>TopScore || (Score==TopScore && NClusters>TopClusters)) {idm=i; idp=j; TopScore=Score; TopClusters=NClusters;}
04023                 }
04024               }
04026             }
04027           }
04028         }
04029 
04030         // Add 2D tracks together across the SM gap, as a result of the match.
04031         // Can also reliably set U/V partners. Set UIDs to -1 to avoid considering
04032         // the old 2D track that is now a subsection of a larger 2D track.
04034         if((1+idm)>0 && (1+idp)>0) {
04035           TrackSegmentCam* Segum = TempMod[0][0][idm];
04036           TrackSegmentCam* Segvm = TempMod[0][1][idm];
04037           TrackSegmentCam* Segup = TempMod[1][0][idp];
04038           TrackSegmentCam* Segvp = TempMod[1][1][idp];
04039 
04040           Segup->AddSegment(Segum); Segup->SetPartner(Segvp); 
04041           Segum->SetUID(-1); Segum->SetPartner(0);
04042           
04043           Segvp->AddSegment(Segvm); Segvp->SetPartner(Segup); 
04044           Segvm->SetUID(-1); Segvm->SetPartner(0);
04045 
04046           Cont=true;
04047         }
04049         
04050       } // End while Cont==true
04051     }
04053 
04054     // Clean up
04055     for(unsigned int i=0; i<2; ++i) {TempMod[i].clear();}
04056   }
04058   
04059 
04060 
04061 
04062   // Find other 2D track matches between views.
04064   Cont=true;
04065   
04066   while(Cont==true) {
04067     Cont=false;
04068     id=-1; TopScore=0.; NClusters=0; TopClusters=0;
04069 
04070     // Loop over the roughly overlapping 2D tracks we stored above
04071     for(unsigned int i=0; i<TempTrackSeg[0].size(); ++i) {
04072       TrackSegmentCam* Segu = TempTrackSeg[0][i];
04073       TrackSegmentCam* Segv = TempTrackSeg[1][i];
04074 
04075       // If have already established U/V partners, but not made the 3D track
04076       // do this now. TempTrack separately stores the U and V components of the
04077       // 3D track, which are currently still track segments, identified by UID==4.
04078       if(Segu->GetPartner()==Segv && Segv->GetPartner()==Segu
04079          && Segu->GetUID()<4 && Segv->GetUID()<4) {
04080 
04081         TempTrack[0].push_back(Segu);
04082         TempTrack[1].push_back(Segv);
04083 
04084         Segu->SetUID(4);
04085         Segv->SetUID(4);
04086       }
04087 
04088       // Otherwise, find the 2D tracks which have the best overlap
04090       else {
04091         if(Segu->GetPartner()==0 && Segv->GetPartner()==0
04092            && Segu->GetUID()>2 && Segv->GetUID()>2) {
04093           
04094           if(Segu->GetBegPlane()<=Segv->GetBegPlane()) {
04095             begplane2=Segu->GetBegPlane();
04096             begplane1=Segv->GetBegPlane();
04097           }
04098 
04099           if(Segv->GetBegPlane()<=Segu->GetBegPlane()) {
04100             begplane2=Segv->GetBegPlane();
04101             begplane1=Segu->GetBegPlane();
04102           }
04103 
04104           if(Segu->GetEndPlane()>=Segv->GetEndPlane()) {
04105             endplane2=Segu->GetEndPlane();
04106             endplane1=Segv->GetEndPlane();
04107           }
04108 
04109           if(Segv->GetEndPlane()>=Segu->GetEndPlane()) {
04110             endplane2=Segv->GetEndPlane();
04111             endplane1=Segu->GetEndPlane();
04112           }
04113 
04114           dplane1=endplane2-begplane2;
04115           dplane2=endplane1-begplane1;
04116                   
04117           Score=double(dplane2*dplane2)/double(dplane1);
04118 
04119           NClusters=int(Segu->GetEntries()+Segv->GetEntries());
04120 
04121           MSG("AlgTrackCamList", Msg::kVerbose) << "NClusters " << NClusters << " TopClusters "<< TopClusters 
04122                                                 << " Score " << Score << " TopScore " << TopScore << " id " << i << endl;
04123           
04124           //      if(Score>TopScore || (Score==TopScore && NClusters>TopClusters)) {id=i; TopScore=Score; TopClusters=NClusters;}
04125           if((Score>TopScore && NClusters>TopClusters-4)
04126              || (Score>0.4*TopScore && NClusters>TopClusters+4) 
04127              || (Score==TopScore && NClusters>TopClusters)) {id=i; TopScore=Score; TopClusters=NClusters;}
04128         }
04129       }
04131     }
04132 
04133     // Set the U/V partners for the 2D tracks having the best overlap
04135     if((1+id)>0) {
04136       MSG("AlgTrackCamList", Msg::kVerbose) << "Final id " << id << endl;
04137 
04138       TrackSegmentCam* Segu = TempTrackSeg[0][id];
04139       TrackSegmentCam* Segv = TempTrackSeg[1][id];
04140       
04141       Segu->SetPartner(Segv); 
04142       Segv->SetPartner(Segu);
04143       Cont=true;
04144     }
04146 
04147   } // End while Cont==true
04149 
04150 
04151 
04152   // Set the cluster flags for those clusters in a 3D track, and print out results
04154   MSG("AlgTrackCamList", Msg::kVerbose) << " ** Tracks Chosen ** " << endl;
04155   
04156   for(int View=0; View<2; ++View) {
04157     MSG("AlgTrackCamList", Msg::kVerbose) << "View: "<<View<<endl;
04158     MSG("AlgTrackCamList", Msg::kVerbose) << "================================" << endl;
04159     
04160     for(unsigned int i=0; i<TempTrack[View].size(); ++i) {
04161       TrackSegmentCam* Seg = TempTrack[View][i];
04162       TrackSegmentCam* SegPartner = Seg->GetPartner();
04163 
04164       MSG("AlgTrackCamList", Msg::kVerbose) << "  Seg    " 
04165                                             << "  begpln: "  << Seg->GetBegPlane()
04166                                             << "  begtpos: " << Seg->GetBegTPos()
04167                                             << "  endpln: "  << Seg->GetEndPlane()
04168                                             << "  endtpos: " << Seg->GetEndTPos() << endl;
04169       MSG("AlgTrackCamList", Msg::kVerbose) << "  Partner" 
04170                                             << "  begpln: "  << SegPartner->GetBegPlane()
04171                                             << "  begtpos: " << SegPartner->GetBegTPos()
04172                                             << "  endpln: "  << SegPartner->GetEndPlane()
04173                                             << "  endtpos: " << SegPartner->GetEndTPos() << endl;
04174       MSG("AlgTrackCamList", Msg::kVerbose) << "================================" << endl;
04175 
04176       for(unsigned int j=0; j<Seg->GetEntries(); ++j) {
04177         Seg->GetCluster(j)->SetTrkFlag(4);      
04178       }
04179     }
04180   }
04182 
04183   return;
04184 }

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

Reimplemented from AlgBase.

Definition at line 5703 of file AlgTrackCamList.cxx.

05704 {
05705 }


Member Data Documentation

vector<HitCam*> AlgTrackCamList::AllHitBank[500] [private]

Definition at line 82 of file AlgTrackCamList.h.

Referenced by FirstUVComparison(), Form3DTracks(), and RunAlg().

vector<ClusterCam*> AlgTrackCamList::ClusterBank[500] [private]
vector<ClusterCam*> AlgTrackCamList::ClusterList[500] [private]

Definition at line 61 of file AlgTrackCamList.h.

Referenced by ClearUp(), and Form3DTracks().

Definition at line 71 of file AlgTrackCamList.h.

Referenced by ClearUp(), Form3DTracks(), and RunAlg().

vector<HitCam*> AlgTrackCamList::HitBank[500] [private]

Definition at line 55 of file AlgTrackCamList.h.

Referenced by ClearUp(), FirstUVComparison(), FormTheClusters(), FormTheHits(), and MatchUV().

vector<HitCam*> AlgTrackCamList::HitsInTracks[2] [private]

Definition at line 64 of file AlgTrackCamList.h.

Referenced by ClearUp(), and NearDetectorTriplets().

double AlgTrackCamList::PECut [private]

Definition at line 80 of file AlgTrackCamList.h.

Referenced by FormTheHits(), IDTrkAndShwClusters(), MatchUV(), and RunAlg().

double AlgTrackCamList::PECut2 [private]

Definition at line 80 of file AlgTrackCamList.h.

Referenced by FormTheHits(), and RunAlg().

Definition at line 69 of file AlgTrackCamList.h.

Referenced by Form2DTracks(), and Join2DTracks().

const double AlgTrackCamList::StripWidth [private]

Definition at line 67 of file AlgTrackCamList.h.

Referenced by ClearUp(), Form3DTracks(), and SecondUVComparison().

Definition at line 74 of file AlgTrackCamList.h.

Referenced by FirstUVComparison(), 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