MFT/AliMFTTrackerMU.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 //====================================================================================================================================================
  17 //
  18 //                          MFT tracker
  19 //
  20 // Class for the creation of the "global muon tracks" built from the tracks reconstructed in the
  21 // muon spectrometer and the clusters of the Muon Forward Tracker
  22 //
  23 // Contact author: antonio.uras@cern.ch
  24 //
  25 //====================================================================================================================================================
  26
  27 #include "TTree.h"
  28 #include "AliLog.h"
  29 #include "AliGeomManager.h"
  30 #include "AliESDEvent.h"
  31 #include "AliESDMuonTrack.h"
  32 #include "AliESDMuonGlobalTrack.h"
  33 #include "AliMFTTrackerMU.h"
  34 #include "TMath.h"
  35 #include "AliRun.h"
  36 #include "AliRunLoader.h"
  37 #include "AliHeader.h"
  38 #include "AliGenEventHeader.h"
  39 #include "TArrayF.h"
  40 #include "AliMFT.h"
  41 #include "AliMUONTrackExtrap.h"
  42 #include "AliMUONTrack.h"
  43 #include "AliMUONESDInterface.h"
  44 #include "AliMuonForwardTrack.h"
  45 #include "AliMUONConstants.h"
  46 #include "TGrid.h"
  47
  48 ClassImp(AliMFTTrackerMU)
  49
  50 const Double_t AliMFTTrackerMU::fRadLengthSi = AliMFTConstants::fRadLengthSi;
  51
  52 //====================================================================================================================================================
  53
  54 AliMFTTrackerMU::AliMFTTrackerMU() :
  55   AliTracker(),
  56   fESD(0),
  57   fMFT(0),
  58   fSegmentation(0),
  59   fNPlanesMFT(0),
  60   fNPlanesMFTAnalyzed(0),
  61   fSigmaClusterCut(2),
  62   fScaleSigmaClusterCut(1.),
  63   fNMaxMissingMFTClusters(0),
  64   fGlobalTrackingDiverged(kFALSE),
  65   fCandidateTracks(0),
  66   fMUONTrack(0),
  67   fCurrentTrack(0),
  68   fFinalBestCandidate(0),
  69   fXExtrapVertex(0),
  70   fYExtrapVertex(0),
  71   fZExtrapVertex(0),
  72   fXExtrapVertexError(0),
  73   fYExtrapVertexError(0),
  74   fBransonCorrection(kFALSE)
  75 {
  76
  77   //--------------------------------------------------------------------
  78   // This is the AliMFTTrackerMU constructor
  79   //--------------------------------------------------------------------
  80
  81   fMFT = (AliMFT*) gAlice->GetDetector("MFT");
  82   fSegmentation = fMFT->GetSegmentation();
  83   SetNPlanesMFT(fSegmentation->GetNPlanes());
  84   AliMUONTrackExtrap::SetField();                 // set the magnetic field for track extrapolations
  85
  86   AliInfo(Form("fMFT = %p, fSegmentation = %p", fMFT, fSegmentation));
  87
  88   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
  89     fMFTClusterArray[iPlane]      = new TClonesArray("AliMFTCluster");
  90     fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
  91     fMFTClusterArrayBack[iPlane]  = new TClonesArray("AliMFTCluster");
  92     fMFTClusterArray[iPlane]      -> SetOwner(kTRUE);
  93     fMFTClusterArrayFront[iPlane] -> SetOwner(kTRUE);
  94     fMFTClusterArrayBack[iPlane]  -> SetOwner(kTRUE);
  95     fMinResearchRadiusAtPlane[iPlane] = 0.;
  96   }
  97
  98   fCandidateTracks = new TClonesArray("AliMuonForwardTrack",50000);
  99
 100 }
 101
 102 //====================================================================================================================================================
 103
 104 AliMFTTrackerMU::~AliMFTTrackerMU() {
 105
 106   // destructor
 107
 108   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 109     fMFTClusterArray[iPlane] -> Delete();
 110     delete fMFTClusterArray[iPlane];
 111     delete fMFTClusterArrayFront[iPlane];
 112     delete fMFTClusterArrayBack[iPlane];
 113   }
 114
 115   delete fCandidateTracks;
 116
 117 }
 118
 119 //====================================================================================================================================================
 120
 121 Int_t AliMFTTrackerMU::LoadClusters(TTree *cTree) {
 122
 123   //--------------------------------------------------------------------
 124   // This function loads the MFT clusters
 125   //--------------------------------------------------------------------
 126
 127   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 128     AliDebug(1, Form("Setting Address for Branch Plane_%02d", iPlane));
 129     cTree->SetBranchAddress(Form("Plane_%02d",iPlane), &fMFTClusterArray[iPlane]);
 130   }
 131
 132   if (!cTree->GetEvent()) return kFALSE;
 133   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 134     AliInfo(Form("plane %02d: nClusters = %d\n", iPlane, fMFTClusterArray[iPlane]->GetEntries()));
 135   }
 136
 137   AddClustersFromUnderlyingEvent();
 138   AddClustersFromPileUpEvents();
 139
 140   SeparateFrontBackClusters();
 141
 142   return 0;
 143
 144 }
 145
 146 //====================================================================================================================================================
 147
 148 void AliMFTTrackerMU::UnloadClusters() {
 149
 150   //--------------------------------------------------------------------
 151   // This function unloads MFT clusters
 152   //--------------------------------------------------------------------
 153
 154   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 155     fMFTClusterArray[iPlane]      -> Clear("C");
 156     fMFTClusterArrayFront[iPlane] -> Clear("C");
 157     fMFTClusterArrayBack[iPlane]  -> Clear("C");
 158   }
 159
 160 }
 161
 162 //====================================================================================================================================================
 163
 164 Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
 165
 166   //--------------------------------------------------------------------
 167   // This functions reconstructs the Muon Forward Tracks
 168   // The clusters must be already loaded !
 169   //--------------------------------------------------------------------
 170
 171   // Tree of AliMuonForwardTrack objects. Created outside the ESD framework for cross-check purposes
 172
 173   TFile *outputFileMuonGlobalTracks = new TFile("MuonGlobalTracks.root", "update");
 174   TTree *outputTreeMuonGlobalTracks = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
 175   TClonesArray *muonForwardTracks = new TClonesArray("AliMuonForwardTrack");
 176   outputTreeMuonGlobalTracks -> Branch("tracks", &muonForwardTracks);
 177
 178   //--------------------------------------------------------------------
 179
 180   fESD = event;
 181
 182   fXExtrapVertex = 0;
 183   fYExtrapVertex = 0;
 184   fZExtrapVertex = 0;
 185   fXExtrapVertexError = AliMFTConstants::fXVertexTolerance;
 186   fYExtrapVertexError = AliMFTConstants::fYVertexTolerance;
 187
 188   // Imposing a fixed primary vertex, in order to keep memory of its position. Taking the primary vertex would imply the risk
 189   // of loosing memory of its position when passing from ESD to AOD, due to possible refitting
 190   const AliESDVertex* esdVert = fESD->GetVertex();
 191   if (esdVert->GetNContributors() > 0 || !strcmp(esdVert->GetTitle(),"vertexer: smearMC")) {
 192     fXExtrapVertex = esdVert->GetX();
 193     fYExtrapVertex = esdVert->GetY();
 194     fZExtrapVertex = esdVert->GetZ();
 195     fXExtrapVertexError = TMath::Max(AliMFTConstants::fXVertexTolerance, esdVert->GetXRes());
 196     fYExtrapVertexError = TMath::Max(AliMFTConstants::fYVertexTolerance, esdVert->GetYRes());
 197     AliInfo(Form("Found ESD vertex from %d contributors (%f +/- %f,  %f +/- %f,  %f)",
 198                  esdVert->GetNContributors(),fXExtrapVertex,fXExtrapVertexError,fYExtrapVertex,fYExtrapVertexError,fZExtrapVertex));
 199   }
 200   else GetVertexFromMC();
 201
 202   //----------- Read ESD MUON tracks -------------------
 203
 204   Int_t nTracksMUON = event->GetNumberOfMuonTracks();
 205
 206   AliInfo(Form("Number of ESD MUON tracks: %d\n", nTracksMUON));
 207
 208   Int_t iTrack=0;
 209   while (iTrack<nTracksMUON) {
 210
 211     fNPlanesMFTAnalyzed = 0;
 212
 213     AliInfo("****************************************************************************************");
 214     AliInfo(Form("***************************   MUON TRACK %3d/%d   ***************************************", iTrack, nTracksMUON));
 215     AliInfo("****************************************************************************************");
 216
 217     fCandidateTracks -> Delete();
 218
 219     fNPlanesMFTAnalyzed = 0;
 220
 221     const AliESDMuonTrack *esdTrack = event->GetMuonTrack(iTrack);
 222     if (fMUONTrack) delete fMUONTrack;
 223     fMUONTrack = new AliMUONTrack();
 224
 225     AliMUONESDInterface::ESDToMUON(*esdTrack, *fMUONTrack, kFALSE);
 226
 227     if (!fMUONTrack->GetTrackParamAtCluster()->First()) {
 228       AliInfo("Skipping track, no parameters available!!!");
 229       iTrack++;
 230       continue;
 231     }
 232
 233     // the track we are going to build, starting from fMUONTrack and adding the MFT clusters
 234     AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
 235     track -> SetMUONTrack(new AliMUONTrack(*fMUONTrack));
 236     track -> SetMCLabel(fMUONTrack->GetMCLabel());
 237     track -> SetMatchTrigger(fMUONTrack->GetMatchTrigger());
 238
 239     // track parameters linearly extrapolated from the first tracking station to the end of the absorber
 240     AliMUONTrackParam trackParamEndOfAbsorber(*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
 241     AliMUONTrackExtrap::ExtrapToZCov(&trackParamEndOfAbsorber, AliMUONConstants::AbsZEnd());   // absorber extends from -90 to -503 cm
 242     Double_t xEndOfAbsorber = trackParamEndOfAbsorber.GetNonBendingCoor();
 243     Double_t yEndOfAbsorber = trackParamEndOfAbsorber.GetBendingCoor();
 244     Double_t rAbsorber      = TMath::Sqrt(xEndOfAbsorber*xEndOfAbsorber + yEndOfAbsorber*yEndOfAbsorber);
 245     track -> SetRAtAbsorberEnd(rAbsorber);
 246
 247     //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
 248
 249     for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) {   /* *** do not reverse the order of this cycle!!!
 250                                                                  *** this reflects the fact that the extrapolation is performed
 251                                                                  *** starting from the last MFT plane back to the origin */
 252
 253       // --------- updating the array of candidates according to the clusters available in the i-th plane ---------
 254
 255       fNPlanesMFTAnalyzed++;
 256
 257       Int_t nCandidates = fCandidateTracks->GetEntriesFast();
 258       for (Int_t iCandidate=0; iCandidate<nCandidates; iCandidate++) {
 259
 260         if (!(fCandidateTracks->UncheckedAt(iCandidate))) continue;
 261         fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iCandidate);
 262
 263         // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array
 264         // (several new tracks can be created for one old track)
 265         if (FindClusterInPlane(iPlane) == kDiverged) {
 266           fGlobalTrackingDiverged = kTRUE;
 267           break;
 268         }
 269
 270         if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
 271           fCandidateTracks->Remove(fCurrentTrack);     // the old track is removed after the check;
 272         }
 273       }
 274       if (fGlobalTrackingDiverged) {
 275         if (fScaleSigmaClusterCut>0) fScaleSigmaClusterCut -= 0.1;
 276         continue;
 277       }
 278
 279       fCandidateTracks->Compress();
 280
 281     }
 282
 283     // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
 284
 285     fGlobalTrackingDiverged = kFALSE;
 286     fScaleSigmaClusterCut = 1.0;
 287
 288     AliDebug(1, "Finished cycle over planes");
 289
 290     iTrack++;
 291
 292     // If we have several final tracks, we must find the best candidate:
 293
 294     Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
 295     AliInfo(Form("nFinalTracks = %d", nFinalTracks));
 296
 297     Int_t nGoodClustersBestCandidate =  0;
 298     Int_t idBestCandidate            =  0;
 299     Double_t bestChi2                = -1.;  // variable defining the best candidate
 300
 301     for (Int_t iFinalCandidate=0; iFinalCandidate<nFinalTracks; iFinalCandidate++) {
 302
 303       if (!(fCandidateTracks->UncheckedAt(iFinalCandidate))) continue;
 304       AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iFinalCandidate);
 305       Int_t nMFTClusters  = finalTrack->GetNMFTClusters();
 306
 307       Double_t chi2 = 0;
 308       for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
 309         AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
 310         chi2 += localCluster->GetLocalChi2();
 311       }
 312       chi2 /= nMFTClusters;
 313
 314       // now comparing the tracks in order to find the best one
 315
 316       if (chi2<bestChi2 || bestChi2<0) {
 317         bestChi2 = chi2;
 318         idBestCandidate = iFinalCandidate;
 319       }
 320
 321     }
 322
 323     if (nFinalTracks) {
 324
 325       AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
 326       newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
 327
 328       new ((*muonForwardTracks)[muonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
 329
 330       //----------------------- Save the information to the AliESDMuonGlobalTrack object
 331
 332       newTrack -> EvalKinem(AliMFTConstants::fZEvalKinem);
 333
 334       AliDebug(2,"Creating a new Muon Global Track");
 335       AliESDMuonGlobalTrack *myESDTrack = event->NewMuonGlobalTrack();
 336       myESDTrack -> SetPxPyPz(newTrack->Px(), newTrack->Py(), newTrack->Pz());
 337
 338       myESDTrack -> SetLabel(newTrack->GetMCLabel());
 339       myESDTrack -> SetChi2OverNdf(newTrack->GetChi2OverNdf());
 340       myESDTrack -> SetCharge(newTrack->GetCharge());
 341       myESDTrack -> SetMatchTrigger(newTrack->GetMatchTrigger());
 342       myESDTrack -> SetNMFTClusters(newTrack->GetNMFTClusters());
 343       myESDTrack -> SetNWrongMFTClustersMC(newTrack->GetNWrongClustersMC());
 344       myESDTrack -> SetFirstTrackingPoint(newTrack->GetMFTCluster(0)->GetX(), newTrack->GetMFTCluster(0)->GetY(), newTrack->GetMFTCluster(0)->GetZ());
 345       myESDTrack -> SetXYAtVertex(newTrack->GetOffsetX(0., AliMFTConstants::fZEvalKinem), newTrack->GetOffsetY(0., AliMFTConstants::fZEvalKinem));
 346       myESDTrack -> SetRAtAbsorberEnd(newTrack->GetRAtAbsorberEnd());
 347       myESDTrack -> SetCovariances(newTrack->GetTrackParamAtMFTCluster(0)->GetCovariances());
 348       myESDTrack -> SetChi2MatchTrigger(esdTrack->GetChi2MatchTrigger());
 349       myESDTrack -> SetMuonClusterMap(esdTrack->GetMuonClusterMap());
 350       myESDTrack -> SetHitsPatternInTrigCh(esdTrack->GetHitsPatternInTrigCh());
 351       myESDTrack -> SetHitsPatternInTrigChTrk(esdTrack->GetHitsPatternInTrigChTrk());
 352       myESDTrack -> Connected(esdTrack->IsConnected());
 353
 354       ULong_t mftClusterPattern = 0;
 355       for (Int_t iCluster=0; iCluster<newTrack->GetNMFTClusters(); iCluster++) {
 356         AliMFTCluster *localCluster = newTrack->GetMFTCluster(iCluster);
 357         mftClusterPattern |= 1 << localCluster->GetPlane();
 358         mftClusterPattern |= IsCorrectMatch(localCluster, newTrack->GetMCLabel()) << (fNMaxPlanes + localCluster->GetPlane());
 359       }
 360       myESDTrack -> SetMFTClusterPattern(mftClusterPattern);
 361
 362       //---------------------------------------------------------------------------------
 363
 364     }
 365
 366     fFinalBestCandidate = NULL;
 367
 368   }
 369
 370   outputTreeMuonGlobalTracks -> Fill();
 371
 372   Int_t myEventID = 0;
 373   while (outputFileMuonGlobalTracks->cd(Form("Event%d",myEventID))) myEventID++;
 374   outputFileMuonGlobalTracks -> mkdir(Form("Event%d",myEventID));
 375   outputFileMuonGlobalTracks -> cd(Form("Event%d",myEventID));
 376   outputTreeMuonGlobalTracks -> Write();
 377   outputFileMuonGlobalTracks -> Close();
 378
 379   muonForwardTracks -> Delete();
 380   delete muonForwardTracks;
 381
 382   return 0;
 383
 384 }
 385
 386 //=========================================================================================================================================
 387
 388 void AliMFTTrackerMU::SeparateFrontBackClusters() {
 389
 390   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 391     fMFTClusterArrayFront[iPlane]->Delete();
 392     fMFTClusterArrayBack[iPlane] ->Delete();
 393     for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
 394       AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
 395       if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
 396         new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
 397       }
 398       else {
 399         new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
 400       }
 401     }
 402   }
 403
 404 }
 405
 406 //==========================================================================================================================================
 407
 408 Int_t AliMFTTrackerMU::FindClusterInPlane(Int_t planeId) {
 409
 410   // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
 411
 412   // propagate track to plane #planeId (both to front and back active sensors)
 413   // look for compatible clusters
 414   // update TrackParam at found cluster (if any) using Kalman Filter
 415
 416   AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
 417
 418   if (planeId == fNPlanesMFT-1) {      // last plane of the telecope
 419     currentParamFront = (*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
 420     currentParamBack  = (*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
 421     currentParamForResearchFront = currentParamFront;
 422     currentParamForResearchBack  = currentParamBack;
 423     if (fBransonCorrection) {
 424       AliMUONTrackExtrap::ExtrapToVertex(&currentParamFront, fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError);
 425       AliMUONTrackExtrap::ExtrapToVertex(&currentParamBack,  fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError);
 426     }
 427     else {
 428       AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, fZExtrapVertex);
 429       AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  fZExtrapVertex);
 430     }
 431     AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError);
 432     AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack,  fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError);
 433   }
 434   else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
 435     currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
 436     currentParamBack  = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
 437     currentParamForResearchFront = currentParamFront;
 438     currentParamForResearchBack  = currentParamBack;
 439     AliMUONTrackExtrap::AddMCSEffect(&currentParamFront,           (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
 440                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
 441     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
 442                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
 443     AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,            (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
 444                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
 445     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
 446                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
 447   }
 448   // for all planes: extrapolation to the Z of the plane
 449   AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront,            -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
 450   AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
 451   AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack,             -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
 452   AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchBack,  -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
 453
 454   //---------------------------------------------------------------------------------------
 455
 456   TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
 457   TMatrixD covBack(5,5);  covBack  = currentParamForResearchBack.GetCovariances();
 458
 459   Double_t squaredError_X_Front = covFront(0,0);
 460   Double_t squaredError_Y_Front = covFront(2,2);
 461   Double_t squaredError_X_Back  = covBack(0,0);
 462   Double_t squaredError_Y_Back  = covBack(2,2);
 463
 464   Double_t corrFact = 1.0;
 465
 466   Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
 467   Double_t researchRadiusBack  = TMath::Sqrt(squaredError_X_Back  + squaredError_Y_Back);
 468   if (0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtPlane[planeId]) {
 469     corrFact = fMinResearchRadiusAtPlane[planeId]/(0.5*(researchRadiusFront+researchRadiusBack));
 470   }
 471
 472   //---------------------------------------------------------------------------------------
 473
 474   Double_t chi2cut = 2.*fScaleSigmaClusterCut*fScaleSigmaClusterCut*fSigmaClusterCut*fSigmaClusterCut;     // depends on the number of variables (here, 2)
 475
 476   // Analyizing the clusters: FRONT ACTIVE ELEMENTS
 477
 478   Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
 479   AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId));
 480
 481   for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
 482
 483     Bool_t isGoodChi2 = kFALSE;
 484
 485     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster);
 486     Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster);     // describes the compatibility between the track and the cluster
 487     if (chi2<chi2cut) isGoodChi2 = kTRUE;
 488
 489     if (isGoodChi2) {
 490       AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 491       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
 492       if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
 493       newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);    // creating new track param and attaching the cluster
 494       AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
 495                        planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
 496       newTrack->SetPlaneExists(planeId);
 497     }
 498     else AliDebug(3, Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 499
 500   }
 501
 502   // Analyizing the clusters: BACK ACTIVE ELEMENTS
 503
 504   Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
 505   AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId));
 506
 507   for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
 508
 509     Bool_t isGoodChi2 = kFALSE;
 510
 511     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster);
 512     Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster);     // describes the compatibility between the track and the cluster
 513     if (chi2<chi2cut) isGoodChi2 = kTRUE;
 514
 515     if (isGoodChi2) {
 516       AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 517       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
 518       if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
 519       newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);    // creating new track param and attaching the cluster
 520       AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
 521                        planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
 522       newTrack->SetPlaneExists(planeId);
 523     }
 524     else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 525
 526   }
 527
 528   //---------------------------------------------------------------------------------------------
 529
 530   return kConverged;
 531
 532 }
 533
 534 //==========================================================================================================================================
 535
 536 Double_t AliMFTTrackerMU::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
 537
 538   // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
 539   // return the corresponding Chi2
 540   // assume the track parameters are given at the Z of the cluster
 541
 542   // Set differences between trackParam and cluster in the bending and non bending directions
 543   Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
 544   Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
 545   AliDebug(3,Form("dX = %f, dY = %f\n", dX, dY));
 546
 547   // Calculate errors and covariances
 548   const TMatrixD& kParamCov = trackParam.GetCovariances();
 549   Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
 550   Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
 551   AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2));
 552   Double_t covXY   = kParamCov(0,2);
 553   Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
 554
 555   // Compute chi2
 556   if (det==0.) return 1.e10;
 557   return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
 558
 559 }
 560
 561 //=========================================================================================================================================
 562
 563 Bool_t AliMFTTrackerMU::IsCorrectMatch(AliMFTCluster *cluster, Int_t labelMC) {
 564
 565   Bool_t result = kFALSE;
 566
 567   // check if the cluster belongs to the correct MC track
 568
 569   for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
 570     if (cluster->GetMCLabel(iTrack)==labelMC) {
 571       result = kTRUE;
 572       break;
 573     }
 574   }
 575
 576   return result;
 577
 578 }
 579
 580 //======================================================================================================================================
 581
 582 void AliMFTTrackerMU::GetVertexFromMC() {
 583
 584   AliRunLoader *runLoader = AliRunLoader::Open("galice.root");
 585   if (!runLoader) {
 586     AliError("no run loader found in file galice.root");
 587     return;
 588   }
 589
 590   runLoader->CdGAFile();
 591   runLoader->LoadgAlice();
 592   runLoader->LoadHeader();
 593   runLoader->GetEvent(gAlice->GetEvNumber());
 594
 595   TArrayF vtx(3);
 596   runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vtx);
 597   AliInfo(Form("Primary vertex from MC found in (%f, %f, %f)\n",vtx[0], vtx[1], vtx[2]));
 598
 599   fXExtrapVertex = gRandom->Gaus(vtx[0], AliMFTConstants::fPrimaryVertexResX);
 600   fYExtrapVertex = gRandom->Gaus(vtx[1], AliMFTConstants::fPrimaryVertexResY);
 601   fZExtrapVertex = gRandom->Gaus(vtx[2], AliMFTConstants::fPrimaryVertexResZ);
 602   fXExtrapVertexError = AliMFTConstants::fXVertexTolerance;
 603   fYExtrapVertexError = AliMFTConstants::fYVertexTolerance;
 604   AliInfo(Form("Set ESD vertex from MC (%f +/- %f,  %f +/- %f,  %f)",
 605                fXExtrapVertex,fXExtrapVertexError,fYExtrapVertex,fYExtrapVertexError,fZExtrapVertex));
 606
 607 }
 608
 609 //======================================================================================================================================
 610
 611 void AliMFTTrackerMU::AddClustersFromUnderlyingEvent() {
 612
 613   AliInfo("Adding clusters from underlying event");
 614
 615   if (!fMFT) return;
 616
 617   TGrid::Connect("alien://");
 618
 619   TFile* fileWithClustersToAdd = TFile::Open(fMFT->GetFileNameForUnderlyingEvent());
 620   if (!fileWithClustersToAdd) return;
 621   if (!(fileWithClustersToAdd->IsOpen())) return;
 622   if (!(fileWithClustersToAdd->cd(Form("Event%d",fMFT->GetUnderlyingEventID())))) return;
 623
 624   TClonesArray *recPointsPerPlaneToAdd[AliMFTConstants::fNMaxPlanes] = {0};
 625   TTree *treeIn = 0;
 626
 627   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) recPointsPerPlaneToAdd[iPlane] = new TClonesArray("AliMFTCluster");
 628
 629   treeIn = (TTree*) gDirectory->Get("TreeR");
 630
 631   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 632     if (!(treeIn->GetBranch(Form("Plane_%02d",iPlane)))) {
 633       for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
 634       return;
 635     }
 636     else treeIn->SetBranchAddress(Form("Plane_%02d",iPlane), &(recPointsPerPlaneToAdd[iPlane]));
 637   }
 638
 639   treeIn -> GetEntry(0);
 640
 641   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 642     printf("plane %d -> before = %d ",iPlane,fMFTClusterArray[iPlane]->GetEntries());
 643     Int_t nClusters = recPointsPerPlaneToAdd[iPlane]->GetEntries();
 644     for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
 645       AliMFTCluster *newCluster = (AliMFTCluster*) recPointsPerPlaneToAdd[iPlane]->At(iCluster);
 646       for (Int_t iTrack=0; iTrack<newCluster->GetNMCTracks(); iTrack++) newCluster->SetMCLabel(iTrack, newCluster->GetMCLabel(iTrack)+AliMFTConstants::fLabelOffsetMC);
 647       new ((*fMFTClusterArray[iPlane])[fMFTClusterArray[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
 648     }
 649     printf("after = %d\n",fMFTClusterArray[iPlane]->GetEntries());
 650   }
 651
 652   for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
 653
 654 }
 655
 656 //======================================================================================================================================
 657
 658 void AliMFTTrackerMU::AddClustersFromPileUpEvents() {
 659
 660   AliInfo("Adding clusters from pile-up event(s)");
 661
 662   if (!fMFT) return;
 663
 664   TGrid::Connect("alien://");
 665
 666   TFile* fileWithClustersToAdd = TFile::Open(fMFT->GetFileNameForPileUpEvents());
 667   if (!fileWithClustersToAdd) return;
 668   if (!(fileWithClustersToAdd->IsOpen())) return;
 669
 670   TClonesArray *recPointsPerPlaneToAdd[AliMFTConstants::fNMaxPlanes] = {0};
 671   TTree *treeIn = 0;
 672
 673   for (Int_t iPileUp=0; iPileUp<AliMFTConstants::fNMaxPileUpEvents; iPileUp++) {
 674
 675     if (!(fileWithClustersToAdd->cd(Form("Event%d",fMFT->GetPileUpEventID(iPileUp))))) continue;
 676
 677     for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) recPointsPerPlaneToAdd[iPlane] = new TClonesArray("AliMFTCluster");
 678
 679     treeIn = (TTree*) gDirectory->Get("TreeR");
 680
 681     for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 682       if (!(treeIn->GetBranch(Form("Plane_%02d",iPlane)))) {
 683         for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
 684         return;
 685       }
 686       else treeIn->SetBranchAddress(Form("Plane_%02d",iPlane), &(recPointsPerPlaneToAdd[iPlane]));
 687     }
 688
 689     treeIn -> GetEntry(0);
 690
 691     for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
 692       printf("plane %d -> before = %d ",iPlane,fMFTClusterArray[iPlane]->GetEntries());
 693       Int_t nClusters = recPointsPerPlaneToAdd[iPlane]->GetEntries();
 694       for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
 695         AliMFTCluster *newCluster = (AliMFTCluster*) recPointsPerPlaneToAdd[iPlane]->At(iCluster);
 696         for (Int_t iTrack=0; iTrack<newCluster->GetNMCTracks(); iTrack++) newCluster->SetMCLabel(iTrack, newCluster->GetMCLabel(iTrack)+AliMFTConstants::fLabelOffsetMC);
 697         new ((*fMFTClusterArray[iPlane])[fMFTClusterArray[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
 698       }
 699       printf("after = %d\n",fMFTClusterArray[iPlane]->GetEntries());
 700     }
 701
 702     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
 703
 704   }
 705
 706 }
 707
 708 //======================================================================================================================================
 709