]>
Commit | Line | Data |
---|---|---|
fa10255f | 1 | /************************************************************************* |
2 | * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /////////////////////////////////////////////////////////////////////////// | |
17 | // Class AliTaskITSUPerf // | |
18 | // Analysis task to produce data and MC histos needed for tracklets // | |
19 | // dNdEta extraction in multiple bins in one go // | |
20 | // Author: ruben.shahoyan@cern.ch // | |
21 | /////////////////////////////////////////////////////////////////////////// | |
22 | ||
23 | #include "TFile.h" | |
24 | #include "TChain.h" | |
25 | #include "TTree.h" | |
26 | #include "TRandom.h" | |
27 | #include "TH1F.h" | |
28 | #include "TH2F.h" | |
29 | #include "TList.h" | |
30 | #include "TObjArray.h" | |
31 | #include "TGeoGlobalMagField.h" | |
32 | ||
33 | #include "AliAnalysisManager.h" | |
34 | ||
35 | #include "AliESDEvent.h" | |
36 | #include "AliESDInputHandler.h" | |
37 | #include "AliESDInputHandlerRP.h" | |
38 | #include "AliCDBPath.h" | |
39 | #include "AliCDBManager.h" | |
40 | #include "AliCDBEntry.h" | |
41 | #include "AliCDBStorage.h" | |
42 | #include "AliGeomManager.h" | |
43 | #include "AliMagF.h" | |
44 | #include "AliMCEventHandler.h" | |
45 | #include "AliMCEvent.h" | |
46 | #include "AliMCParticle.h" | |
47 | #include "AliStack.h" | |
48 | #include "AliGenEventHeader.h" | |
49 | #include "AliCentrality.h" | |
50 | #include "AliLog.h" | |
51 | #include "AliPhysicsSelection.h" | |
52 | #include "AliGenEventHeader.h" | |
53 | #include "AliGenHijingEventHeader.h" | |
54 | #include "AliGenDPMjetEventHeader.h" | |
55 | #include "AliGenCocktailEventHeader.h" | |
56 | #include "AliESDtrackCuts.h" | |
57 | ||
58 | #include "AliTaskITSUPerf.h" | |
59 | #include "../ITS/UPGRADE/AliITSURecoDet.h" | |
60 | #include "../ITS/UPGRADE/AliITSUGeomTGeo.h" | |
61 | #include "../ITS/UPGRADE/AliITSUClusterPix.h" | |
62 | #include "../ITS/UPGRADE/AliITSUTrackCond.h" | |
7bc38dea | 63 | #include "../ITS/UPGRADE/AliITSUAux.h" |
fa10255f | 64 | #include <TGeoManager.h> |
7bc38dea | 65 | #include <TTree.h> |
fa10255f | 66 | |
67 | ||
68 | ClassImp(AliTaskITSUPerf) | |
69 | ||
70 | ||
71 | const char* AliTaskITSUPerf::fgkLabelTypes[AliTaskITSUPerf::kNLabelTypes] = { | |
72 | "ITSokTPCok" | |
73 | ,"ITSokTPCfk" | |
74 | ,"ITSfkTPCok" | |
75 | ,"ITSfkTPCfk" | |
76 | ,"ITSTPCmismatch" | |
77 | ,"ITSTPCnoMatch" | |
78 | }; | |
79 | ||
7bc38dea | 80 | |
81 | typedef struct { | |
82 | Bool_t prim; | |
83 | Bool_t rcbl; | |
84 | Char_t nClITS; | |
85 | Char_t nClTPC; | |
be806700 | 86 | Char_t nClITSMC; |
7bc38dea | 87 | Char_t mcCl[7]; |
88 | Char_t rcCl[7]; | |
89 | Char_t fcCl[7]; | |
be806700 | 90 | Char_t charge; |
b824a5f5 | 91 | Float_t zV; |
7bc38dea | 92 | Float_t ptMC; |
93 | Float_t etaMC; | |
94 | Float_t pt; | |
95 | Float_t eta; | |
b824a5f5 | 96 | Float_t ptTPC; |
97 | Float_t etaTPC; | |
7bc38dea | 98 | Float_t dcaR; |
99 | Float_t dcaZ; | |
100 | Float_t dcaRE; | |
101 | Float_t dcaZE; | |
699d8344 | 102 | Float_t ptE; |
b824a5f5 | 103 | Float_t ptTPCE; |
7bc38dea | 104 | Float_t phi; |
105 | Int_t pdg; | |
106 | Int_t lbl; | |
b824a5f5 | 107 | Int_t lblTPC; |
be806700 | 108 | Int_t spl; |
b824a5f5 | 109 | Int_t evID; |
7bc38dea | 110 | } trInfo_t; |
111 | ||
112 | trInfo_t trackInfo; | |
113 | ||
fa10255f | 114 | //________________________________________________________________________ |
115 | /*//Default constructor | |
116 | AliTaskITSUPerf::AliTaskITSUPerf(const char *name) | |
117 | : AliAnalysisTaskSE(name), | |
118 | */ | |
119 | //________________________________________________________________________ | |
120 | AliTaskITSUPerf::AliTaskITSUPerf(const char *name) | |
121 | : AliAnalysisTaskSE(name) | |
122 | // | |
123 | ,fOutput(0) | |
124 | ,fHistosCentMCLb(0) | |
125 | ,fHistosCent(0) | |
126 | ,fRPTree(0) | |
127 | ,fStack(0) | |
128 | ,fMCEvent(0) | |
129 | ,fVtxSPD(0) | |
130 | ,fVtxTrc(0) | |
131 | ,fESDEvent(0) | |
132 | ,fUseSpecialOutput(kTRUE) | |
133 | ,fUseMC(kTRUE) | |
134 | ,fTrackingCond(0) | |
135 | // | |
136 | ,fGeom(0) | |
137 | ,fITS(0) | |
138 | ,fNSelTracksMC(0) | |
139 | ,fMCStatus(0) | |
097e4cdf | 140 | // |
141 | #ifdef _CLUS_LIST_ | |
142 | ,fClLists(0) | |
143 | #endif | |
144 | // | |
fa10255f | 145 | ,fNPtBins(20) |
146 | ,fNResBins(100) | |
097e4cdf | 147 | ,fMinTPCclusters(60) |
fa10255f | 148 | ,fPtMin(0) |
149 | ,fPtMax(10.) | |
150 | ,fEtaMin(-3.) | |
151 | ,fEtaMax( 3.) | |
152 | ,fZVertexMin(-15.) | |
153 | ,fZVertexMax(15.) | |
154 | // | |
155 | ,fCurrCentBin(0) | |
156 | ,fNCentBins(1) | |
157 | ,fUseCentralityVar(0) | |
60b14b59 | 158 | ,fTPCCut(0) |
7bc38dea | 159 | ,fTree(0) |
fa10255f | 160 | { |
161 | // Constructor | |
162 | ||
163 | DefineOutput(1, TList::Class()); | |
164 | // | |
165 | } | |
166 | ||
167 | //________________________________________________________________________ | |
168 | AliTaskITSUPerf::~AliTaskITSUPerf() | |
169 | { | |
170 | // Destructor | |
171 | // histograms are in the output list and deleted when the output | |
172 | // list is deleted by the TSelector dtor | |
173 | if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { //RRR | |
174 | printf("Deleteing output\n"); | |
175 | delete fOutput; | |
176 | fOutput = 0; | |
177 | } | |
178 | // | |
179 | fHistosCentMCLb.Clear(); | |
180 | fHistosCent.Clear(); | |
181 | // | |
182 | delete fITS; | |
183 | delete fGeom; | |
60b14b59 | 184 | delete fTPCCut; |
fa10255f | 185 | // |
186 | } | |
187 | ||
188 | //________________________________________________________________________ | |
189 | void AliTaskITSUPerf::UserCreateOutputObjects() | |
190 | { | |
191 | // | |
192 | if (fUseSpecialOutput) OpenFile(1); | |
193 | fOutput = new TList(); | |
194 | fOutput->SetOwner(); | |
195 | // | |
196 | // init geometry | |
197 | if (!gGeoManager) { | |
198 | TString geom = "geometry.root"; | |
199 | AliGeomManager::LoadGeometry(geom.Data()); | |
200 | if (!gGeoManager) AliFatal("Failed to load geometry"); | |
201 | } | |
202 | // | |
203 | // create ITS interface | |
204 | fGeom = new AliITSUGeomTGeo(kTRUE,kTRUE); | |
205 | AliITSUClusterPix::SetGeom(fGeom); | |
206 | fITS = new AliITSURecoDet(fGeom,"ITSURecoInterface"); | |
207 | fITS->CreateClusterArrays(); | |
208 | // | |
209 | // Create histograms | |
210 | for (int ib=0;ib<fNCentBins;ib++) { | |
211 | BookHistos(ib); | |
212 | } | |
213 | // | |
214 | int nhist = fOutput->GetEntries(); | |
215 | for (int i=0;i<nhist;i++) { | |
216 | TObject* hst = fOutput->At(i); | |
217 | if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue; | |
218 | ((TH1*)hst)->Sumw2(); | |
219 | } | |
220 | // | |
7bc38dea | 221 | fTree = new TTree("trinfo","trinfo"); |
222 | fTree->Branch("prim",&trackInfo.prim,"prim/O"); | |
223 | fTree->Branch("rcbl",&trackInfo.rcbl,"rcbl/O"); | |
224 | fTree->Branch("nClITS",&trackInfo.nClITS,"nClITS/b"); | |
225 | fTree->Branch("nClTPC",&trackInfo.nClTPC,"nClTPC/b"); | |
be806700 | 226 | fTree->Branch("nClITSMC",&trackInfo.nClITSMC,"nClITSMC/b"); |
7bc38dea | 227 | fTree->Branch("mcCl",&trackInfo.mcCl,"mcCl[7]/b"); |
228 | fTree->Branch("rcCl",&trackInfo.rcCl,"rcCl[7]/b"); | |
229 | fTree->Branch("fcCl",&trackInfo.fcCl,"fcCl[7]/b"); | |
be806700 | 230 | fTree->Branch("charge",&trackInfo.charge,"charge/B"); |
b824a5f5 | 231 | fTree->Branch("zV", &trackInfo.zV,"zV/F"); |
7bc38dea | 232 | fTree->Branch("ptMC", &trackInfo.ptMC,"ptMC/F"); |
233 | fTree->Branch("etaMC", &trackInfo.etaMC,"etaMC/F"); | |
234 | fTree->Branch("pt", &trackInfo.pt,"pt/F"); | |
235 | fTree->Branch("eta", &trackInfo.eta,"eta/F"); | |
b824a5f5 | 236 | fTree->Branch("ptTPC", &trackInfo.ptTPC,"ptTPC/F"); |
237 | fTree->Branch("etaTPC", &trackInfo.etaTPC,"etaTPC/F"); | |
7bc38dea | 238 | fTree->Branch("dcaR", &trackInfo.dcaR,"dcaR/F"); |
239 | fTree->Branch("dcaZ", &trackInfo.dcaZ,"dcaZ/F"); | |
240 | fTree->Branch("dcaRE", &trackInfo.dcaRE,"dcaRE/F"); | |
241 | fTree->Branch("dcaZE", &trackInfo.dcaZE,"dcaZE/F"); | |
699d8344 | 242 | fTree->Branch("ptE", &trackInfo.ptE,"ptE/F"); |
b824a5f5 | 243 | fTree->Branch("ptTPCE", &trackInfo.ptTPCE,"ptTPCE/F"); |
7bc38dea | 244 | fTree->Branch("phi", &trackInfo.phi,"phi/F"); |
b824a5f5 | 245 | fTree->Branch("pdg", &trackInfo.pdg,"pdg/I"); |
7bc38dea | 246 | fTree->Branch("lbl", &trackInfo.lbl,"lbl/I"); |
b824a5f5 | 247 | fTree->Branch("lblTPC", &trackInfo.lblTPC,"lblTPC/I"); |
be806700 | 248 | fTree->Branch("spl", &trackInfo.spl,"spl/I"); |
b824a5f5 | 249 | fTree->Branch("evID", &trackInfo.evID,"evID/I"); |
7bc38dea | 250 | fOutput->Add(fTree); |
251 | // | |
60b14b59 | 252 | fTPCCut = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); |
253 | fTPCCut->SetEtaRange(fEtaMin,fEtaMax); | |
254 | fTPCCut->SetPtRange(fPtMin,fPtMax); | |
255 | fTPCCut->SetPtRange(fPtMin,fPtMax); | |
256 | fTPCCut->SetMinNClustersTPC(fMinTPCclusters); | |
7bc38dea | 257 | printf("Cuts: Eta: %f %f | Pt: %f %f\n",fEtaMin,fEtaMax,fPtMin,fPtMax); |
60b14b59 | 258 | // |
097e4cdf | 259 | #ifdef _CLUS_LIST_ |
260 | fClLists.SetOwner(); | |
261 | #endif | |
fa10255f | 262 | PostData(1, fOutput); |
263 | // | |
264 | } | |
265 | ||
266 | //________________________________________________________________________ | |
267 | void AliTaskITSUPerf::UserExec(Option_t *) | |
268 | { | |
269 | // Main loop | |
270 | static int evID = 0; | |
271 | AliInfo(Form("Event: %d",evID++)); | |
272 | // | |
273 | AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager(); | |
274 | fRPTree = 0; | |
275 | AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler(); | |
b824a5f5 | 276 | if (!handRP) { AliFatal("No RP handler"); return; } |
fa10255f | 277 | // |
278 | fESDEvent = handRP->GetEvent(); | |
279 | if (!fESDEvent) { AliFatal("No AliESDEvent"); return; } | |
280 | // | |
281 | fRPTree = handRP->GetTreeR("ITS"); | |
b824a5f5 | 282 | if (!fRPTree) { AliWarning("Invalid ITS cluster tree"); } |
fa10255f | 283 | // |
284 | AliMCEventHandler* eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler(); | |
285 | if (!eventHandler) { AliFatal("Could not retrieve MC event handler"); return; } | |
286 | fMCEvent = eventHandler->MCEvent(); | |
287 | if (!fMCEvent) { AliFatal("Could not retrieve MC event"); return; } | |
288 | fStack = fMCEvent->Stack(); | |
289 | if (!fStack) { AliFatal("Stack is not available"); return; } | |
290 | // | |
291 | // MC Generator info | |
292 | AliGenEventHeader* mcGenH = 0; | |
293 | mcGenH = fMCEvent->GenEventHeader(); | |
294 | TArrayF vtmc(3); | |
295 | mcGenH->PrimaryVertex(vtmc); | |
296 | for (int i=3;i--;) fVtxMC[i] = vtmc[i]; | |
297 | // | |
298 | printf("MagField: %f\n",fESDEvent->GetMagneticField()); | |
299 | printf("VtxMC : %+e %+e %+e\n",fVtxMC[0],fVtxMC[1],fVtxMC[2]); | |
300 | // | |
301 | fVtxSPD = fESDEvent->GetPrimaryVertexSPD(); | |
302 | // | |
303 | if (fVtxSPD->GetStatus()) printf("VtxSPD: %+e %+e %+e | Nc=%d\n",fVtxSPD->GetX(),fVtxSPD->GetY(),fVtxSPD->GetZ(),fVtxSPD->GetNContributors()); | |
304 | else printf("VtxSPD: N/A\n"); | |
305 | // | |
306 | fVtxTrc = fESDEvent->GetPrimaryVertexTracks(); | |
307 | if (fVtxTrc->GetStatus()) printf("VtxTRC: %+e %+e %+e | Nc=%d\n",fVtxTrc->GetX(),fVtxTrc->GetY(),fVtxTrc->GetZ(),fVtxTrc->GetNContributors()); | |
308 | else printf("VtxTRC: N/A\n"); | |
309 | // | |
310 | BuildMCInfo(); | |
311 | CheckTracks(); | |
312 | // | |
313 | } | |
314 | ||
315 | ||
316 | //________________________________________________________________________ | |
317 | void AliTaskITSUPerf::Terminate(Option_t *) | |
318 | { | |
097e4cdf | 319 | #ifdef _CLUS_LIST_ |
320 | fClLists.Delete(); | |
321 | #endif | |
fa10255f | 322 | Printf("Terminating..."); |
323 | } | |
324 | ||
325 | //_________________________________________________________________________ | |
326 | void AliTaskITSUPerf::BookStandardHistosCentMCLb(Int_t bin, Int_t mcLb) | |
327 | { | |
328 | // book standard set of histos for specific centrality bin and MC status type, | |
329 | // adding them to output list and management array | |
330 | // | |
60b14b59 | 331 | const double kMaxDPt = 0.2; |
332 | const double kMaxDCA = 0.1; | |
fa10255f | 333 | |
334 | TH2* h2=0; | |
335 | // pt resolution | |
336 | h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"dPT2PTvsPT"), | |
337 | Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"#Deltap_{T}/p_{T} vs p_{T}"), | |
338 | fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDPt,kMaxDPt); | |
339 | AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResPTvsPTMC,mcLb,bin) ); | |
340 | // | |
341 | // transverse DCA resolution | |
342 | h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCARvsPT"), | |
343 | Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA R vs p_{T}"), | |
344 | fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA); | |
345 | AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCARvsPTMC,mcLb,bin) ); | |
346 | // | |
347 | // Z DCA resolution | |
348 | h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCAZvsPT"), | |
349 | Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA Z vs p_{T}"), | |
350 | fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA); | |
351 | AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCAZvsPTMC,mcLb,bin) ); | |
352 | // | |
353 | } | |
354 | ||
355 | //_________________________________________________________________________ | |
356 | void AliTaskITSUPerf::BookStandardHistosCent(Int_t bin) | |
357 | { | |
358 | // book standard set of histos for specific centrality bin | |
359 | // | |
360 | TH2* h2=0; | |
361 | // | |
362 | // MC labels combination vs pt | |
097e4cdf | 363 | h2 = new TH2F(Form("B%d_%s_RCBL",bin,"MCLabvsPT"), |
364 | Form("B%d %s Reconstructable",bin,"MCLabComb vs p_{T}"), | |
fa10255f | 365 | fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5); |
097e4cdf | 366 | AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusRcbl,-1,bin) ); |
fa10255f | 367 | for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]); |
368 | // | |
097e4cdf | 369 | h2 = new TH2F(Form("B%d_%s_NotRCBL_prim",bin,"MCLabvsPT"), |
370 | Form("B%d %s Not Reconstructable, Prim.",bin,"MCLabComb vs p_{T}"), | |
371 | fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5); | |
372 | AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusNRcblPrim,-1,bin) ); | |
373 | for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]); | |
374 | // | |
375 | h2 = new TH2F(Form("B%d_%s_NotRCBL_sec",bin,"MCLabvsPT"), | |
376 | Form("B%d %s Not Reconstructable, Sec.",bin,"MCLabComb vs p_{T}"), | |
377 | fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5); | |
378 | AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusNRcblSec,-1,bin) ); | |
379 | for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]); | |
7bc38dea | 380 | // |
381 | h2 = new TH2F(Form("B%d_MCLr_Prim",bin), | |
382 | Form("B%d MCLr Prim vs pt",bin), | |
383 | fNPtBins,fPtMin,fPtMax, 7,-0.5,6.5); | |
384 | AddHisto(&fHistosCent,h2, GetHistoID(kHMCLrPresPrim,-1,bin) ); | |
385 | // | |
386 | h2 = new TH2F(Form("B%d_MCLr_Sec",bin), | |
387 | Form("B%d MCLr Sec vs pt",bin), | |
388 | fNPtBins,fPtMin,fPtMax, 7,-0.5,6.5); | |
389 | AddHisto(&fHistosCent,h2, GetHistoID(kHMCLrPresSec,-1,bin) ); | |
390 | // | |
fa10255f | 391 | } |
392 | ||
393 | //_________________________________________________________________________ | |
394 | void AliTaskITSUPerf::BookHistos(Int_t bin) | |
395 | { | |
396 | // book standard set of histos, adding them to output list and management array | |
397 | // | |
398 | // standard histos | |
399 | for (int ilb=0;ilb<kNLabelTypes;ilb++) { // create similar set of histos for all kind of MC labels | |
400 | BookStandardHistosCentMCLb(bin,ilb); | |
401 | } | |
402 | BookStandardHistosCent(bin); // st.histos regardless of MClabels comb. | |
403 | } | |
404 | ||
405 | //_________________________________________________________________________ | |
406 | void AliTaskITSUPerf::AddHisto(TObjArray* array,TObject* h, Int_t at) | |
407 | { | |
408 | // add single histo to the set | |
409 | if (at>=0) array->AddAtAndExpand(h,at); | |
410 | else array->Add(h); | |
411 | fOutput->Add(h); | |
412 | } | |
413 | ||
414 | //_________________________________________________________________________ | |
415 | Int_t AliTaskITSUPerf::GetHistoID(Int_t htype, Int_t mcStat, Int_t centBin) const | |
416 | { | |
417 | // generate ID for the histo. htype must be <100, mcStat<kNLabelTypes | |
418 | // | |
7bc38dea | 419 | // if (mcStat>=kNLabelTypes) AliFatal(Form("MCStatus ID=%d > max.allowed %d",mcStat,kNLabelTypes)); |
fa10255f | 420 | if (centBin>=fNCentBins) AliFatal(Form("CentBin ID=%d > max.allowed %d",centBin,fNCentBins)); |
421 | // | |
422 | if (centBin>=0) { | |
423 | if (mcStat>=0) { // MCtruthness dependent histos | |
424 | if (htype>=kHNStdHistosCentMC) AliFatal(Form("StdCentMC histo ID=%d > max.allowed %d",htype,kHNStdHistosCentMC)); | |
425 | return htype + kHNStdHistosCentMC*(mcStat + kNLabelTypes*centBin); | |
426 | } | |
427 | else { | |
428 | if (htype>=kHNStdHistosCent) AliFatal(Form("StdCent histo ID=%d > max.allowed %d",htype,kHNStdHistosCent)); | |
429 | return htype + kHNStdHistosCent*centBin; | |
430 | } | |
431 | } | |
432 | // custom histo | |
433 | return htype; | |
434 | // | |
435 | } | |
436 | ||
437 | //_________________________________________________________________________ | |
438 | void AliTaskITSUPerf::CheckTracks() | |
439 | { | |
440 | // build mc truth info for tracks | |
097e4cdf | 441 | Bool_t dump = kFALSE; |
7bc38dea | 442 | Bool_t dump1 = kFALSE; |
b824a5f5 | 443 | static int evid = 0; |
444 | trackInfo.evID = evid++; | |
445 | trackInfo.zV = fESDEvent->GetPrimaryVertex()->GetZ(); | |
fa10255f | 446 | // |
447 | int ntr = fESDEvent->GetNumberOfTracks(); | |
448 | int ntrMC = fStack->GetNtrack(); | |
449 | // | |
450 | double fldBz = fESDEvent->GetMagneticField(); | |
451 | float dcaRZ[2]; | |
452 | // | |
453 | for (int itr=0;itr<ntr;itr++) { | |
454 | AliESDtrack* trc = fESDEvent->GetTrack(itr); | |
455 | // | |
456 | // at the moment we consider only TPC/ITS tracks | |
60b14b59 | 457 | //if (!trc->IsOn(AliESDtrack::kTPCin)) continue; |
458 | if (!fTPCCut->IsSelected(trc)) continue; | |
7bc38dea | 459 | // |
460 | trc->GetDZ(fVtxMC[0],fVtxMC[1],fVtxMC[2], fldBz, dcaRZ); | |
461 | // | |
fa10255f | 462 | Int_t labMC = trc->GetLabel(); |
463 | Int_t labMCabs = TMath::Abs(labMC); | |
464 | Int_t labMCTPC = trc->GetTPCLabel(); | |
465 | Int_t labMCITS = trc->GetITSLabel(); | |
466 | Int_t nClTPC = trc->GetNcls(1); | |
467 | Int_t nClITS = trc->GetNcls(0); | |
468 | // | |
469 | UInt_t& mcStatus = (UInt_t &)fMCStatus[labMCabs]; | |
097e4cdf | 470 | double pt = trc->Pt(); |
471 | double eta = trc->Eta(); | |
472 | // | |
7bc38dea | 473 | TH2* hlr = (TH2*)( (mcStatus&BIT(kMCPrimBit)) ? GetHisto(&fHistosCent,kHMCLrPresPrim) : GetHisto(&fHistosCent,kHMCLrPresSec)); |
474 | if (hlr) { | |
475 | for (int il=0;il<7;il++) { | |
476 | if (mcStatus&(0x1<<(il+kITSHitBits))) hlr->Fill(pt,il); | |
477 | } | |
478 | hlr->Fill(pt,-99); // count tracks | |
479 | } | |
480 | // | |
481 | ||
fa10255f | 482 | Bool_t reject = mcStatus & BIT(kTrCondFail); |
fa10255f | 483 | // |
484 | int mcLabType = GetMCLabType(labMCTPC,labMCITS,nClTPC,nClITS); | |
485 | // | |
60b14b59 | 486 | // if (eta>fEtaMax || eta<fEtaMin) continue; |
487 | // if (nClTPC<fMinTPCclusters) continue; | |
7bc38dea | 488 | // printf("#%3d pt:%5.2f eta:%+5.2f | Lb:%+6d LbTPC:%+6d LbITS:%+6d MCTp:%d | Ntpc:%3d Nits:%3d | Rej:%d\n", |
489 | // itr,pt,eta,labMC,labMCTPC,labMCITS,mcLabType, nClTPC,nClITS,reject); | |
490 | AliMCParticle *part = 0; | |
491 | if (labMCabs>=ntrMC) continue; | |
492 | part = (AliMCParticle*)fMCEvent->GetTrack(labMCabs); | |
493 | double ptMC = part->Pt(); | |
494 | double etaMC = part->Eta(); | |
fa10255f | 495 | // |
7bc38dea | 496 | if (fTree) { |
497 | trackInfo.lbl = labMC; | |
b824a5f5 | 498 | trackInfo.lblTPC = labMCTPC; |
7bc38dea | 499 | trackInfo.prim = (mcStatus&BIT(kMCPrimBit)) != 0; |
500 | trackInfo.rcbl = !reject; | |
501 | trackInfo.nClITS = trc->GetNcls(0); | |
502 | trackInfo.nClTPC = trc->GetNcls(1); | |
503 | trackInfo.ptMC = ptMC; | |
504 | trackInfo.etaMC = etaMC; | |
be806700 | 505 | trackInfo.charge = trc->Charge(); |
7bc38dea | 506 | trackInfo.pt = trc->Pt(); |
be806700 | 507 | trackInfo.eta = trc->Eta(); |
b824a5f5 | 508 | // |
509 | ||
510 | // | |
be806700 | 511 | trackInfo.spl = trc->GetITSModuleIndex(10); |
7bc38dea | 512 | trackInfo.phi = part->Phi(); |
513 | trackInfo.pdg = part->PdgCode(); | |
514 | trackInfo.dcaR = dcaRZ[0]; | |
515 | trackInfo.dcaZ = dcaRZ[1]; | |
516 | trackInfo.dcaRE = TMath::Sqrt(trc->GetSigmaY2()); | |
517 | trackInfo.dcaZE = TMath::Sqrt(trc->GetSigmaZ2()); | |
b824a5f5 | 518 | trackInfo.ptE = TMath::Sqrt(trc->GetSigma1Pt2()); |
519 | ||
520 | const AliExternalTrackParam* ttpc = trc->GetTPCInnerParam(); | |
521 | if (ttpc) { | |
522 | trackInfo.ptTPC = ttpc->Pt(); | |
523 | trackInfo.etaTPC = ttpc->Eta(); | |
524 | trackInfo.ptTPCE = TMath::Sqrt(ttpc->GetSigma1Pt2());//*ttpc->Pt()*ttpc->Pt(); | |
525 | } | |
526 | else { | |
527 | trackInfo.ptTPC = -1; | |
528 | trackInfo.etaTPC = 0; | |
529 | trackInfo.ptTPCE = -1; | |
530 | } | |
7bc38dea | 531 | // |
be806700 | 532 | trackInfo.nClITSMC = 0; |
7bc38dea | 533 | for (int il=0;il<7;il++) { |
be806700 | 534 | trackInfo.mcCl[il] = (mcStatus & (0x1<<(il+kITSHitBits))) != 0; |
535 | if (trackInfo.mcCl[il]) trackInfo.nClITSMC++; | |
b824a5f5 | 536 | trackInfo.rcCl[il] = trc->HasPointOnITSLayer(il); |
7bc38dea | 537 | trackInfo.fcCl[il] = trc->HasSharedPointOnITSLayer(il); |
538 | } | |
539 | fTree->Fill(); | |
540 | } | |
097e4cdf | 541 | if (reject) { |
542 | GetHisto(&fHistosCent,(mcStatus & BIT(kMCPrimBit)) ? kHMatchStatusNRcblPrim:kHMatchStatusNRcblSec)->Fill(pt,mcLabType); // matching status | |
543 | if (dump) { | |
544 | printf("--#%4d Pt:%5.2f Eta:%5.2f |",itr,pt,eta); | |
545 | for (int k=0;k<32;k++) printf("%d", (mcStatus&(0x1<<k)) ? 1:0); | |
546 | printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType); | |
547 | } | |
548 | continue; | |
549 | } | |
550 | GetHisto(&fHistosCent,kHMatchStatusRcbl)->Fill(pt,mcLabType); // matching status | |
fa10255f | 551 | // |
097e4cdf | 552 | TObjArray* clList = 0; |
fa10255f | 553 | // |
7bc38dea | 554 | if (part) { |
555 | /* | |
556 | int pdg = part->PdgCode(); | |
557 | if (TMath::Abs(pdg)!=211) continue; | |
558 | if (nClITS<7) continue; | |
559 | */ | |
097e4cdf | 560 | #ifdef _CLUS_LIST_ |
097e4cdf | 561 | // |
562 | clList = (TObjArray*)fClLists.At(labMCabs); | |
563 | if (clList) { | |
564 | printf("Clusters for track %d MCLabel:%d\n",itr,labMC); | |
565 | for (int ic=0;ic<clList->GetEntriesFast();ic++) { | |
566 | AliITSUClusterPix* cl = (AliITSUClusterPix*) clList->At(ic); | |
567 | cl->Print(); | |
568 | } | |
569 | } | |
097e4cdf | 570 | #endif |
7bc38dea | 571 | } |
097e4cdf | 572 | |
fa10255f | 573 | if (part) { |
574 | // | |
fa10255f | 575 | // |
097e4cdf | 576 | if (dump) { |
577 | printf(" #%4d Pt:%5.2f Eta:%5.2f |",itr,ptMC,etaMC); | |
60b14b59 | 578 | for (int k=0;k<32;k++) printf("%d", (mcStatus&(0x1<<k)) ? 1:0); |
579 | printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType); | |
097e4cdf | 580 | } |
581 | // | |
582 | if ( (mcStatus&BIT(kMCPrimBit)) ) { | |
be806700 | 583 | if (dump1 && TMath::Abs(ptMC-0.5)<0.1) { |
584 | printf("#%4d Pt:%5.2f Eta:%5.2f | ",itr,ptMC,etaMC); | |
585 | for (int k=0;k<7;k++) printf("(%d/%d)", (mcStatus&(0x1<<k)) ? 1:0, trc->HasPointOnITSLayer(k)); | |
7bc38dea | 586 | printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType); |
587 | } | |
fa10255f | 588 | // compare MC vs reco track params |
589 | if (ptMC>0) GetHisto(&fHistosCentMCLb,kHResPTvsPTMC, mcLabType)->Fill(ptMC, (ptMC-pt)/ptMC); | |
590 | // | |
591 | // DCA resolution | |
fa10255f | 592 | GetHisto(&fHistosCentMCLb,kHResDCARvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[0]); |
593 | GetHisto(&fHistosCentMCLb,kHResDCAZvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[1]); | |
594 | // | |
595 | } | |
596 | // | |
597 | } | |
598 | // | |
599 | ||
600 | ||
601 | } | |
602 | // | |
603 | ||
604 | } | |
605 | ||
606 | //_________________________________________________________________________ | |
607 | void AliTaskITSUPerf::BuildMCInfo() | |
608 | { | |
609 | // build mc truth info for tracks | |
610 | // | |
611 | int ntrMC = fStack->GetNtrack(); | |
612 | // | |
613 | if (fMCStatus.GetSize()<ntrMC) fMCStatus.Set(ntrMC); | |
614 | fMCStatus.Reset(); | |
097e4cdf | 615 | #ifdef _CLUS_LIST_ |
616 | fClLists.Delete(); | |
617 | if (fClLists.GetSize()<ntrMC) fClLists.Expand(ntrMC); | |
618 | // | |
619 | // create empty lists for reconstructed tracks | |
620 | int ntr = fESDEvent->GetNumberOfTracks(); | |
621 | for (int itr=0;itr<ntr;itr++) { | |
622 | AliESDtrack* trc = fESDEvent->GetTrack(itr); | |
623 | int lb = TMath::Abs(trc->GetLabel()); | |
624 | if (lb<ntrMC && !fClLists.At(itr)) { | |
625 | fClLists.AddAt(new TObjArray(7),lb); | |
626 | } | |
627 | } | |
628 | #endif | |
fa10255f | 629 | // |
630 | for (int itr=0;itr<ntrMC;itr++) { | |
631 | // | |
632 | AliMCParticle *part = (AliMCParticle*)fMCEvent->GetTrack(itr); | |
633 | if (!part->Charge()) continue; | |
634 | // | |
635 | UInt_t& mcStatus = (UInt_t &)fMCStatus[itr]; | |
636 | mcStatus = 0; | |
637 | // | |
638 | if (fStack->IsPhysicalPrimary(itr)) mcStatus |= BIT(kMCPrimBit); | |
639 | // | |
640 | /* | |
641 | TParticle* tp = part->Particle(); | |
642 | printf("mc%4d IsPri=%d",itr,fStack->IsPhysicalPrimary(itr)); tp->Print(); | |
643 | */ | |
644 | // | |
645 | } | |
646 | // | |
647 | int nTrCond = fTrackingCond ? fTrackingCond->GetEntriesFast() : 0; | |
648 | // fill MC clusters info (if available) | |
649 | if (fRPTree) { | |
650 | // load clusters | |
651 | fITS->LoadClusters(fRPTree); | |
097e4cdf | 652 | fITS->ProcessClusters(); // order in the same way as in reconstruction |
fa10255f | 653 | // |
097e4cdf | 654 | for (int ilr=0;ilr<fITS->GetNLayersActive();ilr++) { |
fa10255f | 655 | AliITSURecoLayer* lr = fITS->GetLayerActive(ilr); |
656 | for (int icl=lr->GetNClusters();icl--;) { | |
657 | AliCluster* cl = lr->GetCluster(icl); | |
658 | // check labels | |
659 | for (int ilb=0;ilb<3;ilb++) { | |
660 | Int_t lab = cl->GetLabel(ilb); | |
661 | if (lab<0||lab>=ntrMC) continue; | |
662 | UInt_t& mcStatus = (UInt_t &)fMCStatus[lab]; | |
663 | mcStatus |= 0x1<<(ilr+kITSHitBits); | |
097e4cdf | 664 | #ifdef _CLUS_LIST_ |
665 | TObjArray *clList = (TObjArray*)fClLists.At(lab); | |
666 | if (clList) clList->AddLast(cl); // create cl.list only for reconstructed tracks | |
667 | #endif | |
fa10255f | 668 | } |
669 | } | |
670 | } | |
671 | // | |
672 | if (nTrCond) { | |
673 | for (int itr=0;itr<ntrMC;itr++) { | |
674 | // | |
675 | UInt_t& mcStatus = (UInt_t &)fMCStatus[itr]; | |
676 | Bool_t fail = kTRUE; | |
677 | UShort_t patt = (mcStatus&0xffff); | |
678 | for (int ic=0;ic<nTrCond;ic++) { | |
679 | AliITSUTrackCond* cnd = (AliITSUTrackCond*) fTrackingCond->At(ic); | |
680 | if (cnd->CheckPattern(patt)) {fail=kFALSE; break;} | |
681 | } | |
682 | if (fail) mcStatus |= BIT(kTrCondFail); | |
7bc38dea | 683 | // |
684 | /* | |
685 | int nlr = AliITSUAux::NumberOfBitsSet(patt); | |
686 | AliMCParticle *part = (AliMCParticle*)fMCEvent->GetTrack(itr); | |
687 | double etaMC = part->Eta(); | |
688 | if (etaMC>=fEtaMin && etaMC<=fEtaMax && nlr) { | |
689 | TH2* hlr = (TH2*)( (mcStatus&BIT(kMCPrimBit)) ? GetHisto(&fHistosCent,kHMCLrPresPrim) : GetHisto(&fHistosCent,kHMCLrPresSec)); | |
690 | if (hlr) { | |
691 | double ptMC = part->Pt(); | |
692 | for (int il=0;il<7;il++) { | |
693 | if (mcStatus&(0x1<<(il+kITSHitBits))) hlr->Fill(ptMC,il); | |
694 | } | |
695 | hlr->Fill(ptMC,-99); // count tracks | |
696 | } | |
697 | } | |
698 | */ | |
fa10255f | 699 | } |
700 | } | |
701 | // | |
702 | } | |
703 | // | |
704 | } | |
705 | ||
706 | ||
707 | //_________________________________________________________________________ | |
708 | Int_t AliTaskITSUPerf::GetMCLabType(Int_t labMCTPC,Int_t labMCITS, Int_t nClTPC, Int_t nClITS) | |
709 | { | |
710 | // build type of track (0: both ITS and TPC correct, 1: ITS corr. TPC fake, 2: ITS fake TPC corr., 3: ITS fake TPC fake, 4: mismatch | |
711 | // | |
712 | int tp = 0; | |
713 | if (nClTPC>0 && nClITS<1) return kITSTPCNoMatch; | |
714 | if (TMath::Abs(labMCTPC)!=TMath::Abs(labMCITS)) tp = kITSTPCMismatch; | |
715 | else { | |
716 | if (labMCITS<0) tp = labMCTPC<0 ? kITS0TPC0 : kITS0TPC1; | |
717 | else tp = labMCTPC<0 ? kITS1TPC0 : kITS1TPC1; | |
718 | } | |
719 | return tp; | |
720 | // | |
721 | } | |
722 | ||
723 | //_________________________________________________________________________ | |
724 | Int_t AliTaskITSUPerf::GetCentralityBin() const | |
725 | { | |
726 | // calculate centrality bin | |
727 | return 0; | |
728 | } | |
729 |