new task to check TRD tracking performances at ESD-MC (only) level.
[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>
21#include <TFile.h>
22#include <TTree.h>
23#include <TROOT.h>
24#include <TChain.h>
25#include <TParticle.h>
26
27#include "AliLog.h"
28#include "AliAnalysisManager.h"
29#include "AliESDEvent.h"
30#include "AliMCEvent.h"
31#include "AliESDInputHandler.h"
32#include "AliMCEventHandler.h"
33
34#include "AliESDtrack.h"
35#include "AliMCParticle.h"
36#include "AliPID.h"
37#include "AliStack.h"
38#include "AliTrackReference.h"
39#include "AliTRDgeometry.h"
40
41#include "AliTRDcheckESD.h"
42
43ClassImp(AliTRDcheckESD)
44
45const Float_t AliTRDcheckESD::xTPC = 290.;
46const Float_t AliTRDcheckESD::xTOF = 365.;
47
48//____________________________________________________________________
49AliTRDcheckESD::AliTRDcheckESD():
50 AliAnalysisTask("ESDchecker", "TRD checker @ ESD level")
51 ,fStatus(0)
52 ,fESD(0x0)
53 ,fMC(0x0)
54 ,fHistos(0x0)
55{
56 //
57 // Default constructor
58 //
59
60 DefineInput(0, TChain::Class());
61 DefineOutput(0, TObjArray::Class());
62}
63
64//____________________________________________________________________
65AliTRDcheckESD::~AliTRDcheckESD()
66{
67 if(fHistos){
68 //fHistos->Delete();
69 delete fHistos;
70 }
71}
72
73//____________________________________________________________________
74void AliTRDcheckESD::ConnectInputData(Option_t *)
75{
76 //
77 // Link the Input Data
78 //
79 TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
80 if(tree) tree->SetBranchStatus("Tracks", 1);
81
82 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
83 fESD = esdH ? esdH->GetEvent() : 0x0;
84
85 if(!HasMC()) return;
86 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
87 fMC = mcH ? mcH->MCEvent() : 0x0;
88}
89
90//____________________________________________________________________
91void AliTRDcheckESD::CreateOutputObjects()
92{
93 //
94 // Create Output Containers (TObjectArray containing 1D histograms)
95 //
96 OpenFile(0, "RECREATE");
97 fHistos = new TObjArray(5);
98 //fHistos->SetOwner(kTRUE);
99
100 TH1 *h = 0x0;
101
102 // clusters per tracklet
103 if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
104 h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
105 h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
106 h->GetYaxis()->SetTitle("entries");
107 } else h->Reset();
108 fHistos->AddAt(h, kNCl);
109
110 // TPC out
111 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
112 h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., 4, -.5, 3.5);
113 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
114 h->GetYaxis()->SetTitle("status bits");
115 h->GetZaxis()->SetTitle("entries");
116 } else h->Reset();
117 fHistos->AddAt(h, kTRDstat);
118}
119
120//____________________________________________________________________
121void AliTRDcheckESD::Exec(Option_t *){
122 //
123 // Run the Analysis
124 //
125 if(!fESD){
126 AliError("ESD not found");
127 return;
128 }
129
130 // Get MC information if available
131 AliStack * fStack = 0x0;
132 if(HasMC() && !fMC){
133 AliWarning("Monte Carlo Event not available");
134 SetMC(kFALSE);
135 } else {
136 if(!(fStack = fMC->Stack())){
137 AliWarning("Cannot get the Monte Carlo Stack");
138 SetMC(kFALSE);
139 }
140 }
141 Bool_t TPCout(0), TRDin(0), TRDout(0), TRDpid(0);
142
143
144 Int_t nTRD = 0, nTPC = 0;
145 //Int_t nTracks = fESD->GetNumberOfTracks();
146 AliESDtrack *esdTrack = 0x0;
147 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
148 TPCout=0;TRDin=0;TRDout=0;TRDpid=0;
149 esdTrack = fESD->GetTrack(itrk);
150 if(esdTrack->GetNcls(1)) nTPC++;
151 if(esdTrack->GetNcls(2)) nTRD++;
152
153 // track status
154 ULong_t status = esdTrack->GetStatus();
155
156 // TRD PID
157 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
158 // pid quality
159 esdTrack->GetTRDpidQuality();
160 // kink index
161 esdTrack->GetKinkIndex(0);
162 // TPC clusters
163 esdTrack->GetNcls(1);
164
165 // look at external track param
166 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
167 Double_t xyz[3];
168 if(op){
169 op->GetXYZ(xyz);
170 op->Global2LocalPosition(xyz, op->GetAlpha());
171 //printf("op @ X[%7.3f]\n", xyz[0]);
172 }
173
174 // read MC info
175 if(!HasMC()) continue;
176
177 Int_t fLabel = esdTrack->GetLabel();
178 if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
179
180 // read MC particle
181 AliMCParticle *mcParticle = 0x0;
182 if(!(mcParticle = fMC->GetTrack(TMath::Abs(fLabel)))){
183 AliWarning(Form("MC particle missing for ESD fLabel %d.", fLabel));
184 continue;
185 }
186
187 AliTrackReference *ref = 0x0;
188 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
189 Int_t iref = 0;
190 while(iref<nRefs){
191 ref = mcParticle->GetTrackReference(iref);
192 if(ref->LocalX() > xTPC) break;
193 ref=0x0; iref++;
194 }
195
196 // read TParticle
197 TParticle *tParticle = mcParticle->Particle();
198 Int_t fPdg = tParticle->GetPdgCode();
199 //tParticle->IsPrimary();
200
201 //printf("[%c] ref[%2d]=", tParticle->IsPrimary() ? 'P' : 'S', iref);
202
203 TPCout=1;
204 if(ref){
205 if(ref->LocalX() > xTOF){
206 //printf(" TOF [");
207 ref = mcParticle->GetTrackReference(iref-1);
208 } else {
209 //printf("%7.2f [", ref->LocalX());
210 TRDin=1;
211 if(esdTrack->GetNcls(2)) TRDout=1;
212 if(esdTrack->GetTRDpidQuality()) TRDpid=1;
213 }
214 } else {
215 //printf(" TPC [");
216 ref = mcParticle->GetTrackReference(iref-1);
217 }
218 Float_t pt = ref->Pt();
219 //printf("%f]\n", pt);
220
221 TH2 *h = (TH2I*)fHistos->At(kTRDstat);
222 if(/*status & AliESDtrack::k*/TPCout) h->Fill(pt, 0);
223 if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, 1);
224 if(/*status & AliESDtrack::k*/TRDout){
225 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
226 h->Fill(pt, 2);
227 }
228 if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, 3);
229 }
230
231 PostData(0, fHistos);
232}
233
234
235//____________________________________________________________________
236void AliTRDcheckESD::Terminate(Option_t *)
237{
238}