PWG4/macros/electrons/makeSysErr.C

   1
   2
   3
   4 void makeSysErr() {
   5
   6   TFile* f0 = new TFile("data/scaled25Oct09/TOTALhistosscaled-LHC09b2-0.root");
   7   TFile* f1 = new TFile("data/scaled25Oct09/TOTALhistosscaled-LHC09b2-1.root");
   8   TFile* f2 = new TFile("data/scaled25Oct09/TOTALhistosscaled-LHC09b2-2.root");
   9   TFile* f3 = new TFile("data/scaled25Oct09/TOTALhistosscaled-LHC09b2-3.root");
  10
  11   double pyscale = (1.E6)*0.5*208*208*100/360; //seconds*lumi*Pb*Pb*acceptance
  12
  13   TH2F* tte0 = (TH2F*)f0->Get("AnaElectron_hPtNPEleTTEScaled");
  14   TH2F* emc0 = (TH2F*)f0->Get("AnaElectron_hPtNPEleEMCALScaled");
  15
  16   TH2F* tte1 = (TH2F*)f1->Get("AnaElectron_hPtNPEleTTEScaled");
  17   TH2F* emc1 = (TH2F*)f1->Get("AnaElectron_hPtNPEleEMCALScaled");
  18
  19   TH2F* tte2 = (TH2F*)f2->Get("AnaElectron_hPtNPEleTTEScaled");
  20   TH2F* emc2 = (TH2F*)f2->Get("AnaElectron_hPtNPEleEMCALScaled");
  21
  22   TH2F* tte3 = (TH2F*)f3->Get("AnaElectron_hPtNPEleTTEScaled");
  23   TH2F* emc3 = (TH2F*)f3->Get("AnaElectron_hPtNPEleEMCALScaled");
  24
  25   tte0->Scale(pyscale);
  26   tte1->Scale(pyscale);
  27   tte2->Scale(pyscale);
  28   tte3->Scale(pyscale);
  29   emc0->Scale(pyscale);
  30   emc1->Scale(pyscale);
  31   emc2->Scale(pyscale);
  32   emc3->Scale(pyscale);
  33
  34   alltte0 = (TH1F*)tte0->ProjectionX("alltte0",1,1);
  35   allemc0 = (TH1F*)emc0->ProjectionX("allemc0",1,1);
  36   alltte0->Rebin(5);
  37   allemc0->Rebin(5);
  38
  39   alltte1 = (TH1F*)tte1->ProjectionX("alltte1",1,1);
  40   allemc1 = (TH1F*)emc1->ProjectionX("allemc1",1,1);
  41   alltte1->Rebin(5);
  42   allemc1->Rebin(5);
  43
  44   alltte2 = (TH1F*)tte2->ProjectionX("alltte2",1,1);
  45   allemc2 = (TH1F*)emc2->ProjectionX("allemc2",1,1);
  46   alltte2->Rebin(5);
  47   allemc2->Rebin(5);
  48
  49   alltte3 = (TH1F*)tte3->ProjectionX("alltte3",1,1);
  50   allemc3 = (TH1F*)emc3->ProjectionX("allemc3",1,1);
  51   alltte3->Rebin(5);
  52   allemc3->Rebin(5);
  53
  54   Double_t tsum = 0.;
  55   Double_t esum = 0.;
  56   for(Int_t i = 1; i <= alltte0->GetNbinsX(); i++) {
  57     Double_t t0 = alltte0->GetBinContent(i);
  58     Double_t t1 = alltte1->GetBinContent(i);
  59     Double_t t2 = alltte2->GetBinContent(i);
  60     Double_t t3 = alltte3->GetBinContent(i);
  61     Double_t e0 = allemc0->GetBinContent(i);
  62     Double_t e1 = allemc1->GetBinContent(i);
  63     Double_t e2 = allemc2->GetBinContent(i);
  64     Double_t e3 = allemc3->GetBinContent(i);
  65
  66     Double_t td01 = (t0 - t1)/(t0+t1+t2+t3);
  67     Double_t td02 = (t0 - t2)/(t0+t1+t2+t3);
  68     Double_t td03 = (t0 - t3)/(t0+t1+t2+t3);
  69
  70     Double_t tave = (TMath::Abs(td01) + TMath::Abs(td02) + TMath::Abs(td03))/2.;
  71     tsum += tave;
  72
  73     Double_t ed01 = (e0 - e1)/(e0+e1+e2+e3);
  74     Double_t ed02 = (e0 - e2)/(e0+e1+e2+e3);
  75     Double_t ed03 = (e0 - e3)/(e0+e1+e2+e3);
  76
  77     Double_t eave = (TMath::Abs(ed01) + TMath::Abs(ed02) + TMath::Abs(ed03))/2.;
  78     printf("%d Average unc = %2.2f\n",i,eave);
  79     esum += eave;
  80   }
  81   Double_t tfinal = tsum/alltte0->GetNbinsX();
  82   Double_t efinal = esum/allemc0->GetNbinsX();
  83   printf("tfinal %f, efinal %f\n",tfinal,efinal);
  84
  85   TCanvas *c1 = new TCanvas();
  86   gPad->SetLogy();
  87   alltte0->SetLineWidth(2);
  88   alltte0->Draw();
  89
  90   TCanvas *c2 = new TCanvas();
  91   gPad->SetLogy();
  92   allemc0->SetLineWidth(2);
  93   allemc0->Draw();
  94
  95   TGraphErrors* terr = new TGraphErrors();
  96   terr->SetName("tteErr");
  97   TGraphErrors* eerr = new TGraphErrors();
  98   eerr->SetName("emcErr");
  99   for(Int_t i = 1; i <= alltte0->GetNbinsX(); i++) {
 100     terr->SetPoint(i-1,alltte0->GetBinCenter(i),alltte0->GetBinContent(i));
 101     terr->SetPointError(i-1,0.,(tfinal)*alltte0->GetBinContent(i));
 102
 103     eerr->SetPoint(i-1,allemc0->GetBinCenter(i),allemc0->GetBinContent(i));
 104     eerr->SetPointError(i-1,0.,(efinal)*allemc0->GetBinContent(i));
 105   }
 106   c1->cd();
 107   terr->SetFillColor(kRed-8);
 108   terr->Draw("3same");
 109   alltte0->Draw("same");
 110
 111   c2->cd();
 112   eerr->SetFillColor(kRed-8);
 113   eerr->Draw("3same");
 114   allemc0->Draw("same");
 115 }