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