PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/convertToMyBins.C

   1 #include <TFile.h>
   2 #include <TH1.h>
   3 #include <TMath.h>
   4 #include <TGraphErrors.h>
   5
   6 #include "my_tools.C"
   7
   8 #include <iostream>
   9
  10 using namespace std;
  11
  12 /*
  13   Info:
  14
  15   Method to convert charged pT histograms to have the same histogram binning
  16   as our histograms have.  They in general go to 100 GeV/c in pT while we only
  17   go to 50. In that range the binning is EXACTLY the same (and the code checks
  18   it).
  19
  20   In charged spectra there is NO NORMALIZATION uncertainty. We add it by hand
  21   to the syst error.
  22
  23   * Need at some point to get the final spectra
  24   * Need at some point to get the final normalization uncertainty
  25
  26
  27   To run:
  28
  29   gSystem->AddIncludePath("-I../macros")
  30   gSystem->AddIncludePath("-I../lib")
  31   gROOT->SetMacroPath(".:../macros")
  32   .L my_tools.C+
  33   .L convertToMyBins.C+
  34
  35
  36   ConvertPPfile("charged_spectra/dNdPt_pp_7TeV.root",    "NOPID_pp_7TeV.root", 0.034);
  37   ConvertPPfile("charged_spectra/dNdPt_pp_2.76TeV.root", "NOPID_pp_2.76TeV.root", 0.034);
  38   ConvertPPfile("charged_spectra/dNdPt_pp_900GeV.root",  "NOPID_pp_900GeV.root", 0.034);
  39
  40
  41
  42
  43   Older methods that can be useful for AA:
  44   The method convertToMyBinsAA is for pp and AA spectra.
  45   The method convertToMyBinsRaa is for RAA.
  46   The method convertToMyBinsOld is for AA data when the centralities had to be merged.
  47
  48  */
  49
  50
  51
  52 TH1D* Convert(TH1D* histIn, Double_t addSyst=0);
  53 TH1D* ConvertGraph(TGraphErrors* graphIn);
  54 void convertToMyBinsAA();
  55 void convertToMyBinsRaa(const Char_t* fileName = "PbPb_RAA_2760GeV_200511.root");
  56 void convertToMyBinsOld(const Char_t* fileName = "dndpt_all_2011-05-15.root");
  57
  58
  59 void  ConvertPPfile(const Char_t* inFileName, const Char_t* outFileName, Double_t normSyst=0.0)
  60 {
  61   TFile* fileDataCharged = FindFileFresh(inFileName);
  62
  63   if(!fileDataCharged)
  64     return;
  65
  66   TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat");
  67   TH1D* hPtCharged = Convert(hPtHelp);
  68   hPtCharged->SetName("hDnDptCharged_PP");
  69   hPtCharged->SetMarkerStyle(29);
  70
  71   hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst");
  72   TH1D* hPtChargedSyst = Convert(hPtHelp, normSyst);
  73   hPtChargedSyst->SetName("hDnDptChargedSyst_PP");
  74   hPtChargedSyst->SetMarkerStyle(29);
  75
  76   TFile* fileOut = new TFile(outFileName, "RECREATE");
  77   hPtCharged->Write();
  78   hPtChargedSyst->Write();
  79   fileOut->Close();
  80 }
  81
  82 //________________________________________________________________________________________
  83 TH1D* Convert(TH1D* histIn, Double_t addSyst)
  84 {
  85   const Int_t nPtBins = 68;
  86   Double_t xBins[nPtBins+1] = {0. ,  0.05, 0.1,  0.15, 0.2,  0.25, 0.3,  0.35, 0.4,  0.45,
  87                                0.5,  0.55, 0.6,  0.65, 0.7,  0.75, 0.8,  0.85, 0.9,  0.95,
  88                                1.0,  1.1 , 1.2,  1.3 , 1.4,  1.5 , 1.6,  1.7 , 1.8,  1.9 ,
  89                                2.0,  2.2 , 2.4,  2.6 , 2.8,  3.0 , 3.2,  3.4 , 3.6,  3.8 ,
  90                                4.0,  4.5 , 5.0,  5.5 , 6.0,  6.5 , 7.0,  8.0 , 9.0,  10.0,
  91                                11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
  92                                26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
  93
  94
  95   TH1D* hist = new TH1D("hist", "", nPtBins, xBins);
  96   hist->SetTitle(histIn->GetTitle());
  97   hist->GetXaxis()->SetTitle(histIn->GetXaxis()->GetTitle());
  98   hist->GetYaxis()->SetTitle(histIn->GetYaxis()->GetTitle());
  99
 100   const Double_t deltapt = 0.0001;
 101
 102   for(Int_t bin = 1; bin <= nPtBins; bin++) {
 103
 104     // check bin size
 105     if(TMath::Abs(hist->GetXaxis()->GetBinLowEdge(bin) -
 106                   histIn->GetXaxis()->GetBinLowEdge(bin)) > deltapt) {
 107       cout << "pt edge low does not agree!!!!!!!!!!!!!!!" << endl;
 108       return 0;
 109     }
 110     if(TMath::Abs(hist->GetXaxis()->GetBinUpEdge(bin) -
 111                   histIn->GetXaxis()->GetBinUpEdge(bin)) > deltapt) {
 112       cout << "pt edge high does not agree!!!!!!!!!!!!!!!" << endl;
 113       return 0;
 114     }
 115     if(TMath::Abs(hist->GetXaxis()->GetBinCenter(bin) -
 116                   histIn->GetXaxis()->GetBinCenter(bin)) > deltapt) {
 117       cout << "pt center does not agree!!!!!!!!!!!!!!!" << endl;
 118       return 0;
 119     }
 120
 121
 122     hist->SetBinContent(bin, histIn->GetBinContent(bin));
 123
 124     Double_t error = histIn->GetBinError(bin);
 125     if(addSyst) {
 126
 127       Double_t extra_error = addSyst*histIn->GetBinContent(bin);
 128       error = TMath::Sqrt(error*error + extra_error*extra_error);
 129     }
 130     hist->SetBinError(bin, error);
 131   }
 132
 133   return hist;
 134 }
 135
 136 //________________________________________________________________________________________
 137 TH1D* ConvertGraph(TGraphErrors* graphIn)
 138 {
 139   const Int_t nPtBins = 68;
 140   Double_t xBins[nPtBins+1] = {0. ,  0.05, 0.1,  0.15, 0.2,  0.25, 0.3,  0.35, 0.4,  0.45,
 141                                0.5,  0.55, 0.6,  0.65, 0.7,  0.75, 0.8,  0.85, 0.9,  0.95,
 142                                1.0,  1.1 , 1.2,  1.3 , 1.4,  1.5 , 1.6,  1.7 , 1.8,  1.9 ,
 143                                2.0,  2.2 , 2.4,  2.6 , 2.8,  3.0 , 3.2,  3.4 , 3.6,  3.8 ,
 144                                4.0,  4.5 , 5.0,  5.5 , 6.0,  6.5 , 7.0,  8.0 , 9.0,  10.0,
 145                                11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
 146                                26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
 147
 148
 149   TH1D* hist = new TH1D("hist", "", nPtBins, xBins);
 150   // hist->SetTitle(graphIn->GetTitle());
 151   // hist->GetXaxis()->SetTitle(graphIn->GetXaxis()->GetTitle());
 152   // hist->GetYaxis()->SetTitle(graphIn->GetYaxis()->GetTitle());
 153
 154   const Double_t deltapt = 0.0001;
 155
 156   // Jacek does not fill the first 3 bins of his graph?
 157   for(Int_t bin = 4; bin <= nPtBins; bin++) {
 158
 159     // check bin size
 160     if(TMath::Abs(hist->GetXaxis()->GetBinCenter(bin) -
 161                   graphIn->GetX()[bin-1]) > deltapt) {
 162       cout << "pt center does not agree!!!!!!!!!!!!!!!" << endl;
 163       return 0;
 164     }
 165
 166
 167     hist->SetBinContent(bin, graphIn->GetY()[bin-1]);
 168     hist->SetBinError(bin, graphIn->GetEY()[bin-1]);
 169   }
 170
 171   return hist;
 172 }
 173
 174 //________________________________________________________________________________________
 175 void convertToMyBinsAA()
 176 {
 177   const Int_t nCent = 7;
 178   const Int_t cent[nCent+1] = {-5, 0, 5, 10, 20, 40, 60, 80};
 179   const Int_t colors[nCent] = {kBlack, kRed+1, kOrange+7, kOrange, kGreen+2, kCyan+2, kBlue+1};
 180
 181   TH1D* hPtCharged[nCent] = {0, 0, 0, 0, 0, 0, 0};
 182   TH1D* hPtChargedSyst[nCent] = {0, 0, 0, 0, 0, 0, 0};
 183
 184   for (Int_t i = 0; i < nCent; i++) {
 185
 186     if(cent[i]==-5) {
 187
 188       TFile* fileDataCharged = FindFileFresh("dNdPt_pp_2.76TeV.root");
 189
 190       if(!fileDataCharged)
 191         return;
 192
 193       TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat");
 194       hPtCharged[i] = Convert(hPtHelp);
 195       hPtCharged[i]->SetName("hDnDptCharged_PP");
 196       hPtCharged[i]->SetMarkerColor(colors[i]);
 197       hPtCharged[i]->SetMarkerStyle(29);
 198       hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst");
 199       hPtChargedSyst[i] = Convert(hPtHelp, 0.034); // 3.4 % norm error
 200       /*
 201         Dear Peter,
 202
 203         for 2.76TeV the systematic error for MB->INEL is 3.40%,, we decided to
 204         drop the uncertainty from VdM scan since we don't really use it...
 205
 206         Michael
 207       */
 208       hPtChargedSyst[i]->SetName("hDnDptChargedSyst_PP");
 209       hPtChargedSyst[i]->SetMarkerColor(colors[i]);
 210       hPtChargedSyst[i]->SetMarkerStyle(29);
 211     } else {
 212
 213       TFile* fileDataCharged = FindFileFresh(Form("dNdPt_PbPb_2.76ATeV_centrality_%d-%d.root", cent[i], cent[i+1]));
 214
 215       if(!fileDataCharged)
 216         return;
 217
 218       TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat");
 219       hPtCharged[i] = Convert(hPtHelp);
 220       hPtCharged[i]->SetName(Form("hDnDptCharged_%d_%d", cent[i], cent[i+1]));
 221       hPtCharged[i]->SetMarkerColor(colors[i]);
 222       hPtCharged[i]->SetMarkerStyle(29);
 223       hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst");
 224       hPtChargedSyst[i] = Convert(hPtHelp);
 225       hPtChargedSyst[i]->SetName(Form("hDnDptChargedSyst_%d_%d", cent[i], cent[i+1]));
 226       hPtChargedSyst[i]->SetMarkerColor(colors[i]);
 227       hPtChargedSyst[i]->SetMarkerStyle(29);
 228     }
 229   }
 230
 231   TFile* fileOut = new TFile("dNdPt_Charged.root", "RECREATE");
 232   for (Int_t i = 0; i < nCent; i++) {
 233     hPtCharged[i]->Write();
 234     hPtChargedSyst[i]->Write();
 235   }
 236   fileOut->Close();
 237 }
 238
 239 //________________________________________________________________________________________
 240 void convertToMyBinsRaa(const Char_t* fileName)
 241 {
 242   TFile* fileData = FindFileFresh(fileName);
 243
 244   if(!fileData)
 245     return;
 246
 247   TGraphErrors* hRaa_0_5 = (TGraphErrors*)fileData->Get("raa_c0_5_stat");
 248   TH1D* hRaaCharged_0_5 = ConvertGraph(hRaa_0_5);
 249   hRaaCharged_0_5->SetName("hRaaCharged_0_5");
 250
 251   TGraphErrors* hRaa_5_10 = (TGraphErrors*)fileData->Get("raa_c5_10_stat");
 252   TH1D* hRaaCharged_5_10 = ConvertGraph(hRaa_5_10);
 253   hRaaCharged_5_10->SetName("hRaaCharged_5_10");
 254
 255   TGraphErrors* hRaa_10_20 = (TGraphErrors*)fileData->Get("raa_c10_20_stat");
 256   TH1D* hRaaCharged_10_20 = ConvertGraph(hRaa_10_20);
 257   hRaaCharged_10_20->SetName("hRaaCharged_10_20");
 258
 259   TGraphErrors* hRaa_20_40 = (TGraphErrors*)fileData->Get("raa_c20_40_stat");
 260   TH1D* hRaaCharged_20_40 = ConvertGraph(hRaa_20_40);
 261   hRaaCharged_20_40->SetName("hRaaCharged_20_40");
 262
 263   TGraphErrors* hRaa_40_60 = (TGraphErrors*)fileData->Get("raa_c40_60_stat");
 264   TH1D* hRaaCharged_40_60 = ConvertGraph(hRaa_40_60);
 265   hRaaCharged_40_60->SetName("hRaaCharged_40_60");
 266
 267   TGraphErrors* hRaa_60_80 = (TGraphErrors*)fileData->Get("raa_c60_80_stat");
 268   TH1D* hRaaCharged_60_80 = ConvertGraph(hRaa_60_80);
 269   hRaaCharged_60_80->SetName("hRaaCharged_60_80");
 270
 271   TGraphErrors* hRaa_0_20 = (TGraphErrors*)fileData->Get("raa_c0_20_stat");
 272   TH1D* hRaaCharged_0_20 = ConvertGraph(hRaa_0_20);
 273   hRaaCharged_0_20->SetName("hRaaCharged_0_20");
 274
 275   TGraphErrors* hRaa_40_80 = (TGraphErrors*)fileData->Get("raa_c40_80_stat");
 276   TH1D* hRaaCharged_40_80 = ConvertGraph(hRaa_40_80);
 277   hRaaCharged_40_80->SetName("hRaaCharged_40_80");
 278
 279
 280   TFile* fileOut = new TFile("Raa_Charged.root", "RECREATE");
 281   hRaaCharged_0_5->Write();
 282   hRaaCharged_5_10->Write();
 283   hRaaCharged_10_20->Write();
 284   hRaaCharged_20_40->Write();
 285   hRaaCharged_40_60->Write();
 286   hRaaCharged_60_80->Write();
 287   hRaaCharged_0_20->Write();
 288   hRaaCharged_40_80->Write();
 289   fileOut->Close();
 290 }
 291
 292 //________________________________________________________________________________________
 293
 294 void convertToMyBinsOld(const Char_t* fileName)
 295 {
 296   TFile* fileData = FindFileFresh(fileName);
 297
 298   if(!fileData)
 299     return;
 300
 301   TH1D* hDnDPt_0_5 = (TH1D*)fileData->Get("dNdPt_PbPb_c0_5");
 302   TH1D* hDnDptCharged_0_5 = Convert(hDnDPt_0_5);
 303   hDnDptCharged_0_5->SetName("hDnDptCharged_0_5");
 304
 305   TH1D* hDnDPt_5_10 = (TH1D*)fileData->Get("dNdPt_PbPb_c5_10");
 306   TH1D* hDnDptCharged_5_10 = Convert(hDnDPt_5_10);
 307   hDnDptCharged_5_10->SetName("hDnDptCharged_5_10");
 308
 309   TH1D* hDnDPt_10_20 = (TH1D*)fileData->Get("dNdPt_PbPb_c10_20");
 310   TH1D* hDnDptCharged_10_20 = Convert(hDnDPt_10_20);
 311   hDnDptCharged_10_20->SetName("hDnDptCharged_10_20");
 312
 313   TH1D* hDnDPt_20_40 = (TH1D*)fileData->Get("dNdPt_PbPb_c20_30");
 314   TH1D* hHelp =  (TH1D*)fileData->Get("dNdPt_PbPb_c30_40");
 315   hDnDPt_20_40->Add(hHelp);
 316   hDnDPt_20_40->Scale(1.0/2.0);
 317   TH1D* hDnDptCharged_20_40 = Convert(hDnDPt_20_40);
 318   hDnDptCharged_20_40->SetName("hDnDptCharged_20_40");
 319
 320   TH1D* hDnDPt_40_60 = (TH1D*)fileData->Get("dNdPt_PbPb_c40_50");
 321   hHelp =  (TH1D*)fileData->Get("dNdPt_PbPb_c50_60");
 322   hDnDPt_40_60->Add(hHelp);
 323   hDnDPt_40_60->Scale(1.0/2.0);
 324   TH1D* hDnDptCharged_40_60 = Convert(hDnDPt_40_60);
 325   hDnDptCharged_40_60->SetName("hDnDptCharged_40_60");
 326
 327   TH1D* hDnDPt_60_80 = (TH1D*)fileData->Get("dNdPt_PbPb_c60_70");
 328   hHelp =  (TH1D*)fileData->Get("dNdPt_PbPb_c70_80");
 329   hDnDPt_60_80->Add(hHelp);
 330   hDnDPt_60_80->Scale(1.0/2.0);
 331   TH1D* hDnDptCharged_60_80 = Convert(hDnDPt_60_80);
 332   hDnDptCharged_60_80->SetName("hDnDptCharged_60_80");
 333
 334
 335   TFile* fileOut = new TFile("dNdPt_Charged.root", "RECREATE");
 336   hDnDptCharged_0_5->Write();
 337   hDnDptCharged_5_10->Write();
 338   hDnDptCharged_10_20->Write();
 339   hDnDptCharged_20_40->Write();
 340   hDnDptCharged_40_60->Write();
 341   hDnDptCharged_60_80->Write();
 342   fileOut->Close();
 343 }