| 79 |
|
datacard << "rate "<<signal->Integral()<<" "<<background1->Integral()<<" "<<background2->Integral()<<"\n"; |
| 80 |
|
datacard << "--------------------------------\n"; |
| 81 |
|
datacard << "lumi lnN " << 1+ PlottingSetup::lumiuncert << " - - luminosity uncertainty\n"; // only affects MC -> signal! |
| 82 |
– |
datacard << "Stat shape 1 - - statistical uncertainty (signal)\n"; |
| 82 |
|
datacard << "Stat shape - 1 - statistical uncertainty (ttbar)\n"; |
| 83 |
|
datacard << "Stat shape - - 1 statistical uncertainty (zjets)\n"; |
| 84 |
< |
datacard << "Sys shape - - 1 systematic uncertainty on ttbar\n"; |
| 85 |
< |
datacard << "Sys shape - 1 - systematic uncertainty on zjets\n"; |
| 84 |
> |
datacard << "Sys shape - 1 - systematic uncertainty on ttbar\n"; |
| 85 |
> |
datacard << "Sys shape - - 1 systematic uncertainty on zjets\n"; |
| 86 |
> |
datacard << "Stat shape 1 - - statistical uncertainty (signal)\n"; |
| 87 |
|
datacard << "JES shape 1 - - uncertainty on Jet Energy Scale\n"; |
| 88 |
|
datacard << "JSU lnN " << 1+uJSU << " - - JZB Scale Uncertainty\n"; |
| 89 |
|
if(uPDF>0) datacard << "PDF lnN " << 1+uPDF << " - - uncertainty from PDFs\n"; |
| 109 |
|
stringstream mSUGRAweight; |
| 110 |
|
float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i); |
| 111 |
|
totalxs+=xs; |
| 112 |
+ |
mSUGRAweight << "(0.5*(2*(process==" << i << "))*" << xs << ")"; |
| 113 |
|
events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str()); |
| 114 |
|
} |
| 115 |
< |
histo->Scale(1.0/totalxs); |
| 115 |
> |
histo->Scale(1.0/totalxs); // to get back to 1 /pb normalization |
| 116 |
|
} |
| 117 |
|
events->Draw("eventNum",addcut.c_str(),"goff"); |
| 118 |
|
float nevents = events->GetSelectedRows(); |
| 148 |
|
} |
| 149 |
|
|
| 150 |
|
void generate_shapes_for_systematic(bool signalonly, bool dataonly,TFile *limfile, TTree *signalevents, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan, string addcut, map < pair<float, float>, map<string, float> > &xsec) { |
| 151 |
+ |
|
| 152 |
|
binning.push_back(8000); |
| 153 |
|
limfile->cd(); |
| 154 |
< |
dout << "Creatig shape templates "; |
| 154 |
> |
dout << "Creating shape templates "; |
| 155 |
|
if(identifier!="") dout << "for "<<identifier; |
| 156 |
|
dout << " ... " << endl; |
| 157 |
|
|
| 158 |
|
TCut limitnJetcut; |
| 159 |
+ |
|
| 160 |
|
if(JES==noJES) limitnJetcut=cutnJets; |
| 161 |
|
else { |
| 162 |
|
if(JES==JESdown) limitnJetcut=cutnJetsJESdown; |
| 163 |
|
if(JES==JESup) limitnJetcut=cutnJetsJESup; |
| 164 |
|
} |
| 165 |
|
|
| 163 |
– |
stringstream mSUGRAweight; |
| 164 |
– |
if(SUSYScanSpace::SUSYscantype==mSUGRA) { |
| 165 |
– |
//for(int i< |
| 166 |
– |
mSUGRAweight << "bla"; |
| 167 |
– |
} else mSUGRAweight << "(1)"; |
| 168 |
– |
|
| 169 |
– |
write_warning(__FUNCTION__,"Not done yet for mSUGRA"); |
| 170 |
– |
|
| 166 |
|
float zjetsestimateuncert=zjetsestimateuncertONPEAK; |
| 167 |
|
float emuncert=emuncertONPEAK; |
| 168 |
|
float emsidebanduncert=emsidebanduncertONPEAK; |
| 215 |
|
cout << " OSSF SB JZB> 0 : " << SBOSSFP->Integral() << " JZB < 0 :" << SBOSSFN->Integral() << endl; |
| 216 |
|
cout << " OSOF SB JZB> 0 : " << SBOSOFP->Integral() << " JZB < 0 :" << SBOSOFN->Integral() << endl; |
| 217 |
|
} |
| 218 |
< |
|
| 218 |
> |
|
| 219 |
|
|
| 220 |
|
flippedLobs->Add(ZOSSFN); |
| 221 |
|
flippedLpred->Add(ZOSSFP); |
| 267 |
|
|
| 268 |
|
for(int i=1;i<=signalStatDn->GetNbinsX();i++) { |
| 269 |
|
float staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i)); |
| 270 |
+ |
cout << "Stat err in bin " << i << " : " << staterr << endl; |
| 271 |
+ |
cout << " prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl; |
| 272 |
+ |
cout << " we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl; |
| 273 |
|
signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr); |
| 274 |
|
signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr); |
| 275 |
|
signal->SetBinError(i,staterr); |
| 665 |
|
|
| 666 |
|
} |
| 667 |
|
|
| 668 |
< |
ShapeDroplet LimitsFromShapes(TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) { |
| 668 |
> |
ShapeDroplet LimitsFromShapes(bool asymptotic, float firstGuess, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) { |
| 669 |
|
map < pair<float, float>, map<string, float> > xsec; |
| 670 |
|
if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile); |
| 671 |
|
|
| 674 |
– |
write_info(__FUNCTION__,"Implementing the shape function"); |
| 672 |
|
time_t timestamp; |
| 673 |
|
tm *now; |
| 674 |
|
timestamp = time(0); |
| 719 |
|
MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1); |
| 720 |
|
if(mceff<0) flipped=1; |
| 721 |
|
doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut); |
| 722 |
< |
float PDFuncert; |
| 722 |
> |
float PDFuncert=0; |
| 723 |
|
int NPdfs = get_npdfs(events); |
| 724 |
|
if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut); |
| 725 |
|
|
| 729 |
– |
|
| 726 |
|
float JZBscale=scaledown; |
| 727 |
|
if(scaleup>scaledown) JZBscale=scaleup; |
| 728 |
|
dout << endl; |
| 731 |
|
dout << " JZB Scale (JSU) : " << JZBscale << endl; |
| 732 |
|
dout << " PDF : " << PDFuncert << endl; |
| 733 |
|
|
| 734 |
< |
prepare_limit_datacard(RunDirectory,limfile,JZBscale,PDFuncert,flipped); |
| 734 |
> |
float SignalIntegral; |
| 735 |
> |
if(flipped) { |
| 736 |
> |
TH1F *flisi = (TH1F*)(limfile->Get("flipped_signal")); |
| 737 |
> |
SignalIntegral= flisi ->Integral(); |
| 738 |
> |
delete flisi; |
| 739 |
> |
} |
| 740 |
> |
else { |
| 741 |
> |
TH1F *flisi = (TH1F*)(limfile->Get("signal")); |
| 742 |
> |
SignalIntegral= flisi ->Integral(); |
| 743 |
> |
delete flisi; |
| 744 |
> |
} |
| 745 |
> |
|
| 746 |
|
TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE"); |
| 747 |
|
|
| 748 |
|
TIter nnextkey(limfile->GetListOfKeys()); |
| 749 |
|
TKey *nkey; |
| 750 |
+ |
|
| 751 |
+ |
if(SignalIntegral==0) SignalIntegral=1; //case when the integral is zero - don't try to normalize, just leave it as it is |
| 752 |
+ |
|
| 753 |
|
while ((nkey = (TKey*)nnextkey())) |
| 754 |
|
{ |
| 755 |
|
TH1F *obj =(TH1F*) nkey->ReadObj(); |
| 757 |
|
if(!flipped&&Contains(obj->GetName(),"flipped")) continue; |
| 758 |
|
if(flipped) obj->SetName(strip_flip_away(obj->GetName())); |
| 759 |
|
final_limfile->cd(); |
| 760 |
+ |
if(Contains(obj->GetName(),"signal")) { |
| 761 |
+ |
obj->Scale(1.0/SignalIntegral); |
| 762 |
+ |
} |
| 763 |
|
obj->Write(); |
| 764 |
|
} |
| 765 |
|
|
| 766 |
+ |
prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert,flipped); |
| 767 |
+ |
|
| 768 |
|
final_limfile->Close(); |
| 769 |
|
limfile->Close(); |
| 770 |
< |
write_info(__FUNCTION__,"Shape root file and datacard have been generated in "+RunDirectory); |
| 771 |
< |
|
| 772 |
< |
int CreatedModelFileExitCode = gSystem->Exec(("bash CreateModel.sh "+RunDirectory+" susydatacard.txt").c_str()); |
| 770 |
> |
dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl; |
| 771 |
> |
stringstream command; |
| 772 |
> |
if(asymptotic) { |
| 773 |
> |
if(firstGuess>0) command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS |
| 774 |
> |
else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS |
| 775 |
> |
} |
| 776 |
> |
else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.5 * firstGuess) << " " << int(2*firstGuess); // ASYMPTOTIC LIMITS |
| 777 |
> |
dout <<"Going to run : " << command.str() << endl; |
| 778 |
> |
int CreatedModelFileExitCode = gSystem->Exec(command.str().c_str()); |
| 779 |
|
cout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl; |
| 780 |
|
assert(CreatedModelFileExitCode==0); |
| 781 |
|
ShapeDroplet alpha; |
| 782 |
|
alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt"); |
| 783 |
|
alpha.PDF=PDFuncert; |
| 784 |
|
alpha.JSU=JZBscale; |
| 785 |
+ |
alpha.SignalIntegral=SignalIntegral; |
| 786 |
|
dout << "Done reading limit information from droplet" << endl; |
| 787 |
|
dout << alpha << endl; |
| 788 |
|
|
| 767 |
– |
write_warning(__FUNCTION__,"Once you're done implementing the limit computation, please uncomment the following line to enable clean up."); |
| 789 |
|
dout << "Everything is saved in " << RunDirectory << endl; |
| 790 |
+ |
dout << "Will transfer model and datacard over for possible post-processing" << endl; |
| 791 |
+ |
dout << " 1) Make sure models directory exists ... " << std::flush; |
| 792 |
+ |
gSystem->Exec("mkdir -p models/"); |
| 793 |
+ |
dout << " ok!" << endl; |
| 794 |
+ |
dout << " 2) Deleting any previous model files with the same name ... " << std::flush; |
| 795 |
+ |
stringstream delcommand; |
| 796 |
+ |
delcommand << "rm models/model_" << name << ".root"; |
| 797 |
+ |
gSystem->Exec(delcommand.str().c_str()); |
| 798 |
+ |
stringstream delcommand2; |
| 799 |
+ |
delcommand2 << "rm models/model_" << name << ".txt"; |
| 800 |
+ |
gSystem->Exec(delcommand2.str().c_str()); |
| 801 |
+ |
dout << " ok!" << endl; |
| 802 |
+ |
dout << " 3) Transfering the new model files ... " << std::flush; |
| 803 |
+ |
stringstream copycommand; |
| 804 |
+ |
copycommand << "cp " << RunDirectory << "/model.root models/model_" << name << ".root"; |
| 805 |
+ |
gSystem->Exec(copycommand.str().c_str()); |
| 806 |
+ |
stringstream copycommand2; |
| 807 |
+ |
copycommand2 << "cp " << RunDirectory << "/susydatacard.txt models/model_" << name << ".txt"; |
| 808 |
+ |
gSystem->Exec(copycommand2.str().c_str()); |
| 809 |
+ |
stringstream copycommand3; |
| 810 |
+ |
copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root"; |
| 811 |
+ |
gSystem->Exec(copycommand3.str().c_str()); |
| 812 |
+ |
dout << " ok!" << endl; |
| 813 |
+ |
dout << " 4) Removing original working directory (" << RunDirectory << ") ... " << std::flush; |
| 814 |
+ |
write_warning(__FUNCTION__,"Watch out : need to uncomment the line below to remove the original working directory again"); |
| 815 |
|
// gSystem->Exec(("rm -r "+RunDirectory).c_str()); |
| 816 |
< |
|
| 817 |
< |
return alpha; |
| 816 |
> |
dout << " ok!" << endl; |
| 817 |
> |
delete limcan; |
| 818 |
> |
return alpha; |
| 819 |
|
} |
| 820 |
|
|
| 821 |
|
void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) { |
