29bf19f2 |
1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
3 | * * |
4 | * Author: The ALICE Off-line Project. * |
5 | * Contributors are mentioned in the code where appropriate. * |
6 | * * |
7 | * Permission to use, copy, modify and distribute this software and its * |
8 | * documentation strictly for non-commercial purposes is hereby granted * |
9 | * without fee, provided that the above copyright notice appears in all * |
10 | * copies and that both the copyright notice and this permission notice * |
11 | * appear in the supporting documentation. The authors make no claims * |
12 | * about the suitability of this software for any purpose. It is * |
13 | * provided "as is" without express or implied warranty. * |
14 | **************************************************************************/ |
15 | |
4ec8e76d |
16 | /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */ |
29bf19f2 |
17 | |
18 | //----------------------------------------------------------------- |
4ec8e76d |
19 | // Base class for handling the pid response // |
20 | // functions of all detectors // |
21 | // and give access to the nsigmas // |
22 | // // |
23 | // Origin: Jens Wiechula, Uni Tuebingen, jens.wiechula@cern.ch // |
29bf19f2 |
24 | //----------------------------------------------------------------- |
25 | |
4ec8e76d |
26 | #include <TList.h> |
27 | #include <TObjArray.h> |
28 | #include <TPRegexp.h> |
29 | #include <TF1.h> |
30 | #include <TSpline.h> |
31 | #include <TFile.h> |
32 | |
33 | #include <AliVEvent.h> |
fd21ec8d |
34 | #include <AliVTrack.h> |
4ec8e76d |
35 | #include <AliLog.h> |
36 | #include <AliPID.h> |
ea235c90 |
37 | #include <AliOADBContainer.h> |
38 | #include <AliTRDPIDReference.h> |
29bf19f2 |
39 | |
40 | #include "AliPIDResponse.h" |
41 | |
42 | ClassImp(AliPIDResponse); |
43 | |
4ec8e76d |
44 | AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) : |
45 | TNamed("PIDResponse","PIDResponse"), |
46 | fITSResponse(isMC), |
47 | fTPCResponse(), |
48 | fTRDResponse(), |
49 | fTOFResponse(), |
fd21ec8d |
50 | fRange(5.), |
51 | fITSPIDmethod(kITSTruncMean), |
4ec8e76d |
52 | fIsMC(isMC), |
53 | fOADBPath(), |
54 | fBeamType("PP"), |
55 | fLHCperiod(), |
56 | fMCperiodTPC(), |
fd21ec8d |
57 | fMCperiodUser(), |
ea235c90 |
58 | fCurrentFile(), |
4ec8e76d |
59 | fRecoPass(0), |
fd21ec8d |
60 | fRecoPassUser(-1), |
4ec8e76d |
61 | fRun(0), |
62 | fOldRun(0), |
63 | fArrPidResponseMaster(0x0), |
64 | fResolutionCorrection(0x0), |
ea235c90 |
65 | fTRDPIDParams(0x0), |
66 | fTRDPIDReference(0x0), |
4ec8e76d |
67 | fTOFTimeZeroType(kBest_T0), |
68 | fTOFres(100.) |
69 | { |
70 | // |
71 | // default ctor |
72 | // |
73 | AliLog::SetClassDebugLevel("AliPIDResponse",10); |
74 | AliLog::SetClassDebugLevel("AliESDpid",10); |
75 | AliLog::SetClassDebugLevel("AliAODpidUtil",10); |
ea235c90 |
76 | |
77 | memset(fTRDslicesForPID,0,sizeof(UInt_t)*2); |
4ec8e76d |
78 | } |
79 | |
80 | //______________________________________________________________________________ |
81 | AliPIDResponse::~AliPIDResponse() |
82 | { |
83 | // |
84 | // dtor |
85 | // |
86 | delete fArrPidResponseMaster; |
ea235c90 |
87 | delete fTRDPIDParams; |
88 | delete fTRDPIDReference; |
4ec8e76d |
89 | } |
90 | |
91 | //______________________________________________________________________________ |
92 | AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) : |
93 | TNamed(other), |
94 | fITSResponse(other.fITSResponse), |
95 | fTPCResponse(other.fTPCResponse), |
96 | fTRDResponse(other.fTRDResponse), |
97 | fTOFResponse(other.fTOFResponse), |
fd21ec8d |
98 | fRange(other.fRange), |
99 | fITSPIDmethod(other.fITSPIDmethod), |
4ec8e76d |
100 | fIsMC(other.fIsMC), |
101 | fOADBPath(other.fOADBPath), |
102 | fBeamType("PP"), |
103 | fLHCperiod(), |
104 | fMCperiodTPC(), |
fd21ec8d |
105 | fMCperiodUser(other.fMCperiodUser), |
ea235c90 |
106 | fCurrentFile(), |
4ec8e76d |
107 | fRecoPass(0), |
fd21ec8d |
108 | fRecoPassUser(other.fRecoPassUser), |
4ec8e76d |
109 | fRun(0), |
110 | fOldRun(0), |
111 | fArrPidResponseMaster(0x0), |
112 | fResolutionCorrection(0x0), |
ea235c90 |
113 | fTRDPIDParams(0x0), |
114 | fTRDPIDReference(0x0), |
4ec8e76d |
115 | fTOFTimeZeroType(AliPIDResponse::kBest_T0), |
116 | fTOFres(100.) |
117 | { |
118 | // |
119 | // copy ctor |
120 | // |
ea235c90 |
121 | memset(fTRDslicesForPID,0,sizeof(UInt_t)*2); |
4ec8e76d |
122 | } |
123 | |
124 | //______________________________________________________________________________ |
125 | AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other) |
126 | { |
127 | // |
128 | // copy ctor |
129 | // |
130 | if(this!=&other) { |
131 | delete fArrPidResponseMaster; |
132 | TNamed::operator=(other); |
133 | fITSResponse=other.fITSResponse; |
134 | fTPCResponse=other.fTPCResponse; |
135 | fTRDResponse=other.fTRDResponse; |
136 | fTOFResponse=other.fTOFResponse; |
fd21ec8d |
137 | fRange=other.fRange; |
138 | fITSPIDmethod=other.fITSPIDmethod; |
4ec8e76d |
139 | fOADBPath=other.fOADBPath; |
140 | fIsMC=other.fIsMC; |
141 | fBeamType="PP"; |
142 | fLHCperiod=""; |
143 | fMCperiodTPC=""; |
fd21ec8d |
144 | fMCperiodUser=other.fMCperiodUser; |
ea235c90 |
145 | fCurrentFile=""; |
4ec8e76d |
146 | fRecoPass=0; |
fd21ec8d |
147 | fRecoPassUser=other.fRecoPassUser; |
4ec8e76d |
148 | fRun=0; |
149 | fOldRun=0; |
150 | fArrPidResponseMaster=0x0; |
151 | fResolutionCorrection=0x0; |
ea235c90 |
152 | fTRDPIDParams=0x0; |
153 | fTRDPIDReference=0x0; |
154 | memset(fTRDslicesForPID,0,sizeof(UInt_t)*2); |
4ec8e76d |
155 | fTOFTimeZeroType=AliPIDResponse::kBest_T0; |
156 | fTOFres=100.; |
157 | } |
158 | return *this; |
159 | } |
160 | |
fd21ec8d |
161 | //______________________________________________________________________________ |
162 | Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const |
163 | { |
164 | // |
165 | // NumberOfSigmas for 'detCode' |
166 | // |
167 | |
168 | switch (detCode){ |
169 | case kDetITS: return NumberOfSigmasITS(track, type); break; |
170 | case kDetTPC: return NumberOfSigmasTPC(track, type); break; |
171 | case kDetTOF: return NumberOfSigmasTOF(track, type); break; |
172 | // case kDetTRD: return ComputeTRDProbability(track, type); break; |
173 | // case kDetPHOS: return ComputePHOSProbability(track, type); break; |
174 | // case kDetEMCAL: return ComputeEMCALProbability(track, type); break; |
175 | // case kDetHMPID: return ComputeHMPIDProbability(track, type); break; |
176 | default: return -999.; |
177 | } |
178 | |
179 | } |
180 | |
181 | //______________________________________________________________________________ |
182 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const |
183 | { |
184 | // |
185 | // Compute PID response of 'detCode' |
186 | // |
187 | |
188 | switch (detCode){ |
189 | case kDetITS: return ComputeITSProbability(track, nSpecies, p); break; |
190 | case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break; |
191 | case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break; |
192 | case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break; |
193 | case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break; |
194 | case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break; |
195 | case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break; |
196 | default: return kDetNoSignal; |
197 | } |
198 | } |
199 | |
200 | //______________________________________________________________________________ |
201 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const |
202 | { |
203 | // |
204 | // Compute PID response for the ITS |
205 | // |
206 | |
207 | // set flat distribution (no decision) |
208 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
209 | |
210 | if ((track->GetStatus()&AliVTrack::kITSin)==0 && |
211 | (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal; |
212 | |
213 | Double_t mom=track->P(); |
214 | Double_t dedx=track->GetITSsignal(); |
215 | Bool_t isSA=kTRUE; |
216 | Double_t momITS=mom; |
217 | ULong_t trStatus=track->GetStatus(); |
218 | if(trStatus&AliVTrack::kTPCin) isSA=kFALSE; |
219 | UChar_t clumap=track->GetITSClusterMap(); |
220 | Int_t nPointsForPid=0; |
221 | for(Int_t i=2; i<6; i++){ |
222 | if(clumap&(1<<i)) ++nPointsForPid; |
223 | } |
224 | |
225 | if(nPointsForPid<3) { // track not to be used for combined PID purposes |
226 | // track->ResetStatus(AliVTrack::kITSpid); |
227 | return kDetNoSignal; |
228 | } |
229 | |
230 | Bool_t mismatch=kTRUE/*, heavy=kTRUE*/; |
231 | for (Int_t j=0; j<AliPID::kSPECIES; j++) { |
232 | Double_t mass=AliPID::ParticleMass(j);//GeV/c^2 |
233 | Double_t bethe=fITSResponse.Bethe(momITS,mass); |
234 | Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA); |
235 | if (TMath::Abs(dedx-bethe) > fRange*sigma) { |
236 | p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; |
237 | } else { |
238 | p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; |
239 | mismatch=kFALSE; |
240 | } |
241 | |
242 | // Check for particles heavier than (AliPID::kSPECIES - 1) |
243 | // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE; |
244 | |
245 | } |
246 | |
247 | if (mismatch){ |
248 | for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES; |
249 | return kDetNoSignal; |
250 | } |
251 | |
252 | |
253 | return kDetPidOk; |
254 | } |
255 | //______________________________________________________________________________ |
256 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const |
257 | { |
258 | // |
259 | // Compute PID response for the TPC |
260 | // |
261 | |
262 | // set flat distribution (no decision) |
263 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
264 | |
265 | // check quality of the track |
266 | if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal; |
267 | |
268 | Double_t mom = track->GetTPCmomentum(); |
269 | |
270 | Double_t dedx=track->GetTPCsignal(); |
271 | Bool_t mismatch=kTRUE/*, heavy=kTRUE*/; |
272 | |
273 | for (Int_t j=0; j<AliPID::kSPECIES; j++) { |
274 | AliPID::EParticleType type=AliPID::EParticleType(j); |
275 | Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type); |
276 | Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type); |
277 | if (TMath::Abs(dedx-bethe) > fRange*sigma) { |
278 | p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; |
279 | } else { |
280 | p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; |
281 | mismatch=kFALSE; |
282 | } |
283 | |
284 | // TODO: Light nuclei, also in TPC pid response |
285 | |
286 | // Check for particles heavier than (AliPID::kSPECIES - 1) |
287 | // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE; |
288 | |
289 | } |
290 | |
291 | if (mismatch){ |
292 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
293 | return kDetNoSignal; |
294 | } |
295 | |
296 | return kDetPidOk; |
297 | } |
298 | //______________________________________________________________________________ |
299 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const |
300 | { |
301 | // |
302 | // Compute PID response for the |
303 | // |
304 | |
305 | // set flat distribution (no decision) |
306 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
307 | |
308 | if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal; |
309 | if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal; |
310 | |
311 | Double_t time[AliPID::kSPECIESN]; |
312 | track->GetIntegratedTimes(time); |
313 | |
314 | Double_t sigma[AliPID::kSPECIES]; |
315 | for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) { |
316 | sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart)); |
317 | } |
318 | |
319 | Bool_t mismatch = kTRUE/*, heavy = kTRUE*/; |
320 | for (Int_t j=0; j<AliPID::kSPECIES; j++) { |
321 | AliPID::EParticleType type=AliPID::EParticleType(j); |
322 | Double_t nsigmas=NumberOfSigmasTOF(track,type); |
323 | |
324 | Double_t sig = sigma[j]; |
325 | if (TMath::Abs(nsigmas) > (fRange+2)) { |
326 | p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig; |
327 | } else |
328 | p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig; |
329 | |
330 | if (TMath::Abs(nsigmas)<5.){ |
331 | Double_t nsigmasTPC=NumberOfSigmasTPC(track,type); |
332 | if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE; |
333 | } |
334 | } |
335 | |
336 | if (mismatch){ |
337 | return kDetMismatch; |
338 | } |
339 | |
340 | // TODO: Light nuclei |
341 | |
342 | return kDetPidOk; |
343 | } |
344 | //______________________________________________________________________________ |
ea235c90 |
345 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const |
fd21ec8d |
346 | { |
347 | // |
348 | // Compute PID response for the |
349 | // |
350 | |
351 | // set flat distribution (no decision) |
352 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
ea235c90 |
353 | if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal; |
354 | |
355 | Float_t mom[6]; |
356 | Double_t dedx[48]; // Allocate space for the maximum number of TRD slices |
357 | Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1; |
358 | AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", fTRDslicesForPID[0], fTRDslicesForPID[1], nslices)); |
359 | for(UInt_t ilayer = 0; ilayer < 6; ilayer++){ |
360 | mom[ilayer] = track->GetTRDmomentum(ilayer); |
361 | for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){ |
362 | dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice); |
363 | } |
364 | } |
365 | fTRDResponse.GetResponse(nslices, dedx, mom, p); |
366 | return kDetPidOk; |
fd21ec8d |
367 | } |
368 | //______________________________________________________________________________ |
369 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability(const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const |
370 | { |
371 | // |
372 | // Compute PID response for the EMCAL |
373 | // |
374 | |
375 | // set flat distribution (no decision) |
376 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
377 | return kDetNoSignal; |
378 | } |
379 | //______________________________________________________________________________ |
380 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const |
381 | { |
382 | // |
383 | // Compute PID response for the PHOS |
384 | // |
385 | |
386 | // set flat distribution (no decision) |
387 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
388 | return kDetNoSignal; |
389 | } |
390 | //______________________________________________________________________________ |
ea235c90 |
391 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const |
fd21ec8d |
392 | { |
393 | // |
394 | // Compute PID response for the HMPID |
395 | // |
396 | |
397 | // set flat distribution (no decision) |
398 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
ea235c90 |
399 | if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal; |
400 | |
401 | track->GetHMPIDpid(p); |
402 | |
403 | return kDetPidOk; |
fd21ec8d |
404 | } |
405 | |
4ec8e76d |
406 | //______________________________________________________________________________ |
407 | void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass) |
408 | { |
409 | // |
410 | // Apply settings for the current event |
411 | // |
412 | fRecoPass=pass; |
413 | |
414 | if (!event) return; |
415 | fRun=event->GetRunNumber(); |
416 | |
417 | if (fRun!=fOldRun){ |
418 | ExecNewRun(); |
419 | fOldRun=fRun; |
420 | } |
421 | |
422 | //TPC resolution parametrisation PbPb |
423 | if ( fResolutionCorrection ){ |
424 | Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event)); |
425 | fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04); |
426 | } |
427 | |
428 | //TOF resolution |
429 | SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType); |
430 | |
431 | } |
432 | |
433 | //______________________________________________________________________________ |
434 | void AliPIDResponse::ExecNewRun() |
435 | { |
436 | // |
437 | // Things to Execute upon a new run |
438 | // |
439 | SetRecoInfo(); |
440 | |
441 | SetITSParametrisation(); |
442 | |
443 | SetTPCPidResponseMaster(); |
444 | SetTPCParametrisation(); |
53d016dc |
445 | |
446 | SetTRDPidResponseMaster(); |
447 | InitializeTRDResponse(); |
4ec8e76d |
448 | |
449 | fTOFResponse.SetTimeResolution(fTOFres); |
450 | } |
451 | |
452 | //_____________________________________________________ |
453 | Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event) |
454 | { |
455 | // |
456 | // Get TPC multiplicity in bins of 150 |
457 | // |
458 | |
459 | const AliVVertex* vertexTPC = event->GetPrimaryVertex(); |
460 | Double_t tpcMulti=0.; |
461 | if(vertexTPC){ |
462 | Double_t vertexContribTPC=vertexTPC->GetNContributors(); |
463 | tpcMulti=vertexContribTPC/150.; |
464 | if (tpcMulti>20.) tpcMulti=20.; |
465 | } |
466 | |
467 | return tpcMulti; |
468 | } |
469 | |
470 | //______________________________________________________________________________ |
471 | void AliPIDResponse::SetRecoInfo() |
472 | { |
473 | // |
474 | // Set reconstruction information |
475 | // |
476 | |
477 | //reset information |
478 | fLHCperiod=""; |
479 | fMCperiodTPC=""; |
480 | |
481 | fBeamType=""; |
482 | |
483 | fBeamType="PP"; |
484 | |
53d016dc |
485 | TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*"); |
4ec8e76d |
486 | //find the period by run number (UGLY, but not stored in ESD and AOD... ) |
487 | if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; } |
488 | else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; } |
489 | else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; } |
490 | else if (fRun>=127719&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; } |
491 | else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; } |
492 | else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; } |
ea235c90 |
493 | else if (fRun>=136851&&fRun<=139517) { |
494 | fLHCperiod="LHC10H"; |
495 | fMCperiodTPC="LHC10H8"; |
496 | if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10"; |
497 | fBeamType="PBPB"; |
498 | } |
4ec8e76d |
499 | else if (fRun>=139699) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; } |
80ab5635 |
500 | |
501 | //exception new pp MC productions from 2011 |
502 | if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2"; |
4ec8e76d |
503 | } |
504 | |
505 | //______________________________________________________________________________ |
506 | void AliPIDResponse::SetITSParametrisation() |
507 | { |
508 | // |
509 | // Set the ITS parametrisation |
510 | // |
511 | } |
512 | |
513 | //______________________________________________________________________________ |
514 | void AliPIDResponse::SetTPCPidResponseMaster() |
515 | { |
516 | // |
517 | // Load the TPC pid response functions from the OADB |
518 | // |
09b50a42 |
519 | //don't load twice for the moment |
520 | if (fArrPidResponseMaster) return; |
521 | |
522 | |
4ec8e76d |
523 | //reset the PID response functions |
524 | delete fArrPidResponseMaster; |
525 | fArrPidResponseMaster=0x0; |
526 | |
527 | TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data())); |
528 | |
ea235c90 |
529 | TFile *f=TFile::Open(fileName.Data()); |
530 | if (f && f->IsOpen() && !f->IsZombie()){ |
531 | fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse")); |
4ec8e76d |
532 | } |
ea235c90 |
533 | delete f; |
4ec8e76d |
534 | |
535 | if (!fArrPidResponseMaster){ |
536 | AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data())); |
537 | return; |
538 | } |
539 | fArrPidResponseMaster->SetOwner(); |
540 | } |
541 | |
542 | //______________________________________________________________________________ |
543 | void AliPIDResponse::SetTPCParametrisation() |
544 | { |
545 | // |
546 | // Change BB parametrisation for current run |
547 | // |
548 | |
549 | if (fLHCperiod.IsNull()) { |
550 | AliFatal("No period set, not changing parametrisation"); |
551 | return; |
552 | } |
553 | |
554 | // |
555 | // Set default parametrisations for data and MC |
556 | // |
557 | |
558 | //data type |
559 | TString datatype="DATA"; |
560 | //in case of mc fRecoPass is per default 1 |
561 | if (fIsMC) { |
562 | datatype="MC"; |
563 | fRecoPass=1; |
564 | } |
565 | |
566 | // |
567 | //reset old splines |
568 | // |
569 | for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){ |
570 | fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0); |
571 | } |
572 | |
573 | // |
574 | //set the new PID splines |
575 | // |
576 | TString period=fLHCperiod; |
577 | if (fArrPidResponseMaster){ |
578 | TObject *grAll=0x0; |
579 | //for MC don't use period information |
580 | // if (fIsMC) period="[A-Z0-9]*"; |
581 | //for MC use MC period information |
582 | if (fIsMC) period=fMCperiodTPC; |
583 | //pattern for the default entry (valid for all particles) |
584 | TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data())); |
585 | |
586 | //loop over entries and filter them |
587 | for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){ |
588 | TObject *responseFunction=fArrPidResponseMaster->At(iresp); |
589 | if (responseFunction==0x0) continue; |
590 | TString responseName=responseFunction->GetName(); |
591 | |
592 | if (!reg.MatchB(responseName)) continue; |
593 | |
594 | TObjArray *arr=reg.MatchS(responseName); |
595 | TString particleName=arr->At(1)->GetName(); |
596 | delete arr; |
597 | if (particleName.IsNull()) continue; |
598 | if (particleName=="ALL") grAll=responseFunction; |
599 | else { |
600 | //find particle id |
601 | for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){ |
602 | TString particle=AliPID::ParticleName(ispec); |
603 | particle.ToUpper(); |
604 | if ( particle == particleName ){ |
605 | fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction); |
606 | fTPCResponse.SetUseDatabase(kTRUE); |
607 | AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName())); |
608 | break; |
609 | } |
610 | } |
611 | } |
612 | } |
613 | |
614 | //set default response function to all particles which don't have a specific one |
615 | if (grAll){ |
616 | for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){ |
617 | if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){ |
618 | fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll); |
619 | AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName())); |
620 | } |
621 | } |
622 | } |
623 | } |
624 | |
625 | // |
626 | // Setup resolution parametrisation |
627 | // |
628 | |
629 | //default |
630 | fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04); |
631 | |
632 | if (fRun>=122195){ |
633 | fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02); |
634 | } |
23425eb2 |
635 | if (fArrPidResponseMaster) |
4ec8e76d |
636 | fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data())); |
637 | |
638 | if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName())); |
639 | } |
640 | |
ea235c90 |
641 | //______________________________________________________________________________ |
642 | void AliPIDResponse::SetTRDPidResponseMaster() |
643 | { |
644 | // |
645 | // Load the TRD pid params and references from the OADB |
646 | // |
647 | if(fTRDPIDParams) return; |
53d016dc |
648 | AliOADBContainer contParams("contParams"); |
649 | |
ea235c90 |
650 | contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams"); |
651 | fTRDPIDParams = (TObjArray *)contParams.GetObject(fRun); |
652 | |
53d016dc |
653 | AliOADBContainer contRefs("contRefs"); |
ea235c90 |
654 | contRefs.InitFromFile(Form("%s/COMMON/PID/dReferencesLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference"); |
655 | fTRDPIDReference = (AliTRDPIDReference *)contRefs.GetObject(fRun); |
656 | } |
657 | |
658 | //______________________________________________________________________________ |
659 | void AliPIDResponse::InitializeTRDResponse(){ |
660 | // |
661 | // Set PID Params and references to the TRD PID response |
662 | // |
663 | fTRDResponse.SetPIDParams(fTRDPIDParams); |
664 | fTRDResponse.Load(fTRDPIDReference); |
665 | if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){ |
666 | fTRDslicesForPID[0] = 0; |
667 | fTRDslicesForPID[1] = 7; |
668 | } |
669 | } |
670 | |
671 | //_________________________________________________________________________ |
672 | Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const { |
673 | // |
674 | // Check whether track is identified as electron under a given electron efficiency hypothesis |
675 | // |
676 | Double_t probs[AliPID::kSPECIES]; |
677 | ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs); |
678 | |
679 | Int_t ntracklets=0; |
680 | Double_t p = 0; |
681 | for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){ |
682 | if(vtrack->GetTRDmomentum(iPl) > 0.){ |
683 | ntracklets++; |
684 | p = vtrack->GetTRDmomentum(iPl); |
685 | } |
686 | } |
687 | |
688 | return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel); |
689 | } |
690 | |