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   if (fMuonTrackReco->GetMCLabel()>=0) fCountRealTracksWithRefMC++;
 495
 496   fLabelMC = fMuonTrackReco->GetMCLabel();
 497
 498   CheckCurrentMuonTrackable();
 499
 500   if (fMuonTrackReco->GetMatchTrigger()) fCountRealTracksWithRefMC_andTrigger++;
 501
 502   // the track we are going to build, starting from fMuonTrackReco and adding the MFT clusters
 503   AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
 504   track -> SetMUONTrack(fMuonTrackReco);
 505   if (fLabelMC>=0 && fStack->Particle(fLabelMC)) track->SetMCTrackRef(fStack->Particle(fLabelMC));
 506   track -> SetMCLabel(fMuonTrackReco->GetMCLabel());
 507   track -> SetMatchTrigger(fMuonTrackReco->GetMatchTrigger());
 508
 509   // track origin
 510   Double_t xVtx=-999., yVtx=-999., zVtx=-999999.;
 511   if (track->GetMCTrackRef()) {
 512     xVtx = track->GetMCTrackRef()->Vx();
 513     yVtx = track->GetMCTrackRef()->Vy();
 514     zVtx = track->GetMCTrackRef()->Vz();
 515   }
 516
 517   // track kinematics
 518   Double_t pt=-999., theta=-999., eta=-999.;
 519   if (track->GetMCTrackRef()) {
 520     pt    = track->GetMCTrackRef()->Pt();
 521     theta = track->GetMCTrackRef()->Theta();
 522     if (theta<0.) theta += TMath::Pi();
 523     eta   = track->GetMCTrackRef()->Eta();
 524   }
 525   else {
 526     AliMUONTrackParam *param = (AliMUONTrackParam*) (fMuonTrackReco->GetTrackParamAtCluster()->First());
 527     pt    = TMath::Sqrt(param->Px()*param->Px() + param->Py()*param->Py());
 528     theta = TMath::ATan(pt/param->Pz());
 529     if (theta<0.) theta += TMath::Pi();
 530     eta   = -1.*TMath::Log(TMath::Tan(0.5*theta));
 531   }
 532   // if the transverse momentum is smaller than the threshold, skip to the next track
 533   if (pt < fLowPtCut) return 3;
 534
 535   // track parameters linearly extrapolated from the first tracking station to the end of the absorber
 536   AliMUONTrackParam trackParamEndOfAbsorber(*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
 537   AliMUONTrackExtrap::ExtrapToZCov(&trackParamEndOfAbsorber, -503.);   // absorber extends from -90 to -503 cm
 538   Double_t xEndOfAbsorber = trackParamEndOfAbsorber.GetNonBendingCoor();
 539   Double_t yEndOfAbsorber = trackParamEndOfAbsorber.GetBendingCoor();
 540   Double_t rAbsorber      = TMath::Sqrt(xEndOfAbsorber*xEndOfAbsorber + yEndOfAbsorber*yEndOfAbsorber);
 541   fHistRadiusEndOfAbsorber -> Fill(rAbsorber);
 542
 543   // if the radial distance of the track at the end of the absorber is smaller than a given radius, skip to the next track
 544   if (rAbsorber < fRAbsorberCut) return 4;
 545
 546   //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
 547
 548   for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) {          // *** do not reverse the order of this cycle!!!
 549                                                                    // *** this reflects the fact that the extrapolation is performed
 550                                                                    // *** starting from the last MFT plane back to the origin
 551
 552     // --------- updating the array of tracks according to the clusters available in the i-th plane ---------
 553
 554     fNPlanesMFTAnalyzed++;
 555
 556     if (fMatchingMode==kRealMatching) {
 557       Int_t nTracksToBeAnalyzed = fCandidateTracks->GetEntriesFast();
 558       for (Int_t iTrack=0; iTrack<nTracksToBeAnalyzed; iTrack++) {
 559         fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
 560         // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array
 561         // (several new tracks can be created for one old track)
 562         FindClusterInPlane(iPlane);
 563         if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
 564           fCandidateTracks->Remove(fCurrentTrack);     // the old track is removed after the check;
 565         }
 566       }
 567       fCandidateTracks->Compress();
 568       if (fIsCurrentMuonTrackable) {
 569         //      fOutputQAFile->cd();
 570         fHistNTracksAfterExtrapolation[iPlane] -> Fill(fCandidateTracks->GetEntriesFast());
 571       }
 572     }
 573
 574     else if (fMatchingMode==kIdealMatching && fIsCurrentMuonTrackable) {
 575       fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
 576       AliDebug(2, Form("plane %02d: fCandidateTracks->GetEntriesFast() = %d   fCandidateTracks->UncheckedAt(0) = %p   fCurrentTrack = %p\n",
 577                        iPlane, fCandidateTracks->GetEntriesFast(), fCandidateTracks->UncheckedAt(0), fCurrentTrack));
 578       AttachGoodClusterInPlane(iPlane);
 579     }
 580
 581   }
 582
 583   // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
 584
 585   AliDebug(1, "Finished cycle over planes");
 586
 587   Double_t momentum = pt * TMath::CosH(eta);
 588   fTxtTrackMomentum = new TLatex(0.10, 0.70, Form("P_{spectro} = %3.1f GeV/c", momentum));
 589
 590   if (fMatchingMode==kIdealMatching) {
 591     AliDebug(1, "Adding track to output tree...\n");
 592     fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
 593     AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
 594     new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
 595     AliDebug(1, "...track added!\n");
 596     fCandidateTracks->Clear();
 597     fCountRealTracksAnalyzedOfEvent++;
 598     fCountRealTracksAnalyzedWithFinalCandidates++;
 599     PrintParticleHistory();
 600     FillPlanesWithTrackHistory();
 601
 602     Double_t chi2AtPlane[fNMaxPlanes] = {0};
 603     Int_t nGoodClusters = 0;
 604     Int_t nMFTClusters  = fFinalBestCandidate->GetNMFTClusters();
 605 //     Int_t nMUONClusters = fFinalBestCandidate->GetNMUONClusters();
 606     Int_t plane = 0;
 607     for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
 608       while (!fFinalBestCandidate->PlaneExists(plane)) plane++;
 609       AliMFTCluster *localCluster = fFinalBestCandidate->GetMFTCluster(iCluster);
 610       chi2AtPlane[plane] = localCluster->GetLocalChi2();
 611       if (IsCorrectMatch(localCluster)) nGoodClusters++;
 612 //       Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster);        // Muon Spectrometer clusters + clusters in the Vertex Telescope
 613 //       Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
 614 //       chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
 615       plane++;
 616     }
 617     for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 618       fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1],
 619                                          0.90*fRPlaneMax[fNPlanesMFT-1],
 620                                          Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
 621     }
 622     fTxtTrackFinalChi2 = new TLatex(0.20, 0.44, Form("#chi^{2}_{final} = %3.1f", chi2AtPlane[0]));
 623
 624     if (fDrawOption) DrawPlanes();
 625     return 5;
 626   }
 627
 628   // If we have several final tracks, we must find the best candidate:
 629
 630   Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
 631   AliDebug(1, Form("nFinalTracks = %d", nFinalTracks));
 632
 633   if (nFinalTracks) fCountRealTracksAnalyzedWithFinalCandidates++;
 634
 635   Double_t theVariable_Best        = -1.;                    // variable defining the best candidate
 636   Bool_t bestCandidateExists       = kFALSE;
 637   Int_t nGoodClustersBestCandidate = 0;
 638   Int_t idBestCandidate            = 0;
 639   Double_t chi2HistoryForBestCandidate[fNMaxPlanes] = {0};  // chi2 on each plane, for the best candidate
 640   Double_t nClustersPerPlane[fNMaxPlanes] = {0};
 641   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 642     chi2HistoryForBestCandidate[iPlane] = -1.;
 643     nClustersPerPlane[iPlane] = fMFTClusterArray[iPlane] -> GetEntries();
 644   }
 645
 646   fTxtFinalCandidates = new TLatex(0.10, 0.78, Form("N_{FinalCandidates} = %d", nFinalTracks));
 647
 648   Int_t nClustersMC = 0;
 649   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) nClustersMC += fIsGoodClusterInPlane[iPlane];
 650
 651   for (Int_t iTrack=0; iTrack<nFinalTracks; iTrack++) {
 652
 653     AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
 654
 655     Double_t chi2AtPlane[fNMaxPlanes] = {0};
 656     Int_t nGoodClusters = 0;
 657     Int_t nMFTClusters  = finalTrack->GetNMFTClusters();
 658 //     Int_t nMUONClusters = finalTrack->GetNMUONClusters();
 659
 660     Int_t plane = 0;
 661     for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
 662       while (!finalTrack->PlaneExists(plane)) plane++;
 663       AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
 664       chi2AtPlane[plane] = localCluster->GetLocalChi2();
 665       if (IsCorrectMatch(localCluster)) nGoodClusters++;
 666 //       Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster);        // Muon Spectrometer clusters + clusters in the Vertex Telescope
 667 //       Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
 668 //       chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
 669       plane++;
 670     }
 671
 672     if (fIsCurrentMuonTrackable) {
 673       //      fOutputQAFile->cd();
 674       fHistNGoodClustersForFinalTracks -> Fill(nGoodClusters);
 675     }
 676
 677     //    fOutputQAFile->cd();
 678
 679     Float_t finalCandidatesInfo[] = {Double_t(fRun),
 680                                      Double_t(fEv),
 681                                      Double_t(fCountRealTracksAnalyzedOfEvent),
 682                                      Double_t(nFinalTracks),
 683                                      Double_t(fLabelMC>=0),
 684                                      xVtx, yVtx, zVtx,
 685                                      Double_t(fMuonTrackReco->GetMatchTrigger()),
 686                                      Double_t(nClustersMC),
 687                                      Double_t(nGoodClusters),
 688                                      pt, theta, eta,
 689                                      chi2AtPlane[0],
 690                                      chi2AtPlane[1],
 691                                      chi2AtPlane[2],
 692                                      chi2AtPlane[3],
 693                                      chi2AtPlane[4],
 694                                      chi2AtPlane[5],
 695                                      chi2AtPlane[6],
 696                                      chi2AtPlane[7],
 697                                      chi2AtPlane[8]};
 698
 699     fNtuFinalCandidates -> Fill(finalCandidatesInfo);
 700
 701     // now comparing the tracks with various criteria, in order to find the best one
 702
 703     Double_t theVariable = 0.;
 704 //     theVariable = chi2AtPlane[0];
 705     for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) theVariable += chi2AtPlane[iCluster];
 706     theVariable /= Double_t(nMFTClusters);
 707
 708
 709     if (theVariable<theVariable_Best || theVariable_Best<0.) {
 710       nGoodClustersBestCandidate = nGoodClusters;
 711       for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = chi2AtPlane[iPlane];
 712       theVariable_Best = theVariable;
 713       fTxtTrackFinalChi2 = new TLatex(0.20, 0.44, Form("#chi^{2}_{final} = %3.1f", chi2HistoryForBestCandidate[0]));
 714       for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 715         fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1],
 716                                            0.90*fRPlaneMax[fNPlanesMFT-1],
 717                                            Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
 718       }
 719       idBestCandidate = iTrack;
 720       bestCandidateExists=kTRUE;
 721     }
 722
 723     // ----------------------------------------------------------
 724
 725   }
 726
 727   if (nFinalTracks) {
 728     fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
 729     AliInfo(Form("fFinalBestCandidate->GetNMFTClusters() = %d\n",  fFinalBestCandidate->GetNMFTClusters()));
 730     PrintParticleHistory();
 731     FillPlanesWithTrackHistory();
 732     AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
 733     newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
 734     new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
 735   }
 736
 737   //  fOutputQAFile->cd();
 738
 739   Float_t finalBestCandidatesInfo[] = {Double_t(fRun),
 740                                        Double_t(fEv),
 741                                        Double_t(fCountRealTracksAnalyzedOfEvent),
 742                                        Double_t(nFinalTracks),
 743                                        Double_t(fLabelMC>=0),
 744                                        xVtx, yVtx, zVtx,
 745                                        Double_t(fMuonTrackReco->GetMatchTrigger()),
 746                                        Double_t(nClustersMC),
 747                                        Double_t(nGoodClustersBestCandidate),
 748                                        pt, theta, eta,
 749                                        chi2HistoryForBestCandidate[0],
 750                                        chi2HistoryForBestCandidate[1],
 751                                        chi2HistoryForBestCandidate[2],
 752                                        chi2HistoryForBestCandidate[3],
 753                                        chi2HistoryForBestCandidate[4],
 754                                        chi2HistoryForBestCandidate[5],
 755                                        chi2HistoryForBestCandidate[6],
 756                                        chi2HistoryForBestCandidate[7],
 757                                        chi2HistoryForBestCandidate[8],
 758                                        nClustersPerPlane[0],
 759                                        nClustersPerPlane[1],
 760                                        nClustersPerPlane[2],
 761                                        nClustersPerPlane[3],
 762                                        nClustersPerPlane[4],
 763                                        nClustersPerPlane[5],
 764                                        nClustersPerPlane[6],
 765                                        nClustersPerPlane[7],
 766                                        nClustersPerPlane[8]};
 767
 768   fNtuFinalBestCandidates -> Fill(finalBestCandidatesInfo);
 769
 770   if (fDrawOption && bestCandidateExists) {
 771     fTxtTrackGoodClusters = new TLatex(0.20, 0.51, Form("N_{GoodClusters} = %d", nGoodClustersBestCandidate));
 772     DrawPlanes();
 773   }
 774
 775   // -------------------------------------------------------------------------------------------
 776
 777   fCandidateTracks->Clear();
 778   fFinalBestCandidate = NULL;
 779
 780   fCountRealTracksAnalyzedOfEvent++;
 781
 782   return 5;
 783
 784 }
 785
 786 //===========================================================================================================================================
 787
 788 void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
 789
 790   AliDebug(2, Form(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId));
 791
 792   // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
 793
 794   // propagate track to plane #planeId (both to front and back active sensors)
 795   // look for compatible clusters
 796   // update TrackParam at found cluster (if any) using Kalman Filter
 797
 798   AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
 799
 800   if (planeId == fNPlanesMFT-1) {      // last plane of the telecope
 801     currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
 802     currentParamBack  = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
 803     currentParamForResearchFront = currentParamFront;
 804     currentParamForResearchBack  = currentParamBack;
 805     Double_t xExtrap = gRandom->Gaus(0,fVertexErrorX);
 806     Double_t yExtrap = gRandom->Gaus(0,fVertexErrorY);
 807     Double_t zExtrap = gRandom->Gaus(0,fVertexErrorZ);
 808     AliMUONTrackExtrap::ExtrapToVertex(&currentParamFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
 809     AliMUONTrackExtrap::ExtrapToVertex(&currentParamBack,  xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
 810     AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
 811     AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack,  xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
 812   }
 813   else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
 814     currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
 815     currentParamBack  = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
 816     currentParamForResearchFront = currentParamFront;
 817     currentParamForResearchBack  = currentParamBack;
 818     AliMUONTrackExtrap::AddMCSEffect(&currentParamFront,           (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
 819                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
 820     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
 821                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
 822     AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,            (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
 823                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
 824     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
 825                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
 826   }
 827   // for all planes: extrapolation to the Z of the plane
 828   AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront,            -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
 829   AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
 830   AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack,             -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
 831   AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchBack,  -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
 832
 833   //---------------------------------------------------------------------------------------
 834
 835   TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
 836   TMatrixD covBack(5,5);  covBack  = currentParamForResearchBack.GetCovariances();
 837
 838   Double_t squaredError_X_Front = covFront(0,0);
 839   Double_t squaredError_Y_Front = covFront(2,2);
 840   Double_t squaredError_X_Back  = covBack(0,0);
 841   Double_t squaredError_Y_Back  = covBack(2,2);
 842
 843   Double_t corrFact = 1.0;
 844
 845   Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
 846   Double_t researchRadiusBack  = TMath::Sqrt(squaredError_X_Back  + squaredError_Y_Back);
 847   if (planeId==fNPlanesMFT-1 && 0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtLastPlane) {
 848     corrFact = fMinResearchRadiusAtLastPlane/(0.5*(researchRadiusFront+researchRadiusBack));
 849   }
 850   if (fIsCurrentMuonTrackable) {
 851     //    fOutputQAFile->cd();
 852     fHistResearchRadius[planeId] -> Fill(0.5*(researchRadiusFront+researchRadiusBack));
 853   }
 854
 855   Double_t position_X_Front = currentParamForResearchFront.GetNonBendingCoor();
 856   Double_t position_Y_Front = currentParamForResearchFront.GetBendingCoor();
 857   Double_t position_X_Back = currentParamForResearchBack.GetNonBendingCoor();
 858   Double_t position_Y_Back = currentParamForResearchBack.GetBendingCoor();
 859   Double_t radialPositionOfTrackFront = TMath::Sqrt(position_X_Front*position_X_Front + position_Y_Front*position_Y_Front);
 860   Double_t radialPositionOfTrackBack  = TMath::Sqrt(position_X_Back*position_X_Back   + position_Y_Back*position_Y_Back);
 861
 862   //---------------------------------------------------------------------------------------
 863
 864   Double_t chi2cut = 2.*fSigmaClusterCut*fSigmaClusterCut;     // depends on the number of variables (here, 2)
 865
 866   // Analyizing the clusters: FRONT ACTIVE ELEMENTS
 867
 868   Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
 869   AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId));
 870
 871   for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
 872
 873     Bool_t isGoodChi2 = kFALSE;
 874
 875     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster);
 876     Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster);     // describes the compatibility between the track and the cluster
 877     if (chi2<chi2cut) isGoodChi2 = kTRUE;
 878
 879     Double_t radialPositionOfClusterFront = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());
 880     if (planeId == fNPlanesMFT-1) {
 881       if (TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront)<fDistanceFromBestClusterAndTrackAtLastPlane ||
 882           fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
 883         fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
 884       }
 885       if (IsCorrectMatch(cluster)) {
 886         fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
 887       }
 888     }
 889
 890     if (fIsCurrentMuonTrackable) {
 891       //      fOutputQAFile->cd();
 892       if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.);     //  chi2/ndf
 893       else                         fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.);     //  chi2/ndf
 894     }
 895
 896     if (isGoodChi2) {
 897       AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 898       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
 899       newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);    // creating new track param and attaching the cluster
 900       AliDebug(1, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
 901                        planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
 902       newTrack->SetPlaneExists(planeId);
 903       AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable));
 904       if (fIsCurrentMuonTrackable) {
 905         Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
 906         AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2()));
 907         Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters();        // Muon Spectrometer clusters + clusters in the Vertex Telescope
 908         Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
 909         //      fOutputQAFile->cd();
 910         if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
 911         else                         fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
 912       }
 913       fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY());
 914     }
 915     else AliDebug(3, Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 916
 917   }
 918
 919   // Analyizing the clusters: BACK ACTIVE ELEMENTS
 920
 921   Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
 922   AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId));
 923
 924   for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
 925
 926     Bool_t isGoodChi2 = kFALSE;
 927
 928     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster);
 929     Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster);     // describes the compatibility between the track and the cluster
 930     if (chi2<chi2cut) isGoodChi2 = kTRUE;
 931
 932     Double_t radialPositionOfClusterBack = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());
 933     if (planeId == fNPlanesMFT-1) {
 934       if (TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack)<fDistanceFromBestClusterAndTrackAtLastPlane ||
 935           fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
 936         fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
 937       }
 938       if (IsCorrectMatch(cluster)) {
 939         fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
 940       }
 941     }
 942
 943     if (fIsCurrentMuonTrackable) {
 944       //      fOutputQAFile->cd();
 945       if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.);     //  chi2/ndf
 946       else                         fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.);     //  chi2/ndf
 947     }
 948
 949     if (isGoodChi2) {
 950       AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 951       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
 952       newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);    // creating new track param and attaching the cluster
 953       AliDebug(1, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
 954                        planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
 955       newTrack->SetPlaneExists(planeId);
 956       AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable));
 957       if (fIsCurrentMuonTrackable) {
 958         Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
 959         AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2()));
 960         Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters();        // Muon Spectrometer clusters + clusters in the Vertex Telescope
 961         Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
 962         //      fOutputQAFile->cd();
 963         if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
 964         else                         fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
 965       }
 966       fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY());
 967     }
 968     else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 969
 970   }
 971
 972   //---------------------------------------------------------------------------------------------
 973
 974   if (planeId == fNPlanesMFT-1) {
 975     if (fIsCurrentMuonTrackable && fDistanceFromGoodClusterAndTrackAtLastPlane>0.) {
 976       //      fOutputQAFile->cd();
 977       fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Fill(TMath::Abs(fDistanceFromBestClusterAndTrackAtLastPlane-
 978                                                                                                        fDistanceFromGoodClusterAndTrackAtLastPlane));
 979       fHistDistanceGoodClusterFromTrackAtLastPlane -> Fill(fDistanceFromGoodClusterAndTrackAtLastPlane);
 980     }
 981   }
 982
 983 }
 984
 985 //==========================================================================================================================================
 986
 987 void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) {
 988
 989   AliDebug(1, Form(">>>> executing AliMuonForwardTrackFinder::AttachGoodClusterInPlane(%d)\n", planeId));
 990
 991   AliMUONTrackParam currentParamFront, currentParamBack;
 992
 993   if (planeId == fNPlanesMFT-1) {      // last plane of the telecope
 994     currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
 995     currentParamBack  = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
 996     AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, 0.);
 997     AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  0.);
 998   }
 999   else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
