From a7f878167a31214a50a16558258c39aa496f712f Mon Sep 17 00:00:00 2001 From: miranov25 Date: Sun, 14 Feb 2021 17:57:16 +0100 Subject: [PATCH 01/15] ATO-525,ATO-496 - addig qptTgl correction to slove dEdx splitting at pt<0.5 GeV First version --- STEER/STEERBase/AliTPCPIDResponse.cxx | 92 +++++++++++++++++++++++---- STEER/STEERBase/AliTPCPIDResponse.h | 25 +++++--- 2 files changed, 94 insertions(+), 23 deletions(-) diff --git a/STEER/STEERBase/AliTPCPIDResponse.cxx b/STEER/STEERBase/AliTPCPIDResponse.cxx index 876af92a0ef..9cd278ba418 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.cxx +++ b/STEER/STEERBase/AliTPCPIDResponse.cxx @@ -78,6 +78,7 @@ AliTPCPIDResponse::AliTPCPIDResponse(): fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios), fOADBContainer(0x0), fPileupCorrection(0x0), + fQPtTglCorrection(0x0), fVoltageMap(72), fLowGainIROCthreshold(-40), fBadIROCthreshhold(-70), @@ -105,6 +106,7 @@ AliTPCPIDResponse::AliTPCPIDResponse(): fOROClongWeight(1.), fPileupCorrectionStrategy(kPileupCorrectionInExpectedSignal), fPileupCorrectionRequested(kFALSE), + fQPtTglCorrectionRequested(kFALSE), fSigmaParametrization(0x0), fMultiplicityNormalization(1), fRecoPassNameUsed(), @@ -243,6 +245,7 @@ AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that): fOROClongWeight(that.fOROClongWeight), fPileupCorrectionStrategy(that.fPileupCorrectionStrategy), fPileupCorrectionRequested(that.fPileupCorrectionRequested), + fQPtTglCorrectionRequested(that.fQPtTglCorrectionRequested), fSigmaParametrization(that.fSigmaParametrization), fMultiplicityNormalization(that.fMultiplicityNormalization), fRecoPassNameUsed(that.fRecoPassNameUsed), @@ -367,7 +370,7 @@ AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that) fOROClongWeight = that.fOROClongWeight; fPileupCorrectionStrategy = that.fPileupCorrectionStrategy; fPileupCorrectionRequested = that.fPileupCorrectionRequested; - fPileupCorrectionRequested = that.fPileupCorrectionRequested; + fQPtTglCorrectionRequested = that.fQPtTglCorrectionRequested; fSigmaParametrization = that.fSigmaParametrization; fMultiplicityNormalization = that.fMultiplicityNormalization; fRecoPassNameUsed = that.fRecoPassNameUsed; @@ -493,7 +496,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, const TSpline3* responseFunction, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection) const + Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const { // Calculates the expected PID signal as the function of // the information stored in the track and the given parameters, @@ -551,7 +554,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, ETPCdEdxSource dedxSource, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection) const + Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const { // Calculates the expected PID signal as the function of // the information stored in the track, for the specified particle type @@ -582,7 +585,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, correctMultiplicity, usePileupCorrection); + return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); } //_________________________________________________________________________ @@ -646,7 +649,7 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, const TSpline3* responseFunction, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection) const + Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const { // Calculates the expected sigma of the PID signal as the function of // the information stored in the track and the given parameters, @@ -671,10 +674,10 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, // 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, correctMultiplicity, usePileupCorrection) * + return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection, useQPtTglCorrection) * fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints); else - return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection)*fRes0[gainScenario]; + return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection, useQPtTglCorrection)*fRes0[gainScenario]; } if (nPoints > 0) { @@ -711,7 +714,7 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, ETPCdEdxSource dedxSource, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection) const + Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const { // Calculates the expected sigma of the PID signal as the function of // the information stored in the track, for the specified particle type @@ -726,7 +729,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, correctMultiplicity, usePileupCorrection); + return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); } @@ -736,7 +739,7 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, ETPCdEdxSource dedxSource, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection) const + Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) 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 @@ -750,8 +753,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, correctMultiplicity, usePileupCorrection); - Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection); + Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); + Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters if (sigma >= 998) return -999; @@ -765,7 +768,7 @@ Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track, ETPCdEdxSource dedxSource, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection /*= kFALSE*/, + Bool_t usePileupCorrection /*= kFALSE*/, Bool_t useQPtTglCorrection, Bool_t ratio/*=kFALSE*/)const { //Calculates the number of sigmas of the PID signal from the expected value @@ -780,7 +783,7 @@ 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, correctMultiplicity, usePileupCorrection); + const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); Double_t delta=-9999.; if (!ratio) delta=dEdx-bethe; @@ -1837,6 +1840,28 @@ Double_t AliTPCPIDResponse::GetPileupCorrectionValue(const AliVTrack* track) con return corrPileup; } + +//_____________________________________________________________________________ +Double_t AliTPCPIDResponse::GetQPtTglCorrectionValue(const AliVTrack* track) const +{ + // + // The correction is an additive value. The corrected dEdx is + // dEdx - corrQPtTgl + // + if (!fQPtTglCorrection) { + return 0.; + } + + //const Double_t trackTgl = TMath::Abs(TMath::SinH(track->GetTPCTgl())); + const Double_t trackTgl = track->GetTgl(); + Double_t corrVals[2] = {track->GetInnerParam()->GetSigned1Pt(), trackTgl}; + const Double_t corrQPtTgl = (1+fQPtTglCorrection->Eval(corrVals)) * 50; + + return corrQPtTgl; +} + + + //______________________________________________________________________________ AliNDLocalRegression* AliTPCPIDResponse::GetPileupCorrectionFromFile(const TString fileName) { @@ -1874,6 +1899,45 @@ AliNDLocalRegression* AliTPCPIDResponse::GetPileupCorrectionFromFile(const TStri return ndLocal; } + +//______________________________________________________________________________ +AliNDLocalRegression* AliTPCPIDResponse::GetQPtTglCorrectionFromFile(const TString fileName) +{ + if (fileName.Contains("alien://") && !gGrid) { + TGrid::Connect("alien"); + } + + TFile* f = TFile::Open(fileName); + if (!f || !f->IsOpen() || f->IsZombie()) { + printf("AliTPCPIDResponse::GetQPtTglCorrectionFromFile: Could not open QPtTgl correction file: %s", fileName.Data()); + delete f; + return 0x0; + } + + // Assume that the only object in the file is the QPtTgl correction + const Int_t numberOfKeys = f->GetListOfKeys()->GetEntries(); + if ( numberOfKeys != 1) { + printf("AliTPCPIDResponse::GetQPtTglCorrectionFromFile: Number of objects in the file %d is not 1", numberOfKeys); + delete f; + return 0x0; + } + + const char* objectName = f->GetListOfKeys()->At(0)->GetName(); + TObject* o = f->Get(objectName); + AliNDLocalRegression* ndLocal = dynamic_cast(o); + if (!ndLocal) { + printf("AliTPCPIDResponse::GetQPtTglCorrectionFromFile: Object in file %s with name %s is not of type AliNDLocalRegression but %s", fileName.Data(), objectName, o->IsA()->GetName()); + delete f; + return 0x0; + } + + ndLocal->SetName("QPtTglCorrection"); + + delete f; + return ndLocal; +} + + //_____________________________________________________________________________ Bool_t AliTPCPIDResponse::InitFromOADB(const Int_t run, const Int_t pass, TString passName, const char* oadbFile/*="$ALICE_PHYSICS/OADB/COMMON/PID/data/TPCPIDResponseOADB.root"*/, diff --git a/STEER/STEERBase/AliTPCPIDResponse.h b/STEER/STEERBase/AliTPCPIDResponse.h index 4459b181d8f..6b49cae2d26 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.h +++ b/STEER/STEERBase/AliTPCPIDResponse.h @@ -173,19 +173,19 @@ class AliTPCPIDResponse: public TNamed { ETPCdEdxSource dedxSource = kdEdxDefault, Bool_t correctEta = kFALSE, Bool_t correctMultiplicity = kFALSE, - Bool_t usePileupCorrection = kFALSE) const; + Bool_t usePileupCorrection = kFALSE, Bool_t useQPtTglCorrection = kFALSE) const; Double_t GetExpectedSigma( const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault, Bool_t correctEta = kFALSE, Bool_t correctMultiplicity = kFALSE, - Bool_t usePileupCorrection = kFALSE) const; + Bool_t usePileupCorrection = kFALSE, Bool_t useQPtTglCorrection = kFALSE) const; Float_t GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault, Bool_t correctEta = kFALSE, Bool_t correctMultiplicity = kFALSE, - Bool_t usePileupCorrection = kFALSE) const; + Bool_t usePileupCorrection = kFALSE, Bool_t useQPtTglCorrection = kFALSE) const; Float_t GetSignalDelta( const AliVTrack* track, AliPID::EParticleType species, @@ -193,7 +193,7 @@ class AliTPCPIDResponse: public TNamed { Bool_t correctEta = kFALSE, Bool_t correctMultiplicity = kFALSE, Bool_t usePileupCorrection = kFALSE, - Bool_t ratio = kFALSE) const; + Bool_t ratio = kFALSE, Bool_t useQPtTglCorrection = kFALSE) const; void SetResponseFunction(TObject* o, AliPID::EParticleType type, @@ -258,12 +258,17 @@ class AliTPCPIDResponse: public TNamed { void SetPileupCorrectionObject(AliNDLocalRegression* correction) { fPileupCorrection = correction; } const AliNDLocalRegression* GetPileupCorrectionObject() const { return fPileupCorrection; } - Bool_t IsPileupCorrectionRequested() const { return fPileupCorrectionRequested; } - Double_t GetPileupCorrectionValue(const AliVTrack* track) const; - static AliNDLocalRegression* GetPileupCorrectionFromFile(const TString fileName); + // + void SetQPtTglCorrectionObject(AliNDLocalRegression* correction) { fQPtTglCorrection = correction; } + const AliNDLocalRegression* GetQPtTglCorrectionObject() const { return fQPtTglCorrection; } + Bool_t IsQPtTglCorrectionRequested() const { return fQPtTglCorrectionRequested; } + Double_t GetQPtTglCorrectionValue(const AliVTrack* track) const; + static AliNDLocalRegression* GetQPtTglCorrectionFromFile(const TString fileName); + + // //===| dEdx type functions |================================================== void SetdEdxType(ETPCdEdxType dEdxType, Int_t dEdxChargeType=0, Int_t dEdxWeightType=0, Double_t dEdxIROCweight=1., Double_t dEdxOROCmedWeight=1., Double_t dEdxOROClongWeight=1.) { fdEdxType=dEdxType; fdEdxChargeType=dEdxChargeType; fdEdxWeightType=dEdxWeightType; fIROCweight=dEdxIROCweight; fOROCmedWeight=dEdxOROCmedWeight; fOROClongWeight=dEdxOROClongWeight; } @@ -318,7 +323,7 @@ class AliTPCPIDResponse: public TNamed { const TSpline3* responseFunction, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection) const; + Bool_t usePileupCorrection, Bool_t useQPtTglCorrection = kFALSE) const; Double_t GetExpectedSigma(const AliVTrack* track, AliPID::EParticleType species, @@ -328,7 +333,7 @@ class AliTPCPIDResponse: public TNamed { const TSpline3* responseFunction, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection) const; + Bool_t usePileupCorrection, Bool_t useQPtTglCorrection = kFALSE) const; // // function for numberical debugging 0 registed splines can be used in the TFormula and tree visualizations // @@ -348,6 +353,7 @@ class AliTPCPIDResponse: public TNamed { TObjArray fResponseFunctions; //! ObjArray of response functions individually for each particle AliOADBContainer* fOADBContainer; //! OADB container with response functions AliNDLocalRegression* fPileupCorrection; // pileup correction object + AliNDLocalRegression* fQPtTglCorrection; // pileup correction object TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers Float_t fLowGainIROCthreshold; //voltage threshold below which the IROC is considered low gain Float_t fBadIROCthreshhold; //voltage threshold for bad IROCS @@ -387,6 +393,7 @@ class AliTPCPIDResponse: public TNamed { ETPCPileupCorrectionStrategy fPileupCorrectionStrategy; // Pileup correction strategy Bool_t fPileupCorrectionRequested; // If pileup correction was configured in the OADB object + Bool_t fQPtTglCorrectionRequested; // If pileup correction was configured in the OADB object //===| New resolution parametrization |======================================= TF1* fSigmaParametrization; // Resolution parametrization From 86bd8aa359301ea2bf9d707012373f6e84661419 Mon Sep 17 00:00:00 2001 From: miranov25 Date: Wed, 17 Feb 2021 17:01:36 +0100 Subject: [PATCH 02/15] ATO-525 - bug fix for ROOT6 - makeCacheTree - have to check upper bound in case above upper bound -statistic get lost --- STAT/AliTreePlayer.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/STAT/AliTreePlayer.cxx b/STAT/AliTreePlayer.cxx index 26b4d8a41bd..2bcd536f730 100644 --- a/STAT/AliTreePlayer.cxx +++ b/STAT/AliTreePlayer.cxx @@ -1647,6 +1647,7 @@ void AliTreePlayer::MakeCacheTreeChunk(TTree * tree, TString varList, TString ou for (Int_t iEntry=firstEntry; iEntrySetEstimate(chunkSize); + if (chunkSize+iEntry>entriesAll) chunkSize=entriesAll-iEntry; // ROOT6 does not handle properly query above limt Int_t entries = tree->Draw(varList.Data(), selection, "goffpara", chunkSize, iEntry); if (entries<=0) break; if (entries > estimate) { From 9b6cb13ea86e4e1ac29168ad0196886317e8ebd7 Mon Sep 17 00:00:00 2001 From: miranov25 Date: Fri, 19 Feb 2021 08:58:57 +0100 Subject: [PATCH 03/15] ATO-525 - bug fix for ROOT6 - makeCacheTree - have to check upper bound in case above upper bound -statistic get lost --- STAT/AliTreePlayer.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/STAT/AliTreePlayer.cxx b/STAT/AliTreePlayer.cxx index 2bcd536f730..c9da85fc4c0 100644 --- a/STAT/AliTreePlayer.cxx +++ b/STAT/AliTreePlayer.cxx @@ -1647,7 +1647,7 @@ void AliTreePlayer::MakeCacheTreeChunk(TTree * tree, TString varList, TString ou for (Int_t iEntry=firstEntry; iEntrySetEstimate(chunkSize); - if (chunkSize+iEntry>entriesAll) chunkSize=entriesAll-iEntry; // ROOT6 does not handle properly query above limt + if (chunkSize+iEntry>=entriesAll) chunkSize=entriesAll-iEntry; // ROOT6 does not handle properly query above limit Int_t entries = tree->Draw(varList.Data(), selection, "goffpara", chunkSize, iEntry); if (entries<=0) break; if (entries > estimate) { From 0ed32806723a56db67e18eeaf71658e3486b55c2 Mon Sep 17 00:00:00 2001 From: miranov25 Date: Mon, 22 Mar 2021 19:42:22 +0100 Subject: [PATCH 04/15] ATO-525 - make histograms - using long for the size print --- STAT/AliTreePlayer.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/STAT/AliTreePlayer.cxx b/STAT/AliTreePlayer.cxx index c9da85fc4c0..1faa56902b3 100644 --- a/STAT/AliTreePlayer.cxx +++ b/STAT/AliTreePlayer.cxx @@ -1172,7 +1172,7 @@ TObjArray * AliTreePlayer::MakeHistograms(TTree * tree, TString hisString, TStr Int_t nExpressions=hisString.CountChar(':')+hisString.CountChar(';')+1; TObjArray * formulaArray = new TObjArray(nExpressions); // array of all expressions - OWNER TString queryString = ""; - Int_t hisSizeFull=0; + Long_t hisSizeFull=0; // // 1.) Analyze formula, book list of TObjString // @@ -1307,7 +1307,7 @@ TObjArray * AliTreePlayer::MakeHistograms(TTree * tree, TString hisString, TStr } } if (verbose&0x1) { - ::Info("AliTreePlayer::MakeHistograms","Total size=%d",hisSizeFull); + ::Info("AliTreePlayer::MakeHistograms","Total size=%lld",hisSizeFull); } } // 2.3 fill histograms From 50701f196e9d5121b3af22a69c7292383ece895a Mon Sep 17 00:00:00 2001 From: miranov25 Date: Mon, 12 Apr 2021 16:44:01 +0200 Subject: [PATCH 05/15] PWGPP-620 - adding sqrt(1+tgl**2) to the energu loss numerical derivative formula --- STEER/ESD/AliESDv0.cxx | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/STEER/ESD/AliESDv0.cxx b/STEER/ESD/AliESDv0.cxx index 295f715e433..c2a9f554c90 100644 --- a/STEER/ESD/AliESDv0.cxx +++ b/STEER/ESD/AliESDv0.cxx @@ -986,15 +986,23 @@ Double_t AliESDv0::GetKFInfoScale(UInt_t p1, UInt_t p2, Int_t type, Double_t d1p // calculate delta of energy loss Double_t bg1=paramP.P() /TDatabasePDG::Instance()->GetParticle(spdg[p1])->Mass(); Double_t bg2=paramN.P() /TDatabasePDG::Instance()->GetParticle(spdg[p2])->Mass(); - Double_t dP1=AliExternalTrackParam::BetheBlochGeant(bg1)/AliExternalTrackParam::BetheBlochGeant(3.); ///relative energu loss - Double_t dP2=AliExternalTrackParam::BetheBlochGeant(bg2)/AliExternalTrackParam::BetheBlochGeant(3.); + Double_t dE1=AliExternalTrackParam::BetheBlochGeant(bg1)/AliExternalTrackParam::BetheBlochGeant(3.); ///relative energu loss + Double_t dE2=AliExternalTrackParam::BetheBlochGeant(bg2)/AliExternalTrackParam::BetheBlochGeant(3.); + dE1*=TMath::Sqrt(1+paramP.GetTgl()*paramP.GetTgl()); + dE2*=TMath::Sqrt(1+paramN.GetTgl()*paramN.GetTgl()); + Double_t mass1=AliPID::ParticleMass(p1); + Double_t mass2=AliPID::ParticleMass(p2); + Double_t E1 = paramP.P()*paramP.P()+mass1*mass1+eLoss*dE1; + Double_t E2 = paramN.P()*paramN.P()+mass2*mass2+eLoss*dE2; + Double_t dP1=TMath::Sqrt(E1*E1-mass1*mass1)-paramP.P(); + Double_t dP2=TMath::Sqrt(E2*E2-mass2*mass2)-paramN.P(); Double_t *pparam1 = (Double_t*)paramP.GetParameter(); Double_t *pparam2 = (Double_t*)paramN.GetParameter(); if (flag&0x1) pparam1[4]+=d1pt; if (flag&0x2) pparam2[4]+=d1pt; - if (flag&0x1) pparam1[4]*=(1+s1pt)*(1.+eLoss*dP1/paramP.P()); /// TODO - use relative energy loss - if (flag&0x2) pparam2[4]*=(1+s1pt)*(1.+eLoss*dP2/paramN.P()); /// + if (flag&0x1) pparam1[4]*=(1+s1pt)*(1.+dP1/paramP.P()); /// TODO - use relative energy loss + if (flag&0x2) pparam2[4]*=(1+s1pt)*(1.+dP2/paramN.P()); /// // AliKFParticle kfp1( paramP, spdg[p1] *TMath::Sign(1,p1) ); From 2640966503c4f5720eaf125ae6aaf2140127253c Mon Sep 17 00:00:00 2001 From: miranov25 Date: Tue, 13 Apr 2021 15:30:09 +0200 Subject: [PATCH 06/15] PWGPP-620 - floating error protection -use 0 as minimal momenta --- STEER/ESD/AliESDv0.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/STEER/ESD/AliESDv0.cxx b/STEER/ESD/AliESDv0.cxx index c2a9f554c90..a66a4b5dd2d 100644 --- a/STEER/ESD/AliESDv0.cxx +++ b/STEER/ESD/AliESDv0.cxx @@ -994,8 +994,8 @@ Double_t AliESDv0::GetKFInfoScale(UInt_t p1, UInt_t p2, Int_t type, Double_t d1p Double_t mass2=AliPID::ParticleMass(p2); Double_t E1 = paramP.P()*paramP.P()+mass1*mass1+eLoss*dE1; Double_t E2 = paramN.P()*paramN.P()+mass2*mass2+eLoss*dE2; - Double_t dP1=TMath::Sqrt(E1*E1-mass1*mass1)-paramP.P(); - Double_t dP2=TMath::Sqrt(E2*E2-mass2*mass2)-paramN.P(); + Double_t dP1=TMath::Sqrt(TMath::Max(E1*E1-mass1*mass1,0.0))-paramP.P(); + Double_t dP2=TMath::Sqrt(TMath::Max(E2*E2-mass2*mass2,0.0))-paramN.P(); Double_t *pparam1 = (Double_t*)paramP.GetParameter(); Double_t *pparam2 = (Double_t*)paramN.GetParameter(); From 233e21943f1d6d598d7ab0348b5e9c4531001464 Mon Sep 17 00:00:00 2001 From: miranov25 Date: Wed, 14 Apr 2021 15:46:22 +0200 Subject: [PATCH 07/15] PWGPP-620 - bug fix - use sqrt --- STEER/ESD/AliESDv0.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/STEER/ESD/AliESDv0.cxx b/STEER/ESD/AliESDv0.cxx index a66a4b5dd2d..758a709e304 100644 --- a/STEER/ESD/AliESDv0.cxx +++ b/STEER/ESD/AliESDv0.cxx @@ -988,15 +988,15 @@ Double_t AliESDv0::GetKFInfoScale(UInt_t p1, UInt_t p2, Int_t type, Double_t d1p Double_t bg2=paramN.P() /TDatabasePDG::Instance()->GetParticle(spdg[p2])->Mass(); Double_t dE1=AliExternalTrackParam::BetheBlochGeant(bg1)/AliExternalTrackParam::BetheBlochGeant(3.); ///relative energu loss Double_t dE2=AliExternalTrackParam::BetheBlochGeant(bg2)/AliExternalTrackParam::BetheBlochGeant(3.); - dE1*=TMath::Sqrt(1+paramP.GetTgl()*paramP.GetTgl()); - dE2*=TMath::Sqrt(1+paramN.GetTgl()*paramN.GetTgl()); + dE1*=TMath::Sqrt(1+(paramP.GetTgl()*paramP.GetTgl())); + dE2*=TMath::Sqrt(1+(paramN.GetTgl()*paramN.GetTgl())); Double_t mass1=AliPID::ParticleMass(p1); Double_t mass2=AliPID::ParticleMass(p2); - Double_t E1 = paramP.P()*paramP.P()+mass1*mass1+eLoss*dE1; - Double_t E2 = paramN.P()*paramN.P()+mass2*mass2+eLoss*dE2; + Double_t E1 = TMath::Sqrt(paramP.P()*paramP.P()+mass1*mass1)+eLoss*dE1; + Double_t E2 = TMath::Sqrt(paramN.P()*paramN.P()+mass2*mass2)+eLoss*dE2; Double_t dP1=TMath::Sqrt(TMath::Max(E1*E1-mass1*mass1,0.0))-paramP.P(); Double_t dP2=TMath::Sqrt(TMath::Max(E2*E2-mass2*mass2,0.0))-paramN.P(); - + //printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",dE1,dE2, E1,E2, dP1,dP2,eLoss); Double_t *pparam1 = (Double_t*)paramP.GetParameter(); Double_t *pparam2 = (Double_t*)paramN.GetParameter(); if (flag&0x1) pparam1[4]+=d1pt; From 8858b3eff22930d4409224699c15ed3e00b98fb6 Mon Sep 17 00:00:00 2001 From: miranov25 Date: Mon, 10 May 2021 15:28:07 +0200 Subject: [PATCH 08/15] PWGPP-620 - bug fix - use proper axis i case of the error query --- STAT/TStatToolkit.cxx | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/STAT/TStatToolkit.cxx b/STAT/TStatToolkit.cxx index 16527eb63ec..11a6684c4f3 100644 --- a/STAT/TStatToolkit.cxx +++ b/STAT/TStatToolkit.cxx @@ -1221,17 +1221,31 @@ TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const ch graphNew->GetXaxis()->SetBinLabel(i+1,xName); graphNew->GetX()[i]+=offset; } - if (tree->GetVar(1)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1))>0){ - for(Int_t i=0;iGetXaxis()->SetBinLabel(i+1,tree->GetHistogram()->GetXaxis()->GetBinLabel(i+1)); + if (tree->GetV3()==NULL) { + if (tree->GetVar(1)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1)) > 0) { + for (Int_t i = 0; i < count; i++) { + graphNew->GetXaxis()->SetBinLabel(i + 1, tree->GetHistogram()->GetXaxis()->GetBinLabel(i + 1)); + } } - } - if (tree->GetVar(0)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1))>0 ){ - for(Int_t i=0;iGetYaxis()->SetBinLabel(i+1,tree->GetHistogram()->GetYaxis()->GetBinLabel(i+1)); + if (tree->GetVar(0)->IsInteger() && strlen(tree->GetHistogram()->GetYaxis()->GetBinLabel(1)) > 0) { + for (Int_t i = 0; i < count; i++) { + graphNew->GetYaxis()->SetBinLabel(i + 1, tree->GetHistogram()->GetYaxis()->GetBinLabel(i + 1)); + } + } + }else{ + if (tree->GetVar(1)->IsInteger() && strlen(tree->GetHistogram()->GetYaxis()->GetBinLabel(1)) > 0) { + for (Int_t i = 0; i < count; i++) { + graphNew->GetXaxis()->SetBinLabel(i + 1, tree->GetHistogram()->GetYaxis()->GetBinLabel(i + 1)); + } + } + if (tree->GetVar(0)->IsInteger() && strlen(tree->GetHistogram()->GetZaxis()->GetBinLabel(1)) > 0) { + for (Int_t i = 0; i < count; i++) { + graphNew->GetYaxis()->SetBinLabel(i + 1, tree->GetHistogram()->GetZaxis()->GetBinLabel(i + 1)); + } } } + graphNew->GetHistogram()->SetTitle(""); graphNew->SetMarkerStyle(mstyle); graphNew->SetMarkerColor(mcolor); graphNew->SetLineColor(mcolor); From 5e263cda9923f250dbf00711aba010a5e132a68c Mon Sep 17 00:00:00 2001 From: miranov25 Date: Mon, 12 Jul 2021 06:48:00 +0200 Subject: [PATCH 09/15] ATO-525 - changing order of arguments --- STEER/STEERBase/AliTPCPIDResponse.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/STEER/STEERBase/AliTPCPIDResponse.h b/STEER/STEERBase/AliTPCPIDResponse.h index 6b49cae2d26..75da02288bf 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.h +++ b/STEER/STEERBase/AliTPCPIDResponse.h @@ -193,7 +193,7 @@ class AliTPCPIDResponse: public TNamed { Bool_t correctEta = kFALSE, Bool_t correctMultiplicity = kFALSE, Bool_t usePileupCorrection = kFALSE, - Bool_t ratio = kFALSE, Bool_t useQPtTglCorrection = kFALSE) const; + Bool_t useQPtTglCorrection = kFALSE, Bool_t ratio = kFALSE) const; void SetResponseFunction(TObject* o, AliPID::EParticleType type, From 49e34a1d0ce96c1fa53f6b7751650b8b154c5e3b Mon Sep 17 00:00:00 2001 From: miranov25 Date: Mon, 12 Jul 2021 09:35:55 +0200 Subject: [PATCH 10/15] ATO-525 - using qPtTglCorrection if requested --- STEER/STEERBase/AliTPCPIDResponse.cxx | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/STEER/STEERBase/AliTPCPIDResponse.cxx b/STEER/STEERBase/AliTPCPIDResponse.cxx index 9cd278ba418..99037c8bd05 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.cxx +++ b/STEER/STEERBase/AliTPCPIDResponse.cxx @@ -526,6 +526,11 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, if (usePileupCorrection && fPileupCorrectionStrategy == kPileupCorrectionInExpectedSignal) { corrPileup = GetPileupCorrectionValue(track); } + //qPtTglCorrection + Double_t corrqPtTgl = 1.; + if (useQPtTglCorrection && fQPtTglCorrectionRequested) { + corrqPtTgl = GetQPtTglCorrectionValue(track); + } if (!correctEta && !correctMultiplicity) { return dEdxSplines + corrPileup; @@ -544,7 +549,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity); } - return dEdxSplines * corrFactorEta * corrFactorMultiplicity + corrPileup; + return dEdxSplines * corrFactorEta * corrFactorMultiplicity * corrqPtTgl + corrPileup; } @@ -1855,7 +1860,7 @@ Double_t AliTPCPIDResponse::GetQPtTglCorrectionValue(const AliVTrack* track) con //const Double_t trackTgl = TMath::Abs(TMath::SinH(track->GetTPCTgl())); const Double_t trackTgl = track->GetTgl(); Double_t corrVals[2] = {track->GetInnerParam()->GetSigned1Pt(), trackTgl}; - const Double_t corrQPtTgl = (1+fQPtTglCorrection->Eval(corrVals)) * 50; + const Double_t corrQPtTgl = (1+fQPtTglCorrection->Eval(corrVals)) ; return corrQPtTgl; } From 033a34f0e10c473f06bb255ad9b1a071b6460c41 Mon Sep 17 00:00:00 2001 From: miranov25 Date: Mon, 12 Jul 2021 13:57:38 +0200 Subject: [PATCH 11/15] ATO-525 - adding static function for the energy loss correction in TPC regions --- STEER/STEERBase/AliTPCPIDResponse.cxx | 106 ++++++++++++++++++++++++++ STEER/STEERBase/AliTPCPIDResponse.h | 16 +++- 2 files changed, 118 insertions(+), 4 deletions(-) diff --git a/STEER/STEERBase/AliTPCPIDResponse.cxx b/STEER/STEERBase/AliTPCPIDResponse.cxx index 99037c8bd05..7c080d1dd60 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.cxx +++ b/STEER/STEERBase/AliTPCPIDResponse.cxx @@ -2447,3 +2447,109 @@ double AliTPCPIDResponse::sigmadEdxPt(double p, double PID, double resol ){ double deltaRel = TMath::Abs(dEdx2-dEdx)/dEdx; return deltaRel; } + + +/// return pt after traversing length in TPC volume assummig AliExternalTrackParam::BetheBlochAleph standard parameterization +/// formula approximation for primary and new to primary tracks +/// \param ptIn - input Pt +/// \param tgl - pz/pt +/// \param mass - particle mass +/// \param type - 0-Argon, 1-Neon +/// \param bz - magnetic field in kG +/// \param fraction - scaling factor for energy loss +/// \return +double AliTPCPIDResponse::GetPtOutHelix(Float_t ptIn, Double_t tgl, Float_t mass, Double_t rIn, Double_t rOut, Int_t type, Float_t bz,Float_t fraction){ + const Float_t kB2C=-0.299792458e-3; + const Float_t kMaxSnp=0.95; + Float_t pin=ptIn*TMath::Sqrt(1+tgl*tgl); + Float_t qPt=1/ptIn; + Float_t bg = pin/mass; + Float_t elossNe = 0.001*0.00156; // page 3 https://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf + Float_t elossAr = 0.001*0.00244; // page 3 https://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf + Float_t dEdx= AliExternalTrackParam::BetheBlochAleph(bg)*fraction; + Double_t e0=TMath::Sqrt(pin*pin+mass*mass); + Float_t eloss=(type==0) ? elossAr:elossNe; + Float_t lengthNorm=TMath::Abs(rOut-rIn)*TMath::Sqrt(1+tgl*tgl); + Float_t C = kB2C*bz*qPt; + if (TMath::Abs(0.5*rIn*C)>1) rIn=2.*kMaxSnp/C; + if (TMath::Abs(0.5*rOut*C)>1) rOut=2.*kMaxSnp/C; + Float_t lIn = 2.*TMath::ASin(0.5*rIn*C)/C; + Float_t lOut = 2.*TMath::ASin((0.5*rOut*C))/C; + lengthNorm=(TMath::Abs(lOut)-TMath::Abs(lIn))*TMath::Sqrt(1+tgl*tgl); + // + Double_t e1=e0-dEdx*eloss*lengthNorm; + if (e1rIn0 && radiusMrIn1 && radiusM Date: Thu, 15 Jul 2021 09:10:23 +0200 Subject: [PATCH 12/15] ATO-525 - adding energy loss correction --- STEER/STEERBase/AliTPCPIDResponse.cxx | 63 +++++++++++++++++++++------ STEER/STEERBase/AliTPCPIDResponse.h | 30 +++++++++---- 2 files changed, 72 insertions(+), 21 deletions(-) diff --git a/STEER/STEERBase/AliTPCPIDResponse.cxx b/STEER/STEERBase/AliTPCPIDResponse.cxx index 7c080d1dd60..ce8a535983f 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.cxx +++ b/STEER/STEERBase/AliTPCPIDResponse.cxx @@ -79,6 +79,9 @@ AliTPCPIDResponse::AliTPCPIDResponse(): fOADBContainer(0x0), fPileupCorrection(0x0), fQPtTglCorrection(0x0), + fEnergyLossCorrectionRequested(kFALSE), // flag fEnergyLossCorrectionRequested + fGasType(0), // type of the gas (0-Argon, 1-Neon) + fEnergyLossFactor(1.4), // energy loss facto - fit parameters in respect to tabulated value fVoltageMap(72), fLowGainIROCthreshold(-40), fBadIROCthreshhold(-70), @@ -496,7 +499,9 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, const TSpline3* responseFunction, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const + Bool_t usePileupCorrection, + Bool_t useQPtTglCorrection, + Bool_t useEnergyLossCorrection) const { // Calculates the expected PID signal as the function of // the information stored in the track and the given parameters, @@ -531,6 +536,12 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, if (useQPtTglCorrection && fQPtTglCorrectionRequested) { corrqPtTgl = GetQPtTglCorrectionValue(track); } + // energy loss correction + Double_t corrELoss = 1.; + if (useEnergyLossCorrection && fEnergyLossCorrectionRequested) { + corrELoss = GetEnergyLossCorrectionValue(track,species); + } + if (!correctEta && !correctMultiplicity) { return dEdxSplines + corrPileup; @@ -549,7 +560,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity); } - return dEdxSplines * corrFactorEta * corrFactorMultiplicity * corrqPtTgl + corrPileup; + return dEdxSplines * corrFactorEta * corrFactorMultiplicity * corrqPtTgl *corrELoss+ corrPileup; } @@ -559,7 +570,9 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, ETPCdEdxSource dedxSource, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const + Bool_t usePileupCorrection, + Bool_t useQPtTglCorrection, + Bool_t useEnergyLossCorrection) const { // Calculates the expected PID signal as the function of // the information stored in the track, for the specified particle type @@ -590,7 +603,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, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); + return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection,useEnergyLossCorrection); } //_________________________________________________________________________ @@ -654,7 +667,9 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, const TSpline3* responseFunction, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const + Bool_t usePileupCorrection, + Bool_t useQPtTglCorrection, + Bool_t useEnergyLossCorrection) const { // Calculates the expected sigma of the PID signal as the function of // the information stored in the track and the given parameters, @@ -679,10 +694,10 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, // 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, correctMultiplicity, usePileupCorrection, useQPtTglCorrection) * + return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection, useQPtTglCorrection,useEnergyLossCorrection) * fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints); else - return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection, useQPtTglCorrection)*fRes0[gainScenario]; + return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection, useQPtTglCorrection,useEnergyLossCorrection)*fRes0[gainScenario]; } if (nPoints > 0) { @@ -719,7 +734,9 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, ETPCdEdxSource dedxSource, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const + Bool_t usePileupCorrection, + Bool_t useQPtTglCorrection, + Bool_t useEnergyLossCorrection) const { // Calculates the expected sigma of the PID signal as the function of // the information stored in the track, for the specified particle type @@ -744,7 +761,9 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, ETPCdEdxSource dedxSource, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const + Bool_t usePileupCorrection, + Bool_t useQPtTglCorrection, + Bool_t useEnergyLossCorrection) 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 @@ -758,8 +777,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, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); - Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); + Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection,useEnergyLossCorrection); + Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection,useEnergyLossCorrection); // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters if (sigma >= 998) return -999; @@ -773,7 +792,9 @@ Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track, ETPCdEdxSource dedxSource, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection /*= kFALSE*/, Bool_t useQPtTglCorrection, + Bool_t usePileupCorrection /*= kFALSE*/, + Bool_t useQPtTglCorrection, + Bool_t useEnergyLossCorrection, Bool_t ratio/*=kFALSE*/)const { //Calculates the number of sigmas of the PID signal from the expected value @@ -788,7 +809,7 @@ 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, correctMultiplicity, usePileupCorrection, useQPtTglCorrection); + const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection,useEnergyLossCorrection); Double_t delta=-9999.; if (!ratio) delta=dEdx-bethe; @@ -2406,6 +2427,22 @@ void AliTPCPIDResponse::GetTF1ParametrizationValues(Double_t values[7], const Al } + // energy loss correction - expected mean energy loss divide energy loss at the TPC inner wall + /// calculate relative energy loss correction in the TPC - in helix approximation for primary tracks + /// \param track + /// \return + Double_t AliTPCPIDResponse::GetEnergyLossCorrectionValue(const AliVTrack* track, Int_t pidCode) const{ + Double_t elossCorr=1; + const Float_t rMin=83; + Float_t mass = AliPID::ParticleMass(pidCode); + // makeFitBGRegion(doCheck, "ldEdxRawDist",period, + // "dEdxMeanToInHelix(1/ldEdxRawDist.qPtMean,ldEdxRawDist.tglMean,AliPID::ParticleMass(pidCenter),83,83+tpcNcr3Dist.binMedian,0,5,10,1.4)","dEdxMeanToInHelix(1/qPtMean,tglMean,AliPID::ParticleMass(pidCenter),83,83+tpcNcr3Dist.binMedian,0,5,10,1.4)<1.1",&info); + elossCorr=dEdxMeanToInHelix(track->Pt(),track->GetTPCTgl(),mass,rMin,rMin+track->GetTPCCrossedRows(),fGasType,track->GetBz(),10,fEnergyLossFactor); + return elossCorr; +} + + + /// This is the empirical ALEPH parameterization of the Bethe-Bloch formula. /// Compared to the AliExternalTrackParam::BetheBlochAleph here it is save version not failing in minimization /// \param bg - beta gamma diff --git a/STEER/STEERBase/AliTPCPIDResponse.h b/STEER/STEERBase/AliTPCPIDResponse.h index 7a60b476a4f..604e423a3ac 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.h +++ b/STEER/STEERBase/AliTPCPIDResponse.h @@ -174,21 +174,24 @@ class AliTPCPIDResponse: public TNamed { Bool_t correctEta = kFALSE, Bool_t correctMultiplicity = kFALSE, Bool_t usePileupCorrection = kFALSE, - Bool_t useQPtTglCorrection = kFALSE) const; + Bool_t useQPtTglCorrection = kFALSE, + Bool_t useEnergyLossCorrection = kFALSE) const; Double_t GetExpectedSigma( const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault, Bool_t correctEta = kFALSE, Bool_t correctMultiplicity = kFALSE, Bool_t usePileupCorrection = kFALSE, - Bool_t useQPtTglCorrection = kFALSE) const; + Bool_t useQPtTglCorrection = kFALSE, + Bool_t useEnergyLossCorrection = kFALSE) const; Float_t GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault, Bool_t correctEta = kFALSE, Bool_t correctMultiplicity = kFALSE, Bool_t usePileupCorrection = kFALSE, - Bool_t useQPtTglCorrection = kFALSE) const; + Bool_t useQPtTglCorrection = kFALSE, + Bool_t useEnergyLossCorrection = kFALSE) const; Float_t GetSignalDelta( const AliVTrack* track, AliPID::EParticleType species, @@ -197,6 +200,7 @@ class AliTPCPIDResponse: public TNamed { Bool_t correctMultiplicity = kFALSE, Bool_t usePileupCorrection = kFALSE, Bool_t useQPtTglCorrection = kFALSE, + Bool_t useEnergyLossCorrection = kFALSE, Bool_t ratio = kFALSE) const; void SetResponseFunction(TObject* o, @@ -259,19 +263,22 @@ class AliTPCPIDResponse: public TNamed { void SetEventPileupProperties(Double_t shift, Double_t pileup, Double_t mult) { fEventPileupProperties[0] = shift; fEventPileupProperties[1] = pileup; fEventPileupProperties[2] = mult; } Float_t GetPileUpProperties(UInt_t i ) {return fEventPileupProperties[i%3];} void SetPileupCorrectionStrategy(ETPCPileupCorrectionStrategy strategy) { fPileupCorrectionStrategy = strategy; } - + // pileup correction void SetPileupCorrectionObject(AliNDLocalRegression* correction) { fPileupCorrection = correction; } const AliNDLocalRegression* GetPileupCorrectionObject() const { return fPileupCorrection; } Bool_t IsPileupCorrectionRequested() const { return fPileupCorrectionRequested; } Double_t GetPileupCorrectionValue(const AliVTrack* track) const; static AliNDLocalRegression* GetPileupCorrectionFromFile(const TString fileName); - // + //qPtTgl correction void SetQPtTglCorrectionObject(AliNDLocalRegression* correction) { fQPtTglCorrection = correction; } const AliNDLocalRegression* GetQPtTglCorrectionObject() const { return fQPtTglCorrection; } Bool_t IsQPtTglCorrectionRequested() const { return fQPtTglCorrectionRequested; } Double_t GetQPtTglCorrectionValue(const AliVTrack* track) const; static AliNDLocalRegression* GetQPtTglCorrectionFromFile(const TString fileName); - + // energy loss correction + void SetEnergyLossCorrection(Int_t gasType, Float_t norm){fGasType=gasType; fEnergyLossFactor=norm;} + Bool_t IsEnergyLossCorrectionRequested() const { return fEnergyLossCorrectionRequested; } + Double_t GetEnergyLossCorrectionValue(const AliVTrack* track, Int_t pidCode) const; // //===| dEdx type functions |================================================== void SetdEdxType(ETPCdEdxType dEdxType, Int_t dEdxChargeType=0, Int_t dEdxWeightType=0, Double_t dEdxIROCweight=1., Double_t dEdxOROCmedWeight=1., Double_t dEdxOROClongWeight=1.) { @@ -331,7 +338,9 @@ class AliTPCPIDResponse: public TNamed { const TSpline3* responseFunction, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection, Bool_t useQPtTglCorrection = kFALSE) const; + Bool_t usePileupCorrection, + Bool_t useQPtTglCorrection = kFALSE, + Bool_t useEnergyLossCorrection = kFALSE) const; Double_t GetExpectedSigma(const AliVTrack* track, AliPID::EParticleType species, @@ -341,7 +350,9 @@ class AliTPCPIDResponse: public TNamed { const TSpline3* responseFunction, Bool_t correctEta, Bool_t correctMultiplicity, - Bool_t usePileupCorrection, Bool_t useQPtTglCorrection = kFALSE) const; + Bool_t usePileupCorrection, + Bool_t useQPtTglCorrection = kFALSE, + Bool_t useEnergyLossCorrection = kFALSE) const; // // function for numberical debugging 0 registed splines can be used in the TFormula and tree visualizations // @@ -362,6 +373,9 @@ class AliTPCPIDResponse: public TNamed { AliOADBContainer* fOADBContainer; //! OADB container with response functions AliNDLocalRegression* fPileupCorrection; // pileup correction object AliNDLocalRegression* fQPtTglCorrection; // pileup correction object + Bool_t fEnergyLossCorrectionRequested; // flag fEnergyLossCorrectionRequested + Int_t fGasType; // type of the gas (0-Argon, 1-Neon) + Float_t fEnergyLossFactor; // energy loss facto - fit parameters in respect to tabulated value TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers Float_t fLowGainIROCthreshold; //voltage threshold below which the IROC is considered low gain Float_t fBadIROCthreshhold; //voltage threshold for bad IROCS From e08aa84b322fb3a6f96a7d19f8a738f4ffcf20cd Mon Sep 17 00:00:00 2001 From: miranov25 Date: Thu, 15 Jul 2021 13:25:00 +0200 Subject: [PATCH 13/15] ATO-525 - fix problem reported by Mesut ... protection against rare floating exception stress test at GSI needed ... --- STEER/STEERBase/AliTPCPIDResponse.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/STEER/STEERBase/AliTPCPIDResponse.cxx b/STEER/STEERBase/AliTPCPIDResponse.cxx index ce8a535983f..fbd9eb3a351 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.cxx +++ b/STEER/STEERBase/AliTPCPIDResponse.cxx @@ -2417,7 +2417,7 @@ void AliTPCPIDResponse::GetTF1ParametrizationValues(Double_t values[7], const Al const Double_t bbAlpeh = AliExternalTrackParam::BetheBlochAleph(bg); // values[0] = fMIP / bbAlpeh; // bbAlelp is already normalized to the MIP - values[0] = 1. / bbAlpeh; + values[0] = (bbAlpeh>0) ? 1. / bbAlpeh:1.; values[1] = track->GetTPCTgl(); values[2] = track->GetTPCmomentum(); values[3] = Double_t(species); From c8238b3d7633ba5daf53287ef0292a8669aac885 Mon Sep 17 00:00:00 2001 From: miranov25 Date: Thu, 15 Jul 2021 13:39:56 +0200 Subject: [PATCH 14/15] ATO-525 - fix problem reported by Mesut ... protection against rare floating exception code crashed if the BG<0.005 - at very low momenta ... for Mesut code crashed because not check of the momentum range stress/crash test at GSI to RUN --- STEER/STEERBase/AliTPCPIDResponse.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/STEER/STEERBase/AliTPCPIDResponse.cxx b/STEER/STEERBase/AliTPCPIDResponse.cxx index fbd9eb3a351..bbc183d92f6 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.cxx +++ b/STEER/STEERBase/AliTPCPIDResponse.cxx @@ -2415,7 +2415,7 @@ void AliTPCPIDResponse::GetTF1ParametrizationValues(Double_t values[7], const Al const Double_t p = track->GetTPCmomentum(); const Double_t bg = p / AliPID::ParticleMassZ(Int_t(species)); const Double_t bbAlpeh = AliExternalTrackParam::BetheBlochAleph(bg); - + // bbAleph not well defined below BG 0.01 - putting to the code protection against outliers // values[0] = fMIP / bbAlpeh; // bbAlelp is already normalized to the MIP values[0] = (bbAlpeh>0) ? 1. / bbAlpeh:1.; values[1] = track->GetTPCTgl(); From 2f249b4d1cff4e54e0b17f2b23ddb4d9b5e4bd9e Mon Sep 17 00:00:00 2001 From: miranov25 Date: Tue, 27 Jul 2021 10:39:59 +0200 Subject: [PATCH 15/15] ATO-525 - using dynamic cast for the tree loading - fix problems with non tree keys --- STAT/AliTreePlayer.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/STAT/AliTreePlayer.cxx b/STAT/AliTreePlayer.cxx index 1faa56902b3..a2eb97c4da9 100644 --- a/STAT/AliTreePlayer.cxx +++ b/STAT/AliTreePlayer.cxx @@ -1734,10 +1734,11 @@ TTree* AliTreePlayer::LoadTrees(const char *inputDataList, const char *chRegExp, for (Int_t iKey = 0; iKey < keys->GetEntries(); iKey++) { if (regExp.Match(keys->At(iKey)->GetName()) == 0) continue; // is selected if (notReg.Match(keys->At(iKey)->GetName()) != 0) continue; // is rejected - TTree *tree = (TTree *) finput->Get(keys->At(iKey)->GetName()); // better to use dynamic cast + TTree *tree = dynamic_cast (finput->Get(keys->At(iKey)->GetName())); // better to use dynamic cast + if (tree==NULL) continue; if (treeBase == NULL) { TFile *finput2 = TFile::Open(fileName.Data(),option.Data()); - treeBase = (TTree *) finput2->Get(keys->At(iKey)->GetName()); + treeBase = dynamic_cast (finput2->Get(keys->At(iKey)->GetName())); } TString fileTitle=tagValue["Title"]; if (fileTitle.Length()){