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