| 1 |
#include <iostream>
|
| 2 |
#include <sstream>
|
| 3 |
|
| 4 |
#include "TFile.h"
|
| 5 |
#include "TTree.h"
|
| 6 |
#include "TH1.h"
|
| 7 |
#include "TH2.h"
|
| 8 |
#include "TMath.h"
|
| 9 |
|
| 10 |
using namespace std;
|
| 11 |
|
| 12 |
int main() {
|
| 13 |
TFile *f = new TFile("/scratch/buchmann/mSUGRA_jzb.root");
|
| 14 |
TTree *events = (TTree*)f->Get("events");
|
| 15 |
|
| 16 |
|
| 17 |
bool beverbose=false;
|
| 18 |
|
| 19 |
|
| 20 |
long nevents=events->GetEntries();
|
| 21 |
// long uppercutoff=nevents;
|
| 22 |
long uppercutoff=nevents;
|
| 23 |
|
| 24 |
cout << "There's a total of " << nevents << " events" << endl;
|
| 25 |
|
| 26 |
float jzb[5];
|
| 27 |
int pfJetGoodNum;
|
| 28 |
float mll;
|
| 29 |
float pt1;
|
| 30 |
float pt2;
|
| 31 |
float pt;
|
| 32 |
int pfJetGoodID[30];
|
| 33 |
|
| 34 |
int id1;
|
| 35 |
int id2;
|
| 36 |
int ch1;
|
| 37 |
int ch2;
|
| 38 |
|
| 39 |
float M0;
|
| 40 |
float M12;
|
| 41 |
|
| 42 |
events->SetBranchAddress("id1",&id1);
|
| 43 |
events->SetBranchAddress("id2",&id2);
|
| 44 |
events->SetBranchAddress("ch1",&ch1);
|
| 45 |
events->SetBranchAddress("ch2",&ch2);
|
| 46 |
events->SetBranchAddress("jzb",&jzb);
|
| 47 |
events->SetBranchAddress("mll",&mll);
|
| 48 |
events->SetBranchAddress("pt1",&pt1);
|
| 49 |
events->SetBranchAddress("pt2",&pt2);
|
| 50 |
events->SetBranchAddress("pt",&pt);
|
| 51 |
events->SetBranchAddress("pfJetGoodNum",&pfJetGoodNum);
|
| 52 |
events->SetBranchAddress("pfJetGoodID",&pfJetGoodID);
|
| 53 |
|
| 54 |
|
| 55 |
events->SetBranchAddress("M0",&M0);
|
| 56 |
events->SetBranchAddress("M12",&M12);
|
| 57 |
|
| 58 |
|
| 59 |
|
| 60 |
|
| 61 |
vector<string> lepcut;
|
| 62 |
//if you adapt these you also need to adapt the definition of the regions below! (e.g. "if((pt1>10&&pt1<15)||(pt2>10&&pt2<15)) categories[0]=1;")
|
| 63 |
lepcut.push_back("10to15");
|
| 64 |
lepcut.push_back("15to20");
|
| 65 |
lepcut.push_back("10to20");
|
| 66 |
lepcut.push_back("two10to15");
|
| 67 |
lepcut.push_back("two10to20");
|
| 68 |
lepcut.push_back("two20to7000");
|
| 69 |
lepcut.push_back("conf1010");
|
| 70 |
lepcut.push_back("conf1510");
|
| 71 |
lepcut.push_back("conf2010");
|
| 72 |
lepcut.push_back("conf1515");
|
| 73 |
lepcut.push_back("conf2015");
|
| 74 |
|
| 75 |
int Mass0min = 0;
|
| 76 |
int Mass0max = 3000;
|
| 77 |
int Mass12min = 0;
|
| 78 |
int Mass12max = 1000;
|
| 79 |
|
| 80 |
vector<float> jzbcut;
|
| 81 |
jzbcut.push_back(0);
|
| 82 |
jzbcut.push_back(50);
|
| 83 |
jzbcut.push_back(100);
|
| 84 |
jzbcut.push_back(150);
|
| 85 |
jzbcut.push_back(200);
|
| 86 |
jzbcut.push_back(250);
|
| 87 |
|
| 88 |
vector<string> predobs;
|
| 89 |
predobs.push_back("pred");
|
| 90 |
predobs.push_back("obs");
|
| 91 |
predobs.push_back("predOF");
|
| 92 |
predobs.push_back("predSF");
|
| 93 |
predobs.push_back("final");
|
| 94 |
|
| 95 |
|
| 96 |
vector<string> mllcut;
|
| 97 |
mllcut.push_back("onpeak");
|
| 98 |
mllcut.push_back("from20");
|
| 99 |
mllcut.push_back("from30");
|
| 100 |
mllcut.push_back("from40");
|
| 101 |
mllcut.push_back("from50");
|
| 102 |
mllcut.push_back("from60");
|
| 103 |
|
| 104 |
|
| 105 |
|
| 106 |
TH2F *histos[predobs.size()][lepcut.size()][mllcut.size()][jzbcut.size()];
|
| 107 |
|
| 108 |
for(int ip=0;ip<predobs.size();ip++) {
|
| 109 |
for(int il=0;il<lepcut.size();il++) {
|
| 110 |
for(int im=0;im<mllcut.size();im++) {
|
| 111 |
for(int ij=0;ij<jzbcut.size();ij++) {
|
| 112 |
stringstream tlist;
|
| 113 |
tlist << "efficiency_" << predobs[ip] << "_lep" << lepcut[il] << "_mll_" << mllcut[im] << "_jzb" << jzbcut[ij];
|
| 114 |
histos[ip][il][im][ij] = new TH2F(tlist.str().c_str(),tlist.str().c_str(),(Mass0max-Mass0min)/20.0,Mass0min,Mass0max,(Mass12max-Mass12min)/20.0,Mass12min,Mass12max);
|
| 115 |
}
|
| 116 |
}
|
| 117 |
}
|
| 118 |
}
|
| 119 |
TH2F *Nevents = new TH2F("Nevents","Nevents",(Mass0max-Mass0min)/20.0,Mass0min,Mass0max,(Mass12max-Mass12min)/20.0,Mass12min,Mass12max);
|
| 120 |
|
| 121 |
int pred=0;
|
| 122 |
int obs=1;
|
| 123 |
int final=2;
|
| 124 |
|
| 125 |
|
| 126 |
for(int i=0;i<uppercutoff;i++) {
|
| 127 |
events->GetEntry(i);
|
| 128 |
if(i%10000 == 0) cout << "\033[1G Working on event " << i << "\033[5G" << flush;
|
| 129 |
int GlobalBin=Nevents->FindBin(M0,M12);
|
| 130 |
Nevents->SetBinContent(GlobalBin,Nevents->GetBinContent(GlobalBin)+1);
|
| 131 |
if(ch1*ch2>0) continue; //same sign ...
|
| 132 |
if(pfJetGoodNum<3) continue; // not enough jets
|
| 133 |
if(!((pfJetGoodNum>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum>=2&&pfJetGoodID[1]!=0))) continue; // jet quality cut
|
| 134 |
float rjzb=(jzb[1]+0.06*pt-3.58612);
|
| 135 |
int rpredobs;
|
| 136 |
int rpredobssign;
|
| 137 |
int categories[lepcut.size()];
|
| 138 |
int mllcategories[mllcut.size()];
|
| 139 |
int jzbregion[jzbcut.size()];
|
| 140 |
int predobscat[predobs.size()];
|
| 141 |
|
| 142 |
if((pt1>10&&pt1<15)||(pt2>10&&pt2<15)) categories[0]=1;
|
| 143 |
else categories[0]=0;
|
| 144 |
if((pt1>15&&pt1<20)||(pt2>15&&pt2<20)) categories[1]=1;
|
| 145 |
else categories[1]=0;
|
| 146 |
if((pt1>10&&pt1<20)||(pt2>10&&pt2<20)) categories[2]=1;
|
| 147 |
else categories[2]=0;
|
| 148 |
if((pt1>10&&pt1<15)&&(pt2>10&&pt2<15)) categories[3]=1;
|
| 149 |
else categories[3]=0;
|
| 150 |
if((pt1>10&&pt1<20)&&(pt2>10&&pt2<20)) categories[4]=1;
|
| 151 |
else categories[4]=0;
|
| 152 |
if((pt1>20)&&(pt2>20)) categories[5]=1;
|
| 153 |
else categories[5]=0;
|
| 154 |
if(pt1>10&&pt2>10) categories[6]=1; //conf1010
|
| 155 |
else categories[6]=0;
|
| 156 |
if(pt1>15&&pt2>10) categories[7]=1; //conf1510
|
| 157 |
else categories[7]=0;
|
| 158 |
if(pt1>20&&pt2>10) categories[8]=1; //conf2010
|
| 159 |
else categories[8]=0;
|
| 160 |
if(pt1>15&&pt2>15) categories[9]=1; //conf1515
|
| 161 |
else categories[9]=0;
|
| 162 |
if(pt1>20&&pt2>15) categories[10]=1; //conf2015
|
| 163 |
else categories[10]=0;
|
| 164 |
|
| 165 |
if(!(rjzb>0&&id1==id2)) {
|
| 166 |
predobscat[0]=1;
|
| 167 |
predobscat[1]=0;
|
| 168 |
} else {
|
| 169 |
predobscat[0]=0;
|
| 170 |
predobscat[1]=1;
|
| 171 |
}
|
| 172 |
if(id1!=id2) predobscat[2]=1;
|
| 173 |
else predobscat[2]=0;
|
| 174 |
if(id1==id2&&rjzb<0) predobscat[3]=1;
|
| 175 |
else predobscat[3]=0;
|
| 176 |
|
| 177 |
if(TMath::Abs(mll-91.2)<20.0) mllcategories[0]=1;
|
| 178 |
else mllcategories[0]=0;
|
| 179 |
if(mll>20.0) mllcategories[1]=1;
|
| 180 |
else mllcategories[1]=0;
|
| 181 |
if(mll>30.0) mllcategories[2]=1;
|
| 182 |
else mllcategories[2]=0;
|
| 183 |
if(mll>40.0) mllcategories[3]=1;
|
| 184 |
else mllcategories[3]=0;
|
| 185 |
if(mll>50.0) mllcategories[4]=1;
|
| 186 |
else mllcategories[4]=0;
|
| 187 |
if(mll>60.0) mllcategories[5]=1;
|
| 188 |
else mllcategories[5]=0;
|
| 189 |
|
| 190 |
|
| 191 |
|
| 192 |
|
| 193 |
for(int ij=0;ij<jzbcut.size();ij++) {
|
| 194 |
if(TMath::Abs(rjzb)>jzbcut[ij]) jzbregion[ij]=1;
|
| 195 |
else jzbregion[ij]=0;
|
| 196 |
}
|
| 197 |
|
| 198 |
if(rjzb>0&&id1==id2) rpredobs=obs;
|
| 199 |
else rpredobs=pred;
|
| 200 |
if(rjzb<0&&id1!=id2) rpredobssign=-1; // this is ttbar JZB<0
|
| 201 |
else rpredobssign=1;
|
| 202 |
|
| 203 |
|
| 204 |
|
| 205 |
|
| 206 |
|
| 207 |
if(beverbose) cout << endl << "____________________________________________ " << endl;
|
| 208 |
for(int ip=0;ip<predobs.size();ip++) {
|
| 209 |
if(!predobscat[ip]) continue;
|
| 210 |
for(int il=0;il<lepcut.size();il++) {
|
| 211 |
if(!categories[il]) continue;
|
| 212 |
for(int im=0;im<mllcut.size();im++) {
|
| 213 |
if(!mllcategories[im]) continue;
|
| 214 |
for(int ij=0;ij<jzbcut.size();ij++) {
|
| 215 |
if(!jzbregion[ij]) continue;
|
| 216 |
int GlobalBin=histos[ip][il][im][ij]->FindBin(M0,M12);
|
| 217 |
float currcontent=histos[ip][il][im][ij]->GetBinContent(GlobalBin);
|
| 218 |
if(ip!=predobs.size()-1) {
|
| 219 |
histos[ip][il][im][ij]->SetBinContent(GlobalBin,currcontent+rpredobssign);
|
| 220 |
if(beverbose) cout << " Written " << rpredobssign << " to " << histos[ip][il][im][ij]->GetName() << endl;
|
| 221 |
}
|
| 222 |
else {
|
| 223 |
if(rpredobs==obs) {
|
| 224 |
histos[ip][il][im][ij]->SetBinContent(GlobalBin,currcontent+rpredobssign);
|
| 225 |
if(beverbose) cout << " Written " << rpredobssign << " to " << histos[ip][il][im][ij]->GetName() << endl;
|
| 226 |
|
| 227 |
}
|
| 228 |
else {
|
| 229 |
histos[ip][il][im][ij]->SetBinContent(GlobalBin,currcontent-rpredobssign);
|
| 230 |
if(beverbose) cout << " Written " << -rpredobssign << " to " << histos[ip][il][im][ij]->GetName() << endl;
|
| 231 |
}
|
| 232 |
}
|
| 233 |
}
|
| 234 |
}
|
| 235 |
}
|
| 236 |
}
|
| 237 |
|
| 238 |
|
| 239 |
|
| 240 |
if(beverbose) cout << "Event " << i << " : " << "id1=" << id1 << " , id2=" << id2 << " , ch1=" << ch1 << " , ch2=" << ch2 << ", mll=" << mll << " , jzb[1]=" << jzb[1] << " , nJets=" << pfJetGoodNum << " , pt1=" << pt1 << " , pt2=" << pt2 << " , M0=" << M0 << " , M12=" << M12 << endl;
|
| 241 |
}
|
| 242 |
|
| 243 |
cout << "Have processed " << Nevents->GetEntries() << " events." << endl;
|
| 244 |
TFile *out = new TFile("NEWalloutputWithMll2complete.root","RECREATE");
|
| 245 |
for(int ip=0;ip<predobs.size();ip++) {
|
| 246 |
for(int il=0;il<lepcut.size();il++) {
|
| 247 |
for(int im=0;im<mllcut.size();im++) {
|
| 248 |
for(int ij=0;ij<jzbcut.size();ij++) {
|
| 249 |
histos[ip][il][im][ij]->Write();
|
| 250 |
}
|
| 251 |
}
|
| 252 |
}
|
| 253 |
}
|
| 254 |
Nevents->Write();
|
| 255 |
out->Close();
|
| 256 |
f->Close();
|
| 257 |
|
| 258 |
return 0;
|
| 259 |
}
|
