LightFlavorRatios/yield_and_ratios/plot_results_systematics.C

0001 void plot_results_systematics()
0002 {
0003
0004   // get systematics from the component histograms
0005   //======================================
0006   std::vector<std::string> fin;
0007   std::vector<std::string> hname;
0008
0009   // daughter track pT cut systematic
0010   fin.push_back("/sphenix/tg/tg01/hf/gregoryottino/lightFlavorPpg16/analysis/LightFlavorRatios/systematic_uncertainty/SystematicDetermination_output_ptmin.root");
0011   hname.push_back("hDifference");
0012
0013   // pT smearing systematic
0014   fin.push_back("/sphenix/tg/tg01/hf/gregoryottino/lightFlavorPpg16/analysis/LightFlavorRatios/systematic_uncertainty/SystematicDetermination_output_ptsmear.root");
0015   hname.push_back("hDifference");
0016
0017   // cuts efficiency systematic
0018   fin.push_back("/sphenix/user/mjpeters/analysis/LightFlavorRatios/swimming_correction/LamdbaKsCutEfficiency_200MeV_hists.root");
0019   hname.push_back("cuteff_syserr_pT");
0020
0021   // The peak extraction systematic
0022   fin.push_back("/sphenix/user/mjpeters/analysis/LightFlavorRatios/yield_and_ratios/fits.root");
0023   hname.push_back("ratio_fit_syserr_vspT");
0024
0025   int nsys = fin.size();
0026   std::cout << " there are " << fin.size() << " systematics input files " << std::endl;
0027
0028   TH1F *hsys[nsys];
0029   int nbinsmax = 0;
0030   int nfilemax = -1;
0031
0032   for(int i=0;i<nsys;++i)
0033     {
0034       TFile *f = new TFile(fin[i].c_str());
0035       if(!f)
0036     {
0037       std::cout << "Failed to open file named " << fin[i] << std::endl;
0038       return;
0039     }
0040       f->GetObject(hname[i].c_str(), hsys[i]);
0041       if(!hsys[i])
0042     {
0043       std::cout << "Failed to get histogram " << hname[i] << " from file " << fin[i] << std::endl;
0044       return;
0045     }
0046       std::cout << " file " << i << " named " << fin[i] << " has histogram " << hname[i] << " with " << hsys[i]->GetNbinsX()  << " bins " << std::endl;
0047       if(hsys[i]->GetNbinsX() > nbinsmax)
0048     {
0049       nbinsmax=hsys[i]->GetNbinsX();
0050       nfilemax = i;
0051     }
0052     }
0053   std::cout << std::endl << " nbinsmax is " << nbinsmax << " for " << fin[nfilemax] << std::endl;
0054
0055   double x[nbinsmax];
0056   for(int j = 0; j<nbinsmax; ++j)
0057     {
0058       x[j] = hsys[nfilemax]->GetBinCenter(j+1);
0059     }
0060
0061   double syscum[nbinsmax];
0062   for(int i = 0; i < nsys; ++i)  // histograms
0063     {
0064       for(int k=0; k < hsys[i]->GetNbinsX();++k)  // bins
0065     {
0066       double xd=hsys[i]->GetBinLowEdge(k+1);
0067       double xu=xd+hsys[i]->GetBinWidth(k+1);
0068       double sys = hsys[i]->GetBinContent(k+1);
0069       // decide which of the fine bins this goes into
0070       for(int j=0; j<nbinsmax;++j)
0071         {
0072           if(x[j] > xd && x[j] < xu)
0073         {
0074           syscum[k] += sys*sys;
0075           std::cout << " file " << i << " file bin " << k << " fine bin " << j << " sys " << sys << " syscum " << syscum[k] << std::endl;
0076         }
0077         }
0078     }
0079     }
0080
0081   // in addition we need to add a systematic for the tracking inefficiency correction to the ratio
0082   // we agreed on 3% flat in pT as a conservative estimate
0083   float track_inefficiency_percent = 0.03;
0084   for(int i=0;i<nbinsmax;++i)
0085     {
0086       syscum[i] += track_inefficiency_percent * track_inefficiency_percent;
0087       //  std::cout << " i " << i << " track_inefficiency_percent " << track_inefficiency_percent << " syscum is now " << syscum[i] << std::endl;
0088     }
0089
0090   // get the ratios
0091   //===========
0092   TFile* f = TFile::Open("/sphenix/user/mjpeters/analysis/LightFlavorRatios/yield_and_ratios/fits.root");
0093   if(!f)
0094     {
0095       std::cout << "Failed to open the file fits.root" << std::endl;
0096       return;
0097     }
0098
0099   TH1F *ratio_hist_pt;
0100   f->GetObject("lambdaKsratio_vspT_lambdafeeddowncorrected_effcorrected_geoacceptancecorrected_cutefficiencycorrected",ratio_hist_pt);
0101   double plotmax = ratio_hist_pt->GetMaximum() * 1.3;
0102
0103   ratio_hist_pt->SetMinimum(0);
0104   ratio_hist_pt->Draw();
0105
0106   // convert ratio histogram to a TGraphMultiErrors
0107   // extract the arrays
0108   int nbins = ratio_hist_pt->GetNbinsX();
0109   float pt[nbins];
0110   float ptelo[nbins];
0111   float ptehi[nbins];
0112   float ratio[nbins];
0113   float ratiostat[nbins];
0114   for (int ipt=0;ipt<nbins;++ipt)
0115   {
0116     pt[ipt] = ratio_hist_pt->GetXaxis()->GetBinCenter(ipt+1);
0117     float ptlo = ratio_hist_pt->GetXaxis()->GetBinLowEdge(ipt+1);
0118     float ptw = ratio_hist_pt->GetXaxis()->GetBinWidth(ipt+1);
0119     ptelo[ipt] = pt[ipt] - ptlo;
0120     ptehi[ipt] = ptlo + ptw - pt[ipt];
0121
0122     ratio[ipt] = ratio_hist_pt->GetBinContent(ipt+1);
0123     ratiostat[ipt] = ratio_hist_pt->GetBinError(ipt+1);
0124   }
0125
0126   // full systematic errors
0127   //=================
0128   double esyslo[nbinsmax];
0129   double esyshi[nbinsmax];
0130   for(int i=0;i<nbinsmax;++i)
0131     {
0132       esyslo[i] = ratio[i]*sqrt(syscum[i]);
0133       esyshi[i] = ratio[i]*sqrt(syscum[i]);
0134       std::cout << " i " << i << " pT " << pt[i] << " ratio " << ratio[i] << " syscum " << syscum[i] << " sqrt(syscum) " << sqrt(syscum[i]) << " esyslo " << esyslo[i] << " esyshi " << esyshi[i] << std::endl;
0135     }
0136
0137   // plot the ratio vs pT
0138   //===============
0139   TCanvas *c=new TCanvas("c","corrected double ratio plot", 50, 50, 800,600);
0140   gPad->SetBottomMargin(0.12);
0141   gPad->SetLeftMargin(0.12);
0142
0143 TGraphMultiErrors *gr = new TGraphMultiErrors("gr","Ratio vs pT", nbins,pt,ratio,ptelo,ptehi,ratiostat,ratiostat);
0144  gr->SetName("ratio_pt_sys");
0145  gr->SetMinimum(0);
0146  gr->SetMaximum(plotmax);
0147  gr->SetMarkerStyle(20);
0148  gr->GetXaxis()->SetRangeUser(0,5.0);
0149
0150  gr->AddYError(nbins,esyslo,esyshi);
0151  gr->SetMarkerStyle(22);
0152  gr->SetMarkerSize(1.5);
0153  gr->SetLineColor(kBlack);
0154  gr->GetAttLine(0)->SetLineColor(kBlack);
0155  gr->GetAttLine(1)->SetLineColor(kBlue);
0156  gr->GetAttFill(1)->SetFillStyle(0);
0157
0158  gr->GetXaxis()->SetTitle("p_{T} (Gev)");
0159  gr->GetYaxis()->SetTitle("#Lambda^{0}/2K_{s}^0 (arbitrary units)");
0160
0161  gr->Draw("APS ; Z ; 5 s=0.5");
0162
0163  c->SaveAs("ratio_allsyst.pdf");
0164
0165  TFile *fout = new TFile("temp_ratio_results_systematic.root","recreate");
0166  gr->Write();
0167  fout->Close();
0168
0169
0170
0171 }