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