]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/MUON/dep/AliAnalysisTaskMuonTrackingEff.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGPP / MUON / dep / AliAnalysisTaskMuonTrackingEff.cxx
CommitLineData
128a8042 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
27de2dfb 16/* $Id$ */
17
128a8042 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.
0ed78680 20//Work on ESD only
128a8042 21//Author: Nicolas LE BRIS - SUBATECH Nantes
79720f5b 22// Modified by Matthieu LENHARDT - SUBATECH Nantes
e1ec553d 23// Modified by Antoine LARDEUX - SUBATECH Nantes
0ed78680 24// Modified by Philippe PILLOT - SUBATECH Nantes
128a8042 25
e1ec553d 26// ROOT includes
0ed78680 27#include <TROOT.h>
79720f5b 28#include <TList.h>
923d3c3c 29#include <THnSparse.h>
e1ec553d 30#include <TObjArray.h>
0ed78680 31#include <TObjString.h>
e1ec553d 32#include <TGeoGlobalMagField.h>
0ed78680 33#include <TVectorD.h>
34#include <TH1F.h>
35#include <TH2F.h>
36#include <TH2D.h>
37#include <TFile.h>
128a8042 38
e1ec553d 39// STEER includes
128a8042 40#include "AliESDEvent.h"
41#include "AliESDMuonTrack.h"
fc7a3fd3 42#include "AliGeomManager.h"
79720f5b 43#include "AliCDBManager.h"
0ed78680 44#include "AliInputEventHandler.h"
45#include "AliCounterCollection.h"
fc7a3fd3 46
e1ec553d 47// ANALYSIS includes
0ed78680 48#include "AliAnalysisDataSlot.h"
49#include "AliAnalysisDataContainer.h"
e1ec553d 50#include "AliAnalysisManager.h"
e1ec553d 51#include "AliCentrality.h"
0ed78680 52#include "AliAnalysisTaskMuonTrackingEff.h"
128a8042 53
54//MUON includes
e1ec553d 55#include "AliMUONCDB.h"
56#include "AliMUONESDInterface.h"
0ed78680 57#include "AliMUONVTrackReconstructor.h"
58#include "AliMUONRecoParam.h"
27f15548 59#include "AliMUONGeometryTransformer.h"
e1ec553d 60#include "AliMUONTrack.h"
61#include "AliMUONTrackParam.h"
62#include "AliMUONTrackExtrap.h"
63#include "AliMUONVCluster.h"
79720f5b 64#include "AliMUONConstants.h"
0ed78680 65#include "AliMUON2DMap.h"
66#include "AliMUONVCalibParam.h"
67#include "AliMUONCalibParamNI.h"
68#include "AliMUONTrackerData.h"
128a8042 69
e1ec553d 70//include MUON/mapping:
71#include "AliMpDEManager.h"
72#include "AliMpSegmentation.h"
73#include "AliMpVSegmentation.h"
74#include "AliMpPad.h"
0ed78680 75#include "AliMpArea.h"
76#include "AliMpDEIterator.h"
77#include "AliMpConstants.h"
78#include "AliMpDDLStore.h"
e1ec553d 79
128a8042 80ClassImp(AliAnalysisTaskMuonTrackingEff)
81
0ed78680 82const Int_t AliAnalysisTaskMuonTrackingEff::fgkNofDE[11] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26, 156};
83const Int_t AliAnalysisTaskMuonTrackingEff::fgkNofBusPath = 888;
84const Int_t AliAnalysisTaskMuonTrackingEff::fgkNofManu = 16828;
128a8042 85
86//________________________________________________________________________
e1ec553d 87AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff() :
88 AliAnalysisTaskSE(),
89 fOCDBLoaded(kFALSE),
90 fOCDBpath(""),
0ed78680 91 fAlignOCDBpath(""),
92 fRecoParamOCDBpath(""),
93 fCentMin(-FLT_MAX),
94 fCentMax(FLT_MAX),
95 fMuonTrackCuts(0x0),
923d3c3c 96 fPtCut(-1.),
97 fUseMCLabel(kFALSE),
0ed78680 98 fEnableDisplay(kFALSE),
fc7a3fd3 99 fTransformer(0x0),
0ed78680 100 fDEPlanes(0x0),
101 fClusters(0x0),
e1ec553d 102 fChamberTDHistList(0x0),
145eb7fe 103 fChamberTTHistList(0x0),
104 fChamberSDHistList(0x0),
105 fExtraHistList(0x0)
128a8042 106{
e1ec553d 107 /// Default constructor
128a8042 108}
128a8042 109
110//________________________________________________________________________
e1ec553d 111AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff(TString name) :
112 AliAnalysisTaskSE(name),
113 fOCDBLoaded(kFALSE),
114 fOCDBpath("raw://"),
0ed78680 115 fAlignOCDBpath(""),
116 fRecoParamOCDBpath(""),
117 fCentMin(-FLT_MAX),
118 fCentMax(FLT_MAX),
119 fMuonTrackCuts(0x0),
923d3c3c 120 fPtCut(-1.),
121 fUseMCLabel(kFALSE),
0ed78680 122 fEnableDisplay(kFALSE),
fc7a3fd3 123 fTransformer(0x0),
0ed78680 124 fDEPlanes(0x0),
125 fClusters(0x0),
e1ec553d 126 fChamberTDHistList(0x0),
145eb7fe 127 fChamberTTHistList(0x0),
128 fChamberSDHistList(0x0),
129 fExtraHistList(0x0)
128a8042 130{
e1ec553d 131 /// Constructor
132
0ed78680 133 // Output slots
134 DefineOutput(1, AliCounterCollection::Class());
e1ec553d 135 DefineOutput(2, TList::Class());
136 DefineOutput(3, TList::Class());
137 DefineOutput(4, TList::Class());
145eb7fe 138 DefineOutput(5, TList::Class());
128a8042 139}
140
e1ec553d 141//________________________________________________________________________
128a8042 142AliAnalysisTaskMuonTrackingEff::~AliAnalysisTaskMuonTrackingEff()
143{
e1ec553d 144 /// Destructor
145 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
0ed78680 146 delete fMuonTrackCuts;
147 delete fClusters;
e1ec553d 148 delete fChamberTDHistList;
149 delete fChamberTTHistList;
145eb7fe 150 delete fChamberSDHistList;
151 delete fExtraHistList;
e1ec553d 152 }
fc7a3fd3 153 delete fTransformer;
0ed78680 154 delete fDEPlanes;
128a8042 155}
156
128a8042 157//________________________________________________________________________
e1ec553d 158void AliAnalysisTaskMuonTrackingEff::NotifyRun()
128a8042 159{
e1ec553d 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();
145eb7fe 167 if (man->IsDefaultStorageSet()) printf("EfficiencyTask: CDB default storage already set!\n");
0ed78680 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 }
145eb7fe 173 if (man->GetRun() > -1) printf("EfficiencyTask: run number already set!\n");
e1ec553d 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
e1ec553d 185 // Mapping
0ed78680 186 if (!AliMpSegmentation::Instance(kFALSE) || !AliMpDDLStore::Instance(kFALSE)) {
187 if (!AliMUONCDB::LoadMapping()) return;
e1ec553d 188 }
189
0ed78680 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 }
145eb7fe 213 }
e1ec553d 214
0ed78680 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
e1ec553d 223 fOCDBLoaded = kTRUE;
128a8042 224}
225
128a8042 226//________________________________________________________________________
e1ec553d 227void AliAnalysisTaskMuonTrackingEff::UserCreateOutputObjects()
128a8042 228{
0ed78680 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
e1ec553d 241 fChamberTDHistList = new TList();
242 fChamberTDHistList->SetOwner();
243 fChamberTTHistList = new TList();
244 fChamberTTHistList->SetOwner();
145eb7fe 245 fChamberSDHistList = new TList();
246 fChamberSDHistList->SetOwner();
247 fExtraHistList = new TList();
248 fExtraHistList->SetOwner();
e1ec553d 249
923d3c3c 250 THnSparse *hn = 0x0;
923d3c3c 251 TString histName, histTitle;
252
253 // centrality bins
e1ec553d 254 Int_t nCentBins = 22;
255 Double_t centRange[2] = {-5., 105.};
923d3c3c 256
257 // prepare binning for THnSparse
258 // 1: Ch or DE Id
259 // 2: centrality
260 // 3: pt
261 // 4: y
0ed78680 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.};
e1ec553d 268
269 for (Int_t iCh = 0; iCh < 10; iCh++)
270 {
271 // histograms per chamber
0ed78680 272 nBins[0] = fgkNofDE[iCh];
273 xMin[0] = 0.; xMax[0] = static_cast<Double_t>(fgkNofDE[iCh]);
e1ec553d 274 histTitle.Form("ChamberNbr %d", iCh+1);
923d3c3c 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);
e1ec553d 284
e1ec553d 285 }
286
287 // global histograms per chamber
923d3c3c 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);
145eb7fe 296
297 //Extra histograms
923d3c3c 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);
0ed78680 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);
e1ec553d 341
342 // post the output data at least once
0ed78680 343 PostData(1, fClusters);
344 PostData(2, fChamberTDHistList);
345 PostData(3, fChamberTTHistList);
346 PostData(4, fChamberSDHistList);
347 PostData(5, fExtraHistList);
128a8042 348}
349
128a8042 350//________________________________________________________________________
e1ec553d 351void AliAnalysisTaskMuonTrackingEff::UserExec(Option_t *)
128a8042 352{
e1ec553d 353 /// Main event loop
354
355 // check the OCDB has been loaded properly
0ed78680 356 if (!fOCDBLoaded) AliFatal("Problem occur while loading OCDB objects");
e1ec553d 357
358 // get the current event
359 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
360 if (!esd) return;
361
362 // get the centrality
0ed78680 363 Double_t cent = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
364 if (cent <= fCentMin || cent > fCentMax) return;
365 static_cast<TH1F*>(fExtraHistList->At(0))->Fill(cent);
e1ec553d 366
367 // loop over tracks
368 AliMUONTrack track;
369 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
0ed78680 370 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
371 AliESDMuonTrack *esdTrack = esd->GetMuonTrack(iTrack);
e1ec553d 372
0ed78680 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;
e1ec553d 377
0ed78680 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);
923d3c3c 383 static_cast<TH1F*>(fExtraHistList->At(3))->Fill(esdTrack->Theta());
384 static_cast<TH1F*>(fExtraHistList->At(4))->Fill(esdTrack->P());
0ed78680 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);
e1ec553d 398
0ed78680 399 // convert to MUON track
e1ec553d 400 AliMUONESDInterface::ESDToMUON(*esdTrack, track);
8917f9d2 401 Double_t trackInfo[6] = {0., cent, pT, y, phi, static_cast<Double_t>(esdTrack->Charge())};
0ed78680 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 }
e1ec553d 462
e1ec553d 463 }
464
465 // post the output data:
0ed78680 466 PostData(1, fClusters);
467 PostData(2, fChamberTDHistList);
468 PostData(3, fChamberTTHistList);
469 PostData(4, fChamberSDHistList);
470 PostData(5, fExtraHistList);
79720f5b 471}
e1ec553d 472
79720f5b 473//________________________________________________________________________
e1ec553d 474void AliAnalysisTaskMuonTrackingEff::Terminate(Option_t *)
79720f5b 475{
e1ec553d 476 /// final plots
79720f5b 477
0ed78680 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;
e1ec553d 492 }
493
0ed78680 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++) {
145eb7fe 502
0ed78680 503 // create the cluster store
504 AliMUON2DMap clustersStore(kTRUE);
145eb7fe 505
0ed78680 506 // loop over fired DEs
507 for (Int_t iDE = 0; iDE < nDEs; iDE++) {
e1ec553d 508
0ed78680 509 // get DE Id
510 TString deKey = static_cast<TObjString*>(deKeys->UncheckedAt(iDE))->GetString();
511 Int_t deId = deKey.Atoi();
e1ec553d 512
0ed78680 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();
145eb7fe 518
0ed78680 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 }
e1ec553d 546
e1ec553d 547 }
548
0ed78680 549 // clean memory
550 delete channelVsManu;
551
79720f5b 552 }
e1ec553d 553
0ed78680 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
e1ec553d 574}
575
576//________________________________________________________________________
0ed78680 577Bool_t AliAnalysisTaskMuonTrackingEff::TagRemovableClusters(AliMUONTrack &track, Bool_t removableChambers[10])
e1ec553d 578{
0ed78680 579 /// Identify clusters/chambers that can be removed from the track
580 /// return kTRUE if the track as it is satisfies the tracking conditions
79720f5b 581
0ed78680 582 for (Int_t i = 0; i < 10; i++) removableChambers[i] = kFALSE;
e1ec553d 583
0ed78680 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++) {
e1ec553d 595
0ed78680 596 AliMUONTrackParam *trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(i));
597 Int_t currentCh = trackParam->GetClusterPtr()->GetChamberId();
598 Int_t currentSt = currentCh/2;
e1ec553d 599
0ed78680 600 nClInCh[currentCh]++;
e1ec553d 601
0ed78680 602 if (currentCh != previousCh) {
603 previousCh = currentCh;
604 nChInSt[currentSt]++;
605 }
e1ec553d 606
0ed78680 607 }
608
609 // tag removable clusters/chambers
610 for (Int_t i = 0; i < nClusters; i++) {
e1ec553d 611
0ed78680 612 AliMUONTrackParam *trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(i));
613 Int_t currentCh = trackParam->GetClusterPtr()->GetChamberId();
614 Int_t currentSt = currentCh/2;
e1ec553d 615
0ed78680 616 // for stations 1, 2 and 3
617 if (currentSt < 3) {
e1ec553d 618
0ed78680 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 }
e1ec553d 630
0ed78680 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 }
e1ec553d 651
e1ec553d 652 }
0ed78680 653
e1ec553d 654 }
0ed78680 655
656 return isValidTrack;
657
128a8042 658}
0246246b 659
e1ec553d 660//________________________________________________________________________
0ed78680 661void AliAnalysisTaskMuonTrackingEff::FindAndRecordMissingClusters(AliMUONTrackParam &param, Int_t chamber,
662 Double_t trackInfo[6])
e1ec553d 663{
0ed78680 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
145eb7fe 745}
746
747//________________________________________________________________________
0ed78680 748void AliAnalysisTaskMuonTrackingEff::Intersect(AliMUONTrackParam &param, Int_t deId, Double_t p[3])
145eb7fe 749{
0ed78680 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
e1ec553d 761}
0246246b 762
e1ec553d 763//________________________________________________________________________
0ed78680 764Bool_t AliAnalysisTaskMuonTrackingEff::OverlapDE(AliMpArea &area, Int_t deId)
e1ec553d 765{
0ed78680 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
e1ec553d 773}
79720f5b 774
e1ec553d 775//________________________________________________________________________
0ed78680 776void AliAnalysisTaskMuonTrackingEff::RecordCluster(Int_t chamber, Int_t deId, AliMpPad pad[2], Double_t trackInfo[6],
777 TString clusterKey, TList *chamberHistList, Bool_t recordChamber)
e1ec553d 778{
0ed78680 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
e1ec553d 799}
800
801//________________________________________________________________________
0ed78680 802Bool_t AliAnalysisTaskMuonTrackingEff::FindPads(Int_t deId, Double_t pos[3], AliMpPad pad[2])
e1ec553d 803{
0ed78680 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
0246246b 821}
79720f5b 822