+ Modifications of the standard tracking algorithm:
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackReconstructorK.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ////////////////////////////////////
19 //
20 // MUON track reconstructor using the kalman method
21 //
22 // This class contains as data:
23 // * the parameters for the track reconstruction
24 //
25 // It contains as methods, among others:
26 // * MakeTracks to build the tracks
27 //
28 ////////////////////////////////////
29
30 #include <stdlib.h>
31 #include <Riostream.h>
32 #include <TDirectory.h>
33 #include <TFile.h>
34 #include <TMatrixD.h>
35
36 #include "AliMUONVTrackReconstructor.h"
37 #include "AliMUONTrackReconstructorK.h"
38 #include "AliMUONData.h"
39 #include "AliMUONConstants.h"
40 #include "AliMUONHitForRec.h"
41 #include "AliMUONObjectPair.h"
42 #include "AliMUONRawCluster.h"
43 #include "AliMUONTrackK.h" 
44 #include "AliLog.h"
45
46 ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
47
48 //__________________________________________________________________________
49 AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(AliMUONData* data, const Option_t* TrackMethod)
50   : AliMUONVTrackReconstructor(data),
51     fTrackMethod(2), //tracking method (2-Kalman 3-Combination-Kalman/Clustering)
52     fMuons(0)
53 {
54   /// Constructor for class AliMUONTrackReconstructorK
55
56   if (strstr(TrackMethod,"Kalman")) {
57     cout << " *** Tracking with the Kalman filter *** " << endl;
58     fTrackMethod = 2;
59   } else if (strstr(TrackMethod,"Combi")) {
60     cout << " *** Combined cluster / track finder ***" << endl;
61     fTrackMethod = 3;
62   } else AliFatal(Form("Tracking method %s not available",TrackMethod));
63   
64   // Memory allocation for the TClonesArray of reconstructed tracks
65   fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
66 }
67
68   //__________________________________________________________________________
69 AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK(void)
70 {
71   /// Destructor for class AliMUONTrackReconstructorK
72   delete fRecTracksPtr;
73 }
74
75   //__________________________________________________________________________
76 void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
77 {
78   /// To add to the list of hits for reconstruction all the raw clusters
79   TTree *treeR;
80   AliMUONHitForRec *hitForRec;
81   AliMUONRawCluster *clus;
82   Int_t iclus, nclus;
83   TClonesArray *rawclusters;
84   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
85
86   treeR = fMUONData->TreeR();
87   if (!treeR) {
88     AliError("TreeR must be loaded");
89     exit(0);
90   }
91   
92   if (fTrackMethod != 3) { //AZ
93     fMUONData->SetTreeAddress("RC"); //AZ
94     fMUONData->GetRawClusters(); // only one entry  
95   }
96
97   // Loop over tracking chambers
98   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
99     rawclusters =fMUONData->RawClusters(ch);
100     nclus = (Int_t) (rawclusters->GetEntries());
101     // Loop over (cathode correlated) raw clusters
102     for (iclus = 0; iclus < nclus; iclus++) {
103       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
104       // new AliMUONHitForRec from raw cluster
105       // and increment number of AliMUONHitForRec's (total and in chamber)
106       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
107       fNHitsForRec++;
108       // more information into HitForRec
109       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
110       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
111       //  original raw cluster
112       hitForRec->SetChamberNumber(ch);
113       hitForRec->SetHitNumber(iclus);
114       // Z coordinate of the raw cluster (cm)
115       hitForRec->SetZ(clus->GetZ(0));
116       StdoutToAliDebug(3,
117                        cout << "Chamber " << ch <<
118                        " raw cluster  " << iclus << " : " << endl;
119                        clus->Print("full");
120                        cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
121                        hitForRec->Print("full");
122                        );
123     } // end of cluster loop
124   } // end of chamber loop
125   SortHitsForRecWithIncreasingChamber(); 
126   
127   AliDebug(1,"End of AddHitsForRecFromRawClusters");
128     if (AliLog::GetGlobalDebugLevel() > 0) {
129       AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
130       for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
131         AliDebug(1, Form("Chamber(0...): %d",ch));
132         AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
133         AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
134         for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
135              hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
136              hit++) {
137           AliDebug(1, Form("HitForRec index(0...): %d",hit));
138           ((*fHitsForRecPtr)[hit])->Dump();
139       }
140     }
141   }
142   
143   return;
144 }
145
146   //__________________________________________________________________________
147 void AliMUONTrackReconstructorK::MakeTracks(void)
148 {
149   /// To make the tracks from the list of segments and points in all stations
150   AliDebug(1,"Enter MakeTracks");
151   // The order may be important for the following Reset's
152   //AZ ResetTracks();
153   MakeTrackCandidates();
154   if (fRecTracksPtr->GetEntriesFast() == 0) return;
155   // Follow tracks in stations(1..) 3, 2 and 1
156   FollowTracks();
157   // Remove double tracks
158   RemoveDoubleTracks();
159   // Propagate tracks to the vertex through absorber
160   ExtrapTracksToVertex();
161   // Fill AliMUONTrack data members
162   FillMUONTrack();
163 }
164
165   //__________________________________________________________________________
166 void AliMUONTrackReconstructorK::MakeTrackCandidates(void)
167 {
168   /// To make initial tracks for Kalman filter from segments in stations(1..)  4 and 5
169   Int_t istat, iseg;
170   TClonesArray *segments;
171   AliMUONObjectPair *segment;
172   AliMUONTrackK *trackK;
173
174   AliDebug(1,"Enter MakeTrackCandidatesK");
175
176   AliMUONTrackK a(this, fHitsForRecPtr);
177   // Loop over stations(1...) 5 and 4
178   for (istat=4; istat>=3; istat--) {
179     // Make segments in the station
180     segments = MakeSegmentsInStation(istat);
181     // Loop over segments in the station
182     for (iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
183       // Transform segments to tracks
184       segment = (AliMUONObjectPair*) ((*segments)[iseg]);
185       trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
186     } // for (iseg=0;...)
187   } // for (istat=4;...)
188   return;
189 }
190
191 //__________________________________________________________________________
192 void AliMUONTrackReconstructorK::FollowTracks(void)
193 {
194   /// Follow tracks using Kalman filter
195   Bool_t ok;
196   Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
197   AliMUONTrackK *trackK;
198   AliMUONHitForRec *hit;
199   AliMUONRawCluster *clus;
200   TClonesArray *rawclusters;
201   clus = 0; rawclusters = 0;
202
203   Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
204   Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
205   Double_t zDipole1 = simpleBPosition + 0.5 * simpleBLength;
206   Double_t zDipole2 = zDipole1 - simpleBLength;
207
208   // Print hits
209   trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);
210
211   if (trackK->DebugLevel() > 0) {
212     for (Int_t i1=0; i1<fNHitsForRec; i1++) {
213       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
214       printf(" Hit # %d %10.4f %10.4f %10.4f",
215              hit->GetChamberNumber(), hit->GetBendingCoor(),
216              hit->GetNonBendingCoor(), hit->GetZ());
217  
218       // from raw clusters
219       rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
220       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
221                                                            GetHitNumber());
222       printf(" %d", clus->GetTrack(1));
223       if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
224       else printf("\n");
225      
226     }
227   } // if (trackK->DebugLevel() > 0)
228
229   icand = -1;
230   Int_t nSeeds;
231   nSeeds = fNRecTracks; // starting number of seeds
232   // Loop over track candidates
233   while (icand < fNRecTracks-1) {
234     icand ++;
235     if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
236     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
237     if (trackK->GetRecover() < 0) continue; // failed track
238
239     // Discard candidate which will produce the double track
240     /*
241     if (icand > 0) {
242       ok = CheckCandidate(icand,nSeeds);
243       if (!ok) {
244         trackK->SetRecover(-1); // mark candidate to be removed
245         continue;
246       }
247     }
248     */
249
250     ok = kTRUE;
251     if (trackK->GetRecover() == 0) 
252       hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
253     else 
254       hit = trackK->GetHitLastOk(); // hit where track stopped
255
256     if (hit) ichamBeg = hit->GetChamberNumber();
257     ichamEnd = 0;
258     // Check propagation direction
259     if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
260     else if (trackK->GetTrackDir() < 0) {
261       ichamEnd = 9; // forward propagation
262       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
263       if (ok) {
264         ichamBeg = ichamEnd;
265         ichamEnd = 6; // backward propagation
266         // Change weight matrix and zero fChi2 for backpropagation
267         trackK->StartBack();
268         trackK->SetTrackDir(1);
269         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
270         ichamBeg = ichamEnd;
271         ichamEnd = 0;
272       }
273     } else {
274       if (trackK->GetBPFlag()) {
275         // backpropagation
276         ichamEnd = 6; // backward propagation
277         // Change weight matrix and zero fChi2 for backpropagation
278         trackK->StartBack();
279         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
280         ichamBeg = ichamEnd;
281         ichamEnd = 0;
282       }
283     }
284
285     if (ok) {
286       trackK->SetTrackDir(1);
287       trackK->SetBPFlag(kFALSE);
288       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
289     }
290     if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed
291
292     // Apply smoother
293     if (trackK->GetRecover() >= 0) {
294       ok = trackK->Smooth();
295       if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
296     }
297
298     // Majority 3 of 4 in first 2 stations
299     if (!ok) continue;
300     chamBits = 0;
301     Double_t chi2max = 0;
302     for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
303       hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
304       chamBits |= BIT(hit->GetChamberNumber());
305       if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
306     }
307     if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) {
308       //trackK->Recover();
309       trackK->SetRecover(-1); //mark candidate to be removed
310       continue;
311     }
312     if (ok) trackK->SetTrackQuality(0); // compute "track quality"
313   } // while
314
315   for (Int_t i=0; i<fNRecTracks; i++) {
316     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
317     if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
318   }
319
320   // Compress TClonesArray
321   fRecTracksPtr->Compress();
322   fNRecTracks = fRecTracksPtr->GetEntriesFast();
323   return;
324 }
325
326 //__________________________________________________________________________
327 Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) const
328 {
329   /// Discards track candidate if it will produce the double track (having
330   /// the same seed segment hits as hits of a good track found before)
331   AliMUONTrackK *track1, *track2;
332   AliMUONHitForRec *hit1, *hit2, *hit;
333   AliMUONObjectPair *segment1, *segment2;
334
335   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
336   hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
337   hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit
338
339   for (Int_t i=0; i<icand; i++) {
340     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
341     //if (track2->GetRecover() < 0) continue;
342     if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
343     segment1 = track1->GetStartSegment();
344     segment2 = track2->GetStartSegment();
345     if (segment1->First()  == segment2->First() &&
346         segment1->Second() == segment2->Second()) {
347       return kFALSE;
348     } else {
349       Int_t nSame = 0;
350       for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
351         hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
352         if (hit == hit1 || hit == hit2) {
353           nSame++;
354           if (nSame == 2) return kFALSE;
355         }
356       } // for (Int_t j=0;
357     }
358   } // for (Int_t i=0;
359   return kTRUE;
360 }
361
362   //__________________________________________________________________________
363 void AliMUONTrackReconstructorK::RemoveDoubleTracks(void)
364 {
365   /// Removes double tracks (sharing more than half of their hits).
366   /// Keeps the track with higher quality
367   AliMUONTrackK *track1, *track2, *trackToKill;
368
369   // Sort tracks according to their quality
370   fRecTracksPtr->Sort();
371
372   // Loop over first track of the pair
373   track1 = (AliMUONTrackK*) fRecTracksPtr->First();
374   Int_t debug = track1->DebugLevel();
375   while (track1) {
376     // Loop over second track of the pair
377     track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
378     while (track2) {
379       // Check whether or not to keep track2
380       if (!track2->KeepTrack(track1)) {
381         if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
382           " " << track2->GetTrackQuality() << endl;
383         trackToKill = track2;
384         track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
385         trackToKill->Kill();
386         fRecTracksPtr->Compress();
387       } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
388     } // track2
389     track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
390   } // track1
391
392   fNRecTracks = fRecTracksPtr->GetEntriesFast();
393   if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
394 }
395
396   //__________________________________________________________________________
397 void AliMUONTrackReconstructorK::ExtrapTracksToVertex(void)
398 {
399   /// Propagates track to the vertex thru absorber
400   /// (using Branson correction for now)
401   Double_t zVertex;
402   zVertex = 0;
403   for (Int_t i=0; i<fNRecTracks; i++) {
404     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson();
405     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
406     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
407     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(1); // with absorber
408   }
409 }
410
411   //__________________________________________________________________________
412 void AliMUONTrackReconstructorK::FillMUONTrack()
413 {
414   // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
415   AliMUONTrackK *track;
416   track = (AliMUONTrackK*) fRecTracksPtr->First();
417   while (track) {
418     track->FillMUONTrack();
419     track = (AliMUONTrackK*) fRecTracksPtr->After(track);
420   } 
421   return;
422 }
423
424   //__________________________________________________________________________
425 void AliMUONTrackReconstructorK::EventDump(void)
426 {
427   /// Dump reconstructed event (track parameters at vertex and at first hit),
428   /// and the particle parameters
429   AliDebug(1,"****** enter EventDump ******");
430   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks)); 
431   
432   AliWarning("AliMUONTrackReconstructorK::EventDump: not implemented"); 
433   
434   // informations about generated particles NO !!!!!!!!
435   
436 //    for (Int_t iPart = 0; iPart < np; iPart++) {
437 //      p = gAlice->Particle(iPart);
438 //      printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
439 //         iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());    
440 //    }
441   return;
442 }
443
444