hbcup/macro/plot_upcsim.C

0001 //
0002 // Make mass, pt etc plots for UPC from the simUPC DSTs
0003 //
0004
0005 void plot_upcsim()
0006 {
0007   const int rebin = 2;
0008   const int NFILES = 9;
0009   vector<std::string> fnames {
0010     "anaupc_cohjpsi.root",
0011     "anaupc_incohjpsi.root",
0012     "anaupc_cohjpsi_mumu.root",
0013     "anaupc_incohjpsi_mumu.root",
0014     "anaupc_qedhiw.root",
0015     "anaupc_qedhiw2.root",
0016     "anaupc_qedhiw_mumu.root",
0017     "anaupc_qedhiw2_mumu.root",
0018     "sHijing/anaupc_hijing.root"
0019     //"sHijing/anaupc_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000019-000000.root"
0020   };
0021
0022   int hijfile = NFILES - 1; // hijing file index, always last
0023
0024   //const double integ_lumi = 4.5e9;   // in inverse barns, pre-sPHENIX
0025   const double integ_lumi = 5.75e9;    // in inverse barns, 2025 AuAu
0026   const double stream_fraction = 1.0;  // fraction of crossings we stream
0027   const double upcsigma[] = { 27.8e-6, 29.5e-6, 27.7e-6, 29.5e-6, 6.374e-3, 17.166e-6, 6.238e-3, 17.190e-6, 7.2*(3.12e3/2e5) };  // in barns
0028   //const double upcsigma[] = { 27.8e-6, 29.5e-6, 27.7e-6, 29.5e-6, 6.374e-3, 17.166e-6, 6.238e-3, 17.190e-6, 7.2 };  // in barns
0029   //double nevents[] = { 1e5, 1e5, 1e6, 0, 0, 0, 0, 0 };  // not used.. we get nevents from datafiles
0030   double nevents[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };  // we get nevents from the datafiles
0031
0032   TFile *tfile[NFILES];
0033   TH1 *h_trig[NFILES];
0034
0035   TH1 *h_mass[NFILES];
0036   TH1 *h_pt[NFILES];
0037   TH1 *h_y[NFILES];
0038
0039   TH1 *h_mass_ls[NFILES];
0040   TH1 *h_pt_ls[NFILES];
0041   TH1 *h_y_ls[NFILES];
0042
0043   double scale[NFILES];   // scaling factor to normalize to same integrated luminosity
0044
0045   for (int ifile=0; ifile<NFILES; ifile++)
0046   {
0047     tfile[ifile] = new TFile(fnames[ifile].c_str(),"READ");
0048
0049     h_mass[ifile] = (TH1*)tfile[ifile]->Get("h_mass");
0050     h_y[ifile] = (TH1*)tfile[ifile]->Get("h_y");
0051     h_pt[ifile] = (TH1*)tfile[ifile]->Get("h_pt");
0052
0053     h_mass_ls[ifile] = (TH1*)tfile[ifile]->Get("h_mass_ls");
0054     h_y_ls[ifile] = (TH1*)tfile[ifile]->Get("h_y_ls");
0055     h_pt_ls[ifile] = (TH1*)tfile[ifile]->Get("h_pt_ls");
0056
0057     h_mass[ifile]->Sumw2();
0058     h_y[ifile]->Sumw2();
0059     h_pt[ifile]->Sumw2();
0060
0061     h_mass[ifile]->SetXTitle("mass (GeV/c^2)");
0062     h_y[ifile]->SetXTitle("y");
0063     h_pt[ifile]->SetXTitle("p_{T} (GeV/c)");
0064
0065     if ( h_mass_ls[ifile] != 0 )
0066     {
0067       h_mass_ls[ifile]->Sumw2();
0068       h_y_ls[ifile]->Sumw2();
0069       h_pt_ls[ifile]->Sumw2();
0070     }
0071
0072     if ( rebin>1 )
0073     {
0074       h_mass[ifile]->Rebin( rebin );
0075
0076       if ( h_mass_ls[ifile] != 0 )
0077       {
0078         h_mass_ls[ifile]->Rebin( rebin );
0079       }
0080     }
0081
0082     h_trig[ifile] = (TH1*)tfile[ifile]->Get("h_trig");
0083     nevents[ifile] = h_trig[ifile]->GetBinContent(1);
0084
0085     // Scaling factor for integrated luminosity expected in sPHENIX Run
0086     scale[ifile] = upcsigma[ifile]*integ_lumi/nevents[ifile];
0087
0088     cout << fnames[ifile] << "\t" << nevents[ifile] << "\t" << scale[ifile] << endl;
0089   }
0090
0091
0092   for (int ifile=0; ifile<NFILES; ifile++)
0093   {
0094     h_mass[ifile]->Scale( scale[ifile] );
0095     h_y[ifile]->Scale( scale[ifile] );
0096     h_pt[ifile]->Scale( scale[ifile] );
0097
0098     if ( h_mass_ls[ifile] != 0 )
0099     {
0100       h_mass_ls[ifile]->Scale( scale[ifile] );
0101       h_y_ls[ifile]->Scale( scale[ifile] );
0102       h_pt_ls[ifile]->Scale( scale[ifile] );
0103     }
0104   }
0105
0106   //=== Make Sums
0107   TH1 *h_tot_mass = (TH1*)h_mass[0]->Clone("h_tot_mass");           // all mass pairs, opp-sign
0108   h_tot_mass->SetXTitle("mass [GeV]");
0109   TH1 *h_tot_mass_ls = (TH1*)h_mass_ls[0]->Clone("h_tot_mass_ls");  // all mass pairs, like-sign
0110   for (int ifile=1; ifile<NFILES; ifile++)
0111   {
0112     h_tot_mass->Add( h_mass[ifile] );
0113
0114     if ( h_mass_ls[ifile] != 0 )
0115     {
0116       h_tot_mass_ls->Add( h_mass_ls[ifile] );
0117     }
0118   }
0119
0120   // J/Psi sums
0121   TH1 *h_jpsi_mass = (TH1*)h_mass[0]->Clone("h_jpsi_mass");           // all mass pairs, opp-sign
0122   h_jpsi_mass->SetXTitle("mass [GeV]");
0123   h_jpsi_mass->SetTitle("J/Psi");
0124   TH1 *h_jpsi_mass_ls = (TH1*)h_mass_ls[0]->Clone("h_jpsi_mass_ls");  // all mass pairs, like-sign
0125   h_jpsi_mass_ls->SetLineColor(4);
0126   for (int ifile=1; ifile<=3; ifile++)
0127   {
0128     h_jpsi_mass->Add( h_mass[ifile] );
0129
0130     if ( h_mass_ls[ifile] != 0 )
0131     {
0132       h_jpsi_mass_ls->Add( h_mass_ls[ifile] );
0133     }
0134   }
0135
0136   // QED sums
0137   TH1 *h_qed_mass = (TH1*)h_mass[4]->Clone("h_qed_mass");           // all mass pairs, opp-sign
0138   h_qed_mass->SetXTitle("mass [GeV]");
0139   h_qed_mass->SetTitle("QED");
0140   TH1 *h_qed_mass_ls = (TH1*)h_mass_ls[4]->Clone("h_qed_mass_ls");  // all mass pairs, like-sign
0141   h_qed_mass_ls->SetLineColor(4);
0142   for (int ifile=5; ifile<=7; ifile++)
0143   {
0144     h_qed_mass->Add( h_mass[ifile] );
0145
0146     if ( h_mass_ls[ifile] != 0 )
0147     {
0148       h_qed_mass_ls->Add( h_mass_ls[ifile] );
0149     }
0150   }
0151
0152   // Hijing (no sum needed)
0153
0154   //=== Mass Plots
0155   int ncv = 0;
0156   TCanvas *ac[100];
0157
0158   //== Total
0159   ac[ncv] = new TCanvas("cmass","mass",800,600);
0160   h_tot_mass->Draw("hist");
0161   h_tot_mass_ls->Draw("histsame");
0162   gPad->SetLogy(1);
0163   ncv++;
0164
0165   //== J/Psi
0166   ac[ncv] = new TCanvas("cmassjpsi","mass jpsi",800,600);
0167   h_jpsi_mass->Draw("hist");
0168   h_jpsi_mass_ls->Draw("histsame");
0169   //gPad->SetLogy(1);
0170   ncv++;
0171
0172   //== J/Psi (e+e- vs mu+mu-)
0173   ac[ncv] = new TCanvas("cmassjpsi_e_vs_mu","mass jpsi e vs mu",800,600);
0174   h_mass[2]->Draw("hist");
0175   for (int ifile=0; ifile<=3; ifile++)
0176   {
0177     h_mass[ifile]->SetLineColor(ifile+1);
0178     h_mass[ifile]->Draw("histsame");
0179   }
0180   TLegend *leg = new TLegend(0.7,0.7,0.9,0.9);
0181   //leg->SetHeader("The Legend Title","C"); // option "C" allows to center the header
0182   leg->AddEntry(h_mass[0],"Coh J/Psi#rightarrowe+e-","l");
0183   leg->AddEntry(h_mass[1],"Incoh J/Psi->e+e-","l");
0184   leg->AddEntry(h_mass[2],"Coh J/Psi->mu+mu-","l");
0185   leg->AddEntry(h_mass[3],"Incoh J/Psi->mu+mu-","l");
0186   leg->Draw();
0187   //gPad->SetLogy(1);
0188   ncv++;
0189
0190   //== QED
0191   ac[ncv] = new TCanvas("cmassqed","mass qed",800,600);
0192   h_qed_mass->Draw("hist");
0193   h_qed_mass_ls->Draw("histsame");
0194   gPad->SetLogy(1);
0195   ncv++;
0196
0197   /*
0198   for (int ifile=NFILES-1; ifile>=0; ifile--)
0199   {
0200     h_mass[ifile]->SetLineColor(ifile+2);
0201     h_mass[ifile]->Draw("histsame");
0202
0203     if ( h_mass_ls[ifile] != 0 )
0204     {
0205       h_mass_ls[ifile]->SetLineColor(ifile+70);
0206       h_mass_ls[ifile]->Draw("histsame");
0207     }
0208   }
0209   ncv++;
0210   */
0211
0212   ac[ncv] = new TCanvas("cmass_ls","mass, like vs unlike",800,600);
0213   h_tot_mass->Draw("hist");
0214   h_tot_mass_ls->SetLineColor(4);
0215   h_tot_mass_ls->Draw("histsame");
0216   ncv++;
0217
0218   ac[ncv] = new TCanvas("cmass_hijing_ls","mass from hijing events, like vs unlike",800,600);
0219   h_mass[hijfile]->Draw("hist");
0220   h_mass_ls[hijfile]->SetLineColor(4);
0221   h_mass_ls[hijfile]->Draw("histsame");
0222   ncv++;
0223
0224   //=== pt
0225   TH1 *h_tot_pt = (TH1*)h_pt[0]->Clone("h_tot_pt");
0226   h_tot_pt->SetXTitle("p_{T} [GeV/c]");
0227   for (int ifile=1; ifile<NFILES; ifile++)
0228   {
0229     h_tot_pt->Add( h_pt[ifile] );
0230   }
0231
0232   ac[ncv] = new TCanvas("cpt","pt",800,600);
0233   h_tot_pt->Draw("hist");
0234   for (int ifile=NFILES-1; ifile>=0; ifile--)
0235   {
0236     h_pt[ifile]->SetLineColor(ifile+2);
0237     h_pt[ifile]->Draw("histsame");
0238   }
0239   ncv++;
0240
0241   //=== y
0242   TH1 *h_tot_y = (TH1*)h_y[0]->Clone("h_tot_y");
0243   h_tot_y->SetXTitle("y");
0244   for (int ifile=1; ifile<NFILES; ifile++)
0245   {
0246     h_tot_y->Add( h_y[ifile] );
0247   }
0248
0249   ac[ncv] = new TCanvas("crap","rapidity",800,600);
0250   h_tot_y->Draw("hist");
0251   for (int ifile=NFILES-1; ifile>=0; ifile--)
0252   {
0253     h_y[ifile]->SetLineColor(ifile+2);
0254     h_y[ifile]->Draw("histsame");
0255   }
0256   ncv++;
0257
0258 }
0259
0260
0261
0262