// Macro to make the required plots with full error information // Desire main function, and two plotting functions which do different // things depending on whether or not want plots for SF or Eff // Ian Connelly 29-12-12 #include #include #include #include #include "TFile.h" #include "TH1D.h" #include "TCanvas.h" #include "TGraphAsymmErrors.h" #include "TString.h" #include "AtlasUtils.C" using namespace std; // Global objects vector vals; vector statUp; vector statDo; vector systUp; vector systDo; TString yaxis; // Function prototypes void getVals(std::string); // Extract values from file int errorCheck(void); // Check vectors are filled void plotEff(void); // Function to plot efficiencies void plotSF(void); // Function to plot scale factor double* createArray(const std::vector&); // Function to create array from vector for TGraphAsymmError void setAbsErrors(vector&, const vector&); // Set error values as absolute (not relative to nominal) void makeTotalError(vector&, const vector&); // Make the total error combining stat and syst // Main function int makePlots(){ cout << "==================================================================" << endl; cout << " Please be aware, currently only uses text files indiscriminantly " << endl; cout << " and does not know if EM or LC. Only option here is to use the " << endl; cout << " correct file when extracting the MC eff information from hists. " << endl; cout << "==================================================================" <> type; if(type != "Eff" && type != "SF"){cerr << "Incorrect type - Use Eff/SF" << endl; return 101;} // Load AtlasStyle gSystem->Load("../AtlasStyle_C.so"); AtlasStyle(); gSystem->Load("libAtlasUtils_C.so"); if(type == "Eff"){yaxis = "Efficiencies"; plotEff();} if(type == "SF") {yaxis = "Scale Factor"; plotSF();} return 0; } void getVals(std::string filepath){ // Clear all vectors vals.clear(); statUp.clear(); statDo.clear(); systUp.clear(); systDo.clear(); // Try to open file ifstream input(filepath.c_str()); std::string line; if(input.is_open()){} else{cerr << "Error in opening file : " << filepath << endl; return;} // Start reading lines while(std::getline(input,line)){ // Initialise array to be filled and set to zero each line float a_vals[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; // Read line std::string systName; TString syst; // Get string before first ampersand (inc. whitespace -> IMPORTANT) sscanf(line.c_str(),"%[^&]", systName.c_str()); syst = systName.c_str(); cout << syst << endl; // Nominal Vals if(syst.EqualTo("Nominal ")){ sscanf(line.c_str(),"%*[^&] & %*f & %*f & %f & %f & %f & %f & %f & %f & %f & %f & %f & %f", &a_vals[0], &a_vals[1], &a_vals[2], &a_vals[3], &a_vals[4], &a_vals[5], &a_vals[6], &a_vals[7], &a_vals[8],&a_vals[9]); vals.push_back(a_vals[0]); vals.push_back(a_vals[1]); vals.push_back(a_vals[2]); vals.push_back(a_vals[3]); vals.push_back(a_vals[4]); vals.push_back(a_vals[5]); vals.push_back(a_vals[6]); vals.push_back(a_vals[7]); vals.push_back(a_vals[8]); vals.push_back(a_vals[9]); } // Statistical error up else if(syst.EqualTo("Total stat (+) ")){ sscanf(line.c_str(),"%*[^&] & %*f & %*f & %f & %f & %f & %f & %f & %f & %f & %f & %f & %f", &a_vals[0], &a_vals[1], &a_vals[2], &a_vals[3], &a_vals[4], &a_vals[5], &a_vals[6], &a_vals[7], &a_vals[8],&a_vals[9]); statUp.push_back(a_vals[0]); statUp.push_back(a_vals[1]); statUp.push_back(a_vals[2]); statUp.push_back(a_vals[3]); statUp.push_back(a_vals[4]); statUp.push_back(a_vals[5]); statUp.push_back(a_vals[6]); statUp.push_back(a_vals[7]); statUp.push_back(a_vals[8]); statUp.push_back(a_vals[9]); } // Statisitical error down else if(syst.EqualTo("Total stat (-) ")){ sscanf(line.c_str(),"%*[^&] & %*f & %*f & %f & %f & %f & %f & %f & %f & %f & %f & %f & %f", &a_vals[0], &a_vals[1], &a_vals[2], &a_vals[3], &a_vals[4], &a_vals[5], &a_vals[6], &a_vals[7], &a_vals[8],&a_vals[9]); statDo.push_back(a_vals[0]); statDo.push_back(a_vals[1]); statDo.push_back(a_vals[2]); statDo.push_back(a_vals[3]); statDo.push_back(a_vals[4]); statDo.push_back(a_vals[5]); statDo.push_back(a_vals[6]); statDo.push_back(a_vals[7]); statDo.push_back(a_vals[8]); statDo.push_back(a_vals[9]); } // Systematics error up else if(syst.EqualTo("Total syst (+) ")){ sscanf(line.c_str(),"%*[^&] & %*f & %*f & %f & %f & %f & %f & %f & %f & %f & %f & %f & %f", &a_vals[0], &a_vals[1], &a_vals[2], &a_vals[3], &a_vals[4], &a_vals[5], &a_vals[6], &a_vals[7], &a_vals[8],&a_vals[9]); systUp.push_back(a_vals[0]); systUp.push_back(a_vals[1]); systUp.push_back(a_vals[2]); systUp.push_back(a_vals[3]); systUp.push_back(a_vals[4]); systUp.push_back(a_vals[5]); systUp.push_back(a_vals[6]); systUp.push_back(a_vals[7]); systUp.push_back(a_vals[8]); systUp.push_back(a_vals[9]); } // Systematics error down else if(syst.EqualTo("Total syst (-) ")){ sscanf(line.c_str(),"%*[^&] & %*f & %*f & %f & %f & %f & %f & %f & %f & %f & %f & %f & %f", &a_vals[0], &a_vals[1], &a_vals[2], &a_vals[3], &a_vals[4], &a_vals[5], &a_vals[6], &a_vals[7], &a_vals[8], &a_vals[9]); systDo.push_back(a_vals[0]); systDo.push_back(a_vals[1]); systDo.push_back(a_vals[2]); systDo.push_back(a_vals[3]); systDo.push_back(a_vals[4]); systDo.push_back(a_vals[5]); systDo.push_back(a_vals[6]); systDo.push_back(a_vals[7]); systDo.push_back(a_vals[8]); systDo.push_back(a_vals[9]); } } // Vectors should be filled! } int errorCheck(){ // Return 0 if all fine // Return !=0 if a vector is not filled if(vals.size() == 0){return 1;} else if(statUp.size() == 0){return 2;} else if(statDo.size() == 0){return 3;} else if(systUp.size() == 0){return 4;} else if(systDo.size() == 0){return 5;} else return 0; } double* createArray( const std::vector< double >& v ){ double* result = new double [v.size()]; memcpy( result, &v.front(), v.size() * sizeof( double ) ); return result; } void setAbsErrors(vector& relError, const vector& nominal){ for(int i = 0; i < relError.size(); i++){ relError.at(i) = relError.at(i)*nominal.at(i)/100; // multiply the central value by the percentage value relError.at(i) = fabs(relError.at(i)); // TGraphAsymError needs absolute values (so ensure for statistical) } return; } void makeTotalError(vector& syst, const vector& stat){ // Use this function before the absolute errors are defined for(int i = 0; i < syst.size(); i++){ syst.at(i) = sqrt(syst.at(i)*syst.at(i) + stat.at(i)*stat.at(i)); } return; } void plotEff(){ gSystem->Load("../AtlasStyle_C.so"); SetAtlasStyle(); gSystem->Load("libAtlasUtils_C.so"); // List of filenames int nFiles = 17; std::string relUncert[] = {"relUncertEff_mv1_60","relUncertEff_mv1_70","relUncertEff_mv1_75","relUncertEff_mv1_80","relUncertEff_sv0_50","relUncertEff_mv1c_50","relUncertEff_mv1c_57","relUncertEff_wcomb_1","relUncertEff_wcomb_2","relUncertEff_wcomb_3","relUncertEff_mv1_30","relUncertEff_mv1_50","relUncertEff_mv1_90","relUncertEff_mv1c_30","relUncertEff_mv1c_70","relUncertEff_mv1c_80","relUncertEff_mv1c_90"}; TString histNames[] = {"MV1_60","MV1_70","MV1_75","MV1_80","SV0_50","MV1c_50","MV1c_57","WCOMB_1","WCOMB_2","WCOMB_3","MV1_30","MV1_50","MV1_90","MV1c_30","MV1c_70","MV1c_80","MV1c_90"}; double effLine[] = {0.6,0.7,0.75,0.8,0.5,0.5,0.57,0.5,0.5,0.5,0.3,0.5,0.9,0.3,0.7,0.8,0.9}; std::string path = "../"; cout << "Is EM or LC : "; TString jtype; cin >> jtype; cout << "Is eµ type (2j,3j) : "; TString emType; cin >> emType; cout << "Is JVF 0.5 ? (y/n) : "; TString jvf; cin >> jvf; if(jvf == "y"){ jtype += "_JVF0.5"; } for(int i = 0; i < nFiles; i++){ // Call getVals - This fills the vectors with values from the text file getVals(path+relUncert[i]+"_em"+emType.Data()+".txt"); // Error check if(errorCheck() != 0){cerr << "Get next file -> Error in vectors being filled : " << errorCheck() << endl; continue;} // Now can use the values and error vectors // Open the MC tt efficiency // TFile* f_ttMC = new TFile("../tt_DL_"+jtype+"_forMCEff.root","OPEN"); TFile* f_ttMC = new TFile("../clusterOutput/results_ttPowheg_UNKNOWN_EM"+emType+"_"+jtype+".root","OPEN"); // ttbar MC efficiency is true b passing the cut divided by the total number of true b TH1D* h_mctotal = (TH1D*)f_ttMC->Get("total_"+histNames[i]+"_b"); TH1D* h_mcpass = (TH1D*)f_ttMC->Get("jets12_pt_"+histNames[i]+"_b"); // Get arrays int nBins = 10; // double arr_ptbins[] = {20000,30000,40000,50000,60000,75000,90000,110000,140000,200000,300000}; double arr_ptmids[] = {25000,35000,45000,55000,67500,82500,100000,125000,170000,250000}; double arr_ptwide[] = {5000,5000,5000,5000,7500,7500,10000,15000,30000,50000}; // Make the total error and store in stat vector makeTotalError(systUp,statUp); makeTotalError(systDo,statDo); // Make sure the statistical errors contain the absolute values setAbsErrors(statUp,vals); setAbsErrors(statDo,vals); // Make sure the total error (stored in syst vectors) contains the absolute values setAbsErrors(systUp,vals); setAbsErrors(systDo,vals); // Get arrays from the vectors double* arr_vals = createArray(vals); double* arr_statUp = createArray(statUp); double* arr_statDo = createArray(statDo); double* arr_systUp = createArray(systUp); double* arr_systDo = createArray(systDo); //cout << arr_vals[3] << " " << vals.at(3) << endl; // Canvas TCanvas* c1 = new TCanvas("canvas","canvas",600,600); c1->cd(); // TGraphAsymmErrors Objects TGraphAsymmErrors tg_MC (h_mcpass,h_mctotal); TGraphAsymmErrors tg_data_stat (nBins, arr_ptmids, arr_vals, arr_ptwide, arr_ptwide, arr_statDo, arr_statUp); TGraphAsymmErrors tg_data_total (nBins, arr_ptmids, arr_vals, arr_ptwide, arr_ptwide, arr_systDo, arr_systUp); tg_MC.SetTitle(histNames[i]); tg_data_stat.SetTitle(histNames[i]); tg_data_total.SetTitle(histNames[i]); tg_MC.SetLineColor(2); tg_data_stat.SetLineColor(1); // tg_data_stat.SetFillColor(2); tg_data_stat.SetMarkerStyle(20); tg_data_stat.SetMarkerSize(1.2); tg_data_stat.SetLineWidth(2.); tg_data_total.SetLineColor(4); tg_data_total.SetFillColor(30); tg_data_total.SetMarkerStyle(20); tg_data_total.SetMarkerSize(1.2); tg_data_total.SetLineWidth(2.); tg_data_total.GetXaxis()->SetTitle("p_{T} [MeV]"); tg_data_total.GetYaxis()->SetTitle(yaxis); tg_data_total.GetYaxis()->SetTitleOffset(1.4); TLine* midline = new TLine(0,effLine[i],300000,effLine[i]); midline->SetLineColor(1); midline->SetLineWidth(1); midline->SetLineStyle(2); cout << midline->GetX1() << " " << midline->GetX2() << " " << midline->GetY1() << " " << midline->GetY2() << endl; // Draw plots total, line, syst, mc tg_data_total.SetMaximum(1.2); tg_data_total.SetMinimum(0); tg_data_total.Draw("2ap"); midline->Draw("same"); tg_data_stat.Draw("Zpsame"); tg_MC.Draw("Zpsame"); // Set labels // c1->SetTitle(histNames[i]); ATLAS_LABEL(0.2,0.89); myText(0.2,0.84,1,"Work in progress"); myText (0.2,0.79,1,"#sqrt{s} = 8 TeV"); myText (0.2,0.74,1,histNames[i]); c1->SetRightMargin(0.1); c1->Print("effs_"+histNames[i]+emType+".eps"); cout << "MC" << endl; tg_MC.Print(); cout <<"Stat" << endl; tg_data_stat.Print(); } return; } void plotSF(){ gSystem->Load("../AtlasStyle_C.so"); SetAtlasStyle(); gSystem->Load("libAtlasUtils_C.so"); // List of filenames double sfLine = 1.0; int nFiles = 17; std::string relUncert[] = {"relUncertSF_mv1_60","relUncertSF_mv1_70","relUncertSF_mv1_75","relUncertSF_mv1_80","relUncertSF_sv0_50","relUncertSF_mv1c_50","relUncertSF_mv1c_57","relUncertSF_wcomb_1","relUncertSF_wcomb_2","relUncertSF_wcomb_3","relUncertSF_mv1_30","relUncertSF_mv1_50","relUncertSF_mv1_90","relUncertSF_mv1c_30","relUncertSF_mv1c_70","relUncertSF_mv1c_80","relUncertSF_mv1c_90"}; TString histNames[] = {"MV1_60","MV1_70","MV1_75","MV1_80","SV0_50","MV1c_50","MV1c_57","WCOMB_1","WCOMB_2","WCOMB_3","MV1_30","MV1_50","MV1_90","MV1c_30","MV1c_70","MV1c_80","MV1c_90"}; std::string path = "../"; cout << "Enter eµ type (2j,3j) : " ; TString emType; cin >> emType; for(int i = 0; i < nFiles; i++){ // Call getVals - This fills the vectors with values from the text file getVals(path+relUncert[i]+"_em"+emType.Data()+".txt"); // Error check if(errorCheck() != 0){cerr << "Get next file -> Error in vectors being filled : " << errorCheck() << endl; continue;} // Now can use the values and error vectors // For SF just use values from the file // Get arrays int nBins = 10; // double arr_ptbins[] = {20000,30000,40000,50000,60000,75000,90000,110000,140000,200000,300000}; double arr_ptmids[] = {25000,35000,45000,55000,67500,82500,100000,125000,170000,250000}; double arr_ptwide[] = {5000,5000,5000,5000,7500,7500,10000,15000,30000,50000}; // Make the total error and store in stat vector makeTotalError(systUp,statUp); makeTotalError(systDo,statDo); // Make sure the statistical errors contain the absolute values setAbsErrors(statUp,vals); setAbsErrors(statDo,vals); // Make sure the total error (stored in syst vectors) contains the absolute values setAbsErrors(systUp,vals); setAbsErrors(systDo,vals); // Get arrays from the vectors double* arr_vals = createArray(vals); double* arr_statUp = createArray(statUp); double* arr_statDo = createArray(statDo); double* arr_systUp = createArray(systUp); double* arr_systDo = createArray(systDo); // Canvas TCanvas* c1 = new TCanvas("canvas","canvas",600,600); c1->cd(); // TGraphAsymmErrors Objects TGraphAsymmErrors tg_data_stat (nBins, arr_ptmids, arr_vals, arr_ptwide, arr_ptwide, arr_statDo, arr_statUp); TGraphAsymmErrors tg_data_total (nBins, arr_ptmids, arr_vals, arr_ptwide, arr_ptwide, arr_systDo, arr_systUp); tg_data_stat.SetTitle(histNames[i]); tg_data_total.SetTitle(histNames[i]); tg_data_stat.SetLineColor(1); tg_data_stat.SetMarkerStyle(20); tg_data_stat.SetMarkerSize(1.2); tg_data_stat.SetLineWidth(2.); tg_data_total.SetLineColor(4); tg_data_total.SetFillColor(30); tg_data_total.SetMarkerStyle(20); tg_data_total.SetMarkerSize(1.2); tg_data_total.SetLineWidth(2.); tg_data_total.GetXaxis()->SetTitle("p_{T} [MeV]"); tg_data_total.GetYaxis()->SetTitle(yaxis); tg_data_total.GetYaxis()->SetTitleOffset(1.4); TLine* midline = new TLine(0,sfLine,300000,sfLine); midline->SetLineColor(1); midline->SetLineWidth(1); midline->SetLineStyle(2); cout << midline->GetX1() << " " << midline->GetX2() << " " << midline->GetY1() << " " << midline->GetY2() << endl; // Draw plots total, line, syst, mc tg_data_total.SetMaximum(1.5); tg_data_total.SetMinimum(0); tg_data_total.Draw("2ap"); midline->Draw("same"); tg_data_stat.Draw("Zpsame"); ATLAS_LABEL(0.2,0.89); myText(0.2,0.84,1,"Work in progress"); myText (0.2,0.79,1,"#sqrt{s} = 8 TeV"); myText (0.2,0.74,1,histNames[i]); c1->SetRightMargin(0.1); c1->Print("sf_"+histNames[i]+emType+".eps"); } return; }