Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-18 08:17:15 

0001 ///*--------------------------------------------------*/
0002 /// eic.cc:
0003 /// Original author: Dr. Ahmed Zafar
0004 /// Date: 2015-2018
0005 ///
0006 ///*--------------------------------------------------*/
0007 /// Modifier: Wenliang (Bill) Li
0008 /// Date: Feb 24 2020
0009 /// Email: wenliang.billlee@gmail.com
0010 ///
0011 /// Comment: Feb 24, 2020: the main function is excuted in main.cc
0012 
0013 #include "eic.h"
0014 
0015 using std::setw;
0016 using std::setprecision;
0017 using std::cout;
0018 using std::cin;
0019 using std::endl;
0020 using std::vector;
0021 using namespace std;
0022 
0023 extern char* DEMPgen_Path;
0024 
0025 //---------------------------------------------------------
0026 // g++ -o pim pimFermi.C `root-config --cflags --glibs`
0027 //---------------------------------------------------------
0028 
0029 //int main() {
0030 // 
0031 //  eic();
0032 //  
0033 //  return 0;
0034 //}
0035 
0036 void eic() {
0037 
0038   Int_t target_direction, kinematics_type;
0039   Double_t EBeam, HBeam;
0040  
0041   cout << "Target Direction (1->Up, 2->Down): "; cin >> target_direction; cout << endl;
0042   cout << "Kinematics type (1->FF, 2->TSSA): ";  cin >> kinematics_type;  cout << endl;
0043   cout << "Enter the number of events: ";        cin >> fNEvents;         cout << endl;
0044   cout << "Enter the file number: ";             cin >> fNFile;           cout << endl;
0045   cout << "Enter the electron beam energy: ";    cin >> EBeam;            cout << endl;
0046   cout << "Enter the hadron beam energy: ";      cin >> HBeam;            cout << endl;
0047  
0048   //    eic(target_direction, kinematics_type, fNEvents);
0049 
0050 }
0051 
0052 /*--------------------------------------------------*/
0053 // 18/01/23 - SJDK- This function is never used since eic() is only called with a json object as the argument. Commented out for now, delete later?
0054 /* void eic(int event_number, int target_direction, int kinematics_type, TString file_name, int fEIC_seed, TString Ejectile, TString RecoilHadron, TString det_location, TString OutputType, double EBeam, double HBeam) {
0055 
0056    TString targetname;
0057    TString charge;
0058 
0059    if( target_direction == 1 ) targetname = "up";
0060    if( target_direction == 2 ) targetname = "down";
0061     
0062    gKinematics_type = kinematics_type;
0063    gfile_name = file_name;
0064 
0065    fNFile = 1;
0066 
0067    fNEvents = event_number;
0068 
0069    fSeed = fEIC_seed;
0070    cout << EBeam << " elec " << HBeam << " RecoilHadrons" << endl; 
0071    fEBeam = EBeam;
0072    fPBeam = HBeam;
0073 
0074    pim* myPim = new pim(fSeed);
0075    myPim->Initilize();
0076    // 09/02/22 - SJDK - Special case for the kaon, if RecoilHadron not specified, default to Lambda
0077    if (Ejectile == "K+"){
0078    if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){
0079    RecoilHadron = "Lambda";
0080    }
0081    else{
0082    RecoilHadron = ExtractParticle(RecoilHadron);
0083    }
0084    Reaction* r1 = new Reaction(Ejectile, RecoilHadron);
0085    r1->process_reaction();
0086    delete r1;
0087    }
0088    else if (Ejectile == "pi+" || Ejectile == "Pion+" ||  Ejectile == "Pi+"){
0089    RecoilHadron = "Neutron";
0090    Ejectile = ExtractParticle(particle);
0091    charge = ExtractCharge(Ejectile);
0092    Reaction* r1 = new Reaction(Ejectile, RecoilHadron);
0093    r1->process_reaction();
0094    delete r1;
0095    }
0096    else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){
0097    RecoilHadron = "Proton";
0098    Ejectile = ExtractParticle(Ejectile);
0099    charge = ExtractCharge(Ejectile);
0100    //Reaction* r1 = new Reaction(Ejectile);
0101    Reaction* r1 = new Reaction(Ejectile, RecoilHadron);
0102    r1->process_reaction();
0103    delete r1;
0104    }
0105    else{
0106    Ejectile = ExtractParticle(Ejectile);
0107    charge = ExtractCharge(Ejectile);
0108    Reaction* r1 = new Reaction(Ejectile);
0109    r1->process_reaction();
0110    delete r1;
0111    }
0112    }
0113 */
0114 
0115 /*--------------------------------------------------*/
0116 /*--------------------------------------------------*/
0117 // SJDK 21/12/22 - Note that this is the one that actually gets used, reads in the .json file
0118 void eic(Json::Value obj) {
0119     
0120   TString targetname;  
0121   TString charge;
0122 
0123   int target_direction = obj["Targ_dir"].asInt();
0124   gKinematics_type     = obj["Kinematics_type"].asInt();
0125 
0126   if( target_direction == 1 ) targetname = "up";
0127   if( target_direction == 2 ) targetname = "down";
0128  
0129   gfile_name = obj["file_name"].asString();
0130  
0131   gPi0_decay = obj["pi0_decay"].asBool();
0132 
0133   fNFile = 1;
0134   fNEvents = obj["n_events"].asUInt64();
0135 
0136   fSeed = obj["generator_seed"].asInt();
0137 
0138   pim* myPim = new pim(fSeed);
0139   myPim->Initilize();
0140  
0141   //    TDatime dsTime;
0142   //    cout << "Start Time:   " << dsTime.GetHour() << ":" << dsTime.GetMinute() << endl;
0143   // 21/12/22 - SJDK - Should do a check if these are defined or not, should crash if not defined or set defaults, see other quantities below
0144   TString ROOTFile = obj["ROOTOut"].asString();
0145   if (ROOTFile == "True" || ROOTFile == "true" || ROOTFile == "TRUE"){
0146     gROOTOut = true;
0147     cout << "ROOT output file enabled." << endl;
0148   }
0149   else{
0150     gROOTOut = false;
0151     cout << "ROOT output file disabled." << endl;
0152   }
0153 
0154   TString Ejectile = obj["ejectile"].asString();
0155   TString RecoilHadron = obj["recoil_hadron"].asString(); // 09/02/22 - SJDK - Added in RecoilHadron type argument for K+
0156   // SJDK - 08/02/22 - This is terrible, need to change this, Ejectile should just be K+
0157   // Add a new flag which, RecoilHadron - where this is specified too, then add conditionals elsewhere based on this
0158   // New conditional, special case for Kaon 
0159   Ejectile = ExtractParticle(Ejectile);
0160   charge = ExtractCharge(Ejectile);
0161   if(RecoilHadron == "Sigma" || RecoilHadron == "sigma"){ // SJDK - 31/01/23 - If RecoilHadron specified as Sigma, interpret this as Sigma0. Also correct for lower case
0162     RecoilHadron = "Sigma0";
0163   }
0164   if (RecoilHadron == "lambda"){ // SJDK - 31/01/23 - Make Lambda selection case insensitive
0165     RecoilHadron = "Lambda"; 
0166   }
0167   if (Ejectile == "K+"){
0168     if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){
0169       RecoilHadron = "Lambda";
0170       cout << "! WARNING !" << endl;
0171       cout << "! WARNING !- K+ production specified but RecoilHadron not recognised, deaulting to Lambda - ! WARNING!" << endl;
0172       cout << "! WARNING !" << endl;
0173     }
0174     else{
0175       RecoilHadron = ExtractParticle(RecoilHadron);
0176     }
0177   }
0178   // SJDK - 19/12/22 - Specify RecoilHadron to neutron/proton for pi+/pi0 production, for pi0 production, may want to adjust? 
0179   else if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){
0180     RecoilHadron = "Neutron";
0181   }
0182   else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){
0183     RecoilHadron = "Proton";
0184   }
0185   else { // SJDK -09/02/22 - Note that in future this could be changed to get different RecoilHadrons in other reactions if desired
0186     RecoilHadron = "";
0187   }
0188 
0189   // SJDK 03/04/23 - Change to how Qsq range is set/chosen, could add as an override variable later too
0190   // Set min/max Qsq values depending upon Ejectile type
0191   if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){
0192     fQsq_Min = 3.5; fQsq_Max = 35.0; // Love Preet changed to 3.5 as SigT parameterization starts from Q2 = 3.5
0193     fW_Min = 2.0; fW_Max = 10.2;
0194     fT_Max_Default = 1.3;
0195   }
0196   else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){
0197     fQsq_Min = 5.0; fQsq_Max = 1000.0;
0198     fW_Min = 2.0; fW_Max = 10.0;
0199     fT_Max_Default = 0.5;
0200   }
0201   else if (Ejectile == "K+"){
0202     fQsq_Min = 1.0; fQsq_Max = 35.0;
0203     fW_Min = 2.0; fW_Max = 10.0;
0204     fT_Max_Default = 2.0;
0205   }
0206   else{
0207     fQsq_Min = 5.0; fQsq_Max = 35.0;
0208     fW_Min = 2.0; fW_Max = 10.0;
0209     fT_Max_Default = 1.3;
0210   }
0211 
0212   // SJDK - 01/06/21
0213   // Set beam energies from .json read in
0214   if (obj.isMember("ebeam")){
0215     fEBeam = obj["ebeam"].asDouble();
0216   }
0217   else{
0218     fEBeam = 5;
0219     cout << "Electron beam energy not specified in .json file, defaulting to 5 GeV." << endl;
0220   }
0221   if (obj.isMember("hbeam")){
0222     fPBeam = obj["hbeam"].asDouble();
0223     fHBeam = obj["hbeam"].asDouble();
0224   }
0225   else{
0226     fPBeam = 100;
0227     fHBeam = 100;
0228     cout << "Ion beam energy not specified in .json file, defaulting to 100 GeV." << endl;
0229   }
0230 
0231   if (obj.isMember("hbeam_part")){
0232     gBeamPart = obj["hbeam_part"].asString();
0233   } 
0234   else {
0235     gBeamPart = "Proton";
0236   }
0237 
0238   // SJDK - 12/01/22
0239   // Set output type as a .json read in
0240   // Should be Pythia6, LUND or HEPMC3
0241   if (obj.isMember("OutputType")){
0242     gOutputType = obj["OutputType"].asString();
0243     if (gOutputType == "Pythia6"){
0244       cout << "Using Pythia6 output format for Fun4All" << endl;
0245     }
0246     else if (gOutputType == "LUND"){
0247       cout << "Using LUND output format" << endl;
0248     }
0249     else if (gOutputType == "HEPMC3"){
0250       cout << "Using HEPMC3 output format for ePIC" << endl;
0251     }
0252     else{
0253       cout << "Output type not recognised!" << endl;
0254       cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl;
0255       gOutputType = "HEPMC3";
0256     }
0257   }
0258   else{
0259     cout << "Output type not specified in .json file!" << endl;
0260     cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl;
0261     gOutputType = "HEPMC3";
0262   }
0263   ///*--------------------------------------------------*/
0264   /// The detector selection is determined here
0265   /// The incidence proton phi angle is 
0266   if (obj.isMember("det_location")){
0267     gDet_location = obj["det_location"].asString();
0268     if (gDet_location == "ip8") {
0269       fProton_incidence_phi = 0.0;
0270     } 
0271     else if (gDet_location == "ip6") {
0272       fProton_incidence_phi = fPi;
0273     }
0274     else {
0275       fProton_incidence_phi = 0.0;
0276       cout << "The interaction point requested is not recognized!" << endl;
0277       cout << "Therefore ip6 is used by default." << endl;
0278     }
0279   }
0280   else{ // 21/12/22 - This could probably be combined with the else statement above in some way
0281     fProton_incidence_phi = 0.0;
0282     cout << "The interaction points was not specified in the .json file!" << endl;
0283     cout << "Therefore ip6 is used by default" << endl;
0284   }
0285 
0286   if (obj.isMember("Ee_Low")){
0287     fScatElec_E_Lo = obj["Ee_Low"].asDouble();
0288   }
0289   else{
0290     fScatElec_E_Lo = 0.5;
0291     cout << "Minumum scattered electron energy not specified in .json file, defaulting to 0.5*EBeam." << endl;
0292   }
0293 
0294   if (obj.isMember("Ee_High")){
0295     fScatElec_E_Hi = obj["Ee_High"].asDouble();
0296   }
0297   else{
0298     fScatElec_E_Hi = 2.5;
0299     cout << "Max scattered electron energy not specified in .json file, defaulting to 2.5*EBeam." << endl;
0300   }
0301 
0302   if (obj.isMember("e_Theta_Low")){
0303     fScatElec_Theta_I = obj["e_Theta_Low"].asDouble() * fDEG2RAD;
0304   }
0305   else{
0306     fScatElec_Theta_I = 60.0 * fDEG2RAD;
0307     cout << "Min scattered electron theta not specified in .json file, defaulting to 60 degrees." << endl;
0308   }
0309 
0310   if (obj.isMember("e_Theta_High")){
0311     fScatElec_Theta_F = obj["e_Theta_High"].asDouble() * fDEG2RAD;
0312   }
0313   else{
0314     fScatElec_Theta_F = 175.0 * fDEG2RAD;
0315     cout << "Max scattered electron theta not specified in .json file, defaulting to 175 degrees." << endl;
0316   }
0317 
0318   if (obj.isMember("EjectileX_Theta_Low")){
0319     fEjectileX_Theta_I = obj["EjectileX_Theta_Low"].asDouble() * fDEG2RAD;
0320   }
0321   else{
0322     fEjectileX_Theta_I = 0.0 * fDEG2RAD;
0323     cout << "Min ejectile X theta not specified in .json file, defaulting to 0 degrees." << endl;
0324   }
0325 
0326   if (obj.isMember("EjectileX_Theta_High")){
0327     fEjectileX_Theta_F = obj["EjectileX_Theta_High"].asDouble() * fDEG2RAD;
0328   }
0329   else{
0330     fEjectileX_Theta_F = 60.0 * fDEG2RAD;
0331     cout << "Max ejectile X theta not specified in .json file, defaulting to 60 degrees." << endl;
0332   }
0333 
0334   // 06/09/23 - SJDK - Added string to check method chosen
0335   TString CalcMethod = obj["calc_method"].asString(); 
0336   if(CalcMethod == "Analytical"){
0337     UseSolve = false;
0338   }
0339   else if (CalcMethod == "Solve"){
0340     UseSolve = true;
0341   }
0342   else{
0343     cout << "! WARNING !" << endl  << "! WARNING !- Calculation method not specified or not recognised, defaulting to Analytical - ! WARNING !" << endl << "! WARNING !" << endl;
0344     UseSolve = false;
0345   }
0346 
0347   // 17/07/25 - SJDK - Add new check to get t_max value
0348   if (obj.isMember("Kin_tMax")){
0349     if( obj["Kin_tMax"].asDouble() < 0 ){
0350       if( abs(obj["Kin_tMax"].asDouble()) > fT_Max_Default){
0351     cout << "! WARNING !" << endl << "! WARNING ! - Max -t entered as a negative number and absolute value exceeds limit for reaction, seting to reaction max! - ! WARNING !"<< endl << "! WARNING !"<< endl;
0352     fT_Max = fT_Max_Default;
0353       }
0354       else{
0355     cout << "! WARNING !" << endl << "! WARNING ! - Max -t entered as a negative number, taking absolute value! - ! WARNING !"<< endl << "! WARNING !"<< endl;
0356     fT_Max = abs(obj["Kin_tMax"].asDouble());
0357       }
0358     }
0359     else if (obj["Kin_tMax"].asDouble() == 0 || obj["Kin_tMax"].asDouble() > fT_Max_Default){
0360       cout << "! WARNING !" << endl << "! WARNING ! - Max -t set to 0 or exceeding paramaterisation limit for reaction, setting to reaction max! - ! WARNING !"<< endl << "! WARNING !"<< endl;
0361       fT_Max = fT_Max_Default;
0362     }
0363     else{
0364       fT_Max = obj["Kin_tMax"].asDouble();
0365     }
0366   }
0367   else{
0368     fT_Max = 1.0;
0369     cout << "! WARNING !" << endl << "! WARNING ! - Max -t not specified, defaulting to 1! - ! WARNING !"<< endl << "! WARNING !"<< endl;
0370   }
0371   SigPar = ReadCrossSectionPar(Ejectile, RecoilHadron);
0372     
0373   if(Ejectile != "pi0"){ // Default case now
0374     Reaction* r1 = new Reaction(Ejectile, RecoilHadron);
0375     r1->process_reaction();
0376     delete r1;
0377   }
0378   else{  // Treat pi0 slightly differently for now
0379     Reaction* r1 = new Reaction(Ejectile);
0380     r1->process_reaction();
0381     delete r1;
0382   }
0383 }
0384 
0385 /*--------------------------------------------------*/
0386 /*--------------------------------------------------*/
0387 
0388 void SetEICSeed (int seed) {
0389   fSeed = seed;
0390 }
0391 
0392 ///*--------------------------------------------------*/
0393 ///*--------------------------------------------------*/
0394 ///  Some utility functions
0395 
0396 ///*--------------------------------------------------*/
0397 /// Extracting the particle type
0398 
0399 TString ExtractParticle(TString particle) {
0400 
0401   /// Make the input particle case insansitive
0402   particle.ToLower();
0403   if (particle.Contains("on")) {
0404     particle.ReplaceAll("on", "");
0405   };
0406     
0407   if (particle.Contains("plus")) {
0408     particle.ReplaceAll("plus", "+");
0409   }
0410 
0411   if (particle.Contains("minus")) {
0412     particle.ReplaceAll("minus", "-");
0413   }
0414 
0415   if (particle.Contains("zero")) {
0416     particle.ReplaceAll("zero", "0");
0417   }
0418 
0419   particle[0] = toupper(particle[0]);
0420   cout << "Particle: " << particle << endl;
0421   return particle;
0422 
0423 }
0424 
0425 ///*--------------------------------------------------*/
0426 /// Extracting the particle charge
0427 
0428 TString ExtractCharge(TString particle) {
0429 
0430   TString charge;
0431 
0432   if (particle.Contains("+") || particle.Contains("plus")) {
0433     charge = "+";
0434   } else if (particle.Contains("-") || particle.Contains("minus")) {
0435     charge = "-";
0436   } else {
0437     charge = "0";
0438   }
0439   return charge;
0440 }
0441 
0442 vector<vector<vector<vector<double>>>> ReadCrossSectionPar(TString EjectileX, TString RecoilHad){
0443   
0444   string sigL_ParamFile, sigT_ParamFile;
0445  
0446   if (EjectileX == "Pi+" && RecoilHad == "Neutron"){
0447     // When pion model parameterised in some way, add Pi+/Neutron case here - 
0448   }
0449   else if (EjectileX == "Pi-" && RecoilHad == "Proton"){
0450     // When pion model parameterised in some way, add Pi-/Proton case here - 
0451   }
0452   else if (EjectileX == "K+" && RecoilHad == "Lambda"){
0453     sigL_ParamFile = Form("%s/src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigL", DEMPgen_Path);
0454     sigT_ParamFile = Form("%s/src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigT", DEMPgen_Path);
0455   }
0456   else if (EjectileX == "K+" && RecoilHad == "Sigma0"){
0457     sigL_ParamFile = Form("%s/src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigL", DEMPgen_Path);
0458     sigT_ParamFile = Form("%s/src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigT", DEMPgen_Path);
0459   }
0460   else if (EjectileX == "Pi0"){
0461     // When pi0 model parameterised, add it here
0462   }
0463   else{
0464     cerr << " !!!!! " << endl << "Warning!" << endl << "Combination of specified ejectile and recoil hadron not found!" << "Cross section parameters cannot be read, check inputs!" << endl << "Warning!" << endl << " !!!!! " << endl;
0465     exit(-1);
0466   }
0467  
0468   //....................................................................................................
0469   // Love's model parameters (Gojko, Stephen and Nishchey helped me to understand this part) 
0470   //....................................................................................................
0471   double ptmp;
0472   std::vector<std::vector<std::vector<std::vector<double>>>> p_vec;
0473   fstream file_vgl; // The parameterization file we will open and loop over
0474 
0475   for (int i = 0; i < 2; i++){
0476     if(i == 0){
0477       file_vgl.open(sigL_ParamFile, ios::in); 
0478     }
0479     if(i == 1){
0480       file_vgl.open(sigT_ParamFile, ios::in); 
0481     }
0482     p_vec.push_back(std::vector<std::vector<std::vector<double>>>());
0483     for(int j=0; j <9; j++){// Loop over all values of W - 2 to 10
0484       p_vec[i].push_back(std::vector<std::vector<double>>());
0485       for(int k=0; k<35; k++){ // Loop over all values of Q2 - 1 to 35 for each w
0486     p_vec[i][j].push_back(std::vector<double>());
0487       
0488     for(int l=0; l<13; l++){ //Loop over all columns at once
0489       file_vgl>>ptmp;
0490       p_vec[i][j][k].push_back(ptmp);
0491     }
0492       }
0493     }
0494     file_vgl.close();// Need to close the file at end of each loop over i 
0495   }       
0496   return p_vec;
0497 }
[ 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!