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

Comparing UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C (file contents):
Revision 1.18 by buchmann, Thu Jun 13 15:46:06 2013 UTC vs.
Revision 1.27 by buchmann, Wed Jun 26 11:08:31 2013 UTC

# Line 1 | Line 1
1   #include <iostream>
2  
3 +
4 +
5   #include <TVirtualIndex.h>
6  
7   #include <RooRealVar.h>
# Line 9 | Line 11
11   #include <RooCategory.h>
12  
13   #include <RooPlot.h>
14 + #include <RooGaussian.h>
15 + #include <RooConstVar.h>
16   #include <RooSimultaneous.h>
17   #include <RooAddPdf.h>
18   #include <RooFitResult.h>
# Line 25 | Line 29
29   #include "RooStats/HypoTestInverterOriginal.h"
30  
31   //#include "ParametrizedEdge.C"
32 + #include "ExtendedMath.C"
33   #include "EdgeModules/RooSUSYTPdf.cxx"
34   #include "EdgeModules/RooSUSYBkgPdf.cxx"
35 + #include "EdgeModules/RooSUSYCompleteBkgPdf.cxx"
36  
37  
38 + #include "md5/md5.h"
39 + #include "md5/md5.C"
40 +
41   using namespace std;
42   using namespace PlottingSetup;
43  
# Line 72 | Line 81 | namespace EdgeFitter {
81    TCut cut;
82    
83    RooDataSet* AllData;
84 <  RooDataSet* SFSample;
85 <  RooDataSet* OFSample;
84 >  RooDataSet* eeSample;
85 >  RooDataSet* mmSample;
86 >  RooDataSet* OFSample;
87 >  RooDataSet* SFSample;
88    
89    bool MarcoDebug=true;
90    
91    float FixedMEdge=-1;
92 <  float FixedMEdgeChi2=-1;
92 >  float FixedMEdgeChi2_H1=-1;
93 >  float FixedMEdgeChi2_H0=-1;
94    
95    bool RejectPointIfNoConvergence=false;
96    
97    string Mode="UndefinedMode";
98    
99 < }
88 <
89 < TLatex* WriteAppleStudel() {
90 <  string sel="Apple Strudel Preliminary";
91 <  TLatex *sele = new TLatex(0.97,0.135,sel.c_str());
92 <  sele->SetNDC(true);
93 <  sele->SetTextColor(TColor::GetColor("#848484"));
94 <  sele->SetTextFont(42);
95 <  sele->SetTextAlign(32);
96 <  sele->SetTextSize(0.03);
97 <  sele->SetTextAngle(270);
98 <  return sele;
99 >  bool AllowTriangle=true;
100   }
101  
102   TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
# Line 277 | Line 278 | float EdgeFitter::calcExclusion(RooWorks
278    
279   }
280  
281 + string GetSkimName(int isample) {
282 +  return removefunnystring(allsamples.collection[isample].filename);
283 + }
284 +
285   TTree* SkimTree(int isample) {
286 <  string TreeName = allsamples.collection[isample].filename;
287 <  TTree* newTree = new TTree("NanoTree",TreeName.c_str());
286 >  string TreeName = GetSkimName(isample);
287 >  cout << "About to skim " << TreeName << " for sample id " << isample << endl;
288 >  TTree* newTree = new TTree(TreeName.c_str(),TreeName.c_str());
289    
290    float mll,edgeWeight;
291    int id1,id2;
# Line 288 | Line 294 | TTree* SkimTree(int isample) {
294    newTree->Branch("mll",&mll,"mll/F");
295    newTree->Branch("id1",&id1,"id1/I");
296    newTree->Branch("id2",&id2,"id2/I");
297 +  newTree->Branch("isample",&isample,"isample/I");
298    
299    float xsweight=1.0;
300    if(allsamples.collection[isample].is_data==false) xsweight=luminosity*allsamples.collection[isample].weight;
301    if(EdgeFitter::MarcoDebug) {
295    dout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
302      dout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
303 +    dout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
304    }
305    
306    float tmll;
# Line 319 | Line 326 | TTree* SkimTree(int isample) {
326     }
327    }
328    
329 <  if(EdgeFitter::MarcoDebug) dout << "     Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
329 >  if(EdgeFitter::MarcoDebug) dout << "   Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
330    return newTree;
331   }
332  
# Line 331 | Line 338 | void EdgeFitter::InitializeVariables(flo
338   }
339  
340   TTree* MergeTrees(vector<TTree*> trees) {
341 +  assert(trees.size()>0);
342    TTree * newtree = (TTree*)trees[0]->CloneTree(0);
343 +  newtree->SetTitle("FullTree");
344 +  newtree->SetName("FullTree");
345    trees[0]->GetListOfClones()->Remove(newtree);
346    trees[0]->ResetBranchAddresses();
347    newtree->ResetBranchAddresses();
# Line 358 | Line 368 | TTree* MergeTrees(vector<TTree*> trees)
368  
369   void EdgeFitter::PrepareDatasets(int is_data) {
370    write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
371 <  TFile *tempout = new TFile("tempout.root","RECREATE");
371 >  
372 >  ensure_directory_exists(PlottingSetup::cbafbasedir+"/EdgeCache/");
373 >  
374 >  stringstream FileName;
375 >  FileName << PlottingSetup::cbafbasedir << "/EdgeCache/" << md5((const char*)cut) << ".root";
376 >  
377 >  
378 >  TFile *fEdgeCache;
379 >  if(PlottingSetup::do_CleanCache) fEdgeCache = new TFile(FileName.str().c_str(),"RECREATE");
380 >  else fEdgeCache = new TFile(FileName.str().c_str(),"UPDATE");
381 >  
382 >  if(fEdgeCache->GetNkeys()==0) {
383 >    ofstream CacheOverview;
384 >    CacheOverview.open((PlottingSetup::cbafbasedir+"/EdgeCache/CacheOverview").c_str(),ios::app);
385 >    CacheOverview << md5((const char*)cut) << ";" << (const char*) cut << endl;
386 >    CacheOverview.close();
387 >  }
388 >  
389    vector<TTree*> SkimmedTrees;
390    TTree *SkimmedTree[(int)allsamples.collection.size()];
391    for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
# Line 368 | Line 395 | void EdgeFitter::PrepareDatasets(int is_
395      if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
396      if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
397      if(EdgeFitter::MarcoDebug) dout << "Considering : " << allsamples.collection[isample].samplename << endl;
398 <    SkimmedTree[isample] = SkimTree(isample);
399 <    tempout->cd();
400 <    SkimmedTree[isample]->Write();
398 >    
399 >    string SkimName = GetSkimName(isample);
400 >    SkimmedTree[isample] = (TTree*)fEdgeCache->Get(SkimName.c_str());
401 >    if(!SkimmedTree[isample]) {
402 >      dout << "Need to generate tree for " << allsamples.collection[isample].filename << endl;
403 >      SkimmedTree[isample] = SkimTree(isample);
404 >      fEdgeCache->cd();
405 >      SkimmedTree[isample]->Write();
406 >    } else {
407 >      dout << "Loaded tree for " << allsamples.collection[isample].filename << " from cache file " << endl;
408 >    }
409 >      
410      SkimmedTrees.push_back(SkimmedTree[isample]);
411    }
412    
# Line 387 | Line 423 | void EdgeFitter::PrepareDatasets(int is_
423    string title="CMS Data";
424    if(is_data!=1) title="CMS SIMULATION";
425    RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
426 <  tempout->cd();
427 <  completetree->Write();
428 <  tempout->Close();
429 <
430 <  
395 <  EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
426 >  fEdgeCache->Close();
427 >  dout << "Edge cache closed." << endl;
428 >    
429 >  EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==0");
430 >  EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==1");
431    EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
432 +  EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
433    EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
434    
435 <  SFSample->SetName("SFSample");
435 >  
436 >  eeSample->SetName("eeSample");
437 >  mmSample->SetName("mmSample");
438    OFSample->SetName("OFSample");
439 +  SFSample->SetName("SFSample");
440    AllData->SetName("AllData");
441    
442 +  eeSample->SetTitle("eeSample");
443 +  mmSample->SetTitle("mmSample");
444 +  OFSample->SetTitle("OFSample");
445 +  SFSample->SetTitle("SFSample");
446 +  AllData->SetTitle("AllData");
447 +  
448    if(EdgeFitter::MarcoDebug) {
449 <    dout << "Number of events in data sample = " << AllData->sumEntries() << endl;
450 <    dout << "Number of events in eemm sample = " << SFSample->sumEntries() << endl;
451 <    dout << "Number of events in em sample = " << OFSample->sumEntries() << endl;
449 >    dout << "Number of (weighted) events in full sample = " << AllData->sumEntries() << " (raw events : " << AllData->numEntries() << ")" <<  endl;
450 >    dout << "Number of (weighted) events in ee sample = " << eeSample->sumEntries() << " (raw events : " << eeSample->numEntries() << ")" <<  endl;
451 >    dout << "Number of (weighted) events in mm sample = " << mmSample->sumEntries() << " (raw events : " << mmSample->numEntries() << ")" <<  endl;
452 >    dout << "Number of (weighted) events in SF sample = " << SFSample->sumEntries() << " (raw events : " << SFSample->numEntries() << ")" <<  endl;
453 >    dout << "Number of (weighted) events in em sample = " << OFSample->sumEntries() << " (raw events : " << OFSample->numEntries() << ")" <<  endl;
454    }
455    
456   }
457  
458   void EdgeFitter::DrawDatasetContent(int is_data) {
459 +  TCanvas* ceedata = new TCanvas("ceedata","ceedata") ;
460 +  
461    RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
462 <  RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("SF sample")) ;
462 >  RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("ee sample")) ;
463    frame1->GetXaxis()->CenterTitle(1);
464    frame1->GetYaxis()->CenterTitle(1);
465 <  SFSample->plotOn(frame1,RooFit::Name("SFdata")) ;
465 >  eeSample->plotOn(frame1,RooFit::Name("eedata"));
466    
467 <  RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
467 >  RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("mm sample")) ;
468    frame2->GetXaxis()->CenterTitle(1);
469    frame2->GetYaxis()->CenterTitle(1);
470 <  OFSample->plotOn(frame2,RooFit::Name("OFdata")) ;
470 >  mmSample->plotOn(frame2,RooFit::Name("mmdata"));
471    
472 <  TCanvas* cSFdata = new TCanvas("cSFdata","cSFdata") ;
473 <  cSFdata->cd() ;
472 >  RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
473 >  frame3->GetXaxis()->CenterTitle(1);
474 >  frame3->GetYaxis()->CenterTitle(1);
475 >  OFSample->plotOn(frame3,RooFit::Name("OFdata"));
476 >  
477 >  TH1F *ee_ref = allsamples.Draw("ee_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1==id2&&id1==0"),is_data, luminosity);
478 >  TH1F *mm_ref = allsamples.Draw("mm_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1==id2&&id1==1"),is_data, luminosity);
479 >  TH1F *of_ref = allsamples.Draw("of_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1!=id2"),is_data, luminosity);
480 >  
481 >  ee_ref->SetFillColor(TColor::GetColor("#3ADF00"));
482 >  mm_ref->SetFillColor(TColor::GetColor("#3ADF00"));
483 >  of_ref->SetFillColor(TColor::GetColor("#3ADF00"));
484 >  
485 >  TH1F *ee_fit = (TH1F*)eeSample->createHistogram("ee_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
486 >  TH1F *mm_fit = (TH1F*)mmSample->createHistogram("mm_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
487 >  TH1F *of_fit = (TH1F*)OFSample->createHistogram("of_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
488 >  
489 >  TLegend *leg = make_legend();
490 >  leg->AddEntry(ee_fit,"RooDataSet content","p");
491 >  leg->AddEntry(ee_ref,"CBAF cross-check","f");
492 >  leg->SetX1(0.5);
493 >  leg->SetY1(0.75);
494 >
495 >  ceedata->cd() ;
496    gPad->SetLeftMargin(0.15);
497 +  frame1->SetMaximum(1.2*frame1->GetMaximum());
498    frame1->GetYaxis()->SetTitleOffset(1.4);
499    frame1->Draw();
500 +  ee_ref->Draw("histo,same");
501 +  ee_fit->Draw("e1,same");
502 +  leg->Draw("same");
503    if(is_data==data) DrawPrelim();
504    else DrawPrelim(PlottingSetup::luminosity,true);
505 <  CompleteSave(cSFdata,"Edge/SF_NoFit");
505 >  CompleteSave(ceedata,"Edge/PreFit_ee");
506    
507 <  TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
433 <  cOFdata->cd() ;
507 >  ceedata->cd() ;
508    gPad->SetLeftMargin(0.15);
509 <  frame2->SetMaximum(frame1->GetMaximum());
509 >  frame2->SetMaximum(1.2*frame2->GetMaximum());
510    frame2->GetYaxis()->SetTitleOffset(1.4);
511    frame2->Draw();
512 +  mm_ref->Draw("histo,same");
513 +  mm_fit->Draw("e1,same");
514 +  leg->Draw("same");
515 +  if(is_data==data) DrawPrelim();
516 +  else DrawPrelim(PlottingSetup::luminosity,true);
517 +  CompleteSave(ceedata,"Edge/PreFit_mm");
518 +  
519 +  TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
520 +  cOFdata->cd() ;
521 +  gPad->SetLeftMargin(0.15);
522 +  frame3->SetMaximum(1.2*frame3->GetMaximum());
523 +  frame3->GetYaxis()->SetTitleOffset(1.4);
524 +  frame3->Draw();
525 +  of_ref->Draw("histo,same");
526 +  of_fit->Draw("e1,same");
527 +  leg->Draw("same");
528    if(is_data==data) DrawPrelim();
529    else DrawPrelim(PlottingSetup::luminosity,true);
530 <  CompleteSave(cOFdata,"Edge/OF_NoFit");
530 >  CompleteSave(cOFdata,"Edge/PreFit_OF");
531 >
532 >  delete ee_fit;
533 >  delete mm_fit;
534 >  delete of_fit;
535 >  
536 >  delete ee_ref;
537 >  delete mm_ref;
538 >  delete of_ref;
539   }
540    
541   string WriteWithError(float central, float error, int digits) {
# Line 478 | Line 576 | Yield EdgeFitter::Get_T_estimate(float j
576  
577   void EdgeFitter::DoFit(int is_data, float jzb_cut) {
578    
481  TLatex *ASP = WriteAppleStudel();
482  
579    stringstream prefix;
580    if(is_data==data) prefix << "data_";
581    if(is_data==mc) prefix << "mc_";
# Line 487 | Line 583 | void EdgeFitter::DoFit(int is_data, floa
583    
584    prefix << EdgeFitter::Mode << "_" << jzb_cut;
585    
586 <
586 >  if(!EdgeFitter::AllowTriangle && EdgeFitter::FixedMEdge>=0) prefix << "__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge);
587 >  
588 >  if(EdgeFitter::AllowTriangle) prefix << "__H1";
589 >  else prefix << "__H0";
590    
591    RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
592    RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
593    RooCategory sample("sample","sample") ;
594 <  sample.defineType("SF");
595 <  //sample.defineType("mm");
594 >  sample.defineType("ee");
595 >  sample.defineType("mm");
596    sample.defineType("OF");
597 +  sample.defineType("SF");
598 +  
599 +  RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::Import("mm",*mmSample),RooFit::Import("OF",*OFSample),RooFit::Import("SF",*SFSample),RooFit::WeightVar(edgeWeight));
600 +  
601 +  //------------------------------------------------------------------------------------------------------------------------------------
602 +  //    _____ _               __
603 +  //   / ____| |             /_ |
604 +  //  | (___ | |_ ___ _ __    | |
605 +  //   \___ \| __/ _ \ '_ \   | |
606 +  //   ____) | ||  __/ |_) |  | |
607 +  //  |_____/ \__\___| .__/   |_|
608 +  //                 | |        
609 +  //                 |_|        
610 +  //  
611 +  // Fit OF to get good initial values for flavor-symmetric background
612    
499  //RooDataSet combData("combData","combined data",mll,Index(sample),Import("SF",*SFSample),Import("mm",*mmSample),Import("OF",*OFSample));
500  RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("SF",*SFSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
613    
614    
615    //First we make a fit to opposite flavor
616 <  RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 10000);
617 <  RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
616 >  RooRealVar fttbarOF("fttbarOF", "fttbarOF", OFSample->sumEntries(), 0, combData.sumEntries());
617 >  RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.77, 0.01, 4.0);
618    RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
619 <  RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
619 >  RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.0272, 0.001, 1.0);
620    RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
621    RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
622    RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
# Line 530 | Line 642 | void EdgeFitter::DoFit(int is_data, floa
642    frameO->Draw();
643    if(is_data==data) DrawPrelim();
644    else DrawPrelim(PlottingSetup::luminosity,true);
645 <  ASP->Draw();
534 <  if(EdgeFitter::FixedMEdge>=0) CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str()+"__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
535 <  else CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str(),false,false);
645 >  CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str(),false,false);
646    delete pof;
647 +  
648 +  TCanvas* pof2 = new TCanvas("pof2","pof2") ;
649 +  pof2->cd();
650 +  gPad->SetLeftMargin(0.15);
651 +  frameO->GetYaxis()->SetTitleOffset(1.4);
652 +  frameO->Draw();
653 +  if(is_data==data) DrawPrelim();
654 +  else DrawPrelim(PlottingSetup::luminosity,true);
655 +  stringstream OFFitParameterSummary;
656 +  OFFitParameterSummary << "#splitline{ftbbarOF = " <<  WriteWithError(fttbarOF.getVal(),fttbarOF.getError(),3)  << "}{";
657 +  OFFitParameterSummary << "#splitline{par1ttbarOF = " <<  WriteWithError(par1ttbarOF.getVal(),par1ttbarOF.getError(),3)  << "}{";
658 +  OFFitParameterSummary << "#splitline{par2ttbarOF = " <<  WriteWithError(par2ttbarOF.getVal(),par2ttbarOF.getError(),3)  << "}{";
659 +  OFFitParameterSummary << "#splitline{par3ttbarOF = " <<  WriteWithError(par3ttbarOF.getVal(),par3ttbarOF.getError(),3)  << "}{";
660 +  OFFitParameterSummary << "par4ttbarOF = " <<  WriteWithError(par4ttbarOF.getVal(),par4ttbarOF.getError(),3);
661 +  OFFitParameterSummary << "}}}}";
662 +  TLatex *OFFitinfobox = new TLatex(0.57,0.75,OFFitParameterSummary.str().c_str());
663 +  OFFitinfobox->SetNDC();
664 +  OFFitinfobox->SetTextSize(0.03);
665 +  OFFitinfobox->Draw();
666 +
667 +  
668 +  CompleteSave(pof2,"Edge/OF__OFFitonly_"+prefix.str()+"__INFO",false,false);
669 +  
670 +  delete pof2;
671      
672    if(resultOF->covQual()!=3) {
673      write_error(__FUNCTION__,"OF fit did not converge!!! Cannot continue!");
674      dout << "covQual is " << resultOF->covQual() << endl;
675 <    EdgeFitter::FixedMEdgeChi2=-1;
675 >    if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
676 >    else EdgeFitter::FixedMEdgeChi2_H0=-1;
677      if(EdgeFitter::RejectPointIfNoConvergence) return;
678    } else {
679      write_info(__FUNCTION__,"OF fit converged");
680    }
681    
682 +  
683 +  //    _____ _               ___  
684 +  //   / ____| |             |__ \
685 +  //  | (___ | |_ ___ _ __      ) |
686 +  //   \___ \| __/ _ \ '_ \    / /
687 +  //   ____) | ||  __/ |_) |  / /_
688 +  //  |_____/ \__\___| .__/  |____|
689 +  //                 | |          
690 +  //                 |_|          
691 +  //    
692 +  // Set up all the models
693 +  
694 +  
695 +  // Step 2a: set up edge position (if fixed), and maximum Z yield
696    float StartingMedge=70;
697    if(EdgeFitter::FixedMEdge>0) StartingMedge=EdgeFitter::FixedMEdge;
549
698  
699    RooDataSet *ZDataSet = (RooDataSet*)EdgeFitter::AllData->reduce("id1==id2 && abs(mll-91.2)<20");
700  
# Line 554 | Line 702 | void EdgeFitter::DoFit(int is_data, floa
702    dout << "maxZ was set to " << maxZ << endl;
703    delete ZDataSet;
704    
705 <  
558 <  // Now same flavor  
559 <  RooRealVar fzSF("fzSF", "fzSF", 39, 39, maxZ);
560 <  RooRealVar meanzSF("meanzSF", "meanzSF", 91.1876, 89, 95);
561 <  //RooRealVar sigmazSF("sigmazSF", "sigmazSF", 0.5);
562 <  RooRealVar sigmazSF("sigmazSF", "sigmazSF", 5, 0.5, 5);
563 <  RooRealVar widthzSF("widthzSF", "widthzSF", 2.94);
564 <  widthzSF.setConstant();
565 <  
566 <  RooRealVar fttbarSF("fttbarSF", "fttbarSF", 100, 0, 100000);
567 <  RooRealVar par1ttbarSF("par1ttbarSF", "par1ttbarSF", 1.02*par1ttbarOF.getVal(), (1.02-0.07)*par1ttbarOF.getVal(), (1.02+0.07)*par1ttbarOF.getVal());
568 <
569 <  RooRealVar fsignalSF("fsignalSF", "fsignalSF", 0, 0, 300);
570 <  RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
571 <  RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
572 <  RooRealVar par3signalSF("par3signalSF", "par3signalSF", StartingMedge, 0, 300);
573 <
574 <  RooVoigtian zSF("zSF", "zSF", mll, meanzSF, widthzSF, sigmazSF);
705 >  RooRealVar JustOne("JustOne","Just One",1) ;
706  
576  if(EdgeFitter::FixedMEdge>0) par3signalSF.setConstant();
707    
708 <  RooSUSYBkgPdf ttbarSF("ttbarSF","ttbarSF", mll , par1ttbarSF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
709 <  RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, sigmazSF, par3signalSF);
708 >  // Step 2b: set up the models!
709 >  write_info(__FUNCTION__,"ee & mm fractions fixed");
710 > /*  RooRealVar eefrac("eefrac","eefrac",eeSample->sumEntries()/(eeSample->sumEntries()+mmSample->sumEntries()),0.1,1.0);
711 >  RooFormulaVar mmfrac("mmfrac","1-eefrac",RooArgSet(eefrac));*/
712 >  RooRealVar eefrac("eefrac","eefrac",eeSample->sumEntries()/(eeSample->sumEntries()+mmSample->sumEntries()));
713 >  RooFormulaVar mmfrac("mmfrac","1-eefrac",RooArgSet(eefrac));
714 >  
715 >  RooRealVar rOneOverEMFrac("rOneOverEMFrac","rOneOverEMFrac",1.02,1.02-0.07,1.02+0.07);
716 >  RooRealVar rmean("rmean","rmean",1.02);
717 >  RooRealVar rwidth("rwidth","rwidth",0.07);
718 >  RooGaussian OneOverEMFrac("OneOverEMFrac","OneOverEMFrac",rOneOverEMFrac,rmean,rwidth);
719 >  
720 >  
721 >  // FIRST: Drell-Yan: different parameters for resolution of DY (and of course yield), but same mean and so on
722 >  RooRealVar widthz("widthz", "widthz", 2.94);
723 >  widthz.setConstant(1);
724 >  RooRealVar meanz("meanz", "meanz", 91.1876, 89, 93);
725 >  RooRealVar fzSF("fzSF","fzSF",39,0,maxZ);
726 >  RooFormulaVar fzee("fzee","eefrac*fzSF",RooArgSet(eefrac,fzSF));
727 >  RooFormulaVar fzmm("fzmm","mmfrac*fzSF",RooArgSet(mmfrac,fzSF));
728 >  RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0.5, 5);
729 >  RooRealVar sigmazmm("sigmazmm", "sigmazmm", 5, 0.5, 5);
730 >  
731 >  RooVoigtian zee("zee", "zee", mll, meanz, widthz, sigmazee);
732 >  RooVoigtian zmm("zmm", "zmm", mll, meanz, widthz, sigmazmm);
733 >  RooAddPdf zSF("zSF","zSF",RooArgList(zee,zmm),RooArgList(fzee,fzmm));
734 >  
735 >  
736 >  //SECOND: Flavor-Symmetry contribution. Only relative fraction is different, all other parameters are the same;
737 >  RooFormulaVar fttbaree("fttbaree","eefrac*fttbarOF*OneOverEMFrac",RooArgSet(eefrac,fttbarOF,OneOverEMFrac));
738 >  RooFormulaVar fttbarmm("fttbarmm","mmfrac*fttbarOF*OneOverEMFrac",RooArgSet(mmfrac,fttbarOF,OneOverEMFrac));
739 >  
740 >  RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
741 >  RooSUSYBkgPdf ttbarmm("ttbarmm","ttbarmm", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
742 >  RooAddPdf ttbarSF("ttbarSF","ttbarSF", RooArgList(ttbaree,ttbarmm), RooArgList(fttbaree,fttbarmm));
743 >  RooAddPdf ttbarSFOne("ttbarSF","ttbarSF", RooArgList(ttbaree,ttbarmm), RooArgList(JustOne,JustOne));
744 >  
745 >  
746 >  //for signal only the resolution and absolute fraction are different for ee / mm
747 >  RooRealVar fsignal("fsignal", "fsignal", 0, 0, 300);
748 >  RooFormulaVar fsignalee("fsignalee","eefrac*fsignal",RooArgSet(eefrac,fsignal));
749 >  RooFormulaVar fsignalmm("fsignalmm","mmfrac*fsignal",RooArgSet(mmfrac,fsignal));
750 >  
751 >  RooRealVar par1signal("par1signal", "par1signal", 45, 20, 100);
752 >  RooRealVar par3signal("par3signal", "par3signal", StartingMedge, 0, 300);
753 >  if(EdgeFitter::FixedMEdge>0) par3signal.setConstant();
754 >  
755 >  RooSUSYTPdf signalee("signalee","signalee", mll , par1signal, sigmazee, par3signal);
756 >  RooSUSYTPdf signalmm("signalmm","signalmm", mll , par1signal, sigmazmm, par3signal);
757 >  RooAddPdf signalSF("signalSF","signalSF",RooArgList(signalee,signalmm),RooArgList(fsignalee,fsignalmm));
758 >  
759 >  if(!EdgeFitter::AllowTriangle) {
760 >    fsignal.setVal(0.0); // kill off the signal if we don't want the triangle
761 >    fsignal.setConstant(1);
762 >    par1signal.setConstant(1);
763 >    par3signal.setConstant(1);
764 >  }
765 >  
766    
767 <  RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
767 >  write_error(__FUNCTION__,"RnD version: fixing most parameters so the fit works very quickly (only one parameters instead of eight) !");
768 > //          fttbarOF.setConstant(1);
769 >            fzSF.setConstant(1);
770 >            meanz.setConstant(1);
771 >            par1ttbarOF.setConstant(1);
772 >            par3ttbarOF.setConstant(1);
773 >            sigmazee.setConstant(1);
774 >            sigmazmm.setConstant(1);
775 >  
776 >  
777 >  //COMPLETE MODEL:
778 >
779 >  RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
780 >  RooAddPdf model_mm("model_mm","model_mm", RooArgList(zmm, ttbarmm, signalmm), RooArgList(fzmm, fttbarmm, fsignalmm));
781 > //  RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
782 >  RooAddPdf model_SF("model_SF","model_SF",RooArgList(model_ee,model_mm),RooArgList(JustOne,JustOne));
783 >  
784 >  
785 >  //    _____ _               ____  
786 >  //   / ____| |             |___ \
787 >  //  | (___ | |_ ___ _ __     __) |
788 >  //   \___ \| __/ _ \ '_ \   |__ <
789 >  //   ____) | ||  __/ |_) |  ___) |
790 >  //  |_____/ \__\___| .__/  |____/
791 >  //                 | |            
792 >  //                 |_|            
793 >  //  
794 >  //  Fitting
795    
796    RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
797 <  simPdf.addPdf(model_SF,"SF") ;
797 >  simPdf.addPdf(model_ee,"ee") ;
798 >  simPdf.addPdf(model_mm,"mm") ;
799    simPdf.addPdf(model_OF,"OF") ;
800  
801    RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended(),RooFit::Minos(true));
# Line 589 | Line 803 | void EdgeFitter::DoFit(int is_data, floa
803    if(result->covQual()!=3) {
804      write_error(__FUNCTION__,"Full fit did not converge!!! Cannot continue!");
805      dout << "covQual is " << result->covQual() << endl;
806 <    EdgeFitter::FixedMEdgeChi2=-1;
806 >    if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
807 >    else EdgeFitter::FixedMEdgeChi2_H0=-1;
808      if(EdgeFitter::RejectPointIfNoConvergence) return;
809    } else {
810      write_info(__FUNCTION__,"Full fit converged");
811    }
812  
813 <  dout << "============================= < OF results > =============================" << endl;
813 >  dout << "============================= < Full results > =============================" << endl;
814    result->Print();
815 <  dout << "============================= < /OF results > =============================" << endl;
815 >  dout << "============================= < /Full results > =============================" << endl;
816  
817 +  //    _____ _               _  _  
818 +  //   / ____| |             | || |  
819 +  //  | (___ | |_ ___ _ __   | || |_
820 +  //   \___ \| __/ _ \ '_ \  |__   _|
821 +  //   ____) | ||  __/ |_) |    | |  
822 +  //  |_____/ \__\___| .__/     |_|  
823 +  //                 | |            
824 +  //                 |_|            
825 +  //
826 +  //   Present results
827    
828    RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
829    frame1->GetXaxis()->CenterTitle(1);
830    frame1->GetYaxis()->CenterTitle(1);
831 <  combData.plotOn(frame1,RooFit::Name("SFdata"),RooFit::Cut("sample==sample::SF")) ;
832 <  simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
833 <  simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("TTbarSFonly"),RooFit::Components("ttbarSF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
834 <  simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("DYSFonly"),RooFit::Components("zSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
835 <  simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("SignalSFonly"),RooFit::Components("signalSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
836 <  
837 <  EdgeFitter::FixedMEdgeChi2 = frame1->chiSquare("FullFit", "SFdata", 3);
838 <  
839 <  
840 <  dout << "Result   : " << endl;
841 <  dout << "f signal : " << fsignalSF.getVal() << " +/- " << fsignalSF.getError() << endl;
842 <  dout << "f ttbar  : " << fttbarSF.getVal() << " +/- " << fttbarSF.getError() << endl;
843 <  dout << "f tt OF  : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
844 <  dout << "f z SF   : " << fzSF.getVal() << " +/- " << fzSF.getError() << endl;
845 <  dout << "#Chi^{2}/NDF     : " << EdgeFitter::FixedMEdgeChi2 << endl;
831 >  combData.plotOn(frame1,RooFit::Name("eedata"),RooFit::Cut("sample==sample::ee")) ;
832 >  simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("FullFitEE"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
833 >  simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("TTbarEEonly"),RooFit::Components("ttbaree"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
834 >  simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("DYEEonly"),RooFit::Components("zee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
835 >  if(EdgeFitter::AllowTriangle) simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("SignalEEonly"),RooFit::Components("signalee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));  
836 >  
837 >  
838 >  RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("mm sample")) ;
839 >  frame2->GetXaxis()->CenterTitle(1);
840 >  frame2->GetYaxis()->CenterTitle(1);
841 >  combData.plotOn(frame2,RooFit::Name("mmdata"),RooFit::Cut("sample==sample::mm")) ;
842 >  simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("FullFitMM"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
843 >  simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("TTbarMMonly"),RooFit::Components("ttbarmm"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
844 >  simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("DYMMonly"),RooFit::Components("zmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
845 >  if(EdgeFitter::AllowTriangle) simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("SignalMMonly"),RooFit::Components("signalmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));  
846 >  
847    
622  // The same plot for the cointrol sample slice
848    RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
849    frame3->GetXaxis()->CenterTitle(1);
850    frame3->GetYaxis()->CenterTitle(1);
851 <  frame3->SetMaximum(frame1->GetMaximum());
852 <  combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
853 <  simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
854 <  simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
855 <  
851 >  combData.plotOn(frame3,RooFit::Name("OFdata"),RooFit::Cut("sample==sample::OF")) ;
852 >  simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("FullFitOF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
853 >  simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("TTbarOFonly"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
854 >  simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("DYOFonly"),RooFit::Components("zOF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
855 >
856 >  RooAbsPdf *ee_result = simPdf.getPdf("ee");
857 >  RooAbsPdf *mm_result = simPdf.getPdf("mm");
858 >  RooAddPdf SF_result("SF_result","SF_result",RooArgList(*ee_result,*mm_result));
859 >  
860 >  
861 >  TCanvas* g = new TCanvas("g","g") ;
862 >  
863 >  RooPlot* frame4 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("SF sample")) ;
864 >  frame4->GetXaxis()->CenterTitle(1);
865 >  frame4->GetYaxis()->CenterTitle(1);
866 >  EdgeFitter::SFSample->plotOn(frame4,RooFit::Name("SFdata")) ;
867 >  SF_result.plotOn(frame4,RooFit::Name("FullFitSF"),RooFit::ProjWData(*EdgeFitter::SFSample),RooFit::LineColor(kBlack));
868 >  SF_result.plotOn(frame4,RooFit::Name("TTbarSFonly"),RooFit::ProjWData(*EdgeFitter::SFSample),RooFit::Components("ttbaree,ttbarmm"),RooFit::LineColor(kBlue),RooFit::LineStyle(kDashed));
869 >  SF_result.plotOn(frame4,RooFit::Name("DYSFonly"),RooFit::ProjWData(*EdgeFitter::SFSample),RooFit::Components("zmm,zee"), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
870 >  if(EdgeFitter::AllowTriangle) SF_result.plotOn(frame4,RooFit::Name("SignalSFonly"),RooFit::Components("signalSF"), RooFit::ProjWData(*EdgeFitter::SFSample), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));  
871    
872 + /// ********************************************************************************************************************************
873    TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
874    c->cd() ;
875    gPad->SetLeftMargin(0.15);
# Line 636 | Line 877 | void EdgeFitter::DoFit(int is_data, floa
877    frame1->Draw();
878    if(is_data==data) DrawPrelim();
879    else DrawPrelim(PlottingSetup::luminosity,true);
880 <  stringstream infotext;
640 <  infotext << "#splitline{Fit results (" << EdgeFitter::Mode << ">" << jzb_cut << "): }{#splitline{";
641 <  infotext << "N(Data) = " << EdgeFitter::SFSample->sumEntries() << "}{#splitline{";
642 <  infotext << "N(Z+Jets) = " << WriteWithError(fzSF.getVal(),fzSF.getError(),3) << "}{#splitline{";
643 <  infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
644 <  infotext << "N(signal) = " << WriteWithError(fsignalSF.getVal(),fsignalSF.getError(),3) << "}{";
645 <  infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";  
646 <
647 <  TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
648 <  infobox->SetNDC();
649 <  infobox->SetTextSize(0.03);
650 <  infobox->Draw();
651 <  ASP->Draw();
652 <  if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_SF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
653 <  else CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
654 <  delete c;
880 >  CompleteSave(c,"Edge/"+prefix.str()+"_ee",false,false);
881    
882 <  TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
883 <  e->cd();
882 >  c->cd() ;
883 >  gPad->SetLeftMargin(0.15);
884 >  frame2->GetYaxis()->SetTitleOffset(1.4);
885 >  frame2->Draw();
886 >  if(is_data==data) DrawPrelim();
887 >  else DrawPrelim(PlottingSetup::luminosity,true);
888 >  CompleteSave(c,"Edge/"+prefix.str()+"_mm",false,false);
889 >  
890 >  c->cd() ;
891    gPad->SetLeftMargin(0.15);
892    frame3->GetYaxis()->SetTitleOffset(1.4);
893    frame3->Draw();
894    if(is_data==data) DrawPrelim();
895    else DrawPrelim(PlottingSetup::luminosity,true);
896 <  ASP->Draw();
897 <  if(EdgeFitter::FixedMEdge>=0) CompleteSave(e,"Edge/"+prefix.str()+"_OF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
898 <  else CompleteSave(e,"Edge/"+prefix.str()+"_OF",false,false);
666 <  delete e;
667 <  
668 <  
669 <  
670 <  
671 < /*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
672 <  f->cd();
896 >  CompleteSave(c,"Edge/"+prefix.str()+"_OF",false,false);
897 >
898 >  c->cd() ;
899    gPad->SetLeftMargin(0.15);
900    frame4->GetYaxis()->SetTitleOffset(1.4);
901    frame4->Draw();
902    if(is_data==data) DrawPrelim();
903    else DrawPrelim(PlottingSetup::luminosity,true);
904 <  CompleteSave(f,"Edge/"+prefix.str()+"_SF");
679 <  delete f;*/
680 <
904 >  CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
905  
906 < /*
907 <  float maxZ=200;
684 <  RooWorkspace* wspace = new RooWorkspace();
685 <  stringstream mllvar;
686 <  mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
687 <  wspace->factory(mllvar.str().c_str());
688 <  wspace->var("mll")->setBins(30);
689 <  wspace->factory("nSig[1.,0.,100.]");
690 <  wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
691 <  wspace->factory("rME[1.12,1.05,1.19]");
692 <  wspace->factory("effUncert[1.]");
693 <  EdgeFitter::prepareLimits(wspace, true);
694 < */
906 >  delete c;
907 > /// ********************************************************************************************************************************
908  
909 < write_warning(__FUNCTION__," A lot missing here to calculate limits");
909 >
910 > //  if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=frame2->chiSquare("FullFit", "EEdata", 3);
911 > //  else EdgeFitter::FixedMEdgeChi2_H0=frame2->chiSquare("FullFit", "EEdata", 3);
912 >  
913 >  
914  
915   }
916  
# Line 713 | Line 930 | void EdgeFitter::DoEdgeFit(string mcjzb,
930    RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
931    
932    
933 +  EdgeFitter::AllowTriangle=false;
934 +  EdgeFitter::DoFit(is_data, jzb_cut);
935 +  
936 +  EdgeFitter::AllowTriangle=true;
937 +  
938    bool ScanMassRange=false;
939    float ScanSteps=5.0;//GeV
940    
941    
720  
942    if(ScanMassRange) {
943      TFile *fscan = new TFile("fscan.root","UPDATE");
944      TGraph *gr = new TGraph();
945 +    TGraph *Rgr = new TGraph();
946      stringstream GrName;
947      GrName << "ScanGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
948 +    stringstream RGrName;
949 +    RGrName << "ScanRatioGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
950      gr->SetName(GrName.str().c_str());
951 +    Rgr->SetName(RGrName.str().c_str());
952  
953      int i=0;
954      for(float tempMedge=10;tempMedge<=300;tempMedge+=ScanSteps) {
955        write_info(__FUNCTION__,"Now testing Medge="+any2string(tempMedge)+" for "+EdgeFitter::Mode+">"+any2string(jzb_cut));
956        EdgeFitter::FixedMEdge=tempMedge;
957        EdgeFitter::DoFit(is_data, jzb_cut);
958 <      if(EdgeFitter::FixedMEdgeChi2<0) continue;
959 <      gr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2);
958 >      if(EdgeFitter::FixedMEdgeChi2_H1<0) continue;
959 >      gr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1);
960 >      Rgr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1/EdgeFitter::FixedMEdgeChi2_H0);
961        i++;
962      }
963      
# Line 750 | Line 976 | void EdgeFitter::DoEdgeFit(string mcjzb,
976      if(is_data) DrawPrelim();
977      else DrawMCPrelim();
978      CompleteSave(ScanCan,ScanCanSave.str());
979 +
980 +    Rgr->GetXaxis()->SetTitle("m_{edge}");
981 +    Rgr->GetXaxis()->CenterTitle();
982 +    Rgr->GetYaxis()->SetTitle("#Chi^{2} / NDF");
983 +    Rgr->GetYaxis()->CenterTitle();
984 +    Rgr->GetYaxis()->SetTitleOffset(0.95);
985 +    Rgr->GetXaxis()->SetTitleOffset(0.9);
986 +    Rgr->SetLineColor(kBlue);
987 +    Rgr->SetTitle("");
988 +    Rgr->Draw("AL");
989 +    ScanCanSave.str("");
990 +    ScanCanSave << "Edge/MEdgeScan_Ratio_"+EdgeFitter::Mode+"_" << jzb_cut;
991 +    if(is_data) DrawPrelim();
992 +    else DrawMCPrelim();
993 +    CompleteSave(ScanCan,ScanCanSave.str());
994      fscan->cd();
995      gr->Write();
996      delete ScanCan;
997      fscan->Close();
998    } else {
999      EdgeFitter::DoFit(is_data, jzb_cut);
1000 +    dout << "Chi^2 (H0) = " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
1001 +    dout << "Chi^2 (H1) = " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
1002    }
1003    
761      
762  
1004    
764  EdgeFitter::DoFit(is_data, jzb_cut);
1005    RooMsgService::instance().setGlobalKillBelow(msglevel);
1006  
767
1007    f->Close();
1008  
1009   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines