Skip to content

ATO-525,ATO-496 - addig qptTgl correction to slove dEdx splitting at pt<0.5 GeV #1319

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 15 commits into from
Closed

ATO-525,ATO-496 - addig qptTgl correction to slove dEdx splitting at pt<0.5 GeV #1319

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
a7f8781
ATO-525,ATO-496 - addig qptTgl correction to slove dEdx splitting at …
miranov25 Feb 14, 2021
86bd8aa
ATO-525 - bug fix for ROOT6 - makeCacheTree - have to check upper bound
miranov25 Feb 17, 2021
9b6cb13
ATO-525 - bug fix for ROOT6 - makeCacheTree - have to check upper bound
miranov25 Feb 19, 2021
0ed3280
ATO-525 - make histograms - using long for the size print
miranov25 Mar 22, 2021
50701f1
PWGPP-620 - adding sqrt(1+tgl**2) to the energu loss numerical deriva…
miranov25 Apr 12, 2021
2640966
PWGPP-620 - floating error protection -use 0 as minimal momenta
miranov25 Apr 13, 2021
233e219
PWGPP-620 - bug fix - use sqrt
miranov25 Apr 14, 2021
8858b3e
PWGPP-620 - bug fix - use proper axis i case of the error query
miranov25 May 10, 2021
5e263cd
ATO-525 - changing order of arguments
miranov25 Jul 12, 2021
49e34a1
ATO-525 - using qPtTglCorrection if requested
miranov25 Jul 12, 2021
033a34f
ATO-525 - adding static function for the energy loss correction in TP…
miranov25 Jul 12, 2021
e140dc5
ATO-525 - adding energy loss correction
miranov25 Jul 15, 2021
e08aa84
ATO-525 - fix problem reported by Mesut ... protection against rare f…
miranov25 Jul 15, 2021
c8238b3
ATO-525 - fix problem reported by Mesut ... protection against rare f…
miranov25 Jul 15, 2021
2f249b4
ATO-525 - using dynamic cast for the tree loading - fix problems wit…
miranov25 Jul 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions STAT/AliTreePlayer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
//
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -1647,6 +1647,7 @@ void AliTreePlayer::MakeCacheTreeChunk(TTree * tree, TString varList, TString ou
for (Int_t iEntry=firstEntry; iEntry<nEntries; iEntry+=chunkSize) {
::Info("Processing chunk","%d",iEntry);
if (estimate < chunkSize) tree->SetEstimate(chunkSize);
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) {
Expand Down
28 changes: 21 additions & 7 deletions STAT/TStatToolkit.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;i<count;i++){
graphNew->GetXaxis()->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;i<count;i++){
graphNew->GetYaxis()->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);
Expand Down
18 changes: 13 additions & 5 deletions STEER/ESD/AliESDv0.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = 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;
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) );
Expand Down
92 changes: 78 additions & 14 deletions STEER/STEERBase/AliTPCPIDResponse.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,7 @@ AliTPCPIDResponse::AliTPCPIDResponse():
fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
fOADBContainer(0x0),
fPileupCorrection(0x0),
fQPtTglCorrection(0x0),
fVoltageMap(72),
fLowGainIROCthreshold(-40),
fBadIROCthreshhold(-70),
Expand Down Expand Up @@ -105,6 +106,7 @@ AliTPCPIDResponse::AliTPCPIDResponse():
fOROClongWeight(1.),
fPileupCorrectionStrategy(kPileupCorrectionInExpectedSignal),
fPileupCorrectionRequested(kFALSE),
fQPtTglCorrectionRequested(kFALSE),
fSigmaParametrization(0x0),
fMultiplicityNormalization(1),
fRecoPassNameUsed(),
Expand Down Expand Up @@ -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),
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);
}

//_________________________________________________________________________
Expand Down Expand Up @@ -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,
Expand All @@ -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) {
Expand Down Expand Up @@ -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
Expand All @@ -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);
}


Expand All @@ -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
Expand All @@ -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;
Expand All @@ -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,
Copy link
Contributor

@wiechula wiechula Jul 12, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Bool_t usePileupCorrection /*= kFALSE*/, Bool_t useQPtTglCorrection,
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
Expand All @@ -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;
Expand Down Expand Up @@ -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)
{
Expand Down Expand Up @@ -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<AliNDLocalRegression*>(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"*/,
Expand Down
Loading