858c3b07d328c2491897db6cef5041ef93518830
[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 "AliMUONRawCluster.h"
42 #include "AliMUONSegment.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, nTRentries;
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     nTRentries = Int_t(treeR->GetEntries());
94     
95     if (!(fMUONData->IsRawClusterBranchesInTree())) {
96       AliError(Form("RawCluster information is not avalaible, nTRentries = %d not equal to 1",nTRentries));
97       exit(0);
98     }
99
100     fMUONData->SetTreeAddress("RC"); //AZ
101     fMUONData->GetRawClusters(); // only one entry  
102   }
103
104   // Loop over tracking chambers
105   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
106     // number of HitsForRec to 0 for the chamber
107     fNHitsForRecPerChamber[ch] = 0;
108     // index of first HitForRec for the chamber
109     if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
110     else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
111     rawclusters =fMUONData->RawClusters(ch);
112     nclus = (Int_t) (rawclusters->GetEntries());
113     // Loop over (cathode correlated) raw clusters
114     for (iclus = 0; iclus < nclus; iclus++) {
115       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
116       // new AliMUONHitForRec from raw cluster
117       // and increment number of AliMUONHitForRec's (total and in chamber)
118       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
119       fNHitsForRec++;
120       (fNHitsForRecPerChamber[ch])++;
121       // more information into HitForRec
122       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
123       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
124       //  original raw cluster
125       hitForRec->SetChamberNumber(ch);
126       hitForRec->SetHitNumber(iclus);
127       // Z coordinate of the raw cluster (cm)
128       hitForRec->SetZ(clus->GetZ(0));
129       
130       StdoutToAliDebug(3,
131                        cout << "Chamber " << ch <<
132                        " raw cluster  " << iclus << " : " << endl;
133                        clus->Print("full");
134                        cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
135                        hitForRec->Print("full");
136                        );
137     } // end of cluster loop
138   } // end of chamber loop
139   SortHitsForRecWithIncreasingChamber(); 
140   
141   AliDebug(1,"End of AddHitsForRecFromRawClusters");
142     if (AliLog::GetGlobalDebugLevel() > 0) {
143       AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
144       for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
145         AliDebug(1, Form("Chamber(0...): %d",ch));
146         AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
147         AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
148         for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
149              hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
150              hit++) {
151           AliDebug(1, Form("HitForRec index(0...): %d",hit));
152           ((*fHitsForRecPtr)[hit])->Dump();
153       }
154     }
155   }
156   
157   return;
158 }
159
160   //__________________________________________________________________________
161 void AliMUONTrackReconstructorK::MakeSegments(void)
162 {
163   /// To make the list of segments in all stations,
164   /// from the list of hits to be reconstructed
165   AliDebug(1,"Enter MakeSegments");
166   //AZ ResetSegments();
167   // Loop over stations
168   for (Int_t st = 3; st < AliMUONConstants::NTrackingCh()/2; st++) MakeSegmentsPerStation(st); 
169   
170   StdoutToAliDebug(3,
171     cout << "end of MakeSegments" << endl;
172     for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) 
173     {
174       cout << "station " << st
175             << "  has " << fNSegments[st] << " segments:"
176             << endl;
177       for (Int_t seg = 0; seg < fNSegments[st]; seg++) 
178       {
179               ((*fSegmentsPtr[st])[seg])->Print();
180       }
181     }
182                    );
183   return;
184 }
185
186   //__________________________________________________________________________
187 void AliMUONTrackReconstructorK::MakeTracks(void)
188 {
189   /// To make the tracks from the list of segments and points in all stations
190   AliDebug(1,"Enter MakeTracks");
191   // The order may be important for the following Reset's
192   //AZ ResetTracks();
193   MakeTrackCandidates();
194   if (fRecTracksPtr->GetEntriesFast() == 0) return;
195   // Follow tracks in stations(1..) 3, 2 and 1
196   FollowTracks();
197   // Remove double tracks
198   RemoveDoubleTracks();
199   // Propagate tracks to the vertex thru absorber
200   GoToVertex();
201   // Fill AliMUONTrack data members
202   FillMUONTrack();
203 }
204
205   //__________________________________________________________________________
206 void AliMUONTrackReconstructorK::MakeTrackCandidates(void)
207 {
208   /// To make initial tracks for Kalman filter from the list of segments
209   Int_t istat, iseg;
210   AliMUONSegment *segment;
211   AliMUONTrackK *trackK;
212
213   AliDebug(1,"Enter MakeTrackCandidatesK");
214
215   AliMUONTrackK a(this, fHitsForRecPtr);
216   // Loop over stations(1...) 5 and 4
217   for (istat=4; istat>=3; istat--) {
218     // Loop over segments in the station
219     for (iseg=0; iseg<fNSegments[istat]; iseg++) {
220       // Transform segments to tracks and evaluate covariance matrix
221       segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
222       trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
223     } // for (iseg=0;...)
224   } // for (istat=4;...)
225   return;
226 }
227
228 //__________________________________________________________________________
229 void AliMUONTrackReconstructorK::FollowTracks(void)
230 {
231   /// Follow tracks using Kalman filter
232   Bool_t ok;
233   Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
234   Double_t zDipole1, zDipole2;
235   AliMUONTrackK *trackK;
236   AliMUONHitForRec *hit;
237   AliMUONRawCluster *clus;
238   TClonesArray *rawclusters;
239   clus = 0; rawclusters = 0;
240
241   zDipole1 = fSimpleBPosition + fSimpleBLength/2;
242   zDipole2 = zDipole1 - fSimpleBLength;
243
244   // Print hits
245   trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);
246
247   if (trackK->DebugLevel() > 0) {
248     for (Int_t i1=0; i1<fNHitsForRec; i1++) {
249       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
250       printf(" Hit # %d %10.4f %10.4f %10.4f",
251              hit->GetChamberNumber(), hit->GetBendingCoor(),
252              hit->GetNonBendingCoor(), hit->GetZ());
253  
254       // from raw clusters
255       rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
256       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
257                                                            GetHitNumber());
258       printf(" %d", clus->GetTrack(1));
259       if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
260       else printf("\n");
261      
262     }
263   } // if (trackK->DebugLevel() > 0)
264
265   icand = -1;
266   Int_t nSeeds;
267   nSeeds = fNRecTracks; // starting number of seeds
268   // Loop over track candidates
269   while (icand < fNRecTracks-1) {
270     icand ++;
271     if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
272     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
273     if (trackK->GetRecover() < 0) continue; // failed track
274
275     // Discard candidate which will produce the double track
276     /*
277     if (icand > 0) {
278       ok = CheckCandidate(icand,nSeeds);
279       if (!ok) {
280         trackK->SetRecover(-1); // mark candidate to be removed
281         continue;
282       }
283     }
284     */
285
286     ok = kTRUE;
287     if (trackK->GetRecover() == 0) 
288       hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
289     else 
290       hit = trackK->GetHitLastOk(); // hit where track stopped
291
292     if (hit) ichamBeg = hit->GetChamberNumber();
293     ichamEnd = 0;
294     // Check propagation direction
295     if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
296     else if (trackK->GetTrackDir() < 0) {
297       ichamEnd = 9; // forward propagation
298       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
299       if (ok) {
300         ichamBeg = ichamEnd;
301         ichamEnd = 6; // backward propagation
302         // Change weight matrix and zero fChi2 for backpropagation
303         trackK->StartBack();
304         trackK->SetTrackDir(1);
305         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
306         ichamBeg = ichamEnd;
307         ichamEnd = 0;
308       }
309     } else {
310       if (trackK->GetBPFlag()) {
311         // backpropagation
312         ichamEnd = 6; // backward propagation
313         // Change weight matrix and zero fChi2 for backpropagation
314         trackK->StartBack();
315         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
316         ichamBeg = ichamEnd;
317         ichamEnd = 0;
318       }
319     }
320
321     if (ok) {
322       trackK->SetTrackDir(1);
323       trackK->SetBPFlag(kFALSE);
324       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
325     }
326     if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed
327
328     // Apply smoother
329     if (trackK->GetRecover() >= 0) {
330       ok = trackK->Smooth();
331       if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
332     }
333
334     // Majority 3 of 4 in first 2 stations
335     if (!ok) continue;
336     chamBits = 0;
337     Double_t chi2max = 0;
338     for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
339       hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
340       chamBits |= BIT(hit->GetChamberNumber());
341       if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
342     }
343     if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) {
344       //trackK->Recover();
345       trackK->SetRecover(-1); //mark candidate to be removed
346       continue;
347     }
348     if (ok) trackK->SetTrackQuality(0); // compute "track quality"
349   } // while
350
351   for (Int_t i=0; i<fNRecTracks; i++) {
352     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
353     if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
354   }
355
356   // Compress TClonesArray
357   fRecTracksPtr->Compress();
358   fNRecTracks = fRecTracksPtr->GetEntriesFast();
359   return;
360 }
361
362 //__________________________________________________________________________
363 Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) const
364 {
365   /// Discards track candidate if it will produce the double track (having
366   /// the same seed segment hits as hits of a good track found before)
367   AliMUONTrackK *track1, *track2;
368   AliMUONHitForRec *hit1, *hit2, *hit;
369
370   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
371   hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
372   hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit
373
374   for (Int_t i=0; i<icand; i++) {
375     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
376     //if (track2->GetRecover() < 0) continue;
377     if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
378
379     if (track1->GetStartSegment() == track2->GetStartSegment()) {
380       return kFALSE;
381     } else {
382       Int_t nSame = 0;
383       for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
384         hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
385         if (hit == hit1 || hit == hit2) {
386           nSame++;
387           if (nSame == 2) return kFALSE;
388         }
389       } // for (Int_t j=0;
390     }
391   } // for (Int_t i=0;
392   return kTRUE;
393 }
394
395   //__________________________________________________________________________
396 void AliMUONTrackReconstructorK::RemoveDoubleTracks(void)
397 {
398   /// Removes double tracks (sharing more than half of their hits).
399   /// Keeps the track with higher quality
400   AliMUONTrackK *track1, *track2, *trackToKill;
401
402   // Sort tracks according to their quality
403   fRecTracksPtr->Sort();
404
405   // Loop over first track of the pair
406   track1 = (AliMUONTrackK*) fRecTracksPtr->First();
407   Int_t debug = track1->DebugLevel();
408   while (track1) {
409     // Loop over second track of the pair
410     track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
411     while (track2) {
412       // Check whether or not to keep track2
413       if (!track2->KeepTrack(track1)) {
414         if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
415           " " << track2->GetTrackQuality() << endl;
416         trackToKill = track2;
417         track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
418         trackToKill->Kill();
419         fRecTracksPtr->Compress();
420       } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
421     } // track2
422     track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
423   } // track1
424
425   fNRecTracks = fRecTracksPtr->GetEntriesFast();
426   if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
427 }
428
429   //__________________________________________________________________________
430 void AliMUONTrackReconstructorK::GoToVertex(void)
431 {
432   /// Propagates track to the vertex thru absorber
433   /// (using Branson correction for now)
434   Double_t zVertex;
435   zVertex = 0;
436   for (Int_t i=0; i<fNRecTracks; i++) {
437     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson();
438     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
439     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
440     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(1); // with absorber
441   }
442 }
443
444   //__________________________________________________________________________
445 void AliMUONTrackReconstructorK::FillMUONTrack()
446 {
447   // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
448   AliMUONTrackK *track;
449   track = (AliMUONTrackK*) fRecTracksPtr->First();
450   while (track) {
451     track->FillMUONTrack();
452     track = (AliMUONTrackK*) fRecTracksPtr->After(track);
453   } 
454   return;
455 }
456
457   //__________________________________________________________________________
458 void AliMUONTrackReconstructorK::EventDump(void)
459 {
460   /// Dump reconstructed event (track parameters at vertex and at first hit),
461   /// and the particle parameters
462   AliDebug(1,"****** enter EventDump ******");
463   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks)); 
464   
465   AliWarning("AliMUONTrackReconstructorK::EventDump: not implemented"); 
466   
467   // informations about generated particles NO !!!!!!!!
468   
469 //    for (Int_t iPart = 0; iPart < np; iPart++) {
470 //      p = gAlice->Particle(iPart);
471 //      printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
472 //         iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());    
473 //    }
474   return;
475 }
476
477