Skip to content

Commit

Permalink
feat: energy loss corrections for proton
Browse files Browse the repository at this point in the history
  • Loading branch information
@c-dilks
c-dilks committed Feb 20, 2024
1 parent 863d390 commit 8366f1b
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 131 deletions.
199 changes: 68 additions & 131 deletions src/iguana/algorithms/clas12/MomentumCorrection.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,10 +44,20 @@ namespace iguana::clas12 {


std::tuple<float,float,float> MomentumCorrection::Transform(float px, float py, float pz, int sec, int pid, float torus) const {
auto cor = torus < 0
? CorrectionInbending(px, py, pz, sec, pid)
: CorrectionOutbending(px, py, pz, sec, pid);
return std::tuple<float,float,float>{cor*px, cor*py, cor*pz};
// energy loss correction
auto e_cor = torus < 0
? EnergyLossInbending(px, py, pz, pid)
: EnergyLossOutbending(px, py, pz, pid);
// momentum correction
auto p_cor = torus < 0
? CorrectionInbending(e_cor * px, e_cor * py, e_cor * pz, sec, pid)
: CorrectionOutbending(e_cor * px, e_cor * py, e_cor * pz, sec, pid);
// return the corrected momentum
return std::tuple<float,float,float>{
e_cor * p_cor * px,
e_cor * p_cor * py,
e_cor * p_cor * pz
};
}


Expand All @@ -64,7 +74,7 @@ namespace iguana::clas12 {
double dp = 0;

// Defining Phi Angle
double Phi = (180/3.1415926)*atan2(Py, Px);
double Phi = (180/M_PI)*atan2(Py, Px);

// (Initial) Shift of the Phi Angle (done to realign sectors whose data is separated when plotted from ±180˚)
if(((sec == 4 || sec == 3) && Phi < 0) || (sec > 4 && Phi < 90)){
Expand Down Expand Up @@ -92,13 +102,11 @@ namespace iguana::clas12 {
phi = PhiLocal - (32/(pp-0.05));
}


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//==================================================================================================================================//
//=======================//=======================// Electron Corrections //=======================//=======================//
//==================================================================================================================================//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

if(pid == 11){
if(sec == 1){
// The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 1] is:
Expand Down Expand Up @@ -126,14 +134,6 @@ namespace iguana::clas12 {
}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//====================================================================================================================================//
//======================//======================// Electron Corrections (End) //======================//======================//
//====================================================================================================================================//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//====================================================================================================================================//
//=========================//=========================// π+ Corrections //=========================//=========================//
Expand Down Expand Up @@ -173,15 +173,6 @@ namespace iguana::clas12 {

}


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//====================================================================================================================================//
//=========================//=========================// π+ Corrections (End) //=========================//=========================//
//====================================================================================================================================//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//====================================================================================================================================//
//==================//==================// π- Corrections (Updated as of 01-13-2023) //==================//==================//
Expand Down Expand Up @@ -225,14 +216,6 @@ namespace iguana::clas12 {
}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//====================================================================================================================================//
//=======================//=======================// π- Corrections (End) //=======================//=======================//
//====================================================================================================================================//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//====================================================================================================================================//
//=======================//=======================// All Proton Corrections //=======================//=======================//
Expand All @@ -259,52 +242,14 @@ namespace iguana::clas12 {
}
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//=====================================================================================================================================//
//=======================//=======================// End of Proton Corrections //=======================//=======================//
//=====================================================================================================================================//
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



return dp/pp;
return dp/pp + 1;
}


/*
* INBENDING USAGE
// The following code is for the Energy Loss Corrections for the proton
double dE_loss = 0;
// Inbending Energy Loss Correction //
if(proth < 27){
dE_loss = exp(-2.739 - 3.932*pro) + 0.002907;
}
if(proth > 27){
dE_loss = exp(-1.2 - 4.228*pro) + 0.007502;
}
double feloss = (pro + dE_loss)/pro;
// Below shows how the corrections are to be applied using the ROOT momentum 4-vector using the above code:
auto fe = dppC(ex, ey, ez, esec, 0) + 1;
auto fpip = dppC(pipx, pipy, pipz, pipsec, 1) + 1;
auto fpim = dppC(pimx, pimy, pimz, pimsec, 2) + 1;
auto fpro = dppC(prox*feloss, proy*feloss, proz*feloss, prosec, 3) + 1;
auto eleC = ROOT::Math::PxPyPzMVector(ex*fe, ey*fe, ez*fe, 0);
auto pipC = ROOT::Math::PxPyPzMVector(pipx*fpip, pipy*fpip, pipz*fpip, 0.13957);
auto pimC = ROOT::Math::PxPyPzMVector(pimx*fpim, pimy*fpim, pimz*fpim, 0.13957);
auto proC = ROOT::Math::PxPyPzMVector(prox*feloss*fpro, proy*feloss*fpro, proz*feloss*fpro, 0.938);
*/




double MomentumCorrection::CorrectionOutbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const {

// skip the correction if it's not defined
if(!( pid == 11 || std::abs(pid) == 211 || pid == 2212 ))
if(!( pid == 11 || std::abs(pid) == 211))
return 1.0;

// Momentum Magnitude
Expand All @@ -314,7 +259,7 @@ namespace iguana::clas12 {
double dp = 0;

// Defining Phi Angle
double Phi = (180/3.1415926)*atan2(Py, Px);
double Phi = (180/M_PI)*atan2(Py, Px);

// (Initial) Shift of the Phi Angle (done to realign sectors whose data is separated when plotted from ±180˚)
if(((sec == 4 || sec == 3) && Phi < 0) || (sec > 4 && Phi < 90)){
Expand Down Expand Up @@ -362,9 +307,6 @@ namespace iguana::clas12 {
dp = ((-6.2488e-08)*phi*phi + (2.23173e-05)*phi + (-0.00227522))*pp*pp + ((1.8372e-05)*phi*phi + (-7.5227e-05)*phi + (0.032636))*pp + ((-6.6566e-05)*phi*phi + (-2.4450e-04)*phi + (-0.072293));
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//======================//======================// Electron Corrections (End) //======================//======================//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//=========================//=========================// π+ Corrections //=========================//=========================//
Expand All @@ -389,77 +331,72 @@ namespace iguana::clas12 {
dp = ((-3.92944e-07)*phi*phi + (1.45848e-05)*phi + (0.00120668))*pp*pp + ((3.7899e-06)*phi*phi + (-1.98219e-04)*phi + (-0.0131312))*pp + ((-3.9961e-06)*phi*phi + (-1.32883e-04)*phi + (0.0294497));
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//=======================//=======================// π+ Corrections (End) //=======================//=======================//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//=======================//=======================// π- Corrections (Start) //=======================//=======================//
//=======================//=======================// π- Corrections //=======================//=======================//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

if(pid == -211){
if(sec == 1){
dp = (7.8044E-06*phi*phi + -9.4703E-05*phi + 4.6696E-03)*pp*pp +(-3.4668E-05*phi*phi + 6.2280E-04*phi + -2.4273E-02)*pp + 2.3566E-05*phi*phi + -5.8519E-04*phi + 3.9226E-02;


}
if(sec == 2){
dp =(-4.6611E-06*phi*phi + -8.1637E-05*phi + 7.5013E-03)*pp*pp +(1.7616E-05*phi*phi + 3.5439E-04*phi + -3.7122E-02)*pp + -1.6286E-05*phi*phi + -2.6545E-04*phi + 4.5659E-02;

}

if(sec == 3){
dp =(4.5270E-06*phi*phi + 2.2578E-04*phi + 5.9214E-03)*pp*pp +(-1.6419E-05*phi*phi + -8.1776E-04*phi + -3.2776E-02)*pp + 1.3734E-05*phi*phi + 6.6125E-04*phi + 4.5784E-02;
}

if(sec == 4){
dp =(-1.3141E-06*phi*phi + 1.9648E-04*phi + 7.6109E-03-0.006)*pp*pp +(8.0912E-06*phi*phi + -8.2672E-04*phi + -4.0495E-02+0.03)*pp + -3.1380E-06*phi*phi + 6.2211E-04*phi + 5.3361E-02-0.04;

}
if(sec == 1){
dp = (7.8044E-06*phi*phi + -9.4703E-05*phi + 4.6696E-03)*pp*pp +(-3.4668E-05*phi*phi + 6.2280E-04*phi + -2.4273E-02)*pp + 2.3566E-05*phi*phi + -5.8519E-04*phi + 3.9226E-02;
}
if(sec == 2){
dp =(-4.6611E-06*phi*phi + -8.1637E-05*phi + 7.5013E-03)*pp*pp +(1.7616E-05*phi*phi + 3.5439E-04*phi + -3.7122E-02)*pp + -1.6286E-05*phi*phi + -2.6545E-04*phi + 4.5659E-02;
}
if(sec == 3){
dp =(4.5270E-06*phi*phi + 2.2578E-04*phi + 5.9214E-03)*pp*pp +(-1.6419E-05*phi*phi + -8.1776E-04*phi + -3.2776E-02)*pp + 1.3734E-05*phi*phi + 6.6125E-04*phi + 4.5784E-02;
}
if(sec == 4){
dp =(-1.3141E-06*phi*phi + 1.9648E-04*phi + 7.6109E-03-0.006)*pp*pp +(8.0912E-06*phi*phi + -8.2672E-04*phi + -4.0495E-02+0.03)*pp + -3.1380E-06*phi*phi + 6.2211E-04*phi + 5.3361E-02-0.04;
}
if(sec == 5){
dp =(-5.4065E-06*phi*phi + -1.6325E-05*phi + 1.2269E-02-0.002)*pp*pp +(1.9512E-05*phi*phi + 1.0228E-04*phi + -6.2351E-02+0.01)*pp + -9.5023E-06*phi*phi + -3.7997E-05*phi + 7.1061E-02-0.02;
}
if(sec == 6){
dp =(-1.1882E-05*phi*phi + 2.0101E-04*phi + 1.1635E-02-0.01)*pp*pp +(5.8488E-05*phi*phi + -6.4709E-04*phi + -5.3833E-02+0.05)*pp + -4.4462E-05*phi*phi + 3.7529E-04*phi + 6.2130E-02-0.06;
}
}

if(sec == 5){
dp =(-5.4065E-06*phi*phi + -1.6325E-05*phi + 1.2269E-02-0.002)*pp*pp +(1.9512E-05*phi*phi + 1.0228E-04*phi + -6.2351E-02+0.01)*pp + -9.5023E-06*phi*phi + -3.7997E-05*phi + 7.1061E-02-0.02;
return dp/pp + 1;
}

}

if(sec == 6){
dp =(-1.1882E-05*phi*phi + 2.0101E-04*phi + 1.1635E-02-0.01)*pp*pp +(5.8488E-05*phi*phi + -6.4709E-04*phi + -5.3833E-02+0.05)*pp + -4.4462E-05*phi*phi + 3.7529E-04*phi + 6.2130E-02-0.06;
double MomentumCorrection::EnergyLossInbending(const float Px, const float Py, const float Pz, const int pid) const {

}
}
// The following code is for the Energy Loss Corrections for the proton
if(pid != 2212)
return 1.0;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//=======================//=======================// π- Corrections (End) //=======================//=======================//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double dE_loss = 0;
auto pro = sqrt(Px*Px + Py*Py + Pz*Pz);
auto proth = atan2(sqrt(Px*Px + Py*Py), Pz)*(180/M_PI);

return dp/pp;
// Inbending Energy Loss Correction //
if(proth < 27){
dE_loss = exp(-2.739 - 3.932*pro) + 0.002907;
}
else {
dE_loss = exp(-1.2 - 4.228*pro) + 0.007502;
}
return (pro + dE_loss)/pro;
}


/*
* OUTBENDING USAGE
double MomentumCorrection::EnergyLossOutbending(const float Px, const float Py, const float Pz, const int pid) const {

// The following code is for the Energy Loss Corrections for the proton
double dE_loss = 0;
// Outbending Energy Loss Correction //
if(proth > 27){
dE_loss = exp(-1.871 - 3.063*pro) + 0.007517;
}
double feloss = (pro + dE_loss)/pro;
// Below shows how the corrections are to be applied using the ROOT momentum 4-vector using the above code:
auto fe = dppC(ex, ey, ez, esec, 0) + 1;
auto fpip = dppC(pipx, pipy, pipz, pipsec, 1) + 1;
auto fpim = dppC(pimx, pimy, pimz, pimsec, 2) + 1;
auto fpro = dppC(prox*feloss, proy*feloss, proz*feloss, prosec, 3) + 1;
// The following code is for the Energy Loss Corrections for the proton
if(pid != 2212)
return 1.0;

auto eleC = ROOT::Math::PxPyPzMVector(ex*fe, ey*fe, ez*fe, 0);
auto pipC = ROOT::Math::PxPyPzMVector(pipx*fpip, pipy*fpip, pipz*fpip, 0.13957);
auto pimC = ROOT::Math::PxPyPzMVector(pimx*fpim, pimy*fpim, pimz*fpim, 0.13957);
auto proC = ROOT::Math::PxPyPzMVector(prox*feloss*fpro, proy*feloss*fpro, proz*feloss*fpro, 0.938);
double dE_loss = 0;
auto pro = sqrt(Px*Px + Py*Py + Pz*Pz);
auto proth = atan2(sqrt(Px*Px + Py*Py), Pz)*(180/M_PI);

*/
// Outbending Energy Loss Correction //
if(proth > 27){
dE_loss = exp(-1.871 - 3.063*pro) + 0.007517;
}
return (pro + dE_loss)/pro;
}


void MomentumCorrection::Stop() {
Expand Down
16 changes: 16 additions & 0 deletions src/iguana/algorithms/clas12/MomentumCorrection.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,22 @@ namespace iguana::clas12 {
/// @returns the correction factor
double CorrectionOutbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const;

/// Energy loss correction for inbending data
/// @param Px @f$p_x@f$
/// @param Py @f$p_y@f$
/// @param Pz @f$p_z@f$
/// @param pid the particle PDG
/// @returns the correction factor
double EnergyLossInbending(const float Px, const float Py, const float Pz, const int pid) const;

/// Energy loss correction for outbending data
/// @param Px @f$p_x@f$
/// @param Py @f$p_y@f$
/// @param Pz @f$p_z@f$
/// @param pid the particle PDG
/// @returns the correction factor
double EnergyLossOutbending(const float Px, const float Py, const float Pz, const int pid) const;

private:

/// `hipo::banklist` index for the particle bank
Expand Down

0 comments on commit 8366f1b

Please sign in to comment.