Compilation warnings
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackFinder.cxx
1 // ROOT includes
2 #include "TObject.h"
3 #include "TClonesArray.h"
4 #include "TObjArray.h"
5 #include "TH2D.h"
6 #include "TH1D.h"
7 #include "TFile.h"
8 #include "TGeoManager.h"
9 #include "TMatrixD.h"
10 #include "TParticle.h"
11 #include "TMath.h"
12 #include "TGraph.h"
13 #include "TEllipse.h"
14 #include "TCanvas.h"
15 #include "TString.h"
16 #include "TLatex.h"
17 #include "TMarker.h"
18 #include "TNtuple.h"
19 #include "TRandom.h"
20 #include "TIterator.h"
21
22 // STEER includes
23 #include "AliLog.h"
24 #include "AliRun.h"
25 #include "AliRunLoader.h"
26 #include "AliLoader.h"
27 #include "AliHeader.h"
28 #include "AliMC.h"
29 #include "AliStack.h"
30 #include "AliMagF.h"
31 #include "AliTracker.h"
32 #include "AliGRPObject.h"
33 #include "AliCDBEntry.h"
34 #include "AliCDBManager.h"
35
36 // MUON includes
37 #include "AliMUONConstants.h"
38 #include "AliMUONTrack.h"
39 #include "AliMUONRecoCheck.h"
40 #include "AliMUONTrackParam.h"
41 #include "AliMUONTrackExtrap.h"
42 #include "AliMUONVTrackStore.h"
43 #include "AliMUONVCluster.h"
44
45 // MFT includes
46 #include "AliMuonForwardTrack.h"
47 #include "AliMFTCluster.h"
48 #include "AliMFT.h"
49 #include "AliMFTSegmentation.h"
50
51 #include "AliMuonForwardTrackFinder.h"
52
53 //====================================================================================================================================================
54 //
55 // Class for the creation of the "global muon tracks" built from the clusters in the 
56 // muon spectrometer and the clusters of the Muon Forward Tracker. QA histograms are also created
57 //
58 // Contact author: antonio.uras@cern.ch
59 //
60 //====================================================================================================================================================
61
62 ClassImp(AliMuonForwardTrackFinder)
63
64 //=====================================================================================================
65
66 AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
67   TObject(),
68   fRun(0),
69   fNEventsToAnalyze(0),
70   fSigmaClusterCut(0),
71   fChi2GlobalCut(0),
72   fSigmaSpectrometerCut(0),
73   fExtrapOriginTransvError(0),
74   fGaussianBlurZVert(0),
75   fNFinalCandidatesCut(0),
76   fReadDir(0),
77   fOutDir(0),
78   fDrawOption(0),
79
80   fDistanceFromGoodClusterAndTrackAtLastPlane(-1),
81   fDistanceFromBestClusterAndTrackAtLastPlane(-1),
82   
83   fRAbsorberCut(0),
84   fLowPtCut(0),
85   fNPlanesMFT(0),
86   fNPlanesMFTAnalyzed(0),
87   fNMaxMissingMFTClusters(0),
88
89   fEv(0),
90   fLabelMC(0),
91
92   fHistPtSpectrometer(0), 
93   fHistPtMuonTrackWithGoodMatch(0),
94   fHistPtMuonTrackWithBadMatch(0),
95   fHistRadiusEndOfAbsorber(0), 
96   fHistNGoodClustersForFinalTracks(0),
97   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane(0),
98   fHistDistanceGoodClusterFromTrackAtLastPlane(0),
99
100   fNtuFinalCandidates1(0),
101   fNtuFinalBestCandidates1(0),
102   fNtuFinalCandidates2(0),
103   fNtuFinalBestCandidates2(0),
104
105   fCanvas(0),
106
107   fTxtMuonHistory(0), 
108   fTxtTrackGoodClusters(0), 
109   fTxtTrackFinalChi2(0), 
110   fTxtFinalCandidates(0), 
111   fTxtDummy(0),
112   fTxtAllClust(0), 
113   fTxtClustGoodChi2(0), 
114   fTxtClustMC(0), 
115   fTxtClustOfTrack(0), 
116   fMrkAllClust(0), 
117   fMrkClustGoodChi2(0), 
118   fMrkClustMC(0), 
119   fMrkClustOfTrack(0),
120
121   fCountRealTracksAnalyzed(0), 
122   fCountRealTracksWithRefMC(0), 
123   fCountRealTracksWithRefMC_andTrigger(0),
124   fCountRealTracksWithRefMC_andTrigger_andGoodPt(0),
125   fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta(0),
126   fCountRealTracksAnalyzedOfEvent(0),
127   fCountRealTracksAnalyzedWithFinalCandidates(0),
128
129   fFileCluster(0),
130   fFileESD(0),
131   fFile_gAlice(0),
132
133   fRunLoader(0),
134   fMFTLoader(0),
135   fMuonRecoCheck(0),
136   fMFTClusterTree(0),
137   fMuonTrackReco(0),
138   fCurrentTrack(0),
139   fIsCurrentMuonTrackable(0),
140   fCandidateTracks(0),
141   fTrackStore(0),
142   fTrackRefStore(0),
143   fNextTrack(0),
144   fStack(0),
145   fMFT(0),
146   fSegmentation(0),
147   fOutputTreeFile(0),
148   fOutputEventTree(0),
149   fMuonForwardTracks(0),
150   fMatchingMode(-1),
151   fMinResearchRadiusAtLastPlane(0),
152   fGRPData(0),
153   fRunInfo(0)
154
155 {
156
157   // Default constructor
158
159   for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) {
160
161     fHistNTracksAfterExtrapolation[iPlane] = 0;
162     fHistChi2Cluster_GoodCluster[iPlane] = 0;
163     fHistChi2Cluster_BadCluster[iPlane] = 0;
164     fHistResearchRadius[iPlane] = 0;      
165     
166     fIsGoodClusterInPlane[iPlane] = kFALSE;
167     
168     fHistChi2Cluster_GoodCluster[iPlane] = 0;
169     fHistChi2Cluster_BadCluster[iPlane]  = 0;
170     
171     fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = 0;
172     fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = 0;
173     
174     fZPlane[iPlane] = 0.;
175     fRPlaneMax[iPlane] = 0.;
176     fRPlaneMin[iPlane] = 0.;
177     
178     for (Int_t i=0; i<4; i++) fGrMFTPlane[iPlane][i] = 0;
179     fCircleExt[iPlane] = 0;
180     fCircleInt[iPlane] = 0;
181     
182     fTxtTrackChi2[iPlane] = 0;
183     
184     fIsClusterCompatible[iPlane] = 0;
185     
186     fMFTClusterArray[iPlane]      = 0;
187     fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
188     fMFTClusterArrayBack[iPlane]  = new TClonesArray("AliMFTCluster");
189
190     fIsPlaneMandatory[iPlane] = kFALSE;
191     
192   }
193
194   //  fNextTrack = 0;
195
196   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane = 0;
197   fHistDistanceGoodClusterFromTrackAtLastPlane = 0;
198
199   fMFTClusterTree = 0;
200   fCandidateTracks = 0;
201
202   fOutputTreeFile = new TFile("MuonGlobalTracks.root", "recreate");
203   fOutputEventTree   = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
204   fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack");
205   fOutputEventTree   -> Branch("tracks", &fMuonForwardTracks);
206
207 }
208
209 //=====================================================================================================
210
211 void AliMuonForwardTrackFinder::Init(Int_t nRun, 
212                                      Char_t *readDir,
213                                      Char_t *outDir,
214                                      Int_t nEventsToAnalyze) {
215   
216   if (fRunLoader) {
217     printf("WARNING: run already initialized!!\n");
218   }
219
220   SetRun(nRun);
221   SetReadDir(readDir);
222   SetOutDir(outDir);
223
224   printf("input  dir = %s\n", fReadDir.Data());
225   printf("output dir = %s\n", fOutDir.Data());
226
227   // -------------------------- initializing files...
228
229   printf("initializing files for run %d...\n", fRun);
230
231   Char_t geoFileName[300];
232   Char_t esdFileName[300];
233   Char_t gAliceName[300];
234   Char_t clusterName[300];
235   
236   sprintf(geoFileName , "%s/geometry.root",      fReadDir.Data());
237   sprintf(esdFileName , "%s/AliESDs.root" ,      fReadDir.Data());
238   sprintf(gAliceName  , "%s/galice.root"  ,      fReadDir.Data());
239   sprintf(clusterName , "%s/MFT.RecPoints.root", fReadDir.Data());
240   
241   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
242   if (!gGeoManager) {
243     TGeoManager::Import(geoFileName);
244     if (!gGeoManager) {
245       printf("getting geometry from file %s failed", geoFileName);
246       return;
247     }
248   }
249   
250   fFileESD = new TFile(esdFileName);
251   if (!fFileESD || !fFileESD->IsOpen()) return;
252   else printf("file %s successfully opened\n", fFileESD->GetName());
253   
254   fMuonRecoCheck = new AliMUONRecoCheck(esdFileName, Form("%s/generated/", fReadDir.Data()));       // Utility class to check reconstruction
255   fFile_gAlice = new TFile(gAliceName);
256   if (!fFile_gAlice || !fFile_gAlice->IsOpen()) return;
257   else printf("file %s successfully opened\n", fFile_gAlice->GetName());
258   
259   fRunLoader = AliRunLoader::Open(gAliceName);
260   gAlice = fRunLoader->GetAliRun();
261   if (!gAlice) fRunLoader->LoadgAlice();
262   fMFT = (AliMFT*) gAlice->GetDetector("MFT"); 
263   fSegmentation = fMFT->GetSegmentation();
264   SetNPlanesMFT(fSegmentation->GetNPlanes());
265
266   if (!SetRunNumber()) return;
267   if (!InitGRP()) return;
268   AliMUONTrackExtrap::SetField();        // set the magnetic field for track extrapolations
269
270   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
271     fZPlane[iPlane]    = fSegmentation->GetPlane(iPlane)->GetZCenter();
272     fRPlaneMax[iPlane] = fSegmentation->GetPlane(iPlane)->GetRMaxSupport();
273     fRPlaneMin[iPlane] = fSegmentation->GetPlane(iPlane)->GetRMinSupport();
274   }
275   
276   // Loading MFT clusters
277   fMFTLoader = fRunLoader->GetDetectorLoader("MFT");
278   fMFTLoader->LoadRecPoints("READ");
279   
280   fMFTClusterTree = fMFTLoader->TreeR();
281
282
283   Int_t nEventsInFile = fMuonRecoCheck->NumberOfEvents();
284   if (!nEventsInFile) {
285     printf("no events available!!!\n");
286     return;
287   }
288   if (nEventsInFile<nEventsToAnalyze || nEventsToAnalyze<0) fNEventsToAnalyze = nEventsInFile;
289   else fNEventsToAnalyze = nEventsToAnalyze;
290
291   fCandidateTracks = new TClonesArray("AliMuonForwardTrack",200000);
292
293   // -------------------------- initializing histograms...
294
295   printf("\ninitializing histograms...\n");
296   BookHistos();
297   SetTitleHistos();
298   printf("... done!\n\n");
299
300   // -------------------------- initializing graphics...
301
302   printf("initializing graphics...\n");
303   BookPlanes();
304   printf("... done!\n\n");
305
306   SetSigmaSpectrometerCut(4.0);
307   SetSigmaClusterCut(4.5);
308   SetChi2GlobalCut(2.0);
309   SetNFinalCandidatesCut(10);
310   SetRAbsorberCut(26.4);
311   SetLowPtCut(0.5);
312   SetExtrapOriginTransvError(1.0);
313
314 }
315
316 //======================================================================================================================================
317
318 Bool_t AliMuonForwardTrackFinder::LoadNextEvent() {
319
320   // load next reconstructed event from the tree
321
322   if (fEv) FillOutputTree();
323
324   if (fEv>=fNEventsToAnalyze) return kFALSE;
325
326   fCountRealTracksAnalyzedOfEvent = 0;
327   
328   printf(" **** analyzing event # %d  \n", fEv);
329   
330   fTrackStore    = fMuonRecoCheck->ReconstructedTracks(fEv);
331   fTrackRefStore = fMuonRecoCheck->ReconstructibleTracks(fEv);
332   
333   fRunLoader->GetEvent(fEv);
334   if (!fMFTLoader->TreeR()->GetEvent()) return kFALSE;
335   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
336     printf("plane %02d: nClusters = %d\n", iPlane, (fMFT->GetRecPointsList(iPlane))->GetEntries());
337     fMFTClusterArray[iPlane] = fMFT->GetRecPointsList(iPlane);
338   }
339   SeparateFrontBackClusters();
340
341   fRunLoader -> LoadKinematics();
342   fStack = fRunLoader->Stack();
343   fNextTrack = fTrackStore->CreateIterator();
344   fMuonForwardTracks->Clear();
345
346   fEv++;
347   
348   return kTRUE;
349
350 }
351
352 //======================================================================================================================================
353
354 Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
355
356   fNPlanesMFTAnalyzed = 0;
357
358   // load next muon track from the reconstructed event
359
360   if (!fCountRealTracksAnalyzed) if (!LoadNextEvent()) return kFALSE;
361
362   while ( !(fMuonTrackReco = static_cast<AliMUONTrack*>(fNextTrack->Next())) ) if (!LoadNextEvent()) return kFALSE;
363
364   printf("**************************************************************************************\n");
365   printf("***************************   MUON TRACK %3d   ***************************************\n", fCountRealTracksAnalyzedOfEvent);
366   printf("**************************************************************************************\n");
367
368   fCountRealTracksAnalyzed++;
369
370   fCandidateTracks -> Clear();
371
372   fLabelMC = -1;
373   fDistanceFromGoodClusterAndTrackAtLastPlane = -1.;
374   fDistanceFromBestClusterAndTrackAtLastPlane = -1.;
375   ResetPlanes();
376
377   TIter nextTrackRef(fTrackRefStore->CreateIterator());
378   AliMUONTrack *trackRef=0;
379   
380   // --------------------------------------- loop on MC generated tracks to find the MC reference...
381   
382   while ( (trackRef = static_cast<AliMUONTrack*>(nextTrackRef())) ) {
383     // number of compatible clusters between trackReco and trackRef
384     Int_t nMatchCluster = fMuonTrackReco->FindCompatibleClusters(*trackRef, fSigmaSpectrometerCut, fIsClusterCompatible);  
385     if ( (fIsClusterCompatible[0] || fIsClusterCompatible[1] || fIsClusterCompatible[2] || fIsClusterCompatible[3]) &&   // before the dipole
386          (fIsClusterCompatible[6] || fIsClusterCompatible[7] || fIsClusterCompatible[8] || fIsClusterCompatible[9]) &&   // after the dipole
387          2*nMatchCluster>fMuonTrackReco->GetNClusters() ) {
388       fMuonTrackReco->SetMCLabel(trackRef->GetUniqueID());   // MC reference has been found for trackReco!
389       break;
390     }
391   }
392   
393   // ------------------------------------- ...done!
394
395   if (fMuonTrackReco->GetMCLabel()>=0) fCountRealTracksWithRefMC++;
396   
397   fLabelMC = fMuonTrackReco->GetMCLabel();
398
399   CheckCurrentMuonTrackable();
400
401   PrintParticleHistory();
402   
403   if (fMuonTrackReco->GetMatchTrigger()) fCountRealTracksWithRefMC_andTrigger++;
404   
405   // the track we are going to build, starting from fMuonTrackReco and adding the MFT clusters
406   AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
407   track -> SetMUONTrack(fMuonTrackReco);
408   if (fLabelMC>=0) track -> SetMCTrackRef(fStack->Particle(fLabelMC));
409   track -> SetMCLabel(fMuonTrackReco->GetMCLabel());
410   track -> SetMatchTrigger(fMuonTrackReco->GetMatchTrigger());
411   
412   // track parameters at the first tracking station in the Muon Spectrometer
413   AliMUONTrackParam *param = (AliMUONTrackParam*) (fMuonTrackReco->GetTrackParamAtCluster()->First());
414   Double_t ptSpectrometer = TMath::Sqrt(param->Px()*param->Px() + param->Py()*param->Py());
415   Double_t thetaSpectrometer = TMath::ATan(ptSpectrometer/param->Pz());
416   if (thetaSpectrometer<0.) thetaSpectrometer += TMath::Pi();
417   Double_t etaSpectrometer = -1.*TMath::Log(TMath::Tan(0.5*thetaSpectrometer));
418   fHistPtSpectrometer -> Fill(ptSpectrometer);
419   
420   // if the transverse momentum in the Muon Spectrometer is smaller than the threshold, skip to the next track
421   if (ptSpectrometer < fLowPtCut) return 3;
422   
423   // track parameters linearly extrapolated from the first tracking station to the end of the absorber
424   AliMUONTrackParam trackParamEndOfAbsorber(*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
425   AliMUONTrackExtrap::ExtrapToZCov(&trackParamEndOfAbsorber, -503.);   // absorber extends from -90 to -503 cm
426   Double_t xEndOfAbsorber = trackParamEndOfAbsorber.GetNonBendingCoor();
427   Double_t yEndOfAbsorber = trackParamEndOfAbsorber.GetBendingCoor();
428   Double_t rAbsorber      = TMath::Sqrt(xEndOfAbsorber*xEndOfAbsorber + yEndOfAbsorber*yEndOfAbsorber);
429   fHistRadiusEndOfAbsorber -> Fill(rAbsorber);
430   
431   // if the radial distance of the track at the end of the absorber is smaller than a radius corresponding to 
432   // 3 degrees as seen from the interaction point, skip to the next track
433   if (rAbsorber < fRAbsorberCut) return 4;
434   
435   //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
436   
437   for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) {          // *** do not reverse the order of this cycle!!! 
438                                                                    // *** this reflects the fact that the extrapolation is performed 
439                                                                    // *** starting from the last MFT plane back to the origin
440     
441     // --------- updating the array of tracks according to the clusters available in the i-th plane ---------
442     
443     fNPlanesMFTAnalyzed++;
444
445     if (fMatchingMode==kRealMatching) {
446       Int_t nTracksToBeAnalyzed = fCandidateTracks->GetEntriesFast();
447       for (Int_t iTrack=0; iTrack<nTracksToBeAnalyzed; iTrack++) {
448         fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
449         // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array 
450         // (several new tracks can be created for one old track)
451         FindClusterInPlane(iPlane);   
452         if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
453           fCandidateTracks->Remove(fCurrentTrack);     // the old track is removed after the check;
454         }
455       }
456       fCandidateTracks->Compress();
457       if (fIsCurrentMuonTrackable) fHistNTracksAfterExtrapolation[iPlane] -> Fill(fCandidateTracks->GetEntriesFast());
458     }
459
460     else if (fMatchingMode==kIdealMatching) {
461       fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
462       printf("plane %02d: fCandidateTracks->GetEntriesFast() = %d   fCandidateTracks->UncheckedAt(0) = %p   fCurrentTrack = %p\n", 
463              iPlane, fCandidateTracks->GetEntriesFast(), fCandidateTracks->UncheckedAt(0), fCurrentTrack);
464       AttachGoodClusterInPlane(iPlane);
465     }
466
467   }      
468   
469   // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
470   
471   if (fMatchingMode==kIdealMatching) {
472     printf("Adding track to output tree...\n");
473     AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
474     new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);  // AU
475     printf("...track added!\n");
476     fCandidateTracks->Clear();
477     fCountRealTracksAnalyzedOfEvent++;
478     return 5;
479   }
480
481   // If we have several final tracks, we must find the best candidate:
482
483   Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
484
485   if (nFinalTracks) fCountRealTracksAnalyzedWithFinalCandidates++;
486   
487   Double_t theVariable_Best        = -1.;                    // variable defining the best candidate
488   Bool_t bestCandidateExists       = kFALSE;
489   Int_t nGoodClustersBestCandidate = 0;
490   Int_t idBestCandidate            = 0;
491   Double_t chi2HistoryForBestCandidate[fMaxNPlanesMFT]={0};    // chi2 on each plane, for the best candidate
492   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = -1.;
493   
494   fTxtFinalCandidates = new TLatex(0.10, 0.78, Form("N_{FinalCandidates} = %d", nFinalTracks));
495   
496   Int_t nClustersMC = 0;
497   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) nClustersMC += fIsGoodClusterInPlane[iPlane];
498
499   for (Int_t iTrack=0; iTrack<nFinalTracks; iTrack++) {
500     
501     AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
502     
503     Double_t chi2AtPlane[fMaxNPlanesMFT]={0};
504     Int_t nGoodClusters = 0;
505     Int_t nMFTClusters  = finalTrack->GetNMFTClusters();
506     Int_t nMUONClusters = finalTrack->GetNMUONClusters();
507
508     Int_t plane = 0;
509     for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
510       while (!finalTrack->PlaneExists(plane)) plane++;
511       AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
512       chi2AtPlane[plane++] = localCluster->GetTrackChi2();
513       if (IsCorrectMatch(localCluster)) nGoodClusters++;
514       Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster);        // Muon Spectrometer clusters + clusters in the Vertex Telescope
515       Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
516       chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
517     }
518     
519     if (fIsCurrentMuonTrackable) fHistNGoodClustersForFinalTracks -> Fill(nGoodClusters);
520
521     fNtuFinalCandidates1 -> Fill(Double_t(fRun),
522                                  Double_t(fEv),
523                                  Double_t(fCountRealTracksAnalyzedOfEvent),
524                                  Double_t(nFinalTracks),
525                                  Double_t(nClustersMC),
526                                  Double_t(nGoodClusters),
527                                  ptSpectrometer,
528                                  thetaSpectrometer,
529                                  etaSpectrometer);
530
531     fNtuFinalCandidates2 -> Fill(chi2AtPlane[0],
532                                  chi2AtPlane[1],
533                                  chi2AtPlane[2],
534                                  chi2AtPlane[3],
535                                  chi2AtPlane[4],
536                                  chi2AtPlane[5],
537                                  chi2AtPlane[6],
538                                  chi2AtPlane[7],
539                                  chi2AtPlane[8],
540                                  chi2AtPlane[9],
541                                  chi2AtPlane[10],
542                                  chi2AtPlane[11],
543                                  chi2AtPlane[12],
544                                  chi2AtPlane[13],
545                                  chi2AtPlane[14]);
546     
547     // now comparing the tracks with various criteria, in order to find the best one
548     
549     Double_t theVariable = 0.;
550     for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) theVariable += chi2AtPlane[iCluster];
551     theVariable /= Double_t(nMFTClusters);
552       
553     if (theVariable<theVariable_Best || theVariable_Best<0.) {
554       nGoodClustersBestCandidate = nGoodClusters;
555       for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = chi2AtPlane[iPlane];
556       theVariable_Best = theVariable;
557       fTxtTrackFinalChi2 = new TLatex(0.20, 0.52, Form("#chi^{2}_{final} = %3.1f", chi2HistoryForBestCandidate[0]));
558       for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
559         fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1], 
560                                            0.90*fRPlaneMax[fNPlanesMFT-1], 
561                                            Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
562       }
563       idBestCandidate = iTrack;
564       bestCandidateExists=kTRUE;
565     }
566
567     // ----------------------------------------------------------
568
569   }
570   
571   if (nFinalTracks) {
572     FillPlanesWithTrackHistory((AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate));
573     AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
574     new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);   // AU
575     fNtuFinalBestCandidates1 -> Fill(Double_t(fRun),
576                                      Double_t(fEv),
577                                      Double_t(fCountRealTracksAnalyzedOfEvent),
578                                      Double_t(nFinalTracks),
579                                      Double_t(nClustersMC),
580                                      Double_t(nGoodClustersBestCandidate),
581                                      ptSpectrometer,
582                                      thetaSpectrometer,
583                                      etaSpectrometer);
584
585     fNtuFinalBestCandidates2 -> Fill(chi2HistoryForBestCandidate[0],
586                                      chi2HistoryForBestCandidate[1],
587                                      chi2HistoryForBestCandidate[2],
588                                      chi2HistoryForBestCandidate[3],
589                                      chi2HistoryForBestCandidate[4],
590                                      chi2HistoryForBestCandidate[5],
591                                      chi2HistoryForBestCandidate[6],
592                                      chi2HistoryForBestCandidate[7],
593                                      chi2HistoryForBestCandidate[8],
594                                      chi2HistoryForBestCandidate[9],
595                                      chi2HistoryForBestCandidate[10],
596                                      chi2HistoryForBestCandidate[11],
597                                      chi2HistoryForBestCandidate[12],
598                                      chi2HistoryForBestCandidate[13],
599                                      chi2HistoryForBestCandidate[14]);
600
601   }
602   
603   if (fDrawOption && bestCandidateExists) {
604     fTxtTrackGoodClusters = new TLatex(0.20, 0.59, Form("N_{GoodClusters} = %d", nGoodClustersBestCandidate));
605     DrawPlanes();
606   }
607
608   if (fIsCurrentMuonTrackable) {
609     if (nGoodClustersBestCandidate==5) fHistPtMuonTrackWithGoodMatch -> Fill(ptSpectrometer);
610     else                               fHistPtMuonTrackWithBadMatch  -> Fill(ptSpectrometer);
611   }
612
613   // -------------------------------------------------------------------------------------------
614
615   fCandidateTracks->Clear();
616   
617   fCountRealTracksAnalyzedOfEvent++;
618
619   return 5;
620   
621 }
622
623 //===========================================================================================================================================
624
625 void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) { 
626   
627   printf(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId);
628
629   // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
630
631   // propagate track to plane #planeId (both to front and back active sensors)
632   // look for compatible clusters
633   // update TrackParam at found cluster (if any) using Kalman Filter
634
635   AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
636
637   if (planeId == fNPlanesMFT-1) {      // last plane of the telecope
638     currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
639     currentParamBack  = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
640     currentParamForResearchFront = currentParamFront;
641     currentParamForResearchBack  = currentParamBack;
642     AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, 0.); 
643     AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  0.); 
644     AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, 0., 0., gRandom->Gaus(0,fGaussianBlurZVert), fExtrapOriginTransvError, fExtrapOriginTransvError); 
645     AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack,  0., 0., gRandom->Gaus(0,fGaussianBlurZVert), fExtrapOriginTransvError, fExtrapOriginTransvError); 
646   }
647   else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
648     currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
649     currentParamBack  = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
650     currentParamForResearchFront = currentParamFront;
651     currentParamForResearchBack  = currentParamBack;
652     AliMUONTrackExtrap::AddMCSEffect(&currentParamFront,           (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
653                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
654     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
655                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
656     AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,            (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
657                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
658     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
659                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
660   }
661   // for all planes: extrapolation to the Z of the plane
662   AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront,            -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());   
663   AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
664   AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack,             -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());   
665   AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchBack,  -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
666
667   //---------------------------------------------------------------------------------------
668
669   TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
670   TMatrixD covBack(5,5);  covBack  = currentParamForResearchBack.GetCovariances();
671   
672   Double_t squaredError_X_Front = covFront(0,0);
673   Double_t squaredError_Y_Front = covFront(2,2);
674   Double_t squaredError_X_Back  = covBack(0,0);
675   Double_t squaredError_Y_Back  = covBack(2,2);
676
677   Double_t corrFact = 1.0;
678
679   Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
680   Double_t researchRadiusBack  = TMath::Sqrt(squaredError_X_Back  + squaredError_Y_Back);
681   if (planeId==fNPlanesMFT-1 && 0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtLastPlane) {
682     corrFact = fMinResearchRadiusAtLastPlane/(0.5*(researchRadiusFront+researchRadiusBack));
683   }
684   if (fIsCurrentMuonTrackable) fHistResearchRadius[planeId] -> Fill(corrFact*0.5*(researchRadiusFront+researchRadiusBack));
685
686   Double_t position_X_Front = currentParamForResearchFront.GetNonBendingCoor();
687   Double_t position_Y_Front = currentParamForResearchFront.GetBendingCoor();
688   Double_t position_X_Back = currentParamForResearchBack.GetNonBendingCoor();
689   Double_t position_Y_Back = currentParamForResearchBack.GetBendingCoor();
690   Double_t radialPositionOfTrackFront = TMath::Sqrt(position_X_Front*position_X_Front + position_Y_Front*position_Y_Front);
691   Double_t radialPositionOfTrackBack  = TMath::Sqrt(position_X_Back*position_X_Back   + position_Y_Back*position_Y_Back);
692
693   //---------------------------------------------------------------------------------------
694
695   Double_t chi2cut = 2.*fSigmaClusterCut*fSigmaClusterCut;     // depends on the number of variables (here, 2)
696   
697   // Analyizing the clusters: FRONT ACTIVE ELEMENTS
698   
699   Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
700   printf("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId);
701   
702   for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
703
704     Bool_t isGoodChi2 = kFALSE;
705
706     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster); 
707     Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster);     // describes the compatibility between the track and the cluster
708     if (chi2<chi2cut) isGoodChi2 = kTRUE;
709
710     Double_t radialPositionOfClusterFront = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());    
711     if (planeId == fNPlanesMFT-1) {
712       if (TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront)<fDistanceFromBestClusterAndTrackAtLastPlane ||
713           fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
714         fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
715       }
716       if (IsCorrectMatch(cluster)) {
717         fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
718       }
719     }
720
721     if (fIsCurrentMuonTrackable) {
722       if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.);     //  chi2/ndf
723       else                         fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.);     //  chi2/ndf
724     }
725
726     if (isGoodChi2) {
727       printf("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
728       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
729       newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);    // creating new track param and attaching the cluster
730       newTrack->SetPlaneExists(planeId);
731       printf("current muon is trackable: %d\n", fIsCurrentMuonTrackable);
732       if (fIsCurrentMuonTrackable) {
733         Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
734         printf("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2());
735         Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters();        // Muon Spectrometer clusters + clusters in the Vertex Telescope
736         Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
737         if (IsCorrectMatch(cluster)) fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
738         else                         fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
739       }
740       fGrMFTPlane[planeId][kClustersGoodChi2] -> SetPoint(fGrMFTPlane[planeId][kClustersGoodChi2]->GetN(), cluster->GetX(), cluster->GetY());
741     }
742     else printf("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
743
744   }
745
746   // Analyizing the clusters: BACK ACTIVE ELEMENTS
747   
748   Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
749   printf("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId);
750   
751   for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
752
753     Bool_t isGoodChi2 = kFALSE;
754
755     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster); 
756     Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster);     // describes the compatibility between the track and the cluster
757     if (chi2<chi2cut) isGoodChi2 = kTRUE;
758
759     Double_t radialPositionOfClusterBack = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());    
760     if (planeId == fNPlanesMFT-1) {
761       if (TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack)<fDistanceFromBestClusterAndTrackAtLastPlane ||
762           fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
763         fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
764       }
765       if (IsCorrectMatch(cluster)) {
766         fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
767       }
768     }
769
770     if (fIsCurrentMuonTrackable) {
771       if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.);     //  chi2/ndf
772       else                         fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.);     //  chi2/ndf
773     }
774
775     if (isGoodChi2) {
776       printf("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
777       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
778       newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);    // creating new track param and attaching the cluster
779       newTrack->SetPlaneExists(planeId);
780       printf("current muon is trackable: %d\n", fIsCurrentMuonTrackable);
781       if (fIsCurrentMuonTrackable) {
782         Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
783         printf("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2());
784         Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters();        // Muon Spectrometer clusters + clusters in the Vertex Telescope
785         Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
786         if (IsCorrectMatch(cluster)) fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
787         else                         fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
788       }
789       fGrMFTPlane[planeId][kClustersGoodChi2] -> SetPoint(fGrMFTPlane[planeId][kClustersGoodChi2]->GetN(), cluster->GetX(), cluster->GetY());
790     }
791     else printf("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
792
793   }
794
795   //---------------------------------------------------------------------------------------------
796
797   if (planeId == fNPlanesMFT-1) {
798     if (fIsCurrentMuonTrackable && fDistanceFromGoodClusterAndTrackAtLastPlane>0.) {
799       fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Fill(TMath::Abs(fDistanceFromBestClusterAndTrackAtLastPlane-
800                                                                                                        fDistanceFromGoodClusterAndTrackAtLastPlane));
801       fHistDistanceGoodClusterFromTrackAtLastPlane -> Fill(fDistanceFromGoodClusterAndTrackAtLastPlane);
802     }
803   }
804   
805 }
806
807 //==========================================================================================================================================
808
809 void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) { 
810   
811   printf(">>>> executing AliMuonForwardTrackFinder::AttachGoodClusterInPlane(%d)\n", planeId);
812
813   AliMUONTrackParam currentParamFront, currentParamBack;
814
815   if (planeId == fNPlanesMFT-1) {      // last plane of the telecope
816     currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
817     currentParamBack  = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
818     AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, 0.); 
819     AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  0.); 
820   }
821   else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
822     printf("fCurrentTrack = %p\n", fCurrentTrack);
823     currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
824     currentParamBack  = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
825     AliMUONTrackExtrap::AddMCSEffect(&currentParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
826                                                           fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
827     AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,  (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
828                                                           fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
829   }
830   // for all planes: linear extrapolation to the Z of the plane
831   AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());   
832   AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack,  -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());   
833
834   Bool_t goodClusterFound = kFALSE;
835   
836   // Analyizing the clusters: FRONT ACTIVE ELEMENTS
837
838   Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
839   
840   printf("nClustersFront = %d\n", nClustersFront);
841   for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
842     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->UncheckedAt(iCluster);
843     printf("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersFront, cluster, fCurrentTrack);
844     if (IsCorrectMatch(cluster)) {
845       fCurrentTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);  // creating new track param and attaching the cluster
846       fCurrentTrack->SetPlaneExists(planeId);
847       goodClusterFound = kTRUE;
848       break;
849     }
850   }
851
852   if (goodClusterFound) return;
853
854   // Analyizing the clusters: BACK ACTIVE ELEMENTS
855
856   Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
857   
858   printf("nClustersBack = %d\n", nClustersBack);
859   for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
860     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->UncheckedAt(iCluster);
861     printf("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersBack, cluster, fCurrentTrack);
862     if (IsCorrectMatch(cluster)) {
863       fCurrentTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);  // creating new track param and attaching the cluster
864       fCurrentTrack->SetPlaneExists(planeId);
865       goodClusterFound = kTRUE;
866       break;
867     }
868   }
869
870 }
871
872 //==========================================================================================================================================
873
874 void AliMuonForwardTrackFinder::CheckCurrentMuonTrackable() {
875
876   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
877     fIsGoodClusterInPlane[iPlane] = kFALSE;
878     Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
879     for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
880       AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster); 
881       for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
882         if (cluster->GetMCLabel(iTrack)==fLabelMC) {
883           fIsGoodClusterInPlane[iPlane] = kTRUE;
884           break;
885         }
886       }
887     }
888   }
889
890   fIsCurrentMuonTrackable = kTRUE;
891   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) fIsCurrentMuonTrackable = (fIsCurrentMuonTrackable&&fIsGoodClusterInPlane[iPlane]);
892
893 }
894
895 //==========================================================================================================================================
896
897 void AliMuonForwardTrackFinder::FillPlanesWithTrackHistory(AliMuonForwardTrack *track) { 
898   
899   // recover track parameters on each planes and look for the corresponding clusters
900
901   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
902     
903     AliMFTCluster *trackCluster = (AliMFTCluster *) track->GetMFTCluster(iPlane);
904
905     fGrMFTPlane[iPlane][kClusterOfTrack] -> SetPoint(fGrMFTPlane[iPlane][kClusterOfTrack]->GetN(), trackCluster->GetX(), trackCluster->GetY());
906     
907     Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
908     
909     for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
910       
911       AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->UncheckedAt(iCluster); 
912
913       fGrMFTPlane[iPlane][kAllClusters] -> SetPoint(fGrMFTPlane[iPlane][kAllClusters]->GetN(), cluster->GetX(), cluster->GetY());
914       
915       if (IsCorrectMatch(cluster)) {
916         fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetPoint(fGrMFTPlane[iPlane][kClusterCorrectMC]->GetN(), cluster->GetX(), cluster->GetY());
917       }
918
919     }
920
921   }
922
923 }
924
925 //======================================================================================================================================
926
927 Bool_t AliMuonForwardTrackFinder::IsCorrectMatch(AliMFTCluster *cluster) {
928
929   Bool_t result = kFALSE;
930
931   // check if the cluster belongs to the correct MC track
932
933   for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
934     if (cluster->GetMCLabel(iTrack)==fLabelMC) {
935       result = kTRUE;
936       break;
937     }
938   }
939
940   printf("returning %d\n", result);
941
942   return result;
943
944 }
945
946 //======================================================================================================================================
947
948 Double_t AliMuonForwardTrackFinder::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
949
950   // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
951   // return the corresponding Chi2
952   // assume the track parameters are given at the Z of the cluster
953   
954   // Set differences between trackParam and cluster in the bending and non bending directions
955   Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
956   Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
957   printf("dX = %f, dY = %f\n", dX, dY);
958   
959   // Calculate errors and covariances
960   const TMatrixD& kParamCov = trackParam.GetCovariances();
961   Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
962   Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
963   printf("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2);
964   Double_t covXY   = kParamCov(0,2);
965   Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
966   
967   // Compute chi2
968   if (det==0.) return 1.e10;
969   return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
970   
971 }
972
973 //=========================================================================================================================================
974
975 void AliMuonForwardTrackFinder::SeparateFrontBackClusters() {
976
977   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
978     fMFTClusterArrayFront[iPlane]->Clear();
979     fMFTClusterArrayBack[iPlane] ->Clear();
980     for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
981       AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
982       if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
983         new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
984       }
985       else {
986         new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
987       }
988     }
989   }
990
991 }
992
993 //=========================================================================================================================================
994
995 Int_t AliMuonForwardTrackFinder::GetNDF(Int_t nClusters) {
996
997   // the same definition as in AliMUONTrack is implemented, since here we just add more clusters to the Muon track
998   
999   Int_t ndf = 2 * nClusters - 5;
1000   return (ndf > 0) ? ndf : 0;
1001
1002 }
1003
1004 //============================================================================================================================================
1005
1006 void AliMuonForwardTrackFinder::BookHistos() {
1007   
1008   const Int_t nMaxNewTracks[]  = {150,     200,   250, 600, 1000};
1009   const Double_t radiusPlane[] = {0.001, 0.010, 0.100, 5.0,  5.0};
1010
1011   fHistPtSpectrometer = new TH1D("hPtSpectrometer", "p_{T} as given by the Muon Spectrometer", 200, 0, 20.); 
1012
1013   fHistPtMuonTrackWithGoodMatch = new TH1D("fHistPtMuonTrackWithGoodMatch", "p_{T} of muon track with good match", 200, 0, 20.); 
1014   fHistPtMuonTrackWithBadMatch  = new TH1D("fHistPtMuonTrackWithBadMatch",  "p_{T} of muon track with bad match",  200, 0, 20.); 
1015
1016   fHistRadiusEndOfAbsorber = new TH1D("hRadiusEndOfAbsorber", "Track radial distance at the end of the absorber",  1000, 0, 100.); 
1017
1018   fHistNGoodClustersForFinalTracks = new TH1D("hNGoodClustersForFinalTracks", "Number of Good Clusters per Final Track", 20, -0.25, 9.75);
1019
1020   fHistDistanceGoodClusterFromTrackAtLastPlane = new TH1D("hDistanceGoodClusterFromTrackAtLastPlane",
1021                                                           "Distance of MC Good Cluster from Track in last MFT plane", 200, 0., 2.);
1022   
1023   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane = 
1024     new TH1D("hDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane",
1025              "Good Cluster distance from track - Best Cluster distance from track in last MFT plane", 200, 0., 2.);
1026   
1027   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1028     
1029     fHistNTracksAfterExtrapolation[iPlane] = new TH1D(Form("hNTracksAfterExtrapolation_pl%02d", iPlane),
1030                                                       Form("Number of Candidates after analysis of MFT plane %02d", iPlane),
1031                                                       nMaxNewTracks[iPlane], -0.5, nMaxNewTracks[iPlane]-0.5);
1032     
1033     fHistResearchRadius[iPlane] = new TH1D(Form("hResearchRadius_pl%02d", iPlane),
1034                                            Form("Research Radius for candidate clusters in MFT plane %02d", iPlane),
1035                                            1000, 0., radiusPlane[iPlane]);
1036     
1037     fHistChi2Cluster_GoodCluster[iPlane] = new TH1D(Form("hChi2Cluster_GoodCluster_pl%02d", iPlane),
1038                                                     Form("#chi^{2}_{clust} for Good clusters in MFT plane %02d", iPlane),
1039                                                     100, 0., 15.);
1040     
1041     fHistChi2Cluster_BadCluster[iPlane] = new TH1D(Form("hChi2Cluster_BadCluster_pl%02d", iPlane),
1042                                                    Form("#chi^{2}_{clust} for Bad clusters in MFT plane %02d", iPlane),
1043                                                    100, 0., 15.);
1044
1045     fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1046                                                                            Form("#chi^{2}/ndf at plane %d for GOOD candidates of trackable muons",iPlane),
1047                                                                            100, 0., 15.);
1048
1049     fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1050                                                                            Form("#chi^{2}/ndf at plane %d for BAD candidates of trackable muons",iPlane),
1051                                                                            100, 0., 15.);
1052
1053   }
1054   
1055   //------------------------------------------
1056   
1057   fHistPtSpectrometer               -> Sumw2();
1058   fHistPtMuonTrackWithGoodMatch     -> Sumw2();
1059   fHistPtMuonTrackWithBadMatch      -> Sumw2();
1060   fHistRadiusEndOfAbsorber          -> Sumw2();
1061   fHistNGoodClustersForFinalTracks  -> Sumw2();
1062
1063   fHistDistanceGoodClusterFromTrackAtLastPlane                                  -> Sumw2();  
1064   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Sumw2();  
1065
1066   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1067     
1068     fHistNTracksAfterExtrapolation[iPlane] -> Sumw2();
1069     fHistResearchRadius[iPlane]            -> Sumw2();
1070     
1071     fHistChi2Cluster_GoodCluster[iPlane]        -> Sumw2();
1072     fHistChi2Cluster_BadCluster[iPlane]         -> Sumw2();
1073
1074     fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
1075     fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> Sumw2();
1076     
1077   }
1078
1079   fNtuFinalCandidates1     = new TNtuple("ntuFinalCandidates1",     "Final Candidates (ALL)", "run:event:muonTrack:nFinalCandidates:nClustersMC:nGoodClusters:ptSpectrometer:thetaSpectrometer:etaSpectrometer");
1080
1081   fNtuFinalBestCandidates1 = new TNtuple("ntuFinalBestCandidates1", "Final Best Candidates",  "run:event:muonTrack:nFinalCandidates:nClustersMC:nGoodClusters:ptSpectrometer:thetaSpectrometer:etaSpectrometer");
1082
1083   fNtuFinalCandidates2     = new TNtuple("ntuFinalCandidates2",     "Final Candidates (ALL)", "chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:chi2AtPlane9:chi2AtPlane10:chi2AtPlane11:chi2AtPlane12:chi2AtPlane13:chi2AtPlane14");
1084
1085   fNtuFinalBestCandidates2 = new TNtuple("ntuFinalBestCandidates2", "Final Best Candidates",  "chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:chi2AtPlane9:chi2AtPlane10:chi2AtPlane11:chi2AtPlane12:chi2AtPlane13:chi2AtPlane14");
1086
1087 }
1088
1089 //============================================================================================================================================
1090
1091 void AliMuonForwardTrackFinder::SetTitleHistos() {
1092
1093   fHistPtSpectrometer              -> SetXTitle("p_{T}  [GeV/c]");
1094   fHistPtMuonTrackWithGoodMatch    -> SetXTitle("p_{T}  [GeV/c]");
1095   fHistPtMuonTrackWithBadMatch     -> SetXTitle("p_{T}  [GeV/c]");
1096   fHistRadiusEndOfAbsorber         -> SetXTitle("R_{abs}  [cm]");
1097   fHistNGoodClustersForFinalTracks -> SetXTitle("N_{GoodClusters}");
1098
1099   fHistDistanceGoodClusterFromTrackAtLastPlane                                  -> SetXTitle("Distance  [cm]");  
1100   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> SetXTitle("Distance  [cm]");  
1101
1102
1103   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1104     
1105     fHistNTracksAfterExtrapolation[iPlane] -> SetXTitle("N_{tracks}");
1106     fHistResearchRadius[iPlane]            -> SetXTitle("Research Radius  [cm]");
1107     
1108     fHistChi2Cluster_GoodCluster[iPlane]         -> SetXTitle("#chi^{2}/ndf");
1109     fHistChi2Cluster_BadCluster[iPlane]          -> SetXTitle("#chi^{2}/ndf");
1110
1111     fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
1112     fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> SetXTitle("#chi^{2}/ndf");
1113     
1114   }
1115
1116 }
1117
1118 //===========================================================================================================================================
1119
1120 void AliMuonForwardTrackFinder::BookPlanes() {
1121
1122   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1123     fGrMFTPlane[iPlane][kAllClusters] = new TGraph();
1124     fGrMFTPlane[iPlane][kAllClusters] -> SetName(Form("fGrMFTPlane_%02d_AllClusters",iPlane));
1125     fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerStyle(20);
1126     //    fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.5);
1127     //    fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.3);
1128     fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.2);
1129   }
1130
1131   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1132     fGrMFTPlane[iPlane][kClustersGoodChi2] = new TGraph();
1133     fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetName(Form("fGrMFTPlane_%02d_ClustersGoodChi2",iPlane));
1134     fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerStyle(20);
1135     //    fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.8);
1136     //    fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.4);
1137     fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.3);
1138     fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerColor(kBlue);
1139   }
1140
1141   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1142     fGrMFTPlane[iPlane][kClusterOfTrack] = new TGraph();
1143     fGrMFTPlane[iPlane][kClusterOfTrack] -> SetName(Form("fGrMFTPlane_%02d_ClustersOfTrack",iPlane));
1144     fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerStyle(25);
1145     //    fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerSize(1.2);
1146     fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerSize(0.9);
1147     fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerColor(kRed);
1148     fGrMFTPlane[iPlane][kClusterOfTrack] -> SetTitle(Form("Plane %d (%3.1f cm)", iPlane, fZPlane[iPlane]));
1149   }
1150
1151   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1152     fGrMFTPlane[iPlane][kClusterCorrectMC] = new TGraph();
1153     fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetName(Form("fGrMFTPlane_%02d_ClustersCorrectMC",iPlane));
1154     fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerStyle(20);
1155     //    fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerSize(0.8);
1156     fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerSize(0.5);
1157     fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerColor(kGreen);
1158   }
1159
1160   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1161     fCircleExt[iPlane] = new TEllipse(0., 0., fRPlaneMax[iPlane], fRPlaneMax[iPlane]);
1162     fCircleInt[iPlane] = new TEllipse(0., 0., fRPlaneMin[iPlane], fRPlaneMin[iPlane]);
1163   }
1164   
1165   fTxtDummy = new TLatex(0.10, 0.67, "Best Candidate:");
1166
1167   //---------------------------------------------------
1168
1169   fMrkAllClust = new TMarker(0.10, 0.32, 20);
1170   fMrkAllClust -> SetMarkerSize(0.5);
1171
1172   fMrkClustGoodChi2 = new TMarker(0.10, 0.26, 20);
1173   fMrkClustGoodChi2 -> SetMarkerSize(0.8);
1174   fMrkClustGoodChi2 -> SetMarkerColor(kBlue);
1175
1176   fMrkClustMC = new TMarker(0.10, 0.20, 20);
1177   fMrkClustMC -> SetMarkerSize(0.8);
1178   fMrkClustMC -> SetMarkerColor(kGreen);
1179
1180   fMrkClustOfTrack = new TMarker(0.10, 0.14, 25);
1181   fMrkClustOfTrack -> SetMarkerSize(1.2);
1182   fMrkClustOfTrack -> SetMarkerColor(kRed);
1183
1184   fTxtAllClust = new TLatex(0.15, 0.30, "All Clusters");
1185   fTxtAllClust -> SetTextSize(0.040);
1186
1187   fTxtClustGoodChi2 = new TLatex(0.15, 0.24, "Clusters involved in the research");
1188   fTxtClustGoodChi2 -> SetTextSize(0.040);
1189
1190   fTxtClustMC = new TLatex(0.15, 0.18, "MC good clusters");
1191   fTxtClustMC -> SetTextSize(0.040);
1192
1193   fTxtClustOfTrack = new TLatex(0.15, 0.12, "Clusters of the best candidate");
1194   fTxtClustOfTrack -> SetTextSize(0.040);
1195
1196 }
1197
1198 //===========================================================================================================================================
1199
1200 void AliMuonForwardTrackFinder::ResetPlanes() {
1201
1202   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1203     for (Int_t iGr=0; iGr<4; iGr++) {
1204       Int_t nOldClusters = fGrMFTPlane[iPlane][iGr]->GetN();
1205       for (Int_t iPoint=nOldClusters-1; iPoint>=0; iPoint--) fGrMFTPlane[iPlane][iGr]->RemovePoint(iPoint);
1206     }
1207   }
1208
1209 }
1210
1211 //===========================================================================================================================================
1212
1213 void AliMuonForwardTrackFinder::PrintParticleHistory() {
1214
1215   TString history = "";
1216   
1217   TParticle *part = fStack->Particle(fLabelMC);
1218   
1219   if (part->GetFirstMother() != -1) {
1220     TParticle *partMother = fStack->Particle(part->GetFirstMother());
1221     if (partMother->GetFirstMother() != -1) history += "...  #rightarrow ";
1222     Char_t newName[100];
1223     PDGNameConverter(partMother->GetName(), newName);
1224     history += Form("%s #rightarrow ", newName);
1225   }
1226   Char_t newName[100];
1227   PDGNameConverter(part->GetName(), newName);
1228   history += Form("%s  at  z = %5.1f cm", newName, part->Vz());
1229   
1230   //  printf("%s", history.Data());
1231   
1232   fTxtMuonHistory = new TLatex(0.10, 0.86, history.Data());
1233   
1234 }
1235
1236 //===========================================================================================================================================
1237
1238 Bool_t AliMuonForwardTrackFinder::IsMother(const Char_t *nameMother) {
1239   
1240   Bool_t result = kFALSE;
1241   
1242   TParticle *part = fStack->Particle(fLabelMC);
1243   
1244   if (part->GetFirstMother() != -1) {
1245     TParticle *partMother = fStack->Particle(part->GetFirstMother());
1246     if (!strcmp(partMother->GetName(), nameMother)) result=kTRUE;
1247   }
1248
1249   return result;
1250
1251 }
1252
1253 //===========================================================================================================================================
1254
1255 void AliMuonForwardTrackFinder::DrawPlanes() {
1256
1257   fCanvas -> Clear();
1258   if (fNPlanesMFT <= 5)       fCanvas -> Divide(3,2);
1259   else if (fNPlanesMFT <= 11) fCanvas -> Divide(4,3);
1260   else if (fNPlanesMFT <= 19) fCanvas -> Divide(5,4);
1261
1262   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1263     
1264     fCanvas->cd(fNPlanesMFT-iPlane+1);
1265     
1266     fGrMFTPlane[iPlane][kClusterOfTrack] -> GetXaxis() -> SetLimits(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1267     fGrMFTPlane[iPlane][kClusterOfTrack] -> GetYaxis() -> SetRangeUser(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1268     fGrMFTPlane[iPlane][kClusterOfTrack] -> GetXaxis() -> SetTitle("X  [cm]");
1269     fGrMFTPlane[iPlane][kClusterOfTrack] -> GetYaxis() -> SetTitle("Y  [cm]");
1270     fGrMFTPlane[iPlane][kClusterOfTrack] -> Draw("ap");
1271
1272     fCircleExt[iPlane] -> Draw("same");
1273     fCircleInt[iPlane] -> Draw("same");
1274     
1275     if (fGrMFTPlane[iPlane][kAllClusters]->GetN())       fGrMFTPlane[iPlane][kAllClusters]      -> Draw("psame");
1276     if (fGrMFTPlane[iPlane][kClustersGoodChi2]->GetN())  fGrMFTPlane[iPlane][kClustersGoodChi2] -> Draw("psame");
1277     if (fGrMFTPlane[iPlane][kClusterOfTrack]->GetN())    fGrMFTPlane[iPlane][kClusterOfTrack]   -> Draw("psame");
1278     if (fGrMFTPlane[iPlane][kClusterCorrectMC]->GetN())  fGrMFTPlane[iPlane][kClusterCorrectMC] -> Draw("psame");
1279
1280     fTxtTrackChi2[iPlane] -> Draw("same");
1281
1282   }
1283
1284   fCanvas -> cd(1);
1285   fTxtMuonHistory       -> Draw();
1286   fTxtDummy             -> Draw("same");
1287   fTxtTrackGoodClusters -> Draw("same");
1288   fTxtTrackFinalChi2    -> Draw("same");
1289   fTxtFinalCandidates   -> Draw("same");
1290
1291   fMrkAllClust      -> Draw("same");
1292   fMrkClustGoodChi2 -> Draw("same");
1293   fMrkClustMC       -> Draw("same");
1294   fMrkClustOfTrack  -> Draw("same");
1295
1296   fTxtAllClust      -> Draw("same");
1297   fTxtClustGoodChi2 -> Draw("same");
1298   fTxtClustMC       -> Draw("same");
1299   fTxtClustOfTrack  -> Draw("same");
1300
1301   //  fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1302   fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1303   if (IsMother("phi")) {
1304     fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1305     fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1306   }
1307   if (IsMother("J/psi")) {
1308     fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1309     fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1310   }
1311     
1312 }
1313
1314 //===========================================================================================================================================
1315
1316 void AliMuonForwardTrackFinder::Terminate() {
1317   
1318   printf("\n");
1319   printf("---------------------------------------------------------------------------------------------------------------\n");
1320   printf("%8d  tracks analyzed\n",                                                     fCountRealTracksAnalyzed);
1321   printf("%8d  tracks with MC ref\n",                                                  fCountRealTracksWithRefMC);
1322   printf("%8d  tracks with MC ref & trigger match\n",                                  fCountRealTracksWithRefMC_andTrigger);
1323   printf("%8d  tracks analyzed with final candidates\n",                               fCountRealTracksAnalyzedWithFinalCandidates);
1324 //   printf("%8d  tracks with MC ref & trigger match & pt>%3.1f GeV/c\n",                 fCountRealTracksWithRefMC_andTrigger_andGoodPt, fLowPtCut);
1325 //   printf("%8d  tracks with MC ref & trigger match & pt>%3.1f GeV/c & correct R_abs\n", fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta, fLowPtCut);
1326   printf("---------------------------------------------------------------------------------------------------------------\n");
1327
1328   WriteOutputTree();
1329   WriteHistos();
1330
1331 }
1332
1333 //==========================================================================================================================================
1334
1335 void AliMuonForwardTrackFinder::FillOutputTree() {
1336
1337   if (!fMuonForwardTracks || !fOutputEventTree) return;
1338
1339   AliDebug(1, Form("Filling output tree %p with %p having %d entries whose 1st entry is %p", 
1340                    fOutputEventTree, fMuonForwardTracks, fMuonForwardTracks->GetEntries(), fMuonForwardTracks->At(0)));
1341   
1342   fOutputEventTree->Fill();
1343   AliDebug(1, Form("\n\nFilled Tree: nEvents = %d!!!!\n\n", Int_t(fOutputEventTree->GetEntries())));
1344
1345 }
1346
1347 //==========================================================================================================================================
1348
1349 void AliMuonForwardTrackFinder::WriteOutputTree() {
1350
1351   if (!fOutputEventTree || !fOutputTreeFile) return;
1352
1353   fOutputTreeFile -> cd();
1354
1355   fOutputEventTree -> Write();
1356   fOutputTreeFile -> Close();
1357
1358 }
1359
1360 //==========================================================================================================================================
1361
1362 void AliMuonForwardTrackFinder::WriteHistos() {
1363
1364   TFile *fileOut = new TFile(Form("%s/MuonGlobalTracking.QA.run%d.root", fOutDir.Data(), fRun), "recreate");
1365
1366   fHistPtSpectrometer              -> Write();
1367   fHistPtMuonTrackWithGoodMatch    -> Write();
1368   fHistPtMuonTrackWithBadMatch     -> Write();
1369   fHistRadiusEndOfAbsorber         -> Write();
1370   fHistNGoodClustersForFinalTracks -> Write();
1371
1372   fHistDistanceGoodClusterFromTrackAtLastPlane                                  -> Write();  
1373   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Write();  
1374
1375   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1376     
1377     fHistNTracksAfterExtrapolation[iPlane] -> Write();
1378     fHistResearchRadius[iPlane]            -> Write();
1379     
1380     fHistChi2Cluster_GoodCluster[iPlane]   -> Write();
1381     fHistChi2Cluster_BadCluster[iPlane]    -> Write();
1382     
1383     fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Write();
1384     fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> Write();
1385
1386   }
1387
1388   fNtuFinalCandidates1     -> Write();
1389   fNtuFinalBestCandidates1 -> Write();
1390   fNtuFinalCandidates2     -> Write();
1391   fNtuFinalBestCandidates2 -> Write();
1392
1393   fileOut -> Close();
1394
1395 }
1396
1397 //===========================================================================================================================================
1398
1399 void AliMuonForwardTrackFinder::PDGNameConverter(const Char_t *nameIn, Char_t *nameOut) {
1400
1401   if      (!strcmp(nameIn, "mu+"))     sprintf(nameOut, "#mu^{+}");
1402   else if (!strcmp(nameIn, "mu-"))     sprintf(nameOut, "#mu^{-}");
1403   else if (!strcmp(nameIn, "pi+"))     sprintf(nameOut, "#pi^{+}");
1404   else if (!strcmp(nameIn, "pi-"))     sprintf(nameOut, "#pi^{-}");
1405   else if (!strcmp(nameIn, "K+"))      sprintf(nameOut, "K^{+}");
1406   else if (!strcmp(nameIn, "K-"))      sprintf(nameOut, "K^{-}");
1407   else if (!strcmp(nameIn, "K*+"))     sprintf(nameOut, "K^{*+}");
1408   else if (!strcmp(nameIn, "K*-"))     sprintf(nameOut, "K^{*-}");
1409   else if (!strcmp(nameIn, "K_S0"))    sprintf(nameOut, "K_{S}^{0}");
1410   else if (!strcmp(nameIn, "K_L0"))    sprintf(nameOut, "K_{L}^{0}");
1411   else if (!strcmp(nameIn, "K0"))      sprintf(nameOut, "K^{0}");
1412   else if (!strcmp(nameIn, "K0_bar"))  sprintf(nameOut, "#bar{K}^{0}");
1413   else if (!strcmp(nameIn, "K*0"))     sprintf(nameOut, "K^{*0}");
1414   else if (!strcmp(nameIn, "K*0_bar")) sprintf(nameOut, "#bar{K}^{*0}");
1415   else if (!strcmp(nameIn, "rho0"))    sprintf(nameOut, "#rho^{0}");
1416   else if (!strcmp(nameIn, "rho+"))    sprintf(nameOut, "#rho^{+}");
1417   else if (!strcmp(nameIn, "rho-"))    sprintf(nameOut, "#rho^{-}");
1418   else if (!strcmp(nameIn, "omega"))   sprintf(nameOut, "#omega");
1419   else if (!strcmp(nameIn, "eta'"))    sprintf(nameOut, "#eta'");
1420   else if (!strcmp(nameIn, "phi"))     sprintf(nameOut, "#phi");
1421
1422   else if (!strcmp(nameIn, "D-"))     sprintf(nameOut, "D^{-}");
1423   else if (!strcmp(nameIn, "D+"))     sprintf(nameOut, "D^{+}");
1424   else if (!strcmp(nameIn, "D0"))     sprintf(nameOut, "D^{0}");
1425   else if (!strcmp(nameIn, "D0_bar")) sprintf(nameOut, "#bar{D}^{0}");
1426   else if (!strcmp(nameIn, "D*-"))    sprintf(nameOut, "D^{*-}");
1427   else if (!strcmp(nameIn, "D*+"))    sprintf(nameOut, "D^{*+}");
1428   else if (!strcmp(nameIn, "D_s+"))   sprintf(nameOut, "D_{s}^{+}");
1429   else if (!strcmp(nameIn, "D*_s+"))  sprintf(nameOut, "D_{s}^{*+}");
1430
1431   else if (!strcmp(nameIn, "B-"))       sprintf(nameOut, "B^{-}");
1432   else if (!strcmp(nameIn, "B+"))       sprintf(nameOut, "B^{+}");
1433   else if (!strcmp(nameIn, "B_s0_bar")) sprintf(nameOut, "#bar{B}_{s}^{0}");
1434
1435   else if (!strcmp(nameIn, "antiproton"))  sprintf(nameOut, "#bar{p}");
1436   else if (!strcmp(nameIn, "proton"))      sprintf(nameOut, "p");
1437   else if (!strcmp(nameIn, "neutron"))     sprintf(nameOut, "n");
1438   else if (!strcmp(nameIn, "Sigma+"))      sprintf(nameOut, "#Sigma^{+}");
1439   else if (!strcmp(nameIn, "Delta+"))      sprintf(nameOut, "#Delta{+}");
1440   else if (!strcmp(nameIn, "Delta--"))     sprintf(nameOut, "#Delta{--}");
1441   else if (!strcmp(nameIn, "Lambda0"))     sprintf(nameOut, "#Lambda_0");
1442   else if (!strcmp(nameIn, "Lambda0_bar")) sprintf(nameOut, "#bar{Lambda}_0");
1443
1444   else sprintf(nameOut, "%s", nameIn);
1445
1446 }
1447
1448 //===========================================================================================================================================
1449
1450 void AliMuonForwardTrackFinder::SetDraw(Bool_t drawOption) { 
1451
1452   fDrawOption = drawOption;
1453
1454   if (!fCanvas) {
1455     fCanvas = new TCanvas("tracking", "tracking", 1200, 800);
1456     fCanvas -> Divide(3,2);
1457   }
1458
1459 }
1460
1461 //===========================================================================================================================================
1462
1463 Bool_t AliMuonForwardTrackFinder::InitGRP() {
1464
1465   //------------------------------------
1466   // Initialization of the GRP entry 
1467   //------------------------------------
1468
1469   AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
1470
1471   if (entry) {
1472
1473     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
1474
1475     if (m) {
1476        AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
1477        m->Print();
1478        fGRPData = new AliGRPObject();
1479        fGRPData->ReadValuesFromMap(m);
1480     }
1481
1482     else {
1483        AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it");
1484        fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
1485        entry->SetOwner(0);
1486     }
1487
1488     //    FIX ME: The unloading of GRP entry is temporarily disabled
1489     //    because ZDC and VZERO are using it in order to initialize
1490     //    their reconstructor objects. In the future one has to think
1491     //    of propagating AliRunInfo to the reconstructors.
1492     //    AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
1493   }
1494
1495   if (!fGRPData) {
1496      AliError("No GRP entry found in OCDB!");
1497      return kFALSE;
1498   }
1499
1500   TString lhcState = fGRPData->GetLHCState();
1501   if (lhcState==AliGRPObject::GetInvalidString()) {
1502     AliError("GRP/GRP/Data entry:  missing value for the LHC state ! Using UNKNOWN");
1503     lhcState = "UNKNOWN";
1504   }
1505
1506   TString beamType = fGRPData->GetBeamType();
1507   if (beamType==AliGRPObject::GetInvalidString()) {
1508     AliError("GRP/GRP/Data entry:  missing value for the beam type ! Using UNKNOWN");
1509     beamType = "UNKNOWN";
1510   }
1511
1512   Float_t beamEnergy = fGRPData->GetBeamEnergy();
1513   if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
1514     AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
1515     beamEnergy = 0;
1516   }
1517
1518   TString runType = fGRPData->GetRunType();
1519   if (runType==AliGRPObject::GetInvalidString()) {
1520     AliError("GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN");
1521     runType = "UNKNOWN";
1522   }
1523
1524   Int_t activeDetectors = fGRPData->GetDetectorMask();
1525   if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
1526     AliError("GRP/GRP/Data entry:  missing value for the detector mask ! Using 1074790399");
1527     activeDetectors = 1074790399;
1528   }
1529   AliDebug(1, Form("activeDetectors = %d", activeDetectors));
1530
1531   fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
1532   fRunInfo->Dump();
1533
1534   // *** Dealing with the magnetic field map
1535
1536   if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
1537     if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
1538       AliInfo("ExpertMode!!! GRP information will be ignored !");
1539       AliInfo("ExpertMode!!! Running with the externally locked B field !");
1540     }
1541     else {
1542       AliInfo("Destroying existing B field instance!");
1543       delete TGeoGlobalMagField::Instance();
1544     }    
1545   }
1546   if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
1547     // Construct the field map out of the information retrieved from GRP.
1548     Bool_t ok = kTRUE;
1549     // L3
1550     Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
1551     if (l3Current == AliGRPObject::GetInvalidFloat()) {
1552       AliError("GRP/GRP/Data entry:  missing value for the L3 current !");
1553       ok = kFALSE;
1554     }
1555     
1556     Char_t l3Polarity = fGRPData->GetL3Polarity();
1557     if (l3Polarity == AliGRPObject::GetInvalidChar()) {
1558       AliError("GRP/GRP/Data entry:  missing value for the L3 polarity !");
1559       ok = kFALSE;
1560     }
1561
1562     // Dipole
1563     Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
1564     if (diCurrent == AliGRPObject::GetInvalidFloat()) {
1565       AliError("GRP/GRP/Data entry:  missing value for the dipole current !");
1566       ok = kFALSE;
1567     }
1568
1569     Char_t diPolarity = fGRPData->GetDipolePolarity();
1570     if (diPolarity == AliGRPObject::GetInvalidChar()) {
1571       AliError("GRP/GRP/Data entry:  missing value for the dipole polarity !");
1572       ok = kFALSE;
1573     }
1574
1575     // read special bits for the polarity convention and map type
1576     Int_t  polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
1577     Bool_t uniformB = fGRPData->IsUniformBMap();
1578
1579     if (ok) { 
1580       AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1), 
1581                                              TMath::Abs(diCurrent) * (diPolarity ? -1:1), 
1582                                              polConvention,uniformB,beamEnergy, beamType.Data());
1583       if (fld) {
1584         TGeoGlobalMagField::Instance()->SetField( fld );
1585         TGeoGlobalMagField::Instance()->Lock();
1586         AliInfo("Running with the B field constructed out of GRP !");
1587       }
1588       else AliFatal("Failed to create a B field map !");
1589     }
1590     else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
1591   }
1592   
1593   return kTRUE;
1594
1595
1596 //====================================================================================================================================================
1597
1598 Bool_t AliMuonForwardTrackFinder::SetRunNumber() {
1599
1600   AliCDBManager *man = AliCDBManager::Instance();
1601
1602   if (!fRunLoader) {
1603     AliError("No run loader found!");
1604     return kFALSE;
1605   }
1606   else {
1607     fRunLoader->LoadHeader();
1608     // read run number from gAlice
1609     if (fRunLoader->GetHeader()) {
1610       man->SetRun(fRunLoader->GetHeader()->GetRun());
1611       fRunLoader->UnloadHeader();
1612     }
1613     else {
1614       AliError("No run-loader header found!");
1615       return kFALSE;
1616     }
1617   }
1618
1619   return kTRUE;
1620
1621 }
1622
1623 //====================================================================================================================================================
1624