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