]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/AliTRDcheckESD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
train ready for Pilot integration with all wagons
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckESD.cxx
CommitLineData
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/////////////////////////////////////////////////////
17//
18// Check basic detector results at ESD level
19// - Geometrical efficiency
20// - Tracking efficiency
21// - PID efficiency
22// - Refit efficiency
23//
24// Author
25// Alex Bercuci <A.Bercuci@gsi.de>
26//
27//////////////////////////////////////////////////////
28
29#include <TClonesArray.h>
30#include <TObjArray.h>
629ae9b8 31#include <TPad.h>
32#include <TLegend.h>
33#include <TF1.h>
1ee39b3a 34#include <TH2I.h>
fb09668f 35#include <TH3S.h>
1ee39b3a 36#include <TGraphErrors.h>
37#include <TGraphAsymmErrors.h>
38#include <TFile.h>
39#include <TTree.h>
40#include <TROOT.h>
41#include <TChain.h>
42#include <TParticle.h>
43
44#include "AliLog.h"
45#include "AliAnalysisManager.h"
46#include "AliESDEvent.h"
47#include "AliESDkink.h"
48#include "AliMCEvent.h"
49#include "AliESDInputHandler.h"
50#include "AliMCEventHandler.h"
51
52#include "AliESDtrack.h"
53#include "AliMCParticle.h"
54#include "AliPID.h"
55#include "AliStack.h"
56#include "AliTrackReference.h"
57
58#include "AliTRDcheckESD.h"
59
60ClassImp(AliTRDcheckESD)
61
62const Float_t AliTRDcheckESD::fgkxTPC = 290.;
63const Float_t AliTRDcheckESD::fgkxTOF = 365.;
629ae9b8 64const UChar_t AliTRDcheckESD::fgkNgraph[kNrefs] ={
650, 4, 2, 10};
9b496b5e 66FILE* AliTRDcheckESD::fgFile = NULL;
1ee39b3a 67
68//____________________________________________________________________
69AliTRDcheckESD::AliTRDcheckESD():
f8f46e4d 70 AliAnalysisTaskSE()
1ee39b3a 71 ,fStatus(0)
629ae9b8 72 ,fNRefFigures(0)
9b496b5e 73 ,fESD(NULL)
74 ,fMC(NULL)
75 ,fHistos(NULL)
76 ,fResults(NULL)
1ee39b3a 77{
78 //
79 // Default constructor
80 //
db99a57a 81 SetNameTitle("checkESD", "Check TRD @ ESD level");
1ee39b3a 82 SetMC(kTRUE);
f8f46e4d 83}
84
f2e89a4c 85//____________________________________________________________________
f8f46e4d 86AliTRDcheckESD::AliTRDcheckESD(char* name):
87 AliAnalysisTaskSE(name)
88 ,fStatus(0)
629ae9b8 89 ,fNRefFigures(0)
f8f46e4d 90 ,fESD(NULL)
91 ,fMC(NULL)
92 ,fHistos(NULL)
93 ,fResults(NULL)
94{
95 //
96 // Default constructor
97 //
98 SetMC(kTRUE);
f2e89a4c 99 SetTitle("Check TRD @ ESD level");
f8f46e4d 100 DefineOutput(1, TObjArray::Class());
1ee39b3a 101}
102
103//____________________________________________________________________
104AliTRDcheckESD::~AliTRDcheckESD()
105{
106// Destructor
107 if(fHistos){
108 //fHistos->Delete();
109 delete fHistos;
110 }
111 if(fResults){
112 fResults->Delete();
113 delete fResults;
114 }
115}
116
117//____________________________________________________________________
f8f46e4d 118void AliTRDcheckESD::UserCreateOutputObjects()
1ee39b3a 119{
120 //
121 // Create Output Containers (TObjectArray containing 1D histograms)
122 //
a96ac33a 123 //OpenFile(0, "RECREATE");
124
1ee39b3a 125 Histos();
126}
127
629ae9b8 128
1ee39b3a 129//____________________________________________________________________
629ae9b8 130Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
1ee39b3a 131{
629ae9b8 132 if(ifig>=fNRefFigures){
133 AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));
134 return kFALSE;
1ee39b3a 135 }
629ae9b8 136 if(!gPad){
137 AliWarning("Please provide a canvas to draw results.");
138 return kFALSE;
1ee39b3a 139 }
140
629ae9b8 141 const Char_t *title[20];
142 TH1 *hF(NULL);
143 if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;
144 TLegend *leg(NULL);
145 TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
146 TObjArray *arr(NULL);
147 switch(ifig){
148 case kNCl: // number of clusters/track
149 ((TH1I*)fResults->At(kNCl))->Draw("c");
150 break;
151 case kTRDstat: // Efficiency
152 if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;
153 leg = new TLegend(.65, .7, .95, .99);
154 leg->SetHeader("TRD Efficiency");
155 leg->SetBorderSize(1); leg->SetFillColor(0);
156 title[0] = "Geometrical (TRDin/TPCout)";
157 title[1] = "Tracking (TRDout/TRDin)";
158 title[2] = "PID (TRDpid/TRDin)";
159 title[3] = "Refit (TRDrefit/TRDin)";
160 hF = new TH1S("hFcheckESD", ";p [GeV/c];Efficiency", 10, 0.1, 10.);
161 hF->SetMaximum(1.3);
162 hF->GetXaxis()->SetMoreLogLabels();
163 hF->Draw("p");
164 for(Int_t ig(0); ig<fgkNgraph[kTRDstat]; ig++){
165 if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
166 g->Draw("pl"); leg->AddEntry(g, title[ig], "pl");
167 //PutTrendValue(name[id], g->GetMean(2));
168 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
1ee39b3a 169 }
629ae9b8 170 leg->Draw(); gPad->SetLogx();
171 break;
172 case kTRDmom: // Energy loss
173 if(!(arr = (TObjArray*)fResults->At(kTRDmom))) return kFALSE;
174 leg = new TLegend(.65, .7, .95, .99);
175 leg->SetHeader("Energy Loss");
176 leg->SetBorderSize(1); leg->SetFillColor(0);
177 title[0] = "Max & 90% quantile";
178 title[1] = "Mean & 60% quantile";
179 hF = new TH1S("hFcheckESD", ";layer;#Delta E", 6, -0.5, 5.5);
180 hF->SetMaximum(1.3);hF->SetMinimum(-.3);
181 hF->Draw("p");
182 for(Int_t ig(0); ig<fgkNgraph[kTRDmom]; ig++){
183 if(!(ga = (TGraphAsymmErrors*)arr->At(ig))) return kFALSE;
184 ga->Draw("pl"); leg->AddEntry(ga, title[ig], "pl");
185 //PutTrendValue(name[id], g->GetMean(2));
186 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
1ee39b3a 187 }
629ae9b8 188 leg->Draw();gPad->SetLogx(kFALSE);
189 break;
190 case kPtRes: // Pt resolution @ vertex
191 if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;
192 gPad->SetMargin(0.1, 0.22, 0.1, 0.023);
193 leg = new TLegend(.78, .1, .99, .98);
194 leg->SetHeader("P_{t} @ DCA");
195 leg->SetBorderSize(1); leg->SetFillColor(0);
196 leg->SetTextAlign(22);
197 leg->SetTextFont(12);
198 leg->SetTextSize(0.03813559);
199 hF = new TH1S("hFcheckESD", ";p_{t} [GeV/c];#Delta p_{t} / p_{t} [%]", 10, 0.2, 10.);
200 hF->SetMaximum(10.);hF->SetMinimum(-3.);
201 hF->GetXaxis()->SetMoreLogLabels();
202 hF->GetXaxis()->SetTitleOffset(1.2);
203 hF->GetYaxis()->CenterTitle();
204 hF->Draw("p");
205 for(Int_t ig(0); ig<fgkNgraph[kPtRes]; ig++){
206 if(!(g = (TGraphErrors*)arr->At(ig))) continue;//return kFALSE;
207 if(!g->GetN()) continue;
208 g->Draw("pl"); leg->AddEntry(g, g->GetTitle(), "pl");
209 //PutTrendValue(name[id], g->GetMean(2));
210 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
211 }
212 leg->Draw();gPad->SetLogx();
213 break;
1ee39b3a 214 }
e2e3cec2 215 return kTRUE;
216}
217
1ee39b3a 218//____________________________________________________________________
f8f46e4d 219void AliTRDcheckESD::UserExec(Option_t *){
1ee39b3a 220 //
221 // Run the Analysis
222 //
f8f46e4d 223 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
224 fMC = MCEvent();
705f8b0a 225
1ee39b3a 226 if(!fESD){
227 AliError("ESD event missing.");
228 return;
229 }
f8f46e4d 230
1ee39b3a 231 // Get MC information if available
9b496b5e 232 AliStack * fStack = NULL;
1ee39b3a 233 if(HasMC()){
234 if(!fMC){
235 AliWarning("MC event missing");
236 SetMC(kFALSE);
237 } else {
238 if(!(fStack = fMC->Stack())){
239 AliWarning("MC stack missing");
240 SetMC(kFALSE);
241 }
242 }
243 }
9b496b5e 244 TH2 *h(NULL);
245
246 AliESDtrack *esdTrack(NULL);
1ee39b3a 247 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
1ee39b3a 248 esdTrack = fESD->GetTrack(itrk);
249
1ee39b3a 250 // track status
fb09668f 251 ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);
1ee39b3a 252 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
253 if(esdTrack->GetKinkIndex(0) > 0) continue;
254
fb09668f 255 //Int_t nTPC(esdTrack->GetNcls(1));
256 Int_t nTRD(esdTrack->GetNcls(2));
257 Double_t pt(esdTrack->Pt());
258 //Double_t eta(esdTrack->Eta());
259 //Double_t phi(esdTrack->Phi());
1ee39b3a 260 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
261 // pid quality
262 //esdTrack->GetTRDntrackletsPID();
fb09668f 263 Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);
1ee39b3a 264
265 // look at external track param
266 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
267 const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
1ee39b3a 268
fb09668f 269 Double_t pt0(0.), eta0(0.), phi0(0.), ptTRD(0.);
9b496b5e 270 // read MC info if available
fb09668f 271 Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);
272 AliMCParticle *mcParticle(NULL);
9b496b5e 273 if(HasMC()){
274 AliTrackReference *ref(NULL);
275 Int_t fLabel(esdTrack->GetLabel());
fb09668f 276 Int_t fIdx(TMath::Abs(fLabel));
277 if(fIdx > fStack->GetNtrack()) continue;
9b496b5e 278
fb09668f 279 // read MC particle
280 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {
9b496b5e 281 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
282 continue;
283 }
fb09668f 284 pt0 = mcParticle->Pt();
285 eta0 = mcParticle->Eta();
286 phi0 = mcParticle->Phi();
287 kPhysPrim = fMC->IsPhysicalPrimary(fIdx);
288
289 // read track references
9b496b5e 290 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
291 if(!nRefs){
fb09668f 292 AliWarning(Form("No TR found for track @ Label[%d].", fLabel));
9b496b5e 293 continue;
294 }
295 Int_t iref = 0;
296 while(iref<nRefs){
297 ref = mcParticle->GetTrackReference(iref);
298 if(ref->LocalX() > fgkxTPC) break;
299 ref=NULL; iref++;
300 }
301 if(ref){
302 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
303 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
304 }
305 } else { // track stopped in TPC
306 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
307 }
fb09668f 308 ptTRD = ref->Pt();kFOUND=kTRUE;
9b496b5e 309 } else { // use reconstructed values
310 if(op){
311 Double_t x(op->GetX());
312 if(x<fgkxTOF && x>fgkxTPC){
fb09668f 313 ptTRD=op->Pt();
9b496b5e 314 kFOUND=kTRUE;
315 }
316 }
1ee39b3a 317
9b496b5e 318 if(!kFOUND && ip){
fb09668f 319 ptTRD=ip->Pt();
9b496b5e 320 kFOUND=kTRUE;
321 }
1ee39b3a 322 }
323
9b496b5e 324 if(kFOUND){
325 h = (TH2I*)fHistos->At(kTRDstat);
fb09668f 326 if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);
327 if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);
328 if(kBarrel && (status & AliESDtrack::kTRDout)){
329 ((TH1*)fHistos->At(kNCl))->Fill(nTRD);
330 h->Fill(ptTRD, kTRDout);
1ee39b3a 331 }
fb09668f 332 if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);
333 if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);
334 }
335 if(HasMC() && kBarrel && (status & AliESDtrack::kTRDout)) {
336 TH3 *h3 = (TH3S*)fHistos->At(kPtRes);
337 Int_t sgn = mcParticle->Charge()>0?1:-1;
338 h3->Fill(pt0, 1.e2*pt/pt0-1.e2, sgn*Pdg2Idx(TMath::Abs(mcParticle->PdgCode())));
1ee39b3a 339 }
1ee39b3a 340 if(ip){
341 h = (TH2I*)fHistos->At(kTRDmom);
fb09668f 342 Float_t pTRD(0.);
1ee39b3a 343 for(Int_t ily=6; ily--;){
fb09668f 344 if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;
345 h->Fill(ip->GetP()-pTRD, ily);
1ee39b3a 346 }
347 }
348 }
f8f46e4d 349 PostData(1, fHistos);
1ee39b3a 350}
351
352//____________________________________________________________________
353TObjArray* AliTRDcheckESD::Histos()
354{
355// Retrieve histograms array if already build or build it
356
357 if(fHistos) return fHistos;
358
359 fHistos = new TObjArray(kNhistos);
360 //fHistos->SetOwner(kTRUE);
629ae9b8 361
9b496b5e 362 TH1 *h = NULL;
1ee39b3a 363
364 // clusters per tracklet
365 if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
fb09668f 366 h = new TH1I("hNCl", "Clusters per TRD track;N_{cl}^{TRD};entries", 100, 0., 200.);
1ee39b3a 367 } else h->Reset();
629ae9b8 368 fHistos->AddAt(h, kNCl); fNRefFigures++;
1ee39b3a 369
370 // status bits histogram
fb09668f 371 const Int_t kNpt(10), kNbits(5);
372 Float_t Pt(0.1), Bits(.5);
373 Float_t binsPt[kNpt+1], binsBits[kNbits+1];
374 for(Int_t i=0;i<kNpt+1; i++,Pt+=(TMath::Exp(i*i*.015)-1.)) binsPt[i]=Pt;
375 for(Int_t i=0; i<kNbits+1; i++,Bits+=1.) binsBits[i]=Bits;
1ee39b3a 376 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
629ae9b8 377 h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);
fb09668f 378 TAxis *ay(h->GetYaxis());
379 ay->SetBinLabel(1, "kTPCout");
380 ay->SetBinLabel(2, "kTRDin");
381 ay->SetBinLabel(3, "kTRDout");
382 ay->SetBinLabel(4, "kTRDpid");
383 ay->SetBinLabel(5, "kTRDrefit");
1ee39b3a 384 } else h->Reset();
385 fHistos->AddAt(h, kTRDstat);
386
387 // energy loss
388 if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
fb09668f 389 h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);
1ee39b3a 390 } else h->Reset();
391 fHistos->AddAt(h, kTRDmom);
629ae9b8 392 if(!HasMC()) return fHistos;
1ee39b3a 393
fb09668f 394 // pt resolution
395 const Int_t kNdpt(100), kNspec(2*AliPID::kSPECIES+1);
396 Float_t DPt(-3.), Spec(-AliPID::kSPECIES-0.5);
397 Float_t binsDPt[kNdpt+1], binsSpec[kNspec+1];
398 for(Int_t i=0; i<kNdpt+1; i++,DPt+=6.e-2) binsDPt[i]=DPt;
399 for(Int_t i=0; i<kNspec+1; i++,Spec+=1.) binsSpec[i]=Spec;
400 if(!(h = (TH3S*)gROOT->FindObject("hPtRes"))){
401 h = new TH3S("hPtRes", "P_{t} resolution @ DCA;p_{t}^{MC} [GeV/c];#Delta p_{t}/p_{t}^{MC} [%];SPECIES", kNpt, binsPt, kNdpt, binsDPt, kNspec, binsSpec);
402 TAxis *az(h->GetZaxis());
403 for(Int_t i(0); i<AliPID::kSPECIES; i++){
404 az->SetBinLabel(5-i, AliPID::ParticleLatexName(i));
405 az->SetBinLabel(7+i, AliPID::ParticleLatexName(i));
406 }
407 } else h->Reset();
408 fHistos->AddAt(h, kPtRes);
409
1ee39b3a 410 return fHistos;
411}
412
413//____________________________________________________________________
414Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
415{
416// Load data from performance file
417
418 if(!TFile::Open(filename)){
419 AliWarning(Form("Couldn't open file %s.", filename));
420 return kFALSE;
421 }
9b496b5e 422 TObjArray *o = NULL;
1ee39b3a 423 if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
424 AliWarning("Missing histogram container.");
425 return kFALSE;
426 }
427 fHistos = (TObjArray*)o->Clone(GetName());
428 gFile->Close();
a96ac33a 429 SETBIT(fStatus, kLoad);
1ee39b3a 430 return kTRUE;
431}
432
433//_______________________________________________________
434Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
435{
436// Dump trending value to default file
437
438 if(!fgFile){
439 fgFile = fopen("TRD.Performance.txt", "at");
440 }
441 fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
442 return kTRUE;
443}
444
445//____________________________________________________________________
446void AliTRDcheckESD::Terminate(Option_t *)
447{
448// Steer post-processing
a96ac33a 449 if(!IsLoad()){
2efff183 450 fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));
a96ac33a 451 if(!fHistos){
452 AliError("Histogram container not found in output");
453 return;
454 }
9b496b5e 455 }
1ee39b3a 456
629ae9b8 457 const Char_t *name[kNrefs] = {
458 "Ncl", "Eff", "Eloss", "PtResDCA"
459 };
460 TObjArray *arr(NULL); TGraph *g(NULL);
461 if(!fResults){
462 fResults = new TObjArray(kNrefs);
463 fResults->SetOwner();
464 fResults->SetName("results");
465 for(Int_t iref(0); iref<kNrefs; iref++){
466 fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);
467 arr->SetName(name[iref]); arr->SetOwner();
468 switch(iref){
469 case kTRDmom:
470 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
471 arr->AddAt(g = new TGraphAsymmErrors(), ig);
472 g->SetLineColor(ig+1);
473 g->SetMarkerColor(ig+1);
474 g->SetMarkerStyle(ig+20);
475 }
476 break;
477 case kPtRes:
478 for(Int_t im(fgkNgraph[iref]/2), ig(im), idx(0); ig--; idx+=ig){
479 arr->AddAt(g = new TGraphErrors(), ig);
480 g->SetLineColor(kRed-idx);
481 g->SetMarkerColor(kRed-idx);
482 g->SetMarkerStyle(20+ig);
483 g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(ig)));
484 arr->AddAt(g = new TGraphErrors(), im+ig);
485 g->SetLineColor(kBlue-idx);
486 g->SetMarkerColor(kBlue-idx);
487 g->SetMarkerStyle(20+ig);
488 g->SetNameTitle(Form("m%d", ig), Form("sys %s", AliPID::ParticleLatexName(ig)));
489 }
490 break;
491 default:
492 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
493 arr->AddAt(g = new TGraphErrors(), ig);
494 g->SetLineColor(ig+1);
495 g->SetMarkerColor(ig+1);
496 g->SetMarkerStyle(ig+20);
497 }
498 break;
499 }
500 }
501 }
502 fNRefFigures = 1;
503
504 // EFFICIENCY
1ee39b3a 505 // geometrical efficiency
629ae9b8 506 TH2I *h2(NULL);
507 if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;
508 arr = (TObjArray*)fResults->At(kTRDstat);
9b496b5e 509 TH1 *h1[2] = {NULL, NULL};
1ee39b3a 510 h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
511 h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
629ae9b8 512 Process(h1, (TGraphErrors*)arr->At(0));
1ee39b3a 513 delete h1[0];delete h1[1];
1ee39b3a 514 // tracking efficiency
515 h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
516 h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
629ae9b8 517 Process(h1, (TGraphErrors*)arr->At(1));
1ee39b3a 518 delete h1[1];
1ee39b3a 519 // PID efficiency
520 h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
629ae9b8 521 Process(h1, (TGraphErrors*)arr->At(2));
1ee39b3a 522 delete h1[1];
1ee39b3a 523 // Refit efficiency
524 h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
629ae9b8 525 Process(h1, (TGraphErrors*)arr->At(3));
1ee39b3a 526 delete h1[1];
629ae9b8 527 fNRefFigures++;
528
529 // ENERGY LOSS
1ee39b3a 530 if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
629ae9b8 531 arr = (TObjArray*)fResults->At(kTRDmom);
532 TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);
1ee39b3a 533 TAxis *ax=h2->GetXaxis();
534 const Int_t nq(4);
535 const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
536 Double_t yq[nq];
537 for(Int_t ily=6; ily--;){
538 h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
539 h1[0]->GetQuantiles(nq,yq,xq);
540 g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
541 g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
542 g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
543 g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
544
545 //printf(" max[%f] mean[%f] q[%f %f %f %f]\n", ax->GetBinCenter(h1[0]->GetMaximumBin()), h1[0]->GetMean(), yq[0], yq[1], yq[2], yq[3]);
546 delete h1[0];
547 }
629ae9b8 548 fNRefFigures++;
549 if(!HasMC()) return;
550
551 // Pt RESOLUTION @ DCA
552 TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};
553 if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;
554 arr = (TObjArray*)fResults->At(kPtRes);
555 TAxis *az(h3->GetZaxis());
556 for(Int_t i(0); i<AliPID::kSPECIES; i++){
557 az->SetRange(5-i, 5-i);
558 gg[1] = (TGraphErrors*)arr->At(4-i);
559 gg[0] = (TGraphErrors*)arr->At(9-i);
560 Process2D((TH2*)h3->Project3D("yx"), gg);
561 //az->SetRange(7+i, 7+i);
562 //Process2D((TH2*)h3->Project3D("yx"), gg);
563 }
564 fNRefFigures++;
1ee39b3a 565}
566
fb09668f 567//____________________________________________________________________
568Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg)
569{
570 switch(pdg){
571 case kElectron: return AliPID::kElectron+1;
572 case kMuonMinus: return AliPID::kMuon+1;
573 case kPiPlus: return AliPID::kPion+1;
574 case kKPlus: return AliPID::kKaon+1;
575 case kProton: return AliPID::kProton+1;
576 }
577 return 0;
578}
579
1ee39b3a 580//____________________________________________________________________
581void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
582{
583// Generic function to process one reference plot
584
585 Int_t n1 = 0, n2 = 0, ip=0;
586 Double_t eff = 0.;
587
588 TAxis *ax = h1[0]->GetXaxis();
589 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
590 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
591 n2 = (Int_t)h1[1]->GetBinContent(ib);
592 eff = n2/Float_t(n1);
593
594 ip=g->GetN();
595 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
596 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
597 }
598}
629ae9b8 599//________________________________________________________
600void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)
601{
602 //
603 // Do the processing
604 //
1ee39b3a 605
629ae9b8 606 Int_t n = 0;
607 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
608 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
609 TF1 f("fg", "gaus", -3.,3.);
610 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
611 Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
612 TH1D *h = h2->ProjectionY("py", ibin, ibin);
613 if(h->GetEntries()<100) continue;
614 //AdjustF1(h, f);
615
616 h->Fit(&f, "QN");
617 Int_t ip = g[0]->GetN();
618 g[0]->SetPoint(ip, x, f.GetParameter(1));
619 g[0]->SetPointError(ip, 0., f.GetParError(1));
620 g[1]->SetPoint(ip, x, f.GetParameter(2));
621 g[1]->SetPointError(ip, 0., f.GetParError(2));
622 }
623 return;
624}
1ee39b3a 625//____________________________________________________________________
626void AliTRDcheckESD::PrintStatus(ULong_t status)
627{
628// Dump track status to stdout
629
630 printf("ITS[i(%d) o(%d) r(%d)] TPC[i(%d) o(%d) r(%d) p(%d)] TRD[i(%d) o(%d) r(%d) p(%d) s(%d)] HMPID[o(%d) p(%d)]\n"
631 ,Bool_t(status & AliESDtrack::kITSin)
632 ,Bool_t(status & AliESDtrack::kITSout)
633 ,Bool_t(status & AliESDtrack::kITSrefit)
634 ,Bool_t(status & AliESDtrack::kTPCin)
635 ,Bool_t(status & AliESDtrack::kTPCout)
636 ,Bool_t(status & AliESDtrack::kTPCrefit)
637 ,Bool_t(status & AliESDtrack::kTPCpid)
638 ,Bool_t(status & AliESDtrack::kTRDin)
639 ,Bool_t(status & AliESDtrack::kTRDout)
640 ,Bool_t(status & AliESDtrack::kTRDrefit)
641 ,Bool_t(status & AliESDtrack::kTRDpid)
642 ,Bool_t(status & AliESDtrack::kTRDStop)
643 ,Bool_t(status & AliESDtrack::kHMPIDout)
644 ,Bool_t(status & AliESDtrack::kHMPIDpid)
645 );
646}
647