ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ZbTools.C
(Generate patch)

Comparing UserCode/cbrown/Development/Plotting/Modules/ZbTools.C (file contents):
Revision 1.6 by buchmann, Fri Aug 17 13:06:08 2012 UTC vs.
Revision 1.12 by buchmann, Mon Nov 19 16:21:52 2012 UTC

# Line 23 | Line 23 | using namespace PlottingSetup;
23   namespace ZbData {
24   vector<float> data_over_mc;
25  
26 < TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120");
27 < TCut LeadingB("Zb30_bTagProbCSVBP[0]>0.679");
28 < TCut EtaB("abs(Zb30_pfJetGoodEta[0])<1.3");
29 < TCut PhiZcut("abs(pfJetDphiZ[0])>2.7");
26 > TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
27 > TCut LeadingB("Zb3010_bTagProbCSVBP[0]>0.898");
28 > TCut EtaB("abs(Zb3010_pfJetGoodEta[0])<1.3");
29 > TCut PhiZcut("abs(Zb3010_pfJetDphiZ[0])>2.7");
30    
31   }
32  
33   using namespace ZbData;
34  
35 +
36 + //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BETWEEN THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
37 + void SpecialCutFlow(string identifier);
38 +  
39 + //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*            AND THIS LINE             /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
40 +
41 +
42   void print_yield(TCut cut) {
43    TH1F *data  = allsamples.Draw("data",    "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,data,luminosity);
44 <  cout << "  Yields for " << (const char*) cut << endl;
45 <  cout << "     data: " << data->Integral() << endl;
44 >  dout << "  Yields for " << (const char*) cut << endl;
45 >  dout << "     data: " << data->Integral() << endl;
46    THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,mc,luminosity);
47    int nhists = mcm.GetHists()->GetEntries();
48 <  cout << "Going to loop over " << nhists << " histograms " << endl;
48 >  dout << "Going to loop over " << nhists << " histograms " << endl;
49    TFile *test = new TFile("test.root","RECREATE");
50    mcm.Write();
51    test->Close();
# Line 46 | Line 53 | void print_yield(TCut cut) {
53      TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
54      cout << "     " << hist->GetName() << ": " << hist->Integral() << endl;
55    }
56 +    
57    cout << endl << endl;
58    
59    delete data;
# Line 103 | Line 111 | Value getfrom2Dmap(TH2F *map, int ixbin,
111    
112   void print_all_b_yields() {
113    cout << "Basic selection with a b jet" << endl;
114 <  print_yield(ZplusBsel&&"Zb30_pfJetGoodNumBtag>0");
114 >  print_yield(ZplusBsel&&"Zb3010_pfJetGoodNumBtag>0");
115    cout << "Leading jet is a b-jet " << endl;
116    print_yield(ZplusBsel&&LeadingB);
117    cout << "|#eta_{b}|<1.3 " << endl;
# Line 113 | Line 121 | void print_all_b_yields() {
121    cout << "10<ptZ<1000 GeV" << endl;
122    print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
123    cout << "#alpha < 0.3" << endl;
124 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb30_pfJetGoodPt[1]/pt)<0.3");
124 >  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.3");
125    cout << "#alpha < 0.2" << endl;
126 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb30_pfJetGoodPt[1]/pt)<0.2");
126 >  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.2");
127    cout << "#alpha < 0.15" << endl;
128 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb30_pfJetGoodPt[1]/pt)<0.15");
128 >  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.15");
129    cout << "#alpha < 0.1" << endl;
130 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb30_pfJetGoodPt[1]/pt)<0.1");
131 <  cout << "#alpha < 0.05" << endl;
132 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb30_pfJetGoodPt[1]/pt)<0.05");
130 >  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.1");
131 > //  cout << "#alpha < 0.05" << endl;
132 > //  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.05");
133   }  
134  
135 + void DrawEvilCutFlow() {
136 +  stringstream MegaCut;
137 +  
138 +  MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut("Zb3010_pfJetGoodNumBtag>0")) << ")";
139 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
140 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
141 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
142 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000")) << ")";
143 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.3")) << ")";
144 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.2")) << ")";
145 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.15")) << ")";
146 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.1")) << ")";
147 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.05")) << ")";
148 +  MegaCut << " ";
149 +  
150 +  TCanvas *can = new TCanvas("can","can");
151 +  can->SetLogy(1);
152 +  TCut basecut("mll>6||mll<10");
153 +  
154 +  
155 +  
156 +  TLegend *leg = allsamples.allbglegend();
157 +  
158 +  
159 +  TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
160 +   THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
161 +
162 +  float runningsum=0;
163 +  for(int i=data->GetNbinsX();i>0;i--) {
164 +    runningsum+=data->GetBinContent(i);
165 +    data->SetBinContent(i,runningsum);
166 +  }
167 +  
168 +  float minimum=data->GetMaximum();
169 +  
170 +  TH1F* h;
171 +  TIter nextOF(mcm.GetHists());
172 +  while ( h = (TH1F*)nextOF() ) {
173 +    float runningsum=0;
174 +    minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
175 +    for(int i=h->GetNbinsX();i>0;i--) {
176 +      runningsum+=h->GetBinContent(i);
177 +      h->SetBinContent(i,runningsum);
178 +      if(runningsum<minimum&&runningsum>0) minimum=runningsum;
179 +    }
180 +  }
181 +  
182 +  
183 +
184 +  data->SetMinimum(0.05*minimum);
185 +  can->cd(1)->SetBottomMargin(0.25);
186 +  data->GetXaxis()->SetBinLabel(1,"Pre-selection");
187 +  data->GetXaxis()->SetBinLabel(2,"Contains b jet");
188 +  data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
189 +  data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
190 +  data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
191 +  data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
192 +  data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
193 +  data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
194 +  data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
195 +  data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
196 +  data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
197 +  data->GetXaxis()->LabelsOption("v");
198 +  
199 +  data->Draw("e1");
200 +  mcm.Draw("same");
201 +  data->Draw("same,e1,axis");
202 +  data->Draw("same,e1");
203 + //  data->Draw("same,TEXT");
204 +  
205 +  leg->Draw();
206 +  
207 +  DrawPrelim();
208 +  
209 +  CompleteSave(can,"CutFlow");
210 + }
211 +
212   void draw_Zb_kin_vars() {
213     draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Z p_{T}","Official/Zpt",1);
214 <
215 <   draw_kin_variable("Zb30_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
216 <   draw_kin_variable("Zb30_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
217 <   draw_kin_variable("Zb30_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/SubLeadingJetPt_Over_ZPt",1);
218 <    draw_kin_variable("Zb30_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<50",20,0,2,"#alpha","Official/alpha_pt_10_to_50",1);
219 <    draw_kin_variable("Zb30_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>50&&pt<100",20,0,2,"#alpha","Official/alpha_50_to_100",1);
220 <    draw_kin_variable("Zb30_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>100&&pt<200",20,0,2,"#alpha","Official/alpha_pt_100_to_200",1);
221 <    draw_kin_variable("Zb30_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>200&&pt<1000",20,0,2,"#alpha","Official/alpha_pt_200_to_1000",1);
214 >   draw_kin_variable("Zb3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
215 >   draw_kin_variable("Zb3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
216 >   draw_kin_variable("Zb3010_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/SubLeadingJetPt_Over_ZPt",1);
217 >   draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<50",   40,0,2,"#alpha","Official/alpha_pt_30_to_50",1);
218 >   draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>50&&pt<75",  40,0,2,"#alpha","Official/alpha_50_to_75",1);
219 >   draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>75&&pt<125", 40,0,2,"#alpha","Official/alpha_pt_75_to_125",1);
220 >   draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>125&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
221 >    draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
222    
223 <  draw_kin_variable("pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
223 >  draw_kin_variable("Zb3010_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
224    
225   /*
226    draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l1}","Lep/pt1",1);
# Line 156 | Line 241 | void draw_Zb_kin_vars() {
241  
242   void draw_mpf_vars() {
243    ZbData::data_over_mc.clear();
244 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
245 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
246 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
247 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
248 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
249 <  
250 <  draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.3",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
251 <  draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.2",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
252 <  draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.15",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
253 <  draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.1",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
254 <  draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.05",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
244 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
245 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
246 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
247 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
248 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
249 >  
250 >  draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.3",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
251 >  draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.2",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
252 >  draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.15",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
253 >  draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.1",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
254 >  draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
255   }  
256  
257 < Value get_Zb_data_over_mc(string variable, TCut cut) {
257 > /*double GetMedian(TH1F *histo) {
258 >    int numBins = histo->GetXaxis()->GetNbins();
259 >    Double_t x[numBins];
260 >    Double_t y[numBins];
261 >    for (int i = 1; i <= numBins; i++) {
262 >        x[i] = histo->GetBinCenter(i);
263 >        y[i] = histo->GetBinContent(i);
264 >    }
265 >    return TMath::Median(numBins, &x, &y);
266 > }*/
267 >
268 > Value get_Zb_data_over_mc(string variable, string variable2, TCut cut, string saveas) {
269   //  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
270 <  TH1F *hdata  = allsamples.Draw("hdata",variable,1000,0,30, "#alpha [GeV]", "events", cut,data,luminosity);
271 <  TH1F *hmc    = allsamples.Draw("hmc" , variable,1000,0,30, "#alpha [GeV]", "events", cut,mc,  luminosity);
270 >  TCanvas *cn = new TCanvas("cn","cn");
271 >  string varname="MPF";
272 >  if(Contains(variable,"pfJetGood")) varname="R_{abs}";
273 >  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
274 >  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
275    
276    float a=hdata->GetMean();
277    float b=hmc->GetMean();
# Line 180 | Line 279 | Value get_Zb_data_over_mc(string variabl
279    float db=hmc->GetMeanError();
280    float factor = a / b;
281    float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
282 +  
283 +  cn->cd();
284 +  hdata->GetYaxis()->SetTitle("events / 0.1");
285 +  hdata->Rebin(4);
286 +  hmc->Rebin(4);
287 +  hdata->GetXaxis()->SetRangeUser(0,2);
288 +  hdata->Draw("e1");
289 +  hmc->Draw("histo,same");
290 +  hdata->Draw("e1,same");
291 +  stringstream summary;
292 +  summary << "#splitline{" << varname << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{ratio: " << factor;
293 +  summary << " #pm " << error << "}}}";
294 +  TText *infobox = write_title(summary.str());
295 +  infobox->SetX(0.75);
296 +  infobox->SetTextSize(0.03);
297 +  infobox->SetY(0.75);
298 +  infobox->Draw();
299 +    
300 +    
301 +  DrawPrelim();
302 +  CompleteSave(cn,"ResponsePlots/"+saveas);
303 +  
304 +  cout << "Data;" << a << ";"<<da<<";MC"<<b<<";"<<db<<endl;
305 +  delete cn;
306    delete hdata;
307    delete hmc;
308 +
309 + //    cout << a << "+/-" << da << "        " << b << "=/-" << db << endl;
310 +
311    return Value(factor,error);
312   }
313    
# Line 190 | Line 316 | Value get_Zb_data_over_mc(string variabl
316   void new_data_mc_agreement_2d() {
317    gStyle->SetOptFit(0);
318    TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
319 <  const int nalphacuts=6;
320 <  float alphacuts[nalphacuts] = {0.0,0.1,0.15,0.2,0.25,0.3};
319 >  const int nalphacuts=5;
320 >  float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3};
321    int PointThree=-1;
322    for(int i=0;i<nalphacuts;i++) {
323 <    if((alphacuts[i]-0.3)<0.01) PointThree=i;
323 >    if((alphacuts[i]-0.3)<0.01) PointThree=i-1;
324    }
325    if(PointThree<0) {
326      write_warning(__FUNCTION__,"alpha=0.3 MUST BE ONE OF THE ALPHA VALUES!");
# Line 202 | Line 328 | void new_data_mc_agreement_2d() {
328    }
329    
330    const int nptcuts=5;
331 <  float ptcuts[nptcuts]={10,50,100,200,1000};
331 >  float ptcuts[nptcuts]={30,50,75,125,1000};
332    
333    float MPF_Results[nalphacuts][nptcuts];
334    float MPF_Errors[nalphacuts][nptcuts];
# Line 210 | Line 336 | void new_data_mc_agreement_2d() {
336    float RABS_Errors[nalphacuts][nptcuts];
337    
338    
339 <  for(int ia=0;ia<nalphacuts-1;ia++) {
339 >  for(int ia=0;ia<nalphacuts;ia++) {
340      for(int ipt=0;ipt<nptcuts-1;ipt++) {
341        stringstream specialcut;
342 <      specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb30_alpha<" << alphacuts[ia+1] << "&&Zb30_alpha>" << alphacuts[ia] << "))";
343 <      Value MPF_data_over_mc = get_Zb_data_over_mc("mpf",TCut(basecut && specialcut.str().c_str()));
344 <      Value RABS_data_over_mc = get_Zb_data_over_mc("Zb30_pfJetGoodPt[0]/pt",TCut(basecut && specialcut.str().c_str()));
342 >      specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<" << alphacuts[ia] << "))";
343 >      cout << specialcut.str() << endl;
344 >      Value MPF_data_over_mc = get_Zb_data_over_mc("mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]));
345 >      Value RABS_data_over_mc = get_Zb_data_over_mc("Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]));
346        
347        MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
348        MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
# Line 223 | Line 350 | void new_data_mc_agreement_2d() {
350        RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
351        RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
352        
353 <      cout << "alpha cut at " << alphacuts[ia] << " and pt range " << ptcuts[ipt] << " , " << ptcuts[ipt+1] << " \t : \t" << specialcut.str() << " \t \t --> " << MPF_data_over_mc << " (MPF) " << RABS_data_over_mc << " (RABS)" << endl;
353 >      cout << "alpha cut at " << alphacuts[ia+1] << " and pt range " << ptcuts[ipt] << " , " << ptcuts[ipt+1] << " \t : \t" << specialcut.str() << " \t \t --> " << MPF_data_over_mc << " (MPF) " << RABS_data_over_mc << " (RABS)" << endl;
354      }
355    }
356    
# Line 244 | Line 371 | void new_data_mc_agreement_2d() {
371      TGraphErrors *mpf_gr =new TGraphErrors();
372      TGraphErrors *rabs_gr =new TGraphErrors();
373      for(int ia=0;ia<nalphacuts;ia++) {
374 <      mpf_gr->SetPoint(ia,alphacuts[ia],MPF_Results[ia][ipt]);
375 <      mpf_gr->SetPointError(ia,0,MPF_Errors[ia][ipt]);
376 <      rabs_gr->SetPoint(ia,alphacuts[ia],RABS_Results[ia][ipt]);
377 <      rabs_gr->SetPointError(ia,0,RABS_Errors[ia][ipt]);
374 >      cout << "Setting a point for an alpha cut at " << alphacuts[ia] << " with value " << MPF_Results[ia][ipt] << " (MPF) and " << RABS_Results[ia][ipt] << " (RABS)" << endl;
375 >      mpf_gr->SetPoint(ia,alphacuts[ia],MPF_Results[ia][ipt]);     // avoiding empty bins (nan's)
376 >      mpf_gr->SetPointError(ia,0,MPF_Errors[ia][ipt]);             // avoiding empty bins (nan's)
377 >      rabs_gr->SetPoint(ia,alphacuts[ia],RABS_Results[ia][ipt]); // avoiding empty bins (nan's)
378 >      rabs_gr->SetPointError(ia,0,RABS_Errors[ia][ipt]);         // avoiding empty bins (nan's)
379      }
380 +    cout << "About to add a point in pt (" << ptcuts[ipt]+ptcuts[ipt+1] << " ) to face value : " << MPF_Results[PointThree][ipt] << endl;
381      FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1]),MPF_Results[PointThree][ipt]);
382      FaceValueAtPoint3->SetPointError(ipt,0,MPF_Errors[PointThree][ipt]);
383      
# Line 258 | Line 387 | void new_data_mc_agreement_2d() {
387      mpf_gr->Fit("pol1","QE"); // quiet, use minos
388      mpf_gr->GetXaxis()->SetTitle("#alpha");
389      mpf_gr->GetXaxis()->CenterTitle();
390 <    mpf_gr->GetYaxis()->SetTitle("data/mc");
390 >    mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
391      mpf_gr->GetYaxis()->CenterTitle();
392      TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
393      float mpf_a=mpf_pol->GetParameter(0);
# Line 291 | Line 420 | void new_data_mc_agreement_2d() {
420      rabs_gr->Draw("AP");
421      rabs_gr->GetXaxis()->SetTitle("#alpha");
422      rabs_gr->GetXaxis()->CenterTitle();
423 <    rabs_gr->GetYaxis()->SetTitle("data/mc");
423 >    rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
424      rabs_gr->GetYaxis()->CenterTitle();
425      
426      TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
# Line 328 | Line 457 | void new_data_mc_agreement_2d() {
457    MPF_FinalGraph->Draw("AP");
458    MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
459    MPF_FinalGraph->GetXaxis()->CenterTitle();
460 <  MPF_FinalGraph->GetYaxis()->SetTitle("data/mc");
460 >  MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
461    MPF_FinalGraph->GetYaxis()->CenterTitle();
462    MPF_FinalGraph->Draw("AP");
463    
# Line 355 | Line 484 | void new_data_mc_agreement_2d() {
484    RABS_FinalGraph->Draw("AP");
485    RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
486    RABS_FinalGraph->GetXaxis()->CenterTitle();
487 <  RABS_FinalGraph->GetYaxis()->SetTitle("data/mc");
487 >  RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
488    RABS_FinalGraph->GetYaxis()->CenterTitle();
489    RABS_FinalGraph->Draw("AP");
490    TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
# Line 375 | Line 504 | void new_data_mc_agreement_2d() {
504    filename2 << "Extrapolation/FINAL_RESULT_RABS";
505    CompleteSave(can,filename2.str());
506    
507 <  
507 >  can->cd();
508    TCanvas *FaceValueCanvas = new TCanvas("FaceValueCanvas","FaceValueCanvas");
509    FaceValueAtPoint3->SetMarkerStyle(21);
510    FaceValueAtPoint3->SetMarkerStyle(21);
# Line 400 | Line 529 | void new_data_mc_agreement_2d() {
529    float FaceValResultError=faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0));
530    
531    filename2.str("");
532 <  filename2 << "Extrapolation/FINAL_RESULT_RABS";
532 >  filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
533    CompleteSave(can,filename2.str());
534    
535    cout << "FINAL RESULTS: " << endl;
# Line 414 | Line 543 | void new_data_mc_agreement_2d() {
543   }
544  
545  
546 + void DoPUStudy(string identifier) {
547 +  
548 +  float numVtxcuts[7]={0,10,15,20};
549 +  
550 +  TCanvas *can = new TCanvas("can","can");
551 +  TH1F *malpha[8];
552 +  TH1F *dalpha[8];
553 +  cout << identifier << ";";
554 +  cout << endl;
555 +  for(int i=1;i<5;i++) {
556 +    stringstream sSpecialCut;
557 +    if(i<4) {
558 +      sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
559 +      cout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
560 +    } else {
561 +      sSpecialCut << "numVtx>" << numVtxcuts[i-1];
562 +      cout << numVtxcuts[i-1] << ";";
563 +    }
564 +    
565 +    TCut SpecialCut(sSpecialCut.str().c_str());
566 +    
567 +    
568 +    malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
569 +    dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
570 +    
571 +    
572 +    float a=dalpha[i]->GetMean();
573 +    float b=malpha[i]->GetMean();
574 +    float da=dalpha[i]->GetMeanError();
575 +    float db=malpha[i]->GetMeanError();
576 +    float factor = a / b;
577 +    float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
578 +  
579 +    cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
580 +    malpha[i]->SetLineColor(i+1);
581 +    if(i==0) malpha[i]->Draw("histo");
582 +    else malpha[i]->Draw("histo,same");
583 +  }
584 +  
585 +  CompleteSave(can,"AlphaOverview_"+identifier);
586 +  delete can;
587 +  
588 + /*  for(int i=1;i<8;i++) {
589 +    delete malpha[i];
590 +    delete dalpha[i];
591 +  }
592 +  */
593 +  
594 + }
595 +
596 + void ScenarioComparison() {
597 +  
598 +  TGraphAsymmErrors *gr[5];
599 +  TF1 *fit[5];
600 +  bool dofit=false;
601 +  
602 +  TCanvas *can = new TCanvas("can","can",500,500);
603 +  
604 +  TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
605 +  leg->SetBorderSize(0);
606 +  leg->SetFillColor(kWhite);
607 +  
608 +  for(int i=0;i<5;i++) {
609 +    gr[i] = new TGraphAsymmErrors();
610 +    double x[4]   = {5,12.5,17.5,20};
611 +    double dxl[4] = {5,2.5,2.5,2.5};
612 +    double dxh[4] = {5,2.5,2.5,10};
613 +    
614 +    string name="";
615 +    string color="";
616 +    
617 +    if(i==0) {
618 +      double y[4]   = {1.008,1.013,1.003,1.007};
619 +      double dyl[4] = {0.011,0.009,0.013,0.020};
620 +      double dyh[4] = {0.011,0.009,0.013,0.020};
621 +      name="30/30";
622 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
623 +      color="#a01d99";
624 +    }
625 +    if(i==1) {
626 +      double y[4]   = {1.008+0.001,1.013+0.001,1.003+0.001,1.007+0.001};
627 +      double dyl[4] = {0.011,0.009,0.013,0.020};
628 +      double dyh[4] = {0.011,0.009,0.013,0.020};
629 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
630 +      name="30/10";
631 +      color="#2464bc";
632 +    }
633 +    if(i==2) {
634 +      double y[4]   = {1.014,1.026,1.014,1.059};
635 +      double dyl[4] = {0.011,0.010,0.014,0.023};
636 +      double dyh[4] = {0.011,0.010,0.014,0.023};
637 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
638 +      name="15/15";
639 +      color="#f49230";
640 +    }
641 +    if(i==3) {
642 +      double y[4]   = {1.015,1.020,1.008,1.050};
643 +      double dyl[4] = {0.011,0.010,0.014,0.022};
644 +      double dyh[4] = {0.011,0.010,0.014,0.022};
645 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
646 +      name="20/15";
647 +      color="#c72f3c";
648 +    }
649 +    if(i==4) {
650 +      double y[4]   = {1.015+0.001,1.020+0.001,1.008+0.001,1.050+0.001};
651 +      double dyl[4] = {0.011,0.010,0.014,0.022};
652 +      double dyh[4] = {0.011,0.010,0.014,0.022};
653 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
654 +      name="20/20";
655 +      color="#7cb433";
656 +    }
657 +    
658 +    
659 +    gr[i]->SetTitle();
660 +    gr[i]->GetXaxis()->SetTitle("N(Vtx)");
661 +    gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
662 +    gr[i]->GetXaxis()->CenterTitle();
663 +    gr[i]->GetYaxis()->CenterTitle();
664 +    gr[i]->SetName(name.c_str());
665 +    gr[i]->SetMarkerSize(0.00001);
666 +    gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
667 +    gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
668 +    gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
669 +    if(dofit) {
670 +      gr[i]->Fit("pol1");
671 +      fit[i] = (TF1*)gr[i]->GetFunction("pol1");
672 +      fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
673 +    }
674 +    leg->AddEntry(gr[i],name.c_str(),"l");
675 +    if(i==0) gr[i]->Draw("AP*");
676 +    else gr[i]->Draw("P");
677 +  }
678 +  leg->Draw();
679 +  DrawPrelim();
680 +  
681 +  CompleteSave(can,"ScenarioComparison");
682 +  delete can;
683 +  
684 + }
685 +
686 + void ScenarioComparisonInclusive() {
687 +  
688 +  TGraphAsymmErrors *gr[5];
689 +  TF1 *fit[5];
690 +  bool dofit=true;
691 +  
692 +  TCanvas *can = new TCanvas("can","can",500,500);
693 +  
694 +  TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
695 +  leg->SetBorderSize(0);
696 +  leg->SetFillColor(kWhite);
697 +  
698 +  for(int i=0;i<5;i++) {
699 +    gr[i] = new TGraphAsymmErrors();
700 +    double x[4]   = {5,12.5,17.5,20};
701 +    double dxl[4] = {5,2.5,2.5,2.5};
702 +    double dxh[4] = {5,2.5,2.5,10};
703 +    
704 +    string name="";
705 +    string color="";
706 +    
707 +    if(i==0) {
708 +      double y[4]   = {1.01103,1.01496,1.01486,1.02785};
709 +      double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
710 +      double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
711 +      name="30/30";
712 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
713 +      color="#a01d99";
714 +    }
715 +    if(i==1) {
716 +      double y[4]   = {1.01103+0.001,1.01496+0.001,1.01486+0.001,1.02785+0.001};
717 +      double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
718 +      double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
719 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
720 +      name="30/10";
721 +      color="#2464bc";
722 +    }
723 +    if(i==2) {
724 +      double y[4]   = {1.01575,1.02879,1.0325,1.05456};
725 +      double dyl[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
726 +      double dyh[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
727 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
728 +      name="15/15";
729 +      color="#f49230";
730 +    }
731 +    if(i==3) {
732 +      double y[4]   = {1.01221,1.01922,1.02396,1.04677};
733 +      double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
734 +      double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
735 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
736 +      name="20/15";
737 +      color="#c72f3c";
738 +    }
739 +    if(i==4) {
740 +      double y[4]   = {1.01221+0.001,1.01922+0.001,1.02396+0.001,1.04677+0.001};
741 +      double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
742 +      double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
743 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
744 +      name="20/20";
745 +      color="#7cb433";
746 +    }
747 +    
748 +    
749 +    gr[i]->SetTitle();
750 +    gr[i]->GetXaxis()->SetTitle("N(Vtx)");
751 +    gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
752 +    gr[i]->GetXaxis()->CenterTitle();
753 +    gr[i]->GetYaxis()->CenterTitle();
754 +    gr[i]->SetName(name.c_str());
755 +    gr[i]->SetMarkerSize(0.00001);
756 +    gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
757 +    gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
758 +    gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
759 +    if(dofit) {
760 +      gr[i]->Fit("pol1");
761 +      fit[i] = (TF1*)gr[i]->GetFunction("pol1");
762 +      fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
763 +    }
764 +    leg->AddEntry(gr[i],name.c_str(),"l");
765 +    if(i==0) gr[i]->Draw("AP*");
766 +    else gr[i]->Draw("P");
767 +  }
768 +  leg->Draw();
769 +  DrawPrelim();
770 +  
771 +  CompleteSave(can,"ScenarioComparisonInclusive");
772 +  delete can;
773 +  
774 + }
775 +
776 + void compare_selection(string identifier) {
777 +  bool recognized_scenario=false;
778 +  
779 +  cout << "Running with identifier " << identifier << endl;
780 +  if(identifier=="Zb3010") {
781 +      LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
782 +      EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
783 +      PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
784 +      recognized_scenario=true;
785 +  }
786 +  
787 +  if(identifier=="Zb3010") {
788 +      LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
789 +      EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
790 +      PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
791 +      recognized_scenario=true;
792 +  }
793 +  
794 +  if(identifier=="Zb40") {
795 +      LeadingB=TCut ("Zb40_bTagProbCSVBP[0]>0.898");
796 +      EtaB=TCut("abs(Zb40_pfJetGoodEta[0])<1.3");
797 +      PhiZcut=TCut("abs(Zb40_pfJetDphiZ[0])>2.7");
798 +      recognized_scenario=true;
799 +  }
800 +  
801 +  if(identifier=="Zb301030") {
802 +      LeadingB=TCut ("Zb301030_bTagProbCSVBP[0]>0.898");
803 +      EtaB=TCut("abs(Zb301030_pfJetGoodEta[0])<1.3");
804 +      PhiZcut=TCut("abs(Zb301030_pfJetDphiZ[0])>2.7");
805 +      recognized_scenario=true;
806 +  }
807 +  
808 +  if(identifier=="Zb1530") {
809 +      LeadingB=TCut ("Zb1530_bTagProbCSVBP[0]>0.898");
810 +      EtaB=TCut("abs(Zb1530_pfJetGoodEta[0])<1.3");
811 +      PhiZcut=TCut("abs(Zb1530_pfJetDphiZ[0])>2.7");
812 +      recognized_scenario=true;
813 +  }
814 +  
815 +  if(identifier=="Zb301010") {
816 +      LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
817 +      EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
818 +      PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
819 +      recognized_scenario=true;
820 +  }
821 +  
822 +  if(identifier=="Zb1510") {
823 +      LeadingB=TCut ("Zb1510_bTagProbCSVBP[0]>0.898");
824 +      EtaB=TCut("abs(Zb1510_pfJetGoodEta[0])<1.3");
825 +      PhiZcut=TCut("abs(Zb1510_pfJetDphiZ[0])>2.7");
826 +      recognized_scenario=true;
827 +  }
828 +  
829 +  if(identifier=="Zb301010") {
830 +      LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
831 +      EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
832 +      PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
833 +      recognized_scenario=true;
834 +  }
835 +  
836 +  if(identifier=="ZbMikko") {
837 +      LeadingB=TCut ("ZbMikko_bTagProbCSVBP[0]>0.898 && ZbMikko__IsClean");
838 +      EtaB=TCut("abs(ZbMikko_pfJetGoodEta[0])<1.3 && ZbMikko__IsClean");
839 +      PhiZcut=TCut("abs(ZbMikko_pfJetDphiZ[0])>2.7 && ZbMikko__IsClean");
840 +      recognized_scenario=true;
841 +  }
842 +  
843 +  if(identifier=="Zb3010_SecEta3") {
844 +      LeadingB=TCut ("Zb3010_SecEta3_bTagProbCSVBP[0]>0.898");
845 +      EtaB=TCut("abs(Zb3010_SecEta3_pfJetGoodEta[0])<1.3");
846 +      PhiZcut=TCut("abs(Zb3010_SecEta3_pfJetDphiZ[0])>2.7");
847 +      recognized_scenario=true;
848 +  }
849 +  
850 +  if(identifier=="Zb3010_SecEta5") {
851 +      LeadingB=TCut ("Zb3010_SecEta5_bTagProbCSVBP[0]>0.898");
852 +      EtaB=TCut("abs(Zb3010_SecEta5_pfJetGoodEta[0])<1.3");
853 +      PhiZcut=TCut("abs(Zb3010_SecEta5_pfJetDphiZ[0])>2.7");
854 +      recognized_scenario=true;
855 +  }
856 +  
857 +  if(identifier=="Zb3010_p5Clean") {
858 +      LeadingB=TCut ("Zb3010_p5Clean_bTagProbCSVBP[0]>0.898 && Zb3010_p5Clean_IsClean");
859 +      EtaB=TCut("abs(Zb3010_p5Clean_pfJetGoodEta[0])<1.3 && Zb3010_p5Clean_IsClean");
860 +      PhiZcut=TCut("abs(Zb3010_p5Clean_pfJetDphiZ[0])>2.7 && Zb3010_p5Clean_IsClean");
861 +      recognized_scenario=true;
862 +  }
863 +  
864 +  if(identifier=="Zb3010CHS") {
865 +      LeadingB=TCut ("Zb3010CHS_bTagProbCSVBP[0]>0.898");
866 +      EtaB=TCut("abs(Zb3010CHS_pfJetGoodEta[0])<1.3");
867 +      PhiZcut=TCut("abs(Zb3010CHS_pfJetDphiZ[0])>2.7");
868 +      recognized_scenario=true;
869 +  }
870 +  
871 +  assert(recognized_scenario);
872 +  
873 +  cout << "Cuts: " << endl;
874 +  cout << "   " << (const char*) LeadingB << endl;
875 +  cout << "   " << (const char*) EtaB << endl;
876 +  cout << "   " << (const char*) PhiZcut << endl;
877 +  
878 +  cout << endl << endl;
879 +  
880 + //    ScenarioComparison();
881 + //    ScenarioComparisonInclusive();
882 +  
883 +  
884 + //     DoPUStudy(identifier);
885 +  
886 + //   SpecialCutFlow(identifier);
887 + //    draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
888 + //   draw_kin_variable(identifier+"_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",40,0,200,"Leading Jet Pt (B)","Official/"+identifier+"/LeadingJetPt__"+identifier,1);
889 + //   draw_kin_variable(identifier+"_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",40,0,200,"Sub-Leading Jet Pt [GeV]","Official/"+identifier+"/Jet2Pt__"+identifier,1);
890 + //   draw_kin_variable(identifier+"_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/"+identifier+"/SubLeadingJetPt_Over_ZPt__"+identifier,1);
891 +  draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<50"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_30_to_50__"+identifier,1);
892 +  draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),20,0,2,"#alpha","Official/"+identifier+"/alpha_50_to_75__"+identifier,1);
893 +  draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>75&&pt<125"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_75_to_125__"+identifier,1);
894 +  draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>125&&pt<1000"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_125_to_1000__"+identifier,1);
895 +  draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut,20,0,2,"#alpha","Official/"+identifier+"/alpha_full__"+identifier,1);
896 + //   draw_kin_variable(identifier+"_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","Official/"+identifier+"/DeltaPhiZBlead__"+identifier,0);
897 +  
898 + //   draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"Z+b MPF","Official/"+identifier+"/mpf_alpha_smaller_0p3___"+identifier,0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
899 + //   draw_kin_variable(identifier+"_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"p_{T} b-jet / p_{T} Z","Official/"+identifier+"/ptb_over_ptz___alpha_smaller_0p3______"+identifier,0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
900 +
901 +  
902 + //   */
903 +  
904 + }
905  
906 < void do_basic_ZB_analysis() {
906 > void do_basic_ZB_analysis(int do_inclusive) {
907 >  
908 >  //https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagPOG#2011_Data_and_MC
909 >  //https://twiki.cern.ch/twiki/pub/CMS/BtagPOG/SFb-ttbar_payload.txt
910 >  
911 >  if(do_inclusive) {
912 >    LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>=-10000");
913 >    cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)*((id1==id2&&id1==0)*0.95+(id1==id2&&id1==1)*0.88+(id1!=id2)*0.92)))");
914 >    write_info(__FUNCTION__,"Inclusive selection - no btag requirement, no btag efficiency/mistag correction (weight)");
915 >  } else {
916 >    LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
917 >    cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)*((id1==id2&&id1==0)*0.95+(id1==id2&&id1==1)*0.88+(id1!=id2)*0.92)*Zb3010_BTagWgt))");
918 >    write_info(__FUNCTION__,"Exclusive selection - switched on btag requirement and btag efficiency/mistag correction (weight)");
919 >  }
920 >  
921 >  bool doquick=true;
922 >  
923    TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
924    
925 < //  print_all_b_yields();
926 < //  draw_mpf_vars();
927 <  draw_Zb_kin_vars();
925 > //     compare_selection("Zb3010_p5Clean");
926 > //     compare_selection("Zb3010CHS");
927 > //     compare_selection("Zb3010");
928 > //     compare_selection("Zb301030");
929 > //     compare_selection("Zb1530");
930 > //     compare_selection("Zb301010");
931 > //     compare_selection("Zb3010_SecEta3");
932 > //     compare_selection("Zb3010_SecEta5");
933 > //     compare_selection("Zb1510");
934 > //     compare_selection("Zb301010");
935 > //     compare_selection("ZbMikko");
936 > //     compare_selection("Zb3010");
937 > //     compare_selection("Zb40");
938 >
939 >  if(!doquick) {
940 >    print_all_b_yields();
941 >    draw_mpf_vars();
942 >    //  DrawEvilCutFlow();
943 >        
944 >    draw_Zb_kin_vars();
945 >        
946 >   }
947    
948    new_data_mc_agreement_2d();
949    
950 +  write_error(__FUNCTION__,"Need to implement switch for inclusive result; inclusive result needs to have btag efficiency reweighting deactivated!");
951    delete zbcanvas;
952   }
953 +
954 +
955 +
956 +
957 + //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BELOW THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
958 +
959 +
960 + const char* concatenate(string bla1, string bla2) {
961 +  stringstream bla;
962 +  bla << bla1 << bla2;
963 +  return bla.str().c_str();
964 + }
965 +
966 + void SpecialCutFlow(string identifier) {
967 +  stringstream MegaCut;
968 +  
969 +  
970 +  
971 +  MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut(concatenate(identifier,"_pfJetGoodNumBtag>0"))) << ")";
972 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
973 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
974 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
975 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000").c_str())) << ")";
976 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
977 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
978 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
979 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
980 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
981 +  MegaCut << " ";
982 +  
983 +  TCanvas *can = new TCanvas("can","can");
984 +  can->SetLogy(1);
985 +  TCut basecut(ZplusBsel);
986 +  
987 +  
988 +  
989 +  TLegend *leg = allsamples.allbglegend();
990 +  
991 +  
992 +  TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
993 +   THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
994 +
995 +  float runningsum=0;
996 +  for(int i=data->GetNbinsX();i>0;i--) {
997 +    runningsum+=data->GetBinContent(i);
998 +    data->SetBinContent(i,runningsum);
999 +  }
1000 +  
1001 +  float minimum=data->GetMaximum();
1002 +  
1003 +  TH1F* h;
1004 +  TIter nextOF(mcm.GetHists());
1005 +  
1006 +  float mcsum=0;
1007 +  float zb=0;
1008 +  while ( h = (TH1F*)nextOF() ) {
1009 +    float runningsum=0;
1010 +    minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
1011 +    for(int i=h->GetNbinsX();i>0;i--) {
1012 +      runningsum+=h->GetBinContent(i);
1013 +      h->SetBinContent(i,runningsum);
1014 +      if(runningsum<minimum&&runningsum>0) minimum=runningsum;
1015 +    }
1016 +    if(Contains(h->GetName(),"Z_b_")) {
1017 +      zb=h->GetBinContent(h->GetNbinsX());
1018 +    }
1019 +    mcsum+=h->GetBinContent(h->GetNbinsX());
1020 +  }
1021 +  
1022 +  
1023 +
1024 +  data->SetMinimum(0.05*minimum);
1025 +  can->cd(1)->SetBottomMargin(0.25);
1026 +  data->GetXaxis()->SetBinLabel(1,"Pre-selection");
1027 +  data->GetXaxis()->SetBinLabel(2,"Contains b jet");
1028 +  data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
1029 +  data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
1030 +  data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
1031 +  data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
1032 +  data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
1033 +  data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
1034 +  data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
1035 +  data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
1036 +  data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
1037 +  data->GetXaxis()->LabelsOption("v");
1038 +  
1039 +  stringstream purityinfo;
1040 +  purityinfo << "Purity: " << std::setprecision(3) << 100*zb/mcsum << " %";
1041 +  stringstream neventsinfo;
1042 +  neventsinfo << "Nevents: " << data->GetBinContent(data->GetNbinsX());
1043 +  TH1F *crap = new TH1F("crap","",1,0,1);
1044 +  crap->SetLineColor(kWhite);
1045 +  leg->AddEntry(crap,purityinfo.str().c_str(),"l");
1046 +  leg->AddEntry(crap,neventsinfo.str().c_str(),"l");
1047 +  
1048 +  data->Draw("e1");
1049 +  mcm.Draw("same");
1050 +  data->Draw("same,e1,axis");
1051 +  data->Draw("same,e1");
1052 + //  data->Draw("same,TEXT");
1053 +  
1054 +  leg->Draw();
1055 +  
1056 +  DrawPrelim();
1057 +  
1058 +  CompleteSave(can,"CutFlow___"+identifier);
1059 +  
1060 +  delete crap;
1061 +  delete data;
1062 +  
1063 +  
1064 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines