1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */
18 //-----------------------------------------------------------------
19 // Base class for handling the pid response //
20 // functions of all detectors //
21 // and give access to the nsigmas //
23 // Origin: Jens Wiechula, Uni Tuebingen, jens.wiechula@cern.ch //
24 //-----------------------------------------------------------------
27 #include <TObjArray.h>
34 #include <AliVEvent.h>
35 #include <AliVTrack.h>
38 #include <AliOADBContainer.h>
39 #include <AliTRDPIDResponseObject.h>
40 #include <AliTOFPIDParams.h>
42 #include "AliPIDResponse.h"
44 #include "AliCentrality.h"
46 ClassImp(AliPIDResponse);
48 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
49 TNamed("PIDResponse","PIDResponse"),
56 fITSPIDmethod(kITSTruncMean),
68 fArrPidResponseMaster(0x0),
69 fResolutionCorrection(0x0),
70 fTRDPIDResponseObject(0x0),
81 AliLog::SetClassDebugLevel("AliPIDResponse",0);
82 AliLog::SetClassDebugLevel("AliESDpid",0);
83 AliLog::SetClassDebugLevel("AliAODpidUtil",0);
85 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
88 //______________________________________________________________________________
89 AliPIDResponse::~AliPIDResponse()
94 delete fArrPidResponseMaster;
95 delete fTRDPIDResponseObject;
96 if (fTOFPIDParams) delete fTOFPIDParams;
99 //______________________________________________________________________________
100 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
102 fITSResponse(other.fITSResponse),
103 fTPCResponse(other.fTPCResponse),
104 fTRDResponse(other.fTRDResponse),
105 fTOFResponse(other.fTOFResponse),
106 fEMCALResponse(other.fEMCALResponse),
107 fRange(other.fRange),
108 fITSPIDmethod(other.fITSPIDmethod),
110 fOADBPath(other.fOADBPath),
114 fMCperiodUser(other.fMCperiodUser),
117 fRecoPassUser(other.fRecoPassUser),
120 fArrPidResponseMaster(0x0),
121 fResolutionCorrection(0x0),
122 fTRDPIDResponseObject(0x0),
125 fEMCALPIDParams(0x0),
127 fCurrCentrality(0.0),
128 fTuneMConData(kFALSE)
133 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
136 //______________________________________________________________________________
137 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
143 delete fArrPidResponseMaster;
144 TNamed::operator=(other);
145 fITSResponse=other.fITSResponse;
146 fTPCResponse=other.fTPCResponse;
147 fTRDResponse=other.fTRDResponse;
148 fTOFResponse=other.fTOFResponse;
149 fEMCALResponse=other.fEMCALResponse;
151 fITSPIDmethod=other.fITSPIDmethod;
152 fOADBPath=other.fOADBPath;
157 fMCperiodUser=other.fMCperiodUser;
160 fRecoPassUser=other.fRecoPassUser;
163 fArrPidResponseMaster=0x0;
164 fResolutionCorrection=0x0;
165 fTRDPIDResponseObject=0x0;
167 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
170 fCurrentEvent=other.fCurrentEvent;
175 //______________________________________________________________________________
176 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
179 // NumberOfSigmas for 'detCode'
183 case kDetITS: return NumberOfSigmasITS(track, type); break;
184 case kDetTPC: return NumberOfSigmasTPC(track, type); break;
185 case kDetTOF: return NumberOfSigmasTOF(track, type); break;
186 // case kDetTRD: return ComputeTRDProbability(track, type); break;
187 // case kDetPHOS: return ComputePHOSProbability(track, type); break;
188 // case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
189 // case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
190 default: return -999.;
195 //______________________________________________________________________________
196 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
198 AliVCluster *matchedClus = NULL;
203 Double_t fClsE = -1.;
205 Int_t nMatchClus = -1;
209 nMatchClus = track->GetEMCALcluster();
214 charge = track->Charge();
216 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
220 // matched cluster is EMCAL
221 if(matchedClus->IsEMCAL()){
223 fClsE = matchedClus->E();
227 // NSigma value really meaningful only for electrons!
228 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
237 //______________________________________________________________________________
238 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
240 AliVCluster *matchedClus = NULL;
245 Double_t fClsE = -1.;
247 // initialize eop and shower shape parameters
249 for(Int_t i = 0; i < 4; i++){
250 showershape[i] = -1.;
253 Int_t nMatchClus = -1;
257 nMatchClus = track->GetEMCALcluster();
262 charge = track->Charge();
264 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
268 // matched cluster is EMCAL
269 if(matchedClus->IsEMCAL()){
271 fClsE = matchedClus->E();
274 // fill used EMCAL variables here
276 showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
277 showershape[1] = matchedClus->GetM02(); // long axis
278 showershape[2] = matchedClus->GetM20(); // short axis
279 showershape[3] = matchedClus->GetDispersion(); // dispersion
281 // NSigma value really meaningful only for electrons!
282 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
290 //______________________________________________________________________________
291 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
294 // Compute PID response of 'detCode'
298 case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
299 case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
300 case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
301 case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
302 case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
303 case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
304 case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
305 default: return kDetNoSignal;
309 //______________________________________________________________________________
310 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
313 // Compute PID response for the ITS
316 // set flat distribution (no decision)
317 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
319 if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
320 (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
322 Double_t mom=track->P();
323 Double_t dedx=track->GetITSsignal();
326 ULong_t trStatus=track->GetStatus();
327 if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
328 UChar_t clumap=track->GetITSClusterMap();
329 Int_t nPointsForPid=0;
330 for(Int_t i=2; i<6; i++){
331 if(clumap&(1<<i)) ++nPointsForPid;
334 if(nPointsForPid<3) { // track not to be used for combined PID purposes
335 // track->ResetStatus(AliVTrack::kITSpid);
339 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
340 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
341 Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
342 Double_t bethe=fITSResponse.Bethe(momITS,mass);
343 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
344 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
345 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
347 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
351 // Check for particles heavier than (AliPID::kSPECIES - 1)
352 // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
357 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
364 //______________________________________________________________________________
365 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
368 // Compute PID response for the TPC
371 // set flat distribution (no decision)
372 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
374 // check quality of the track
375 if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
377 Double_t mom = track->GetTPCmomentum();
379 Double_t dedx=track->GetTPCsignal();
380 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
382 if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
384 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
385 AliPID::EParticleType type=AliPID::EParticleType(j);
386 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
387 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
388 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
389 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
391 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
395 // TODO: Light nuclei, also in TPC pid response
397 // Check for particles heavier than (AliPID::kSPECIES - 1)
398 // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
403 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
409 //______________________________________________________________________________
410 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
413 // Compute PID response for the
416 Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
418 // set flat distribution (no decision)
419 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
421 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
422 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
424 Double_t time[AliPID::kSPECIESN];
425 track->GetIntegratedTimes(time);
427 // Double_t sigma[AliPID::kSPECIES];
428 // for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
429 // sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
432 Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
433 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
434 AliPID::EParticleType type=AliPID::EParticleType(j);
435 Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
437 // Double_t sig = sigma[j];
438 Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
439 Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
440 if (TMath::Abs(nsigmas) > (fRange+2)) {
441 if(nsigmas < fTOFtail)
442 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
444 p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
446 if(nsigmas < fTOFtail)
447 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
449 p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
452 /* OLD Gaussian shape
453 if (TMath::Abs(nsigmas) > (fRange+2)) {
454 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
456 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
459 if (TMath::Abs(nsigmas)<5.){
460 Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
461 if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
469 // TODO: Light nuclei
473 //______________________________________________________________________________
474 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
477 // Compute PID response for the
480 // set flat distribution (no decision)
481 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
482 if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
485 Double_t dedx[48]; // Allocate space for the maximum number of TRD slices
486 Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
487 AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
488 for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
489 mom[ilayer] = track->GetTRDmomentum(ilayer);
490 for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
491 dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
494 fTRDResponse.GetResponse(nslices, dedx, mom, p);
497 //______________________________________________________________________________
498 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
501 // Compute PID response for the EMCAL
504 AliVCluster *matchedClus = NULL;
509 Double_t fClsE = -1.;
511 Int_t nMatchClus = -1;
515 nMatchClus = track->GetEMCALcluster();
521 charge = track->Charge();
523 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
527 // matched cluster is EMCAL
528 if(matchedClus->IsEMCAL()){
530 fClsE = matchedClus->E();
534 // compute the probabilities
535 if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){
537 // in case everything is OK
545 // in all other cases set flat distribution (no decision)
546 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
550 //______________________________________________________________________________
551 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
554 // Compute PID response for the PHOS
557 // set flat distribution (no decision)
558 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
561 //______________________________________________________________________________
562 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
565 // Compute PID response for the HMPID
568 // set flat distribution (no decision)
569 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
570 if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
572 track->GetHMPIDpid(p);
577 //______________________________________________________________________________
578 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
581 // Apply settings for the current event
588 fRun=event->GetRunNumber();
595 //TPC resolution parametrisation PbPb
596 if ( fResolutionCorrection ){
597 Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
598 fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
602 SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
605 // Get and set centrality
606 AliCentrality *centrality = event->GetCentrality();
608 fCurrCentrality = centrality->GetCentralityPercentile("V0M");
611 fCurrCentrality = -1;
615 //______________________________________________________________________________
616 void AliPIDResponse::ExecNewRun()
619 // Things to Execute upon a new run
623 SetITSParametrisation();
625 SetTPCPidResponseMaster();
626 SetTPCParametrisation();
628 SetTRDPidResponseMaster();
629 InitializeTRDResponse();
631 SetEMCALPidResponseMaster();
632 InitializeEMCALResponse();
634 SetTOFPidResponseMaster();
635 InitializeTOFResponse();
638 //_____________________________________________________
639 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
642 // Get TPC multiplicity in bins of 150
645 const AliVVertex* vertexTPC = event->GetPrimaryVertex();
646 Double_t tpcMulti=0.;
648 Double_t vertexContribTPC=vertexTPC->GetNContributors();
649 tpcMulti=vertexContribTPC/150.;
650 if (tpcMulti>20.) tpcMulti=20.;
656 //______________________________________________________________________________
657 void AliPIDResponse::SetRecoInfo()
660 // Set reconstruction information
672 TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
673 TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
675 //find the period by run number (UGLY, but not stored in ESD and AOD... )
676 if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
677 else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
678 else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
679 else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
680 else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
681 else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
682 else if (fRun>=136851&&fRun<=139517) {
684 fMCperiodTPC="LHC10H8";
685 if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
688 else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
689 //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
690 else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
691 else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
692 else if (fRun>=166529) {
694 fMCperiodTPC="LHC11A10";
696 if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
700 //exception new pp MC productions from 2011
701 if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
702 // exception for 11f1
703 if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
706 //______________________________________________________________________________
707 void AliPIDResponse::SetITSParametrisation()
710 // Set the ITS parametrisation
714 //______________________________________________________________________________
715 void AliPIDResponse::SetTPCPidResponseMaster()
718 // Load the TPC pid response functions from the OADB
720 //don't load twice for the moment
721 if (fArrPidResponseMaster) return;
724 //reset the PID response functions
725 delete fArrPidResponseMaster;
726 fArrPidResponseMaster=0x0;
728 TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
730 TFile *f=TFile::Open(fileName.Data());
731 if (f && f->IsOpen() && !f->IsZombie()){
732 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
736 if (!fArrPidResponseMaster){
737 AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
740 fArrPidResponseMaster->SetOwner();
743 //______________________________________________________________________________
744 void AliPIDResponse::SetTPCParametrisation()
747 // Change BB parametrisation for current run
750 if (fLHCperiod.IsNull()) {
751 AliFatal("No period set, not changing parametrisation");
756 // Set default parametrisations for data and MC
760 TString datatype="DATA";
761 //in case of mc fRecoPass is per default 1
763 if(!fTuneMConData) datatype="MC";
770 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
771 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
775 TString period=fLHCperiod;
776 if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
778 AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
781 //set the new PID splines
783 if (fArrPidResponseMaster){
784 Int_t recopass = fRecoPass;
785 if(fTuneMConData) recopass = fRecoPassUser;
787 //for MC don't use period information
788 // if (fIsMC) period="[A-Z0-9]*";
789 //for MC use MC period information
790 //pattern for the default entry (valid for all particles)
791 TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
793 //loop over entries and filter them
794 for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
795 TObject *responseFunction=fArrPidResponseMaster->At(iresp);
796 if (responseFunction==0x0) continue;
797 TString responseName=responseFunction->GetName();
799 if (!reg.MatchB(responseName)) continue;
801 TObjArray *arr=reg.MatchS(responseName);
802 TString particleName=arr->At(1)->GetName();
804 if (particleName.IsNull()) continue;
805 if (particleName=="ALL") grAll=responseFunction;
808 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
809 TString particle=AliPID::ParticleName(ispec);
811 if ( particle == particleName ){
812 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
813 fTPCResponse.SetUseDatabase(kTRUE);
814 AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
822 //set default response function to all particles which don't have a specific one
824 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
825 if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
826 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
827 AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
835 AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
838 // Setup resolution parametrisation
842 fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
845 fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
847 if (fArrPidResponseMaster)
848 fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
850 if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
853 //______________________________________________________________________________
854 void AliPIDResponse::SetTRDPidResponseMaster()
857 // Load the TRD pid params and references from the OADB
859 if(fTRDPIDResponseObject) return;
860 AliOADBContainer contParams("contParams");
862 Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
864 AliError("Failed initializing PID Response Object from OADB");
866 AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
867 fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
868 if(!fTRDPIDResponseObject){
869 AliError(Form("TRD Response not found in run %d", fRun));
873 AliOADBContainer contRefs("contRefs");
874 Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
876 AliInfo("Failed Loading References for TRD");
878 AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
879 fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
880 if(!fTRDPIDReference){
881 AliError(Form("TRD References not found in OADB Container for run %d", fRun));
887 //______________________________________________________________________________
888 void AliPIDResponse::InitializeTRDResponse(){
890 // Set PID Params and references to the TRD PID response
892 fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
896 void AliPIDResponse::SetTRDPIDmethod(AliTRDPIDResponse::ETRDPIDMethod method){
898 fTRDResponse.SetPIDmethod(method);
899 if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
900 // backward compatibility for setting with 8 slices
901 fTRDslicesForPID[0] = 0;
902 fTRDslicesForPID[1] = 7;
905 if(method==AliTRDPIDResponse::kLQ1D){
906 fTRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
907 fTRDslicesForPID[1] = 0;
909 if(method==AliTRDPIDResponse::kLQ2D){
910 fTRDslicesForPID[0] = 1;
911 fTRDslicesForPID[1] = 7;
914 AliDebug(1,Form("Slice Range set to %d - %d",fTRDslicesForPID[0],fTRDslicesForPID[1]));
917 //______________________________________________________________________________
918 void AliPIDResponse::SetTOFPidResponseMaster()
921 // Load the TOF pid params from the OADB
923 TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
924 if (oadbf->IsOpen()) {
925 AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
926 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
927 if (fTOFPIDParams) delete fTOFPIDParams;
928 fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
932 AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
938 //______________________________________________________________________________
939 void AliPIDResponse::InitializeTOFResponse(){
941 // Set PID Params to the TOF PID response
943 for (Int_t i=0;i<4;i++) {
944 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
946 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
951 //_________________________________________________________________________
952 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
954 // Check whether track is identified as electron under a given electron efficiency hypothesis
956 Double_t probs[AliPID::kSPECIES];
957 ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
959 Int_t ntracklets = vtrack->GetTRDntrackletsPID();
960 // Take mean of the TRD momenta in the given tracklets
961 Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
963 for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
964 if(vtrack->GetTRDmomentum(iPl) > 0.){
965 trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
968 p = TMath::Mean(nmomenta, trdmomenta);
970 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
973 //______________________________________________________________________________
974 void AliPIDResponse::SetEMCALPidResponseMaster()
977 // Load the EMCAL pid response functions from the OADB
979 TObjArray* fEMCALPIDParamsRun = NULL;
980 TObjArray* fEMCALPIDParamsPass = NULL;
982 if(fEMCALPIDParams) return;
983 AliOADBContainer contParams("contParams");
985 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
987 AliError("Failed initializing PID Params from OADB");
990 AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
992 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
993 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
994 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
996 if(!fEMCALPIDParams){
997 AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
998 AliInfo("Will take the standard LHC11d instead ...");
1000 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1001 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1002 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1004 if(!fEMCALPIDParams){
1005 AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1011 //______________________________________________________________________________
1012 void AliPIDResponse::InitializeEMCALResponse(){
1014 // Set PID Params to the EMCAL PID response
1016 fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1019 //_________________________________________________________________________
1020 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1022 // Set TOF response function
1023 // Input option for event_time used
1026 Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1027 if(t0spread < 10) t0spread = 80;
1029 // T0 from TOF algorithm
1031 Bool_t flagT0TOF=kFALSE;
1032 Bool_t flagT0T0=kFALSE;
1033 Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1034 Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1035 Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1038 Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1039 Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1040 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1041 estimatedT0event[i]=0.0;
1042 estimatedT0resolution[i]=0.0;
1043 startTimeMask[i] = 0;
1046 Float_t resT0A=75,resT0C=65,resT0AC=55;
1047 if(vevent->GetT0TOF()){ // check if T0 detector information is available
1052 AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1054 if (tofHeader) { // read global info and T0-TOF
1055 fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1056 t0spread = tofHeader->GetT0spread(); // read t0 sprad
1057 if(t0spread < 10) t0spread = 80;
1060 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1061 startTime[i]=tofHeader->GetDefaultEventTimeVal();
1062 startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1063 if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1066 TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1067 TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1068 TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1069 for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1070 Int_t icurrent = (Int_t)ibin->GetAt(j);
1071 startTime[icurrent]=t0Bin->GetAt(j);
1072 startTimeRes[icurrent]=t0ResBin->GetAt(j);
1073 if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1077 // for cut of 3 sigma on t0 spread
1078 Float_t t0cut = 3 * t0spread;
1079 if(t0cut < 500) t0cut = 500;
1081 if(option == kFILL_T0){ // T0-FILL is used
1082 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1083 estimatedT0event[i]=0.0;
1084 estimatedT0resolution[i]=t0spread;
1086 fTOFResponse.SetT0event(estimatedT0event);
1087 fTOFResponse.SetT0resolution(estimatedT0resolution);
1090 if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1092 fTOFResponse.SetT0event(startTime);
1093 fTOFResponse.SetT0resolution(startTimeRes);
1094 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1095 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1096 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1100 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1101 estimatedT0event[i]=0.0;
1102 estimatedT0resolution[i]=t0spread;
1103 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1105 fTOFResponse.SetT0event(estimatedT0event);
1106 fTOFResponse.SetT0resolution(estimatedT0resolution);
1109 else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1110 Float_t t0AC=-10000;
1114 t0AC= vevent->GetT0TOF()[0];
1115 t0A= vevent->GetT0TOF()[1];
1116 t0C= vevent->GetT0TOF()[2];
1119 Float_t t0t0Best = 0;
1120 Float_t t0t0BestRes = 9999;
1122 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1124 t0t0BestRes = resT0AC;
1127 else if(TMath::Abs(t0C) < t0cut){
1129 t0t0BestRes = resT0C;
1132 else if(TMath::Abs(t0A) < t0cut){
1134 t0t0BestRes = resT0A;
1138 if(flagT0TOF){ // if T0-TOF info is available
1139 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1140 if(t0t0BestRes < 999){
1141 if(startTimeRes[i] < t0spread){
1142 Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1143 Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1144 estimatedT0event[i]=t0best / wtot;
1145 estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1146 startTimeMask[i] = t0used+1;
1149 estimatedT0event[i]=t0t0Best;
1150 estimatedT0resolution[i]=t0t0BestRes;
1151 startTimeMask[i] = t0used;
1155 estimatedT0event[i]=startTime[i];
1156 estimatedT0resolution[i]=startTimeRes[i];
1157 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1159 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1161 fTOFResponse.SetT0event(estimatedT0event);
1162 fTOFResponse.SetT0resolution(estimatedT0resolution);
1164 else{ // if no T0-TOF info is available
1165 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1166 fTOFResponse.SetT0binMask(i,t0used);
1167 if(t0t0BestRes < 999){
1168 estimatedT0event[i]=t0t0Best;
1169 estimatedT0resolution[i]=t0t0BestRes;
1172 estimatedT0event[i]=0.0;
1173 estimatedT0resolution[i]=t0spread;
1176 fTOFResponse.SetT0event(estimatedT0event);
1177 fTOFResponse.SetT0resolution(estimatedT0resolution);
1181 else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1182 Float_t t0AC=-10000;
1186 t0AC= vevent->GetT0TOF()[0];
1187 t0A= vevent->GetT0TOF()[1];
1188 t0C= vevent->GetT0TOF()[2];
1191 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1192 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1193 estimatedT0event[i]=t0AC;
1194 estimatedT0resolution[i]=resT0AC;
1195 fTOFResponse.SetT0binMask(i,6);
1198 else if(TMath::Abs(t0C) < t0cut){
1199 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1200 estimatedT0event[i]=t0C;
1201 estimatedT0resolution[i]=resT0C;
1202 fTOFResponse.SetT0binMask(i,4);
1205 else if(TMath::Abs(t0A) < t0cut){
1206 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1207 estimatedT0event[i]=t0A;
1208 estimatedT0resolution[i]=resT0A;
1209 fTOFResponse.SetT0binMask(i,2);
1213 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1214 estimatedT0event[i]=0.0;
1215 estimatedT0resolution[i]=t0spread;
1216 fTOFResponse.SetT0binMask(i,0);
1219 fTOFResponse.SetT0event(estimatedT0event);
1220 fTOFResponse.SetT0resolution(estimatedT0resolution);
1222 delete [] startTime;
1223 delete [] startTimeRes;
1224 delete [] startTimeMask;
1225 delete [] estimatedT0event;
1226 delete [] estimatedT0resolution;