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