Bug in raw histo filling corrected
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDpidRefMaker.cxx
CommitLineData
1ee39b3a 1#include "TPDGCode.h"
2#include "TFile.h"
3#include "TTree.h"
4#include "TEventList.h"
5#include "TObjArray.h"
6#include "TH2.h"
7
8#include "AliLog.h"
9#include "AliESDtrack.h"
10#include "AliTrackReference.h"
11
12#include "AliTRDReconstructor.h"
13#include "AliTRDtrackV1.h"
14#include "AliTRDseedV1.h"
15#include "AliTRDpidRefMaker.h"
16#include "info/AliTRDv0Info.h"
17
18
19// Defines and implements basic functionality for building reference data for TRD PID.
20//
21// Here is the list of functionality provided by this class
22// 1.
23// 2.
24//
25// Authors:
26// Alex Bercuci <A.Bercuci@gsi.de>
27// Alex Wilk <wilka@uni-muenster.de>
28// Markus Fasel <mfasel@gsi.de>
29// Markus Heide <mheide@uni-muenster.de>
30//
31
32ClassImp(AliTRDpidRefMaker)
33ClassImp(AliTRDpidRefMaker::AliTRDpidRefData)
34ClassImp(AliTRDpidRefMaker::AliTRDpidRefDataArray)
35
36//________________________________________________________________________
37AliTRDpidRefMaker::AliTRDpidRefDataArray::AliTRDpidRefDataArray() :
38 fNtracklets(0)
b91fdd71 39 ,fData(NULL)
1ee39b3a 40{
41 // Constructor of data array
42 fData = new AliTRDpidRefData[AliTRDgeometry::kNlayer];
43}
44
45//________________________________________________________________________
46AliTRDpidRefMaker::AliTRDpidRefDataArray::~AliTRDpidRefDataArray()
47{
48 // Destructor
49 delete [] fData;
50}
51
52//________________________________________________________________________
53void AliTRDpidRefMaker::AliTRDpidRefDataArray::PushBack(Int_t ly, Int_t p, Float_t *dedx)
54{
55// Add PID data to the end of the array
56 fData[fNtracklets].fPLbin= (ly<<4) | (p&0xf);
57 memcpy(fData[fNtracklets].fdEdx, dedx, 8*sizeof(Float_t));
58 fNtracklets++;
59}
60
61//________________________________________________________________________
62void AliTRDpidRefMaker::AliTRDpidRefDataArray::Reset()
63{
64// Reset content
65
66 if(!fNtracklets) return;
67 while(fNtracklets--){
68 fData[fNtracklets].fPLbin = 0xff;
69 memset(fData[fNtracklets].fdEdx, 0, 8*sizeof(Float_t));
70 }
71 fNtracklets=0;
72}
73
74//________________________________________________________________________
705f8b0a 75AliTRDpidRefMaker::AliTRDpidRefMaker()
76 :AliTRDrecoTask()
77 ,fReconstructor(NULL)
78 ,fV0s(NULL)
79 ,fData(NULL)
80 ,fInfo(NULL)
81 ,fPIDdataArray(NULL)
82 ,fRefPID(kMC)
83 ,fRefP(kMC)
84 ,fPIDbin(0xff)
85 ,fFreq(1.)
86 ,fP(-1.)
87 ,fPthreshold(0.)
88{
89 //
90 // Default constructor
91 //
92 SetNameTitle("PIDrefMaker", "PID Reference Maker");
93}
94
95//________________________________________________________________________
1ee39b3a 96AliTRDpidRefMaker::AliTRDpidRefMaker(const char *name, const char *title)
97 :AliTRDrecoTask(name, title)
b91fdd71 98 ,fReconstructor(NULL)
99 ,fV0s(NULL)
100 ,fData(NULL)
101 ,fInfo(NULL)
102 ,fPIDdataArray(NULL)
1ee39b3a 103 ,fRefPID(kMC)
104 ,fRefP(kMC)
105 ,fPIDbin(0xff)
106 ,fFreq(1.)
107 ,fP(-1.)
108 ,fPthreshold(0.5)
109{
110 //
111 // Default constructor
112 //
113
114 fReconstructor = new AliTRDReconstructor();
115 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
116 memset(fdEdx, 0, 10*sizeof(Float_t));
117 memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
118
705f8b0a 119 DefineInput(2, TObjArray::Class()); // v0 list
120 DefineInput(3, TObjArray::Class()); // pid info list
121 DefineOutput(2, TTree::Class());
1ee39b3a 122}
123
124
125//________________________________________________________________________
126AliTRDpidRefMaker::~AliTRDpidRefMaker()
127{
128 if(fPIDdataArray) delete fPIDdataArray;
129 if(fReconstructor) delete fReconstructor;
130}
131
132//________________________________________________________________________
133Bool_t AliTRDpidRefMaker::CheckQuality(AliTRDseedV1* /*trklt*/)
134{
135// Place holder for checking tracklet quality for PID.
136 return kTRUE;
137}
138
139
140//________________________________________________________________________
141Float_t* AliTRDpidRefMaker::CookdEdx(AliTRDseedV1 *trklt)
142{
143 trklt->CookdEdx(AliTRDpidUtil::kNNslices);
144 memcpy(fdEdx, trklt->GetdEdx(), AliTRDpidUtil::kNNslices*sizeof(Float_t));
145 return fdEdx;
146}
147
148
149//________________________________________________________________________
f8f46e4d 150void AliTRDpidRefMaker::UserCreateOutputObjects()
1ee39b3a 151{
152 // Create histograms
153 // Called once
154
1ee39b3a 155 fContainer = new TObjArray();
156 fContainer->SetName(Form("Moni%s", GetName()));
157
158 TH2 *h2 = new TH2I("hPDG","Particle abundance", AliPID::kSPECIES, -0.5, 4.5, AliTRDCalPID::kNMom, -0.5, AliTRDCalPID::kNMom-0.5);
159 TAxis *ax = h2->GetXaxis();
160 ax->SetNdivisions(505);
161 ax->SetTitle("Particle species");
162 for(Int_t is=AliPID::kSPECIES; is--;) ax->SetBinLabel(is+1, AliPID::ParticleShortName(is));
163 h2->GetYaxis()->SetTitle("P bins");
164 h2->GetYaxis()->SetNdivisions(511);
165
166 fContainer->AddAt(h2, 0);
167
168 fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
169 fData->Branch("s", &fPIDbin, "s/b");
170 fData->Branch("data", &fPIDdataArray);
171}
172
173//________________________________________________________________________
f8f46e4d 174void AliTRDpidRefMaker::UserExec(Option_t *)
1ee39b3a 175{
176 // Main loop
177 // Called for each event
178
87c9fc00 179 if(!(fTracks = dynamic_cast<TObjArray*>(GetInputData(1)))) return;
180 if(!(fV0s = dynamic_cast<TObjArray*>(GetInputData(2)))) return;
181 if(!(fInfo = dynamic_cast<TObjArray*>(GetInputData(3)))) return;
705f8b0a 182
87c9fc00 183 AliDebug(1, Form("Analyse N[%d] tracks", fTracks->GetEntriesFast()));
b91fdd71 184 AliTRDtrackInfo *track = NULL;
185 AliTRDtrackV1 *trackTRD = NULL;
186 AliTrackReference *ref = NULL;
187 const AliTRDtrackInfo::AliESDinfo *infoESD = NULL;
1ee39b3a 188 for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
189 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
190 if(!track->HasESDtrack()) continue;
191 trackTRD = track->GetTrack();
192 infoESD = track->GetESDinfo();
193 Double32_t *infoPID = infoESD->GetSliceIter();
194 Int_t n = infoESD->GetNSlices() - AliTRDgeometry::kNlayer;
46760dda 195 if(n==0){
196 AliWarning(Form("dEdx info missing in ESD track %d", itrk));
197 continue;
198 }
1ee39b3a 199 Double32_t *p = &infoPID[n];
4226db3e 200 AliDebug(4, Form("n[%d] p[GeV/c]{%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f}", n, p[0], p[1], p[2], p[3], p[4], p[5]));
1ee39b3a 201
202 ULong_t status = track->GetStatus();
203 if(!(status&AliESDtrack::kTPCout)) continue;
204
205 // fill the pid information
206 SetRefPID(fRefPID, track, fPID);
207
208 // prepare PID data array
209 if(!fPIDdataArray){
210 fPIDdataArray = new AliTRDpidRefDataArray();
211 } else fPIDdataArray->Reset();
212
213 // fill PID information
214 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
215
216 // fill P & dE/dx information
217 if(HasFriends()){ // from TRD track
218 if(!trackTRD) continue;
b91fdd71 219 AliTRDseedV1 *trackletTRD(NULL);
1ee39b3a 220 trackTRD -> SetReconstructor(fReconstructor);
221 if(!(trackletTRD = trackTRD -> GetTracklet(ily))) continue;
222 if(!CheckQuality(trackletTRD)) continue;
223 if(!CookdEdx(trackletTRD)) continue;
224
225 // fill momentum information
226 fP = 0.;
227 switch(fRefP){
228 case kMC:
229 if(!(ref = track->GetTrackRef(trackletTRD))) continue;
230 fP = ref->P();
231 break;
232 case kRec:
233 fP = trackletTRD->GetMomentum();
234 break;
235 default: continue;
236 }
237 } else { // from ESD track
238 // fill momentum information
239 switch(fRefP){
240 case kMC:
241 if(!(ref = track->GetTrackRef(ily))) continue;
242 fP = ref->P();
243 break;
244 case kRec:
245 fP = p[ily];
246 break;
247 default: continue;
248 }
249 Double32_t *it = &infoPID[ily*AliTRDCalPID::kNSlicesNN];
250 for(Int_t is=AliTRDCalPID::kNSlicesNN; is--; it++) fdEdx[is] = (*it);
251 }
252
253 // momentum threshold
254 if(fP < fPthreshold) continue;
255
256 // store information
257 fPIDdataArray->PushBack(ily, AliTRDpidUtil::GetMomentumBin(fP), fdEdx);
258 }
259
260 Fill();
261 }
262
f8f46e4d 263 PostData(1, fContainer);
264 PostData(2, fData);
1ee39b3a 265}
266
267
268//________________________________________________________________________
269void AliTRDpidRefMaker::Fill()
270{
271// Fill data tree
272
273 if(!fPIDdataArray->fNtracklets) return;
274 fPIDbin = TMath::LocMax(AliPID::kSPECIES, fPID); // get particle type
87c9fc00 275 // Fill data tree
1ee39b3a 276 fData->Fill();
277
278
279 // fill monitor
280 for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
281 Int_t pBin = fPIDdataArray->fData[ily].fPLbin & 0xf;
282 ((TH2*)fContainer->At(0))->Fill(fPIDbin, pBin);
283 }
284}
285
286//________________________________________________________________________
287void AliTRDpidRefMaker::LinkPIDdata()
288{
289// Link data tree to data members
290 fData->SetBranchAddress("s", &fPIDbin);
291 fData->SetBranchAddress("data", &fPIDdataArray);
292}
293
294//________________________________________________________________________
295void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Float_t *pid)
296{
297// Fill the reference PID values "pid" from "source" object
298// according to the option "select". Possible options are
299// - kV0 - v0 based PID
300// - kMC - MC truth [default]
301// - kRec - outside detectors
302
303
304 memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
305 switch(select){
306 case kV0:
307 {
308 AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
309 if(!track){
310 AliError("No trackInfo found");
311 return;
312 }
313
314 //Get V0 PID decisions from the AliTRDv0Info for all particle species (implemented so far : electrons from conversions, pions from K0s and protons from Lambdas) :
315 AliTRDv0Info v0;
316 for(Int_t is=AliPID::kSPECIES; is--;) fPID[is] = v0.GetV0PID(is, track);
317 }
318 break;
319 case kMC:
320 if(!HasMCdata()){
321 AliError("Could not retrive reference PID from MC");
322 return;
323 }
324 {
325 AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
326 switch(track->GetPDG()){
327 case kElectron:
328 case kPositron:
329 fPID[AliPID::kElectron] = 1.;
330 break;
331 case kMuonPlus:
332 case kMuonMinus:
333 fPID[AliPID::kMuon] = 1.;
334 break;
335 case kPiPlus:
336 case kPiMinus:
337 fPID[AliPID::kPion] = 1.;
338 break;
339 case kKPlus:
340 case kKMinus:
341 fPID[AliPID::kKaon] = 1.;
342 break;
343 case kProton:
344 case kProtonBar:
345 fPID[AliPID::kProton] = 1.;
346 break;
347 }
348 }
349 break;
350 case kRec:
351 {
352 AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
353 AliTRDtrackV1 *trackTRD = track->GetTrack();
354 trackTRD -> SetReconstructor(fReconstructor);
355 //fReconstructor -> SetOption("nn");
356 trackTRD -> CookPID();
357 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
358 pid[iPart] = trackTRD -> GetPID(iPart);
359 AliDebug(4, Form("PDG is (in V0info) %d %f", iPart, pid[iPart]));
360 }
361 }
362 break;
363 default:
364 AliWarning("PID reference source not implemented");
365 return;
366 }
367 AliDebug(4, Form("Ref PID [%] : %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f]"
368 ,AliPID::ParticleShortName(0), 1.e2*fPID[0]
369 ,AliPID::ParticleShortName(1), 1.e2*fPID[1]
370 ,AliPID::ParticleShortName(2), 1.e2*fPID[2]
371 ,AliPID::ParticleShortName(3), 1.e2*fPID[3]
372 ,AliPID::ParticleShortName(4), 1.e2*fPID[4]
373 ));
374}
375
376//________________________________________________________________________
377void AliTRDpidRefMaker::SetAbundance(Float_t train)
378{
379// Split data sample between trainning and testing
380
381 if(train<0. || train >1.){
382 AliWarning("The input data should be in the interval [0, 1]");
383 return;
384 }
385
386 fFreq = train;
387}
388