// Original Author: Loic Quertenmont #include "Analysis_Global.h" #include "Analysis_Global.h" #include "Analysis_CommonFunction.h" #include "Analysis_PlotFunction.h" #include "Analysis_PlotStructure.h" #include "Analysis_Samples.h" #include "tdrstyle.C" #include "Math/QuantFuncMathCore.h" #include "TMath.h" #include "TGraphAsymmErrors.h" #include "TPaletteAxis.h" #include "TColor.h" using namespace std; /////////////////////////// FUNCTION DECLARATION ///////////////////////////// void MassPrediction(string InputPattern, unsigned int CutIndex, string HistoSuffix="Mass", bool showMC=true, string Data="Data8TeV"); void PredictionAndControlPlot(string InputPattern, string Data, unsigned int CutIndex, unsigned int CutIndex_Flip); void CutFlow(string InputPattern, unsigned int CutIndex=0); void SelectionPlot (string InputPattern, unsigned int CutIndex, unsigned int CutIndexTight); void Make2DPlot_Core(string ResultPattern, unsigned int CutIndex); void SignalMassPlot(string InputPattern, unsigned int CutIndex); void GetSystematicOnPrediction(string InputPattern, string DataName="Data8TeV"); void MakeExpLimitpLot(string Input, string Output); void CosmicBackgroundSystematic(string InputPattern, string DataType="8TeV"); void CheckPrediction(string InputPattern, string HistoSuffix="_Flip", string DataType="Data8TeV"); void CheckPredictionBin(string InputPattern, string HistoSuffix="_Flip", string DataType="Data8TeV", string bin=""); void CollisionBackgroundSystematicFromFlip(string InputPattern, string DataType="Data8TeV"); void Make2DPlot_Special(string ResultPattern,string ResultPattern2); //, unsigned int CutIndex); void CompareRecoAndGenPt(string ResultPattern); void CheckPUDistribution(string InputPattern, unsigned int CutIndex); std::vector<stSample> samples; /////////////////////////// CODE PARAMETERS ///////////////////////////// void Analysis_Step5() { 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); //gStyle->SetTextFont(43); GetSampleDefinition(samples); string InputPattern; unsigned int CutIndex; unsigned int CutIndex_Flip=1; unsigned int CutIndexTight; std::vector<string> Legends; std::vector<string> Inputs; // CheckPUDistribution("Results/Type0/", 0); // CheckPUDistribution("Results/Type4/", 0); // return; // Collect the functions that make plots that are in paper Make2DPlot_Special("Results/Type0/", "Results/Type5/"); InputPattern = "Results/Type0/"; CutIndex = 4; CutIndexTight = 84; MassPrediction(InputPattern, CutIndex, "Mass", false, "8TeV_LooseNoSMMC"); MassPrediction(InputPattern, CutIndex, "Mass", false, "7TeV_LooseNoSMMC"); MassPrediction(InputPattern, CutIndexTight, "Mass", false, "8TeV_TightNoSMMC"); MassPrediction(InputPattern, CutIndexTight, "Mass", false, "7TeV_TightNoSMMC"); InputPattern = "Results/Type2/"; CutIndex = 16; CutIndexTight = 905; CutIndex_Flip=16; Make2DPlot_Core(InputPattern, 0); MassPrediction(InputPattern, CutIndex, "Mass", false, "8TeV_LooseNoSMMC"); MassPrediction(InputPattern, CutIndex, "Mass", false, "7TeV_LooseNoSMMC"); MassPrediction(InputPattern, CutIndexTight, "Mass", false, "8TeV_TightNoSMMC"); MassPrediction(InputPattern, CutIndexTight, "Mass", false, "7TeV_TightNoSMMC"); PredictionAndControlPlot(InputPattern, "Data7TeV", CutIndex, CutIndex_Flip); PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); InputPattern = "Results/Type3/"; CutIndex = 96; CutIndexTight = 96; CutIndex_Flip=54; PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); InputPattern = "Results/Type4/"; CutIndex = 21; CutIndexTight = 263; CutIndex_Flip=21; PredictionAndControlPlot(InputPattern, "Data7TeV", CutIndex, CutIndex_Flip); PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); InputPattern = "Results/Type5/"; CutIndex = 48; CutIndexTight = 48; CutIndex_Flip=2; InitdEdx("dedxRASmi"); PredictionAndControlPlot(InputPattern, "Data7TeV", CutIndex, CutIndex_Flip); PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); return; Make2DPlot_Special("Results/Type0/", "Results/Type5/"); CompareRecoAndGenPt("Results/Type0/"); InputPattern = "Results/Type0/"; CutIndex = 4; CutIndexTight = 84; //set of cuts from the array, 0 means no cut Make2DPlot_Core(InputPattern, 0); MassPrediction(InputPattern, CutIndex, "Mass", true, "8TeV_Loose"); MassPrediction(InputPattern, CutIndex, "Mass", true, "7TeV_Loose"); MassPrediction(InputPattern, CutIndex, "Mass", false, "8TeV_LooseNoSMMC"); MassPrediction(InputPattern, CutIndex, "Mass", false, "7TeV_LooseNoSMMC"); MassPrediction(InputPattern, CutIndexTight, "Mass", true, "8TeV_Tight"); MassPrediction(InputPattern, CutIndexTight, "Mass", true, "7TeV_Tight"); MassPrediction(InputPattern, CutIndexTight, "Mass", false, "8TeV_TightNoSMMC"); MassPrediction(InputPattern, CutIndexTight, "Mass", false, "7TeV_TightNoSMMC"); PredictionAndControlPlot(InputPattern, "Data7TeV", CutIndex, CutIndex_Flip); PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); CutFlow(InputPattern, CutIndex); SelectionPlot(InputPattern, CutIndex, CutIndexTight); InputPattern = "Results/Type2/"; CutIndex = 16; CutIndexTight = 905; CutIndex_Flip=16; Make2DPlot_Core(InputPattern, 0); MassPrediction(InputPattern, CutIndex, "Mass", true, "8TeV_Loose"); MassPrediction(InputPattern, CutIndex, "Mass", true, "7TeV_Loose"); MassPrediction(InputPattern, CutIndex, "Mass", false, "8TeV_LooseNoSMMC"); MassPrediction(InputPattern, CutIndex, "Mass", false, "7TeV_LooseNoSMMC"); MassPrediction(InputPattern, CutIndexTight, "Mass", true, "8TeV_Tight"); MassPrediction(InputPattern, CutIndexTight, "Mass", true, "7TeV_Tight"); MassPrediction(InputPattern, CutIndexTight, "Mass", false, "8TeV_TightNoSMMC"); MassPrediction(InputPattern, CutIndexTight, "Mass", false, "7TeV_TightNoSMMC"); MassPrediction(InputPattern, 1, "Mass_Flip", true, "8TeV_Loose"); MassPrediction(InputPattern, 1, "Mass_Flip", true, "7TeV_Loose"); MassPrediction(InputPattern, CutIndex_Flip, "Mass_Flip", true, "8TeV_Tight"); MassPrediction(InputPattern, CutIndex_Flip, "Mass_Flip", true, "7TeV_Tight"); MassPrediction(InputPattern, CutIndex_Flip, "Mass_Flip", false, "8TeV_TightNoSMMC"); MassPrediction(InputPattern, CutIndex_Flip, "Mass_Flip", false, "7TeV_TightNoSMMC"); PredictionAndControlPlot(InputPattern, "Data7TeV", CutIndex, CutIndex_Flip); PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); CutFlow(InputPattern, CutIndex); SelectionPlot(InputPattern, CutIndex, CutIndexTight); GetSystematicOnPrediction(InputPattern, "Data7TeV"); GetSystematicOnPrediction(InputPattern, "Data8TeV"); CheckPrediction(InputPattern, "_Flip", "Data7TeV"); CheckPrediction(InputPattern, "_Flip", "Data8TeV"); InputPattern = "Results/Type3/"; CutIndex = 96; CutIndexTight = 96; CutIndex_Flip=54; PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); Make2DPlot_Core(InputPattern, 0); CutFlow(InputPattern, CutIndex); SelectionPlot(InputPattern, CutIndex, CutIndexTight); CosmicBackgroundSystematic(InputPattern, "8TeV"); CheckPrediction(InputPattern, "", "Data8TeV"); CheckPrediction(InputPattern, "_Flip", "Data8TeV"); CollisionBackgroundSystematicFromFlip(InputPattern, "Data8TeV"); //CheckPredictionBin(InputPattern, "_Flip", "Data8TeV", "0"); //CheckPredictionBin(InputPattern, "_Flip", "Data8TeV", "1"); //CheckPredictionBin(InputPattern, "_Flip", "Data8TeV", "2"); //CheckPredictionBin(InputPattern, "_Flip", "Data8TeV", "3"); //CheckPredictionBin(InputPattern, "_Flip", "Data8TeV", "4"); //CheckPredictionBin(InputPattern, "_Flip", "Data8TeV", "5"); InputPattern = "Results/Type4/"; CutIndex = 21; CutIndexTight = 263; CutIndex_Flip=21; Make2DPlot_Core(InputPattern, 0); PredictionAndControlPlot(InputPattern, "Data7TeV", CutIndex, CutIndex_Flip); PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); CheckPrediction(InputPattern, "", "Data7TeV"); CheckPrediction(InputPattern, "_Flip", "Data7TeV"); CheckPrediction(InputPattern, "", "Data8TeV"); CheckPrediction(InputPattern, "_Flip", "Data8TeV"); CollisionBackgroundSystematicFromFlip(InputPattern, "Data7TeV"); CollisionBackgroundSystematicFromFlip(InputPattern, "Data8TeV"); SelectionPlot(InputPattern, CutIndex, CutIndexTight); CutFlow(InputPattern, CutIndexTight); // CompareRecoAndGenPt(InputPattern); InputPattern = "Results/Type5/"; CutIndex = 48; CutIndexTight = 48; CutIndex_Flip=2; InitdEdx("dedxRASmi"); Make2DPlot_Core(InputPattern, 0); PredictionAndControlPlot(InputPattern, "Data7TeV", CutIndex, CutIndex_Flip); PredictionAndControlPlot(InputPattern, "Data8TeV", CutIndex, CutIndex_Flip); SelectionPlot(InputPattern, CutIndex, CutIndexTight); CutFlow(InputPattern); //This function has not yet been reviewed after july's update // MakeExpLimitpLot("Results_1toys_lp/dedxASmi/combined/Eta15/PtMin35/Type0/EXCLUSION/Stop200.info","tmp1.png"); } ////////////////////////////////////////////////////////////////////////////////////////////////////////// // General code for validation plots and final paper plots production // Make the plot of the mass distibution: Observed, data-driven prediciton and signal expectation void MassPrediction(string InputPattern, unsigned int CutIndex, string HistoSuffix, bool showMC, string DataName){ if(DataName.find("7TeV")!=string::npos){SQRTS=7.0;}else{SQRTS=8.0;} bool IsTkOnly = (InputPattern.find("Type0",0)<std::string::npos); double SystError = 0.20; TH1D *Pred8TeV=NULL, *Data8TeV=NULL, *Pred7TeV=NULL, *Data7TeV=NULL, *MCPred=NULL, *MC=NULL, *Signal=NULL; TFile* InputFile = new TFile((InputPattern + "/Histos.root").c_str()); if(!InputFile || InputFile->IsZombie() || !InputFile->IsOpen() || InputFile->TestBit(TFile::kRecovered) )return; string SName,SLeg; //README: Comments or uncomment lines below in order to decide what you want to see on your plot if(DataName.find("8TeV")!=string::npos){ SName="Gluino_8TeV_M1000_f10"; SLeg="Gluino (M = 1000 GeV/#font[12]{c}^{2})"; if(!IsTkOnly){SName="GMStau_8TeV_M308"; SLeg="Stau (M = 308 GeV/#font[12]{c}^{2})";} Pred8TeV = ((TH2D*)GetObjectFromPath(InputFile, string("Data8TeV/Pred_") + HistoSuffix ))->ProjectionY("TmpPredMass" ,CutIndex+1,CutIndex+1,"o"); Data8TeV = ((TH2D*)GetObjectFromPath(InputFile, string("Data8TeV/" ) + HistoSuffix ))->ProjectionY("TmpDataMass" ,CutIndex+1,CutIndex+1,"o"); // Pred7TeV = ((TH2D*)GetObjectFromPath(InputFile, string("Data7TeV/Pred_") + HistoSuffix ))->ProjectionY("TmpPred7TeVMass" ,CutIndex+1,CutIndex+1,"o"); // Data7TeV = ((TH2D*)GetObjectFromPath(InputFile, string("Data7TeV/" ) + HistoSuffix ))->ProjectionY("TmpData7TeVMass" ,CutIndex+1,CutIndex+1,"o"); if(showMC)MCPred = ((TH2D*)GetObjectFromPath(InputFile, string("MCTr_8TeV/Pred_" ) + HistoSuffix ))->ProjectionY("TmpMCPred" ,CutIndex+1,CutIndex+1,"o"); if(showMC)MC = ((TH2D*)GetObjectFromPath(InputFile, string("MCTr_8TeV/" ) + HistoSuffix ))->ProjectionY("TmpMCMass" ,CutIndex+1,CutIndex+1,"o"); Signal = ((TH2D*)GetObjectFromPath(InputFile, string(SName+"/" ) + HistoSuffix ))->ProjectionY("TmpSignalMass" ,CutIndex+1,CutIndex+1,"o"); }else{ SName="Gluino_7TeV_M1000_f10"; SLeg="Gluino (M = 1000 GeV/#font[12]{c}^{2})"; if(!IsTkOnly){SName="GMStau_7TeV_M308"; SLeg="Stau (M = 308 GeV/#font[12]{c}^{2})";} Pred8TeV = ((TH2D*)GetObjectFromPath(InputFile, string("Data7TeV/Pred_") + HistoSuffix ))->ProjectionY("TmpPredMass" ,CutIndex+1,CutIndex+1,"o"); Data8TeV = ((TH2D*)GetObjectFromPath(InputFile, string("Data7TeV/" ) + HistoSuffix ))->ProjectionY("TmpDataMass" ,CutIndex+1,CutIndex+1,"o"); // Pred7TeV = ((TH2D*)GetObjectFromPath(InputFile, string("Data7TeV/Pred_") + HistoSuffix ))->ProjectionY("TmpPred7TeVMass" ,CutIndex+1,CutIndex+1,"o"); // Data7TeV = ((TH2D*)GetObjectFromPath(InputFile, string("Data7TeV/" ) + HistoSuffix ))->ProjectionY("TmpData7TeVMass" ,CutIndex+1,CutIndex+1,"o"); if(showMC)MCPred = ((TH2D*)GetObjectFromPath(InputFile, string("MCTr_7TeV/Pred_" ) + HistoSuffix ))->ProjectionY("TmpMCPred" ,CutIndex+1,CutIndex+1,"o"); if(showMC)MC = ((TH2D*)GetObjectFromPath(InputFile, string("MCTr_7TeV/" ) + HistoSuffix ))->ProjectionY("TmpMCMass" ,CutIndex+1,CutIndex+1,"o"); Signal = ((TH2D*)GetObjectFromPath(InputFile, string(SName+"/" ) + HistoSuffix ))->ProjectionY("TmpSignalMass" ,CutIndex+1,CutIndex+1,"o"); } //rescale MC samples to prediction if(Data8TeV && Pred8TeV){ //if(Data7TeV && Pred7TeV)Data7TeV->Scale(Pred8TeV->Integral()/Pred7TeV->Integral()); //if(Pred7TeV) Pred7TeV->Scale(Pred8TeV->Integral()/Pred7TeV->Integral()); if(MC) MC ->Scale(Pred8TeV->Integral()/MCPred->Integral()); if(MCPred) MCPred->Scale(Pred8TeV->Integral()/MCPred->Integral()); } if(Data7TeV && Pred7TeV){ if(MC) MC ->Scale(Pred7TeV->Integral()/MCPred->Integral()); if(MCPred) MCPred->Scale(Pred7TeV->Integral()/MCPred->Integral()); } //compute integral for few mass window if(Data8TeV && Pred8TeV){ for(double M=0;M<=1000;M+=100){ if(M>400 && (int)M%200!=0)continue; double D = Data8TeV->Integral( Data8TeV->GetXaxis()->FindBin(M), Data8TeV->GetXaxis()->FindBin(2000.0)); double P = Pred8TeV->Integral( Pred8TeV->GetXaxis()->FindBin(M), Pred8TeV->GetXaxis()->FindBin(2000.0)); double Perr = 0; for(int i=Pred8TeV->GetXaxis()->FindBin(M);i<Pred8TeV->GetXaxis()->FindBin(2000.0);i++){ Perr += pow(Pred8TeV->GetBinError(i),2); } Perr = sqrt(Perr); printf("%4.0f<M<2000 --> Obs=%9.3f Data-Pred = %9.3f +- %8.3f(syst+stat) %9.3f (syst) %9.3f (stat)\n", M, D, P, sqrt(Perr*Perr + pow(P*SystError,2)), P*SystError, Perr); } } //Rebin the histograms and find who is the highest double Max = 1.0; double Min=0.01; if(Data8TeV){Data8TeV->Rebin(4); Max=std::max(Max, Data8TeV->GetMaximum());} if(Pred8TeV){Pred8TeV->Rebin(4); Max=std::max(Max, Pred8TeV->GetMaximum());} if(Data7TeV){Data7TeV->Rebin(4); Max=std::max(Max, Data7TeV->GetMaximum());} if(Pred7TeV){Pred7TeV->Rebin(4); Max=std::max(Max, Pred7TeV->GetMaximum());} if(Signal){Signal->Rebin(4); Max=std::max(Max, Signal->GetMaximum());} if(MC) {MC ->Rebin(4); Max=std::max(Max, MC ->GetMaximum());} if(MCPred){MCPred->Rebin(4); Max=std::max(Max, MCPred->GetMaximum());} Max*=5.0; //compute error bands associated to the predictions TH1D *Pred8TeVErr=NULL, *Pred7TeVErr=NULL, *PredMCErr=NULL; if(Pred8TeV)Pred8TeVErr = (TH1D*) Pred8TeV->Clone("Pred8TeVErr"); if(Pred7TeV)Pred7TeVErr = (TH1D*) Pred7TeV->Clone("Pred7TeVErr"); if(MCPred)PredMCErr = (TH1D*) MCPred->Clone("PredMCErr"); if(Pred8TeV){for(unsigned int i=0;i<(unsigned int)Pred8TeV->GetNbinsX();i++){ double error = sqrt(pow(Pred8TeVErr->GetBinError(i),2) + pow(Pred8TeVErr->GetBinContent(i)*SystError,2)); Pred8TeVErr->SetBinError(i,error); if(Pred8TeVErr->GetBinContent(i)<Min && i>5){for(unsigned int j=i+1;j<(unsigned int)Pred8TeVErr->GetNbinsX();j++)Pred8TeVErr->SetBinContent(j,0);} }} if(Pred7TeV){for(unsigned int i=0;i<(unsigned int)Pred7TeV->GetNbinsX();i++){ double error = sqrt(pow(Pred7TeVErr->GetBinError(i),2) + pow(Pred7TeVErr->GetBinContent(i)*SystError,2)); Pred7TeVErr->SetBinError(i,error); if(Pred7TeVErr->GetBinContent(i)<Min && i>5){for(unsigned int j=i+1;j<(unsigned int)Pred7TeVErr->GetNbinsX();j++)Pred7TeVErr->SetBinContent(j,0);} }} if(MCPred){for(unsigned int i=0;i<(unsigned int)MCPred->GetNbinsX();i++){ double error = sqrt(pow(PredMCErr->GetBinError(i),2) + pow(PredMCErr->GetBinContent(i)*SystError,2)); PredMCErr->SetBinError(i,error); if(PredMCErr->GetBinContent(i)<Min && i>5){for(unsigned int j=i+1;j<(unsigned int)PredMCErr->GetNbinsX();j++)PredMCErr->SetBinContent(j,0);} }} TH1D *Pred8TeVR=NULL, *Pred7TeVR=NULL, *PredMCR=NULL; if(Pred8TeV){ Pred8TeVR = (TH1D*)Pred8TeV->Clone("Pred8TeVR"); Pred8TeVR->Reset(); for(unsigned int i=0;i<(unsigned int)Pred8TeV->GetNbinsX();i++){ double Perr=0; double P = Pred8TeV->IntegralAndError(i,Pred8TeV->GetNbinsX()+1, Perr); if(P<=0)continue; double Derr=0; double D = Data8TeV->IntegralAndError(i,Data8TeV->GetNbinsX()+1, Derr); Perr = sqrt(Perr*Perr + pow(P*SystError,2)); Pred8TeVR->SetBinContent(i, D/P); Pred8TeVR->SetBinError (i, sqrt(pow(Derr*P,2) +pow(Perr*D,2))/pow(P,2)); }} if(Pred7TeV){ Pred7TeVR = (TH1D*)Pred7TeV->Clone("Pred7TeVR"); Pred7TeVR->Reset(); for(unsigned int i=0;i<(unsigned int)Pred7TeV->GetNbinsX();i++){ double Perr=0; double P = Pred7TeV->IntegralAndError(i,Pred7TeV->GetNbinsX()+1, Perr); if(P<=0)continue; double Derr=0; double D = Data7TeV->IntegralAndError(i,Data7TeV->GetNbinsX()+1, Derr); Perr = sqrt(Perr*Perr + pow(P*SystError,2)); Pred7TeVR->SetBinContent(i, D/P); Pred7TeVR->SetBinError (i, sqrt(pow(Derr*P,2) +pow(Perr*D,2))/pow(P,2)); }} if(MCPred){ PredMCR = (TH1D*)Pred8TeV->Clone("PredMCR"); PredMCR->Reset(); for(unsigned int i=0;i<(unsigned int)MCPred->GetNbinsX();i++){ double Perr=0; double P = MCPred->IntegralAndError(i,MCPred->GetNbinsX()+1, Perr); if(P<=0)continue; double Derr=0; double D = MC ->IntegralAndError(i,MC ->GetNbinsX()+1, Derr); Perr = sqrt(Perr*Perr + pow(P*SystError,2)); PredMCR->SetBinContent(i, D/P); PredMCR->SetBinError (i, sqrt(pow(Derr*P,2) +pow(Perr*D,2))/pow(P,2)); }} //Prepare the canvas for drawing and draw everything on it std::vector<string> legend; TLegend* leg; TCanvas* c1 = new TCanvas("c1","c1,",600,600); char YAxisLegend[1024]; sprintf(YAxisLegend,"Tracks / %2.0f GeV/#font[12]{c}^{2}",(Data8TeV!=NULL?Data8TeV:Data7TeV)->GetXaxis()->GetBinWidth(1)); //Loop twice to make plots with and without ratio box for(int r=0; r<2; r++) { if(r==1) { TPad* t1 = new TPad("t1","t1", 0.0, 0.20, 1.0, 1.0); t1->Draw(); t1->cd(); t1->SetLogy(true); t1->SetTopMargin(0.06); } TH1D* frame = new TH1D("frame", "frame", 1,0,1400); frame->GetXaxis()->SetNdivisions(505); frame->SetTitle(""); frame->SetStats(kFALSE); frame->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); frame->GetYaxis()->SetTitle(YAxisLegend); frame->GetYaxis()->SetTitleOffset(1.40); frame->SetMaximum(Max); frame->SetMinimum(Min); frame->SetAxisRange(0,1400,"X"); frame->Draw("AXIS"); if(Signal){ Signal->SetMarkerStyle(21); Signal->SetMarkerColor(8); Signal->SetMarkerSize(1.5); Signal->SetLineColor(1); Signal->SetFillColor(8); Signal->Draw("same HIST"); } if(MCPred){ // PredMCErr->SetLineStyle(0); PredMCErr->SetLineColor(38); PredMCErr->SetFillColor(38); PredMCErr->SetFillStyle(1001); PredMCErr->SetMarkerStyle(23); PredMCErr->SetMarkerColor(9); PredMCErr->SetMarkerSize(1.0); PredMCErr->Draw("same E5"); MCPred->SetMarkerStyle(29); MCPred->SetMarkerColor(9); MCPred->SetMarkerSize(1.5); MCPred->SetLineColor(9); MCPred->SetFillColor(0); MCPred->Draw("same HIST P"); } if(MC){ MC->SetFillStyle(3004); MC->SetLineColor(22); MC->SetFillColor(11); MC->SetMarkerStyle(0); MC->Draw("same HIST E1"); } if(Pred7TeV){ Pred7TeVErr->SetLineColor(7); Pred7TeVErr->SetFillColor(7); Pred7TeVErr->SetFillStyle(1001); Pred7TeVErr->SetMarkerStyle(22); Pred7TeVErr->SetMarkerColor(7); Pred7TeVErr->SetMarkerSize(1.0); Pred7TeVErr->Draw("same E5"); Pred7TeV->SetMarkerStyle(26); Pred7TeV->SetMarkerColor(2); Pred7TeV->SetMarkerSize(1.5); Pred7TeV->SetLineColor(2); Pred7TeV->SetFillColor(0); Pred7TeV->Draw("same HIST P"); } if(Data7TeV){ Data7TeV->SetBinContent(Data7TeV->GetNbinsX(), Data7TeV->GetBinContent(Data7TeV->GetNbinsX()) + Data7TeV->GetBinContent(Data7TeV->GetNbinsX()+1)); Data7TeV->SetMarkerStyle(24); Data7TeV->SetMarkerColor(1); Data7TeV->SetMarkerSize(1.0); Data7TeV->SetLineColor(1); Data7TeV->SetFillColor(0); Data7TeV->Draw("E1 same"); } if(Pred8TeV){ Pred8TeVErr->SetLineColor(5); Pred8TeVErr->SetFillColor(5); Pred8TeVErr->SetFillStyle(1001); Pred8TeVErr->SetMarkerStyle(22); Pred8TeVErr->SetMarkerColor(5); Pred8TeVErr->SetMarkerSize(1.0); Pred8TeVErr->Draw("same E5"); Pred8TeV->SetMarkerStyle(22); Pred8TeV->SetMarkerColor(2); Pred8TeV->SetMarkerSize(1.5); Pred8TeV->SetLineColor(2); Pred8TeV->SetFillColor(0); Pred8TeV->Draw("same HIST P"); } if(Data8TeV){ Data8TeV->SetBinContent(Data8TeV->GetNbinsX(), Data8TeV->GetBinContent(Data8TeV->GetNbinsX()) + Data8TeV->GetBinContent(Data8TeV->GetNbinsX()+1)); Data8TeV->SetMarkerStyle(20); Data8TeV->SetMarkerColor(1); Data8TeV->SetMarkerSize(1.0); Data8TeV->SetLineColor(1); Data8TeV->SetFillColor(0); Data8TeV->Draw("E1 same"); } //Fill the legend leg = new TLegend(0.7,0.93,0.35,showMC?0.66:0.75); // leg->SetHeader(LegendFromType(InputPattern).c_str()); leg->SetFillStyle(0); leg->SetBorderSize(0); leg->SetTextFont(43); leg->SetTextSize(20); if(Data8TeV){leg->AddEntry(Data8TeV, "Observed" ,"P");} if(Pred8TeV){TH1D* PredLeg8TeV = (TH1D*)Pred8TeV->Clone("RescLeg12"); PredLeg8TeV->SetFillColor(Pred8TeVErr->GetFillColor()); PredLeg8TeV->SetFillStyle(Pred8TeVErr->GetFillStyle()); leg->AddEntry(PredLeg8TeV, "Data-based SM prediction" ,"PF"); } if(Data7TeV){leg->AddEntry(Data7TeV, "Observed (2011)" ,"P");} if(Pred7TeV){TH1D* PredLeg7TeV = (TH1D*)Pred7TeV->Clone("RescLeg11"); PredLeg7TeV->SetFillColor(Pred7TeVErr->GetFillColor()); PredLeg7TeV->SetFillStyle(Pred7TeVErr->GetFillStyle()); leg->AddEntry(PredLeg7TeV, "Data-based SM prediction (2011)" ,"PF"); } if(MC ){leg->AddEntry(MC, "Simulation" ,"LF");} if(MCPred){TH1D* MCPredLeg = (TH1D*) MCPred->Clone("RescMCLeg"); MCPredLeg->SetFillColor(PredMCErr->GetFillColor()); MCPredLeg->SetFillStyle(PredMCErr->GetFillStyle()); leg->AddEntry(MCPredLeg, "SM prediction (MC)" ,"PF"); } if(Signal)leg->AddEntry(Signal, SLeg.c_str() ,"F"); leg->Draw(); //Redraw axis ticks gPad->RedrawAxis(); //add CMS label and save DrawPreliminary(LegendFromType(InputPattern), SQRTS, IntegratedLuminosityFromE(SQRTS), false); c1->SetLogy(true); c1->cd(); if(r==0) SaveCanvas(c1, InputPattern, string("RescaleNoRatio_") + HistoSuffix + "_" + DataName); else { c1->cd(); TPad* t2 = new TPad("t2","t2", 0.0, 0.0, 1.0, 0.2); t2->Draw(); t2->cd(); t2->SetGridy(true); t2->SetPad(0,0.0,1.0,0.2); t2->SetTopMargin(0); t2->SetBottomMargin(0.5); TH1D* frameR = new TH1D("frameR", "frameR", 1,0, 1400); frameR->GetXaxis()->SetNdivisions(505); frameR->SetTitle(""); frameR->SetStats(kFALSE); frameR->GetXaxis()->SetTitle(""); frameR->GetXaxis()->SetTitle("Mass (GeV/#font[12]{c}^{2})"); frameR->GetYaxis()->SetTitle("Obs / Pred"); frameR->SetMaximum(2.00); frameR->SetMinimum(0.00); frameR->GetYaxis()->SetLabelFont(43); //give the font size in pixel (instead of fraction) frameR->GetYaxis()->SetLabelSize(15); //font size frameR->GetYaxis()->SetTitleFont(43); //give the font size in pixel (instead of fraction) frameR->GetYaxis()->SetTitleSize(15); //font size frameR->GetYaxis()->SetNdivisions(505); frameR->GetXaxis()->SetLabelFont(43); //give the font size in pixel (instead of fraction) frameR->GetXaxis()->SetLabelSize(15); //font size frameR->GetXaxis()->SetTitleFont(43); //give the font size in pixel (instead of fraction) frameR->GetXaxis()->SetTitleSize(15); //font size frameR->GetXaxis()->SetTitleOffset(3.75); frameR->Draw("AXIS"); if(MCPred){ PredMCR->SetMarkerStyle(29); PredMCR->SetMarkerColor(9); PredMCR->SetMarkerSize(1.5); PredMCR->SetLineColor(9); PredMCR->SetFillColor(0); PredMCR->Draw("same E1"); } if(Pred7TeV){ Pred7TeVR->SetMarkerStyle(26); Pred7TeVR->SetMarkerColor(2); Pred7TeVR->SetMarkerSize(1.5); Pred7TeVR->SetLineColor(2); Pred7TeVR->SetFillColor(0); Pred7TeVR->Draw("same E1"); } if(Pred8TeV){ Pred8TeVR->SetMarkerStyle(22); Pred8TeVR->SetMarkerColor(2); Pred8TeVR->SetMarkerSize(1.5); Pred8TeVR->SetLineColor(2); Pred8TeVR->SetFillColor(0); Pred8TeVR->Draw("same E1"); } TLine* LineAtOne = new TLine(0,1,1400,1); LineAtOne->SetLineStyle(3); LineAtOne->Draw(); c1->cd(); SaveCanvas(c1, InputPattern, string("Rescale_") + HistoSuffix + "_" + DataName); } } delete c1; InputFile->Close(); } // make some control plots to show that ABCD method can be used void PredictionAndControlPlot(string InputPattern, string Data, unsigned int CutIndex, unsigned int CutIndex_Flip){ if(Data.find("7TeV")!=string::npos){SQRTS=7.0;}else{SQRTS=8.0;} TCanvas* c1; TObject** Histos = new TObject*[10]; std::vector<string> legend; TypeMode = TypeFromPattern(InputPattern); string LegendTitle = LegendFromType(InputPattern);; TFile* InputFile = new TFile((InputPattern + "/Histos.root").c_str()); TH1D* CtrlPt_S1_Is = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S1_Is" ); CtrlPt_S1_Is ->Rebin(5); TH1D* CtrlPt_S1_Im = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S1_Im" ); CtrlPt_S1_Im ->Rebin(1); TH1D* CtrlPt_S1_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S1_TOF"); CtrlPt_S1_TOF->Rebin(1); TH1D* CtrlPt_S2_Is = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S2_Is" ); CtrlPt_S2_Is ->Rebin(5); TH1D* CtrlPt_S2_Im = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S2_Im" ); CtrlPt_S2_Im ->Rebin(1); TH1D* CtrlPt_S2_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S2_TOF"); CtrlPt_S2_TOF->Rebin(1); TH1D* CtrlPt_S3_Is = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S3_Is" ); CtrlPt_S3_Is ->Rebin(5); TH1D* CtrlPt_S3_Im = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S3_Im" ); CtrlPt_S3_Im ->Rebin(1); TH1D* CtrlPt_S3_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S3_TOF"); CtrlPt_S3_TOF->Rebin(1); TH1D* CtrlPt_S4_Is = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S4_Is" ); CtrlPt_S4_Is ->Rebin(5); TH1D* CtrlPt_S4_Im = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S4_Im" ); CtrlPt_S4_Im ->Rebin(1); TH1D* CtrlPt_S4_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S4_TOF"); CtrlPt_S4_TOF->Rebin(1); TH1D* CtrlIs_S1_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlIs_S1_TOF"); CtrlIs_S1_TOF->Rebin(1); TH1D* CtrlIs_S2_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlIs_S2_TOF"); CtrlIs_S2_TOF->Rebin(1); TH1D* CtrlIs_S3_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlIs_S3_TOF"); CtrlIs_S3_TOF->Rebin(1); TH1D* CtrlIs_S4_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlIs_S4_TOF"); CtrlIs_S4_TOF->Rebin(1); TH1D* CtrlIm_S1_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlIm_S1_TOF"); CtrlIm_S1_TOF->Rebin(1); TH1D* CtrlIm_S2_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlIm_S2_TOF"); CtrlIm_S2_TOF->Rebin(1); TH1D* CtrlIm_S3_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlIm_S3_TOF"); CtrlIm_S3_TOF->Rebin(1); TH1D* CtrlIm_S4_TOF = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlIm_S4_TOF"); CtrlIm_S4_TOF->Rebin(1); TH1D* CtrlPt_S1_TOF_Binned[MaxPredBins]; TH1D* CtrlPt_S2_TOF_Binned[MaxPredBins]; TH1D* CtrlPt_S3_TOF_Binned[MaxPredBins]; TH1D* CtrlPt_S4_TOF_Binned[MaxPredBins]; if(TypeMode==3) PredBins=6; else PredBins=0; for(int i=0; i<PredBins; i++) { char Suffix[1024]; sprintf(Suffix,"_%i",i); string Bin=Suffix; CtrlPt_S1_TOF_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S1_TOF_Binned"+Bin); CtrlPt_S1_TOF_Binned[i]->Rebin(1); CtrlPt_S2_TOF_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S2_TOF_Binned"+Bin); CtrlPt_S2_TOF_Binned[i]->Rebin(1); CtrlPt_S3_TOF_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S3_TOF_Binned"+Bin); CtrlPt_S3_TOF_Binned[i]->Rebin(1); CtrlPt_S4_TOF_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, Data+"/CtrlPt_S4_TOF_Binned"+Bin); CtrlPt_S4_TOF_Binned[i]->Rebin(1); } std::vector<std::string> PtLimitsNames; if(TypeMode!=3) { PtLimitsNames.push_back(" 45 < p_{T} < 60 GeV/#font[12]{c}"); PtLimitsNames.push_back(" 60 < p_{T} < 80 GeV/#font[12]{c}"); PtLimitsNames.push_back(" 80 < p_{T} < 100 GeV/#font[12]{c}"); PtLimitsNames.push_back("100 < p_{T} GeV/#font[12]{c}"); } else { PtLimitsNames.push_back(" 80 < p_{T} < 120 GeV/#font[12]{c}"); PtLimitsNames.push_back("120 < p_{T} < 170 GeV/#font[12]{c}"); PtLimitsNames.push_back("170 < p_{T} < 240 GeV/#font[12]{c}"); PtLimitsNames.push_back("240 < p_{T} GeV/#font[12]{c}"); } c1 = new TCanvas("c1","c1,",600,600); legend.clear(); if(CtrlPt_S1_Is->Integral()>0)CtrlPt_S1_Is->Scale(1/CtrlPt_S1_Is->Integral()); if(CtrlPt_S2_Is->Integral()>0)CtrlPt_S2_Is->Scale(1/CtrlPt_S2_Is->Integral()); if(CtrlPt_S3_Is->Integral()>0)CtrlPt_S3_Is->Scale(1/CtrlPt_S3_Is->Integral()); if(CtrlPt_S4_Is->Integral()>0)CtrlPt_S4_Is->Scale(1/CtrlPt_S4_Is->Integral()); Histos[0] = CtrlPt_S1_Is; legend.push_back(PtLimitsNames[0]); Histos[1] = CtrlPt_S2_Is; legend.push_back(PtLimitsNames[1]); Histos[2] = CtrlPt_S3_Is; legend.push_back(PtLimitsNames[2]); Histos[3] = CtrlPt_S4_Is; legend.push_back(PtLimitsNames[3]); char YAxisTitle[100]; sprintf(YAxisTitle,"Fraction of tracks / %0.3f",((TH1D*)Histos[0])->GetBinWidth(1)); DrawSuperposedHistos((TH1**)Histos, legend, "E1", dEdxS_Legend, YAxisTitle,TypeMode!=5?0:0.7,TypeMode!=5?0.5:1.0, 0,0); DrawLegend(Histos,legend,"", "P", 0.77, 0.92, 0.43, 0.05); c1->SetLogy(true); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS), false); SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Pt_IsSpectrum"); delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); if(CtrlPt_S1_Im->Integral()>0)CtrlPt_S1_Im->Scale(1/CtrlPt_S1_Im->Integral()); if(CtrlPt_S2_Im->Integral()>0)CtrlPt_S2_Im->Scale(1/CtrlPt_S2_Im->Integral()); if(CtrlPt_S3_Im->Integral()>0)CtrlPt_S3_Im->Scale(1/CtrlPt_S3_Im->Integral()); if(CtrlPt_S4_Im->Integral()>0)CtrlPt_S4_Im->Scale(1/CtrlPt_S4_Im->Integral()); Histos[0] = CtrlPt_S1_Im; legend.push_back(PtLimitsNames[0]); Histos[1] = CtrlPt_S2_Im; legend.push_back(PtLimitsNames[1]); Histos[2] = CtrlPt_S3_Im; legend.push_back(PtLimitsNames[2]); Histos[3] = CtrlPt_S4_Im; legend.push_back(PtLimitsNames[3]); DrawSuperposedHistos((TH1**)Histos, legend, "E1", dEdxM_Legend, "arbitrary units", 3.0,5, 0,0); DrawLegend(Histos,legend,"","P", 0.77, 0.92, 0.43, 0.05); c1->SetLogy(true); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Pt_ImSpectrum"); delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); if(CtrlPt_S1_TOF->Integral()>0)CtrlPt_S1_TOF->Scale(1/CtrlPt_S1_TOF->Integral()); if(CtrlPt_S2_TOF->Integral()>0)CtrlPt_S2_TOF->Scale(1/CtrlPt_S2_TOF->Integral()); if(CtrlPt_S3_TOF->Integral()>0)CtrlPt_S3_TOF->Scale(1/CtrlPt_S3_TOF->Integral()); if(CtrlPt_S4_TOF->Integral()>0)CtrlPt_S4_TOF->Scale(1/CtrlPt_S4_TOF->Integral()); Histos[0] = CtrlPt_S1_TOF; legend.push_back(PtLimitsNames[0]); Histos[1] = CtrlPt_S2_TOF; legend.push_back(PtLimitsNames[1]); Histos[2] = CtrlPt_S3_TOF; legend.push_back(PtLimitsNames[2]); Histos[3] = CtrlPt_S4_TOF; legend.push_back(PtLimitsNames[3]); sprintf(YAxisTitle,"Fraction of tracks / %0.3f",((TH1D*)Histos[0])->GetBinWidth(1)); // DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta", "fraction of tracks", 0.5,2.2, 1E-6,1); DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta", YAxisTitle, 1.0, 1.4, 0, 0); DrawLegend(Histos,legend, "" ,"P", 0.80, 0.92, 0.43, 0.05); c1->SetLogy(true); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS), false); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Pt_TOFSpectrum"); c1->SetLogy(false); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Pt_TOFSpectrumNoLog"); delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); c1->SetLogy(true); if(CtrlIs_S1_TOF->Integral()>0)CtrlIs_S1_TOF->Scale(1/CtrlIs_S1_TOF->Integral()); if(CtrlIs_S2_TOF->Integral()>0)CtrlIs_S2_TOF->Scale(1/CtrlIs_S2_TOF->Integral()); if(CtrlIs_S3_TOF->Integral()>0)CtrlIs_S3_TOF->Scale(1/CtrlIs_S3_TOF->Integral()); if(CtrlIs_S4_TOF->Integral()>0)CtrlIs_S4_TOF->Scale(1/CtrlIs_S4_TOF->Integral()); Histos[0] = CtrlIs_S1_TOF; legend.push_back("0.00 < I_{as} < 0.05"); Histos[1] = CtrlIs_S2_TOF; legend.push_back("0.05 < I_{as} < 0.10"); Histos[2] = CtrlIs_S3_TOF; legend.push_back("0.10 < I_{as} < 0.20"); Histos[3] = CtrlIs_S4_TOF; legend.push_back("0.20 < I_{as}"); sprintf(YAxisTitle,"Fraction of tracks / %0.3f",((TH1D*)Histos[0])->GetBinWidth(1)); DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta", YAxisTitle, 1,1.4, 0.000005,0.6); CtrlIs_S4_TOF->Draw("E1 same"); //redraw this histogram to make sure it is on top of the other ones DrawLegend(Histos,legend, "","P", 0.77, 0.92, 0.43, 0.05); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS), false); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Is_TOFSpectrumLog"); c1->SetLogy(false); // DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); // DrawLegend(Histos,legend, "","P"); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Is_TOFSpectrum"); delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); if(CtrlIm_S1_TOF->Integral()>0)CtrlIm_S1_TOF->Scale(1/CtrlIm_S1_TOF->Integral()); if(CtrlIm_S2_TOF->Integral()>0)CtrlIm_S2_TOF->Scale(1/CtrlIm_S2_TOF->Integral()); if(CtrlIm_S3_TOF->Integral()>0)CtrlIm_S3_TOF->Scale(1/CtrlIm_S3_TOF->Integral()); if(CtrlIm_S4_TOF->Integral()>0)CtrlIm_S4_TOF->Scale(1/CtrlIm_S4_TOF->Integral()); Histos[0] = CtrlIm_S1_TOF; legend.push_back("3.5<I_{as}<3.8"); Histos[1] = CtrlIm_S2_TOF; legend.push_back("3.8<I_{as}<4.1"); Histos[2] = CtrlIm_S3_TOF; legend.push_back("4.1<I_{as}<4.4"); DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta", "arbitrary units", 1,1.7, 0,0); DrawLegend(Histos,legend,"","P", 0.77, 0.92, 0.43, 0.05); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Im_TOFSpectrum"); delete c1; for(int i=0; i<PredBins; i++) { char Suffix[1024]; sprintf(Suffix,"_%i",i); string Bin=Suffix; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); if(CtrlPt_S1_TOF_Binned[i]->Integral()>0)CtrlPt_S1_TOF_Binned[i]->Scale(1/CtrlPt_S1_TOF_Binned[i]->Integral()); if(CtrlPt_S2_TOF_Binned[i]->Integral()>0)CtrlPt_S2_TOF_Binned[i]->Scale(1/CtrlPt_S2_TOF_Binned[i]->Integral()); if(CtrlPt_S3_TOF_Binned[i]->Integral()>0)CtrlPt_S3_TOF_Binned[i]->Scale(1/CtrlPt_S3_TOF_Binned[i]->Integral()); if(CtrlPt_S4_TOF_Binned[i]->Integral()>0)CtrlPt_S4_TOF_Binned[i]->Scale(1/CtrlPt_S4_TOF_Binned[i]->Integral()); Histos[0] = CtrlPt_S1_TOF_Binned[i]; legend.push_back(PtLimitsNames[0]); Histos[1] = CtrlPt_S2_TOF_Binned[i]; legend.push_back(PtLimitsNames[1]); Histos[2] = CtrlPt_S3_TOF_Binned[i]; legend.push_back(PtLimitsNames[2]); Histos[3] = CtrlPt_S4_TOF_Binned[i]; legend.push_back(PtLimitsNames[3]); DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta", "arbitrary units", 0,2, 0,0); DrawLegend(Histos,legend,LegendTitle,"P", 0.77, 0.92, 0.43, 0.05); c1->SetLogy(true); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Pt_TOFSpectrum_Binned"+Bin); c1->SetLogy(false); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Control_")+Data+"_Pt_TOFSpectrumNoLog_Binned"+Bin); delete c1; } if(TypeMode<3) {//These plots only made for analyses using mass distribution //Show P, I and TOF distribution in the signal region (observed and predicted) TH2D* Pred_P = (TH2D*)GetObjectFromPath(InputFile, Data+"/Pred_P"); TH2D* Pred_I = (TH2D*)GetObjectFromPath(InputFile, Data+"/Pred_I"); TH2D* Pred_TOF = (TH2D*)GetObjectFromPath(InputFile, Data+"/Pred_TOF"); TH2D* Data_I = (TH2D*)GetObjectFromPath(InputFile, Data+"/RegionD_I"); TH2D* Data_P = (TH2D*)GetObjectFromPath(InputFile, Data+"/RegionD_P"); TH2D* Data_TOF = (TH2D*)GetObjectFromPath(InputFile, Data+"/RegionD_TOF"); c1 = new TCanvas("c1","c1,",600,600); legend.clear(); c1->SetLogy(true); Histos[0] = (TH1D*)(Data_P->ProjectionY("PA",CutIndex+1,CutIndex+1,"o")); legend.push_back("Observed"); Histos[1] = (TH1D*)(Pred_P->ProjectionY("PB",CutIndex+1,CutIndex+1,"o")); legend.push_back("Predicted"); ((TH1D*)Histos[0])->Scale(1/std::max(((TH1D*)Histos[0])->Integral(),1.0)); ((TH1D*)Histos[1])->Scale(1/std::max(((TH1D*)Histos[1])->Integral(),1.0)); ((TH1D*)Histos[0])->Rebin(10); ((TH1D*)Histos[1])->Rebin(10); DrawSuperposedHistos((TH1**)Histos, legend, "Hist E1", "p (Gev/c)", "u.a.", 0,1500, 0,0); DrawLegend(Histos,legend,"","P"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_PSpectrum"); delete Histos[0]; delete Histos[1]; delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); c1->SetLogy(true); Histos[0] = (TH1D*)(Data_I->ProjectionY("IA",CutIndex+1,CutIndex+1,"o")); legend.push_back("Observed"); Histos[1] = (TH1D*)(Pred_I->ProjectionY("IB",CutIndex+1,CutIndex+1,"o")); legend.push_back("Predicted"); ((TH1D*)Histos[0])->Scale(1/std::max(((TH1D*)Histos[0])->Integral(),1.0)); ((TH1D*)Histos[1])->Scale(1/std::max(((TH1D*)Histos[1])->Integral(),1.0)); ((TH1D*)Histos[0])->Rebin(2); ((TH1D*)Histos[1])->Rebin(2); DrawSuperposedHistos((TH1**)Histos, legend, "Hist E1", dEdxM_Legend, "u.a.", 0,6, 0,0); DrawLegend(Histos,legend,"", "P"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_ISpectrum"); delete Histos[0]; delete Histos[1]; delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); c1->SetLogy(true); Histos[0] = (TH1D*)(Data_TOF->ProjectionY("TA",CutIndex+1,CutIndex+1,"o")); legend.push_back("Observed"); Histos[1] = (TH1D*)(Pred_TOF->ProjectionY("TB",CutIndex+1,CutIndex+1,"o")); legend.push_back("Predicted"); ((TH1D*)Histos[0])->Scale(1/std::max(((TH1D*)Histos[0])->Integral(),1.0)); ((TH1D*)Histos[1])->Scale(1/std::max(((TH1D*)Histos[1])->Integral(),1.0)); ((TH1D*)Histos[0])->Rebin(2); ((TH1D*)Histos[1])->Rebin(2); DrawSuperposedHistos((TH1**)Histos, legend, "Hist E1", "1/#beta", "u.a.", 0,0, 0,0); DrawLegend(Histos,legend,"", "P"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_TOFSpectrum"); delete Histos[0]; delete Histos[1]; delete c1; //Show P, I and TOF distribution in the region with TOF < 1(observed and predicted) TH2D* Pred_P_Flip = (TH2D*)GetObjectFromPath(InputFile, Data+"/Pred_P_Flip"); TH2D* Pred_I_Flip = (TH2D*)GetObjectFromPath(InputFile, Data+"/Pred_I_Flip"); TH2D* Pred_TOF_Flip = (TH2D*)GetObjectFromPath(InputFile, Data+"/Pred_TOF_Flip"); TH2D* Data_I_Flip = (TH2D*)GetObjectFromPath(InputFile, Data+"/RegionD_I_Flip"); TH2D* Data_P_Flip = (TH2D*)GetObjectFromPath(InputFile, Data+"/RegionD_P_Flip"); TH2D* Data_TOF_Flip = (TH2D*)GetObjectFromPath(InputFile, Data+"/RegionD_TOF_Flip"); c1 = new TCanvas("c1","c1,",600,600); legend.clear(); c1->SetLogy(true); Histos[0] = (TH1D*)(Data_P_Flip ->ProjectionY("PA_Flip",CutIndex_Flip+1,CutIndex_Flip+1,"o")); legend.push_back("Observed"); Histos[1] = (TH1D*)(Pred_P_Flip ->ProjectionY("PB_Flip",CutIndex_Flip+1,CutIndex_Flip+1,"o")); legend.push_back("Predicted"); ((TH1D*)Histos[0])->Scale(1/std::max(((TH1D*)Histos[0])->Integral(),1.0)); ((TH1D*)Histos[1])->Scale(1/std::max(((TH1D*)Histos[1])->Integral(),1.0)); ((TH1D*)Histos[0])->Rebin(10); ((TH1D*)Histos[1])->Rebin(10); DrawSuperposedHistos((TH1**)Histos, legend, "Hist E1", "p (Gev/c)", "u.a.", 0,1500, 0,0); DrawLegend(Histos,legend,"","P"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_PSpectrum_Flip"); delete Histos[0]; delete Histos[1]; delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); c1->SetLogy(true); Histos[0] = (TH1D*)(Data_I_Flip ->ProjectionY("IA_Flip",CutIndex_Flip+1,CutIndex_Flip+1,"o")); legend.push_back("Observed"); Histos[1] = (TH1D*)(Pred_I_Flip ->ProjectionY("IB_Flip",CutIndex_Flip+1,CutIndex_Flip+1,"o")); legend.push_back("Predicted"); ((TH1D*)Histos[0])->Scale(1/std::max(((TH1D*)Histos[0])->Integral(),1.0)); ((TH1D*)Histos[1])->Scale(1/std::max(((TH1D*)Histos[1])->Integral(),1.0)); ((TH1D*)Histos[0])->Rebin(2); ((TH1D*)Histos[1])->Rebin(2); DrawSuperposedHistos((TH1**)Histos, legend, "Hist E1", dEdxM_Legend, "u.a.", 0,6, 0,0); DrawLegend(Histos,legend,"","P"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_ISpectrum_Flip"); delete Histos[0]; delete Histos[1]; delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); c1->SetLogy(true); Histos[0] = (TH1D*)(Data_TOF_Flip ->ProjectionY("TA_Flip",CutIndex_Flip+1,CutIndex_Flip+1,"o")); legend.push_back("Observed"); Histos[1] = (TH1D*)(Pred_TOF_Flip ->ProjectionY("TB_Flip",CutIndex_Flip+1,CutIndex_Flip+1,"o")); legend.push_back("Predicted"); ((TH1D*)Histos[0])->Scale(1/std::max(((TH1D*)Histos[0])->Integral(),1.0)); ((TH1D*)Histos[1])->Scale(1/std::max(((TH1D*)Histos[1])->Integral(),1.0)); ((TH1D*)Histos[0])->Rebin(2); ((TH1D*)Histos[1])->Rebin(2); DrawSuperposedHistos((TH1**)Histos, legend, "Hist E1", "1/#beta", "u.a.", 0,0, 0,0); DrawLegend(Histos,legend, "" ,"P"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); if(TypeMode>=2)SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_TOFSpectrum_Flip"); delete Histos[0]; delete Histos[1]; delete c1; } if (TypeMode==4) { TH1D* HCuts_I = (TH1D*)GetObjectFromPath(InputFile, "HCuts_I"); TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, "HCuts_TOF"); TH2D* Pred_Ias = (TH2D*)GetObjectFromPath(InputFile, Data+"/RegionH_Ias"); TH2D* Data_Ias = (TH2D*)GetObjectFromPath(InputFile, Data+"/RegionD_Ias"); TH1D* Data_C = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_C"); TH1D* Data_G = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_G"); TH1D* Data_H = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_H"); TH1D* Data_D = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_D"); TH1D* Data_P = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_P"); cout << " Ias cut is " << HCuts_I->GetBinContent(CutIndex+1) << " , 1/beta cut is " << HCuts_TOF->GetBinContent(CutIndex+1) << endl; cout << " C = " << Data_C->GetBinContent(CutIndex+1) << endl; cout << " G = " << Data_G->GetBinContent(CutIndex+1) << endl; cout << " H = " << Data_H->GetBinContent(CutIndex+1) << endl; cout << " D = " << Data_D->GetBinContent(CutIndex+1) << endl; cout << " Prediction is " << Data_P->GetBinContent(CutIndex+1) << " +/- " << Data_P->GetBinError(CutIndex+1) << endl; if (Data_P->GetBinContent(CutIndex+1) > Data_D->GetBinContent(CutIndex+1)){ cout << " Data is fewer by " << ((Data_P->GetBinContent(CutIndex+1)- Data_D->GetBinContent(CutIndex+1))/Data_P->GetBinContent(CutIndex+1))*100.0 << "%" << endl; } else { cout << " Prediction is fewer by " << ((Data_D->GetBinContent(CutIndex+1)-Data_P->GetBinContent(CutIndex+1))/Data_P->GetBinContent(CutIndex+1))*100.0 << "%" << endl; } double factor = Data_C->GetBinContent(CutIndex+1)/Data_G->GetBinContent(CutIndex+1) ; // cout << " factor = " << factor << endl; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); c1->SetLogy(true); Histos[0] = (TH1D*)(Data_Ias->ProjectionY("IasD",CutIndex+1,CutIndex+1,"o")); legend.push_back("Observed"); Histos[1] = (TH1D*)(Pred_Ias->ProjectionY("IasH",CutIndex+1,CutIndex+1,"o")); legend.push_back("Prediction"); ((TH1D*)Histos[1])->Scale(factor); DrawSuperposedHistos((TH1**)Histos, legend, "E1", "I_{as} ", "Tracks", 0,0, 0.05, 400000); DrawLegend(Histos,legend,"","P"); c1->SetLogy(true); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_IasSpectrum"); delete Histos[0]; delete Histos[1]; delete c1; } if(TypeMode!=2){ for(int S=0;S<2;S++){ if(TypeMode==0 && S>0)continue; string suffix = ""; if(S==1)suffix = "_Flip"; TH1D* H_D = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_D"+suffix); TH1D* H_P = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_P"+suffix); TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, "HCuts_TOF"+suffix); TH1D* HCuts_Pt = (TH1D*)GetObjectFromPath(InputFile, "HCuts_Pt"+suffix); TH1D* HCuts_I = (TH1D*)GetObjectFromPath(InputFile, "HCuts_I"+suffix); std::map<std::pair<float, float>, TGraphErrors*> mapPred; std::map<std::pair<float, float>, TGraphErrors*> mapObs; std::map<std::pair<float, float>, TGraphErrors*> mapRatio; double max = 1; double min = 1000; double xmin=100; double xmax=0; double binZero=-8, binOne=-1; for(int CutIndex_=1;CutIndex_<H_P->GetNbinsX();CutIndex_++){ float PtCut = HCuts_Pt ->GetBinContent(CutIndex_+1); //float ICut = HCuts_I ->GetBinContent(CutIndex_+1); float TCut = HCuts_TOF->GetBinContent(CutIndex_+1); //if(std::make_pair(PtCut, -1)!=P1 && std::make_pair(PtCut, -1)!=P2 && std::make_pair(PtCut, -1)!=P3 && //std::make_pair(-1, TCut)!=P1 &&std::make_pair(-1, TCut)!=P2 && std::make_pair(-1, TCut)!=P3) //cout << PtCut << endl; if(TypeMode==3) TCut=-1; if(mapPred.find(std::make_pair(PtCut, TCut))!=mapPred.end())continue; int N =0;for(int i=CutIndex_;i<H_P->GetNbinsX();i++)if(float(HCuts_Pt->GetBinContent(i+1))==PtCut && (float(HCuts_TOF->GetBinContent(i+1))==TCut || TypeMode==3))N++; mapPred [std::make_pair(PtCut, TCut)] = new TGraphErrors(N); mapObs [std::make_pair(PtCut, TCut)] = new TGraphErrors(N); mapRatio[std::make_pair(PtCut, TCut)] = new TGraphErrors(N); int N_i = 0; for(int i=CutIndex_;i<H_P->GetNbinsX();i++){ if(float(HCuts_Pt->GetBinContent(i+1))!=PtCut || (float(HCuts_TOF->GetBinContent(i+1))!=TCut && TypeMode!=3))continue; if(isnan((float)H_P->GetBinContent(i+1)))continue; if(S!=1 && H_P->GetBinContent(i+1)<=0)continue; if(TypeMode!=3) { xmin = std::min(xmin, HCuts_I->GetBinContent(i+1)); xmax = std::max(xmax, HCuts_I->GetBinContent(i+1)); if(N_i==0) binZero=HCuts_I->GetBinContent(i+1); if(N_i==1) binOne=HCuts_I->GetBinContent(i+1); } else { xmin = std::min(xmin, HCuts_TOF->GetBinContent(i+1)); xmax = std::max(xmax, HCuts_TOF->GetBinContent(i+1)); if(N_i==0) binZero=HCuts_TOF->GetBinContent(i+1); if(N_i==1) binOne=HCuts_TOF->GetBinContent(i+1); } max = std::max(max, H_P->GetBinContent(i+1)); // min = std::min(min, std::max(H_P->GetBinContent(CutIndex_+i+1), 0.1) ); double P = H_P->GetBinContent(i+1); double Perr = H_P->GetBinError(i+1); if(S==1 && P<=0){P = 3/2.0; Perr=3/2.0; max = std::max(max, 3.0);} if(TypeMode!=3) mapPred[std::make_pair(PtCut, TCut)]->SetPoint (N_i, HCuts_I->GetBinContent(i+1), P); else mapPred[std::make_pair(PtCut, TCut)]->SetPoint (N_i, HCuts_TOF->GetBinContent(i+1), P); mapPred[std::make_pair(PtCut, TCut)]->SetPointError(N_i, 0 , sqrt(pow(Perr,2) + pow(0.2*H_P->GetBinContent(i+1),2) ) ); mapPred[std::make_pair(PtCut, TCut)]->Set(N_i+1); max = std::max(max, H_D->GetBinContent(CutIndex_+i+1)); min = std::min(min, std::max(H_D->GetBinContent(CutIndex_+i+1), 0.1)); if(TypeMode!=3) mapObs [std::make_pair(PtCut, TCut)]->SetPoint (N_i, HCuts_I->GetBinContent(i+1), H_D->GetBinContent(i+1)); else mapObs [std::make_pair(PtCut, TCut)]->SetPoint (N_i, HCuts_TOF->GetBinContent(i+1), H_D->GetBinContent(i+1)); mapObs [std::make_pair(PtCut, TCut)]->SetPointError(N_i, 0 , H_D->GetBinError (i+1)); mapObs [std::make_pair(PtCut, TCut)]->Set(N_i+1); if(TypeMode!=3) mapRatio[std::make_pair(PtCut, TCut)]->SetPoint(N_i, HCuts_I->GetBinContent(i+1), (H_D->GetBinContent(i+1)<=0 || P<=0) ? 0 : H_D->GetBinContent(i+1)/P); else {mapRatio[std::make_pair(PtCut, TCut)]->SetPoint(N_i, HCuts_TOF->GetBinContent(i+1), (H_D->GetBinContent(i+1)<=0 || P<=0) ? 0 : H_D->GetBinContent(i+1) / P);} mapRatio[std::make_pair(PtCut, TCut)]->SetPointError(N_i, 0, H_D->GetBinContent(i+1)<=0 ? 0 : sqrt(pow(Perr * H_D->GetBinContent(i+1),2) + pow(P * H_D->GetBinError (i+1),2) ) / pow(P,2) ); mapRatio[std::make_pair(PtCut, TCut)]->Set(N_i+1); N_i++; } } //Save twice once with a ratio and once without for(int r=0; r<2; r++) { c1 = new TCanvas("c1","c1,",600,600); TPad* t1; if(r==1) { t1 = new TPad("t1","t1", 0.0, 0.20, 1.0, 1.0); t1->Draw(); t1->cd(); t1->SetLogy(true); t1->SetTopMargin(0.06); } TH1D* frame = new TH1D("frame", "frame", 1,std::max(xmin,0.02),xmax); frame->GetXaxis()->SetNdivisions(505); frame->SetTitle(""); frame->SetStats(kFALSE); if(TypeMode==3) frame->GetXaxis()->SetTitle("1/#beta threshold"); else frame->GetXaxis()->SetTitle((dEdxS_Legend + " threshold").c_str()); //Cumulative plot so not per anything //char YAxisTitle[100]; //sprintf(YAxisTitle,"Tracks/%0.3f",fabs(binOne-binZero)); frame->GetYaxis()->SetTitle("Tracks"); frame->GetYaxis()->SetTitleOffset(1.50); frame->SetMaximum(max*10); frame->SetMinimum(min*0.1); frame->Draw("AXIS"); TLegend* LEG; if(S==0) LEG = new TLegend(0.31,0.72,0.77,0.92); else LEG = new TLegend(0.15,0.72,0.6,0.92); LEG->SetFillColor(0); LEG->SetFillStyle(0); LEG->SetBorderSize(0); LEG->SetTextFont(43); LEG->SetTextSize(20); std::pair<float, float> P1, P2, P3; string L1, L2, L3; string Data1, Data2, Data3; if(TypeMode==0){ P1 = std::make_pair(55.0, -1.0); L1 = "Predicted (p_{T} > 55 GeV/#font[12]{c})"; P2 = std::make_pair(65.0, -1.0); L2 = "Predicted (p_{T} > 65 GeV/#font[12]{c})"; P3 = std::make_pair(75.0, -1.0); L3 = "Predicted (p_{T} > 75 GeV/#font[12]{c})"; Data1 = "Observed (p_{T} > 55 GeV/#font[12]{c})"; Data2 = "Observed (p_{T} > 65 GeV/#font[12]{c})"; Data3 = "Observed (p_{T} > 75 GeV/#font[12]{c})"; }else if(TypeMode==3){ P1 = std::make_pair(110.0, -1.0); L1 = "Predicted (p_{T} > 110 GeV/#font[12]{c})"; P2 = std::make_pair(170.0, -1.0); L2 = "Predicted (p_{T} > 170 GeV/#font[12]{c})"; P3 = std::make_pair(230.0, -1.0); L3 = "Predicted (p_{T} > 230 GeV/#font[12]{c})"; Data1 = "Observed (p_{T} > 110 GeV/#font[12]{c})"; Data2 = "Observed (p_{T} > 170 GeV/#font[12]{c})"; Data3 = "Observed (p_{T} > 230 GeV/#font[12]{c})"; }else if(TypeMode==4 && S==0){ P1 = std::make_pair( -1.0 , 1.075); L1 = "Predicted (1/#beta > 1.075)"; P2 = std::make_pair( -1.0 , 1.100); L2 = "Predicted (1/#beta > 1.100)"; P3 = std::make_pair( -1.0 , 1.125); L3 = "Predicted (1/#beta > 1.125)"; Data1 = "Observed (1/#beta > 1.075)"; Data2 = "Observed (1/#beta > 1.100)"; Data3 = "Observed (1/#beta > 1.125)"; }else if(TypeMode==4 && S==1){ P1 = std::make_pair( -1.0 , 0.925); L1 = "Predicted (1/#beta < 0.925)"; P2 = std::make_pair( -1.0 , 0.900); L2 = "Predicted (1/#beta < 0.900)"; P3 = std::make_pair( -1.0 , 0.875); L3 = "Predicted (1/#beta < 0.875)"; Data1 = "Observed (1/#beta < 0.925)"; Data2 = "Observed (1/#beta < 0.900)"; Data3 = "Observed (1/#beta < 0.875)"; }else if(TypeMode==5){ P1 = std::make_pair( 75.0, -1.0); L1 = "Predicted (p_{T} > 75 GeV/#font[12]{c})"; P2 = std::make_pair(100.0, -1.0); L2 = "Predicted (p_{T} > 100 GeV/#font[12]{c})"; P3 = std::make_pair(125.0, -1.0); L3 = "Predicted (p_{T} > 125 GeV/#font[12]{c})"; Data1 = "Observed (p_{T} > 75 GeV/#font[12]{c})"; Data2 = "Observed (p_{T} > 100 GeV/#font[12]{c})"; Data3 = "Observed (p_{T} > 125 GeV/#font[12]{c})"; } c1->cd(); if(r==1) t1->cd(); mapObs[P1]->SetMarkerColor(2); mapObs[P1]->SetMarkerStyle(20); mapObs[P1]->Draw("P"); mapPred[P1]->SetLineColor(2); mapRatio[P1]->SetLineColor(2); mapPred[P1]->SetLineWidth(2.0); mapRatio[P1]->SetLineWidth(2.0); mapPred[P1]->SetFillColor(2); mapRatio[P1]->SetFillColor(2); mapPred[P1]->SetFillStyle(3004); mapRatio[P1]->SetFillStyle(3004); mapPred[P1]->Draw("3L"); //I (Chris) think the plot looks better with only two lines as it keeps the lines from bleeding into one another so removing middle line. /* mapObs[P2]->SetMarkerColor(4); mapObs[P2]->SetMarkerStyle(20); mapObs[P2]->Draw("P"); mapPred[P2]->SetLineColor(4); mapRatio[P2]->SetLineColor(4); mapPred[P2]->SetLineWidth(2.0); mapRatio[P2]->SetLineWidth(2.0); mapPred[P2]->SetFillColor(4); mapRatio[P2]->SetFillColor(4); mapPred[P2]->SetFillStyle(3005); mapRatio[P2]->SetFillStyle(3005); mapPred[P2]->Draw("3L"); */ mapObs[P3]->SetMarkerColor(8); mapObs[P3]->SetMarkerStyle(20); mapObs[P3]->Draw("P"); mapPred[P3]->SetLineColor(8); mapRatio[P3]->SetLineColor(8); mapPred[P3]->SetLineWidth(2.0); mapRatio[P3]->SetLineWidth(2.0); mapPred[P3]->SetFillColor(8); mapRatio[P3]->SetFillColor(8); mapPred[P3]->SetFillStyle(3007); mapRatio[P3]->SetFillStyle(3007); mapPred[P3]->Draw("3L"); //TH1D* obsLeg = (TH1D*)mapObs [P1]->Clone("ObsLeg"); //obsLeg->SetMarkerColor(1); LEG->AddEntry(mapObs[P1], Data1.c_str(),"P"); LEG->AddEntry(mapPred[P1], L1.c_str(),"FL"); //LEG->AddEntry(mapPred[P2], L2.c_str(),"FL"); //LEG->AddEntry(mapObs[P2], Data2.c_str(),"P"); LEG->AddEntry(mapObs[P3], Data3.c_str(),"P"); LEG->AddEntry(mapPred[P3], L3.c_str(),"FL"); LEG->Draw("same"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS), false); if(r==0) { c1->SetLogy(1); SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_NPredVsNObs"+suffix); } else { c1->cd(); TPad* t2; t2 = new TPad("t2","t2", 0.0, 0.0, 1.0, 0.2); t2->Draw(); t2->cd(); t2->SetGridy(true); t2->SetPad(0,0.0,1.0,0.2); t2->SetTopMargin(0); t2->SetBottomMargin(0.5); TH1D* frameR = new TH1D("frameR", "frameR", 1,std::max(xmin,0.02),xmax); frameR->GetXaxis()->SetNdivisions(505); frameR->SetTitle(""); frameR->SetStats(kFALSE); frameR->GetXaxis()->SetTitle("");//dEdxS_Legend.c_str()); frameR->GetYaxis()->SetTitle("Obs / Pred"); frameR->GetYaxis()->SetLabelFont(43); //give the font size in pixel (instead of fraction) frameR->GetYaxis()->SetLabelSize(15); //font size frameR->GetYaxis()->SetTitleFont(43); //give the font size in pixel (instead of fraction) frameR->GetYaxis()->SetTitleSize(15); //font size frameR->GetYaxis()->SetNdivisions(505); frameR->GetXaxis()->SetLabelFont(43); //give the font size in pixel (instead of fraction) frameR->GetXaxis()->SetLabelSize(15); //font size frameR->GetXaxis()->SetTitleFont(43); //give the font size in pixel (instead of fraction) frameR->GetXaxis()->SetTitleSize(15); //font size frameR->GetXaxis()->SetTitleOffset(3.75); frameR->SetMaximum(1.40); frameR->SetMinimum(0.60); frameR->Draw("AXIS"); t2->cd(); mapRatio[P1]->Draw("3C"); //mapRatio[P2]->Draw("3C"); mapRatio[P3]->Draw("3C"); TLine* LineAtOne = new TLine(xmin,1,xmax,1); LineAtOne->SetLineStyle(3); LineAtOne->Draw(); c1->cd(); SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_NPredVsNObsWRatio"+suffix); delete t2; } delete c1; } } } //README: Draw a map of the prediction/observation for the various selection //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_D = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_D"); //TH1D* H_P = (TH1D*)GetObjectFromPath(InputFile, Data+"/H_P"); //c1 = new TCanvas("c1","c1,",600,600); legend.clear(); //TH2D* DataVsPred = new TH2D("DataVsPred","DataVsPred",17,30,200, 8,0.05,0.5); //TH2D* DataMap = new TH2D("DataMap" ,"DataMap" ,17,30,200, 8,0.05,0.5); //TH2D* PredMap = new TH2D("PredMap" ,"PredMap" ,17,30,200, 8,0.05,0.5); //for(unsigned int CutIndex_=0;CutIndex_<H_P->GetNbinsX();CutIndex_++){ // double P = H_P->GetBinContent(CutIndex_+1); // double D = H_D->GetBinContent(CutIndex_+1); // double Err = sqrt( pow(H_P->GetBinError(CutIndex_+1),2) + std::max(D,1.0) ); // double NSigma = (D-P)/Err; // DataMap->SetBinContent(DataVsPred->GetXaxis()->FindBin(HCuts_Pt->GetBinContent(CutIndex_+1)), DataVsPred->GetYaxis()->FindBin(HCuts_I->GetBinContent(CutIndex_+1)), D); // PredMap->SetBinContent(DataVsPred->GetXaxis()->FindBin(HCuts_Pt->GetBinContent(CutIndex_+1)), DataVsPred->GetYaxis()->FindBin(HCuts_I->GetBinContent(CutIndex_+1)), P); // if(isnan(P) || P<=0)continue; //Is <=0 only when prediction failed or is not meaningful (i.e. WP=(0,0,0) ) // DataVsPred->SetBinContent(DataVsPred->GetXaxis()->FindBin(HCuts_Pt->GetBinContent(CutIndex_+1)), DataVsPred->GetYaxis()->FindBin(HCuts_I->GetBinContent(CutIndex_+1)), NSigma); //} //DataVsPred->SetMinimum(-3); //DataVsPred->SetMaximum(3); //Histos[0] = DataVsPred; legend.push_back("Observed"); //DrawSuperposedHistos((TH1**)Histos, legend, "COLZ", "PtCut", "ICut", 0,0, 0,0); //DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); //SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_PredVsObs"); //delete c1; //PredMap->SetMinimum(1E-2); //DataMap->SetMinimum(1E-2); //PredMap->SetMaximum(std::max(PredMap->GetMaximum(),DataMap->GetMaximum())); //DataMap->SetMaximum(std::max(PredMap->GetMaximum(),DataMap->GetMaximum())); //c1 = new TCanvas("c1","c1,",600,600); legend.clear(); //c1->SetLogz(true); //Histos[0] = PredMap; legend.push_back("Observed"); //DrawSuperposedHistos((TH1**)Histos, legend, "COLZ", "PtCut", "ICut", 0,0, 0,0); //DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); //SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_Pred"); //delete c1; //c1 = new TCanvas("c1","c1,",600,600); legend.clear(); //c1->SetLogz(true); //Histos[0] = DataMap; legend.push_back("Observed"); //DrawSuperposedHistos((TH1**)Histos, legend, "COLZ", "PtCut", "ICut", 0,0, 0,0); //DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); //SaveCanvas(c1,InputPattern,string("Prediction_")+Data+"_Data"); //delete c1; InputFile->Close(); } // print the event flow table for all signal point and Data and MCTr void CutFlow(string InputPattern, unsigned int CutIndex){ TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); if(!InputFile)std::cout << "FileProblem\n"; TH1D* HCuts_Pt = (TH1D*)GetObjectFromPath(InputFile, "HCuts_Pt"); TH1D* HCuts_I = (TH1D*)GetObjectFromPath(InputFile, "HCuts_I"); TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, "HCuts_TOF"); char Buffer[1024]; sprintf(Buffer,"%s/CutFlow_%03i_Pt%03.0f_I%05.3f_TOF%04.3f.txt",InputPattern.c_str(),CutIndex,HCuts_Pt->GetBinContent(CutIndex+1),HCuts_I->GetBinContent(CutIndex+1),HCuts_TOF->GetBinContent(CutIndex+1)); FILE* pFile = fopen(Buffer,"w"); stPlots plots; if(stPlots_InitFromFile(InputFile, plots,"Data8TeV")){ stPlots_Dump(plots, pFile, CutIndex); stPlots_Clear(&plots);} if(stPlots_InitFromFile(InputFile, plots,"Data7TeV")){ stPlots_Dump(plots, pFile, CutIndex); stPlots_Clear(&plots); } if(stPlots_InitFromFile(InputFile, plots,"MCTr_7TeV")){ stPlots_Dump(plots, pFile, CutIndex); stPlots_Clear(&plots); } if(stPlots_InitFromFile(InputFile, plots,"MCTr_8TeV")){ stPlots_Dump(plots, pFile, CutIndex); stPlots_Clear(&plots); } if(stPlots_InitFromFile(InputFile, plots,"Cosmic8TeV")){ stPlots_Dump(plots, pFile, CutIndex); stPlots_Clear(&plots); } for(unsigned int s=0;s<samples.size();s++){ if(samples[s].Type!=2 || !samples[s].MakePlot)continue; if(stPlots_InitFromFile(InputFile, plots, samples[s].Name)){ stPlots_Dump(plots, pFile, CutIndex); stPlots_Clear(&plots); } } fclose(pFile); InputFile->Close(); } // make all plots of the preselection and selection variables as well as some plots showing 2D planes void SelectionPlot(string InputPattern, unsigned int CutIndex, unsigned int CutIndexTight){ string LegendTitle = LegendFromType(InputPattern);; TypeMode = TypeFromPattern(InputPattern); TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); stPlots Data8TeVPlots, Data7TeVPlots, MCTr8TeVPlots, MCTr7TeVPlots, Cosmic7TeVPlots, Cosmic8TeVPlots, SignPlots[samples.size()]; TypeMode = TypeFromPattern(InputPattern); stPlots_InitFromFile(InputFile, Data8TeVPlots,"Data8TeV"); stPlots_InitFromFile(InputFile, Data7TeVPlots,"Data7TeV"); stPlots_InitFromFile(InputFile, MCTr8TeVPlots ,"MCTr_8TeV"); stPlots_InitFromFile(InputFile, MCTr7TeVPlots ,"MCTr_7TeV"); if(TypeMode==3) { stPlots_InitFromFile(InputFile, Cosmic8TeVPlots,"Cosmic8TeV"); //stPlots_InitFromFile(InputFile, Cosmic7TeVPlots,"Cosmic8TeV"); } for(unsigned int s=0;s<samples.size();s++){ if (samples[s].Name!="Gluino_7TeV_M300_f10" && samples[s].Name!="Gluino_7TeV_M600_f10" && samples[s].Name!="Gluino_7TeV_M800_f10" && samples[s].Name!="Gluino_8TeV_M1200_f100" && samples[s].Name!="Gluino_8TeV_M300_f10" && samples[s].Name!="Gluino_8TeV_M600_f10" && samples[s].Name!="Gluino_8TeV_M800_f10" && samples[s].Name!="GMStau_7TeV_M247" && samples[s].Name!="GMStau_7TeV_M370" && samples[s].Name!="GMStau_7TeV_M494" && samples[s].Name!="GMStau_8TeV_M247" && samples[s].Name!="GMStau_8TeV_M370" && samples[s].Name!="GMStau_8TeV_M494" && samples[s].Name!="DY_7TeV_M100_Q1o3" && samples[s].Name!="DY_7TeV_M600_Q1o3" && samples[s].Name!="DY_7TeV_M100_Q2o3" && samples[s].Name!="DY_7TeV_M600_Q2o3" && samples[s].Name!="DY_8TeV_M100_Q1o3" && samples[s].Name!="DY_8TeV_M600_Q1o3" && samples[s].Name!="DY_8TeV_M100_Q2o3" && samples[s].Name!="DY_8TeV_M600_Q2o3" && samples[s].Name!="DY_8TeV_M400_Q1" && samples[s].Name!="DY_8TeV_M400_Q3" && samples[s].Name!="DY_8TeV_M400_Q5" && samples[s].Name!="DY_7TeV_M400_Q1" && samples[s].Name!="DY_7TeV_M400_Q3" && samples[s].Name!="DY_7TeV_M400_Q5" && samples[s].Name!="Gluino_8TeV_M500_f100" && samples[s].Name!="Stop_8TeV_M500") continue; if(!stPlots_InitFromFile(InputFile, SignPlots[s],samples[s].Name)){printf("Missing sample %s\n",samples[s].Name.c_str());continue;} if(samples[s].Name=="Gluino_8TeV_M600_f10" || samples[s].Name=="DY_8TeV_M600_Q2o3")stPlots_Draw(SignPlots[s], InputPattern + "/Selection_" + samples[s].Name, LegendTitle, CutIndex); } SQRTS=8; stPlots_Draw(Data8TeVPlots, InputPattern + "/Selection_Data8TeV", LegendTitle, CutIndex); if(TypeMode!=3) {SQRTS=7; stPlots_Draw(Data7TeVPlots, InputPattern + "/Selection_Data7TeV", LegendTitle, CutIndex);} SQRTS=8; stPlots_Draw(MCTr8TeVPlots , InputPattern + "/Selection_MCTr_8TeV" , LegendTitle, CutIndex); if(TypeMode!=3) {SQRTS=7; stPlots_Draw(MCTr7TeVPlots , InputPattern + "/Selection_MCTr_7TeV" , LegendTitle, CutIndex);} if(TypeMode==3) { stPlots_Draw(Cosmic8TeVPlots, InputPattern + "/Selection_Cosmic8TeV", LegendTitle, CutIndex); //stPlots_Draw(Cosmic11Plots, InputPattern + "/Selection_Cosmic11", LegendTitle, CutIndex); } if(TypeMode!=3) stPlots_DrawComparison(InputPattern + "/Selection_Comp_Data" , LegendTitle, CutIndex, CutIndexTight, &Data8TeVPlots, &Data7TeVPlots, &MCTr8TeVPlots, &MCTr7TeVPlots); if(TypeMode<=2){ SQRTS=7; stPlots_DrawComparison(InputPattern + "/Selection_Comp_7TeV_Gluino", LegendTitle, CutIndex, CutIndexTight, &Data7TeVPlots, &MCTr7TeVPlots, &SignPlots[JobIdToIndex("Gluino_7TeV_M300_f10",samples)], &SignPlots[JobIdToIndex("Gluino_7TeV_M600_f10",samples)], &SignPlots[JobIdToIndex("Gluino_7TeV_M800_f10",samples)]);} if(TypeMode<=2){ SQRTS=8; stPlots_DrawComparison(InputPattern + "/Selection_Comp_8TeV_Gluino", LegendTitle, CutIndex, CutIndexTight, &Data8TeVPlots, &MCTr8TeVPlots, &SignPlots[JobIdToIndex("Gluino_8TeV_M300_f10",samples)], &SignPlots[JobIdToIndex("Gluino_8TeV_M600_f10",samples)], &SignPlots[JobIdToIndex("Gluino_8TeV_M800_f10",samples)]);} if(TypeMode==3){ //SQRTS=7; stPlots_DrawComparison(InputPattern + "/Selection_Comp_Cosmic_7TeV", LegendTitle, CutIndex, CutIndexTight, &Data7TeVPlots,&SignPlots[JobIdToIndex("Gluino_7TeV_M800_f10",samples)]); SQRTS=8; stPlots_DrawComparison(InputPattern + "/Selection_Comp_8TeV_Cosmic", LegendTitle, CutIndex, CutIndexTight, &Data8TeVPlots, &MCTr8TeVPlots, &Cosmic8TeVPlots, &SignPlots[JobIdToIndex("Gluino_8TeV_M1200_f100",samples)]); SQRTS=8; stPlots_DrawComparison(InputPattern + "/Selection_Comp_Signal_8TeV", LegendTitle, CutIndex, CutIndexTight, &Data8TeVPlots, &SignPlots[JobIdToIndex("Gluino_8TeV_M500_f100",samples)], &SignPlots[JobIdToIndex("Stop_8TeV_M500",samples)], &SignPlots[JobIdToIndex("GMStau_8TeV_M494",samples)]); //SQRTS=78; stPlots_DrawComparison(InputPattern + "/Selection_Comp_Cosmic_78TeV", LegendTitle, CutIndex, CutIndexTight, &Data8TeVPlots, &Data7TeVPlots, &Cosmic8TeVPlots, &SignPlots[JobIdToIndex("Gluino_7TeV_M800_f10",samples)], &SignPlots[JobIdToIndex("Gluino_8TeV_M800_f10",samples)]); } if(TypeMode==0 || TypeMode==4){ SQRTS=8; stPlots_DrawComparison(InputPattern + "/Selection_Comp_8TeV_DY_QG" , LegendTitle, CutIndex, CutIndexTight, &Data8TeVPlots, &MCTr8TeVPlots, &SignPlots[JobIdToIndex("DY_8TeV_M400_Q1",samples)], &SignPlots[JobIdToIndex("DY_8TeV_M400_Q3",samples)], &SignPlots[JobIdToIndex("DY_8TeV_M400_Q5",samples)]); SQRTS=7; stPlots_DrawComparison(InputPattern + "/Selection_Comp_7TeV_DY_QG" , LegendTitle, CutIndex, CutIndexTight, &Data7TeVPlots, &MCTr7TeVPlots, &SignPlots[JobIdToIndex("DY_7TeV_M400_Q1",samples)], &SignPlots[JobIdToIndex("DY_7TeV_M400_Q3",samples)], &SignPlots[JobIdToIndex("DY_7TeV_M400_Q5",samples)]); } if(TypeMode<3){ SQRTS=7; stPlots_DrawComparison(InputPattern + "/Selection_Comp_7TeV_DY" , LegendTitle, CutIndex, CutIndexTight, &Data7TeVPlots, &MCTr7TeVPlots, &SignPlots[JobIdToIndex("DY_7TeV_M100_Q2o3",samples)], &SignPlots[JobIdToIndex("DY_7TeV_M600_Q2o3",samples)]); SQRTS=8; stPlots_DrawComparison(InputPattern + "/Selection_Comp_8TeV_DY" , LegendTitle, CutIndex, CutIndexTight, &Data8TeVPlots, &MCTr8TeVPlots, &SignPlots[JobIdToIndex("DY_8TeV_M100_Q2o3",samples)], &SignPlots[JobIdToIndex("DY_8TeV_M600_Q2o3",samples)]); } if(TypeMode==5 || TypeMode==3){ if(TypeMode==5) {SQRTS=7; stPlots_DrawComparison(InputPattern + "/Selection_Comp_7TeV_DY" , LegendTitle, CutIndex, CutIndexTight, &Data7TeVPlots, &MCTr7TeVPlots, &SignPlots[JobIdToIndex("DY_7TeV_M100_Q1o3",samples)], &SignPlots[JobIdToIndex("DY_7TeV_M100_Q2o3",samples)], &SignPlots[JobIdToIndex("DY_7TeV_M600_Q2o3",samples)]);} SQRTS=8; stPlots_DrawComparison(InputPattern + "/Selection_Comp_8TeV_DY" , LegendTitle, CutIndex, CutIndexTight, &Data8TeVPlots, &MCTr8TeVPlots, &SignPlots[JobIdToIndex("DY_8TeV_M100_Q1o3",samples)], &SignPlots[JobIdToIndex("DY_8TeV_M100_Q2o3",samples)], &SignPlots[JobIdToIndex("DY_8TeV_M600_Q2o3",samples)]); } stPlots_Clear(&Data8TeVPlots); if(TypeMode!=3) stPlots_Clear(&Data7TeVPlots); stPlots_Clear(&MCTr8TeVPlots); if(TypeMode!=3) stPlots_Clear(&MCTr7TeVPlots); for(unsigned int s=0;s<samples.size();s++){ if (samples[s].Name!="Gluino_7TeV_M300_f10" && samples[s].Name!="Gluino_7TeV_M600_f10" && samples[s].Name!="Gluino_7TeV_M800_f10" && "Gluino_8TeV_M300_f10" && samples[s].Name!="Gluino_8TeV_M600_f10" && samples[s].Name!="Gluino_8TeV_M800_f10" && samples[s].Name!="GMStau_7TeV_M247" && samples[s].Name!="GMStau_7TeV_M370" && samples[s].Name!="GMStau_7TeV_M494" && samples[s].Name!="GMStau_8TeV_M247" && samples[s].Name!="GMStau_8TeV_M370" && samples[s].Name!="GMStau_8TeV_M494" && samples[s].Name!="DY_7TeV_M100_Q1o3" && samples[s].Name!="DY_7TeV_M600_Q1o3" && samples[s].Name!="DY_7TeV_M100_Q2o3" && samples[s].Name!="DY_7TeV_M600_Q2o3" && samples[s].Name!="DY_8TeV_M100_Q1o3" && samples[s].Name!="DY_8TeV_M600_Q1o3" && samples[s].Name!="DY_8TeV_M100_Q2o3" && samples[s].Name!="DY_8TeV_M600_Q2o3" && samples[s].Name!="DY_8TeV_M400_Q1" && samples[s].Name!="DY_8TeV_M400_Q3" && samples[s].Name!="DY_8TeV_M400_Q5") continue; if(!stPlots_InitFromFile(InputFile, SignPlots[s],samples[s].Name)) continue; stPlots_Clear(&SignPlots[s]); } InputFile->Close(); } // Determine the systematic uncertainty by computing datadriven prediction from different paths (only works with 3D ABCD method) void GetSystematicOnPrediction(string InputPattern, string DataName){ if(DataName.find("7TeV")!=string::npos){SQRTS=7.0;}else{SQRTS=8.0;} TypeMode = TypeFromPattern(InputPattern); if(TypeMode!=2)return; string LegendTitle = LegendFromType(InputPattern);; 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_A = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_A"); TH1D* H_B = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_B"); TH1D* H_C = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_C"); //TH1D* H_D = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_D"); TH1D* H_E = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_E"); TH1D* H_F = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_F"); TH1D* H_G = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_G"); TH1D* H_H = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_H"); //TH1D* H_P = (TH1D*)GetObjectFromPath(InputFile, DataName+"/H_P"); int ArrN[8] = {0,0,0,0,0,0,0,0}; double ArrPred[5][8][20]; double ArrErr[5][8][20]; int ArrPredN[5][8]; for(unsigned int i=0;i<6;i++){for(unsigned int j=0;j<8;j++){ArrPredN[i][j]=0;}} //double ArrMean [8][20]; double ArrSigma[8][20]; double ArrDist [8][20]; double ArrSum [8][20]; double ArrSyst [8][20]; double ArrStat [8][20]; double ArrStatB[8][20]; double ArrPt [8][20]; double ArrI [8][20]; double ArrT [8][20]; std::vector<int> Index; std::vector<int> Plot; //variation on TOF cut 50, 0.05 1.05->1.2 Index.push_back(16); Plot.push_back(0); Index.push_back(17); Plot.push_back(0); Index.push_back(18); Plot.push_back(0); Index.push_back(19); Plot.push_back(0); Index.push_back(20); Plot.push_back(0); Index.push_back(21); Plot.push_back(0); Index.push_back(22); Plot.push_back(0); //variation on I cut 50, 0.05->0.225 1.05 Index.push_back(16); Plot.push_back(1); Index.push_back(30); Plot.push_back(1); Index.push_back(44); Plot.push_back(1); Index.push_back(58); Plot.push_back(1); Index.push_back(72); Plot.push_back(1); Index.push_back(86); Plot.push_back(1); Index.push_back(100); Plot.push_back(1); Index.push_back(114); Plot.push_back(1); //variation on Pt cut 50->90 0.30 1.05 Index.push_back(156); Plot.push_back(2); Index.push_back(366); Plot.push_back(2); Index.push_back(576); Plot.push_back(2); Index.push_back(786); Plot.push_back(2); Index.push_back( 996); Plot.push_back(2); Index.push_back(1206); Plot.push_back(2); Index.push_back(1416); Plot.push_back(2); Index.push_back(1626); Plot.push_back(2); //variation on Pt cut 50->90 0.1 1.1 Index.push_back(46); Plot.push_back(3); Index.push_back(256); Plot.push_back(3); Index.push_back(466); Plot.push_back(3); Index.push_back(676); Plot.push_back(3); Index.push_back(886); Plot.push_back(3); Index.push_back(1096); Plot.push_back(3); Index.push_back(1306); Plot.push_back(3); Index.push_back(1516); Plot.push_back(3); //variation on Pt cut 50->90 0.15 1.05 // Index.push_back(72); Plot.push_back(4); // Index.push_back(492); Plot.push_back(4); // Index.push_back(912); Plot.push_back(4); // Index.push_back(1332); Plot.push_back(4); // Index.push_back(1752); Plot.push_back(4); // Index.push_back(2172); Plot.push_back(4); // Index.push_back(2592); Plot.push_back(4); // Index.push_back(2802); Plot.push_back(4); //variation on Pt cut 50->90 0.10 1.20 Index.push_back(50); Plot.push_back(4); Index.push_back(260); Plot.push_back(4); Index.push_back(470); Plot.push_back(4); Index.push_back(680); Plot.push_back(4); Index.push_back(890); Plot.push_back(4); Index.push_back(1100); Plot.push_back(4); Index.push_back(1310); Plot.push_back(4); Index.push_back(1520); Plot.push_back(4); //Not used Index.push_back(82 + 4); Plot.push_back(5); Index.push_back(154+ 4); Plot.push_back(5); Index.push_back(226+ 4); Plot.push_back(5); Index.push_back(298+ 4); Plot.push_back(5); Index.push_back(370+ 4); Plot.push_back(5); Index.push_back(442+ 4); Plot.push_back(5); Index.push_back(514+ 4); Plot.push_back(5); Index.push_back(586+ 4); Plot.push_back(5); Index.push_back(658+ 4); Plot.push_back(5); Index.push_back(730+ 4); Plot.push_back(5); Index.push_back(802+ 4); Plot.push_back(5); //variation on Pt cut 50->115 0.05 1.05 Index.push_back(16); Plot.push_back(6); Index.push_back(436); Plot.push_back(6); Index.push_back(856); Plot.push_back(6); Index.push_back(1276); Plot.push_back(6); Index.push_back(1696); Plot.push_back(6); Index.push_back(2116); Plot.push_back(6); Index.push_back(2536); Plot.push_back(6); Index.push_back(2746); Plot.push_back(6); for(unsigned int i=0;i<Index.size();i++){ int CutIndex = Index[i]; const double& A=H_A->GetBinContent(CutIndex+1); const double& B=H_B->GetBinContent(CutIndex+1); const double& C=H_C->GetBinContent(CutIndex+1); //const double& D=H_D->GetBinContent(CutIndex+1); const double& E=H_E->GetBinContent(CutIndex+1); const double& F=H_F->GetBinContent(CutIndex+1); const double& G=H_G->GetBinContent(CutIndex+1); const double& H=H_H->GetBinContent(CutIndex+1); double Pred[8]; double Err [8]; double N = 0; double Sigma = 0; double Mean = 0; for(unsigned int p=0;p<7;p++){ Pred[p] = -1; Err [p] = -1; if(p==0){ if(A<25 || F<25 || G<25 || E<25)continue; Pred[p] = (A*F*G)/(E*E); Err [p] = Pred [p] * sqrt( 1/A + 1/F + 1/G + 4/E); }else if(p==1){ if(A<25 || H<25 || E<25)continue; Pred[p] = ((A*H)/E); Err [p] = Pred[p] * sqrt( 1/A+1/H+1/E ); }else if (p==2){ if(B<25 || G<25 || E<25)continue; Pred[p] = ((B*G)/E); Err [p] = Pred[p] * sqrt( 1/B+ 1/G+ 1/E ); }else if (p==3){ if(F<25 || C<25 || E<25)continue; Pred[p] = ((F*C)/E); Err [p] = Pred[p] * sqrt( 1/F + 1/C + 1/E ); } if(Pred[p]>=0){ N++; Mean += Pred[p]/pow(Err [p],2); Sigma += 1 /pow(Err [p],2); } ArrPred [p][Plot[i]][ArrN[Plot[i]]] = Pred[p]; ArrErr [p][Plot[i]][ArrN[Plot[i]]] = Err [p]; if(Pred[p]>=0)ArrPredN[p][Plot[i]]++; } Mean = Mean/Sigma; Sigma = sqrt(Sigma); double Dist = fabs(Pred[0] - Mean); double Sum=0, Stat=0, Syst=0, StatB=0; for(unsigned int p=0;p<7;p++){ if(Pred[p]>=0){ Sum += pow(Pred[p]-Mean,2); Stat += pow(Err [p],2); StatB += Err [p]; } } Sum = sqrt(Sum/(N-1)); Stat = sqrt(Stat)/N; StatB= StatB/N; Syst = sqrt(Sum*Sum - Stat*Stat); //printf("pT>%6.2f I> %6.2f TOF>%6.2f : ", HCuts_Pt ->GetBinContent(CutIndex+1), HCuts_I ->GetBinContent(CutIndex+1), HCuts_TOF->GetBinContent(CutIndex+1)); //printf("A =%6.2E, B=%6.2E, C=%6.2E, D=%6.2E E =%6.2E, F=%6.2E, G=%6.2E, H=%6.2E\n", A,B,C,D, E, F, G, H); printf("--> N = %1.0f Mean = %8.2E Sigma=%8.2E Dist=%8.2E Sum=%8.2E Stat=%8.2E Syst=%8.2E\n", N, Mean, Sigma/Mean, Dist/Mean, Sum/Mean, Stat/Mean, Syst/Mean); if(N>0){ //ArrMean [Plot[i]][ArrN[Plot[i]]] = Mean; ArrSigma [Plot[i]][ArrN[Plot[i]]] = Sigma/Mean; ArrDist [Plot[i]][ArrN[Plot[i]]] = Dist/Mean; ArrSum [Plot[i]][ArrN[Plot[i]]] = Sum/Mean; ArrSyst [Plot[i]][ArrN[Plot[i]]] = Syst/Mean; ArrStat [Plot[i]][ArrN[Plot[i]]] = Stat/Mean; ArrStatB [Plot[i]][ArrN[Plot[i]]] = StatB/Mean; ArrPt [Plot[i]][ArrN[Plot[i]]] = HCuts_Pt ->GetBinContent(CutIndex+1); ; ArrI [Plot[i]][ArrN[Plot[i]]] = HCuts_I ->GetBinContent(CutIndex+1); ; ArrT [Plot[i]][ArrN[Plot[i]]] = HCuts_TOF->GetBinContent(CutIndex+1); ; ArrN[Plot[i]]++; } } TGraphErrors* graph_T0 = new TGraphErrors(ArrPredN[0][0],ArrT [0],ArrPred[0][0],0,ArrErr[0][0]); graph_T0->SetLineColor(1); graph_T0->SetMarkerColor(1); graph_T0->SetMarkerStyle(20); TGraphErrors* graph_T1 = new TGraphErrors(ArrPredN[1][0],ArrT [0],ArrPred[1][0],0,ArrErr[1][0]); graph_T1->SetLineColor(2); graph_T1->SetMarkerColor(2); graph_T1->SetMarkerStyle(21); TGraphErrors* graph_T2 = new TGraphErrors(ArrPredN[2][0],ArrT [0],ArrPred[2][0],0,ArrErr[2][0]); graph_T2->SetLineColor(4); graph_T2->SetMarkerColor(4); graph_T2->SetMarkerStyle(22); TGraphErrors* graph_T3 = new TGraphErrors(ArrPredN[3][0],ArrT [0],ArrPred[3][0],0,ArrErr[3][0]); graph_T3->SetLineColor(8); graph_T3->SetMarkerColor(8); graph_T3->SetMarkerStyle(23); TGraphErrors* graph_I0 = new TGraphErrors(ArrPredN[0][1],ArrI [1],ArrPred[0][1],0,ArrErr[0][1]); graph_I0->SetLineColor(1); graph_I0->SetMarkerColor(1); graph_I0->SetMarkerStyle(20); TGraphErrors* graph_I1 = new TGraphErrors(ArrPredN[1][1],ArrI [1],ArrPred[1][1],0,ArrErr[1][1]); graph_I1->SetLineColor(2); graph_I1->SetMarkerColor(2); graph_I1->SetMarkerStyle(21); TGraphErrors* graph_I2 = new TGraphErrors(ArrPredN[2][1],ArrI [1],ArrPred[2][1],0,ArrErr[2][1]); graph_I2->SetLineColor(4); graph_I2->SetMarkerColor(4); graph_I2->SetMarkerStyle(22); TGraphErrors* graph_I3 = new TGraphErrors(ArrPredN[3][1],ArrI [1],ArrPred[3][1],0,ArrErr[3][1]); graph_I3->SetLineColor(8); graph_I3->SetMarkerColor(8); graph_I3->SetMarkerStyle(23); TGraphErrors* graph_P0 = new TGraphErrors(ArrPredN[0][6],ArrPt[6],ArrPred[0][6],0,ArrErr[0][6]); graph_P0->SetLineColor(1); graph_P0->SetMarkerColor(1); graph_P0->SetMarkerStyle(20); TGraphErrors* graph_P1 = new TGraphErrors(ArrPredN[1][6],ArrPt[6],ArrPred[1][6],0,ArrErr[1][6]); graph_P1->SetLineColor(2); graph_P1->SetMarkerColor(2); graph_P1->SetMarkerStyle(21); TGraphErrors* graph_P2 = new TGraphErrors(ArrPredN[2][6],ArrPt[6],ArrPred[2][6],0,ArrErr[2][6]); graph_P2->SetLineColor(4); graph_P2->SetMarkerColor(4); graph_P2->SetMarkerStyle(22); TGraphErrors* graph_P3 = new TGraphErrors(ArrPredN[3][6],ArrPt[6],ArrPred[3][6],0,ArrErr[3][6]); graph_P3->SetLineColor(8); graph_P3->SetMarkerColor(8); graph_P3->SetMarkerStyle(23); TLegend* LEG = NULL; LEG = new TLegend(0.50,0.65,0.80,0.90); LEG->SetFillColor(0); LEG->SetBorderSize(0); LEG->AddEntry(graph_T0, "D=AFG/EE" ,"LP"); LEG->AddEntry(graph_T1, "D=AH/E" ,"LP"); LEG->AddEntry(graph_T2, "D=BG/E" ,"LP"); LEG->AddEntry(graph_T3, "D=FC/E" ,"LP"); TCanvas* c1; c1 = new TCanvas("c1", "c1",600,600); c1->SetLogy(true); TMultiGraph* MGTOF = new TMultiGraph(); MGTOF->Add(graph_T0 ,"LP"); MGTOF->Add(graph_T1 ,"LP"); MGTOF->Add(graph_T2 ,"LP"); MGTOF->Add(graph_T3 ,"LP"); MGTOF->Draw("A"); MGTOF->SetTitle(""); MGTOF->GetXaxis()->SetTitle("1/#beta selection"); MGTOF->GetYaxis()->SetTitle("Number of expected backgrounds"); MGTOF->GetYaxis()->SetTitleOffset(1.70); MGTOF->GetYaxis()->SetRangeUser(500,1E6); LEG->Draw(); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Systematics_")+DataName+"_TOF_Value"); delete c1; c1 = new TCanvas("c1", "c1",600,600); TMultiGraph* MGI = new TMultiGraph(); c1->SetLogy(true); MGI->Add(graph_I0 ,"LP"); MGI->Add(graph_I1 ,"LP"); MGI->Add(graph_I2 ,"LP"); MGI->Add(graph_I3 ,"LP"); MGI->Draw("A"); MGI->SetTitle(""); MGI->GetXaxis()->SetTitle("I_{as} selection"); MGI->GetYaxis()->SetTitle("Number of expected backgrounds"); MGI->GetYaxis()->SetTitleOffset(1.70); MGI->GetYaxis()->SetRangeUser(500,1E6); LEG->Draw(); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Systematics_")+DataName+"_I_Value"); delete c1; c1 = new TCanvas("c1", "c1",600,600); c1->SetLogy(true); TMultiGraph* MGP = new TMultiGraph(); MGP->Add(graph_P0 ,"LP"); MGP->Add(graph_P1 ,"LP"); MGP->Add(graph_P2 ,"LP"); MGP->Add(graph_P3 ,"LP"); MGP->Draw("A"); MGP->SetTitle(""); MGP->GetXaxis()->SetTitle("p_{T} selection"); MGP->GetYaxis()->SetTitle("Number of expected backgrounds"); MGP->GetYaxis()->SetTitleOffset(1.70); MGP->GetYaxis()->SetRangeUser(500,1E6); LEG->Draw(); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,InputPattern,string("Systematics_")+DataName+"_P_Value"); delete c1; for(unsigned int p=0;p<3;p++){ string Title; string Name; if(p==0){ Title = "1/#beta selection"; Name="TOF_"; } if(p==1){ Title = "dEdx selection"; Name="I_"; } if(p==2){ Title = "p_{T} selection"; Name="pT_"; } c1 = new TCanvas("c1","c1", 600, 600); TGraph* graph_s; if(p==0)graph_s = new TGraph(ArrN[p],ArrT [p],ArrSigma[p]); if(p==1)graph_s = new TGraph(ArrN[p],ArrI [p],ArrSigma[p]); if(p==2)graph_s = new TGraph(ArrN[p],ArrPt[p],ArrSigma[p]); graph_s->SetTitle(""); graph_s->GetYaxis()->SetTitle("Prediction #sigma/#mu"); graph_s->GetYaxis()->SetTitleOffset(1.70); graph_s->GetXaxis()->SetTitle(Title.c_str()); graph_s->Draw("AC*"); SaveCanvas(c1,InputPattern,string(string("Systematics_")+DataName+"_")+Name+"Sigma"); delete c1; c1 = new TCanvas("c1","c1", 600, 600); TGraph* graph_d; if(p==0)graph_d = new TGraph(ArrN[p],ArrT [p],ArrDist[p]); if(p==1)graph_d = new TGraph(ArrN[p],ArrI [p],ArrDist[p]); if(p==2)graph_d = new TGraph(ArrN[p],ArrPt[p],ArrDist[p]); graph_d->SetTitle(""); graph_d->GetYaxis()->SetTitle("Prediction Dist/#mu"); graph_d->GetYaxis()->SetTitleOffset(1.70); graph_d->GetXaxis()->SetTitle(Title.c_str()); graph_d->Draw("AC*"); SaveCanvas(c1,InputPattern,string(string("Systematics_")+DataName+"_")+Name+"Dist"); delete c1; c1 = new TCanvas("c1","c1", 600, 600); TGraph* graph_sum; if(p==0)graph_sum = new TGraph(ArrN[p],ArrT [p],ArrSum[p]); if(p==1)graph_sum = new TGraph(ArrN[p],ArrI [p],ArrSum[p]); if(p==2)graph_sum = new TGraph(ArrN[p+2],ArrPt[p+2],ArrSum[p+2]); graph_sum->SetTitle(""); graph_sum->GetYaxis()->SetTitle("Prediction #sigma_{Stat+Syst}/#mu"); graph_sum->GetYaxis()->SetTitleOffset(1.70); graph_sum->GetXaxis()->SetTitle(Title.c_str()); graph_sum->Draw("AC*"); graph_sum->GetYaxis()->SetRangeUser(0,0.5); if(p==2){ TGraph* graph_sum2 = new TGraph(ArrN[p+1],ArrPt[p+1],ArrSum[p+1]); graph_sum2->SetLineColor(2); graph_sum2->SetMarkerColor(2); graph_sum2->Draw("C*"); TGraph* graph_sum3 = new TGraph(ArrN[p+0],ArrPt[p+0],ArrSum[p+0]); graph_sum3->SetLineColor(4); graph_sum3->SetMarkerColor(4); graph_sum3->Draw("C*"); //TGraph* graph_sum4 = new TGraph(ArrN[p+3],ArrPt[p+3],ArrSum[p+3]); //graph_sum4->SetLineColor(8); //graph_sum4->SetMarkerColor(8); //graph_sum4->Draw("C*"); LEG = new TLegend(0.50,0.65,0.80,0.90); LEG->SetFillColor(0); LEG->SetBorderSize(0); LEG->AddEntry(graph_sum2, "I_{as}>0.10 & 1/#beta>1.10", "L"); LEG->AddEntry(graph_sum3, "I_{as}>0.30 & 1/#beta>1.05", "L"); LEG->AddEntry(graph_sum, "I_{as}>0.10 & 1/#beta>1.20", "L"); LEG->Draw(); } SaveCanvas(c1,InputPattern,string(string("Systematics_")+DataName+"_")+Name+"Sum"); delete c1; c1 = new TCanvas("c1","c1", 600, 600); TGraph* graph_stat; if(p==0)graph_stat = new TGraph(ArrN[p],ArrT [p],ArrStat[p]); if(p==1)graph_stat = new TGraph(ArrN[p],ArrI [p],ArrStat[p]); if(p==2)graph_stat = new TGraph(ArrN[p+2],ArrPt[p+2],ArrStat[p+2]); graph_stat->SetTitle(""); graph_stat->GetYaxis()->SetTitle("Prediction #sigma_{Stat}/#mu"); graph_stat->GetYaxis()->SetTitleOffset(1.70); graph_stat->GetXaxis()->SetTitle(Title.c_str()); graph_stat->Draw("AC*"); graph_stat->GetYaxis()->SetRangeUser(0,0.25); if(p==2){ TGraph* graph_stat2 = new TGraph(ArrN[p+1],ArrPt[p+1],ArrStat[p+1]); graph_stat2->SetLineColor(2); graph_stat2->SetMarkerColor(2); graph_stat2->Draw("C*"); TGraph* graph_stat3 = new TGraph(ArrN[p+0],ArrPt[p+0],ArrStat[p+0]); graph_stat3->SetLineColor(4); graph_stat3->SetMarkerColor(4); graph_stat3->Draw("C*"); //TGraph* graph_stat4 = new TGraph(ArrN[p+3],ArrPt[p+3],ArrStat[p+3]); //graph_stat4->SetLineColor(8); //graph_stat4->SetMarkerColor(8); //graph_stat4->Draw("C*"); LEG = new TLegend(0.50,0.65,0.80,0.90); LEG->SetFillColor(0); LEG->SetBorderSize(0); LEG->AddEntry(graph_stat2, "I_{as}>0.10 & 1/#beta>1.10", "L"); LEG->AddEntry(graph_stat3, "I_{as}>0.30 & 1/#beta>1.05", "L"); LEG->AddEntry(graph_stat, "I_{as}>0.10 & 1/#beta>1.20", "L"); LEG->Draw(); } SaveCanvas(c1,InputPattern,string(string("Systematics_")+DataName+"_")+Name+"Stat"); delete c1; c1 = new TCanvas("c1","c1", 600, 600); TGraph* graph_statB; if(p==0)graph_statB = new TGraph(ArrN[p],ArrT [p],ArrStat[p]); if(p==1)graph_statB = new TGraph(ArrN[p],ArrI [p],ArrStat[p]); if(p==2)graph_statB = new TGraph(ArrN[p+2],ArrPt[p+2],ArrStatB[p+2]); graph_statB->SetTitle(""); graph_statB->GetYaxis()->SetTitle("Prediction #sigma_{Stat}/#mu"); graph_statB->GetYaxis()->SetTitleOffset(1.70); graph_statB->GetXaxis()->SetTitle(Title.c_str()); graph_statB->Draw("AC*"); graph_statB->GetYaxis()->SetRangeUser(0,0.25); if(p==2){ TGraph* graph_statB2 = new TGraph(ArrN[p+1],ArrPt[p+1],ArrStatB[p+1]); graph_statB2->SetLineColor(2); graph_statB2->SetMarkerColor(2); graph_statB2->Draw("C*"); TGraph* graph_statB3 = new TGraph(ArrN[p+0],ArrPt[p+0],ArrStatB[p+0]); graph_statB3->SetLineColor(4); graph_statB3->SetMarkerColor(4); graph_statB3->Draw("C*"); //TGraph* graph_statB4 = new TGraph(ArrN[p+3],ArrPt[p+3],ArrStat[p+3]); //graph_statB4->SetLineColor(8); //graph_statB4->SetMarkerColor(8); //graph_statB4->Draw("C*"); LEG = new TLegend(0.50,0.65,0.80,0.90); LEG->SetFillColor(0); LEG->SetBorderSize(0); LEG->AddEntry(graph_statB2, "I_{as}>0.10 & 1/#beta>1.10", "L"); LEG->AddEntry(graph_statB3, "I_{as}>0.30 & 1/#beta>1.05", "L"); LEG->AddEntry(graph_statB, "I_{as}>0.10 & 1/#beta>1.20", "L"); LEG->Draw(); } SaveCanvas(c1,InputPattern,string(string("Systematics_")+DataName+"_")+Name+"StatB"); delete c1; c1 = new TCanvas("c1","c1", 600, 600); TGraph* graph_syst; if(p==0)graph_syst = new TGraph(ArrN[p],ArrT [p],ArrSyst[p]); if(p==1)graph_syst = new TGraph(ArrN[p],ArrI [p],ArrSyst[p]); if(p==2)graph_syst = new TGraph(ArrN[p+2],ArrPt[p+2],ArrSyst[p+2]); graph_syst->SetTitle(""); graph_syst->GetYaxis()->SetTitle("Prediction #sigma_{Syst}/#mu"); graph_syst->GetYaxis()->SetTitleOffset(1.70); graph_syst->GetXaxis()->SetTitle(Title.c_str()); graph_syst->Draw("AC*"); graph_syst->GetXaxis()->SetRangeUser(40,100); graph_syst->GetYaxis()->SetRangeUser(0,0.25); if(p==2){ TGraph* graph_syst2 = new TGraph(ArrN[p+1],ArrPt[p+1],ArrSyst[p+1]); graph_syst2->SetLineColor(2); graph_syst2->SetMarkerColor(2); graph_syst2->Draw("C*"); TGraph* graph_syst3 = new TGraph(ArrN[p+0],ArrPt[p+0],ArrSyst[p+0]); graph_syst3->SetLineColor(4); graph_syst3->SetMarkerColor(4); graph_syst3->Draw("C*"); //TGraph* graph_syst4 = new TGraph(ArrN[p+3],ArrPt[p+3],ArrSyst[p+3]); //graph_syst4->SetLineColor(8); //graph_syst4->SetMarkerColor(8); //graph_syst4->Draw("C*"); LEG = new TLegend(0.50,0.65,0.80,0.90); LEG->SetFillColor(0); LEG->SetBorderSize(0); LEG->AddEntry(graph_syst2, "I_{as}>0.10 & 1/#beta>1.10", "L"); LEG->AddEntry(graph_syst3, "I_{as}>0.30 & 1/#beta>1.05", "L"); LEG->AddEntry(graph_syst, "I_{as}>0.10 & 1/#beta>1.20", "L"); LEG->Draw(); } SaveCanvas(c1,InputPattern,string(string("Systematics_")+DataName+"_")+Name+"Syst"); delete c1; } InputFile->Close(); } ////////////////////////////////////////////////////////////////////////////////////////////////////////// // WARNING... ALL THE FUNCTIONS BELOW HAVE NOT BEEN REVIEWED SINCE JULY UPDATE... THEY ARE PROBABLY STILL WORKING FINE, BUT CAN NOT BE SURE TILL SOMEONE TRY // PLEASE TAKE SOME TIME TO REVIEW THE FUNCTION BELOW IF YOU NEED TO USE THEM, AND MOVE THEM ABOVE THIS WARNING WHEN VALIDATED void SignalMassPlot(string InputPattern, unsigned int CutIndex){ string SavePath = InputPattern + "MassPlots/"; MakeDirectories(SavePath); string Input = InputPattern + "Histos.root"; TFile* InputFile = new TFile(Input.c_str()); for(unsigned int s=0;s<samples.size();s++){ if(samples[s].Type!=2 || !samples[s].MakePlot)continue; TH1D* Mass = GetCutIndexSliceFromTH2((TH2D*)GetObjectFromPath(InputFile, samples[s].Name + "/Mass" ), CutIndex, "SignalMass"); Mass->Scale(1.0/Mass->Integral()); char YAxisLegend[1024]; sprintf(YAxisLegend,"#tracks / %2.0f GeV/c^{2}",Mass->GetXaxis()->GetBinWidth(1)); TCanvas* c1 = new TCanvas("c1","c1", 600, 600); Mass->SetAxisRange(0,1250,"X"); // Mass->SetAxisRange(Min,Max,"Y"); Mass->SetTitle(""); // Mass->SetStats(kFALSE); Mass->GetXaxis()->SetTitle("m (GeV/c^{2})"); Mass->GetYaxis()->SetTitle(YAxisLegend); Mass->SetLineWidth(2); Mass->SetLineColor(Color[0]); Mass->SetMarkerColor(Color[0]); Mass->SetMarkerStyle(Marker[0]); Mass->Draw("HIST E1"); c1->SetLogy(true); SaveCanvas(c1,SavePath,samples[s].Name); delete c1; } } void Make2DPlot_Core(string InputPattern, unsigned int CutIndex){ TCanvas* c1; TLegend* leg; string Input = InputPattern + "Histos.root"; string outpath = InputPattern; MakeDirectories(outpath); TypeMode = TypeFromPattern(InputPattern); string LegendTitle = LegendFromType(InputPattern);; string S1 = "Gluino_8TeV_M300_f10"; //double Q1=1; string S2 = "Gluino_8TeV_M600_f10"; //double Q2=1; string S3 = "Gluino_8TeV_M900_f10"; //double Q3=1; string Da = "Data8TeV"; string outName = "2DPlots"; if(TypeMode==3) { S1 = "Gluino_8TeV_M500_f100"; S2 = "Gluino_8TeV_M1000_f100"; S3 = "Gluino_8TeV_M1500_f100"; } if(TypeMode==4){ S1 = "DY_8TeV_M400_Q1"; //Q1=1; S2 = "DY_8TeV_M400_Q2"; //Q2=2; S3 = "DY_8TeV_M400_Q3"; //Q3=3; } if(TypeMode==5){ S1 = "DY_8TeV_M300_Q1o3"; //Q1=1/3.0; S2 = "DY_8TeV_M300_Q2o3"; //Q2=2/3.0; S3 = "DY_8TeV_M600_Q2o3"; //Q3=2/3.0; } int S1i = JobIdToIndex(S1,samples); if(S1i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } int S2i = JobIdToIndex(S2,samples); if(S2i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } int S3i = JobIdToIndex(S3,samples); if(S3i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); TH1D* Signal1Mass = GetCutIndexSliceFromTH2((TH2D*)GetObjectFromPath(InputFile, S1+"/Mass" ), CutIndex, "S1Mass" ); TH2D* Signal1PtIs = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S1+"/AS_PtIs" ), CutIndex, "S1PtIs_zy"); TH2D* Signal1PIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S1+"/AS_PIm" ), CutIndex, "S1PIm_zy" ); TH2D* Signal1TOFIs= GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S1+"/AS_TOFIs"), CutIndex, "S1TIs_zy" ); TH2D* Signal1TOFIm= GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S1+"/AS_TOFIm"), CutIndex, "S1TIm_zy" ); TH2D* Signal1PtTOF= (TH2D*)GetObjectFromPath(InputFile, S1+"/BS_PtTOF" ); Signal1PtTOF->RebinX(2); Signal1PtTOF->RebinY(2); TH1D* Signal2Mass = GetCutIndexSliceFromTH2((TH2D*)GetObjectFromPath(InputFile, S2+"/Mass" ), CutIndex, "S2Mass" ); TH2D* Signal2PtIs = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S2+"/AS_PtIs" ), CutIndex, "S2PtIs_zy"); TH2D* Signal2PIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S2+"/AS_PIm" ), CutIndex, "S2PIm_zy" ); TH2D* Signal2TOFIs= GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S2+"/AS_TOFIs"), CutIndex, "S2TIs_zy" ); TH2D* Signal2TOFIm= GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S2+"/AS_TOFIm"), CutIndex, "S2TIm_zy" ); TH2D* Signal2PtTOF= (TH2D*)GetObjectFromPath(InputFile, S2+"/BS_PtTOF" ); Signal2PtTOF->RebinX(2); Signal2PtTOF->RebinY(2); TH1D* Signal3Mass = GetCutIndexSliceFromTH2((TH2D*)GetObjectFromPath(InputFile, S3+"/Mass" ), CutIndex, "S3Mass" ); TH2D* Signal3PtIs = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S3+"/AS_PtIs" ), CutIndex, "S3PtIs_zy"); TH2D* Signal3PIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S3+"/AS_PIm" ), CutIndex, "S3PIm_zy" ); TH2D* Signal3TOFIs= GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S3+"/AS_TOFIs"), CutIndex, "S3TIs_zy" ); TH2D* Signal3TOFIm= GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S3+"/AS_TOFIm"), CutIndex, "S3TIm_zy" ); TH2D* Signal3PtTOF= (TH2D*)GetObjectFromPath(InputFile, S3+"/BS_PtTOF" ); Signal3PtTOF->RebinX(2); Signal3PtTOF->RebinY(2); TH2D* Data_PtIs = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, Da+"/AS_PtIs" ), CutIndex); TH2D* Data_PIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, Da+"/AS_PIm" ), CutIndex); TH2D* Data_TOFIs = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, Da+"/AS_TOFIs"), CutIndex); TH2D* Data_TOFIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, Da+"/AS_TOFIm"), CutIndex); TH2D* Data_PtTOF = (TH2D*)GetObjectFromPath(InputFile, Da+"/BS_PtTOF" ); Data_PtTOF->RebinX(2); Data_PtTOF->RebinY(2); // TH2D* Data_PIm_075 = (TH2D*)Data_PIm->Clone(); Data_PIm_075->Reset(); // TH2D* Data_PIm_150 = (TH2D*)Data_PIm->Clone(); Data_PIm_150->Reset(); // TH2D* Data_PIm_300 = (TH2D*)Data_PIm->Clone(); Data_PIm_300->Reset(); // TH2D* Data_PIm_450 = (TH2D*)Data_PIm->Clone(); Data_PIm_450->Reset(); // TH2D* Data_PIm_All = (TH2D*)Data_PIm->Clone(); Data_PIm_All->Reset(); // for(unsigned int i=0;i<(unsigned int)Data_PIm->GetNbinsX();i++){ // for(unsigned int j=0;j<(unsigned int)Data_PIm->GetNbinsY();j++){ // if(Data_PIm->GetBinContent(i,j)<=0)continue; // double M = GetMass(Data_PIm->GetXaxis ()->GetBinCenter(i), Data_PIm->GetYaxis ()->GetBinCenter(j), false); // if(isnan(M))continue; // if (M<100){ Data_PIm_075->SetBinContent(i,j, Data_PIm->GetBinContent(i,j) ); } // else if(M<200){ Data_PIm_150->SetBinContent(i,j, Data_PIm->GetBinContent(i,j) ); } // else if(M<300){ Data_PIm_300->SetBinContent(i,j, Data_PIm->GetBinContent(i,j) ); } // else if(M<395){ Data_PIm_450->SetBinContent(i,j, Data_PIm->GetBinContent(i,j) ); } // else { Data_PIm_All->SetBinContent(i,j, Data_PIm->GetBinContent(i,j) ); } // }} Signal1Mass = (TH1D*) Signal1Mass->Rebin(2); Signal2Mass = (TH1D*) Signal2Mass->Rebin(2); Signal3Mass = (TH1D*) Signal3Mass->Rebin(2); double Min = 1E-3; double Max = 1E4; char YAxisLegend[1024]; sprintf(YAxisLegend,"#tracks / %2.0f GeV/c^{2}",Signal1Mass->GetXaxis()->GetBinWidth(1)); c1 = new TCanvas("c1","c1", 600, 600); Signal1Mass->SetAxisRange(0,1250,"X"); Signal1Mass->SetAxisRange(Min,Max,"Y"); Signal1Mass->SetTitle(""); Signal1Mass->SetStats(kFALSE); Signal1Mass->GetXaxis()->SetTitle("m (GeV/c^{2})"); Signal1Mass->GetYaxis()->SetTitle(YAxisLegend); Signal1Mass->SetLineWidth(2); Signal1Mass->SetLineColor(Color[0]); Signal1Mass->SetMarkerColor(Color[0]); Signal1Mass->SetMarkerStyle(Marker[0]); Signal1Mass->Draw("HIST E1"); Signal2Mass->Draw("HIST E1 same"); Signal2Mass->SetLineColor(Color[1]); Signal2Mass->SetMarkerColor(Color[1]); Signal2Mass->SetMarkerStyle(Marker[1]); Signal2Mass->SetLineWidth(2); Signal3Mass->SetLineWidth(2); Signal3Mass->SetLineColor(Color[2]); Signal3Mass->SetMarkerColor(Color[2]); Signal3Mass->SetMarkerStyle(Marker[2]); Signal3Mass->Draw("HIST E1 same"); c1->SetLogy(true); TLine* line300 = new TLine(300, Min, 300, Max); line300->SetLineWidth(2); line300->SetLineColor(Color[0]); line300->SetLineStyle(2); line300->Draw("same"); TLine* line500 = new TLine(500, Min, 500, Max); line500->SetLineWidth(2); line500->SetLineColor(Color[1]); line500->SetLineStyle(2); line500->Draw("same"); TLine* line800 = new TLine(800, Min, 800, Max); line800->SetLineWidth(2); line800->SetLineColor(Color[2]); line800->SetLineStyle(2); line800->Draw("same"); leg = new TLegend(0.80,0.93,0.80 - 0.20,0.93 - 6*0.03); leg->SetFillColor(0); leg->SetBorderSize(0); leg->AddEntry(Signal1Mass, samples[S1i].Legend.c_str() ,"P"); leg->AddEntry(Signal2Mass, samples[S2i].Legend.c_str() ,"P"); leg->AddEntry(Signal3Mass, samples[S3i].Legend.c_str() ,"P"); leg->Draw(); DrawPreliminary((LegendTitle + " - Simulation").c_str(), SQRTS, ""); SaveCanvas(c1, outpath, outName + "_Mass"); delete c1; c1 = new TCanvas("c1","c1", 600, 600); c1->SetLogz(true); Data_PtIs->SetTitle(""); Data_PtIs->SetStats(kFALSE); Data_PtIs->GetXaxis()->SetTitle("p (GeV/c)"); Data_PtIs->GetYaxis()->SetTitle(dEdxS_Legend.c_str()); Data_PtIs->SetAxisRange(0,1750,"X"); Data_PtIs->SetMarkerSize (0.2); Data_PtIs->SetMarkerColor(Color[4]); Data_PtIs->SetFillColor(Color[4]); Data_PtIs->Draw("COLZ"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1, outpath, outName + "_Data_PtIs", true); delete c1; c1 = new TCanvas("c1","c1", 600, 600); c1->SetLogz(true); Data_PIm->SetTitle(""); Data_PIm->SetStats(kFALSE); Data_PIm->GetXaxis()->SetTitle("p (GeV/c)"); Data_PIm->GetYaxis()->SetTitle(dEdxM_Legend.c_str()); Data_PIm->SetAxisRange(0,1750,"X"); Data_PIm->SetAxisRange(0,TypeMode==5?3:15,"Y"); Data_PIm->SetMarkerSize (0.2); Data_PIm->SetMarkerColor(Color[4]); Data_PIm->SetFillColor(Color[4]); Data_PIm->Draw("COLZ"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1, outpath, outName + "_Data_PIm", true); delete c1; c1 = new TCanvas("c1","c1", 600, 600); c1->SetLogz(true); Data_TOFIs->SetTitle(""); Data_TOFIs->SetStats(kFALSE); Data_TOFIs->GetXaxis()->SetTitle("1/#beta_{TOF}"); Data_TOFIs->GetYaxis()->SetTitle(dEdxS_Legend.c_str()); Data_TOFIs->SetAxisRange(0,1750,"X"); Data_TOFIs->SetMarkerSize (0.2); Data_TOFIs->SetMarkerColor(Color[4]); Data_TOFIs->SetFillColor(Color[4]); Data_TOFIs->Draw("COLZ"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1, outpath, outName + "_Data_TOFIs", true); delete c1; c1 = new TCanvas("c1","c1", 600, 600); c1->SetLogz(true); Data_TOFIm->SetTitle(""); Data_TOFIm->SetStats(kFALSE); Data_TOFIm->GetXaxis()->SetTitle("1/#beta_{TOF}"); Data_TOFIm->GetYaxis()->SetTitle(dEdxM_Legend.c_str()); Data_TOFIm->SetAxisRange(0,15,"Y"); Data_TOFIm->SetMarkerSize (0.2); Data_TOFIm->SetMarkerColor(Color[4]); Data_TOFIm->SetFillColor(Color[4]); Data_TOFIm->Draw("COLZ"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1, outpath, outName + "_Data_TOFIm", true); delete c1; c1 = new TCanvas("c1","c1", 600, 600); c1->SetLogz(true); Data_PtTOF->SetTitle(""); Data_PtTOF->SetStats(kFALSE); Data_PtTOF->GetXaxis()->SetTitle("p_{T} (GeV/c)"); Data_PtTOF->GetYaxis()->SetTitle("1/#beta"); Data_PtTOF->SetAxisRange(0,1750,"X"); Data_PtTOF->SetMarkerSize (0.2); Data_PtTOF->SetMarkerColor(Color[4]); Data_PtTOF->SetFillColor(Color[4]); Data_PtTOF->Draw("COLZ"); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1, outpath, outName + "_Data_PtTOF", false); delete c1; c1 = new TCanvas("c1","c1", 600, 600); Signal3PtIs->SetTitle(""); Signal3PtIs->SetStats(kFALSE); Signal3PtIs->GetXaxis()->SetTitle("p (GeV/c)"); Signal3PtIs->GetYaxis()->SetTitle(dEdxS_Legend.c_str()); Signal3PtIs->SetAxisRange(0,1750,"X"); Signal3PtIs->Scale(1/Signal3PtIs->Integral()); Signal3PtIs->SetMarkerSize (0.2); Signal3PtIs->SetMarkerColor(Color[2]); Signal3PtIs->SetFillColor(Color[2]); Signal3PtIs->Draw("BOX"); Signal2PtIs->Scale(1/Signal2PtIs->Integral()); Signal2PtIs->SetMarkerSize (0.2); Signal2PtIs->SetMarkerColor(Color[1]); Signal2PtIs->SetFillColor(Color[1]); Signal2PtIs->Draw("BOX same"); Signal1PtIs->Scale(1/Signal1PtIs->Integral()); Signal1PtIs->SetMarkerSize (0.2); Signal1PtIs->SetMarkerColor(Color[0]); Signal1PtIs->SetFillColor(Color[0]); Signal1PtIs->Draw("BOX same"); leg = new TLegend(0.80,0.93,0.80 - 0.40,0.93 - 6*0.03); leg->SetFillColor(0); leg->SetBorderSize(0); leg->AddEntry(Signal1PtIs, samples[S1i].Legend.c_str() ,"F"); leg->AddEntry(Signal2PtIs, samples[S2i].Legend.c_str() ,"F"); leg->AddEntry(Signal3PtIs, samples[S3i].Legend.c_str() ,"F"); leg->Draw(); DrawPreliminary(LegendTitle + " - Simulation", SQRTS, ""); SaveCanvas(c1, outpath, outName + "_PtIs", true); delete c1; c1 = new TCanvas("c1","c1", 600, 600); // c1->SetLogz(true); Signal1PIm->SetTitle(""); Signal1PIm->SetStats(kFALSE); Signal1PIm->GetXaxis()->SetTitle("p (GeV/c)"); Signal1PIm->GetYaxis()->SetTitle(dEdxM_Legend.c_str()); Signal1PIm->SetAxisRange(0,1750,"X"); Signal1PIm->SetAxisRange(0,TypeMode==5?3:15,"Y"); Signal1PIm->Scale(10/Signal1PIm->Integral()); Signal1PIm->SetMarkerSize (0.2); Signal1PIm->SetMarkerColor(Color[2]); Signal1PIm->SetFillColor(Color[2]); Signal1PIm->Draw("BOX"); Signal2PIm->Scale(10/Signal2PIm->Integral()); Signal2PIm->SetMarkerSize (0.2); Signal2PIm->SetMarkerColor(Color[1]); Signal2PIm->SetFillColor(Color[1]); Signal2PIm->Draw("BOX same"); Signal3PIm->Scale(10/Signal3PIm->Integral()); Signal3PIm->SetMarkerSize (0.2); Signal3PIm->SetMarkerColor(Color[0]); Signal3PIm->SetFillColor(Color[0]); Signal3PIm->Draw("BOX same"); // TF1* MassLine800 = GetMassLineQ(samples[S3i].Mass, Q3, true); // MassLine800->SetLineColor(kGray+3); // MassLine800->SetLineWidth(2); // MassLine800->Draw("same"); // TF1* MassLine500 = GetMassLineQ(samples[S2i].Mass, Q2, true); // MassLine500->SetLineColor(kBlue-7); // MassLine500->SetLineWidth(2); // MassLine500->Draw("same"); // TF1* MassLine300 = GetMassLineQ(samples[S1i].Mass, Q1, true); // MassLine300->SetLineColor(kRed-7); // MassLine300->SetLineWidth(2); // MassLine300->Draw("same"); leg = new TLegend(0.80,0.93,0.80 - 0.40,0.93 - 6*0.03); leg->SetFillColor(0); leg->SetBorderSize(0); leg->AddEntry(Signal1PIm, samples[S1i].Legend.c_str() ,"F"); leg->AddEntry(Signal2PIm, samples[S2i].Legend.c_str() ,"F"); leg->AddEntry(Signal3PIm, samples[S3i].Legend.c_str() ,"F"); leg->Draw(); DrawPreliminary((LegendTitle + " - Simulation").c_str(), SQRTS, ""); SaveCanvas(c1, outpath, outName + "_PIm", true); delete c1; c1 = new TCanvas("c1","c1", 600, 600); Signal3TOFIs->SetTitle(""); Signal3TOFIs->SetStats(kFALSE); Signal3TOFIs->GetXaxis()->SetTitle("1/#beta_{TOF}"); Signal3TOFIs->GetYaxis()->SetTitle(dEdxS_Legend.c_str()); Signal3TOFIs->SetAxisRange(0,1750,"X"); Signal3TOFIs->Scale(1/Signal3TOFIs->Integral()); Signal3TOFIs->SetMarkerSize (0.2); Signal3TOFIs->SetMarkerColor(Color[2]); Signal3TOFIs->SetFillColor(Color[2]); Signal3TOFIs->Draw("BOX"); Signal2TOFIs->Scale(1/Signal2TOFIs->Integral()); Signal2TOFIs->SetMarkerSize (0.2); Signal2TOFIs->SetMarkerColor(Color[1]); Signal2TOFIs->SetFillColor(Color[1]); Signal2TOFIs->Draw("BOX same"); Signal1TOFIs->Scale(1/Signal1TOFIs->Integral()); Signal1TOFIs->SetMarkerSize (0.2); Signal1TOFIs->SetMarkerColor(Color[0]); Signal1TOFIs->SetFillColor(Color[0]); Signal1TOFIs->Draw("BOX same"); leg = new TLegend(0.80,0.93,0.80 - 0.40,0.93 - 6*0.03); leg->SetFillColor(0); leg->SetBorderSize(0); leg->AddEntry(Signal1TOFIs, samples[S1i].Legend.c_str() ,"F"); leg->AddEntry(Signal2TOFIs, samples[S2i].Legend.c_str() ,"F"); leg->AddEntry(Signal3TOFIs, samples[S3i].Legend.c_str() ,"F"); leg->Draw(); DrawPreliminary((LegendTitle + " - Simulation").c_str(), SQRTS, ""); SaveCanvas(c1, outpath, outName + "_TOFIs", true); delete c1; c1 = new TCanvas("c1","c1", 600, 600); Signal3TOFIm->SetTitle(""); Signal3TOFIm->SetStats(kFALSE); Signal3TOFIm->GetXaxis()->SetTitle("1/#beta_{TOF}"); Signal3TOFIm->GetYaxis()->SetTitle(dEdxM_Legend.c_str()); Signal3TOFIm->SetAxisRange(0,1750,"X"); Signal3TOFIm->SetAxisRange(0,15,"Y"); Signal3TOFIm->Scale(1/Signal3TOFIm->Integral()); Signal3TOFIm->SetMarkerSize (0.2); Signal3TOFIm->SetMarkerColor(Color[2]); Signal3TOFIm->SetFillColor(Color[2]); Signal3TOFIm->Draw("BOX"); Signal2TOFIm->Scale(1/Signal2TOFIm->Integral()); Signal2TOFIm->SetMarkerSize (0.2); Signal2TOFIm->SetMarkerColor(Color[1]); Signal2TOFIm->SetFillColor(Color[1]); Signal2TOFIm->Draw("BOX same"); Signal1TOFIm->Scale(1/Signal1TOFIm->Integral()); Signal1TOFIm->SetMarkerSize (0.2); Signal1TOFIm->SetMarkerColor(Color[0]); Signal1TOFIm->SetFillColor(Color[0]); Signal1TOFIm->Draw("BOX same"); leg = new TLegend(0.80,0.93,0.80 - 0.40,0.93 - 6*0.03); leg->SetFillColor(0); leg->SetBorderSize(0); leg->AddEntry(Signal1TOFIm, samples[S1i].Legend.c_str() ,"F"); leg->AddEntry(Signal2TOFIm, samples[S2i].Legend.c_str() ,"F"); leg->AddEntry(Signal3TOFIm, samples[S3i].Legend.c_str() ,"F"); leg->Draw(); DrawPreliminary(LegendTitle + " - Simulation", SQRTS, ""); SaveCanvas(c1, outpath, outName + "_TOFIm", true); delete c1; c1 = new TCanvas("c1","c1", 600, 600); Signal3PtTOF->SetTitle(""); Signal3PtTOF->SetStats(kFALSE); Signal3PtTOF->GetXaxis()->SetTitle("p_{T} (GeV/c)"); Signal3PtTOF->GetYaxis()->SetTitle("1/#beta"); Signal3PtTOF->SetAxisRange(0,1750,"X"); Signal3PtTOF->Scale(1/Signal3PtTOF->Integral()); Signal3PtTOF->SetMarkerSize (0.2); Signal3PtTOF->SetMarkerColor(Color[2]); Signal3PtTOF->SetFillColor(Color[2]); Signal3PtTOF->Draw("BOX"); Signal2PtTOF->Scale(1/Signal2PtTOF->Integral()); Signal2PtTOF->SetMarkerSize (0.2); Signal2PtTOF->SetMarkerColor(Color[1]); Signal2PtTOF->SetFillColor(Color[1]); Signal2PtTOF->Draw("BOX same"); Signal1PtTOF->Scale(1/Signal1PtTOF->Integral()); Signal1PtTOF->SetMarkerSize (0.2); Signal1PtTOF->SetMarkerColor(Color[0]); Signal1PtTOF->SetFillColor(Color[0]); Signal1PtTOF->Draw("BOX same"); leg = new TLegend(0.80,0.93,0.80 - 0.40,0.93 - 6*0.03); leg->SetFillColor(0); leg->SetBorderSize(0); leg->AddEntry(Signal1PtTOF, samples[S1i].Legend.c_str() ,"F"); leg->AddEntry(Signal2PtTOF, samples[S2i].Legend.c_str() ,"F"); leg->AddEntry(Signal3PtTOF, samples[S3i].Legend.c_str() ,"F"); leg->Draw(); DrawPreliminary(LegendTitle + " - Simulation", SQRTS, ""); SaveCanvas(c1, outpath, outName + "_PtTOF", false); delete c1; // c1 = new TCanvas("c1","c1", 600, 600); // Data_PIm_075->SetTitle(""); // Data_PIm_075->SetStats(kFALSE); // Data_PIm_075->GetXaxis()->SetTitle("p (GeV/c)"); // Data_PIm_075->GetYaxis()->SetTitle(dEdxM_Legend.c_str()); // Data_PIm_075->SetAxisRange(3,10,"Y"); // Data_PIm_075->SetAxisRange(0,2000,"X"); // Data_PIm_075->SetMarkerSize (0.6); // Data_PIm_075->SetMarkerColor(Color[4]); // Data_PIm_075->SetMarkerStyle(Marker[4]); // Data_PIm_075->SetFillColor(Color[4]); // Data_PIm_075->Draw(""); // Data_PIm_150->SetMarkerSize (0.8); // Data_PIm_150->SetMarkerColor(Color[3]); // Data_PIm_150->SetMarkerStyle(Marker[3]); // Data_PIm_150->SetFillColor(Color[3]); // Data_PIm_150->Draw("same"); // Data_PIm_300->SetMarkerSize (1.0); // Data_PIm_300->SetMarkerColor(Color[2]); // Data_PIm_300->SetMarkerStyle(Marker[2]); // Data_PIm_300->SetFillColor(Color[2]); // Data_PIm_300->Draw("same"); // Data_PIm_450->SetMarkerSize (1.2); // Data_PIm_450->SetMarkerColor(Color[1]); // Data_PIm_450->SetMarkerStyle(Marker[1]); // Data_PIm_450->SetFillColor(Color[1]); // Data_PIm_450->Draw("same"); // Data_PIm_All->SetMarkerSize (1.4); // Data_PIm_All->SetMarkerColor(Color[0]); // Data_PIm_All->SetFillColor(Color[0]); // Data_PIm_All->Draw("same"); // for(double m=100;m<1000;m+=100){ // TF1* MassLine = GetMassLine(m, false); // MassLine->SetLineColor(1); // MassLine->SetLineWidth(1); // MassLine->Draw("same"); // } // leg = new TLegend(0.80,0.93,0.80 - 0.30,0.93 - 6*0.03); // leg->SetHeader(LegendFromType(InputPattern).c_str()); // leg->SetFillColor(0); // leg->SetBorderSize(0); // leg->AddEntry(Data_PIm_075, "M < 100 GeV","P"); // leg->AddEntry(Data_PIm_150, "100 < M < 200 GeV","P"); // leg->AddEntry(Data_PIm_300, "200 < M < 300 GeV","P"); // leg->AddEntry(Data_PIm_450, "300 < M < 400 GeV","P"); // leg->AddEntry(Data_PIm_All, "400 < M GeV" ,"P"); // leg->Draw(); // DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); // SaveCanvas(c1, outpath, "Data_PIm_Colored", true); // delete c1; } void MakeExpLimitpLot(string Input, string Output){ TH2D* ExpLimitPlot = new TH2D("ExpLimitPlot","ExpLimitPlot", 10,37.5,87.5,13,0.0875,0.4175); std::vector<double> VectPt; std::vector<double> VectI; std::vector<double> VectExpLim; FILE* pFile = fopen(Input.c_str(),"r"); if(!pFile){ printf("Not Found: %s\n",Input.c_str()); return; } unsigned int Index; double Pt, I, TOF, MassMin, MassMax; double NData, NPred, NPredErr, SignalEff; double ExpLimit; char Model[256], Tmp[2048]; while ( ! feof (pFile) ){ fscanf(pFile,"%s Testing CutIndex= %d (Pt>%lf I>%lf TOF>%lf) %lf<M<%lf Ndata=%lE NPred=%lE+-%lE SignalEff=%lf --> %lE expected",Model,&Index,&Pt,&I,&TOF,&MassMin,&MassMax,&NData,&NPred,&NPredErr,&SignalEff,&ExpLimit); fgets(Tmp, 256 , pFile); // if(Pt<80 && I<0.38)printf("%s Testing CutIndex= %d (Pt>%f I>%f TOF>%f) %f<M<%f Ndata=%E NPred=%E+-%E SignalEff=%f --> %E expected %s",Model,Index,Pt,I,TOF,MassMin,MassMax,NData,NPred,NPredErr,SignalEff,ExpLimit, Tmp); // ExpLimitPlot->SetBinContent(PtMap[Pt],IsMap[I],ExpLimit); ExpLimitPlot->Fill(Pt,I,ExpLimit); } fclose(pFile); TCanvas* c1 = new TCanvas("c1","c1",600,600); c1->SetLogz(true); ExpLimitPlot->SetTitle(""); ExpLimitPlot->SetStats(kFALSE); ExpLimitPlot->GetXaxis()->SetTitle("Pt Cut"); ExpLimitPlot->GetYaxis()->SetTitle("I Cut"); ExpLimitPlot->GetXaxis()->SetTitleOffset(1.1); ExpLimitPlot->GetYaxis()->SetTitleOffset(1.70); ExpLimitPlot->GetXaxis()->SetNdivisions(505); ExpLimitPlot->GetYaxis()->SetNdivisions(505); ExpLimitPlot->SetMaximum(1E0); ExpLimitPlot->SetMinimum(1E-2); ExpLimitPlot->Draw("COLZ"); c1->SaveAs(Output.c_str()); delete c1; return; } void CosmicBackgroundSystematic(string InputPattern, string DataType){ if(DataType.find("7TeV")!=string::npos){SQRTS=7.0;}else{SQRTS=8.0;} string SavePath = InputPattern; MakeDirectories(SavePath); TCanvas* c1; TH1** Histos = new TH1*[10]; std::vector<string> legend; string LegendTitle = LegendFromType(InputPattern); TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); TH1D* HCuts_Pt = (TH1D*)GetObjectFromPath(InputFile, "HCuts_Pt"); TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, "HCuts_TOF"); TH2D* H_D_DzSidebands = ((TH2D*)GetObjectFromPath(InputFile, "Data" + DataType + "/H_D_DzSidebands")); TH2D* H_D_DzSidebands_Cosmic = (TH2D*)GetObjectFromPath(InputFile, "Cosmic8TeV/H_D_DzSidebands"); TH1D* H_D_Cosmic = (TH1D*)GetObjectFromPath(InputFile, "Cosmic8TeV/H_D"); //TH1D* H_D_Data = (TH1D*)GetObjectFromPath(InputFile, "Data8TeV/H_D"); std::vector<int> Index; std::vector<int> Plot; for(int CutIndex=0; CutIndex<HCuts_Pt->GetNbinsX(); CutIndex++) { if(fabs(HCuts_TOF->GetBinContent(CutIndex+1)-1.25)<0.001) {Index.push_back(CutIndex); Plot.push_back(0);} if(fabs(HCuts_TOF->GetBinContent(CutIndex+1)-1.3)<0.001) {Index.push_back(CutIndex); Plot.push_back(1);} if(fabs(HCuts_TOF->GetBinContent(CutIndex+1)-1.375)<0.001) {Index.push_back(CutIndex); Plot.push_back(2);} } const int TimeRegions=3; TH1F *Pred[TimeRegions*DzRegions]; TH1F *StatSyst[TimeRegions]; TH1F *Stat[TimeRegions]; TH1F *Syst[TimeRegions]; string Preds[TimeRegions] = {"125", "130", "135"}; string RegionNames[DzRegions]={"Region0","Region1","Region2","Region3","Region4", "Region5"}; string LegendNames[DzRegions]={"dz < 20 cm","20 cm < dz < 30 cm","30 cm < dz < 50 cm","50 cm < dz < 70 cm","70 cm < dz < 120 cm", "dz > 120 cm"}; for(int i=0; i<TimeRegions; i++) { StatSyst[i] = new TH1F(("StatSyst_TOF" + Preds[i]).c_str(), "StatSyst_TOF100", 9, 95, 365); Stat[i] = new TH1F(("Stat_TOF" + Preds[i]).c_str(), "Stat_TOF100", 9, 95, 365); Syst[i] = new TH1F(("Syst_TOF" + Preds[i]).c_str(), "Syst_TOF110", 9, 95, 365); for(int Region=0; Region<DzRegions; Region++) { string Name="Pred_TOF" + Preds[i] + "_"+ RegionNames[Region]; Pred[i*DzRegions+Region] = new TH1F(Name.c_str(), Name.c_str(), 9, 95, 365); } } const double alpha = 1 - 0.6827; //cout << endl << endl; for(unsigned int i=0; i<Index.size(); i++) { int CutIndex=Index[i]; double D_Cosmic = H_D_Cosmic->GetBinContent(CutIndex+1); double D_Cosmic_Var = pow(ROOT::Math::gamma_quantile_c(alpha/2,D_Cosmic+1,1) - D_Cosmic,2); int Bin=Pred[0]->FindBin(HCuts_Pt->GetBinContent(CutIndex+1)); double Sigma = 0; double Mean = 0; int N=0; for(int Region=2; Region<DzRegions; Region++) { double D_Sideband = H_D_DzSidebands->GetBinContent(CutIndex+1, Region+1); double D_Sideband_Cosmic = H_D_DzSidebands_Cosmic->GetBinContent(CutIndex+1, Region+1); double NPred = D_Sideband * D_Cosmic / D_Sideband_Cosmic; double D_Sideband_Cosmic_Var = pow(ROOT::Math::gamma_quantile_c(alpha/2,D_Sideband_Cosmic+1,1) - D_Sideband_Cosmic,2); double D_Sideband_Var = pow(ROOT::Math::gamma_quantile_c(alpha/2,D_Sideband+1,1) - D_Sideband,2); double NPredErr = sqrt( (pow(D_Cosmic/D_Sideband_Cosmic,2)*D_Sideband_Var) + (pow(D_Sideband/D_Sideband_Cosmic,2)*D_Cosmic_Var) + (pow((D_Cosmic*(D_Sideband)/(D_Sideband_Cosmic*D_Sideband_Cosmic)),2)*D_Sideband_Cosmic_Var) ); Pred[Plot[i]*DzRegions + Region]->SetBinContent(Bin, NPred); Pred[Plot[i]*DzRegions + Region]->SetBinError(Bin, NPredErr); Mean+=NPred/pow(NPredErr,2); Sigma+=1/pow(NPredErr,2); N++; if(fabs(HCuts_TOF->GetBinContent(CutIndex+1)-1.375)<0.001 && fabs(HCuts_Pt->GetBinContent(CutIndex+1)-230)<0.001) { cout << endl << "D Sideband " << D_Sideband << " D_Cosmic " << D_Cosmic << " D_Sideband_Cosmic " << D_Sideband_Cosmic << endl; cout << "For Dz region " << LegendNames[Region] << " NPred " << NPred << " +- " << NPredErr << endl; } } Mean = Mean/Sigma; Sigma = sqrt(Sigma); double SUM=0, STAT=0, SYST=0; if(fabs(HCuts_TOF->GetBinContent(CutIndex+1)-1.375)<0.001 && fabs(HCuts_Pt->GetBinContent(CutIndex+1)-230)<0.001) cout << endl << "Stat before " << STAT << endl; for(int p=0;p<DzRegions;p++){ SUM += pow(Pred[Plot[i]*DzRegions + p]->GetBinContent(Bin)-Mean,2); STAT += pow(Pred[Plot[i]*DzRegions + p]->GetBinError(Bin),2); if(fabs(HCuts_TOF->GetBinContent(CutIndex+1)-1.375)<0.001 && fabs(HCuts_Pt->GetBinContent(CutIndex+1)-230)<0.001) cout << "STAT with " << p << " Preds " << STAT << endl; } if(fabs(HCuts_TOF->GetBinContent(CutIndex+1)-1.375)<0.001 && fabs(HCuts_Pt->GetBinContent(CutIndex+1)-230)<0.001) cout << "Sum of stat " << STAT << endl; SUM = sqrt(SUM/(N-1)); STAT = sqrt(STAT/(N)); SYST = sqrt(SUM*SUM - STAT*STAT); if(fabs(HCuts_TOF->GetBinContent(CutIndex+1)-1.375)<0.001 && fabs(HCuts_Pt->GetBinContent(CutIndex+1)-230)<0.001) { cout << "Average " << STAT << endl; cout << "Mean " << Mean << " Sigma " << Sigma << " Stat " << STAT/Mean << " StatSyst " << SUM/Mean << "Syst " << SYST/Mean << endl;} Stat[Plot[i]]->SetBinContent(Bin, STAT/Mean); StatSyst[Plot[i]]->SetBinContent(Bin, SUM/Mean); Syst[Plot[i]]->SetBinContent(Bin, SYST/Mean); } //cout << endl << endl; for(int i=0; i<TimeRegions; i++) { c1 = new TCanvas("c1","c1,",600,600); legend.clear(); for(int Region=2; Region<DzRegions; Region++) { Histos[Region-2] = Pred[i*DzRegions + Region]; legend.push_back(LegendNames[Region]); } DrawSuperposedHistos((TH1**)Histos, legend, "E1", "Pt Cut", "Predicted", 0,0, 0,0); DrawLegend((TObject**)Histos,legend,LegendTitle,"P",0.8, 0.9, 0.4, 0.05); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,(DataType + "CosmicPrediction_TOF" + Preds[i]).c_str()); delete c1; } c1 = new TCanvas("c1","c1,",600,600); legend.clear(); for(int i=0; i<TimeRegions; i++) { Histos[i] = StatSyst[i]; legend.push_back(Preds[i]); } DrawSuperposedHistos((TH1**)Histos, legend, "", "Pt Cut", "Stat+Syst Rel. Error", 0,0, 0,1.4); DrawLegend((TObject**)Histos,legend,LegendTitle,"P"); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,DataType +"CosmicStatSyst"); delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); for(int i=0; i<TimeRegions; i++) { Histos[i] = Stat[i]; legend.push_back(Preds[i]); } DrawSuperposedHistos((TH1**)Histos, legend, "", "Pt Cut", "Stat Rel. Error", 0,0, 0,1.4); DrawLegend((TObject**)Histos,legend,LegendTitle,"P"); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,DataType +"CosmicStat"); delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); for(int i=0; i<TimeRegions; i++) { Histos[i] = Syst[i]; legend.push_back(Preds[i]); } DrawSuperposedHistos((TH1**)Histos, legend, "", "Pt Cut", "Syst Rel. Error", 0,0, 0,1.4); DrawLegend((TObject**)Histos,legend,LegendTitle,"P"); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,DataType +"CosmicSyst"); delete c1; } void CheckPrediction(string InputPattern, string HistoSuffix, string DataType){ TypeMode = TypeFromPattern(InputPattern); if(TypeMode==0)return; if(DataType.find("7TeV")!=string::npos){SQRTS=7.0;}else{SQRTS=8.0;} std::vector<string> legend; string LegendTitle = LegendFromType(InputPattern); string SavePath = InputPattern; MakeDirectories(SavePath); TCanvas* c1; TH1** Histos = new TH1*[100]; TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); TH1D* HCuts_Pt = (TH1D*)GetObjectFromPath(InputFile, string("HCuts_Pt") + HistoSuffix); TH1D* HCuts_I = (TH1D*)GetObjectFromPath(InputFile, string("HCuts_I") + HistoSuffix); TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, string("HCuts_TOF") + HistoSuffix); TH1D* H_D = (TH1D*)GetObjectFromPath(InputFile, string(DataType+"/H_D") + HistoSuffix); TH1D* H_P = (TH1D*)GetObjectFromPath(InputFile, string(DataType+"/H_P") + HistoSuffix); std::vector<int> Index; std::vector<int> Plot; std::vector<double> TOFCuts; double TOFCutMax=0, TOFCutMin=9999; map<std::pair<double, double>,int> CutMap; int countPlots=0; for(int CutIndex=1; CutIndex<HCuts_TOF->GetNbinsX(); CutIndex++) { TOFCuts.push_back(HCuts_TOF->GetBinContent(CutIndex+1)); if(HCuts_TOF->GetBinContent(CutIndex+1)<TOFCutMin) TOFCutMin=HCuts_TOF->GetBinContent(CutIndex+1); if(HCuts_TOF->GetBinContent(CutIndex+1)>TOFCutMax) TOFCutMax=HCuts_TOF->GetBinContent(CutIndex+1); std::pair<double, double> key(HCuts_I->GetBinContent(CutIndex+1), HCuts_Pt->GetBinContent(CutIndex+1)); //New combination of TOF and I cuts if(CutMap.find(key)==CutMap.end()) { CutMap[key]=countPlots; countPlots++; } } std::vector<TH1D*> Pred; std::vector<TH1D*> Data; std::vector<TH1D*> Ratio; for(int i=0; i<countPlots; i++) { char DataName[1024]; sprintf(DataName,"Data_%i",i); TH1D* TempData = new TH1D(DataName, DataName, TOFCuts.size(), TOFCutMin, TOFCutMax); Data.push_back(TempData); char PredName[1024]; sprintf(PredName,"Pred_%i",i); TH1D* TempPred = new TH1D(PredName, PredName, TOFCuts.size(), TOFCutMin, TOFCutMax); Pred.push_back(TempPred); char RatioName[1024]; sprintf(RatioName,"Ratio_%i",i); TH1D* TempRatio = new TH1D(RatioName, RatioName, TOFCuts.size(), TOFCutMin, TOFCutMax); Ratio.push_back(TempRatio); } for(int CutIndex=1; CutIndex<HCuts_TOF->GetNbinsX(); CutIndex++) { std::pair<double, double> key(HCuts_I->GetBinContent(CutIndex+1), HCuts_Pt->GetBinContent(CutIndex+1)); int plot = CutMap.find(key)->second; int bin = Data[plot]->FindBin(HCuts_TOF->GetBinContent(CutIndex+1)); double D = H_D->GetBinContent(CutIndex+1); Data[plot]->SetBinContent(bin, D); double P = H_P->GetBinContent(CutIndex+1); double Perr = H_P->GetBinError(CutIndex+1); Pred[plot]->SetBinContent(bin, P); Pred[plot]->SetBinError(bin, Perr); Ratio[plot]->SetBinContent(bin, D/P); Ratio[plot]->SetBinError(bin, sqrt( D/(P*P) + pow(D*Perr/(P*P),2) )); } for(int i=0; i<countPlots; i++) { map<std::pair<double, double>,int>::iterator it; double ICut=-1, PtCut=-1; for ( it=CutMap.begin() ; it != CutMap.end(); it++ ) if((*it).second==i) { ICut = (*it).first.first; PtCut = (*it).first.second; } c1 = new TCanvas("c1","c1,",600,600); legend.clear(); Histos[0] = Data[i]; legend.push_back("Observed"); Histos[1] = Pred[i]; legend.push_back("Prediction"); if (TypeMode!=4) { DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta Cut", "Tracks", 0, 0, 1, 100000); } else { DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta Cut", "Tracks", 0, 0, 1, 2200000); } DrawLegend((TObject**)Histos,legend,LegendTitle,"P", 0.5); c1->SetLogy(true); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); char Title[1024]; if(ICut>-1 && PtCut>-1 && PtCut>=100) sprintf(Title,"Pred%s_I%0.2f_Pt%3.0f_",HistoSuffix.c_str(), ICut, PtCut); else if(ICut>-1 && PtCut>-1) sprintf(Title,"Pred%s_I%0.2f_Pt%2.0f_",HistoSuffix.c_str(), ICut, PtCut); else if(PtCut>-1 && PtCut>=100) sprintf(Title,"Pred%s_Pt%3.0f_",HistoSuffix.c_str(), PtCut); else if(PtCut>-1 && PtCut>=10) sprintf(Title,"Pred%s_Pt%2.0f_",HistoSuffix.c_str(), PtCut); else if(ICut>-1) sprintf(Title,"Pred%s_I%0.2f_",HistoSuffix.c_str(),ICut); SaveCanvas(c1,SavePath,Title + DataType); delete c1; delete Histos[0]; delete Histos[1]; } legend.clear(); for(int i=0; i<countPlots; i++) { map<std::pair<double, double>,int>::iterator it; double ICut=-1, PtCut=-1; for ( it=CutMap.begin() ; it != CutMap.end(); it++ ) if((*it).second==i) { ICut = (*it).first.first; PtCut = (*it).first.second; } char LegendName[1024]; if(ICut>-1 && PtCut>-1) sprintf(LegendName,"I>%0.2f Pt>%3.0f",ICut, PtCut); else if(PtCut>-1 && PtCut>=100) sprintf(LegendName,"Pt>%3.0f",PtCut); else if(PtCut>-1 && PtCut>=10) sprintf(LegendName,"Pt>%2.0f",PtCut); else if(ICut>-1) sprintf(LegendName,"I>%0.2f",ICut); Histos[i] = Ratio[i]; legend.push_back(LegendName); } c1 = new TCanvas("c1","c1,",600,600); DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta Cut", "Data/MC", 0, 0, 0,0); DrawLegend((TObject**)Histos,legend,LegendTitle,"P"); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,"Pred_Ratio_" + DataType + HistoSuffix); delete c1; } void CollisionBackgroundSystematicFromFlip(string InputPattern, string DataType){ if(DataType.find("7TeV")!=string::npos){SQRTS=7.0;}else{SQRTS=8.0;} TypeMode = TypeFromPattern(InputPattern); string SavePath = InputPattern; MakeDirectories(SavePath); TCanvas* c1; std::vector<string> legend; TGraphErrors** Graphs = new TGraphErrors*[10]; string LegendTitle = LegendFromType(InputPattern); TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); TH1D* HCuts_Pt = (TH1D*)GetObjectFromPath(InputFile, "HCuts_Pt"); TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, "HCuts_TOF"); TH1D* HCuts_I = (TH1D*)GetObjectFromPath(InputFile, "HCuts_I"); TH1D* HCuts_Pt_Flip = (TH1D*)GetObjectFromPath(InputFile, "HCuts_Pt_Flip"); TH1D* HCuts_TOF_Flip = (TH1D*)GetObjectFromPath(InputFile, "HCuts_TOF_Flip"); TH1D* HCuts_I_Flip = (TH1D*)GetObjectFromPath(InputFile, "HCuts_I_Flip"); TH1D* H_A = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_A"); //TH1D* H_B = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_B"); TH1D* H_C = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_C"); //TH1D* H_D = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_D"); TH1D* H_E = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_E"); TH1D* H_F = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_F"); TH1D* H_G = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_G"); TH1D* H_H = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_H"); TH1D* H_A_Flip = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_A_Flip"); TH1D* H_B_Flip = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_B_Flip"); TH1D* H_C_Flip = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_C_Flip"); TH1D* H_D_Flip = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_D_Flip"); TH1D* H_E_Flip = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_E_Flip"); TH1D* H_F_Flip = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_F_Flip"); TH1D* H_G_Flip = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_G_Flip"); TH1D* H_H_Flip = (TH1D*)GetObjectFromPath(InputFile, DataType + "/H_H_Flip"); TH1D* H_B_Binned[MaxPredBins]; TH1D* H_F_Binned[MaxPredBins]; TH1D* H_H_Binned[MaxPredBins]; TH1D* H_B_Flip_Binned[MaxPredBins]; TH1D* H_D_Flip_Binned[MaxPredBins]; TH1D* H_F_Flip_Binned[MaxPredBins]; TH1D* H_H_Flip_Binned[MaxPredBins]; if(TypeMode==3) { for(int i=0; i<MaxPredBins; i++) { char Suffix[1024]; sprintf(Suffix,"_%i",i); string Bin=Suffix; H_B_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (DataType + "/H_B_Binned" + Bin).c_str()); H_F_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (DataType + "/H_F_Binned" + Bin).c_str()); H_H_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (DataType + "/H_H_Binned" + Bin).c_str()); H_B_Flip_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (DataType + "/H_B_Binned_Flip" + Bin).c_str()); H_D_Flip_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (DataType + "/H_D_Binned_Flip" + Bin).c_str()); H_F_Flip_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (DataType + "/H_F_Binned_Flip" + Bin).c_str()); H_H_Flip_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (DataType + "/H_H_Binned_Flip" + Bin).c_str()); } } vector<double> InversebetaCuts; if (TypeMode!=4) { InversebetaCuts.push_back(1.10); InversebetaCuts.push_back(1.20); } else { InversebetaCuts.push_back(1.05); InversebetaCuts.push_back(1.15); } std::vector<int> Index; std::vector<int> Index_Flip; std::vector<int> Plot; for(int CutIndex_Flip=0; CutIndex_Flip<HCuts_Pt_Flip->GetNbinsX(); CutIndex_Flip++) { if(fabs(HCuts_TOF_Flip->GetBinContent(CutIndex_Flip+1)-( 2-InversebetaCuts[0]))<0.001 || fabs(HCuts_TOF_Flip->GetBinContent(CutIndex_Flip+1)-( 2-InversebetaCuts[1]))<0.001) { for(int CutIndex=0; CutIndex<HCuts_Pt->GetNbinsX(); CutIndex++) { if(fabs(HCuts_TOF_Flip->GetBinContent(CutIndex_Flip+1)-(2-HCuts_TOF->GetBinContent(CutIndex+1)))<0.0001 && fabs(HCuts_Pt_Flip->GetBinContent(CutIndex_Flip+1)-HCuts_Pt->GetBinContent(CutIndex+1))<0.0001 && fabs(HCuts_I_Flip->GetBinContent(CutIndex_Flip+1)-HCuts_I->GetBinContent(CutIndex+1))<0.0001) { if ( (TypeMode == 4) && (fabs(HCuts_I->GetBinContent(CutIndex+1)) > 0.36)) { continue; } else { Index.push_back(CutIndex); Index_Flip.push_back(CutIndex_Flip); if(fabs(HCuts_TOF_Flip->GetBinContent(CutIndex_Flip+1)-( 2-InversebetaCuts[0]))<0.001) Plot.push_back(0); if(fabs(HCuts_TOF_Flip->GetBinContent(CutIndex_Flip+1)-( 2-InversebetaCuts[1]))<0.001) Plot.push_back(1); } } } } } const int TimeRegions=2; int NCuts = Index.size()/TimeRegions; double Pred[TimeRegions][3][NCuts]; double PredErr[TimeRegions][3][NCuts]; double StatSyst[TimeRegions][NCuts]; double Stat[TimeRegions][NCuts]; double Syst[TimeRegions][NCuts]; double PtCut[TimeRegions][NCuts]; double IasCut[TimeRegions][NCuts]; double PredN[TimeRegions]={0}; string Preds[TimeRegions] = {"110", "120"}; if (TypeMode == 4) { Preds[0] = "105"; Preds[1] = "115"; } string PredsLegend[TimeRegions] = {"1/#beta>1.1", "1/#beta>1.2"}; if (TypeMode == 4) { PredsLegend[0] = "1/#beta>1.05"; PredsLegend[1] = "1/#beta>1.15"; } string LegendNames[3]={"BH/F", "BH'/F'", "BD'/B'"}; if (TypeMode == 4) { LegendNames[0] = "CH/G"; LegendNames[1] = "CH'/G'"; LegendNames[2] = "CD'/C'"; } string outfile = SavePath + "BkgUncertainty_" + DataType.substr(4) +".txt"; ofstream fout(outfile.c_str()); if ( ! fout.good() ) { cout << "unable to create file " << outfile << endl; return; } char record[400]; sprintf(record, " This is %5s data", DataType.substr(4).c_str()); fout << record << endl << endl; sprintf(record, " %-10s%-10s%-14s%-14s%-14s%-14s%-14s%-14s%-14s%-14s%-14s", "Ias", "1/beta", "Pred1", "Pred2", "Pred3", "DeltaPred1", "DeltaPred2", "DeltaPred3", "STAT", "STATSYST", "SYST"); fout << record << endl ; for(unsigned int i=0; i<Index.size(); i++) { int CutIndex=Index[i]; int Point=PredN[Plot[i]]; double Sigma = 0; double Mean = 0; int N=0; double A = H_A->GetBinContent(Index[i]); //double B = H_B->GetBinContent(Index[i]); double C = H_C->GetBinContent(Index[i]); //double D = H_D->GetBinContent(Index[i]); double E = H_E->GetBinContent(Index[i]); double F = H_F->GetBinContent(Index[i]); double G = H_G->GetBinContent(Index[i]); double H = H_H->GetBinContent(Index[i]); double A_Flip = H_A_Flip->GetBinContent(Index_Flip[i]); double B_Flip = H_B_Flip->GetBinContent(Index_Flip[i]); double C_Flip = H_C_Flip->GetBinContent(Index_Flip[i]); double D_Flip = H_D_Flip->GetBinContent(Index_Flip[i]); double E_Flip = H_E_Flip->GetBinContent(Index_Flip[i]); double F_Flip = H_F_Flip->GetBinContent(Index_Flip[i]); double G_Flip = H_G_Flip->GetBinContent(Index_Flip[i]); double H_Flip = H_H_Flip->GetBinContent(Index_Flip[i]); double NPred[3]={0}; double NPredErr[3]={0}; //Calculate the three predictions and their statistical errors //The statisitical uncertainty due to the shared region is not included as it is correlated between the predictions if(TypeMode==2) { NPred[0] = (A*F*G)/(E*E); NPredErr[0] = sqrt( ((pow(F*G,2)* A + pow(A*G,2)*F + pow(A*F,2)*G)/pow(E,4)) + (pow((2*A*F*G)/pow(E,3),2)*E)); NPred[1] = (A*F_Flip*G_Flip)/(E_Flip*E_Flip); NPredErr[1] = sqrt( ((pow(A_Flip*G_Flip,2)* F_Flip + pow(A*F_Flip,2)*G_Flip)/pow(E_Flip,4)) + (pow((2*A_Flip*F*G_Flip)/pow(E_Flip,3),2)*E_Flip)); NPred[2] = (A*B_Flip*C_Flip)/(A_Flip*A_Flip); NPredErr[2] = sqrt( ((pow(A*C_Flip,2)* B_Flip + pow(A*B_Flip,2)*C_Flip)/pow(A_Flip,4)) + (pow((2*A*B_Flip*C_Flip)/pow(A_Flip,3),2)*A_Flip)); } if(TypeMode==3) { for(int j=0; j<MaxPredBins; j++) { //double A_Binned = H_A_Binned[j]->GetBinContent(Index[i]); double B_Binned = H_B_Binned[j]->GetBinContent(Index[i]); //double C_Binned = H_C_Binned[j]->GetBinContent(Index[i]); //double D_Binned = H_D_Binned[j]->GetBinContent(Index[i]); //double E_Binned = H_E_Binned[j]->GetBinContent(Index[i]); double F_Binned = H_F_Binned[j]->GetBinContent(Index[i]); //double G_Binned = H_G_Binned[j]->GetBinContent(Index[i]); double H_Binned = H_H_Binned[j]->GetBinContent(Index[i]); //double A_Flip_Binned = H_A_Flip_Binned[j]->GetBinContent(Index_Flip[i]); double B_Flip_Binned = H_B_Flip_Binned[j]->GetBinContent(Index_Flip[i]); //double C_Flip_Binned = H_C_Flip_Binned[j]->GetBinContent(Index_Flip[i]); double D_Flip_Binned = H_D_Flip_Binned[j]->GetBinContent(Index_Flip[i]); //double E_Flip_Binned = H_E_Flip_Binned[j]->GetBinContent(Index_Flip[i]); double F_Flip_Binned = H_F_Flip_Binned[j]->GetBinContent(Index_Flip[i]); //double G_Flip_Binned = H_G_Flip_Binned[j]->GetBinContent(Index_Flip[i]); double H_Flip_Binned = H_H_Flip_Binned[j]->GetBinContent(Index_Flip[i]); NPred[0]+=H_Binned*B_Binned/F_Binned; NPredErr[0] += (pow(B_Binned/F_Binned,2)*H_Binned) + (pow((B_Binned*(H_Binned)/(F_Binned*F_Binned)),2)*F_Binned); NPred[1]+=H_Flip_Binned*B_Binned/F_Flip_Binned; NPredErr[1] += (pow(B_Binned/F_Flip_Binned,2)*H_Flip_Binned) + (pow((B_Binned*(H_Flip_Binned)/(F_Flip_Binned*F_Flip_Binned)),2)*F_Flip_Binned); NPred[2]+=D_Flip_Binned*B_Binned/B_Flip_Binned; NPredErr[2] += (pow(B_Binned/B_Flip_Binned,2)*D_Flip_Binned) + (pow((B_Binned*(D_Flip_Binned)/(B_Flip_Binned*B_Flip_Binned)),2)*B_Flip_Binned); } } if(TypeMode==4) { NPred[0] = ((C*H)/G); NPredErr[0] = (pow(C/G,2)*H) + (pow((H*(C)/(G*G)),2)*G); NPred[1] = ((C*H_Flip)/G_Flip); NPredErr[1] = (pow(C/G_Flip,2)*H_Flip) + (pow((H_Flip*(C)/(G_Flip*G_Flip)),2)*G_Flip); NPred[2] = ((C*D_Flip)/C_Flip); NPredErr[2] = (pow(C/C_Flip,2)*D_Flip) + (pow((D_Flip*(C)/(C_Flip*C_Flip)),2)*C_Flip); } for(int Region=0; Region<3; Region++) { NPredErr[Region]=sqrt(NPredErr[Region]); Pred[Plot[i]][Region][Point] = NPred[Region]; PredErr[Plot[i]][Region][Point] = NPredErr[Region]; Mean+=NPred[Region]/pow(NPredErr[Region],2); Sigma+=1/pow(NPredErr[Region],2); N++; } Mean = Mean/Sigma; Sigma = sqrt(Sigma); double SUM=0, STAT=0, SYST=0; for(int p=0;p<3;p++){ SUM += pow(Pred[Plot[i]][p][Point]-Mean,2); STAT += pow(PredErr[Plot[i]][p][Point],2); } SUM = sqrt(SUM/(N-1)); STAT = sqrt(STAT/(N)); //Use N here as averaging the various statistical uncertaties, correlation due to the shared is accounted for above. if(SUM*SUM > STAT*STAT) SYST = sqrt(SUM*SUM - STAT*STAT); else SYST=0; Stat[Plot[i]][Point]=STAT/Mean; StatSyst[Plot[i]][Point]=SUM/Mean; Syst[Plot[i]][Point]=SYST/Mean; PtCut[Plot[i]][Point]=HCuts_Pt->GetBinContent(CutIndex); IasCut[Plot[i]][Point]=HCuts_I->GetBinContent(CutIndex); PredN[Plot[i]]++; sprintf(record, "%-5i%-10.3f%-10.3f%-14.3f%-14.3f%-14.3f%-14.3f%-14.3f%-14.3f%-14.3f%-14.3f%-14.3f", i, HCuts_I->GetBinContent(CutIndex+1), HCuts_TOF->GetBinContent(CutIndex+1), Pred[Plot[i]][0][Point], Pred[Plot[i]][1][Point], Pred[Plot[i]][2][Point], PredErr[Plot[i]][0][Point], PredErr[Plot[i]][1][Point], PredErr[Plot[i]][2][Point],Stat[Plot[i]][Point], StatSyst[Plot[i]][Point], Syst[Plot[i]][Point]); fout << record << endl ; } fout.close(); TMultiGraph* PredGraphs; for(int i=0; i<TimeRegions; i++) { c1 = new TCanvas("c1","c1,",600,600); legend.clear(); PredGraphs = new TMultiGraph(); for(int Region=0; Region<3; Region++) { if (TypeMode!=4) { Graphs[Region] = new TGraphErrors(PredN[Plot[i]],PtCut[Plot[i]],Pred[Plot[i]][Region],0,PredErr[Plot[i]][Region]); legend.push_back(LegendNames[Region]); } else { std::cout << PredN[Plot[i]] << " " << *IasCut[Plot[i]] << " " << *Pred[Plot[i]][Region] << " " << *PredErr[Plot[i]][Region] << std::endl; Graphs[Region] = new TGraphErrors(PredN[Plot[1]],IasCut[Plot[i]],Pred[Plot[i]][Region],0,PredErr[Plot[i]][Region]); legend.push_back(LegendNames[Region]); } Graphs[Region]->SetLineColor(Color[Region]); Graphs[Region]->SetMarkerColor(Color[Region]); Graphs[Region]->SetMarkerStyle(GraphStyle[Region]); PredGraphs->Add(Graphs[Region],"LP"); } PredGraphs->Draw("A"); PredGraphs->SetTitle(""); PredGraphs->GetXaxis()->SetTitle("P_T cut"); if (TypeMode==4) PredGraphs->GetXaxis()->SetTitle("I_{as} cut"); PredGraphs->GetYaxis()->SetTitle("Number of expected backgrounds"); PredGraphs->GetYaxis()->SetTitleOffset(1.70); //if(i==0) PredGraphs->GetYaxis()->SetRangeUser(0,50000); //else PredGraphs->GetYaxis()->SetRangeUser(0,2400); //if (TypeMode == 4) //{ double yup = Pred[i][0][0]; double ydown = 1000000.0000; for(int Region=0; Region<3; Region++) { double Predmin = Pred[i][Region][NCuts-1]; if (Predmin < ydown) ydown = Predmin; } PredGraphs->GetYaxis()->SetRangeUser(ydown*0.4, yup*1.4); c1->SetLogy(true); //} DrawLegend((TObject**)Graphs,legend,LegendTitle,"P",0.8, 0.9, 0.4, 0.05); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,(DataType + "CollisionPrediction_TOF" + Preds[i]).c_str()); delete c1; delete PredGraphs; } c1 = new TCanvas("c1","c1,",600,600); legend.clear(); PredGraphs = new TMultiGraph(); for(int i=0; i<TimeRegions; i++) { if (TypeMode!=4) { Graphs[i] = new TGraphErrors(PredN[Plot[i]],PtCut[Plot[i]],StatSyst[Plot[i]],0,0); legend.push_back(PredsLegend[i]); } else { Graphs[i] = new TGraphErrors(PredN[Plot[i]],IasCut[Plot[i]],StatSyst[Plot[i]],0,0); legend.push_back(PredsLegend[i]); } Graphs[i]->SetLineColor(Color[i]); Graphs[i]->SetMarkerColor(Color[i]); Graphs[i]->SetMarkerStyle(GraphStyle[i]); PredGraphs->Add(Graphs[i],"LP"); } PredGraphs->Draw("A"); PredGraphs->SetTitle(""); PredGraphs->GetXaxis()->SetTitle("P_T cut"); if (TypeMode==4) PredGraphs->GetXaxis()->SetTitle("I_{as} cut"); PredGraphs->GetYaxis()->SetTitle("Relative Statistical + Systematic Uncertainty"); PredGraphs->GetYaxis()->SetTitleOffset(1.70); PredGraphs->GetYaxis()->SetRangeUser(0,0.40); c1->SetLogy(0); DrawLegend((TObject**)Graphs,legend,LegendTitle,"P",0.8, 0.9, 0.4, 0.05); PredGraphs->GetYaxis()->SetRangeUser(0, 0.4); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,DataType + "CollisionStatSyst"); delete c1; delete PredGraphs; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); PredGraphs = new TMultiGraph(); for(int i=0; i<TimeRegions; i++) { if (TypeMode!=4) { Graphs[i] = new TGraphErrors(PredN[Plot[i]],PtCut[Plot[i]],Stat[Plot[i]],0,0); legend.push_back(PredsLegend[i]); } else { Graphs[i] = new TGraphErrors(PredN[Plot[i]],IasCut[Plot[i]],Stat[Plot[i]],0,0); legend.push_back(PredsLegend[i]); } Graphs[i]->SetLineColor(Color[i]); Graphs[i]->SetMarkerColor(Color[i]); Graphs[i]->SetMarkerStyle(GraphStyle[i]); PredGraphs->Add(Graphs[i],"LP"); } PredGraphs->Draw("A"); PredGraphs->SetTitle(""); PredGraphs->GetXaxis()->SetTitle("P_T cut"); if (TypeMode==4) PredGraphs->GetXaxis()->SetTitle("I_{as} cut"); PredGraphs->GetYaxis()->SetTitle("Relative Statistical Uncertainty"); PredGraphs->GetYaxis()->SetTitleOffset(1.70); PredGraphs->GetYaxis()->SetRangeUser(0,0.40); c1->SetLogy(0); DrawLegend((TObject**)Graphs,legend,LegendTitle,"P",0.8, 0.9, 0.4, 0.05); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,DataType + "CollisionStat"); delete c1; delete PredGraphs; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); PredGraphs = new TMultiGraph(); for(int i=0; i<TimeRegions; i++) { if (TypeMode!=4) { Graphs[i] = new TGraphErrors(PredN[Plot[i]],PtCut[Plot[i]],Syst[Plot[i]],0,0); legend.push_back(PredsLegend[i]); } else { Graphs[i] = new TGraphErrors(PredN[Plot[i]],IasCut[Plot[i]],Syst[Plot[i]],0,0); legend.push_back(PredsLegend[i]); } Graphs[i]->SetLineColor(Color[i]); Graphs[i]->SetMarkerColor(Color[i]); Graphs[i]->SetMarkerStyle(GraphStyle[i]); PredGraphs->Add(Graphs[i],"LP"); } PredGraphs->Draw("A"); PredGraphs->SetTitle(""); PredGraphs->GetXaxis()->SetTitle("P_T cut"); if (TypeMode==4) PredGraphs->GetXaxis()->SetTitle("I_{as} cut"); PredGraphs->GetYaxis()->SetTitle("Relative Systemtic Uncertainty"); PredGraphs->GetYaxis()->SetTitleOffset(1.70); PredGraphs->GetYaxis()->SetRangeUser(0,0.40); c1->SetLogy(0); DrawLegend((TObject**)Graphs,legend,LegendTitle,"P",0.8, 0.9, 0.4, 0.05); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,DataType + "CollisionSyst"); delete c1; delete PredGraphs; /* c1 = new TCanvas("c1","c1,",600,600); legend.clear(); for(int i=0; i<TimeRegions; i++) { Histos[i] = Stat[i]; legend.push_back(PredsLegend[i]); } DrawSuperposedHistos((TH1**)Histos, legend, "", "Pt Cut", "Stat Rel. Error", 0,0, 0,1.4); DrawLegend((TObject**)Histos,legend,LegendTitle,"P"); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,DataType + "CollisionStat"); delete c1; c1 = new TCanvas("c1","c1,",600,600); legend.clear(); for(int i=0; i<TimeRegions; i++) { Histos[i] = Syst[i]; legend.push_back(PredsLegend[i]); } DrawSuperposedHistos((TH1**)Histos, legend, "", "Pt Cut", "Syst Rel. Error", 0,0, 0,1.4); DrawLegend((TObject**)Histos,legend,LegendTitle,"P"); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,DataType +"CollisionSyst"); delete c1; */ } void CheckPredictionBin(string InputPattern, string HistoSuffix, string DataType, string bin){ TypeMode = TypeFromPattern(InputPattern); if(TypeMode==0)return; std::vector<string> legend; string LegendTitle = LegendFromType(InputPattern); string SavePath = InputPattern; MakeDirectories(SavePath); TCanvas* c1; TH1** Histos = new TH1*[100]; TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); TH1D* HCuts_Pt = (TH1D*)GetObjectFromPath(InputFile, string("HCuts_Pt") + HistoSuffix); TH1D* HCuts_I = (TH1D*)GetObjectFromPath(InputFile, string("HCuts_I") + HistoSuffix); TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, string("HCuts_TOF") + HistoSuffix); TH1D* H_D = (TH1D*)GetObjectFromPath(InputFile, string(DataType+"/H_D_Binned" + HistoSuffix + "_" + bin)); TH1D* H_P = (TH1D*)GetObjectFromPath(InputFile, string(DataType+"/H_P_Binned" + HistoSuffix + "_" + bin)); std::vector<int> Index; std::vector<int> Plot; std::vector<double> TOFCuts; double TOFCutMax=0, TOFCutMin=9999; map<std::pair<double, double>,int> CutMap; int countPlots=0; for(int CutIndex=1; CutIndex<HCuts_TOF->GetNbinsX(); CutIndex++) { TOFCuts.push_back(HCuts_TOF->GetBinContent(CutIndex+1)); if(HCuts_TOF->GetBinContent(CutIndex+1)<TOFCutMin) TOFCutMin=HCuts_TOF->GetBinContent(CutIndex+1); if(HCuts_TOF->GetBinContent(CutIndex+1)>TOFCutMax) TOFCutMax=HCuts_TOF->GetBinContent(CutIndex+1); std::pair<double, double> key(HCuts_I->GetBinContent(CutIndex+1), HCuts_Pt->GetBinContent(CutIndex+1)); //New combination of TOF and I cuts if(CutMap.find(key)==CutMap.end()) { CutMap[key]=countPlots; countPlots++; } } std::vector<TH1D*> Pred; std::vector<TH1D*> Data; std::vector<TH1D*> Ratio; for(int i=0; i<countPlots; i++) { char DataName[1024]; sprintf(DataName,"Data_%i",i); TH1D* TempData = new TH1D(DataName, DataName, TOFCuts.size(), TOFCutMin, TOFCutMax); Data.push_back(TempData); char PredName[1024]; sprintf(PredName,"Pred_%i",i); TH1D* TempPred = new TH1D(PredName, PredName, TOFCuts.size(), TOFCutMin, TOFCutMax); Pred.push_back(TempPred); char RatioName[1024]; sprintf(RatioName,"Ratio_%i",i); TH1D* TempRatio = new TH1D(RatioName, RatioName, TOFCuts.size(), TOFCutMin, TOFCutMax); Ratio.push_back(TempRatio); } for(int CutIndex=1; CutIndex<HCuts_TOF->GetNbinsX(); CutIndex++) { std::pair<double, double> key(HCuts_I->GetBinContent(CutIndex+1), HCuts_Pt->GetBinContent(CutIndex+1)); int plot = CutMap.find(key)->second; int histo_bin = Data[plot]->FindBin(HCuts_TOF->GetBinContent(CutIndex+1)); double D = H_D->GetBinContent(CutIndex+1); Data[plot]->SetBinContent(histo_bin, D); double P = H_P->GetBinContent(CutIndex+1); double Perr = H_P->GetBinError(CutIndex+1); Pred[plot]->SetBinContent(histo_bin, P); Pred[plot]->SetBinError(histo_bin, Perr); Ratio[plot]->SetBinContent(histo_bin, D/P); Ratio[plot]->SetBinError(histo_bin, sqrt( D/(P*P) + pow(D*Perr/(P*P),2) )); } for(int i=0; i<countPlots; i++) { map<std::pair<double, double>,int>::iterator it; double ICut=-1, PtCut=-1; for ( it=CutMap.begin() ; it != CutMap.end(); it++ ) if((*it).second==i) { ICut = (*it).first.first; PtCut = (*it).first.second; } c1 = new TCanvas("c1","c1,",600,600); legend.clear(); Histos[0] = Data[i]; legend.push_back("Obs"); Histos[1] = Pred[i]; legend.push_back("Pred"); DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta Cut", "Tracks", 0, 0, 0,0); DrawLegend((TObject**)Histos,legend,LegendTitle,"P", 0.5); c1->SetLogy(true); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); char Title[1024]; if(ICut>-1 && PtCut>-1) sprintf(Title,"Pred%s_I%0.2f_Pt%3.0f_",HistoSuffix.c_str(), ICut, PtCut); else if(PtCut>-1) sprintf(Title,"Pred%s_Pt%3.0f_",HistoSuffix.c_str(),PtCut); else if(ICut>-1) sprintf(Title,"Pred%s_I%0.2f_",HistoSuffix.c_str(),ICut); SaveCanvas(c1,SavePath,Title + DataType + "_Binned_" + bin); delete c1; delete Histos[0]; delete Histos[1]; } legend.clear(); for(int i=0; i<countPlots; i++) { map<std::pair<double, double>,int>::iterator it; double ICut=-1, PtCut=-1; for ( it=CutMap.begin() ; it != CutMap.end(); it++ ) if((*it).second==i) { ICut = (*it).first.first; PtCut = (*it).first.second; } char LegendName[1024]; if(ICut>-1 && PtCut>-1) sprintf(LegendName,"I>%0.2f Pt>%3.0f",ICut, PtCut); else if(PtCut>-1) sprintf(LegendName,"Pt>%3.0f",PtCut); else if(ICut>-1) sprintf(LegendName,"I>%0.2f",ICut); Histos[i] = Ratio[i]; legend.push_back(LegendName); } c1 = new TCanvas("c1","c1,",600,600); DrawSuperposedHistos((TH1**)Histos, legend, "E1", "1/#beta Cut", "Data/MC", 0, 0, 0,0); DrawLegend((TObject**)Histos,legend,LegendTitle,"P"); c1->SetLogy(false); DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); SaveCanvas(c1,SavePath,"Pred_Ratio_" + DataType + HistoSuffix); delete c1; } void Make2DPlot_Special(string InputPattern, string InputPattern2){//, unsigned int CutIndex){ TCanvas* c1; TLegend* leg; string Input = InputPattern + "Histos.root"; string outpath = InputPattern; MakeDirectories(outpath); TypeMode = TypeFromPattern(InputPattern); string LegendTitle = LegendFromType(InputPattern);; string S1 = "DY_8TeV_M400_Q2o3"; //double Q1=1; string S2 = "DY_8TeV_M400_Q1"; //double Q2=1; string S3 = "DY_8TeV_M400_Q3"; //double Q3=1; string Da = "Data8TeV"; string outName = "2DPlotsS"; int S1i = JobIdToIndex(S1,samples); if(S1i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } int S2i = JobIdToIndex(S2,samples); if(S2i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } int S3i = JobIdToIndex(S3,samples); if(S3i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); // TH2D* Signal1PIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S1+"/AS_PIm" ), CutIndex, "S1PIm_zy" ); // TH2D* Signal2PIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S2+"/AS_PIm" ), CutIndex, "S2PIm_zy" ); // TH2D* Signal3PIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, S3+"/AS_PIm" ), CutIndex, "S3PIm_zy" ); // TH2D* Data_PIm = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile, Da+"/AS_PIm" ), CutIndex); TH2D* Signal1PIm = (TH2D*)GetObjectFromPath(InputFile, S1+"/BS_PImHD" ); TH2D* Signal2PIm = (TH2D*)GetObjectFromPath(InputFile, S2+"/BS_PImHD" ); TH2D* Signal3PIm = (TH2D*)GetObjectFromPath(InputFile, S3+"/BS_PImHD" ); TH2D* Data_PIm = (TH2D*)GetObjectFromPath(InputFile, Da+"/BS_PImHD" ); TFile* InputFile2 = new TFile((InputPattern2 + "Histos.root").c_str()); // TH2D* Signal1PIm2 = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile2, S1+"/AS_PIm" ), CutIndex, "S1PIm_zy" ); //// TH2D* Signal2PIm2 = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile2, S2+"/AS_PIm" ), CutIndex, "S2PIm_zy" ); //// TH2D* Signal3PIm2 = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile2, S3+"/AS_PIm" ), CutIndex, "S3PIm_zy" ); // TH2D* Data_PIm2 = GetCutIndexSliceFromTH3((TH3D*)GetObjectFromPath(InputFile2, Da+"/AS_PIm" ), CutIndex); TH2D* Signal1PIm2 = (TH2D*)GetObjectFromPath(InputFile2, S1+"/BS_PImHD"); TH2D* Data_PIm2 = (TH2D*)GetObjectFromPath(InputFile2, Da+"/BS_PImHD"); Signal1PIm->Add(Signal1PIm2); // Signal2PIm->Add(Signal2PIm2); // Signal3PIm->Add(Signal3PIm2); Data_PIm ->Add(Data_PIm2); gStyle->SetPalette(53); c1 = new TCanvas("c1","c1", 600, 600); c1->SetLeftMargin (0.14); c1->SetRightMargin (0.20); c1->SetLogz(true); Data_PIm->SetTitle(""); Data_PIm->SetStats(kFALSE); Data_PIm->GetXaxis()->SetTitle("p (GeV/#font[12]{c})"); Data_PIm->GetYaxis()->SetTitle(dEdxM_Legend.c_str()); Data_PIm->GetYaxis()->SetTitleOffset(1.40); Data_PIm->GetZaxis()->SetTitleOffset(1.60); Data_PIm->GetZaxis()->SetTitle("Tracks / [ 2.4 (GeV/#font[12]{c}) #times 0.03 (MeV/cm) ]"); Data_PIm->GetZaxis()->CenterTitle(true); printf("%f - %f\n", Data_PIm->GetXaxis()->GetBinWidth(3), Data_PIm->GetYaxis()->GetBinWidth(3) ); Data_PIm->SetAxisRange(50,1750,"X"); Data_PIm->SetAxisRange(0,20,"Y"); Data_PIm->SetMarkerSize (0.2); Data_PIm->SetMarkerColor(Color[1]); Data_PIm->SetFillColor(Color[1]); Data_PIm->Draw("COLZ"); // DrawPreliminary(LegendTitle, SQRTS, IntegratedLuminosityFromE(SQRTS)); // SaveCanvas(c1, outpath, outName + "_Data_PIm", true); // delete c1; // c1 = new TCanvas("c1","c1", 600, 600); // c1->SetLogz(true); Signal3PIm->SetTitle(""); Signal3PIm->SetStats(kFALSE); Signal3PIm->GetXaxis()->SetTitle("p (GeV/#font[12]{c})"); Signal3PIm->GetYaxis()->SetTitle(dEdxM_Legend.c_str()); Signal3PIm->SetAxisRange(50,1750,"X"); Signal3PIm->GetYaxis()->SetTitleOffset(1.40); Signal3PIm->GetZaxis()->SetTitleOffset(1.60); Signal3PIm->SetAxisRange(50,1750,"X"); Signal3PIm->SetAxisRange(0,20,"Y"); // Signal3PIm->SetAxisRange(0,TypeMode==5?3:15,"Y"); Signal3PIm->Scale(5000/Signal3PIm->Integral()); Signal3PIm->SetMarkerStyle (1); Signal3PIm->SetMarkerColor(Color[4]); Signal3PIm->SetFillColor(Color[4]); Signal3PIm->Draw("SCAT same"); Signal2PIm->Scale(5000/Signal2PIm->Integral()); Signal2PIm->SetMarkerStyle (1); Signal2PIm->SetMarkerColor(Color[3]); Signal2PIm->SetFillColor(Color[3]); Signal2PIm->Draw("SCAT same"); Signal1PIm->Scale(5000/Signal1PIm->Integral()); Signal1PIm->SetMarkerStyle(1); Signal1PIm->SetMarkerColor(Color[2]); Signal1PIm->SetFillColor(Color[2]); Signal1PIm->Draw("SCAT same"); // TF1* MassLine800 = GetMassLineQ(samples[S3i].Mass, Q3, true); // MassLine800->SetLineColor(kGray+3); // MassLine800->SetLineWidth(2); // MassLine800->Draw("same"); // TF1* MassLine500 = GetMassLineQ(samples[S2i].Mass, Q2, true); // MassLine500->SetLineColor(kBlue-7); // MassLine500->SetLineWidth(2); // MassLine500->Draw("same"); // TF1* MassLine300 = GetMassLineQ(samples[S1i].Mass, Q1, true); // MassLine300->SetLineColor(kRed-7); // MassLine300->SetLineWidth(2); // MassLine300->Draw("same"); TBox* box = new TBox(50.0, 2.8, 1200.0, 3.0); box->SetFillColor(1); box->SetFillStyle(3004); box->Draw("same"); leg = new TLegend(0.80,0.92,0.80 - 0.41,0.92 - 6*0.03); leg->SetTextFont(43); //give the font size in pixel (instead of fraction) leg->SetTextSize(18); //font size leg->SetFillColor(0); leg->SetFillStyle(0); leg->SetBorderSize(0); leg->AddEntry(Data_PIm, "Data (#sqrt{s} = 8 TeV)" ,"F"); leg->AddEntry(Signal3PIm, samples[S3i].Legend.c_str() ,"F"); leg->AddEntry(Signal2PIm, samples[S2i].Legend.c_str() ,"F"); leg->AddEntry(Signal1PIm, samples[S1i].Legend.c_str() ,"F"); leg->AddEntry(box, "Excluded" ,"F"); leg->Draw(); DrawPreliminary("", SQRTS, IntegratedLuminosityFromE(SQRTS), false); SaveCanvas(c1, outpath, outName + "_PIm", false); delete c1; delete leg; c1 = new TCanvas("c1","c1", 600, 600); c1->SetLeftMargin (0.14); c1->SetRightMargin (0.20); c1->SetLogz(true); Data_PIm->Draw("COLZ"); box->Draw("same"); leg = new TLegend(0.80,0.92,0.80 - 0.20,0.92 - 0.04); leg->SetTextFont(43); leg->SetTextSize(18); leg->SetFillColor(0); leg->SetFillStyle(0); leg->SetBorderSize(0); //leg->AddEntry(Data_PIm, "Data (#sqrt{s}=8 TeV)" ,"F"); leg->AddEntry(box, "Excluded" ,"F"); leg->Draw(); DrawPreliminary("", SQRTS, IntegratedLuminosityFromE(SQRTS), false); SaveCanvas(c1, outpath, outName + "_PImData", false); delete c1; delete leg; c1 = new TCanvas("c1","c1", 600, 600); c1->SetLeftMargin (0.14); c1->SetRightMargin (0.16); c1->SetLogz(true); int backPalette[] = {1, 1}; gStyle->SetPalette(2, backPalette); Data_PIm->SetFillColor(1); Data_PIm->Draw("COL"); Signal3PIm->Draw("SCAT same"); Signal2PIm->Draw("SCAT same"); Signal1PIm->Draw("SCAT same"); box->Draw("same"); leg = new TLegend(0.80,0.92,0.80 - 0.41,0.92 - 5*0.04); leg->SetTextFont(43); leg->SetTextSize(18); leg->SetFillColor(0); leg->SetFillStyle(0); leg->SetBorderSize(0); leg->AddEntry(Data_PIm, "Data (#sqrt{s} = 8 TeV)" ,"F"); leg->AddEntry(Signal3PIm, samples[S3i].Legend.c_str() ,"F"); leg->AddEntry(Signal2PIm, samples[S2i].Legend.c_str() ,"F"); leg->AddEntry(Signal1PIm, samples[S1i].Legend.c_str() ,"F"); leg->AddEntry(box, "Excluded" ,"F"); leg->Draw(); DrawPreliminary("", SQRTS, IntegratedLuminosityFromE(SQRTS), false); SaveCanvas(c1, outpath, outName + "_PImSignal", false); delete c1; delete leg; gStyle->SetPalette(1); } void CompareRecoAndGenPt(string InputPattern){ TCanvas* c1; TLegend* leg; string outpath = InputPattern; MakeDirectories(outpath); TypeMode = TypeFromPattern(InputPattern); string LegendTitle = LegendFromType(InputPattern);; string S1 = "DY_7TeV_M400_Q2o3"; //double Q1=1/3.0; string S2 = "DY_7TeV_M400_Q1"; //double Q2=1; string S3 = "DY_7TeV_M400_Q2"; //double Q3=1; int S1i = JobIdToIndex(S1,samples); if(S1i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } int S2i = JobIdToIndex(S2,samples); if(S2i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } int S3i = JobIdToIndex(S3,samples); if(S3i<0){ printf("There is no signal corresponding to the JobId Given\n"); return; } TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); TH2D* Signal1genrecopt = (TH2D*)GetObjectFromPath(InputFile, S1+"/genrecopT"); TH2D* Signal2genrecopt = (TH2D*)GetObjectFromPath(InputFile, S2+"/genrecopT"); TH2D* Signal3genrecopt = (TH2D*)GetObjectFromPath(InputFile, S3+"/genrecopT"); c1 = new TCanvas("c1","c1", 600, 600); Signal1genrecopt->SetTitle(""); Signal1genrecopt->SetStats(kFALSE); Signal1genrecopt->GetXaxis()->SetTitle("gen p_{T} (GeV/c)"); Signal1genrecopt->GetYaxis()->SetTitle("reco p_{T} (GeV/c)"); Signal1genrecopt->GetYaxis()->SetTitleOffset(1.70); Signal1genrecopt->SetAxisRange(0,1200,"X"); Signal1genrecopt->SetAxisRange(0,1200,"Y"); Signal1genrecopt->GetXaxis()->SetNdivisions(206, 0); Signal1genrecopt->GetYaxis()->SetNdivisions(206, 0); Signal1genrecopt->SetMarkerSize (0.2); Signal1genrecopt->SetMarkerColor(Color[2]); Signal1genrecopt->SetFillColor(Color[2]); Signal1genrecopt->Draw("BOX"); Signal2genrecopt->SetMarkerSize (0.2); Signal2genrecopt->SetMarkerColor(Color[1]); Signal2genrecopt->SetFillColor(Color[1]); Signal2genrecopt->Draw("BOX same"); Signal3genrecopt->SetMarkerSize (0.2); Signal3genrecopt->SetMarkerColor(Color[0]); Signal3genrecopt->SetFillColor(Color[0]); Signal3genrecopt->Draw("BOX same"); TLine* diagonal = new TLine(0, 0, 1200, 1200); diagonal->SetLineWidth(2); diagonal->SetLineColor(Color[0]); diagonal->SetLineStyle(2); diagonal->Draw("same"); TLine* diagonalf = new TLine(0, 0, 800, 1200); diagonalf->SetLineWidth(2); diagonalf->SetLineColor(Color[0]); diagonalf->SetLineStyle(2); diagonalf->Draw("same"); TLine* diagonalm = new TLine(0, 0, 1200, 600); diagonalm->SetLineWidth(2); diagonalm->SetLineColor(Color[0]); diagonalm->SetLineStyle(2); diagonalm->Draw("same"); leg = new TLegend(0.80,0.93,0.80 - 0.40,0.93 - 6*0.03); leg->SetFillColor(0); leg->SetBorderSize(0); leg->AddEntry(Signal1genrecopt, samples[S1i].Legend.c_str() ,"F"); leg->AddEntry(Signal2genrecopt, samples[S2i].Legend.c_str() ,"F"); leg->AddEntry(Signal3genrecopt, samples[S3i].Legend.c_str() ,"F"); leg->Draw(); DrawPreliminary("Simulation", SQRTS, ""); SaveCanvas(c1, outpath, "SIM_Validation_Pt", false); delete c1; return; } void CheckPUDistribution(string InputPattern, unsigned int CutIndex){ TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); if(!InputFile)std::cout << "FileProblem\n"; TCanvas* c1 = new TCanvas("c1","c1",600, 600); TH1F* frame = new TH1F("frame", "frame", 1,0,55); frame->GetXaxis()->SetNdivisions(505); frame->SetTitle(""); frame->SetStats(kFALSE); frame->GetXaxis()->SetTitle("#Vertices"); //frame->GetYaxis()->SetTitle(YAxisLegend); frame->GetYaxis()->SetTitleOffset(1.50); frame->SetMaximum(1); frame->SetMinimum(1E-3); frame->SetAxisRange(0,45,"X"); frame->Draw("AXIS"); for(unsigned int s=0;s<samples.size();s++){ if(samples[s].Type!=2)continue; stPlots plots; if(stPlots_InitFromFile(InputFile, plots, samples[s].Name)){ TH1F* plot = (TH1F*)plots.BS_NVertex_NoEventWeight->Clone((samples[s].Name + "_plot").c_str()); plot->Scale(1.0/plot->Integral()); if(samples[s].Pileup=="S10"){ plot->SetLineColor(2); plot->SetMarkerColor(2);} if(samples[s].Pileup=="S3") { plot->SetLineColor(8); plot->SetMarkerColor(8);} if(samples[s].Pileup=="S4") { plot->SetLineColor(4); plot->SetMarkerColor(4);} plot->Draw("E1 same"); printf("%30s --> %5s <--> %f\n", samples[s].Name.c_str(), samples[s].Pileup.c_str(), plot->GetMean()); stPlots_Clear(&plots); } } TLegend* leg = new TLegend(0.79,0.93,0.35,0.75); leg->SetHeader("Pileup Scenario:"); leg->SetFillStyle(0); leg->SetBorderSize(0); TLegendEntry* entry = NULL; entry = leg->AddEntry("", "S4" ,"P"); entry->SetMarkerColor(4); entry = leg->AddEntry("", "S10" ,"P"); entry->SetMarkerColor(2); leg->Draw("same"); c1->SetLogy(true); SaveCanvas(c1, "test/", "pileup"); delete c1; InputFile->Close(); } /* void CheckPUDistribution(string InputPattern, unsigned int CutIndex){ TFile* InputFile = new TFile((InputPattern + "Histos.root").c_str()); if(!InputFile)std::cout << "FileProblem\n"; for(unsigned int s=0;s<samples.size();s++){ if(samples[s].Type!=2)continue; if(samples[s].Name.find("DY")==string::npos || samples[s].Name.find("o3")!=string::npos)continue; stPlots plots; if(stPlots_InitFromFile(InputFile, plots, samples[s].Name)){ TCanvas* c1 = new TCanvas("c1","c1",600, 600); TH2F* plot = (TH2F*)plots.BS_PImHD->Clone((samples[s].Name + "_plot").c_str()); plot->Rebin2D(4,4); plot->Draw("COLZ"); c1->SetLogz(true); SaveCanvas(c1, "test/", samples[s].Name); delete c1; stPlots_Clear(&plots); } } InputFile->Close(); } */