EICMCAnalysis/analysis_eicsmear/hera_dis_countbins.C

0001 int
0002 hera_dis_countbins()
0003 {
0004   TFile *fin = new TFile("hera_dis_histo.root","OPEN");
0005   THnSparse *hfull = (THnSparse*)fin->Get("hn_dis");
0006
0007   TH2F* hxQ2 = (TH2F*)hfull->Projection(1,0);
0008
0009   /* beam energies */
0010   float ebeam_e = 27.5;
0011   float ebeam_p = 820;
0012
0013   /* Pythia generated events, cross section, and luminosity */
0014   float pythia_ngen = 999999;
0015   float pythia_xsec = 0.80120995103272474; // in microbarn
0016   float convert_microbarn_to_femtobarn = 1e9;
0017   float pythia_lumi = pythia_ngen / ( pythia_xsec * convert_microbarn_to_femtobarn );
0018
0019   /* target luminosity and scaling factor */
0020   float target_lumi = 0.5; // in inverse femtobarn
0021   float lumi_scaling = target_lumi / pythia_lumi;
0022
0023   cout << "Pythia luminosity:  " << pythia_lumi << " fb^-1" << endl;
0024   cout << "Target luminosity:  " << target_lumi << " fb^-1" << endl;
0025   cout << "Luminosity scaling: " << lumi_scaling << endl;
0026
0027   /* create tree to store information */
0028   TFile *fout = new TFile("hera_dis_tree.root", "RECREATE");
0029
0030   TTree *tcount = new TTree("tcount", "A tree with counts in kinematics bins");
0031   float t_pbeam_lepton = 0;
0032   float t_pbeam_proton = 0;
0033   float t_s = 0;
0034   float t_x = 0;
0035   float t_Q2 = 0;
0036   float t_y = 0;
0037   float t_N = 0;
0038   float t_stdev_N = 0;
0039   tcount->Branch("pbeam_lepton", &t_pbeam_lepton, "pbeam_lepton/F");
0040   tcount->Branch("pbeam_proton", &t_pbeam_proton, "pbeam_proton/F");
0041   tcount->Branch("s", &t_s, "s/F");
0042   tcount->Branch("x", &t_x, "x/F");
0043   tcount->Branch("Q2", &t_Q2, "Q2/F");
0044   tcount->Branch("y", &t_y, "y/F");
0045   tcount->Branch("N", &t_N, "N/F");
0046   tcount->Branch("stdev_N", &t_stdev_N, "stdev_N/F");
0047
0048   /* copy beam parameters */
0049   t_pbeam_lepton = ebeam_e;
0050   t_pbeam_proton = ebeam_p;
0051
0052   /* center of mass energy */
0053   t_s = 4 * t_pbeam_lepton * t_pbeam_proton;
0054
0055   /* loop over all bins */
0056   for ( int bin_x = 1; bin_x <= hxQ2->GetNbinsX(); bin_x++ )
0057     {
0058
0059       for ( int bin_y = 1; bin_y <= hxQ2->GetNbinsY(); bin_y++ )
0060     {
0061
0062
0063       t_x = TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) ) )
0064                      + ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) + hxQ2->GetXaxis()->GetBinWidth(bin_x) ) ) ) );
0065
0066       t_Q2 = TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_y) ) )
0067                       + ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_y) + hxQ2->GetYaxis()->GetBinWidth(bin_y) ) ) ) );
0068
0069       t_y = t_Q2 / ( t_x * t_s );
0070
0071       t_N = hxQ2->GetBinContent( bin_x, bin_y ) * lumi_scaling;
0072
0073       /* skip kinematics bins wth y > 0.95 and y < 1e-2 */
0074       if ( t_y > 0.95 || t_y < 1e-2 )
0075         continue;
0076
0077       t_stdev_N = 1./(sqrt(t_N));
0078
0079       tcount->Fill();
0080
0081       std::cout.precision(2);
0082
0083       cout << "lepton = " << std::fixed << t_pbeam_lepton
0084            << " x proton = " << std::fixed << t_pbeam_proton
0085            << " , sqrt(s) = " << std::fixed << sqrt( t_s )
0086            << " , x = " << std::scientific << t_x
0087            << " , Q2 = " << std::scientific << t_Q2
0088            << " , y = " << std::fixed << t_y
0089            << " , N = " << std::scientific << t_N
0090            << endl;
0091
0092
0093     }
0094
0095     }
0096
0097   TCanvas *c1 = new TCanvas();
0098   hxQ2->Draw("colz");
0099   c1->SetLogx(1);
0100   c1->SetLogy(1);
0101   c1->SetLogz(1);
0102
0103   TCanvas *c2 = new TCanvas();
0104   tcount->Draw("N:Q2", "x > 1e-2 && x < 2e-2");
0105
0106   tcount->Write();
0107   fout->Close();
0108
0109   return 0;
0110 }