1000     AliDebug(2, Form("fCurrentTrack = %p\n", fCurrentTrack));
1001     currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
1002     currentParamBack  = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
1003     AliMUONTrackExtrap::AddMCSEffect(&currentParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
1004                                                           fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
1005     AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,  (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
1006                                                           fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
1007   }
1008   // for all planes: linear extrapolation to the Z of the plane
1009   AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
1010   AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack,  -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
1011
1012   Bool_t goodClusterFound = kFALSE;
1013
1014   // Analyizing the clusters: FRONT ACTIVE ELEMENTS
1015
1016   Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
1017
1018   AliDebug(1, Form("nClustersFront = %d\n", nClustersFront));
1019   for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
1020     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->UncheckedAt(iCluster);
1021     AliDebug(2, Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersFront, cluster, fCurrentTrack));
1022     if (IsCorrectMatch(cluster)) {
1023       fCurrentTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);  // creating new track param and attaching the cluster
1024       fCurrentTrack->SetPlaneExists(planeId);
1025       goodClusterFound = kTRUE;
1026       break;
1027     }
1028   }
1029
1030   if (goodClusterFound) return;
1031
1032   // Analyizing the clusters: BACK ACTIVE ELEMENTS
1033
1034   Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
1035
1036   AliDebug(1, Form("nClustersBack = %d\n", nClustersBack));
1037   for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
1038     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->UncheckedAt(iCluster);
1039     AliDebug(2,Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersBack, cluster, fCurrentTrack));
1040     if (IsCorrectMatch(cluster)) {
1041       fCurrentTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);  // creating new track param and attaching the cluster
1042       fCurrentTrack->SetPlaneExists(planeId);
1043       goodClusterFound = kTRUE;
1044       break;
1045     }
1046   }
1047
1048 }
1049
1050 //==========================================================================================================================================
1051
1052 void AliMuonForwardTrackFinder::CheckCurrentMuonTrackable() {
1053
1054   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1055     fIsGoodClusterInPlane[iPlane] = kFALSE;
1056     Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
1057     for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1058       AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
1059       for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
1060         if (cluster->GetMCLabel(iTrack)==fLabelMC) {
1061           fIsGoodClusterInPlane[iPlane] = kTRUE;
1062           break;
1063         }
1064       }
1065     }
1066   }
1067
1068   fIsCurrentMuonTrackable = kTRUE;
1069   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) fIsCurrentMuonTrackable = (fIsCurrentMuonTrackable&&fIsGoodClusterInPlane[iPlane]);
1070
1071 }
1072
1073 //==========================================================================================================================================
1074
1075 void AliMuonForwardTrackFinder::FillPlanesWithTrackHistory() {
1076
1077   // Fill planes with the clusters
1078
1079   Int_t cluster = 0;
1080   AliDebug(2, Form("fFinalBestCandidate->GetNMFTClusters() = %d\n",  fFinalBestCandidate->GetNMFTClusters()));
1081   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1082     if (fFinalBestCandidate->PlaneExists(iPlane)) {
1083       AliMFTCluster *trackCluster = fFinalBestCandidate->GetMFTCluster(cluster++);
1084       fGrMFTPlane[kClusterOfTrack][iPlane] -> SetPoint(fGrMFTPlane[kClusterOfTrack][iPlane]->GetN(), trackCluster->GetX(), trackCluster->GetY());
1085     }
1086     Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
1087     for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1088       AliMFTCluster *myCluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->UncheckedAt(iCluster);
1089       fGrMFTPlane[kAllClusters][iPlane] -> SetPoint(fGrMFTPlane[kAllClusters][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY());
1090       if (IsCorrectMatch(myCluster)) {
1091         fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetPoint(fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY());
1092       }
1093     }
1094   }
1095
1096 }
1097
1098 //======================================================================================================================================
1099
1100 Bool_t AliMuonForwardTrackFinder::IsCorrectMatch(AliMFTCluster *cluster) {
1101
1102   Bool_t result = kFALSE;
1103
1104   // check if the cluster belongs to the correct MC track
1105
1106   for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
1107     if (cluster->GetMCLabel(iTrack)==fLabelMC) {
1108       result = kTRUE;
1109       break;
1110     }
1111   }
1112
1113   AliDebug(2,Form("returning %d\n", result));
1114
1115   return result;
1116
1117 }
1118
1119 //======================================================================================================================================
1120
1121 Double_t AliMuonForwardTrackFinder::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
1122
1123   // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
1124   // return the corresponding Chi2
1125   // assume the track parameters are given at the Z of the cluster
1126
1127   // Set differences between trackParam and cluster in the bending and non bending directions
1128   Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
1129   Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
1130   AliDebug(3,Form("dX = %f, dY = %f\n", dX, dY));
1131
1132   // Calculate errors and covariances
1133   const TMatrixD& kParamCov = trackParam.GetCovariances();
1134   Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
1135   Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
1136   AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2));
1137   Double_t covXY   = kParamCov(0,2);
1138   Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
1139
1140   // Compute chi2
1141   if (det==0.) return 1.e10;
1142   return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
1143
1144 }
1145
1146 //=========================================================================================================================================
1147
1148 void AliMuonForwardTrackFinder::SeparateFrontBackClusters() {
1149
1150   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1151     fMFTClusterArrayFront[iPlane]->Clear();
1152     fMFTClusterArrayBack[iPlane] ->Clear();
1153     for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
1154       AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
1155       if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
1156         new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
1157       }
1158       else {
1159         new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
1160       }
1161     }
1162   }
1163
1164 }
1165
1166 //=========================================================================================================================================
1167
1168 Int_t AliMuonForwardTrackFinder::GetNDF(Int_t nClusters) {
1169
1170   // the same definition as in AliMUONTrack is implemented, since here we just add more clusters to the Muon track
1171
1172   Int_t ndf = 2 * nClusters - 5;
1173   return (ndf > 0) ? ndf : 0;
1174
1175 }
1176
1177 //============================================================================================================================================
1178
1179 void AliMuonForwardTrackFinder::BookHistos() {
1180
1181   const Int_t nMaxNewTracks[]  = {150,     200,   250, 600, 1000};
1182   const Double_t radiusPlane[] = {0.010, 0.010, 0.050, 0.5,  1.5};
1183
1184   fHistRadiusEndOfAbsorber = new TH1D("hRadiusEndOfAbsorber", "Track radial distance at the end of the absorber",  1000, 0, 100.);
1185
1186   fHistNGoodClustersForFinalTracks = new TH1D("hNGoodClustersForFinalTracks", "Number of Good Clusters per Final Track", 20, -0.25, 9.75);
1187
1188   fHistDistanceGoodClusterFromTrackAtLastPlane = new TH1D("hDistanceGoodClusterFromTrackAtLastPlane",
1189                                                           "Distance of MC Good Cluster from Track in last MFT plane", 200, 0., 2.);
1190
1191   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane =
1192     new TH1D("hDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane",
1193              "Good Cluster distance from track - Best Cluster distance from track in last MFT plane", 200, 0., 2.);
1194
1195   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1196
1197     fHistNTracksAfterExtrapolation[iPlane] = new TH1D(Form("hNTracksAfterExtrapolation_pl%02d", iPlane),
1198                                                       Form("Number of Candidates after analysis of MFT plane %02d", iPlane),
1199                                                       nMaxNewTracks[iPlane], -0.5, nMaxNewTracks[iPlane]-0.5);
1200
1201     fHistResearchRadius[iPlane] = new TH1D(Form("hResearchRadius_pl%02d", iPlane),
1202                                            Form("Research Radius for candidate clusters in MFT plane %02d", iPlane),
1203                                            1000, 0., radiusPlane[iPlane]);
1204
1205     fHistChi2Cluster_GoodCluster[iPlane] = new TH1D(Form("hChi2Cluster_GoodCluster_pl%02d", iPlane),
1206                                                     Form("#chi^{2}_{clust} for Good clusters in MFT plane %02d", iPlane),
1207                                                     100, 0., 15.);
1208
1209     fHistChi2Cluster_BadCluster[iPlane] = new TH1D(Form("hChi2Cluster_BadCluster_pl%02d", iPlane),
1210                                                    Form("#chi^{2}_{clust} for Bad clusters in MFT plane %02d", iPlane),
1211                                                    100, 0., 15.);
1212
1213     fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1214                                                                                  Form("#chi^{2}/ndf at plane %d for GOOD candidates of trackable muons",iPlane),
1215                                                                                  100, 0., 15.);
1216
1217     fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1218                                                                                 Form("#chi^{2}/ndf at plane %d for BAD candidates of trackable muons",iPlane),
1219                                                                                 100, 0., 15.);
1220
1221   }
1222
1223   //------------------------------------------
1224
1225   fHistRadiusEndOfAbsorber          -> Sumw2();
1226   fHistNGoodClustersForFinalTracks  -> Sumw2();
1227
1228   fHistDistanceGoodClusterFromTrackAtLastPlane                                  -> Sumw2();
1229   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Sumw2();
1230
1231   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1232
1233     fHistNTracksAfterExtrapolation[iPlane] -> Sumw2();
1234     fHistResearchRadius[iPlane]            -> Sumw2();
1235
1236     fHistChi2Cluster_GoodCluster[iPlane]        -> Sumw2();
1237     fHistChi2Cluster_BadCluster[iPlane]         -> Sumw2();
1238
1239     fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
1240     fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> Sumw2();
1241
1242   }
1243
1244   fNtuFinalCandidates     = new TNtuple("ntuFinalCandidates",     "Final Candidates (ALL)", "run:event:muonTrack:nFinalCandidates:MCTrackRefExists:xVtx:yVtx:zVtx:triggerMatch:nClustersMC:nGoodClusters:pt:theta:eta:chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8");
1245
1246   fNtuFinalBestCandidates = new TNtuple("ntuFinalBestCandidates", "Final Best Candidates",  "run:event:muonTrack:nFinalCandidates:MCTrackRefExists:xVtx:yVtx:zVtx:triggerMatch:nClustersMC:nGoodClusters:pt:theta:eta:chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:nClustersAtPlane0:nClustersAtPlane1:nClustersAtPlane2:nClustersAtPlane3:nClustersAtPlane4:nClustersAtPlane5:nClustersAtPlane6:nClustersAtPlane7:nClustersAtPlane8");
1247
1248 }
1249
1250 //============================================================================================================================================
1251
1252 void AliMuonForwardTrackFinder::SetTitleHistos() {
1253
1254   fHistRadiusEndOfAbsorber         -> SetXTitle("R_{abs}  [cm]");
1255   fHistNGoodClustersForFinalTracks -> SetXTitle("N_{GoodClusters}");
1256
1257   fHistDistanceGoodClusterFromTrackAtLastPlane                                  -> SetXTitle("Distance  [cm]");
1258   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> SetXTitle("Distance  [cm]");
1259
1260
1261   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1262
1263     fHistNTracksAfterExtrapolation[iPlane] -> SetXTitle("N_{tracks}");
1264     fHistResearchRadius[iPlane]            -> SetXTitle("Research Radius  [cm]");
1265
1266     fHistChi2Cluster_GoodCluster[iPlane]         -> SetXTitle("#chi^{2}/ndf");
1267     fHistChi2Cluster_BadCluster[iPlane]          -> SetXTitle("#chi^{2}/ndf");
1268
1269     fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
1270     fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> SetXTitle("#chi^{2}/ndf");
1271
1272   }
1273
1274 }
1275
1276 //===========================================================================================================================================
1277
1278 void AliMuonForwardTrackFinder::BookPlanes() {
1279
1280   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1281     fGrMFTPlane[kAllClusters][iPlane] = new TGraph();
1282     fGrMFTPlane[kAllClusters][iPlane] -> SetName(Form("fGrMFTPlane_%02d_AllClusters",iPlane));
1283     fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerStyle(20);
1284     //    fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.5);
1285     //    fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.3);
1286     fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.2);
1287   }
1288
1289   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1290     fGrMFTPlane[kClustersGoodChi2][iPlane] = new TGraph();
1291     fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersGoodChi2",iPlane));
1292     fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerStyle(20);
1293     //    fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.8);
1294     //    fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.4);
1295     fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.3);
1296     fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerColor(kBlue);
1297   }
1298
1299   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1300     fGrMFTPlane[kClusterOfTrack][iPlane] = new TGraph();
1301     fGrMFTPlane[kClusterOfTrack][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersOfTrack",iPlane));
1302     fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerStyle(25);
1303     //    fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(1.2);
1304     fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(0.9);
1305     fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerColor(kRed);
1306     fGrMFTPlane[kClusterOfTrack][iPlane] -> SetTitle(Form("Plane %d (%3.1f cm)", iPlane, fZPlane[iPlane]));
1307   }
1308
1309   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1310     fGrMFTPlane[kClusterCorrectMC][iPlane] = new TGraph();
1311     fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersCorrectMC",iPlane));
1312     fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerStyle(20);
1313     //    fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.8);
1314     fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.5);
1315     fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerColor(kGreen);
1316   }
1317
1318   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1319     fCircleExt[iPlane] = new TEllipse(0., 0., fRPlaneMax[iPlane], fRPlaneMax[iPlane]);
1320     fCircleInt[iPlane] = new TEllipse(0., 0., fRPlaneMin[iPlane], fRPlaneMin[iPlane]);
1321   }
1322
1323   fTxtDummy = new TLatex(0.10, 0.59, "Best Candidate:");
1324
1325   //---------------------------------------------------
1326
1327   fMrkAllClust = new TMarker(0.10, 0.32, 20);
1328   fMrkAllClust -> SetMarkerSize(0.5);
1329
1330   fMrkClustGoodChi2 = new TMarker(0.10, 0.26, 20);
1331   fMrkClustGoodChi2 -> SetMarkerSize(0.8);
1332   fMrkClustGoodChi2 -> SetMarkerColor(kBlue);
1333
1334   fMrkClustMC = new TMarker(0.10, 0.20, 20);
1335   fMrkClustMC -> SetMarkerSize(0.8);
1336   fMrkClustMC -> SetMarkerColor(kGreen);
1337
1338   fMrkClustOfTrack = new TMarker(0.10, 0.14, 25);
1339   fMrkClustOfTrack -> SetMarkerSize(1.2);
1340   fMrkClustOfTrack -> SetMarkerColor(kRed);
1341
1342   fTxtAllClust = new TLatex(0.15, 0.30, "All Clusters");
1343   fTxtAllClust -> SetTextSize(0.040);
1344
1345   fTxtClustGoodChi2 = new TLatex(0.15, 0.24, "Clusters involved in the research");
1346   fTxtClustGoodChi2 -> SetTextSize(0.040);
1347
1348   fTxtClustMC = new TLatex(0.15, 0.18, "MC good clusters");
1349   fTxtClustMC -> SetTextSize(0.040);
1350
1351   fTxtClustOfTrack = new TLatex(0.15, 0.12, "Clusters of the best candidate");
1352   fTxtClustOfTrack -> SetTextSize(0.040);
1353
1354 }
1355
1356 //===========================================================================================================================================
1357
1358 void AliMuonForwardTrackFinder::ResetPlanes() {
1359
1360   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1361     for (Int_t iGr=0; iGr<4; iGr++) {
1362       Int_t nOldClusters = fGrMFTPlane[iGr][iPlane]->GetN();
1363       for (Int_t iPoint=nOldClusters-1; iPoint>=0; iPoint--) fGrMFTPlane[iGr][iPlane]->RemovePoint(iPoint);
1364     }
1365   }
1366
1367 }
1368
1369 //===========================================================================================================================================
1370
1371 void AliMuonForwardTrackFinder::PrintParticleHistory() {
1372
1373   AliDebug(1, "Entering");
1374
1375   TString history = "";
1376
1377   TParticle *part = 0;
1378   if (fLabelMC>=0) part = fStack->Particle(fLabelMC);
1379
1380   AliDebug(1, Form("fStack->Particle(fLabelMC) = %p", part));
1381
1382   if (part) {
1383     if (part->GetFirstMother() != -1) {
1384       TParticle *partMother = fStack->Particle(part->GetFirstMother());
1385       AliDebug(1, Form("fStack->Particle(part->GetFirstMother() = %p", partMother));
1386       if (partMother) {
1387         Char_t newName[100];
1388         if (partMother->GetFirstMother() != -1) history += "...  #rightarrow ";
1389         PDGNameConverter(partMother->GetName(), newName);
1390         history += Form("%s #rightarrow ", newName);
1391       }
1392     }
1393     Char_t newName[100];
1394     PDGNameConverter(part->GetName(), newName);
1395     history += Form("%s  at  z = %5.1f cm", newName, part->Vz());
1396     //  printf("%s", history.Data());
1397   }
1398   else history += "NO AVAILABLE HISTORY";
1399
1400   fTxtMuonHistory = new TLatex(0.10, 0.86, history.Data());
1401
1402   // Filling particle history in the fFinalBestCandidate
1403
1404   if (part) {
1405     for (Int_t iParent=0; iParent<AliMuonForwardTrack::fgkNParentsMax; iParent++) {
1406       if (part->GetFirstMother() == -1) break;
1407       if (!(fStack->Particle(part->GetFirstMother()))) break;
1408       AliDebug(1, Form("fStack->Particle(part->GetFirstMother() = %p", fStack->Particle(part->GetFirstMother())));
1409       fFinalBestCandidate->SetParentMCLabel(iParent, part->GetFirstMother());
1410       fFinalBestCandidate->SetParentPDGCode(iParent, fStack->Particle(part->GetFirstMother())->GetPdgCode());
1411       part = fStack->Particle(part->GetFirstMother());
1412     }
1413   }
1414
1415 }
1416
1417 //===========================================================================================================================================
1418
1419 Bool_t AliMuonForwardTrackFinder::IsMother(Char_t *nameMother) {
1420
1421   Bool_t result = kFALSE;
1422
1423   TParticle *part = 0;
1424   if (fLabelMC>=0) part = fStack->Particle(fLabelMC);
1425
1426   if (part) {
1427     if (part->GetFirstMother() != -1) {
1428       TParticle *partMother = fStack->Particle(part->GetFirstMother());
1429       if (partMother) {
1430         if (!strcmp(partMother->GetName(), nameMother)) result=kTRUE;
1431       }
1432     }
1433   }
1434
1435   return result;
1436
1437 }
1438
1439 //===========================================================================================================================================
1440
1441 void AliMuonForwardTrackFinder::DrawPlanes() {
1442
1443   fCanvas -> Clear();
1444   if (fNPlanesMFT <= 5)       fCanvas -> Divide(3,2);
1445   else if (fNPlanesMFT <= 11) fCanvas -> Divide(4,3);
1446   else if (fNPlanesMFT <= 19) fCanvas -> Divide(5,4);
1447
1448   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1449
1450     fCanvas->cd(fNPlanesMFT-iPlane+1);
1451
1452     fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetLimits(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1453     fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetRangeUser(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1454     fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetTitle("X  [cm]");
1455     fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetTitle("Y  [cm]");
1456     fGrMFTPlane[kClusterOfTrack][iPlane] -> Draw("ap");
1457
1458     fCircleExt[iPlane] -> Draw("same");
1459     fCircleInt[iPlane] -> Draw("same");
1460
1461     if (fGrMFTPlane[kAllClusters][iPlane]->GetN())       fGrMFTPlane[kAllClusters][iPlane]      -> Draw("psame");
1462     if (fGrMFTPlane[kClustersGoodChi2][iPlane]->GetN())  fGrMFTPlane[kClustersGoodChi2][iPlane] -> Draw("psame");
1463     if (fGrMFTPlane[kClusterOfTrack][iPlane]->GetN())    fGrMFTPlane[kClusterOfTrack][iPlane]   -> Draw("psame");
1464     if (fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN())  fGrMFTPlane[kClusterCorrectMC][iPlane] -> Draw("psame");
1465
1466     fTxtTrackChi2[iPlane] -> Draw("same");
1467
1468   }
1469
1470   fCanvas -> cd(1);
1471   fTxtMuonHistory       -> Draw();
1472   fTxtDummy             -> Draw("same");
1473   if (fMatchingMode==kRealMatching) fTxtTrackGoodClusters -> Draw("same");
1474   fTxtTrackFinalChi2    -> Draw("same");
1475   fTxtTrackMomentum     -> Draw("same");
1476   if (fMatchingMode==kRealMatching) fTxtFinalCandidates   -> Draw("same");
1477
1478   fMrkAllClust      -> Draw("same");
1479   fMrkClustGoodChi2 -> Draw("same");
1480   fMrkClustMC       -> Draw("same");
1481   fMrkClustOfTrack  -> Draw("same");
1482
1483   fTxtAllClust      -> Draw("same");
1484   fTxtClustGoodChi2 -> Draw("same");
1485   fTxtClustMC       -> Draw("same");
1486   fTxtClustOfTrack  -> Draw("same");
1487
1488   //  fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1489   fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1490   if (IsMother("phi")) {
1491     fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1492     fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1493   }
1494   if (IsMother("J/psi")) {
1495     fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1496     fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1497   }
1498
1499 }
1500
1501 //===========================================================================================================================================
1502
1503 void AliMuonForwardTrackFinder::Terminate() {
1504
1505   AliInfo("");
1506   AliInfo("---------------------------------------------------------------------------------------------------------------");
1507   AliInfo(Form("%8d  tracks analyzed",                                                     fCountRealTracksAnalyzed));
1508   AliInfo(Form("%8d  tracks with MC ref",                                                  fCountRealTracksWithRefMC));
1509   AliInfo(Form("%8d  tracks with MC ref & trigger match",                                  fCountRealTracksWithRefMC_andTrigger));
1510   if (fMatchingMode==kRealMatching) {
1511     AliInfo(Form("%8d  tracks analyzed with final candidates",                             fCountRealTracksAnalyzedWithFinalCandidates));
1512   }
1513   else {
1514     AliInfo(Form("%8d  tracks matched with their MC clusters",                             fCountRealTracksAnalyzedWithFinalCandidates));
1515   }
1516 //   printf("%8d  tracks with MC ref & trigger match & pt>%3.1f GeV/c",                 fCountRealTracksWithRefMC_andTrigger_andGoodPt, fLowPtCut);
1517 //   printf("%8d  tracks with MC ref & trigger match & pt>%3.1f GeV/c & correct R_abs", fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta, fLowPtCut);
1518   AliInfo("---------------------------------------------------------------------------------------------------------------");
1519
1520   WriteOutputTree();
1521   WriteHistos();
1522
1523 }
1524
1525 //==========================================================================================================================================
1526
1527 void AliMuonForwardTrackFinder::FillOutputTree() {
1528
1529   if (!fMuonForwardTracks || !fOutputEventTree) return;
1530
1531   AliDebug(1, Form("Filling output tree %p with %p having %d entries whose 1st entry is %p",
1532                    fOutputEventTree, fMuonForwardTracks, fMuonForwardTracks->GetEntries(), fMuonForwardTracks->At(0)));
1533
1534   //  fOutputTreeFile->cd();
1535   fOutputEventTree->Fill();
1536   AliDebug(1, Form("\nFilled Tree: nEvents = %d!!!!\n", Int_t(fOutputEventTree->GetEntries())));
1537
1538 }
1539
1540 //==========================================================================================================================================
1541
1542 void AliMuonForwardTrackFinder::WriteOutputTree() {
1543
1544   if (!fOutputEventTree || !fOutputTreeFile) return;
1545
1546   fOutputTreeFile -> cd();
1547
1548   fOutputEventTree -> Write();
1549   fOutputTreeFile -> Close();
1550
1551 }
1552
1553 //==========================================================================================================================================
1554
1555 void AliMuonForwardTrackFinder::WriteHistos() {
1556
1557   fOutputQAFile = new TFile(Form("MuonGlobalTracking.QA.run%d.root", fRun), "recreate");
1558   fOutputQAFile -> cd();
1559
1560   fHistRadiusEndOfAbsorber         -> Write();
1561   fHistNGoodClustersForFinalTracks -> Write();
1562
1563   fHistDistanceGoodClusterFromTrackAtLastPlane                                  -> Write();
1564   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Write();
1565
1566   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1567
1568     fHistNTracksAfterExtrapolation[iPlane] -> Write();
1569     fHistResearchRadius[iPlane]            -> Write();
1570
1571     fHistChi2Cluster_GoodCluster[iPlane]   -> Write();
1572     fHistChi2Cluster_BadCluster[iPlane]    -> Write();
1573
1574     fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Write();
1575     fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> Write();
1576
1577   }
1578
1579   fNtuFinalCandidates     -> Write();
1580   fNtuFinalBestCandidates -> Write();
1581
1582   fOutputQAFile -> Close();
1583
1584 }
1585
1586 //===========================================================================================================================================
1587
1588 void AliMuonForwardTrackFinder::PDGNameConverter(const Char_t *nameIn, Char_t *nameOut) {
1589
1590   if      (!strcmp(nameIn, "mu+"))     sprintf(nameOut, "#mu^{+}");
1591   else if (!strcmp(nameIn, "mu-"))     sprintf(nameOut, "#mu^{-}");
1592   else if (!strcmp(nameIn, "pi+"))     sprintf(nameOut, "#pi^{+}");
1593   else if (!strcmp(nameIn, "pi-"))     sprintf(nameOut, "#pi^{-}");
1594   else if (!strcmp(nameIn, "K+"))      sprintf(nameOut, "K^{+}");
1595   else if (!strcmp(nameIn, "K-"))      sprintf(nameOut, "K^{-}");
1596   else if (!strcmp(nameIn, "K*+"))     sprintf(nameOut, "K^{*+}");
1597   else if (!strcmp(nameIn, "K*-"))     sprintf(nameOut, "K^{*-}");
1598   else if (!strcmp(nameIn, "K_S0"))    sprintf(nameOut, "K_{S}^{0}");
1599   else if (!strcmp(nameIn, "K_L0"))    sprintf(nameOut, "K_{L}^{0}");
1600   else if (!strcmp(nameIn, "K0"))      sprintf(nameOut, "K^{0}");
1601   else if (!strcmp(nameIn, "K0_bar"))  sprintf(nameOut, "#bar{K}^{0}");
1602   else if (!strcmp(nameIn, "K*0"))     sprintf(nameOut, "K^{*0}");
1603   else if (!strcmp(nameIn, "K*0_bar")) sprintf(nameOut, "#bar{K}^{*0}");
1604   else if (!strcmp(nameIn, "rho0"))    sprintf(nameOut, "#rho^{0}");
1605   else if (!strcmp(nameIn, "rho+"))    sprintf(nameOut, "#rho^{+}");
1606   else if (!strcmp(nameIn, "rho-"))    sprintf(nameOut, "#rho^{-}");
1607   else if (!strcmp(nameIn, "omega"))   sprintf(nameOut, "#omega");
1608   else if (!strcmp(nameIn, "eta'"))    sprintf(nameOut, "#eta'");
1609   else if (!strcmp(nameIn, "phi"))     sprintf(nameOut, "#phi");
1610
1611   else if (!strcmp(nameIn, "D-"))     sprintf(nameOut, "D^{-}");
1612   else if (!strcmp(nameIn, "D+"))     sprintf(nameOut, "D^{+}");
1613   else if (!strcmp(nameIn, "D0"))     sprintf(nameOut, "D^{0}");
1614   else if (!strcmp(nameIn, "D0_bar")) sprintf(nameOut, "#bar{D}^{0}");
1615   else if (!strcmp(nameIn, "D*-"))    sprintf(nameOut, "D^{*-}");
1616   else if (!strcmp(nameIn, "D*+"))    sprintf(nameOut, "D^{*+}");
1617   else if (!strcmp(nameIn, "D_s+"))   sprintf(nameOut, "D_{s}^{+}");
1618   else if (!strcmp(nameIn, "D*_s+"))  sprintf(nameOut, "D_{s}^{*+}");
1619
1620   else if (!strcmp(nameIn, "B-"))       sprintf(nameOut, "B^{-}");
1621   else if (!strcmp(nameIn, "B+"))       sprintf(nameOut, "B^{+}");
1622   else if (!strcmp(nameIn, "B_s0_bar")) sprintf(nameOut, "#bar{B}_{s}^{0}");
1623
1624   else if (!strcmp(nameIn, "antiproton"))  sprintf(nameOut, "#bar{p}");
1625   else if (!strcmp(nameIn, "proton"))      sprintf(nameOut, "p");
1626   else if (!strcmp(nameIn, "neutron"))     sprintf(nameOut, "n");
1627   else if (!strcmp(nameIn, "Sigma+"))      sprintf(nameOut, "#Sigma^{+}");
1628   else if (!strcmp(nameIn, "Delta+"))      sprintf(nameOut, "#Delta{+}");
1629   else if (!strcmp(nameIn, "Delta--"))     sprintf(nameOut, "#Delta{--}");
1630   else if (!strcmp(nameIn, "Lambda0"))     sprintf(nameOut, "#Lambda_0");
1631   else if (!strcmp(nameIn, "Lambda0_bar")) sprintf(nameOut, "#bar{Lambda}_0");
1632
1633   else sprintf(nameOut, "%s", nameIn);
1634
1635 }
1636
1637 //===========================================================================================================================================
1638
1639 void AliMuonForwardTrackFinder::SetDraw(Bool_t drawOption) {
1640
1641   fDrawOption = drawOption;
1642
1643   if (!fCanvas) {
1644     fCanvas = new TCanvas("tracking", "tracking", 1200, 800);
1645     fCanvas -> Divide(3,2);
1646   }
1647
1648 }
1649
1650 //===========================================================================================================================================
1651
1652 Bool_t AliMuonForwardTrackFinder::InitGRP() {
1653
1654   //------------------------------------
1655   // Initialization of the GRP entry
1656   //------------------------------------
1657
1658   AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
1659
1660   if (entry) {
1661
1662     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
1663
1664     if (m) {
1665        AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
1666        m->Print();
1667        fGRPData = new AliGRPObject();
1668        fGRPData->ReadValuesFromMap(m);
1669     }
1670
1671     else {
1672        AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it");
1673        fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
1674        entry->SetOwner(0);
1675     }
1676
1677     //    FIX ME: The unloading of GRP entry is temporarily disabled
1678     //    because ZDC and VZERO are using it in order to initialize
1679     //    their reconstructor objects. In the future one has to think
1680     //    of propagating AliRunInfo to the reconstructors.
1681     //    AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
1682   }
1683
1684   if (!fGRPData) {
1685      AliError("No GRP entry found in OCDB!");
1686      return kFALSE;
1687   }
1688
1689   TString lhcState = fGRPData->GetLHCState();
1690   if (lhcState==AliGRPObject::GetInvalidString()) {
1691     AliError("GRP/GRP/Data entry:  missing value for the LHC state ! Using UNKNOWN");
1692     lhcState = "UNKNOWN";
1693   }
1694
1695   TString beamType = fGRPData->GetBeamType();
1696   if (beamType==AliGRPObject::GetInvalidString()) {
1697     AliError("GRP/GRP/Data entry:  missing value for the beam type ! Using UNKNOWN");
1698     beamType = "UNKNOWN";
1699   }
1700
1701   Float_t beamEnergy = fGRPData->GetBeamEnergy();
1702   if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
1703     AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
1704     beamEnergy = 0;
1705   }
1706
1707   TString runType = fGRPData->GetRunType();
1708   if (runType==AliGRPObject::GetInvalidString()) {
1709     AliError("GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN");
1710     runType = "UNKNOWN";
1711   }
1712
1713   Int_t activeDetectors = fGRPData->GetDetectorMask();
1714   if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
1715     AliError("GRP/GRP/Data entry:  missing value for the detector mask ! Using 1074790399");
1716     activeDetectors = 1074790399;
1717   }
1718   AliDebug(1, Form("activeDetectors = %d", activeDetectors));
1719
1720   fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
1721   fRunInfo->Dump();
1722
1723   // *** Dealing with the magnetic field map
1724
1725   if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
1726     if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
1727       AliInfo("ExpertMode!!! GRP information will be ignored !");
1728       AliInfo("ExpertMode!!! Running with the externally locked B field !");
1729     }
1730     else {
1731       AliInfo("Destroying existing B field instance!");
1732       delete TGeoGlobalMagField::Instance();
1733     }
1734   }
1735   if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
1736     // Construct the field map out of the information retrieved from GRP.
1737     Bool_t ok = kTRUE;
1738     // L3
1739     Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
1740     if (l3Current == AliGRPObject::GetInvalidFloat()) {
1741       AliError("GRP/GRP/Data entry:  missing value for the L3 current !");
1742       ok = kFALSE;
1743     }
1744
1745     Char_t l3Polarity = fGRPData->GetL3Polarity();
1746     if (l3Polarity == AliGRPObject::GetInvalidChar()) {
1747       AliError("GRP/GRP/Data entry:  missing value for the L3 polarity !");
1748       ok = kFALSE;
1749     }
1750
1751     // Dipole
1752     Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
1753     if (diCurrent == AliGRPObject::GetInvalidFloat()) {
1754       AliError("GRP/GRP/Data entry:  missing value for the dipole current !");
1755       ok = kFALSE;
1756     }
1757
1758     Char_t diPolarity = fGRPData->GetDipolePolarity();
1759     if (diPolarity == AliGRPObject::GetInvalidChar()) {
1760       AliError("GRP/GRP/Data entry:  missing value for the dipole polarity !");
1761       ok = kFALSE;
1762     }
1763
1764     // read special bits for the polarity convention and map type
1765     Int_t  polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
1766     Bool_t uniformB = fGRPData->IsUniformBMap();
1767
1768     if (ok) {
1769       AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1),
1770                                              TMath::Abs(diCurrent) * (diPolarity ? -1:1),
1771                                              polConvention,uniformB,beamEnergy, beamType.Data());
1772       if (fld) {
1773         TGeoGlobalMagField::Instance()->SetField( fld );
1774         TGeoGlobalMagField::Instance()->Lock();
1775         AliInfo("Running with the B field constructed out of GRP !");
1776       }
1777       else AliFatal("Failed to create a B field map !");
1778     }
1779     else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
1780   }
1781
1782   return kTRUE;
1783 }
1784
1785 //====================================================================================================================================================
1786
1787 Bool_t AliMuonForwardTrackFinder::SetRunNumber() {
1788
1789   AliCDBManager *man = AliCDBManager::Instance();
1790
1791   if (!fRunLoader) {
1792     AliError("No run loader found!");
1793     return kFALSE;
1794   }
1795   else {
1796     fRunLoader->LoadHeader();
1797     // read run number from gAlice
1798     if (fRunLoader->GetHeader()) {
1799       man->SetRun(fRunLoader->GetHeader()->GetRun());
1800       fRunLoader->UnloadHeader();
1801     }
1802     else {
1803       AliError("No run-loader header found!");
1804       return kFALSE;
1805     }
1806   }
1807
1808   return kTRUE;
1809
1810 }
1811
1812 //====================================================================================================================================================
1813