Skip to content
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

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

Open
wants to merge 15 commits into
base: master
Choose a base branch
from
+306 −41
Open

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 1 commit
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
106 changes: 106 additions & 0 deletions STEER/STEERBase/AliTPCPIDResponse.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 (e1<mass) return 0;
return TMath::Sqrt(e1*e1-mass*mass)/TMath::Sqrt(1+tgl*tgl);
}
Comment on lines +2498 to +2520
Copy link
Contributor

Choose a reason for hiding this comment

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

  • Why mix float and double precision?
  • Is double needed? If not I propose to make all float otherwise all double.
  • nearly all values can be const.
  • Why is lengthNorm overwritten? The first value is not used.


/// calcucate dEdx ratio in region 0 and region 0
/// \param ptIn - pt at the inner wall (rIn0) in (GeV)
/// \param tgl - pz/pt
/// \param mass - particle mass
/// \param rIn0 - rIn0 (e.g rIn0 of IROC)
/// \param rOut0 - rOut0 (e.g rOut of OROC)
/// \param rIn1 - rIn1 (e.g rIn0 of OROC medium)
/// \param rOut1 - rIn1 (e.g rOut of OROC medium)
/// \param type - 0-Argon, 1-Neon
/// \param bz - magnetic field (kGaus)
/// \param nStep - number of steps for Euler integration
/// \return
double AliTPCPIDResponse::dEdxRationR0R1Helix(Float_t ptIn, Float_t tgl, Float_t mass, Double_t rIn0,Float_t rOut0, Double_t rIn1,Float_t rOut1, Int_t type, Float_t bz, Int_t nStep){
Copy link
Contributor

Choose a reason for hiding this comment

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

Also in this function, why mit float and double. Make all variables const which don't change.

Float_t pin=ptIn*TMath::Sqrt(1+tgl*tgl);
Double_t rMin=TMath::Min(rIn0,rIn1);
Double_t rMax=TMath::Max(rOut0,rOut1);
Float_t rStep=(rMax-rMin)/nStep;
//
Float_t ptOut=ptIn;
Float_t dEdxSum0=0, dEdxSumN0=0;
Float_t dEdxSum1=0, dEdxSumN1=0;
for (Int_t iStep=0; iStep<nStep; iStep++) {
Float_t radiusM=rMin+rStep*(iStep*0.5);
ptOut = AliTPCPIDResponse::GetPtOutHelix(ptOut, tgl, mass, rMin+rStep*iStep, rMin+rStep*(iStep+1), type, bz);
if (ptOut <= 0) return 1;
Float_t pOut = ptOut * TMath::Sqrt(1 + tgl * tgl);
Float_t dEdxOut = AliExternalTrackParam::BetheBlochAleph(pOut / mass);
if (radiusM>rIn0 && radiusM<rOut0){
dEdxSum0+=dEdxOut;
dEdxSumN0++;
}
if (radiusM>rIn1 && radiusM<rOut1){
dEdxSum1+=dEdxOut;
dEdxSumN1++;
}
}
return (dEdxSum0/dEdxSumN0)/(dEdxSum1/dEdxSumN1);
}



/// return ration of mean dEdx alogn trajectroy to the dEdx at the beginning of trajectory
/// \param ptIn - pt at the input (GeV)
/// \param tgl - pz/pt
/// \param mass - particle mass
/// \param rIn - inner radius
/// \param rOut - outer radius
/// \param type - 0-Argon, 1-Neon
/// \param bz - magnetic field
/// \param nStep - Euler integration approximation with n steps
/// \param useInnerPt - default kFLASE - used for debugging purposes
/// \return
double AliTPCPIDResponse::dEdxMeanToInHelix(Float_t ptIn, Float_t tgl, Float_t mass, Double_t rIn,Float_t rOut, Int_t type, Float_t bz, Int_t nStep, Float_t scale, Bool_t useInnerPt){
Copy link
Contributor

Choose a reason for hiding this comment

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

Also in this function, why mit float and double. Make all variables const which don't change.

Float_t pin=ptIn*TMath::Sqrt(1+tgl*tgl);
Float_t dEdxIn=AliExternalTrackParam::BetheBlochAleph(pin/mass);
Float_t rStep=(rOut-rIn)/nStep;
Float_t ptOut=ptIn;
Float_t dEdxSum=0;
Float_t dEdxSumN=0;
for (Int_t iStep=0; iStep<nStep; iStep++) {
Float_t pt=(useInnerPt) ? ptIn:ptOut;
ptOut = AliTPCPIDResponse::GetPtOutHelix(pt, tgl, mass, rIn+rStep*iStep, rIn+rStep*(iStep+1), type, bz,scale);
if (ptOut <= 0) continue;
Float_t pOut = ptOut * TMath::Sqrt(1 + tgl * tgl);
Float_t dEdxOut = AliExternalTrackParam::BetheBlochAleph(pOut / mass);
dEdxSum+=dEdxOut;
dEdxSumN++;
}
if (dEdxSumN==0) return 1;
return (dEdxSum/dEdxSumN)/dEdxIn;
}
16 changes: 12 additions & 4 deletions STEER/STEERBase/AliTPCPIDResponse.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,27 +173,31 @@ class AliTPCPIDResponse: public TNamed {
ETPCdEdxSource dedxSource = kdEdxDefault,
Bool_t correctEta = kFALSE,
Bool_t correctMultiplicity = kFALSE,
Bool_t usePileupCorrection = kFALSE, Bool_t useQPtTglCorrection = 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, Bool_t useQPtTglCorrection = 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, Bool_t useQPtTglCorrection = kFALSE) const;
Bool_t usePileupCorrection = kFALSE,
Bool_t useQPtTglCorrection = kFALSE) const;

Float_t GetSignalDelta( 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, Bool_t ratio = kFALSE) const;
Bool_t useQPtTglCorrection = kFALSE,
Bool_t ratio = kFALSE) const;

void SetResponseFunction(TObject* o,
AliPID::EParticleType type,
Expand Down Expand Up @@ -316,6 +320,10 @@ class AliTPCPIDResponse: public TNamed {
//
static Double_t BetheBlochAleph(Double_t bg, Double_t kp1, Double_t kp2, Double_t kp3, Double_t kp4, Double_t kp5);
static double sigmadEdxPt(double p, double PID, double resol=0.01 );
static double GetPtOutHelix(Float_t ptIn, Double_t tgl, Float_t mass, Double_t rIn=84, Double_t rOut=245, Int_t type=0, Float_t bz=5,Float_t fraction=1);
static double dEdxRationR0R1Helix(Float_t ptIn, Float_t tgl, Float_t mass, Double_t rIn0=93,Float_t rOut0=245, Double_t rIn1=93,Float_t rOut1=245, Int_t type=0, Float_t bz=5, Int_t nStep=20);
static double dEdxMeanToInHelix(Float_t ptIn, Float_t tgl, Float_t mass, Double_t rIn=93,Float_t rOut=245, Int_t type=0, Float_t bz=5, Int_t nStep=10, Float_t scale=1, Bool_t useInnerPt=kFALSE);

protected:
Double_t GetExpectedSignal(const AliVTrack* track,
AliPID::EParticleType species,
Expand Down