From 87da02055b83ef70791da14c3bdb7b62d5b0d159 Mon Sep 17 00:00:00 2001 From: morsch Date: Tue, 16 Apr 2013 06:32:40 +0000 Subject: [PATCH] Attached you can find TPC and EMCAL PID response files for the OADB and a larger patch: o preparation of multiplicity treatment for TPC PID o updated TPC splines for (nearly) all periods o first splines for LHC13b2_fix_1 (no eta maps, yet, to be updated...) o update in multiplicity treatment for EMCAL PID Jens Wiechula --- ANALYSIS/AliAnalysisTaskPIDResponse.cxx | 14 +- ANALYSIS/AliAnalysisTaskPIDResponse.h | 6 +- ANALYSIS/macros/AddTaskPIDResponse.C | 12 +- STEER/AOD/AliAODpidUtil.cxx | 6 +- STEER/ESD/AliESDpid.cxx | 6 +- STEER/STEERBase/AliEMCALPIDResponse.cxx | 55 +++- STEER/STEERBase/AliEMCALPIDResponse.h | 5 +- STEER/STEERBase/AliPIDResponse.cxx | 164 ++++++++-- STEER/STEERBase/AliPIDResponse.h | 18 +- STEER/STEERBase/AliTPCPIDResponse.cxx | 418 +++++++++++++++++++++--- STEER/STEERBase/AliTPCPIDResponse.h | 64 +++- 11 files changed, 671 insertions(+), 97 deletions(-) diff --git a/ANALYSIS/AliAnalysisTaskPIDResponse.cxx b/ANALYSIS/AliAnalysisTaskPIDResponse.cxx index 816e5c4135f..49af37ba7f9 100644 --- a/ANALYSIS/AliAnalysisTaskPIDResponse.cxx +++ b/ANALYSIS/AliAnalysisTaskPIDResponse.cxx @@ -46,7 +46,8 @@ fRecoPass(0), fIsTunedOnData(kFALSE), fTunedOnDataMask(0), fRecoPassTuned(0), -fUseTPCEtaCorrection(kFALSE)//TODO: In future, default kTRUE +fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE +fUseTPCMultiplicityCorrection(kFALSE)//TODO: In future, default kTRUE { // // Dummy constructor @@ -67,7 +68,8 @@ fRecoPass(0), fIsTunedOnData(kFALSE), fTunedOnDataMask(0), fRecoPassTuned(0), -fUseTPCEtaCorrection(kFALSE)//TODO: In future, default kTRUE +fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE +fUseTPCMultiplicityCorrection(kFALSE)//TODO: In future, default kTRUE { // // Default constructor @@ -120,7 +122,7 @@ void AliAnalysisTaskPIDResponse::UserCreateOutputObjects() } } delete arr; - } + } } //______________________________________________________________________________ @@ -135,9 +137,11 @@ void AliAnalysisTaskPIDResponse::UserExec(Option_t */*option*/) if (fRun!=fOldRun){ SetRecoInfo(); fOldRun=fRun; + + fPIDResponse->SetUseTPCEtaCorrection(fUseTPCEtaCorrection); + fPIDResponse->SetUseTPCMultiplicityCorrection(fUseTPCMultiplicityCorrection); } - fPIDResponse->SetUseTPCEtaCorrection(fUseTPCEtaCorrection); fPIDResponse->InitialiseEvent(event,fRecoPass); AliESDpid *pidresp = dynamic_cast(fPIDResponse); if(pidresp && AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){ @@ -178,6 +182,8 @@ void AliAnalysisTaskPIDResponse::SetRecoInfo() fPIDResponse->SetCurrentFile(fileName.Data()); } + fPIDResponse->SetCurrentAliRootRev(prodInfo.GetAlirootSvnVersion()); + if (prodInfo.IsMC() == kTRUE) fIsMC=kTRUE; // protection if user didn't use macro switch if ( (prodInfo.IsMC() == kFALSE) && (fIsMC == kFALSE) ) { // reco pass is needed only for data fRecoPass = prodInfo.GetRecoPass(); diff --git a/ANALYSIS/AliAnalysisTaskPIDResponse.h b/ANALYSIS/AliAnalysisTaskPIDResponse.h index 2e95651b0de..d842ffde2bb 100644 --- a/ANALYSIS/AliAnalysisTaskPIDResponse.h +++ b/ANALYSIS/AliAnalysisTaskPIDResponse.h @@ -47,6 +47,9 @@ public: void SetUseTPCEtaCorrection(Bool_t useTPCEtaCorrection) { fUseTPCEtaCorrection = useTPCEtaCorrection; }; Bool_t UseTPCEtaCorrection() const { return fUseTPCEtaCorrection; }; + + void SetUseTPCMultiplicityCorrection(Bool_t useMultiplicityCorrection = kTRUE) { fUseTPCMultiplicityCorrection = useMultiplicityCorrection; }; + Bool_t UseTPCMultiplicityCorrection() const { return fUseTPCMultiplicityCorrection; }; void SetSpecialDetectorResponse(const char* det) { fSpecialDetResponse=det; } @@ -65,7 +68,8 @@ private: Int_t fTunedOnDataMask; // mask to activate tuning on data on specific detectors Int_t fRecoPassTuned; // Reco pass tuned on data for MC - Bool_t fUseTPCEtaCorrection; // Use TPC eta correction + Bool_t fUseTPCEtaCorrection; // Use TPC eta correction + Bool_t fUseTPCMultiplicityCorrection; // Use TPC multiplicity correction // void SetRecoInfo(); diff --git a/ANALYSIS/macros/AddTaskPIDResponse.C b/ANALYSIS/macros/AddTaskPIDResponse.C index a8fb7a9d869..b407a6fec35 100644 --- a/ANALYSIS/macros/AddTaskPIDResponse.C +++ b/ANALYSIS/macros/AddTaskPIDResponse.C @@ -1,7 +1,8 @@ AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE, Bool_t tuneOnData=kFALSE, Int_t recoPass=2, Bool_t cachePID=kFALSE, TString detResponse="", - Bool_t useTPCEtaCorrection = kTRUE) + Bool_t useTPCEtaCorrection = kTRUE, + Bool_t useTPCMultiplicityCorrection = kFALSE) { // Macro to connect a centrality selection task to an existing analysis manager. AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); @@ -22,10 +23,6 @@ AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE, AliPIDResponseInputHandler *pidResponseIH = new AliPIDResponseInputHandler(); multiInputHandler->AddInputEventHandler(pidResponseIH); - if (autoMCesd && - multiInputHandler->GetFirstInputEventHandler()->IsA()==AliESDInputHandler::Class() && - multiInputHandler->GetFirstMCEventHandler() - ) isMC=kTRUE; pidResponseIH->SetIsMC(isMC); return 0x0; @@ -36,10 +33,6 @@ AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE, printf("PIDResponse: Initialising AliAnalysisTaskPIDResponse\n"); printf("========================================================================================\n"); - if ( autoMCesd && (inputHandler->IsA() == AliESDInputHandler::Class()) ) { - isMC=mgr->GetMCtruthEventHandler()!=0x0; - } - AliAnalysisTaskPIDResponse *pidTask = new AliAnalysisTaskPIDResponse("PIDResponseTask"); // pidTask->SelectCollisionCandidates(AliVEvent::kMB); pidTask->SetIsMC(isMC); @@ -52,6 +45,7 @@ AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE, pidTask->SetCachePID(cachePID); pidTask->SetSpecialDetectorResponse(detResponse); pidTask->SetUseTPCEtaCorrection(useTPCEtaCorrection); + pidTask->SetUseTPCMultiplicityCorrection(useTPCMultiplicityCorrection); mgr->AddTask(pidTask); // AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("PIDResponseQA", diff --git a/STEER/AOD/AliAODpidUtil.cxx b/STEER/AOD/AliAODpidUtil.cxx index f5d31fb1bc5..45cda466436 100644 --- a/STEER/AOD/AliAODpidUtil.cxx +++ b/STEER/AOD/AliAODpidUtil.cxx @@ -91,8 +91,10 @@ Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const { if(kGood){ //TODO maybe introduce different dEdxSources? - Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection()); - Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection()); + Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(), + this->UseTPCMultiplicityCorrection()); + Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(), + this->UseTPCMultiplicityCorrection()); dedx = gRandom->Gaus(bethe,sigma); // if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5; diff --git a/STEER/ESD/AliESDpid.cxx b/STEER/ESD/AliESDpid.cxx index 7c39219f931..79043951de6 100644 --- a/STEER/ESD/AliESDpid.cxx +++ b/STEER/ESD/AliESDpid.cxx @@ -123,8 +123,10 @@ Float_t AliESDpid::GetTPCsignalTunedOnData(const AliVTrack *t) const { if(kGood){ //TODO maybe introduce different dEdxSources? - Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection()); - Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection()); + Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(), + this->UseTPCMultiplicityCorrection()); + Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(), + this->UseTPCMultiplicityCorrection()); dedx = gRandom->Gaus(bethe,sigma); // if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5; } diff --git a/STEER/STEERBase/AliEMCALPIDResponse.cxx b/STEER/STEERBase/AliEMCALPIDResponse.cxx index 7efcb58b4e6..91c06189a99 100644 --- a/STEER/STEERBase/AliEMCALPIDResponse.cxx +++ b/STEER/STEERBase/AliEMCALPIDResponse.cxx @@ -61,6 +61,7 @@ ClassImp(AliEMCALPIDResponse) AliEMCALPIDResponse::AliEMCALPIDResponse(): TObject(), fNorm(NULL), + fCurrCentrality(-1.), fkPIDParams(NULL) { // @@ -74,6 +75,7 @@ AliEMCALPIDResponse::AliEMCALPIDResponse(): AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other): TObject(other), fNorm(other.fNorm), + fCurrCentrality(other.fCurrCentrality), fkPIDParams(other.fkPIDParams) { // @@ -93,6 +95,7 @@ AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& // Make copy TObject::operator=(other); fNorm = other.fNorm; + fCurrCentrality = other.fCurrCentrality; fkPIDParams = other.fkPIDParams; @@ -264,6 +267,12 @@ const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int // 6 = probLow (not used for electrons) // 7 = probHigh (not used for electrons) // + // for PbPb the parametrization is done centrality dependent (marked by TString "Centrality") + // so first the correct centrality bin has to be found + + // **** Centrality bins (hard coded for the moment) + const Int_t nCent = 7; + Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90}; if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL; if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons @@ -271,19 +280,53 @@ const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int TObjArray * particlePar = dynamic_cast(fkPIDParams->At(nParticle)); if(!particlePar) return NULL; - TIter parIter(particlePar); const TVectorD *parameters = NULL; Double_t momMin = 0.; Double_t momMax = 0.; - while((parameters = static_cast(parIter()))){ + // is the centrality dependent parametrization used + TString arrayName = particlePar->GetName(); + + // centrality dependent parametrization + if(arrayName.Contains("Centrality")){ + + for(Int_t iCent = 0; iCent < nCent; iCent++){ - momMin = (*parameters)[0]; - momMax = (*parameters)[1]; + if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){ + + TObjArray * centPar = dynamic_cast(particlePar->At(iCent)); + if(!centPar) return NULL; + + TIter centIter(centPar); + parameters = NULL; + momMin = 0.; + momMax = 0.; + + while((parameters = static_cast(centIter()))){ + + momMin = (*parameters)[0]; + momMax = (*parameters)[1]; + + if( fPt > momMin && fPt < momMax ) return parameters; + + } + } + } + } - if( fPt > momMin && fPt < momMax ) return parameters; + // NO centrality dependent parametrization + else{ - } + TIter parIter(particlePar); + while((parameters = static_cast(parIter()))){ + + momMin = (*parameters)[0]; + momMax = (*parameters)[1]; + + if( fPt > momMin && fPt < momMax ) return parameters; + + } + } AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt)); return parameters; diff --git a/STEER/STEERBase/AliEMCALPIDResponse.h b/STEER/STEERBase/AliEMCALPIDResponse.h index 57c17fb8abb..ee5b9300df6 100644 --- a/STEER/STEERBase/AliEMCALPIDResponse.h +++ b/STEER/STEERBase/AliEMCALPIDResponse.h @@ -35,6 +35,7 @@ public : //Setters void SetPIDParams(const TObjArray * params) { fkPIDParams = params; } + void SetCentrality(Float_t currentCentrality) { fCurrCentrality = currentCentrality;} // EMCAL probability @@ -46,11 +47,13 @@ private: TF1 *fNorm; // Gauss function for normalizing NON electron probabilities + Double_t fCurrCentrality; // current (in the current event) centrality percentile + const TObjArray *fkPIDParams; // PID Params const TVectorD* GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const; - ClassDef(AliEMCALPIDResponse, 1) + ClassDef(AliEMCALPIDResponse, 2) }; #endif // #ifdef AliEMCALPIDResponse_cxx diff --git a/STEER/STEERBase/AliPIDResponse.cxx b/STEER/STEERBase/AliPIDResponse.cxx index 4806cc91851..68756167a4b 100644 --- a/STEER/STEERBase/AliPIDResponse.cxx +++ b/STEER/STEERBase/AliPIDResponse.cxx @@ -71,6 +71,7 @@ fLHCperiod(), fMCperiodTPC(), fMCperiodUser(), fCurrentFile(), +fCurrentAliRootRev(-1), fRecoPass(0), fRecoPassUser(-1), fRun(-1), @@ -81,7 +82,8 @@ fResT0AC(55.), fArrPidResponseMaster(NULL), fResolutionCorrection(NULL), fOADBvoltageMaps(NULL), -fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE +fUseTPCEtaCorrection(kFALSE), +fUseTPCMultiplicityCorrection(kFALSE), fTRDPIDResponseObject(NULL), fTOFtail(1.1), fTOFPIDParams(NULL), @@ -132,6 +134,7 @@ fLHCperiod(), fMCperiodTPC(), fMCperiodUser(other.fMCperiodUser), fCurrentFile(), +fCurrentAliRootRev(other.fCurrentAliRootRev), fRecoPass(0), fRecoPassUser(other.fRecoPassUser), fRun(-1), @@ -143,6 +146,7 @@ fArrPidResponseMaster(NULL), fResolutionCorrection(NULL), fOADBvoltageMaps(NULL), fUseTPCEtaCorrection(other.fUseTPCEtaCorrection), +fUseTPCMultiplicityCorrection(other.fUseTPCMultiplicityCorrection), fTRDPIDResponseObject(NULL), fTOFtail(1.1), fTOFPIDParams(NULL), @@ -184,6 +188,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other) fMCperiodTPC=""; fMCperiodUser=other.fMCperiodUser; fCurrentFile=""; + fCurrentAliRootRev=other.fCurrentAliRootRev; fRecoPass=0; fRecoPassUser=other.fRecoPassUser; fRun=-1; @@ -195,12 +200,14 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other) fResolutionCorrection=NULL; fOADBvoltageMaps=NULL; fUseTPCEtaCorrection=other.fUseTPCEtaCorrection; + fUseTPCMultiplicityCorrection=other.fUseTPCMultiplicityCorrection; fTRDPIDResponseObject=NULL; fEMCALPIDParams=NULL; fTOFtail=1.1; fTOFPIDParams=NULL; fHMPIDPIDParams=NULL; fCurrentEvent=other.fCurrentEvent; + } return *this; } @@ -271,8 +278,7 @@ Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack, //get number of sigmas according the selected TPC gain configuration scenario const AliVTrack *track=static_cast(vtrack); -// return 0.; - Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection); + Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection); return nSigma; } @@ -544,6 +550,20 @@ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run) fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04); } + // Set up TPC multiplicity for PbPb + //TODO Will NOT give the desired number for AODs -> Needs new variable/function in future. + // Fatal, if AOD event and correction enabled + //printf("DETECTED class: %s (%d)\n\n\n\n", event->IsA()->GetName(), fUseTPCMultiplicityCorrection);//TODO + if (fUseTPCMultiplicityCorrection && strcmp(event->IsA()->GetName(), "AliESDEvent") != 0) { + AliFatal("TPC multiplicity correction is enabled, but will NOT work for AOD events, only for ESD => Disabled multiplicity correction!"); + fUseTPCMultiplicityCorrection = kFALSE; + } + + if (fUseTPCMultiplicityCorrection) + fTPCResponse.SetCurrentEventMultiplicity(event->GetNumberOfTracks()); + else + fTPCResponse.SetCurrentEventMultiplicity(0); + //TOF resolution SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod()); @@ -556,6 +576,10 @@ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run) else{ fCurrCentrality = -1; } + + // Set centrality percentile for EMCAL + fEMCALResponse.SetCentrality(fCurrCentrality); + } //______________________________________________________________________________ @@ -661,8 +685,18 @@ void AliPIDResponse::SetRecoInfo() // 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 && fRun < 194480) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; } - if (fRun >= 194480) { fLHCperiod="LHC13B"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; } + if (fRun >= 186636 && fRun < 194480) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; } + + // New parametrisation for 2013 pPb runs + if (fRun >= 194480) { + fLHCperiod="LHC13B"; + fBeamType="PPB"; + + if (fCurrentAliRootRev >= 61605) + fMCperiodTPC="LHC13B2_FIX"; + else + fMCperiodTPC="LHC12G"; + } //exception new pp MC productions from 2011 if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; } @@ -881,13 +915,13 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod; if (fIsMC) { - if (!fTuneMConData) { + if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) { period=fMCperiodTPC; dataType="MC"; } fRecoPass = 1; - if (!fTuneMConData && fMCperiodTPC.IsNull()) { + if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) && fMCperiodTPC.IsNull()) { AliFatal("MC detected, but no MC period set -> Not changing eta maps!"); return; } @@ -912,13 +946,14 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass)); if (statusCont) { AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction"); + fUseTPCEtaCorrection = kFALSE; } else { AliInfo(Form("Loading TPC eta correction map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data())); TH2D* etaMap = 0x0; - if (fIsMC && !fTuneMConData) { + if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) { TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass); etaMap = dynamic_cast(etaMapsCont.GetDefaultObject(searchMap.Data())); if (!etaMap) { @@ -933,6 +968,7 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac if (!etaMap) { AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun)); + fUseTPCEtaCorrection = kFALSE; } else { TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY); @@ -941,6 +977,7 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) { AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun)); fTPCResponse.SetEtaCorrMap(0x0); + fUseTPCEtaCorrection = kFALSE; } else { AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s", @@ -951,10 +988,17 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac } else { AliError(Form("Failed to set TPC eta correction map for run %d (map was loaded, but couldn't be refined) -> Disabled eta correction!!!", fRun)); + fUseTPCEtaCorrection = kFALSE; } } } + // If there was some problem loading the eta maps, it makes no sense to load the sigma maps (that require eta corrected data) + if (fUseTPCEtaCorrection == kFALSE) { + AliError("Failed to load TPC eta correction map required by sigma maps -> Using old parametrisation for sigma"); + return; + } + // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta)) AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass)); @@ -968,7 +1012,7 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac TObjArray* etaSigmaPars = 0x0; - if (fIsMC && !fTuneMConData) { + if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) { TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass); etaSigmaPars = dynamic_cast(etaSigmaMapsCont.GetDefaultObject(searchMap.Data())); if (!etaSigmaPars) { @@ -1091,16 +1135,16 @@ void AliPIDResponse::SetTPCParametrisation() TString datatype="DATA"; //in case of mc fRecoPass is per default 1 if (fIsMC) { - if(!fTuneMConData) datatype="MC"; + if(!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) datatype="MC"; fRecoPass=1; } // period TString period=fLHCperiod; - if (fIsMC && !fTuneMConData) period=fMCperiodTPC; + if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) period=fMCperiodTPC; Int_t recopass = fRecoPass; - if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) recopass = fRecoPassUser; + if(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) recopass = fRecoPassUser; AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data())); Bool_t found=kFALSE; @@ -1223,8 +1267,86 @@ void AliPIDResponse::SetTPCParametrisation() AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data())); } + + // + // Setup multiplicity correction + // + if (fUseTPCMultiplicityCorrection && (fBeamType.CompareTo("PBPB") == 0 || fBeamType.CompareTo("AA") == 0)) { + AliInfo("Multiplicity correction enabled!"); + + //TODO After testing, load parameters from outside + if (period.Contains("LHC11A10")) {//LHC11A10A + AliInfo("Using multiplicity correction parameters for 11a10!"); + fTPCResponse.SetParameterMultiplicityCorrection(0, 6.90133e-06); + fTPCResponse.SetParameterMultiplicityCorrection(1, -1.22123e-03); + fTPCResponse.SetParameterMultiplicityCorrection(2, 1.80220e-02); + fTPCResponse.SetParameterMultiplicityCorrection(3, 0.1); + fTPCResponse.SetParameterMultiplicityCorrection(4, 6.45306e-03); + + fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -2.85505e-07); + fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, -1.31911e-06); + fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5); + + fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -4.29665e-05); + fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 1.37023e-02); + fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -6.36337e-01); + fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.13479e-02); + } + else if (period.Contains("LHC10H") && recopass == 2) { + AliInfo("Using multiplicity correction parameters for 10h.pass2!"); + fTPCResponse.SetParameterMultiplicityCorrection(0, 3.21636e-07); + fTPCResponse.SetParameterMultiplicityCorrection(1, -6.65876e-04); + fTPCResponse.SetParameterMultiplicityCorrection(2, 1.28786e-03); + fTPCResponse.SetParameterMultiplicityCorrection(3, 1.47677e-02); + fTPCResponse.SetParameterMultiplicityCorrection(4, 0); + + fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, 7.23591e-08); + fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 2.7469e-06); + fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5); + + fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -1.22590e-05); + fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 6.88888e-03); + fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -3.20788e-01); + fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.07345e-02); + } + else { + AliError(Form("Multiplicity correction is enabled, but no multiplicity correction parameters have been found for period %s.pass%d -> Mulitplicity correction DISABLED!", period.Data(), recopass)); + fUseTPCMultiplicityCorrection = kFALSE; + fTPCResponse.ResetMultiplicityCorrectionFunctions(); + } + } + else { + // Just set parameters such that overall correction factor is 1, i.e. no correction. + // This is just a reasonable choice for the parameters for safety reasons. Disabling + // the multiplicity correction will anyhow skip the calculation of the corresponding + // correction factor inside THIS class. Nevertheless, experts can access the TPCPIDResponse + // directly and use it for calculations - which should still give valid results, even if + // the multiplicity correction is explicitely enabled in such expert calls. + + AliInfo(Form("Multiplicity correction %sdisabled (%s)!", fUseTPCMultiplicityCorrection ? "automatically " : "", + fUseTPCMultiplicityCorrection ? "no PbPb or AA" : "requested by user")); + + fUseTPCMultiplicityCorrection = kFALSE; + fTPCResponse.ResetMultiplicityCorrectionFunctions(); + } + + /* + //TODO NOW start + for (Int_t i = 0; i <= 4 + 1; i++) { + printf("parMultCorr: %d, %e\n", i, fTPCResponse.GetMultiplicityCorrectionFunction()->GetParameter(i)); + } + for (Int_t j = 0; j <= 2 + 1; j++) { + printf("parMultCorrTanTheta: %d, %e\n", j, fTPCResponse.GetMultiplicityCorrectionFunctionTanTheta()->GetParameter(j)); + } + for (Int_t j = 0; j <= 3 + 1; j++) { + printf("parMultSigmaCorr: %d, %e\n", j, fTPCResponse.GetMultiplicitySigmaCorrectionFunction()->GetParameter(j)); + } + + //TODO NOW end + */ + // - // Setup resolution parametrisation + // Setup old resolution parametrisation // //default @@ -1805,9 +1927,10 @@ Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID: // the following call is needed in order to fill the transient data member // fTPCsignalTuned which is used in the TPCPIDResponse to judge // if using tuned on data - if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) this->GetTPCsignalTunedOnData(track); + if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) + this->GetTPCsignalTunedOnData(track); - return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection); + return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection); } //______________________________________________________________________________ @@ -1886,10 +2009,10 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTPC(const AliVPartic // the following call is needed in order to fill the transient data member // fTPCsignalTuned which is used in the TPCPIDResponse to judge // if using tuned on data - if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) + if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) this->GetTPCsignalTunedOnData(track); - val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, ratio); + val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, ratio); return GetTPCPIDStatus(track); } @@ -1902,6 +2025,7 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTOF(const AliVPartic // AliVTrack *track=(AliVTrack*)vtrack; val=GetSignalDeltaTOFold(track, type, ratio); + return GetTOFPIDStatus(track); } @@ -2005,7 +2129,7 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability (const A Double_t dedx=track->GetTPCsignal(); Bool_t mismatch=kTRUE/*, heavy=kTRUE*/; - if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) dedx = this->GetTPCsignalTunedOnData(track); + if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) dedx = this->GetTPCsignalTunedOnData(track); Double_t bethe = 0.; Double_t sigma = 0.; @@ -2013,8 +2137,8 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability (const A for (Int_t j=0; j fRange*sigma) { p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; diff --git a/STEER/STEERBase/AliPIDResponse.h b/STEER/STEERBase/AliPIDResponse.h index 97cbe5d85da..fa00d13fdb4 100644 --- a/STEER/STEERBase/AliPIDResponse.h +++ b/STEER/STEERBase/AliPIDResponse.h @@ -130,6 +130,9 @@ public: void InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run=-1); void SetCurrentFile(const char* file) { fCurrentFile=file; } + + void SetCurrentAliRootRev(Int_t alirootRev) { fCurrentAliRootRev = alirootRev; } + Int_t GetCurrentAliRootRev() const { return fCurrentAliRootRev; } // cache PID in the track void SetCachePID(Bool_t cache) { fCachePID=cache; } @@ -150,6 +153,9 @@ public: void SetUseTPCEtaCorrection(Bool_t useEtaCorrection = kTRUE) { fUseTPCEtaCorrection = useEtaCorrection; }; Bool_t UseTPCEtaCorrection() const { return fUseTPCEtaCorrection; }; + void SetUseTPCMultiplicityCorrection(Bool_t useMultiplicityCorrection = kTRUE) { fUseTPCMultiplicityCorrection = useMultiplicityCorrection; }; + Bool_t UseTPCMultiplicityCorrection() const { return fUseTPCMultiplicityCorrection; }; + // TOF setting void SetTOFtail(Float_t tail=1.1){if(tail > 0) fTOFtail=tail; else printf("TOF tail should be greater than 0 (nothing done)\n");}; void SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option); @@ -199,6 +205,7 @@ private: TString fMCperiodTPC; //! corresponding MC period to use for the TPC splines TString fMCperiodUser; // MC prodution requested by the user TString fCurrentFile; //! name of currently processed file + Int_t fCurrentAliRootRev; //! Aliroot rev. used to reconstruct the data Int_t fRecoPass; //! reconstruction pass Int_t fRecoPassUser; // reconstruction pass explicitly set by the user Int_t fRun; //! current run number @@ -207,10 +214,11 @@ private: Float_t fResT0C; //! T0C resolution in current run Float_t fResT0AC; //! T0A.and.T0C resolution in current run - TObjArray *fArrPidResponseMaster; //! TPC pid splines - TF1 *fResolutionCorrection; //! TPC resolution correction - AliOADBContainer* fOADBvoltageMaps; //! container with the voltage maps - Bool_t fUseTPCEtaCorrection; // Use TPC eta correction + TObjArray *fArrPidResponseMaster; //! TPC pid splines + TF1 *fResolutionCorrection; //! TPC resolution correction + AliOADBContainer* fOADBvoltageMaps; //! container with the voltage maps + Bool_t fUseTPCEtaCorrection; // Use TPC eta correction + Bool_t fUseTPCMultiplicityCorrection; // Use TPC multiplicity correction AliTRDPIDResponseObject *fTRDPIDResponseObject; //! TRD PID Response Object @@ -305,7 +313,7 @@ private: EDetPidStatus GetPHOSPIDStatus(const AliVTrack *track) const; EDetPidStatus GetEMCALPIDStatus(const AliVTrack *track) const; - ClassDef(AliPIDResponse, 11); //PID response handling + ClassDef(AliPIDResponse, 12); //PID response handling }; #endif diff --git a/STEER/STEERBase/AliTPCPIDResponse.cxx b/STEER/STEERBase/AliTPCPIDResponse.cxx index bbf370f9fa8..367d8cc2e91 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.cxx +++ b/STEER/STEERBase/AliTPCPIDResponse.cxx @@ -71,12 +71,26 @@ AliTPCPIDResponse::AliTPCPIDResponse(): fMagField(0.), fhEtaCorr(0x0), fhEtaSigmaPar1(0x0), - fSigmaPar0(0.0) + fSigmaPar0(0.0), + fCurrentEventMultiplicity(0), + fCorrFuncMultiplicity(0x0), + fCorrFuncMultiplicityTanTheta(0x0), + fCorrFuncSigmaMultiplicity(0x0) { // // The default constructor // for (Int_t i=0; iSetDirectory(0); } + + // Copy multiplicity correction functions + if (that.fCorrFuncMultiplicity) { + fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity)); + } + + if (that.fCorrFuncMultiplicityTanTheta) { + fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta)); + } + + if (that.fCorrFuncSigmaMultiplicity) { + fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity)); + } } //_________________________________________________________________________ @@ -185,6 +225,7 @@ AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that) fBadOROCthreshhold=that.fBadOROCthreshhold; fMaxBadLengthFraction=that.fMaxBadLengthFraction; fMagField=that.fMagField; + fCurrentEventMultiplicity=that.fCurrentEventMultiplicity; for (Int_t i=0; iEval(mom/mass) * chargeFactor; - - if (!correctEta) + + if (!correctEta && !correctMultiplicity) return dEdxSplines; - //TODO Alternatively take current track dEdx - //return dEdxSplines * GetEtaCorrection(track, dEdx); - return dEdxSplines * GetEtaCorrection(track, dEdxSplines); + Double_t corrFactorEta = 1.0; + Double_t corrFactorMultiplicity = 1.0; + + if (correctEta) { + corrFactorEta = GetEtaCorrection(track, dEdxSplines); + //TODO Alternatively take current track dEdx + //corrFactorEta = GetEtaCorrection(track, dEdx); + } + + if (correctMultiplicity) + corrFactorMultiplicity = GetMultiplicityCorrection(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity); + + return dEdxSplines * corrFactorEta * corrFactorMultiplicity; } @@ -359,13 +430,15 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource, - Bool_t correctEta) const + Bool_t correctEta, + Bool_t correctMultiplicity) const { // Calculates the expected PID signal as the function of // the information stored in the track, for the specified particle type // // At the moment, these signals are just the results of calling the - // Bethe-Bloch formula plus, if desired, taking into account the eta dependence. + // Bethe-Bloch formula plus, if desired, taking into account the eta dependence + // and the multiplicity dependence (for PbPb). // This can be improved. By taking into account the number of // assigned clusters and/or the track dip angle, for example. // @@ -389,7 +462,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, } // Charge factor already taken into account inside the following function call - return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta); + return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); } //_________________________________________________________________________ @@ -451,7 +524,8 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, Double_t dEdx, Int_t nPoints, const TSpline3* responseFunction, - Bool_t correctEta) const + Bool_t correctEta, + Bool_t correctMultiplicity) const { // Calculates the expected sigma of the PID signal as the function of // the information stored in the track and the given parameters, @@ -460,20 +534,37 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, if (!responseFunction) return 999; - + + //TODO Check whether it makes sense to set correctMultiplicity to kTRUE while correctEta might be kFALSE // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation if (!fhEtaSigmaPar1 || !correctEta) { if (nPoints != 0) - return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE) * + return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity) * fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints); else - return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE)*fRes0[gainScenario]; + return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity)*fRes0[gainScenario]; } if (nPoints > 0) { + // Use eta correction (+ eta-dependent sigma) Double_t sigmaPar1 = GetSigmaPar1(track, species, dEdx, responseFunction); - return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE)*sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints); + if (correctMultiplicity) { + // In addition, take into account multiplicity dependence of mean and sigma of dEdx + Double_t dEdxExpectedEtaCorrected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); + + // GetMultiplicityCorrection and GetMultiplicitySigmaCorrection both need the eta corrected dEdxExpected + Double_t multiplicityCorrFactor = GetMultiplicityCorrection(track, dEdxExpectedEtaCorrected, fCurrentEventMultiplicity); + Double_t multiplicitySigmaCorrFactor = GetMultiplicitySigmaCorrection(dEdxExpectedEtaCorrected, fCurrentEventMultiplicity); + + // multiplicityCorrFactor to correct dEdxExpected for multiplicity. In addition: Correction factor for sigma + return (dEdxExpectedEtaCorrected * multiplicityCorrFactor) + * (sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints) * multiplicitySigmaCorrFactor); + } + else { + return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE)* + sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints); + } } else { // One should never have/take tracks with 0 dEdx clusters! @@ -486,7 +577,8 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource, - Bool_t correctEta) const + Bool_t correctEta, + Bool_t correctMultiplicity) const { // Calculates the expected sigma of the PID signal as the function of // the information stored in the track, for the specified particle type @@ -501,7 +593,7 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) return 999; //TODO: better handling! - return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta); + return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity); } @@ -509,7 +601,8 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource, - Bool_t correctEta) const + Bool_t correctEta, + Bool_t correctMultiplicity) const { //Calculates the number of sigmas of the PID signal from the expected value //for a given particle species in the presence of multiple gain scenarios @@ -523,8 +616,8 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) return -999; //TODO: Better handling! - Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta); - Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta); + Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); + Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity); // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters if (sigma >= 998) return -999; @@ -537,7 +630,8 @@ Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource, Bool_t correctEta, - Bool_t ratio/*=kFALSE*/) const + Bool_t correctMultiplicity, + Bool_t ratio/*=kFALSE*/)const { //Calculates the number of sigmas of the PID signal from the expected value //for a given particle species in the presence of multiple gain scenarios @@ -551,12 +645,12 @@ Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track, if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) return -9999.; //TODO: Better handling! - const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta); + const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); Double_t delta=-9999.; if (!ratio) delta=dEdx-bethe; else if (bethe>1.e-20) delta=dEdx/bethe; - + return delta; } @@ -677,12 +771,7 @@ Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, Double_t dE if (tpcSignal < 1.) return 1.; - // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()). - // However, this value is not available for AODs and, thus, not for AliVTrack. - // Fortunately, the following formula allows to approximate the local tanTheta with the - // global theta angle -> This is for by far most of the tracks the same, but gives at - // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok. - Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0); + Double_t tanTheta = GetTrackTanTheta(track); Int_t binX = fhEtaCorr->GetXaxis()->FindBin(tanTheta); Int_t binY = fhEtaCorr->GetYaxis()->FindBin(1. / tpcSignal); @@ -718,7 +807,8 @@ Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EPa if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) return 1.; - Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE); + // For the eta correction, do NOT take the multiplicity corrected value of dEdx + Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); //TODO Alternatively take current track dEdx //return GetEtaCorrection(track, dEdx); @@ -756,7 +846,8 @@ Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, Ali Double_t etaCorr = 0.; if (species < AliPID::kUnknown) { - Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE); + // For the eta correction, do NOT take the multiplicity corrected value of dEdx + Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); etaCorr = GetEtaCorrection(track, dEdxSplines); } else { @@ -770,7 +861,6 @@ Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, Ali } - //_________________________________________________________________________ Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, const TSpline3* responseFunction) const { @@ -791,8 +881,9 @@ Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EPartic // sigma parameter! // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!) // of such a particle, which by assumption then has this dEdx value - - Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE); + + // For the eta correction, do NOT take the multiplicity corrected value of dEdx + Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); if (dEdxExpected < 1.) return 999; @@ -824,7 +915,7 @@ Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EPartic Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const { // - // Get eta correction for the given track. + // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track. // if (!fhEtaSigmaPar1) @@ -865,6 +956,7 @@ Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap) return kTRUE; } + //_________________________________________________________________________ Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0) { @@ -893,6 +985,256 @@ Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0 } +//_________________________________________________________________________ +Double_t AliTPCPIDResponse::GetTrackTanTheta(const AliVTrack *track) const +{ + // Extract the tanTheta from the information available in the AliVTrack + + // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()). + // However, this value is not available for AODs and, thus, not for AliVTrack. + // Fortunately, the following formula allows to approximate the local tanTheta with the + // global theta angle -> This is for by far most of the tracks the same, but gives at + // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok. + + + // Alternatively (in average, the effect is found to be negligable!): + // Take local tanTheta from TPC inner wall, if available (currently only for ESDs available) + /*const AliExternalTrackParam* innerParam = track->GetInnerParam(); + if (innerParam) { + return innerParam->GetTgl(); + }*/ + + return TMath::Tan(-track->Theta() + TMath::Pi() / 2.0); +} + + +//_________________________________________________________________________ +Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, const Double_t dEdxExpected, const Int_t multiplicity) const +{ + // Calculate the multiplicity correction factor for this track for the given multiplicity. + // The parameter dEdxExpected should take into account the eta correction already! + + // Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta + // => Use eta corrected value for dEdexExpected. + + if (dEdxExpected <= 0 || multiplicity <= 0) + return 1.0; + + const Double_t dEdxExpectedInv = 1. / dEdxExpected; + Double_t relSlope = fCorrFuncMultiplicity->Eval(dEdxExpectedInv); + + const Double_t tanTheta = GetTrackTanTheta(track); + relSlope += fCorrFuncMultiplicityTanTheta->Eval(tanTheta); + + return (1. + relSlope * multiplicity); +} + + +//_________________________________________________________________________ +Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const +{ + // + // Get multiplicity correction for the given track (w.r.t. the mulitplicity of the current event) + // + + //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse + + // No correction in case of multiplicity <= 0 + if (fCurrentEventMultiplicity <= 0) + return 1.; + + Double_t dEdx = -1; + Int_t nPoints = -1; + ETPCgainScenario gainScenario = kGainScenarioInvalid; + TSpline3* responseFunction = 0x0; + + if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) + return 1.; + + //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid? + // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) + Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); + + return GetMultiplicityCorrection(track, dEdxExpected, fCurrentEventMultiplicity); +} + + +//_________________________________________________________________________ +Double_t AliTPCPIDResponse::GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const +{ + // + // Get multiplicity corrected dEdx for the given track. For the correction, the expected dEdx of + // the specified species will be used. If the species is set to AliPID::kUnknown, the + // dEdx of the track is used instead. + // WARNING: In the latter case, the correction might not be as good as if the + // expected dEdx is used, which is the way the correction factor is designed + // for. + // In any case, one has to decide carefully to which expected signal one wants to + // compare the corrected value - to the corrected or uncorrected. + // Anyhow, a safer way of looking e.g. at the n-sigma is to call + // the corresponding function GetNumberOfSigmas! + // + + //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse + + Double_t dEdx = -1; + Int_t nPoints = -1; + ETPCgainScenario gainScenario = kGainScenarioInvalid; + TSpline3* responseFunction = 0x0; + + // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case + // it is not used anyway, so this causes no trouble. + if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) + return -1.; + + + // No correction in case of multiplicity <= 0 + if (fCurrentEventMultiplicity <= 0) + return dEdx; + + Double_t multiplicityCorr = 0.; + + // TODO Normally, one should use the eta corrected values in BOTH of the following cases. Does it make sense to use the uncorrected dEdx values? + if (species < AliPID::kUnknown) { + // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course). + // However, one needs the eta corrected value! + Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); + multiplicityCorr = GetMultiplicityCorrection(track, dEdxSplines, fCurrentEventMultiplicity); + } + else { + // One needs the eta corrected value to determine the multiplicity correction factor! + Double_t etaCorr = GetEtaCorrection(track, dEdx); + multiplicityCorr = GetMultiplicityCorrection(track, dEdx * etaCorr, fCurrentEventMultiplicity); + } + + if (multiplicityCorr <= 0) + return -1.; + + return dEdx / multiplicityCorr; +} + + +//_________________________________________________________________________ +Double_t AliTPCPIDResponse::GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, + ETPCdEdxSource dedxSource) const +{ + // + // Get multiplicity and eta corrected dEdx for the given track. For the correction, + // the expected dEdx of the specified species will be used. If the species is set + // to AliPID::kUnknown, the dEdx of the track is used instead. + // WARNING: In the latter case, the correction might not be as good as if the + // expected dEdx is used, which is the way the correction factor is designed + // for. + // In any case, one has to decide carefully to which expected signal one wants to + // compare the corrected value - to the corrected or uncorrected. + // Anyhow, a safer way of looking e.g. at the n-sigma is to call + // the corresponding function GetNumberOfSigmas! + // + + Double_t dEdx = -1; + Int_t nPoints = -1; + ETPCgainScenario gainScenario = kGainScenarioInvalid; + TSpline3* responseFunction = 0x0; + + // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case + // it is not used anyway, so this causes no trouble. + if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) + return -1.; + + Double_t multiplicityCorr = 0.; + Double_t etaCorr = 0.; + + if (species < AliPID::kUnknown) { + // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) + Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); + etaCorr = GetEtaCorrection(track, dEdxSplines); + multiplicityCorr = GetMultiplicityCorrection(track, dEdxSplines * etaCorr, fCurrentEventMultiplicity); + } + else { + etaCorr = GetEtaCorrection(track, dEdx); + multiplicityCorr = GetMultiplicityCorrection(track, dEdx * etaCorr, fCurrentEventMultiplicity); + } + + if (multiplicityCorr <= 0 || etaCorr <= 0) + return -1.; + + return dEdx / multiplicityCorr / etaCorr; +} + + +//_________________________________________________________________________ +Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const Double_t dEdxExpected, const Int_t multiplicity) const +{ + // Calculate the multiplicity sigma correction factor for the corresponding expected dEdx and for the given multiplicity. + // The parameter dEdxExpected should take into account the eta correction already! + + // Multiplicity dependence of sigma depends on the real dEdx at zero multiplicity, + // i.e. the eta (only) corrected dEdxExpected value has to be used + // since all maps etc. have been created for ~zero multiplicity + + if (dEdxExpected <= 0 || multiplicity <= 0) + return 1.0; + + Double_t relSigmaSlope = fCorrFuncSigmaMultiplicity->Eval(1. / dEdxExpected); + + return (1. + relSigmaSlope * multiplicity); +} + + +//_________________________________________________________________________ +Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const +{ + // + // Get multiplicity sigma correction for the given track (w.r.t. the mulitplicity of the current event) + // + + //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse + + // No correction in case of multiplicity <= 0 + if (fCurrentEventMultiplicity <= 0) + return 1.; + + Double_t dEdx = -1; + Int_t nPoints = -1; + ETPCgainScenario gainScenario = kGainScenarioInvalid; + TSpline3* responseFunction = 0x0; + + if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) + return 1.; + + //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid? + // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) + Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); + + return GetMultiplicitySigmaCorrection(dEdxExpected, fCurrentEventMultiplicity); +} + + +//_________________________________________________________________________ +void AliTPCPIDResponse::ResetMultiplicityCorrectionFunctions() +{ + // Default values: No correction, i.e. overall correction factor should be one + + // This function should always return 0.0, if multcorr disabled + fCorrFuncMultiplicity->SetParameter(0, 0.); + fCorrFuncMultiplicity->SetParameter(1, 0.); + fCorrFuncMultiplicity->SetParameter(2, 0.); + fCorrFuncMultiplicity->SetParameter(3, 0.); + fCorrFuncMultiplicity->SetParameter(4, 0.); + + // This function should always return 0., if multCorr disabled + fCorrFuncMultiplicityTanTheta->SetParameter(0, 0.); + fCorrFuncMultiplicityTanTheta->SetParameter(1, 0.); + fCorrFuncMultiplicityTanTheta->SetParameter(2, 0.); + + // This function should always return 0.0, if mutlcorr disabled + fCorrFuncSigmaMultiplicity->SetParameter(0, 0.); + fCorrFuncSigmaMultiplicity->SetParameter(1, 0.); + fCorrFuncSigmaMultiplicity->SetParameter(2, 0.); + fCorrFuncSigmaMultiplicity->SetParameter(3, 0.); +} + + //_________________________________________________________________________ Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner, Double_t* trackPositionOuter, diff --git a/STEER/STEERBase/AliTPCPIDResponse.h b/STEER/STEERBase/AliTPCPIDResponse.h index 6dc5a9863b6..9bf2639a6b7 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.h +++ b/STEER/STEERBase/AliTPCPIDResponse.h @@ -20,6 +20,7 @@ #include #include #include +#include #include "AliPID.h" #include "AliVTrack.h" @@ -89,8 +90,10 @@ public: const TH2D* GetEtaCorrMap() const { return fhEtaCorr; }; Bool_t SetEtaCorrMap(TH2D* hMap); - Double_t GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const; + Double_t GetTrackTanTheta(const AliVTrack *track) const; + Double_t GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const; + Double_t GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const; const TH2D* GetSigmaPar1Map() const { return fhEtaSigmaPar1; }; @@ -99,25 +102,56 @@ public: Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const; + + const TF1* GetMultiplicityCorrectionFunction() const { return fCorrFuncMultiplicity; }; + void SetParameterMultiplicityCorrection(Int_t parIndex, Double_t parValue) + { if (fCorrFuncMultiplicity) fCorrFuncMultiplicity->SetParameter(parIndex, parValue); }; + + const TF1* GetMultiplicityCorrectionFunctionTanTheta() const { return fCorrFuncMultiplicityTanTheta; }; + void SetParameterMultiplicityCorrectionTanTheta(Int_t parIndex, Double_t parValue) + { if (fCorrFuncMultiplicityTanTheta) fCorrFuncMultiplicityTanTheta->SetParameter(parIndex, parValue); }; + + const TF1* GetMultiplicitySigmaCorrectionFunction() const { return fCorrFuncSigmaMultiplicity; }; + void SetParameterMultiplicitySigmaCorrection(Int_t parIndex, Double_t parValue) + { if (fCorrFuncSigmaMultiplicity) fCorrFuncSigmaMultiplicity->SetParameter(parIndex, parValue); }; + + void ResetMultiplicityCorrectionFunctions(); + + void SetCurrentEventMultiplicity(Int_t value) { fCurrentEventMultiplicity = value; }; + Int_t GetCurrentEventMultiplicity() const { return fCurrentEventMultiplicity; }; + + Double_t GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const; + + Double_t GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const; + + Double_t GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const; + + Double_t GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, + ETPCdEdxSource dedxSource = kdEdxDefault) const; //NEW void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario ); Double_t GetExpectedSignal( const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault, - Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE + Bool_t correctEta = kFALSE, + Bool_t correctMultiplicity = kFALSE) const; Double_t GetExpectedSigma( const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault, - Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE + Bool_t correctEta = kFALSE, + Bool_t correctMultiplicity = kFALSE) const; Float_t GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault, - Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE + Bool_t correctEta = kFALSE, + Bool_t correctMultiplicity = kFALSE) const; Float_t GetSignalDelta( const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault, - Bool_t correctEta = kFALSE, Bool_t ratio=kFALSE) const; + Bool_t correctEta = kFALSE, + Bool_t correctMultiplicity = kFALSE, + Bool_t ratio = kFALSE) const; void SetResponseFunction(TObject* o, AliPID::EParticleType type, @@ -155,7 +189,8 @@ public: AliPID::EParticleType n=AliPID::kKaon) const { // // Deprecated function (for backward compatibility). Please use - // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource ) + // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource, + // Bool_t correctEta, Bool_t correctMultiplicity) // instead!TODO // @@ -175,7 +210,8 @@ protected: AliPID::EParticleType species, Double_t dEdx, const TSpline3* responseFunction, - Bool_t correctEta) const; + Bool_t correctEta, + Bool_t correctMultiplicity) const; Double_t GetExpectedSigma(const AliVTrack* track, AliPID::EParticleType species, @@ -183,10 +219,15 @@ protected: Double_t dEdx, Int_t nPoints, const TSpline3* responseFunction, - Bool_t correctEta) const; + Bool_t correctEta, + Bool_t correctMultiplicity) const; Double_t GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const; + Double_t GetMultiplicityCorrection(const AliVTrack *track, const Double_t dEdxExpected, const Int_t multiplicity) const; + + Double_t GetMultiplicitySigmaCorrection(const Double_t dEdxExpected, const Int_t multiplicity) const; + Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, const TSpline3* responseFunction) const; @@ -221,8 +262,13 @@ private: TH2D* fhEtaSigmaPar1; //! Map for parameter 1 of the dEdx sigma parametrisation Double_t fSigmaPar0; // Parameter 0 of the dEdx sigma parametrisation + + Int_t fCurrentEventMultiplicity; // Multiplicity of the current event + TF1* fCorrFuncMultiplicity; //! Function to correct for the multiplicity dependence of the TPC dEdx + TF1* fCorrFuncMultiplicityTanTheta; //! Function to correct the additional tanTheta dependence of the multiplicity dependence of the TPC dEdx + TF1* fCorrFuncSigmaMultiplicity; //! Function to correct for the multiplicity dependence of the TPC dEdx resolution - ClassDef(AliTPCPIDResponse,5) // TPC PID class + ClassDef(AliTPCPIDResponse,6) // TPC PID class }; #endif -- 2.43.0