From ae80f1cdb069aa87a72900e6780b013da0e4a36d Mon Sep 17 00:00:00 2001 From: jonzarling Date: Wed, 9 Dec 2020 16:52:45 -0500 Subject: [PATCH 1/4] Add non-default option to generate with uniform t-distribution, minor histogram updates --- .../Simulation/genEtaRegge/genEtaRegge.cc | 62 +++++++++++++++++-- 1 file changed, 58 insertions(+), 4 deletions(-) diff --git a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc index dde138678..5890d0b23 100644 --- a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc +++ b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc @@ -94,11 +94,19 @@ TH1D *thrown_dalitzZ; TH1D *thrown_Egamma; TH2D *thrown_dalitzXY; TH2D *thrown_theta_vs_p; +TH2D *thrown_theta_vs_p_eta; TH1D *cobrems_vs_E; char input_file_name[250]="eta548.in"; char output_file_name[250]="eta_gen.hddm"; +// Non-default option to generate uniform t-distribution from tmin to tmax +/// (calculating cross section at fixed t_uniform_eval) +bool gen_uniform_t=false; +double t_uniform_eval=-1; // Only used if gen_uniform is true. Currently, +float t_min_uniform=0; // takes min(t_min_uniform,t_0) so as to avoid unphysical t_values +float t_max_uniform=-100; // takes max(t_max_uniform,t_max) so as to avoid unphysical t_values + void Usage(void){ printf("genEtaRegge: generator for eta production based on Regge trajectory formalism.\n"); printf(" Usage: genEtaRegge \n"); @@ -106,6 +114,9 @@ void Usage(void){ printf(" -O (default: eta_gen.hddm)\n"); printf(" -I (default: eta548.in)\n"); printf(" -R (default: 10000)\n"); + printf(" -U (generate uniform t-distribution instead of sloped)\n"); + printf(" -tmin (min t value for uniform dist., if -U specified, default=physical min.)\n"); + printf(" -tmax (max t value for uniform dist., if -U specified, default=2.0)\n"); printf(" -h (Print this message and exit.)\n"); printf("Coupling constants, photon beam energy range, and eta decay products are\n"); printf("specified in the file.\n"); @@ -124,7 +135,25 @@ void ParseCommandLineArguments(int narg, char* argv[]) } for(int i=1; i=5) { + char *check_str = (char *)"-tmin", *matches=NULL; + matches=strstr(ptr,check_str); + if(matches) { + sscanf(&ptr[5],"%f",&t_min_uniform); + t_min_uniform=-abs(t_min_uniform); // Make sure value is negative, no matter what was supplied + } + + check_str = (char *)"-tmax"; matches=NULL; + matches=strstr(ptr,check_str); + if(matches) { + sscanf(&ptr[5],"%f",&t_max_uniform); + t_max_uniform=-abs(t_max_uniform); // Make sure value is negative, no matter what was supplied + } + } + + // For all other command line options (which are exactly one character long) if(ptr[0] == '-'){ switch(ptr[1]){ case 'h': Usage(); break; @@ -146,6 +175,9 @@ void ParseCommandLineArguments(int narg, char* argv[]) case 'd': debug=true; break; + case 'U': + gen_uniform_t=true; + break; default: break; } @@ -386,7 +418,7 @@ void WriteEvent(unsigned int eventNumber,TLorentzVector &beam, float vert[3], // Create some diagnostic histograms void CreateHistograms(string beamConfigFile){ - thrown_t=new TH1D("thrown_t","Thrown -t distribution",1000,0.,2.0); + thrown_t=new TH1D("thrown_t","Thrown -t distribution",1000,0.,3.0); thrown_t->SetXTitle("-t [GeV^{2}]"); thrown_dalitzZ=new TH1D("thrown_dalitzZ","thrown dalitz Z",110,-0.05,1.05); thrown_Egamma=new TH1D("thrown_Egamma","Thrown E_{#gamma} distribution", @@ -399,6 +431,11 @@ void CreateHistograms(string beamConfigFile){ thrown_theta_vs_p->SetXTitle("p [GeV/c]"); thrown_theta_vs_p->SetYTitle("#theta [degrees]"); + thrown_theta_vs_p_eta=new TH2D("thrown_theta_vs_p_eta","#eta #theta_{LAB} vs. p", + 120,0,12.,180,0.,180.); + thrown_theta_vs_p_eta->SetXTitle("p [GeV/c]"); + thrown_theta_vs_p_eta->SetYTitle("#theta [degrees]"); + BeamProperties beamProp(beamConfigFile); cobrems_vs_E = (TH1D*)beamProp.GetFlux(); } @@ -441,6 +478,8 @@ void GraphCrossSection(double &xsec_max){ t_old=t; } TGraph *Gxsec=new TGraph(10000,t_array,xsec_array); + TString xsec_title = "#eta Cross Section at E_{#gamma}="+to_string(Egamma)+" GeV;-t [GeV^{2}];d#sigma/dt [#mub/GeV^{2}]"; + Gxsec->SetTitle(xsec_title); Gxsec->Write("Cross section"); cout << "Total cross section at " << Egamma << " GeV = "<< sum @@ -455,6 +494,7 @@ int main(int narg, char *argv[]) { ParseCommandLineArguments(narg, argv); + // open ROOT file string rootfilename="eta_gen.root"; TFile *rootfile=new TFile(rootfilename.c_str(),"RECREATE", @@ -525,6 +565,8 @@ int main(int narg, char *argv[]) cout << "number of decay particles = " << num_decay_particles << endl; + if(gen_uniform_t) cout << "GENERATING DATA WITH UNIFORM T-DIST FROM " << t_min_uniform << " TO " << t_max_uniform << endl; + bool use_evtgen = false; #ifdef HAVE_EVTGEN // check to see if we should use EvtGen @@ -846,8 +888,19 @@ int main(int narg, char *argv[]) sin_theta_over_2=sin(0.5*theta_cm); t=t0-4.*p_gamma*p_eta*sin_theta_over_2*sin_theta_over_2; - xsec=CrossSection(s,t,p_gamma,p_eta,theta_cm); - + xsec=CrossSection(s,t,p_gamma,p_eta,theta_cm); + + // If generating a sample uniform in t, we need to fix t and re-calculate theta_cm based on it. Others do not depend on t. + if(gen_uniform_t) { + //Cross section at fixed t value (t_tmp) + double t_tmp = t_uniform_eval; + double theta_cm_tmp = 2.*asin(0.5*sqrt( (t0-t_tmp)/(p_gamma*p_eta) ) ); + xsec=CrossSection(s,t_tmp,p_gamma,p_eta,theta_cm_tmp); + // Make t uniform, calculate theta_cm based off of it + t=myrand->Uniform(t_max_uniform, min( float(t0),t_min_uniform) ); // If t_min_uniform provided is unphysical, then use physical t_min. + theta_cm=2.*asin(0.5*sqrt( (t0-t)/(p_gamma*p_eta) ) ); + } + // Generate a test value for the cross section xsec_test=myrand->Uniform(xsec_max); } @@ -875,6 +928,7 @@ int main(int narg, char *argv[]) //proton4.Print(); thrown_theta_vs_p->Fill(proton4.P(),180./M_PI*proton4.Theta()); + thrown_theta_vs_p_eta->Fill(eta4.P(),180./M_PI*eta4.Theta()); // Other diagnostic histograms thrown_t->Fill(-t); From 01b6344cf5ed60eb4b929aa1846103bea12e5757 Mon Sep 17 00:00:00 2001 From: jonzarling Date: Wed, 9 Dec 2020 17:38:53 -0500 Subject: [PATCH 2/4] Add u-channel support, check that tmax is physical --- src/programs/Simulation/genEtaRegge/genEtaRegge.cc | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc index 5890d0b23..f595e3f72 100644 --- a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc +++ b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc @@ -105,7 +105,7 @@ char output_file_name[250]="eta_gen.hddm"; bool gen_uniform_t=false; double t_uniform_eval=-1; // Only used if gen_uniform is true. Currently, float t_min_uniform=0; // takes min(t_min_uniform,t_0) so as to avoid unphysical t_values -float t_max_uniform=-100; // takes max(t_max_uniform,t_max) so as to avoid unphysical t_values +float t_max_uniform=-3; // takes max(t_max_uniform,t_max) so as to avoid unphysical t_values void Usage(void){ printf("genEtaRegge: generator for eta production based on Regge trajectory formalism.\n"); @@ -418,7 +418,7 @@ void WriteEvent(unsigned int eventNumber,TLorentzVector &beam, float vert[3], // Create some diagnostic histograms void CreateHistograms(string beamConfigFile){ - thrown_t=new TH1D("thrown_t","Thrown -t distribution",1000,0.,3.0); + thrown_t=new TH1D("thrown_t","Thrown -t distribution",1000,0.,abs(t_max_uniform)); thrown_t->SetXTitle("-t [GeV^{2}]"); thrown_dalitzZ=new TH1D("thrown_dalitzZ","thrown dalitz Z",110,-0.05,1.05); thrown_Egamma=new TH1D("thrown_Egamma","Thrown E_{#gamma} distribution", @@ -565,6 +565,11 @@ int main(int narg, char *argv[]) cout << "number of decay particles = " << num_decay_particles << endl; + if( abs(t_max_uniform)>21. ) { // Max t at GlueX endpoint energy is about 20 GeV^2. Reset value to protect against inefficient accept/reject. + t_max_uniform=-21; + cout << "tmax provided is larger than physically allowed t at GlueX highest E, resetting to physical max........" << endl; + } + if(gen_uniform_t) cout << "GENERATING DATA WITH UNIFORM T-DIST FROM " << t_min_uniform << " TO " << t_max_uniform << endl; bool use_evtgen = false; @@ -899,6 +904,7 @@ int main(int narg, char *argv[]) // Make t uniform, calculate theta_cm based off of it t=myrand->Uniform(t_max_uniform, min( float(t0),t_min_uniform) ); // If t_min_uniform provided is unphysical, then use physical t_min. theta_cm=2.*asin(0.5*sqrt( (t0-t)/(p_gamma*p_eta) ) ); + if( std::isnan(theta_cm)==true ) xsec=-1.; // Lazy person's way of skipping unphysical theta_cm. Breaking do/while to accept event will never be satisfied for this case. } // Generate a test value for the cross section From cd05020ecaf25910f5c8c561d6e63a1616ac5682 Mon Sep 17 00:00:00 2001 From: jonzarling Date: Thu, 10 Dec 2020 11:24:33 -0500 Subject: [PATCH 3/4] Add u-channel support, check that tmax is physical --- src/programs/Simulation/genEtaRegge/genEtaRegge.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc index f595e3f72..293d15bfb 100644 --- a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc +++ b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc @@ -116,7 +116,7 @@ void Usage(void){ printf(" -R (default: 10000)\n"); printf(" -U (generate uniform t-distribution instead of sloped)\n"); printf(" -tmin (min t value for uniform dist., if -U specified, default=physical min.)\n"); - printf(" -tmax (max t value for uniform dist., if -U specified, default=2.0)\n"); + printf(" -tmax (max t value for uniform dist., if -U specified, default=3.0)\n"); printf(" -h (Print this message and exit.)\n"); printf("Coupling constants, photon beam energy range, and eta decay products are\n"); printf("specified in the file.\n"); From 0c727201951a228c2a36d0f6275c698980edcfdf Mon Sep 17 00:00:00 2001 From: jonzarling Date: Thu, 10 Dec 2020 11:27:50 -0500 Subject: [PATCH 4/4] Make usage text agree with default tmax --- src/programs/Simulation/genEtaRegge/genEtaRegge.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc index 293d15bfb..ab5c807ad 100644 --- a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc +++ b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc @@ -124,6 +124,7 @@ void Usage(void){ exit(0); } + //----------- // ParseCommandLineArguments //-----------