]>
Commit | Line | Data |
---|---|---|
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 | ||
32 | ClassImp(AliTRDpidRefMaker) | |
33 | ClassImp(AliTRDpidRefMaker::AliTRDpidRefData) | |
34 | ClassImp(AliTRDpidRefMaker::AliTRDpidRefDataArray) | |
35 | ||
36 | //________________________________________________________________________ | |
37 | AliTRDpidRefMaker::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 | //________________________________________________________________________ | |
46 | AliTRDpidRefMaker::AliTRDpidRefDataArray::~AliTRDpidRefDataArray() | |
47 | { | |
48 | // Destructor | |
49 | delete [] fData; | |
50 | } | |
51 | ||
52 | //________________________________________________________________________ | |
53 | void 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 | //________________________________________________________________________ | |
62 | void 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 | ||
705f8b0a | 74 | //________________________________________________________________________ |
75 | AliTRDpidRefMaker::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 | ||
1ee39b3a | 95 | //________________________________________________________________________ |
96 | AliTRDpidRefMaker::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 | //________________________________________________________________________ | |
126 | AliTRDpidRefMaker::~AliTRDpidRefMaker() | |
127 | { | |
128 | if(fPIDdataArray) delete fPIDdataArray; | |
129 | if(fReconstructor) delete fReconstructor; | |
130 | } | |
131 | ||
132 | //________________________________________________________________________ | |
133 | Bool_t AliTRDpidRefMaker::CheckQuality(AliTRDseedV1* /*trklt*/) | |
134 | { | |
135 | // Place holder for checking tracklet quality for PID. | |
136 | return kTRUE; | |
137 | } | |
138 | ||
139 | ||
140 | //________________________________________________________________________ | |
141 | Float_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 | //________________________________________________________________________ | |
150 | void AliTRDpidRefMaker::ConnectInputData(Option_t *opt) | |
151 | { | |
152 | AliTRDrecoTask::ConnectInputData(opt); | |
705f8b0a | 153 | fV0s = dynamic_cast<TObjArray*>(GetInputData(2)); |
154 | fInfo = dynamic_cast<TObjArray*>(GetInputData(3)); | |
b91fdd71 | 155 | } |
1ee39b3a | 156 | |
157 | //________________________________________________________________________ | |
f8f46e4d | 158 | void AliTRDpidRefMaker::UserCreateOutputObjects() |
1ee39b3a | 159 | { |
160 | // Create histograms | |
161 | // Called once | |
162 | ||
f8f46e4d | 163 | OpenFile(1, "RECREATE"); |
1ee39b3a | 164 | fContainer = new TObjArray(); |
165 | fContainer->SetName(Form("Moni%s", GetName())); | |
166 | ||
167 | TH2 *h2 = new TH2I("hPDG","Particle abundance", AliPID::kSPECIES, -0.5, 4.5, AliTRDCalPID::kNMom, -0.5, AliTRDCalPID::kNMom-0.5); | |
168 | TAxis *ax = h2->GetXaxis(); | |
169 | ax->SetNdivisions(505); | |
170 | ax->SetTitle("Particle species"); | |
171 | for(Int_t is=AliPID::kSPECIES; is--;) ax->SetBinLabel(is+1, AliPID::ParticleShortName(is)); | |
172 | h2->GetYaxis()->SetTitle("P bins"); | |
173 | h2->GetYaxis()->SetNdivisions(511); | |
174 | ||
175 | fContainer->AddAt(h2, 0); | |
176 | ||
177 | fData = new TTree(GetName(), Form("Reference data for %s", GetName())); | |
178 | fData->Branch("s", &fPIDbin, "s/b"); | |
179 | fData->Branch("data", &fPIDdataArray); | |
180 | } | |
181 | ||
182 | //________________________________________________________________________ | |
f8f46e4d | 183 | void AliTRDpidRefMaker::UserExec(Option_t *) |
1ee39b3a | 184 | { |
185 | // Main loop | |
186 | // Called for each event | |
187 | ||
705f8b0a | 188 | AliInfo(Form("Analyse N[%d] tracks", fTracks->GetEntriesFast())); |
189 | ||
b91fdd71 | 190 | AliTRDtrackInfo *track = NULL; |
191 | AliTRDtrackV1 *trackTRD = NULL; | |
192 | AliTrackReference *ref = NULL; | |
193 | const AliTRDtrackInfo::AliESDinfo *infoESD = NULL; | |
1ee39b3a | 194 | for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){ |
195 | track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk); | |
196 | if(!track->HasESDtrack()) continue; | |
197 | trackTRD = track->GetTrack(); | |
198 | infoESD = track->GetESDinfo(); | |
199 | Double32_t *infoPID = infoESD->GetSliceIter(); | |
200 | Int_t n = infoESD->GetNSlices() - AliTRDgeometry::kNlayer; | |
201 | Double32_t *p = &infoPID[n]; | |
4226db3e | 202 | 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 | 203 | |
204 | ULong_t status = track->GetStatus(); | |
205 | if(!(status&AliESDtrack::kTPCout)) continue; | |
206 | ||
207 | // fill the pid information | |
208 | SetRefPID(fRefPID, track, fPID); | |
209 | ||
210 | // prepare PID data array | |
211 | if(!fPIDdataArray){ | |
212 | fPIDdataArray = new AliTRDpidRefDataArray(); | |
213 | } else fPIDdataArray->Reset(); | |
214 | ||
215 | // fill PID information | |
216 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){ | |
217 | ||
218 | // fill P & dE/dx information | |
219 | if(HasFriends()){ // from TRD track | |
220 | if(!trackTRD) continue; | |
b91fdd71 | 221 | AliTRDseedV1 *trackletTRD(NULL); |
1ee39b3a | 222 | trackTRD -> SetReconstructor(fReconstructor); |
223 | if(!(trackletTRD = trackTRD -> GetTracklet(ily))) continue; | |
224 | if(!CheckQuality(trackletTRD)) continue; | |
225 | if(!CookdEdx(trackletTRD)) continue; | |
226 | ||
227 | // fill momentum information | |
228 | fP = 0.; | |
229 | switch(fRefP){ | |
230 | case kMC: | |
231 | if(!(ref = track->GetTrackRef(trackletTRD))) continue; | |
232 | fP = ref->P(); | |
233 | break; | |
234 | case kRec: | |
235 | fP = trackletTRD->GetMomentum(); | |
236 | break; | |
237 | default: continue; | |
238 | } | |
239 | } else { // from ESD track | |
240 | // fill momentum information | |
241 | switch(fRefP){ | |
242 | case kMC: | |
243 | if(!(ref = track->GetTrackRef(ily))) continue; | |
244 | fP = ref->P(); | |
245 | break; | |
246 | case kRec: | |
247 | fP = p[ily]; | |
248 | break; | |
249 | default: continue; | |
250 | } | |
251 | Double32_t *it = &infoPID[ily*AliTRDCalPID::kNSlicesNN]; | |
252 | for(Int_t is=AliTRDCalPID::kNSlicesNN; is--; it++) fdEdx[is] = (*it); | |
253 | } | |
254 | ||
255 | // momentum threshold | |
256 | if(fP < fPthreshold) continue; | |
257 | ||
258 | // store information | |
259 | fPIDdataArray->PushBack(ily, AliTRDpidUtil::GetMomentumBin(fP), fdEdx); | |
260 | } | |
261 | ||
262 | Fill(); | |
263 | } | |
264 | ||
f8f46e4d | 265 | PostData(1, fContainer); |
266 | PostData(2, fData); | |
1ee39b3a | 267 | } |
268 | ||
269 | ||
270 | //________________________________________________________________________ | |
271 | void AliTRDpidRefMaker::Fill() | |
272 | { | |
273 | // Fill data tree | |
274 | ||
275 | if(!fPIDdataArray->fNtracklets) return; | |
276 | fPIDbin = TMath::LocMax(AliPID::kSPECIES, fPID); // get particle type | |
277 | // Fill data tree | |
705f8b0a | 278 | AliInfo(Form("fData[%p]", (void*)fData)); |
1ee39b3a | 279 | fData->Fill(); |
280 | ||
281 | ||
282 | // fill monitor | |
283 | for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){ | |
284 | Int_t pBin = fPIDdataArray->fData[ily].fPLbin & 0xf; | |
285 | ((TH2*)fContainer->At(0))->Fill(fPIDbin, pBin); | |
286 | } | |
287 | } | |
288 | ||
289 | //________________________________________________________________________ | |
290 | void AliTRDpidRefMaker::LinkPIDdata() | |
291 | { | |
292 | // Link data tree to data members | |
293 | fData->SetBranchAddress("s", &fPIDbin); | |
294 | fData->SetBranchAddress("data", &fPIDdataArray); | |
295 | } | |
296 | ||
297 | //________________________________________________________________________ | |
298 | void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Float_t *pid) | |
299 | { | |
300 | // Fill the reference PID values "pid" from "source" object | |
301 | // according to the option "select". Possible options are | |
302 | // - kV0 - v0 based PID | |
303 | // - kMC - MC truth [default] | |
304 | // - kRec - outside detectors | |
305 | ||
306 | ||
307 | memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t)); | |
308 | switch(select){ | |
309 | case kV0: | |
310 | { | |
311 | AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source); | |
312 | if(!track){ | |
313 | AliError("No trackInfo found"); | |
314 | return; | |
315 | } | |
316 | ||
317 | //Get V0 PID decisions from the AliTRDv0Info for all particle species (implemented so far : electrons from conversions, pions from K0s and protons from Lambdas) : | |
318 | AliTRDv0Info v0; | |
319 | for(Int_t is=AliPID::kSPECIES; is--;) fPID[is] = v0.GetV0PID(is, track); | |
320 | } | |
321 | break; | |
322 | case kMC: | |
323 | if(!HasMCdata()){ | |
324 | AliError("Could not retrive reference PID from MC"); | |
325 | return; | |
326 | } | |
327 | { | |
328 | AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source); | |
329 | switch(track->GetPDG()){ | |
330 | case kElectron: | |
331 | case kPositron: | |
332 | fPID[AliPID::kElectron] = 1.; | |
333 | break; | |
334 | case kMuonPlus: | |
335 | case kMuonMinus: | |
336 | fPID[AliPID::kMuon] = 1.; | |
337 | break; | |
338 | case kPiPlus: | |
339 | case kPiMinus: | |
340 | fPID[AliPID::kPion] = 1.; | |
341 | break; | |
342 | case kKPlus: | |
343 | case kKMinus: | |
344 | fPID[AliPID::kKaon] = 1.; | |
345 | break; | |
346 | case kProton: | |
347 | case kProtonBar: | |
348 | fPID[AliPID::kProton] = 1.; | |
349 | break; | |
350 | } | |
351 | } | |
352 | break; | |
353 | case kRec: | |
354 | { | |
355 | AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source); | |
356 | AliTRDtrackV1 *trackTRD = track->GetTrack(); | |
357 | trackTRD -> SetReconstructor(fReconstructor); | |
358 | //fReconstructor -> SetOption("nn"); | |
359 | trackTRD -> CookPID(); | |
360 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
361 | pid[iPart] = trackTRD -> GetPID(iPart); | |
362 | AliDebug(4, Form("PDG is (in V0info) %d %f", iPart, pid[iPart])); | |
363 | } | |
364 | } | |
365 | break; | |
366 | default: | |
367 | AliWarning("PID reference source not implemented"); | |
368 | return; | |
369 | } | |
370 | AliDebug(4, Form("Ref PID [%] : %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f]" | |
371 | ,AliPID::ParticleShortName(0), 1.e2*fPID[0] | |
372 | ,AliPID::ParticleShortName(1), 1.e2*fPID[1] | |
373 | ,AliPID::ParticleShortName(2), 1.e2*fPID[2] | |
374 | ,AliPID::ParticleShortName(3), 1.e2*fPID[3] | |
375 | ,AliPID::ParticleShortName(4), 1.e2*fPID[4] | |
376 | )); | |
377 | } | |
378 | ||
379 | //________________________________________________________________________ | |
380 | void AliTRDpidRefMaker::SetAbundance(Float_t train) | |
381 | { | |
382 | // Split data sample between trainning and testing | |
383 | ||
384 | if(train<0. || train >1.){ | |
385 | AliWarning("The input data should be in the interval [0, 1]"); | |
386 | return; | |
387 | } | |
388 | ||
389 | fFreq = train; | |
390 | } | |
391 |