cascadeAnalysis/systematics/CalculateSystematic.C

0001 void CalculateSystematic()
0002 {
0003   //TFile* inputFile = TFile::Open("cutVarianceSystematics.root", "READ");
0004   //TFile* inputFile = TFile::Open("geoAccEffSystematics.root", "READ");
0005   //TFile* inputFile = TFile::Open("trackAccSystematics.root", "READ");
0006   TFile* inputFile = TFile::Open("fitFunctionSystematics.root", "READ");
0007   //TFile* inputFile = TFile::Open("cutEffSystematics.root", "READ");
0008
0009   const std::vector<std::string> kinVars  = {"pT", "eta", "phi", "rap"};
0010   const std::vector<std::string> variants = {"base", "base_min", "base_max"};
0011   //const std::vector<std::string> variants = {"base", "base_max"};
0012
0013   // Load all histograms and detach from file
0014   std::map<std::string, TH1F*> hMap;
0015   for (const auto& var : variants)
0016   {
0017     for (const auto& kin : kinVars)
0018     {
0019       for (const auto& par : {"KS0", "Xi"})
0020       {
0021         std::string suffix  = (var == "base") ? "" : var.substr(4); // "" / "_min" / "_max"
0022         std::string name    = "h_" + kin + "_base_" + par + suffix;
0023         TH1F* h             = (TH1F*)inputFile->Get(name.c_str());
0024         h->SetDirectory(0);
0025         hMap[var + "_" + kin + "_" + par] = h;
0026       }
0027     }
0028   }
0029   inputFile->Close();
0030
0031   // Compute Xi/(2*KS0) ratios for each variant and kinematic variable
0032   std::map<std::string, TH1F*> ratioMap;
0033   for (const auto& var : variants)
0034   {
0035     for (const auto& kin : kinVars)
0036     {
0037       TH1F* hXi  = hMap.at(var + "_" + kin + "_Xi");
0038       TH1F* hKS0 = hMap.at(var + "_" + kin + "_KS0");
0039
0040       TH1F* hRatio = (TH1F*)hXi->Clone(("ratio_" + var + "_" + kin).c_str());
0041       hRatio->SetDirectory(0);
0042       hRatio->Divide(hXi, hKS0, 1, 2); // Xi / (2*KS0) via scale factor
0043       ratioMap[var + "_" + kin] = hRatio;
0044     }
0045   }
0046
0047   // Calculate systematic uncertainty as max bin-by-bin deviation of min/max from base
0048   //TFile* outputFile = TFile::Open("cutVarianceSystematicUncertainties.root", "RECREATE");
0049   //TFile* outputFile = TFile::Open("geoAccEffSystematicUncertainties.root", "RECREATE");
0050   //TFile* outputFile = TFile::Open("trackAccSystematicUncertainties.root", "RECREATE");
0051   //TFile* outputFile = TFile::Open("fitFunctionSystematicUncertainties.root", "RECREATE");
0052   //TFile* outputFile = TFile::Open("cutEffSystematicUncertainties.root", "RECREATE");
0053
0054   for (const auto& kin : kinVars)
0055   {
0056     TH1F* hBase = ratioMap.at("base_"     + kin);
0057     TH1F* hMin  = ratioMap.at("base_min_" + kin);
0058     TH1F* hMax  = ratioMap.at("base_max_" + kin);
0059
0060     TH1F* hSyst = (TH1F*)hBase->Clone(("h_syst_" + kin).c_str());
0061     hSyst->SetDirectory(0);
0062     hSyst->Reset();
0063
0064     for (int b = 1; b <= hBase->GetNbinsX(); b++)
0065     {
0066       double base    = hBase->GetBinContent(b);
0067       double baseErr = hBase->GetBinError(b);
0068       double diffMin = std::abs(hMin->GetBinContent(b) - base);
0069       double diffMax = std::abs(hMax->GetBinContent(b) - base);
0070       double diffErr = diffMax > diffMin ? hMax->GetBinError(b) : hMin->GetBinError(b);
0071       //double diffErr = hMax->GetBinError(b);
0072       double diff = std::max(diffMin, diffMax);
0073       //double diff = diffMax;
0074       //hSyst->SetBinContent(b, std::max(diffMin, diffMax));
0075       double barlowSys = diff;
0076       //if (diff*diff > (baseErr*baseErr + diffErr*diffErr))
0077       //{
0078       //  barlowSys = std::sqrt(diff*diff - (baseErr*baseErr + diffErr*diffErr));
0079       //}
0080       //else
0081       //{
0082       //  barlowSys = 0.0;
0083       //}
0084       //std::cout << "diff: " << diff << std::endl;
0085       //std::cout << "baseErr: " << baseErr << std::endl;
0086       //std::cout << "diffErr: " << diffErr << std::endl;
0087       //std::cout << "barlowSys: " << barlowSys << std::endl;
0088       hSyst->SetBinContent(b, barlowSys);
0089       std::cout << "Kin: " << kin << ", bin: " << b << ", barlowSys: " << 100*barlowSys/base << std::endl;
0090     }
0091
0092     //outputFile->cd();
0093     //hBase->Write(("h_ratio_base_" + kin).c_str());
0094     //hSyst->Write(("h_syst_"       + kin).c_str());
0095   }
0096
0097   //outputFile->Close();
0098   //std::cout << "Systematics saved to cutVarianceSystematicUncertainties.root" << std::endl;
0099   //std::cout << "Systematics saved to geoAccEffSystematicUncertainties.root" << std::endl;
0100   //std::cout << "Systematics saved to trackAccSystematicUncertainties.root" << std::endl;
0101   //std::cout << "Systematics saved to fitFunctionSystematicUncertainties.root" << std::endl;
0102   //std::cout << "Systematics saved to cutEffSystematicUncertainties.root" << std::endl;
0103 }