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