PWG3/muondep/AliAnalysisTaskMuonTrackingEff.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 /* $Id$ */
  17
  18 //Class to calculate the intrinsic efficiency of the detection elements of the
  19 //MUON tracking chambers in function of the position in the detection element.
  20 //WOrk on ESD only
  21 //Author:  Nicolas LE BRIS - SUBATECH Nantes
  22 // Modified by Matthieu LENHARDT - SUBATECH Nantes
  23 // Modified by Antoine LARDEUX - SUBATECH Nantes
  24
  25 // ROOT includes
  26 #include <TList.h>
  27 #include <TH3F.h>
  28 #include <THnSparse.h>
  29 #include <TObjArray.h>
  30 #include <TGeoGlobalMagField.h>
  31 #include <TVector3.h>
  32
  33 // STEER includes
  34 #include "AliESDEvent.h"
  35 #include "AliESDMuonTrack.h"
  36 #include "AliGeomManager.h"
  37 #include "AliCDBManager.h"
  38 #include "AliESDVZERO.h"
  39
  40 // ANALYSIS includes
  41 #include "AliInputEventHandler.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliAnalysisTaskMuonTrackingEff.h"
  44 #include "AliCentrality.h"
  45 #include "AliVVertex.h"
  46
  47 //MUON includes
  48 #include "AliMUONCDB.h"
  49 #include "AliMUONESDInterface.h"
  50 #include "AliMUONGeometryTransformer.h"
  51 #include "AliMUONTrack.h"
  52 #include "AliMUONTrackParam.h"
  53 #include "AliMUONTrackExtrap.h"
  54 #include "AliMUONVCluster.h"
  55 #include "AliMUONConstants.h"
  56
  57 //include MUON/mapping:
  58 #include "AliMpDEManager.h"
  59 #include "AliMpSegmentation.h"
  60 #include "AliMpVSegmentation.h"
  61 #include "AliMpPad.h"
  62
  63 ClassImp(AliAnalysisTaskMuonTrackingEff)
  64
  65 const Int_t AliAnalysisTaskMuonTrackingEff::fgkNbrOfDetectionElt[10] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26};
  66 const Int_t AliAnalysisTaskMuonTrackingEff::fgkOffset = 100;
  67
  68 //________________________________________________________________________
  69 AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff() :
  70   AliAnalysisTaskSE(),
  71   fOCDBLoaded(kFALSE),
  72   fOCDBpath(""),
  73   fMatchTrig(kFALSE),
  74   fApplyAccCut(kFALSE),
  75   fPDCACut(-1.),
  76   fChi2Cut(-1.),
  77   fPtCut(-1.),
  78   fUseMCLabel(kFALSE),
  79   fCurrentCentrality(0.),
  80   fCurrentTrack(0x0),
  81   fTransformer(0x0),
  82   fDetEltTDHistList(0x0),
  83   fDetEltTTHistList(0x0),
  84   fDetEltSDHistList(0x0),
  85   fChamberTDHistList(0x0),
  86   fChamberTTHistList(0x0),
  87   fChamberSDHistList(0x0),
  88   fExtraHistList(0x0)
  89 {
  90   /// Default constructor
  91 }
  92
  93 //________________________________________________________________________
  94 AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff(TString name) :
  95   AliAnalysisTaskSE(name),
  96   fOCDBLoaded(kFALSE),
  97   fOCDBpath("raw://"),
  98   fMatchTrig(kFALSE),
  99   fApplyAccCut(kFALSE),
 100   fPDCACut(-1.),
 101   fChi2Cut(-1.),
 102   fPtCut(-1.),
 103   fUseMCLabel(kFALSE),
 104   fCurrentCentrality(100.),
 105   fCurrentTrack(0x0),
 106   fTransformer(0x0),
 107   fDetEltTDHistList(0x0),
 108   fDetEltTTHistList(0x0),
 109   fDetEltSDHistList(0x0),
 110   fChamberTDHistList(0x0),
 111   fChamberTTHistList(0x0),
 112   fChamberSDHistList(0x0),
 113   fExtraHistList(0x0)
 114 {
 115   /// Constructor
 116
 117   // Output slots 0 to 5 writes into a TClonesArray:
 118   DefineOutput(1, TList::Class());
 119   DefineOutput(2, TList::Class());
 120   DefineOutput(3, TList::Class());
 121   DefineOutput(4, TList::Class());
 122   DefineOutput(5, TList::Class());
 123   DefineOutput(6, TList::Class());
 124   DefineOutput(7, TList::Class());
 125 }
 126
 127 //________________________________________________________________________
 128 AliAnalysisTaskMuonTrackingEff::~AliAnalysisTaskMuonTrackingEff()
 129 {
 130   /// Destructor
 131   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 132     delete fDetEltTDHistList;
 133     delete fDetEltTTHistList;
 134     delete fDetEltSDHistList;
 135     delete fChamberTDHistList;
 136     delete fChamberTTHistList;
 137     delete fChamberSDHistList;
 138     delete fExtraHistList;
 139   }
 140   delete fTransformer;
 141 }
 142
 143 //________________________________________________________________________
 144 void AliAnalysisTaskMuonTrackingEff::NotifyRun()
 145 {
 146   /// Load the OCDB and the Geometry
 147
 148   // Load it only once
 149   if (fOCDBLoaded) return;
 150
 151   // OCDB
 152   AliCDBManager* man = AliCDBManager::Instance();
 153   if (man->IsDefaultStorageSet()) printf("EfficiencyTask: CDB default storage already set!\n");
 154   else man->SetDefaultStorage(fOCDBpath.Data());
 155   if (man->GetRun() > -1) printf("EfficiencyTask: run number already set!\n");
 156   else man->SetRun(fCurrentRunNumber);
 157
 158   // Geometry
 159   if (!AliGeomManager::GetGeometry()) {
 160     AliGeomManager::LoadGeometry();
 161     if (!AliGeomManager::GetGeometry()) return;
 162     if (!AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return;
 163   }
 164   fTransformer = new AliMUONGeometryTransformer();
 165   fTransformer->LoadGeometryData();
 166
 167   // Magnetic field for track extrapolation
 168   if (!TGeoGlobalMagField::Instance()->GetField()) {
 169     if (!AliMUONCDB::LoadField()) return;
 170   }
 171
 172   // Mapping
 173   if (!AliMpSegmentation::Instance(kFALSE)) {
 174     if (!AliMUONCDB::LoadMapping(kTRUE)) return;
 175   }
 176
 177   // RecoParam for refitting
 178   if (!AliMUONESDInterface::GetTracker()) {
 179     AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
 180     if (!recoParam) return;
 181     AliMUONESDInterface::ResetTracker(recoParam);
 182   }
 183
 184   fOCDBLoaded = kTRUE;
 185 }
 186
 187 //________________________________________________________________________
 188 void AliAnalysisTaskMuonTrackingEff::UserCreateOutputObjects()
 189 {
 190   /// Define efficiency histograms
 191
 192   fDetEltTDHistList  = new TList();
 193   fDetEltTDHistList->SetOwner();
 194   fDetEltTTHistList  = new TList();
 195   fDetEltTTHistList->SetOwner();
 196   fDetEltSDHistList  = new TList();
 197   fDetEltSDHistList->SetOwner();
 198   fChamberTDHistList = new TList();
 199   fChamberTDHistList->SetOwner();
 200   fChamberTTHistList = new TList();
 201   fChamberTTHistList->SetOwner();
 202   fChamberSDHistList = new TList();
 203   fChamberSDHistList->SetOwner();
 204   fExtraHistList = new TList();
 205   fExtraHistList->SetOwner();
 206
 207   THnSparse *hn = 0x0;
 208   TH3F *h3 = 0x0;
 209   TString histName, histTitle;
 210
 211   // centrality bins
 212   Int_t nCentBins = 22;
 213   Double_t centRange[2] = {-5., 105.};
 214
 215   // prepare binning for THnSparse
 216   // 1: Ch or DE Id
 217   // 2: centrality
 218   // 3: pt
 219   // 4: y
 220   // 5: sign
 221   const Int_t nDims = 5;
 222   Int_t nBins[nDims] = {0., nCentBins, 20, 15, 2};
 223   Double_t xMin[nDims] = {0., centRange[0], 0., -4., -2.};
 224   Double_t xMax[nDims] = {0., centRange[1], 20., -2.5, 2.};
 225
 226   // global index of DE in the lists
 227   Int_t iDEGlobal = 0;
 228
 229   for (Int_t iCh = 0; iCh < 10; iCh++)
 230   {
 231     // histograms per chamber
 232     nBins[0] = fgkNbrOfDetectionElt[iCh];
 233     xMin[0] = 0.; xMax[0] = static_cast<Double_t>(fgkNbrOfDetectionElt[iCh]);
 234     histTitle.Form("ChamberNbr %d", iCh+1);
 235     histName.Form("TD_ChamberNbr%d", iCh+1);
 236     hn = new THnSparseT<TArrayF>(histName, histTitle, nDims, nBins, xMin, xMax);
 237     fChamberTDHistList->AddAt(hn, iCh);
 238     histName.Form("TT_ChamberNbr%d",iCh+1);
 239     hn = new THnSparseT<TArrayF>(histName, histTitle, nDims, nBins, xMin, xMax);
 240     fChamberTTHistList->AddAt(hn, iCh);
 241     histName.Form("SD_ChamberNbr%d", iCh+1);
 242     hn = new THnSparseT<TArrayF>(histName, histTitle, nDims, nBins, xMin, xMax);
 243     fChamberSDHistList->AddAt(hn, iCh);
 244
 245     // histograms per DE
 246     for (Int_t iDE = 0; iDE < fgkNbrOfDetectionElt[iCh]; iDE++)
 247     {
 248       Int_t deId = FromLocalId2DetElt(iCh, iDE);
 249       histTitle.Form("detEltNbr %d",deId);
 250       if(iCh < 4)
 251       {// chambers 1 -> 4
 252         histName.Form("TD_detEltNbr%d",deId);
 253         h3 = new TH3F(histName, histTitle, 12, -10.0 , 110.0, 12, -10.0, 110.0, nCentBins, centRange[0], centRange[1]);
 254         fDetEltTDHistList->AddAt(h3, iDEGlobal);
 255         histName.Form("TT_detEltNbr%d",deId);
 256         h3 = new TH3F(histName, histTitle, 12, -10.0 , 110.0, 12, -10.0, 110.0, nCentBins, centRange[0], centRange[1]);
 257         fDetEltTTHistList->AddAt(h3, iDEGlobal);
 258         histName.Form("SD_detEltNbr%d",deId);
 259         h3 = new TH3F(histName, histTitle, 12, -10.0 , 110.0, 12, -10.0, 110.0, nCentBins, centRange[0], centRange[1]);
 260         fDetEltSDHistList->AddAt(h3, iDEGlobal);
 261       }
 262       else
 263       {// chambers 5 -> 10
 264         histName.Form("TD_detEltNbr%d",deId);
 265         h3 = new TH3F(histName, histTitle, 28, -140.0, 140.0, 8, -40.0, 40.0, nCentBins, centRange[0], centRange[1]);
 266         fDetEltTDHistList->AddAt(h3, iDEGlobal);
 267         histName.Form("TT_detEltNbr%d",deId);
 268         h3 = new TH3F(histName, histTitle, 28, -140.0, 140.0, 8, -40.0, 40.0, nCentBins, centRange[0], centRange[1]);
 269         fDetEltTTHistList->AddAt(h3, iDEGlobal);
 270         histName.Form("SD_detEltNbr%d",deId);
 271         h3 = new TH3F(histName, histTitle, 28, -140.0, 140.0, 8, -40.0, 40.0, nCentBins, centRange[0], centRange[1]);
 272         fDetEltSDHistList->AddAt(h3, iDEGlobal);
 273       }
 274       iDEGlobal++;
 275     }
 276   }
 277
 278   // global histograms per chamber
 279   nBins[0] = 10;
 280   xMin[0] = 0.5; xMax[0] = 10.5;
 281   hn = new THnSparseT<TArrayF>("TD_Chambers 11", "Chambers 11", nDims, nBins, xMin, xMax);
 282   fChamberTDHistList->AddAt(hn, 10);
 283   hn = new THnSparseT<TArrayF>("TT_Chambers 11", "Chambers 11", nDims, nBins, xMin, xMax);
 284   fChamberTTHistList->AddAt(hn, 10);
 285   hn = new THnSparseT<TArrayF>("SD_Chambers 11", "Chambers 11", nDims, nBins, xMin, xMax);
 286   fChamberSDHistList->AddAt(hn, 10);
 287
 288   //Extra histograms
 289   TH1F *fHistCent = new TH1F("fHistCent", "centrality distribution", nCentBins, centRange[0], centRange[1]);
 290   fExtraHistList->AddAt(fHistCent,0);
 291   TH1F *fHistPt = new TH1F("fHistPt", "pt distribution", 250, 0., 50.);
 292   fExtraHistList->AddAt(fHistPt,1);
 293   TH1F *fHistY = new TH1F("fHistY", "y distribution", 60, -4., -2.5);
 294   fExtraHistList->AddAt(fHistY,2);
 295   TH1F *fHistTheta = new TH1F("fHistTheta", "theta distribution", 120, 2.8, 3.2);
 296   fExtraHistList->AddAt(fHistTheta,3);
 297   TH1F *fHistP = new TH1F("fHistP", "momentum distribution", 250, 0., 500.);
 298   fExtraHistList->AddAt(fHistP,4);
 299
 300   // post the output data at least once
 301   PostData(1, fDetEltTDHistList);
 302   PostData(2, fDetEltTTHistList);
 303   PostData(3, fDetEltSDHistList);
 304   PostData(4, fChamberTDHistList);
 305   PostData(5, fChamberTTHistList);
 306   PostData(6, fChamberSDHistList);
 307   PostData(7, fExtraHistList);
 308 }
 309
 310 //________________________________________________________________________
 311 void AliAnalysisTaskMuonTrackingEff::UserExec(Option_t *)
 312 {
 313   /// Main event loop
 314
 315   // check the OCDB has been loaded properly
 316   if (!fOCDBLoaded) return;
 317
 318   // get the current event
 319   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 320   if (!esd) return;
 321
 322   // get the centrality
 323   fCurrentCentrality = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
 324   static_cast<TH1F*>(fExtraHistList->At(0))->Fill(fCurrentCentrality);
 325
 326   // loop over tracks
 327   AliMUONTrack track;
 328   Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
 329   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++)
 330   {
 331     AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
 332
 333     if(!esdTrack->ContainTrackerData()) continue;
 334
 335     if(fMatchTrig && !esdTrack->ContainTriggerData()) continue;
 336
 337     Double_t thetaTrackAbsEnd = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 338     Double_t eta = esdTrack->Eta();
 339     if(fApplyAccCut && !(thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 10. && eta >= -4. && eta <= -2.5)) continue;
 340
 341     if (fPDCACut > 0.) {
 342       const AliVVertex* primaryVertex = esd->GetPrimaryVertexSPD();
 343       TVector3 trackDcaAtVz(esdTrack->GetNonBendingCoorAtDCA(), esdTrack->GetBendingCoorAtDCA(), primaryVertex->GetZ());
 344       TVector3 vertex(primaryVertex->GetX(), primaryVertex->GetY(), primaryVertex->GetZ());
 345       TVector3 meanDca(-0.46, -0.92, 0.); // LHC10h1
 346       TVector3 dcaAtVz = trackDcaAtVz - vertex - meanDca;
 347       Double_t correctedDca = dcaAtVz.Mag(); // it should also be equal to dcaAtVz.Pt().
 348       Double_t pMean = 0.5 * (esdTrack->P() + esdTrack->PUncorrected());
 349       Double_t cutVariable = pMean * correctedDca;
 350       Double_t cutValue = (thetaTrackAbsEnd > 3.) ? 63. : 120.;
 351       cutValue = TMath::Sqrt(cutValue*cutValue + 0.4*0.4*esdTrack->P()*esdTrack->P());
 352       if ( cutVariable > fPDCACut*cutValue ) continue;
 353     }
 354
 355     if (fChi2Cut > 0. && esdTrack->GetNormalizedChi2() > fChi2Cut) continue;
 356
 357     if (fPtCut > 0. && esdTrack->Pt() < fPtCut) continue;
 358
 359     if (fUseMCLabel && esdTrack->GetLabel() < 0) continue;
 360
 361     fCurrentTrack = esdTrack;
 362     static_cast<TH1F*>(fExtraHistList->At(1))->Fill(esdTrack->Pt());
 363     static_cast<TH1F*>(fExtraHistList->At(2))->Fill(esdTrack->Y());
 364     static_cast<TH1F*>(fExtraHistList->At(3))->Fill(esdTrack->Theta());
 365     static_cast<TH1F*>(fExtraHistList->At(4))->Fill(esdTrack->P());
 366
 367     AliMUONESDInterface::ESDToMUON(*esdTrack, track);
 368
 369     TrackParamLoop(track.GetTrackParamAtCluster());
 370   }
 371
 372   // post the output data:
 373   PostData(1, fDetEltTDHistList);
 374   PostData(2, fDetEltTTHistList);
 375   PostData(3, fDetEltSDHistList);
 376   PostData(4, fChamberTDHistList);
 377   PostData(5, fChamberTTHistList);
 378   PostData(6, fChamberSDHistList);
 379   PostData(7, fExtraHistList);
 380 }
 381
 382 //________________________________________________________________________
 383 void AliAnalysisTaskMuonTrackingEff::Terminate(Option_t *)
 384 {
 385   /// final plots
 386 }
 387
 388 //________________________________________________________________________
 389 void AliAnalysisTaskMuonTrackingEff::TrackParamLoop(const TObjArray* trackParams)
 390 {
 391   /// Loop on all the track params and fill the histos
 392
 393   Bool_t trackFilter[10];
 394   memset(trackFilter, kFALSE, 10*sizeof(Bool_t));
 395   Bool_t chamberResponse[10];
 396   memset(chamberResponse, kFALSE, 10*sizeof(Bool_t));
 397
 398   // check if the chamber responds
 399   Int_t nTrackParams = (Int_t) trackParams->GetEntriesFast();
 400   for (Int_t iTrackParam = 0; iTrackParam < nTrackParams; ++iTrackParam)
 401   {
 402     Int_t chamberId = static_cast<AliMUONTrackParam*>(trackParams->UncheckedAt(iTrackParam))->GetClusterPtr()->GetChamberId();
 403     trackFilter[chamberId] = kTRUE;
 404     chamberResponse[chamberId] = kTRUE;
 405   }
 406
 407   // To make sure the calculation of the efficiency of a given chamber (DE) is not biased by the tracking algorithm
 408   // we must make sure the track would have been reconstructed whatever this chamber (DE) has responded or not.
 409   // If the track is valid for a given chamber, the following code set trackFilter[chamberId] to kTRUE.
 410   for (Int_t station = 0; station < 4; ++station)
 411   {
 412     Int_t filter;
 413     Int_t ch1 = 2*station;
 414     Int_t ch2 = 2*station + 1;
 415     Int_t ch3 = 2*station + 2;
 416     Int_t ch4 = 2*station + 3;
 417     if (station < 3 )
 418     {
 419       filter           = trackFilter[ch1];
 420       trackFilter[ch1] = trackFilter[ch2];
 421       trackFilter[ch2] = filter;
 422     }
 423     else
 424     {
 425       if (chamberResponse[ch3] && chamberResponse[ch4])
 426       {
 427         filter           = trackFilter[ch1];
 428         trackFilter[ch1] = trackFilter[ch2];
 429         trackFilter[ch2] = filter;
 430       }
 431       else
 432       {
 433         trackFilter[ch1] = kFALSE;
 434         trackFilter[ch2] = kFALSE;
 435       }
 436
 437       if (chamberResponse[ch1] && chamberResponse[ch2])
 438       {
 439         filter           = trackFilter[ch3];
 440         trackFilter[ch3] = trackFilter[ch4];
 441         trackFilter[ch4] = filter;
 442       }
 443       else
 444       {
 445         trackFilter[ch3] = kFALSE;
 446         trackFilter[ch4] = kFALSE;
 447       }
 448     }
 449   }
 450
 451   // loop over track parameters
 452   Int_t oldChamber = -1;
 453   for (Int_t iTrackParam = 0; iTrackParam < nTrackParams; ++iTrackParam)
 454   {
 455     AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(trackParams->UncheckedAt(iTrackParam));
 456     AliMUONVCluster* cluster = trackParam->GetClusterPtr();
 457
 458     Int_t newChamber = cluster->GetChamberId();
 459
 460     Int_t detElt = cluster->GetDetElemId();
 461
 462     ///track position in the global coordinate system
 463     Double_t posXG = trackParam->GetNonBendingCoor();
 464     Double_t posYG = trackParam->GetBendingCoor();
 465     Double_t posZG = trackParam->GetZ();
 466
 467     ///track position in the coordinate system of the DE
 468     Double_t posXL, posYL, posZL;
 469     fTransformer->Global2Local(detElt, posXG, posYG, posZG, posXL, posYL, posZL);
 470
 471     // fill histograms if the track is valid for this chamber
 472     if(trackFilter[newChamber])
 473     {
 474
 475       // fill histograms of the cluster positions on the detection element of the TRACKS DETECTED (TD)
 476       FillTDHistos(newChamber, detElt, posXL, posYL);
 477
 478       // fill histograms of the cluster positions on the detection element of ALL THE TRACKS (TT)
 479       FillTTHistos(newChamber, detElt, posXL, posYL);
 480
 481     } else {
 482
 483      FillSDHistos(newChamber, detElt, posXL, posYL);
 484
 485     }
 486
 487     // look for missing cluster(s) if any
 488     if (newChamber != oldChamber)
 489     {
 490       if (newChamber > oldChamber + 1)
 491       {
 492         Int_t nbrMissChamber = newChamber - (oldChamber + 1);
 493
 494         // find the DE(s) that should have been fired and fill the corresponding histograms
 495         FindAndFillMissedDetElt(trackParam, trackFilter, oldChamber+1, nbrMissChamber);
 496       }
 497
 498       // in case the last chamber has not responded
 499       if ( iTrackParam == nTrackParams-1 && newChamber != 9) FindAndFillMissedDetElt(trackParam, trackFilter, 9, 1);
 500     }
 501
 502     oldChamber = newChamber;
 503   }
 504 }
 505
 506 //________________________________________________________________________
 507 void AliAnalysisTaskMuonTrackingEff::FindAndFillMissedDetElt(const AliMUONTrackParam* trackParam,
 508                                                              const Bool_t* trackFilter,
 509                                                              Int_t firstMissCh, Int_t nbrMissCh)
 510 {
 511   /// Find which detection elements should have been hit but were missed, and fill the TT histos appropriately
 512
 513   // copy track parameters for extrapolation
 514   AliMUONTrackParam extrapTrackParam(*trackParam);
 515
 516   // loop over missing chambers
 517   for (Int_t iCh = 0; iCh < nbrMissCh; ++iCh)
 518   {
 519     Int_t chamber = firstMissCh + iCh;
 520
 521     // skip this chamber if the track is not valid for it
 522     if(!trackFilter[chamber]) continue;
 523
 524     Int_t nbrOfDetElt =  AliMpDEManager::GetNofDEInChamber(chamber, kTRUE);
 525
 526     Double_t pos1[6] = {0, 0, 0, 0, 0, 0};
 527     Double_t pos2[6] = {0, 0, 0, 0, 0, 0};
 528     Double_t posMiss[2] = {0, 0};
 529
 530     // track position at the chamber z
 531     pos1[2] = AliMUONConstants::DefaultChamberZ(chamber);
 532     AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, pos1[2]);
 533     pos1[0] = extrapTrackParam.GetNonBendingCoor();
 534     pos1[1] = extrapTrackParam.GetBendingCoor();
 535
 536     // track position at the chamber z + dz (where dz = distance between the 2 chamber in the station)
 537     pos2[2] = AliMUONConstants::DefaultChamberZ(chamber) + AliMUONConstants::DzCh();
 538     AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, pos2[2]);
 539     pos2[0] = extrapTrackParam.GetNonBendingCoor();
 540     pos2[1] = extrapTrackParam.GetBendingCoor();
 541
 542     // loop over all the detection element of the chamber
 543     for (Int_t iDE = 0; iDE < nbrOfDetElt; iDE++)
 544     {
 545       Int_t deId = (chamber + 1)*fgkOffset + iDE;
 546
 547       // track positions (at chamber z and chamber z + dz) in the local coordinate system of the DE
 548       fTransformer->Global2Local(deId, pos1[0], pos1[1], pos1[2], pos1[3], pos1[4], pos1[5]);
 549       fTransformer->Global2Local(deId, pos2[0], pos2[1], pos2[2], pos2[3], pos2[4], pos2[5]);
 550
 551       // track position at z=0 in the local coordinate system of the DE
 552       CoordinatesOfMissingCluster(pos1[3], pos1[4], pos1[5], pos2[3], pos2[4], pos2[5], posMiss[0], posMiss[1]);
 553
 554       // check if the track cross this DE and fill the corresponding histogram
 555       if (CoordinatesInDetElt(deId, posMiss[0], posMiss[1])) FillTTHistos(chamber, deId, posMiss[0], posMiss[1]);
 556     }
 557   }
 558 }
 559
 560 //________________________________________________________________________
 561 void AliAnalysisTaskMuonTrackingEff::CoordinatesOfMissingCluster(Double_t x1, Double_t y1, Double_t z1,
 562                                                                  Double_t x2, Double_t y2, Double_t z2,
 563                                                                  Double_t& x, Double_t& y) const
 564 {
 565   /// Compute the coordinates of the missing cluster. They are defined by the intersection between
 566   /// the straigth line joining two extrapolated points (1 and 2) and the detection element plane.
 567   /// In the local coordinates, this means Z=0 in the parametric equation of the line.
 568   Double_t t = - z1 / (z2 - z1);
 569   x = t * (x2 - x1) + x1;
 570   y = t * (y2 - y1) + y1;
 571 }
 572
 573 //________________________________________________________________________
 574 Bool_t AliAnalysisTaskMuonTrackingEff::CoordinatesInDetElt(Int_t DeId, Double_t x, Double_t y) const
 575 {
 576   /// Return kTRUE if the coordinates are in the Detection Element.
 577   /// This is done by checking if a pad correspond to the (x, y) position.
 578   const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(DeId, AliMp::kCath0);
 579   const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(DeId, AliMp::kCath1);
 580   if (!seg1 || !seg2) return kFALSE;
 581   AliMpPad pad1 = seg1->PadByPosition(x, y, kFALSE);
 582   AliMpPad pad2 = seg2->PadByPosition(x, y, kFALSE);
 583   return (pad1.IsValid() && pad2.IsValid());
 584 }
 585
 586 //________________________________________________________________________
 587 void AliAnalysisTaskMuonTrackingEff::FillTDHistos(Int_t chamber, Int_t detElt, Double_t posXL, Double_t posYL)
 588 {
 589   /// Fill the histo for detected tracks
 590   static_cast<TH3F*>(fDetEltTDHistList->At(FromDetElt2iDet(chamber, detElt)))->Fill(posXL, posYL, fCurrentCentrality);
 591   Double_t x[5] = {0., fCurrentCentrality, fCurrentTrack->Pt(), fCurrentTrack->Y(), fCurrentTrack->Charge()};
 592   x[0] = static_cast<Double_t>(FromDetElt2LocalId(chamber, detElt));
 593   static_cast<THnSparse*>(fChamberTDHistList->At(chamber))->Fill(x);
 594   x[0] = static_cast<Double_t>(chamber+1);
 595   static_cast<THnSparse*>(fChamberTDHistList->At(10))->Fill(x);
 596 }
 597
 598 //________________________________________________________________________
 599 void AliAnalysisTaskMuonTrackingEff::FillTTHistos(Int_t chamber, Int_t detElt, Double_t posXL, Double_t posYL)
 600 {
 601   /// Fill the histo for all tracks
 602   static_cast<TH3F*>(fDetEltTTHistList->At(FromDetElt2iDet(chamber, detElt)))->Fill(posXL, posYL, fCurrentCentrality);
 603   Double_t x[5] = {0., fCurrentCentrality, fCurrentTrack->Pt(), fCurrentTrack->Y(), fCurrentTrack->Charge()};
 604   x[0] = static_cast<Double_t>(FromDetElt2LocalId(chamber, detElt));
 605   static_cast<THnSparse*>(fChamberTTHistList->At(chamber))->Fill(x);
 606   x[0] = static_cast<Double_t>(chamber+1);
 607   static_cast<THnSparse*>(fChamberTTHistList->At(10))->Fill(x);
 608 }
 609
 610 //________________________________________________________________________
 611 void AliAnalysisTaskMuonTrackingEff::FillSDHistos(Int_t chamber, Int_t detElt, Double_t posXL, Double_t posYL)
 612 {
 613   /// Fill the histo for single detected tracks
 614   static_cast<TH3F*>(fDetEltSDHistList->At(FromDetElt2iDet(chamber, detElt)))->Fill(posXL, posYL, fCurrentCentrality);
 615   Double_t x[5] = {0., fCurrentCentrality, fCurrentTrack->Pt(), fCurrentTrack->Y(), fCurrentTrack->Charge()};
 616   x[0] = static_cast<Double_t>(FromDetElt2LocalId(chamber, detElt));
 617   static_cast<THnSparse*>(fChamberSDHistList->At(chamber))->Fill(x);
 618   x[0] = static_cast<Double_t>(chamber+1);
 619   static_cast<THnSparse*>(fChamberSDHistList->At(10))->Fill(x);
 620 }
 621
 622 //________________________________________________________________________
 623 Int_t AliAnalysisTaskMuonTrackingEff::FromDetElt2iDet(Int_t chamber, Int_t detElt) const
 624 {
 625   /// Connexion between the detection element Id and its position in the list of histograms
 626   Int_t iDet = FromDetElt2LocalId(chamber, detElt);
 627   for (Int_t iCh = chamber-1; iCh >=0; iCh--) iDet += fgkNbrOfDetectionElt[iCh];
 628   return iDet;
 629 }
 630
 631 //________________________________________________________________________
 632 Int_t AliAnalysisTaskMuonTrackingEff::FromDetElt2LocalId(Int_t chamber, Int_t detElt) const
 633 {
 634   /// Connexion between the detection element Id and its number in the chamber
 635   return detElt - fgkOffset*(chamber+1);
 636 }
 637
 638 //________________________________________________________________________
 639 Int_t AliAnalysisTaskMuonTrackingEff::FromLocalId2DetElt(Int_t chamber, Int_t iDet) const
 640 {
 641   /// Connexion between the number of the detection element in the chamber and its Id
 642   return iDet + fgkOffset*(chamber+1);
 643 }
 644