Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-17 08:07:02 

0001 ///*--------------------------------------------------*/
0002 /// main.cxx:
0003 ///*--------------------------------------------------*/
0004 /// Old Generator name: DEMP_EvGen
0005 /// Original author: Rory Evens
0006 /// Date: 2015-2018
0007 ///
0008 ///*--------------------------------------------------*/
0009 /// New Generator name: DEMPgen
0010 /// Modifier: Wenliang (Bill) Li
0011 /// Date: Feb 24 2020
0012 /// Email: wenliang.billlee@gmail.com
0013 ///
0014 /// Comment: March 04, 2020: 
0015 
0016 #pragma clang diagnostic ignored "-Wdeprecated-declarations"
0017 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0018 
0019 //C++lib includes
0020 #include <stdio.h>
0021 #include <stdlib.h>
0022 #include <iostream>
0023 #include <fstream>
0024 #include <iomanip>
0025 #include "json/json.h"
0026 #include "json/json-forwards.h"
0027 
0028 //Root includes
0029 #include "TROOT.h"
0030 #include "TRandom.h"
0031 #include "TFile.h"
0032 #include "TH1.h"
0033 #include "TCanvas.h"
0034 #include "TApplication.h"
0035 #include "TMath.h"
0036 #include "TVector3.h"
0037 #include "TSystem.h"
0038 
0039 //Project includes
0040 #include "Particle.hxx"
0041 #include "CustomRand.hxx"
0042 #include "ScatteredParticleGen.hxx"
0043 #include "ProductGen.hxx"
0044 #include "TreeBuilder.hxx"
0045 #include "Constants.hxx"
0046 #include "DEMPEvent.hxx"
0047 #include "SigmaCalc.hxx"
0048 #include "TargetGen.hxx"
0049 #include "MatterEffects.hxx"
0050 #include "FSI.hxx"
0051 
0052 #include "eic_evgen/eic.h"
0053 
0054 using namespace std;
0055 using namespace constants;
0056 
0057 TFile * AsymmFile;
0058 char* DEMPgen_Path;
0059 
0060 Json::Value obj; //Declared here for global access
0061 
0062 string get_date(void);
0063 
0064 int Gen_seed;
0065 
0066 int main(int argc, char** argv){
0067 
0068   cout << endl << "DEMPgen Copyright (C) 2024 Z. Ahmed, R.S. Evans. I. Goel. G.M. Huber, S.J.D. Kay, W.B. Li, L. Preet, A. Usman." << endl;
0069   cout << "This program comes with ABSOLUTELY NO WARRANTY." << endl; 
0070   cout << "This is free software, and you are welcome to redistribute it under certain conditions; contact authors for details." << endl;
0071   cout << "See - https://github.com/JeffersonLab/DEMPgen - for author contact details." << endl << endl;
0072 
0073   // Grab DEMPgen_path
0074   // Check if the environment variable was found
0075   DEMPgen_Path = getenv("DEMPgen");
0076   if (DEMPgen_Path == nullptr) {
0077     cerr << "!!!!! ERROR !!!!! DEMPgen environment variable not set !!!!! ERROR !!!!!" << endl;
0078     cerr << "!!!!! ERROR !!!!! Source the setup.sh (or .csh) script and rerun !!!!! ERROR !!!!!" << endl;
0079     exit(0);
0080   }
0081   
0082   // Parsing of config file.
0083   //ifstream ifs("../Config.json");
0084   ifstream ifs(argv[1]);
0085   Json::Reader reader;
0086   reader.parse(ifs, obj);
0087  
0088   unsigned long long int nEvents = obj["n_events"].asUInt64();
0089   cout << "Generating "<< nEvents << " events."<<endl;
0090 
0091   TString file_name = obj["file_name"].asString();
0092 
0093   if (file_name == "") {
0094     file_name.Form("%i", nEvents);  
0095     file_name = file_name + "_" + get_date();
0096   } 
0097   
0098   int gen_seed = obj["generator_seed"].asInt();
0099   TString particle = obj["particle"].asString();
0100 
0101   if (obj["experiment"].asString() == "eic") {
0102  
0103     cout << "EIC is used " << endl;
0104  
0105     int target_direction = obj["Targ_dir"].asInt();
0106     int kinematics_type = obj["Kinematics_type"].asInt();
0107     double EBeam = obj["ebeam"].asDouble();
0108     double HBeam = obj["hbeam"].asDouble();
0109     TString hadron = obj["hadron"].asString(); // SJDK 08/02/22 - Add the hadron type as an argument
0110     //      bool = obj["pi0_particle"].asBool()
0111     eic(obj);
0112  
0113   }
0114 
0115   else if (obj["experiment"].asString() == "solid") {
0116  
0117     gROOT->ProcessLine( "gErrorIgnoreLevel = 3001;");
0118 
0119     Gen_seed = gen_seed;
0120 
0121     cout << "SoLID is used " << endl;
0122     cout << "ROOT based error printouts supressed, comment line 117 of src/main.cxx and recompile to re-enable" << endl;
0123     
0124     if (obj["ionisation"].asBool())
0125       cout << "Ionisation Enabled" << endl;
0126     if (obj["bremsstrahlung"].asBool())
0127       cout << "Bremsstrahlung Enabled" << endl;
0128     if (obj["fermi_momentum"].asBool())
0129       cout << "Fermi Momentum Enabled" << endl;
0130     if (obj["multiple_scattering"].asBool())
0131       cout << "Multiple Scattering Enabled" << endl;
0132     if (obj["final_state_interaction"].asBool())
0133       cout << "FSI Enabled" << endl;
0134    
0135     MatterEffects* ME = new MatterEffects();
0136 
0137     // Check an asymmetries file exists to load in for later use. Throw error if not.
0138     if(gSystem->AccessPathName(Form("%s/data/input/Asymmetries.root", DEMPgen_Path)) == kTRUE){
0139       cerr << Form("!!!!! ERROR !!!!! %s/data/input/Asymmetries.root not found  !!!!! ERROR !!!!!", DEMPgen_Path) << endl;
0140       cerr << "!!!!! ERROR !!!!! Check path! !!!!! ERROR !!!!!" << endl;
0141       exit(1);
0142     }
0143     else{
0144       // Load in the file containing asymmetry information
0145       AsymmFile = new TFile(Form("%s/data/input/Asymmetries.root", DEMPgen_Path));
0146     }
0147     
0148     // Initilization of DEMPEvent objects for different reference frames
0149     DEMPEvent* VertEvent = new DEMPEvent("Vert");
0150     // VertEvent contains particles with kinematic properties at the vertex.
0151     // VertEvent is viewed in the laboratory reference frame.
0152     // This event is not to be modified after initial generation of the event.
0153     DEMPEvent* RestEvent = new DEMPEvent("RF");
0154     // Event in rest frame of the target neutron.
0155     DEMPEvent* CofMEvent = new DEMPEvent("CofM");
0156     // Event in the rest frame of the system's center of mass.
0157     DEMPEvent* TConEvent = new DEMPEvent("TCon"); 
0158     // Event rotated to follow the Trento Conventions
0159     DEMPEvent* LCorEvent = new DEMPEvent("Lab");
0160     // LCorEvent is the event after all corrections and effects have been applied,
0161     // in the laboratory reference frame.
0162     
0163     SigmaCalc* Sig = new SigmaCalc(VertEvent, CofMEvent, RestEvent, TConEvent);
0164     
0165     // Retrieval of pointers to particles in VertEvent.
0166     // For clarity: these are the same objects as are referenced by VertEvent,
0167     // not copies. Operations on these objects affect the original. 
0168     Particle* VertBeamElec = VertEvent->BeamElec;
0169     Particle* VertTargNeut = VertEvent->TargNeut;
0170     Particle* VertScatElec = VertEvent->ScatElec;
0171     Particle* VertProdPion = VertEvent->ProdPion;
0172     Particle* VertProdProt = VertEvent->ProdProt;
0173     
0174     Particle* Photon = VertEvent->VirtPhot;
0175     Photon->SetName("VirtPhot");
0176 
0177     Particle* FSIProt = new Particle(proton_mass_mev, "FSIProt", pid_prot);
0178    
0179     Particle* InTotal = new Particle();
0180     Particle* OutTotal = new Particle();
0181     
0182     VertBeamElec->SetThetaPhiE(0, 0, obj["beam_energy"].asDouble());
0183    
0184     double elecERange[2] = {obj["scat_elec_Emin"].asDouble(),
0185       obj["scat_elec_Emax"].asDouble()};
0186     double elecThetaRange[2] = {obj["scat_elec_thetamin"].asDouble()/DEG,
0187       obj["scat_elec_thetamax"].asDouble()/DEG};
0188     double elecPhiRange[2] = {0, 360/DEG};
0189    
0190     TargetGen * NeutGen = new TargetGen(neutron_mass_mev, obj["fermi_momentum"].asBool());
0191     
0192     ScatteredParticleGen * ElecGen =
0193       new ScatteredParticleGen(electron_mass_mev,
0194                    elecERange,
0195                    elecThetaRange,
0196                    elecPhiRange);
0197 
0198 
0199     FSI* FSIobj = new FSI();
0200 
0201     ProductGen * ProtonPionGen = new ProductGen(Photon,
0202                         VertTargNeut);
0203     
0204     int nSuccess = 0;
0205     int nFail = 0;
0206     int nNeg = 0;
0207     int nCut = 0;
0208    
0209     int FSIfail = 0;
0210    
0211     int event_status = 0;
0212    
0213     file_name = Form("%s/data/output/Solid_DEMP_%s.root", DEMPgen_Path, file_name.Data());
0214 
0215     TreeBuilder * Output = new TreeBuilder(file_name.Data(), "t1");
0216 
0217     Output->AddEvent(VertEvent);
0218     //Output->AddEvent(CofMEvent);
0219     //Output->AddEvent(RestEvent);
0220     Output->AddEvent(LCorEvent);
0221    
0222     //if (obj["final_state_interaction"].asBool())
0223     Output->AddParticle(FSIProt);
0224    
0225     Output->AddParticle(Photon);
0226 
0227     // These parameters are calculated using multiple reference frames (DEMPEvent objects),
0228     // and need to be added to the output seperately.
0229     double sigma_l;
0230     double sigma_t;
0231     double sigma_tt;
0232     double sigma_lt;
0233    
0234     double sigma_uu;
0235     double sigma_ut;
0236    
0237     double sigma_k0;
0238     double sigma_k1;
0239     double sigma_k2;
0240     double sigma_k3;
0241     double sigma_k4;
0242    
0243     double sigma_k5 = 0;
0244    
0245     double sigma;
0246    
0247     double weight;
0248     double epsilon;
0249    
0250     double targetthickness, airthickness, targwindowthickness;
0251    
0252     targwindowthickness = 0.0120*Window_Density / Window_Thickness;
0253    
0254     Output -> AddDouble(&sigma_l,"sigma_l");
0255     Output -> AddDouble(&sigma_t,"sigma_t");
0256     Output -> AddDouble(&sigma_tt,"sigma_tt");
0257     Output -> AddDouble(&sigma_lt,"sigma_lt");
0258     Output -> AddDouble(&sigma_uu,"sigma_uu");
0259     Output -> AddDouble(&sigma_ut,"sigma_ut");
0260    
0261     Output -> AddDouble(&sigma_k0,"AutPhiMinusPhiS");
0262     Output -> AddDouble(&sigma_k1,"AutPhiS");
0263     Output -> AddDouble(&sigma_k2,"Aut2PhiMinusPhiS");
0264     Output -> AddDouble(&sigma_k3,"AutPhiPlusPhiS");
0265     Output -> AddDouble(&sigma_k4,"Aut3PhiMinusPhiS");
0266     Output -> AddDouble(&sigma_k5,"Aut2PhiPlusPhiS");
0267    
0268     Output -> AddDouble(&sigma, "sigma");
0269    
0270     Output -> AddDouble(&weight,"weight");
0271     Output -> AddDouble(&epsilon, "epsilon");
0272    
0273     //if(obj["final_state_interaction"].asBool()){
0274     Output -> AddDouble(FSIobj->PhaseShiftWeight, "PhaseShiftWeight");
0275     Output -> AddDouble(FSIobj->WilliamsWeight, "WilliamsWeight");
0276     Output -> AddDouble(FSIobj->DedrickWeight, "DedrickWeight");
0277     Output -> AddDouble(FSIobj->CatchenWeight, "CatchenWeight");
0278     //}
0279    
0280     // Event global variables to add to output
0281     Output -> AddDouble(VertEvent->qsq_GeV, "Vert_Qsq_GeV");
0282     Output -> AddDouble(VertEvent->w_GeV, "Vert_w_GeV");
0283     Output -> AddDouble(VertEvent->t_GeV, "Vert_t_GeV");
0284     Output -> AddDouble(VertEvent->t_para_GeV, "Vert_t_para_GeV");
0285     Output -> AddDouble(VertEvent->t_prime_GeV, "Vert_t_prime_GeV");
0286     Output -> AddDouble(VertEvent->negt, "Vert_negt_GeV");
0287     Output -> AddDouble(VertEvent->x_B, "Vert_x_B");
0288     Output -> AddDouble(VertEvent->Phi_deg, "Vert_Phi");
0289     Output -> AddDouble(VertEvent->Phi_s_deg, "Vert_Phi_s");
0290     Output -> AddDouble(VertEvent->Theta_deg, "Vert_Theta");
0291    
0292     Output -> AddDouble(LCorEvent->qsq_GeV, "Lab_Qsq_GeV");
0293     Output -> AddDouble(LCorEvent->w_GeV, "Lab_w_GeV");
0294     Output -> AddDouble(LCorEvent->t_GeV, "Lab_t_GeV");
0295     Output -> AddDouble(LCorEvent->t_para_GeV, "Lab_t_para_GeV");
0296     Output -> AddDouble(LCorEvent->t_prime_GeV, "Lab_t_prime_GeV");
0297     Output -> AddDouble(LCorEvent->negt, "Lab_negt_GeV");
0298     Output -> AddDouble(LCorEvent->x_B, "Lab_x_B");
0299     Output -> AddDouble(LCorEvent->Phi_deg, "Lab_Phi");
0300     Output -> AddDouble(LCorEvent->Phi_s_deg, "Lab_Phi_s");
0301     Output -> AddDouble(LCorEvent->Theta_deg, "Lab_Theta");
0302    
0303     Output -> AddDouble(VertEvent->Vertex_x, "Vertex_x");
0304     Output -> AddDouble(VertEvent->Vertex_y, "Vertex_y");
0305     Output -> AddDouble(VertEvent->Vertex_z, "Vertex_z");
0306 
0307     // Main loop of the generator (SoLID module)
0308     for (int i=0; i<nEvents; i++){
0309    
0310       if (i%100 == 0)
0311     cout <<"\r"<< (double)i/nEvents * 100 << "%";
0312    
0313       // Generate vertex location
0314       *VertEvent->Vertex_x = gRandom->Uniform(-0.25, 0.25);
0315       *VertEvent->Vertex_y = gRandom->Uniform(-0.25,0.25);
0316       *VertEvent->Vertex_z = gRandom->Uniform(-370,-330);
0317    
0318       // Reset inc. electron to beam energy
0319    
0320       VertBeamElec->SetThetaPhiE(0, 0, obj["beam_energy"].asDouble());
0321    
0322       // Perform matter effects on inc electron
0323       // These effects occur before the reaction, so affect the vertex values
0324       targetthickness = ((*VertEvent->Vertex_z+370.0) * Helium_Density)/
0325     (ME->X0(Helium_Z, Helium_A));
0326       if (obj["ionisation"].asBool()){
0327     ME->IonLoss(VertBeamElec, Window_A, Window_Z, Window_Density, targwindowthickness);
0328     ME->IonLoss(VertBeamElec, Helium_A, Helium_Z, Helium_Density, targetthickness);
0329       }
0330       if (obj["bremsstrahlung"].asBool()){
0331     ME->BremLoss(VertBeamElec,
0332              targetthickness*ME->b(Helium_Z)
0333              /ME->X0(Helium_Z, Helium_A));
0334     ME->BremLoss(VertBeamElec,
0335              targwindowthickness*ME->b(Window_Z)
0336              /ME->X0(Window_Z, Window_A));
0337       }
0338       if (obj["multiple_scattering"].asBool()){
0339     ME->MultiScatter(VertBeamElec,
0340              targetthickness / ME->X0(Helium_Z, Helium_A));
0341     ME->MultiScatter(VertBeamElec,
0342              targwindowthickness / ME->X0(Window_Z, Window_A));
0343       }
0344    
0345       // Generate target and scattered electron
0346       *VertTargNeut = *NeutGen->GetParticle();
0347       *VertScatElec = *ElecGen->GetParticle();
0348       *Photon = *VertBeamElec - *VertScatElec;
0349       
0350       // Solve for remaining particles
0351       event_status = ProtonPionGen->Solve();
0352       if (event_status == 0)
0353     nSuccess ++;
0354       if (event_status == 1){
0355     nFail ++;
0356     continue;
0357       }
0358       *VertProdPion = *ProtonPionGen->ProdPion();
0359       *VertProdProt = *ProtonPionGen->ProdProton();
0360       
0361       VertEvent->Update();
0362       // VertEvent and its components are not to be modified beyond this point.
0363    
0364       // Copy VertEvent and boost to CofM frame.
0365       *CofMEvent = *VertEvent;
0366       CofMEvent->Boost(-VertEvent->CoM());
0367       CofMEvent->Update();
0368    
0369       // Copy VertEvent and boost to rest frame
0370       *RestEvent = *VertEvent;
0371       RestEvent->Boost(-(VertEvent->TargNeut->Vect()*(1/VertEvent->TargNeut->E())));
0372       RestEvent->Update();
0373    
0374       // Copy RestEvent and rotate to Trento coords.
0375       *TConEvent = *RestEvent;
0376       TConEvent->Rotate(RestEvent->VirtPhot->Theta(),-RestEvent->ScatElec->Phi());
0377       TConEvent->Update();
0378    
0379       // Copy VertEvent onto LCorEvent, revert beamelec back to beam value
0380       *LCorEvent = *VertEvent;
0381       LCorEvent->BeamElec->SetThetaPhiE(0, 0, obj["beam_energy"].asDouble());
0382       LCorEvent->Update();
0383    
0384       // Get cross-sections
0385       sigma_l = Sig->sigma_l();
0386       sigma_t = Sig->sigma_t();
0387       sigma_lt = Sig->sigma_lt();
0388       sigma_tt = Sig->sigma_tt();
0389       sigma_uu = Sig->sigma_uu();
0390       sigma_ut = Sig->sigma_ut();
0391    
0392       sigma_k0 = Sig->Sigma_k(0);
0393       sigma_k1 = Sig->Sigma_k(1);
0394       sigma_k2 = Sig->Sigma_k(2);
0395       sigma_k3 = Sig->Sigma_k(3);
0396       sigma_k4 = Sig->Sigma_k(4);
0397    
0398       sigma = Sig->sigma();
0399 
0400       epsilon = Sig->epsilon();
0401 
0402       // Cuts
0403    
0404       if (obj["Qsq_cut"].asBool()){
0405     if (*VertEvent->qsq_GeV<obj["Qsq_min"].asDouble()){
0406       nSuccess --;
0407       nCut++;
0408       continue;
0409     }
0410       }
0411    
0412       if (obj["w_cut"].asBool()){
0413     if (*VertEvent->w_GeV<obj["w_min"].asDouble()){
0414       nSuccess --;
0415       nCut++;
0416       continue;
0417     }
0418       }
0419    
0420       if (obj["t_cut"].asBool()){
0421     if (*VertEvent->t_GeV<obj["t_min"].asDouble()){
0422       nSuccess --;
0423       nCut++;
0424       continue;
0425     }
0426       }
0427    
0428       if (sigma<=0){
0429     //cout << "NEGATIVE EVENT" << endl;
0430     nNeg++;
0431       }
0432    
0433       if (obj["weight_cut"].asBool()){
0434     if (sigma<=0) {
0435       nSuccess --;
0436       // cout << "NEGATIVE EVENT" << endl;
0437       // cout << "Sigma_l: \t" << sigma_l << endl;
0438       // cout << "Sigma_t: \t" << sigma_t << endl;
0439       // cout << "Sigma_lt: \t" << sigma_lt << endl;
0440       // cout << "Sigma_tt: \t" << sigma_tt << endl;
0441       // cout << "Sigma_uu: \t" << sigma_uu << endl;
0442       // cout << "Sigma_ut: \t" << sigma_ut << endl;
0443       // cout << "Sigma: \t" << sigma << endl;
0444       continue;
0445     }
0446       }
0447       
0448       weight = Sig->weight(nEvents);
0449    
0450       // Final State Interaction.
0451        
0452       if (obj["final_state_interaction"].asBool()){
0453     *FSIobj->VertInPion = *VertProdPion;
0454     FSIfail = FSIobj->Generate();
0455     if (FSIfail == 1){
0456       cout << "FSI Generation Failure!" << endl;
0457       // No instances of this so far.
0458       continue;
0459     }
0460     *FSIProt = *FSIobj->VertOutProt;
0461     *LCorEvent->ProdPion = *FSIobj->VertOutPion;
0462     //if (FSIfail == 0) cout << "FSI Generation Successful" << endl;
0463     FSIobj -> CalculateWeights();
0464       }
0465    
0466       // Final Cons Law Check
0467    
0468       *InTotal = (*(VertEvent->BeamElec)+
0469           *(VertEvent->TargNeut)
0470           );
0471       *OutTotal = (*(VertEvent->ScatElec)+
0472            *(LCorEvent->ProdPion)+
0473            *(VertEvent->ProdProt)
0474            );
0475       if (obj["final_state_interaction"].asBool()){
0476     *InTotal += *(FSIobj->VertTargProt);
0477     *OutTotal += *(FSIobj->VertOutProt);
0478       }
0479    
0480       if ((*InTotal-*OutTotal).Px() > 1.0)
0481     cout << "Px Violation" << endl;
0482       if ((*InTotal-*OutTotal).Py() > 1.0)
0483     cout << "Px Violation" << endl;
0484       if ((*InTotal-*OutTotal).Pz() > 1.0)
0485     cout << "Px Violation" << endl;
0486       if ((*InTotal-*OutTotal).E() > 1.0)
0487     cout << "Px Violation" << endl;
0488    
0489       //Matter Effects~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0490       
0491       targetthickness = ((-330.0 - *VertEvent->Vertex_z) * Helium_Density)/
0492     (ME->X0(Helium_Z, Helium_A));
0493       // Assuming path along z only inside the target.
0494       // Shape of target is not really being considered here.
0495       // A similar approximation is made for through-air distances.
0496       // Exact geometry of the detector is not considered, only large angle
0497       // or small angle.
0498       // If more detailed analysis is desired, geant can do it better.
0499       // Turn off effects and use GEMC.
0500    
0501       //cout << targetthickness << endl;
0502    
0503       if (obj["ionisation"].asBool()){
0504    
0505     //Scattered Electron
0506     ME->IonLoss(LCorEvent->ScatElec, Window_A, Window_Z,
0507             Window_Density, targwindowthickness);
0508     ME->IonLoss(LCorEvent->ScatElec, Helium_A, Helium_Z,
0509             Helium_Density, targetthickness);
0510     airthickness = 450 * Air_Density;
0511     if (LCorEvent->ScatElec->Theta() < 16) // Large Angle
0512       airthickness = 950 * Air_Density;
0513     ME->IonLoss(LCorEvent->ScatElec, Air_A, Air_Z,
0514             Air_Density, airthickness);
0515    
0516     //Proton
0517     airthickness = 450 * Air_Density;
0518     if (LCorEvent->ProdProt->Theta() < 16) // Large Angle
0519       airthickness = 950 * Air_Density;
0520    
0521     ME->IonLoss(LCorEvent->ProdProt, Helium_A, Helium_Z,
0522             Helium_Density, targetthickness);
0523    
0524     ME->IonLoss(LCorEvent->ProdProt, Window_A, Window_Z,
0525             Window_Density, targwindowthickness);
0526    
0527     ME->IonLoss(LCorEvent->ProdProt, Air_A, Air_Z,
0528             Air_Density, airthickness);
0529    
0530     //Pion
0531     airthickness = 450 * Air_Density;
0532     if (LCorEvent->ProdPion->Theta() < 16) // Large Angle
0533       airthickness = 950 * Air_Density;
0534    
0535     ME->IonLoss(LCorEvent->ProdPion, Helium_A, Helium_Z,
0536             Helium_Density, targetthickness);
0537    
0538     ME->IonLoss(LCorEvent->ProdPion, Window_A, Window_Z,
0539             Window_Density, targwindowthickness);
0540    
0541     ME->IonLoss(LCorEvent->ProdPion, Air_A, Air_Z,
0542             Air_Density, airthickness);
0543    
0544       }
0545    
0546       if (obj["bremsstrahlung"].asBool()){
0547    
0548     //Scattered Electron
0549     ME->BremLoss(LCorEvent->ScatElec,
0550              targetthickness*ME->b(Helium_Z)
0551              /ME->X0(Helium_Z, Helium_A));
0552     ME->BremLoss(LCorEvent->ScatElec,
0553              targwindowthickness*ME->b(Window_Z)
0554              /ME->X0(Window_Z, Window_A));
0555     airthickness = 450 * Air_Density / ME->X0(Air_Z, Air_A);
0556     if (LCorEvent->ScatElec->Theta() < 16) // Large Angle
0557       airthickness = 950 * Air_Density / ME->X0(Air_Z, Air_A);
0558     ME->BremLoss(LCorEvent->ScatElec, airthickness*ME->b(Air_Z));
0559    
0560       }
0561    
0562       if (obj["multiple_scattering"].asBool()){
0563     //Scattered Electron
0564     ME->MultiScatter(LCorEvent->ScatElec,
0565              targetthickness / ME->X0(Helium_Z, Helium_A));
0566     ME->MultiScatter(LCorEvent->ScatElec,
0567              targwindowthickness / ME->X0(Window_Z, Window_A));
0568     airthickness = 450 * Air_Density / ME->X0(Air_Z,Air_A);
0569     if (LCorEvent->ScatElec->Theta() < 16) // Large Angle
0570       airthickness = 950 * Air_Density / ME->X0(Air_Z,Air_A);
0571     ME->MultiScatter(LCorEvent->ScatElec, airthickness);
0572    
0573     //Proton
0574     ME->MultiScatter(LCorEvent->ProdProt,
0575              targetthickness / ME->X0(Helium_Z, Helium_A));
0576     ME->MultiScatter(LCorEvent->ProdProt,
0577              targwindowthickness / ME->X0(Window_Z, Window_A));
0578     airthickness = 450 * Air_Density / ME->X0(Air_Z,Air_A);
0579     if (LCorEvent->ProdProt->Theta() < 16) // Large Angle
0580       airthickness = 950 * Air_Density / ME->X0(Air_Z,Air_A);
0581     ME->MultiScatter(LCorEvent->ProdProt, airthickness);
0582    
0583     //Pion
0584     ME->MultiScatter(LCorEvent->ProdPion,
0585              targetthickness / ME->X0(Helium_Z, Helium_A));
0586     ME->MultiScatter(LCorEvent->ProdPion,
0587              targwindowthickness / ME->X0(Window_Z, Window_A));
0588     airthickness = 450 * Air_Density / ME->X0(Air_Z,Air_A);
0589     if (LCorEvent->ProdPion->Theta() < 16) // Large Angle
0590       airthickness = 950 * Air_Density / ME->X0(Air_Z,Air_A);
0591     ME->MultiScatter(LCorEvent->ProdPion, airthickness);
0592    
0593       }
0594       
0595       Output->Fill();
0596    
0597     }
0598    
0599     cout << endl;
0600    
0601     if (nSuccess>0)
0602       Output->Save();
0603    
0604     cout << "Successful Events: \t" << nSuccess << endl;
0605     cout << "Failed Events: \t\t" << nFail << endl;
0606     cout << "Negative Events: \t\t" << nNeg << endl;
0607     cout << "Cut Events: \t\t" << nCut << endl;
0608 
0609     // Checks against old event generator:
0610    
0611     if(nEvents <0){
0612    
0613       int tests = -nEvents;
0614       int N = 10000;
0615    
0616       TFile * Check = new TFile("RootFiles/DEMP_Ee_11_Events_10000_File_0.root");
0617       TTree * t1 = (TTree*)Check->Get("t1");
0618    
0619       cout << "Running Debug/Check Values" << endl;
0620    
0621       Double_t Epsilon, Qsq_GeV, T_GeV, W_GeV, x, y, z, WSq_GeV, EventWeight, PhaseShiftWeight, PhaseSpaceWeight;
0622       Double_t Qsq_Corrected_GeV, T_Corrected_GeV, W_Corrected_GeV;
0623    
0624       double WilliamsWeight, DedrickWeight, CatchenWeight;
0625    
0626       Double_t ScatElec_Energy_Col_GeV,ScatElec_MomX_Col_GeV,ScatElec_MomY_Col_GeV,ScatElec_MomZ_Col_GeV,ScatElec_Mom_Col_GeV;
0627       Double_t ScatElec_Phi_Col,ScatElec_Theta_Col;
0628    
0629       Double_t ScatElec_Corrected_Energy_Col_GeV,ScatElec_Corrected_MomX_Col_GeV,ScatElec_Corrected_MomY_Col_GeV,ScatElec_Corrected_MomZ_Col_GeV,ScatElec_Corrected_Mom_Col_GeV;
0630       Double_t ScatElec_Corrected_Phi_Col,ScatElec_Corrected_Theta_Col;
0631    
0632       Double_t Pion_FSI_Energy_Col_GeV,Pion_FSI_MomX_Col_GeV,Pion_FSI_MomY_Col_GeV,Pion_FSI_MomZ_Col_GeV, Pion_FSI_Mom_Col_GeV;
0633       Double_t Pion_FSI_Phi_Col, Pion_FSI_Theta_Col;
0634    
0635       Double_t Pion_Energy_Col_GeV,Pion_MomX_Col_GeV,Pion_MomY_Col_GeV,Pion_MomZ_Col_GeV, Pion_Mom_Col_GeV;
0636       Double_t Pion_Phi_Col, Pion_Theta_Col;
0637    
0638       Double_t Pion_Corrected_Energy_Col_GeV,Pion_Corrected_MomX_Col_GeV,Pion_Corrected_MomY_Col_GeV,Pion_Corrected_MomZ_Col_GeV, Pion_Corrected_Mom_Col_GeV;
0639       Double_t Pion_Corrected_Phi_Col, Pion_Corrected_Theta_Col;
0640    
0641       Double_t RecoilProton_Energy_Col_GeV, RecoilProton_MomX_Col_GeV, RecoilProton_MomY_Col_GeV, RecoilProton_MomZ_Col_GeV, RecoilProton_Mom_Col_GeV;
0642       Double_t RecoilProton_Phi_Col, RecoilProton_Theta_Col;
0643    
0644       Double_t RecoilProton_Corrected_Energy_Col_GeV, RecoilProton_Corrected_MomX_Col_GeV, RecoilProton_Corrected_MomY_Col_GeV, RecoilProton_Corrected_MomZ_Col_GeV;
0645       Double_t RecoilProton_Corrected_Phi_Col, RecoilProton_Corrected_Theta_Col, RecoilProton_Corrected_Mom_Col_GeV;
0646    
0647       Double_t AsymPhiMinusPhi_S, AsymPhi_S, Asym2PhiMinusPhi_S, AsymPhiPlusPhi_S, Asym3PhiMinusPhi_S, Asym2PhiPlusPhi_S;
0648    
0649       double T_Para_GeV, T_Prime_GeV;
0650    
0651       t1->SetBranchAddress("Epsilon", &Epsilon );
0652       t1->SetBranchAddress("Qsq_GeV", &Qsq_GeV );
0653       t1->SetBranchAddress("T_GeV", &T_GeV );
0654       t1->SetBranchAddress("W_GeV", &W_GeV );
0655       t1->SetBranchAddress("T_Para_GeV", &T_Para_GeV);
0656       //t1->SetBranchAddress("T_Prime_GeV", &T_Prime_GeV);
0657    
0658       t1->SetBranchAddress("Qsq_Corrected_GeV", &Qsq_Corrected_GeV );
0659       t1->SetBranchAddress("T_Corrected_GeV", &T_Corrected_GeV );
0660       t1->SetBranchAddress("W_Corrected_GeV", &W_Corrected_GeV );
0661    
0662       t1->SetBranchAddress("ScatElec_Energy_Col_GeV", &ScatElec_Energy_Col_GeV );
0663       t1->SetBranchAddress("ScatElec_MomX_Col_GeV", &ScatElec_MomX_Col_GeV );
0664       t1->SetBranchAddress("ScatElec_MomY_Col_GeV", &ScatElec_MomY_Col_GeV );
0665       t1->SetBranchAddress("ScatElec_MomZ_Col_GeV", &ScatElec_MomZ_Col_GeV );
0666       t1->SetBranchAddress("ScatElec_Mom_Col_GeV", &ScatElec_Mom_Col_GeV );
0667       t1->SetBranchAddress("ScatElec_Theta_Col", &ScatElec_Theta_Col );
0668       t1->SetBranchAddress("ScatElec_Phi_Col", &ScatElec_Phi_Col );
0669    
0670       t1->SetBranchAddress("ScatElec_Corrected_Energy_Col_GeV", &ScatElec_Corrected_Energy_Col_GeV );
0671       t1->SetBranchAddress("ScatElec_Corrected_MomX_Col_GeV", &ScatElec_Corrected_MomX_Col_GeV );
0672       t1->SetBranchAddress("ScatElec_Corrected_MomY_Col_GeV", &ScatElec_Corrected_MomY_Col_GeV );
0673       t1->SetBranchAddress("ScatElec_Corrected_MomZ_Col_GeV", &ScatElec_Corrected_MomZ_Col_GeV );
0674       t1->SetBranchAddress("ScatElec_Corrected_Mom_Col_GeV", &ScatElec_Corrected_Mom_Col_GeV );
0675       t1->SetBranchAddress("ScatElec_Corrected_Theta_Col", &ScatElec_Corrected_Theta_Col );
0676       t1->SetBranchAddress("ScatElec_Corrected_Phi_Col", &ScatElec_Corrected_Phi_Col );
0677    
0678       // t1->SetBranchAddress("Pion_FSI_Energy_Col_GeV", &Pion_FSI_Energy_Col_GeV );
0679       // t1->SetBranchAddress("Pion_FSI_MomX_Col_GeV", &Pion_FSI_MomX_Col_GeV );
0680       // t1->SetBranchAddress("Pion_FSI_MomY_Col_GeV", &Pion_FSI_MomY_Col_GeV );
0681       // t1->SetBranchAddress("Pion_FSI_MomZ_Col_GeV", &Pion_FSI_MomZ_Col_GeV );
0682       // t1->SetBranchAddress("Pion_FSI_Mom_Col_GeV", &Pion_FSI_Mom_Col_GeV );
0683       // t1->SetBranchAddress("Pion_FSI_Theta_Col", &Pion_FSI_Theta_Col );
0684       // t1->SetBranchAddress("Pion_FSI_Phi_Col", &Pion_FSI_Phi_Col );
0685    
0686       t1->SetBranchAddress("Pion_Energy_Col_GeV", &Pion_Energy_Col_GeV );
0687       t1->SetBranchAddress("Pion_MomX_Col_GeV", &Pion_MomX_Col_GeV );
0688       t1->SetBranchAddress("Pion_MomY_Col_GeV", &Pion_MomY_Col_GeV );
0689       t1->SetBranchAddress("Pion_MomZ_Col_GeV", &Pion_MomZ_Col_GeV );
0690       t1->SetBranchAddress("Pion_Mom_Col_GeV", &Pion_Mom_Col_GeV );
0691       t1->SetBranchAddress("Pion_Theta_Col", &Pion_Theta_Col );
0692       t1->SetBranchAddress("Pion_Phi_Col", &Pion_Phi_Col );
0693    
0694       t1->SetBranchAddress("Pion_Corrected_Energy_Col_GeV", &Pion_Corrected_Energy_Col_GeV );
0695       t1->SetBranchAddress("Pion_Corrected_MomX_Col_GeV", &Pion_Corrected_MomX_Col_GeV );
0696       t1->SetBranchAddress("Pion_Corrected_MomY_Col_GeV", &Pion_Corrected_MomY_Col_GeV );
0697       t1->SetBranchAddress("Pion_Corrected_MomZ_Col_GeV", &Pion_Corrected_MomZ_Col_GeV );
0698       t1->SetBranchAddress("Pion_Corrected_Mom_Col_GeV", &Pion_Corrected_Mom_Col_GeV );
0699       t1->SetBranchAddress("Pion_Corrected_Theta_Col", &Pion_Corrected_Theta_Col );
0700       t1->SetBranchAddress("Pion_Corrected_Phi_Col", &Pion_Corrected_Phi_Col );
0701    
0702       t1->SetBranchAddress("RecoilProton_Energy_Col_GeV", &RecoilProton_Energy_Col_GeV );
0703       t1->SetBranchAddress("RecoilProton_MomX_Col_GeV", &RecoilProton_MomX_Col_GeV );
0704       t1->SetBranchAddress("RecoilProton_MomY_Col_GeV", &RecoilProton_MomY_Col_GeV );
0705       t1->SetBranchAddress("RecoilProton_MomZ_Col_GeV", &RecoilProton_MomZ_Col_GeV );
0706       t1->SetBranchAddress("RecoilProton_Mom_Col_GeV", &RecoilProton_Mom_Col_GeV );
0707       t1->SetBranchAddress("RecoilProton_Theta_Col", &RecoilProton_Theta_Col );
0708       t1->SetBranchAddress("RecoilProton_Phi_Col", &RecoilProton_Phi_Col );
0709    
0710       t1->SetBranchAddress("RecoilProton_Corrected_Energy_Col_GeV", &RecoilProton_Corrected_Energy_Col_GeV );
0711       t1->SetBranchAddress("RecoilProton_Corrected_MomX_Col_GeV", &RecoilProton_Corrected_MomX_Col_GeV );
0712       t1->SetBranchAddress("RecoilProton_Corrected_MomY_Col_GeV", &RecoilProton_Corrected_MomY_Col_GeV );
0713       t1->SetBranchAddress("RecoilProton_Corrected_MomZ_Col_GeV", &RecoilProton_Corrected_MomZ_Col_GeV );
0714       t1->SetBranchAddress("RecoilProton_Corrected_Mom_Col_GeV", &RecoilProton_Corrected_Mom_Col_GeV );
0715       t1->SetBranchAddress("RecoilProton_Corrected_Theta_Col", &RecoilProton_Corrected_Theta_Col );
0716       t1->SetBranchAddress("RecoilProton_Corrected_Phi_Col", &RecoilProton_Corrected_Phi_Col );
0717    
0718       t1->SetBranchAddress("EventWeight", &EventWeight );
0719       t1->SetBranchAddress("WilliamsWeight", &WilliamsWeight );
0720       t1->SetBranchAddress("DedrickWeight", &DedrickWeight );
0721       t1->SetBranchAddress("CatchenWeight", &CatchenWeight );
0722    
0723       double ZASig_T, ZASig_L, ZASig_LT, ZASig_TT, ZASigma_Lab, ZASigma_UU, RorySigma_UT;
0724       double Theta, Phi, Phi_S;
0725    
0726       t1->SetBranchAddress("ZASig_T", &ZASig_T);
0727       t1->SetBranchAddress("ZASig_L", &ZASig_L);
0728       t1->SetBranchAddress("ZASig_LT", &ZASig_LT);
0729       t1->SetBranchAddress("ZASig_TT", &ZASig_TT);
0730       t1->SetBranchAddress("ZASigma_Lab", &ZASigma_Lab);
0731       t1->SetBranchAddress("ZASigma_UU", &ZASigma_UU);
0732       t1->SetBranchAddress("RorySigma_UT", &RorySigma_UT);
0733    
0734       t1->SetBranchAddress("AsymPhiMinusPhi_S", &AsymPhiMinusPhi_S);
0735       t1->SetBranchAddress("AsymPhi_S", &AsymPhi_S);
0736       t1->SetBranchAddress("Asym2PhiMinusPhi_S", &Asym2PhiMinusPhi_S);
0737       t1->SetBranchAddress("AsymPhiPlusPhi_S", &AsymPhiPlusPhi_S);
0738       t1->SetBranchAddress("Asym3PhiMinusPhi_S", &Asym3PhiMinusPhi_S);
0739       t1->SetBranchAddress("Asym2PhiPlusPhi_S", &Asym2PhiPlusPhi_S);
0740    
0741       t1->SetBranchAddress("Phi", &Phi);
0742       t1->SetBranchAddress("PhiS", &Phi_S);
0743       t1->SetBranchAddress("Photon_Theta_Col", &Theta);
0744    
0745       double Jacobian_CM, Jacobian_CM_RF, Jacobian_CM_Col;
0746       double Flux_Factor_RF, Flux_Factor_Col, A;
0747    
0748       t1->SetBranchAddress("Jacobian_CM",&Jacobian_CM);
0749       t1->SetBranchAddress("Jacobian_CM_RF",&Jacobian_CM_RF);
0750       t1->SetBranchAddress("Jacobian_CM_Col",&Jacobian_CM_Col);
0751       t1->SetBranchAddress("A", &A);
0752    
0753       t1->SetBranchAddress("Flux_Factor_RF",&Flux_Factor_RF);
0754       t1->SetBranchAddress("Flux_Factor_Col",&Flux_Factor_Col);
0755    
0756       bool ALERT = false;
0757       for (int i=0; i<tests; i++){
0758    
0759     t1->GetEntry(i);
0760    
0761     VertScatElec->SetThetaPhiE(ScatElec_Theta_Col/DEG, ScatElec_Phi_Col/DEG,
0762                    ScatElec_Energy_Col_GeV * 1000);
0763     *VertTargNeut = *NeutGen->GetParticle();
0764    
0765     /// Setting Photon
0766     *Photon = *VertBeamElec - *VertScatElec;
0767    
0768     ProtonPionGen->Solve(Pion_Theta_Col/DEG,Pion_Phi_Col/DEG);   
0769 
0770     /// Setting Pion
0771     *VertProdPion = *ProtonPionGen->ProdPion();
0772     /// Setting Proton
0773     *VertProdProt = *ProtonPionGen->ProdProton();
0774    
0775     VertEvent->Update();
0776    
0777     *CofMEvent = *VertEvent;
0778     CofMEvent->Boost(-VertEvent->CoM());
0779     CofMEvent->Update();
0780    
0781     *RestEvent = *VertEvent;
0782     RestEvent->Boost(-(VertEvent->TargNeut->Vect()*(1/VertEvent->TargNeut->E())));
0783     RestEvent->Update();
0784    
0785     *TConEvent = *RestEvent;
0786     TConEvent->Rotate(RestEvent->VirtPhot->Theta(),-RestEvent->ScatElec->Phi());
0787     TConEvent->Update();
0788    
0789     double qsq_GeV = *VertEvent->qsq_GeV;
0790     double w_GeV = *VertEvent->w_GeV;
0791     double t_GeV = *VertEvent->t_GeV;
0792     double t_prime_GeV = *VertEvent->t_prime_GeV;
0793     double t_para_GeV = *VertEvent->t_para_GeV;
0794     double phi = *TConEvent->Phi;
0795     if (phi<0) phi+=2*TMath::Pi();
0796     double phi_s = *TConEvent->Phi_s;
0797     if (phi_s<0) phi_s+=2*TMath::Pi();
0798     double theta = *RestEvent->Theta;
0799     if (theta <0) theta+=2*TMath::Pi();
0800       
0801     sigma_l = Sig->sigma_l();
0802     sigma_t = Sig->sigma_t();
0803     sigma_lt = Sig->sigma_lt();
0804     sigma_tt = Sig->sigma_tt();
0805     sigma_uu = Sig->sigma_uu();
0806     sigma_ut = Sig->sigma_ut();
0807    
0808     sigma_k0 = Sig->Sigma_k(0);
0809     sigma_k1 = Sig->Sigma_k(1);
0810     sigma_k2 = Sig->Sigma_k(2);
0811     sigma_k3 = Sig->Sigma_k(3);
0812     sigma_k4 = Sig->Sigma_k(4);
0813    
0814     sigma = Sig->sigma();
0815    
0816     if (obj["final_state_interaction"].asBool()){
0817       *FSIobj->VertInPion = *VertProdPion;
0818       FSIobj->VertOutPion -> SetPxPyPzE(Pion_Corrected_MomX_Col_GeV*1000,
0819                         Pion_Corrected_MomY_Col_GeV*1000,
0820                         Pion_Corrected_MomZ_Col_GeV*1000,
0821                         Pion_Corrected_Energy_Col_GeV*1000);
0822       //FSIfail = FSIobj->Generate();
0823       //if (FSIfail == 1){
0824       //  cout << "FSI Generation Failure!" << endl;
0825       //  continue;
0826       //}
0827       FSIobj->GenerateNoRandom();
0828       *FSIProt = *FSIobj->VertOutProt;
0829       *LCorEvent->ProdPion = *FSIobj->VertOutPion;
0830       FSIobj -> CalculateWeights();
0831     }
0832    
0833     int printw = 20;
0834    
0835     ALERT = false;
0836     if(TMath::Abs((sigma-ZASigma_Lab)/sigma)>0.01) {ALERT = true; cout << "SIGMA:\t"<<TMath::Abs((sigma-ZASigma_Lab)/sigma)<<endl;}
0837     if(TMath::Abs((sigma_l-ZASig_L)/sigma_l)>0.01) {ALERT = true; cout << "SIGMA_L:\t"<<TMath::Abs((sigma_l-ZASig_L)/sigma_l)<<endl;}
0838     if(TMath::Abs((sigma_t-ZASig_T)/sigma_t)>0.01) {ALERT = true; cout << "SIGMA_T:\t"<<TMath::Abs((sigma_t-ZASig_T)/sigma_t)<<endl;}
0839     if(TMath::Abs((sigma_lt-ZASig_LT)/sigma_lt)>0.01) {ALERT = true; cout << "SIGMA_LT:\t"<<TMath::Abs((sigma_lt-ZASig_LT)/sigma_lt)<<endl;}
0840     if(TMath::Abs((sigma_tt-ZASig_TT)/sigma_tt)>0.01) {ALERT = true; cout << "SIGMA_TT:\t"<<TMath::Abs((sigma_tt-ZASig_TT)/sigma_tt)<<endl;}
0841     if(TMath::Abs((sigma_uu-ZASigma_UU)/sigma_uu)>0.01) {ALERT = true; cout << "SIGMA_UU:\t"<<TMath::Abs((sigma_uu-ZASigma_UU)/sigma_uu)<<endl;}
0842     if(TMath::Abs((sigma_ut-RorySigma_UT)/sigma_ut)>0.01) {ALERT = true; cout << "SIGMA_UT:\t"<<TMath::Abs((sigma_ut-RorySigma_UT)/sigma_ut)<<endl;}
0843     if(TMath::Abs((Sig->weight(N)-EventWeight)/EventWeight)>0.01) {ALERT = true; cout << "WEIGHT:\t"<<TMath::Abs((Sig->weight(N)-EventWeight)/EventWeight)<<endl;}
0844    
0845     if(ALERT){
0846    
0847       cout<<left<<setw(printw)<<"Event:"<<left<<setw(printw)<<i<<endl;
0848       cout<<left<<setw(printw)<<"Kinematics:"<<endl;
0849    
0850       cout<<left<<setw(printw)<< "ElecE:"<<left<<setw(printw) << VertScatElec->E()/1000<<left<<setw(printw)<<ScatElec_Energy_Col_GeV << endl;
0851       cout<<left<<setw(printw)<< "ElecTh:"<<left<<setw(printw)<<VertScatElec->Theta()*DEG<<left<<setw(printw)<<ScatElec_Theta_Col << endl;
0852       cout<<left<<setw(printw)<< "ElecPhi:"<<left<<setw(printw)<<VertScatElec->Phi()*DEG <<left<<setw(printw)<<ScatElec_Phi_Col << endl;
0853    
0854       cout<<left<<setw(printw)<< "PionE:"<<left<<setw(printw)<<VertProdPion->E()/1000<<left<<setw(printw)<<Pion_Energy_Col_GeV << endl;
0855       cout<<left<<setw(printw)<< "PionTh:"<<left<<setw(printw)<<VertProdPion->Theta()*DEG<<left<<setw(printw)<<Pion_Theta_Col<< endl;
0856       cout<<left<<setw(printw)<< "PionPhi:"<<left<<setw(printw)<<VertProdPion->Phi()*DEG<<left<<setw(printw)<<Pion_Phi_Col<< endl;
0857    
0858       cout<<left<<setw(printw)<< "ProtE:"<<left<<setw(printw)<<VertProdProt->E()/1000<<left<<setw(printw)<<RecoilProton_Energy_Col_GeV << endl;
0859       cout<<left<<setw(printw)<< "ProtTh:"<<left<<setw(printw)<<VertProdProt->Theta()*DEG<<left<<setw(printw)<<RecoilProton_Theta_Col<< endl;
0860       cout<<left<<setw(printw)<< "ProtPhi:"<<left<<setw(printw)<<VertProdProt->Phi()*DEG<<left<<setw(printw)<<RecoilProton_Phi_Col << endl;
0861    
0862       cout<<left<<setw(printw)<< "Qsq:"<<left<<setw(printw) << qsq_GeV<<left<<setw(printw)<<Qsq_GeV << endl;
0863       cout<<left<<setw(printw)<< "W:"<<left<<setw(printw) << w_GeV<<left<<setw(printw)<<W_GeV << endl;
0864       cout<<left<<setw(printw)<< "t:"<<left<<setw(printw) << t_GeV<<left<<setw(printw)<<T_GeV << endl;
0865       cout<<left<<setw(printw)<< "t_para:"<<left<<setw(printw) << t_para_GeV<<left<<setw(printw)<<T_Para_GeV << endl;
0866       cout<<left<<setw(printw)<< "t':"<<left<<setw(printw) << t_prime_GeV<<left<<setw(printw)<<T_Prime_GeV << endl;
0867       cout<<left<<setw(printw)<<"Theta:"<<left<<setw(printw)<<theta*DEG<<left<<setw(printw)<<Theta<<endl;
0868       cout<<left<<setw(printw)<<"Phi"<<left<<setw(printw)<<phi*DEG<<left<<setw(printw)<<Phi<<endl;
0869       cout<<left<<setw(printw)<<"PhiS"<<left<<setw(printw)<<phi_s*DEG<<left<<setw(printw)<<Phi_S<<endl;
0870       cout<<left<<setw(printw)<< endl;
0871    
0872       cout<<left<<setw(printw)<<"Cross Sections: "<<endl;
0873       cout<<left<<setw(printw)<<"sigma:"<<left<<setw(printw)<< sigma<<left<<setw(printw)<<ZASigma_Lab << endl;
0874       cout<<left<<setw(printw)<<"sigma_l:"<<left<<setw(printw)<< sigma_l<<left<<setw(printw)<<ZASig_L << endl;
0875       cout<<left<<setw(printw)<<"sigma_t:"<<left<<setw(printw)<< sigma_t<<left<<setw(printw)<<ZASig_T << endl;
0876       cout<<left<<setw(printw)<<"sigma_lt:"<<left<<setw(printw)<< sigma_lt<<left<<setw(printw)<<ZASig_LT << endl;
0877       cout<<left<<setw(printw)<<"sigma_tt:"<<left<<setw(printw)<< sigma_tt<<left<<setw(printw)<<ZASig_TT << endl;
0878       cout<<left<<setw(printw)<<"sigma_uu:"<<left<<setw(printw)<< sigma_uu<<left<<setw(printw)<<ZASigma_UU << endl;
0879       cout<<left<<setw(printw)<<"sigma_ut:"<<left<<setw(printw)<< sigma_ut<<left<<setw(printw)<<RorySigma_UT << endl;
0880       cout<<left<<setw(printw)<<"AsyP-Ps:"<<left<<setw(printw)<< sigma_k0<<left<<setw(printw)<<AsymPhiMinusPhi_S << endl;
0881       cout<<left<<setw(printw)<<"AsyPs:"<<left<<setw(printw)<< sigma_k1<<left<<setw(printw)<<AsymPhi_S << endl;
0882       cout<<left<<setw(printw)<<"Asy2P-Ps:"<<left<<setw(printw)<< sigma_k2<<left<<setw(printw)<<Asym2PhiMinusPhi_S << endl;
0883       cout<<left<<setw(printw)<<"AsyP+Ps:"<<left<<setw(printw)<< sigma_k3<<left<<setw(printw)<<AsymPhiPlusPhi_S << endl;
0884       cout<<left<<setw(printw)<<"Asy3P-Ps:"<<left<<setw(printw)<< sigma_k4<<left<<setw(printw)<<Asym3PhiMinusPhi_S << endl;
0885       cout<<left<<setw(printw)<<"Asy2P_Ps:"<<left<<setw(printw)<< sigma_k5<<left<<setw(printw)<<Asym2PhiPlusPhi_S << endl;
0886       cout<<left<<setw(printw)<<"Epsilon:"<<left<<setw(printw)<< Sig->epsilon()<<left<<setw(printw)<<Epsilon << endl;
0887       cout<<left<<setw(printw)<<"Jacobian CM:"<<left<<setw(printw)<<Sig->jacobian_cm()<<left<<setw(printw)<<Jacobian_CM<<endl;
0888       cout<<left<<setw(printw)<<"Jacobian CM-Col:"<<left<<setw(printw)<<Sig->jacobian_cm_col()<<left<<setw(printw)<<Jacobian_CM_Col<<endl;
0889       cout<<left<<setw(printw)<<"Flux Factor Col:"<<left<<setw(printw)<<Sig->fluxfactor_col()<<left<<setw(printw)<<Flux_Factor_Col<<endl;
0890       cout<<left<<setw(printw)<<"A:"<<left<<setw(printw)<<Sig->jacobian_A()<<left<<setw(printw)<<A<<endl;
0891       cout<<left<<setw(printw)<<"Event Weight:"<<left<<setw(printw)<<Sig->weight(N)<<left<<setw(printw)<<EventWeight<<left<<setw(printw)<< endl;
0892       //cout<<left<<setw(printw)<<EventWeight/Sig->weight(100000) << endl;
0893       cout<<left<<setw(printw)<< endl;
0894    
0895       cout << (ZASigma_UU+RorySigma_UT)*Flux_Factor_Col*A*Jacobian_CM_Col << endl;
0896    
0897       cout<<left<<setw(printw)<< endl;
0898    
0899       // Sig->Asyms->at(0)->PrintPars();
0900    
0901       cout<<left<<setw(printw)<< endl;
0902       cout<<left<<setw(printw)<<"WilliamWeight:"<<left<<setw(printw)<<*FSIobj->WilliamsWeight<<left<<setw(printw)<<WilliamsWeight<<left<<setw(printw)<<endl;
0903       cout<<left<<setw(printw)<<"DedrickWeight:"<<left<<setw(printw)<<*FSIobj->DedrickWeight<<left<<setw(printw)<<DedrickWeight<<left<<setw(printw)<<endl;
0904       cout<<left<<setw(printw)<<"CatchenWeight:"<<left<<setw(printw)<<*FSIobj->CatchenWeight<<left<<setw(printw)<<CatchenWeight<<left<<setw(printw)<<endl;
0905    
0906    
0907       cout << "--------------------------------------------------------------------------------------------"<<endl;
0908     }
0909       }
0910     }
0911  
0912   }
0913 
0914   else {
0915  
0916     cerr << endl;
0917     cerr << "/*--------------------------------------------------*/" << endl;
0918     cerr << "/*--------------------------------------------------*/" << endl;
0919     cerr << "/*                      Error!                      */" << endl;
0920     cerr << "/*            Not a valid experiment!               */" << endl;
0921     cerr << "/*            Option use EIC or solid!              */" << endl;
0922     cerr << "/*--------------------------------------------------*/" << endl;
0923     cerr << "/*--------------------------------------------------*/" << endl;
0924     cerr << endl;
0925  
0926     exit(2);
0927  
0928   }
0929  
0930   return 0;
0931 }
0932 
0933 string get_date(void) {
0934  
0935   int MAX_DATE = 22;
0936 
0937   time_t now;
0938   char the_date[MAX_DATE];
0939 
0940   the_date[0] = '\0';
0941 
0942   now = time(NULL);
0943 
0944   if (now != -1)
0945     {
0946       strftime(the_date, MAX_DATE, "%Y_%h_%d_%H_%M_%S", localtime(&now));
0947     }
0948 
0949   return string(the_date);
0950 }
[ 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!