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