EICMCAnalysis/analysis_eicsmear/plot_g1_projection.C

0001 #include "calculate_g1.C"
0002
0003 int
0004 plot_g1_projection(TString fin_name)
0005 {
0006   TFile *fin = new TFile(fin_name,"OPEN");
0007
0008   TH2F *hxQ2 = (TH2F*)fin->Get("hxQ2");
0009   TTree *tin = (TTree*)fin->Get("tcount");
0010
0011   /* beam energies */
0012   float ebeam_e = 18;
0013   float ebeam_p = 275;
0014
0015   /* Retrieve Pythia generated events luminosity information */
0016   TObjString* gen_lumi_ifb_s = (TObjString*)fin->Get("luminosity");
0017   float gen_lumi_ifb = (TString((gen_lumi_ifb_s)->GetName())).Atof();
0018
0019   cout << "Pythia luminosity:  " << gen_lumi_ifb << " fb^-1" << endl;
0020
0021   /* collect all x-bins */
0022   //set<float> s_binc_x;
0023
0024   /* Prepare TPaveText for plots */
0025   TString str_ebeam = TString::Format("%.0f GeV x %.0f GeV", ebeam_e, ebeam_p );
0026   TString str_lumin = TString::Format("L = %.4f fb^{-1}", gen_lumi_ifb );
0027
0028   TPaveText *pt_ebeam_lumi_ul = new TPaveText(1e-5,1e3,1e-3,1e4,"none");
0029   pt_ebeam_lumi_ul->SetFillStyle(0);
0030   pt_ebeam_lumi_ul->SetLineColor(0);
0031   pt_ebeam_lumi_ul->AddText(str_ebeam);
0032   pt_ebeam_lumi_ul->AddText(str_lumin);
0033
0034   TPaveText *pt_ebeam_lumi_ll = new TPaveText(1e2,45,1e3,50,"none");
0035   pt_ebeam_lumi_ll->SetFillStyle(0);
0036   pt_ebeam_lumi_ll->SetLineColor(0);
0037   pt_ebeam_lumi_ll->AddText(str_ebeam);
0038   pt_ebeam_lumi_ll->AddText(str_lumin);
0039
0040   /* Prepare inelasticity (y) cut lines for plots */
0041   TF1 *f_y095 = new TF1("f_y095", "4*x*[0]*[1]*[2]", 1e-5, 1);
0042   f_y095->SetParameter( 0, ebeam_e);
0043   f_y095->SetParameter( 1, ebeam_p);
0044   f_y095->SetParameter( 2, 0.95);
0045   TF1 *f_y001 = (TF1*)f_y095->Clone("f_y01");
0046   f_y001->SetParameter(2 , 0.01);
0047
0048   /* make x-Q2 plot */
0049   TCanvas *c1 = new TCanvas();
0050   c1->SetRightMargin(0.12);
0051   hxQ2->Draw("colz");
0052   c1->SetLogx(1);
0053   c1->SetLogy(1);
0054   c1->SetLogz(1);
0055
0056   f_y095->Draw("same");
0057   f_y001->Draw("same");
0058
0059   TPaveText *pt_N = new TPaveText(1e-5,1e2,1e-3,3e2,"none");
0060   pt_N->SetFillStyle(0);
0061   pt_N->SetLineColor(0);
0062   pt_N->AddText("z = # events");
0063
0064   pt_N->Draw();
0065   pt_ebeam_lumi_ul->Draw();
0066   gPad->RedrawAxis();
0067
0068
0069   /* plot g1 statistical uncertainty for different x-Q2 bins */
0070   TH2F* hxQ2_g1_sigma = hxQ2->Clone("x_Q2_g1_sigma");
0071   hxQ2_g1_sigma->GetZaxis()->SetNdivisions(505);
0072   hxQ2_g1_sigma->Reset();
0073
0074   hxQ2_g1_sigma->GetXaxis()->SetTitle("x");
0075   hxQ2_g1_sigma->GetYaxis()->SetTitle("Q^{2} (GeV^{2})");
0076
0077   /* Filling happens in the step below- draw afterwards */
0078
0079   //  tcount->Draw("Q2:x >> x_Q2_g1_sigma","stdev_g1 * (N > 0 && x >= 5e-5 && stdev_g1<1)");
0080
0081   /* normalize to number of entries in each bin */
0082   // hxQ2_g1_sigma->Divide(hxQ2);
0083
0084
0085   /* plot g1 vs Q2 for various x */
0086   TCanvas *c2 = new TCanvas("g1","",700,800);
0087   c2->SetLogx();
0088
0089   TH1F* hframe_g1 = (TH1F*)hxQ2->ProjectionY()->Clone("h1_g1");
0090   hframe_g1->Reset();
0091   hframe_g1->GetXaxis()->CenterTitle();
0092   hframe_g1->GetYaxis()->CenterTitle();
0093   hframe_g1->GetXaxis()->SetTitle("Q^{2} (GeV^{2})");
0094   //hframe_g1->GetYaxis()->SetTitle("g_{1}(x,Q^{2}) + const(x)");
0095   hframe_g1->GetYaxis()->SetTitle("const(x)");
0096   hframe_g1->GetXaxis()->SetRangeUser(0.99,1500);
0097   hframe_g1->GetYaxis()->SetRangeUser(2,50);
0098   hframe_g1->GetYaxis()->SetNdivisions(505);
0099
0100   hframe_g1->Draw("");
0101
0102   /* collect all x-bins */
0103   set<float> s_binc_x;
0104   for ( int bin_x = 1; bin_x <= hxQ2->GetNbinsX(); bin_x++ )
0105     {
0106       t_x = TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) ) )
0107                      + ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) + hxQ2->GetXaxis()->GetBinWidth(bin_x) ) ) ) );
0108       s_binc_x.insert(t_x);
0109     }
0110
0111   /* Prepare legend to record offsets */
0112   TPaveText *leg_offset_x = new TPaveText(0,0,1,1,"none");
0113   leg_offset_x->SetFillStyle(0);
0114   leg_offset_x->SetLineColor(0);
0115   //  leg_offset_x->SetMargin(0);
0116   leg_offset_x->SetTextSize(0.1);
0117
0118
0119
0120   /* draw graphs */
0121   TCanvas *ctmp = new TCanvas();
0122   float offset = 48;
0123
0124   cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
0125   cout << "Plot offsets: " << endl;
0126
0127   for ( set<float>::iterator itx = s_binc_x.begin();
0128     itx != s_binc_x.end(); itx++ )
0129     {
0130       /* Skip x < 5e-5 */
0131       if ( *itx < 5e-5 )
0132     continue;
0133
0134       /* print values */
0135       cout << "offset = " << (int)offset << " for x = " << *itx << endl;
0136
0137       /* Add to legend */
0138       TString str_offset = TString::Format("const(x = %.2g) = %d", *itx, (int)offset );
0139       leg_offset_x->AddText(str_offset);
0140
0141       /* Create graph */
0142       ctmp->cd();
0143
0144       unsigned npoints = tcount->GetEntries( TString::Format("x > 0.99*%f && x < 1.01*%f", *itx, *itx ) );
0145       tcount->Draw( TString::Format("%f:Q2:stdev_g1", offset),
0146             TString::Format("x > 0.99*%f && x < 1.01*%f", *itx, *itx ) );
0147
0148       TGraphErrors* gnew = new TGraphErrors( npoints, tcount->GetV2(), tcount->GetV1(), 0, tcount->GetV3() );
0149       gnew->SetMarkerColor(kRed);
0150
0151       c2->cd();
0152       gnew->Draw("PLsame");
0153
0154       offset -= 2;
0155
0156       /* Fill TH2 histogram */
0157       double x_j = 0;
0158       double Q2_j = 0;
0159       double g1sig_j = 0;
0160       double dummy_j = 0;
0161       for ( unsigned j = 0; j < gnew->GetN(); j++ )
0162     {
0163       gnew->GetPoint(j, Q2_j, dummy_j);
0164       g1sig_j = gnew->GetErrorY(j);
0165       x_j = *itx;
0166       hxQ2_g1_sigma->Fill(x_j, Q2_j, g1sig_j);
0167     }
0168     }
0169
0170   pt_ebeam_lumi_ll->Draw();
0171   gPad->RedrawAxis();
0172
0173   c2->Print("plots/g1_projection_eic_sphenix.eps");
0174
0175   TCanvas *c2_legend = new TCanvas("g1_legend","",200,800);
0176   leg_offset_x->Draw();
0177   c2_legend->Print("plots/g1_projection_eic_sphenix_legend.eps");
0178
0179   delete ctmp;
0180
0181
0182   /* Draw TH2 with g1 uncertainty defined above */
0183   TCanvas *c5 = new TCanvas();
0184   c5->SetRightMargin(0.12);
0185
0186   hxQ2_g1_sigma->Draw("colz");
0187   hxQ2_g1_sigma->GetZaxis()->SetRangeUser(1e-3,2e1);
0188   c5->SetLogx(1);
0189   c5->SetLogy(1);
0190   c5->SetLogz(1);
0191
0192   f_y095->Draw("same");
0193   f_y001->Draw("same");
0194
0195   TPaveText *pt_g1_stdev = new TPaveText(1e-5,1e2,1e-3,3e2,"none");
0196   pt_g1_stdev->SetFillStyle(0);
0197   pt_g1_stdev->SetLineColor(0);
0198   pt_g1_stdev->AddText("z = #sigma_{g_{1}}");
0199
0200   pt_g1_stdev->Draw();
0201   pt_ebeam_lumi_ul->Draw();
0202   gPad->RedrawAxis();
0203
0204   c5->Print("plots/g1_projection_eic_sphenix_2D.eps");
0205
0206   return 0;
0207 }