ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/various_assignments/QuickScan/QuickScan.C
Revision: 1.3
Committed: Mon Mar 5 13:44:49 2012 UTC (13 years, 2 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +17 -2 lines
Log Message:
Added categories for the different pt configurations

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 buchmann 1.3 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 buchmann 1.1
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 buchmann 1.3 int categories[lepcut.size()];
138 buchmann 1.1 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 buchmann 1.3 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 buchmann 1.1
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 buchmann 1.2 if(ip!=predobs.size()-1) {
219     histos[ip][il][im][ij]->SetBinContent(GlobalBin,currcontent+rpredobssign);
220 buchmann 1.1 if(beverbose) cout << " Written " << rpredobssign << " to " << histos[ip][il][im][ij]->GetName() << endl;
221 buchmann 1.2 }
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 buchmann 1.1 }
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 buchmann 1.3 TFile *out = new TFile("NEWalloutputWithMll2complete.root","RECREATE");
245 buchmann 1.1 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 buchmann 1.2 }