DEMPgen/src/main.cxx

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