Adding GRP/CTP/CTPtiming OCDB entry in AliRoot, from the file
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDcheckESD.cxx
CommitLineData
9a281e83 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
ddd1a033 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//////////////////////////////////////////////////////
9a281e83 28
29#include <TClonesArray.h>
30#include <TObjArray.h>
31#include <TObject.h>
32#include <TH2I.h>
48797466 33#include <TGraphErrors.h>
9a281e83 34#include <TFile.h>
35#include <TTree.h>
36#include <TROOT.h>
37#include <TChain.h>
38#include <TParticle.h>
39
40#include "AliLog.h"
41#include "AliAnalysisManager.h"
42#include "AliESDEvent.h"
d12237d6 43#include "AliESDkink.h"
9a281e83 44#include "AliMCEvent.h"
45#include "AliESDInputHandler.h"
46#include "AliMCEventHandler.h"
47
48#include "AliESDtrack.h"
49#include "AliMCParticle.h"
50#include "AliPID.h"
51#include "AliStack.h"
52#include "AliTrackReference.h"
9a281e83 53
54#include "AliTRDcheckESD.h"
55
56ClassImp(AliTRDcheckESD)
57
965e895b 58const Float_t AliTRDcheckESD::fgkxTPC = 290.;
59const Float_t AliTRDcheckESD::fgkxTOF = 365.;
59589fc9 60FILE* AliTRDcheckESD::fgFile = 0x0;
9a281e83 61
62//____________________________________________________________________
63AliTRDcheckESD::AliTRDcheckESD():
6da3eee3 64 AliAnalysisTask("checkESD", "ESD checker for TRD info")
9a281e83 65 ,fStatus(0)
66 ,fESD(0x0)
67 ,fMC(0x0)
68 ,fHistos(0x0)
69{
70 //
71 // Default constructor
72 //
49c7f0a5 73 SetMC(kTRUE);
9a281e83 74 DefineInput(0, TChain::Class());
75 DefineOutput(0, TObjArray::Class());
76}
77
78//____________________________________________________________________
79AliTRDcheckESD::~AliTRDcheckESD()
80{
ddd1a033 81// Destructor
9a281e83 82 if(fHistos){
83 //fHistos->Delete();
84 delete fHistos;
85 }
86}
87
88//____________________________________________________________________
89void AliTRDcheckESD::ConnectInputData(Option_t *)
90{
91 //
92 // Link the Input Data
93 //
94 TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
95 if(tree) tree->SetBranchStatus("Tracks", 1);
96
97 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
98 fESD = esdH ? esdH->GetEvent() : 0x0;
99
100 if(!HasMC()) return;
101 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
102 fMC = mcH ? mcH->MCEvent() : 0x0;
103}
104
105//____________________________________________________________________
106void AliTRDcheckESD::CreateOutputObjects()
107{
108 //
109 // Create Output Containers (TObjectArray containing 1D histograms)
110 //
111 OpenFile(0, "RECREATE");
c2fd5953 112 Histos();
9a281e83 113}
114
965e895b 115//____________________________________________________________________
ea4502f6 116TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
965e895b 117{
ddd1a033 118// Retrieve graph with "id"
119// Possible options are :
120// "b" - build graph if none found
121// "c" - clear existing graph
122
ea4502f6 123 Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
124 kCLEAR = strstr(opt, "c"); // clear existing graph
965e895b 125
126 const Char_t *name[] = {
127 "Geo", "Trk", "Pid", "Ref"
128 };
129 const Char_t *title[] = {
130 "TRD geometrical efficiency (TRDin/TPCout)"
131 ,"TRD tracking efficiency (TRDout/TRDin)"
132 ,"TRD PID efficiency (TRDpid/TRDin)"
133 ,"TRD refit efficiency (TRDrefit/TRDin)"
134 };
135 const Int_t ngr = sizeof(name)/sizeof(Char_t*);
ddd1a033 136 if(ngr != kNgraphs){
965e895b 137 AliWarning("No of graphs defined different from definition");
138 return 0x0;
139 }
140
141 TObjArray *res = 0x0;
142 if(!(res = (TObjArray*)fHistos->At(kResults)) ||
143 (id < 0 || id >= ngr)){
ea4502f6 144 AliWarning("Graph array missing.");
965e895b 145 return 0x0;
146 }
147
148 TGraphErrors *g = 0x0;
149 if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
59589fc9 150 if(kCLEAR){
151 for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
152 } else {
5e17c62a 153 PutTrendValue(name[id], g->GetMean(2));
154 PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
59589fc9 155 }
965e895b 156 } else {
157 if(kBUILD){
158 g = new TGraphErrors();
159 g->SetNameTitle(name[id], title[id]);
160 res->AddAt(g, id);
161 }
162 }
163 return g;
164}
165
9a281e83 166//____________________________________________________________________
167void AliTRDcheckESD::Exec(Option_t *){
168 //
169 // Run the Analysis
170 //
171 if(!fESD){
48797466 172 AliError("ESD event missing.");
9a281e83 173 return;
174 }
175
176 // Get MC information if available
177 AliStack * fStack = 0x0;
2b104c26 178 if(HasMC()){
179 if(!fMC){
180 AliWarning("MC event missing");
9a281e83 181 SetMC(kFALSE);
2b104c26 182 } else {
183 if(!(fStack = fMC->Stack())){
184 AliWarning("MC stack missing");
185 SetMC(kFALSE);
186 }
9a281e83 187 }
188 }
ddd1a033 189 Bool_t bTRDin(0), bTRDout(0), bTRDpid(0);
9a281e83 190
9a281e83 191 AliESDtrack *esdTrack = 0x0;
192 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
ddd1a033 193 bTRDin=0;bTRDout=0;bTRDpid=0;
9a281e83 194 esdTrack = fESD->GetTrack(itrk);
d12237d6 195
196// if(esdTrack->GetNcls(1)) nTPC++;
197// if(esdTrack->GetNcls(2)) nTRD++;
9a281e83 198
199 // track status
d12237d6 200 ULong_t status = esdTrack->GetStatus();
2b104c26 201 //PrintStatus(status);
d12237d6 202
203 // define TPC out tracks
204 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
205 if(esdTrack->GetKinkIndex(0) > 0) continue;
9a281e83 206
207 // TRD PID
208 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
209 // pid quality
ed15ef4f 210 //esdTrack->GetTRDntrackletsPID();
9a281e83 211
212 // look at external track param
213 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
214 Double_t xyz[3];
215 if(op){
216 op->GetXYZ(xyz);
217 op->Global2LocalPosition(xyz, op->GetAlpha());
218 //printf("op @ X[%7.3f]\n", xyz[0]);
219 }
220
221 // read MC info
222 if(!HasMC()) continue;
223
224 Int_t fLabel = esdTrack->GetLabel();
225 if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
226
227 // read MC particle
228 AliMCParticle *mcParticle = 0x0;
7aad0c47 229 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
48797466 230 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
9a281e83 231 continue;
232 }
233
234 AliTrackReference *ref = 0x0;
235 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
48797466 236 if(!nRefs){
237 AliWarning(Form("Track refs missing. Label[%d].", fLabel));
238 continue;
239 }
9a281e83 240 Int_t iref = 0;
241 while(iref<nRefs){
242 ref = mcParticle->GetTrackReference(iref);
965e895b 243 if(ref->LocalX() > fgkxTPC) break;
9a281e83 244 ref=0x0; iref++;
245 }
246
247 // read TParticle
965e895b 248 //TParticle *tParticle = mcParticle->Particle();
48797466 249 //Int_t fPdg = tParticle->GetPdgCode();
d12237d6 250 // reject secondaries
251 //if(!tParticle->IsPrimary()) continue;
9a281e83 252
d12237d6 253 if(ref){
965e895b 254 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
48797466 255 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
9a281e83 256 } else {
ddd1a033 257 bTRDin=1;
258 if(esdTrack->GetNcls(2)) bTRDout=1;
259 if(esdTrack->GetTRDntrackletsPID()>=4) bTRDpid=1;
9a281e83 260 }
d12237d6 261 } else { // track stopped in TPC
48797466 262 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
263 }
d12237d6 264 // get the MC pt !!
9a281e83 265 Float_t pt = ref->Pt();
d12237d6 266
9a281e83 267 TH2 *h = (TH2I*)fHistos->At(kTRDstat);
d12237d6 268 if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
ddd1a033 269 if(/*status & AliESDtrack::k*/bTRDin) h->Fill(pt, kTRDin);
270 if(/*status & AliESDtrack::k*/bTRDout){
9a281e83 271 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
d12237d6 272 h->Fill(pt, kTRDout);
9a281e83 273 }
ddd1a033 274 if(/*status & AliESDtrack::k*/bTRDpid) h->Fill(pt, kTRDpid);
965e895b 275 if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
9a281e83 276 }
9a281e83 277 PostData(0, fHistos);
278}
279
c2fd5953 280//____________________________________________________________________
281TObjArray* AliTRDcheckESD::Histos()
282{
ddd1a033 283// Retrieve histograms array if already build or build it
284
c2fd5953 285 if(fHistos) return fHistos;
286
287 fHistos = new TObjArray(kNhistos);
288 //fHistos->SetOwner(kTRUE);
289
290 TH1 *h = 0x0;
291
292 // clusters per tracklet
293 if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
294 h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
295 h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
296 h->GetYaxis()->SetTitle("entries");
297 } else h->Reset();
298 fHistos->AddAt(h, kNCl);
299
300 // status bits histogram
301 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
302 h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
303 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
304 h->GetYaxis()->SetTitle("status bits");
305 h->GetZaxis()->SetTitle("entries");
306 } else h->Reset();
307 fHistos->AddAt(h, kTRDstat);
308
309 // results array
310 TObjArray *res = new TObjArray();
311 res->SetName("Results");
312 fHistos->AddAt(res, kResults);
313 return fHistos;
314}
315
965e895b 316//____________________________________________________________________
317Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
318{
ddd1a033 319// Load data from performance file
320
965e895b 321 if(!TFile::Open(filename)){
322 AliWarning(Form("Couldn't open file %s.", filename));
323 return kFALSE;
324 }
325 TObjArray *o = 0x0;
326 if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
327 AliWarning("Missing histogram container.");
328 return kFALSE;
329 }
330 fHistos = (TObjArray*)o->Clone(GetName());
331 gFile->Close();
332 return kTRUE;
333}
9a281e83 334
59589fc9 335//_______________________________________________________
5e17c62a 336Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
59589fc9 337{
ddd1a033 338// Dump trending value to default file
339
59589fc9 340 if(!fgFile){
341 fgFile = fopen("TRD.Performance.txt", "at");
342 }
5e17c62a 343 fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
59589fc9 344 return kTRUE;
345}
346
9a281e83 347//____________________________________________________________________
348void AliTRDcheckESD::Terminate(Option_t *)
349{
ddd1a033 350// Steer post-processing
351
965e895b 352 TObjArray *res = 0x0;
353 if(!(res = (TObjArray*)fHistos->At(kResults))){
354 AliWarning("Graph container missing.");
355 return;
356 }
ddd1a033 357 if(!res->GetEntriesFast()) res->Expand(kNgraphs);
48797466 358
359 // geometrical efficiency
965e895b 360 TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
361 TH1 *h1[2] = {0x0, 0x0};
dad3aa6b 362 h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
363 h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
965e895b 364 Process(h1, GetGraph(0));
ea4502f6 365 delete h1[0];delete h1[1];
48797466 366
367 // tracking efficiency
dad3aa6b 368 h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
369 h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
965e895b 370 Process(h1, GetGraph(1));
ea4502f6 371 delete h1[1];
48797466 372
373 // PID efficiency
dad3aa6b 374 h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
965e895b 375 Process(h1, GetGraph(2));
ea4502f6 376 delete h1[1];
965e895b 377
378 // Refit efficiency
dad3aa6b 379 h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
965e895b 380 Process(h1, GetGraph(3));
ea4502f6 381 delete h1[1];
965e895b 382}
383
384//____________________________________________________________________
385void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
386{
ddd1a033 387// Generic function to process one reference plot
388
965e895b 389 Int_t n1 = 0, n2 = 0, ip=0;
390 Double_t eff = 0.;
391
392 TAxis *ax = h1[0]->GetXaxis();
48797466 393 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
394 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
395 n2 = (Int_t)h1[1]->GetBinContent(ib);
396 eff = n2/Float_t(n1);
397
398 ip=g->GetN();
399 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
400 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
401 }
2b104c26 402}
403
404//____________________________________________________________________
405void AliTRDcheckESD::PrintStatus(ULong_t status)
406{
ddd1a033 407// Dump track status to stdout
408
2b104c26 409 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"
410 ,Bool_t(status & AliESDtrack::kITSin)
411 ,Bool_t(status & AliESDtrack::kITSout)
412 ,Bool_t(status & AliESDtrack::kITSrefit)
413 ,Bool_t(status & AliESDtrack::kTPCin)
414 ,Bool_t(status & AliESDtrack::kTPCout)
415 ,Bool_t(status & AliESDtrack::kTPCrefit)
416 ,Bool_t(status & AliESDtrack::kTPCpid)
417 ,Bool_t(status & AliESDtrack::kTRDin)
418 ,Bool_t(status & AliESDtrack::kTRDout)
419 ,Bool_t(status & AliESDtrack::kTRDrefit)
420 ,Bool_t(status & AliESDtrack::kTRDpid)
421 ,Bool_t(status & AliESDtrack::kTRDStop)
422 ,Bool_t(status & AliESDtrack::kHMPIDout)
423 ,Bool_t(status & AliESDtrack::kHMPIDpid)
424 );
9a281e83 425}
2b104c26 426