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.

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

feat: momentum corrections #104

Merged
merged 9 commits into from
Feb 21, 2024
Merged

feat: momentum corrections #104

merged 9 commits into from
Feb 21, 2024

Conversation

c-dilks
Copy link
Member

@c-dilks c-dilks commented Feb 10, 2024

Adapted from https://clasweb.jlab.org/wiki/index.php/CLAS12_Momentum_Corrections#tab=Correction_Code from @RichCap

Diffs from @RichCap's code to this PR's (click to expand):

Inbending corrections
1c1
< auto dppC = [&](float Px, float Py, float Pz, int sec, int ivec){
---
> double MomentumCorrection::CorrectionInbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const {
3,6c3,5
<   // ivec = 0 --> Electron Corrections
<   // ivec = 1 --> π+ Corrections
<   // ivec = 2 --> π- Corrections
<   // ivec = 3 --> Proton Corrections (NOT UPDATED YET)
---
>   // skip the correction if it's not defined
>   if(!( pid == particle::electron || pid == particle::pi_plus || pid == particle::pi_minus || pid == particle::proton ))
>     return 1.0;
15c14
<   double Phi = (180/3.1415926)*atan2(Py, Px);
---
>   double Phi = (180/M_PI)*atan2(Py, Px);
29c28
<   if(ivec == 0){
---
>   if(pid == particle::electron){
34c33
<   if(ivec == 1 || ivec == 3){
---
>   if(pid == particle::pi_plus || pid == particle::proton){
39c38
<   if(ivec == 2){
---
>   if(pid == particle::pi_minus){
43d41
<
49,50c47
<
<   if(ivec == 0){
---
>   if(pid == particle::electron){
79,86d75
<   //======================//======================//     Electron Corrections (End)     //======================//======================//
<   //====================================================================================================================================//
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<
<
<
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<   //====================================================================================================================================//
90c79
<   if(vec == 1){
---
>   if(pid == particle::pi_plus){
124,132d112
<
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<   //====================================================================================================================================//
<   //=========================//=========================//  π+ Corrections (End)  //=========================//=========================//
<   //====================================================================================================================================//
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<
<
<
138,141c118
<
<
<
<   if(ivec == 2){
---
>   if(pid == particle::pi_minus){
181,188d157
<   //=======================//=======================//      π- Corrections (End)      //=======================//=======================//
<   //====================================================================================================================================//
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<
<
<
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<   //====================================================================================================================================//
192c161
<   if(ivec == 3){
---
>   if(pid == particle::proton){
194c163
<       dp = ((1 + TMath::Sign(1, (pp - 1.4)))/2)*((4.4034e-03)*pp + (-0.01703)) + ((1 + TMath::Sign(1, -(pp - 1.4)))/2)*((-0.10898)*(pp - 1.4)*(pp - 1.4) + (-0.09574)*(pp - 1.4) + ((4.4034e-03)*1.4 + (-0.01703)));
---
>       dp = ((1 + std::copysign(1, (pp - 1.4)))/2)*((4.4034e-03)*pp + (-0.01703)) + ((1 + std::copysign(1, -(pp - 1.4)))/2)*((-0.10898)*(pp - 1.4)*(pp - 1.4) + (-0.09574)*(pp - 1.4) + ((4.4034e-03)*1.4 + (-0.01703)));
197c166
<       dp = ((1 + TMath::Sign(1, (pp - 1.5)))/2)*((0.01318)*pp + (-0.03403)) + ((1 + TMath::Sign(1, -(pp - 1.5)))/2)*((-0.09829)*(pp - 1.5)*(pp - 1.5) + (-0.0986)*(pp - 1.5) + ((0.01318)*1.5 + (-0.03403)));
---
>       dp = ((1 + std::copysign(1, (pp - 1.5)))/2)*((0.01318)*pp + (-0.03403)) + ((1 + std::copysign(1, -(pp - 1.5)))/2)*((-0.09829)*(pp - 1.5)*(pp - 1.5) + (-0.0986)*(pp - 1.5) + ((0.01318)*1.5 + (-0.03403)));
200c169
<       dp = ((1 + TMath::Sign(1, (pp - 1.05)))/2)*((-4.7052e-03)*pp + (1.2410e-03)) + ((1 + TMath::Sign(1, -(pp - 1.05)))/2)*((-0.22721)*(pp - 1.05)*(pp - 1.05) + (-0.09702)*(pp - 1.05) + ((-4.7052e-03)*1.05 + (1.2410e-03)));
---
>       dp = ((1 + std::copysign(1, (pp - 1.05)))/2)*((-4.7052e-03)*pp + (1.2410e-03)) + ((1 + std::copysign(1, -(pp - 1.05)))/2)*((-0.22721)*(pp - 1.05)*(pp - 1.05) + (-0.09702)*(pp - 1.05) + ((-4.7052e-03)*1.05 + (1.2410e-03)));
203c172
<       dp = ((1 + TMath::Sign(1, (pp - 1.4)))/2)*((-1.0900e-03)*pp + (-4.0573e-03)) + ((1 + TMath::Sign(1, -(pp - 1.4)))/2)*((-0.09236)*(pp - 1.4)*(pp - 1.4) + (-0.073)*(pp - 1.4) + ((-1.0900e-03)*1.4 + (-4.0573e-03)));
---
>       dp = ((1 + std::copysign(1, (pp - 1.4)))/2)*((-1.0900e-03)*pp + (-4.0573e-03)) + ((1 + std::copysign(1, -(pp - 1.4)))/2)*((-0.09236)*(pp - 1.4)*(pp - 1.4) + (-0.073)*(pp - 1.4) + ((-1.0900e-03)*1.4 + (-4.0573e-03)));
206c175
<       dp = ((1 + TMath::Sign(1, (pp - 1.5)))/2)*((7.3965e-03)*pp + (-0.02428)) + ((1 + TMath::Sign(1, -(pp - 1.5)))/2)*((-0.09539)*(pp - 1.5)*(pp - 1.5) + (-0.09263)*(pp - 1.5) + ((7.3965e-03)*1.5 + (-0.02428)));
---
>       dp = ((1 + std::copysign(1, (pp - 1.5)))/2)*((7.3965e-03)*pp + (-0.02428)) + ((1 + std::copysign(1, -(pp - 1.5)))/2)*((-0.09539)*(pp - 1.5)*(pp - 1.5) + (-0.09263)*(pp - 1.5) + ((7.3965e-03)*1.5 + (-0.02428)));
209c178
<       dp = ((1 + TMath::Sign(1, (pp - 1.15)))/2)*((-7.6214e-03)*pp + (8.1014e-03)) + ((1 + TMath::Sign(1, -(pp - 1.15)))/2)*((-0.12718)*(pp - 1.15)*(pp - 1.15) + (-0.06626)*(pp - 1.15) + ((-7.6214e-03)*1.15 + (8.1014e-03)));
---
>       dp = ((1 + std::copysign(1, (pp - 1.15)))/2)*((-7.6214e-03)*pp + (8.1014e-03)) + ((1 + std::copysign(1, -(pp - 1.15)))/2)*((-0.12718)*(pp - 1.15)*(pp - 1.15) + (-0.06626)*(pp - 1.15) + ((-7.6214e-03)*1.15 + (8.1014e-03)));
213,219c182,183
<   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<   //=====================================================================================================================================//
<   //=======================//=======================//    End of Proton Corrections    //=======================//=======================//
<   //=====================================================================================================================================//
<   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<
<
---
>   return dp/pp + 1;
> }
221,222d184
<   return dp/pp;
< };
Outbending corrections
1,21c1,5
< auto dppC = [&](float Px, float Py, float Pz, int sec, int ivec, int corEl, int corPip, int corPim, int corPro){
<   // ivec = 0 --> Electron Corrections
<   // ivec = 1 --> π+ Pion Corrections
<   // ivec = 2 --> π- Pion Corrections (Not available yet)
<   // ivec = 3 --> Proton Corrections  (Not available yet)
<
<   // corEl ==> Gives the version of the electron correction (On/Off)
<   // corEl == 0 --> No Correction
<   // corEl != 0 --> Final Extended Version of Corrections
<
<   // corPip ==> Gives the version of the π+ Pion correction (On/Off)
<   // corPip == 0 --> No Correction
<   // corPip != 0 --> Final (Extended) Version of Corrections
<
<   // corPim ==> Gives the version of the π- Pion correction (On/Off)
<   // corPim == 0 --> No Correction
<   // corPim != 0 --> (Not Available Yet)
<
<   // corPro ==> Gives the version of the Proton correction (On/Off)
<   // corPro == 0 --> No Correction
<   // corPro != 0 --> (Not Available Yet)
---
> 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 == particle::electron || pid == particle::pi_plus || pid == particle::pi_minus ))
>     return 1.0;
30c14
<   double Phi = (180/3.1415926)*atan2(Py, Px);
---
>   double Phi = (180/M_PI)*atan2(Py, Px);
43c27
<   if(ivec == 0){
---
>   if(pid == particle::electron){
47c31
<   if(ivec == 1 || ivec == 3){
---
>   if(pid == particle::pi_plus || pid == particle::proton){
51c35
<   if(ivec == 2){
---
>   if(pid == particle::pi_minus){
55,73d38
<   ////////////////////////////////////////////////////////////////////////////////////////////////
<   //===============//===============//     No Corrections     //===============//===============//
<   ////////////////////////////////////////////////////////////////////////////////////////////////
<   if(corEl == 0  && ivec == 0){ // No Electron Correction
<     dp = 0;
<   }
<   if(corPip == 0 && ivec == 1){ // No π+ Correction
<     dp = 0;
<   }
<   if(corPim == 0 && ivec == 2){ // No π- Correction
<     dp = 0;
<   }
<   if(corPro == 0 && ivec == 3){ // No Proton Correction
<     dp = 0;
<   }
<   //////////////////////////////////////////////////////////////////////////////////////////////////
<   //==============//==============//     No Corrections (End)     //==============//==============//
<   //////////////////////////////////////////////////////////////////////////////////////////////////
<
77c42
<   if(ivec == 0 && corEl != 0){
---
>   if(pid == particle::electron){
97,99d61
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<   //======================//======================//     Electron Corrections (End)     //======================//======================//
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
104c66
<   if(ivec == 1 && corPip != 0){
---
>   if(pid == particle::pi_plus){
124,127d85
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<   //=======================//=======================//      π+ Corrections (End)      //=======================//=======================//
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<
130c88
<   //=======================//=======================//      π- Corrections (Start)      //=======================//=======================//
---
>   //=======================//=======================//      π- Corrections            //=======================//=======================//
132,133c90
<
<   if(ivec == 2){
---
>   if(pid == particle::pi_minus){
136,137d92
<
<
141d95
<
143d96
<
147d99
<
150d101
<
152d102
<
155d104
<
157d105
<
160d107
<
164,166c111,112
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
<   //=======================//=======================//      π- Corrections (End)      //=======================//=======================//
<   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
---
>   return dp/pp + 1;
> }
168,169d113
<   return dp/pp;
< };

@c-dilks c-dilks linked an issue Feb 12, 2024 that may be closed by this pull request
@c-dilks c-dilks force-pushed the momentum-corrections branch from 9d99971 to a09ef90 Compare February 20, 2024 22:17
@c-dilks c-dilks marked this pull request as ready for review February 20, 2024 22:36
@c-dilks c-dilks force-pushed the momentum-corrections branch from c163255 to d5af10b Compare February 20, 2024 22:37
@c-dilks
Copy link
Member Author

c-dilks commented Feb 20, 2024

Suppressing gavalian/hipo#49 in UBSAN

@c-dilks c-dilks merged commit 3516e5e into main Feb 21, 2024
55 checks passed
@c-dilks c-dilks deleted the momentum-corrections branch February 21, 2024 00:17
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

add momentum corrections
1 participant