Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-04-04 08:09:26 

0001 #ifndef __DISCRETIAZATIONCORRECTION_H__
0002 #define __DISCRETIAZATIONCORRECTION_H__
0003 
0004 //////////////////////////////////////////////////////////////////////////
0005 //      Discretization Correction for two point correlators     //
0006 //                                  //
0007 //      Skaydi                          //
0008 //      1 Oct 2025                      //
0009 //                                  //
0010 //  discretization correction and jet size corrections      //
0011 //  see attached documentation for derivation and discussion    //
0012 //  Useful to compare calo only to truth/track          //
0013 //  Works for 2pt and integrated 3pt, No idea on full 3     //
0014 //////////////////////////////////////////////////////////////////////////
0015 
0016 #include <math.h>
0017 #include <TH1.h>
0018 #include <TF1.h>
0019 #include <map>
0020 #include <vector>
0021 #include <string>
0022 #include <utility>
0023 
0024 
0025 class DiscretizationCorrection
0026 {   
0027     public:
0028         DiscretizationCorrection(bool Discretize=false, float binsize=0., bool SizeCorr=false, float jetradius=0.)
0029         {
0030             if(Discretize){
0031                 setIsDiscretizationCorrection(Discretize);
0032                 setLatticeSize(binsize);
0033             }
0034             if(SizeCorr){
0035                 setIsSizeCorrection(SizeCorr);
0036                 setJetRadius(jetradius);
0037             }
0038             fillr2s();
0039         }
0040         ~DiscretizationCorrection(){};
0041         
0042         float getLatticeSize(){return this->m_latticesize;}
0043         float getJetRadius(){return this->m_jetradius;}
0044         bool getIsSizeCorrection() { return this->isSizeCorrection;}
0045         bool getIsDiscretizationCorrection(){ return this->isDiscretizationCorrection;}
0046 
0047         void setLatticeSize(float binsize){
0048             this->m_latticesize=binsize;
0049             return;
0050         }
0051         void setJetRadius(float jetradius){
0052             this->m_jetradius=jetradius;
0053             return;
0054         }
0055         void setIsSizeCorrection(bool size){
0056             this->isSizeCorrection=size;
0057             return;
0058         }
0059         void setIsDiscretizationCorrection(bool discr){
0060             this->isDiscretizationCorrection=discr;
0061             return;
0062         }
0063     
0064         float GetCorrectionFactor(float RL){
0065             if(Correction_factors.find(RL) != Correction_factors.end()){
0066                 return Correction_factors[RL];
0067             }
0068             else{
0069                 bool isLess = true, isBigger=false;
0070                 float LessVal= 0., BigVal=0., Val=0., diff=0.;
0071                 for(auto c:Correction_factors)
0072                 {
0073                     if(isBigger && ! isLess ) break;
0074                     if(c.first < RL){
0075                         LessVal=c.second;
0076                         isLess=true;
0077                         isBigger=false;
0078                         diff=RL-c.first;
0079                         continue;
0080                     }
0081                     else if(c.first == RL) 
0082                         return c.second;
0083                     else if(c.first > RL){
0084                         BigVal=c.second;
0085                         isLess=false;
0086                         isBigger=true;
0087                         float diff2=c.first-RL;
0088                         float dr=diff2+diff;
0089                         Val=LessVal*diff/dr+BigVal*diff2/dr;
0090                         return Val;
0091                     }
0092                     else continue;
0093                 }
0094             }
0095         }
0096 
0097         void CaluculateCorrectionFactor(float maxRL=this->sPHENIX_maxRL, int order=0)
0098         {
0099             std::string Plot_title="Efficiency correction factor";
0100             if(isSizeCorrection){
0101                 PlotTitle+=" for R="+std::to_string(m_jetradius)+" jets"
0102             }
0103             else{
0104                 PlotTitle+=" for whole calorimeter";
0105             }
0106             if(isDiscretizationCorrection){
0107                 PlotTitle+=" with bin width #delta = " +std::to_string(m_latticesize)+" in #eta and #varphi ";
0108             }
0109             else
0110             {
0111                 PlotTitle+=" on a continuous lattice ";
0112             }
0113             PlotTitle+="; #Delta R_{L}^{2}; #varepsilon";
0114             Correction_factor_hist=new TH1F("h_corr_fact", PlotTitle.c_str(), (int) std::pow(maxRL/m_latticesize,2), 0, std::pow(maxRL,2)); 
0115             this->correction_order=order;
0116             if(isDiscretizationCorrection)
0117             {
0118                 int nSteps=std::pow(maxRL, 2)/std::pow(m_latticesize,2);
0119                 for(int i=0; i<nSteps; i++)
0120                 {
0121                     float eps=1.;
0122                     float rl=std::sqrt( (float) i );
0123                     if(correction_order==0){
0124                         if(isSizeCorrection)
0125                         {
0126                             eps=all_order_continuous_jet(rl)*std::pow(all_order_discrete_jet(rl), -1);
0127                         }
0128                         eps=eps*all_order_discrete_full(rl)*std::pow(all_order_continuous_full(rl), -1);    
0129                         else{
0130                             if(isSizeCorrection)
0131                             {
0132                                 eps=order_by_order_cont_jet(rl)*std::pow(all_order_discrete_jet(rl), -1);
0133                             }
0134                             eps=eps*all_order_discrete_full(rl)*std::pow(all_order_continuous_full(rl), -1);
0135                         }
0136                     }
0137                     Correction_factors[rl]=eps;
0138                     Correction_factor_hist->Fill(i, eps);
0139                 }
0140             }
0141             else
0142             {
0143                 //no clue how best to do this, this should just be a function???
0144                 continue;
0145             }
0146             return;
0147 
0148         }           
0149             
0150         
0151     private:
0152         
0153         float m_latticesize     =0.;
0154         float m_jetradius       =0.;
0155         bool isSizeCorrection       =false;
0156         bool isDiscretizationCorrection =false;
0157         float sPHENIX_maxRL         =std::sqrt(std::pow(M_PI,2)+std::pow(2.2, 2));
0158 
0159         int correction_order = 3; //how many powers of R_L to expand the elliptic integrals to
0160         
0161         std::map<float, float> Correction_factors;
0162         TH1F* Correction_factor_hist;
0163         float order_by_order_cont_jet(float RL)
0164         {
0165             float ep = R_L*pow(2*M_PI*std::pow(m_jetradius, 2), -1);
0166             float ic = std::sqrt(4*std::pow(m_jetradius, 2)-std::pow(R_L, 2)) + m_jetradius*std::arccos(R_L*std::pow(2*m_jetradius, -1));
0167             if(correction_order >=3 ){
0168                     ic+= -R_L/3.*std::pow(std::arccos(R_L*std::pow(2*m_jetradius, -1)), 3);
0169             }
0170             if(correction_order >= 5 ){
0171                 ic+=R_L*std::pow(120.*m_jetradius, -1)*(4*std::pow(m_jetradius, 2) - std::pow(R_L, 2))*std::pow(std::arccos(R_L*std::pow(2*m_jetradius, -1)), 5);
0172             }
0173             if(correction_order > 5)
0174             {
0175                 for(int i=0; i<=correction_order; i++)
0176                 {
0177                     if(i%2 != 1) continue;
0178                     ic+=0; //getting this expansion  is going to be tricky 
0179                 }
0180             }
0181             float c = ep*ic;
0182             return c;
0183         }
0184         float all_order_continuous_jet(float RL)
0185         {
0186             //This is the solution of the integral for geometric pairs at distance RL where both points lie within the jet
0187             float I_inc=R_L*std::pow(2*M_PI*std::pow(m_jetradius, 2), -1)*std::sqrt(4*std::pow(m_jetradius, 2)-std::pow(RL,2)) * all_order_continous_full(RL);
0188             float elipint=std::ellint_2f(std::arccos(RL*std::pow(2*m_jetradius, -1)), std::pow(RL*std::pow(m_jetradius, -1), 2)) *all_order_continuous_full(RL);
0189             I_inc+=elipint;
0190             return I_inc;
0191         }
0192         float all_order_continuous_full(float RL)
0193         {
0194             //this is the total pairs in the phase space defined by the points lying within the jet and their pairs at distance R_l irrespective of the pair being in the jet radius
0195             float I_Total=std::pow(2*M_PI, 2)*RL*m_jetradius;
0196             return I_Total;
0197         }
0198         float all_order_discrete_jet(float RL) 
0199         {
0200             //pairs with both in jet
0201             float c=0;
0202             int i=floor(std::pow((RL+m_jetradius)/m_latticesize, 2));
0203             for(int j=0; j<=i; j++) 
0204             {
0205                 if(j < std::pow(m_jetradius/m_latticesize, 2)) continue;
0206                 c+=this->r2s.at(j);
0207             }
0208             float ct=all_order_discrete_full(RL) - c;
0209             return ct;
0210         }
0211         float all_order_discrete_full(float RL)
0212         {
0213             //pairs with one in jet
0214             float c=0;
0215             int i=floor(std::pow(RL/m_latticesize, 2));
0216             int counter=0;
0217             for(auto n=this->r2s.begin(); n!=this->r2s.end(); ++n)
0218             {
0219                 c+=*n;
0220                 counter++;
0221                 if(counter >=i) break;
0222             }
0223             c=c*r2s.at(i);
0224             return c;
0225         }
0226         int sum_of_squares(int r) 
0227         {
0228             int r2=0.;
0229             std::vector<int> divs=getDivisors(r);
0230             std::pair<int, int> d1_3=divisors_mod_1_3(divs);
0231             r2=4*(d1_3.first-d1_3.second);
0232             return r2;
0233         }
0234 
0235         std::pair<int,int> divisors_mod_1_3(std::vector<int> divisors)
0236         {
0237             std::pair<int, int> d1_3 {0,0};
0238             for(auto i:divisors)
0239             {
0240                 if(i%4==1) d1_3.first++;
0241                 else if (i%4==3) d1_3.second++;
0242                 else continue;
0243             }
0244             return d1_3;
0245         }
0246 
0247         std::vector<int> getDivisors(int r)
0248         {
0249             std::vector<int> divisors; 
0250             for(int i=1; i<=r; i++)
0251             {
0252                 if(r % i ==0 ) divisors.push_back(i);
0253                 else continue;
0254             }
0255             return divisors;
0256         }
0257         void fillr2s()
0258         {
0259             float max=0.;
0260             if(isSizeCorrection) max=m_jetradius/m_latticesize;
0261             else max=std::sqrt(std::pow(2.2, 2) +std::pow(2*M_PI, 2))/m_latticesize; //set to sPHENIX phase space
0262             for(int i=0; i<std::pow(max, 2); i++)
0263             {
0264                 this->r2s.push_back(sum_of_squares(r));
0265             }
0266             return;
0267         }
0268         std::vector<int> r2s;
0269         TF1 *Continous_jet_f1, *Continuous_full_f1; 
0270         TF1 *Discrete_jet_f1,  *Discete_full_f1;
0271         TF1 *Corretion_funct;
0272         
0273 }
0274 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!