fix the error calculations in generating function cumulants when weights are used
[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
16
17#include <TClonesArray.h>
18#include <TObjArray.h>
19#include <TObject.h>
20#include <TH2I.h>
48797466 21#include <TGraphErrors.h>
9a281e83 22#include <TFile.h>
23#include <TTree.h>
24#include <TROOT.h>
25#include <TChain.h>
26#include <TParticle.h>
27
28#include "AliLog.h"
29#include "AliAnalysisManager.h"
30#include "AliESDEvent.h"
d12237d6 31#include "AliESDkink.h"
9a281e83 32#include "AliMCEvent.h"
33#include "AliESDInputHandler.h"
34#include "AliMCEventHandler.h"
35
36#include "AliESDtrack.h"
37#include "AliMCParticle.h"
38#include "AliPID.h"
39#include "AliStack.h"
40#include "AliTrackReference.h"
41#include "AliTRDgeometry.h"
42
43#include "AliTRDcheckESD.h"
44
45ClassImp(AliTRDcheckESD)
46
47const Float_t AliTRDcheckESD::xTPC = 290.;
48const Float_t AliTRDcheckESD::xTOF = 365.;
49
50//____________________________________________________________________
51AliTRDcheckESD::AliTRDcheckESD():
52 AliAnalysisTask("ESDchecker", "TRD checker @ ESD level")
53 ,fStatus(0)
54 ,fESD(0x0)
55 ,fMC(0x0)
56 ,fHistos(0x0)
57{
58 //
59 // Default constructor
60 //
61
62 DefineInput(0, TChain::Class());
63 DefineOutput(0, TObjArray::Class());
64}
65
66//____________________________________________________________________
67AliTRDcheckESD::~AliTRDcheckESD()
68{
69 if(fHistos){
70 //fHistos->Delete();
71 delete fHistos;
72 }
73}
74
75//____________________________________________________________________
76void AliTRDcheckESD::ConnectInputData(Option_t *)
77{
78 //
79 // Link the Input Data
80 //
81 TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
82 if(tree) tree->SetBranchStatus("Tracks", 1);
83
84 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
85 fESD = esdH ? esdH->GetEvent() : 0x0;
86
87 if(!HasMC()) return;
88 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
89 fMC = mcH ? mcH->MCEvent() : 0x0;
90}
91
92//____________________________________________________________________
93void AliTRDcheckESD::CreateOutputObjects()
94{
95 //
96 // Create Output Containers (TObjectArray containing 1D histograms)
97 //
98 OpenFile(0, "RECREATE");
99 fHistos = new TObjArray(5);
100 //fHistos->SetOwner(kTRUE);
101
102 TH1 *h = 0x0;
103
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");
109 } else h->Reset();
110 fHistos->AddAt(h, kNCl);
111
d12237d6 112 // status bits histogram
9a281e83 113 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
d12237d6 114 h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
9a281e83 115 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
116 h->GetYaxis()->SetTitle("status bits");
117 h->GetZaxis()->SetTitle("entries");
118 } else h->Reset();
119 fHistos->AddAt(h, kTRDstat);
48797466 120
d12237d6 121 // results array
48797466 122 TObjArray *res = new TObjArray();
123 res->SetName("Results");
124 fHistos->AddAt(res, kResults);
9a281e83 125}
126
127//____________________________________________________________________
128void AliTRDcheckESD::Exec(Option_t *){
129 //
130 // Run the Analysis
131 //
132 if(!fESD){
48797466 133 AliError("ESD event missing.");
9a281e83 134 return;
135 }
136
137 // Get MC information if available
138 AliStack * fStack = 0x0;
139 if(HasMC() && !fMC){
48797466 140 AliWarning("MC event missing");
9a281e83 141 SetMC(kFALSE);
142 } else {
143 if(!(fStack = fMC->Stack())){
48797466 144 AliWarning("MC stack missing");
9a281e83 145 SetMC(kFALSE);
146 }
147 }
d12237d6 148 Bool_t TRDin(0), TRDout(0), TRDpid(0);
9a281e83 149
9a281e83 150 AliESDtrack *esdTrack = 0x0;
151 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
d12237d6 152 TRDin=0;TRDout=0;TRDpid=0;
9a281e83 153 esdTrack = fESD->GetTrack(itrk);
d12237d6 154
155// if(esdTrack->GetNcls(1)) nTPC++;
156// if(esdTrack->GetNcls(2)) nTRD++;
9a281e83 157
158 // track status
d12237d6 159 ULong_t status = esdTrack->GetStatus();
160
161 // define TPC out tracks
162 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
163 if(esdTrack->GetKinkIndex(0) > 0) continue;
9a281e83 164
165 // TRD PID
166 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
167 // pid quality
168 esdTrack->GetTRDpidQuality();
9a281e83 169
170 // look at external track param
171 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
172 Double_t xyz[3];
173 if(op){
174 op->GetXYZ(xyz);
175 op->Global2LocalPosition(xyz, op->GetAlpha());
176 //printf("op @ X[%7.3f]\n", xyz[0]);
177 }
178
179 // read MC info
180 if(!HasMC()) continue;
181
182 Int_t fLabel = esdTrack->GetLabel();
183 if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
184
185 // read MC particle
186 AliMCParticle *mcParticle = 0x0;
187 if(!(mcParticle = fMC->GetTrack(TMath::Abs(fLabel)))){
48797466 188 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
9a281e83 189 continue;
190 }
191
192 AliTrackReference *ref = 0x0;
193 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
48797466 194 if(!nRefs){
195 AliWarning(Form("Track refs missing. Label[%d].", fLabel));
196 continue;
197 }
9a281e83 198 Int_t iref = 0;
199 while(iref<nRefs){
200 ref = mcParticle->GetTrackReference(iref);
201 if(ref->LocalX() > xTPC) break;
202 ref=0x0; iref++;
203 }
204
205 // read TParticle
206 TParticle *tParticle = mcParticle->Particle();
48797466 207 //Int_t fPdg = tParticle->GetPdgCode();
d12237d6 208 // reject secondaries
209 //if(!tParticle->IsPrimary()) continue;
9a281e83 210
d12237d6 211 if(ref){
212 if(ref->LocalX() > xTOF){ // track skipping TRD fiducial volume
48797466 213 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
9a281e83 214 } else {
d12237d6 215 TRDin=1;
216 if(esdTrack->GetNcls(2)) TRDout=1;
217 if(esdTrack->GetTRDpidQuality()) TRDpid=1;
9a281e83 218 }
d12237d6 219 } else { // track stopped in TPC
48797466 220 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
221 }
d12237d6 222 // get the MC pt !!
9a281e83 223 Float_t pt = ref->Pt();
d12237d6 224
9a281e83 225 TH2 *h = (TH2I*)fHistos->At(kTRDstat);
d12237d6 226 if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
227 if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, kTRDin);
9a281e83 228 if(/*status & AliESDtrack::k*/TRDout){
229 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
d12237d6 230 h->Fill(pt, kTRDout);
9a281e83 231 }
d12237d6 232 if(status & AliESDtrack::kTRDpid) h->Fill(pt, kTRDpid);
9a281e83 233 }
9a281e83 234 PostData(0, fHistos);
235}
236
237
238//____________________________________________________________________
239void AliTRDcheckESD::Terminate(Option_t *)
240{
48797466 241 TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
242 TAxis *ax = h2->GetXaxis();
243 TObjArray *res = (TObjArray*)fHistos->At(kResults);
244
245 TH1 *h1[2] = {0x0, 0x0};
246 TGraphErrors *g = 0x0;
d12237d6 247 res->Expand(3);
48797466 248 Int_t n1 = 0, n2 = 0, ip=0;
249 Double_t eff = 0.;
250
251 // geometrical efficiency
d12237d6 252 h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
253 h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
48797466 254 res->Add(g = new TGraphErrors());
255 g->SetNameTitle("geom", "TRD geometrical efficiency (TRDin/TPCout)");
256 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
257 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
258 n2 = (Int_t)h1[1]->GetBinContent(ib);
259 eff = n2/Float_t(n1);
260
261 ip=g->GetN();
262 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
263 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
264 }
265
266 // tracking efficiency
d12237d6 267 h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
268 h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
48797466 269 res->Add(g = new TGraphErrors());
270 g->SetNameTitle("tracking", "TRD tracking efficiency (TRDout/TRDin)");
271 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
272 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
273 n2 = (Int_t)h1[1]->GetBinContent(ib);
274 eff = n2/Float_t(n1);
275
276 ip=g->GetN();
277 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
278 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
279 }
280
281 // PID efficiency
d12237d6 282 h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
283 h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
48797466 284 res->Add(g = new TGraphErrors());
285 g->SetNameTitle("PID", "TRD PID efficiency (TRDpid/TRDin)");
286 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
287 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
288 n2 = (Int_t)h1[1]->GetBinContent(ib);
289 eff = n2/Float_t(n1);
290
291 ip=g->GetN();
292 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
293 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
294 }
9a281e83 295}
d12237d6 296
297// if(TRDin && !TRDout){
298// printf("[%c] x[%2d]=%7.2f pt[%7.3f] id[%d] kink[%d] label[%d] in[%d] out[%d] pdg[%d]\n", tParticle->IsPrimary() ? 'P' : 'S', iref, ref->LocalX(), pt, esdTrack->GetID(), esdTrack->GetKinkIndex(0), fLabel, TRDin, TRDout, tParticle->GetPdgCode());
299// printf("\tstatus ITS[%d %d] TPC[%d %d %d] \n"
300// ,Bool_t(status & AliESDtrack::kITSin)
301// ,Bool_t(status & AliESDtrack::kITSout)
302// ,Bool_t(status & AliESDtrack::kTPCin)
303// ,Bool_t(status & AliESDtrack::kTPCout)
304// ,Bool_t(status & AliESDtrack::kTPCrefit)
305// );
306// printf("clusters ITS[%d] TPC[%d] TRD[%d]\n", esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2));
307// if(esdTrack->GetKinkIndex(0) != 0){
308// AliESDkink *kink = fESD->GetKink(TMath::Abs(esdTrack->GetKinkIndex(0)));
309// if(kink) printf("index[%d %d] label[%d %d]\n", kink->GetIndex(0), kink->GetIndex(1), kink->GetLabel(0), kink->GetLabel(1));
310// }
311// }