Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-06-01 07:07:13 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Christopher Dilks, Connor Pecar, Duane Byer, Sanghwa Park, Brian Page
0003 
0004 /* NOTE:
0005  * if you make changes, MAINTAIN DOCUMENTATION IN ../doc/kinematics.md
0006  */
0007 
0008 #include "Kinematics.h"
0009 ClassImp(Kinematics)
0010 
0011 Kinematics::Kinematics(
0012     Double_t enEleBeam, /*GeV*/
0013     Double_t enIonBeam, /*GeV*/
0014     Double_t crossAng /*mrad*/    
0015     )
0016 {
0017   srand(time(NULL));
0018   // importing from local python script for ML predictions
0019   // requires tensorflow, energyflow packages installed
0020 #ifdef SIDIS_MLPRED
0021   efnpackage = py::module_::import("testEFlowimport");
0022   pfnimport = efnpackage.attr("eflowPredict");
0023 #endif
0024   // set ion mass
0025   IonMass = ProtonMass();
0026   
0027   // revise crossing angle
0028   crossAng *= 1e-3; // mrad -> rad
0029   crossAng = -1*TMath::Abs(crossAng); // take -1*abs(crossAng) to enforce the correct sign
0030 
0031   // TEST SETTINGS - for debugging calculations ///////////////
0032   // set main frame, used for calculations where there is ambiguity which frame is the correct frame to use
0033   mainFrame = fHeadOn; // fLab, fHeadOn
0034   // set method for determining `vecQ` 4-momentum components for certain recon methods (JB,DA,(e)Sigma)
0035   qComponentsMethod = qQuadratic; // qQuadratic, qHadronic, qElectronic
0036   /////////////////////////////////////////////////////////////
0037   
0038   // set beam 4-momenta
0039   // - electron beam points toward negative z
0040   // - ion beam points toward positive z, negative x
0041   // - crossing angle about y-axis is therefore negative, in right-handed coord. system
0042   // - north is +x, east is +z
0043   Double_t momEleBeam = EMtoP(enEleBeam,ElectronMass());
0044   Double_t momIonBeam = EMtoP(enIonBeam,IonMass);
0045   vecEleBeam.SetPxPyPzE(
0046       0,
0047       0,
0048       -momEleBeam,
0049       enEleBeam
0050       );
0051   vecIonBeam.SetPxPyPzE(
0052       momIonBeam * TMath::Sin(crossAng),
0053       0,
0054       momIonBeam * TMath::Cos(crossAng),
0055       enIonBeam
0056       );
0057   s = (vecEleBeam+vecIonBeam).M2();
0058 
0059   // calculate transformations for head-on frame boost
0060   // - boost lab frame -> c.o.m. frame of proton and ion Beams
0061   BvecBoost = vecEleBeam + vecIonBeam;
0062   Bboost = -1*BvecBoost.BoostVector();
0063   // - boost c.o.m. frame of beams -> back to a frame with energies (nearly) the Original beam energies
0064   OvecBoost.SetXYZT( 0.0, 0.0, BvecBoost[2], BvecBoost[3] );
0065   Oboost = OvecBoost.BoostVector();
0066   // - boost beams to c.o.m. frame of beams
0067   this->BoostToBeamComFrame(vecEleBeam,BvecEleBeam);
0068   this->BoostToBeamComFrame(vecIonBeam,BvecIonBeam);
0069   // - rotation of beams about y to remove x-components
0070   rotAboutY = -TMath::ATan2( BvecIonBeam.Px(), BvecIonBeam.Pz() );
0071   BvecEleBeam.RotateY(rotAboutY);
0072   BvecIonBeam.RotateY(rotAboutY);
0073   // - rotation of beams about x to remove y-components
0074   rotAboutX = TMath::ATan2( BvecIonBeam.Py(), BvecIonBeam.Pz() );
0075 
0076   // default transverse spin (needed for phiS calculation)
0077   tSpin = 1; // +1=spin-up, -1=spin-down
0078   lSpin = 1;
0079 
0080   // default proton polarization
0081   polT = 0.80;
0082   polL = 0.;
0083   polBeam = 0.;
0084 
0085   // random number generator (for asymmetry injection
0086   RNG = std::make_unique<TRandomMixMax>(91874); // (TODO: fixed seed?)
0087 
0088   // reset counters
0089   countPIDsmeared=countPIDtrue=0;
0090 };
0091 
0092 
0093 // ---------------------------------
0094 // calculates 4-momenta components of q and W (`vecQ` and `vecW`) as well as
0095 // derived invariants `W` and `nu`
0096 // - use the quadratic method
0097 // - requires `Q2`, `y`, `Pxh`, `Pyh`
0098 void Kinematics::GetQWNu_quadratic(){
0099 
0100   double f,px,py,pz,pE;
0101   switch(mainFrame) {
0102     case fLab:
0103       f = y*(vecIonBeam.Dot(vecEleBeam));
0104       pz = vecIonBeam.Pz();
0105       py = vecIonBeam.Py();
0106       px = vecIonBeam.Px();
0107       pE = vecIonBeam.E();
0108       break;
0109     case fHeadOn:
0110       f = y*(HvecIonBeam.Dot(HvecEleBeam));
0111       pz = HvecIonBeam.Pz();
0112       py = HvecIonBeam.Py();
0113       px = HvecIonBeam.Px();
0114       pE = HvecIonBeam.E();
0115       break;
0116   };
0117   double hx = Pxh - px;
0118   double hy = Pyh - py;
0119 
0120   double a = 1.0 - (pE*pE)/(pz*pz);
0121   double b = (2*pE/(pz*pz))*(px*hx + py*hy + f);
0122   double c = Q2 - hx*hx - hy*hy - (1/(pz*pz))*pow( (f+px*hx+py*hy) ,2.0);
0123   double disc = b*b - 4*a*c; // discriminant
0124 
0125   double qz1, qz2, qE1, qE2, qE, qz;
0126   if(disc>=0 && TMath::Abs(a)>1e-6) {
0127     qE1 = (-1*b+sqrt(b*b-4*a*c))/(2*a);
0128     qE2 = (-1*b-sqrt(b*b-4*a*c))/(2*a);
0129     qz1 = (-1*f + pE*qE1 - px*hx - py*hy)/(pz);
0130     qz2 = (-1*f + pE*qE2 - px*hx - py*hy)/(pz);
0131 
0132     if(fabs(qE1) < fabs(qE2)){
0133       qE = qE1;
0134       qz = qz1;
0135     }
0136     else{
0137       qE = qE2;
0138       qz = qz2;
0139     }
0140 
0141     vecQ.SetPxPyPzE(hx, hy, qz, qE);
0142     if(mainFrame==fHeadOn) 
0143       this->TransformBackToLabFrame(vecQ,vecQ); // need to be in lab frame for downstream `CalculateHadronKinematics`
0144     vecW = vecIonBeam + vecQ;
0145     W = vecW.M();
0146     Nu = vecIonBeam.Dot(vecQ)/IonMass;
0147 
0148   } else {
0149     // this happens a lot more often if mainFrame==fLab
0150     cerr << "ERROR: in Kinematics::GetQWNu_quadratic, ";
0151     if(isnan(disc))             cerr << "discriminant is NaN";
0152     else if(disc<0)             cerr << "negative discriminant";
0153     else if(TMath::Abs(a)<1e-6) cerr << "zero denominator";
0154     else                        cerr << "unknown reason";
0155     cerr << "; skipping event" << endl;
0156     // cerr << "       p=(" << px << "," << py << "," << pz << "," << pE << ") " << endl;
0157     // cerr << "       a=" << a << endl;
0158     // cerr << "       b=" << b << endl;
0159     // cerr << "       c=" << c << endl;
0160     // cerr << "       disc=" << disc << endl;
0161     // cerr << "       hx=" << hx << endl;
0162     // cerr << "       hy=" << hy << endl;
0163     // cerr << "       Pxh=" << Pxh << endl;
0164     // cerr << "       Pyh=" << Pyh << endl;
0165     // cerr << "       f=" << f << endl;
0166     // cerr << "       y=" << y << endl;
0167     // cerr << "       Q2=" << Q2 << endl;
0168     reconOK = false;
0169   }
0170 };
0171 
0172 
0173 // calculates 4-momenta components of q and W (`vecQ` and `vecW`) as well as
0174 // derived invariants `W` and `nu`
0175 // - use the hadronic sum 4-momentum `hadronSumVec`
0176 void Kinematics::GetQWNu_hadronic(){
0177   switch(mainFrame) {
0178     case fLab:
0179       vecQ = hadronSumVec - vecIonBeam;
0180       vecW = hadronSumVec;
0181       break;
0182     case fHeadOn:
0183       vecQ = hadronSumVec - HvecIonBeam;
0184       vecW = hadronSumVec;
0185       this->TransformBackToLabFrame(vecQ,vecQ); // need to be in lab frame for downstream `CalculateHadronKinematics`
0186       this->TransformBackToLabFrame(vecW,vecW);
0187       break;
0188   }
0189   W = vecW.M();
0190   Nu = vecIonBeam.Dot(vecQ)/IonMass;
0191 };
0192 
0193 
0194 // calculates 4-momenta components of q and W (`vecQ` and `vecW`) as well as
0195 // derived invariants `W` and `nu`
0196 // - use the scattered electron 4-momentum `vecElectron`
0197 void Kinematics::GetQWNu_electronic(){
0198   vecQ = vecEleBeam - vecElectron;
0199   vecW = vecEleBeam + vecIonBeam - vecElectron;
0200   W = vecW.M();
0201   Nu = vecIonBeam.Dot(vecQ) / IonMass;
0202 };
0203 
0204 void Kinematics::GetQWNu_ML(){
0205   hfsinfo.clear();
0206   float pidadj = 0;
0207   if(nHFS >= 2){
0208     std::vector<float> partHold;
0209     for(int i = 0; i < nHFS; i++){
0210       double pidsgn=(hfspid[i]/abs(hfspid[i]));
0211       if(abs(hfspid[i])==211) pidadj = 0.4*pidsgn;
0212       if(abs(hfspid[i])==22) pidadj = 0.2*pidsgn;
0213       if(abs(hfspid[i])==321) pidadj = 0.6*pidsgn;
0214       if(abs(hfspid[i])==2212) pidadj = 0.8*pidsgn;
0215       if(abs(hfspid[i])==11) pidadj = 1.0*pidsgn;      
0216       partHold.push_back(hfspx[i]);
0217       partHold.push_back(hfspy[i]);
0218       partHold.push_back(hfspz[i]);
0219       partHold.push_back(hfsE[i]);
0220       partHold.push_back(pidadj);
0221       hfsinfo.push_back(partHold);
0222       partHold.clear();
0223     }
0224     double Q2ele, Q2DA, Q2JB;
0225     double xele, xDA, xJB;
0226     TLorentzVector vecQEle;
0227     globalinfo.clear();
0228     this->CalculateDISbyElectron();
0229     vecQEle.SetPxPyPzE(vecQ.Px(), vecQ.Py(), vecQ.Pz(), vecQ.E());
0230     Q2ele = Q2;
0231     xele = x;
0232     this->CalculateDISbyDA();
0233     Q2DA = Q2;
0234     xDA = x;
0235     this->CalculateDISbyJB();
0236     Q2JB = Q2;
0237     xJB = x;
0238     if( Q2DA > 0 && Q2DA < 1e4){
0239       globalinfo.push_back(log10(Q2DA));
0240     }
0241     else{
0242       globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0)));
0243     }
0244     if( Q2ele > 0 && Q2ele < 1e4){
0245       globalinfo.push_back(log10(Q2ele));
0246     }
0247     else{
0248       globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0)));
0249     }
0250     if( Q2JB > 0 && Q2JB < 1e4){
0251       globalinfo.push_back(log10(Q2JB));
0252     }
0253     else{
0254       globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0)));
0255     }        
0256     if(xDA>0 && xDA < 1){
0257       globalinfo.push_back(-1*log10(xDA));
0258     }
0259     else{
0260       globalinfo.push_back(-1*log10( (float) (rand()) / (float) (RAND_MAX/1.0)  ));
0261     }
0262     if(xele>0 && xele < 1){
0263       globalinfo.push_back(-1*log10(xele));
0264     }
0265     else{
0266       globalinfo.push_back(-1*log10( (float) (rand()) / (float) (RAND_MAX/1.0)  ));
0267     }
0268     if(xJB>0 && xJB < 1){
0269       globalinfo.push_back(-1*log10(xJB));
0270     }
0271     else{
0272       globalinfo.push_back( -1*log10((float) (rand()) / (float) (RAND_MAX/1.0) ) );
0273     }
0274     globalinfo.push_back(vecQEle.Px());
0275     globalinfo.push_back(vecQEle.Py());
0276     globalinfo.push_back(vecQEle.Pz());
0277     globalinfo.push_back(vecQEle.E());
0278 #ifdef SIDIS_MLPRED
0279     py::object nnoutput = pfnimport(hfsinfo, globalinfo);
0280     std::vector<float> nnvecq = nnoutput.cast<std::vector<float>>();
0281     vecQ.SetPxPyPzE(nnvecq[0],nnvecq[1],nnvecq[2],nnvecq[3]);
0282 #endif
0283   }
0284   else{
0285     this->CalculateDISbyElectron();
0286   }
0287   vecW = vecEleBeam + vecIonBeam - vecElectron; 
0288   W = vecW.M();
0289   Nu = vecIonBeam.Dot(vecQ) / IonMass;
0290   
0291 }
0292 
0293 // ------------------------------------------------------
0294 
0295 
0296 // function to call different reconstruction methods
0297 Bool_t Kinematics::CalculateDIS(TString recmethod){
0298 
0299   reconOK = true;
0300 
0301   // transform to the head-on frame; not needed by all reconstruction methods,
0302   // but best to make sure this is done up front
0303   this->TransformToHeadOnFrame(vecEleBeam,HvecEleBeam);
0304   this->TransformToHeadOnFrame(vecIonBeam,HvecIonBeam);
0305   this->TransformToHeadOnFrame(vecElectron,HvecElectron);
0306   
0307   // calculate primary DIS variables, including Q2,x,y,W,nu
0308   if     (recmethod.CompareTo( "Ele", TString::kIgnoreCase)==0)    { this->CalculateDISbyElectron(); }
0309   else if(recmethod.CompareTo( "DA", TString::kIgnoreCase)==0)     { this->CalculateDISbyDA(); }
0310   else if(recmethod.CompareTo( "JB", TString::kIgnoreCase)==0)     { this->CalculateDISbyJB(); }
0311   else if(recmethod.CompareTo( "Mixed", TString::kIgnoreCase)==0)  { this->CalculateDISbyMixed(); }
0312   else if(recmethod.CompareTo( "Sigma", TString::kIgnoreCase)==0)  { this->CalculateDISbySigma(); }
0313   else if(recmethod.CompareTo( "eSigma", TString::kIgnoreCase)==0) { this->CalculateDISbyeSigma(); }
0314   else if(recmethod.CompareTo( "ML", TString::kIgnoreCase)==0) { this->CalculateDISbyML(); }
0315   else {
0316     cerr << "ERROR: unknown reconstruction method" << endl;
0317     return false;
0318   };
0319 
0320   // calculate SIDIS boost vectors
0321   // - lab frame -> C.o.m. frame of virtual photon and ion
0322   CvecBoost = vecQ + vecIonBeam;
0323   Cboost = -1*CvecBoost.BoostVector();
0324   // - lab frame -> Ion rest frame
0325   IvecBoost = vecIonBeam;
0326   Iboost = -1*IvecBoost.BoostVector();
0327 
0328   // calculate depolarization
0329   // - calculate epsilon, the ratio of longitudinal and transverse photon flux [hep-ph/0611265]
0330   // - these calculations are Lorentz invariant
0331   gamma = 2*ProtonMass()*x / TMath::Sqrt(Q2);
0332   epsilon = ( 1 - y - TMath::Power(gamma*y,2)/4 ) /
0333     ( 1 - y + y*y/2 + TMath::Power(gamma*y,2)/4 );
0334   // - factors A,B,C,V,W (see [hep-ph/0611265] using notation from [1408.5721])
0335   depolA = y*y / (2 - 2*epsilon);
0336   depolB = depolA * epsilon;
0337   depolC = depolA * TMath::Sqrt(1-epsilon*epsilon);
0338   depolV = depolA * TMath::Sqrt(2*epsilon*(1+epsilon));
0339   depolW = depolA * TMath::Sqrt(2*epsilon*(1-epsilon));
0340   // - factor ratios (see [1807.10606] eq. 2.3)
0341   if(depolA==0) depolP1=depolP2=depolP3=depolP4=0;
0342   else {
0343     depolP1 = depolB / depolA;
0344     depolP2 = depolC / depolA;
0345     depolP3 = depolV / depolA;
0346     depolP4 = depolW / depolA;
0347   };
0348 
0349   return reconOK;
0350 };
0351 
0352 // calculate DIS kinematics using scattered electron
0353 // - needs `vecElectron` set
0354 void Kinematics::CalculateDISbyElectron() {
0355   this->GetQWNu_electronic(); // set `vecQ`, `vecW`, `W`, `Nu`
0356   Q2 = -1 * vecQ.M2();
0357   x = Q2 / ( 2 * vecQ.Dot(vecIonBeam) );
0358   y = vecIonBeam.Dot(vecQ) / vecIonBeam.Dot(vecEleBeam);
0359 };
0360 
0361 // calculate DIS kinematics using JB method
0362 void Kinematics::CalculateDISbyJB(){
0363   switch(mainFrame) {
0364     case fLab:    y = sigmah/(2*vecEleBeam.E()); break;
0365     case fHeadOn: y = sigmah/(2*HvecEleBeam.E()); break;
0366   };
0367   Q2 = (Pxh*Pxh + Pyh*Pyh)/(1-y);
0368   x = Q2/(s*y);
0369   switch(qComponentsMethod) {
0370     case qQuadratic:  this->GetQWNu_quadratic(); break;
0371     case qHadronic:   this->GetQWNu_hadronic(); break;
0372     case qElectronic: this->GetQWNu_electronic(); break;
0373   }
0374 };
0375 
0376 // calculate DIS kinematics using DA method
0377 // requires 'vecElectron' set
0378 void Kinematics::CalculateDISbyDA(){
0379   float thetah = acos( (Pxh*Pxh+Pyh*Pyh - sigmah*sigmah)/(Pxh*Pxh+Pyh*Pyh+sigmah*sigmah) );
0380   float thetae;
0381   switch(mainFrame) {
0382     case fLab:
0383       thetae = vecElectron.Theta();
0384       Q2 = 4.0*vecEleBeam.E()*vecEleBeam.E()*sin(thetah)*(1+cos(thetae))/(sin(thetah)+sin(thetae)-sin(thetah+thetae));
0385       break;
0386     case fHeadOn:
0387       thetae = HvecElectron.Theta();
0388       Q2 = 4.0*HvecEleBeam.E()*HvecEleBeam.E()*sin(thetah)*(1+cos(thetae))/(sin(thetah)+sin(thetae)-sin(thetah+thetae));
0389       break;
0390   };
0391   y = (sin(thetae)*(1-cos(thetah)))/(sin(thetah)+sin(thetae)-sin(thetah+thetae));
0392   x = Q2/(s*y);
0393   switch(qComponentsMethod) {
0394     case qQuadratic:  this->GetQWNu_quadratic(); break;
0395     case qHadronic:   this->GetQWNu_hadronic(); break;
0396     case qElectronic: this->GetQWNu_electronic(); break;
0397   }
0398 };
0399 
0400 // calculate DIS kinematics using mixed method
0401 // requires 'vecElectron' set
0402 void Kinematics::CalculateDISbyMixed(){
0403   this->GetQWNu_electronic();
0404   Q2 = -1*vecQ.M2();
0405   switch(mainFrame) {
0406     case fLab:    y = sigmah/(2*vecEleBeam.E()); break;
0407     case fHeadOn: y = sigmah/(2*HvecEleBeam.E()); break;
0408   }
0409   x = Q2/(s*y);
0410 };
0411 
0412 // calculate DIS kinematics using Sigma method
0413 // requires 'vecElectron' set
0414 void Kinematics::CalculateDISbySigma(){
0415   switch(mainFrame) {
0416     case fLab:
0417       y = sigmah/(sigmah + vecElectron.E()*(1-cos(vecElectron.Theta())));
0418       Q2 = (vecElectron.Px()*vecElectron.Px() + vecElectron.Py()*vecElectron.Py())/(1-y);
0419       break;
0420     case fHeadOn:
0421       y = sigmah/(sigmah + HvecElectron.E()*(1-cos(HvecElectron.Theta())));
0422       Q2 = (HvecElectron.Px()*HvecElectron.Px() + HvecElectron.Py()*HvecElectron.Py())/(1-y);
0423       break;
0424   }
0425   x = Q2/(s*y);
0426   switch(qComponentsMethod) {
0427     case qQuadratic:  this->GetQWNu_quadratic(); break;
0428     case qHadronic:   this->GetQWNu_hadronic(); break;
0429     case qElectronic: this->GetQWNu_electronic(); break;
0430   }
0431 };
0432 
0433 // calculate DIS kinematics using eSigma method
0434 // requires 'vecElectron' set
0435 void Kinematics::CalculateDISbyeSigma(){
0436   this->CalculateDISbySigma();
0437   Double_t xsigma = x;
0438   Double_t ysigma = y; 
0439   vecQ = vecEleBeam - vecElectron;
0440   Q2 = -1 * vecQ.M2();
0441   y = Q2/(s*xsigma);
0442   x = xsigma;
0443   switch(qComponentsMethod) {
0444     case qQuadratic:  this->GetQWNu_quadratic(); break;
0445     case qHadronic:   this->GetQWNu_hadronic(); break;
0446     case qElectronic: this->GetQWNu_electronic(); break;
0447   }
0448 };
0449 
0450 void Kinematics::CalculateDISbyML() {
0451   this->GetQWNu_ML(); // set `vecQ`, `vecW`, `W`, `Nu`   
0452   Q2 = -1 * vecQ.M2();
0453   x = Q2 / ( 2 * vecQ.Dot(vecIonBeam) );
0454   y = vecIonBeam.Dot(vecQ) / vecIonBeam.Dot(vecEleBeam);
0455 };
0456 
0457 
0458 // calculate hadron kinematics
0459 // - calculate DIS kinematics first, so we have `vecQ`, etc.
0460 // - needs `vecHadron` set
0461 void Kinematics::CalculateHadronKinematics() {
0462   // hadron momentum
0463   pLab = vecHadron.P();
0464   pTlab = vecHadron.Pt();
0465   phiLab = vecHadron.Phi();
0466   etaLab = vecHadron.Eta();
0467   // hadron z
0468   z = vecIonBeam.Dot(vecHadron) / vecIonBeam.Dot(vecQ);
0469   // missing mass
0470   mX = (vecW-vecHadron).M(); // missing mass
0471   // boosts
0472   this->BoostToComFrame(vecHadron,CvecHadron);
0473   this->BoostToComFrame(vecQ,CvecQ);
0474   this->BoostToIonFrame(vecHadron,IvecHadron);
0475   this->BoostToIonFrame(vecQ,IvecQ);
0476   this->BoostToIonFrame(vecElectron,IvecElectron);
0477   // feynman-x: calculated in photon+ion c.o.m. frame
0478   xF = 2 * CvecHadron.Vect().Dot(CvecQ.Vect()) /
0479       (W * CvecQ.Vect().Mag());
0480   // phiH: calculated in ion rest frame
0481   phiH = AdjAngle(PlaneAngle(
0482       IvecQ.Vect(), IvecElectron.Vect(),
0483       IvecQ.Vect(), IvecHadron.Vect()
0484       ));
0485   // phiS: calculated in ion rest frame
0486   tSpin = RNG->Uniform() < 0.5 ? 1 : -1;
0487   lSpin = RNG->Uniform() < 0.5 ? 1 : -1;
0488   vecSpin.SetXYZT(0,1,0,0); // Pauli-Lubanski pseudovector, in lab frame
0489   this->BoostToIonFrame(vecSpin,IvecSpin); // boost to ion rest frame
0490   phiS = AdjAngle(PlaneAngle(
0491       IvecQ.Vect(), IvecElectron.Vect(),
0492       IvecQ.Vect(), IvecSpin.Vect()
0493       ));
0494   // pT, in perp frame (transverse to q): calculated in ion rest frame
0495   pT = Reject(
0496       IvecHadron.Vect(),
0497       IvecQ.Vect()
0498       ).Mag();
0499   // qT
0500   qT = pT / z;
0501 };
0502 
0503 // validate transformations to the head-on frame
0504 void Kinematics::ValidateHeadOnFrame() {
0505   this->BoostToIonFrame(vecEleBeam,IvecEleBeam);
0506   this->BoostToIonFrame(vecIonBeam,IvecIonBeam);
0507   this->TransformToHeadOnFrame(vecIonBeam,HvecIonBeam);
0508   this->TransformToHeadOnFrame(vecEleBeam,HvecEleBeam);
0509   this->TransformToHeadOnFrame(vecIonBeam,HvecIonBeam);
0510   this->TransformToHeadOnFrame(vecElectron,HvecElectron);
0511   this->TransformToHeadOnFrame(vecHadron,HvecHadron);
0512   printf("\nVALIDATION:\n");
0513   printf("lab E:     "); vecEleBeam.Print();
0514   printf("lab I:     "); vecIonBeam.Print();
0515   printf("ion RF  E: "); IvecEleBeam.Print();
0516   printf("ion RF  I: "); IvecIonBeam.Print();
0517   printf("head-on E: "); HvecEleBeam.Print();
0518   printf("head-on I: "); HvecIonBeam.Print();
0519   printf("---\n");
0520   printf("lab electron:     "); vecElectron.Print();
0521   printf("head-on electron: "); HvecElectron.Print();
0522   printf("difference:       "); (vecElectron-HvecElectron).Print();
0523   printf("---\n");
0524   printf("lab hadron:     "); vecHadron.Print();
0525   printf("head-on hadron: "); HvecHadron.Print();
0526   printf("difference:     "); (vecHadron-HvecHadron).Print();
0527 };
0528 
0529 
0530 // add a 4-momentum to the hadronic final state
0531 void Kinematics::AddToHFS(TLorentzVector p4_) {
0532   TLorentzVector p4 = p4_;
0533   if(mainFrame==fHeadOn) this->TransformToHeadOnFrame(p4,p4);
0534   sigmah += (p4.E() - p4.Pz());
0535   Pxh += p4.Px();
0536   Pyh += p4.Py();
0537   hadronSumVec += p4;
0538   countHadrons++;
0539 };
0540 
0541 // add  information from a reconstructed particle to HFSTree
0542 // branches containing full HFS
0543 void Kinematics::AddToHFSTree(TLorentzVector p4, int pid) {
0544   hfspx.push_back(p4.Px());
0545   hfspy.push_back(p4.Py());
0546   hfspz.push_back(p4.Pz());
0547   hfsE.push_back(p4.E());
0548   hfspid.push_back(pid);
0549 
0550   nHFS++;
0551 };
0552 
0553 // add a specific track and its matched true four momentum
0554 // to HFSTree for final kinematic calculations (not full HFS)
0555 void Kinematics::AddTrackToHFSTree(TLorentzVector p4, int pid) {
0556   selectedHadPx.push_back(p4.Px());
0557   selectedHadPy.push_back(p4.Py());
0558   selectedHadPz.push_back(p4.Pz());
0559   selectedHadE.push_back(p4.E());
0560   selectedHadPID.push_back(pid);  
0561 };
0562 
0563 // subtract electron from hadronic final state variables
0564 void Kinematics::SubtractElectronFromHFS() {
0565   if(!isnan(vecElectron.E())){
0566     switch(mainFrame) {
0567       case fLab:
0568         sigmah -= (vecElectron.E() - vecElectron.Pz());
0569         Pxh -= vecElectron.Px();
0570         Pyh -= vecElectron.Py();
0571         hadronSumVec -= vecElectron;
0572         break;
0573       case fHeadOn:
0574         this->TransformToHeadOnFrame(vecElectron,HvecElectron);
0575         sigmah -= (HvecElectron.E() - HvecElectron.Pz());
0576         Pxh -= HvecElectron.Px();
0577         Pyh -= HvecElectron.Py();
0578         hadronSumVec -= HvecElectron;
0579         break;
0580     }
0581     countHadrons--;    
0582   } else {
0583     cerr << "ERROR: electron energy is NaN" << endl;
0584     // TODO: kill event
0585   }
0586 };
0587 
0588 
0589 // reset some variables for the hadronic final state
0590 void Kinematics::ResetHFS() {  
0591   sigmah = Pxh = Pyh = 0;
0592   hadronSumVec.SetPxPyPzE(0,0,0,0);
0593   countHadrons = 0;
0594   nHFS = 0;
0595   
0596   hfspx.clear();
0597   hfspy.clear();
0598   hfspz.clear();
0599   hfsE.clear();
0600   hfspid.clear();
0601   
0602   selectedHadPx.clear();
0603   selectedHadPy.clear();
0604   selectedHadPz.clear();
0605   selectedHadE.clear();
0606   selectedHadPID.clear();    
0607 };
0608 
0609 
0610 // DELPHES-only methods //////////////////////////////////////////////////////
0611 #ifndef EXCLUDE_DELPHES
0612 // calculates reconstructed hadronic final state variables from DELPHES tree branches
0613 // expects 'vecElectron' set
0614 // - calculates `sigmah`, `Pxh`, and `Pyh` in the lab and head-on frames
0615 void Kinematics::GetHFS(
0616     TObjArrayIter itTrack,
0617     TObjArrayIter itEFlowTrack,
0618     TObjArrayIter itEFlowPhoton,
0619     TObjArrayIter itEFlowNeutralHadron,
0620     TObjArrayIter itpfRICHTrack,
0621     TObjArrayIter itDIRCepidTrack,   TObjArrayIter itDIRChpidTrack,
0622     TObjArrayIter itBTOFepidTrack,   TObjArrayIter itBTOFhpidTrack,
0623     TObjArrayIter itdualRICHagTrack, TObjArrayIter itdualRICHcfTrack
0624     ) {
0625 
0626   // resets
0627   this->ResetHFS();
0628   itTrack.Reset();
0629   itEFlowTrack.Reset();
0630   itEFlowPhoton.Reset();
0631   itEFlowNeutralHadron.Reset();
0632 
0633   // track loop
0634   while(Track *track = (Track*)itTrack() ){
0635     TLorentzVector  trackp4 = track->P4();
0636     if(!isnan(trackp4.E())){
0637       if( std::abs(track->Eta) < 4.0  ){
0638 
0639         int pid = GetTrackPID( // get smeared PID
0640             track,
0641             itpfRICHTrack,
0642             itDIRCepidTrack, itDIRChpidTrack,
0643             itBTOFepidTrack, itBTOFhpidTrack,
0644             itdualRICHagTrack, itdualRICHcfTrack
0645             );
0646 
0647         if(pid != -1){ // if smeared PID determined, set mass of `trackp4` accordingly
0648           trackp4.SetPtEtaPhiM(trackp4.Pt(),trackp4.Eta(),trackp4.Phi(),correctMass(pid));
0649           countPIDsmeared++;
0650         }
0651         else { // otherwise, assume true PID
0652           countPIDtrue++;
0653           //continue; // drop events if PID not smeared // TODO [critical]: this is more realistic, but resolutions are much worse
0654         }
0655 
0656         this->AddToHFS(trackp4);
0657     this->AddToHFSTree(trackp4,pid);
0658       }
0659     }    
0660   }
0661 
0662   // eflow high |eta| track loop
0663   while(Track *eflowTrack = (Track*)itEFlowTrack() ){
0664     TLorentzVector eflowTrackp4 = eflowTrack->P4();
0665     int pid = GetTrackPID( // get smeared PID                                                                     
0666             eflowTrack,
0667             itpfRICHTrack,
0668             itDIRCepidTrack, itDIRChpidTrack,
0669             itBTOFepidTrack, itBTOFhpidTrack,
0670             itdualRICHagTrack, itdualRICHcfTrack
0671             );
0672 
0673     if(!isnan(eflowTrackp4.E())){
0674       if(std::abs(eflowTrack->Eta) >= 4.0){
0675         this->AddToHFS(eflowTrackp4);
0676     this->AddToHFSTree(eflowTrackp4, pid);
0677       }
0678     }
0679   }
0680   
0681   // eflow photon loop
0682   while(Tower* towerPhoton = (Tower*)itEFlowPhoton() ){
0683     TLorentzVector  towerPhotonp4 = towerPhoton->P4();
0684     if(!isnan(towerPhotonp4.E())){
0685       if( std::abs(towerPhoton->Eta) < 4.0  ){
0686         this->AddToHFS(towerPhotonp4);
0687     this->AddToHFSTree(towerPhotonp4,22);
0688       }
0689     }
0690   }
0691 
0692   // eflow neutral hadron loop
0693   while(Tower* towerNeutralHadron = (Tower*)itEFlowNeutralHadron() ){
0694     TLorentzVector  towerNeutralHadronp4 = towerNeutralHadron->P4();
0695     if(!isnan(towerNeutralHadronp4.E())){
0696       if( std::abs(towerNeutralHadron->Eta) < 4.0 ){
0697         this->AddToHFS(towerNeutralHadronp4);
0698     this->AddToHFSTree(towerNeutralHadronp4,-1);
0699       }
0700     }
0701   }
0702   
0703   // remove electron from hadronic final state
0704   this->SubtractElectronFromHFS();
0705 };
0706 
0707 
0708 // calculates generated truth hadronic final state variables from DELPHES tree branches
0709 void Kinematics::GetTrueHFS(TObjArrayIter itParticle){
0710 
0711   // resets
0712   this->ResetHFS();
0713   itParticle.Reset();
0714 
0715   // truth loop
0716   while(GenParticle *partTrue = (GenParticle*)itParticle() ) {
0717     if(partTrue->Status == 1) this->AddToHFS(partTrue->P4());
0718   }
0719 
0720   // remove electron from hadronic final state
0721   this->SubtractElectronFromHFS();
0722 };
0723 
0724 
0725 // get PID information from PID systems tracks
0726 int Kinematics::GetTrackPID(
0727     Track *track,
0728     TObjArrayIter itpfRICHTrack,
0729     TObjArrayIter itDIRCepidTrack, TObjArrayIter itDIRChpidTrack,
0730     TObjArrayIter itBTOFepidTrack, TObjArrayIter itBTOFhpidTrack,
0731     TObjArrayIter itdualRICHagTrack, TObjArrayIter itdualRICHcfTrack
0732     ) {
0733 
0734   itpfRICHTrack.Reset();
0735   itDIRCepidTrack.Reset();   itDIRChpidTrack.Reset();
0736   itBTOFepidTrack.Reset();   itBTOFhpidTrack.Reset();
0737   itdualRICHagTrack.Reset(); itdualRICHcfTrack.Reset();
0738   GenParticle *trackParticle = (GenParticle*)track->Particle.GetObject();
0739   GenParticle *detectorParticle;
0740 
0741   // TODO: make this less repetitive:
0742 
0743   while(Track *detectorTrack = (Track*)itpfRICHTrack() ){
0744     detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0745     if( detectorParticle == trackParticle ) return detectorTrack->PID;
0746   }
0747 
0748   while(Track *detectorTrack = (Track*)itDIRCepidTrack() ){
0749     detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0750     if( detectorParticle == trackParticle ) return detectorTrack->PID;
0751   }
0752   while(Track *detectorTrack = (Track*)itDIRChpidTrack() ){
0753     detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0754     if( detectorParticle == trackParticle ) return detectorTrack->PID;
0755   }
0756 
0757   while(Track *detectorTrack = (Track*)itBTOFepidTrack() ){
0758     detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0759     if( detectorParticle == trackParticle ) return detectorTrack->PID;
0760   }
0761   while(Track *detectorTrack = (Track*)itBTOFhpidTrack() ){
0762     detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0763     if( detectorParticle == trackParticle ) return detectorTrack->PID;
0764   }
0765 
0766   while(Track *detectorTrack = (Track*)itdualRICHagTrack() ){
0767     detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0768     if( detectorParticle == trackParticle ) return detectorTrack->PID;
0769   }
0770   while(Track *detectorTrack = (Track*)itdualRICHcfTrack() ){
0771     detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0772     if( detectorParticle == trackParticle ) return detectorTrack->PID;
0773   }
0774 
0775   return -1; // not found
0776 }
0777 
0778 #endif // ifndef EXCLUDE_DELPHES
0779 // end DELPHES-only methods //////////////////////////////////////////////////////
0780 
0781 
0782 // BOOSTS
0783 /////////////////
0784 
0785 // boost from Lab frame `Lvec` to photon+ion C.o.m. frame `Cvec`
0786 void Kinematics::BoostToComFrame(TLorentzVector Lvec, TLorentzVector &Cvec) {
0787   Cvec=Lvec;
0788   Cvec.Boost(Cboost);
0789 };
0790 
0791 // boost from Lab frame `Lvec` to Ion rest frame `Ivec`
0792 void Kinematics::BoostToIonFrame(TLorentzVector Lvec, TLorentzVector &Ivec) {
0793   Ivec=Lvec;
0794   Ivec.Boost(Iboost);
0795 };
0796 
0797 // boost from Lab frame `Lvec` to ion+electron Beam c.o.m. frame `Bvec`
0798 void Kinematics::BoostToBeamComFrame(TLorentzVector Lvec, TLorentzVector &Bvec) {
0799   Bvec=Lvec;
0800   Bvec.Boost(Bboost);
0801 };
0802 
0803 // transform from Lab frame `Lvec` to Head-on frame `Hvec`
0804 void Kinematics::TransformToHeadOnFrame(TLorentzVector Lvec, TLorentzVector &Hvec) {
0805   this->BoostToBeamComFrame(Lvec,Hvec); // boost to c.o.m. frame of beams
0806   Hvec.RotateY(rotAboutY); // remove x-component of beams
0807   Hvec.RotateX(rotAboutX); // remove y-component of beams
0808   Hvec.Boost(Oboost); // return to frame where beam energies are (nearly) the original
0809 };
0810 
0811 // transform from Head-on frame `Hvec` back to Lab frame `Lvec`
0812 void Kinematics::TransformBackToLabFrame(TLorentzVector Hvec, TLorentzVector &Lvec) {
0813   Lvec=Hvec;
0814   Lvec.Boost(-1*Oboost); // revert boosts and rotations
0815   Lvec.RotateX(-rotAboutX);
0816   Lvec.RotateY(-rotAboutY);
0817   Lvec.Boost(-1*Bboost); // boost from c.o.m. frame of beams back to lab frame
0818 };
0819 
0820 
0821 // test a fake asymmetry, for fit code validation
0822 // - assigns `tSpin` based on desired fake asymmetry
0823 void Kinematics::InjectFakeAsymmetry() {
0824   // modulations, including depolarization factors [1807.10606 eq. 2.2-3]
0825   moduVal[0] = TMath::Sin(phiH-phiS); // sivers
0826   moduVal[1] = TMath::Sin(phiH+phiS) * depolP1; // transversity*collins
0827   // fake amplitudes
0828   ampVal[0] = 0.1;
0829   ampVal[1] = 0.1;
0830   // fake dependence on SIDIS kinematics (linear in x)
0831   asymInject = 0;
0832   asymInject +=  ampVal[0]/0.2 * x * moduVal[0];
0833   asymInject += -ampVal[1]/0.2 * x * moduVal[1];
0834   //asymInject = ampVal[0]*moduVal[0] + ampVal[1]*moduVal[1]; // constant
0835   // apply polarization
0836   asymInject *= polT;
0837   // generate random number in [0,1]
0838   RN = RNG->Uniform();
0839   tSpin = (RN<0.5*(1+asymInject)) ? 1 : -1;
0840 };
0841 
0842 
0843 Kinematics::~Kinematics() {
0844 };
0845
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.5 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!