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.14 by buchmann, Mon Nov 26 12:06:26 2012 UTC vs.
Revision 1.16 by buchmann, Wed Nov 28 17:51:14 2012 UTC

# Line 30 | Line 30 | namespace ZbData {
30    
31   }
32  
33 + class ZbPointResult {
34 + public:
35 +  ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value Data, Value MC);
36 +  float PtLow;
37 +  float PtHigh;
38 +  float AlphaLow;
39 +  float AlphaHigh;
40 +  string variable;
41 +  Value Data;
42 +  Value MC;
43 +  Value Ratio;
44 + };
45 +
46 + ZbPointResult::ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value data, Value mc) :
47 +  variable(var), PtLow(ptlow), PtHigh(pthigh), AlphaLow(alphalow), AlphaHigh(alphahigh), Data(data), MC(mc) {
48 +    Ratio = data/mc;
49 + }
50 +  
51 + std::ostream &operator<<(std::ostream &ostr, ZbPointResult &b)
52 + {
53 +  ostr << b.PtLow << ";" << b.PtHigh << ";" << b.AlphaLow << ";" << b.AlphaHigh << ";" << b.Data.getValue() << ";" << b.Data.getError() << ";" << b.MC.getValue() << ";" << b.MC.getError() << ";" << b.Ratio.getValue() << ";" << b.Ratio.getError();
54 +  return ostr;
55 + }
56 +
57 +
58 + class ZbResultContainer {
59 + public:
60 +  string name;  
61 +  vector<ZbPointResult> MPF_Rabs_location;
62 +  Value MPF_Result_FaceValue;
63 +  Value Rabs_Result_FaceValue;
64 +  Value MPF_Result_Extrapolated;
65 +  Value Rabs_Result_Extrapolated;
66 +  
67 +  ZbResultContainer(string name);
68 +  ZbResultContainer(const ZbResultContainer &old);
69 + };
70 +
71 +
72 + ZbResultContainer::ZbResultContainer(string _name) {
73 +  name=_name;
74 +  cout<< "Result Container for \"" << name << "\" has been initialized" << endl;
75 +  
76 + }
77 +
78 + std::ostream &operator<<(std::ostream &ostr, ZbResultContainer &b)
79 + {
80 +  ostr << "Results for " << b.name << endl;
81 +  ostr << " MPF and Rabs locations: " << endl;
82 +  for(int i=0;i<b.MPF_Rabs_location.size();i++) ostr << "   " << b.MPF_Rabs_location[i] << endl;
83 +  ostr << " Results: " << endl;
84 +  ostr << "   Face Value: " << endl;
85 +  ostr << "     MPF " << b.MPF_Result_FaceValue << endl;
86 +  ostr << "     Rabs " << b.Rabs_Result_FaceValue << endl;
87 +  ostr << "   Extrapolated: " << endl;
88 +  ostr << "     MPF " << b.MPF_Result_Extrapolated << endl;
89 +  ostr << "             Rabs " << b.Rabs_Result_Extrapolated << endl;
90 +  
91 +  return ostr;
92 + }
93 +
94 +
95 +
96   using namespace ZbData;
97  
98  
# Line 101 | Line 164 | void draw_kin_variable(string variable,
164    delete ckin;
165   }
166    
104 – Value getfrom2Dmap(TH2F *map, int ixbin, int iybin) {
105 –  Value sum(0,0);
106 –  for(int iy=1;iy<=iybin;iy++) {
107 –    sum=sum+Value(map->GetBinContent(ixbin,iy),map->GetBinError(ixbin,iy));
108 –  }
109 –  return sum;
110 – }
111 –  
167   void print_all_b_yields() {
168    cout << "Basic selection with a b jet" << endl;
169    print_yield(ZplusBsel&&"Zb3010_pfJetGoodNumBtag>0");
# Line 134 | Line 189 | void print_all_b_yields() {
189   //  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.05");
190   }  
191  
192 < void DrawEvilCutFlow() {
192 > void TEST_DrawEvilCutFlow() {
193    stringstream MegaCut;
194    
195    MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut("Zb3010_pfJetGoodNumBtag>0")) << ")";
# Line 267 | Line 322 | void draw_mpf_vars() {
322      return TMath::Median(numBins, &x, &y);
323   }*/
324  
325 < Value get_Zb_data_over_mc(string variable, string variable2, TCut cut, string saveas) {
325 > Value get_Zb_data_over_mc(string ContainerName, ZbResultContainer &res, string variable, string variable2, TCut cut, string saveas, float pt1, float pt2, float a1, float a2) {
326   //  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
327    TCanvas *cn = new TCanvas("cn","cn");
328    string varname="MPF";
# Line 276 | Line 331 | Value get_Zb_data_over_mc(string variabl
331    TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
332    hmc->SetLineColor(kBlue);
333    
279 –  
280 –  
281 –  
334    float a=hdata->GetMean();
335    float b=hmc->GetMean();
336    float da=hdata->GetMeanError();
# Line 311 | Line 363 | Value get_Zb_data_over_mc(string variabl
363    gaus_data->SetLineColor(kRed);
364    
365    hdata->GetXaxis()->SetRangeUser(0,2);
366 +  hdata->GetYaxis()->SetTitleOffset(1.3);
367    hdata->Draw("e1");
368    hmc->Draw("histo,same");
369    hdata->Draw("e1,same");
370    gaus_data->Draw("same");
371    gaus_mc->Draw("same");
372    stringstream summary;
373 <  summary << "#splitline{" << varname << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{";
373 >  summary << "#splitline{#bf{" << varname << "}}{#splitline{" << pt1 << " < p_{T}^{Z} < " << pt2 << ", " << a1 << " #leq #alpha < " << a2 << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{";
374    summary << "#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{#splitline{ratio: " << factor;
375    summary << " #pm " << error << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
376    summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
377    summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "+/-" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
378 <  summary << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) << "+/-" << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) * sqrt((gaus_data->GetParError(1) * gaus_data->GetParError(1))/(gaus_data->GetParameter(1)*gaus_data->GetParameter(1)) + (gaus_mc->GetParError(1) * gaus_mc->GetParError(1))/(gaus_mc->GetParameter(1)*gaus_mc->GetParameter(1)) ) << "}}}}}}}";
378 >  summary << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) << "+/-" << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) * sqrt((gaus_data->GetParError(1) * gaus_data->GetParError(1))/(gaus_data->GetParameter(1)*gaus_data->GetParameter(1)) + (gaus_mc->GetParError(1) * gaus_mc->GetParError(1))/(gaus_mc->GetParameter(1)*gaus_mc->GetParameter(1)) ) << "}}}}}}}}";
379 >  
380    TText *infobox = write_title(summary.str());
381    infobox->SetX(0.75);
382    infobox->SetTextSize(0.03);
# Line 331 | Line 385 | Value get_Zb_data_over_mc(string variabl
385      
386      
387    DrawPrelim();
388 <  CompleteSave(cn,"ResponsePlots/"+saveas);
388 >  CompleteSave(cn,ContainerName+"/ResponsePlots/"+saveas);
389    
390 +  res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
391    cout << "Data;" << a << ";"<<da<<";MC"<<b<<";"<<db<<endl;
392    delete cn;
393    delete hdata;
394    delete hmc;
395  
341 – //    cout << a << "+/-" << da << "        " << b << "=/-" << db << endl;
342 –
396    return Value(factor,error);
397   }
398    
399 + ZbResultContainer new_data_mc_agreement_2d(string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt) {
400 +  ZbResultContainer results(ContainerName);
401    
402 <
403 < void new_data_mc_agreement_2d() {
402 >  cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)))")*BTagWgt;
403 >  LeadingB=BTagCut;
404    
405    gStyle->SetOptFit(0);
406    TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
407 <  const int nalphacuts=6;
408 <  float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
407 > //  const int nalphacuts=6;
408 > //  float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
409 >  
410 > //  const int nptcuts=5;
411 > //  float ptcuts[nptcuts]={30,50,75,125,1000};
412 >  
413 >  write_info(__FUNCTION__,"JUST TESTING - PLEASE REMOVE THIS LINE");const int nalphacuts=2;float alphacuts[nalphacuts] = {0.1,0.3};const int nptcuts=3;float ptcuts[nptcuts]={30,75,1000};
414    
355 –  const int nptcuts=5;
356 –  float ptcuts[nptcuts]={30,50,75,125,1000};
415    
416    float MPF_Results[nalphacuts][nptcuts];
417    float MPF_Errors[nalphacuts][nptcuts];
# Line 364 | Line 422 | void new_data_mc_agreement_2d() {
422    for(int ia=0;ia<nalphacuts;ia++) {
423      for(int ipt=0;ipt<nptcuts-1;ipt++) {
424        stringstream specialcut;
425 <      if(ia>0) specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<" << alphacuts[ia] << " && Zb3010_alpha>" << alphacuts[ia-1] << "))";
426 <      else     specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<" << alphacuts[ia] << "))";
425 >      float alow,ahigh;
426 >      if(ia>0) {
427 >        specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << " && " << AlphaVariable << ">" << alphacuts[ia-1] << "))";
428 >        alow=alphacuts[ia-1];
429 >        ahigh=alphacuts[ia];
430 >      } else {
431 >        specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << "))";
432 >        alow=0;
433 >        ahigh=alphacuts[ia];
434 >      }
435        cout << specialcut.str() << endl;
436 <      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]));
437 <      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]));
436 >      Value MPF_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
437 >      Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
438        
439        MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
440        MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
# Line 425 | Line 491 | void new_data_mc_agreement_2d() {
491      }
492      
493      stringstream specialcut;
494 <    specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<0.3))";
495 <    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])+"__INCLUSIVEalpha_0.3");
496 <    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])+"__INCLUSIVEalpha_0.3");
494 >    specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<0.3))";
495 >    Value MPF_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3",ptcuts[ipt],ptcuts[ipt+1],0,0.3);
496 >    Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3",ptcuts[ipt],ptcuts[ipt+1],0,0.3);
497      cout << "About to add a point in pt (" << (ptcuts[ipt]+ptcuts[ipt+1])/2 << " ) to face value : " << MPF_data_over_mc << " (MPF) , " <<  RABS_data_over_mc << "(RABS)" << endl;
498      MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
499      MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
# Line 456 | Line 522 | void new_data_mc_agreement_2d() {
522      histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
523      histo->GetXaxis()->CenterTitle();
524      histo->GetYaxis()->CenterTitle();
525 +    histo->GetYaxis()->SetTitleOffset(1.3);
526      histo->SetStats(0);
527      histo->GetXaxis()->SetRangeUser(0,0.4);
528      histo->GetYaxis()->SetRangeUser(min,max);
# Line 490 | Line 557 | void new_data_mc_agreement_2d() {
557      string filenamebkp=filename.str();
558      filename << "__MPF";
559      DrawPrelim();
560 <    CompleteSave(can,filename.str());
560 >    CompleteSave(can,ContainerName+"/"+filename.str());
561      cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
562      
563      filename.str("");
# Line 544 | Line 611 | void new_data_mc_agreement_2d() {
611      RABS_ExtraPolatedResults[ipt]=rabs_a;
612      filename << filenamebkp << "__RABS";
613      DrawPrelim();
614 <    CompleteSave(can,filename.str());
614 >    CompleteSave(can,ContainerName+"/"+filename.str());
615      cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
616    }
617    
# Line 578 | Line 645 | void new_data_mc_agreement_2d() {
645    
646    stringstream filename2;
647    filename2 << "Extrapolation/FINAL_RESULT_MPF";
648 <  CompleteSave(can,filename2.str());
648 >  CompleteSave(can,ContainerName+"/"+filename2.str());
649    
650    
651    RABS_FinalGraph->SetMarkerStyle(21);
# Line 606 | Line 673 | void new_data_mc_agreement_2d() {
673    
674    filename2.str("");
675    filename2 << "Extrapolation/FINAL_RESULT_RABS";
676 <  CompleteSave(can,filename2.str());
676 >  CompleteSave(can,ContainerName+"/"+filename2.str());
677    
678    can->cd();
679    MPF_FaceValueAtPoint3->SetMarkerStyle(21);
# Line 638 | Line 705 | void new_data_mc_agreement_2d() {
705    
706    filename2.str("");
707    filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
708 <  CompleteSave(can,filename2.str());
708 >  CompleteSave(can,ContainerName+"/"+filename2.str());
709    
710    can->cd();
711    RABS_FaceValueAtPoint3->SetMarkerStyle(21);
# Line 670 | Line 737 | void new_data_mc_agreement_2d() {
737    
738    filename2.str("");
739    filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
740 <  CompleteSave(can,filename2.str());
740 >  CompleteSave(can,ContainerName+"/"+filename2.str());
741    
742    cout << "FINAL RESULTS: " << endl;
743    cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
# Line 680 | Line 747 | void new_data_mc_agreement_2d() {
747    
748    delete can;
749    
750 <  
750 >  results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
751 >  results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
752 >  results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
753 >  results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
754 >  cout << results << endl;
755 >  return results;
756   }
757  
758  
759 < void DoPUStudy(string identifier) {
759 > void TEST_DoPUStudy(string identifier) {
760    
761    float numVtxcuts[7]={0,10,15,20};
762    
# Line 734 | Line 806 | void DoPUStudy(string identifier) {
806    
807   }
808  
809 < void ScenarioComparison() {
809 > void TEST_ScenarioComparison() {
810   //very dumb way to do this.
811    TGraphAsymmErrors *gr[5];
812    TF1 *fit[5];
# Line 824 | Line 896 | void ScenarioComparison() {
896    
897   }
898  
899 < void ScenarioComparisonInclusive() {
899 > void TEST_ScenarioComparisonInclusive() {
900    //dumbest way ever to do this. but ok we only need to do it once.
901    TGraphAsymmErrors *gr[5];
902    TF1 *fit[5];
# Line 914 | Line 986 | void ScenarioComparisonInclusive() {
986    
987   }
988  
989 < void compare_selection(string identifier) {
989 > void TEST_compare_selection(string identifier) {
990    bool recognized_scenario=false;
991    
992    cout << "Running with identifier " << identifier << endl;
# Line 1068 | Line 1140 | void GetNumberEventsInsideOutsideAlphaWi
1140    delete mc_OUT;
1141   }
1142  
1143 < void GetNumberEventsInsideOutsideAlphaWindow() {
1143 > void TEST_GetNumberEventsInsideOutsideAlphaWindow() {
1144    TCanvas *ca = new TCanvas("ca","ca");
1145   //   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1146   //   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
# Line 1107 | Line 1179 | void GetNumberEventsInsideOutsideAlphaWi
1179  
1180  
1181      
1182 < void do_basic_ZB_analysis(int do_inclusive) {
1182 > void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
1183    
1184    //https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagPOG#2011_Data_and_MC
1185    //https://twiki.cern.ch/twiki/pub/CMS/BtagPOG/SFb-ttbar_payload.txt
1186    
1187    cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1116 –  if(do_inclusive) {
1117 –    LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>=-10000");
1118 –    cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)))");
1119 –    write_info(__FUNCTION__,"Inclusive selection - no btag requirement, no btag efficiency/mistag correction (weight)");
1120 –  } else {
1121 –    LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
1122 –    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))");
1123 –    write_info(__FUNCTION__,"Exclusive selection - switched on btag requirement and btag efficiency/mistag correction (weight)");
1124 –  }
1125 –  
1188    bool doquick=true;
1189    
1190    TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
# Line 1151 | Line 1213 | void do_basic_ZB_analysis(int do_inclusi
1213     }
1214    
1215   //  GetNumberEventsInsideOutsideAlphaWindow();
1216 <  new_data_mc_agreement_2d();
1216 >  
1217 >  dout <<"THIS IS WHERE IT GETS REALLY INTERESTING!!!!" << endl;
1218 >  if(DoCompleteAnalysis) {
1219 >    /* need to run the following scenarios:
1220 >     * - normal
1221 >     * - inclusive
1222 >     * - medium WP
1223 >     * - alpha up
1224 >     * - alpha down
1225 >     */
1226 >    
1227 >    ZbResultContainer inclusive = new_data_mc_agreement_2d("Zb3010_alpha","inclusive",TCut("Zb3010_bTagProbCSVBP[0]>-10000"),TCut("1.0"));
1228 >    ZbResultContainer Reference = new_data_mc_agreement_2d("Zb3010_alpha","reference",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1229 >    ZbResultContainer effmistagUp = new_data_mc_agreement_2d("Zb3010_alpha","effmistagUp",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp"));
1230 >    ZbResultContainer effmistagDown = new_data_mc_agreement_2d("Zb3010_alpha","effmistagDown",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTDown"));
1231 >    ZbResultContainer medium = new_data_mc_agreement_2d("Zb3010_alpha","medium",TCut("Zb3010_bTagProbCSVBP[0]>0.679"),TCut("ZbCHS3010_BTagWgtM"));
1232 >    ZbResultContainer loose = new_data_mc_agreement_2d("Zb3010_alpha","loose",TCut("Zb3010_bTagProbCSVBP[0]>0.244"),TCut("ZbCHS3010_BTagWgtL"));
1233 >    ZbResultContainer AlphaUp = new_data_mc_agreement_2d("Zb3010_alphaUp","AlphaUp",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1234 >    ZbResultContainer AlphaDown = new_data_mc_agreement_2d("Zb3010_alphaDown","AlphaDown",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1235 >    
1236 >    //and now let's compute some systematics!
1237 >    //1 alpha variation
1238 >    float SysMPF_Alpha_down = Reference.MPF_Result_FaceValue.getValue()-AlphaDown.MPF_Result_FaceValue.getValue();
1239 >    float SysMPF_Alpha_up = Reference.MPF_Result_FaceValue.getValue()-AlphaUp.MPF_Result_FaceValue.getValue();
1240 >    
1241 >    //2 btagging efficiency/mistag correction
1242 >    float SysMPF_EFF_down = Reference.MPF_Result_FaceValue.getValue()-effmistagDown.MPF_Result_FaceValue.getValue();
1243 >    float SysMPF_EFF_up = Reference.MPF_Result_FaceValue.getValue()-effmistagUp.MPF_Result_FaceValue.getValue();
1244 >    
1245 >    //3 purity
1246 >    zbcanvas->cd();
1247 >    write_info(__FUNCTION__,"Don't know the exact purity at the moment, still needs to be determined!");
1248 >    TGraphErrors *gSysMPF_Purity = new TGraphErrors();
1249 >    gSysMPF_Purity->SetPoint(0,0.05,inclusive.MPF_Result_FaceValue.getValue());
1250 >    gSysMPF_Purity->SetPointError(0,0.0,inclusive.MPF_Result_FaceValue.getError());
1251 >    gSysMPF_Purity->SetPointError(1,0.20,loose.MPF_Result_FaceValue.getError());//x value?
1252 >    gSysMPF_Purity->SetPointError(1,0.0,loose.MPF_Result_FaceValue.getError());
1253 >    gSysMPF_Purity->SetPointError(2,0.54,medium.MPF_Result_FaceValue.getError());//x value?
1254 >    gSysMPF_Purity->SetPointError(2,0.0,medium.MPF_Result_FaceValue.getError());
1255 >    gSysMPF_Purity->SetPointError(3,0.82,Reference.MPF_Result_FaceValue.getError());//x value?
1256 >    gSysMPF_Purity->SetPointError(3,0.0,Reference.MPF_Result_FaceValue.getError());
1257 >    
1258 >    gSysMPF_Purity->GetXaxis()->SetTitle("Purity");
1259 >    gSysMPF_Purity->GetXaxis()->CenterTitle();
1260 >    gSysMPF_Purity->GetYaxis()->SetTitle("C_{abs}");
1261 >    gSysMPF_Purity->GetYaxis()->CenterTitle();
1262 >    
1263 >    gSysMPF_Purity->Fit("pol1");
1264 >    gSysMPF_Purity->Draw();
1265 >    CompleteSave(zbcanvas,"Systematics/Purity");
1266 >    
1267 >    TF1 *PurityFit = (TF1*)gSysMPF_Purity->GetFunction("pol1");
1268 >    float ExtrapolatedResult=PurityFit->GetParameter(0)+1.0*PurityFit->GetParameter(1);
1269 >    
1270 >    float SysMPF_Purity = abs(ExtrapolatedResult-Reference.MPF_Result_FaceValue.getValue());
1271 >    
1272 >    float sys_alpha  = abs(SysMPF_Alpha_down)>abs(SysMPF_Alpha_up)?abs(SysMPF_Alpha_down):abs(SysMPF_Alpha_up);
1273 >    float sys_eff    = abs(SysMPF_EFF_down)>abs(SysMPF_EFF_up)?abs(SysMPF_EFF_down):abs(SysMPF_EFF_up);
1274 >    float sys_purity = abs(SysMPF_Purity);
1275 >    float sys        = sqrt(sys_alpha*sys_alpha+sys_eff*sys_eff+sys_purity*sys_purity);
1276 >    
1277 >    dout << "MPF METHOD:" << endl;
1278 >    dout << "   Systematics : " << endl;
1279 >    dout << "      Alpha variation      : " << SysMPF_Alpha_down << "   ,  " << SysMPF_Alpha_up << "   ( " << sys_alpha << ") " << endl;
1280 >    dout << "      eff/mistag variation : " << SysMPF_EFF_down << "   ,  " << SysMPF_EFF_up << "   ( " << sys_eff << ") " << endl;
1281 >    dout << "      purity               : " << SysMPF_Purity << "   ( " << sys_purity << " )" << endl;
1282 >    dout << "      total                : " << sys << endl;
1283 >    dout << endl << endl;
1284 >    dout << "   FINAL RESULT : " << endl;
1285 >    dout << "     " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << sys << " (sys) " << endl;
1286 >    
1287 >    LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
1288 >    
1289 >  } else {
1290 >    ZbResultContainer Reference = new_data_mc_agreement_2d("Zb3010_alpha","reference",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgt"));
1291 >  }
1292 >  
1293    
1294    delete zbcanvas;
1295   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines