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