Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hera_dis_countbins.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hera_dis_countbins.C
1 int
2 hera_dis_countbins()
3 {
4  TFile *fin = new TFile("hera_dis_histo.root","OPEN");
5  THnSparse *hfull = (THnSparse*)fin->Get("hn_dis");
6 
7  TH2F* hxQ2 = (TH2F*)hfull->Projection(1,0);
8 
9  /* beam energies */
10  float ebeam_e = 27.5;
11  float ebeam_p = 820;
12 
13  /* Pythia generated events, cross section, and luminosity */
14  float pythia_ngen = 999999;
15  float pythia_xsec = 0.80120995103272474; // in microbarn
16  float convert_microbarn_to_femtobarn = 1e9;
17  float pythia_lumi = pythia_ngen / ( pythia_xsec * convert_microbarn_to_femtobarn );
18 
19  /* target luminosity and scaling factor */
20  float target_lumi = 0.5; // in inverse femtobarn
21  float lumi_scaling = target_lumi / pythia_lumi;
22 
23  cout << "Pythia luminosity: " << pythia_lumi << " fb^-1" << endl;
24  cout << "Target luminosity: " << target_lumi << " fb^-1" << endl;
25  cout << "Luminosity scaling: " << lumi_scaling << endl;
26 
27  /* create tree to store information */
28  TFile *fout = new TFile("hera_dis_tree.root", "RECREATE");
29 
30  TTree *tcount = new TTree("tcount", "A tree with counts in kinematics bins");
31  float t_pbeam_lepton = 0;
32  float t_pbeam_proton = 0;
33  float t_s = 0;
34  float t_x = 0;
35  float t_Q2 = 0;
36  float t_y = 0;
37  float t_N = 0;
38  float t_stdev_N = 0;
39  tcount->Branch("pbeam_lepton", &t_pbeam_lepton, "pbeam_lepton/F");
40  tcount->Branch("pbeam_proton", &t_pbeam_proton, "pbeam_proton/F");
41  tcount->Branch("s", &t_s, "s/F");
42  tcount->Branch("x", &t_x, "x/F");
43  tcount->Branch("Q2", &t_Q2, "Q2/F");
44  tcount->Branch("y", &t_y, "y/F");
45  tcount->Branch("N", &t_N, "N/F");
46  tcount->Branch("stdev_N", &t_stdev_N, "stdev_N/F");
47 
48  /* copy beam parameters */
49  t_pbeam_lepton = ebeam_e;
50  t_pbeam_proton = ebeam_p;
51 
52  /* center of mass energy */
53  t_s = 4 * t_pbeam_lepton * t_pbeam_proton;
54 
55  /* loop over all bins */
56  for ( int bin_x = 1; bin_x <= hxQ2->GetNbinsX(); bin_x++ )
57  {
58 
59  for ( int bin_y = 1; bin_y <= hxQ2->GetNbinsY(); bin_y++ )
60  {
61 
62 
63  t_x = TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) ) )
64  + ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) + hxQ2->GetXaxis()->GetBinWidth(bin_x) ) ) ) );
65 
66  t_Q2 = TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_y) ) )
67  + ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_y) + hxQ2->GetYaxis()->GetBinWidth(bin_y) ) ) ) );
68 
69  t_y = t_Q2 / ( t_x * t_s );
70 
71  t_N = hxQ2->GetBinContent( bin_x, bin_y ) * lumi_scaling;
72 
73  /* skip kinematics bins wth y > 0.95 and y < 1e-2 */
74  if ( t_y > 0.95 || t_y < 1e-2 )
75  continue;
76 
77  t_stdev_N = 1./(sqrt(t_N));
78 
79  tcount->Fill();
80 
81  std::cout.precision(2);
82 
83  cout << "lepton = " << std::fixed << t_pbeam_lepton
84  << " x proton = " << std::fixed << t_pbeam_proton
85  << " , sqrt(s) = " << std::fixed << sqrt( t_s )
86  << " , x = " << std::scientific << t_x
87  << " , Q2 = " << std::scientific << t_Q2
88  << " , y = " << std::fixed << t_y
89  << " , N = " << std::scientific << t_N
90  << endl;
91 
92 
93  }
94 
95  }
96 
97  TCanvas *c1 = new TCanvas();
98  hxQ2->Draw("colz");
99  c1->SetLogx(1);
100  c1->SetLogy(1);
101  c1->SetLogz(1);
102 
103  TCanvas *c2 = new TCanvas();
104  tcount->Draw("N:Q2", "x > 1e-2 && x < 2e-2");
105 
106  tcount->Write();
107  fout->Close();
108 
109  return 0;
110 }