// Original Author: Loic Quertenmont #include "Analysis_Global.h" #include "Analysis_CommonFunction.h" #include "Analysis_PlotFunction.h" #include "Analysis_Samples.h" #include "tdrstyle.C" #include "TGraphAsymmErrors.h" using namespace std; class stAllInfo{ public: double Mass, MassMean, MassSigma, MassCut; double XSec_Th, XSec_Err, XSec_Exp, XSec_ExpUp, XSec_ExpDown, XSec_Exp2Up, XSec_Exp2Down, XSec_Obs; double Eff, Eff_SYSTP, Eff_SYSTI, Eff_SYSTM, Eff_SYSTT, Eff_SYSTPU, TotalUnc; double EffE, EffE_SYSTP, EffE_SYSTI, EffE_SYSTM, EffE_SYSTT, EffE_SYSTPU; double Significance; double XSec_5Sigma; double Index, WP_Pt, WP_I, WP_TOF; float NData, NPred, NPredErr, NSign; double LInt; stAllInfo(string path=""){ //Default Values Mass = 0; MassMean = 0; MassSigma = 0; MassCut = 0; Index = 0; WP_Pt = 0; WP_I = 0; WP_TOF = 0; XSec_Th = 0; XSec_Err = 0; XSec_Exp = 1E50; XSec_ExpUp = 1E50; XSec_ExpDown = 1E50; XSec_Exp2Up = 1E50; XSec_Exp2Down = 1E50; XSec_Obs = 1E50; Significance = 1E50; XSec_5Sigma = 1E50; Eff = 0; Eff_SYSTP = 0; Eff_SYSTI = 0; Eff_SYSTM = 0; Eff_SYSTT = 0; Eff_SYSTPU = 0; EffE = 0; EffE_SYSTP = 0; EffE_SYSTI = 0; EffE_SYSTM = 0; EffE_SYSTT = 0; EffE_SYSTPU = 0; NData = 0; NPred = 0; NPredErr = 0; NSign = 0; LInt = 0; if(path=="")return; FILE* pFile = fopen(path.c_str(),"r"); if(!pFile){return; printf("Can't open %s\n",path.c_str()); return;} fscanf(pFile,"Mass : %lf\n",&Mass); fscanf(pFile,"MassMean : %lf\n",&MassMean); fscanf(pFile,"MassSigma : %lf\n",&MassSigma); fscanf(pFile,"MassCut : %lf\n",&MassCut); fscanf(pFile,"Index : %lf\n",&Index); fscanf(pFile,"WP_Pt : %lf\n",&WP_Pt); fscanf(pFile,"WP_I : %lf\n",&WP_I); fscanf(pFile,"WP_TOF : %lf\n",&WP_TOF); fscanf(pFile,"Eff : %lf +- %lf\n",&Eff , &EffE ); fscanf(pFile,"Eff_SystP : %lf +- %lf\n",&Eff_SYSTP , &EffE_SYSTP ); fscanf(pFile,"Eff_SystI : %lf +- %lf\n",&Eff_SYSTI , &EffE_SYSTI ); fscanf(pFile,"Eff_SystM : %lf +- %lf\n",&Eff_SYSTM , &EffE_SYSTM ); fscanf(pFile,"Eff_SystT : %lf +- %lf\n",&Eff_SYSTT , &EffE_SYSTT ); fscanf(pFile,"Eff_SystPU : %lf +- %lf\n",&Eff_SYSTPU, &EffE_SYSTPU); fscanf(pFile,"TotalUnc : %lf\n",&TotalUnc); fscanf(pFile,"Signif : %lf\n",&Significance); fscanf(pFile,"XSec_Th : %lf\n",&XSec_Th); fscanf(pFile,"XSec_Exp : %lf\n",&XSec_Exp); fscanf(pFile,"XSec_ExpUp : %lf\n",&XSec_ExpUp); fscanf(pFile,"XSec_ExpDown : %lf\n",&XSec_ExpDown); fscanf(pFile,"XSec_Exp2Up : %lf\n",&XSec_Exp2Up); fscanf(pFile,"XSec_Exp2Down: %lf\n",&XSec_Exp2Down); fscanf(pFile,"XSec_Obs : %lf\n",&XSec_Obs); fscanf(pFile,"NData : %E\n" ,&NData); fscanf(pFile,"NPred : %E\n" ,&NPred); fscanf(pFile,"NPredErr : %E\n" ,&NPredErr); fscanf(pFile,"NSign : %E\n" ,&NSign); fscanf(pFile,"LInt : %lf\n",&LInt); fscanf(pFile,"XSec_5Sigma : %lf\n",&XSec_5Sigma); fclose(pFile); } void Save(string path=""){ FILE* pFile = fopen(path.c_str(),"w"); if(!pFile)printf("Can't open file : %s\n",path.c_str()); fprintf(pFile,"Mass : %f\n",Mass); fprintf(pFile,"MassMean : %f\n",MassMean); fprintf(pFile,"MassSigma : %f\n",MassSigma); fprintf(pFile,"MassCut : %f\n",MassCut); fprintf(pFile,"Index : %f\n",Index); fprintf(pFile,"WP_Pt : %f\n",WP_Pt); fprintf(pFile,"WP_I : %f\n",WP_I); fprintf(pFile,"WP_TOF : %f\n",WP_TOF); fprintf(pFile,"Eff : %f +- %f\n",Eff , EffE); fprintf(pFile,"Eff_SystP : %f +- %f\n",Eff_SYSTP , EffE_SYSTP); fprintf(pFile,"Eff_SystI : %f +- %f\n",Eff_SYSTI , EffE_SYSTI); fprintf(pFile,"Eff_SystM : %f +- %f\n",Eff_SYSTM , EffE_SYSTM); fprintf(pFile,"Eff_SystT : %f +- %f\n",Eff_SYSTT , EffE_SYSTT); fprintf(pFile,"Eff_SystPU : %f +- %f\n",Eff_SYSTPU, EffE_SYSTPU); fprintf(pFile,"TotalUnc : %f\n",TotalUnc); fprintf(pFile,"Signif : %f\n",Significance); fprintf(pFile,"XSec_Th : %f\n",XSec_Th); fprintf(pFile,"XSec_Exp : %f\n",XSec_Exp); fprintf(pFile,"XSec_ExpUp : %f\n",XSec_ExpUp); fprintf(pFile,"XSec_ExpDown : %f\n",XSec_ExpDown); fprintf(pFile,"XSec_Exp2Up : %f\n",XSec_Exp2Up); fprintf(pFile,"XSec_Exp2Down: %f\n",XSec_Exp2Down); fprintf(pFile,"XSec_Obs : %f\n",XSec_Obs); fprintf(pFile,"NData : %+6.2E\n",NData); fprintf(pFile,"NPred : %+6.2E\n",NPred); fprintf(pFile,"NPredErr : %+6.2E\n",NPredErr); fprintf(pFile,"NSign : %+6.2E\n",NSign); fprintf(pFile,"LInt : %f\n",LInt); fprintf(pFile,"XSec_5Sigma : %f\n",XSec_5Sigma); fclose(pFile); } }; string EXCLUSIONDIR = "EXCLUSION"; //Background prediction rescale and uncertainty double RescaleFactor = 1.0; double RescaleError = 0.2; //final Plot y-axis range double PlotMinScale = 0.0001; double PlotMaxScale = 700; //Easy flag to skip running time consuming Cls expected limits. True runs the limit, false does not bool FullExpLimit=true; void Optimize(string InputPattern, string Data, string signal, bool shape, bool cutFromFile); double GetSignalMeanHSCPPerEvent(string InputPattern, unsigned int CutIndex, double MinRange, double MaxRange); TGraph* MakePlot(FILE* pFile, FILE* talkFile, string InputPattern, string ModelName, int XSectionType, std::vector<stSample>& modelSamples, double& LInt); TGraph* CheckSignalUncertainty(FILE* pFile, FILE* talkFile, string InputPattern, string ModelName, std::vector<stSample>& modelSample); void DrawModelLimitWithBand(string InputPattern); void DrawRatioBands(string InputPattern); void printSummary(FILE* pFile, FILE* talkFile, string InputPattern, string ModelName, std::vector<stSample>& modelSamples); void printSummaryPaper(FILE* pFile, FILE* talkFile, string InputPattern, string ModelName, std::vector<stSample>& modelSamples); string toLatex(double value); void makeDataCard(string outpath, string rootPath, string ChannelName, string SignalName, double Obs, double Pred, double PredRelErr, double Sign, double SignalStat, double SignalUnc, bool Shape); void saveHistoForLimit(TH1* histo, string Name, string Id); void saveVariationHistoForLimit(TH1* histo, TH1* vardown, string Name, string variationName); void testShapeBasedAnalysis(string InputPattern, string signal); double computeSignificance(string datacard, bool expected, string& signal, string massStr, float Strength); bool runCombine(bool fastOptimization, bool getXsection, bool getSignificance, string& InputPattern, string& signal, unsigned int CutIndex, bool Shape, bool Temporary, stAllInfo& result, TH1* MassData, TH1* MassPred, TH1* MassSign, TH1* MassSignP, TH1* MassSignI, TH1* MassSignM, TH1* MassSignT, TH1* MassSignPU); bool Combine(string InputPattern, string signal7, string signal8); bool useSample(int TypeMode, string sample); double MinRange = 0; double MaxRange = 1999; int CurrentSampleIndex; std::vector<stSample> samples; std::vector<std::string> modelVector; std::map<std::string, std::vector<stSample> > modelMap; string SHAPESTRING=""; void Analysis_Step6(string MODE="COMPILE", string InputPattern="", string signal=""){ setTDRStyle(); gStyle->SetPadTopMargin (0.06); gStyle->SetPadBottomMargin(0.12); gStyle->SetPadRightMargin (0.16); gStyle->SetPadLeftMargin (0.14); gStyle->SetTitleSize(0.05, "XYZ"); gStyle->SetTitleXOffset(1.1); gStyle->SetTitleYOffset(1.4); gStyle->SetPalette(1); gStyle->SetNdivisions(505,"X"); gStyle->SetNdivisions(550,"Y"); //gStyle->SetTextFont(43); if(MODE=="COMPILE")return; string Data; if(MODE.find("SHAPE")!=string::npos){SHAPESTRING="SHAPE";}else{SHAPESTRING="";} if(MODE.find("COMPUTELIMIT")!=string::npos || MODE.find("OPTIMIZE")!=string::npos){ if(signal.find("7TeV")!=string::npos){Data = "Data7TeV"; SQRTS=7.0; EXCLUSIONDIR+="7TeV"; } if(signal.find("8TeV")!=string::npos){Data = "Data8TeV"; SQRTS=8.0; EXCLUSIONDIR+="8TeV"; } printf("EXCLUSIONDIR = %s\nData = %s\n",EXCLUSIONDIR.c_str(), Data.c_str()); if(MODE.find("COMPUTELIMIT")!=string::npos){Optimize(InputPattern, Data, signal, SHAPESTRING!="", true); return;} if(MODE.find("OPTIMIZE")!=string::npos){ Optimize(InputPattern, Data, signal, SHAPESTRING!="", false); return;} //testShapeBasedAnalysis(InputPattern,signal); //use the second part if you want to run shape based analyssi on optimal point form c&c } if(MODE.find("COMBINE")!=string::npos){ printf("COMBINE!!!\n"); string signal7TeV = signal; if(signal7TeV.find("_8TeV")!=string::npos) signal7TeV = signal7TeV.replace(signal7TeV.find("_8TeV"),5, "_7TeV"); string signal8TeV = signal; if(signal8TeV.find("_7TeV")!=string::npos) signal8TeV = signal8TeV.replace(signal8TeV.find("_7TeV"),5, "_8TeV"); string EXCLUSIONDIR_SAVE = EXCLUSIONDIR; //2011 Limits Data = "Data7TeV"; SQRTS=7.0; EXCLUSIONDIR=EXCLUSIONDIR_SAVE+"7TeV"; Optimize(InputPattern, Data, signal7TeV, SHAPESTRING!="", true); //2012 Limits Data = "Data8TeV"; SQRTS=8.0; EXCLUSIONDIR=EXCLUSIONDIR_SAVE+"8TeV"; Optimize(InputPattern, Data, signal8TeV, SHAPESTRING!="", true); //Combined Limits EXCLUSIONDIR=EXCLUSIONDIR_SAVE+"COMB"; SQRTS=78.0; Combine(InputPattern, signal7TeV, signal8TeV); return; } if(MODE.find("7TeV")!=string::npos){Data = "Data7TeV"; SQRTS=7.0; EXCLUSIONDIR+="7TeV"; } if(MODE.find("8TeV")!=string::npos){Data = "Data8TeV"; SQRTS=8.0; EXCLUSIONDIR+="8TeV"; } printf("EXCLUSIONDIR = %s\nData = %s\n",EXCLUSIONDIR.c_str(), Data.c_str()); string TkPattern = "Results/Type0/"; string MuPattern = "Results/Type2/"; string MOPattern = "Results/Type3/"; string HQPattern = "Results/Type4/"; string LQPattern = "Results/Type5/"; bool Combine = (MODE.find("COMB")!=string::npos); if(Combine){EXCLUSIONDIR+="COMB"; SQRTS=78.0;} if(Combine) {PlotMinScale=1E-6; PlotMaxScale=300;} string outpath = string("Results/"+SHAPESTRING+EXCLUSIONDIR+"/"); MakeDirectories(outpath); //determine the list of models that are considered GetSampleDefinition(samples); if(SQRTS!=78.0) keepOnlySamplesAt7and8TeVX(samples, SQRTS); for(unsigned int s=0;s<samples.size();s++){ if(samples[s].Type!=2)continue; //printf("Name-->Model >> %30s --> %s\n",samples[s].Name.c_str(), samples[s].ModelName().c_str()); if(SQRTS== 7.0 && samples[s].Name.find("_7TeV")==string::npos){continue;} if(SQRTS== 8.0 && samples[s].Name.find("_8TeV")==string::npos){continue;} // if(SQRTS==78.0){if(samples[s].Name.find("_7TeV")==string::npos){continue;}else{samples[s].Name.replace(samples[s].Name.find("_7TeV"),5, ""); } } if(SQRTS==78.0){if(samples[s].Name.find("_8TeV")==string::npos){continue;}else{samples[s].Name.replace(samples[s].Name.find("_8TeV"),5, ""); } } modelMap[samples[s].ModelName()].push_back(samples[s]); if(modelMap[samples[s].ModelName()].size()==1)modelVector.push_back(samples[s].ModelName()); } //unti we have all the samples at both 7 and 8TeV, add the 7TeV models // if(SQRTS== 8.0){ // for(unsigned int s=0;s<samples.size();s++){ // if(samples[s].Type!=2)continue; // // if(modelMap.find(samples[s].ModelName())==modelMap.end()){ // modelMap[samples[s].ModelName()].push_back(samples[s]); // if(modelMap[samples[s].ModelName()].size()==1)modelVector.push_back(samples[s].ModelName()); // } // } // } //based on the modelMap //DrawRatioBands(TkPattern); //DrawRatioBands(MuPattern); //DrawRatioBands(MOPattern); //DrawRatioBands(LQPattern); //draw the cross section limit for all model DrawModelLimitWithBand(TkPattern); DrawModelLimitWithBand(MuPattern); DrawModelLimitWithBand(MOPattern); DrawModelLimitWithBand(LQPattern); //make plots of the observed limit for all signal model (and mass point) and save the result in a latex table TCanvas* c1; TLegend* LEG; double LInt = 0; FILE* pFile = fopen((outpath+string("Analysis_Step6_Result") + ".txt").c_str(),"w"); FILE* talkFile = fopen((outpath + "TalkPlots" + ".txt").c_str(),"w"); fprintf(pFile , "\\documentclass{article}\n"); fprintf(pFile , "\\begin{document}\n\n"); fprintf(talkFile, "\\documentclass{article}\n"); fprintf(talkFile, "\\usepackage{rotating}\n"); fprintf(talkFile, "\\begin{document}\n\n"); fprintf(talkFile, "\\begin{tiny}\n\n"); fprintf(pFile , "%% %50s\n", "TkOnly"); fprintf(pFile , "\\begin{table}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile, "\\begin{sidewaystable}\n \\centering\n \\begin{tabular}{|l|cccccccc|}\n \\hline\n"); fprintf(talkFile,"Sample & Mass(GeV) & Pt(GeV) & $I_{as}$ & TOF & Mass Cut (GeV) & N pred & N observed & Eff & Signif \\\\\n"); fprintf(talkFile, "\\hline\n"); TGraph** TkGraphs = new TGraph*[modelVector.size()]; for(unsigned int k=0; k<modelVector.size(); k++){ TkGraphs[k] = MakePlot(pFile,talkFile,TkPattern,modelVector[k], 2, modelMap[modelVector[k]], LInt); } fprintf(pFile ," \\end{tabular}\n\\end{table}\n\n"); fprintf(talkFile," \\end{tabular}\n\\end{sidewaystable}\n\n"); fprintf(pFile , "%% %50s\n", "TkMuon"); fprintf(pFile , "\\begin{table}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile, "\\begin{sidewaystable}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile,"Sample & Mass(GeV) & Pt(GeV) & $I_{as}$ & $#beta^{-1]$ & Mass Cut (GeV) & N pred & N observed & Eff \\\\\n"); fprintf(talkFile, "\\hline\n"); TGraph** MuGraphs = new TGraph*[modelVector.size()]; for(unsigned int k=0; k<modelVector.size(); k++){ bool isNeutral = false;if(modelVector[k].find("GluinoN")!=string::npos || modelVector[k].find("StopN")!=string::npos)isNeutral = true; if(isNeutral) continue;//skip charged suppressed models MuGraphs[k] = MakePlot(pFile,talkFile,MuPattern,modelVector[k], 2, modelMap[modelVector[k]], LInt); } fprintf(pFile ," \\end{tabular}\n\\end{table}\n\n"); fprintf(talkFile," \\end{tabular}\n\\end{sidewaystable}\n\n"); fprintf(pFile , "%% %50s\n", "MuOnly"); fprintf(pFile , "\\begin{table}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile, "\\begin{sidewaystable}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile,"Sample & Mass(GeV) & Pt(GeV) & $I_{as}$ & $#beta^{-1]$ & Mass Cut (GeV) & N pred & N observed & Eff \\\\\n"); fprintf(talkFile, "\\hline\n"); TGraph** MOGraphs = new TGraph*[modelVector.size()]; for(unsigned int k=0; k<modelVector.size(); k++){ bool isNeutral = false;if(modelVector[k].find("GluinoN")!=string::npos || modelVector[k].find("StopN")!=string::npos)isNeutral = true; if(isNeutral) continue;//skip charged suppressed models MOGraphs[k] = MakePlot(pFile,talkFile,MOPattern,modelVector[k], 2, modelMap[modelVector[k]], LInt); } fprintf(pFile ," \\end{tabular}\n\\end{table}\n\n"); fprintf(talkFile," \\end{tabular}\n\\end{sidewaystable}\n\n"); fprintf(pFile , "%% %50s\n", "multiple charge"); fprintf(pFile , "\\begin{table}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile, "\\begin{sidewaystable}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile,"Sample & Mass(GeV) & Pt(GeV) & $I_{as}$ & $#beta^{-1]$ & Mass Cut (GeV) & N pred & N observed & Eff \\\\\n"); fprintf(talkFile, "\\hline\n"); TGraph** HQGraphs = new TGraph*[modelVector.size()]; for(unsigned int k=0; k<modelVector.size(); k++){ bool ismHSCP = false; if(modelVector[k].find("DY_Q")!=string::npos && modelVector[k].find("o3")==string::npos ) ismHSCP = true; if(!ismHSCP) continue; // only mHSCP models HQGraphs[k] = MakePlot(pFile,talkFile,HQPattern,modelVector[k], 2, modelMap[modelVector[k]], LInt); } fprintf(pFile ," \\end{tabular}\n\\end{table}\n\n"); fprintf(talkFile," \\end{tabular}\n\\end{sidewaystable}\n\n"); fprintf(pFile , "%% %50s\n", "fractionnally charge"); fprintf(pFile , "\\begin{table}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile, "\\begin{sidewaystable}\n \\centering\n \\begin{tabular}{|l|cccccc|}\n \\hline\n"); fprintf(talkFile,"Sample & Mass(GeV) & Pt(GeV) & $I_{as}$ & $#beta^{-1]$ & Mass Cut (GeV) & N pred & N observed & Eff \\\\\n"); fprintf(talkFile, "\\hline\n"); TGraph** LQGraphs = new TGraph*[modelVector.size()]; for(unsigned int k=0; k<modelVector.size(); k++){ bool isFractional = false;if(modelVector[k].find("1o3")!=string::npos || modelVector[k].find("2o3")!=string::npos ||modelVector[k].find("Q1")!=string::npos)isFractional = true; if(!isFractional) continue;//skip q>=1 charged suppressed models LQGraphs[k] = MakePlot(pFile,talkFile,LQPattern,modelVector[k], 2, modelMap[modelVector[k]], LInt); } fprintf(pFile ," \\end{tabular}\n\\end{table}\n\n"); fprintf(talkFile," \\end{tabular}\n\\end{sidewaystable}\n\n"); fprintf(pFile ,"\\end{document}\n\n"); fprintf(talkFile,"\\end{document}\n\n"); if(SQRTS==8.0){ fprintf(pFile,"%%TKONLY\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummary(pFile, talkFile, TkPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"%%TKTOF\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummary(pFile, talkFile, MuPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"%%MUONLY\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummary(pFile, talkFile, MOPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"%%Q>1\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummary(pFile, talkFile, HQPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"%%Q<1\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummary(pFile, talkFile, LQPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"\n\n\n\n\n\n"); fprintf(pFile,"%%TKONLY\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummaryPaper(pFile, talkFile, TkPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"%%TKTOF\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummaryPaper(pFile, talkFile, MuPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"%%MUONLY\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummaryPaper(pFile, talkFile, MOPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"%%Q>1\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummaryPaper(pFile, talkFile, HQPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); fprintf(pFile,"%%Q<1\n"); fprintf(pFile,"Sample & Mass & Cut & \\multicolumn{4}{c|}{$\\sqrt{s}=7TeV$} & \\multicolumn{4}{c|}{$\\sqrt{s}=8TeV$} & \\multicolumn{2}{c|}{$\\sqrt{s}=7+8TeV$} \\\\\\hline\n"); fprintf(pFile," & (GeV) & (GeV) & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & Eff & $\\sigma_{TH}$ & $\\sigma_{obs}$ & $\\sigma_{pred}$ & $\\mu_{obs}$ & $\\mu_{pred}$ \\\\\\hline\n"); for(unsigned int k=0; k<modelVector.size(); k++){printSummaryPaper(pFile, talkFile, LQPattern , modelVector[k], modelMap[modelVector[k]]); } fprintf(pFile,"\\hline\n\n\n"); } //print a table with all uncertainty on signal efficiency c1 = new TCanvas("c1", "c1",600,600); c1->SetLeftMargin(0.15); TMultiGraph* TkSystGraphs = new TMultiGraph(); LEG = new TLegend(0.20,0.75,0.80,0.90); LEG->SetNColumns(2) ; LEG->SetFillColor(0); LEG->SetFillStyle(0); LEG->SetBorderSize(0); fprintf(pFile ,"\n\n %20s \n\n", LegendFromType(TkPattern).c_str()); fprintf(pFile , "%20s Eff --> PScale | DeDxScale | PUScale | TotalUncertainty \n","Model"); fprintf(talkFile, "\\hline\n%20s & Eff & PScale & DeDxScale & PUScale & TotalUncertainty \\\\\n","Model"); int Graphs=0; for(unsigned int k=0; k<modelVector.size(); k++){ TGraph* Uncertainty = CheckSignalUncertainty(pFile,talkFile,TkPattern, modelVector[k], modelMap[modelVector[k]]); if(Uncertainty!=NULL && useSample(0, modelVector[k])) { Uncertainty->SetLineColor(Color[Graphs]); Uncertainty->SetMarkerColor(Color[Graphs]); Uncertainty->SetMarkerStyle(20); Uncertainty->SetLineWidth(2); TkSystGraphs->Add(Uncertainty,"LP"); LEG->AddEntry(Uncertainty, modelVector[k].c_str() ,"L"); Graphs++; } } if(Graphs>0) { TkSystGraphs->Draw("A"); TkSystGraphs->SetTitle(""); TkSystGraphs->GetXaxis()->SetTitle("Mass (GeV)"); TkSystGraphs->GetYaxis()->SetTitle("Relative Uncertainty"); TkSystGraphs->GetYaxis()->SetTitleOffset(1.40); TkSystGraphs->GetYaxis()->SetRangeUser(0.0, 0.60); TkSystGraphs->GetYaxis()->SetNdivisions(520, "X"); LEG->Draw(); c1->SetLogy(false); c1->SetGridy(false); DrawPreliminary(LegendFromType(TkPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,"Results/"+SHAPESTRING+EXCLUSIONDIR+"/", "TkUncertainty"); delete c1; delete TkSystGraphs; delete LEG; } c1 = new TCanvas("c1", "c1",600,600); c1->SetLeftMargin(0.15); TMultiGraph* MuSystGraphs = new TMultiGraph(); LEG = new TLegend(0.20,0.75,0.80,0.90); LEG->SetNColumns(2) ; LEG->SetFillColor(0); LEG->SetFillStyle(0); LEG->SetBorderSize(0); fprintf(pFile ,"\n\n %20s \n\n", LegendFromType(MuPattern).c_str()); fprintf(pFile, "%20s Eff --> PScale | DeDxScale | PUScale | TOFScale | TotalUncertainty \n","Model"); fprintf(talkFile, "\\hline\n%20s & Eff & PScale & DeDxScale & PUScale & TOFScale & TotalUncertainty \\\\\n","Model"); Graphs=0; for(unsigned int k=0; k<modelVector.size(); k++){ TGraph* Uncertainty = CheckSignalUncertainty(pFile,talkFile,MuPattern, modelVector[k], modelMap[modelVector[k]]); if(Uncertainty!=NULL && useSample(2, modelVector[k])) { Uncertainty->SetLineColor(Color[Graphs]); Uncertainty->SetMarkerColor(Color[Graphs]); Uncertainty->SetMarkerStyle(20); Uncertainty->SetLineWidth(2); MuSystGraphs->Add(Uncertainty,"LP"); LEG->AddEntry(Uncertainty, modelVector[k].c_str() ,"L"); Graphs++; } } if(Graphs>0) { MuSystGraphs->Draw("A"); MuSystGraphs->SetTitle(""); MuSystGraphs->GetXaxis()->SetTitle("Mass (GeV)"); MuSystGraphs->GetYaxis()->SetTitle("Relative Uncertainty"); MuSystGraphs->GetYaxis()->SetTitleOffset(1.40); MuSystGraphs->GetYaxis()->SetRangeUser(0.0, 0.6); MuSystGraphs->GetYaxis()->SetNdivisions(520, "X"); LEG->Draw(); c1->SetLogy(false); c1->SetGridy(false); DrawPreliminary(LegendFromType(MuPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,"Results/"+SHAPESTRING+EXCLUSIONDIR+"/", "MuUncertainty"); delete c1; delete MuSystGraphs; delete LEG; } fprintf(pFile ,"\n\n %20s \n\n", LegendFromType(MOPattern).c_str()); fprintf(pFile, "%20s Eff --> PScale | DeDxScale | PUScale | TOFScale | TotalUncertainty \n","Model"); fprintf(talkFile, "\\hline\n%20s & Eff & PScale & DeDxScale & PUScale & TOFScale & TotalUncertainty \\\\\n","Model"); c1 = new TCanvas("c1", "c1",600,600); c1->SetLeftMargin(0.15); TMultiGraph* MOSystGraphs = new TMultiGraph(); LEG = new TLegend(0.20,0.75,0.80,0.90); LEG->SetNColumns(2) ; LEG->SetFillColor(0); LEG->SetFillStyle(0); LEG->SetBorderSize(0); Graphs=0; for(unsigned int k=0; k<modelVector.size(); k++){ TGraph* Uncertainty = CheckSignalUncertainty(pFile,talkFile,MOPattern, modelVector[k], modelMap[modelVector[k]]); if(Uncertainty!=NULL && useSample(3, modelVector[k])) { Uncertainty->SetLineColor(Color[Graphs]); Uncertainty->SetMarkerColor(Color[Graphs]); Uncertainty->SetMarkerStyle(20); Uncertainty->SetLineWidth(2); MOSystGraphs->Add(Uncertainty,"LP"); LEG->AddEntry(Uncertainty, modelVector[k].c_str() ,"L"); Graphs++; } } if(Graphs>0) { MOSystGraphs->Draw("A"); MOSystGraphs->SetTitle(""); MOSystGraphs->GetXaxis()->SetTitle("Mass (GeV)"); MOSystGraphs->GetYaxis()->SetTitle("Relative Uncertainty"); MOSystGraphs->GetYaxis()->SetTitleOffset(1.40); MOSystGraphs->GetYaxis()->SetRangeUser(0.0, 0.6); MOSystGraphs->GetYaxis()->SetNdivisions(520, "X"); LEG->Draw(); c1->SetLogy(false); c1->SetGridy(false); DrawPreliminary(LegendFromType(MOPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,"Results/"+SHAPESTRING+EXCLUSIONDIR+"/", "MOUncertainty"); delete c1; delete MOSystGraphs; delete LEG; } c1 = new TCanvas("c1", "c1",600,600); c1->SetLeftMargin(0.15); TMultiGraph* LQSystGraphs = new TMultiGraph(); LEG = new TLegend(0.20,0.75,0.80,0.90); LEG->SetNColumns(2) ; LEG->SetFillColor(0); LEG->SetFillStyle(0); LEG->SetBorderSize(0); fprintf(pFile ,"\n\n %20s \n\n", LegendFromType(LQPattern).c_str()); fprintf(pFile , "%20s Eff --> PScale | DeDxScale | PUScale | TotalUncertainty \n","Model"); fprintf(talkFile, "\\hline\n%20s & Eff & PScale & DeDxScale & PUScale & TotalUncertainty \\\\\n","Model"); Graphs=0; for(unsigned int k=0; k<modelVector.size(); k++){ TGraph* Uncertainty = CheckSignalUncertainty(pFile,talkFile,LQPattern, modelVector[k], modelMap[modelVector[k]]); if(Uncertainty!=NULL && useSample(5, modelVector[k])) { Uncertainty->SetLineColor(Color[Graphs]); Uncertainty->SetMarkerColor(Color[Graphs]); Uncertainty->SetMarkerStyle(20); Uncertainty->SetLineWidth(2); LQSystGraphs->Add(Uncertainty,"LP"); LEG->AddEntry(Uncertainty, modelVector[k].c_str() ,"L"); Graphs++; } } if(Graphs>0) { LQSystGraphs->Draw("A"); LQSystGraphs->SetTitle(""); LQSystGraphs->GetXaxis()->SetTitle("Mass (GeV)"); LQSystGraphs->GetYaxis()->SetTitle("Relative Uncertainty"); LQSystGraphs->GetYaxis()->SetTitleOffset(1.40); LQSystGraphs->GetYaxis()->SetRangeUser(0.0, 0.6); LQSystGraphs->GetYaxis()->SetNdivisions(520, "X"); LEG->Draw(); c1->SetLogy(false); c1->SetGridy(false); DrawPreliminary(LegendFromType(LQPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,"Results/"+SHAPESTRING+EXCLUSIONDIR+"/", "LQUncertainty"); delete c1; delete LQSystGraphs; delete LEG; } c1 = new TCanvas("c1", "c1",600,600); c1->SetLeftMargin(0.15); TMultiGraph* HQSystGraphs = new TMultiGraph(); LEG = new TLegend(0.20,0.75,0.80,0.90); LEG->SetNColumns(2) ; LEG->SetFillColor(0); LEG->SetFillStyle(0); LEG->SetBorderSize(0); fprintf(pFile ,"\n\n %20s \n\n", LegendFromType(HQPattern).c_str()); fprintf(pFile , "%20s Eff --> PScale | DeDxScale | PUScale | TotalUncertainty \n","Model"); fprintf(talkFile, "\\hline\n%20s & Eff & PScale & DeDxScale & PUScale & TotalUncertainty \\\\\n","Model"); Graphs=0; for(unsigned int k=0; k<modelVector.size(); k++){ TGraph* Uncertainty = CheckSignalUncertainty(pFile,talkFile,HQPattern, modelVector[k], modelMap[modelVector[k]]); if(Uncertainty!=NULL && useSample(4, modelVector[k])) { Uncertainty->SetLineColor(Color[Graphs]); Uncertainty->SetMarkerColor(Color[Graphs]); Uncertainty->SetMarkerStyle(20); Uncertainty->SetLineWidth(2); HQSystGraphs->Add(Uncertainty,"LP"); LEG->AddEntry(Uncertainty, modelVector[k].c_str() ,"L"); Graphs++; } } if(Graphs>0) { HQSystGraphs->Draw("A"); HQSystGraphs->SetTitle(""); HQSystGraphs->GetXaxis()->SetTitle("Mass (GeV)"); HQSystGraphs->GetYaxis()->SetTitle("Relative Uncertainty"); HQSystGraphs->GetYaxis()->SetTitleOffset(1.40); HQSystGraphs->GetYaxis()->SetRangeUser(0.0, 0.6); HQSystGraphs->GetYaxis()->SetNdivisions(520, "X"); LEG->Draw(); c1->SetLogy(false); c1->SetGridy(false); DrawPreliminary(LegendFromType(HQPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,"Results/"+SHAPESTRING+EXCLUSIONDIR+"/", "HQUncertainty"); delete c1; delete HQSystGraphs; delete LEG; } std::cout<<"TESTA\n"; //Get Theoretical xsection and error bands TGraph** ThXSec = new TGraph*[modelVector.size()]; TCutG ** ThXSecErr = new TCutG* [modelVector.size()]; for(unsigned int k=0; k<modelVector.size(); k++){ if(modelVector[k].find("Gluino")!=string::npos){ if(SQRTS==7){ ThXSec [k] = new TGraph(sizeof(THXSEC7TeV_Gluino_Mass)/sizeof(double),THXSEC7TeV_Gluino_Mass,THXSEC7TeV_Gluino_Cen); ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr",sizeof(THXSEC7TeV_Gluino_Mass)/sizeof(double),THXSEC7TeV_Gluino_Mass,THXSEC7TeV_Gluino_Low,THXSEC7TeV_Gluino_High, PlotMinScale, PlotMaxScale); }else if(SQRTS==8){ ThXSec [k] = new TGraph(sizeof(THXSEC8TeV_Gluino_Mass)/sizeof(double),THXSEC8TeV_Gluino_Mass,THXSEC8TeV_Gluino_Cen); ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr",sizeof(THXSEC8TeV_Gluino_Mass)/sizeof(double),THXSEC8TeV_Gluino_Mass,THXSEC8TeV_Gluino_Low,THXSEC8TeV_Gluino_High, PlotMinScale, PlotMaxScale); } else { const int NMass=sizeof(THXSEC8TeV_Gluino_Mass)/sizeof(double); double ones[NMass]; for(int i=0; i<NMass; i++) ones[i]=1; ThXSec [k] = new TGraph(NMass,THXSEC8TeV_Gluino_Mass,ones); } }else if(modelVector[k].find("Stop" )!=string::npos){ if(SQRTS==7){ ThXSec [k] = new TGraph(sizeof(THXSEC7TeV_Stop_Mass)/sizeof(double),THXSEC7TeV_Stop_Mass,THXSEC7TeV_Stop_Cen); ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr",sizeof(THXSEC7TeV_Stop_Mass)/sizeof(double),THXSEC7TeV_Stop_Mass,THXSEC7TeV_Stop_Low,THXSEC7TeV_Stop_High, PlotMinScale, PlotMaxScale); }else if(SQRTS==8){ ThXSec [k] = new TGraph(sizeof(THXSEC8TeV_Stop_Mass)/sizeof(double),THXSEC8TeV_Stop_Mass,THXSEC8TeV_Stop_Cen); ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr",sizeof(THXSEC8TeV_Stop_Mass)/sizeof(double),THXSEC8TeV_Stop_Mass,THXSEC8TeV_Stop_Low,THXSEC8TeV_Stop_High, PlotMinScale, PlotMaxScale); } else { const int NMass=sizeof(THXSEC8TeV_Stop_Mass)/sizeof(double); double ones[NMass]; for(int i=0; i<NMass; i++) ones[i]=1; ThXSec [k] = new TGraph(NMass,THXSEC8TeV_Stop_Mass,ones); } }else if(modelVector[k].find("GMStau" )!=string::npos){ if(SQRTS==7){ ThXSec [k] = MakePlot(NULL, NULL, TkPattern,modelVector[k], 0, modelMap[modelVector[k]], LInt); ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr", sizeof(THXSEC7TeV_GMStau_Mass)/sizeof(double),THXSEC7TeV_GMStau_Mass,THXSEC7TeV_GMStau_Low,THXSEC7TeV_GMStau_High, PlotMinScale, PlotMaxScale); }else if(SQRTS==8){ ThXSec [k] = MakePlot(NULL, NULL, TkPattern,modelVector[k], 0, modelMap[modelVector[k]], LInt); ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr", sizeof(THXSEC8TeV_GMStau_Mass)/sizeof(double),THXSEC8TeV_GMStau_Mass,THXSEC8TeV_GMStau_Low,THXSEC8TeV_GMStau_High, PlotMinScale, PlotMaxScale); } else { const int NMass=sizeof(THXSEC8TeV_GMStau_Mass)/sizeof(double); double ones[NMass]; for(int i=0; i<NMass; i++) ones[i]=1; ThXSec [k] = new TGraph(NMass,THXSEC8TeV_GMStau_Mass,ones); } }else if(modelVector[k].find("PPStau" )!=string::npos){ if(SQRTS==7){ ThXSec [k] = MakePlot(NULL, NULL, TkPattern,modelVector[k], 0, modelMap[modelVector[k]], LInt); ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr", sizeof(THXSEC7TeV_PPStau_Mass)/sizeof(double),THXSEC7TeV_PPStau_Mass,THXSEC7TeV_PPStau_Low,THXSEC7TeV_PPStau_High, PlotMinScale, PlotMaxScale); }else if(SQRTS==8){ ThXSec [k] = MakePlot(NULL, NULL, TkPattern,modelVector[k], 0, modelMap[modelVector[k]], LInt); ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr", sizeof(THXSEC8TeV_PPStau_Mass)/sizeof(double),THXSEC8TeV_PPStau_Mass,THXSEC8TeV_PPStau_Low,THXSEC8TeV_PPStau_High, PlotMinScale, PlotMaxScale); } else { const int NMass=sizeof(THXSEC8TeV_PPStau_Mass)/sizeof(double); double ones[NMass]; for(int i=0; i<NMass; i++) ones[i]=1; ThXSec [k] = new TGraph(NMass,THXSEC8TeV_PPStau_Mass,ones); } }else{ if(modelVector[k].find("o3")!=string::npos){ ThXSec [k] = MakePlot(NULL, NULL, LQPattern,modelVector[k], 0, modelMap[modelVector[k]], LInt); }else if(modelVector[k].find("DY_Q")!=string::npos){ ThXSec [k] = MakePlot(NULL, NULL, HQPattern,modelVector[k], 0, modelMap[modelVector[k]], LInt); }else{ ThXSec [k] = MakePlot(NULL, NULL, TkPattern,modelVector[k], 0, modelMap[modelVector[k]], LInt); } double* XSecErrLow = new double[ThXSec[k]->GetN()]; double* XSecErrHigh = new double[ThXSec[k]->GetN()]; //assume 10% error on xsection for(int i=0;i<ThXSec[k]->GetN();i++){ XSecErrLow[i] = ThXSec[k]->GetY()[i]*0.90; XSecErrHigh[i] = ThXSec[k]->GetY()[i]*1.10; } ThXSecErr[k] = GetErrorBand(modelVector[k]+"ThErr", ThXSec[k]->GetN(),ThXSec[k]->GetX(),XSecErrLow,XSecErrHigh, PlotMinScale, PlotMaxScale); } } std::cout<<"TESTB\n"; //Print the excluded mass range fprintf(pFile,"\n\n\n-----------------------\n Mass range excluded \n-------------------------\n"); fprintf(pFile,"-----------------------\n0%% TK-ONLY \n-------------------------\n"); for(unsigned int k=0; k<modelVector.size(); k++){ if(TkGraphs[k]->GetN()==0) continue; if(TkGraphs[k]->GetX()[TkGraphs[k]->GetN()-1]<0) continue; fprintf(pFile,"%20s --> Excluded mass below %8.3fGeV\n", modelVector[k].c_str(), FindIntersectionBetweenTwoGraphs(TkGraphs[k], ThXSec[k], TkGraphs[k]->GetX()[0], TkGraphs[k]->GetX()[TkGraphs[k]->GetN()-1], 1, 0.00)); } fprintf(pFile,"-----------------------\n0%% MU+TOF \n-------------------------\n"); for(unsigned int k=0; k<modelVector.size(); k++){ bool isNeutral = false;if(modelVector[k].find("GluinoN")!=string::npos || modelVector[k].find("StopN")!=string::npos)isNeutral = true; if(isNeutral) continue;//skip charged suppressed models if(MuGraphs[k]->GetN()==0) continue; if(MuGraphs[k]->GetX()[MuGraphs[k]->GetN()-1]<0) continue; double minMass=-1, maxMass=-1; FindRangeBetweenTwoGraphs(MuGraphs[k], ThXSec[k], MuGraphs[k]->GetX()[0], MuGraphs[k]->GetX()[MuGraphs[k]->GetN()-1], 1, 0.00, minMass, maxMass); fprintf(pFile,"%20s --> Excluded mass range %8.3f - %8.3fGeV\n", modelVector[k].c_str(), minMass, maxMass); //fprintf(pFile,"%20s --> Excluded mass below %8.3fGeV\n", modelVector[k].c_str(), FindIntersectionBetweenTwoGraphs(MuGraphs[k], ThXSec[k], MuGraphs[k]->GetX()[0], MuGraphs[k]->GetX()[MuGraphs[k]->GetN()-1], 1, 0.00)); } fprintf(pFile,"-----------------------\n0%% MU+Only \n-------------------------\n"); for(unsigned int k=0; k<modelVector.size(); k++){ bool isNeutral = false;if(modelVector[k].find("GluinoN")!=string::npos || modelVector[k].find("StopN")!=string::npos)isNeutral = true; if(isNeutral) continue;//skip charged suppressed models if(MOGraphs[k]->GetN()==0) continue; if(MOGraphs[k]->GetX()[MOGraphs[k]->GetN()-1]<0) continue; fprintf(pFile,"%20s --> Excluded mass below %8.3fGeV\n", modelVector[k].c_str(), FindIntersectionBetweenTwoGraphs(MOGraphs[k], ThXSec[k], MOGraphs[k]->GetX()[0], MOGraphs[k]->GetX()[MOGraphs[k]->GetN()-1], 1, 0.00)); } fprintf(pFile,"-----------------------\n0%% Q>1 \n-------------------------\n"); for(unsigned int k=0; k<modelVector.size(); k++){ bool ismHSCP = false; if(modelVector[k].find("DY_Q")!=string::npos && modelVector[k].find("o3")==string::npos ) ismHSCP = true; if(!ismHSCP) continue; if(HQGraphs[k]->GetN()==0) continue; if(HQGraphs[k]->GetX()[HQGraphs[k]->GetN()-1]<0) continue; double minMass=-1, maxMass=-1; FindRangeBetweenTwoGraphs(HQGraphs[k], ThXSec[k], HQGraphs[k]->GetX()[0], HQGraphs[k]->GetX()[HQGraphs[k]->GetN()-1], 1, 0.00, minMass, maxMass); fprintf(pFile,"%20s --> Excluded mass range %8.3f - %8.3fGeV\n", modelVector[k].c_str(), minMass, maxMass); //fprintf(pFile,"%20s --> Excluded mass below %8.3fGeV\n", modelVector[k].c_str(), FindIntersectionBetweenTwoGraphs(HQGraphs[k], ThXSec[k], HQGraphs[k]->GetX()[0], HQGraphs[k]->GetX()[HQGraphs[k]->GetN()-1], 1, 0.00)); } fprintf(pFile,"-----------------------\n0%% Q<1 \n-------------------------\n"); for(unsigned int k=0; k<modelVector.size(); k++){ bool isFractional = false;if(modelVector[k].find("1o3")!=string::npos || modelVector[k].find("2o3")!=string::npos || modelVector[k].find("DY_Q1")!=string::npos)isFractional = true; if(!isFractional) continue;//skip non fractional charge models if(LQGraphs[k]->GetN()==0) continue; if(LQGraphs[k]->GetX()[LQGraphs[k]->GetN()-1]<0) continue; fprintf(pFile,"%20s --> Excluded mass below %8.3fGeV\n", modelVector[k].c_str(), FindIntersectionBetweenTwoGraphs(LQGraphs[k], ThXSec[k], LQGraphs[k]->GetX()[0], LQGraphs[k]->GetX()[LQGraphs[k]->GetN()-1], 1, 0.00)); } fclose(pFile); std::cout<<"TESTC\n"; //Make the final plot with all curves in it // I don't like much this part because it is dependent of what is in Analysis_Samples.h in an hardcoded way std::map<string, TGraph*> TkGraphMap; std::map<string, TGraph*> MuGraphMap; std::map<string, TGraph*> LQGraphMap; std::map<string, TGraph*> HQGraphMap; std::map<string, TGraph*> MOGraphMap; std::map<string, TGraph*> ThGraphMap; std::map<string, TCutG* > ThErrorMap; for(unsigned int k=0; k<modelVector.size(); k++){ TkGraphMap[modelVector[k]] = TkGraphs [k]; MuGraphMap[modelVector[k]] = MuGraphs [k]; LQGraphMap[modelVector[k]] = LQGraphs [k]; HQGraphMap[modelVector[k]] = HQGraphs [k]; MOGraphMap[modelVector[k]] = MOGraphs [k]; ThGraphMap[modelVector[k]] = ThXSec [k]; ThErrorMap[modelVector[k]] = ThXSecErr[k]; } ThGraphMap["Gluino_f10" ]->SetLineColor(4); ThGraphMap["Gluino_f10" ]->SetMarkerColor(4); ThGraphMap["Gluino_f10" ]->SetLineWidth(1); ThGraphMap["Gluino_f10" ]->SetLineStyle(1); ThGraphMap["Gluino_f10" ]->SetMarkerStyle(1); MuGraphMap["Gluino_f10" ]->SetLineColor(4); MuGraphMap["Gluino_f10" ]->SetMarkerColor(4); MuGraphMap["Gluino_f10" ]->SetLineWidth(2); MuGraphMap["Gluino_f10" ]->SetLineStyle(1); MuGraphMap["Gluino_f10" ]->SetMarkerStyle(22); MuGraphMap["Gluino_f50" ]->SetLineColor(4); MuGraphMap["Gluino_f50" ]->SetMarkerColor(4); MuGraphMap["Gluino_f50" ]->SetLineWidth(2); MuGraphMap["Gluino_f50" ]->SetLineStyle(1); MuGraphMap["Gluino_f50" ]->SetMarkerStyle(23); TkGraphMap["Gluino_f10" ]->SetLineColor(4); TkGraphMap["Gluino_f10" ]->SetMarkerColor(4); TkGraphMap["Gluino_f10" ]->SetLineWidth(2); TkGraphMap["Gluino_f10" ]->SetLineStyle(1); TkGraphMap["Gluino_f10" ]->SetMarkerStyle(22); TkGraphMap["Gluino_f50" ]->SetLineColor(4); TkGraphMap["Gluino_f50" ]->SetMarkerColor(4); TkGraphMap["Gluino_f50" ]->SetLineWidth(2); TkGraphMap["Gluino_f50" ]->SetLineStyle(1); TkGraphMap["Gluino_f50" ]->SetMarkerStyle(23); TkGraphMap["GluinoN_f10" ]->SetLineColor(4); TkGraphMap["GluinoN_f10" ]->SetMarkerColor(4); TkGraphMap["GluinoN_f10" ]->SetLineWidth(2); TkGraphMap["GluinoN_f10" ]->SetLineStyle(1); TkGraphMap["GluinoN_f10" ]->SetMarkerStyle(26); MOGraphMap["Gluino_f10" ]->SetLineColor(4); MOGraphMap["Gluino_f10" ]->SetMarkerColor(4); MOGraphMap["Gluino_f10" ]->SetLineWidth(2); MOGraphMap["Gluino_f10" ]->SetLineStyle(1); MOGraphMap["Gluino_f10" ]->SetMarkerStyle(22); MOGraphMap["Gluino_f50" ]->SetLineColor(4); MOGraphMap["Gluino_f50" ]->SetMarkerColor(4); MOGraphMap["Gluino_f50" ]->SetLineWidth(2); MOGraphMap["Gluino_f50" ]->SetLineStyle(1); MOGraphMap["Gluino_f50" ]->SetMarkerStyle(23); MOGraphMap["Gluino_f100" ]->SetLineColor(4); MOGraphMap["Gluino_f100" ]->SetMarkerColor(4); MOGraphMap["Gluino_f100" ]->SetLineWidth(2); MOGraphMap["Gluino_f100" ]->SetLineStyle(1); MOGraphMap["Gluino_f100" ]->SetMarkerStyle(26); ThGraphMap["Stop" ]->SetLineColor(2); ThGraphMap["Stop" ]->SetMarkerColor(2); ThGraphMap["Stop" ]->SetLineWidth(1); ThGraphMap["Stop" ]->SetLineStyle(2); ThGraphMap["Stop" ]->SetMarkerStyle(1); MuGraphMap["Stop" ]->SetLineColor(2); MuGraphMap["Stop" ]->SetMarkerColor(2); MuGraphMap["Stop" ]->SetLineWidth(2); MuGraphMap["Stop" ]->SetLineStyle(1); MuGraphMap["Stop" ]->SetMarkerStyle(21); TkGraphMap["Stop" ]->SetLineColor(2); TkGraphMap["Stop" ]->SetMarkerColor(2); TkGraphMap["Stop" ]->SetLineWidth(2); TkGraphMap["Stop" ]->SetLineStyle(1); TkGraphMap["Stop" ]->SetMarkerStyle(21); TkGraphMap["StopN" ]->SetLineColor(2); TkGraphMap["StopN" ]->SetMarkerColor(2); TkGraphMap["StopN" ]->SetLineWidth(2); TkGraphMap["StopN" ]->SetLineStyle(1); TkGraphMap["StopN" ]->SetMarkerStyle(25); MOGraphMap["Stop" ]->SetLineColor(2); MOGraphMap["Stop" ]->SetMarkerColor(2); MOGraphMap["Stop" ]->SetLineWidth(2); MOGraphMap["Stop" ]->SetLineStyle(1); MOGraphMap["Stop" ]->SetMarkerStyle(21); ThGraphMap["GMStau" ]->SetLineColor(1); ThGraphMap["GMStau" ]->SetMarkerColor(1); ThGraphMap["GMStau" ]->SetLineWidth(1); ThGraphMap["GMStau" ]->SetLineStyle(3); ThGraphMap["GMStau" ]->SetMarkerStyle(1); ThGraphMap["PPStau" ]->SetLineColor(6); ThGraphMap["PPStau" ]->SetMarkerColor(6); ThGraphMap["PPStau" ]->SetLineWidth(1); ThGraphMap["PPStau" ]->SetLineStyle(4); ThGraphMap["PPStau" ]->SetMarkerStyle(1); // ThGraphMap["DCRho08HyperK"]->SetLineColor(4); ThGraphMap["DCRho08HyperK"]->SetMarkerColor(4); ThGraphMap["DCRho08HyperK"]->SetLineWidth(1); ThGraphMap["DCRho08HyperK"]->SetLineStyle(3); ThGraphMap["DCRho08HyperK"]->SetMarkerStyle(1); //ThGraphMap["DCRho12HyperK"]->SetLineColor(2); ThGraphMap["DCRho12HyperK"]->SetMarkerColor(2); ThGraphMap["DCRho12HyperK"]->SetLineWidth(1); ThGraphMap["DCRho12HyperK"]->SetLineStyle(2); ThGraphMap["DCRho12HyperK"]->SetMarkerStyle(1); //ThGraphMap["DCRho16HyperK"]->SetLineColor(1); ThGraphMap["DCRho16HyperK"]->SetMarkerColor(1); ThGraphMap["DCRho16HyperK"]->SetLineWidth(1); ThGraphMap["DCRho16HyperK"]->SetLineStyle(1); ThGraphMap["DCRho16HyperK"]->SetMarkerStyle(1); MuGraphMap["GMStau" ]->SetLineColor(1); MuGraphMap["GMStau" ]->SetMarkerColor(1); MuGraphMap["GMStau" ]->SetLineWidth(2); MuGraphMap["GMStau" ]->SetLineStyle(1); MuGraphMap["GMStau" ]->SetMarkerStyle(20); MuGraphMap["PPStau" ]->SetLineColor(6); MuGraphMap["PPStau" ]->SetMarkerColor(6); MuGraphMap["PPStau" ]->SetLineWidth(2); MuGraphMap["PPStau" ]->SetLineStyle(1); MuGraphMap["PPStau" ]->SetMarkerStyle(20); // MuGraphMap["DCRho08HyperK"]->SetLineColor(4); MuGraphMap["DCRho08HyperK"]->SetMarkerColor(4); MuGraphMap["DCRho08HyperK"]->SetLineWidth(2); MuGraphMap["DCRho08HyperK"]->SetLineStyle(1); MuGraphMap["DCRho08HyperK"]->SetMarkerStyle(22); //MuGraphMap["DCRho12HyperK"]->SetLineColor(2); MuGraphMap["DCRho12HyperK"]->SetMarkerColor(2); MuGraphMap["DCRho12HyperK"]->SetLineWidth(2); MuGraphMap["DCRho12HyperK"]->SetLineStyle(1); MuGraphMap["DCRho12HyperK"]->SetMarkerStyle(23); //MuGraphMap["DCRho16HyperK"]->SetLineColor(1); MuGraphMap["DCRho16HyperK"]->SetMarkerColor(1); MuGraphMap["DCRho16HyperK"]->SetLineWidth(2); MuGraphMap["DCRho16HyperK"]->SetLineStyle(1); MuGraphMap["DCRho16HyperK"]->SetMarkerStyle(26); TkGraphMap["GMStau" ]->SetLineColor(1); TkGraphMap["GMStau" ]->SetMarkerColor(1); TkGraphMap["GMStau" ]->SetLineWidth(2); TkGraphMap["GMStau" ]->SetLineStyle(1); TkGraphMap["GMStau" ]->SetMarkerStyle(20); TkGraphMap["PPStau" ]->SetLineColor(6); TkGraphMap["PPStau" ]->SetMarkerColor(6); TkGraphMap["PPStau" ]->SetLineWidth(2); TkGraphMap["PPStau" ]->SetLineStyle(1); TkGraphMap["PPStau" ]->SetMarkerStyle(20); // TkGraphMap["DCRho08HyperK"]->SetLineColor(4); TkGraphMap["DCRho08HyperK"]->SetMarkerColor(4); TkGraphMap["DCRho08HyperK"]->SetLineWidth(2); TkGraphMap["DCRho08HyperK"]->SetLineStyle(1); TkGraphMap["DCRho08HyperK"]->SetMarkerStyle(22); //TkGraphMap["DCRho12HyperK"]->SetLineColor(2); TkGraphMap["DCRho12HyperK"]->SetMarkerColor(2); TkGraphMap["DCRho12HyperK"]->SetLineWidth(2); TkGraphMap["DCRho12HyperK"]->SetLineStyle(1); TkGraphMap["DCRho12HyperK"]->SetMarkerStyle(23); // TkGraphMap["DCRho16HyperK"]->SetLineColor(1); TkGraphMap["DCRho16HyperK"]->SetMarkerColor(1); TkGraphMap["DCRho16HyperK"]->SetLineWidth(2); TkGraphMap["DCRho16HyperK"]->SetLineStyle(1); TkGraphMap["DCRho16HyperK"]->SetMarkerStyle(26); ThGraphMap["DY_Q1o3" ]->SetLineColor(41); ThGraphMap["DY_Q1o3" ]->SetMarkerColor(41); ThGraphMap["DY_Q1o3" ]->SetLineWidth(1); ThGraphMap["DY_Q1o3" ]->SetLineStyle(9); ThGraphMap["DY_Q1o3" ]->SetMarkerStyle(1); TkGraphMap["DY_Q1o3" ]->SetLineColor(41); TkGraphMap["DY_Q1o3" ]->SetMarkerColor(41); TkGraphMap["DY_Q1o3" ]->SetLineWidth(2); TkGraphMap["DY_Q1o3" ]->SetLineStyle(1); TkGraphMap["DY_Q1o3" ]->SetMarkerStyle(33); MOGraphMap["DY_Q1o3" ]->SetLineColor(41); MOGraphMap["DY_Q1o3" ]->SetMarkerColor(41); MOGraphMap["DY_Q1o3" ]->SetLineWidth(2); MOGraphMap["DY_Q1o3" ]->SetLineStyle(1); MOGraphMap["DY_Q1o3" ]->SetMarkerStyle(33); LQGraphMap["DY_Q1o3" ]->SetLineColor(41); LQGraphMap["DY_Q1o3" ]->SetMarkerColor(41); LQGraphMap["DY_Q1o3" ]->SetLineWidth(2); LQGraphMap["DY_Q1o3" ]->SetLineStyle(1); LQGraphMap["DY_Q1o3" ]->SetMarkerStyle(33); ThGraphMap["DY_Q2o3" ]->SetLineColor(43); ThGraphMap["DY_Q2o3" ]->SetMarkerColor(43); ThGraphMap["DY_Q2o3" ]->SetLineWidth(1); ThGraphMap["DY_Q2o3" ]->SetLineStyle(10); ThGraphMap["DY_Q2o3" ]->SetMarkerStyle(1); TkGraphMap["DY_Q2o3" ]->SetLineColor(43); TkGraphMap["DY_Q2o3" ]->SetMarkerColor(43); TkGraphMap["DY_Q2o3" ]->SetLineWidth(2); TkGraphMap["DY_Q2o3" ]->SetLineStyle(1); TkGraphMap["DY_Q2o3" ]->SetMarkerStyle(34); MuGraphMap["DY_Q2o3" ]->SetLineColor(43); MuGraphMap["DY_Q2o3" ]->SetMarkerColor(43); MuGraphMap["DY_Q2o3" ]->SetLineWidth(2); MuGraphMap["DY_Q2o3" ]->SetLineStyle(1); MuGraphMap["DY_Q2o3" ]->SetMarkerStyle(34); MOGraphMap["DY_Q2o3" ]->SetLineColor(43); MOGraphMap["DY_Q2o3" ]->SetMarkerColor(43); MOGraphMap["DY_Q2o3" ]->SetLineWidth(2); MOGraphMap["DY_Q2o3" ]->SetLineStyle(1); MOGraphMap["DY_Q2o3" ]->SetMarkerStyle(34); LQGraphMap["DY_Q2o3" ]->SetLineColor(43); LQGraphMap["DY_Q2o3" ]->SetMarkerColor(43); LQGraphMap["DY_Q2o3" ]->SetLineWidth(2); LQGraphMap["DY_Q2o3" ]->SetLineStyle(1); LQGraphMap["DY_Q2o3" ]->SetMarkerStyle(34); LQGraphMap["DY_Q1" ]->SetLineColor(46); LQGraphMap["DY_Q1" ]->SetMarkerColor(46); LQGraphMap["DY_Q1" ]->SetLineWidth(2); LQGraphMap["DY_Q1" ]->SetLineStyle(1); LQGraphMap["DY_Q1" ]->SetMarkerStyle(20); ThGraphMap["DY_Q1" ]->SetLineColor(46); ThGraphMap["DY_Q1" ]->SetMarkerColor(46); ThGraphMap["DY_Q1" ]->SetLineWidth(1); ThGraphMap["DY_Q1" ]->SetLineStyle(1); ThGraphMap["DY_Q1" ]->SetMarkerStyle(1); MuGraphMap["DY_Q1" ]->SetLineColor(46); MuGraphMap["DY_Q1" ]->SetMarkerColor(46); MuGraphMap["DY_Q1" ]->SetLineWidth(2); MuGraphMap["DY_Q1" ]->SetLineStyle(1); MuGraphMap["DY_Q1" ]->SetMarkerStyle(20); HQGraphMap["DY_Q1" ]->SetLineColor(46); HQGraphMap["DY_Q1" ]->SetMarkerColor(46); HQGraphMap["DY_Q1" ]->SetLineWidth(2); HQGraphMap["DY_Q1" ]->SetLineStyle(1); HQGraphMap["DY_Q1" ]->SetMarkerStyle(20); ThGraphMap["DY_Q2" ]->SetLineColor(2 ); ThGraphMap["DY_Q2" ]->SetMarkerColor(2 ); ThGraphMap["DY_Q2" ]->SetLineWidth(1); ThGraphMap["DY_Q2" ]->SetLineStyle(2); ThGraphMap["DY_Q2" ]->SetMarkerStyle(1); HQGraphMap["DY_Q2" ]->SetLineColor(2 ); HQGraphMap["DY_Q2" ]->SetMarkerColor(2 ); HQGraphMap["DY_Q2" ]->SetLineWidth(2); HQGraphMap["DY_Q2" ]->SetLineStyle(1); HQGraphMap["DY_Q2" ]->SetMarkerStyle(21); ThGraphMap["DY_Q3" ]->SetLineColor(1 ); ThGraphMap["DY_Q3" ]->SetMarkerColor(1 ); ThGraphMap["DY_Q3" ]->SetLineWidth(1); ThGraphMap["DY_Q3" ]->SetLineStyle(9); ThGraphMap["DY_Q3" ]->SetMarkerStyle(1); HQGraphMap["DY_Q3" ]->SetLineColor(1 ); HQGraphMap["DY_Q3" ]->SetMarkerColor(1 ); HQGraphMap["DY_Q3" ]->SetLineWidth(2); HQGraphMap["DY_Q3" ]->SetLineStyle(1); HQGraphMap["DY_Q3" ]->SetMarkerStyle(22); ThGraphMap["DY_Q4" ]->SetLineColor(6 ); ThGraphMap["DY_Q4" ]->SetMarkerColor(6 ); ThGraphMap["DY_Q4" ]->SetLineWidth(1); ThGraphMap["DY_Q4" ]->SetLineStyle(3); ThGraphMap["DY_Q4" ]->SetMarkerStyle(1); HQGraphMap["DY_Q4" ]->SetLineColor(6 ); HQGraphMap["DY_Q4" ]->SetMarkerColor(6 ); HQGraphMap["DY_Q4" ]->SetLineWidth(2); HQGraphMap["DY_Q4" ]->SetLineStyle(1); HQGraphMap["DY_Q4" ]->SetMarkerStyle(23); ThGraphMap["DY_Q5" ]->SetLineColor(4 ); ThGraphMap["DY_Q5" ]->SetMarkerColor(4 ); ThGraphMap["DY_Q5" ]->SetLineWidth(1); ThGraphMap["DY_Q5" ]->SetLineStyle(8); ThGraphMap["DY_Q5" ]->SetMarkerStyle(1); HQGraphMap["DY_Q5" ]->SetLineColor(4 ); HQGraphMap["DY_Q5" ]->SetMarkerColor(4 ); HQGraphMap["DY_Q5" ]->SetLineWidth(2); HQGraphMap["DY_Q5" ]->SetLineStyle(1); HQGraphMap["DY_Q5" ]->SetMarkerStyle(29); ThGraphMap["DY_Q6" ]->SetLineColor(9 ); ThGraphMap["DY_Q6" ]->SetMarkerColor(9 ); ThGraphMap["DY_Q6" ]->SetLineWidth(1); ThGraphMap["DY_Q6" ]->SetLineStyle(6); ThGraphMap["DY_Q6" ]->SetMarkerStyle(1); HQGraphMap["DY_Q6" ]->SetLineColor(9 ); HQGraphMap["DY_Q6" ]->SetMarkerColor(9 ); HQGraphMap["DY_Q6" ]->SetLineWidth(2); HQGraphMap["DY_Q6" ]->SetLineStyle(1); HQGraphMap["DY_Q6" ]->SetMarkerStyle(33); ThGraphMap["DY_Q7" ]->SetLineColor(12); ThGraphMap["DY_Q7" ]->SetMarkerColor(12); ThGraphMap["DY_Q7" ]->SetLineWidth(1); ThGraphMap["DY_Q7" ]->SetLineStyle(7); ThGraphMap["DY_Q7" ]->SetMarkerStyle(1); HQGraphMap["DY_Q7" ]->SetLineColor(12); HQGraphMap["DY_Q7" ]->SetMarkerColor(12); HQGraphMap["DY_Q7" ]->SetLineWidth(2); HQGraphMap["DY_Q7" ]->SetLineStyle(1); HQGraphMap["DY_Q7" ]->SetMarkerStyle(34); ThGraphMap["DY_Q8" ]->SetLineColor(14); ThGraphMap["DY_Q8" ]->SetMarkerColor(14); ThGraphMap["DY_Q8" ]->SetLineWidth(1); ThGraphMap["DY_Q8" ]->SetLineStyle(10); ThGraphMap["DY_Q8" ]->SetMarkerStyle(1); HQGraphMap["DY_Q8" ]->SetLineColor(14); HQGraphMap["DY_Q8" ]->SetMarkerColor(14); HQGraphMap["DY_Q8" ]->SetLineWidth(2); HQGraphMap["DY_Q8" ]->SetLineStyle(1); HQGraphMap["DY_Q8" ]->SetMarkerStyle(24); std::cout<<"TESTD\n"; ThGraphMap["GMStau" ]->SetMinimum(PlotMinScale); c1 = new TCanvas("c1", "c1",600,600); TMultiGraph* MGMu = new TMultiGraph(); if(!Combine) { MGMu->Add(ThGraphMap["Gluino_f10" ] ,"L"); MGMu->Add(ThGraphMap["Stop" ] ,"L"); MGMu->Add(ThGraphMap["GMStau" ] ,"L"); MGMu->Add(ThGraphMap["PPStau" ] ,"L"); MGMu->Add(ThGraphMap["DY_Q2o3" ] ,"L"); MGMu->Add(ThGraphMap["DY_Q1" ] ,"L"); } MGMu->Add(MuGraphMap["Gluino_f10" ] ,"LP"); MGMu->Add(MuGraphMap["Gluino_f50" ] ,"LP"); MGMu->Add(MuGraphMap["Stop" ] ,"LP"); MGMu->Add(MuGraphMap["GMStau" ] ,"LP"); MGMu->Add(MuGraphMap["PPStau" ] ,"LP"); MGMu->Add(MuGraphMap["DY_Q2o3" ] ,"LP"); MGMu->Add(MuGraphMap["DY_Q1" ] ,"LP"); MGMu->Draw("A"); if(!Combine) { ThErrorMap["Gluino_f10"]->Draw("f"); ThErrorMap["Stop" ]->Draw("f"); ThErrorMap["GMStau" ]->Draw("f"); ThErrorMap["PPStau" ]->Draw("f"); ThErrorMap["DY_Q2o3" ]->Draw("f"); ThErrorMap["DY_Q1" ]->Draw("f"); }else{ TLine* LineAtOne = new TLine(50,1,1550,1); LineAtOne->SetLineStyle(3); LineAtOne->Draw(); } MGMu->Draw("same"); MGMu->SetTitle(""); MGMu->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MGMu->GetYaxis()->SetTitle(Combine?"95% CL limit on #sigma/#sigma_{th}":"95% CL limit on #sigma (pb)"); MGMu->GetYaxis()->SetTitleOffset(1.40); MGMu->GetYaxis()->SetRangeUser(PlotMinScale,PlotMaxScale); MGMu->GetXaxis()->SetRangeUser(50,1550); DrawPreliminary(LegendFromType(MuPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS), false); TLegend* LEGMu = !Combine ? new TLegend(0.50,0.92-7*0.043,0.83,0.92) : new TLegend(0.50,0.15,0.83,0.15+7*0.043); LEGMu->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGMu->SetTextSize(18); //font size LEGMu->SetFillColor(0); LEGMu->SetFillStyle(0); LEGMu->SetBorderSize(0); LEGMu->AddEntry(MuGraphMap["Gluino_f50"] , "gluino; 50% #tilde{g}g" ,"LP"); LEGMu->AddEntry(MuGraphMap["Gluino_f10"] , "gluino; 10% #tilde{g}g" ,"LP"); LEGMu->AddEntry(MuGraphMap["Stop" ] , "stop" ,"LP"); LEGMu->AddEntry(MuGraphMap["PPStau" ] , "stau; dir. prod." ,"LP"); LEGMu->AddEntry(MuGraphMap["GMStau" ] , "stau" ,"LP"); LEGMu->AddEntry(MuGraphMap["DY_Q2o3" ], "|Q| = 2e/3" ,"LP"); LEGMu->AddEntry(MuGraphMap["DY_Q1" ], "|Q| = 1e" ,"LP"); TLegend* LEGTh = new TLegend(0.15,0.92-(1+6)*0.043,0.50,0.92); LEGTh->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGTh->SetTextSize(18); //font size if(!Combine) { LEGTh->SetHeader("Theoretical Prediction"); LEGTh->SetFillColor(0); //LEGTh->SetFillStyle(0); LEGTh->SetBorderSize(0); TGraph* GlThLeg = (TGraph*) ThGraphMap["Gluino_f10"]->Clone("GluinoThLeg"); GlThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGTh->AddEntry(GlThLeg, "gluino (NLO+NLL)" ,"LF"); TGraph* StThLeg = (TGraph*) ThGraphMap["Stop" ]->Clone("StopThLeg"); StThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGTh->AddEntry(StThLeg ,"stop (NLO+NLL)" ,"LF"); TGraph* PPStauThLeg = (TGraph*) ThGraphMap["PPStau" ]->Clone("PPStauThLeg"); PPStauThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGTh->AddEntry(PPStauThLeg ,"stau, dir. prod. (NLO)" ,"LF"); TGraph* StauThLeg = (TGraph*) ThGraphMap["GMStau" ]->Clone("StauThLeg"); StauThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGTh->AddEntry(StauThLeg ,"stau (NLO)" ,"LF"); TGraph* DYQ2o3ThLeg = (TGraph*) ThGraphMap["DY_Q2o3" ]->Clone("DYQ2o3ThLeg"); DYQ2o3ThLeg->SetFillColor(ThErrorMap["DY_Q2o3"]->GetFillColor()); LEGTh->AddEntry(DYQ2o3ThLeg ,"|Q| = 2e/3 (LO)" ,"LF"); TGraph* DYQ1ThLeg = (TGraph*) ThGraphMap["DY_Q1" ]->Clone("DYQ1ThLeg"); DYQ1ThLeg->SetFillColor(ThErrorMap["DY_Q1"]->GetFillColor()); LEGTh->AddEntry(DYQ1ThLeg ,"|Q| = 1e (LO)" ,"LF"); LEGTh->Draw(); } LEGMu->Draw(); c1->SetLogy(true); SaveCanvas(c1, outpath, string("MuExclusionLog")); delete c1; c1 = new TCanvas("c1", "c1",600,600); TMultiGraph* MGTk = new TMultiGraph(); if(!Combine) { MGTk->Add(ThGraphMap["Gluino_f10" ] ,"L"); MGTk->Add(ThGraphMap["Stop" ] ,"L"); MGTk->Add(ThGraphMap["GMStau" ] ,"L"); MGTk->Add(ThGraphMap["PPStau" ] ,"L"); MGTk->Add(ThGraphMap["DY_Q2o3" ] ,"L"); } MGTk->Add(TkGraphMap["Gluino_f10" ] ,"LP"); MGTk->Add(TkGraphMap["Gluino_f50" ] ,"LP"); MGTk->Add(TkGraphMap["GluinoN_f10"] ,"LP"); MGTk->Add(TkGraphMap["Stop" ] ,"LP"); MGTk->Add(TkGraphMap["StopN" ] ,"LP"); MGTk->Add(TkGraphMap["GMStau" ] ,"LP"); MGTk->Add(TkGraphMap["PPStau" ] ,"LP"); MGTk->Add(TkGraphMap["DY_Q2o3" ] ,"LP"); MGTk->Draw("A"); if(!Combine) { ThErrorMap["Gluino_f10"]->Draw("f"); ThErrorMap["Stop" ]->Draw("f"); ThErrorMap["GMStau" ]->Draw("f"); ThErrorMap["PPStau" ]->Draw("f"); ThErrorMap["DY_Q2o3" ]->Draw("f"); }else{ TLine* LineAtOne = new TLine(50,1,1550,1); LineAtOne->SetLineStyle(3); LineAtOne->Draw(); } MGTk->Draw("same"); MGTk->SetTitle(""); MGTk->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MGTk->GetYaxis()->SetTitle(Combine?"95% CL limit on #sigma/#sigma_{th}":"95% CL limit on #sigma (pb)"); MGTk->GetYaxis()->SetTitleOffset(1.40); MGTk->GetYaxis()->SetRangeUser(PlotMinScale,PlotMaxScale); //if(Combine) MGTk->GetYaxis()->SetRangeUser(PlotMinScale,50); //else MGTk->GetYaxis()->SetRangeUser(PlotMinScale,700); MGTk->GetXaxis()->SetRangeUser(50,1550); DrawPreliminary(LegendFromType(TkPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS), false); TLegend* LEGTk = !Combine ? new TLegend(0.50,0.92-8*0.043,0.83,0.92) : new TLegend(0.50,0.15,0.83,0.15+8*0.043); LEGTk->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGTk->SetTextSize(18); //font size LEGTk->SetFillColor(0); LEGTk->SetFillStyle(0); LEGTk->SetBorderSize(0); LEGTk->AddEntry(TkGraphMap["Gluino_f50" ], "gluino; 50% #tilde{g}g" ,"LP"); LEGTk->AddEntry(TkGraphMap["Gluino_f10" ], "gluino; 10% #tilde{g}g" ,"LP"); LEGTk->AddEntry(TkGraphMap["GluinoN_f10"], "gluino; 10% #tilde{g}g; CS" ,"LP"); LEGTk->AddEntry(TkGraphMap["Stop" ], "stop" ,"LP"); LEGTk->AddEntry(TkGraphMap["StopN" ], "stop; CS" ,"LP"); LEGTk->AddEntry(TkGraphMap["PPStau" ], "stau; dir. prod." ,"LP"); LEGTk->AddEntry(TkGraphMap["GMStau" ], "stau" ,"LP"); LEGTk->AddEntry(TkGraphMap["DY_Q2o3" ], "|Q| = 2e/3" ,"LP"); TLegend* LEGThTk = new TLegend(0.15,0.92-(1+6)*0.043,0.50,0.92); LEGThTk->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGThTk->SetTextSize(18); //font size if(!Combine) { LEGThTk->SetHeader("Theoretical Prediction"); LEGThTk->SetFillColor(0); //LEGThTk->SetFillStyle(0); LEGThTk->SetBorderSize(0); TGraph* GlThLeg = (TGraph*) ThGraphMap["Gluino_f10"]->Clone("GluinoThLeg"); GlThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGThTk->AddEntry(GlThLeg, "gluino (NLO+NLL)" ,"LF"); TGraph* StThLeg = (TGraph*) ThGraphMap["Stop" ]->Clone("StopThLeg"); StThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGThTk->AddEntry(StThLeg ,"stop (NLO+NLL)" ,"LF"); TGraph* PPStauThLeg = (TGraph*) ThGraphMap["PPStau" ]->Clone("PPStauThLeg"); PPStauThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGThTk->AddEntry(PPStauThLeg ,"stau; dir. prod. (NLO)" ,"LF"); TGraph* StauThLeg = (TGraph*) ThGraphMap["GMStau" ]->Clone("StauThLeg"); StauThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGThTk->AddEntry(StauThLeg ,"stau (NLO)" ,"LF"); TGraph* DYQ2o3ThLeg = (TGraph*) ThGraphMap["DY_Q2o3" ]->Clone("DYQ2o3ThLeg"); DYQ2o3ThLeg->SetFillColor(ThErrorMap["DY_Q2o3"]->GetFillColor()); LEGThTk->AddEntry(DYQ2o3ThLeg ,"|Q| = 2e/3 (LO)" ,"LF"); //TGraph* DYQ1ThLeg = (TGraph*) ThGraphMap["DY_Q1" ]->Clone("DYQ1ThLeg"); //DYQ1ThLeg->SetFillColor(ThErrorMap["DY_Q1"]->GetFillColor()); //LEGThTk->AddEntry(DYQ2o3ThLeg ,"|Q| = 1e (LO)" ,"LF"); LEGThTk->Draw(); } if(!Combine) LEGThTk->Draw(); LEGTk->Draw(); c1->SetLogy(true); SaveCanvas(c1, outpath, string("TkExclusionLog")); delete c1; ////////////////////////////////////// plot requested by G. Landsberg if(Combine){ c1 = new TCanvas("c1", "c1",600,600); MGTk = new TMultiGraph(); MuGraphMap["GMStau" ]->SetLineColor(2); MuGraphMap["GMStau" ]->SetMarkerColor(2); MuGraphMap["GMStau" ]->SetLineWidth(2); MuGraphMap["GMStau" ]->SetLineStyle(1); MuGraphMap["GMStau" ]->SetMarkerStyle(22); TkGraphMap["GMStau" ]->SetLineColor(1); TkGraphMap["GMStau" ]->SetMarkerColor(1); TkGraphMap["GMStau" ]->SetLineWidth(2); TkGraphMap["GMStau" ]->SetLineStyle(1); TkGraphMap["GMStau" ]->SetMarkerStyle(23); MuGraphMap["PPStau" ]->SetLineColor(2); MuGraphMap["PPStau" ]->SetMarkerColor(2); MuGraphMap["PPStau" ]->SetLineWidth(2); MuGraphMap["PPStau" ]->SetLineStyle(2); MuGraphMap["PPStau" ]->SetMarkerStyle(26); TkGraphMap["PPStau" ]->SetLineColor(1); TkGraphMap["PPStau" ]->SetMarkerColor(1); TkGraphMap["PPStau" ]->SetLineWidth(2); TkGraphMap["PPStau" ]->SetLineStyle(2); TkGraphMap["PPStau" ]->SetMarkerStyle(32); MGTk->Add(MuGraphMap["PPStau" ] ,"LP"); MGTk->Add(TkGraphMap["PPStau" ] ,"LP"); MGTk->Add(MuGraphMap["GMStau" ] ,"LP"); MGTk->Add(TkGraphMap["GMStau" ] ,"LP"); MGTk->Draw("A"); TLine* LineAtOne = new TLine(90,1,570,1); LineAtOne->SetLineStyle(3); LineAtOne->Draw(); //60.6,533.4 MGTk->Draw("same"); MGTk->SetTitle(""); MGTk->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MGTk->GetYaxis()->SetTitle(Combine?"95% CL limit on #sigma/#sigma_{th}":"95% CL limit on #sigma (pb)"); MGTk->GetYaxis()->SetTitleOffset(1.40); MGTk->GetYaxis()->SetRangeUser(1e-3,20); //if(Combine) MGTk->GetYaxis()->SetRangeUser(PlotMinScale,50); //else MGTk->GetYaxis()->SetRangeUser(PlotMinScale,700); MGTk->GetXaxis()->SetRangeUser(90,570); DrawPreliminary("", SQRTS, IntegratedLuminosityFromE(SQRTS), false); LEGTk = !Combine ? new TLegend(0.50,0.92-3*0.043,0.83,0.92) : new TLegend(0.45,0.15+4*0.043,0.80,0.15+7*0.043); LEGTk->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGTk->SetTextSize(18); //font size LEGTk->SetFillColor(0); LEGTk->SetFillStyle(0); LEGTk->SetBorderSize(0); LEGTk->SetHeader("Stau prod. (direct+indirect)"); LEGTk->AddEntry(MuGraphMap["GMStau" ], LegendFromType(MuPattern).c_str() ,"LP"); LEGTk->AddEntry(TkGraphMap["GMStau" ], LegendFromType(TkPattern).c_str() ,"LP"); LEGTk->Draw(); LEGTh = !Combine ? new TLegend(0.50,0.92-3*0.043,0.83,0.92) : new TLegend(0.45,0.15+0*0.043,0.80,0.15+3*0.043); LEGTh->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGTh->SetTextSize(18); //font size LEGTh->SetFillColor(0); LEGTh->SetFillStyle(0); LEGTh->SetBorderSize(0); LEGTh->SetHeader("Stau prod. (direct)"); LEGTh->AddEntry(MuGraphMap["PPStau" ], LegendFromType(MuPattern).c_str() ,"LP"); LEGTh->AddEntry(TkGraphMap["PPStau" ], LegendFromType(TkPattern).c_str() ,"LP"); LEGTh->Draw(); c1->SetLogy(true); SaveCanvas(c1, outpath, string("StauExclusionLog")); delete c1; } ////////////////////////////////////// /* c1 = new TCanvas("c1", "c1",600,600); TMultiGraph* MGDCMu = new TMultiGraph(); if(!Combine) { // MGDCMu->Add(ThGraphMap["DCRho08HyperK"] ,"L"); //MGDCMu->Add(ThGraphMap["DCRho12HyperK"] ,"L"); //MGDCMu->Add(ThGraphMap["DCRho16HyperK"] ,"L"); } // MGDCMu->Add(MuGraphMap["DCRho08HyperK"] ,"LP"); //MGDCMu->Add(MuGraphMap["DCRho12HyperK"] ,"LP"); //MGDCMu->Add(MuGraphMap["DCRho16HyperK"] ,"LP"); MGDCMu->Draw("A"); MGDCMu->Draw("same"); MGDCMu->SetTitle(""); MGDCMu->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MGDCMu->GetYaxis()->SetTitle(Combine?"95% CL limit on #sigma/#sigma_{th}":"95% CL limit on #sigma (pb)"); MGDCMu->GetYaxis()->SetTitleOffset(1.70); MGDCMu->GetYaxis()->SetRangeUser(PlotMinScale,PlotMaxScale); MGDCMu->GetXaxis()->SetRangeUser(50,1550); DrawPreliminary("Tracker + TOF", SQRTS, LInt); TLegend* LEGDCMu = new TLegend(0.50,0.65,0.80,0.9); LEGDCMu->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGDCMu->SetTextSize(18); //font size // LEGDCMu->SetHeader("Tracker + TOF"); LEGDCMu->SetFillColor(0); LEGDCMu->SetFillStyle(0); LEGDCMu->SetBorderSize(0); // LEGDCMu->AddEntry(MuGraphMap["DCRho08HyperK"] , "Hyper-K, #tilde{#rho} = 0.8 TeV" ,"LP"); LEGDCMu->AddEntry(MuGraphMap["DCRho12HyperK"] , "Hyper-K, #tilde{#rho} = 1.2 TeV" ,"LP"); LEGDCMu->AddEntry(MuGraphMap["DCRho16HyperK"] , "Hyper-K, #tilde{#rho} = 1.6 TeV" ,"LP"); TLegend* LEGDCTh = new TLegend(0.15,0.7,0.49,0.9); LEGDCTh->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGDCTh->SetTextSize(18); //font size if(!Combine) { LEGDCTh->SetHeader("Theoretical Prediction"); LEGDCTh->SetFillColor(0); LEGDCTh->SetFillStyle(0); LEGDCTh->SetBorderSize(0); // TGraph* DCRho08HyperKThLeg = (TGraph*) ThGraphMap["DCRho08HyperK"]->Clone("DCRho08HyperKThLeg"); // DCRho08HyperKThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); // LEGDCTh->AddEntry(DCRho08HyperKThLeg ,"Hyper-K, #tilde{#rho} = 0.8 TeV (LO)" ,"L"); TGraph* DCRho12HyperKThLeg = (TGraph*) ThGraphMap["DCRho12HyperK"]->Clone("DCRho12HyperKThLeg"); DCRho12HyperKThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGDCTh->AddEntry(DCRho12HyperKThLeg ,"Hyper-K, #tilde{#rho} = 1.2 TeV (LO)" ,"L"); TGraph* DCRho16HyperKThLeg = (TGraph*) ThGraphMap["DCRho16HyperK"]->Clone("DCRho16HyperKThLeg"); DCRho16HyperKThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); LEGDCTh->AddEntry(DCRho16HyperKThLeg ,"Hyper-K, #tilde{#rho} = 1.6 TeV (LO)" ,"L"); LEGDCTh->Draw(); } LEGDCMu->Draw(); c1->SetLogy(true); SaveCanvas(c1, outpath, string("MuDCExclusionLog")); delete c1; c1 = new TCanvas("c1", "c1",600,600); TMultiGraph* MGDCTk = new TMultiGraph(); if(!Combine) { // MGDCTk->Add(ThGraphMap["DCRho08HyperK"] ,"L"); MGDCTk->Add(ThGraphMap["DCRho12HyperK"] ,"L"); MGDCTk->Add(ThGraphMap["DCRho16HyperK"] ,"L"); } // MGDCTk->Add(TkGraphMap["DCRho08HyperK"] ,"LP"); MGDCTk->Add(TkGraphMap["DCRho12HyperK"] ,"LP"); MGDCTk->Add(TkGraphMap["DCRho16HyperK"] ,"LP"); MGDCTk->Draw("A"); MGDCTk->Draw("same"); MGDCTk->SetTitle(""); MGDCTk->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MGDCTk->GetYaxis()->SetTitle(Combine?"95% CL limit on #sigma/#sigma_{th}":"95% CL limit on #sigma (pb)"); MGDCTk->GetYaxis()->SetTitleOffset(1.70); MGDCTk->GetYaxis()->SetRangeUser(PlotMinScale,PlotMaxScale); MGDCTk->GetXaxis()->SetRangeUser(50,1550); DrawPreliminary("Tracker - Only", SQRTS, LInt); TLegend* LEGDCTk = new TLegend(0.50,0.65,0.80,0.90); LEGDCTk->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGDCTk->SetTextSize(18); //font size // LEGDCTk->SetHeader("Tracker - Only"); LEGDCTk->SetFillColor(0); LEGDCTk->SetFillStyle(0); LEGDCTk->SetBorderSize(0); LEGDCTk->AddEntry(TkGraphMap["DCRho08HyperK"] , "Hyper-K, #tilde{#rho} = 0.8 TeV" ,"LP"); // LEGDCTk->AddEntry(TkGraphMap["DCRho12HyperK"] , "Hyper-K, #tilde{#rho} = 1.2 TeV" ,"LP"); LEGDCTk->AddEntry(TkGraphMap["DCRho16HyperK"] , "Hyper-K, #tilde{#rho} = 1.6 TeV" ,"LP"); LEGDCTk->Draw(); LEGDCTh->Draw(); c1->SetLogy(true); SaveCanvas(c1, outpath, string("TkDCExclusionLog")); delete c1; */ c1 = new TCanvas("c1", "c1",600,600); TMultiGraph* MGMO = new TMultiGraph(); if(!Combine) { MGMO->Add(ThGraphMap["Gluino_f10" ] ,"L"); MGMO->Add(ThGraphMap["Stop" ] ,"L"); //MGMO->Add(ThGraphMap["DY_Q1o3" ] ,"L"); //MGMO->Add(ThGraphMap["DY_Q2o3" ] ,"L"); } //TGraph* DYQ1o3ThLeg = (TGraph*) ThGraphMap["DY_Q1o3" ]->Clone("DYQ1o3ThLeg"); //DYQ1o3ThLeg->SetFillColor(ThErrorMap["DY_Q1o3"]->GetFillColor()); //LQLEGTh->AddEntry(DYQ1o3ThLeg ,"|Q| = e/3 (LO)" ,"LF"); //TGraph* DYQ2o3ThLeg = (TGraph*) ThGraphMap["DY_Q2o3" ]->Clone("DYQ2o3ThLeg"); //DYQ2o3ThLeg->SetFillColor(ThErrorMap["DY_Q2o3"]->GetFillColor()); //LQLEGTh->AddEntry(DYQ2o3ThLeg ,"|Q| = 2e/3 (LO)" ,"LF"); //LQLEGTh->Draw(); MGMO->Add(MOGraphMap["Gluino_f10" ] ,"LP"); MGMO->Add(MOGraphMap["Gluino_f50" ] ,"LP"); MGMO->Add(MOGraphMap["Gluino_f100"] ,"LP"); MGMO->Add(MOGraphMap["Stop" ] ,"LP"); //MGMO->Add(MOGraphMap["DY_Q1o3" ] ,"LP"); //MGMO->Add(MOGraphMap["DY_Q2o3" ] ,"LP"); MGMO->Draw("A"); if(!Combine) { ThErrorMap["Gluino_f10"]->Draw("f"); ThErrorMap["Stop" ]->Draw("f"); //ThErrorMap["DY_Q1o3" ]->Draw("f"); //ThErrorMap["DY_Q2o3" ]->Draw("f"); }else{ TLine* LineAtOne = new TLine(50,1,1550,1); LineAtOne->SetLineStyle(3); LineAtOne->Draw(); } MGMO->Draw("same"); MGMO->SetTitle(""); MGMO->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MGMO->GetYaxis()->SetTitle(Combine?"95% CL limit on #sigma/#sigma_{th}":"95% CL limit on #sigma (pb)"); MGMO->GetYaxis()->SetTitleOffset(1.40); MGMO->GetYaxis()->SetRangeUser(PlotMinScale,PlotMaxScale); MGMO->GetXaxis()->SetRangeUser(50,1550); DrawPreliminary(LegendFromType(MOPattern).c_str(), 8.0, IntegratedLuminosityFromE(8.0), false); TLegend* LEGMO = !Combine ? new TLegend(0.50,0.92-4*0.043,0.83,0.92) : new TLegend(0.55,0.25,0.80,0.25+4*0.043); LEGMO->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGMO->SetTextSize(18); //font size LEGMO->SetFillColor(0); LEGMO->SetFillStyle(0); LEGMO->SetBorderSize(0); LEGMO->AddEntry(MOGraphMap["Gluino_f100" ], "gluino; 100% #tilde{g}g" ,"LP"); LEGMO->AddEntry(MOGraphMap["Gluino_f50" ], "gluino; 50% #tilde{g}g" ,"LP"); LEGMO->AddEntry(MOGraphMap["Gluino_f10" ], "gluino; 10% #tilde{g}g" ,"LP"); LEGMO->AddEntry(MOGraphMap["Stop" ], "stop" ,"LP"); //LEGMO->AddEntry(TkGraphMap["DY_Q1o3" ], "|Q| = e/3" ,"LP"); //LEGMO->AddEntry(TkGraphMap["DY_Q2o3" ], "|Q| = 2e/3" ,"LP"); LEGMO->Draw(); TLegend* MOLEGTh = new TLegend(0.15,0.92-(1+2)*0.043,0.50,0.92); MOLEGTh->SetTextFont(43); //give the font size in pixel (instead of fraction) MOLEGTh->SetTextSize(18); //font size if(!Combine) { MOLEGTh->SetHeader("Theoretical Prediction"); MOLEGTh->SetFillColor(0); //MOLEGTh->SetFillStyle(0); MOLEGTh->SetBorderSize(0); TGraph* GlThLeg = (TGraph*) ThGraphMap["Gluino_f10"]->Clone("GluinoThLeg"); GlThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); MOLEGTh->AddEntry(GlThLeg, "gluino (NLO+NLL)" ,"LF"); TGraph* StThLeg = (TGraph*) ThGraphMap["Stop" ]->Clone("StopThLeg"); StThLeg->SetFillColor(ThErrorMap["Gluino_f10"]->GetFillColor()); MOLEGTh->AddEntry(StThLeg ,"stop (NLO+NLL)" ,"LF"); //TGraph* DYQ1o3ThLeg = (TGraph*) ThGraphMap["DY_Q1o3" ]->Clone("DYQ1o3ThLeg"); //DYQ1o3ThLeg->SetFillColor(ThErrorMap["DY_Q1o3"]->GetFillColor()); //MOLEGTh->AddEntry(DYQ1o3ThLeg ,"|Q| = e/3 (LO)" ,"LF"); //TGraph* DYQ2o3ThLeg = (TGraph*) ThGraphMap["DY_Q2o3" ]->Clone("DYQ2o3ThLeg"); //DYQ2o3ThLeg->SetFillColor(ThErrorMap["DY_Q2o3"]->GetFillColor()); //MOLEGTh->AddEntry(DYQ2o3ThLeg ,"|Q| = 2e/3 (LO)" ,"LF"); MOLEGTh->Draw(); } c1->SetLogy(true); SaveCanvas(c1, outpath, string("MOExclusionLog")); delete c1; /////////////////////////////// LQ Analysis TLegend* LQLEGTh = new TLegend(0.15,0.92-(1+2)*0.043,0.50,0.92); LQLEGTh->SetTextFont(43); //give the font size in pixel (instead of fraction) LQLEGTh->SetTextSize(18); //font size c1 = new TCanvas("c1", "c1",600,600); if(!Combine) { LQLEGTh->SetHeader("Theoretical Prediction"); LQLEGTh->SetFillColor(0); //LQLEGTh->SetFillStyle(0); LQLEGTh->SetBorderSize(0); TGraph* DYQ1o3ThLeg = (TGraph*) ThGraphMap["DY_Q1o3" ]->Clone("DYQ1o3ThLeg"); DYQ1o3ThLeg->SetFillColor(ThErrorMap["DY_Q1o3"]->GetFillColor()); LQLEGTh->AddEntry(DYQ1o3ThLeg ,"|Q| = e/3 (LO)" ,"LF"); TGraph* DYQ2o3ThLeg = (TGraph*) ThGraphMap["DY_Q2o3" ]->Clone("DYQ2o3ThLeg"); DYQ2o3ThLeg->SetFillColor(ThErrorMap["DY_Q2o3"]->GetFillColor()); LQLEGTh->AddEntry(DYQ2o3ThLeg ,"|Q| = 2e/3 (LO)" ,"LF"); LQLEGTh->Draw(); } TMultiGraph* MGLQ = new TMultiGraph(); if(!Combine) { MGLQ->Add(ThGraphMap["DY_Q1o3" ] ,"L"); MGLQ->Add(ThGraphMap["DY_Q2o3" ] ,"L"); } MGLQ->Add(LQGraphMap["DY_Q1o3" ] ,"LP"); MGLQ->Add(LQGraphMap["DY_Q2o3" ] ,"LP"); MGLQ->Draw("A"); if(!Combine) { ThErrorMap["DY_Q1o3" ]->Draw("f"); ThErrorMap["DY_Q2o3" ]->Draw("f"); }else{ TLine* LineAtOne = new TLine(75,1,625,1); LineAtOne->SetLineStyle(3); LineAtOne->Draw(); } MGLQ->Draw("same"); MGLQ->SetTitle(""); MGLQ->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MGLQ->GetYaxis()->SetTitle(Combine?"95% CL limit on #sigma/#sigma_{th}":"95% CL limit on #sigma (pb)"); MGLQ->GetYaxis()->SetTitleOffset(1.40); MGLQ->GetYaxis()->SetRangeUser(PlotMinScale,PlotMaxScale); MGLQ->GetXaxis()->SetRangeUser(75,625); DrawPreliminary(LegendFromType(LQPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS), false); TLegend* LEGLQ = !Combine ? new TLegend(0.50,0.92-2*0.043,0.83,0.92) : new TLegend(0.20,0.88-2*0.043,0.50,0.88); LEGLQ->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGLQ->SetTextSize(18); //font size // LEGLQ->SetHeader("Q<1"); LEGLQ->SetFillColor(0); LEGLQ->SetFillStyle(0); LEGLQ->SetBorderSize(0); LEGLQ->AddEntry(TkGraphMap["DY_Q1o3" ], "|Q| = e/3" ,"LP"); LEGLQ->AddEntry(TkGraphMap["DY_Q2o3" ], "|Q| = 2e/3" ,"LP"); if(!Combine) LQLEGTh->Draw(); LEGLQ->Draw(); c1->SetLogy(true); SaveCanvas(c1, outpath, string("LQExclusionLog")); delete c1; c1 = new TCanvas("c1", "c1",600,600); TMultiGraph* MGHQ = new TMultiGraph(); if(!Combine) { MGHQ->Add(ThGraphMap["DY_Q1" ] ,"L"); MGHQ->Add(ThGraphMap["DY_Q2" ] ,"L"); MGHQ->Add(ThGraphMap["DY_Q3" ] ,"L"); MGHQ->Add(ThGraphMap["DY_Q4" ] ,"L"); MGHQ->Add(ThGraphMap["DY_Q5" ] ,"L"); MGHQ->Add(ThGraphMap["DY_Q6" ] ,"L"); MGHQ->Add(ThGraphMap["DY_Q7" ] ,"L"); MGHQ->Add(ThGraphMap["DY_Q8" ] ,"L"); } MGHQ->Add(HQGraphMap["DY_Q1" ] ,"LP"); MGHQ->Add(HQGraphMap["DY_Q2" ] ,"LP"); MGHQ->Add(HQGraphMap["DY_Q3" ] ,"LP"); MGHQ->Add(HQGraphMap["DY_Q4" ] ,"LP"); MGHQ->Add(HQGraphMap["DY_Q5" ] ,"LP"); MGHQ->Add(HQGraphMap["DY_Q6" ] ,"LP"); MGHQ->Add(HQGraphMap["DY_Q7" ] ,"LP"); MGHQ->Add(HQGraphMap["DY_Q8" ] ,"LP"); MGHQ->Draw("A"); if(!Combine) { ThErrorMap["DY_Q1"]->Draw("f"); ThErrorMap["DY_Q2"]->Draw("f"); ThErrorMap["DY_Q3"]->Draw("f"); ThErrorMap["DY_Q4"]->Draw("f"); ThErrorMap["DY_Q5"]->Draw("f"); ThErrorMap["DY_Q6"]->Draw("f"); ThErrorMap["DY_Q7"]->Draw("f"); ThErrorMap["DY_Q8"]->Draw("f"); }else{ TLine* LineAtOne = new TLine(50,1,1050,1); LineAtOne->SetLineStyle(3); LineAtOne->Draw(); } MGHQ->Draw("same"); MGHQ->SetTitle(""); MGHQ->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MGHQ->GetYaxis()->SetTitle(Combine?"95% CL limit on #sigma/#sigma_{th}":"95% CL limit on #sigma (pb)"); MGHQ->GetYaxis()->SetTitleOffset(1.40); MGHQ->GetYaxis()->SetRangeUser(PlotMinScale,PlotMaxScale); //MGHQ->GetYaxis()->SetRangeUser(PlotMinScale,100); MGHQ->GetXaxis()->SetRangeUser(50,1050); DrawPreliminary(LegendFromType(HQPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS), false); TLegend* LEGHQ = !Combine ? new TLegend(0.62,0.92-0.043-8*0.043,0.83,0.92-0.043) : new TLegend(0.55,0.35,0.80,0.35+6*0.043); // TLegend* LEGHQ = !Combine ? new TLegend(0.62,0.92-5*0.043,0.83,0.92) : new TLegend(0.55,0.35,0.80,0.35+6*0.043); LEGHQ->SetTextFont(43); //give the font size in pixel (instead of fraction) LEGHQ->SetTextSize(18); //font size LEGHQ->SetFillColor(0); //LEGHQ->SetFillStyle(0); LEGHQ->SetBorderSize(0); LEGHQ->AddEntry(HQGraphMap["DY_Q1"] , "|Q| = 1e " ,"LP"); LEGHQ->AddEntry(HQGraphMap["DY_Q2"] , "|Q| = 2e " ,"LP"); LEGHQ->AddEntry(HQGraphMap["DY_Q3"] , "|Q| = 3e " ,"LP"); LEGHQ->AddEntry(HQGraphMap["DY_Q4"] , "|Q| = 4e " ,"LP"); LEGHQ->AddEntry(HQGraphMap["DY_Q5"] , "|Q| = 5e " ,"LP"); LEGHQ->AddEntry(HQGraphMap["DY_Q6"] , "|Q| = 6e " ,"LP"); LEGHQ->AddEntry(HQGraphMap["DY_Q7"] , "|Q| = 7e " ,"LP"); LEGHQ->AddEntry(HQGraphMap["DY_Q8"] , "|Q| = 8e " ,"LP"); TLegend* HQLEGTh = new TLegend(0.35,0.92-(1+8)*0.043,0.57,0.92); // TLegend* HQLEGTh = new TLegend(0.3,0.92-(1+5)*0.043,0.57,0.92); HQLEGTh->SetTextFont(43); //give the font size in pixel (instead of fraction) HQLEGTh->SetTextSize(18); //font size if(!Combine){ HQLEGTh->SetHeader("Theoretical Prediction"); HQLEGTh->SetFillColor(0); HQLEGTh->SetFillStyle(0); HQLEGTh->SetBorderSize(0); TGraph* Q1ThLeg = (TGraph*) ThGraphMap["DY_Q1"]->Clone("HSCPQ1ThLeg"); Q1ThLeg->SetFillColor(ThErrorMap["DY_Q1"]->GetFillColor()); HQLEGTh->AddEntry(Q1ThLeg, "|Q| = 1e (LO)" ,"LF"); TGraph* Q2ThLeg = (TGraph*) ThGraphMap["DY_Q2"]->Clone("HSCPQ2ThLeg"); Q2ThLeg->SetFillColor(ThErrorMap["DY_Q2"]->GetFillColor()); HQLEGTh->AddEntry(Q2ThLeg, "|Q| = 2e (LO)" ,"LF"); TGraph* Q3ThLeg = (TGraph*) ThGraphMap["DY_Q3"]->Clone("HSCPQ3ThLeg"); Q3ThLeg->SetFillColor(ThErrorMap["DY_Q3"]->GetFillColor()); HQLEGTh->AddEntry(Q3ThLeg, "|Q| = 3e (LO)" ,"LF"); TGraph* Q4ThLeg = (TGraph*) ThGraphMap["DY_Q4"]->Clone("HSCPQ4ThLeg"); Q4ThLeg->SetFillColor(ThErrorMap["DY_Q4"]->GetFillColor()); HQLEGTh->AddEntry(Q4ThLeg, "|Q| = 4e (LO)" ,"LF"); TGraph* Q5ThLeg = (TGraph*) ThGraphMap["DY_Q5"]->Clone("HSCPQ5ThLeg"); Q5ThLeg->SetFillColor(ThErrorMap["DY_Q5"]->GetFillColor()); HQLEGTh->AddEntry(Q5ThLeg, "|Q| = 5e (LO)" ,"LF"); TGraph* Q6ThLeg = (TGraph*) ThGraphMap["DY_Q6"]->Clone("HSCPQ6ThLeg"); Q6ThLeg->SetFillColor(ThErrorMap["DY_Q6"]->GetFillColor()); HQLEGTh->AddEntry(Q6ThLeg, "|Q| = 6e (LO)" ,"LF"); TGraph* Q7ThLeg = (TGraph*) ThGraphMap["DY_Q7"]->Clone("HSCPQ7ThLeg"); Q7ThLeg->SetFillColor(ThErrorMap["DY_Q7"]->GetFillColor()); HQLEGTh->AddEntry(Q7ThLeg, "|Q| = 7e (LO)" ,"LF"); TGraph* Q8ThLeg = (TGraph*) ThGraphMap["DY_Q8"]->Clone("HSCPQ8ThLeg"); Q8ThLeg->SetFillColor(ThErrorMap["DY_Q8"]->GetFillColor()); HQLEGTh->AddEntry(Q8ThLeg, "|Q| = 8e (LO)" ,"LF"); HQLEGTh->Draw(); } LEGHQ->Draw(); c1->SetLogy(true); SaveCanvas(c1, outpath, string("HQExclusionLog")); delete c1; return; } TGraph* CheckSignalUncertainty(FILE* pFile, FILE* talkFile, string InputPattern, string ModelName, std::vector<stSample>& modelSample){ int TypeMode = TypeFromPattern(InputPattern); string prefix = "BUG"; switch(TypeMode){ case 0: prefix = "Tk"; break; case 2: prefix = "Mu"; break; case 3: prefix = "Mo"; break; case 4: prefix = "HQ"; break; case 5: prefix = "LQ"; break; } unsigned int N = 0; double* Mass = new double [modelSample.size()]; double* MassErr = new double [modelSample.size()]; double* SystP = new double [modelSample.size()]; double* SystErrP = new double [modelSample.size()]; double* SystI = new double [modelSample.size()]; double* SystErrI = new double [modelSample.size()]; double* SystPU = new double [modelSample.size()]; double* SystErrPU = new double [modelSample.size()]; double* SystT = new double [modelSample.size()]; double* SystErrT = new double [modelSample.size()]; double* SystTr = new double [modelSample.size()]; double* SystErrTr = new double [modelSample.size()]; double* SystRe = new double [modelSample.size()]; double* SystErrRe = new double [modelSample.size()]; double* SystMB = new double [modelSample.size()]; double* SystErrMB = new double [modelSample.size()]; double* SystTotal = new double [modelSample.size()]; double* SystErrTotal = new double [modelSample.size()]; double* SystTotal2 = new double [modelSample.size()]; for(unsigned int s=0;s<modelSample.size();s++){ if(modelSample[s].Type!=2)continue; bool IsNeutral = (modelSample[s].ModelName().find("N")!=std::string::npos); if(TypeMode!=0 && IsNeutral)continue; stAllInfo tmp(InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR + "/"+modelSample[s].Name+".txt"); if(tmp.Eff==0) continue; Mass[N] = tmp.Mass; MassErr[N] = 0.0; SystP[N] = (tmp.Eff_SYSTP - tmp.Eff)/tmp.Eff; SystErrP[N] = (tmp.EffE_SYSTP)/tmp.Eff; SystI[N] = (tmp.Eff_SYSTI - tmp.Eff)/tmp.Eff; SystErrI[N] = (tmp.EffE_SYSTI)/tmp.Eff; if(TypeMode==5){ if(modelSample[s].ModelName().find("1o3")!=string::npos) SystI[N]=-0.25; SystErrP[N] = 0; if(modelSample[s].ModelName().find("2o3")!=string::npos) SystI[N]=-0.10; SystErrI[N] = 0; } SystPU[N] = (tmp.Eff_SYSTPU - tmp.Eff)/tmp.Eff; SystErrPU[N] = (tmp.EffE_SYSTPU)/tmp.Eff; SystT[N] = (tmp.Eff_SYSTT - tmp.Eff)/tmp.Eff; SystErrT[N] = (tmp.EffE_SYSTT )/tmp.Eff; SystRe[N] = -0.02; SystErrRe[N] = 0.0; SystMB[N]=0.; SystErrMB[N]=0.0; // if((modelSample[s].ModelName().find("Q2")!=string::npos && modelSample[s].ModelName().find("Q2o3")==string::npos) || modelSample[s].ModelName().find("Q3")!=string::npos || modelSample[s].ModelName().find("Q4")!=string::npos || modelSample[s].ModelName().find("Q5")!=string::npos) SystMB[N]=-0.2; if(modelSample[s].ModelName().find("DY_Q")!=string::npos && modelSample[s].ModelName().find("o3")==string::npos) SystMB[N]=-0.2; if(SQRTS==7) { if(modelSample[s].ModelName().find("1o3")!=string::npos) SystTr[N] = -1*sqrt(0.15*0.15 + 0.02*0.02 + 0.05*0.05); else if(modelSample[s].ModelName().find("2o3")!=string::npos) SystTr[N] = -1*sqrt(0.03*0.03 + 0.02*0.02 + 0.05*0.05); else if(IsNeutral) SystTr[N] = -0.05; else SystTr[N] = -1*sqrt(0.05*0.05 + 0.02*0.02 + 0.02*0.02); SystErrTr[N] = 0.0; }else { if(IsNeutral) SystTr[N] = -0.01; else if(modelSample[s].ModelName().find("1o3")!=string::npos) SystTr[N] = -1*sqrt(0.15*0.15 + 0.04*0.04 + 0.05*0.05); else if(modelSample[s].ModelName().find("2o3")!=string::npos) SystTr[N] = -1*sqrt(0.03*0.03 + 0.04*0.04 + 0.05*0.05); else SystTr[N] = -1*sqrt(0.05*0.05 + 0.04*0.04 + 0.01*0.01); SystErrTr[N] = 0.0; } // double Ptemp=max(SystP[N], 0.0), Itemp=max(SystI[N], 0.0), PUtemp=max(SystPU[N], 0.0), Ttemp=max(SystT[N], 0.0); double Ptemp=SystP[N], Itemp=SystI[N], PUtemp=SystPU[N], Ttemp=SystT[N]; SystTotal[N] = -1*sqrt(Ptemp*Ptemp + Itemp*Itemp + PUtemp*PUtemp + Ttemp*Ttemp + SystTr[N]*SystTr[N] + SystRe[N]*SystRe[N] + SystMB[N]*SystMB[N]); // SystErrTotal[N] = sqrt(pow(SystErrP[N],2) + pow(SystErrI[N],2) + pow(SystErrPU[N],2) + pow(SystErrT[N],2) + pow(SystErrTr[N],2) + pow(SystErrRe[N],2) + pow(SystErrMB[N],2) ); SystErrTotal[N] = (tmp.EffE)/tmp.Eff; SystTotal2[N] = -1*SystTotal[N]; if(TypeMode==0 || TypeMode==5)fprintf(pFile, "%30s %7.3f --> %7.3f | %7.3f | %7.3f | %7.3f" ,modelSample[N].Name.c_str(), tmp.Eff, SystP[N], SystI[N], SystPU[N] , SystTotal[N]); else fprintf(pFile, "%30s %7.3f --> %7.3f | %7.3f | %7.3f | %7.3f | %7.3f",modelSample[N].Name.c_str(), tmp.Eff, SystP[N], SystI[N], SystPU[N], SystT[N], SystTotal[N]); //Check if they are variated samples for this point for(unsigned int sv=0;sv<samples.size();sv++){ if(samples[sv].Type!=3)continue; if(samples[sv].Name.find(modelSample[s].Name)!=0)continue; //we expect to have the beginning of the name being identical if(samples[sv].Name.length()!=modelSample[s].Name.length()+3) continue; // variated samples are exactly 3char longer stAllInfo tmpVaried(InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR + "/"+samples[sv].Name+".txt"); if(tmp.Mass<=0) continue; // fprintf(pFile, "%20s%10s %7.3f --> RelDiff=%7.3f\n", " Variation --> ",samples[sv].Name.c_str()+modelSample[s].Name.length()+1, tmpVaried.Eff, (tmpVaried.Eff-tmp.Eff)/tmp.Eff); fprintf(pFile, " | %7.3f (%2s)", (tmpVaried.Eff - tmp.Eff)/tmp.Eff, samples[sv].Name.c_str()+modelSample[s].Name.length()+1); } fprintf(pFile, "\n"); N++; } TGraph* graphSystP = NULL; TGraph* graphSystI = NULL; TGraph* graphSystPU = NULL; TGraph* graphSystT = NULL; TGraph* graphSystTr = NULL; TGraph* graphSystRe = NULL; TGraph* graphSystMB = NULL; TGraph* graphSystTotal = NULL; TGraph* graphSystTotal2 = NULL; if(N>0) { TCanvas* c2 = new TCanvas("c2", "c2",600,600); c2->SetLeftMargin(0.15); // graphSystP = new TGraphErrors(N,Mass,SystP, MassErr, SystErrP); // graphSystI = new TGraphErrors(N,Mass,SystI, MassErr,SystErrI); // graphSystPU = new TGraphErrors(N,Mass,SystPU, MassErr,SystErrPU); // graphSystT = new TGraphErrors(N,Mass,SystT, MassErr,SystErrT); // graphSystTr = new TGraphErrors(N,Mass,SystTr, MassErr,SystErrTr); // graphSystRe = new TGraphErrors(N,Mass,SystRe, MassErr,SystErrRe); // graphSystMB = new TGraphErrors(N,Mass,SystMB, MassErr,SystErrMB); // graphSystTotal = new TGraphErrors(N,Mass,SystTotal, MassErr, SystErrTotal); graphSystP = new TGraph(N,Mass,SystP);//, MassErr, SystErrP); graphSystI = new TGraph(N,Mass,SystI);//, MassErr,SystErrI); graphSystPU = new TGraph(N,Mass,SystPU);//, MassErr,SystErrPU); graphSystT = new TGraph(N,Mass,SystT);//, MassErr,SystErrT); graphSystTr = new TGraph(N,Mass,SystTr);//, MassErr,SystErrTr); graphSystRe = new TGraph(N,Mass,SystRe);//, MassErr,SystErrRe); graphSystMB = new TGraph(N,Mass,SystMB);//, MassErr,SystErrMB); graphSystTotal = new TGraphErrors(N,Mass,SystTotal, MassErr, SystErrTotal); graphSystTotal2 = new TGraph(N,Mass,SystTotal2);//, MassErr, SystErrTotal); TMultiGraph* SystGraphs = new TMultiGraph(); graphSystTotal->GetYaxis()->SetTitle("CrossSection ( pb )"); graphSystTotal->SetLineColor(Color[0]); graphSystTotal->SetMarkerColor(Color[0]); graphSystTotal->SetMarkerStyle(20); graphSystTotal->SetLineWidth(2); graphSystP->SetLineColor(Color[1]); graphSystP->SetMarkerColor(Color[1]); graphSystP->SetMarkerStyle(Marker[1]); graphSystP->SetLineWidth(2); graphSystI->SetLineColor(Color[2]); graphSystI->SetMarkerColor(Color[2]); graphSystI->SetMarkerStyle(Marker[2]); graphSystI->SetLineWidth(2); graphSystPU->SetLineColor(Color[3]); graphSystPU->SetMarkerColor(Color[3]); graphSystPU->SetMarkerStyle(Marker[3]);graphSystPU->SetLineWidth(2); graphSystT->SetLineColor(Color[4]); graphSystT->SetMarkerColor(Color[4]); graphSystT->SetMarkerStyle(Marker[4]); graphSystT->SetLineWidth(2); graphSystTr->SetLineColor(Color[5]); graphSystTr->SetMarkerColor(Color[5]); graphSystTr->SetMarkerStyle(Marker[5]);graphSystTr->SetLineWidth(2); graphSystRe->SetLineColor(Color[6]); graphSystRe->SetMarkerColor(Color[6]); graphSystRe->SetMarkerStyle(Marker[6]);graphSystRe->SetLineWidth(2); graphSystMB->SetLineColor(Color[7]); graphSystMB->SetMarkerColor(Color[7]); graphSystMB->SetMarkerStyle(Marker[7]);graphSystMB->SetLineWidth(2); SystGraphs->Add(graphSystP,"LP0"); SystGraphs->Add(graphSystTr,"LP"); SystGraphs->Add(graphSystRe,"LP"); if(TypeMode!=3)SystGraphs->Add(graphSystI,"LP0"); SystGraphs->Add(graphSystPU,"LP0"); if(TypeMode!=0 && TypeMode!=5)SystGraphs->Add(graphSystT,"LP0"); if(TypeMode==4) SystGraphs->Add(graphSystMB,"LP"); SystGraphs->Add(graphSystTotal,"P"); SystGraphs->Draw("A"); SystGraphs->SetTitle(""); SystGraphs->GetXaxis()->SetTitle("Mass (GeV)"); SystGraphs->GetYaxis()->SetTitle("Relative Change in Efficiency"); SystGraphs->GetYaxis()->SetTitleOffset(1.40); SystGraphs->GetYaxis()->SetRangeUser(-0.55, 0.35); SystGraphs->GetYaxis()->SetNdivisions(520, "X"); TLegend* LEG = new TLegend(0.20,0.9,0.80,0.75); LEG->SetFillColor(0); LEG->SetFillStyle(0); LEG->SetBorderSize(0); LEG->AddEntry(graphSystTr, "Trigger" ,"L"); LEG->AddEntry(graphSystRe, "Reconstruction" ,"L"); if(TypeMode==4)LEG->AddEntry(graphSystMB, "MB" ,"L"); LEG->AddEntry(graphSystP, "P" ,"L"); if(TypeMode!=3)LEG->AddEntry(graphSystI, "dE/dx" ,"L"); LEG->AddEntry(graphSystPU, "Pile Up" ,"L"); if(TypeMode!=0 && TypeMode!=5)LEG->AddEntry(graphSystT, "1/#beta" ,"L"); LEG->AddEntry(graphSystTotal, "Total" ,"P"); LEG->SetNColumns(2); LEG->Draw(); c2->SetLogy(false); c2->SetGridy(false); DrawPreliminary(LegendFromType(InputPattern).c_str(), SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c2,"Results/"+SHAPESTRING+EXCLUSIONDIR+"/", string(prefix+ ModelName + "Uncertainty")); delete c2; //delete SystGraphs; } return graphSystTotal2; } TGraph* MakePlot(FILE* pFile, FILE* talkFile, string InputPattern, string ModelName, int XSectionType, std::vector<stSample>& modelSamples, double& LInt){ std::vector<int> signalPoints; for(unsigned int i=0;i<modelSamples.size();i++) if(XSectionType==0 || stAllInfo(InputPattern+""+SHAPESTRING+EXCLUSIONDIR+"/" + modelSamples[i].Name +".txt").XSec_Exp<1E10) { //Skip 100GeV for DY Q=7 and Q=8 if(XSectionType>0 && (ModelName.find("DY_Q7")!=string::npos || ModelName.find("DY_Q8")!=string::npos) && stAllInfo(InputPattern+""+SHAPESTRING+EXCLUSIONDIR+"/" + modelSamples[i].Name +".txt").Mass==100.0 )continue; signalPoints.push_back(i); } unsigned int N = signalPoints.size(); double* Mass = new double [signalPoints.size()]; double* XSecTh = new double [signalPoints.size()]; double* XSecObs = new double [signalPoints.size()]; double* XSecExp = new double [signalPoints.size()]; stAllInfo* Infos = new stAllInfo[signalPoints.size()]; bool FileFound=false; // int I=0; for(unsigned int i=0;i<signalPoints.size();i++){ Infos [i] = stAllInfo(InputPattern+""+SHAPESTRING+EXCLUSIONDIR+"/" + modelSamples[signalPoints[i]].Name +".txt"); if(Infos[i].Mass!=0) FileFound=true; // if(XSectionType>0 && Infos[i].XSec_Exp>1E10)continue; Mass [i] = Infos[i].Mass; XSecTh [i] = Infos[i].XSec_Th; XSecObs [i] = Infos[i].XSec_Obs; XSecExp [i] = Infos[i].XSec_Exp; LInt = std::max(LInt, Infos[i].LInt); //printf("%i %s\n", (int)FileFound, (InputPattern+""+SHAPESTRING+EXCLUSIONDIR+"/" + modelSamples[signalPoints[i]].Name +".txt").c_str()); // I++; } // N=I; if(XSectionType>0 && FileFound){ //for(unsigned int i=0;i<N;i++)printf("%-18s %5.0f --> Pt>%+6.1f & I>%+5.3f & TOF>%+4.3f & M>%3.0f--> NData=%2.0f NPred=%6.1E+-%6.1E NSign=%6.1E (Eff=%3.2f) Local Significance %3.2f\n",ModelName.c_str(),Infos[i].Mass,Infos[i].WP_Pt,Infos[i].WP_I,Infos[i].WP_TOF,Infos[i].MassCut, Infos[i].NData, Infos[i].NPred, Infos[i].NPredErr, Infos[i].NSign, Infos[i].Eff, Infos[i].Significance); for(unsigned int i=0;i<signalPoints.size();i++){ //for(unsigned int i=0;i<N;i++){ if(Infos[i].WP_TOF==-1){fprintf(pFile,"%-20s & %4.0f & %6.0f & %5.3f & / & %4.0f & %6.3f $\\pm$ %6.3f & %2.0f & %4.3f & %6.1E & %6.1E & %6.1E & %3.2f \\\\\n", ModelName.c_str(), Infos[i].Mass, Infos[i].WP_Pt,Infos[i].WP_I,Infos[i].MassCut, Infos[i].NPred, Infos[i].NPredErr, Infos[i].NData, Infos[i].Eff, Infos[i].XSec_Th,Infos[i].XSec_Exp, Infos[i].XSec_Obs, Infos[i].Significance); }else{ fprintf(pFile,"%-20s & %4.0f & %6.0f & %5.3f & %4.3f & %4.0f & %6.3f $\\pm$ %6.3f & %2.0f & %4.3f & %6.1E & %6.1E & %6.1E & %3.2f \\\\\n", ModelName.c_str(), Infos[i].Mass, Infos[i].WP_Pt,Infos[i].WP_I,Infos[i].WP_TOF,Infos[i].MassCut, Infos[i].NPred, Infos[i].NPredErr, Infos[i].NData, Infos[i].Eff, Infos[i].XSec_Th,Infos[i].XSec_Exp, Infos[i].XSec_Obs, Infos[i].Significance); } bool IsNeutral = (ModelName.find("N",0)<std::string::npos); if(Infos[i].WP_TOF==-1 && (ModelName=="GMSB Stau" || (int)Infos[i].Mass%200==0)) { fprintf(talkFile,"%-20s & %4.0f & %6.0f & %5.3f & / & %4.0f & %6.3f $\\pm$ %6.3f & %2.0f & %4.3f & %3.2f \\\\\n", ModelName.c_str(), Infos[i].Mass, Infos[i].WP_Pt,Infos[i].WP_I,Infos[i].MassCut, Infos[i].NPred, Infos[i].NPredErr, Infos[i].NData, Infos[i].Eff, Infos[i].Significance); fprintf(talkFile, "\\hline\n"); } if(Infos[i].WP_TOF!=-1 && !IsNeutral) { fprintf(talkFile,"%-20s & %4.0f & %6.0f & %5.3f & %4.3f & %4.0f & %6.3f $\\pm$ %6.3f & %2.0f & %4.3f %3.2f \\\\\n", ModelName.c_str(), Infos[i].Mass, Infos[i].WP_Pt,Infos[i].WP_I,Infos[i].WP_TOF,Infos[i].MassCut, Infos[i].NPred, Infos[i].NPredErr, Infos[i].NData, Infos[i].Eff, Infos[i].Significance); fprintf(talkFile, "\\hline\n"); } } } TGraph* graph = NULL; if(XSectionType==0)graph = new TGraph(N,Mass,XSecTh); if(XSectionType==1)graph = new TGraph(N,Mass,XSecExp); if(XSectionType==2)graph = new TGraph(N,Mass,XSecObs); graph->SetTitle(""); graph->GetYaxis()->SetTitle("CrossSection ( pb )"); graph->GetYaxis()->SetTitleOffset(1.40); return graph; } void printSummary(FILE* pFile, FILE* talkFile, string InputPattern, string ModelName, std::vector<stSample>& modelSamples){ TypeMode = TypeFromPattern(InputPattern); for(unsigned int i=0;i<modelSamples.size();i++){ string signal7TeV = modelSamples[i].Name; if(signal7TeV.find("_8TeV")!=string::npos) signal7TeV = signal7TeV.replace(signal7TeV.find("_8TeV"),5, "_7TeV"); string signal8TeV = modelSamples[i].Name; if(signal8TeV.find("_7TeV")!=string::npos) signal8TeV = signal8TeV.replace(signal8TeV.find("_7TeV"),5, "_8TeV"); string signal = signal8TeV; if(signal .find("_8TeV")!=string::npos) signal = signal .replace(signal .find("_8TeV"),5, ""); stAllInfo Infos7(InputPattern+""+SHAPESTRING+"EXCLUSION7TeV"+"/" + signal7TeV +".txt"); stAllInfo Infos8(InputPattern+""+SHAPESTRING+"EXCLUSION8TeV"+"/" + signal8TeV +".txt"); stAllInfo InfosC(InputPattern+""+SHAPESTRING+"EXCLUSIONCOMB"+"/" + signal +".txt"); if(Infos7.Mass<=0 && Infos8.Mass<=0 && InfosC.Mass<=0)continue; if(Infos7.Eff<=0 && Infos8.Eff<=0)continue; double Mass = std::max(Infos7.Mass, Infos8.Mass); TString ModelNameTS = ModelName.c_str(); ModelNameTS.ReplaceAll("_"," "); ModelNameTS.ReplaceAll("8TeV",""); ModelNameTS.ReplaceAll("7TeV",""); if(ModelNameTS.Contains("Stop") && ((int)(Mass)/100)%2!=0)continue; if(ModelNameTS.Contains("Gluino") && ((int)(Mass)/100)%2!=1)continue; if(ModelNameTS.Contains("DY") && ((int)(Mass)/100)%2!=0)continue; if(ModelNameTS.Contains("DC") )continue; char massCut[255]; if(TypeMode<3){sprintf(massCut,"$>%.0f$",Infos8.MassCut);}else{sprintf(massCut," - ");} char Results7[255]; if(Infos7.Mass>0 && TypeMode!=3){sprintf(Results7, "%6.2f & %6.2E & %6.2E & %6.2E", Infos7.Eff, Infos7.XSec_Th,Infos7.XSec_Obs, Infos7.XSec_Exp);}else{sprintf(Results7, " - & - & - & - ");} char Results8[255]; if(Infos8.Mass>0){sprintf(Results8, "%6.2f & %6.2E & %6.2E & %6.2E", Infos8.Eff, Infos8.XSec_Th,Infos8.XSec_Obs, Infos8.XSec_Exp);}else{sprintf(Results8, " - & - & - & - ");} char ResultsC[255]; if(InfosC.Mass>0 && TypeMode!=3){sprintf(ResultsC, "%6.2E & %6.2E", InfosC.XSec_Obs, InfosC.XSec_Exp);} else if(InfosC.Mass>0){sprintf(ResultsC, "%6.2E & %6.2E", Infos8.XSec_Obs/Infos8.XSec_Th, Infos8.XSec_Exp/Infos8.XSec_Th);} else{sprintf(ResultsC, " - & - ");} // char Results7[255]; if(Infos7.Mass>0){sprintf(Results7, "%10s & %10s & %10s & %10s", toLatexRounded(Infos7.Eff).c_str(), toLatexRounded(Infos7.XSec_Th).c_str(),toLatexRounded(Infos7.XSec_Obs).c_str(), toLatexRounded(Infos7.XSec_Exp).c_str());}else{sprintf(Results7, "%10s & %10s & %10s & %10s", "", "", "", "");} // char Results8[255]; if(Infos8.Mass>0){sprintf(Results8, "%10s & %10s & %10s & %10s", toLatexRounded(Infos8.Eff).c_str(), toLatexRounded(Infos8.XSec_Th).c_str(),toLatexRounded(Infos8.XSec_Obs).c_str(), toLatexRounded(Infos8.XSec_Exp).c_str());}else{sprintf(Results8, "%10s & %10s & %10s & %10s", "", "", "", "");} // char ResultsC[255]; if(InfosC.Mass>0){sprintf(ResultsC, "%10s & %10s", toLatexRounded(InfosC.XSec_Obs).c_str(), toLatexRounded(InfosC.XSec_Exp).c_str());}else{sprintf(ResultsC, "%10s & %10s", "", "");} fprintf(pFile,"%-20s & %4.0f & %-7s & %s & %s & %s\\\\\n", ModelNameTS.Data(), Mass, massCut, Results7, Results8, ResultsC); } } void printSummaryPaper(FILE* pFile, FILE* talkFile, string InputPattern, string ModelName, std::vector<stSample>& modelSamples){ for(unsigned int i=0;i<modelSamples.size();i++){ TypeMode = TypeFromPattern(InputPattern); string signal7TeV = modelSamples[i].Name; if(signal7TeV.find("_8TeV")!=string::npos) signal7TeV = signal7TeV.replace(signal7TeV.find("_8TeV"),5, "_7TeV"); string signal8TeV = modelSamples[i].Name; if(signal8TeV.find("_7TeV")!=string::npos) signal8TeV = signal8TeV.replace(signal8TeV.find("_7TeV"),5, "_8TeV"); string signal = signal8TeV; if(signal .find("_8TeV")!=string::npos) signal = signal .replace(signal .find("_8TeV"),5, ""); stAllInfo Infos7(InputPattern+""+SHAPESTRING+"EXCLUSION7TeV"+"/" + signal7TeV +".txt"); stAllInfo Infos8(InputPattern+""+SHAPESTRING+"EXCLUSION8TeV"+"/" + signal8TeV +".txt"); stAllInfo InfosC(InputPattern+""+SHAPESTRING+"EXCLUSIONCOMB"+"/" + signal +".txt"); if(Infos7.Mass<=0 && Infos8.Mass<=0 && InfosC.Mass<=0)continue; if(Infos7.Eff<=0 && Infos8.Eff<=0)continue; double Mass = std::max(Infos7.Mass, Infos8.Mass); TString ModelNameTS = ModelName.c_str(); ModelNameTS.ReplaceAll("8TeV",""); ModelNameTS.ReplaceAll("7TeV",""); if((ModelNameTS.Contains("Gluino_f10") && !ModelNameTS.Contains("f100") && ((int)(Mass)/100)%4==3 && TypeMode==0) || (ModelNameTS.Contains("GluinoN_f10") && !ModelNameTS.Contains("f100") && ((int)(Mass)/100)%4==3 && TypeMode==0) || (ModelNameTS.Contains("Gluino_f50") && ((int)(Mass)/100)%4==3 && TypeMode==3) || (ModelNameTS.Contains("Gluino_f100") && ((int)(Mass)/100)%4==3 && TypeMode==3) || (ModelNameTS.Contains("Stop") && ((int)(Mass)/100)%3==2 && TypeMode==0) || (ModelNameTS.Contains("GMStau") && (Mass==126 || Mass==308 || Mass==494) && TypeMode==2) || (ModelNameTS.Contains("PPStau") && (Mass==126 || Mass==308 || Mass==494) && TypeMode==2) || (ModelNameTS.Contains("DY") && ModelNameTS.Contains("Q1") && !ModelNameTS.Contains("o3") && ((int)(Mass)/100)%3==2 && TypeMode==2) || (ModelNameTS.Contains("DY") && ((int)(Mass)/100)%3==2 && TypeMode==4) || (ModelNameTS.Contains("DY") && ((int)(Mass)/100)%1==0 && TypeMode==5)) { fprintf(pFile,"%s\\\\\n", ModelName.c_str()); char massCut[255]; if(TypeMode<3){sprintf(massCut,"$>%.0f$",Infos8.MassCut);}else{sprintf(massCut," - ");} char Results7[255]; if(Infos7.Mass>0 && TypeMode!=3){sprintf(Results7, "%s & %s & %6.2f", toLatex(Infos7.XSec_Exp).c_str(), toLatex(Infos7.XSec_Obs).c_str(), Infos7.Eff);}else{sprintf(Results7, " - & - & - ");} char Results8[255]; if(Infos8.Mass>0){sprintf(Results8, "%s & %s & %6.2f", toLatex(Infos8.XSec_Exp).c_str(), toLatex(Infos8.XSec_Obs).c_str(), Infos8.Eff);}else{sprintf(Results8, " - & - & - & - ");} char ResultsC[255]; if(InfosC.Mass>0 && TypeMode!=3){sprintf(ResultsC, "%s & %s", toLatex(InfosC.XSec_Exp).c_str(), toLatex(InfosC.XSec_Obs).c_str());} else if(InfosC.Mass>0){sprintf(ResultsC, "%s & %s", toLatex(Infos8.XSec_Exp/Infos8.XSec_Th).c_str(), toLatex(Infos8.XSec_Obs/Infos8.XSec_Th).c_str());} else{sprintf(ResultsC, " - & - ");} // char Results7[255]; if(Infos7.Mass>0){sprintf(Results7, "%10s & %10s & %10s & %10s", toLatexRounded(Infos7.Eff).c_str(), toLatexRounded(Infos7.XSec_Th).c_str(),toLatexRounded(Infos7.XSec_Obs).c_str(), toLatexRounded(Infos7.XSec_Exp).c_str());}else{sprintf(Results7, "%10s & %10s & %10s & %10s", "", "", "", "");} // char Results8[255]; if(Infos8.Mass>0){sprintf(Results8, "%10s & %10s & %10s & %10s", toLatexRounded(Infos8.Eff).c_str(), toLatexRounded(Infos8.XSec_Th).c_str(),toLatexRounded(Infos8.XSec_Obs).c_str(), toLatexRounded(Infos8.XSec_Exp).c_str());}else{sprintf(Results8, "%10s & %10s & %10s & %10s", "", "", "", "");} // char ResultsC[255]; if(InfosC.Mass>0){sprintf(ResultsC, "%10s & %10s", toLatexRounded(InfosC.XSec_Obs).c_str(), toLatexRounded(InfosC.XSec_Exp).c_str());}else{sprintf(ResultsC, "%10s & %10s", "", "");} fprintf(pFile," %4.0f & %-7s & %s & %s & %s\\\\\n", Mass, massCut, Results7, Results8, ResultsC); } } } double GetSignalMeanHSCPPerEvent(string InputPattern, unsigned int CutIndex, double MinRange_, double MaxRange_){ TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); TH2D* Mass = (TH2D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/Mass"); TH2D* MaxEventMass = (TH2D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/MaxEventMass"); TH1D* NTracksPassingSelection = Mass->ProjectionY("NTracksPassingSelection",CutIndex+1,CutIndex+1); TH1D* NEventsPassingSelection = MaxEventMass->ProjectionY("NEventsPassingSelection",CutIndex+1,CutIndex+1); double NTracks = NTracksPassingSelection->Integral(NTracksPassingSelection->GetXaxis()->FindBin(MinRange_), NTracksPassingSelection->GetXaxis()->FindBin(MaxRange_)); double NEvents = NEventsPassingSelection->Integral(NEventsPassingSelection->GetXaxis()->FindBin(MinRange_), NEventsPassingSelection->GetXaxis()->FindBin(MaxRange_)); double toReturn = (float)std::max(1.0,NTracks/ NEvents); delete Mass; delete MaxEventMass; delete NTracksPassingSelection; delete NEventsPassingSelection; delete InputFile; return toReturn; } void DrawModelLimitWithBand(string InputPattern){ int TypeMode = TypeFromPattern(InputPattern); string prefix = "BUG"; switch(TypeMode){ case 0: prefix = "Tk"; break; case 2: prefix = "Mu"; break; case 3: prefix = "Mo"; break; case 4: prefix = "HQ"; break; case 5: prefix = "LQ"; break; } double LInt = 0; for(unsigned int k=0; k<modelVector.size(); k++){ bool isNeutral = false;if(modelVector[k].find("GluinoN")!=string::npos || modelVector[k].find("StopN")!=string::npos)isNeutral = true; if(TypeMode!=0 && isNeutral) continue; unsigned int N = modelMap[modelVector[k]].size(); stAllInfo Infos;double Mass[N], XSecTh[N], XSecExp[N],XSecObs[N], XSecExpUp[N],XSecExpDown[N],XSecExp2Up[N],XSecExp2Down[N]; for(unsigned int i=0;i<N;i++){ Infos = stAllInfo(InputPattern+""+SHAPESTRING+EXCLUSIONDIR+"/" + modelMap[modelVector[k]][i].Name +".txt"); Mass [i]=Infos.Mass; XSecTh [i]=Infos.XSec_Th; XSecObs [i]=Infos.XSec_Obs; XSecExp [i]=Infos.XSec_Exp; XSecExpUp [i]=Infos.XSec_ExpUp; XSecExpDown [i]=Infos.XSec_ExpDown; XSecExp2Up [i]=Infos.XSec_Exp2Up; XSecExp2Down[i]=Infos.XSec_Exp2Down; LInt =std::max(LInt, Infos.LInt); } TGraph* graphtheory = new TGraph(N,Mass,XSecTh); TGraph* graphobs = new TGraph(N,Mass,XSecObs); TGraph* graphexp = new TGraph(N,Mass,XSecExp); TCutG* ExpErr = GetErrorBand("ExpErr" ,N,Mass,XSecExpDown ,XSecExpUp , PlotMinScale, PlotMaxScale); TCutG* Exp2SigmaErr = GetErrorBand("Exp2SigmaErr",N,Mass,XSecExp2Down,XSecExp2Up, PlotMinScale, PlotMaxScale); graphtheory->SetLineStyle(3); graphtheory->SetFillColor(kBlue); graphexp->SetLineStyle(4); graphexp->SetLineColor(kRed); graphexp->SetMarkerStyle(); graphexp->SetMarkerSize(0.); Exp2SigmaErr->SetFillColor(kYellow); Exp2SigmaErr->SetLineColor(kWhite); ExpErr->SetFillColor(kGreen); ExpErr->SetLineColor(kWhite); graphobs->SetLineColor(kBlack); graphobs->SetLineWidth(2); graphobs->SetMarkerColor(kBlack); graphobs->SetMarkerStyle(23); TCanvas* c1 = new TCanvas("c1", "c1",600,600); TMultiGraph* MG = new TMultiGraph(); MG->Add(graphexp ,"LP"); MG->Add(graphobs ,"LP"); MG->Add(graphtheory ,"L"); MG->Draw("A"); Exp2SigmaErr->Draw("f"); ExpErr ->Draw("f"); MG ->Draw("same"); MG->SetTitle(""); MG->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MG->GetYaxis()->SetTitle("#sigma (pb)"); MG->GetYaxis()->SetTitleOffset(1.70); MG->GetYaxis()->SetRangeUser(PlotMinScale,PlotMaxScale); DrawPreliminary(SQRTS, LInt); TLegend* LEG = EXCLUSIONDIR.find("COMB")==string::npos ? new TLegend(0.45,0.58,0.65,0.90) : new TLegend(0.45,0.10,0.65,0.42); //TLegend* LEG = new TLegend(0.40,0.65,0.8,0.90); string headerstr = "95% CL Limits ("; headerstr += LegendFromType(InputPattern) + string(")"); LEG->SetHeader(headerstr.c_str()); LEG->SetFillColor(0); LEG->SetBorderSize(0); LEG->AddEntry(graphtheory, modelMap[modelVector[k]][0].ModelLegend().c_str() ,"L"); LEG->AddEntry(graphexp , "Expected" ,"L"); LEG->AddEntry(ExpErr , "Expected #pm 1#sigma" ,"F"); LEG->AddEntry(Exp2SigmaErr, "Expected #pm 2#sigma ","F"); LEG->AddEntry(graphobs , "Observed" ,"LP"); LEG->Draw(); c1->SetLogy(true); SaveCanvas(c1,"Results/"+SHAPESTRING+EXCLUSIONDIR+"/", string(prefix+ modelVector[k] + "ExclusionLog")); delete c1; } } // This code make the Expected Limit error band divided by expected limit plot for all signal models // I don't like much this function... I started to rewrite it, but more work is still needed to improve it. // I don't think two loops are needed, neither all these arrays... void DrawRatioBands(string InputPattern) { int TypeMode = TypeFromPattern(InputPattern); string prefix = "BUG"; switch(TypeMode){ case 0: prefix = "Tk"; break; case 2: prefix = "Mu"; break; case 3: prefix = "Mo"; break; case 4: prefix = "HQ"; break; case 5: prefix = "LQ"; break; } TCanvas* c2 = new TCanvas("c2", "c2",600,800); TGraph** graphAtheory = new TGraph*[modelVector.size()]; TGraph** graphAobs = new TGraph*[modelVector.size()]; TGraph** graphAexp = new TGraph*[modelVector.size()]; TCutG** ExpAErr = new TCutG* [modelVector.size()]; TCutG** Exp2SigmaAErr = new TCutG* [modelVector.size()]; TPad** padA = new TPad* [modelVector.size()]; double step, top; double LInt = 0; top= 1.0/(modelVector.size()+2); step=(1.0-2.*top)/(modelVector.size()); for(unsigned int k=0; k<modelVector.size(); k++){ bool isNeutral = false;if(modelVector[k].find("GluinoN")!=string::npos || modelVector[k].find("StopN")!=string::npos)isNeutral = true; if(TypeMode!=0 && isNeutral) continue; TPad* pad; if(k<(modelVector.size()-1)){ pad = new TPad(Form("pad%i",k),Form("ExpErr%i",k),0.1,1-top-(k+1)*step,0.9,1-top-step*k);//lower left x, y, topright x, y pad->SetBottomMargin(0.); }else { pad = new TPad(Form("pad%i",k),Form("ExpErr%i",k),0.1,0.0,0.9,1-top-step*(k));//lower left x, y, topright x, y pad->SetBottomMargin(top/(step+top)); } pad->SetLeftMargin(0.1); pad->SetRightMargin(0.); pad->SetTopMargin(0.); padA[k] = pad; padA[k]->Draw(); } for(unsigned int k=0; k<modelVector.size(); k++){ bool isNeutral = false;if(modelVector[k].find("GluinoN")!=string::npos || modelVector[k].find("StopN")!=string::npos)isNeutral = true; if(TypeMode>0 && isNeutral) continue; TMultiGraph* MG = new TMultiGraph(); unsigned int N = modelMap[modelVector[k]].size(); stAllInfo Infos;double Mass[N], XSecTh[N], XSecExp[N],XSecObs[N], XSecExpUp[N],XSecExpDown[N],XSecExp2Up[N],XSecExp2Down[N]; for(unsigned int i=0;i<N;i++){ Infos = stAllInfo(InputPattern+""+SHAPESTRING+EXCLUSIONDIR+"/" + modelMap[modelVector[k]][i].Name +".txt"); Mass [i]=Infos.Mass; XSecTh [i]=Infos.XSec_Th; XSecObs [i]=Infos.XSec_Obs /Infos.XSec_Exp; XSecExp [i]=Infos.XSec_Exp /Infos.XSec_Exp; XSecExpUp [i]=Infos.XSec_ExpUp /Infos.XSec_Exp; XSecExpDown [i]=Infos.XSec_ExpDown /Infos.XSec_Exp; XSecExp2Up [i]=Infos.XSec_Exp2Up /Infos.XSec_Exp; XSecExp2Down[i]=Infos.XSec_Exp2Down/Infos.XSec_Exp; LInt =std::max(LInt, Infos.LInt); } TGraph* graphtheory = new TGraph(N,Mass,XSecTh); TGraph* graphobs = new TGraph(N,Mass,XSecObs); TGraph* graphexp = new TGraph(N,Mass,XSecExp); TCutG* ExpErr = GetErrorBand(Form("ExpErr%i",k) ,N,Mass,XSecExpDown ,XSecExpUp, 0.0, 3.0); TCutG* Exp2SigmaErr = GetErrorBand(Form("Exp2SigmaErr%i",k),N,Mass,XSecExp2Down,XSecExp2Up, 0.0, 3.0); graphAtheory [k] = graphtheory; graphAobs [k] = graphobs; graphAexp [k] = graphexp; ExpAErr [k] = ExpErr; Exp2SigmaAErr[k] = Exp2SigmaErr; graphAtheory [k]->SetLineStyle(3); graphAexp [k]->SetLineStyle(4); graphAexp [k]->SetLineColor(kRed); graphAexp [k]->SetMarkerStyle(); graphAexp [k]->SetMarkerSize(0.); Exp2SigmaAErr[k]->SetFillColor(kYellow); Exp2SigmaAErr[k]->SetLineColor(kWhite); ExpAErr [k]->SetFillColor(kGreen); ExpAErr [k]->SetLineColor(kWhite); graphAobs [k]->SetLineColor(kBlack); graphAobs [k]->SetLineWidth(2); graphAobs [k]->SetMarkerColor(kBlack); graphAobs [k]->SetMarkerStyle(23); padA[k]->cd(); int masst[2] = {0,1250}; int xsect[2] = {2, 1}; TGraph* graph = new TGraph(2,masst,xsect); //fake graph to set xaxis right graph->SetMarkerSize(0.); MG->Add(graph ,"P"); MG->Add(graphAobs[k] ,"LP"); MG->Draw("A"); if(k==0){ TLegend* LEG; LEG = new TLegend(0.13,0.01,0.32,0.99); string headerstr; headerstr = LegendFromType(InputPattern); LEG->SetHeader(headerstr.c_str()); LEG->SetFillStyle(0); LEG->SetBorderSize(0); LEG->AddEntry(ExpAErr[0], "Expected #pm 1#sigma","F"); LEG->SetMargin(0.1); LEG->Draw(); } if(k==1){ TLegend* LEG; LEG = new TLegend(0.13,0.01,0.32,0.99); string headerstr; LEG->SetFillStyle(0); LEG->SetBorderSize(0); LEG->AddEntry(Exp2SigmaAErr[0], "Expected #pm 2#sigma","F"); LEG->AddEntry(graphAobs[0],"Observed" ,"LP"); LEG->SetMargin(0.1); LEG->Draw(); } Exp2SigmaAErr[k]->Draw("f"); ExpAErr[k] ->Draw("f"); MG->Draw("same"); MG->SetTitle(""); if(k==modelVector.size()-1) { MG->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); MG->GetXaxis()->SetTitleSize(0.2); MG->GetXaxis()->SetLabelSize(0.2); } TPaveText *pt; if(TypeMode==0) { if(k!=modelVector.size()-1) pt = new TPaveText(0.45, 0.6, 0.95, 0.87,"LBNDC"); else pt = new TPaveText(0.45, 0.82, 0.95, 0.935,"LBNDC"); } else { if(k!=modelVector.size()-1) pt = new TPaveText(0.55, 0.6, 0.95, 0.87,"LBNDC"); else pt = new TPaveText(0.55, 0.82, 0.95, 0.935,"LBNDC"); } pt->SetBorderSize(0); pt->SetLineWidth(0); pt->SetFillStyle(kWhite); TText *text = pt->AddText(modelMap[modelVector[k]][0].ModelLegend().c_str()); text ->SetTextAlign(12); text ->SetTextSize(0.3); if(k==modelVector.size()-1) text ->SetTextSize(0.5*text ->GetTextSize()); pt->Draw(); MG->GetXaxis()->SetRangeUser(0,1250); MG->GetXaxis()->SetNdivisions(506,"Z"); MG->GetYaxis()->SetRangeUser(0.001,2.99); MG->GetYaxis()->SetNdivisions(303, "Z"); MG->GetYaxis()->SetLabelSize(0.3); if(k==modelVector.size()-1){ MG->GetYaxis()->SetLabelSize(0.15); } } c2->cd(); DrawPreliminary(SQRTS, LInt); TPaveText *pt = new TPaveText(0.1, 0., 0.15, 0.7,"NDC"); string tmp = "95% CL Limits (Relative to Expected Limit)"; TText *text = pt->AddText(tmp.c_str()); text ->SetTextAlign(12); text ->SetTextAngle(90); text ->SetTextSize(0.04); pt->SetBorderSize(0); pt->SetFillColor(0); pt->Draw(); SaveCanvas(c2,"Results/"+SHAPESTRING+EXCLUSIONDIR+"/", string(prefix+"LimitsRatio")); delete c2; } //will run on all possible selection and try to identify which is the best one for this sample void Optimize(string InputPattern, string Data, string signal, bool shape, bool cutFromFile){ printf("Optimize selection for %s in %s\n",signal.c_str(), InputPattern.c_str());fflush(stdout); //get the typeMode from pattern TypeMode = TypeFromPattern(InputPattern); if (TypeMode == 3) RescaleError = 0.0; //Done in step 4 if (TypeMode == 4) RescaleError = 0.20; //Identify the signal sample GetSampleDefinition(samples); CurrentSampleIndex = JobIdToIndex(signal,samples); if(CurrentSampleIndex<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } if(Data.find("7TeV")!=string::npos){SQRTS=7.0;} //IntegratedLuminosity = IntegratedLuminosityFromE(SQRTS); } if(Data.find("8TeV")!=string::npos){SQRTS=8.0;} //IntegratedLuminosity = IntegratedLuminosityFromE(SQRTS); } //For muon only don't run on neutral samples as near zero efficiency can make jobs take very long time if((signal.find("Gluino")!=string::npos || signal.find("Stop")!=string::npos) && signal.find("N")!=string::npos && TypeMode==3) return; //Load all input histograms TFile*InputFile = new TFile((InputPattern + "Histos.root").c_str()); TH1D* HCuts_Pt = (TH1D*)GetObjectFromPath(InputFile, "HCuts_Pt"); TH1D* HCuts_I = (TH1D*)GetObjectFromPath(InputFile, "HCuts_I"); TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, "HCuts_TOF"); TH1D* H_Lumi = (TH1D*)GetObjectFromPath(InputFile, Data+"/IntLumi"); TH1D* H_A = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_A"); TH1D* H_B = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_B"); TH1D* H_C = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_C"); TH1D* H_D = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_D"); TH1D* H_E = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_E"); TH1D* H_F = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_F"); TH1D* H_G = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_G"); TH1D* H_H = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_H"); TH1D* H_P = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_P"); TH1D* H_S = (TH1D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/TOF"); TH2D* MassData = (TH2D*)GetObjectFromPath(InputFile, Data+"/Mass"); TH2D* MassPred = (TH2D*)GetObjectFromPath(InputFile, Data+"/Pred_Mass"); TH2D* MassSign = (TH2D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/Mass" ); if(!MassSign){printf("The sample %s is not present in the root file, returns\n", signal.c_str());return;} TH2D* MassSignP = (TH2D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/Mass_SystP"); TH2D* MassSignI = (TH2D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/Mass_SystI"); TH2D* MassSignM = (TH2D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/Mass_SystM"); TH2D* MassSignT = (TH2D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/Mass_SystT"); TH2D* MassSignPU = (TH2D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/Mass_SystPU" ); TH1D* TotalE = (TH1D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/TotalE" ); TH1D* TotalEPU = (TH1D*)GetObjectFromPath(InputFile, samples[CurrentSampleIndex].Name + "/TotalEPU" ); if(TypeMode==3) { //Need to add in systematic uncertainty of NPred, can't do it later because it comes from two different sources that have different uncertainties TH1D* H_P_Coll = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_P_Coll"); TH1D* H_P_Cosmic = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_P_Cosmic"); for(int i=0; i<H_P->GetNbinsX()+2; i++) { H_P->SetBinError(i, sqrt(H_P->GetBinError(i)*H_P->GetBinError(i) + H_P_Coll->GetBinContent(i)*0.2*H_P_Coll->GetBinContent(i)*0.2 + H_P_Cosmic->GetBinContent(i)*0.8*H_P_Cosmic->GetBinContent(i)*0.8)); } } //If Take the cuts From File --> Load the actual cut index int OptimCutIndex = -1; //int OptimMassWindow; if(cutFromFile){ FILE* pFile = fopen("Analysis_Cuts.txt","r"); if(!pFile){printf("Can't open %s\n","Analysis_Cuts.txt"); return;} while(true){ char line[4096]; string Name_; int TypeMode_; double cutPt_; double cutI_; double cutTOF_; int massWindow_; if(!fgets(line, 4096, pFile))break; char* pch=strtok(line,","); int Arg=0; string tmp; while (pch!=NULL){ switch(Arg){ case 0: tmp = pch; Name_ = tmp.substr(tmp.find("\"")+1,tmp.rfind("\"")-tmp.find("\"")-1); break; case 1: sscanf(pch, "%d", &TypeMode_); break; case 2: sscanf(pch, "%lf", &cutPt_ ); break; case 3: sscanf(pch, "%lf", &cutI_ ); break; case 4: sscanf(pch, "%lf", &cutTOF_); break; case 5: sscanf(pch, "%d", &massWindow_);break; default:break; } pch=strtok(NULL,",");Arg++; } //printf("%s %i %f %f %f %i\n",Name_.c_str(), TypeMode_, cutPt_, cutI_, cutTOF_, massWindow_); if(TypeMode_!=TypeMode)continue; //Not reading the cut line for the right TypeMode string signalNameWithoutEnergy = signal; char str7TeV[]="_7TeV"; char str8TeV[]="_8TeV"; if(signalNameWithoutEnergy.find(str7TeV)!=string::npos)signalNameWithoutEnergy.erase(signalNameWithoutEnergy.find(str7TeV), string(str7TeV).length()); if(signalNameWithoutEnergy.find(str8TeV)!=string::npos)signalNameWithoutEnergy.erase(signalNameWithoutEnergy.find(str8TeV), string(str8TeV).length()); //printf("%s vs %s\n",Name_.c_str(), signalNameWithoutEnergy.c_str()); if(Name_!=signalNameWithoutEnergy )continue; //Not reading the cut line for the right signal sample //We are looking at the right sample --> Now loop over all cuts and identify the cut index of the optimal cut double MinDistance = 10000; for(int CutIndex=0;CutIndex<HCuts_Pt->GetNbinsX();CutIndex++){ double cutDistance = fabs(cutPt_ - HCuts_Pt ->GetBinContent(CutIndex+1)) + fabs(cutI_ - HCuts_I ->GetBinContent(CutIndex+1)) + fabs(cutTOF_ - HCuts_TOF ->GetBinContent(CutIndex+1)); if(cutDistance<MinDistance){MinDistance=cutDistance; OptimCutIndex=CutIndex; }//OptimMassWindow=massWindow_;} } printf("Closest cut index to the cuts provided: %i\n",OptimCutIndex); break; } fclose(pFile); if(OptimCutIndex<0){printf("DID NOT FIND THE CUT TO USE FOR THIS SAMPLE %s\n",signal.c_str());return;} } //normalise the signal samples to XSection * IntLuminosity double LInt = H_Lumi->GetBinContent(1); double norm = samples[CurrentSampleIndex].XSec*LInt/TotalE ->Integral(); //normalize the samples to the actual lumi used for limits double normPU= samples[CurrentSampleIndex].XSec*LInt/(TotalEPU->Integral()>0?TotalEPU->Integral():TotalE->Integral()); MassSign ->Scale(norm); MassSignP ->Scale(norm); MassSignI ->Scale(norm); MassSignM ->Scale(norm); MassSignT ->Scale(norm); MassSignPU->Scale(normPU); //Compute mass range for the cut&count search double Mean=-1,Width=-1; if(!shape && TypeMode<=2){ TH1D* tmpMassSignProj = MassSign->ProjectionY("MassSignProj0",1,1); Mean = tmpMassSignProj->GetMean(); Width = tmpMassSignProj->GetRMS(); MinRange = std::max(0.0, Mean-2*Width); MinRange = tmpMassSignProj->GetXaxis()->GetBinLowEdge(tmpMassSignProj->GetXaxis()->FindBin(MinRange)); //Round to a bin value to avoid counting prpoblem due to the binning. delete tmpMassSignProj; }else{ MinRange = 0; } //prepare output directory and log file string outpath = InputPattern + "/"+SHAPESTRING+EXCLUSIONDIR+"/"; MakeDirectories(outpath); FILE* pFile = fopen((outpath+"/"+signal+".info").c_str(),"w"); if(!pFile)printf("Can't open file : %s\n",(outpath+"/"+signal+".info").c_str()); stAllInfo result; stAllInfo toReturn; //loop on all possible selections and determine which is the best one for(int CutIndex=0;CutIndex<MassData->GetNbinsX();CutIndex++){ //if(CutIndex>25)break; //for debugging purposes //if we don't want to optimize but take instead the cuts from a file, we can skip all other cuts if(OptimCutIndex>=0 && CutIndex!=OptimCutIndex)continue; //make sure the pT cut is high enough to get some statistic for both ABCD and mass shape if(HCuts_Pt ->GetBinContent(CutIndex+1) < 50 && TypeMode!=4){printf("Skip cut=%i because of too lose pT cut\n", CutIndex); continue; } //make sure we have a reliable prediction of the number of events if(OptimCutIndex<0 && H_E->GetBinContent(CutIndex+1) >0 && (H_A->GetBinContent(CutIndex+1)<25 || H_F->GetBinContent(CutIndex+1)<25 || H_G->GetBinContent(CutIndex+1)<25)){printf("Skip cut=%i because of unreliable ABCD prediction\n", CutIndex); continue;} //Skip events where Prediction (AFG/EE) is not reliable if(OptimCutIndex<0 && H_E->GetBinContent(CutIndex+1)==0 && H_A->GetBinContent(CutIndex+1)>0 && (H_C->GetBinContent(CutIndex+1)<25 || H_B->GetBinContent(CutIndex+1)<25)){printf("Skip cut=%i because of unreliable ABCD prediction\n", CutIndex); continue;} //Skip events where Prediction (CB/A) is not reliable if(OptimCutIndex<0 && H_F->GetBinContent(CutIndex+1)>0 && H_A->GetBinContent(CutIndex+1)==0 && (H_B->GetBinContent(CutIndex+1)<25 || H_H->GetBinContent(CutIndex+1)<25)){printf("Skip cut=%i because of unreliable ABCD prediction\n", CutIndex); continue;} //Skip events where Prediction (CB/A) is not reliable if(OptimCutIndex<0 && H_G->GetBinContent(CutIndex+1)>0 && H_F->GetBinContent(CutIndex+1)==0 && (H_C->GetBinContent(CutIndex+1)<25 || H_H->GetBinContent(CutIndex+1)<25)){printf("Skip cut=%i because of unreliable ABCD prediction\n", CutIndex); continue;} //Skip events where Prediction (CH/G) is not reliable //make sure we have a reliable prediction of the shape if(TypeMode<=2){ double N_P = H_P->GetBinContent(CutIndex+1); if(H_E->GetBinContent(CutIndex+1) >0 && (H_A->GetBinContent(CutIndex+1)<0.25*N_P || H_F->GetBinContent(CutIndex+1)<0.25*N_P || H_G->GetBinContent(CutIndex+1)<0.25*N_P)){printf("Skip cut=%i because of unreliable mass prediction\n", CutIndex); continue;} //Skip events where Mass Prediction is not reliable if(H_E->GetBinContent(CutIndex+1)==0 && (H_C->GetBinContent(CutIndex+1)<0.25*N_P || H_B->GetBinContent(CutIndex+1)<0.25*N_P)){printf("Skip cut=%i because of unreliable mass prediction\n", CutIndex); continue;} //Skip events where Mass Prediction is not reliable } //prepare outputs result structure result = toReturn; result.MassMean = Mean; result.MassSigma = Width; result.MassCut = TypeMode<=2?MinRange:0; result.Mass = samples[JobIdToIndex(signal,samples)].Mass; result.XSec_Th = samples[JobIdToIndex(signal,samples)].XSec; result.XSec_Err = samples[JobIdToIndex(signal,samples)].XSec * 0.15; toReturn = result; result.Index = CutIndex; result.WP_Pt = HCuts_Pt ->GetBinContent(CutIndex+1); result.WP_I = HCuts_I ->GetBinContent(CutIndex+1); result.WP_TOF = HCuts_TOF->GetBinContent(CutIndex+1); result.LInt = LInt; //compute the limit for this point and check it run sucessfully (some point may be skipped because of lack of statistics or other reasons) //best expected limit //if(TypeMode<=2){if(!runCombine(true, true, false, InputPattern, signal, CutIndex, shape, true, result, MassData, MassPred, MassSign, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU))continue; //}else {if(!runCombine(true, true, false, InputPattern, signal, CutIndex, shape, true, result, MassData, MassPred, MassSign, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU))continue; //} //no need to precompute the reach when not optimizing the cuts //if(OptimCutIndex<0){ //best significance --> is actually best reach if(TypeMode<=2){if(!runCombine(true, false, true, InputPattern, signal, CutIndex, shape, true, result, MassData, MassPred, MassSign, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU)){printf("runCombine did not converge\n"); continue;} //}else {if(!runCombine(true, false, true, InputPattern, signal, CutIndex, shape, true, result, H_D, H_P, H_S, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU)){printf("runCombine did not converge\n"); continue;} }else {if(!runCombine(true, false, true, InputPattern, signal, CutIndex, shape, true, result, H_D, H_P, MassSign, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU)){printf("runCombine did not converge\n"); continue;} } //}else{ //result.XSec_5Sigma=0.0001;//Dummy number --> will be recomputed later on... but it must be >0 //} //report the result for this point in the log file fprintf(pFile ,"%10s: Testing CutIndex=%4i (Pt>%6.2f I>%6.3f TOF>%6.3f) %3.0f<M<inf Ndata=%+6.2E NPred=%6.3E+-%6.3E SignalEff=%6.3f ExpLimit=%6.3E (%6.3E) Reach=%6.3E",signal.c_str(),CutIndex,HCuts_Pt ->GetBinContent(CutIndex+1), HCuts_I ->GetBinContent(CutIndex+1), HCuts_TOF->GetBinContent(CutIndex+1), MinRange,result.NData,result.NPred, result.NPredErr,result.Eff,result.XSec_Exp, result.XSec_Obs, result.XSec_5Sigma);fflush(stdout); fprintf(stdout ,"%10s: Testing CutIndex=%4i (Pt>%6.2f I>%6.3f TOF>%6.3f) %3.0f<M<inf Ndata=%+6.2E NPred=%6.3E+-%6.3E SignalEff=%6.3f ExpLimit=%6.3E (%6.3E) Reach=%6.3E",signal.c_str(),CutIndex,HCuts_Pt ->GetBinContent(CutIndex+1), HCuts_I ->GetBinContent(CutIndex+1), HCuts_TOF->GetBinContent(CutIndex+1), MinRange,result.NData,result.NPred, result.NPredErr,result.Eff,result.XSec_Exp, result.XSec_Obs, result.XSec_5Sigma);fflush(stdout); // if(result.XSec_Exp<toReturn.XSec_Exp){ if(OptimCutIndex>=0 || (result.XSec_5Sigma>0 && result.XSec_5Sigma<toReturn.XSec_5Sigma)){ toReturn=result; fprintf(pFile ," BestSelection\n");fflush(stdout); fprintf(stdout ," BestSelection\n");fflush(stdout); }else{ fprintf(pFile ,"\n");fflush(stdout); fprintf(stdout ,"\n");fflush(stdout); } }//end of selection cut loop fclose(pFile); //recompute the limit for the final point and save the output in the final directory (also save some plots for the shape based analysis) if(TypeMode<=2){runCombine(false, true, true, InputPattern, signal, toReturn.Index, shape, false, toReturn, MassData, MassPred, MassSign, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU); //}else {runCombine(false, true, true, InputPattern, signal, toReturn.Index, shape, false, toReturn, H_D, H_P, H_S, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU); }else {runCombine(false, true, true, InputPattern, signal, toReturn.Index, shape, false, toReturn, H_D, H_P, MassSign, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU); } //all done, save the result to file toReturn.Save(outpath+"/"+signal+".txt"); } // produce the Higgs combine stat tool datacard void makeDataCard(string outpath, string rootPath, string ChannelName, string SignalName, double Obs, double Pred, double PredRelErr, double Sign, double SignStat, double SignalUnc, bool Shape){ double LumiUnc = (SQRTS==7?1.022:1.044); FILE* pFile = fopen(outpath.c_str(), "w"); fprintf(pFile, "imax 1\n"); fprintf(pFile, "jmax *\n"); fprintf(pFile, "kmax *\n"); if(Shape){ fprintf(pFile, "-------------------------------\n"); fprintf(pFile, "shapes * * %s $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATIC\n",rootPath.c_str()); } fprintf(pFile, "-------------------------------\n"); fprintf(pFile, "bin %s\n",ChannelName.c_str()); fprintf(pFile, "Observation %f\n",Obs); fprintf(pFile, "-------------------------------\n"); fprintf(pFile, "bin %s %s\n",ChannelName.c_str(), ChannelName.c_str()); fprintf(pFile, "process %s pred\n",SignalName.c_str()); fprintf(pFile, "process 0 1\n"); fprintf(pFile, "rate %f %f\n",Sign,std::max(1E-4, Pred) ); //if Pred<1E-4 we have troubles when merging datacards fprintf(pFile, "-------------------------------\n"); fprintf(pFile, "%35s %6s %5.3f 1.0 \n","Lumi" , "lnN", LumiUnc); fprintf(pFile, "%35s %6s - %5.3f\n",(ChannelName+"systP").c_str(), "lnN", PredRelErr); fprintf(pFile, "%35s %6s %5.3f - \n",(ChannelName+"systS").c_str(), "lnN", SignalUnc); if(!Shape){ fprintf(pFile, "%35s %6s %5.3f - \n",(ChannelName+"statS").c_str(), "lnN", std::min(SignStat,2.0)); }else{ fprintf(pFile, "%35s %6s 1.000 - \n",(ChannelName+"statS").c_str(), "shapeN2"); fprintf(pFile, "%35s %6s - 1 \n",(ChannelName+"statP").c_str(), "shapeN2"); fprintf(pFile, "%35s %6s 1.000 - \n","mom" , "shapeN2"); fprintf(pFile, "%35s %6s 1.000 - \n","ias" , "shapeN2"); fprintf(pFile, "%35s %6s 1.000 - \n","ih" , "shapeN2"); fprintf(pFile, "%35s %6s 1.000 - \n","tof" , "shapeN2"); fprintf(pFile, "%35s %6s 1.000 - \n","pu" , "shapeN2"); } fclose(pFile); } //save histogram in root file (and it's statistical variation if it's not observed date histogram) void saveHistoForLimit(TH1* histo, string Name, string Id){ histo ->Write( Name .c_str()); if(Name=="data_obs")return; TH1* statup = (TH1*)histo->Clone((Name+"_stat"+Id+"Up").c_str()); TH1* statdown = (TH1*)histo->Clone((Name+"_stat"+Id+"Down").c_str()); for(int ibin=1; ibin<=statup->GetXaxis()->GetNbins(); ibin++){ statup ->SetBinContent(ibin,statup ->GetBinContent(ibin) + histo->GetBinError(ibin)); statdown->SetBinContent(ibin,statdown->GetBinContent(ibin) - histo->GetBinError(ibin)); } statup ->Write((Name+"_stat"+Id+"Up" ).c_str()); statdown->Write((Name+"_stat"+Id+"Down").c_str()); delete statup; delete statdown; } //save histogram with variation in root file (save it as symetrical up and down variation) void saveVariationHistoForLimit(TH1* histo, TH1* vardown, string Name, string variationName){ TH1* varup = (TH1*)histo ->Clone((Name+"_"+variationName+"Up" ).c_str()); varup->Add(vardown,-1); //varup=x varup->Add(histo,1); varup ->Write((Name+"_"+variationName+"Up" ).c_str()); vardown->Write((Name+"_"+variationName+"Down").c_str()); } //function for debugging only, should be remove soon //this function just compute the result of the shape based analysis using the optimal point from the cut&count analysis void testShapeBasedAnalysis(string InputPattern, string signal){ CurrentSampleIndex = JobIdToIndex(signal, samples); if(CurrentSampleIndex<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } int s = CurrentSampleIndex; string outpath = InputPattern + "/"+SHAPESTRING+EXCLUSIONDIR+"/"; MakeDirectories(outpath); //Get Optimal cut from cut&count optimization stAllInfo result = stAllInfo(InputPattern+"/"+EXCLUSIONDIR+"/"+signal+".txt"); //load all intput histograms TFile* InputFile = new TFile((InputPattern+"/Histos.root").c_str()); TH1D* H_Lumi = (TH1D*)GetObjectFromPath(InputFile, "Data7TeV/IntLumi"); TH2D* MassData = (TH2D*)GetObjectFromPath(InputFile, "Data7TeV/Mass"); TH2D* MassPred = (TH2D*)GetObjectFromPath(InputFile, "Data7TeV/Pred_Mass"); TH2D* MassSign = (TH2D*)GetObjectFromPath(InputFile, samples[s].Name+"/Mass"); TH2D* MassSignP = (TH2D*)GetObjectFromPath(InputFile, samples[s].Name+"/Mass_SystP"); TH2D* MassSignI = (TH2D*)GetObjectFromPath(InputFile, samples[s].Name+"/Mass_SystI"); TH2D* MassSignM = (TH2D*)GetObjectFromPath(InputFile, samples[s].Name+"/Mass_SystM"); TH2D* MassSignT = (TH2D*)GetObjectFromPath(InputFile, samples[s].Name+"/Mass_SystT"); TH2D* MassSignPU = (TH2D*)GetObjectFromPath(InputFile, samples[s].Name+"/Mass_SystPU"); TH1D* TotalE = (TH1D*)GetObjectFromPath(InputFile, samples[s].Name+"/TotalE" ); TH1D* TotalEPU = (TH1D*)GetObjectFromPath(InputFile, samples[s].Name+"/TotalEPU" ); //normalise the signal samples to XSection * IntLuminosity double LInt = H_Lumi->GetBinContent(1); double norm = samples[CurrentSampleIndex].XSec*LInt/TotalE ->Integral(); //normalize the samples to the actual lumi used for limits double normPU= samples[CurrentSampleIndex].XSec*LInt/TotalEPU->Integral(); MassSign ->Scale(norm); MassSignP ->Scale(norm); MassSignI ->Scale(norm); MassSignM ->Scale(norm); MassSignT ->Scale(norm); MassSignPU->Scale(normPU); bool Shape = true; //find range if(Shape){ MinRange = 0; }else{ MinRange = std::max(0.0, result.MassMean-2*result.MassSigma); MinRange = MassSign->GetYaxis()->GetBinLowEdge(MassSign->GetYaxis()->FindBin(MinRange)); //Round to a bin value to avoid counting prpoblem due to the binning. } //compute shape based limits and save it's output runCombine(false, true, true, InputPattern, signal, result.Index, Shape, false, result, MassData, MassPred, MassSign, MassSignP, MassSignI, MassSignM, MassSignT, MassSignPU); //all done, save the results to file result.Save(InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/"+signal+".txt"); } //compute the significance using a ProfileLikelihood (assuming datacard is already produced) double computeSignificance(string datacard, bool expected, string& signal, string massStr, float Strength){ double toReturn = -1; char strengthStr[255]; sprintf(strengthStr,"--expectSignal=%f",Strength); string CodeToExecute = "cd /tmp/;"; if(expected)CodeToExecute += "combine -M ProfileLikelihood -n " + signal + " -m " + massStr + " --significance -t 100 " + strengthStr + " " + datacard + " &> shape_" + signal + ".log;"; else CodeToExecute += "combine -M ProfileLikelihood -n " + signal + " -m " + massStr + " --significance " + datacard + " &> shape_" + signal + ".log;"; CodeToExecute += "cd $OLDPWD;"; system(CodeToExecute.c_str()); char line[4096]; FILE* sFile = fopen((string("/tmp/shape_")+signal + ".log").c_str(), "r"); if(!sFile)std::cout<<"FILE NOT OPEN:"<< (string("/tmp/shape_")+signal + ".log").c_str() << endl; int LineIndex=0; int GarbageI; double GarbageD; while(fgets(line, 4096, sFile)){LineIndex++; if(!expected && LineIndex==3){sscanf(line,"Significance: %lf",&toReturn); break;} // if( expected && LineIndex==7){sscanf(line,"median expected limit: r < %lf @ 95%%CL (%i toyMC)",&toReturn,&GarbageI); break;} if( expected && LineIndex==6){sscanf(line,"mean expected limit: r < %lf +/- %lf @ 95%%CL (%i toyMC)",&toReturn, &GarbageD, &GarbageI); break;} continue; }fclose(sFile); return toReturn; } //run the higgs combine stat tool using predicted mass shape distribution (possibly do shape based analysis and/or cut on mass) OR 1D histogram output from ABCD (only do cut and count without mass cut) bool runCombine(bool fastOptimization, bool getXsection, bool getSignificance, string& InputPattern, string& signal, unsigned int CutIndex, bool Shape, bool Temporary, stAllInfo& result, TH1* MassData, TH1* MassPred, TH1* MassSign, TH1* MassSignP, TH1* MassSignI, TH1* MassSignM, TH1* MassSignT, TH1* MassSignPU){ TH1D *MassDataProj=NULL, *MassPredProj=NULL, *MassSignProj=NULL, *MassSignProjP=NULL, *MassSignProjI=NULL, *MassSignProjM=NULL, *MassSignProjT=NULL, *MassSignProjPU=NULL; double NData, NPredErr, NPred, NSign, NSignP, NSignI, NSignM, NSignT, NSignPU; double NSignErr, NSignPErr, NSignIErr, NSignMErr, NSignTErr, NSignPUErr; double signalsMeanHSCPPerEvent = GetSignalMeanHSCPPerEvent(InputPattern,CutIndex, MinRange, MaxRange); //IF 2D histograms --> we get all the information from there (and we can do shape based analysis AND/OR cut on mass) if(MassData->InheritsFrom("TH2")){ //make the projection of all the 2D input histogram to get the shape for this single point MassDataProj = ((TH2D*)MassData )->ProjectionY("MassDataProj" ,CutIndex+1,CutIndex+1); MassPredProj = ((TH2D*)MassPred )->ProjectionY("MassPredProj" ,CutIndex+1,CutIndex+1); MassSignProj = ((TH2D*)MassSign )->ProjectionY("MassSignProj" ,CutIndex+1,CutIndex+1); MassSignProjP = ((TH2D*)MassSignP )->ProjectionY("MassSignProP" ,CutIndex+1,CutIndex+1); MassSignProjI = ((TH2D*)MassSignI )->ProjectionY("MassSignProI" ,CutIndex+1,CutIndex+1); MassSignProjM = ((TH2D*)MassSignM )->ProjectionY("MassSignProM" ,CutIndex+1,CutIndex+1); MassSignProjT = ((TH2D*)MassSignT )->ProjectionY("MassSignProT" ,CutIndex+1,CutIndex+1); MassSignProjPU = ((TH2D*)MassSignPU)->ProjectionY("MassSignProPU" ,CutIndex+1,CutIndex+1); //count events in the allowed range (infinite for shape based and restricted for cut&count) NData = MassDataProj->Integral(MassDataProj->GetXaxis()->FindBin(MinRange), MassDataProj->GetXaxis()->FindBin(MaxRange)); NPred = MassPredProj->Integral(MassPredProj->GetXaxis()->FindBin(MinRange), MassPredProj->GetXaxis()->FindBin(MaxRange)); NPredErr = pow(NPred*RescaleError,2); for(int i=MassPredProj->GetXaxis()->FindBin(MinRange); i<=MassPredProj->GetXaxis()->FindBin(MaxRange) ;i++){NPredErr+=pow(MassPredProj->GetBinError(i),2);}NPredErr=sqrt(NPredErr); NSign = (MassSignProj ->IntegralAndError(MassSignProj ->GetXaxis()->FindBin(MinRange), MassSignProj ->GetXaxis()->FindBin(MaxRange), NSignErr )) / signalsMeanHSCPPerEvent; NSignErr /= signalsMeanHSCPPerEvent; NSignP = (MassSignProjP ->IntegralAndError(MassSignProjP ->GetXaxis()->FindBin(MinRange), MassSignProjP ->GetXaxis()->FindBin(MaxRange), NSignPErr )) / signalsMeanHSCPPerEvent; NSignPErr /= signalsMeanHSCPPerEvent; NSignI = (MassSignProjI ->IntegralAndError(MassSignProjI ->GetXaxis()->FindBin(MinRange), MassSignProjI ->GetXaxis()->FindBin(MaxRange), NSignIErr )) / signalsMeanHSCPPerEvent; NSignIErr /= signalsMeanHSCPPerEvent; NSignM = (MassSignProjM ->IntegralAndError(MassSignProjM ->GetXaxis()->FindBin(MinRange), MassSignProjM ->GetXaxis()->FindBin(MaxRange), NSignMErr )) / signalsMeanHSCPPerEvent; NSignMErr /= signalsMeanHSCPPerEvent; NSignT = (MassSignProjT ->IntegralAndError(MassSignProjT ->GetXaxis()->FindBin(MinRange), MassSignProjT ->GetXaxis()->FindBin(MaxRange), NSignTErr )) / signalsMeanHSCPPerEvent; NSignTErr /= signalsMeanHSCPPerEvent; NSignPU = (MassSignProjPU->IntegralAndError(MassSignProjPU->GetXaxis()->FindBin(MinRange), MassSignProjPU->GetXaxis()->FindBin(MaxRange), NSignPUErr)) / signalsMeanHSCPPerEvent; NSignPUErr/= signalsMeanHSCPPerEvent; //IF 1D histograms --> we get all the information from the ABCD method output }else{ Shape=false; //can not do shape based if we don't get the shapes NData = MassData ->GetBinContent(CutIndex+1); NPredErr = MassPred ->GetBinError (CutIndex+1); NPred = MassPred ->GetBinContent(CutIndex+1); //NSign = MassSign ->GetBinContent(CutIndex+1) / signalsMeanHSCPPerEvent; //NSignErr = MassSign ->GetBinError (CutIndex+1) / signalsMeanHSCPPerEvent; MassSignProj = ((TH2D*)MassSign )->ProjectionY("MassSignPro" ,CutIndex+1,CutIndex+1); MassSignProjP = ((TH2D*)MassSignP )->ProjectionY("MassSignProP" ,CutIndex+1,CutIndex+1); MassSignProjI = ((TH2D*)MassSignI )->ProjectionY("MassSignProI" ,CutIndex+1,CutIndex+1); MassSignProjM = ((TH2D*)MassSignM )->ProjectionY("MassSignProM" ,CutIndex+1,CutIndex+1); MassSignProjT = ((TH2D*)MassSignT )->ProjectionY("MassSignProT" ,CutIndex+1,CutIndex+1); MassSignProjPU = ((TH2D*)MassSignPU)->ProjectionY("MassSignProPU" ,CutIndex+1,CutIndex+1); NSign = MassSignProj ->IntegralAndError(0, MassSignProj ->GetNbinsX()+1, NSignErr) / signalsMeanHSCPPerEvent; NSignErr /= signalsMeanHSCPPerEvent; NSignP = MassSignProjP ->IntegralAndError(0, MassSignProjP ->GetNbinsX()+1, NSignPErr) / signalsMeanHSCPPerEvent; NSignPErr /= signalsMeanHSCPPerEvent; NSignI = MassSignProjI ->IntegralAndError(0, MassSignProjI ->GetNbinsX()+1, NSignIErr) / signalsMeanHSCPPerEvent; NSignIErr /= signalsMeanHSCPPerEvent; NSignM = MassSignProjM ->IntegralAndError(0, MassSignProjM ->GetNbinsX()+1, NSignMErr) / signalsMeanHSCPPerEvent; NSignMErr /= signalsMeanHSCPPerEvent; NSignT = MassSignProjT ->IntegralAndError(0, MassSignProjT ->GetNbinsX()+1, NSignTErr) / signalsMeanHSCPPerEvent; NSignTErr /= signalsMeanHSCPPerEvent; NSignPU = MassSignProjPU->IntegralAndError(0, MassSignProjPU->GetNbinsX()+1, NSignPUErr) / signalsMeanHSCPPerEvent; NSignPUErr/= signalsMeanHSCPPerEvent; NPredErr = sqrt(pow((NPred* RescaleError), 2) + pow(NPredErr,2)); // incorporate background uncertainty } //skip pathological selection point if(isnan((float)NPred))return false; if(NPred<=0){return false;} //Is <=0 only when prediction failed or is not meaningful (i.e. WP=(0,0,0) ) if(!Shape && NPred>1000){return false;} //When NPred is too big, expected limits just take an infinite time! //compute all efficiencies (not really needed anymore, but it's nice to look at these numbers afterward) double Eff = NSign / (result.XSec_Th*result.LInt); double EffP = NSignP / (result.XSec_Th*result.LInt); double EffI = NSignI / (result.XSec_Th*result.LInt); double EffM = NSignM / (result.XSec_Th*result.LInt); double EffT = NSignT / (result.XSec_Th*result.LInt); double EffPU = NSignPU / (result.XSec_Th*result.LInt); double EffErr = NSignErr / (result.XSec_Th*result.LInt); double EffPErr = NSignPErr / (result.XSec_Th*result.LInt); double EffIErr = NSignIErr / (result.XSec_Th*result.LInt); double EffMErr = NSignMErr / (result.XSec_Th*result.LInt); double EffTErr = NSignTErr / (result.XSec_Th*result.LInt); double EffPUErr = NSignPUErr / (result.XSec_Th*result.LInt); if(Eff==0)return false; // if(Eff<=1E-5)return false; // if Eff<0.001% -> limit will hardly converge and we are probably not interested by this point anyway //no way that this point is optimal bool pointMayBeOptimal = (fastOptimization && !getXsection && getSignificance && ((NPred-3*NPredErr)<=result.NPred || Eff>=result.Eff)); //save these info to the result structure result.Eff = Eff; result.Eff_SYSTP = EffP; result.Eff_SYSTI = EffI; result.Eff_SYSTM = EffM; result.Eff_SYSTT = EffT; result.Eff_SYSTPU= EffPU; result.EffE = EffErr; result.EffE_SYSTP = EffPErr; result.EffE_SYSTI = EffIErr; result.EffE_SYSTM = EffMErr; result.EffE_SYSTT = EffTErr; result.EffE_SYSTPU= EffPUErr; result.NData = NData; result.NPred = NPred; result.NPredErr = NPredErr; result.NSign = NSign; // NSign/=(result.XSec_Th*100.0); //normalize xsection to 10fb NSign/=(100.0); //normalize xsection to 10fb //for shape based analysis we need to save all histograms into a root file char CutIndexStr[255];sprintf(CutIndexStr, "SQRTS%02.0fCut%03.0f",SQRTS, result.Index); if(Shape){ //prepare the histograms and variation //scale to 10fb xsection and to observed events instead of observed tracks MassSignProj ->Scale(1.0/(result.XSec_Th*signalsMeanHSCPPerEvent*100)); MassSignProjP ->Scale(1.0/(result.XSec_Th*signalsMeanHSCPPerEvent*100)); MassSignProjI ->Scale(1.0/(result.XSec_Th*signalsMeanHSCPPerEvent*100)); MassSignProjM ->Scale(1.0/(result.XSec_Th*signalsMeanHSCPPerEvent*100)); MassSignProjT ->Scale(1.0/(result.XSec_Th*signalsMeanHSCPPerEvent*100)); MassSignProjPU->Scale(1.0/(result.XSec_Th*signalsMeanHSCPPerEvent*100)); //Rebin --> keep CPU time reasonable and error small MassDataProj ->Rebin(2); MassPredProj ->Rebin(2); MassSignProj ->Rebin(2); MassSignProjP ->Rebin(2); MassSignProjI ->Rebin(2); MassSignProjM ->Rebin(2); MassSignProjT ->Rebin(2); MassSignProjPU->Rebin(2); //make histo that will contains the shapes for limit string shapeFilePath = "/tmp/shape_"+signal+".root"; TFile* out = new TFile(shapeFilePath.c_str(),"RECREATE"); out->cd(); out->mkdir(CutIndexStr); out->cd(CutIndexStr); //save histo into the file saveHistoForLimit(MassDataProj, "data_obs",""); saveHistoForLimit(MassPredProj, "pred", "P"); saveHistoForLimit(MassSignProj, signal, "S"); saveVariationHistoForLimit(MassSignProj, MassSignProjP , signal, "mom"); saveVariationHistoForLimit(MassSignProj, MassSignProjI , signal, "ias"); saveVariationHistoForLimit(MassSignProj, MassSignProjM , signal, "ih"); saveVariationHistoForLimit(MassSignProj, MassSignProjT , signal, "tof"); saveVariationHistoForLimit(MassSignProj, MassSignProjPU, signal, "pu"); //close the output file out->Close(); } //Need to set what the systematic uncertainty on signal is bool IsNeutral = (signal.find("N")!=std::string::npos); //Reconstruction Efficiency uncertainty double UncEffP=(EffP-Eff)/Eff; double UncEffI=(EffI-Eff)/Eff; double UncEffPU=(EffPU-Eff)/Eff; double UncEffT=(EffT-Eff)/Eff; double UncEffRe = -0.02; double UncEffTr = -10; double UncEffMB = 0.0; //Reset Reco and dEdx uncertainty for fractional as dedicated samples for this if(TypeMode==5){ if(signal.find("1o3")!=string::npos) {UncEffI= -0.25; UncEffRe=0.;} if(signal.find("2o3")!=string::npos) {UncEffI= -0.10; UncEffRe=0.;} } //Reset MB for mCHAMP // if((signal.find("Q2")!=string::npos && signal.find("Q2o3")==string::npos) || signal.find("Q3")!=string::npos || signal.find("Q4")!=string::npos || signal.find("Q5")!=string::npos) UncEffMB=-0.2; if(signal.find("DY_Q")!=string::npos && signal.find("o3")==string::npos) UncEffMB=-0.2; //Trigger efficiency uncertainty if(SQRTS==7) { //Numbers here are 0.05 for muon reconstruction uncertainty, 0.02 for muon trigger synchronization, 0.05 for MET charge suppresses //and 0.02 for MET non charge suppressed if(signal.find("1o3")!=string::npos) UncEffTr = -1*sqrt(0.15*0.15 + 0.02*0.02 + 0.05*0.05); else if(signal.find("2o3")!=string::npos) UncEffTr = -1*sqrt(0.03*0.03 + 0.02*0.02 + 0.05*0.05); else if(IsNeutral) UncEffTr = -0.05; else UncEffTr = -1*sqrt(0.05*0.05 + 0.02*0.02 + 0.02*0.02); } else { //Numbers here are 0.05 for muon reconstruction uncertainty, 0.04 for muon trigger synchronization, 0.01 for MET (charge suppresses or not) if(IsNeutral) UncEffTr = -0.01; else if(signal.find("1o3")!=string::npos) UncEffTr = -1*sqrt(0.15*0.15 + 0.04*0.04 + 0.05*0.05); else if(signal.find("2o3")!=string::npos) UncEffTr = -1*sqrt(0.03*0.03 + 0.04*0.04 + 0.05*0.05); else UncEffTr = -1*sqrt(0.05*0.05 + 0.04*0.04 + 0.01*0.01); } printf("uncertainties = %f %f %f %f %f %f %f\n", UncEffP, UncEffI, UncEffPU, UncEffT, UncEffRe, UncEffTr, UncEffMB); if(isnan((float)UncEffPU))UncEffPU=0.0; double SignalUnc = 1 + sqrt(UncEffP*UncEffP + UncEffI*UncEffI + UncEffPU*UncEffPU + UncEffT*UncEffT + UncEffTr*UncEffTr + UncEffRe*UncEffRe + UncEffMB*UncEffMB); result.TotalUnc = SignalUnc-1; printf("SIGNAL UNCERTAINTY = %f\n",SignalUnc); //build the combine datacard, the same code is used both for cut&count and shape base char TypeStr[255]; sprintf(TypeStr,"Type%i", TypeMode); string JobName = TypeStr+signal; string datacardPath = "/tmp/shape_"+JobName+".dat"; makeDataCard(datacardPath,string("shape_")+JobName+".root", CutIndexStr,signal, NData, NPred, 1.0+(Shape?RescaleError:NPredErr/NPred), NSign, 1.0+fabs(EffErr/Eff), SignalUnc, Shape); char massStr[255]; sprintf(massStr,"%.0f",result.Mass); string test = massStr + signal; if(getSignificance && Temporary){ if(NPred<0.001) NPred=0.001; double SignifValue=0.0; double PrevSignifValue=0; double Strength=0.1*(3*sqrt(NPred)/NSign); if(result.XSec_5Sigma>0 && result.XSec_5Sigma<1E48)Strength=result.XSec_5Sigma/result.XSec_Th; // double SignifValue=0.0;double Strength=0.0005; if(result.XSec_5Sigma>0 && result.XSec_5Sigma<1E50)Strength=result.XSec_5Sigma/result.XSec_Th; double previousXSec_5Sigma=result.XSec_5Sigma; result.XSec_5Sigma = -1; //find signal strength needed to get a 5sigma significance unsigned int l=0; double CountDecrease=0; for(l=0;l<10 && pointMayBeOptimal;l++){ PrevSignifValue = SignifValue; SignifValue = computeSignificance(datacardPath, true, (JobName), massStr, Strength); printf("SIGNAL STRENGTH = %E --> SIGNIFICANCE=%E\n",Strength,SignifValue);fflush(stdout); if(SignifValue<=PrevSignifValue || SignifValue<=0){CountDecrease++;}else{CountDecrease=0;} if(CountDecrease>=2){result.XSec_5Sigma = 1E49; break;} //we found the signal strength that lead to a significance close enough to the 5sigma to stop the loop //OR we know that this point is not going to be a good one --> can do a coarse approximation since the begining // if(fabs(SignifValue-5)<0.75 || (fastOptimization && Strength>=previousXSec_5Sigma && SignifValue<5)){ if(fabs(SignifValue-5)<1.00 || (fastOptimization && previousXSec_5Sigma>=0 && Strength>=previousXSec_5Sigma/(result.XSec_Th/100.0) && SignifValue<5 && SignifValue>=0)){ if(fabs(SignifValue-5)<1.00){printf("Full Converge\n"); }else{printf("Fast Converge\n");} result.XSec_5Sigma = Strength * (5/SignifValue) * (result.XSec_Th/100.0);//xsection in pb printf("XSection for 5sigma discovery = %f = %f * %f * %f\n",result.XSec_5Sigma, Strength, (5/SignifValue), (result.XSec_Th/100.0)); break; } //Not yet at the right significance, change the strength to get close if(isinf((float)SignifValue)){Strength/=5; continue;} //strength is already way too high if(SignifValue<=0 ){Strength*=10; continue;} //strength is already way too low if(SignifValue>5 ){Strength*=std::max( 0.1,(4.95/SignifValue)); continue;} //5/significance could be use but converge faster with 4.9 if(SignifValue<5 ){Strength*=std::min(10.0,(5.05/SignifValue)); continue;} //5/significance could be use, but it converges faster with 5.1 break; } } if(getXsection){ //prepare and run the script that will run the external "combine" tool from the Higgs group //If very low background range too small, set limit at 0.001. Only affects scanning range not final limit if(NPred<0.001) NPred=0.001; char rangeStr[255];sprintf(rangeStr," --rMin %f --rMax %f ", 0.0f, 2*(3*sqrt(NPred)/NSign) ); printf("%f/%f --> %s\n",NSign,NPred,rangeStr); string CodeToExecute = "cd /tmp/;"; CodeToExecute += "combine -M Asymptotic -n " + JobName + " -m " + massStr + rangeStr + " shape_" + JobName+".dat &> shape_" + JobName + ".log;"; CodeToExecute += "cd $OLDPWD;cp /tmp/shape_" + JobName + ".* " + InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/." + ";"; system(CodeToExecute.c_str()); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. TFile* file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+".Asymptotic.mH"+massStr+".root").c_str()); if(!file || file->IsZombie())return false; TTree* tree = (TTree*)file->Get("limit"); if(!tree)return false; double Tmass, Tlimit, TlimitErr; float TquantExp; tree->GetBranch("mh" )->SetAddress(&Tmass ); tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("limitErr" )->SetAddress(&TlimitErr); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.025f){ result.XSec_Exp2Down = Tlimit*(result.XSec_Th/100.0); }else if(TquantExp==0.160f){ result.XSec_ExpDown = Tlimit*(result.XSec_Th/100.0); }else if(TquantExp==0.500f){ result.XSec_Exp = Tlimit*(result.XSec_Th/100.0); }else if(TquantExp==0.840f){ result.XSec_ExpUp = Tlimit*(result.XSec_Th/100.0); }else if(TquantExp==0.975f){ result.XSec_Exp2Up = Tlimit*(result.XSec_Th/100.0); }else if(TquantExp==-1 ){ result.XSec_Obs = Tlimit*(result.XSec_Th/100.0); //will be overwritten afterward }else{printf("Quantil %f unused by the analysis --> check the code\n", TquantExp); } } file->Close(); //RUN FULL HYBRID CLS LIMIT (just for observed limit so far, because it is very slow for expected limits --> should be updated --> FIXME) CodeToExecute = "cd /tmp/;"; CodeToExecute += "combine -M HybridNew -n " + JobName + " -m " + massStr + rangeStr + " shape_" + JobName+".dat > shape_" + JobName + ".log;"; CodeToExecute += "cd $OLDPWD; cp /tmp/shape_" + JobName + ".* " + InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/." + ";"; system(CodeToExecute.c_str()); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+".HybridNew.mH"+massStr+".root").c_str()); if(!file || file->IsZombie())return false; tree = (TTree*)file->Get("limit"); if(!tree)return false; tree->GetBranch("mh" )->SetAddress(&Tmass ); tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("limitErr" )->SetAddress(&TlimitErr); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==-1 ){ result.XSec_Obs = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f unused by the analysis --> check the code\n", TquantExp); } } file->Close(); if(FullExpLimit) { //Create grid to find expected limits in Cls //Number of different signal strengths to try int gridPoints=30; //Normalize to 10/fb double Down2 = 0.5*result.XSec_Exp2Down*100/result.XSec_Th; double Up2 = 1.3*result.XSec_Exp2Up*100/result.XSec_Th; double Step=(Up2-Down2)/gridPoints; CodeToExecute = "cd /tmp/;"; for (int i=0; i<gridPoints+1; i++) { char Seed[1024]; sprintf(Seed,"%i",i); char PointStr[1024]; double Point = Down2 + i*Step; sprintf(PointStr,"%6.8f",Point); //Don't include mass string here or else it won't work CodeToExecute += "combine shape_" + JobName + ".dat -M HybridNew --freq --fork 1 -T 500 --clsAcc 0 -n " + JobName + " --saveHybridResult --saveToys -s " + Seed + " -i 8 --rMax 1E20 --singlePoint " + PointStr + " >> shape_" + JobName + "Exp.log;"; } CodeToExecute += "hadd -f higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root higgsCombine"+JobName+".HybridNew.mH120.*.root >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.5 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.16 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.84 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.025 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.975 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "cd $OLDPWD; cp /tmp/shape_" + JobName + "Exp.* " + InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/." + ";"; system(CodeToExecute.c_str()); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.500.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.500.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.500f){ result.XSec_Exp = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.5 --> check the code\n", TquantExp); } } file->Close(); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.160.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.160.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.160f){ result.XSec_ExpDown = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.16 --> check the code\n", TquantExp); } } file->Close(); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.840.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.840.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.840f){ result.XSec_ExpUp = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.84 --> check the code\n", TquantExp); } } file->Close(); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.025.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.025.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.025f){ result.XSec_Exp2Down = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.025 --> check the code\n", TquantExp); } } file->Close(); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.975.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.975.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.975f){ result.XSec_Exp2Up = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.975 --> check the code\n", TquantExp); } } file->Close(); } } if(!Temporary && getSignificance){ result.Significance = computeSignificance(datacardPath, false, (JobName), massStr, 1.0); } //makePlots (only for shape based analysis) if(Shape && !Temporary){ TCanvas* c1 = new TCanvas("c1", "c1",1200,600); c1->Divide(2,1); (c1->cd(1))->SetLogy(true); double Max = 2.0 * std::max(std::max(MassDataProj->GetMaximum(), MassPredProj->GetMaximum()), MassSignProj->GetMaximum()); double Min = 0.01; MassSignProj->SetStats(kFALSE); MassSignProj->SetMaximum(Max); MassSignProj->SetMinimum(Min); MassSignProj->GetXaxis()->SetRangeUser(0,1400); MassSignProj->GetXaxis()->SetNdivisions(505,"X"); MassSignProj->SetMarkerStyle(21); MassSignProj->SetMarkerColor(kBlue-9); MassSignProj->SetMarkerSize(1.5); MassSignProj->SetLineColor(kBlue-9); MassSignProj->SetFillColor(kBlue-9); MassSignProj->Draw("HIST"); TH1D* MassPredSignProj = (TH1D*)MassPredProj->Clone("predWithSign"); MassPredSignProj->Add(MassSignProj); MassPredSignProj->SetMarkerStyle(0); MassPredSignProj->SetLineColor(kBlue-10); MassPredSignProj->SetLineStyle(1); MassPredSignProj->SetLineWidth(4); MassPredSignProj->SetFillColor(0); MassPredSignProj->SetFillStyle(0); MassPredSignProj->Draw("same HIST C"); MassPredProj->SetBinContent(MassPredProj->GetNbinsX(), MassPredProj->GetBinContent(MassPredProj->GetNbinsX()) + MassPredProj->GetBinContent(MassPredProj->GetNbinsX()+1)); MassPredProj->SetMarkerStyle(22); MassPredProj->SetMarkerColor(2); MassPredProj->SetMarkerSize(1.5); MassPredProj->SetLineColor(2); MassPredProj->SetFillColor(8); MassPredProj->SetFillStyle(3001); MassPredProj->Draw("same E3"); MassDataProj->SetBinContent(MassDataProj->GetNbinsX(), MassDataProj->GetBinContent(MassDataProj->GetNbinsX()) + MassDataProj->GetBinContent(MassDataProj->GetNbinsX()+1)); MassDataProj->SetMarkerStyle(20); MassDataProj->SetMarkerColor(1); MassDataProj->SetMarkerSize(1.0); MassDataProj->SetLineColor(1); MassDataProj->SetFillColor(0); MassDataProj->Draw("same E1"); TLegend* leg = new TLegend(0.40,0.75,0.80,0.93); leg->SetHeader(NULL); leg->SetFillColor(0); leg->SetFillStyle(0); leg->SetBorderSize(0); leg->AddEntry(MassDataProj,"Data", "P"); leg->AddEntry(MassPredProj,"Prediction", "FP"); leg->AddEntry(MassSignProj,signal.c_str(), "F"); leg->AddEntry(MassPredSignProj,(string("Prediction + ")+signal).c_str(), "L"); leg->Draw(); (c1->cd(2))->SetLogy(true); // (c1->cd(2))->SetGridy(true); TH1* MassSignProjRatio = (TH1D*)MassSignProj ->Clone(signal.c_str() ); MassSignProjRatio ->SetLineColor(1); MassSignProjRatio ->SetMarkerColor(1); MassSignProjRatio ->SetMarkerStyle(0); TH1* MassSignProjPRatio = (TH1D*)MassSignProjP->Clone("mom"); MassSignProjPRatio->SetLineColor(2); MassSignProjPRatio->SetMarkerColor(2); MassSignProjPRatio->SetMarkerStyle(20); TH1* MassSignProjIRatio = (TH1D*)MassSignProjI->Clone("Ias"); MassSignProjIRatio->SetLineColor(4); MassSignProjIRatio->SetMarkerColor(4); MassSignProjIRatio->SetMarkerStyle(21); TH1* MassSignProjMRatio = (TH1D*)MassSignProjM->Clone("Ih"); MassSignProjMRatio->SetLineColor(3); MassSignProjMRatio->SetMarkerColor(3); MassSignProjMRatio->SetMarkerStyle(22); TH1* MassSignProjTRatio = (TH1D*)MassSignProjT->Clone("TOF"); MassSignProjTRatio->SetLineColor(8); MassSignProjTRatio->SetMarkerColor(8); MassSignProjTRatio->SetMarkerStyle(23); TH1* MassSignProjLRatio = (TH1D*)MassSignProjPU->Clone("pu"); MassSignProjLRatio->SetLineColor(6); MassSignProjLRatio->SetMarkerColor(6); MassSignProjLRatio->SetMarkerStyle(33); //MassSignProjPRatio->Divide(MassSignProjPRatio, MassSignProjRatio,1,1, "B"); //MassSignProjIRatio->Divide(MassSignProjIRatio, MassSignProjRatio,1,1, "B"); //MassSignProjMRatio->Divide(MassSignProjMRatio, MassSignProjRatio,1,1, "B"); //MassSignProjTRatio->Divide(MassSignProjTRatio, MassSignProjRatio,1,1, "B"); //MassSignProjLRatio->Divide(MassSignProjLRatio, MassSignProjRatio,1,1, "B"); //MassSignProjRatio ->Divide(MassSignProjRatio , MassSignProjRatio,1,1, "B"); MassSignProjRatio->SetStats(kFALSE); MassSignProjRatio->SetFillColor(0); MassSignProjRatio->SetLineWidth(2); MassSignProjRatio->GetXaxis()->SetRangeUser(0,1400); MassSignProjRatio->GetXaxis()->SetNdivisions(505,"X"); MassSignProjRatio->GetYaxis()->SetNdivisions(505,"X"); MassSignProjRatio->SetMaximum(MassSignProjPRatio->GetMaximum()*2.0); MassSignProjRatio->SetMinimum(Min); //MassSignProjRatio->SetMaximum(2);//Max); //MassSignProjRatio->SetMinimum(0);//Min); //MassSignProjRatio->Reset(); //use this histogram as a framework only MassSignProjRatio->Draw("HIST"); //TLine l1(0,1,1400,1);l1.SetLineColor(1); l1.SetLineWidth(2); l1.Draw("same"); MassSignProjPRatio->Draw("same E1"); MassSignProjIRatio->Draw("same E1"); MassSignProjMRatio->Draw("same E1"); MassSignProjTRatio->Draw("same E1"); MassSignProjLRatio->Draw("same E1"); TLegend* leg2 = new TLegend(0.45,0.65,0.85,0.93); leg2->SetHeader("Variations"); leg2->SetFillColor(0); leg2->SetFillStyle(0); leg2->SetBorderSize(0); leg2->AddEntry(MassSignProjRatio,signal.c_str(), "L"); leg2->AddEntry(MassSignProjPRatio,"momentum", "P"); leg2->AddEntry(MassSignProjIRatio,"Ias", "P"); leg2->AddEntry(MassSignProjMRatio,"Ih", "P"); leg2->AddEntry(MassSignProjTRatio,"tof", "P"); leg2->AddEntry(MassSignProjLRatio,"pu", "P"); leg2->Draw(); SaveCanvas(c1, InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/shape", signal, true); delete leg2; delete leg; delete c1; delete MassPredSignProj; delete MassSignProjRatio; delete MassSignProjPRatio; delete MassSignProjIRatio; delete MassSignProjMRatio; delete MassSignProjTRatio; delete MassSignProjLRatio; } //all done, clean everything and return true delete MassDataProj; delete MassPredProj; delete MassSignProj; delete MassSignProjP; delete MassSignProjI; delete MassSignProjM; delete MassSignProjT; delete MassSignProjPU; return true; } bool Combine(string InputPattern, string signal7, string signal8){ // CurrentSampleIndex = JobIdToIndex(signal, samples); if(CurrentSampleIndex<0){ printf("There is no signal corresponding to the JobId Given\n"); return false; } // int s = CurrentSampleIndex; int TypeMode = TypeFromPattern(InputPattern); string outpath = InputPattern + "/"+SHAPESTRING+EXCLUSIONDIR+"/"; MakeDirectories(outpath); //Get Optimal cut from sample11 stAllInfo result11 = stAllInfo(InputPattern+"/EXCLUSION7TeV/"+signal7+".txt"); //Get Optimal cut from sample12 stAllInfo result12 = stAllInfo(InputPattern+"/EXCLUSION8TeV/"+signal8+".txt"); stAllInfo result = result12; char massStr[255]; sprintf(massStr,"%.0f",result.Mass); string signal = signal7; if(signal.find("_7TeV")!=string::npos){signal.replace(signal.find("_7TeV"),5, "");} char TypeStr[100] ;sprintf(TypeStr,"Type%i", TypeMode); string JobName = TypeStr+signal; FILE* pFileTmp = NULL; bool is7TeVPresent = true; pFileTmp = fopen((InputPattern+"/EXCLUSION7TeV/shape_"+(TypeStr+signal7)+".dat").c_str(), "r"); if(!pFileTmp){is7TeVPresent=false;}else{fclose(pFileTmp);} if(TypeMode==3) is7TeVPresent=false; bool is8TeVPresent = true; pFileTmp = fopen((InputPattern+"/EXCLUSION8TeV/shape_"+(TypeStr+signal8)+".dat").c_str(), "r"); if(!pFileTmp){is8TeVPresent=false;}else{fclose(pFileTmp);} string CodeToExecute = "combineCards.py "; if(is7TeVPresent)CodeToExecute+=" " + InputPattern+"/EXCLUSION7TeV/shape_"+(TypeStr+signal7)+".dat "; if(is8TeVPresent)CodeToExecute+=" " + InputPattern+"/EXCLUSION8TeV/shape_"+(TypeStr+signal8)+".dat "; CodeToExecute+=" > " + outpath+"shape_"+JobName+".dat "; system(CodeToExecute.c_str()); printf("%s \n",CodeToExecute.c_str()); result.XSec_Th = 1.0; //Muon only uses just 2012 if(TypeMode==3) { result.XSec_Obs=result12.XSec_Obs/result12.XSec_Th; result.XSec_Exp=result12.XSec_Exp/result12.XSec_Th; result.XSec_ExpUp=result12.XSec_ExpUp/result12.XSec_Th; result.XSec_ExpDown=result12.XSec_ExpDown/result12.XSec_Th; result.XSec_Exp2Up=result12.XSec_Exp2Up/result12.XSec_Th; result.XSec_Exp2Down=result12.XSec_Exp2Down/result12.XSec_Th; result.Save(InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/"+signal+".txt"); return true; } result.LInt = result11.LInt + result12.LInt ; result.NSign = result11.NSign + result12.NSign; result.NData = result11.NData + result12.NData; result.NPred = result11.NPred + result12.NPred; result.NPredErr = sqrt(pow(result11.NPredErr,2) + pow(result12.NPredErr,2)); //compute combined significance CodeToExecute = "cp " + outpath+"shape_"+JobName+".dat /tmp/.;"; system(CodeToExecute.c_str()); result.Significance = computeSignificance(string("shape_")+JobName+".dat", false, signal, massStr, 1.0); printf("Combined Significance = %f (%s)\n", result.Significance, (outpath+"shape_"+JobName+".dat").c_str()); double NPred = result.NPred; double NSign = result.NSign / 100.0; //ALL CODE BELOW IS A BIT DIFFERENT THAN THE ONE USED IN runCombined, BECAUSE HERE WE KEEP THE RESULTS ON LIMIT IN TERMS OF SIGNAL STRENGTH (r=SigmaObs/SigmaTH) if(true){ //prepare and run the script that will run the external "combine" tool from the Higgs group //If very low background range too small, set limit at 0.001. Only affects scanning range not final limit if(NPred<0.001) NPred=0.001; char rangeStr[255];sprintf(rangeStr," --rMin %f --rMax %f ", 0.0f, 2*(3*sqrt(NPred)/NSign) ); printf("%f/%f --> %s\n",NSign,NPred,rangeStr); string CodeToExecute = "cp " + outpath+"shape_"+JobName+".dat /tmp/.;"; CodeToExecute += "cd /tmp/;"; CodeToExecute += "combine -M Asymptotic -n " + JobName + " -m " + massStr + rangeStr + " shape_" + JobName+".dat &> shape_" + JobName + ".log;"; CodeToExecute += "cd $OLDPWD;cp /tmp/shape_" + JobName + ".* " + InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/." + ";"; system(CodeToExecute.c_str()); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. TFile* file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+".Asymptotic.mH"+massStr+".root").c_str()); if(!file || file->IsZombie())return false; TTree* tree = (TTree*)file->Get("limit"); if(!tree)return false; double Tmass, Tlimit, TlimitErr; float TquantExp; tree->GetBranch("mh" )->SetAddress(&Tmass ); tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("limitErr" )->SetAddress(&TlimitErr); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.025f){ result.XSec_Exp2Down = Tlimit/100.0; }else if(TquantExp==0.160f){ result.XSec_ExpDown = Tlimit/100.0; }else if(TquantExp==0.500f){ result.XSec_Exp = Tlimit/100.0; }else if(TquantExp==0.840f){ result.XSec_ExpUp = Tlimit/100.0; }else if(TquantExp==0.975f){ result.XSec_Exp2Up = Tlimit/100.0; }else if(TquantExp==-1 ){ result.XSec_Obs = Tlimit/100.0; }else{printf("Quantil %f unused by the analysis --> check the code\n", TquantExp); } } file->Close(); //RUN FULL HYBRID CLS LIMIT (just for observed limit so far, because it is very slow for expected limits --> should be updated --> FIXME) CodeToExecute = "cd /tmp/;"; CodeToExecute += "combine -M HybridNew -n " + JobName + " -m " + massStr + rangeStr + " shape_" + JobName+".dat &> shape_" + JobName + ".log;"; CodeToExecute += "cd $OLDPWD;cp /tmp/shape_" + JobName + ".* " + InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/." + ";"; system(CodeToExecute.c_str()); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+".HybridNew.mH"+massStr+".root").c_str()); if(!file || file->IsZombie())return false; tree = (TTree*)file->Get("limit"); if(!tree)return false; tree->GetBranch("mh" )->SetAddress(&Tmass ); tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("limitErr" )->SetAddress(&TlimitErr); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); // if(TquantExp==0.025f){ result.XSec_Exp2Down = Tlimit*(result.XSec_Th/100.0); // }else if(TquantExp==0.160f){ result.XSec_ExpDown = Tlimit*(result.XSec_Th/100.0); // }else if(TquantExp==0.500f){ result.XSec_Exp = Tlimit*(result.XSec_Th/100.0); // }else if(TquantExp==0.840f){ result.XSec_ExpUp = Tlimit*(result.XSec_Th/100.0); // }else if(TquantExp==0.975f){ result.XSec_Exp2Up = Tlimit*(result.XSec_Th/100.0); // }else if(TquantExp==-1 ){ result.XSec_Obs = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f unused by the analysis --> check the code\n", TquantExp); } } file->Close(); if(FullExpLimit) { //Create grid to find expected limits in Cls //Number of different signal strengths to try int gridPoints=30; //Normalize to 10/fb double Down2 = 0.5*result.XSec_Exp2Down*100/result.XSec_Th; double Up2 = 1.3*result.XSec_Exp2Up*100/result.XSec_Th; double Step=(Up2-Down2)/gridPoints; CodeToExecute = "cd /tmp/;"; CodeToExecute += "echo 'Finding Expected Limits' > shape_" + JobName + "Exp.log;"; for (int i=0; i<gridPoints+1; i++) { char Seed[1024]; sprintf(Seed,"%i",i); char PointStr[1024]; double Point = Down2 + i*Step; sprintf(PointStr,"%6.8f",Point); //Don't include mass string here or else it won't work CodeToExecute += "combine shape_" + JobName + ".dat -M HybridNew --freq --fork 1 -T 500 --clsAcc 0 -n " + JobName + " --saveHybridResult --saveToys -s " + Seed + " -i 8 --rMax 1E20 --singlePoint " + PointStr + " >> shape_" + JobName + "Exp.log;"; } CodeToExecute += "hadd -f higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root higgsCombine"+JobName+".HybridNew.mH120.*.root >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.5 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.16 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.84 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.025 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "combine shape_" + JobName+".dat -M HybridNew --grid=higgsCombine"+JobName+".HybridNew.mH"+massStr+"grid.root -n " + JobName + "Expected --expectedFromGrid 0.975 >> shape_" + JobName + "Exp.log;"; CodeToExecute += "cd $OLDPWD; cp /tmp/shape_" + JobName + "Exp.* " + InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/." + ";"; system(CodeToExecute.c_str()); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.500.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.500.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.500f){ result.XSec_Exp = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.5 --> check the code\n", TquantExp); } } file->Close(); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.160.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.160.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.160f){ result.XSec_ExpDown = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.16 --> check the code\n", TquantExp); } } file->Close(); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.840.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.840.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.840f){ result.XSec_ExpUp = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.84 --> check the code\n", TquantExp); } } file->Close(); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.025.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.025.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.025f){ result.XSec_ExpDown = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.025 --> check the code\n", TquantExp); } } file->Close(); //if all went well, the combine tool created a new file containing the result of the limit in the form of a TTree //we can open this TTree and access the values for the expected limit, uncertainty bands, and observed limits. file = TFile::Open((string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.975.root").c_str()); if(!file || file->IsZombie()){printf("Can't fine file %s", (string("/tmp/")+"higgsCombine"+JobName+"Expected.HybridNew.mH120.quant0.975.root").c_str()); return false;} tree = (TTree*)file->Get("limit"); if(!tree){printf("Can't find tree named limit"); return false;} tree->GetBranch("limit" )->SetAddress(&Tlimit ); tree->GetBranch("quantileExpected")->SetAddress(&TquantExp); for(int ientry=0;ientry<tree->GetEntriesFast();ientry++){ tree->GetEntry(ientry); if(TquantExp==0.975f){ result.XSec_ExpDown = Tlimit*(result.XSec_Th/100.0); }else{printf("Quantil %f should be 0.975 --> check the code\n", TquantExp); } } file->Close(); } } //all done, save the results to file result.Save(InputPattern+"/"+SHAPESTRING+EXCLUSIONDIR+"/"+signal+".txt"); return true; } bool useSample(int TypeMode, string sample) { if(TypeMode==0 && (sample=="Gluino_f10" || sample=="GluinoN_f10" || sample=="StopN" || sample=="Stop" || sample=="DY_Q2o3" || sample=="GMStau" || sample=="PPStau")) return true; if(TypeMode==2 && (sample=="Gluino_f10" || sample=="Gluino_f50" || sample=="Stop" || sample=="GMStau" || sample=="PPStau" || sample=="DY_Q2o3" || sample=="DY_Q1")) return true; if(TypeMode==3 && (sample=="Stop" || sample=="Gluino_f10" || sample=="Gluino_f50" || sample=="Gluino_f100")) return true; if(TypeMode==4) return true; if(TypeMode==5) return true; return false; } string toLatex(double value) { char toReturn[255]; if(value<0.00095) { int exponent=0; while (value<1) { exponent++; value*=10; } sprintf(toReturn,"$ %6.1f \\times 10^{-%i}$",value, exponent); } else sprintf(toReturn,"%.2g",value); string stringReturn = toReturn; return stringReturn; }