4 #include "TEventList.h"
9 #include "AliESDtrack.h"
10 #include "AliTrackReference.h"
12 #include "AliTRDReconstructor.h"
13 #include "AliTRDtrackV1.h"
14 #include "AliTRDseedV1.h"
15 #include "AliTRDpidRefMaker.h"
16 #include "info/AliTRDv0Info.h"
17 #include "info/AliTRDpidInfo.h"
20 // Defines and implements basic functionality for building reference data for TRD PID.
22 // Here is the list of functionality provided by this class
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>
33 ClassImp(AliTRDpidRefMaker)
35 //________________________________________________________________________
36 AliTRDpidRefMaker::AliTRDpidRefMaker()
50 // Default constructor
52 SetNameTitle("PIDrefMaker", "PID Reference Maker");
55 //________________________________________________________________________
56 AliTRDpidRefMaker::AliTRDpidRefMaker(const char *name, const char *title)
57 :AliTRDrecoTask(name, title)
70 // Default constructor
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));
78 DefineInput(2, TObjArray::Class()); // v0 list
79 DefineInput(3, TObjArray::Class()); // pid info list
80 DefineOutput(2, TTree::Class());
84 //________________________________________________________________________
85 AliTRDpidRefMaker::~AliTRDpidRefMaker()
87 if(fPIDdataArray) delete fPIDdataArray;
88 if(fReconstructor) delete fReconstructor;
91 //________________________________________________________________________
92 Bool_t AliTRDpidRefMaker::CheckQuality(AliTRDseedV1* /*trklt*/)
94 // Place holder for checking tracklet quality for PID.
99 //________________________________________________________________________
100 Float_t* AliTRDpidRefMaker::CookdEdx(AliTRDseedV1 *trklt)
102 trklt->CookdEdx(AliTRDpidUtil::kNNslices);
103 memcpy(fdEdx, trklt->GetdEdx(), AliTRDpidUtil::kNNslices*sizeof(Float_t));
108 //________________________________________________________________________
109 void AliTRDpidRefMaker::UserCreateOutputObjects()
114 fContainer = new TObjArray();
115 fContainer->SetName(Form("Moni%s", GetName()));
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);
124 fContainer->AddAt(h2, 0);
127 fData = new TTree("RefPID", "Reference data for PID");
128 fData->Branch("data", &fPIDdataArray);
131 //________________________________________________________________________
132 void AliTRDpidRefMaker::UserExec(Option_t *)
135 // Called for each event
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;
141 AliDebug(1, Form("Entries: Tracks[%d] V0[%d] PID[%d]", fTracks->GetEntriesFast(), fV0s->GetEntriesFast(), fInfo->GetEntriesFast()));
142 AliTRDtrackInfo *track = NULL;
143 AliTRDtrackV1 *trackTRD = NULL;
144 AliTrackReference *ref = NULL;
145 const AliTRDtrackInfo::AliESDinfo *infoESD = NULL;
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;
154 AliWarning(Form("dEdx info missing in ESD track %d", itrk));
157 Double32_t *p = &infoPID[n];
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]));
160 ULong_t status = track->GetStatus();
161 if(!(status&AliESDtrack::kTRDpid)) continue;
163 // fill the pid information
164 SetRefPID(fRefPID, track, fPID);
166 Int_t idx(TMath::LocMax(AliPID::kSPECIES, fPID));
167 if(fPID[idx]<1.e-5) continue;
169 // prepare PID data array
171 fPIDdataArray = new AliTRDpidInfo();
172 } else fPIDdataArray->Reset();
173 fPIDdataArray->SetPID(idx);
175 // fill PID information
176 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
178 // fill P & dE/dx information
179 if(HasFriends()){ // from TRD track
180 if(!trackTRD) continue;
181 AliTRDseedV1 *trackletTRD(NULL);
182 trackTRD -> SetReconstructor(fReconstructor);
183 if(!(trackletTRD = trackTRD -> GetTracklet(ily))) continue;
184 if(!CheckQuality(trackletTRD)) continue;
185 if(!CookdEdx(trackletTRD)) continue;
187 // fill momentum information
191 if(!(ref = track->GetTrackRef(trackletTRD))) continue;
195 fP = trackletTRD->GetMomentum();
199 } else { // from ESD track
200 // fill momentum information
203 if(!(ref = track->GetTrackRef(ily))) continue;
211 Double32_t *it = &infoPID[ily*AliTRDCalPID::kNSlicesNN];
212 for(Int_t is=AliTRDCalPID::kNSlicesNN; is--; it++) fdEdx[is] = (*it);
215 // momentum threshold
216 if(fP < fPthreshold) continue;
219 fPIDdataArray->PushBack(ily, AliTRDpidUtil::GetMomentumBin(fP), fdEdx);
225 PostData(1, fContainer);
230 //________________________________________________________________________
231 void AliTRDpidRefMaker::Fill()
235 if(!fPIDdataArray->GetNtracklets()) return;
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);
247 //________________________________________________________________________
248 void AliTRDpidRefMaker::LinkPIDdata()
250 // Link data tree to data members
251 fData->SetBranchAddress("data", &fPIDdataArray);
254 //________________________________________________________________________
255 void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, AliTRDtrackInfo *track, Float_t *pid)
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
264 AliError("No trackInfo found");
267 memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
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) :
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);
283 AliError("Could not retrive reference PID from MC");
286 switch(track->GetPDG()){
289 fPID[AliPID::kElectron] = 1.;
293 fPID[AliPID::kMuon] = 1.;
297 fPID[AliPID::kPion] = 1.;
301 fPID[AliPID::kKaon] = 1.;
305 fPID[AliPID::kProton] = 1.;
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]));
322 AliWarning("PID reference source not implemented");
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]
334 //________________________________________________________________________
335 void AliTRDpidRefMaker::SetAbundance(Float_t train)
337 // Split data sample between trainning and testing
339 if(train<0. || train >1.){
340 AliWarning("The input data should be in the interval [0, 1]");