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.4 by buchmann, Fri Aug 10 11:47:53 2012 UTC vs.
Revision 1.10 by buchmann, Wed Nov 14 14:52:54 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("Zb30_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;
60   }
61  
# Line 100 | Line 109 | Value getfrom2Dmap(TH2F *map, int ixbin,
109    return sum;
110   }
111    
103  
104 void data_mc_agreement_2d() {
105  
106  vector<int> emptyvector;
107  float binX[5]={10,50,100,200,1000};
108  float binY[5]={0.0,0.1,0.15,0.2,0.3};
109  vector<float>binningx;
110  vector<float>binningy;
111  binningx.assign(binX,binX+5);
112  binningy.assign(binY,binY+5);
113  
114  TCut cut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
115  
116  TCanvas *agrcan = new TCanvas("agrcan","agrcan");
117  
118  cout << "Drawing 2d for data " << endl;
119  TH2F *alphazptplane_data = allsamples.Draw("alphazptplane_data","alpha:pt",binningx,binningy,"Z p_{T}","#alpha",cut,data,luminosity,emptyvector,false);
120  cout << "Drawing 2d for MC " << endl;
121  TH2F *alphazptplane_mc   = allsamples.Draw("alphazptplane_mc",  "alpha:pt",binningx,binningy,"Z p_{T}","#alpha",cut,mc,  luminosity,emptyvector,false);
122  
123  
124  TFile *f = new TFile("test.root","RECREATE");
125  alphazptplane_data->Write();
126  f->Close();
127  
128  for(int ix=1;ix<binningx.size();ix++) {
129    cout << "pt in " << binningx[ix-1] << " , " << binningx[ix] << endl;
130    TGraphErrors gr;
131    for(int iy=1;iy<binningy.size();iy++) {
132      Value mc   = Value((getfrom2Dmap(alphazptplane_mc,ix,iy)).getValue(),0.0);// not taking MC statistics into account
133      Value data = getfrom2Dmap(alphazptplane_data,ix,iy);
134      cout << " For an alpha cut at " << binningy[iy] << " we get : mc (" << mc << ") , data (" << data << ") --> data/mc = " << data/mc << endl;
135      gr.SetPoint(iy-1,binningy[iy],(data/mc).getValue());
136      gr.SetPointError(iy-1,0,(data/mc).getError());
137    }
138
139    gr.Draw("AP*");
140    gr.GetYaxis()->SetTitle("data / mc");
141    gr.GetYaxis()->CenterTitle();
142    gr.GetXaxis()->SetTitle("#alpha");
143    gr.GetXaxis()->CenterTitle();
144    gr.GetXaxis()->SetRangeUser(0,0.4);
145    gr.GetYaxis()->SetRangeUser(0.5,1.5);
146    gr.Draw("AP*");
147    DrawPrelim();
148    stringstream ptbin;
149    ptbin << binningx[ix-1] << " < p_{T}^{Z} < " << binningx[ix] << " GeV";
150    TText *ptt = write_title(ptbin.str());
151    ptt->SetX(0.7);
152    ptt->SetY(0.7);
153    ptt->SetTextSize(0.03);
154    ptt->Draw();
155    TLine *l = new TLine(0.0,1.0,0.4,1.0);
156    l->SetLineStyle(2);
157    l->SetLineColor(kBlue);
158    l->Draw();
159    
160    stringstream saveas;
161    saveas << "Zplusb_data_over_mc___" << binningx[ix-1] << "to" << binningx[ix];
162    DrawPrelim();
163    CompleteSave(agrcan,saveas.str());
164 //    agrcan->SaveAs("tester.png");
165  }
166  
167  delete agrcan;
168  
169 }
170
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 180 | 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_alpha<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_alpha<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_alpha<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_alpha<0.1");
131 <  cout << "#alpha < 0.05" << endl;
132 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<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("Zb30_pfBJetDphiZ",ZplusBsel&&LeadingB,100,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",1);
214 <  draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Z p_{T}","Zpt",1);
215 <  draw_kin_variable("Zb30_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Sub-Leading Jet Pt","Jet2Pt",1);
216 <  draw_kin_variable("Zb30_pfJetGoodPtlBtag[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Leading Jet Pt (B)","LeadingJetPt",1);
217 <  draw_kin_variable("Zb30_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"#alpha","alpha",1);
213 >   draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Z p_{T}","Official/Zpt",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("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);
227 >  draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,0,200,"p_{T}^{e1}","Lep/pt1_e",1);
228 >  draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,0,200,"p_{T}^{#mu1}","Lep/pt1_m",1);
229 >  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l2}","Lep/pt2",1);
230 >  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,0,200,"p_{T}^{e2}","Lep/pt2_e",1);
231 >  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,0,200,"p_{T}^{#mu2}","Lep/pt2_m",1);
232 >  draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta1",1);
233 >  draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,-3.1,3.1,"#eta^{e1}","Lep/eta1_e",1);
234 >  draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,-3.1,3.1,"#eta^{m1}","Lep/eta1_m",1);
235 >  draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta2",1);
236 >  draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,-3.1,3.1,"#eta^{e2}","Lep/eta2_e",1);
237 >  draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,-3.1,3.1,"#eta^{#mu2}","Lep/eta2_m",1);
238 >  draw_kin_variable("numVtx",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",30,0,30,"N(Vtx)","Lep/nVtx",1);
239 >  */
240   }
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) {
# Line 227 | Line 267 | Value get_Zb_data_over_mc(string variabl
267    float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
268    delete hdata;
269    delete hmc;
270 +
271 + //    cout << a << "+/-" << da << "        " << b << "=/-" << db << endl;
272 +
273    return Value(factor,error);
274   }
275    
# Line 235 | Line 278 | Value get_Zb_data_over_mc(string variabl
278   void new_data_mc_agreement_2d() {
279    gStyle->SetOptFit(0);
280    TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
238  
281    const int nalphacuts=5;
282    float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3};
283 +  int PointThree=-1;
284 +  for(int i=0;i<nalphacuts;i++) {
285 +    if((alphacuts[i]-0.3)<0.01) PointThree=i-1;
286 +  }
287 +  if(PointThree<0) {
288 +    write_warning(__FUNCTION__,"alpha=0.3 MUST BE ONE OF THE ALPHA VALUES!");
289 +    assert(PointThree<0);
290 +  }
291 +  
292    const int nptcuts=5;
293 <  float ptcuts[nptcuts]={10,50,100,200,1000};
293 >  float ptcuts[nptcuts]={30,50,75,125,1000};
294    
295    float MPF_Results[nalphacuts][nptcuts];
296    float MPF_Errors[nalphacuts][nptcuts];
# Line 250 | Line 301 | void new_data_mc_agreement_2d() {
301    for(int ia=0;ia<nalphacuts;ia++) {
302      for(int ipt=0;ipt<nptcuts-1;ipt++) {
303        stringstream specialcut;
304 <      specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (alpha<" << alphacuts[ia] << "))";
304 >      specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<" << alphacuts[ia] << "))";
305 >      cout << specialcut.str() << endl;
306        Value MPF_data_over_mc = get_Zb_data_over_mc("mpf",TCut(basecut && specialcut.str().c_str()));
307 <      Value RABS_data_over_mc = get_Zb_data_over_mc("Zb30_pfJetGoodPt[0]/pt",TCut(basecut && specialcut.str().c_str()));
307 >      Value RABS_data_over_mc = get_Zb_data_over_mc("Zb3010_pfJetGoodPt[0]/pt",TCut(basecut && specialcut.str().c_str()));
308        
309        MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
310        MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
# Line 260 | Line 312 | void new_data_mc_agreement_2d() {
312        RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
313        RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
314        
315 <      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;
315 >      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;
316      }
317    }
318    
# Line 271 | Line 323 | void new_data_mc_agreement_2d() {
323    
324    TGraphErrors *MPF_FinalGraph = new TGraphErrors();
325    TGraphErrors *RABS_FinalGraph = new TGraphErrors();
326 +  TGraphErrors *FaceValueAtPoint3 = new TGraphErrors();
327    
328    for(int ipt=0;ipt<nptcuts-1;ipt++) {
329      
# Line 280 | Line 333 | void new_data_mc_agreement_2d() {
333      TGraphErrors *mpf_gr =new TGraphErrors();
334      TGraphErrors *rabs_gr =new TGraphErrors();
335      for(int ia=0;ia<nalphacuts;ia++) {
336 <      mpf_gr->SetPoint(ia,alphacuts[ia],MPF_Results[ia][ipt]);
337 <      mpf_gr->SetPointError(ia,0,MPF_Errors[ia][ipt]);
338 <      rabs_gr->SetPoint(ia,alphacuts[ia],RABS_Results[ia][ipt]);
339 <      rabs_gr->SetPointError(ia,0,RABS_Errors[ia][ipt]);
336 >      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;
337 >      mpf_gr->SetPoint(ia,alphacuts[ia],MPF_Results[ia][ipt]);     // avoiding empty bins (nan's)
338 >      mpf_gr->SetPointError(ia,0,MPF_Errors[ia][ipt]);             // avoiding empty bins (nan's)
339 >      rabs_gr->SetPoint(ia,alphacuts[ia],RABS_Results[ia][ipt]); // avoiding empty bins (nan's)
340 >      rabs_gr->SetPointError(ia,0,RABS_Errors[ia][ipt]);         // avoiding empty bins (nan's)
341      }
342 +    cout << "About to add a point in pt (" << ptcuts[ipt]+ptcuts[ipt+1] << " ) to face value : " << MPF_Results[PointThree][ipt] << endl;
343 +    FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1]),MPF_Results[PointThree][ipt]);
344 +    FaceValueAtPoint3->SetPointError(ipt,0,MPF_Errors[PointThree][ipt]);
345      
346      mpf_gr->SetMarkerStyle(21);
347      mpf_gr->SetMarkerColor(kBlue);
# Line 292 | Line 349 | void new_data_mc_agreement_2d() {
349      mpf_gr->Fit("pol1","QE"); // quiet, use minos
350      mpf_gr->GetXaxis()->SetTitle("#alpha");
351      mpf_gr->GetXaxis()->CenterTitle();
352 <    mpf_gr->GetYaxis()->SetTitle("data/mc");
352 >    mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
353      mpf_gr->GetYaxis()->CenterTitle();
354      TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
355      float mpf_a=mpf_pol->GetParameter(0);
# Line 325 | Line 382 | void new_data_mc_agreement_2d() {
382      rabs_gr->Draw("AP");
383      rabs_gr->GetXaxis()->SetTitle("#alpha");
384      rabs_gr->GetXaxis()->CenterTitle();
385 <    rabs_gr->GetYaxis()->SetTitle("data/mc");
385 >    rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
386      rabs_gr->GetYaxis()->CenterTitle();
387      
388      TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
# Line 362 | Line 419 | void new_data_mc_agreement_2d() {
419    MPF_FinalGraph->Draw("AP");
420    MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
421    MPF_FinalGraph->GetXaxis()->CenterTitle();
422 <  MPF_FinalGraph->GetYaxis()->SetTitle("data/mc");
422 >  MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
423    MPF_FinalGraph->GetYaxis()->CenterTitle();
424    MPF_FinalGraph->Draw("AP");
425    
# Line 374 | Line 431 | void new_data_mc_agreement_2d() {
431    rmpf_final_mark->SetTextSize(0.03);
432    rmpf_final_mark->Draw();
433    DrawPrelim();
434 +  float MPFResult=1/mpf_pol0->GetParameter(0);
435 +  float MPFResultError=mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0));
436    
437    stringstream filename2;
438    filename2 << "Extrapolation/FINAL_RESULT_MPF";
# Line 387 | Line 446 | void new_data_mc_agreement_2d() {
446    RABS_FinalGraph->Draw("AP");
447    RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
448    RABS_FinalGraph->GetXaxis()->CenterTitle();
449 <  RABS_FinalGraph->GetYaxis()->SetTitle("data/mc");
449 >  RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
450    RABS_FinalGraph->GetYaxis()->CenterTitle();
451    RABS_FinalGraph->Draw("AP");
452    TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
# Line 400 | Line 459 | void new_data_mc_agreement_2d() {
459    rabs_final_mark->Draw();
460    DrawPrelim();
461    
462 +  float RabsResult=1/rabs_pol0->GetParameter(0);
463 +  float RabsResultError=rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0));
464 +  
465    filename2.str("");
466    filename2 << "Extrapolation/FINAL_RESULT_RABS";
467    CompleteSave(can,filename2.str());
468    
469 +  can->cd();
470 +  TCanvas *FaceValueCanvas = new TCanvas("FaceValueCanvas","FaceValueCanvas");
471 +  FaceValueAtPoint3->SetMarkerStyle(21);
472 +  FaceValueAtPoint3->SetMarkerStyle(21);
473 +  FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
474 +  FaceValueAtPoint3->Draw("AP");
475 +  FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
476 +  FaceValueAtPoint3->GetXaxis()->CenterTitle();
477 +  FaceValueAtPoint3->GetYaxis()->SetTitle("data/mc for #alpha<0.3");
478 +  FaceValueAtPoint3->GetYaxis()->CenterTitle();
479 +  FaceValueAtPoint3->Draw("AP");
480 +  TF1 *faceval_pol0=(TF1*)FaceValueAtPoint3->GetFunction("pol0");
481 +  stringstream faceval_result;
482 +  faceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0)) << "}{ (MPF at #alpha<0.3)}";;
483 +  TText *faceval_final_mark = write_title(faceval_result.str());
484 +  faceval_final_mark->SetX(0.75);
485 +  faceval_final_mark->SetY(0.75);
486 +  faceval_final_mark->SetTextSize(0.03);
487 +  faceval_final_mark->Draw();
488 +  DrawPrelim();
489 +  
490 +  float FaceValResult=1/faceval_pol0->GetParameter(0);
491 +  float FaceValResultError=faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0));
492 +  
493 +  filename2.str("");
494 +  filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
495 +  CompleteSave(can,filename2.str());
496 +  
497 +  cout << "FINAL RESULTS: " << endl;
498 +  cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
499 +  cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
500 +  cout << "Face value at 0.3 :" << FaceValResult << " +/- " << FaceValResultError << endl;
501    
502    delete can;
503    
# Line 411 | Line 505 | void new_data_mc_agreement_2d() {
505   }
506  
507  
508 + void DoPUStudy(string identifier) {
509 +  
510 +  float numVtxcuts[7]={0,10,15,20};
511 +  
512 +  TCanvas *can = new TCanvas("can","can");
513 +  TH1F *malpha[8];
514 +  TH1F *dalpha[8];
515 +  cout << identifier << ";";
516 +  cout << endl;
517 +  for(int i=1;i<5;i++) {
518 +    stringstream sSpecialCut;
519 +    if(i<4) {
520 +      sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
521 +      cout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
522 +    } else {
523 +      sSpecialCut << "numVtx>" << numVtxcuts[i-1];
524 +      cout << numVtxcuts[i-1] << ";";
525 +    }
526 +    
527 +    TCut SpecialCut(sSpecialCut.str().c_str());
528 +    
529 +    
530 +    malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
531 +    dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
532 +    
533 +    
534 +    float a=dalpha[i]->GetMean();
535 +    float b=malpha[i]->GetMean();
536 +    float da=dalpha[i]->GetMeanError();
537 +    float db=malpha[i]->GetMeanError();
538 +    float factor = a / b;
539 +    float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
540 +  
541 +    cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
542 +    malpha[i]->SetLineColor(i+1);
543 +    if(i==0) malpha[i]->Draw("histo");
544 +    else malpha[i]->Draw("histo,same");
545 +  }
546 +  
547 +  CompleteSave(can,"AlphaOverview_"+identifier);
548 +  delete can;
549 +  
550 + /*  for(int i=1;i<8;i++) {
551 +    delete malpha[i];
552 +    delete dalpha[i];
553 +  }
554 +  */
555 +  
556 + }
557 +
558 + void ScenarioComparison() {
559 +  
560 +  TGraphAsymmErrors *gr[5];
561 +  TF1 *fit[5];
562 +  bool dofit=false;
563 +  
564 +  TCanvas *can = new TCanvas("can","can",500,500);
565 +  
566 +  TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
567 +  leg->SetBorderSize(0);
568 +  leg->SetFillColor(kWhite);
569 +  
570 +  for(int i=0;i<5;i++) {
571 +    gr[i] = new TGraphAsymmErrors();
572 +    double x[4]   = {5,12.5,17.5,20};
573 +    double dxl[4] = {5,2.5,2.5,2.5};
574 +    double dxh[4] = {5,2.5,2.5,10};
575 +    
576 +    string name="";
577 +    string color="";
578 +    
579 +    if(i==0) {
580 +      double y[4]   = {1.008,1.013,1.003,1.007};
581 +      double dyl[4] = {0.011,0.009,0.013,0.020};
582 +      double dyh[4] = {0.011,0.009,0.013,0.020};
583 +      name="30/30";
584 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
585 +      color="#a01d99";
586 +    }
587 +    if(i==1) {
588 +      double y[4]   = {1.008+0.001,1.013+0.001,1.003+0.001,1.007+0.001};
589 +      double dyl[4] = {0.011,0.009,0.013,0.020};
590 +      double dyh[4] = {0.011,0.009,0.013,0.020};
591 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
592 +      name="30/10";
593 +      color="#2464bc";
594 +    }
595 +    if(i==2) {
596 +      double y[4]   = {1.014,1.026,1.014,1.059};
597 +      double dyl[4] = {0.011,0.010,0.014,0.023};
598 +      double dyh[4] = {0.011,0.010,0.014,0.023};
599 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
600 +      name="15/15";
601 +      color="#f49230";
602 +    }
603 +    if(i==3) {
604 +      double y[4]   = {1.015,1.020,1.008,1.050};
605 +      double dyl[4] = {0.011,0.010,0.014,0.022};
606 +      double dyh[4] = {0.011,0.010,0.014,0.022};
607 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
608 +      name="20/15";
609 +      color="#c72f3c";
610 +    }
611 +    if(i==4) {
612 +      double y[4]   = {1.015+0.001,1.020+0.001,1.008+0.001,1.050+0.001};
613 +      double dyl[4] = {0.011,0.010,0.014,0.022};
614 +      double dyh[4] = {0.011,0.010,0.014,0.022};
615 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
616 +      name="20/20";
617 +      color="#7cb433";
618 +    }
619 +    
620 +    
621 +    gr[i]->SetTitle();
622 +    gr[i]->GetXaxis()->SetTitle("N(Vtx)");
623 +    gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
624 +    gr[i]->GetXaxis()->CenterTitle();
625 +    gr[i]->GetYaxis()->CenterTitle();
626 +    gr[i]->SetName(name.c_str());
627 +    gr[i]->SetMarkerSize(0.00001);
628 +    gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
629 +    gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
630 +    gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
631 +    if(dofit) {
632 +      gr[i]->Fit("pol1");
633 +      fit[i] = (TF1*)gr[i]->GetFunction("pol1");
634 +      fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
635 +    }
636 +    leg->AddEntry(gr[i],name.c_str(),"l");
637 +    if(i==0) gr[i]->Draw("AP*");
638 +    else gr[i]->Draw("P");
639 +  }
640 +  leg->Draw();
641 +  DrawPrelim();
642 +  
643 +  CompleteSave(can,"ScenarioComparison");
644 +  delete can;
645 +  
646 + }
647 +
648 + void ScenarioComparisonInclusive() {
649 +  
650 +  TGraphAsymmErrors *gr[5];
651 +  TF1 *fit[5];
652 +  bool dofit=true;
653 +  
654 +  TCanvas *can = new TCanvas("can","can",500,500);
655 +  
656 +  TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
657 +  leg->SetBorderSize(0);
658 +  leg->SetFillColor(kWhite);
659 +  
660 +  for(int i=0;i<5;i++) {
661 +    gr[i] = new TGraphAsymmErrors();
662 +    double x[4]   = {5,12.5,17.5,20};
663 +    double dxl[4] = {5,2.5,2.5,2.5};
664 +    double dxh[4] = {5,2.5,2.5,10};
665 +    
666 +    string name="";
667 +    string color="";
668 +    
669 +    if(i==0) {
670 +      double y[4]   = {1.01103,1.01496,1.01486,1.02785};
671 +      double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
672 +      double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
673 +      name="30/30";
674 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
675 +      color="#a01d99";
676 +    }
677 +    if(i==1) {
678 +      double y[4]   = {1.01103+0.001,1.01496+0.001,1.01486+0.001,1.02785+0.001};
679 +      double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
680 +      double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
681 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
682 +      name="30/10";
683 +      color="#2464bc";
684 +    }
685 +    if(i==2) {
686 +      double y[4]   = {1.01575,1.02879,1.0325,1.05456};
687 +      double dyl[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
688 +      double dyh[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
689 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
690 +      name="15/15";
691 +      color="#f49230";
692 +    }
693 +    if(i==3) {
694 +      double y[4]   = {1.01221,1.01922,1.02396,1.04677};
695 +      double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
696 +      double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
697 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
698 +      name="20/15";
699 +      color="#c72f3c";
700 +    }
701 +    if(i==4) {
702 +      double y[4]   = {1.01221+0.001,1.01922+0.001,1.02396+0.001,1.04677+0.001};
703 +      double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
704 +      double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
705 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
706 +      name="20/20";
707 +      color="#7cb433";
708 +    }
709 +    
710 +    
711 +    gr[i]->SetTitle();
712 +    gr[i]->GetXaxis()->SetTitle("N(Vtx)");
713 +    gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
714 +    gr[i]->GetXaxis()->CenterTitle();
715 +    gr[i]->GetYaxis()->CenterTitle();
716 +    gr[i]->SetName(name.c_str());
717 +    gr[i]->SetMarkerSize(0.00001);
718 +    gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
719 +    gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
720 +    gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
721 +    if(dofit) {
722 +      gr[i]->Fit("pol1");
723 +      fit[i] = (TF1*)gr[i]->GetFunction("pol1");
724 +      fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
725 +    }
726 +    leg->AddEntry(gr[i],name.c_str(),"l");
727 +    if(i==0) gr[i]->Draw("AP*");
728 +    else gr[i]->Draw("P");
729 +  }
730 +  leg->Draw();
731 +  DrawPrelim();
732 +  
733 +  CompleteSave(can,"ScenarioComparisonInclusive");
734 +  delete can;
735 +  
736 + }
737 +
738 + void compare_selection(string identifier) {
739 +  bool recognized_scenario=false;
740 +  
741 +  cout << "Running with identifier " << identifier << endl;
742 +  if(identifier=="Zb3010") {
743 +      LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
744 +      EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
745 +      PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
746 +      recognized_scenario=true;
747 +  }
748 +  
749 +  if(identifier=="Zb3010") {
750 +      LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
751 +      EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
752 +      PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
753 +      recognized_scenario=true;
754 +  }
755 +  
756 +  if(identifier=="Zb40") {
757 +      LeadingB=TCut ("Zb40_bTagProbCSVBP[0]>0.898");
758 +      EtaB=TCut("abs(Zb40_pfJetGoodEta[0])<1.3");
759 +      PhiZcut=TCut("abs(Zb40_pfJetDphiZ[0])>2.7");
760 +      recognized_scenario=true;
761 +  }
762 +  
763 +  if(identifier=="Zb301030") {
764 +      LeadingB=TCut ("Zb301030_bTagProbCSVBP[0]>0.898");
765 +      EtaB=TCut("abs(Zb301030_pfJetGoodEta[0])<1.3");
766 +      PhiZcut=TCut("abs(Zb301030_pfJetDphiZ[0])>2.7");
767 +      recognized_scenario=true;
768 +  }
769 +  
770 +  if(identifier=="Zb1530") {
771 +      LeadingB=TCut ("Zb1530_bTagProbCSVBP[0]>0.898");
772 +      EtaB=TCut("abs(Zb1530_pfJetGoodEta[0])<1.3");
773 +      PhiZcut=TCut("abs(Zb1530_pfJetDphiZ[0])>2.7");
774 +      recognized_scenario=true;
775 +  }
776 +  
777 +  if(identifier=="Zb301010") {
778 +      LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
779 +      EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
780 +      PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
781 +      recognized_scenario=true;
782 +  }
783 +  
784 +  if(identifier=="Zb1510") {
785 +      LeadingB=TCut ("Zb1510_bTagProbCSVBP[0]>0.898");
786 +      EtaB=TCut("abs(Zb1510_pfJetGoodEta[0])<1.3");
787 +      PhiZcut=TCut("abs(Zb1510_pfJetDphiZ[0])>2.7");
788 +      recognized_scenario=true;
789 +  }
790 +  
791 +  if(identifier=="Zb301010") {
792 +      LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
793 +      EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
794 +      PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
795 +      recognized_scenario=true;
796 +  }
797 +  
798 +  if(identifier=="ZbMikko") {
799 +      LeadingB=TCut ("ZbMikko_bTagProbCSVBP[0]>0.898 && ZbMikko__IsClean");
800 +      EtaB=TCut("abs(ZbMikko_pfJetGoodEta[0])<1.3 && ZbMikko__IsClean");
801 +      PhiZcut=TCut("abs(ZbMikko_pfJetDphiZ[0])>2.7 && ZbMikko__IsClean");
802 +      recognized_scenario=true;
803 +  }
804 +  
805 +  if(identifier=="Zb3010_SecEta3") {
806 +      LeadingB=TCut ("Zb3010_SecEta3_bTagProbCSVBP[0]>0.898");
807 +      EtaB=TCut("abs(Zb3010_SecEta3_pfJetGoodEta[0])<1.3");
808 +      PhiZcut=TCut("abs(Zb3010_SecEta3_pfJetDphiZ[0])>2.7");
809 +      recognized_scenario=true;
810 +  }
811 +  
812 +  if(identifier=="Zb3010_SecEta5") {
813 +      LeadingB=TCut ("Zb3010_SecEta5_bTagProbCSVBP[0]>0.898");
814 +      EtaB=TCut("abs(Zb3010_SecEta5_pfJetGoodEta[0])<1.3");
815 +      PhiZcut=TCut("abs(Zb3010_SecEta5_pfJetDphiZ[0])>2.7");
816 +      recognized_scenario=true;
817 +  }
818 +  
819 +  if(identifier=="Zb3010_p5Clean") {
820 +      LeadingB=TCut ("Zb3010_p5Clean_bTagProbCSVBP[0]>0.898 && Zb3010_p5Clean_IsClean");
821 +      EtaB=TCut("abs(Zb3010_p5Clean_pfJetGoodEta[0])<1.3 && Zb3010_p5Clean_IsClean");
822 +      PhiZcut=TCut("abs(Zb3010_p5Clean_pfJetDphiZ[0])>2.7 && Zb3010_p5Clean_IsClean");
823 +      recognized_scenario=true;
824 +  }
825 +  
826 +  if(identifier=="Zb3010CHS") {
827 +      LeadingB=TCut ("Zb3010CHS_bTagProbCSVBP[0]>0.898");
828 +      EtaB=TCut("abs(Zb3010CHS_pfJetGoodEta[0])<1.3");
829 +      PhiZcut=TCut("abs(Zb3010CHS_pfJetDphiZ[0])>2.7");
830 +      recognized_scenario=true;
831 +  }
832 +  
833 +  assert(recognized_scenario);
834 +  
835 +  cout << "Cuts: " << endl;
836 +  cout << "   " << (const char*) LeadingB << endl;
837 +  cout << "   " << (const char*) EtaB << endl;
838 +  cout << "   " << (const char*) PhiZcut << endl;
839 +  
840 +  cout << endl << endl;
841 +  
842 + //    ScenarioComparison();
843 + //    ScenarioComparisonInclusive();
844 +  
845 +  
846 + //     DoPUStudy(identifier);
847 +  
848 + //   SpecialCutFlow(identifier);
849 + //    draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
850 + //   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);
851 + //   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);
852 + //   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);
853 +  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);
854 +  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);
855 +  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);
856 +  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);
857 +  draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut,20,0,2,"#alpha","Official/"+identifier+"/alpha_full__"+identifier,1);
858 + //   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);
859 +  
860 + //   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");
861 + //   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");
862 +
863 +  
864 + //   */
865 +  
866 + }
867  
868   void do_basic_ZB_analysis() {
869 +  
870 +  
871    TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
417  write_warning(__FUNCTION__,"TESTING NEW DATA MC AGREEMENT");
872    
873 <  /*
874 < //  print_all_b_yields();
875 < //  draw_mpf_vars();
873 > //     compare_selection("Zb3010_p5Clean");
874 > //     compare_selection("Zb3010CHS");
875 > //     compare_selection("Zb3010");
876 > //     compare_selection("Zb301030");
877 > //     compare_selection("Zb1530");
878 > //     compare_selection("Zb301010");
879 > //     compare_selection("Zb3010_SecEta3");
880 > //     compare_selection("Zb3010_SecEta5");
881 > //     compare_selection("Zb1510");
882 > //     compare_selection("Zb301010");
883 > //     compare_selection("ZbMikko");
884 > //     compare_selection("Zb3010");
885 > //     compare_selection("Zb40");
886 >
887 >  
888 >  print_all_b_yields();
889 >  draw_mpf_vars();
890 > //  DrawEvilCutFlow();
891 >  
892    draw_Zb_kin_vars();
893    
424  data_mc_agreement_2d();
425  */
894    new_data_mc_agreement_2d();
895    
896 +  
897    delete zbcanvas;
898   }
899 +
900 +
901 +
902 +
903 + //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BELOW THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
904 +
905 +
906 + const char* concatenate(string bla1, string bla2) {
907 +  stringstream bla;
908 +  bla << bla1 << bla2;
909 +  return bla.str().c_str();
910 + }
911 +
912 + void SpecialCutFlow(string identifier) {
913 +  stringstream MegaCut;
914 +  
915 +  
916 +  
917 +  MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut(concatenate(identifier,"_pfJetGoodNumBtag>0"))) << ")";
918 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
919 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
920 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
921 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000").c_str())) << ")";
922 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
923 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
924 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
925 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
926 +  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
927 +  MegaCut << " ";
928 +  
929 +  TCanvas *can = new TCanvas("can","can");
930 +  can->SetLogy(1);
931 +  TCut basecut(ZplusBsel);
932 +  
933 +  
934 +  
935 +  TLegend *leg = allsamples.allbglegend();
936 +  
937 +  
938 +  TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
939 +   THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
940 +
941 +  float runningsum=0;
942 +  for(int i=data->GetNbinsX();i>0;i--) {
943 +    runningsum+=data->GetBinContent(i);
944 +    data->SetBinContent(i,runningsum);
945 +  }
946 +  
947 +  float minimum=data->GetMaximum();
948 +  
949 +  TH1F* h;
950 +  TIter nextOF(mcm.GetHists());
951 +  
952 +  float mcsum=0;
953 +  float zb=0;
954 +  while ( h = (TH1F*)nextOF() ) {
955 +    float runningsum=0;
956 +    minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
957 +    for(int i=h->GetNbinsX();i>0;i--) {
958 +      runningsum+=h->GetBinContent(i);
959 +      h->SetBinContent(i,runningsum);
960 +      if(runningsum<minimum&&runningsum>0) minimum=runningsum;
961 +    }
962 +    if(Contains(h->GetName(),"Z_b_")) {
963 +      zb=h->GetBinContent(h->GetNbinsX());
964 +    }
965 +    mcsum+=h->GetBinContent(h->GetNbinsX());
966 +  }
967 +  
968 +  
969 +
970 +  data->SetMinimum(0.05*minimum);
971 +  can->cd(1)->SetBottomMargin(0.25);
972 +  data->GetXaxis()->SetBinLabel(1,"Pre-selection");
973 +  data->GetXaxis()->SetBinLabel(2,"Contains b jet");
974 +  data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
975 +  data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
976 +  data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
977 +  data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
978 +  data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
979 +  data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
980 +  data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
981 +  data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
982 +  data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
983 +  data->GetXaxis()->LabelsOption("v");
984 +  
985 +  stringstream purityinfo;
986 +  purityinfo << "Purity: " << std::setprecision(3) << 100*zb/mcsum << " %";
987 +  stringstream neventsinfo;
988 +  neventsinfo << "Nevents: " << data->GetBinContent(data->GetNbinsX());
989 +  TH1F *crap = new TH1F("crap","",1,0,1);
990 +  crap->SetLineColor(kWhite);
991 +  leg->AddEntry(crap,purityinfo.str().c_str(),"l");
992 +  leg->AddEntry(crap,neventsinfo.str().c_str(),"l");
993 +  
994 +  data->Draw("e1");
995 +  mcm.Draw("same");
996 +  data->Draw("same,e1,axis");
997 +  data->Draw("same,e1");
998 + //  data->Draw("same,TEXT");
999 +  
1000 +  leg->Draw();
1001 +  
1002 +  DrawPrelim();
1003 +  
1004 +  CompleteSave(can,"CutFlow___"+identifier);
1005 +  
1006 +  delete crap;
1007 +  delete data;
1008 +  
1009 +  
1010 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines