MFT/AliMuonForwardTrackFinder.cxx

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