ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/various_assignments/QuickScan/QuickScan.C
Revision: 1.2
Committed: Tue Feb 28 10:24:14 2012 UTC (13 years, 2 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.1: +15 -2 lines
Log Message:
Fixed an issue with the quick prediction

File Contents

# User Rev Content
1 buchmann 1.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    
70     int Mass0min = 0;
71     int Mass0max = 3000;
72     int Mass12min = 0;
73     int Mass12max = 1000;
74    
75     vector<float> jzbcut;
76     jzbcut.push_back(0);
77     jzbcut.push_back(50);
78     jzbcut.push_back(100);
79     jzbcut.push_back(150);
80     jzbcut.push_back(200);
81     jzbcut.push_back(250);
82    
83     vector<string> predobs;
84     predobs.push_back("pred");
85     predobs.push_back("obs");
86     predobs.push_back("predOF");
87     predobs.push_back("predSF");
88     predobs.push_back("final");
89    
90    
91     vector<string> mllcut;
92     mllcut.push_back("onpeak");
93     mllcut.push_back("from20");
94     mllcut.push_back("from30");
95     mllcut.push_back("from40");
96     mllcut.push_back("from50");
97     mllcut.push_back("from60");
98    
99    
100    
101     TH2F *histos[predobs.size()][lepcut.size()][mllcut.size()][jzbcut.size()];
102    
103     for(int ip=0;ip<predobs.size();ip++) {
104     for(int il=0;il<lepcut.size();il++) {
105     for(int im=0;im<mllcut.size();im++) {
106     for(int ij=0;ij<jzbcut.size();ij++) {
107     stringstream tlist;
108     tlist << "efficiency_" << predobs[ip] << "_lep" << lepcut[il] << "_mll_" << mllcut[im] << "_jzb" << jzbcut[ij];
109     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);
110     }
111     }
112     }
113     }
114     TH2F *Nevents = new TH2F("Nevents","Nevents",(Mass0max-Mass0min)/20.0,Mass0min,Mass0max,(Mass12max-Mass12min)/20.0,Mass12min,Mass12max);
115    
116     int pred=0;
117     int obs=1;
118     int final=2;
119    
120    
121     for(int i=0;i<uppercutoff;i++) {
122     events->GetEntry(i);
123     if(i%10000 == 0) cout << "\033[1G Working on event " << i << "\033[5G" << flush;
124     int GlobalBin=Nevents->FindBin(M0,M12);
125     Nevents->SetBinContent(GlobalBin,Nevents->GetBinContent(GlobalBin)+1);
126     if(ch1*ch2>0) continue; //same sign ...
127     if(pfJetGoodNum<3) continue; // not enough jets
128     if(!((pfJetGoodNum>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum>=2&&pfJetGoodID[1]!=0))) continue; // jet quality cut
129     float rjzb=(jzb[1]+0.06*pt-3.58612);
130     int rpredobs;
131     int rpredobssign;
132     int categories[6];
133     int mllcategories[mllcut.size()];
134     int jzbregion[jzbcut.size()];
135     int predobscat[predobs.size()];
136    
137     if((pt1>10&&pt1<15)||(pt2>10&&pt2<15)) categories[0]=1;
138     else categories[0]=0;
139     if((pt1>15&&pt1<20)||(pt2>15&&pt2<20)) categories[1]=1;
140     else categories[1]=0;
141     if((pt1>10&&pt1<20)||(pt2>10&&pt2<20)) categories[2]=1;
142     else categories[2]=0;
143     if((pt1>10&&pt1<15)&&(pt2>10&&pt2<15)) categories[3]=1;
144     else categories[3]=0;
145     if((pt1>10&&pt1<20)&&(pt2>10&&pt2<20)) categories[4]=1;
146     else categories[4]=0;
147     if((pt1>20)&&(pt2>20)) categories[5]=1;
148     else categories[5]=0;
149    
150     if(!(rjzb>0&&id1==id2)) {
151     predobscat[0]=1;
152     predobscat[1]=0;
153     } else {
154     predobscat[0]=0;
155     predobscat[1]=1;
156     }
157     if(id1!=id2) predobscat[2]=1;
158     else predobscat[2]=0;
159     if(id1==id2&&rjzb<0) predobscat[3]=1;
160     else predobscat[3]=0;
161    
162     if(TMath::Abs(mll-91.2)<20.0) mllcategories[0]=1;
163     else mllcategories[0]=0;
164     if(mll>20.0) mllcategories[1]=1;
165     else mllcategories[1]=0;
166     if(mll>30.0) mllcategories[2]=1;
167     else mllcategories[2]=0;
168     if(mll>40.0) mllcategories[3]=1;
169     else mllcategories[3]=0;
170     if(mll>50.0) mllcategories[4]=1;
171     else mllcategories[4]=0;
172     if(mll>60.0) mllcategories[5]=1;
173     else mllcategories[5]=0;
174    
175    
176    
177    
178     for(int ij=0;ij<jzbcut.size();ij++) {
179     if(TMath::Abs(rjzb)>jzbcut[ij]) jzbregion[ij]=1;
180     else jzbregion[ij]=0;
181     }
182    
183     if(rjzb>0&&id1==id2) rpredobs=obs;
184     else rpredobs=pred;
185     if(rjzb<0&&id1!=id2) rpredobssign=-1; // this is ttbar JZB<0
186     else rpredobssign=1;
187    
188    
189    
190    
191    
192     if(beverbose) cout << endl << "____________________________________________ " << endl;
193     for(int ip=0;ip<predobs.size();ip++) {
194     if(!predobscat[ip]) continue;
195     for(int il=0;il<lepcut.size();il++) {
196     if(!categories[il]) continue;
197     for(int im=0;im<mllcut.size();im++) {
198     if(!mllcategories[im]) continue;
199     for(int ij=0;ij<jzbcut.size();ij++) {
200     if(!jzbregion[ij]) continue;
201     int GlobalBin=histos[ip][il][im][ij]->FindBin(M0,M12);
202     float currcontent=histos[ip][il][im][ij]->GetBinContent(GlobalBin);
203 buchmann 1.2 if(ip!=predobs.size()-1) {
204     histos[ip][il][im][ij]->SetBinContent(GlobalBin,currcontent+rpredobssign);
205 buchmann 1.1 if(beverbose) cout << " Written " << rpredobssign << " to " << histos[ip][il][im][ij]->GetName() << endl;
206 buchmann 1.2 }
207     else {
208     if(rpredobs==obs) {
209     histos[ip][il][im][ij]->SetBinContent(GlobalBin,currcontent+rpredobssign);
210     if(beverbose) cout << " Written " << rpredobssign << " to " << histos[ip][il][im][ij]->GetName() << endl;
211    
212     }
213     else {
214     histos[ip][il][im][ij]->SetBinContent(GlobalBin,currcontent-rpredobssign);
215     if(beverbose) cout << " Written " << -rpredobssign << " to " << histos[ip][il][im][ij]->GetName() << endl;
216     }
217     }
218 buchmann 1.1 }
219     }
220     }
221     }
222    
223    
224    
225     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;
226     }
227    
228     cout << "Have processed " << Nevents->GetEntries() << " events." << endl;
229     TFile *out = new TFile("NEWalloutputWithMll2.root","RECREATE");
230     for(int ip=0;ip<predobs.size();ip++) {
231     for(int il=0;il<lepcut.size();il++) {
232     for(int im=0;im<mllcut.size();im++) {
233     for(int ij=0;ij<jzbcut.size();ij++) {
234     histos[ip][il][im][ij]->Write();
235     }
236     }
237     }
238     }
239     Nevents->Write();
240     out->Close();
241     f->Close();
242    
243     return 0;
244 buchmann 1.2 }