HFDiJetMomImbalance/macros/Draw_BDiJetImbalance_ncoll.C

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // Draw the Di b-jet <xj> vs Ncoll
0004 //
0005 // Useful references:
0006 // CMS di b-jet: http://cds.cern.ch/record/2202805/files/HIN-16-005-pas.pdf
0007 //
0008 ////////////////////////////////////////////////////////////////////////////////
0009 //
0010 // Darren McGlinchey
0011 // 29 Jun 2017
0012 //
0013 ////////////////////////////////////////////////////////////////////////////////
0014
0015 #include <TFile.h>
0016 #include <TChain.h>
0017 #include <TTree.h>
0018 #include <TH1.h>
0019 #include <TCanvas.h>
0020 #include <TLatex.h>
0021 #include <TStyle.h>
0022 #include <TMath.h>
0023 #include <TLegend.h>
0024 #include <TRandom3.h>
0025 #include <TString.h>
0026 #include <TSystem.h>
0027 #include <TGraphErrors.h>
0028
0029 #include <iostream>
0030
0031 using namespace std;
0032
0033 void calc_mxj(TH1D* htmp, double &mxj, double &mxj_e)
0034 {
0035
0036   double sum = 0;
0037   double w = 0;
0038   double err = 0;
0039   for (int i = 1; i <= htmp->GetNbinsX(); i++)
0040   {
0041     double c = htmp->GetBinContent(i);
0042     if (c > 0)
0043     {
0044       sum += c * htmp->GetBinCenter(i);
0045       w += c;
0046       err += TMath::Power(c * htmp->GetBinError(i), 2);
0047     }
0048   }
0049   double mean = sum / w;
0050   err = TMath::Sqrt(err) / w;
0051
0052   mxj = mean;
0053   mxj_e = err;
0054 }
0055
0056 void Draw_BDiJetImbalance_ncoll()
0057 {
0058   //==========================================================================//
0059   // SET RUNNING CONDITIONS
0060   //==========================================================================//
0061   bool print_plots = true;
0062
0063   const char* ppFile = "dibjet_imbalance_histos_pp.root";
0064   const char* auauFile = "dibjet_imbalance_histos_AuAu0-10.root";
0065
0066   const int NCENT = 4;
0067   const char* centlim[] = {"0-10%", "10-20%", "20-40%", "40-60%", "60-92%"};
0068   double ncoll[] = {962, 603, 296, 94, 15};
0069   double lcent[] = {549, 344, 338, 107, 27}; // nn luminosity (for scaling)
0070
0071
0072   //==========================================================================//
0073   // DECLARE VARIABLES
0074   //==========================================================================//
0075
0076   TH1D* hxj_pp;
0077   TH1D* hxj_auau[NCENT];
0078
0079   TGraphErrors *gmxj = new TGraphErrors();
0080   gmxj->SetMarkerStyle(kFullCircle);
0081   gmxj->SetMarkerColor(kBlack);
0082   gmxj->SetLineColor(kBlack);
0083
0084
0085   double mxj_pp;
0086   double mxj_e_pp;
0087
0088   double mxj_auau[NCENT];
0089   double mxj_e_auau[NCENT];
0090
0091   //-- other
0092   TFile *fin;
0093   char name[500];
0094
0095   //==========================================================================//
0096   // READ RESULTS FROM FILE
0097   //==========================================================================//
0098   cout << endl;
0099   cout << "--> Reading results from file" << endl;
0100
0101   cout << "----> Reading pp from " << ppFile << endl;
0102   fin = TFile::Open(ppFile);
0103   if (!fin)
0104   {
0105     cout << "ERROR!! Unable to read " << ppFile << endl;
0106     return;
0107   }
0108
0109   sprintf(name, "hdijet_xj_fix");
0110   hxj_pp = (TH1D*) fin->Get(name);
0111   if ( !hxj_pp )
0112   {
0113     cout << "ERRRO! Unable to find " << name << " in " << ppFile << endl;
0114     return;
0115   }
0116   hxj_pp->SetDirectory(0);
0117   hxj_pp->SetName("hxj_pp");
0118
0119   fin->Close();
0120   delete fin;
0121
0122
0123
0124   cout << "----> Reading auau from " << auauFile << endl;
0125   fin = TFile::Open(auauFile);
0126   if (!fin)
0127   {
0128     cout << "ERROR!! Unable to read " << auauFile << endl;
0129     return;
0130   }
0131
0132   sprintf(name, "hdijet_xj_fix");
0133   hxj_auau[0] = (TH1D*) fin->Get(name);
0134   if ( !hxj_auau[0] )
0135   {
0136     cout << "ERRRO! Unable to find " << name << " in " << auauFile << endl;
0137     return;
0138   }
0139   hxj_auau[0]->SetDirectory(0);
0140   hxj_auau[0]->SetName("hxj_auau_0");
0141
0142   fin->Close();
0143   delete fin;
0144
0145   //==========================================================================//
0146   // SCALE UNCERTAINTIES FOR OTHER CENTRALITY BINS
0147   //==========================================================================//
0148   cout << endl;
0149   cout << "--> Scaling AuAu centrality bins" << endl;
0150
0151   for (int icent = 1; icent < NCENT; icent++)
0152   {
0153     sprintf(name, "hxj_auau_%i", icent);
0154     hxj_auau[icent] = (TH1D*) hxj_auau[0]->Clone(name);
0155
0156     double sf = 1./TMath::Sqrt(lcent[icent] / lcent[0]);
0157     for (int ix = 1; ix <= hxj_auau[icent]->GetNbinsX(); ix++)
0158     {
0159       double be = hxj_auau[icent]->GetBinError(ix);
0160
0161       hxj_auau[icent]->SetBinError(ix, be * sf);
0162     } // ix
0163   }
0164
0165
0166   //==========================================================================//
0167   // CALCULATE <XJ>
0168   //==========================================================================//
0169   cout << endl;
0170   cout << "--> Calculating <x_j>" << endl;
0171
0172   //-- pp
0173   calc_mxj(hxj_pp, mxj_pp, mxj_e_pp);
0174
0175   gmxj->SetPoint(0, 1, mxj_pp);
0176   gmxj->SetPointError(0, 0, mxj_e_pp);
0177
0178   //-- AuAu
0179   for (int i = 0; i < NCENT; i++)
0180   {
0181     calc_mxj(hxj_auau[i], mxj_auau[i], mxj_e_auau[i]);
0182
0183     gmxj->SetPoint(i + 1, ncoll[i], mxj_auau[i]);
0184     gmxj->SetPointError(i + 1, 0, mxj_e_auau[i]);
0185
0186   } // i
0187
0188
0189   //==========================================================================//
0190   // PLOT OBJECTS
0191   //==========================================================================//
0192
0193   TH1D* haxis = new TH1D("haxis", ";N_{coll};#LTx_{j}#GT", 1200, -100, 1100);
0194   // haxis->SetMinimum(0.551);
0195   // haxis->SetMaximum(0.749);
0196   haxis->SetMinimum(0.601);
0197   haxis->SetMaximum(0.779);
0198
0199   TLatex lbl;
0200   lbl.SetTextAlign(22);
0201
0202   //==========================================================================//
0203   // PLOT
0204   //==========================================================================//
0205
0206   TCanvas *cxj = new TCanvas("cxj", "xj", 800, 600);
0207
0208   cxj->cd();
0209   // gPad->SetLogx();
0210   haxis->Draw();
0211
0212   gmxj->Draw("P");
0213
0214   // cent labels
0215   lbl.DrawLatex(1, mxj_pp + 0.02, "p+p");
0216
0217   for (int i = 0; i < NCENT; i++)
0218     lbl.DrawLatex(ncoll[i], mxj_auau[i] + 0.02 + 0.01*(i%2), centlim[i]);
0219
0220   TLegend *leg = new TLegend(0.10, 0.20, 0.70, 0.55);
0221   leg->SetFillStyle(0);
0222   leg->SetBorderSize(0);
0223   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation", "");
0224   leg->AddEntry("", "#sqrt{s_{_{NN}}} = 200 GeV", "");
0225   leg->AddEntry("", "Di b-jets (Pythia8, CTEQ6L)", "");
0226   leg->AddEntry("", "anti-k_{T}, R = 0.4, |#eta|<0.7", "");
0227   leg->AddEntry("", "p+p: 200 pb^{-1}, 60% Eff., 40% Pur.", "");
0228   leg->AddEntry("", "Au+Au: 240B Events, 40% Eff., 40% Pur.", "");
0229   // leg->AddEntry("", "200 pb^{-1} p+p, 240B Au+Au", "");
0230   leg->Draw("same");
0231
0232   //==========================================================================//
0233   // PRINT PLOTS
0234   //==========================================================================//
0235   if ( print_plots )
0236   {
0237     cout << endl;
0238     cout << "--> Printing results" << endl;
0239
0240     cout << "p+p <xj> = " << mxj_pp << " +/- " << mxj_e_pp << endl;
0241     for (int i = 0; i < NCENT; i++)
0242     {
0243       cout << "Au+Au " << centlim[i] << " <xj> = " << mxj_auau[i] << " +/- " << mxj_e_auau[i] << endl;
0244     }
0245
0246     cxj->Print("dibjet_imbalance_ncoll.pdf");
0247     cxj->Print("dibjet_imbalance_ncoll.C");
0248     cxj->Print("dibjet_imbalance_ncoll.root");
0249   }
0250
0251
0252
0253
0254
0255 }