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