1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 #include <TClonesArray.h>
18 #include <TObjArray.h>
21 #include <TGraphErrors.h>
26 #include <TParticle.h>
29 #include "AliAnalysisManager.h"
30 #include "AliESDEvent.h"
31 #include "AliESDkink.h"
32 #include "AliMCEvent.h"
33 #include "AliESDInputHandler.h"
34 #include "AliMCEventHandler.h"
36 #include "AliESDtrack.h"
37 #include "AliMCParticle.h"
40 #include "AliTrackReference.h"
42 #include "AliTRDcheckESD.h"
44 ClassImp(AliTRDcheckESD)
46 const Int_t AliTRDcheckESD::fgkNgraphs = 4;
47 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
48 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
50 //____________________________________________________________________
51 AliTRDcheckESD::AliTRDcheckESD():
52 AliAnalysisTask("ESDchecker", "TRD checker @ ESD level")
59 // Default constructor
62 DefineInput(0, TChain::Class());
63 DefineOutput(0, TObjArray::Class());
66 //____________________________________________________________________
67 AliTRDcheckESD::~AliTRDcheckESD()
75 //____________________________________________________________________
76 void AliTRDcheckESD::ConnectInputData(Option_t *)
79 // Link the Input Data
81 TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
82 if(tree) tree->SetBranchStatus("Tracks", 1);
84 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
85 fESD = esdH ? esdH->GetEvent() : 0x0;
88 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
89 fMC = mcH ? mcH->MCEvent() : 0x0;
92 //____________________________________________________________________
93 void AliTRDcheckESD::CreateOutputObjects()
96 // Create Output Containers (TObjectArray containing 1D histograms)
98 OpenFile(0, "RECREATE");
99 fHistos = new TObjArray(5);
100 //fHistos->SetOwner(kTRUE);
104 // clusters per tracklet
105 if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
106 h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
107 h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
108 h->GetYaxis()->SetTitle("entries");
110 fHistos->AddAt(h, kNCl);
112 // status bits histogram
113 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
114 h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
115 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
116 h->GetYaxis()->SetTitle("status bits");
117 h->GetZaxis()->SetTitle("entries");
119 fHistos->AddAt(h, kTRDstat);
122 TObjArray *res = new TObjArray();
123 res->SetName("Results");
124 fHistos->AddAt(res, kResults);
127 //____________________________________________________________________
128 TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t*)
130 Bool_t kBUILD = 1, // build graph if none found
131 kCLEAR = 1; // clear existing graph
133 const Char_t *name[] = {
134 "Geo", "Trk", "Pid", "Ref"
136 const Char_t *title[] = {
137 "TRD geometrical efficiency (TRDin/TPCout)"
138 ,"TRD tracking efficiency (TRDout/TRDin)"
139 ,"TRD PID efficiency (TRDpid/TRDin)"
140 ,"TRD refit efficiency (TRDrefit/TRDin)"
142 const Int_t ngr = sizeof(name)/sizeof(Char_t*);
143 if(ngr != fgkNgraphs){
144 AliWarning("No of graphs defined different from definition");
148 TObjArray *res = 0x0;
149 if(!(res = (TObjArray*)fHistos->At(kResults)) ||
150 (id < 0 || id >= ngr)){
151 AliWarning("Graph missing.");
155 TGraphErrors *g = 0x0;
156 if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
157 if(kCLEAR) for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
160 g = new TGraphErrors();
161 g->SetNameTitle(name[id], title[id]);
168 //____________________________________________________________________
169 void AliTRDcheckESD::Exec(Option_t *){
174 AliError("ESD event missing.");
178 // Get MC information if available
179 AliStack * fStack = 0x0;
181 AliWarning("MC event missing");
184 if(!(fStack = fMC->Stack())){
185 AliWarning("MC stack missing");
189 Bool_t TRDin(0), TRDout(0), TRDpid(0);
191 AliESDtrack *esdTrack = 0x0;
192 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
193 TRDin=0;TRDout=0;TRDpid=0;
194 esdTrack = fESD->GetTrack(itrk);
196 // if(esdTrack->GetNcls(1)) nTPC++;
197 // if(esdTrack->GetNcls(2)) nTRD++;
200 ULong_t status = esdTrack->GetStatus();
202 // define TPC out tracks
203 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
204 if(esdTrack->GetKinkIndex(0) > 0) continue;
207 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
209 esdTrack->GetTRDpidQuality();
211 // look at external track param
212 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
216 op->Global2LocalPosition(xyz, op->GetAlpha());
217 //printf("op @ X[%7.3f]\n", xyz[0]);
221 if(!HasMC()) continue;
223 Int_t fLabel = esdTrack->GetLabel();
224 if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
227 AliMCParticle *mcParticle = 0x0;
228 if(!(mcParticle = fMC->GetTrack(TMath::Abs(fLabel)))){
229 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
233 AliTrackReference *ref = 0x0;
234 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
236 AliWarning(Form("Track refs missing. Label[%d].", fLabel));
241 ref = mcParticle->GetTrackReference(iref);
242 if(ref->LocalX() > fgkxTPC) break;
247 //TParticle *tParticle = mcParticle->Particle();
248 //Int_t fPdg = tParticle->GetPdgCode();
249 // reject secondaries
250 //if(!tParticle->IsPrimary()) continue;
253 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
254 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
257 if(esdTrack->GetNcls(2)) TRDout=1;
258 if(esdTrack->GetTRDpidQuality()>=4) TRDpid=1;
260 } else { // track stopped in TPC
261 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
264 Float_t pt = ref->Pt();
266 TH2 *h = (TH2I*)fHistos->At(kTRDstat);
267 if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
268 if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, kTRDin);
269 if(/*status & AliESDtrack::k*/TRDout){
270 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
271 h->Fill(pt, kTRDout);
273 if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, kTRDpid);
274 if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
276 PostData(0, fHistos);
279 //____________________________________________________________________
280 Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
282 if(!TFile::Open(filename)){
283 AliWarning(Form("Couldn't open file %s.", filename));
287 if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
288 AliWarning("Missing histogram container.");
291 fHistos = (TObjArray*)o->Clone(GetName());
296 //____________________________________________________________________
297 void AliTRDcheckESD::Terminate(Option_t *)
299 TObjArray *res = 0x0;
300 if(!(res = (TObjArray*)fHistos->At(kResults))){
301 AliWarning("Graph container missing.");
304 if(!res->GetEntriesFast()) res->Expand(fgkNgraphs);
306 // geometrical efficiency
307 TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
308 TH1 *h1[2] = {0x0, 0x0};
309 h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
310 h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
311 Process(h1, GetGraph(0));
313 // tracking efficiency
314 h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
315 h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
316 Process(h1, GetGraph(1));
319 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
320 h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
321 Process(h1, GetGraph(2));
324 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
325 h1[1] = h2->ProjectionX("px1", kTRDref, kTRDref);
326 Process(h1, GetGraph(3));
329 //____________________________________________________________________
330 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
332 Int_t n1 = 0, n2 = 0, ip=0;
335 TAxis *ax = h1[0]->GetXaxis();
336 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
337 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
338 n2 = (Int_t)h1[1]->GetBinContent(ib);
339 eff = n2/Float_t(n1);
342 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
343 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);