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.5 by buchmann, Wed Jun 27 08:53:23 2012 UTC vs.
Revision 1.15 by buchmann, Thu Jun 13 11:41:33 2013 UTC

# Line 1 | Line 1
1   #include <iostream>
2  
3 + #include <TVirtualIndex.h>
4 +
5   #include <RooRealVar.h>
6   #include <RooArgSet.h>
7   #include <RooDataSet.h>
# Line 23 | Line 25
25   #include "RooStats/HypoTestInverterOriginal.h"
26  
27   //#include "ParametrizedEdge.C"
28 < #include "/shome/pablom/RooFit/Pdfs/RooSUSYTPdf.cxx"
29 < #include "/shome/pablom/RooFit/Pdfs/RooSUSYBkgPdf.cxx"
28 > #include "EdgeModules/RooSUSYTPdf.cxx"
29 > #include "EdgeModules/RooSUSYBkgPdf.cxx"
30  
31  
32   using namespace std;
# Line 51 | Line 53 | namespace EdgeFitter {
53    
54    void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree*);
55    void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree*);
56 <  void getMedianLimit(RooWorkspace *ws,vector<RooDataSet*> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
56 >  void getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
57    void InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut);
58    void PrepareDatasets(int);
59 +  void DrawDatasetContent(int);
60    void DoFit(int is_data, float jzb_cut);
61    string RandomStorageFile();
62    Yield Get_Z_estimate(float,int);
63    Yield Get_T_estimate(float,int);
64 <  float calcExclusion(RooWorkspace *ws, RooDataSet *data = NULL);
65 <  void prepareLimits(RooWorkspace *ws);
63 <  vector<RooDataSet*> generateToys(RooWorkspace *ws, int nToys);
64 >  float calcExclusion(RooWorkspace *ws, RooDataSet data, bool calcExclusion);
65 >  vector<RooDataSet> generateToys(RooWorkspace *ws, int nToys);
66    void prepareLimits(RooWorkspace *ws, bool ComputeBands);
67    TGraph* prepareLM(float mass, float nEv);
68    
# Line 70 | Line 72 | namespace EdgeFitter {
72    TCut cut;
73    
74    RooDataSet* AllData;
75 <  RooDataSet* eeSample;
76 <  RooDataSet* mmSample;
77 <  RooDataSet* emSample;
75 >  RooDataSet* SFSample;
76 >  RooDataSet* OFSample;
77 >  
78 >  bool MarcoDebug=true;
79 >  
80 >  float FixedMEdge=-1;
81 >  float FixedMEdgeChi2=-1;
82 >  
83 >  bool RejectPointIfNoConvergence=false;
84 >  
85 >  string Mode="UndefinedMode";
86    
77  bool MarcoDebug;
87   }
88  
89   TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
# Line 93 | Line 102 | TGraph* EdgeFitter::prepareLM(float mass
102    return lm;
103   }
104  
105 < vector<RooDataSet*> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
105 > vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
106 >  ws->ls();
107    ws->var("nSig")->setVal(0.);
108    ws->var("nSig")->setConstant(true);
109    RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
110 <  vector<RooDataSet*> theToys;
110 >  vector<RooDataSet> theToys;
111    
112    RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
113    mcEE.generate(nToys,44,true);
# Line 106 | Line 116 | vector<RooDataSet*> EdgeFitter::generate
116    RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
117    mcOSOF.generate(nToys,44,true);
118    
119 <  RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
119 >  RooRealVar mll("m_{ll}","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
120    RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
121    RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
122    RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
# Line 120 | Line 130 | vector<RooDataSet*> EdgeFitter::generate
130      stringstream toyname;
131      toyname << "theToy_" << i;
132      write_warning(__FUNCTION__,"Problem while adding toys");
133 < //    RooDataSet *toyData = RooDataSet(toyname.str(),toyname.str(),observables,RooFit::Index(ws->cat("cat")),RooFit::Import("OSOF",*toyOSOF),RooFit::Import("EE",*toyEE),RooFit::Import("MM",*toyMM));
134 < //    theToys.push_back(toyData);
133 >    RooDataSet toyData = RooDataSet(toyname.str().c_str(),toyname.str().c_str(),observables,RooFit::Index(const_cast<RooCategory&>(*ws->cat("cat"))),RooFit::Import("OSOF",*toyOSOF),RooFit::Import("EE",*toyEE),RooFit::Import("MM",*toyMM));
134 >    theToys.push_back(toyData);
135    }
136    ws->var("nSig")->setVal(17.0);
137    ws->var("nSig")->setConstant(false);
138    return theToys;
139   }
140  
141 < void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet*> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
141 > void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
142    TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
143    vector<float> theLimits;
144 <  for(int itoy=0;itoy<theToys.size();itoy++) {
145 <    float theLimit = calcExclusion(ws,theToys[itoy]);
144 >  for(int itoy=0;itoy<(int)theToys.size();itoy++) {
145 >    float theLimit = calcExclusion(ws,theToys[itoy],false);
146      if(theLimit > 0 ) gauLimit->Fill(theLimit);
147    }
148    const Int_t nQ = 4;
# Line 150 | Line 160 | void EdgeFitter::getMedianLimit(RooWorks
160  
161   void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
162    if(ComputeBands) {
163 <    vector<RooDataSet*> theToys = EdgeFitter::generateToys(ws,50);
163 >    vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
164      vector<float> medVals;
165      vector<float> medLimits;
166      vector<float> sigmaLimitsDown;
# Line 179 | Line 189 | void EdgeFitter::prepareLimits(RooWorksp
189        theVals.push_back((float)i);
190        ws->var("m0")->setVal((float)i);
191        ws->var("m0")->setConstant(true);
192 <      theLimits.push_back(calcExclusion(ws));
192 > //      theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
193 >      write_error(__FUNCTION__,"Error while casting roo data set");
194      }
195      
196 <    for(int i=0;i<theLimits.size();i++) {
196 >    for(int i=0;i<(int)theLimits.size();i++) {
197        if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
198          cout << i << " : " << theVals[i] << endl;
199          theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
# Line 194 | Line 205 | void EdgeFitter::prepareLimits(RooWorksp
205   }
206    
207  
208 < float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet *data) {
208 > float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
209 >  int numberOfToys=50;
210    RooRealVar mu("mu","nSig",0,10000,"");
211    RooArgSet poi = RooArgSet(mu);
212    RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
# Line 203 | Line 215 | float EdgeFitter::calcExclusion(RooWorks
215    model->SetWorkspace(*ws);
216    model->SetPdf("combModel");
217    model->SetParametersOfInterest(poi);
218 <  if(!data) data = (RooDataSet*)ws->data("data_obs");
218 > //  if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
219  
220 <  RooStats::ProfileLikelihoodCalculator plc(*data, *model);
220 >  RooStats::ProfileLikelihoodCalculator plc(data, *model);
221    plc.SetNullParameters(*nullParams);
222    plc.SetTestSize(0.05);
223 +  
224    RooStats::LikelihoodInterval* interval = plc.GetInterval();
225    RooStats::HypoTestResult *htr = plc.GetHypoTest();
226    double theLimit = interval->UpperLimit( mu );
227 <  cout << "Significance " << htr->Significance() << endl;
227 > //  double significance = htr->Significance();
228    
229    ws->defineSet("obs","nB");
230    ws->defineSet("poi","nSig");
# Line 235 | Line 248 | float EdgeFitter::calcExclusion(RooWorks
248    slrts.SetAltParameters(*sb_model.GetSnapshot());
249    RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
250    
251 <  RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(*data, sb_model, b_model,0,0);
251 >  RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
252    hc.SetTestStatistic(2);
253 <  hc.SetNumberOfToys(50);
253 >  hc.SetNumberOfToys(numberOfToys);
254    
255    RooStats::HypoTestInverterOriginal hcInv =  RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
256    hcInv.SetTestSize(0.05);
# Line 253 | Line 266 | float EdgeFitter::calcExclusion(RooWorks
266   }
267  
268   TTree* SkimTree(int isample) {
269 <  TTree* newTree = allsamples.collection[isample].events->CloneTree(0);
269 >  string TreeName = allsamples.collection[isample].filename;
270 >  TTree* newTree = new TTree("NanoTree",TreeName.c_str());
271 >  
272 >  float mll,edgeWeight;
273 >  int id1,id2;
274 >  
275 >  newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
276 >  newTree->Branch("mll",&mll,"mll/F");
277 >  newTree->Branch("id1",&id1,"id1/I");
278 >  newTree->Branch("id2",&id2,"id2/I");
279 >  
280    float xsweight=1.0;
281    if(allsamples.collection[isample].is_data==false) xsweight=luminosity*allsamples.collection[isample].weight;
282    if(EdgeFitter::MarcoDebug) {
283      cout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
284      cout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
285    }
286 +  
287 +  float tmll;
288 +  int tid1,tid2;
289 +  allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
290 +  allsamples.collection[isample].events->SetBranchAddress("id1",&tid1);
291 +  allsamples.collection[isample].events->SetBranchAddress("id2",&tid2);
292 +  
293    TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
294 +  TTreeFormula *Weight = new TTreeFormula("Weight", cutWeight, allsamples.collection[isample].events);
295 +  
296    float wgt=1.0;
265  allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
297    for (Int_t entry = 0 ; entry < allsamples.collection[isample].events->GetEntries() ; entry++) {
298     allsamples.collection[isample].events->LoadTree(entry);
299     if (select->EvalInstance()) {
300       allsamples.collection[isample].events->GetEntry(entry);
301 <     wgt=wgt*xsweight;
301 >     mll=tmll;
302 >     id1=tid1;
303 >     id2=tid2;
304 >     wgt=Weight->EvalInstance();
305 >     edgeWeight=wgt*xsweight;
306       newTree->Fill();
307     }
308    }
# Line 282 | Line 317 | void EdgeFitter::InitializeVariables(flo
317    jzbmax=_jzbmax;
318    cut=_cut;
319   }
320 <  
320 >
321 > TTree* MergeTrees(vector<TTree*> trees) {
322 >  cout << "Quick overview of tree addresses : " << endl;
323 >  for(int i=0;i<(int)trees.size();i++) {
324 >    cout << "Tree " << i << " : " << trees[i] << endl;
325 >  }
326 >  cout << "Address of first tree : " << trees[0] << endl;
327 >  TTree * newtree = (TTree*)trees[0]->CloneTree(0);
328 >  trees[0]->GetListOfClones()->Remove(newtree);
329 >  trees[0]->ResetBranchAddresses();
330 >  newtree->ResetBranchAddresses();
331 >  
332 >  for(int itree=0;itree<trees.size();itree++) {
333 >    newtree->CopyAddresses(trees[itree]);
334 >    Long64_t nentries = trees[itree]->GetEntries();
335 >    for (Long64_t iEntry=0;iEntry<nentries;iEntry++) {
336 >      trees[itree]->GetEntry(iEntry);
337 >      newtree->Fill();
338 >    }
339 >    trees[itree]->ResetBranchAddresses(); // Disconnect from new tree
340 >    if (newtree->GetTreeIndex()) {
341 >      newtree->GetTreeIndex()->Append(trees[itree]->GetTreeIndex(),kTRUE);
342 >    }
343 >    if (newtree && newtree->GetTreeIndex()) {
344 >      newtree->GetTreeIndex()->Append(0,kFALSE); // Force the sorting
345 >    }
346 >  }
347 >  return newtree;
348 > }
349 >    
350 >      
351 >
352   void EdgeFitter::PrepareDatasets(int is_data) {
287  TTree *completetree;
353    write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
354 <  bool hashit=0;
355 <  for(int isample=0;isample<allsamples.collection.size();isample++) {
354 >  TFile *tempout = new TFile("tempout.root","RECREATE");
355 >  vector<TTree*> SkimmedTrees;
356 >  TTree *SkimmedTree[(int)allsamples.collection.size()];
357 >  for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
358      if(!allsamples.collection[isample].is_active) continue;
359      if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
360      if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
361      if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
362      if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
363      if(EdgeFitter::MarcoDebug) cout << "Considering : " << allsamples.collection[isample].samplename << endl;
364 <    if(!hashit) {
365 <      hashit=true;
366 <      completetree = SkimTree(isample)->CloneTree();
367 <    } else {
301 <      completetree->CopyEntries(SkimTree(isample));
302 <    }
303 <    if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
364 >    SkimmedTree[isample] = SkimTree(isample);
365 >    tempout->cd();
366 >    SkimmedTree[isample]->Write();
367 >    SkimmedTrees.push_back(SkimmedTree[isample]);
368    }
369    
370 <  RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
370 >  TTree *completetree  = MergeTrees(SkimmedTrees);
371 >  
372 >  if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
373 >  
374 >  RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
375    RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
376    RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
377 <  RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
378 <  RooRealVar weight("weight","weight",0,1000,"");
311 <  RooArgSet observables(mll,jzb,id1,id2,weight);
377 >  RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
378 >  RooArgSet observables(mll,id1,id2,edgeWeight);
379    
380    string title="CMS Data";
381    if(is_data!=1) title="CMS SIMULATION";
382 <  RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","weight");
382 >  RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
383 >  tempout->cd();
384    completetree->Write();
385 < //  delete completetree;
385 >  tempout->Close();
386 >
387    
388 <  EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2");
389 <  EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2");
321 <  EdgeFitter::emSample = (RooDataSet*)LAllData.reduce("id1!=id2");
388 >  EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
389 >  EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
390    EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
391    
392 <  eeSample->SetName("eeSample");
393 <  mmSample->SetName("mmSample");
326 <  emSample->SetName("emSample");
392 >  SFSample->SetName("SFSample");
393 >  OFSample->SetName("OFSample");
394    AllData->SetName("AllData");
395    
396    if(EdgeFitter::MarcoDebug) {
397 <    cout << "Number of events in data sample = " << AllData->numEntries() << endl;
398 <    cout << "Number of events in ee sample = " << eeSample->numEntries() << endl;
399 <    cout << "Number of events in mm sample = " << mmSample->numEntries() << endl;
333 <    cout << "Number of events in em sample = " << emSample->numEntries() << endl;
397 >    cout << "Number of events in data sample = " << AllData->sumEntries() << endl;
398 >    cout << "Number of events in eemm sample = " << SFSample->sumEntries() << endl;
399 >    cout << "Number of events in em sample = " << OFSample->sumEntries() << endl;
400    }
401 +  
402 + }
403 +
404 + void EdgeFitter::DrawDatasetContent(int is_data) {
405 +  RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
406 +  RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("SF sample")) ;
407 +  frame1->GetXaxis()->CenterTitle(1);
408 +  frame1->GetYaxis()->CenterTitle(1);
409 +  SFSample->plotOn(frame1,RooFit::Name("SFdata")) ;
410 +  
411 +  RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
412 +  frame2->GetXaxis()->CenterTitle(1);
413 +  frame2->GetYaxis()->CenterTitle(1);
414 +  OFSample->plotOn(frame2,RooFit::Name("OFdata")) ;
415 +  
416 +  TCanvas* cSFdata = new TCanvas("cSFdata","cSFdata") ;
417 +  cSFdata->cd() ;
418 +  gPad->SetLeftMargin(0.15);
419 +  frame1->GetYaxis()->SetTitleOffset(1.4);
420 +  frame1->Draw();
421 +  if(is_data==data) DrawPrelim();
422 +  else DrawPrelim(PlottingSetup::luminosity,true);
423 +  CompleteSave(cSFdata,"Edge/SF_NoFit");
424 +  
425 +  TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
426 +  cOFdata->cd() ;
427 +  gPad->SetLeftMargin(0.15);
428 +  frame2->GetYaxis()->SetTitleOffset(1.4);
429 +  frame2->Draw();
430 +  if(is_data==data) DrawPrelim();
431 +  else DrawPrelim(PlottingSetup::luminosity,true);
432 +  CompleteSave(cOFdata,"Edge/OF_NoFit");
433 + }
434 +  
435 + string WriteWithError(float central, float error, int digits) {
436 +  float ref=central;
437 +  if(ref<0) ref=-central;
438 +  int HighestSigDigit = 0;
439 +  if(ref>1) HighestSigDigit = int(log(ref)/log(10))+1;
440 +  else HighestSigDigit = int(log(ref)/(log(10)));
441 +  
442 +  float divider=pow(10.0,(double(HighestSigDigit-digits)));
443 +  
444 +  stringstream result;
445 +  result << divider*int(central/divider+0.5) << " #pm " << divider*int(error/divider+0.5);
446 +  return result.str();
447   }
448  
449 +
450   string EdgeFitter::RandomStorageFile() {
451    TRandom3 *r = new TRandom3(0);
452    int rho = (int)r->Uniform(1,10000000);
# Line 356 | Line 469 | Yield EdgeFitter::Get_T_estimate(float j
469   }
470  
471   void EdgeFitter::DoFit(int is_data, float jzb_cut) {
472 <  RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
472 >  RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
473 >  RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
474    RooCategory sample("sample","sample") ;
475 <  sample.defineType("ee");
475 >  sample.defineType("SF");
476    //sample.defineType("mm");
477 <  sample.defineType("em");
478 <  //RooDataSet combData("combData","combined data",mll,Index(sample),Import("ee",*eeSample),Import("mm",*mmSample),Import("em",*emSample));
479 <  RooDataSet combData("combData","combined data",mll,RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::Import("em",*emSample));
477 >  sample.defineType("OF");
478 >  
479 >  //RooDataSet combData("combData","combined data",mll,Index(sample),Import("SF",*SFSample),Import("mm",*mmSample),Import("OF",*OFSample));
480 >  RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("SF",*SFSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
481 >  
482    
367
368
483    //First we make a fit to opposite flavor
484 <  RooRealVar fttbarem("fttbarem", "fttbarem", 100, 0, 10000);
485 <  RooRealVar par1ttbarem("par1ttbarem", "par1ttbarem", 1.6, 0.01, 4.0);
486 <  RooRealVar par2ttbarem("par2ttbarem", "par2ttbarem", 1.0);
487 <  RooRealVar par3ttbarem("par3ttbarem", "par3ttbarem", 0.028, 0.001, 1.0);
488 <  RooRealVar par4ttbarem("par4ttbarem", "par4ttbarem", 2.0);
489 <  RooSUSYBkgPdf ttbarem("ttbarem","ttbarem", mll , par1ttbarem, par2ttbarem, par3ttbarem, par4ttbarem);
490 <  RooAddPdf model_em("model_em","model_em", ttbarem, fttbarem);
484 >  RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 10000);
485 >  RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
486 >  RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
487 >  RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
488 >  RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
489 >  RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
490 >  RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
491    RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
492 <  simPdfOF.addPdf(model_em,"em");
493 <  RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save());
494 <  resultOF->Print();
492 >  simPdfOF.addPdf(model_OF,"OF");
493 >  RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save(),RooFit::Extended(),RooFit::Minos(true));
494 >  //resultOF->Print();
495 >  
496 >  if(resultOF->covQual()!=3) {
497 >    write_error(__FUNCTION__,"OF fit did not converge!!! Cannot continue!");
498 >    cout << "covQual is " << resultOF->covQual() << endl;
499 >    EdgeFitter::FixedMEdgeChi2=-1;
500 >    if(EdgeFitter::RejectPointIfNoConvergence) return;
501 >  } else {
502 >    write_info(__FUNCTION__,"OF fit converged");
503 >  }
504  
505 <  RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarem");
505 >  RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarOF");
506    float resultOFpar1 = resultOFpar1_->getVal();
507 <  //RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarem");
507 >  //RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarOF");
508    //float resultOFpar2 = resultOFpar2_->getVal();
509    //cout << "caca2.txt" << endl;
510  
511 <  RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarem");
511 >  RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarOF");
512    float resultOFpar3 = resultOFpar3_->getVal();
513  
514 <  //RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarem");
514 >  //RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarOF");
515    //float resultOFpar4 = resultOFpar4_->getVal();
516    //cout << "caca4.txt" << endl;
517 +  
518 +  float StartingMedge=70;
519 +  if(EdgeFitter::FixedMEdge>0) StartingMedge=EdgeFitter::FixedMEdge;
520  
521  
522    // Now same flavor  
523 <  RooRealVar fzee("fzee", "fzee", 5, 0, 100000);
524 <  RooRealVar meanzee("meanzee", "meanzee", 91.1876, 89, 95);
525 <  //RooRealVar sigmazee("sigmazee", "sigmazee", 0.5);
526 <  RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0, 100);
527 <  RooRealVar widthzee("widthzee", "widthzee", 2.94);
528 <  
529 <  RooRealVar fttbaree("fttbaree", "fttbaree", 100, 0, 100000);
530 <  RooRealVar par1ttbaree("par1ttbaree", "par1ttbaree", resultOFpar1, 0, 100);
531 <  RooRealVar par2ttbaree("par2ttbaree", "par2ttbaree", 1.0);
532 <  RooRealVar par3ttbaree("par3ttbaree", "par3ttbaree", resultOFpar3, 0, 100);
533 <  RooRealVar par4ttbaree("par4ttbaree", "par4ttbaree", 2.0);
523 >  RooRealVar fzSF("fzSF", "fzSF", 5, 0, 100000);
524 >  RooRealVar meanzSF("meanzSF", "meanzSF", 91.1876, 89, 95);
525 >  //RooRealVar sigmazSF("sigmazSF", "sigmazSF", 0.5);
526 >  RooRealVar sigmazSF("sigmazSF", "sigmazSF", 5, 0.5, 5);
527 >  RooRealVar widthzSF("widthzSF", "widthzSF", 2.94);
528 >  
529 >  RooRealVar fttbarSF("fttbarSF", "fttbarSF", 100, 0, 100000);
530 >  RooRealVar par1ttbarSF("par1ttbarSF", "par1ttbarSF", 1.02*resultOFpar1, 0, 100);
531 >  RooRealVar par2ttbarSF("par2ttbarSF", "par2ttbarSF", 1.0);
532 >  RooRealVar par3ttbarSF("par3ttbarSF", "par3ttbarSF", resultOFpar3, 0, 100);
533 >  RooRealVar par4ttbarSF("par4ttbarSF", "par4ttbarSF", 2.0);
534  
535 <  RooRealVar fsignalee("fsignalee", "fsignalee", 10, 0, 400);
536 <  RooRealVar par1signalee("par1signalee", "par1signalee", 45, 20, 100);
537 <  RooRealVar par2signalee("par2signalee", "par2signalee", 2, 1, 10);
538 <  RooRealVar par3signalee("par3signalee", "par3signalee", 45, 0, 200);
539 <
540 <  RooVoigtian zee("zee", "zee", mll, meanzee, widthzee, sigmazee);
541 <
542 <  
543 <  RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbaree, par2ttbaree, par3ttbaree, par4ttbaree);
544 <  //RooSUSYTPdf signalee("signalee","signalee", mll , par1signalee, par2signalee, par3signalee);
545 <  RooSUSYTPdf signalee("signalee","signalee", mll , par1signalee, sigmazee, par3signalee);
546 <
547 <  //RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
548 <  RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
549 <  RooAddPdf model_emu("model_emu","model_emu", RooArgList(ttbaree), RooArgList(fttbaree));
535 >  RooRealVar fsignalSF("fsignalSF", "fsignalSF", 10, 0, 300);
536 >  RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
537 >  RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
538 >  RooRealVar par3signalSF("par3signalSF", "par3signalSF", StartingMedge, 0, 300);
539 >
540 >  RooVoigtian zSF("zSF", "zSF", mll, meanzSF, widthzSF, sigmazSF);
541 >
542 >  if(EdgeFitter::FixedMEdge>0) par3signalSF.setConstant();
543 >  
544 > /*  par1ttbarOF.setConstant(1);
545 >  par2ttbarOF.setConstant(1);
546 >  par3ttbarOF.setConstant(1);
547 >  par4ttbarOF.setConstant(1);
548 >  fttbarOF.setConstant(1);*/
549 >  
550 >  RooSUSYBkgPdf ttbarSF("ttbarSF","ttbarSF", mll , par1ttbarSF, par2ttbarSF, par3ttbarSF, par4ttbarSF);
551 >  //RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, par2signalSF, par3signalSF);
552 >  RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, sigmazSF, par3signalSF);
553 >  
554 > /*  par1ttbarSF.setConstant(true);
555 >  par2ttbarSF.setConstant(true);
556 >  par3ttbarSF.setConstant(true);
557 >  par4ttbarSF.setConstant(true);*/
558 >  
559 >
560 >  //RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
561 >  RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
562 > //  RooAddPdf model_OF("model_OF","model_OF", RooArgList(ttbarSF), RooArgList(fttbarSF));
563  
564    
565    RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
566 <  simPdf.addPdf(model_ee,"ee") ;
567 <  simPdf.addPdf(model_emu,"em") ;
566 >  simPdf.addPdf(model_SF,"SF") ;
567 >  simPdf.addPdf(model_OF,"OF") ;
568  
569 <  RooFitResult *result = simPdf.fitTo(combData, RooFit::Save());
431 <  result->Print();
569 >  RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended(),RooFit::Minos(true));
570    
571 <  RooPlot* frame1 = mll.frame(RooFit::Bins(25),RooFit::Title("EE sample")) ;
572 <  combData.plotOn(frame1,RooFit::Cut("sample==sample::ee")) ;
573 <  simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
574 <  simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("ttbaree"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
575 <  simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("zee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
576 <  simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("signalee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
577 <
578 <  RooPlot* frame2 = mll.frame(RooFit::Bins(25),RooFit::Title("MM sample")) ;
441 <  combData.plotOn(frame2,RooFit::Cut("sample==sample::mm")) ;
442 <  simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
443 <  simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("ttbarmm"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
444 <  simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("zmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
445 <  simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("signalmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
571 >  if(result->covQual()!=3) {
572 >    write_error(__FUNCTION__,"Full fit did not converge!!! Cannot continue!");
573 >    cout << "covQual is " << result->covQual() << endl;
574 >    EdgeFitter::FixedMEdgeChi2=-1;
575 >    if(EdgeFitter::RejectPointIfNoConvergence) return;
576 >  } else {
577 >    write_info(__FUNCTION__,"Full fit converged");
578 >  }
579  
580 + //  result->Print();
581 +  
582 +  RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
583 +  frame1->GetXaxis()->CenterTitle(1);
584 +  frame1->GetYaxis()->CenterTitle(1);
585 +  combData.plotOn(frame1,RooFit::Name("SFdata"),RooFit::Cut("sample==sample::SF")) ;
586 +  simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
587 +  simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("TTbarSFonly"),RooFit::Components("ttbarSF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
588 +  simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("DYSFonly"),RooFit::Components("zSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
589 +  simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("SignalSFonly"),RooFit::Components("signalSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
590 +  
591 +  EdgeFitter::FixedMEdgeChi2 = frame1->chiSquare("FullFit", "SFdata", 3);
592 +  
593    
594    cout << "Result   : " << endl;
595 <  cout << "f signal : " << fsignalee.getVal() << " +/- " << fsignalee.getError() << endl;
596 <  cout << "f ttbar  : " << fttbaree.getVal() << " +/- " << fttbaree.getError() << endl;
597 <  cout << "f tt em  : " << fttbarem.getVal() << " +/- " << fttbarem.getError() << endl;
598 <  cout << "f z ee   : " << fzee.getVal() << " +/- " << fzee.getError() << endl;
595 >  cout << "f signal : " << fsignalSF.getVal() << " +/- " << fsignalSF.getError() << endl;
596 >  cout << "f ttbar  : " << fttbarSF.getVal() << " +/- " << fttbarSF.getError() << endl;
597 >  cout << "f tt OF  : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
598 >  cout << "f z SF   : " << fzSF.getVal() << " +/- " << fzSF.getError() << endl;
599 >  cout << "#Chi^{2}/NDF     : " << EdgeFitter::FixedMEdgeChi2 << endl;
600    
601    // The same plot for the cointrol sample slice
602 <  RooPlot* frame3 = mll.frame(RooFit::Bins(25),RooFit::Title("EM sample")) ;
603 <  combData.plotOn(frame3,RooFit::Cut("sample==sample::em")) ;
604 <  simPdfOF.plotOn(frame3,RooFit::Slice(sample,"em"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
605 <  simPdfOF.plotOn(frame3,RooFit::Slice(sample,"em"),RooFit::Components("ttbarem"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
602 >  RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
603 >  frame3->GetXaxis()->CenterTitle(1);
604 >  frame3->GetYaxis()->CenterTitle(1);
605 >  frame3->SetMaximum(frame1->GetMaximum());
606 >  combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
607 >  simPdfOF.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
608 >  simPdfOF.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
609 >  
610    
611    stringstream prefix;
612    if(is_data==data) prefix << "data_";
613    if(is_data==mc) prefix << "mc_";
614    if(is_data==mcwithsignal) prefix << "mcwithS_";
615    
616 <  prefix << "JZB_" << jzb_cut;
616 >  prefix << EdgeFitter::Mode << "_" << jzb_cut;
617    
467 /*  cout << "fsignalee : " << fsignalee << endl;
468  cout << "fttbaree : " << fttbaree << endl;
469  cout << "fzee : " << fzee << endl;*/
618    
619    
620    TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
# Line 476 | Line 624 | void EdgeFitter::DoFit(int is_data, floa
624    frame1->Draw();
625    if(is_data==data) DrawPrelim();
626    else DrawPrelim(PlottingSetup::luminosity,true);
627 <  CompleteSave(c,"Edge/"+prefix.str()+"eemm");
627 >  stringstream infotext;
628 >  infotext << "#splitline{Fit results (" << EdgeFitter::Mode << ">" << jzb_cut << "): }{#splitline{";
629 >  infotext << "N(Data) = " << EdgeFitter::SFSample->sumEntries() << "}{#splitline{";
630 >  infotext << "N(Z+Jets) = " << WriteWithError(fzSF.getVal(),fzSF.getError(),3) << "}{#splitline{";
631 >  infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
632 >  infotext << "N(signal) = " << WriteWithError(fsignalSF.getVal(),fsignalSF.getError(),3) << "}{";
633 >  infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";  
634 >
635 >  TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
636 >  infobox->SetNDC();
637 >  infobox->SetTextSize(0.03);
638 >  infobox->Draw();
639 >  if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_SF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
640 >  else CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
641    delete c;
642    
482  TCanvas* d = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
483  d->cd() ;
484  gPad->SetLeftMargin(0.15);
485  frame2->GetYaxis()->SetTitleOffset(1.4);
486  frame2->Draw();
487  if(is_data==data) DrawPrelim();
488  else DrawPrelim(PlottingSetup::luminosity,true);
489  CompleteSave(d,"Edge/"+prefix.str()+"mm");
490  delete d;
491  //c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw();
492  
643    TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
644    e->cd();
645    gPad->SetLeftMargin(0.15);
# Line 497 | Line 647 | void EdgeFitter::DoFit(int is_data, floa
647    frame3->Draw();
648    if(is_data==data) DrawPrelim();
649    else DrawPrelim(PlottingSetup::luminosity,true);
650 <  CompleteSave(e,"Edge/"+prefix.str()+"emu");
650 >  if(EdgeFitter::FixedMEdge>=0) CompleteSave(e,"Edge/"+prefix.str()+"_OF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
651 >  else CompleteSave(e,"Edge/"+prefix.str()+"_OF",false,false);
652    delete e;
653    
654 +  
655 +  
656 +  
657   /*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
658    f->cd();
659    gPad->SetLeftMargin(0.15);
# Line 507 | Line 661 | void EdgeFitter::DoFit(int is_data, floa
661    frame4->Draw();
662    if(is_data==data) DrawPrelim();
663    else DrawPrelim(PlottingSetup::luminosity,true);
664 <  CompleteSave(f,"Edge/"+prefix.str()+"eemm");
664 >  CompleteSave(f,"Edge/"+prefix.str()+"_SF");
665    delete f;*/
666 <  
666 >
667 >
668 > /*
669 >  float maxZ=200;
670 >  RooWorkspace* wspace = new RooWorkspace();
671 >  stringstream mllvar;
672 >  mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
673 >  wspace->factory(mllvar.str().c_str());
674 >  wspace->var("mll")->setBins(30);
675 >  wspace->factory("nSig[1.,0.,100.]");
676 >  wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
677 >  wspace->factory("rME[1.12,1.05,1.19]");
678 >  wspace->factory("effUncert[1.]");
679 >  EdgeFitter::prepareLimits(wspace, true);
680 > */
681 >
682 > write_warning(__FUNCTION__," A lot missing here to calculate limits");
683 >
684   }
685  
686   void EdgeFitter::DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree *signalevents=0) {
# Line 522 | Line 693 | void EdgeFitter::DoEdgeFit(string mcjzb,
693    
694    EdgeFitter::PrepareDatasets(is_data);
695    
696 +  EdgeFitter::DrawDatasetContent(is_data);
697 +  
698    RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
699    RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
700 +  
701 +  
702 +  bool ScanMassRange=false;
703 +  
704 +  
705 +  
706 +  if(ScanMassRange) {
707 +    TFile *fscan = new TFile("fscan.root","UPDATE");
708 +    TGraph *gr = new TGraph();
709 +    stringstream GrName;
710 +    GrName << "ScanGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
711 +    gr->SetName(GrName.str().c_str());
712 +
713 +    int i=0;
714 +    for(float tempMedge=10;tempMedge<=300;tempMedge+=5.0) {
715 +      write_info(__FUNCTION__,"Now testing Medge="+any2string(tempMedge)+" for "+EdgeFitter::Mode+">"+any2string(jzb_cut));
716 +      EdgeFitter::FixedMEdge=tempMedge;
717 +      EdgeFitter::DoFit(is_data, jzb_cut);
718 +      if(EdgeFitter::FixedMEdgeChi2<0) continue;
719 +      gr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2);
720 +      i++;
721 +    }
722 +    
723 +    TCanvas *ScanCan = new TCanvas("ScanCan","ScanCan",500,500);
724 +    gr->GetXaxis()->SetTitle("m_{edge}");
725 +    gr->GetXaxis()->CenterTitle();
726 +    gr->GetYaxis()->SetTitle("#Chi^{2} / NDF");
727 +    gr->GetYaxis()->CenterTitle();
728 +    gr->GetYaxis()->SetTitleOffset(0.95);
729 +    gr->GetXaxis()->SetTitleOffset(0.9);
730 +    gr->SetLineColor(kBlue);
731 +    gr->SetTitle("");
732 +    gr->Draw("AL");
733 +    stringstream ScanCanSave;
734 +    ScanCanSave << "Edge/MEdgeScan_"+EdgeFitter::Mode+"_" << jzb_cut;
735 +    if(is_data) DrawPrelim();
736 +    else DrawMCPrelim();
737 +    CompleteSave(ScanCan,ScanCanSave.str());
738 +    fscan->cd();
739 +    gr->Write();
740 +    delete ScanCan;
741 +    fscan->Close();
742 +  } else {
743 +    EdgeFitter::DoFit(is_data, jzb_cut);
744 +  }
745 +  
746 +      
747 +  
748 +  
749    EdgeFitter::DoFit(is_data, jzb_cut);
750    RooMsgService::instance().setGlobalKillBelow(msglevel);
751  
# Line 533 | Line 755 | void EdgeFitter::DoEdgeFit(string mcjzb,
755   }
756  
757   void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
758 <  for(int icut=0;icut<jzb_cut.size();icut++) {
758 >  
759 >  EdgeFitter::Mode="JZB";
760 >  if(mcjzb=="met[4]") EdgeFitter::Mode="MET";
761 >  
762 >  for(int icut=0;icut<(int)jzb_cut.size();icut++) {
763      stringstream addcut;
764      if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
765      if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines