/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */ //----------------------------------------------------------------- // Base class for handling the pid response // // functions of all detectors // // and give access to the nsigmas // // // // Origin: Jens Wiechula, Uni Tuebingen, jens.wiechula@cern.ch // //----------------------------------------------------------------- #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "AliPIDResponse.h" #include "AliDetectorPID.h" #include "AliCentrality.h" ClassImp(AliPIDResponse); AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) : TNamed("PIDResponse","PIDResponse"), fITSResponse(isMC), fTPCResponse(), fTRDResponse(), fTOFResponse(), fEMCALResponse(), fRange(5.), fITSPIDmethod(kITSTruncMean), fIsMC(isMC), fCachePID(kTRUE), fOADBPath(), fCustomTPCpidResponse(), fBeamType("PP"), fLHCperiod(), fMCperiodTPC(), fMCperiodUser(), fCurrentFile(), fRecoPass(0), fRecoPassUser(-1), fRun(0), fOldRun(0), fResT0A(75.), fResT0C(65.), fResT0AC(55.), fArrPidResponseMaster(NULL), fResolutionCorrection(NULL), fOADBvoltageMaps(NULL), fTRDPIDResponseObject(NULL), fTOFtail(1.1), fTOFPIDParams(NULL), fEMCALPIDParams(NULL), fCurrentEvent(NULL), fCurrCentrality(0.0), fTuneMConData(kFALSE) { // // default ctor // AliLog::SetClassDebugLevel("AliPIDResponse",0); AliLog::SetClassDebugLevel("AliESDpid",0); AliLog::SetClassDebugLevel("AliAODpidUtil",0); } //______________________________________________________________________________ AliPIDResponse::~AliPIDResponse() { // // dtor // delete fArrPidResponseMaster; delete fTRDPIDResponseObject; delete fTOFPIDParams; } //______________________________________________________________________________ AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) : TNamed(other), fITSResponse(other.fITSResponse), fTPCResponse(other.fTPCResponse), fTRDResponse(other.fTRDResponse), fTOFResponse(other.fTOFResponse), fEMCALResponse(other.fEMCALResponse), fRange(other.fRange), fITSPIDmethod(other.fITSPIDmethod), fIsMC(other.fIsMC), fCachePID(other.fCachePID), fOADBPath(other.fOADBPath), fCustomTPCpidResponse(other.fCustomTPCpidResponse), fBeamType("PP"), fLHCperiod(), fMCperiodTPC(), fMCperiodUser(other.fMCperiodUser), fCurrentFile(), fRecoPass(0), fRecoPassUser(other.fRecoPassUser), fRun(0), fOldRun(0), fResT0A(75.), fResT0C(65.), fResT0AC(55.), fArrPidResponseMaster(NULL), fResolutionCorrection(NULL), fOADBvoltageMaps(NULL), fTRDPIDResponseObject(NULL), fTOFtail(1.1), fTOFPIDParams(NULL), fEMCALPIDParams(NULL), fCurrentEvent(NULL), fCurrCentrality(0.0), fTuneMConData(kFALSE) { // // copy ctor // } //______________________________________________________________________________ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other) { // // copy ctor // if(this!=&other) { delete fArrPidResponseMaster; TNamed::operator=(other); fITSResponse=other.fITSResponse; fTPCResponse=other.fTPCResponse; fTRDResponse=other.fTRDResponse; fTOFResponse=other.fTOFResponse; fEMCALResponse=other.fEMCALResponse; fRange=other.fRange; fITSPIDmethod=other.fITSPIDmethod; fOADBPath=other.fOADBPath; fCustomTPCpidResponse=other.fCustomTPCpidResponse; fIsMC=other.fIsMC; fCachePID=other.fCachePID; fBeamType="PP"; fLHCperiod=""; fMCperiodTPC=""; fMCperiodUser=other.fMCperiodUser; fCurrentFile=""; fRecoPass=0; fRecoPassUser=other.fRecoPassUser; fRun=0; fOldRun=0; fResT0A=75.; fResT0C=65.; fResT0AC=55.; fArrPidResponseMaster=NULL; fResolutionCorrection=NULL; fOADBvoltageMaps=NULL; fTRDPIDResponseObject=NULL; fEMCALPIDParams=NULL; fTOFtail=1.1; fTOFPIDParams=NULL; fCurrentEvent=other.fCurrentEvent; } return *this; } //______________________________________________________________________________ Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const { // // NumberOfSigmas for 'detCode' // switch (detCode){ case kDetITS: return NumberOfSigmasITS(track, type); break; case kDetTPC: return NumberOfSigmasTPC(track, type); break; case kDetTOF: return NumberOfSigmasTOF(track, type); break; case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break; default: return -999.; } } //______________________________________________________________________________ Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const { // // NumberOfSigmas for 'detCode' // return NumberOfSigmas((EDetCode)(1<GetDetectorPID(); const EDetector detector=kITS; if ( detPID && detPID->HasNumberOfSigmas(detector)){ return detPID->GetNumberOfSigmas(detector, type); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetNumberOfSigmas(detector, type); } return GetNumberOfSigmasITS(track, type); } //______________________________________________________________________________ Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const { // // Calculate the number of sigmas in the TPC // AliVTrack *track=(AliVTrack*)vtrack; // look for cached value first const AliDetectorPID *detPID=track->GetDetectorPID(); const EDetector detector=kTPC; if ( detPID && detPID->HasNumberOfSigmas(detector)){ return detPID->GetNumberOfSigmas(detector, type); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetNumberOfSigmas(detector, type); } return GetNumberOfSigmasTPC(track, type); } //______________________________________________________________________________ Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack, AliPID::EParticleType type, AliTPCPIDResponse::ETPCdEdxSource dedxSource) { //get number of sigmas according the selected TPC gain configuration scenario const AliVTrack *track=static_cast(vtrack); Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource); return nSigma; } //______________________________________________________________________________ Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const { // // Calculate the number of sigmas in the TOF // AliVTrack *track=(AliVTrack*)vtrack; // look for cached value first const AliDetectorPID *detPID=track->GetDetectorPID(); const EDetector detector=kTOF; if ( detPID && detPID->HasNumberOfSigmas(detector)){ return detPID->GetNumberOfSigmas(detector, type); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetNumberOfSigmas(detector, type); } return GetNumberOfSigmasTOF(track, type); } //______________________________________________________________________________ Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const { // // Calculate the number of sigmas in the EMCAL // AliVTrack *track=(AliVTrack*)vtrack; // look for cached value first const AliDetectorPID *detPID=track->GetDetectorPID(); const EDetector detector=kEMCAL; if ( detPID && detPID->HasNumberOfSigmas(detector)){ return detPID->GetNumberOfSigmas(detector, type); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetNumberOfSigmas(detector, type); } return GetNumberOfSigmasEMCAL(track, type); } //______________________________________________________________________________ Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const { // // emcal nsigma with eop and showershape // AliVTrack *track=(AliVTrack*)vtrack; AliVCluster *matchedClus = NULL; Double_t mom = -1.; Double_t pt = -1.; Double_t EovP = -1.; Double_t fClsE = -1.; // initialize eop and shower shape parameters eop = -1.; for(Int_t i = 0; i < 4; i++){ showershape[i] = -1.; } Int_t nMatchClus = -1; Int_t charge = 0; // Track matching nMatchClus = track->GetEMCALcluster(); if(nMatchClus > -1){ mom = track->P(); pt = track->Pt(); charge = track->Charge(); matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus); if(matchedClus){ // matched cluster is EMCAL if(matchedClus->IsEMCAL()){ fClsE = matchedClus->E(); EovP = fClsE/mom; // fill used EMCAL variables here eop = EovP; // E/p showershape[0] = matchedClus->GetNCells(); // number of cells in cluster showershape[1] = matchedClus->GetM02(); // long axis showershape[2] = matchedClus->GetM20(); // short axis showershape[3] = matchedClus->GetDispersion(); // dispersion // look for cached value first const AliDetectorPID *detPID=track->GetDetectorPID(); const EDetector detector=kEMCAL; if ( detPID && detPID->HasNumberOfSigmas(detector)){ return detPID->GetNumberOfSigmas(detector, type); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetNumberOfSigmas(detector, type); } // NSigma value really meaningful only for electrons! return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); } } } return -999; } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const { // // Compute PID response of 'detCode' // switch (detCode){ case kDetITS: return ComputeITSProbability(track, nSpecies, p); break; case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break; case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break; case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break; case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break; case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break; case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break; default: return kDetNoSignal; } } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const { // // Compute PID response of 'detCode' // return ComputePIDProbability((EDetCode)(1<GetDetectorPID(); const EDetector detector=kITS; if ( detPID && detPID->HasRawProbabilitiy(detector)){ return detPID->GetRawProbability(detector, p, nSpecies); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetRawProbability(detector, p, nSpecies); } return GetComputeITSProbability(track, nSpecies, p); } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const { // // Compute PID response for the TPC // // look for cached value first const AliDetectorPID *detPID=track->GetDetectorPID(); const EDetector detector=kTPC; if ( detPID && detPID->HasRawProbabilitiy(detector)){ return detPID->GetRawProbability(detector, p, nSpecies); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetRawProbability(detector, p, nSpecies); } return GetComputeTPCProbability(track, nSpecies, p); } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const { // // Compute PID response for the // const AliDetectorPID *detPID=track->GetDetectorPID(); const EDetector detector=kTOF; if ( detPID && detPID->HasRawProbabilitiy(detector)){ return detPID->GetRawProbability(detector, p, nSpecies); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetRawProbability(detector, p, nSpecies); } return GetComputeTOFProbability(track, nSpecies, p); } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const { // // Compute PID response for the // const AliDetectorPID *detPID=track->GetDetectorPID(); const EDetector detector=kTRD; // chacke only for the default method (1d at the moment) if (PIDmethod!=AliTRDPIDResponse::kLQ1D) return GetComputeTRDProbability(track, nSpecies, p, PIDmethod); if ( detPID && detPID->HasRawProbabilitiy(detector)){ return detPID->GetRawProbability(detector, p, nSpecies); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetRawProbability(detector, p, nSpecies); } return GetComputeTRDProbability(track, nSpecies, p); } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const { // // Compute PID response for the EMCAL // const AliDetectorPID *detPID=track->GetDetectorPID(); const EDetector detector=kEMCAL; if ( detPID && detPID->HasRawProbabilitiy(detector)){ return detPID->GetRawProbability(detector, p, nSpecies); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetRawProbability(detector, p, nSpecies); } return GetComputeEMCALProbability(track, nSpecies, p); } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const { // // Compute PID response for the PHOS // // look for cached value first // if (track->GetDetectorPID()){ // return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies); // } // set flat distribution (no decision) for (Int_t j=0; jGetDetectorPID(); const EDetector detector=kHMPID; if ( detPID && detPID->HasRawProbabilitiy(detector)){ return detPID->GetRawProbability(detector, p, nSpecies); } else if (fCachePID) { FillTrackDetectorPID(track, detector); detPID=track->GetDetectorPID(); return detPID->GetRawProbability(detector, p, nSpecies); } return GetComputeHMPIDProbability(track, nSpecies, p); } //______________________________________________________________________________ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run) { // // Apply settings for the current event // fRecoPass=pass; fCurrentEvent=NULL; if (!event) return; fCurrentEvent=event; if (run>0) fRun=run; else fRun=event->GetRunNumber(); if (fRun!=fOldRun){ ExecNewRun(); fOldRun=fRun; } //TPC resolution parametrisation PbPb if ( fResolutionCorrection ){ Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event)); fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04); } //TOF resolution SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod()); // Get and set centrality AliCentrality *centrality = event->GetCentrality(); if(centrality){ fCurrCentrality = centrality->GetCentralityPercentile("V0M"); } else{ fCurrCentrality = -1; } } //______________________________________________________________________________ void AliPIDResponse::ExecNewRun() { // // Things to Execute upon a new run // SetRecoInfo(); SetITSParametrisation(); SetTPCPidResponseMaster(); SetTPCParametrisation(); SetTRDPidResponseMaster(); InitializeTRDResponse(); SetEMCALPidResponseMaster(); InitializeEMCALResponse(); SetTOFPidResponseMaster(); InitializeTOFResponse(); if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField()); } //______________________________________________________________________________ Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event) { // // Get TPC multiplicity in bins of 150 // const AliVVertex* vertexTPC = event->GetPrimaryVertex(); Double_t tpcMulti=0.; if(vertexTPC){ Double_t vertexContribTPC=vertexTPC->GetNContributors(); tpcMulti=vertexContribTPC/150.; if (tpcMulti>20.) tpcMulti=20.; } return tpcMulti; } //______________________________________________________________________________ void AliPIDResponse::SetRecoInfo() { // // Set reconstruction information // //reset information fLHCperiod=""; fMCperiodTPC=""; fBeamType=""; fBeamType="PP"; TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*"); TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*"); //find the period by run number (UGLY, but not stored in ESD and AOD... ) if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; } else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; } else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; } else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; } else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; } else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; } else if (fRun>=136851&&fRun<=139846) { fLHCperiod="LHC10H"; fMCperiodTPC="LHC10H8"; if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10"; fBeamType="PBPB"; } else if (fRun>=139847&&fRun<=146974) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; } //TODO: periods 11B (146975-150721), 11C (150722-155837) are not yet treated assume 11d for the moment else if (fRun>=146975&&fRun<=155837) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; } else if (fRun>=155838&&fRun<=159649) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; } // also for 11e (159650-162750),f(162751-165771) use 11d else if (fRun>=159650&&fRun<=162750) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; } else if (fRun>=162751&&fRun<=165771) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; } else if (fRun>=165772 && fRun<=170718) { fLHCperiod="LHC11H"; fMCperiodTPC="LHC11A10"; fBeamType="PBPB"; if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17"; } if (fRun>=170719 && fRun<=177311) { fLHCperiod="LHC12A"; fBeamType="PP"; /*fMCperiodTPC="";*/ } // for the moment use LHC12b parameters up to LHC12e if (fRun>=177312 /*&& fRun<=179356*/) { fLHCperiod="LHC12B"; fBeamType="PP"; /*fMCperiodTPC="";*/ } // if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; /*fMCperiodTPC="";*/ } // if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; /*fMCperiodTPC="";*/ } // if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; /*fMCperiodTPC="";*/ } // if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; /*fMCperiodTPC="";*/ } // if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; /*fMCperiodTPC="";*/ } // if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ } // for the moment use 12g parametrisation for all full gain runs (LHC12f+) if (fRun >= 186636 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ } //exception new pp MC productions from 2011 if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2"; // exception for 11f1 if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1"; } //______________________________________________________________________________ void AliPIDResponse::SetITSParametrisation() { // // Set the ITS parametrisation // } //______________________________________________________________________________ void AliPIDResponse::SetTPCPidResponseMaster() { // // Load the TPC pid response functions from the OADB // Load the TPC voltage maps from OADB // //don't load twice for the moment if (fArrPidResponseMaster) return; //reset the PID response functions delete fArrPidResponseMaster; fArrPidResponseMaster=NULL; TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data())); TFile *f=NULL; if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse; TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data())); f=TFile::Open(fileNamePIDresponse.Data()); if (f && f->IsOpen() && !f->IsZombie()){ fArrPidResponseMaster=dynamic_cast(f->Get("TPCPIDResponse")); } delete f; TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data())); f=TFile::Open(fileNameVoltageMaps.Data()); if (f && f->IsOpen() && !f->IsZombie()){ fOADBvoltageMaps=dynamic_cast(f->Get("TPCvoltageSettings")); } delete f; if (!fArrPidResponseMaster){ AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data())); return; } fArrPidResponseMaster->SetOwner(); if (!fOADBvoltageMaps) { AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data())); } fArrPidResponseMaster->SetOwner(); } //______________________________________________________________________________ void AliPIDResponse::SetTPCParametrisation() { // // Change BB parametrisation for current run // // //reset old splines // fTPCResponse.ResetSplines(); if (fLHCperiod.IsNull()) { AliError("No period set, not changing parametrisation"); return; } // // Set default parametrisations for data and MC // //data type TString datatype="DATA"; //in case of mc fRecoPass is per default 1 if (fIsMC) { if(!fTuneMConData) datatype="MC"; fRecoPass=1; } // period TString period=fLHCperiod; if (fIsMC && !fTuneMConData) period=fMCperiodTPC; AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data())); Bool_t found=kFALSE; // //set the new PID splines // if (fArrPidResponseMaster){ Int_t recopass = fRecoPass; if(fTuneMConData) recopass = fRecoPassUser; //for MC don't use period information //if (fIsMC) period="[A-Z0-9]*"; //for MC use MC period information //pattern for the default entry (valid for all particles) TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data())); //find particle id ang gain scenario for (Int_t igainScenario=0; igainScenarioGetEntriesFast();++iresp) { TObject *responseFunction=fArrPidResponseMaster->At(iresp); if (responseFunction==NULL) continue; TString responseName=responseFunction->GetName(); if (!reg.MatchB(responseName)) continue; TObjArray *arr=reg.MatchS(responseName); if (!arr) continue; TObject* tmp=NULL; tmp=arr->At(1); if (!tmp) continue; TString particleName=tmp->GetName(); tmp=arr->At(3); if (!tmp) continue; TString gainScenarioName=tmp->GetName(); delete arr; if (particleName.IsNull()) continue; if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction; else { for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec) { TString particle=AliPID::ParticleName(ispec); particle.ToUpper(); //std::cout<GetName())); found=kTRUE; // overwrite default with proton spline (for light nuclei) if (ispec==AliPID::kProton) grAll=responseFunction; break; } } } } if (grAll) { for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec) { if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec, (AliTPCPIDResponse::ETPCgainScenario)igainScenario)) { fTPCResponse.SetResponseFunction( grAll, (AliPID::EParticleType)ispec, (AliTPCPIDResponse::ETPCgainScenario)igainScenario ); fTPCResponse.SetUseDatabase(kTRUE); AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName())); found=kTRUE; } } } } } else AliInfo("no fArrPidResponseMaster"); if (!found){ AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data())); } // // Setup resolution parametrisation // //default fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04); if (fRun>=122195){ fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02); } if (fRun>=186636){ // if (fRun>=188356){ fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05); } if (fArrPidResponseMaster) fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data())); if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName())); //read in the voltage map TVectorF* gsm = 0x0; if (fOADBvoltageMaps) gsm=dynamic_cast(fOADBvoltageMaps->GetObject(fRun)); if (gsm) { fTPCResponse.SetVoltageMap(*gsm); TString vals; AliInfo(Form("Reading the voltage map for run %d\n",fRun)); vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);} AliInfo(vals.Data()); vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);} AliInfo(vals.Data()); vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);} AliInfo(vals.Data()); vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);} AliInfo(vals.Data()); } else AliInfo("no voltage map, ideal default assumed"); } //______________________________________________________________________________ void AliPIDResponse::SetTRDPidResponseMaster() { // // Load the TRD pid params and references from the OADB // if(fTRDPIDResponseObject) return; AliOADBContainer contParams("contParams"); Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject"); if(statusResponse){ AliError("Failed initializing PID Response Object from OADB"); } else { AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data())); fTRDPIDResponseObject = dynamic_cast(contParams.GetObject(fRun)); if(!fTRDPIDResponseObject){ AliError(Form("TRD Response not found in run %d", fRun)); } } } //______________________________________________________________________________ void AliPIDResponse::InitializeTRDResponse(){ // // Set PID Params and references to the TRD PID response // fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject); } //______________________________________________________________________________ void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{ if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){ // backward compatibility for setting with 8 slices TRDslicesForPID[0] = 0; TRDslicesForPID[1] = 7; } else{ if(method==AliTRDPIDResponse::kLQ1D){ TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx TRDslicesForPID[1] = 0; } if(method==AliTRDPIDResponse::kLQ2D){ TRDslicesForPID[0] = 1; TRDslicesForPID[1] = 7; } } AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1])); } //______________________________________________________________________________ void AliPIDResponse::SetTOFPidResponseMaster() { // // Load the TOF pid params from the OADB // if (fTOFPIDParams) delete fTOFPIDParams; fTOFPIDParams=NULL; TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data())); if (oadbf && oadbf->IsOpen()) { AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data())); AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb"); if (oadbc) fTOFPIDParams = dynamic_cast(oadbc->GetObject(fRun,"TOFparams")); oadbf->Close(); delete oadbc; } delete oadbf; if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved"); } //______________________________________________________________________________ void AliPIDResponse::InitializeTOFResponse(){ // // Set PID Params to the TOF PID response // AliInfo("TOF PID Params loaded from OADB"); AliInfo(Form(" TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution())); AliInfo(Form(" StartTime method %d",fTOFPIDParams->GetStartTimeMethod())); AliInfo(Form(" TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f", fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3))); for (Int_t i=0;i<4;i++) { fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i)); } fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution()); AliInfo("TZERO resolution loaded from ESDrun/AODheader"); Float_t t0Spread[4]; for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i); AliInfo(Form(" TZERO spreads from data: (A+C)/2 %f A %f C %f (A'-C')/2: %f",t0Spread[0],t0Spread[1],t0Spread[2],t0Spread[3])); Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3]; Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3]; if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) { fResT0AC=t0Spread[3]; fResT0A=TMath::Sqrt(a); fResT0C=TMath::Sqrt(c); } else { AliInfo(" TZERO spreads not present or inconsistent, loading default"); fResT0A=75.; fResT0C=65.; fResT0AC=55.; } AliInfo(Form(" TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC)); } //______________________________________________________________________________ Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const { // // Check whether track is identified as electron under a given electron efficiency hypothesis // Double_t probs[AliPID::kSPECIES]; ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod); Int_t ntracklets = vtrack->GetTRDntrackletsPID(); // Take mean of the TRD momenta in the given tracklets Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes]; Int_t nmomenta = 0; for(Int_t iPl=0;iPlGetTRDmomentum(iPl) > 0.){ trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); } } p = TMath::Mean(nmomenta, trdmomenta); return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod); } //______________________________________________________________________________ void AliPIDResponse::SetEMCALPidResponseMaster() { // // Load the EMCAL pid response functions from the OADB // TObjArray* fEMCALPIDParamsRun = NULL; TObjArray* fEMCALPIDParamsPass = NULL; if(fEMCALPIDParams) return; AliOADBContainer contParams("contParams"); Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams"); if(statusPars){ AliError("Failed initializing PID Params from OADB"); } else { AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data())); fEMCALPIDParamsRun = dynamic_cast(contParams.GetObject(fRun)); if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass))); if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles"))); if(!fEMCALPIDParams){ AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass)); AliInfo("Will take the standard LHC11d instead ..."); fEMCALPIDParamsRun = dynamic_cast(contParams.GetObject(156477)); if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast(fEMCALPIDParamsRun->FindObject(Form("pass%d",1))); if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles"))); if(!fEMCALPIDParams){ AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data())); } } } } //______________________________________________________________________________ void AliPIDResponse::InitializeEMCALResponse(){ // // Set PID Params to the EMCAL PID response // fEMCALResponse.SetPIDParams(fEMCALPIDParams); } //______________________________________________________________________________ void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const { // // create detector PID information and setup the transient pointer in the track // // check if detector number is inside accepted range if (detector == kNdetectors) return; // get detector pid AliDetectorPID *detPID=const_cast(track->GetDetectorPID()); if (!detPID) { detPID=new AliDetectorPID; (const_cast(track))->SetDetectorPID(detPID); } //check if values exist if (detPID->HasRawProbabilitiy(detector) && detPID->HasNumberOfSigmas(detector)) return; //TODO: which particles to include? See also the loops below... Double_t values[AliPID::kSPECIESC]={0}; //nsigmas for (Int_t ipart=0; ipartSetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC); //probabilities EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values); detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status); } //______________________________________________________________________________ void AliPIDResponse::FillTrackDetectorPID() { // // create detector PID information and setup the transient pointer in the track // if (!fCurrentEvent) return; for (Int_t itrack=0; itrackGetNumberOfTracks(); ++itrack){ AliVTrack *track=dynamic_cast(fCurrentEvent->GetTrack(itrack)); if (!track) continue; for (Int_t idet=0; idetGetEventTimeSpread(); if(t0spread < 10) t0spread = 80; // T0 from TOF algorithm Bool_t flagT0TOF=kFALSE; Bool_t flagT0T0=kFALSE; Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()]; Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()]; Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()]; // T0-TOF arrays Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()]; Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()]; for(Int_t i=0;iGetT0TOF()){ // check if T0 detector information is available flagT0T0=kTRUE; } AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader(); if (tofHeader) { // read global info and T0-TOF fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution()); t0spread = tofHeader->GetT0spread(); // read t0 sprad if(t0spread < 10) t0spread = 80; flagT0TOF=kTRUE; for(Int_t i=0;iGetDefaultEventTimeVal(); startTimeRes[i]=tofHeader->GetDefaultEventTimeRes(); if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread; } TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues(); TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues(); TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes(); for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins Int_t icurrent = (Int_t)ibin->GetAt(j); startTime[icurrent]=t0Bin->GetAt(j); startTimeRes[icurrent]=t0ResBin->GetAt(j); if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread; } } // for cut of 3 sigma on t0 spread Float_t t0cut = 3 * t0spread; if(t0cut < 500) t0cut = 500; if(option == kFILL_T0){ // T0-FILL is used for(Int_t i=0;iGetT0TOF()[1]; t0C= vevent->GetT0TOF()[2]; // t0AC= vevent->GetT0TOF()[0]; t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C; resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C); t0AC /= resT0AC*resT0AC; } Float_t t0t0Best = 0; Float_t t0t0BestRes = 9999; Int_t t0used=0; if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){ t0t0Best = t0AC; t0t0BestRes = resT0AC; t0used=6; } else if(TMath::Abs(t0C) < t0cut){ t0t0Best = t0C; t0t0BestRes = resT0C; t0used=4; } else if(TMath::Abs(t0A) < t0cut){ t0t0Best = t0A; t0t0BestRes = resT0A; t0used=2; } if(flagT0TOF){ // if T0-TOF info is available for(Int_t i=0;iGetT0TOF()[1]; t0C= vevent->GetT0TOF()[2]; // t0AC= vevent->GetT0TOF()[0]; t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C; resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C); t0AC /= resT0AC*resT0AC; } if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){ for(Int_t i=0;iGetITSsignal(); if (dEdx<=0) return -999.; UChar_t clumap=track->GetITSClusterMap(); Int_t nPointsForPid=0; for(Int_t i=2; i<6; i++){ if(clumap&(1<P(); //check for ITS standalone tracks Bool_t isSA=kTRUE; if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE; //TODO: in case of the electron, use the SA parametrisation, // this needs to be changed if ITS provides a parametrisation // for electrons also for ITS+TPC tracks return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron)); } //______________________________________________________________________________ Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const { // // Calculate the number of sigmas in the TPC // AliVTrack *track=(AliVTrack*)vtrack; Double_t mom = track->GetTPCmomentum(); Double_t sig = track->GetTPCsignal(); if(fTuneMConData) sig = this->GetTPCsignalTunedOnData(track); UInt_t sigN = track->GetTPCsignalN(); Double_t nSigma = -999.; if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type); return nSigma; } //______________________________________________________________________________ Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const { // // Calculate the number of sigmas in the EMCAL // AliVTrack *track=(AliVTrack*)vtrack; AliVCluster *matchedClus = NULL; Double_t mom = -1.; Double_t pt = -1.; Double_t EovP = -1.; Double_t fClsE = -1.; Int_t nMatchClus = -1; Int_t charge = 0; // Track matching nMatchClus = track->GetEMCALcluster(); if(nMatchClus > -1){ mom = track->P(); pt = track->Pt(); charge = track->Charge(); matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus); if(matchedClus){ // matched cluster is EMCAL if(matchedClus->IsEMCAL()){ fClsE = matchedClus->E(); EovP = fClsE/mom; // NSigma value really meaningful only for electrons! return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); } } } return -999; } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const { // // Compute PID response of 'detCode' // switch (detCode){ case kITS: return GetComputeITSProbability(track, nSpecies, p); break; case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break; case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break; case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break; case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break; case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break; case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break; default: return kDetNoSignal; } } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const { // // Compute PID response for the ITS // // look for cached value first // only the non SA tracks are cached if (track->GetDetectorPID()){ return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies); } // set flat distribution (no decision) for (Int_t j=0; jGetStatus()&AliVTrack::kITSin)==0 && (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal; //check for ITS standalone tracks Bool_t isSA=kTRUE; if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE; Double_t mom=track->P(); Double_t dedx=track->GetITSsignal(); Double_t momITS=mom; UChar_t clumap=track->GetITSClusterMap(); Int_t nPointsForPid=0; for(Int_t i=2; i<6; i++){ if(clumap&(1<ResetStatus(AliVTrack::kITSpid); return kDetNoSignal; } Bool_t mismatch=kTRUE/*, heavy=kTRUE*/; for (Int_t j=0; j fRange*sigma) { p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; } else { p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; mismatch=kFALSE; } // Check for particles heavier than (AliPID::kSPECIES - 1) // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE; } if (mismatch){ for (Int_t j=0; jGetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal; Double_t mom = track->GetTPCmomentum(); Double_t dedx=track->GetTPCsignal(); Bool_t mismatch=kTRUE/*, heavy=kTRUE*/; if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track); for (Int_t j=0; jGetTPCsignalN(),type); if (TMath::Abs(dedx-bethe) > fRange*sigma) { p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; } else { p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; mismatch=kFALSE; } } if (mismatch){ for (Int_t j=0; jGetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal; if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal; Bool_t mismatch = kTRUE/*, heavy = kTRUE*/; for (Int_t j=0; jP(),expTime,AliPID::ParticleMassZ(type)); if (TMath::Abs(nsigmas) > (fRange+2)) { if(nsigmas < fTOFtail) p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig; else p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig; } else{ if(nsigmas < fTOFtail) p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig; else p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig; } if (TMath::Abs(nsigmas)<5.){ Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type); if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE; } } if (mismatch){ return kDetMismatch; } // TODO: Light nuclei return kDetPidOk; } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const { // // Compute PID response for the // UInt_t TRDslicesForPID[2]; SetTRDSlices(TRDslicesForPID,PIDmethod); // set flat distribution (no decision) for (Int_t j=0; jGetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal; Float_t mom[6]={0.}; Double_t dedx[48]={0.}; // Allocate space for the maximum number of TRD slices Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1; AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", TRDslicesForPID[0], TRDslicesForPID[1], nslices)); for(UInt_t ilayer = 0; ilayer < 6; ilayer++){ mom[ilayer] = track->GetTRDmomentum(ilayer); for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){ dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice); } } fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod); return kDetPidOk; } //______________________________________________________________________________ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const { // // Compute PID response for the EMCAL // for (Int_t j=0; jGetEMCALcluster(); if(nMatchClus > -1){ mom = track->P(); pt = track->Pt(); charge = track->Charge(); matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus); if(matchedClus){ // matched cluster is EMCAL if(matchedClus->IsEMCAL()){ fClsE = matchedClus->E(); EovP = fClsE/mom; // compute the probabilities if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){ // in case everything is OK return kDetPidOk; } } } } // in all other cases set flat distribution (no decision) for (Int_t j=0; jGetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal; track->GetHMPIDpid(p); return kDetPidOk; }