PWGPP/MUON/dep/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 // Modified by Philippe PILLOT - SUBATECH Nantes
  25
  26 // ROOT includes
  27 #include <TROOT.h>
  28 #include <TList.h>
  29 #include <THnSparse.h>
  30 #include <TObjArray.h>
  31 #include <TObjString.h>
  32 #include <TGeoGlobalMagField.h>
  33 #include <TVectorD.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TH2D.h>
  37 #include <TFile.h>
  38
  39 // STEER includes
  40 #include "AliESDEvent.h"
  41 #include "AliESDMuonTrack.h"
  42 #include "AliGeomManager.h"
  43 #include "AliCDBManager.h"
  44 #include "AliInputEventHandler.h"
  45 #include "AliCounterCollection.h"
  46
  47 // ANALYSIS includes
  48 #include "AliAnalysisDataSlot.h"
  49 #include "AliAnalysisDataContainer.h"
  50 #include "AliAnalysisManager.h"
  51 #include "AliCentrality.h"
  52 #include "AliAnalysisTaskMuonTrackingEff.h"
  53
  54 //MUON includes
  55 #include "AliMUONCDB.h"
  56 #include "AliMUONESDInterface.h"
  57 #include "AliMUONVTrackReconstructor.h"
  58 #include "AliMUONRecoParam.h"
  59 #include "AliMUONGeometryTransformer.h"
  60 #include "AliMUONTrack.h"
  61 #include "AliMUONTrackParam.h"
  62 #include "AliMUONTrackExtrap.h"
  63 #include "AliMUONVCluster.h"
  64 #include "AliMUONConstants.h"
  65 #include "AliMUON2DMap.h"
  66 #include "AliMUONVCalibParam.h"
  67 #include "AliMUONCalibParamNI.h"
  68 #include "AliMUONTrackerData.h"
  69
  70 //include MUON/mapping:
  71 #include "AliMpDEManager.h"
  72 #include "AliMpSegmentation.h"
  73 #include "AliMpVSegmentation.h"
  74 #include "AliMpPad.h"
  75 #include "AliMpArea.h"
  76 #include "AliMpDEIterator.h"
  77 #include "AliMpConstants.h"
  78 #include "AliMpDDLStore.h"
  79
  80 ClassImp(AliAnalysisTaskMuonTrackingEff)
  81
  82 const Int_t AliAnalysisTaskMuonTrackingEff::fgkNofDE[11] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26, 156};
  83 const Int_t AliAnalysisTaskMuonTrackingEff::fgkNofBusPath = 888;
  84 const Int_t AliAnalysisTaskMuonTrackingEff::fgkNofManu = 16828;
  85
  86 //________________________________________________________________________
  87 AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff() :
  88   AliAnalysisTaskSE(),
  89   fOCDBLoaded(kFALSE),
  90   fOCDBpath(""),
  91   fAlignOCDBpath(""),
  92   fRecoParamOCDBpath(""),
  93   fCentMin(-FLT_MAX),
  94   fCentMax(FLT_MAX),
  95   fMuonTrackCuts(0x0),
  96   fPtCut(-1.),
  97   fUseMCLabel(kFALSE),
  98   fEnableDisplay(kFALSE),
  99   fTransformer(0x0),
 100   fDEPlanes(0x0),
 101   fClusters(0x0),
 102   fChamberTDHistList(0x0),
 103   fChamberTTHistList(0x0),
 104   fChamberSDHistList(0x0),
 105   fExtraHistList(0x0)
 106 {
 107   /// Default constructor
 108 }
 109
 110 //________________________________________________________________________
 111 AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff(TString name) :
 112   AliAnalysisTaskSE(name),
 113   fOCDBLoaded(kFALSE),
 114   fOCDBpath("raw://"),
 115   fAlignOCDBpath(""),
 116   fRecoParamOCDBpath(""),
 117   fCentMin(-FLT_MAX),
 118   fCentMax(FLT_MAX),
 119   fMuonTrackCuts(0x0),
 120   fPtCut(-1.),
 121   fUseMCLabel(kFALSE),
 122   fEnableDisplay(kFALSE),
 123   fTransformer(0x0),
 124   fDEPlanes(0x0),
 125   fClusters(0x0),
 126   fChamberTDHistList(0x0),
 127   fChamberTTHistList(0x0),
 128   fChamberSDHistList(0x0),
 129   fExtraHistList(0x0)
 130 {
 131   /// Constructor
 132
 133   // Output slots
 134   DefineOutput(1, AliCounterCollection::Class());
 135   DefineOutput(2, TList::Class());
 136   DefineOutput(3, TList::Class());
 137   DefineOutput(4, TList::Class());
 138   DefineOutput(5, TList::Class());
 139 }
 140
 141 //________________________________________________________________________
 142 AliAnalysisTaskMuonTrackingEff::~AliAnalysisTaskMuonTrackingEff()
 143 {
 144   /// Destructor
 145   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 146     delete fMuonTrackCuts;
 147     delete fClusters;
 148     delete fChamberTDHistList;
 149     delete fChamberTTHistList;
 150     delete fChamberSDHistList;
 151     delete fExtraHistList;
 152   }
 153   delete fTransformer;
 154   delete fDEPlanes;
 155 }
 156
 157 //________________________________________________________________________
 158 void AliAnalysisTaskMuonTrackingEff::NotifyRun()
 159 {
 160   /// Load the OCDB and the Geometry
 161
 162   // Load it only once
 163   if (fOCDBLoaded) return;
 164
 165   // OCDB
 166   AliCDBManager* man = AliCDBManager::Instance();
 167   if (man->IsDefaultStorageSet()) printf("EfficiencyTask: CDB default storage already set!\n");
 168   else {
 169     man->SetDefaultStorage(fOCDBpath.Data());
 170     if (!fAlignOCDBpath.IsNull()) man->SetSpecificStorage("MUON/Align/Data",fAlignOCDBpath.Data());
 171     if (!fRecoParamOCDBpath.IsNull()) man->SetSpecificStorage("MUON/Calib/RecoParam",fRecoParamOCDBpath.Data());
 172   }
 173   if (man->GetRun() > -1) printf("EfficiencyTask: run number already set!\n");
 174   else man->SetRun(fCurrentRunNumber);
 175
 176   // Geometry
 177   if (!AliGeomManager::GetGeometry()) {
 178     AliGeomManager::LoadGeometry();
 179     if (!AliGeomManager::GetGeometry()) return;
 180     if (!AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return;
 181   }
 182   fTransformer = new AliMUONGeometryTransformer();
 183   fTransformer->LoadGeometryData();
 184
 185   // Mapping
 186   if (!AliMpSegmentation::Instance(kFALSE) || !AliMpDDLStore::Instance(kFALSE)) {
 187     if (!AliMUONCDB::LoadMapping()) return;
 188   }
 189
 190   // vectors (x0, y0, z0, a, b, c) defining the plane of each DE in the global frame
 191   Double_t pl0[3] = {0., 0., 0.};
 192   Double_t pl1[3] = {0., 0., 1.};
 193   Double_t pg0[3], pg1[3];
 194   fDEPlanes = new TObjArray(1026);
 195   fDEPlanes->SetOwner(kTRUE);
 196   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
 197     AliMpDEIterator it;
 198     it.First(i);
 199     while (!it.IsDone()) {
 200       Int_t deId = it.CurrentDEId();
 201       fTransformer->Local2Global(deId, pl0[0], pl0[1], pl0[2], pg0[0], pg0[1], pg0[2]);
 202       fTransformer->Local2Global(deId, pl1[0], pl1[1], pl1[2], pg1[0], pg1[1], pg1[2]);
 203       TVectorD *plane = new TVectorD(6);
 204       (*plane)[0] = pg0[0];
 205       (*plane)[1] = pg0[1];
 206       (*plane)[2] = pg0[2];
 207       (*plane)[3] = pg1[0] - pg0[0];
 208       (*plane)[4] = pg1[1] - pg0[1];
 209       (*plane)[5] = pg1[2] - pg0[2];
 210       fDEPlanes->AddAt(plane, deId);
 211       it.Next();
 212     }
 213   }
 214
 215   // Prepare the tracker (will load RecoParam and magnetic field if not already done)
 216   if (!AliMUONESDInterface::GetTracker()) AliMUONESDInterface::ResetTracker();
 217
 218   // get the trackCuts for this run
 219   if (!fMuonTrackCuts) AliFatal("You must specify the requested selections (AliMuonTrackCut obj is missing)");
 220   fMuonTrackCuts->SetRun(fInputHandler);
 221   fMuonTrackCuts->Print();
 222
 223   fOCDBLoaded = kTRUE;
 224 }
 225
 226 //________________________________________________________________________
 227 void AliAnalysisTaskMuonTrackingEff::UserCreateOutputObjects()
 228 {
 229   /// Define output objects
 230
 231   // count detected, accepted and expected clusters
 232   fClusters = new AliCounterCollection(GetOutputSlot(1)->GetContainer()->GetName());
 233   fClusters->AddRubric("Cluster", "Detected/Accepted/Expected");
 234   fClusters->AddRubric("Chamber", AliMpConstants::NofTrackingChambers());
 235   fClusters->AddRubric("DE", fgkNofDE[10]);
 236   fClusters->AddRubric("BusPatch", fgkNofBusPath);
 237   fClusters->AddRubric("Manu", fgkNofManu);
 238   fClusters->AddRubric("Channel", AliMpConstants::ManuNofChannels());
 239   fClusters->Init();
 240
 241   fChamberTDHistList = new TList();
 242   fChamberTDHistList->SetOwner();
 243   fChamberTTHistList = new TList();
 244   fChamberTTHistList->SetOwner();
 245   fChamberSDHistList = new TList();
 246   fChamberSDHistList->SetOwner();
 247   fExtraHistList = new TList();
 248   fExtraHistList->SetOwner();
 249
 250   THnSparse *hn = 0x0;
 251   TString histName, histTitle;
 252
 253   // centrality bins
 254   Int_t nCentBins = 22;
 255   Double_t centRange[2] = {-5., 105.};
 256
 257   // prepare binning for THnSparse
 258   // 1: Ch or DE Id
 259   // 2: centrality
 260   // 3: pt
 261   // 4: y
 262   // 5: phi
 263   // 6: sign
 264   const Int_t nDims = 6;
 265   Int_t nBins[nDims] = {0, nCentBins, 20, 15, 15, 2};
 266   Double_t xMin[nDims] = {0., centRange[0], 0., -4., 0., -2.};
 267   Double_t xMax[nDims] = {0., centRange[1], 20., -2.5, TMath::TwoPi(), 2.};
 268
 269   for (Int_t iCh = 0; iCh < 10; iCh++)
 270   {
 271     // histograms per chamber
 272     nBins[0] = fgkNofDE[iCh];
 273     xMin[0] = 0.; xMax[0] = static_cast<Double_t>(fgkNofDE[iCh]);
 274     histTitle.Form("ChamberNbr %d", iCh+1);
 275     histName.Form("TD_ChamberNbr%d", iCh+1);
 276     hn = new THnSparseT<TArrayF>(histName, histTitle, nDims, nBins, xMin, xMax);
 277     fChamberTDHistList->AddAt(hn, iCh);
 278     histName.Form("TT_ChamberNbr%d",iCh+1);
 279     hn = new THnSparseT<TArrayF>(histName, histTitle, nDims, nBins, xMin, xMax);
 280     fChamberTTHistList->AddAt(hn, iCh);
 281     histName.Form("SD_ChamberNbr%d", iCh+1);
 282     hn = new THnSparseT<TArrayF>(histName, histTitle, nDims, nBins, xMin, xMax);
 283     fChamberSDHistList->AddAt(hn, iCh);
 284
 285   }
 286
 287   // global histograms per chamber
 288   nBins[0] = 10;
 289   xMin[0] = 0.5; xMax[0] = 10.5;
 290   hn = new THnSparseT<TArrayF>("TD_Chambers 11", "Chambers 11", nDims, nBins, xMin, xMax);
 291   fChamberTDHistList->AddAt(hn, 10);
 292   hn = new THnSparseT<TArrayF>("TT_Chambers 11", "Chambers 11", nDims, nBins, xMin, xMax);
 293   fChamberTTHistList->AddAt(hn, 10);
 294   hn = new THnSparseT<TArrayF>("SD_Chambers 11", "Chambers 11", nDims, nBins, xMin, xMax);
 295   fChamberSDHistList->AddAt(hn, 10);
 296
 297   //Extra histograms
 298   TH1F *fHistCent = new TH1F("fHistCent", "centrality distribution", nCentBins, centRange[0], centRange[1]);
 299   fExtraHistList->AddAt(fHistCent,0);
 300   TH1F *fHistPt = new TH1F("fHistPt", "pt distribution", 250, 0., 50.);
 301   fExtraHistList->AddAt(fHistPt,1);
 302   TH1F *fHistY = new TH1F("fHistY", "y distribution", 60, -4., -2.5);
 303   fExtraHistList->AddAt(fHistY,2);
 304   TH1F *fHistTheta = new TH1F("fHistTheta", "theta distribution", 120, 2.8, 3.2);
 305   fExtraHistList->AddAt(fHistTheta,3);
 306   TH1F *fHistP = new TH1F("fHistP", "momentum distribution", 250, 0., 500.);
 307   fExtraHistList->AddAt(fHistP,4);
 308   TH1F *fHistZ = new TH1F("fHistZ", "Z distribution", 200, -100., 100.);
 309   fExtraHistList->AddAt(fHistZ,5);
 310   TH1F *fHistPhi = new TH1F("fHistPhi", "phi distribution", 60, 0., TMath::TwoPi());
 311   fExtraHistList->AddAt(fHistPhi,6);
 312   TH1F *fHistPtRap2p5To2p75 = new TH1F("fHistPtRap2p5To2p75", "2.5 < y < 2.75", 250, 0., 50.);
 313   fExtraHistList->AddAt(fHistPtRap2p5To2p75,7);
 314   TH1F *fHistPtRap2p75To3p0 = new TH1F("fHistPtRap2p75To3p0", "2.75 < y < 3.0", 250, 0., 50.);
 315   fExtraHistList->AddAt(fHistPtRap2p75To3p0,8);
 316   TH1F *fHistPtRap3p0To3p25 = new TH1F("fHistPtRap3p0To3p25", "3.0 < y < 3.25", 250, 0., 50.);
 317   fExtraHistList->AddAt(fHistPtRap3p0To3p25,9);
 318   TH1F *fHistPtRap3p25To3p5 = new TH1F("fHistPtRap3p25To3p5", "3.25 < y < 3.5", 250, 0., 50.);
 319   fExtraHistList->AddAt(fHistPtRap3p25To3p5,10);
 320   TH1F *fHistPtRap3p5To3p75 = new TH1F("fHistPtRap3p5To3p75", "3.5 < y < 3.75", 250, 0., 50.);
 321   fExtraHistList->AddAt(fHistPtRap3p5To3p75,11);
 322   TH1F *fHistPtRap3p75To4p0 = new TH1F("fHistPtRap3p75To4p0", "3.75 < y < 4.0", 250, 0., 50.);
 323   fExtraHistList->AddAt(fHistPtRap3p75To4p0,12);
 324   TH1F *fHistRapPt1p0To2p0 = new TH1F("fHistRapPt1p0To2p0", "1.0 < pT < 2.0", 60, -4., -2.5);
 325   fExtraHistList->AddAt(fHistRapPt1p0To2p0,13);
 326   TH1F *fHistRapPt2p0To5p0 = new TH1F("fHistRapPt2p0To5p0", "2.0 < pT < 5.0", 60, -4., -2.5);
 327   fExtraHistList->AddAt(fHistRapPt2p0To5p0,14);
 328   TH1F *fHistRapPt5p0To8p0 = new TH1F("fHistRapPt5p0To8p0", "5.0 < pT < 8.0", 60, -4., -2.5);
 329   fExtraHistList->AddAt(fHistRapPt5p0To8p0,15);
 330   TH1F *fHistRapPt2p0To4p0 = new TH1F("fHistRapPt2p0To4p0", "2.0 < pT < 4.0", 60, -4., -2.5);
 331   fExtraHistList->AddAt(fHistRapPt2p0To4p0,16);
 332   TH1F *fHistRapPt4p0To8p0 = new TH1F("fHistRapPt4p0To8p0", "4.0 < pT < 8.0", 60, -4., -2.5);
 333   fExtraHistList->AddAt(fHistRapPt4p0To8p0,17);
 334
 335   TH2D *hDXYOverDXYMax = new TH2D("hDXYOverDXYMax", "DXY / DXYMax;DX / DXMax;DY / DYMax", 100, -1., 1., 100, -1., 1.);
 336   fExtraHistList->AddAt(hDXYOverDXYMax,18);
 337   TH2D *hDXOverDXMax = new TH2D("hDXOverDXMax", "DX / DXMax vs pXZ;pXZ;DX / DXMax", 50, 0., 500., 100, -1., 1.);
 338   fExtraHistList->AddAt(hDXOverDXMax,19);
 339   TH2D *hDYOverDYMax = new TH2D("hDYOverDYMax", "DY / DYMax vs pYZ;pYZ;DY / DYMax", 50, 0., 500., 100, -1., 1.);
 340   fExtraHistList->AddAt(hDYOverDYMax,20);
 341
 342   // post the output data at least once
 343   PostData(1, fClusters);
 344   PostData(2, fChamberTDHistList);
 345   PostData(3, fChamberTTHistList);
 346   PostData(4, fChamberSDHistList);
 347   PostData(5, fExtraHistList);
 348 }
 349
 350 //________________________________________________________________________
 351 void AliAnalysisTaskMuonTrackingEff::UserExec(Option_t *)
 352 {
 353   /// Main event loop
 354
 355   // check the OCDB has been loaded properly
 356   if (!fOCDBLoaded) AliFatal("Problem occur while loading OCDB objects");
 357
 358   // get the current event
 359   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 360   if (!esd) return;
 361
 362   // get the centrality
 363   Double_t cent = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
 364   if (cent <= fCentMin || cent > fCentMax) return;
 365   static_cast<TH1F*>(fExtraHistList->At(0))->Fill(cent);
 366
 367   // loop over tracks
 368   AliMUONTrack track;
 369   Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
 370   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 371     AliESDMuonTrack *esdTrack = esd->GetMuonTrack(iTrack);
 372
 373     // apply track selections
 374     Double_t pT = esdTrack->Pt();
 375     if (!esdTrack->ContainTrackerData() || !fMuonTrackCuts->IsSelected(esdTrack) || (fPtCut > 0. && pT < fPtCut) ||
 376         (fUseMCLabel && (esdTrack->GetLabel() < 0 || esdTrack->TestBit(BIT(22))))) continue;
 377
 378     // fill histograms
 379     Double_t y = esdTrack->Y();
 380     Double_t phi = esdTrack->Phi();
 381     static_cast<TH1F*>(fExtraHistList->At(1))->Fill(pT);
 382     static_cast<TH1F*>(fExtraHistList->At(2))->Fill(y);
 383     static_cast<TH1F*>(fExtraHistList->At(3))->Fill(esdTrack->Theta());
 384     static_cast<TH1F*>(fExtraHistList->At(4))->Fill(esdTrack->P());
 385     static_cast<TH1F*>(fExtraHistList->At(5))->Fill(esdTrack->GetZ());
 386     static_cast<TH1F*>(fExtraHistList->At(6))->Fill(phi);
 387     if ( (y < -2.5) && (y > -2.75) ) static_cast<TH1F*>(fExtraHistList->At(7))->Fill(pT);
 388     if ( (y < -2.75) && (y > -3.0) ) static_cast<TH1F*>(fExtraHistList->At(8))->Fill(pT);
 389     if ( (y < -3.0) && (y > -3.25) ) static_cast<TH1F*>(fExtraHistList->At(9))->Fill(pT);
 390     if ( (y < -3.25) && (y > -3.5) ) static_cast<TH1F*>(fExtraHistList->At(10))->Fill(pT);
 391     if ( (y < -3.5) && (y > -3.75) ) static_cast<TH1F*>(fExtraHistList->At(11))->Fill(pT);
 392     if ( (y < -3.75) && (y > -4.0) ) static_cast<TH1F*>(fExtraHistList->At(12))->Fill(pT);
 393     if ( (pT > 1.0) && (pT < 2.0) ) static_cast<TH1F*>(fExtraHistList->At(13))->Fill(y);
 394     if ( (pT > 2.0) && (pT < 5.0) ) static_cast<TH1F*>(fExtraHistList->At(14))->Fill(y);
 395     if ( (pT > 5.0) && (pT < 8.0) ) static_cast<TH1F*>(fExtraHistList->At(15))->Fill(y);
 396     if ( (pT > 2.0) && (pT < 4.0) ) static_cast<TH1F*>(fExtraHistList->At(16))->Fill(y);
 397     if ( (pT > 4.0) && (pT < 8.0) ) static_cast<TH1F*>(fExtraHistList->At(17))->Fill(y);
 398
 399     // convert to MUON track
 400     AliMUONESDInterface::ESDToMUON(*esdTrack, track);
 401     Double_t trackInfo[6] = {0., cent, pT, y, phi, static_cast<Double_t>(esdTrack->Charge())};
 402
 403     // tag the removable clusters/chambers, i.e. not needed to fulfill the tracking conditions
 404     Bool_t removableChambers[10];
 405     Bool_t isValidTrack = TagRemovableClusters(track, removableChambers);
 406
 407     // loop over clusters
 408     Int_t previousCh = -1;
 409     TObjArray *trackParams = track.GetTrackParamAtCluster();
 410     AliMUONTrackParam *trackParam = static_cast<AliMUONTrackParam*>(trackParams->First());
 411     while (trackParam) {
 412
 413       AliMUONVCluster* cluster = trackParam->GetClusterPtr();
 414       Int_t currentCh = cluster->GetChamberId();
 415       Int_t currentDE = cluster->GetDetElemId();
 416
 417       // find the pads at the position of the cluster
 418       Double_t pos[3] = {cluster->GetX(), cluster->GetY(), cluster->GetZ()};
 419       AliMpPad pad[2];
 420       FindPads(currentDE, pos, pad);
 421
 422       AliMUONTrackParam *nextTrackParam = static_cast<AliMUONTrackParam*>(trackParams->After(trackParam));
 423       Int_t nextCh = nextTrackParam ? nextTrackParam->GetClusterPtr()->GetChamberId() : 10;
 424
 425       // record all clusters/chambers
 426       RecordCluster(currentCh, currentDE, pad, trackInfo, "Detected", fChamberSDHistList, (currentCh != nextCh));
 427
 428       // record removable clusters/chambers
 429       if (trackParam->IsRemovable()) {
 430         Bool_t recordChamber = (removableChambers[currentCh] && currentCh != nextCh);
 431         RecordCluster(currentCh, currentDE, pad, trackInfo, "Accepted", fChamberTDHistList, recordChamber);
 432         RecordCluster(currentCh, currentDE, pad, trackInfo, "Expected", fChamberTTHistList, recordChamber);
 433       }
 434
 435       // record missing clusters/chambers prior to current one
 436       while (previousCh < currentCh-1 && (isValidTrack || previousCh < 5))
 437         FindAndRecordMissingClusters(*trackParam, ++previousCh, trackInfo);
 438
 439       // stop if we reached station 4 or 5 and the track is not valid
 440       if (!isValidTrack && currentCh > 5) break;
 441
 442       // record missing cluster on the same chamber
 443       if (currentCh != previousCh && currentCh != nextCh)
 444         FindAndRecordMissingClusters(*trackParam, currentCh, trackInfo);
 445
 446       if (nextTrackParam) {
 447
 448         // prepare next step
 449         previousCh = currentCh;
 450         trackParam = nextTrackParam;
 451
 452       } else {
 453
 454         // record missing clusters/chambers next to the last chamber
 455         while (++currentCh < 10 && (isValidTrack || currentCh < 6))
 456           FindAndRecordMissingClusters(*trackParam, currentCh, trackInfo);
 457
 458         break;
 459       }
 460
 461     }
 462
 463   }
 464
 465   // post the output data:
 466   PostData(1, fClusters);
 467   PostData(2, fChamberTDHistList);
 468   PostData(3, fChamberTTHistList);
 469   PostData(4, fChamberSDHistList);
 470   PostData(5, fExtraHistList);
 471 }
 472
 473 //________________________________________________________________________
 474 void AliAnalysisTaskMuonTrackingEff::Terminate(Option_t *)
 475 {
 476   /// final plots
 477
 478   if (!fEnableDisplay) return;
 479
 480   fClusters = static_cast<AliCounterCollection*>(GetOutputData(1));
 481   if (!fClusters) return;
 482
 483   // load mapping locally if not already done
 484   AliCDBManager* man = AliCDBManager::Instance();
 485   if (!man->IsDefaultStorageSet()) {
 486     if (gROOT->IsBatch()) man->SetDefaultStorage("alien://folder=/alice/simulation/2008/v4-15-Release/Ideal");
 487     else man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 488   }
 489   if (man->GetRun() < 0) man->SetRun(0);
 490   if (!AliMpSegmentation::Instance(kFALSE) || !AliMpDDLStore::Instance(kFALSE)) {
 491     if (!AliMUONCDB::LoadMapping()) return;
 492   }
 493
 494   TString clusterKey[3] = {"Detected", "Accepted", "Expected"};
 495
 496   // list of fired DEs
 497   TObjArray *deKeys = fClusters->GetKeyWords("DE").Tokenize(",");
 498   Int_t nDEs = deKeys->GetEntriesFast();
 499
 500   // loop over cluster types
 501   for (Int_t iKey = 0; iKey < 3; iKey++) {
 502
 503     // create the cluster store
 504     AliMUON2DMap clustersStore(kTRUE);
 505
 506     // loop over fired DEs
 507     for (Int_t iDE = 0; iDE < nDEs; iDE++) {
 508
 509       // get DE Id
 510       TString deKey = static_cast<TObjString*>(deKeys->UncheckedAt(iDE))->GetString();
 511       Int_t deId = deKey.Atoi();
 512
 513       // get numbers of clusters in this DE for each manu/channel combination
 514       TH2D *channelVsManu = fClusters->Get("channel", "Manu",
 515                                            Form("Cluster:%s/DE:%s", clusterKey[iKey].Data(), deKey.Data()));
 516       Int_t nManus = channelVsManu->GetNbinsX();
 517       Int_t nChannels = channelVsManu->GetNbinsY();
 518
 519       // loop over fired manus
 520       for (Int_t iManu = 1; iManu <= nManus; iManu++) {
 521
 522         // get manu Id
 523         TString manuKey = channelVsManu->GetXaxis()->GetBinLabel(iManu);
 524         Int_t manuId = manuKey.Atoi();
 525
 526         // loop over fired channels
 527         for (Int_t iChannel = 1; iChannel <= nChannels; iChannel++) {
 528
 529           // get channel Id
 530           TString channelKey = channelVsManu->GetYaxis()->GetBinLabel(iChannel);
 531           Int_t channelId = channelKey.Atoi();
 532
 533           // get the number of clusters in this pad
 534           Int_t nClusters = static_cast<Int_t>(channelVsManu->GetBinContent(iManu, iChannel));
 535           if (nClusters < 1) continue;
 536
 537           // register the clusters
 538           AliMUONVCalibParam* c = static_cast<AliMUONVCalibParam*>(clustersStore.FindObject(deId, manuId));
 539           if (!c) {
 540             c = new AliMUONCalibParamNI(1, AliMpConstants::ManuNofChannels(), deId, manuId);
 541             clustersStore.Add(c);
 542           }
 543           c->SetValueAsInt(channelId, 0, nClusters);
 544
 545         }
 546
 547       }
 548
 549       // clean memory
 550       delete channelVsManu;
 551
 552     }
 553
 554     // create the tracker data
 555     TString suffix = GetName();
 556     suffix.ReplaceAll("MuonTrackingEfficiency","");
 557     AliMUONTrackerData clustersData(Form("%sClusters%s", clusterKey[iKey].Data(), suffix.Data()),
 558                                     Form("%s clusters %s", clusterKey[iKey].Data(), suffix.Data()), 1, kFALSE);
 559     clustersData.SetDimensionName(0, "count");
 560     clustersData.Add(clustersStore);
 561
 562     // save it to a file
 563     TFile *outFile = TFile::Open("DisplayResults.root", "UPDATE");
 564     if (outFile && outFile->IsOpen()) {
 565       clustersData.Write(0x0, TObject::kOverwrite);
 566       outFile->Close();
 567     }
 568
 569   }
 570
 571   // clean memory
 572   delete deKeys;
 573
 574 }
 575
 576 //________________________________________________________________________
 577 Bool_t AliAnalysisTaskMuonTrackingEff::TagRemovableClusters(AliMUONTrack &track, Bool_t removableChambers[10])
 578 {
 579   /// Identify clusters/chambers that can be removed from the track
 580   /// return kTRUE if the track as it is satisfies the tracking conditions
 581
 582   for (Int_t i = 0; i < 10; i++) removableChambers[i] = kFALSE;
 583
 584   // check if track is valid as it is
 585   UInt_t requestedStationMask = AliMUONESDInterface::GetTracker()->GetRecoParam()->RequestedStationMask();
 586   Bool_t request2ChInSameSt45 = !AliMUONESDInterface::GetTracker()->GetRecoParam()->MakeMoreTrackCandidates();
 587   Bool_t isValidTrack = track.IsValid(requestedStationMask, request2ChInSameSt45);
 588
 589   // count the number of clusters per chamber and the number of chambers hit per station
 590   Int_t nClInCh[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 591   Int_t nChInSt[5] = {0, 0, 0, 0, 0};
 592   Int_t previousCh = -1;
 593   Int_t nClusters = track.GetNClusters();
 594   for (Int_t i = 0; i < nClusters; i++) {
 595
 596     AliMUONTrackParam *trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(i));
 597     Int_t currentCh = trackParam->GetClusterPtr()->GetChamberId();
 598     Int_t currentSt = currentCh/2;
 599
 600     nClInCh[currentCh]++;
 601
 602     if (currentCh != previousCh) {
 603       previousCh = currentCh;
 604       nChInSt[currentSt]++;
 605     }
 606
 607   }
 608
 609   // tag removable clusters/chambers
 610   for (Int_t i = 0; i < nClusters; i++) {
 611
 612     AliMUONTrackParam *trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(i));
 613     Int_t currentCh = trackParam->GetClusterPtr()->GetChamberId();
 614     Int_t currentSt = currentCh/2;
 615
 616     // for stations 1, 2 and 3
 617     if (currentSt < 3) {
 618
 619       if (!((1 << currentSt) & requestedStationMask) || // station not required
 620           nChInSt[currentSt] == 2) { // or both chamber hit in the station
 621
 622         removableChambers[currentCh] = kTRUE;
 623         trackParam->SetRemovable(kTRUE);
 624
 625       } else if (nClInCh[currentCh] == 2) { // 2 clusters in the chamber
 626
 627         trackParam->SetRemovable(kTRUE);
 628
 629       }
 630
 631     } else { // for stations 4 and 5
 632
 633       // if the track is already not valid we can certainly not remove more cluster on station 4 or 5
 634       if (!isValidTrack) continue;
 635
 636       if (((request2ChInSameSt45 && // tracking starts with 2 chamber hits in the same station:
 637             nChInSt[4-currentSt/4] == 2) || // --> 2 hits in the other stations
 638            (!request2ChInSameSt45 && // or tracking starts with 2 chamber hits in stations 4&5:
 639             nChInSt[3]+nChInSt[4] >= 3)) && // --> at least 3 hits in stations 4&5
 640           (!((1 << currentSt) & requestedStationMask) || // + this station not requested
 641            nChInSt[currentSt] == 2)) { // or 2 hits in it
 642
 643         removableChambers[currentCh] = kTRUE;
 644         trackParam->SetRemovable(kTRUE);
 645
 646       } else if (nClInCh[currentCh] == 2) { // 2 clusters in the chamber
 647
 648         trackParam->SetRemovable(kTRUE);
 649
 650       }
 651
 652     }
 653
 654   }
 655
 656   return isValidTrack;
 657
 658 }
 659
 660 //________________________________________________________________________
 661 void AliAnalysisTaskMuonTrackingEff::FindAndRecordMissingClusters(AliMUONTrackParam &param, Int_t chamber,
 662                                                                   Double_t trackInfo[6])
 663 {
 664   /// Find which detection elements should have been hit and record the missing clusters
 665
 666   static const Double_t maxDZ = 0.01; // max distance between extrapolated track and DE to stop to extrapolate
 667   static const Double_t maxDevX = 0.01; // max X-deviation per dZ(cm) between the linear and the correct extrapolation
 668   static const Double_t maxDevY = 0.1; // max Y-deviation per dZ(cm) between the linear and the correct extrapolation
 669   static const Double_t minDXY = 1.; // min half-size of track area to intersect with DE to account for bad DE area
 670
 671   Bool_t missingChamber = (chamber != param.GetClusterPtr()->GetChamberId());
 672   Int_t startDE = param.GetClusterPtr()->GetDetElemId();
 673   Double_t pos[3], maxDX, maxDY, dZ = 0.;
 674   AliMUONTrackParam param1, param2;
 675   AliMpPad pad[2];
 676
 677   // extrapolate the track to the missing chamber if needed
 678   param1.SetParameters(param.GetParameters());
 679   param1.SetZ(param.GetZ());
 680   if (missingChamber && !AliMUONTrackExtrap::ExtrapToZ(&param1, AliMUONConstants::DefaultChamberZ(chamber))) return;
 681
 682   Double_t pX = param1.Px();
 683   Double_t pY = param1.Py();
 684   Double_t pZ = param1.Pz();
 685   Double_t pXZ = TMath::Sqrt(pX*pX + pZ*pZ);
 686   Double_t pYZ = TMath::Sqrt(pY*pY + pZ*pZ);
 687
 688   // loop over DEs
 689   AliMpDEIterator it;
 690   it.First(chamber);
 691   while (!it.IsDone()) {
 692     Int_t deId = it.CurrentDEId();
 693
 694     // skip current cluster
 695     if (deId == startDE) {
 696       it.Next();
 697       continue;
 698     }
 699
 700     // reset parameters
 701     param2.SetParameters(param1.GetParameters());
 702     param2.SetZ(param1.GetZ());
 703
 704     // check if the track can cross this DE
 705     Int_t nStep = 0;
 706     Bool_t crossDE = kFALSE;
 707     Bool_t extrapOk = kTRUE;
 708     do {
 709
 710       // plots to check that the correct extrapolation is within the area opened around the linear extrapolation
 711       if (nStep >= 1) {
 712         Double_t dX = pos[0]-param2.GetNonBendingCoor();
 713         Double_t dY = pos[1]-param2.GetBendingCoor();
 714         Double_t dXOverDXMax = (dZ > 0.) ? dX*pXZ/(maxDevX*dZ) : 0.;
 715         Double_t dYOverDYMax = (dZ > 0.) ? dY*pYZ/(maxDevY*dZ) : 0.;
 716         static_cast<TH2D*>(fExtraHistList->At(18))->Fill(dXOverDXMax, dYOverDYMax);
 717         static_cast<TH2D*>(fExtraHistList->At(19))->Fill(pXZ, dXOverDXMax);
 718         static_cast<TH2D*>(fExtraHistList->At(20))->Fill(pYZ, dYOverDYMax);
 719         if (dXOverDXMax >= 1.) printf("st = %d; pXZ = %f; dZ = %f; dX = %f; dX*PXZ/(dZ*maxDevX) = %f\n",
 720                                       chamber/2, pXZ, dZ, dX, dXOverDXMax);
 721         if (dYOverDYMax >= 1.) printf("st = %d; pYZ = %f; dZ = %f; dY = %f; dY*PYZ/(dZ*maxDevY) = %f\n",
 722                                       chamber/2, pYZ, dZ, dY, dYOverDYMax);
 723       }
 724       nStep++;
 725
 726       // build the area in which the track can eventually cross this DE and check if it overlaps with the DE area
 727       Intersect(param2, deId, pos);
 728       dZ = TMath::Abs(pos[2]-param2.GetZ());
 729       maxDX = minDXY + maxDevX*dZ/pXZ;
 730       maxDY = minDXY + maxDevY*dZ/pYZ;
 731       AliMpArea area(pos[0], pos[1], maxDX, maxDY);
 732       crossDE = OverlapDE(area, deId);
 733
 734     } while(crossDE && dZ > maxDZ && (extrapOk = AliMUONTrackExtrap::ExtrapToZ(&param2, pos[2])));
 735
 736     // find the pads (if any) at the position of the missing cluster and register it
 737     if (crossDE && extrapOk && FindPads(deId, pos, pad)) {
 738       RecordCluster(chamber, deId, pad, trackInfo, "Expected", fChamberTTHistList, missingChamber);
 739       missingChamber = kFALSE;
 740     }
 741
 742     it.Next();
 743   }
 744
 745 }
 746
 747 //________________________________________________________________________
 748 void AliAnalysisTaskMuonTrackingEff::Intersect(AliMUONTrackParam &param, Int_t deId, Double_t p[3])
 749 {
 750   /// Find the intersection point between the track (assuming straight line) and the DE in the global frame
 751
 752   Double_t pos[3] = {param.GetNonBendingCoor(), param.GetBendingCoor(), param.GetZ()};
 753   Double_t slope[2] = {param.GetNonBendingSlope(), param.GetBendingSlope()};
 754   TVectorD &plane = *(static_cast<TVectorD*>(fDEPlanes->UncheckedAt(deId)));
 755
 756   p[2] = (plane[3]*(slope[0]*pos[2]-pos[0]+plane[0]) + plane[4]*(slope[1]*pos[2]-pos[1]+plane[1]) + plane[5]*plane[2]) /
 757   (plane[3]*slope[0] + plane[4]*slope[1] + plane[5]);
 758   p[0] = slope[0]*(p[2]-pos[2]) + pos[0];
 759   p[1] = slope[1]*(p[2]-pos[2]) + pos[1];
 760
 761 }
 762
 763 //________________________________________________________________________
 764 Bool_t AliAnalysisTaskMuonTrackingEff::OverlapDE(AliMpArea &area, Int_t deId)
 765 {
 766   /// Check whether (global) area overlaps with the given DE
 767
 768   AliMpArea* globalDEArea = fTransformer->GetDEArea(deId);
 769   if (!globalDEArea) return kFALSE;
 770
 771   return area.Overlap(*globalDEArea);
 772
 773 }
 774
 775 //________________________________________________________________________
 776 void AliAnalysisTaskMuonTrackingEff::RecordCluster(Int_t chamber, Int_t deId, AliMpPad pad[2], Double_t trackInfo[6],
 777                                                    TString clusterKey, TList *chamberHistList, Bool_t recordChamber)
 778 {
 779   /// Register the cluster in the given stores
 780
 781   // register the pads
 782   for (Int_t iCath = 0; iCath < 2; iCath++) if (pad[iCath].IsValid()) {
 783     Int_t manuId = pad[iCath].GetManuId();
 784     Int_t busPatchId = AliMpDDLStore::Instance()->GetBusPatchId(deId,manuId);
 785     fClusters->Count(Form("Cluster:%s/Chamber:%d/DE:%d/BusPatch:%d/Manu:%d/Channel:%d",
 786                      clusterKey.Data(), chamber, deId, busPatchId, manuId, pad[iCath].GetManuChannel()));
 787   }
 788
 789   // register the DE
 790   trackInfo[0] = static_cast<Double_t>(deId%100);
 791   static_cast<THnSparse*>(chamberHistList->At(chamber))->Fill(trackInfo);
 792
 793   // register the chamber
 794   if (recordChamber) {
 795     trackInfo[0] = static_cast<Double_t>(chamber+1);
 796     static_cast<THnSparse*>(chamberHistList->At(10))->Fill(trackInfo);
 797   }
 798
 799 }
 800
 801 //________________________________________________________________________
 802 Bool_t AliAnalysisTaskMuonTrackingEff::FindPads(Int_t deId, Double_t pos[3], AliMpPad pad[2])
 803 {
 804   /// Look for pads at the cluster's location
 805
 806   static const AliMpPad emptyPad;
 807
 808   // compute the cluster position in the DE frame
 809   Double_t localPos[3];
 810   fTransformer->Global2Local(deId, pos[0], pos[1], pos[2], localPos[0], localPos[1], localPos[2]);
 811
 812   // find pads at this position
 813   for (Int_t iCath = 0; iCath < 2; iCath++) {
 814     const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->GetMpSegmentation(deId, AliMp::GetCathodType(iCath));
 815     if (seg) pad[iCath] = seg->PadByPosition(localPos[0], localPos[1], kFALSE);
 816     else pad[iCath] = emptyPad;
 817   }
 818
 819   return (pad[0].IsValid() || pad[1].IsValid());
 820
 821 }
 822