Small update of the calibration code for the case of runs without TRD (Raphaelle)
[u/mrichter/AliRoot.git] / STEER / AliPIDResponse.cxx
CommitLineData
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>
29bf19f2 37
38#include "AliPIDResponse.h"
39
40ClassImp(AliPIDResponse);
41
4ec8e76d 42AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
43TNamed("PIDResponse","PIDResponse"),
44fITSResponse(isMC),
45fTPCResponse(),
46fTRDResponse(),
47fTOFResponse(),
fd21ec8d 48fRange(5.),
49fITSPIDmethod(kITSTruncMean),
4ec8e76d 50fIsMC(isMC),
51fOADBPath(),
52fBeamType("PP"),
53fLHCperiod(),
54fMCperiodTPC(),
fd21ec8d 55fMCperiodUser(),
4ec8e76d 56fRecoPass(0),
fd21ec8d 57fRecoPassUser(-1),
4ec8e76d 58fRun(0),
59fOldRun(0),
60fArrPidResponseMaster(0x0),
61fResolutionCorrection(0x0),
62fTOFTimeZeroType(kBest_T0),
63fTOFres(100.)
64{
65 //
66 // default ctor
67 //
68 AliLog::SetClassDebugLevel("AliPIDResponse",10);
69 AliLog::SetClassDebugLevel("AliESDpid",10);
70 AliLog::SetClassDebugLevel("AliAODpidUtil",10);
71}
72
73//______________________________________________________________________________
74AliPIDResponse::~AliPIDResponse()
75{
76 //
77 // dtor
78 //
79 delete fArrPidResponseMaster;
80}
81
82//______________________________________________________________________________
83AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
84TNamed(other),
85fITSResponse(other.fITSResponse),
86fTPCResponse(other.fTPCResponse),
87fTRDResponse(other.fTRDResponse),
88fTOFResponse(other.fTOFResponse),
fd21ec8d 89fRange(other.fRange),
90fITSPIDmethod(other.fITSPIDmethod),
4ec8e76d 91fIsMC(other.fIsMC),
92fOADBPath(other.fOADBPath),
93fBeamType("PP"),
94fLHCperiod(),
95fMCperiodTPC(),
fd21ec8d 96fMCperiodUser(other.fMCperiodUser),
4ec8e76d 97fRecoPass(0),
fd21ec8d 98fRecoPassUser(other.fRecoPassUser),
4ec8e76d 99fRun(0),
100fOldRun(0),
101fArrPidResponseMaster(0x0),
102fResolutionCorrection(0x0),
103fTOFTimeZeroType(AliPIDResponse::kBest_T0),
104fTOFres(100.)
105{
106 //
107 // copy ctor
108 //
109}
110
111//______________________________________________________________________________
112AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
113{
114 //
115 // copy ctor
116 //
117 if(this!=&other) {
118 delete fArrPidResponseMaster;
119 TNamed::operator=(other);
120 fITSResponse=other.fITSResponse;
121 fTPCResponse=other.fTPCResponse;
122 fTRDResponse=other.fTRDResponse;
123 fTOFResponse=other.fTOFResponse;
fd21ec8d 124 fRange=other.fRange;
125 fITSPIDmethod=other.fITSPIDmethod;
4ec8e76d 126 fOADBPath=other.fOADBPath;
127 fIsMC=other.fIsMC;
128 fBeamType="PP";
129 fLHCperiod="";
130 fMCperiodTPC="";
fd21ec8d 131 fMCperiodUser=other.fMCperiodUser;
4ec8e76d 132 fRecoPass=0;
fd21ec8d 133 fRecoPassUser=other.fRecoPassUser;
4ec8e76d 134 fRun=0;
135 fOldRun=0;
136 fArrPidResponseMaster=0x0;
137 fResolutionCorrection=0x0;
138 fTOFTimeZeroType=AliPIDResponse::kBest_T0;
139 fTOFres=100.;
140 }
141 return *this;
142}
143
144//______________________________________________________________________________
fd21ec8d 145Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
146{
147 //
148 // NumberOfSigmas for 'detCode'
149 //
150
151 switch (detCode){
152 case kDetITS: return NumberOfSigmasITS(track, type); break;
153 case kDetTPC: return NumberOfSigmasTPC(track, type); break;
154 case kDetTOF: return NumberOfSigmasTOF(track, type); break;
155// case kDetTRD: return ComputeTRDProbability(track, type); break;
156// case kDetPHOS: return ComputePHOSProbability(track, type); break;
157// case kDetEMCAL: return ComputeEMCALProbability(track, type); break;
158// case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
159 default: return -999.;
160 }
161
162}
163
164//______________________________________________________________________________
165AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
166{
167 //
168 // Compute PID response of 'detCode'
169 //
170
171 switch (detCode){
172 case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
173 case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
174 case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
175 case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
176 case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
177 case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
178 case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
179 default: return kDetNoSignal;
180 }
181}
182
183//______________________________________________________________________________
184AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
185{
186 //
187 // Compute PID response for the ITS
188 //
189
190 // set flat distribution (no decision)
191 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
192
193 if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
194 (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
195
196 Double_t mom=track->P();
197 Double_t dedx=track->GetITSsignal();
198 Bool_t isSA=kTRUE;
199 Double_t momITS=mom;
200 ULong_t trStatus=track->GetStatus();
201 if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
202 UChar_t clumap=track->GetITSClusterMap();
203 Int_t nPointsForPid=0;
204 for(Int_t i=2; i<6; i++){
205 if(clumap&(1<<i)) ++nPointsForPid;
206 }
207
208 if(nPointsForPid<3) { // track not to be used for combined PID purposes
209 // track->ResetStatus(AliVTrack::kITSpid);
210 return kDetNoSignal;
211 }
212
213 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
214 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
215 Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
216 Double_t bethe=fITSResponse.Bethe(momITS,mass);
217 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
218 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
219 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
220 } else {
221 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
222 mismatch=kFALSE;
223 }
224
225 // Check for particles heavier than (AliPID::kSPECIES - 1)
226 // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
227
228 }
229
230 if (mismatch){
231 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
232 return kDetNoSignal;
233 }
234
235
236 return kDetPidOk;
237}
238//______________________________________________________________________________
239AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
240{
241 //
242 // Compute PID response for the TPC
243 //
244
245 // set flat distribution (no decision)
246 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
247
248 // check quality of the track
249 if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
250
251 Double_t mom = track->GetTPCmomentum();
252
253 Double_t dedx=track->GetTPCsignal();
254 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
255
256 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
257 AliPID::EParticleType type=AliPID::EParticleType(j);
258 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
259 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
260 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
261 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
262 } else {
263 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
264 mismatch=kFALSE;
265 }
266
267 // TODO: Light nuclei, also in TPC pid response
268
269 // Check for particles heavier than (AliPID::kSPECIES - 1)
270// if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
271
272 }
273
274 if (mismatch){
275 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
276 return kDetNoSignal;
277 }
278
279 return kDetPidOk;
280}
281//______________________________________________________________________________
282AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
283{
284 //
285 // Compute PID response for the
286 //
287
288 // set flat distribution (no decision)
289 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
290
291 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
292 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
293
294 Double_t time[AliPID::kSPECIESN];
295 track->GetIntegratedTimes(time);
296
297 Double_t sigma[AliPID::kSPECIES];
298 for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
299 sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
300 }
301
302 Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
303 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
304 AliPID::EParticleType type=AliPID::EParticleType(j);
305 Double_t nsigmas=NumberOfSigmasTOF(track,type);
306
307 Double_t sig = sigma[j];
308 if (TMath::Abs(nsigmas) > (fRange+2)) {
309 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
310 } else
311 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
312
313 if (TMath::Abs(nsigmas)<5.){
314 Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
315 if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
316 }
317 }
318
319 if (mismatch){
320 return kDetMismatch;
321 }
322
323 // TODO: Light nuclei
324
325 return kDetPidOk;
326}
327//______________________________________________________________________________
328AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
329{
330 //
331 // Compute PID response for the
332 //
333
334 // set flat distribution (no decision)
335 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
336 return kDetNoSignal;
337}
338//______________________________________________________________________________
339AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability(const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
340{
341 //
342 // Compute PID response for the EMCAL
343 //
344
345 // set flat distribution (no decision)
346 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
347 return kDetNoSignal;
348}
349//______________________________________________________________________________
350AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
351{
352 //
353 // Compute PID response for the PHOS
354 //
355
356 // set flat distribution (no decision)
357 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
358 return kDetNoSignal;
359}
360//______________________________________________________________________________
361AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
362{
363 //
364 // Compute PID response for the HMPID
365 //
366
367 // set flat distribution (no decision)
368 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
369 return kDetNoSignal;
370}
371
372//______________________________________________________________________________
4ec8e76d 373void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
374{
375 //
376 // Apply settings for the current event
377 //
378 fRecoPass=pass;
379
380 if (!event) return;
381 fRun=event->GetRunNumber();
382
383 if (fRun!=fOldRun){
384 ExecNewRun();
385 fOldRun=fRun;
386 }
387
388 //TPC resolution parametrisation PbPb
389 if ( fResolutionCorrection ){
390 Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
391 fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
392 }
393
394 //TOF resolution
395 SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
396
397}
398
399//______________________________________________________________________________
400void AliPIDResponse::ExecNewRun()
401{
402 //
403 // Things to Execute upon a new run
404 //
405 SetRecoInfo();
406
407 SetITSParametrisation();
408
409 SetTPCPidResponseMaster();
410 SetTPCParametrisation();
411
412 fTOFResponse.SetTimeResolution(fTOFres);
413}
414
415//_____________________________________________________
416Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
417{
418 //
419 // Get TPC multiplicity in bins of 150
420 //
421
422 const AliVVertex* vertexTPC = event->GetPrimaryVertex();
423 Double_t tpcMulti=0.;
424 if(vertexTPC){
425 Double_t vertexContribTPC=vertexTPC->GetNContributors();
426 tpcMulti=vertexContribTPC/150.;
427 if (tpcMulti>20.) tpcMulti=20.;
428 }
429
430 return tpcMulti;
431}
432
433//______________________________________________________________________________
434void AliPIDResponse::SetRecoInfo()
435{
436 //
437 // Set reconstruction information
438 //
439
440 //reset information
441 fLHCperiod="";
442 fMCperiodTPC="";
443
444 fBeamType="";
445
446 fBeamType="PP";
447
448 //find the period by run number (UGLY, but not stored in ESD and AOD... )
449 if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
450 else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
451 else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
452 else if (fRun>=127719&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
453 else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
454 else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
455 else if (fRun>=136851&&fRun<=139517) { fLHCperiod="LHC10H"; fMCperiodTPC="LHC10H8"; fBeamType="PBPB"; }
456 else if (fRun>=139699) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
457
458}
459
460//______________________________________________________________________________
461void AliPIDResponse::SetITSParametrisation()
462{
463 //
464 // Set the ITS parametrisation
465 //
466}
467
468//______________________________________________________________________________
469void AliPIDResponse::SetTPCPidResponseMaster()
470{
471 //
472 // Load the TPC pid response functions from the OADB
473 //
09b50a42 474 //don't load twice for the moment
475 if (fArrPidResponseMaster) return;
476
477
4ec8e76d 478 //reset the PID response functions
479 delete fArrPidResponseMaster;
480 fArrPidResponseMaster=0x0;
481
482 TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
483
484 TFile f(fileName.Data());
485 if (f.IsOpen() && !f.IsZombie()){
486 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f.Get("TPCPIDResponse"));
487 f.Close();
488 }
489
490 if (!fArrPidResponseMaster){
491 AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
492 return;
493 }
494 fArrPidResponseMaster->SetOwner();
495}
496
497//______________________________________________________________________________
498void AliPIDResponse::SetTPCParametrisation()
499{
500 //
501 // Change BB parametrisation for current run
502 //
503
504 if (fLHCperiod.IsNull()) {
505 AliFatal("No period set, not changing parametrisation");
506 return;
507 }
508
509 //
510 // Set default parametrisations for data and MC
511 //
512
513 //data type
514 TString datatype="DATA";
515 //in case of mc fRecoPass is per default 1
516 if (fIsMC) {
517 datatype="MC";
518 fRecoPass=1;
519 }
520
521 //
522 //reset old splines
523 //
524 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
525 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
526 }
527
528 //
529 //set the new PID splines
530 //
531 TString period=fLHCperiod;
532 if (fArrPidResponseMaster){
533 TObject *grAll=0x0;
534 //for MC don't use period information
535// if (fIsMC) period="[A-Z0-9]*";
536 //for MC use MC period information
537 if (fIsMC) period=fMCperiodTPC;
538//pattern for the default entry (valid for all particles)
539 TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
540
541 //loop over entries and filter them
542 for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
543 TObject *responseFunction=fArrPidResponseMaster->At(iresp);
544 if (responseFunction==0x0) continue;
545 TString responseName=responseFunction->GetName();
546
547 if (!reg.MatchB(responseName)) continue;
548
549 TObjArray *arr=reg.MatchS(responseName);
550 TString particleName=arr->At(1)->GetName();
551 delete arr;
552 if (particleName.IsNull()) continue;
553 if (particleName=="ALL") grAll=responseFunction;
554 else {
555 //find particle id
556 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
557 TString particle=AliPID::ParticleName(ispec);
558 particle.ToUpper();
559 if ( particle == particleName ){
560 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
561 fTPCResponse.SetUseDatabase(kTRUE);
562 AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
563 break;
564 }
565 }
566 }
567 }
568
569 //set default response function to all particles which don't have a specific one
570 if (grAll){
571 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
572 if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
573 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
574 AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
575 }
576 }
577 }
578 }
579
580 //
581 // Setup resolution parametrisation
582 //
583
584 //default
585 fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
586
587 if (fRun>=122195){
588 fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
589 }
23425eb2 590 if (fArrPidResponseMaster)
4ec8e76d 591 fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
592
593 if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
594}
595