Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
eic_sphenix_sidis_countbins.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file eic_sphenix_sidis_countbins.C
1 int
2 eic_sphenix_sidis_countbins()
3 {
4  TFile *fin = new TFile("output/eic_sphenix_dis_histo_1M.root","OPEN");
5  THnSparse *hfull = (THnSparse*)fin->Get("hn_sidis_pion_plus_accept");
6  THnSparse *hfull_fullaccept = (THnSparse*)fin->Get("hn_sidis_pion_plus");
7 
8  TH2F* hxQ2 = (TH2F*)hfull->Projection(1,0);
9  TH2F* hxQ2_fullaccept = (TH2F*)hfull_fullaccept->Projection(1,0);
10 
11  /* beam energies */
12  float ebeam_e = 20;
13  float ebeam_p = 250;
14 
15  /* Pythia generated events, cross section, and luminosity */
16  float pythia_ngen = 1e6;
17  float pythia_xsec = 0.60785660255324614; // in microbarn
18  float convert_microbarn_to_femtobarn = 1e9;
19  float pythia_lumi = pythia_ngen / ( pythia_xsec * convert_microbarn_to_femtobarn );
20 
21  /* target luminosity and scaling factor */
22  float target_lumi = pythia_lumi; // in inverse femtobarn
23  float lumi_scaling = target_lumi / pythia_lumi;
24 
25  cout << "Pythia luminosity: " << pythia_lumi << " fb^-1" << endl;
26  cout << "Target luminosity: " << target_lumi << " fb^-1" << endl;
27  cout << "Luminosity scaling: " << lumi_scaling << endl;
28 
29  /* Minimum number of bin entries to accept bin */
30  int Nmin = 1;
31 
32  /* create tree to store information */
33  TTree *tcount = new TTree("tcount", "A tree with counts in kinematics bins");
34  float t_pbeam_lepton = 0;
35  float t_pbeam_proton = 0;
36  float t_s = 0;
37  float t_x = 0;
38  float t_Q2 = 0;
39  float t_y = 0;
40  float t_z = 0;
41  float t_pT = 0;
42  float t_N = 0;
43  float t_stdev_N = 0;
44  tcount->Branch("pbeam_lepton", &t_pbeam_lepton, "pbeam_lepton/F");
45  tcount->Branch("pbeam_proton", &t_pbeam_proton, "pbeam_proton/F");
46  tcount->Branch("s", &t_s, "s/F");
47  tcount->Branch("x", &t_x, "x/F");
48  tcount->Branch("Q2", &t_Q2, "Q2/F");
49  tcount->Branch("y", &t_y, "y/F");
50  tcount->Branch("z", &t_z, "z/F");
51  tcount->Branch("pT", &t_pT, "pT/F");
52  tcount->Branch("N", &t_N, "N/F");
53  tcount->Branch("stdev_N", &t_stdev_N, "stdev_N/F");
54 
55  /* copy beam parameters */
56  t_pbeam_lepton = ebeam_e;
57  t_pbeam_proton = ebeam_p;
58 
59  /* center of mass energy */
60  t_s = 4 * ebeam_e * ebeam_p;
61 
62 
63  /* Get all axis to loop over */
64  TAxis *ax_x = hfull->GetAxis(0);
65  TAxis *ax_Q2 = hfull->GetAxis(1);
66  TAxis *ax_z = hfull->GetAxis(2);
67  TAxis *ax_pT = hfull->GetAxis(3);
68 
69  ax_x->Print();
70  ax_Q2->Print();
71  ax_z->Print();
72  ax_pT->Print();
73 
74  /* collect all x-bins */
75  // set<float> s_binc_x;
76 
77  /* loop over all bins */
78  for ( int bin_x = 1; bin_x <= ax_x->GetNbins(); bin_x++ )
79  {
80  for ( int bin_Q2 = 1; bin_Q2 <= ax_Q2->GetNbins(); bin_Q2++ )
81  {
82  t_x = TMath::Power(10, 0.5 * ( ( TMath::Log10( ax_x->GetBinLowEdge(bin_x) ) )
83  + ( TMath::Log10( ax_x->GetBinLowEdge(bin_x) + ax_x->GetBinWidth(bin_x) ) ) ) );
84 
85  t_Q2 = TMath::Power(10, 0.5 * ( ( TMath::Log10( ax_Q2->GetBinLowEdge(bin_Q2) ) )
86  + ( TMath::Log10( ax_Q2->GetBinLowEdge(bin_Q2) + ax_Q2->GetBinWidth(bin_Q2) ) ) ) );
87 
88  /* Calculate inelasticity y */
89  t_y = t_Q2 / ( t_x * t_s );
90 
91  /* skip kinematics bins wth y > 0.95 and y < 1e-2 */
92  if ( t_y > 0.95 || t_y < 1e-2 )
93  continue;
94 
95  for ( int bin_z = 1; bin_z <= ax_z->GetNbins(); bin_z++ )
96  {
97 
98  t_z = ax_z->GetBinCenter( bin_z );
99 
100  for ( int bin_pT = 1; bin_pT <= ax_pT->GetNbins(); bin_pT++ )
101  {
102 
103  t_pT = ax_pT->GetBinCenter( bin_pT );
104 
105  /* Get bin entries */
106  int binloc[4] = { bin_x, bin_Q2, bin_z, bin_pT };
107  t_N = hfull->GetBinContent( binloc ) * lumi_scaling;
108 
109  /* skip bins with small number of entries */
110  if ( t_N < Nmin )
111  continue;
112 
113  t_stdev_N = 1./(sqrt(t_N));
114 
115  tcount->Fill();
116  //s_binc_x.insert(t_x);
117 
118  /* print values */
119  std::cout.precision(2);
120 
121  bool toscreen = false;
122  if ( toscreen )
123  {
124  cout // <<"lepton = " << std::fixed << t_pbeam_lepton
125  // << " x proton = " << std::fixed << t_pbeam_proton
126  // << " , sqrt(s) = " << std::fixed << sqrt( t_s )
127  << " , x = " << std::scientific << t_x
128  << " , Q2 = " << std::scientific << t_Q2
129  << " , y = " << std::fixed << t_y
130  << " , z = " << std::fixed << t_z
131  << " , pT = " << std::fixed << t_pT
132  << " , N = " << std::scientific << t_N
133  << endl;
134  }
135  }
136  }
137  }
138  }
139 
140  /* Prepare TPaveText for plots */
141  TString str_ebeam = TString::Format("%.0f GeV x %.0f GeV", ebeam_e, ebeam_p );
142  TString str_lumin = TString::Format("L = %.4f fb^{-1}", target_lumi );
143 
144  TPaveText *pt_ebeam_lumi_ul = new TPaveText(1e-5,1e3,1e-3,1e4,"none");
145  pt_ebeam_lumi_ul->SetFillStyle(0);
146  pt_ebeam_lumi_ul->SetLineColor(0);
147  pt_ebeam_lumi_ul->AddText(str_ebeam);
148  pt_ebeam_lumi_ul->AddText(str_lumin);
149 
150  TPaveText *pt_ebeam_lumi_ll = new TPaveText(1e2,45,1e3,50,"none");
151  pt_ebeam_lumi_ll->SetFillStyle(0);
152  pt_ebeam_lumi_ll->SetLineColor(0);
153  pt_ebeam_lumi_ll->AddText(str_ebeam);
154  pt_ebeam_lumi_ll->AddText(str_lumin);
155 
156  /* Prepare inelasticity (y) cut lines for plots */
157  TF1 *f_y095 = new TF1("f_y095", "4*x*[0]*[1]*[2]", 1e-5, 1);
158  f_y095->SetParameter( 0, ebeam_e);
159  f_y095->SetParameter( 1, ebeam_p);
160  f_y095->SetParameter( 2, 0.95);
161  TF1 *f_y001 = (TF1*)f_y095->Clone("f_y01");
162  f_y001->SetParameter(2 , 0.01);
163 
164  /* make x-Q2 plot */
165  TCanvas *c1 = new TCanvas();
166  c1->SetRightMargin(0.12);
167  hxQ2->Draw("colz");
168  c1->SetLogx(1);
169  c1->SetLogy(1);
170  c1->SetLogz(1);
171 
172  f_y095->Draw("same");
173  f_y001->Draw("same");
174 
175  pt_ebeam_lumi_ul->Draw();
176  gPad->RedrawAxis();
177 
178  /* make x-Q2 plot for 'perfect' acceptance */
179  TCanvas *c4 = new TCanvas();
180  c4->SetRightMargin(0.12);
181  hxQ2_fullaccept->Draw("colz");
182  c4->SetLogx(1);
183  c4->SetLogy(1);
184  c4->SetLogz(1);
185 
186  f_y095->Draw("same");
187  f_y001->Draw("same");
188 
189  pt_ebeam_lumi_ul->Draw();
190  gPad->RedrawAxis();
191 
192  /* make x-Q2 acceptance fraction pot */
193  TCanvas *c3 = new TCanvas();
194  c3->SetRightMargin(0.12);
195  TH2F* hxQ2_acceptance_ratio = hxQ2->Clone("x_Q2_acceptance_ratio");
196  hxQ2_acceptance_ratio->GetZaxis()->SetNdivisions(505);
197  hxQ2_acceptance_ratio->Divide(hxQ2_fullaccept);
198  hxQ2_acceptance_ratio->Draw("colz");
199  c3->SetLogx(1);
200  c3->SetLogy(1);
201 
202  f_y095->Draw("same");
203  f_y001->Draw("same");
204 
205  pt_ebeam_lumi_ul->Draw();
206  gPad->RedrawAxis();
207 
208  /* create tree to store information */
209  TFile *fout = new TFile("output/eic_sphenix_sidis_tree.root", "RECREATE");
210 
211  /* Write tree to output */
212  tcount->Write();
213  fout->Close();
214 
215  return 0;
216 }