EICMCAnalysis/analysis_eicsmear/calculate_g1_projection_per_bin.C

0001 #include "calculate_g1.C"
0002
0003 int
0004 calculate_g1_projection_per_bin(TString fin_name)
0005 {
0006   TFile *fin = new TFile(fin_name,"OPEN");
0007   THnSparse *hfull = (THnSparse*)fin->Get("hn_dis_electron_accept");
0008
0009   TH2F* hxQ2 = (TH2F*)hfull->Projection(1,0);
0010   hxQ2->SetName("hxQ2");
0011
0012   /* beam energies */
0013   float ebeam_e = 18;
0014   float ebeam_p = 275;
0015
0016   /* Retrieve Pythia generated events luminosity information */
0017   /* @TODO: Get this from string in file */
0018   TObjString* gen_crossSection_s = (TObjString*)fin->Get("crossSection");
0019   TObjString* gen_nEvents_s = (TObjString*)fin->Get("nEvents");
0020   TObjString* gen_lumi_ifb_s = (TObjString*)fin->Get("luminosity");
0021
0022   float gen_crossSection_mb = (TString((gen_crossSection_s)->GetName())).Atof(); // microbarn (10^-6) from pythiaeRHIC
0023   float gen_nEvents = (TString((gen_nEvents_s)->GetName())).Atof();
0024   float gen_lumi_ifb = (TString((gen_lumi_ifb_s)->GetName())).Atof();
0025
0026   float target_lumi_ifb = 10.;
0027
0028   cout << "Pythia cross section (microbarn):  " << gen_crossSection_mb << " mb" << endl;
0029   //cout << "Pythia cross section (femtobarn):  " << gen_crossSection_fb << " fb" << endl;
0030   cout << "Pythia generated (scaled) luminosity:  " << gen_lumi_ifb << " fb^-1" << endl;
0031   cout << "Target luminosity:  " << target_lumi_ifb << " fb^-1" << endl;
0032
0033   /* create tree to store information */
0034   TTree *tcount = new TTree("tcount", "A tree with counts in kinematics bins");
0035   float t_pbeam_lepton = 0;
0036   float t_pbeam_proton = 0;
0037   float t_lumi_ifb = 0;
0038   float t_s = 0;
0039   float t_x = 0;
0040   float t_Q2 = 0;
0041   float t_y = 0;
0042   float t_N = 0;
0043   float t_stdev_N = 0;
0044   float t_stdev_g1 = 0;
0045   float t_stdev_A1 = 0;
0046   tcount->Branch("pbeam_lepton", &t_pbeam_lepton, "pbeam_lepton/F");
0047   tcount->Branch("pbeam_proton", &t_pbeam_proton, "pbeam_proton/F");
0048   tcount->Branch("luminosity", &t_lumi_ifb, "luminosity/F"); // in inverse femtobarn
0049   tcount->Branch("s", &t_s, "s/F");
0050   tcount->Branch("x", &t_x, "x/F");
0051   tcount->Branch("Q2", &t_Q2, "Q2/F");
0052   tcount->Branch("y", &t_y, "y/F");
0053   tcount->Branch("N", &t_N, "N/F");
0054   tcount->Branch("stdev_N", &t_stdev_N, "stdev_N/F");
0055   tcount->Branch("stdev_g1", &t_stdev_g1, "stdev_g1/F");
0056   tcount->Branch("stdev_A1", &t_stdev_A1, "stdev_A1/F");
0057
0058   /* copy beam parameters */
0059   t_pbeam_lepton = ebeam_e;
0060   t_pbeam_proton = ebeam_p;
0061   t_lumi_ifb     = gen_lumi_ifb;
0062
0063   /* center of mass energy */
0064   t_s = 4 * ebeam_e * ebeam_p;
0065
0066   /* collect all x-bins */
0067   set<float> s_binc_x;
0068
0069   /* loop over all bins */
0070   for ( int bin_x = 1; bin_x <= hxQ2->GetNbinsX(); bin_x++ )
0071     {
0072       for ( int bin_y = 1; bin_y <= hxQ2->GetNbinsY(); bin_y++ )
0073     {
0074       t_x = TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) ) )
0075                      + ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) + hxQ2->GetXaxis()->GetBinWidth(bin_x) ) ) ) );
0076
0077       t_Q2 = TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_y) ) )
0078                       + ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_y) + hxQ2->GetYaxis()->GetBinWidth(bin_y) ) ) ) );
0079
0080       t_y = t_Q2 / ( t_x * t_s );
0081
0082       t_N = (float)hxQ2->GetBinContent( bin_x, bin_y ) * (target_lumi_ifb / gen_lumi_ifb);
0083
0084       /* skip bins with no entries */
0085       if ( t_N < 1 )
0086         continue;
0087
0088       t_stdev_N = 1./(sqrt(t_N));
0089       t_stdev_g1 = get_g1_sigma( t_x, t_Q2, t_y, t_N, 0.000511 );
0090       t_stdev_A1 = get_A1_sigma( t_x, t_Q2, t_y, t_N, 0.000511 );
0091
0092       tcount->Fill();
0093       s_binc_x.insert(t_x);
0094
0095       /* print values */
0096       std::cout.precision(3);
0097
0098       cout << "lepton = " << std::fixed << t_pbeam_lepton
0099            << " x proton = " << std::fixed << t_pbeam_proton
0100            << " , L = " << std::fixed << t_lumi_ifb << " fb^-1"
0101            << " , sqrt(s) = " << std::fixed << sqrt( t_s ) << " GeV"
0102            << " , x = " << std::scientific << t_x
0103            << " , Q2 = " << std::scientific << t_Q2
0104            << " , y = " << std::fixed << t_y
0105            << " , N = " << std::scientific << t_N
0106            << " , g1 = " << std::fixed << -1
0107            << " , g1_stdev = " << std::fixed << t_stdev_g1
0108            << " , A1_stdev = " << std::fixed << t_stdev_A1
0109            << endl;
0110     }
0111     }
0112
0113   /* create tree to store information */
0114   TFile *fout = new TFile("output/eic_sphenix_dis_tree.root", "RECREATE");
0115   fout->cd();
0116
0117   /* save input histogram in output */
0118   hxQ2->Write();
0119
0120   /* Write tree to output */
0121   tcount->Write();
0122
0123   /* copy string */
0124   gen_lumi_ifb_s->Write("luminosity");
0125
0126   /* close file */
0127   fout->Close();
0128
0129   return 0;
0130 }