]>
Commit | Line | Data |
---|---|---|
1ee39b3a | 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: AliTRDefficiency.cxx 27496 2008-07-22 08:35:45Z cblume $ */ | |
17 | ||
18 | //////////////////////////////////////////////////////////////////////////// | |
19 | // // | |
20 | // Reconstruction QA // | |
21 | // // | |
22 | // Authors: // | |
23 | // Markus Fasel <M.Fasel@gsi.de> // | |
24 | // // | |
25 | //////////////////////////////////////////////////////////////////////////// | |
26 | ||
eb05d549 | 27 | #include <TROOT.h> |
668a0654 | 28 | #include <TFile.h> |
eb05d549 | 29 | #include <TStyle.h> |
1ee39b3a | 30 | #include <TClonesArray.h> |
31 | #include <TObjArray.h> | |
32 | #include <TProfile.h> | |
e2e3cec2 | 33 | #include <TPad.h> |
eb05d549 | 34 | #include <TCanvas.h> |
a310e49b | 35 | #include <TLegend.h> |
eb05d549 | 36 | #include <THnSparse.h> |
37 | #include <TH2.h> | |
38 | #include <TH3.h> | |
39 | #include <THStack.h> | |
1ee39b3a | 40 | #include "TTreeStream.h" |
41 | ||
42 | #include "AliPID.h" | |
43 | #include "AliESDtrack.h" | |
44 | #include "AliTrackReference.h" | |
45 | #include "AliExternalTrackParam.h" | |
46 | #include "AliTracker.h" | |
47 | #include "AliAnalysisManager.h" | |
48 | ||
eb05d549 | 49 | #include "AliTRDgeometry.h" |
50 | #include "AliTRDtrackV1.h" | |
1ee39b3a | 51 | #include "Cal/AliTRDCalPID.h" |
52 | #include "AliTRDefficiency.h" | |
53 | #include "info/AliTRDtrackInfo.h" | |
54 | ||
55 | ClassImp(AliTRDefficiency) | |
56 | ||
57 | //____________________________________________________________________ | |
58 | AliTRDefficiency::AliTRDefficiency() | |
f2e89a4c | 59 | :AliTRDrecoTask() |
db99a57a | 60 | ,fMissed(NULL) |
eb05d549 | 61 | ,fProj(NULL) |
f8f46e4d | 62 | { |
63 | // | |
64 | // Default constructor | |
65 | // | |
f2e89a4c | 66 | SetNameTitle("TRDefficiency", "TRD barrel tracking efficiency checker"); |
f8f46e4d | 67 | } |
68 | ||
f2e89a4c | 69 | //____________________________________________________________________ |
f8f46e4d | 70 | AliTRDefficiency::AliTRDefficiency(char* name) |
71 | :AliTRDrecoTask(name, "TRD barrel tracking efficiency checker") | |
db99a57a | 72 | ,fMissed(NULL) |
eb05d549 | 73 | ,fProj(NULL) |
1ee39b3a | 74 | { |
75 | // | |
76 | // Default constructor | |
77 | // | |
78 | } | |
79 | ||
80 | //____________________________________________________________________ | |
81 | AliTRDefficiency::~AliTRDefficiency() | |
82 | { | |
83 | // Destructor | |
84 | if(fMissed){ | |
85 | fMissed->Delete(); | |
86 | delete fMissed; | |
87 | } | |
88 | } | |
89 | ||
eb05d549 | 90 | // //____________________________________________________________________ |
91 | // void AliTRDefficiency::UserCreateOutputObjects() | |
92 | // { | |
93 | // // | |
94 | // // Create output objects | |
95 | // // | |
96 | // | |
97 | // const Int_t nbins = AliTRDCalPID::kNMom; | |
98 | // Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.}; | |
99 | // | |
100 | // TH1 *h = NULL; | |
101 | // fContainer = new TObjArray(); fContainer->SetOwner(); | |
102 | // for(Int_t is=0; is<AliPID::kSPECIES; is++){ | |
103 | // fContainer->Add(h = new TProfile(Form("h%s", AliTRDCalPID::GetPartSymb(is)), AliPID::ParticleShortName(is), nbins, xbins)); | |
104 | // h->SetLineColor(AliTRDCalPID::GetPartColor(is)); | |
105 | // h->SetMarkerColor(AliTRDCalPID::GetPartColor(is)); | |
106 | // h->SetMarkerStyle(24); | |
107 | // } | |
108 | // fContainer->Add(h = new TProfile("h", "", nbins, xbins)); | |
109 | // h->SetMarkerStyle(7); | |
110 | // PostData(1, fContainer); | |
111 | // } | |
112 | ||
1ee39b3a | 113 | //____________________________________________________________________ |
eb05d549 | 114 | TH1* AliTRDefficiency::PlotBasicEff(const AliTRDtrackV1 *track) |
1ee39b3a | 115 | { |
eb05d549 | 116 | // plot TRD efficiency based on ESD info |
1ee39b3a | 117 | |
eb05d549 | 118 | if(!fkESD){ |
119 | AliDebug(4, "No ESD info."); | |
120 | return NULL; | |
121 | } | |
1ee39b3a | 122 | |
eb05d549 | 123 | THnSparse *H(NULL); |
124 | if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hEFF")))){ | |
125 | AliWarning("No output container defined."); | |
126 | return NULL; | |
1ee39b3a | 127 | } |
eb05d549 | 128 | if(track) fkTrack = track; |
129 | ||
668a0654 | 130 | Double_t val[7]; memset(val, 0, 7*sizeof(Double_t)); |
eb05d549 | 131 | ULong_t status(fkESD->GetStatus()); |
132 | val[0] =((status&AliESDtrack::kTRDin)?1:0) + | |
133 | ((status&AliESDtrack::kTRDStop)?1:0) + | |
134 | ((status&AliESDtrack::kTRDout)?2:0); | |
ef09211d | 135 | val[1] = fkESD->Phi(); |
136 | val[2] = fkESD->Eta(); | |
137 | val[3] = GetPtBin(fPt>0.?fPt:fkESD->Pt()); | |
eb05d549 | 138 | val[4] = 0.; |
139 | if(fkMC){ | |
140 | if(fkMC->GetLabel() == fkMC->GetTRDlabel()) val[4] = 0.; | |
141 | else if(fkMC->GetLabel() == -fkMC->GetTRDlabel()) val[4] = 1.; | |
142 | else val[4] = -1.; | |
143 | } | |
0f2c4c4f | 144 | val[5] = fkTrack?fkTrack->GetNumberOfTracklets():0; |
668a0654 | 145 | // down scale PID resolution |
0f2c4c4f | 146 | Int_t spc(fSpecies); if(spc==-6) spc=0; if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3; |
2b33fd9c | 147 | val[6] = spc; |
0f2c4c4f | 148 | |
149 | if(fkTrack) AliDebug(2, Form("%3d[%s] tracklets[%d] label[%2d %+2d] species[%d %d] in[%c] out[%c] stop[%c]", | |
150 | fkESD->GetId(), fkTrack->IsTPCseeded()?"TPC":"ITS", Int_t(val[5]), fkMC->GetLabel(), fkMC->GetTRDlabel(), fSpecies, spc, | |
151 | (status&AliESDtrack::kTRDin)?'x':'o', (status&AliESDtrack::kTRDout)?'x':'o', (status&AliESDtrack::kTRDStop)?'x':'o')); | |
eb05d549 | 152 | H->Fill(val); |
153 | return NULL; | |
154 | } | |
155 | ||
156 | // //____________________________________________________________________ | |
157 | // TH1* AliTRDefficiency::PlotMC(const AliTRDtrackV1 *track) | |
158 | // { | |
159 | // // plot TRD efficiency based on MC info | |
160 | // | |
161 | // if(!HasMC()) return NULL; | |
162 | // if(!fkESD){ | |
163 | // AliDebug(4, "No ESD info."); | |
164 | // return NULL; | |
165 | // } | |
166 | // if(!fkMC){ | |
167 | // AliDebug(4, "No MC info."); | |
168 | // return NULL; | |
169 | // } | |
170 | // | |
171 | // THnSparse *H(NULL); | |
172 | // if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hMC")))){ | |
173 | // AliWarning("No output container defined."); | |
174 | // return NULL; | |
175 | // } | |
176 | // if(track) fkTrack = track; | |
177 | // Double_t val[11]; memset(val, 0, 11*sizeof(Double_t)); | |
178 | // ULong_t status(fkESD->GetStatus()); | |
179 | // val[0] =((status&AliESDtrack::kTRDin)?1:0) + | |
180 | // ((status&AliESDtrack::kTRDStop)?1:0) + | |
181 | // ((status&AliESDtrack::kTRDout)?2:0); | |
182 | // val[1] = fkESD->Phi(); | |
183 | // val[2] = fkESD->Eta(); | |
184 | // val[3] = DebugLevel()>=1?GetPtBin(fkESD->Pt()):GetPtBinSignificant(fkESD->Pt()); | |
185 | // if(fkTrack){ // read track status in debug mode with friends | |
186 | // val[4] = fkTrack->GetStatusTRD(-1); | |
187 | // for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++) val[5+ily]=fkTrack->GetStatusTRD(ily); | |
188 | // } | |
189 | // H->Fill(val); | |
190 | // | |
191 | // } | |
1ee39b3a | 192 | |
193 | //____________________________________________________________________ | |
eb05d549 | 194 | void AliTRDefficiency::LocalUserExec(Option_t *) |
1ee39b3a | 195 | { |
196 | // | |
eb05d549 | 197 | // Do it obsolete |
1ee39b3a | 198 | // |
199 | ||
eb05d549 | 200 | Int_t labelsacc[10000]; |
1ee39b3a | 201 | memset(labelsacc, 0, sizeof(Int_t) * 10000); |
202 | ||
5935a6da | 203 | fTracks = dynamic_cast<TObjArray *>(GetInputData(1)); |
204 | if(!fTracks) return; | |
205 | if(!fTracks->GetEntriesFast()) return; | |
206 | else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName())); | |
1ee39b3a | 207 | if(!fMissed){ |
208 | fMissed = new TClonesArray("AliTRDtrackInfo", 10); | |
209 | fMissed->SetOwner(); | |
210 | } | |
211 | ||
212 | Float_t mom; | |
213 | Int_t selection[10000], nselect = 0; | |
214 | ULong_t status; Int_t pidx; | |
215 | Int_t nTRD = 0, nTPC = 0, nMiss = 0; | |
db99a57a | 216 | AliTRDtrackInfo *track = NULL; |
217 | AliTrackReference *ref = NULL; | |
218 | AliExternalTrackParam *esd = NULL; | |
1ee39b3a | 219 | for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){ |
220 | track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk); | |
221 | ||
222 | if(!track->HasESDtrack()) continue; | |
223 | status = track->GetStatus(); | |
224 | ||
225 | // missing TPC propagation - interesting for SA | |
226 | if(!(status&AliESDtrack::kTPCout)) continue; | |
227 | ||
228 | // missing MC info. | |
229 | if(HasMCdata() && track->GetNTrackRefs() <= 1) continue; | |
230 | ||
231 | nTPC++; | |
232 | selection[nselect++]=itrk; | |
233 | ref = track->GetTrackRef(0); | |
234 | esd = track->GetESDinfo()->GetOuterParam(); | |
235 | mom = ref ? ref->P(): esd->P(); | |
236 | pidx = AliTRDCalPID::GetPartIndex(track->GetPDG()); | |
237 | pidx = TMath::Max(pidx, 0); | |
8f7b1226 | 238 | AliDebug(4, Form("PID: %d", pidx)); |
1ee39b3a | 239 | |
240 | //Int_t n = track->GetNumberOfClusters(); | |
241 | // where are this tracklets ??? | |
242 | //if(ncls0 > ncls1) printf("%3d ESD[%3d] TRD[%3d|%3d]\n", itrk, ncls0, ncls1, n); | |
243 | if(track->GetNumberOfClustersRefit()){ | |
244 | ((TProfile*)fContainer->At(pidx))->Fill(mom, 1.); | |
245 | labelsacc[nTRD] = track->GetLabel(); | |
246 | nTRD++; | |
247 | continue; | |
248 | } | |
249 | ||
250 | ||
251 | ||
252 | Float_t xmed, xleng; | |
253 | Int_t iref = 1; Bool_t found = kFALSE; | |
254 | while((ref = track->GetTrackRef(iref))){ | |
255 | xmed = .5*(ref->LocalX() + track->GetTrackRef(iref-1)->LocalX()); | |
256 | xleng= (ref->LocalX() - track->GetTrackRef(iref-1)->LocalX()); | |
257 | if(TMath::Abs(xmed - 298.5) < .5 && | |
258 | TMath::Abs(xleng - 3.7) < .1){ | |
259 | found = kTRUE; | |
260 | break; | |
261 | } | |
262 | iref++; | |
263 | } | |
264 | if(!found){ | |
265 | nTPC--; | |
266 | // track missing first layer. Maybe interesting for SA. | |
267 | continue; | |
268 | } | |
269 | nselect--; | |
270 | new ((*fMissed)[nMiss]) AliTRDtrackInfo(*track); | |
271 | nMiss++; | |
272 | } | |
273 | AliDebug(2, Form("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast())); | |
274 | ||
275 | ||
276 | // Find double tracks | |
277 | Float_t threshold = 10.; | |
db99a57a | 278 | AliTrackReference *refMiss = NULL; |
279 | AliExternalTrackParam *op = NULL; | |
280 | AliTRDtrackInfo *tt = NULL; | |
1ee39b3a | 281 | for(Int_t imiss=0; imiss<nMiss; imiss++){ |
282 | //printf("Searching missing %d ...\n", imiss); | |
283 | ||
284 | // get outer param of missed | |
285 | tt = (AliTRDtrackInfo*)fMissed->UncheckedAt(imiss); | |
286 | op = tt->GetESDinfo()->GetOuterParam(); | |
287 | Double_t alpha = op->GetAlpha(), cosa = TMath::Cos(alpha), sina = TMath::Sin(alpha); | |
288 | ||
289 | Double_t xyz[3], x0, y0, z0, x, y, z, dx, dy, dz, d; | |
290 | ||
291 | Bool_t bFOUND = kFALSE; | |
292 | for(Int_t iselect=0; iselect<nselect; iselect++){ | |
293 | track = (AliTRDtrackInfo*)fTracks->UncheckedAt(selection[iselect]); | |
294 | ||
295 | // check first MC ... if available | |
296 | d = 0; | |
297 | for(Int_t iref=0; iref<track->GetNTrackRefs(); iref++){ | |
298 | if(!(ref = track->GetTrackRef(iref))) continue; | |
299 | if((refMiss = tt->GetTrackRef(iref))){ | |
300 | dy = ref->LocalY() - refMiss->LocalY(); | |
301 | dz = ref->Z() - refMiss->Z(); | |
302 | } else { | |
303 | // compare missOP with refTrackRef in LTC | |
304 | x0 = ref->LocalX(); | |
305 | op->GetYAt(x0, AliTracker::GetBz(), y0); | |
306 | op->GetZAt(x0, AliTracker::GetBz(), z0); | |
307 | dy = y0 - ref->LocalY(); | |
308 | dz = z0 - ref->Z(); | |
309 | } | |
310 | d += (dy*dy + dz*dz); | |
311 | } | |
312 | //printf("\td[%d] = %f N[%d]\n", selection[iselect], d, track->GetNTrackRefs()); | |
313 | if((track->GetNTrackRefs())){ | |
314 | d /= track->GetNTrackRefs(); | |
315 | if(d < threshold){ | |
316 | //printf("\t\tFound %2d in ref[%3d] : d[%f]\n", imiss, selection[iselect], d/track->GetNTrackRefs()); | |
317 | bFOUND = kTRUE; break; | |
318 | } | |
319 | } | |
320 | ||
321 | // process outer param ... always available | |
322 | // compare missOP with OP in GTC | |
323 | esd = track->GetESDinfo()->GetOuterParam(); | |
324 | esd->GetXYZ(xyz); | |
325 | x0 = esd->GetX(); | |
326 | op->GetYAt(x0, AliTracker::GetBz(), y0); | |
327 | op->GetZAt(x0, AliTracker::GetBz(), z0); | |
328 | x = x0*cosa - y0*sina; | |
329 | y = x0*sina + y0*cosa; | |
330 | z = z0; | |
331 | dx=xyz[0]-x; | |
332 | dy=xyz[1]-y; | |
333 | dz=xyz[2]-z; | |
334 | d = dx*dx+dy*dy+dz*dz; | |
335 | //printf("\td[%d] = %f op\n", selection[iselect], d); | |
336 | if(d < threshold){ | |
337 | //printf("\t\tFound %2d in op[%3d] : d[%f] dx[%5.2f] dy[%5.2f] dz[%5.2f]\n", imiss, selection[iselect], d, dx, dy, dz); | |
338 | bFOUND = kTRUE; break; | |
339 | } | |
340 | } | |
341 | if(bFOUND) nTPC--; | |
342 | else{ | |
343 | ref = tt->GetTrackRef(0); | |
344 | mom = ref ? ref->P(): op->P(); | |
345 | pidx = AliTRDCalPID::GetPartIndex(tt->GetPDG()); | |
346 | pidx = TMath::Max(pidx, 0); | |
347 | ((TProfile*)fContainer->At(pidx))->Fill(mom, 0.); | |
348 | AliDebug(2, Form(" NOT bFOUND Id[%d] Mom[%f]\n", tt->GetTrackId(), mom)); | |
349 | } | |
350 | } | |
351 | ||
352 | AliDebug(2, Form("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast())); | |
353 | ||
354 | //fMissed->Delete(); | |
355 | // check for double countings | |
356 | Int_t indices[10000]; memset(indices, 0, sizeof(Int_t) * 10000); | |
357 | TMath::Sort(nTRD, labelsacc, indices); | |
358 | if(DebugLevel() > 2){ | |
359 | for(Int_t itk = 0; itk < nTRD - 1; itk++) | |
360 | if(labelsacc[indices[itk]] ==labelsacc[indices[itk + 1]]) printf("Double counted MC track: %d\n", labelsacc[indices[itk]]); | |
361 | } | |
1ee39b3a | 362 | } |
363 | ||
1ee39b3a | 364 | //____________________________________________________________________ |
365 | Bool_t AliTRDefficiency::GetRefFigure(Int_t ifig) | |
366 | { | |
367 | // Steer reference figures | |
368 | ||
e2e3cec2 | 369 | if(!gPad){ |
370 | AliWarning("Please provide a canvas to draw results."); | |
371 | return kFALSE; | |
372 | } | |
373 | gPad->SetLogx(); | |
374 | ||
a310e49b | 375 | TLegend *leg(NULL); |
376 | Bool_t bFIRST(kTRUE); | |
377 | TProfile *h(NULL); | |
1ee39b3a | 378 | switch(ifig){ |
379 | case 0: | |
380 | h = (TProfile*)fContainer->At(AliPID::kSPECIES); | |
381 | for(Int_t is=0; is<AliPID::kSPECIES; is++){ | |
382 | h->Add((TProfile*)fContainer->At(is)); | |
383 | } | |
a310e49b | 384 | h->SetMarkerStyle(24); |
1ee39b3a | 385 | h->SetMarkerColor(kBlack); |
386 | h->SetLineColor(kBlack); | |
e2e3cec2 | 387 | h->SetTitle("TRD Efficiency integrated"); |
a310e49b | 388 | h->SetXTitle("p [GeV/c]"); |
e2e3cec2 | 389 | h->GetXaxis()->SetMoreLogLabels(); |
a310e49b | 390 | h->SetYTitle("Efficiency"); |
391 | h->GetYaxis()->CenterTitle(); | |
1ee39b3a | 392 | h->Draw("e1"); |
393 | break; | |
394 | case 1: | |
395 | bFIRST = kTRUE; | |
396 | for(Int_t is=0; is<AliPID::kSPECIES; is++){ | |
397 | if(!(h = (TProfile*)fContainer->At(is))) continue; | |
a310e49b | 398 | h->SetMarkerStyle(24); |
1ee39b3a | 399 | if(bFIRST){ |
400 | h->Draw("e1"); | |
a310e49b | 401 | h->SetXTitle("p [GeV/c]"); |
e2e3cec2 | 402 | h->GetXaxis()->SetMoreLogLabels(); |
a310e49b | 403 | h->SetYTitle("Efficiency"); |
404 | h->GetYaxis()->CenterTitle(); | |
405 | h->GetYaxis()->SetRangeUser(0.8, 1.05); | |
406 | leg=new TLegend(.7, .2, .98, .6); | |
407 | leg->SetHeader("Species"); | |
408 | leg->SetBorderSize(0); | |
409 | leg->SetFillStyle(0); | |
410 | leg->AddEntry(h, h->GetTitle(), "pl"); | |
411 | } else { | |
412 | leg->AddEntry(h, h->GetTitle(), "pl"); | |
413 | h->Draw("same e1"); | |
414 | } | |
1ee39b3a | 415 | bFIRST = kFALSE; |
416 | } | |
7fe4e88b | 417 | if(leg) leg->Draw(); |
1ee39b3a | 418 | break; |
419 | } | |
420 | return kTRUE; | |
421 | } | |
422 | ||
eb05d549 | 423 | //________________________________________________________ |
424 | TObjArray* AliTRDefficiency::Histos() | |
425 | { | |
426 | // | |
427 | // Define histograms | |
428 | // | |
429 | ||
430 | if(fContainer) return fContainer; | |
431 | ||
432 | fContainer = new TObjArray(1); fContainer->SetOwner(kTRUE); | |
433 | THnSparse *H(NULL); | |
434 | TString st; | |
435 | ||
436 | //++++++++++++++++++++++ | |
437 | // cluster to detector | |
438 | if(!(H = (THnSparseI*)gROOT->FindObject("hEFF"))){ | |
668a0654 | 439 | const Int_t mdim(7); |
eb05d549 | 440 | Int_t nlabel(1); |
0f2c4c4f | 441 | const Char_t *eTitle[mdim] = {"status", "#phi [rad]", "eta", "p_{t} [bin]", "label", "N_{trklt}", "chg*spec*rc"}; |
ef09211d | 442 | const Int_t eNbins[mdim] = { 5, 144, 45, fNpt-1, nlabel, AliTRDgeometry::kNlayer-2, 5}; |
443 | const Double_t eMin[mdim] = { -0.5, -TMath::Pi(), -.9, -0.5, -0.5, 1.5, -2.5}, | |
444 | eMax[mdim] = { 4.5, TMath::Pi(), .9, fNpt-1.5, nlabel-0.5, 5.5, 2.5}; | |
eb05d549 | 445 | st = "basic efficiency;"; |
446 | // define minimum info to be saved in non debug mode | |
447 | Int_t ndim=DebugLevel()>=1?mdim:(HasMCdata()?5:4); | |
448 | for(Int_t idim(0); idim<ndim; idim++){ st += eTitle[idim]; st+=";";} | |
449 | H = new THnSparseI("hEFF", st.Data(), ndim, eNbins, eMin, eMax); | |
450 | /* TAxis *ax(H->GetAxis(0)); const Char_t *lTRDflag[] = {"!TRDin", "TRDin", "TRDin&TRDStop", "TRDin&TRDout", "TRDin&TRDout&TRDStop"}; | |
451 | for(Int_t ibin(1); ibin<=ax->GetNbins(); ibin++) ax->SetBinLabel(ibin, lTRDflag[ibin-1]);*/ | |
452 | } else H->Reset(); | |
453 | fContainer->AddAt(H, 0); | |
454 | ||
455 | return fContainer; | |
456 | } | |
1ee39b3a | 457 | |
458 | //____________________________________________________________________ | |
459 | Bool_t AliTRDefficiency::PostProcess() | |
460 | { | |
eb05d549 | 461 | // Fit, Project, Combine, Extract values from the containers filled during execution |
462 | ||
463 | if (!fContainer) { | |
464 | AliError("ERROR: list not available"); | |
465 | return kFALSE; | |
466 | } | |
467 | if(!fProj){ | |
468 | AliInfo("Building array of projections ..."); | |
0f2c4c4f | 469 | fProj = new TObjArray(200); fProj->SetOwner(kTRUE); |
eb05d549 | 470 | } |
0f2c4c4f | 471 | // set pt/p segmentation. guess from data |
472 | THnSparse *H(NULL); | |
473 | if(!(H = (THnSparse*)fContainer->FindObject("hEFF"))){ | |
474 | AliError("Missing/Wrong data @ hEFF."); | |
475 | return kFALSE; | |
476 | } | |
477 | fNpt=H->GetAxis(3)->GetNbins()+1; | |
478 | if(!MakeMomSegmentation()) return kFALSE; | |
479 | ||
eb05d549 | 480 | if(!MakeProjectionBasicEff()) return kFALSE; |
481 | return kTRUE; | |
482 | } | |
483 | ||
484 | //____________________________________________________________________ | |
485 | Bool_t AliTRDefficiency::MakeProjectionBasicEff() | |
486 | { | |
487 | // Make basic efficiency plots | |
488 | ||
489 | if(!fContainer || !fProj){ | |
490 | AliError("Missing data container."); | |
491 | return kFALSE; | |
492 | } | |
493 | THnSparse *H(NULL); | |
494 | if(!(H = (THnSparse*)fContainer->FindObject("hEFF"))){ | |
495 | AliError("Missing/Wrong data @ hEFF."); | |
496 | return kFALSE; | |
497 | } | |
498 | Int_t ndim(H->GetNdimensions()); //Bool_t debug(ndim>Int_t(kNdimCl)); | |
499 | TAxis *aa[11], *al(NULL); memset(aa, 0, sizeof(TAxis*) * 11); | |
500 | for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id); | |
501 | if(H->GetNdimensions() > 4) al = H->GetAxis(4); | |
502 | Int_t nlab=al?3:1; | |
503 | ||
504 | // define rebinning strategy | |
505 | //const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5}; | |
506 | AliTRDrecoProjection hp[15]; TObjArray php(15); | |
507 | const Char_t *stat[] = {"!TRDin", "TRDin", "TRDin&TRDStop", "TRDin&TRDout", "TRDin&TRDout&TRDStop"}; | |
0f2c4c4f | 508 | const Char_t *lab[] = {"MC-Bad", "MC-Good", "MC-Accept"}; |
eb05d549 | 509 | Int_t ih(0); |
510 | for(Int_t ilab(0); ilab<nlab; ilab++){ | |
511 | for(Int_t istat(0); istat<5; istat++){ | |
eb05d549 | 512 | hp[ih].Build(Form("HEff%d%d", ilab, istat), |
0f2c4c4f | 513 | Form("Efficiency :: Stat[#bf{%s}] %s", stat[istat], nlab>1?lab[ilab]:""), 2, 1, 3, aa); |
eb05d549 | 514 | //hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
515 | php.AddLast(&hp[ih++]); //np[isel]++; | |
516 | } | |
517 | } | |
518 | AliInfo(Form("Build %3d 3D projections.", ih)); | |
519 | ||
668a0654 | 520 | AliTRDrecoProjection *pr0(NULL), *pr1(NULL); |
eb05d549 | 521 | Int_t istatus, ilab(0), coord[11]; memset(coord, 0, sizeof(Int_t) * 11); Double_t v = 0.; |
522 | for (Long64_t ib(0); ib < H->GetNbins(); ib++) { | |
523 | v = H->GetBinContent(ib, coord); if(v<1.) continue; | |
524 | istatus = coord[0]-1; | |
0f2c4c4f | 525 | ilab=al?0:coord[4]; |
eb05d549 | 526 | Int_t isel = ilab*5+istatus; |
668a0654 | 527 | if(isel>=ih){ |
528 | AliError(Form("Wrong selection %d [%3d]", isel, ih)); | |
529 | return kFALSE; | |
530 | } | |
531 | if(!(pr0=(AliTRDrecoProjection*)php.At(isel))) { | |
532 | AliError(Form("Missing projection %d", isel)); | |
533 | return kFALSE; | |
534 | } | |
535 | if(strcmp(pr0->H()->GetName(), Form("HEff%d%d", ilab, istatus))!=0){ | |
536 | AliError(Form("Projection mismatch :: request[HEff%d%d] found[%s]", ilab, istatus, pr0->H()->GetName())); | |
537 | return kFALSE; | |
538 | } | |
eb05d549 | 539 | for(Int_t jh(0); jh<1/*np[isel]*/; jh++) ((AliTRDrecoProjection*)php.At(isel+jh))->Increment(coord, v); |
540 | } | |
668a0654 | 541 | if(HasDump3D()){ |
542 | TDirectory *cwd = gDirectory; | |
543 | TFile::Open(Form("EffDump_%s.root", H->GetName()), "RECREATE"); | |
544 | for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){ | |
545 | if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue; | |
546 | if(!pr0->H()) continue; | |
547 | TH3 *h3=(TH3*)pr0->H()->Clone(); | |
548 | h3->Write(); | |
549 | } | |
550 | gFile->Close(); | |
551 | cwd->cd(); | |
552 | } | |
553 | ||
0f2c4c4f | 554 | Int_t jh(0); |
eb05d549 | 555 | for(; ih--; ){ |
556 | if(!hp[ih].H()) continue; | |
5047978d | 557 | for(Int_t ipt(0); ipt<=fNpt; ipt++) fProj->AddAt(hp[ih].Projection2Dbin(ipt), jh++); |
eb05d549 | 558 | } |
559 | ||
0f2c4c4f | 560 | /* AliTRDrecoProjection prLab; TH2 *hLab[3] = {0}; TH1 *hpLab[3] = {0}; |
eb05d549 | 561 | for(ilab=0; ilab<nlab; ilab++){ |
562 | if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, 3)))) continue; | |
563 | prLab=(*pr0); | |
564 | prLab.SetNameTitle(Form("HEffLb%d", ilab), "Sum over status"); | |
565 | prLab.H()->SetNameTitle(Form("HEffLb%d", ilab), Form("Efficiency :: #bf{%s} Propagated Tracks", lab[ilab])); | |
566 | if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, 4)))) continue; | |
567 | prLab+=(*pr1); | |
4b92200c | 568 | prLab.Projection2D(1, 10, -1, kFALSE); |
eb05d549 | 569 | if((hLab[ilab] = (TH2*)gDirectory->Get(Form("%sEn", prLab.H()->GetName())))) fProj->AddAt(hLab[ilab], jh++); |
570 | if((hpLab[ilab] = prLab.H()->Project3D("z"))) fProj->AddAt(hpLab[ilab], jh++); | |
0f2c4c4f | 571 | }*/ |
eb05d549 | 572 | |
573 | for(Int_t istat(0); istat<5; istat++) { | |
574 | if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, istat)))) { | |
0f2c4c4f | 575 | // sum over MC labels if available |
eb05d549 | 576 | for(ilab=1; ilab<nlab; ilab++){ |
577 | if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, istat)))) continue; | |
578 | (*pr0)+=(*pr1); | |
579 | } | |
580 | pr0->H()->SetNameTitle(Form("HEff%d", istat), Form("Efficiency :: Stat[#bf{%s}]", stat[istat])); | |
5047978d | 581 | for(Int_t ipt(0); ipt<=fNpt; ipt++) fProj->AddAt(pr0->Projection2Dbin(ipt), jh++); |
eb05d549 | 582 | |
583 | if(istat>1 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff01"))) (*pr1)+=(*pr0); | |
584 | if(istat>2 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff02"))) (*pr1)+=(*pr0); | |
585 | if(istat>3 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff03"))) (*pr1)+=(*pr0); | |
586 | } | |
587 | } | |
0f2c4c4f | 588 | // Project 2D tracks |
589 | const char suffix[] = {'A', 'T', 'P'}; | |
590 | const char *sname[] = {"All", "Trk", "Prp"}; | |
591 | for(Int_t istat(0); istat<3; istat++){ | |
592 | if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, istat+1)))) continue; | |
593 | pr0->H()->SetNameTitle(Form("HEff%c", suffix[istat]), Form("Efficiency :: %s Tracks", sname[istat])); | |
5047978d | 594 | for(Int_t ipt(-1); ipt<=fNpt; ipt++) fProj->AddAt(pr0->Projection2Dbin(ipt), jh++); |
eb05d549 | 595 | } |
0f2c4c4f | 596 | |
597 | // Efficiency | |
598 | TH2F *h2T(NULL), *h2P(NULL); | |
599 | for(Int_t ipt(-1); ipt<=fNpt; ipt++){ | |
600 | if(!(h2T = (TH2F*)fProj->FindObject(Form("HEffT%d_2D", ipt)))) continue; | |
601 | if(!(h2P = (TH2F*)fProj->FindObject(Form("HEffP%d_2D", ipt)))) continue; | |
602 | h2P->Divide(h2T); | |
603 | PutTrendValue(ipt<0?"pt":(Form("pt%d", ipt)), GetMeanStat(h2P, 0.01, ">")); | |
eb05d549 | 604 | } |
0f2c4c4f | 605 | /* // process MC label |
eb05d549 | 606 | if(hEff[2]){ |
607 | for(ilab=0; ilab<nlab; ilab++){ | |
608 | if(!hLab[ilab]) continue; | |
668a0654 | 609 | // remove fakes |
610 | TH2 *hEff1 = (TH2*)hLab[ilab]->Clone(Form("%sN", hLab[ilab]->GetName())); | |
611 | for(Int_t ix(1); ix<=hLab[ilab]->GetNbinsX(); ix++){ | |
612 | for(Int_t iy(1); iy<=hLab[ilab]->GetNbinsY(); iy++){ | |
613 | if(hLab[ilab]->GetBinContent(ix, iy)<5) hEff1->SetBinContent(ix, iy, 0.); | |
614 | }} | |
615 | hEff1->Divide(hEff[2]); | |
616 | fProj->AddAt(hEff1, jh++); | |
eb05d549 | 617 | } |
618 | } | |
619 | if(hpEff[2]){ | |
620 | for(ilab=0; ilab<nlab; ilab++){ | |
621 | if(!hpLab[ilab]) continue; | |
668a0654 | 622 | TH1 *hpEff1 = (TH1*)hpLab[ilab]->Clone(Form("%sN", hpLab[ilab]->GetName())); |
623 | hpEff1->Divide(hpEff[2]); | |
624 | fProj->AddAt(hpEff1, jh++); | |
eb05d549 | 625 | } |
0f2c4c4f | 626 | }*/ |
eb05d549 | 627 | AliInfo(Form("Done %3d 2D projections.", jh)); |
1ee39b3a | 628 | return kTRUE; |
629 | } | |
eb05d549 | 630 | |
631 | //____________________________________________________________________ | |
632 | void AliTRDefficiency::MakeSummary() | |
633 | { | |
634 | // Build summary plots | |
635 | if(!fProj){ | |
636 | AliError("Missing results"); | |
637 | return; | |
638 | } | |
639 | TVirtualPad *p(NULL); TCanvas *cOut(NULL); | |
640 | TH2 *h2(NULL); | |
641 | gStyle->SetPalette(1); | |
642 | ||
4b92200c | 643 | Int_t nbins(DebugLevel()==0?3:20); |
eb05d549 | 644 | //calculate true pt bin |
645 | Float_t ptBins[23]; ptBins[0] = 0.; | |
646 | if(nbins==3){ // significant bins | |
647 | ptBins[1] = 0.5; | |
648 | ptBins[2] = 0.8; | |
649 | ptBins[3] = 1.5; | |
650 | ptBins[4] = 5.; | |
651 | ptBins[5] = 10.; | |
652 | } else if(nbins==20){ // logarithmic bins | |
653 | ptBins[1] = 0.5; | |
654 | Float_t dpt(0.002); | |
655 | for(Int_t ipt(1); ipt<21; ipt++) ptBins[ipt+1] = ptBins[ipt]+(TMath::Exp(ipt*ipt*dpt)-1.); | |
656 | ptBins[22] = 10.; | |
657 | } else { | |
658 | AliError(Form("Unknown no.[%d] of bins in the p_t spectrum", nbins)); | |
659 | return;// kFALSE; | |
660 | } | |
661 | ||
33056e04 | 662 | const Char_t cid[]={'T','P'}; |
eb05d549 | 663 | cOut = new TCanvas(Form("%s_Eff", GetName()), "TRD Efficiency", 1536, 1536); cOut->Divide(2,2,1.e-5,1.e-5); |
664 | // tracking eff :: eta-phi dependence | |
665 | for(Int_t it(0); it<2; it++){ | |
668a0654 | 666 | if(!(h2 = (TH2*)fProj->FindObject(Form("HEff%cEnN", cid[it])))) { |
667 | AliError(Form("Missing \"HEff%cEnN\".", cid[it])); | |
eb05d549 | 668 | continue; |
669 | } | |
33056e04 | 670 | h2->SetContour(10); h2->Scale(1.e2); SetRangeZ(h2, 80, 100, 5); |
eb05d549 | 671 | h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle(); |
672 | p=cOut->cd(it+1);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712); | |
673 | h2->Draw("colz"); | |
33056e04 | 674 | MakeDetectorPlot(); |
eb05d549 | 675 | } |
676 | if(!(h2 = (TH2*)fProj->FindObject("HEff0En"))) { | |
677 | AliError("Missing \"HEff0En\"."); | |
678 | return; | |
679 | } | |
680 | p=cOut->cd(3);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712); | |
33056e04 | 681 | h2->Draw("colz"); MakeDetectorPlot(); |
eb05d549 | 682 | // tracking eff :: pt dependence |
683 | TH1 *h[2] = {0}; | |
668a0654 | 684 | if(!(h[0] = (TH1*)fProj->FindObject("HEffP_zN"))){ |
685 | AliError("Missing \"HEffP_zN\"."); | |
eb05d549 | 686 | return; |
687 | } | |
668a0654 | 688 | if(!(h[1] = (TH1*)fProj->FindObject("HEffT_zN"))){ |
689 | AliError("Missing \"HEffT_zN\"."); | |
eb05d549 | 690 | return; |
691 | } | |
692 | TH1 *h1[3] = {0}; | |
693 | Color_t color[] = {kGreen, kBlue, kRed}; | |
694 | for(Int_t il=0;il<3;il++){ | |
695 | h1[il]=new TH1F(Form("h1Eff%d", il), "", nbins+2, ptBins); | |
696 | h1[il]->SetFillColor(color[il]); | |
697 | h1[il]->SetFillStyle(il==2?3002:1001); | |
698 | h1[il]->SetLineColor(color[il]); | |
668a0654 | 699 | h1[il]->SetMarkerStyle(4); |
700 | h1[il]->SetMarkerColor(color[il]); | |
eb05d549 | 701 | h1[il]->SetLineWidth(1); |
702 | } | |
703 | for(Int_t ip(0);ip<=(nbins+1);ip++){ | |
704 | h1[0]->SetBinContent(ip+1, 1.e2*h[0]->GetBinContent(ip)); // propagated | |
705 | h1[1]->SetBinContent(ip+1, 1.e2*(h[1]->GetBinContent(ip) - h[0]->GetBinContent(ip))); // stopped | |
706 | h1[2]->SetBinContent(ip+1, 1.e2*(1 - h[1]->GetBinContent(ip))); // missed | |
707 | } | |
33056e04 | 708 | const Char_t *labEff[] = {"Propagated", "Stopped", "Missed"}; |
eb05d549 | 709 | THStack *hs = new THStack("hEff","Tracking Efficiency;p_{t} [GeV/c];Efficiency [%]"); |
710 | TLegend *leg = new TLegend(0.671371,0.1313559,0.9576613,0.2923729,NULL,"brNDC"); | |
711 | leg->SetBorderSize(0); leg->SetFillColor(kWhite); leg->SetFillStyle(1001); | |
712 | for(Int_t ic(0); ic<3;ic++){ hs->Add(h1[ic]);leg->AddEntry(h1[ic], labEff[ic], "f");} | |
713 | p=cOut->cd(4); p->SetLeftMargin(0.08266129); p->SetRightMargin(0.0141129);p->SetTopMargin(0.006355932);p->SetLogx(); | |
668a0654 | 714 | hs->Draw(); //hs->Draw("same nostack,e1p"); |
715 | leg->Draw(); | |
eb05d549 | 716 | hs->GetXaxis()->SetRangeUser(0.6,10.); |
717 | hs->GetXaxis()->SetMoreLogLabels(); | |
718 | hs->GetXaxis()->SetTitleOffset(1.2); | |
719 | hs->GetYaxis()->SetNdivisions(513); | |
668a0654 | 720 | hs->SetMinimum(75.); |
eb05d549 | 721 | hs->GetYaxis()->CenterTitle(); |
722 | cOut->SaveAs(Form("%s.gif", cOut->GetName())); | |
723 | ||
724 | cOut = new TCanvas(Form("%s_MC", GetName()), "TRD Label", 1536, 1536); cOut->Divide(2,2,1.e-5,1.e-5); | |
725 | for(Int_t ipad(0); ipad<3; ipad++){ | |
726 | p=cOut->cd(ipad+1);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712); | |
668a0654 | 727 | if(!(h2 = (TH2*)fProj->FindObject(Form("HEffLb%dEnN", ipad)))) continue; |
eb05d549 | 728 | h2->SetContour(10); |
668a0654 | 729 | h2->Scale(1.e2); SetRangeZ(h2, ipad==1?50:0., ipad==1?90.:50., ipad==1?0.01:0.01); |
eb05d549 | 730 | h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle(); |
731 | h2->Draw("colz"); | |
33056e04 | 732 | MakeDetectorPlot(); |
eb05d549 | 733 | } |
33056e04 | 734 | color[0] = kRed; color[1] = kGreen; color[2] = kBlue; |
eb05d549 | 735 | for(Int_t il=0;il<3;il++){ |
668a0654 | 736 | if(!(h[il] = (TH1D*)fProj->FindObject(Form("HEffLb%d_zN", il)))) continue; |
eb05d549 | 737 | h1[il]=new TH1F(Form("h1Lab%d", il), "", nbins+2, ptBins); |
738 | for(Int_t ip(0);ip<=(nbins+1);ip++) h1[il]->SetBinContent(ip+1, 1.e2*h[il]->GetBinContent(ip)); | |
33056e04 | 739 | h1[il]->SetFillColor(color[il]); |
eb05d549 | 740 | h1[il]->SetFillStyle(il==2?3002:1001); |
33056e04 | 741 | h1[il]->SetLineColor(color[il]); |
eb05d549 | 742 | h1[il]->SetLineWidth(1); |
743 | } | |
33056e04 | 744 | const Char_t *labMC[] = {"TRD != ESD [bad]", "TRD == ESD [good]", "TRD == -ESD [accept]"}; |
eb05d549 | 745 | leg = new TLegend(0.671371,0.1313559,0.9576613,0.2923729,NULL,"brNDC"); |
746 | leg->SetBorderSize(0); leg->SetFillColor(kWhite); leg->SetFillStyle(1001); | |
747 | hs = new THStack("hLab","TRD Label;p_{t} [GeV/c];Efficiency [%]"); | |
748 | hs->Add(h1[1]);leg->AddEntry(h1[1], labMC[1], "f"); // good | |
749 | hs->Add(h1[2]);leg->AddEntry(h1[2], labMC[2], "f"); // accept | |
750 | hs->Add(h1[0]);leg->AddEntry(h1[0], labMC[0], "f"); // bad | |
751 | p=cOut->cd(4); p->SetLeftMargin(0.08266129); p->SetRightMargin(0.0141129);p->SetTopMargin(0.006355932); p->SetLogx(); | |
668a0654 | 752 | hs->Draw(/*"nostack,e1p"*/); leg->Draw(); |
eb05d549 | 753 | cOut->Modified();cOut->Update(); |
754 | hs->GetXaxis()->SetRangeUser(0.6,10.); | |
755 | hs->GetXaxis()->SetMoreLogLabels(); | |
756 | hs->GetXaxis()->SetTitleOffset(1.2); | |
757 | hs->GetYaxis()->SetNdivisions(513); | |
668a0654 | 758 | hs->SetMinimum(50.); |
eb05d549 | 759 | hs->GetYaxis()->CenterTitle(); |
760 | cOut->SaveAs(Form("%s.gif", cOut->GetName())); | |
761 | } |