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