EICMCAnalysis/analysis_eicsmear/eic_sphenix_sidis_countbins.C

0001 int
0002 eic_sphenix_sidis_countbins()
0003 {
0004   TFile *fin = new TFile("output/eic_sphenix_dis_histo_1M.root","OPEN");
0005   THnSparse *hfull = (THnSparse*)fin->Get("hn_sidis_pion_plus_accept");
0006   THnSparse *hfull_fullaccept = (THnSparse*)fin->Get("hn_sidis_pion_plus");
0007
0008   TH2F* hxQ2 = (TH2F*)hfull->Projection(1,0);
0009   TH2F* hxQ2_fullaccept = (TH2F*)hfull_fullaccept->Projection(1,0);
0010
0011   /* beam energies */
0012   float ebeam_e = 20;
0013   float ebeam_p = 250;
0014
0015   /* Pythia generated events, cross section, and luminosity */
0016   float pythia_ngen = 1e6;
0017   float pythia_xsec = 0.60785660255324614; // in microbarn
0018   float convert_microbarn_to_femtobarn = 1e9;
0019   float pythia_lumi = pythia_ngen / ( pythia_xsec * convert_microbarn_to_femtobarn );
0020
0021   /* target luminosity and scaling factor */
0022   float target_lumi = pythia_lumi; // in inverse femtobarn
0023   float lumi_scaling = target_lumi / pythia_lumi;
0024
0025   cout << "Pythia luminosity:  " << pythia_lumi << " fb^-1" << endl;
0026   cout << "Target luminosity:  " << target_lumi << " fb^-1" << endl;
0027   cout << "Luminosity scaling: " << lumi_scaling << endl;
0028
0029   /* Minimum number of bin entries to accept bin */
0030   int Nmin = 1;
0031
0032   /* create tree to store information */
0033   TTree *tcount = new TTree("tcount", "A tree with counts in kinematics bins");
0034   float t_pbeam_lepton = 0;
0035   float t_pbeam_proton = 0;
0036   float t_s = 0;
0037   float t_x = 0;
0038   float t_Q2 = 0;
0039   float t_y = 0;
0040   float t_z = 0;
0041   float t_pT = 0;
0042   float t_N = 0;
0043   float t_stdev_N = 0;
0044   tcount->Branch("pbeam_lepton", &t_pbeam_lepton, "pbeam_lepton/F");
0045   tcount->Branch("pbeam_proton", &t_pbeam_proton, "pbeam_proton/F");
0046   tcount->Branch("s", &t_s, "s/F");
0047   tcount->Branch("x", &t_x, "x/F");
0048   tcount->Branch("Q2", &t_Q2, "Q2/F");
0049   tcount->Branch("y", &t_y, "y/F");
0050   tcount->Branch("z", &t_z, "z/F");
0051   tcount->Branch("pT", &t_pT, "pT/F");
0052   tcount->Branch("N", &t_N, "N/F");
0053   tcount->Branch("stdev_N", &t_stdev_N, "stdev_N/F");
0054
0055   /* copy beam parameters */
0056   t_pbeam_lepton = ebeam_e;
0057   t_pbeam_proton = ebeam_p;
0058
0059   /* center of mass energy */
0060   t_s = 4 * ebeam_e * ebeam_p;
0061
0062
0063   /* Get all axis to loop over */
0064   TAxis *ax_x = hfull->GetAxis(0);
0065   TAxis *ax_Q2 = hfull->GetAxis(1);
0066   TAxis *ax_z = hfull->GetAxis(2);
0067   TAxis *ax_pT = hfull->GetAxis(3);
0068
0069   ax_x->Print();
0070   ax_Q2->Print();
0071   ax_z->Print();
0072   ax_pT->Print();
0073
0074   /* collect all x-bins */
0075   //  set<float> s_binc_x;
0076
0077   /* loop over all bins */
0078   for ( int bin_x = 1; bin_x <= ax_x->GetNbins(); bin_x++ )
0079     {
0080       for ( int bin_Q2 = 1; bin_Q2 <= ax_Q2->GetNbins(); bin_Q2++ )
0081         {
0082           t_x = TMath::Power(10, 0.5 * ( ( TMath::Log10( ax_x->GetBinLowEdge(bin_x) ) )
0083                                          + ( TMath::Log10( ax_x->GetBinLowEdge(bin_x) + ax_x->GetBinWidth(bin_x) ) ) ) );
0084
0085           t_Q2 = TMath::Power(10, 0.5 * ( ( TMath::Log10( ax_Q2->GetBinLowEdge(bin_Q2) ) )
0086                                           + ( TMath::Log10( ax_Q2->GetBinLowEdge(bin_Q2) + ax_Q2->GetBinWidth(bin_Q2) ) ) ) );
0087
0088           /* Calculate inelasticity y */
0089           t_y = t_Q2 / ( t_x * t_s );
0090
0091           /* skip kinematics bins wth y > 0.95 and y < 1e-2 */
0092           if ( t_y > 0.95 || t_y < 1e-2 )
0093             continue;
0094
0095           for ( int bin_z = 1; bin_z <= ax_z->GetNbins(); bin_z++ )
0096             {
0097
0098               t_z = ax_z->GetBinCenter( bin_z );
0099
0100               for ( int bin_pT = 1; bin_pT <= ax_pT->GetNbins(); bin_pT++ )
0101                 {
0102
0103                   t_pT = ax_pT->GetBinCenter( bin_pT );
0104
0105                   /* Get bin entries */
0106           int binloc[4] = { bin_x, bin_Q2, bin_z, bin_pT };
0107                   t_N = hfull->GetBinContent( binloc ) * lumi_scaling;
0108
0109                   /* skip bins with small number of entries */
0110                   if ( t_N < Nmin )
0111                     continue;
0112
0113                   t_stdev_N = 1./(sqrt(t_N));
0114
0115                   tcount->Fill();
0116                   //s_binc_x.insert(t_x);
0117
0118                   /* print values */
0119                   std::cout.precision(2);
0120
0121           bool toscreen = false;
0122           if ( toscreen )
0123             {
0124               cout    //     <<"lepton = " << std::fixed << t_pbeam_lepton
0125             //     << " x proton = " << std::fixed << t_pbeam_proton
0126             //     << " , sqrt(s) = " << std::fixed << sqrt( t_s )
0127             << " , x = " << std::scientific << t_x
0128             << " , Q2 = " << std::scientific << t_Q2
0129             << " , y = " << std::fixed << t_y
0130             << " , z = " << std::fixed << t_z
0131             << " , pT = " << std::fixed << t_pT
0132             << " , N = " << std::scientific << t_N
0133             << endl;
0134             }
0135                 }
0136             }
0137         }
0138     }
0139
0140   /* Prepare TPaveText for plots */
0141   TString str_ebeam = TString::Format("%.0f GeV x %.0f GeV", ebeam_e, ebeam_p );
0142   TString str_lumin = TString::Format("L = %.4f fb^{-1}", target_lumi );
0143
0144   TPaveText *pt_ebeam_lumi_ul = new TPaveText(1e-5,1e3,1e-3,1e4,"none");
0145   pt_ebeam_lumi_ul->SetFillStyle(0);
0146   pt_ebeam_lumi_ul->SetLineColor(0);
0147   pt_ebeam_lumi_ul->AddText(str_ebeam);
0148   pt_ebeam_lumi_ul->AddText(str_lumin);
0149
0150   TPaveText *pt_ebeam_lumi_ll = new TPaveText(1e2,45,1e3,50,"none");
0151   pt_ebeam_lumi_ll->SetFillStyle(0);
0152   pt_ebeam_lumi_ll->SetLineColor(0);
0153   pt_ebeam_lumi_ll->AddText(str_ebeam);
0154   pt_ebeam_lumi_ll->AddText(str_lumin);
0155
0156   /* Prepare inelasticity (y) cut lines for plots */
0157   TF1 *f_y095 = new TF1("f_y095", "4*x*[0]*[1]*[2]", 1e-5, 1);
0158   f_y095->SetParameter( 0, ebeam_e);
0159   f_y095->SetParameter( 1, ebeam_p);
0160   f_y095->SetParameter( 2, 0.95);
0161   TF1 *f_y001 = (TF1*)f_y095->Clone("f_y01");
0162   f_y001->SetParameter(2 , 0.01);
0163
0164   /* make x-Q2 plot */
0165   TCanvas *c1 = new TCanvas();
0166   c1->SetRightMargin(0.12);
0167   hxQ2->Draw("colz");
0168   c1->SetLogx(1);
0169   c1->SetLogy(1);
0170   c1->SetLogz(1);
0171
0172   f_y095->Draw("same");
0173   f_y001->Draw("same");
0174
0175   pt_ebeam_lumi_ul->Draw();
0176   gPad->RedrawAxis();
0177
0178   /* make x-Q2 plot for 'perfect' acceptance */
0179   TCanvas *c4 = new TCanvas();
0180   c4->SetRightMargin(0.12);
0181   hxQ2_fullaccept->Draw("colz");
0182   c4->SetLogx(1);
0183   c4->SetLogy(1);
0184   c4->SetLogz(1);
0185
0186   f_y095->Draw("same");
0187   f_y001->Draw("same");
0188
0189   pt_ebeam_lumi_ul->Draw();
0190   gPad->RedrawAxis();
0191
0192   /* make x-Q2 acceptance fraction pot */
0193   TCanvas *c3 = new TCanvas();
0194   c3->SetRightMargin(0.12);
0195   TH2F* hxQ2_acceptance_ratio = hxQ2->Clone("x_Q2_acceptance_ratio");
0196   hxQ2_acceptance_ratio->GetZaxis()->SetNdivisions(505);
0197   hxQ2_acceptance_ratio->Divide(hxQ2_fullaccept);
0198   hxQ2_acceptance_ratio->Draw("colz");
0199   c3->SetLogx(1);
0200   c3->SetLogy(1);
0201
0202   f_y095->Draw("same");
0203   f_y001->Draw("same");
0204
0205   pt_ebeam_lumi_ul->Draw();
0206   gPad->RedrawAxis();
0207
0208   /* create tree to store information */
0209   TFile *fout = new TFile("output/eic_sphenix_sidis_tree.root", "RECREATE");
0210
0211   /* Write tree to output */
0212   tcount->Write();
0213   fout->Close();
0214
0215   return 0;
0216 }