src/eic_evgen/eic.cc

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