ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C
Revision: 1.21
Committed: Fri Jun 14 15:24:32 2013 UTC (11 years, 10 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.20: +82 -27 lines
Log Message: 
Commiting TEST (stating that again: TEST) version of code allowing H1 and H0 testing (which is done) and R(SF/OF), which is still experimental; we now have a gaussian distribution centered at 1.02 with witdth 0.07 which is multiplied with the integral from OF to get the flavor-symmetric contribution in SF; also added additional plot to illustrate CHi2(H1)/Chi2(H0) for mass scans

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2    
3 buchmann 1.21
4    
5 buchmann 1.10 #include <TVirtualIndex.h>
6    
7 buchmann 1.2 #include <RooRealVar.h>
8     #include <RooArgSet.h>
9     #include <RooDataSet.h>
10 buchmann 1.4 #include <RooMCStudy.h>
11     #include <RooCategory.h>
12    
13 buchmann 1.5 #include <RooPlot.h>
14 buchmann 1.21 #include <RooGaussian.h>
15     #include <RooConstVar.h>
16 buchmann 1.5 #include <RooSimultaneous.h>
17     #include <RooAddPdf.h>
18     #include <RooFitResult.h>
19     #include <RooVoigtian.h>
20     #include <RooMsgService.h>
21    
22 buchmann 1.3 #include <RooStats/ModelConfig.h>
23     #include "RooStats/ProfileLikelihoodCalculator.h"
24     #include "RooStats/LikelihoodInterval.h"
25     #include "RooStats/HypoTestResult.h"
26     #include "RooStats/SimpleLikelihoodRatioTestStat.h"
27     #include "RooStats/ProfileLikelihoodTestStat.h"
28     #include "RooStats/HybridCalculatorOriginal.h"
29     #include "RooStats/HypoTestInverterOriginal.h"
30    
31 buchmann 1.5 //#include "ParametrizedEdge.C"
32 buchmann 1.8 #include "EdgeModules/RooSUSYTPdf.cxx"
33     #include "EdgeModules/RooSUSYBkgPdf.cxx"
34 buchmann 1.5
35 buchmann 1.2
36 buchmann 1.1 using namespace std;
37     using namespace PlottingSetup;
38    
39    
40 buchmann 1.2
41    
42 buchmann 1.1 ShapeDroplet LimitsFromEdge(float low_fullCLs, float high_CLs, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
43     write_error(__FUNCTION__,"Not implemented edge limits yet (returning crap)");
44     ShapeDroplet beta;beta.observed=-12345;beta.SignalIntegral=1;return beta;
45     }
46    
47    
48     void PrepareEdgeShapes(string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
49     write_error(__FUNCTION__,"Not implemented edge shape storage yet");
50     }
51    
52    
53 buchmann 1.2 ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
54    
55    
56     namespace EdgeFitter {
57    
58     void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree*);
59     void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree*);
60 buchmann 1.7 void getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
61 buchmann 1.2 void InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut);
62     void PrepareDatasets(int);
63 buchmann 1.15 void DrawDatasetContent(int);
64 buchmann 1.5 void DoFit(int is_data, float jzb_cut);
65 buchmann 1.2 string RandomStorageFile();
66     Yield Get_Z_estimate(float,int);
67     Yield Get_T_estimate(float,int);
68 buchmann 1.7 float calcExclusion(RooWorkspace *ws, RooDataSet data, bool calcExclusion);
69     vector<RooDataSet> generateToys(RooWorkspace *ws, int nToys);
70 buchmann 1.4 void prepareLimits(RooWorkspace *ws, bool ComputeBands);
71     TGraph* prepareLM(float mass, float nEv);
72 buchmann 1.2
73     float jzbmax;
74     float mllmin;
75     float mllmax;
76     TCut cut;
77    
78     RooDataSet* AllData;
79 buchmann 1.10 RooDataSet* SFSample;
80     RooDataSet* OFSample;
81 buchmann 1.2
82 buchmann 1.6 bool MarcoDebug=true;
83 buchmann 1.11
84     float FixedMEdge=-1;
85 buchmann 1.21 float FixedMEdgeChi2_H1=-1;
86     float FixedMEdgeChi2_H0=-1;
87 buchmann 1.11
88 buchmann 1.14 bool RejectPointIfNoConvergence=false;
89    
90     string Mode="UndefinedMode";
91    
92 buchmann 1.21 bool AllowTriangle=true;
93 buchmann 1.2 }
94    
95 buchmann 1.18 TLatex* WriteAppleStudel() {
96     string sel="Apple Strudel Preliminary";
97     TLatex *sele = new TLatex(0.97,0.135,sel.c_str());
98     sele->SetNDC(true);
99     sele->SetTextColor(TColor::GetColor("#848484"));
100     sele->SetTextFont(42);
101     sele->SetTextAlign(32);
102     sele->SetTextSize(0.03);
103     sele->SetTextAngle(270);
104     return sele;
105     }
106    
107 buchmann 1.4 TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
108     float massLM[1];
109     massLM[0]=mass;
110     float accEffLM[1];
111     accEffLM[0]=nEv/PlottingSetup::luminosity;
112     TGraph *lm = new TGraph(1, massLM, accEffLM);
113     lm->GetXaxis()->SetNoExponent(true);
114     lm->GetXaxis()->SetTitle("m_{cut} [GeV]");
115     lm->GetYaxis()->SetTitle("#sigma #times A [pb] 95% CL UL");
116     lm->GetXaxis()->SetLimits(0.,300.);
117     lm->GetYaxis()->SetRangeUser(0.,0.08);
118     lm->SetMarkerStyle(34);
119     lm->SetMarkerColor(kRed);
120     return lm;
121     }
122    
123 buchmann 1.7 vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
124 buchmann 1.10 ws->ls();
125 buchmann 1.4 ws->var("nSig")->setVal(0.);
126     ws->var("nSig")->setConstant(true);
127     RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
128 buchmann 1.7 vector<RooDataSet> theToys;
129 buchmann 1.4
130     RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
131     mcEE.generate(nToys,44,true);
132     RooMCStudy mcMM(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"MM"));
133     mcMM.generate(nToys,44,true);
134     RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
135     mcOSOF.generate(nToys,44,true);
136    
137 buchmann 1.10 RooRealVar mll("m_{ll}","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
138 buchmann 1.4 RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
139     RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
140     RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
141     RooRealVar weight("weight","weight",0,1000,"");
142     RooArgSet observables(mll,jzb,id1,id2,weight);
143    
144     for(int i=0;i<nToys;i++) {
145     RooDataSet* toyEE = (RooDataSet*)mcEE.genData(i);
146     RooDataSet* toyMM = (RooDataSet*)mcMM.genData(i);
147     RooDataSet* toyOSOF = (RooDataSet*)mcOSOF.genData(i);
148     stringstream toyname;
149     toyname << "theToy_" << i;
150     write_warning(__FUNCTION__,"Problem while adding toys");
151 buchmann 1.7 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));
152     theToys.push_back(toyData);
153 buchmann 1.4 }
154     ws->var("nSig")->setVal(17.0);
155     ws->var("nSig")->setConstant(false);
156     return theToys;
157     }
158    
159 buchmann 1.7 void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
160 buchmann 1.4 TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
161     vector<float> theLimits;
162 buchmann 1.9 for(int itoy=0;itoy<(int)theToys.size();itoy++) {
163 buchmann 1.7 float theLimit = calcExclusion(ws,theToys[itoy],false);
164 buchmann 1.4 if(theLimit > 0 ) gauLimit->Fill(theLimit);
165     }
166     const Int_t nQ = 4;
167     Double_t yQ[nQ] = {0.,0.,0.,0.};
168     Double_t xQ[nQ] = {0.022750132,0.1586552539,0.84134474609999998,0.977249868};
169     gauLimit->GetQuantiles(nQ,yQ,xQ);
170     median = gauLimit->GetMean();
171     // median = median1(gauLimit);
172     twoSigmaDown = abs(yQ[0]-median);
173     sigmaDown = abs(yQ[1]-median);
174     sigmaUp = abs(yQ[2]-median);
175     twoSigmaUp = abs(yQ[3]-median);
176 buchmann 1.18 dout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
177 buchmann 1.4 }
178    
179     void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
180     if(ComputeBands) {
181 buchmann 1.7 vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
182 buchmann 1.4 vector<float> medVals;
183     vector<float> medLimits;
184     vector<float> sigmaLimitsDown;
185     vector<float> sigmaLimitsUp;
186     vector<float> twoSigmaLimitsDown;
187     vector<float> twoSigmaLimitsUp;
188     for(int i=20;i<=320;i+=40) {
189     ws->var("nSig")->setVal(10.0);
190     medVals.push_back((float)i);
191     ws->var("m0")->setVal((float)i);
192     ws->var("m0")->setConstant(true);
193     float Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp;
194     EdgeFitter::getMedianLimit(ws,theToys,Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp);
195     medLimits.push_back(Smedian);
196     sigmaLimitsDown.push_back(SsigmaDown);
197     sigmaLimitsUp.push_back(SsigmaUp);
198     twoSigmaLimitsDown.push_back(StwoSigmaDown);
199     twoSigmaLimitsUp.push_back(StwoSigmaUp);
200     }
201     write_warning(__FUNCTION__,"Still need to store limits");
202     } else {
203     vector<float> theVals;
204     vector<float> theLimits;
205     for(int i=20;i<300;i+=5) {
206     ws->var("nSig")->setVal(0.0);
207     theVals.push_back((float)i);
208     ws->var("m0")->setVal((float)i);
209     ws->var("m0")->setConstant(true);
210 buchmann 1.7 // theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
211     write_error(__FUNCTION__,"Error while casting roo data set");
212 buchmann 1.4 }
213    
214 buchmann 1.7 for(int i=0;i<(int)theLimits.size();i++) {
215 buchmann 1.4 if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
216 buchmann 1.18 dout << i << " : " << theVals[i] << endl;
217 buchmann 1.4 theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
218     write_warning(__FUNCTION__,"Need to store limits");
219     }
220     write_warning(__FUNCTION__,"Need to store limits");
221     }
222     }
223     }
224    
225    
226 buchmann 1.7 float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
227     int numberOfToys=50;
228 buchmann 1.3 RooRealVar mu("mu","nSig",0,10000,"");
229     RooArgSet poi = RooArgSet(mu);
230     RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
231     nullParams->setRealValue("nSig",0);
232     RooStats::ModelConfig *model = new RooStats::ModelConfig();
233     model->SetWorkspace(*ws);
234     model->SetPdf("combModel");
235     model->SetParametersOfInterest(poi);
236 buchmann 1.7 // if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
237 buchmann 1.3
238 buchmann 1.7 RooStats::ProfileLikelihoodCalculator plc(data, *model);
239 buchmann 1.3 plc.SetNullParameters(*nullParams);
240     plc.SetTestSize(0.05);
241 buchmann 1.7
242 buchmann 1.3 RooStats::LikelihoodInterval* interval = plc.GetInterval();
243     RooStats::HypoTestResult *htr = plc.GetHypoTest();
244     double theLimit = interval->UpperLimit( mu );
245 buchmann 1.7 // double significance = htr->Significance();
246 buchmann 1.3
247     ws->defineSet("obs","nB");
248     ws->defineSet("poi","nSig");
249    
250     RooStats::ModelConfig b_model = RooStats::ModelConfig("B_model", ws);
251     b_model.SetPdf(*ws->pdf("combModel"));
252     b_model.SetObservables(*ws->set("obs"));
253     b_model.SetParametersOfInterest(*ws->set("poi"));
254     ws->var("nSig")->setVal(0.0); //# important!
255     b_model.SetSnapshot(*ws->set("poi"));
256    
257     RooStats::ModelConfig sb_model = RooStats::ModelConfig("S+B_model", ws);
258     sb_model.SetPdf(*ws->pdf("combModel"));
259     sb_model.SetObservables(*ws->set("obs"));
260     sb_model.SetParametersOfInterest(*ws->set("poi"));
261     ws->var("nSig")->setVal(64.0); //# important!
262     sb_model.SetSnapshot(*ws->set("poi"));
263    
264     RooStats::SimpleLikelihoodRatioTestStat slrts = RooStats::SimpleLikelihoodRatioTestStat(*b_model.GetPdf(),*sb_model.GetPdf());
265     slrts.SetNullParameters(*b_model.GetSnapshot());
266     slrts.SetAltParameters(*sb_model.GetSnapshot());
267     RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
268    
269 buchmann 1.7 RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
270 buchmann 1.3 hc.SetTestStatistic(2);
271 buchmann 1.7 hc.SetNumberOfToys(numberOfToys);
272 buchmann 1.3
273     RooStats::HypoTestInverterOriginal hcInv = RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
274     hcInv.SetTestSize(0.05);
275     hcInv.UseCLs(true);
276     hcInv.RunFixedScan(5,theLimit-0.5,theLimit+0.5);;
277     RooStats::HypoTestInverterResult* hcInt = hcInv.GetInterval();
278     float ulError = hcInt->UpperLimitEstimatedError();
279     theLimit = hcInt->UpperLimit();
280 buchmann 1.18 dout << "Found upper limit : " << theLimit << "+/-" << ulError << endl;
281 buchmann 1.3
282     return theLimit/PlottingSetup::luminosity;
283    
284     }
285    
286 buchmann 1.2 TTree* SkimTree(int isample) {
287 buchmann 1.20 string TreeName = removefunnystring(allsamples.collection[isample].filename);
288     cout << "About to skim " << TreeName << " for sample id " << isample << endl;
289 buchmann 1.15 TTree* newTree = new TTree("NanoTree",TreeName.c_str());
290    
291     float mll,edgeWeight;
292     int id1,id2;
293    
294     newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
295     newTree->Branch("mll",&mll,"mll/F");
296     newTree->Branch("id1",&id1,"id1/I");
297     newTree->Branch("id2",&id2,"id2/I");
298 buchmann 1.20 newTree->Branch("isample",&isample,"isample/I");
299 buchmann 1.15
300 buchmann 1.2 float xsweight=1.0;
301     if(allsamples.collection[isample].is_data==false) xsweight=luminosity*allsamples.collection[isample].weight;
302     if(EdgeFitter::MarcoDebug) {
303 buchmann 1.20 dout << " Going to reduce it with cut " << EdgeFitter::cut << endl;
304 buchmann 1.18 dout << " Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
305 buchmann 1.2 }
306 buchmann 1.15
307 buchmann 1.6 float tmll;
308 buchmann 1.15 int tid1,tid2;
309 buchmann 1.6 allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
310 buchmann 1.15 allsamples.collection[isample].events->SetBranchAddress("id1",&tid1);
311     allsamples.collection[isample].events->SetBranchAddress("id2",&tid2);
312 buchmann 1.6
313 buchmann 1.2 TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
314 buchmann 1.10 TTreeFormula *Weight = new TTreeFormula("Weight", cutWeight, allsamples.collection[isample].events);
315 buchmann 1.15
316 buchmann 1.2 float wgt=1.0;
317     for (Int_t entry = 0 ; entry < allsamples.collection[isample].events->GetEntries() ; entry++) {
318     allsamples.collection[isample].events->LoadTree(entry);
319     if (select->EvalInstance()) {
320     allsamples.collection[isample].events->GetEntry(entry);
321 buchmann 1.15 mll=tmll;
322     id1=tid1;
323     id2=tid2;
324 buchmann 1.10 wgt=Weight->EvalInstance();
325 buchmann 1.6 edgeWeight=wgt*xsweight;
326 buchmann 1.2 newTree->Fill();
327     }
328     }
329    
330 buchmann 1.20 if(EdgeFitter::MarcoDebug) dout << " Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
331 buchmann 1.2 return newTree;
332     }
333    
334     void EdgeFitter::InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut) {
335     mllmin=_mllmin;
336     mllmax=_mllmax;
337     jzbmax=_jzbmax;
338     cut=_cut;
339     }
340 buchmann 1.10
341     TTree* MergeTrees(vector<TTree*> trees) {
342 buchmann 1.15 TTree * newtree = (TTree*)trees[0]->CloneTree(0);
343 buchmann 1.20 newtree->SetTitle("FullTree");
344     newtree->SetName("FullTree");
345 buchmann 1.10 trees[0]->GetListOfClones()->Remove(newtree);
346     trees[0]->ResetBranchAddresses();
347     newtree->ResetBranchAddresses();
348    
349 buchmann 1.15 for(int itree=0;itree<trees.size();itree++) {
350 buchmann 1.10 newtree->CopyAddresses(trees[itree]);
351     Long64_t nentries = trees[itree]->GetEntries();
352     for (Long64_t iEntry=0;iEntry<nentries;iEntry++) {
353     trees[itree]->GetEntry(iEntry);
354     newtree->Fill();
355     }
356     trees[itree]->ResetBranchAddresses(); // Disconnect from new tree
357     if (newtree->GetTreeIndex()) {
358     newtree->GetTreeIndex()->Append(trees[itree]->GetTreeIndex(),kTRUE);
359     }
360     if (newtree && newtree->GetTreeIndex()) {
361     newtree->GetTreeIndex()->Append(0,kFALSE); // Force the sorting
362     }
363     }
364     return newtree;
365     }
366    
367    
368    
369 buchmann 1.2 void EdgeFitter::PrepareDatasets(int is_data) {
370 buchmann 1.5 write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
371 buchmann 1.15 TFile *tempout = new TFile("tempout.root","RECREATE");
372 buchmann 1.10 vector<TTree*> SkimmedTrees;
373     TTree *SkimmedTree[(int)allsamples.collection.size()];
374 buchmann 1.7 for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
375 buchmann 1.2 if(!allsamples.collection[isample].is_active) continue;
376     if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
377     if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
378     if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
379     if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
380 buchmann 1.18 if(EdgeFitter::MarcoDebug) dout << "Considering : " << allsamples.collection[isample].samplename << endl;
381 buchmann 1.15 SkimmedTree[isample] = SkimTree(isample);
382     tempout->cd();
383     SkimmedTree[isample]->Write();
384     SkimmedTrees.push_back(SkimmedTree[isample]);
385 buchmann 1.2 }
386 buchmann 1.10
387     TTree *completetree = MergeTrees(SkimmedTrees);
388    
389 buchmann 1.18 if(EdgeFitter::MarcoDebug) dout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
390 buchmann 1.2
391 buchmann 1.10 RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
392 buchmann 1.2 RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
393     RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
394 buchmann 1.6 RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
395     RooArgSet observables(mll,id1,id2,edgeWeight);
396 buchmann 1.2
397     string title="CMS Data";
398     if(is_data!=1) title="CMS SIMULATION";
399 buchmann 1.6 RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
400 buchmann 1.15 tempout->cd();
401 buchmann 1.2 completetree->Write();
402 buchmann 1.15 tempout->Close();
403 buchmann 1.2
404 buchmann 1.10 EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
405     EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
406 buchmann 1.2 EdgeFitter::AllData = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
407    
408 buchmann 1.10 SFSample->SetName("SFSample");
409     OFSample->SetName("OFSample");
410 buchmann 1.2 AllData->SetName("AllData");
411    
412     if(EdgeFitter::MarcoDebug) {
413 buchmann 1.20 dout << "Number of (weighted) events in full sample = " << AllData->sumEntries() << " (raw events : " << AllData->numEntries() << ")" << endl;
414     dout << "Number of (weighted) events in eemm sample = " << SFSample->sumEntries() << " (raw events : " << SFSample->numEntries() << ")" << endl;
415     dout << "Number of (weighted) events in em sample = " << OFSample->sumEntries() << " (raw events : " << OFSample->numEntries() << ")" << endl;
416 buchmann 1.2 }
417 buchmann 1.6
418 buchmann 1.2 }
419    
420 buchmann 1.15 void EdgeFitter::DrawDatasetContent(int is_data) {
421 buchmann 1.20 TCanvas* cSFdata = new TCanvas("cSFdata","cSFdata") ;
422    
423 buchmann 1.15 RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
424     RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("SF sample")) ;
425     frame1->GetXaxis()->CenterTitle(1);
426     frame1->GetYaxis()->CenterTitle(1);
427     SFSample->plotOn(frame1,RooFit::Name("SFdata")) ;
428    
429     RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
430     frame2->GetXaxis()->CenterTitle(1);
431     frame2->GetYaxis()->CenterTitle(1);
432     OFSample->plotOn(frame2,RooFit::Name("OFdata")) ;
433    
434 buchmann 1.20 // THStack sf_ref = allsamples.DrawStack("sf_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1==id2"),is_data, luminosity);
435     // THStack of_ref = allsamples.DrawStack("of_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1!=id2"),is_data, luminosity);
436    
437     TH1F *sf_ref = allsamples.Draw("sf_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1==id2"),is_data, luminosity);
438     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);
439    
440     TH1F *sf_fit = (TH1F*)SFSample->createHistogram("sf_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
441     TH1F *of_fit = (TH1F*)OFSample->createHistogram("of_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
442    
443     TLegend *leg = allsamples.allbglegend();
444     leg->AddEntry(sf_fit,"RooDataSet content","p");
445     leg->SetX1(0.58);
446    
447 buchmann 1.15 cSFdata->cd() ;
448     gPad->SetLeftMargin(0.15);
449     frame1->GetYaxis()->SetTitleOffset(1.4);
450     frame1->Draw();
451 buchmann 1.20 sf_ref->Draw("histo,same");
452     sf_fit->Draw("e1,same");
453     leg->Draw("same");
454 buchmann 1.15 if(is_data==data) DrawPrelim();
455     else DrawPrelim(PlottingSetup::luminosity,true);
456     CompleteSave(cSFdata,"Edge/SF_NoFit");
457    
458     TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
459     cOFdata->cd() ;
460     gPad->SetLeftMargin(0.15);
461 buchmann 1.16 frame2->SetMaximum(frame1->GetMaximum());
462 buchmann 1.15 frame2->GetYaxis()->SetTitleOffset(1.4);
463     frame2->Draw();
464 buchmann 1.20 of_ref->Draw("histo,same");
465     of_fit->Draw("e1,same");
466     leg->Draw("same");
467 buchmann 1.15 if(is_data==data) DrawPrelim();
468     else DrawPrelim(PlottingSetup::luminosity,true);
469     CompleteSave(cOFdata,"Edge/OF_NoFit");
470 buchmann 1.20
471     delete sf_fit;
472     delete of_fit;
473    
474     delete sf_ref;
475     delete of_ref;
476    
477 buchmann 1.15 }
478    
479 buchmann 1.10 string WriteWithError(float central, float error, int digits) {
480     float ref=central;
481     if(ref<0) ref=-central;
482     int HighestSigDigit = 0;
483     if(ref>1) HighestSigDigit = int(log(ref)/log(10))+1;
484     else HighestSigDigit = int(log(ref)/(log(10)));
485    
486     float divider=pow(10.0,(double(HighestSigDigit-digits)));
487    
488     stringstream result;
489     result << divider*int(central/divider+0.5) << " #pm " << divider*int(error/divider+0.5);
490     return result.str();
491     }
492    
493    
494 buchmann 1.2 string EdgeFitter::RandomStorageFile() {
495     TRandom3 *r = new TRandom3(0);
496     int rho = (int)r->Uniform(1,10000000);
497     stringstream RandomFile;
498     RandomFile << PlottingSetup::cbafbasedir << "/exchange/TempEdgeFile_" << rho << ".root";
499     delete r;
500     return RandomFile.str();
501     }
502    
503     Yield EdgeFitter::Get_Z_estimate(float jzb_cut, int icut) {
504     if(MarcoDebug) write_error(__FUNCTION__,"Returning random Z yield");
505     Yield a(123,45,67); return a;
506     return PlottingSetup::allresults.predictions[icut].Zbkg;
507     }
508    
509     Yield EdgeFitter::Get_T_estimate(float jzb_cut, int icut) {
510     if(MarcoDebug) write_error(__FUNCTION__,"Returning random TTbar yield");
511     Yield a(1234,56,78); return a;
512     return PlottingSetup::allresults.predictions[icut].Flavorsym;
513     }
514    
515 buchmann 1.5 void EdgeFitter::DoFit(int is_data, float jzb_cut) {
516 buchmann 1.16
517     stringstream prefix;
518     if(is_data==data) prefix << "data_";
519     if(is_data==mc) prefix << "mc_";
520     if(is_data==mcwithsignal) prefix << "mcwithS_";
521    
522     prefix << EdgeFitter::Mode << "_" << jzb_cut;
523    
524 buchmann 1.21 if(EdgeFitter::AllowTriangle) prefix << "__H1";
525     else prefix << "__H0";
526 buchmann 1.16
527 buchmann 1.10 RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
528 buchmann 1.6 RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
529 buchmann 1.5 RooCategory sample("sample","sample") ;
530 buchmann 1.10 sample.defineType("SF");
531 buchmann 1.5 //sample.defineType("mm");
532 buchmann 1.10 sample.defineType("OF");
533 buchmann 1.6
534 buchmann 1.10 //RooDataSet combData("combData","combined data",mll,Index(sample),Import("SF",*SFSample),Import("mm",*mmSample),Import("OF",*OFSample));
535     RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("SF",*SFSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
536 buchmann 1.6
537 buchmann 1.5
538     //First we make a fit to opposite flavor
539 buchmann 1.21 RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 1.5*combData.sumEntries());
540 buchmann 1.10 RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
541     RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
542     RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
543     RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
544     RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
545     RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
546 buchmann 1.5 RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
547 buchmann 1.10 simPdfOF.addPdf(model_OF,"OF");
548 buchmann 1.14 RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save(),RooFit::Extended(),RooFit::Minos(true));
549 buchmann 1.18 dout << "============================= < OF results > =============================" << endl;
550 buchmann 1.17 resultOF->Print();
551 buchmann 1.18 dout << "============================= < /OF results > =============================" << endl;
552 buchmann 1.14
553 buchmann 1.16
554     RooPlot* frameO = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample"));
555     frameO->GetXaxis()->CenterTitle(1);
556     frameO->GetYaxis()->CenterTitle(1);
557     combData.plotOn(frameO,RooFit::Name("OFdata"),RooFit::Cut("sample==sample::OF")) ;
558     simPdfOF.plotOn(frameO,RooFit::Slice(sample,"OF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
559     simPdfOF.plotOn(frameO,RooFit::Slice(sample,"OF"),RooFit::Name("TTbarOFonly"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
560    
561     TCanvas* pof = new TCanvas("pof","pof") ;
562     pof->cd() ;
563     gPad->SetLeftMargin(0.15);
564     frameO->GetYaxis()->SetTitleOffset(1.4);
565     frameO->Draw();
566     if(is_data==data) DrawPrelim();
567     else DrawPrelim(PlottingSetup::luminosity,true);
568 buchmann 1.17 if(EdgeFitter::FixedMEdge>=0) CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str()+"__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
569     else CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str(),false,false);
570 buchmann 1.16 delete pof;
571    
572 buchmann 1.14 if(resultOF->covQual()!=3) {
573     write_error(__FUNCTION__,"OF fit did not converge!!! Cannot continue!");
574 buchmann 1.18 dout << "covQual is " << resultOF->covQual() << endl;
575 buchmann 1.21 if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
576     else EdgeFitter::FixedMEdgeChi2_H0=-1;
577 buchmann 1.14 if(EdgeFitter::RejectPointIfNoConvergence) return;
578     } else {
579     write_info(__FUNCTION__,"OF fit converged");
580     }
581 buchmann 1.11
582     float StartingMedge=70;
583     if(EdgeFitter::FixedMEdge>0) StartingMedge=EdgeFitter::FixedMEdge;
584 buchmann 1.5
585    
586 buchmann 1.18 RooDataSet *ZDataSet = (RooDataSet*)EdgeFitter::AllData->reduce("id1==id2 && abs(mll-91.2)<20");
587    
588 buchmann 1.17 float maxZ = ZDataSet->sumEntries();
589 buchmann 1.18 dout << "maxZ was set to " << maxZ << endl;
590 buchmann 1.17 delete ZDataSet;
591    
592    
593 buchmann 1.5 // Now same flavor
594 buchmann 1.17 RooRealVar fzSF("fzSF", "fzSF", 39, 39, maxZ);
595 buchmann 1.10 RooRealVar meanzSF("meanzSF", "meanzSF", 91.1876, 89, 95);
596     //RooRealVar sigmazSF("sigmazSF", "sigmazSF", 0.5);
597     RooRealVar sigmazSF("sigmazSF", "sigmazSF", 5, 0.5, 5);
598     RooRealVar widthzSF("widthzSF", "widthzSF", 2.94);
599 buchmann 1.21 widthzSF.setConstant(1);
600    
601     RooRealVar fttbarSFx("fttbarSFx","fttbarSFx",0.2,1.8);
602     RooRealVar mRsfof("mRsfof","mRsfof",1.02);
603     RooRealVar wRsfof("wRsfof","wRsfof",0.07);
604    
605     RooGaussian RSFOF("RSFOF","RSFOF",fttbarSFx,mRsfof,wRsfof);
606     RooProdPdf fttbarSF("fttbarSF","fttbarOF*RSFOF",RooArgList(RSFOF,fttbarOF));
607 buchmann 1.10
608 buchmann 1.17 RooRealVar par1ttbarSF("par1ttbarSF", "par1ttbarSF", 1.02*par1ttbarOF.getVal(), (1.02-0.07)*par1ttbarOF.getVal(), (1.02+0.07)*par1ttbarOF.getVal());
609 buchmann 1.5
610 buchmann 1.18 RooRealVar fsignalSF("fsignalSF", "fsignalSF", 0, 0, 300);
611 buchmann 1.10 RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
612     RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
613 buchmann 1.11 RooRealVar par3signalSF("par3signalSF", "par3signalSF", StartingMedge, 0, 300);
614 buchmann 1.5
615 buchmann 1.10 RooVoigtian zSF("zSF", "zSF", mll, meanzSF, widthzSF, sigmazSF);
616 buchmann 1.5
617 buchmann 1.11 if(EdgeFitter::FixedMEdge>0) par3signalSF.setConstant();
618 buchmann 1.5
619 buchmann 1.17 RooSUSYBkgPdf ttbarSF("ttbarSF","ttbarSF", mll , par1ttbarSF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
620 buchmann 1.10 RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, sigmazSF, par3signalSF);
621 buchmann 1.21
622     RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
623    
624 buchmann 1.10
625 buchmann 1.21 if(!EdgeFitter::AllowTriangle) {
626     fsignalSF.setVal(0.0); // kill off the signal if we don't want the triangle
627     fsignalSF.setConstant(1);
628     par2signalSF.setConstant(1);
629     par3signalSF.setConstant(1);
630     }
631 buchmann 1.5
632     RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
633 buchmann 1.10 simPdf.addPdf(model_SF,"SF") ;
634 buchmann 1.14 simPdf.addPdf(model_OF,"OF") ;
635 buchmann 1.5
636 buchmann 1.21
637 buchmann 1.14 RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended(),RooFit::Minos(true));
638    
639     if(result->covQual()!=3) {
640     write_error(__FUNCTION__,"Full fit did not converge!!! Cannot continue!");
641 buchmann 1.18 dout << "covQual is " << result->covQual() << endl;
642 buchmann 1.21 if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
643     else EdgeFitter::FixedMEdgeChi2_H0=-1;
644 buchmann 1.14 if(EdgeFitter::RejectPointIfNoConvergence) return;
645     } else {
646     write_info(__FUNCTION__,"Full fit converged");
647     }
648    
649 buchmann 1.21 dout << "============================= < Full results > =============================" << endl;
650 buchmann 1.17 result->Print();
651 buchmann 1.21 dout << "============================= < /Full results > =============================" << endl;
652 buchmann 1.17
653 buchmann 1.5
654 buchmann 1.10 RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
655     frame1->GetXaxis()->CenterTitle(1);
656 buchmann 1.14 frame1->GetYaxis()->CenterTitle(1);
657 buchmann 1.11 combData.plotOn(frame1,RooFit::Name("SFdata"),RooFit::Cut("sample==sample::SF")) ;
658     simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
659     simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("TTbarSFonly"),RooFit::Components("ttbarSF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
660     simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("DYSFonly"),RooFit::Components("zSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
661 buchmann 1.21 if(EdgeFitter::AllowTriangle) simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("SignalSFonly"),RooFit::Components("signalSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
662 buchmann 1.5
663 buchmann 1.21 if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=frame1->chiSquare("FullFit", "SFdata", 3);
664     else EdgeFitter::FixedMEdgeChi2_H0=frame1->chiSquare("FullFit", "SFdata", 3);
665    
666 buchmann 1.18 dout << "Result : " << endl;
667 buchmann 1.21 if(EdgeFitter::AllowTriangle) dout << "f signal : " << fsignalSF.getVal() << " +/- " << fsignalSF.getError() << endl;
668     // dout << "f ttbar : " << fttbarSF.getVal() << " +/- " << fttbarSF.getError() << endl;
669     dout << "f ttbar : " << fttbarSF.getVal() << " +/- NO ERROR CUZ ITS A PDF "<< endl;
670 buchmann 1.18 dout << "f tt OF : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
671     dout << "f z SF : " << fzSF.getVal() << " +/- " << fzSF.getError() << endl;
672 buchmann 1.21 if(EdgeFitter::AllowTriangle) dout << "#Chi^{2}/NDF : " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
673     else dout << "#Chi^{2}/NDF : " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
674 buchmann 1.5
675     // The same plot for the cointrol sample slice
676 buchmann 1.10 RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
677     frame3->GetXaxis()->CenterTitle(1);
678 buchmann 1.14 frame3->GetYaxis()->CenterTitle(1);
679 buchmann 1.10 frame3->SetMaximum(frame1->GetMaximum());
680     combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
681 buchmann 1.16 simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
682     simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
683 buchmann 1.5
684    
685     TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
686     c->cd() ;
687     gPad->SetLeftMargin(0.15);
688     frame1->GetYaxis()->SetTitleOffset(1.4);
689     frame1->Draw();
690     if(is_data==data) DrawPrelim();
691     else DrawPrelim(PlottingSetup::luminosity,true);
692 buchmann 1.10 stringstream infotext;
693 buchmann 1.14 infotext << "#splitline{Fit results (" << EdgeFitter::Mode << ">" << jzb_cut << "): }{#splitline{";
694 buchmann 1.15 infotext << "N(Data) = " << EdgeFitter::SFSample->sumEntries() << "}{#splitline{";
695 buchmann 1.10 infotext << "N(Z+Jets) = " << WriteWithError(fzSF.getVal(),fzSF.getError(),3) << "}{#splitline{";
696 buchmann 1.21 //infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
697     infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),0,3) << "}{#splitline{";
698     write_warning(any2string(__LINE__),"Don't have the error yet, need to complete this");
699     if(EdgeFitter::AllowTriangle) {
700     infotext << "N(signal) = " << WriteWithError(fsignalSF.getVal(),fsignalSF.getError(),3) << "}{";
701     infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";
702     } else infotext << "}{}}}}}";
703 buchmann 1.10
704     TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
705     infobox->SetNDC();
706     infobox->SetTextSize(0.03);
707     infobox->Draw();
708 buchmann 1.14 if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_SF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
709     else CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
710 buchmann 1.5 delete c;
711    
712     TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
713     e->cd();
714     gPad->SetLeftMargin(0.15);
715     frame3->GetYaxis()->SetTitleOffset(1.4);
716     frame3->Draw();
717     if(is_data==data) DrawPrelim();
718     else DrawPrelim(PlottingSetup::luminosity,true);
719 buchmann 1.14 if(EdgeFitter::FixedMEdge>=0) CompleteSave(e,"Edge/"+prefix.str()+"_OF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
720     else CompleteSave(e,"Edge/"+prefix.str()+"_OF",false,false);
721 buchmann 1.5 delete e;
722    
723 buchmann 1.6
724    
725    
726 buchmann 1.5 /* TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
727     f->cd();
728     gPad->SetLeftMargin(0.15);
729     frame4->GetYaxis()->SetTitleOffset(1.4);
730     frame4->Draw();
731     if(is_data==data) DrawPrelim();
732     else DrawPrelim(PlottingSetup::luminosity,true);
733 buchmann 1.10 CompleteSave(f,"Edge/"+prefix.str()+"_SF");
734 buchmann 1.5 delete f;*/
735 buchmann 1.7
736 buchmann 1.10
737     /*
738     float maxZ=200;
739     RooWorkspace* wspace = new RooWorkspace();
740     stringstream mllvar;
741     mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
742     wspace->factory(mllvar.str().c_str());
743     wspace->var("mll")->setBins(30);
744     wspace->factory("nSig[1.,0.,100.]");
745     wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
746     wspace->factory("rME[1.12,1.05,1.19]");
747     wspace->factory("effUncert[1.]");
748     EdgeFitter::prepareLimits(wspace, true);
749     */
750    
751     write_warning(__FUNCTION__," A lot missing here to calculate limits");
752    
753 buchmann 1.5 }
754    
755 buchmann 1.2 void EdgeFitter::DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree *signalevents=0) {
756    
757 buchmann 1.5 TCut _cut(cut&&PlottingSetup::basiccut&&PlottingSetup::passtrig);
758    
759     TFile *f = new TFile("workingfile.root","RECREATE");
760    
761     EdgeFitter::InitializeVariables(PlottingSetup::iMllLow,PlottingSetup::iMllHigh,PlottingSetup::jzbHigh,_cut);
762 buchmann 1.2
763     EdgeFitter::PrepareDatasets(is_data);
764 buchmann 1.5
765 buchmann 1.15 EdgeFitter::DrawDatasetContent(is_data);
766    
767 buchmann 1.5 RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
768     RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
769 buchmann 1.11
770    
771 buchmann 1.21
772     EdgeFitter::AllowTriangle=false;
773     EdgeFitter::DoFit(is_data, jzb_cut);
774    
775     write_info(__FUNCTION__,"TAKING SHORTCUT");return;
776    
777     EdgeFitter::AllowTriangle=true;
778    
779 buchmann 1.18 bool ScanMassRange=false;
780 buchmann 1.17 float ScanSteps=5.0;//GeV
781 buchmann 1.12
782 buchmann 1.11
783     if(ScanMassRange) {
784 buchmann 1.12 TFile *fscan = new TFile("fscan.root","UPDATE");
785     TGraph *gr = new TGraph();
786 buchmann 1.21 TGraph *Rgr = new TGraph();
787 buchmann 1.12 stringstream GrName;
788 buchmann 1.14 GrName << "ScanGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
789 buchmann 1.21 stringstream RGrName;
790     RGrName << "ScanRatioGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
791 buchmann 1.12 gr->SetName(GrName.str().c_str());
792 buchmann 1.21 Rgr->SetName(RGrName.str().c_str());
793 buchmann 1.12
794 buchmann 1.11 int i=0;
795 buchmann 1.17 for(float tempMedge=10;tempMedge<=300;tempMedge+=ScanSteps) {
796 buchmann 1.14 write_info(__FUNCTION__,"Now testing Medge="+any2string(tempMedge)+" for "+EdgeFitter::Mode+">"+any2string(jzb_cut));
797 buchmann 1.11 EdgeFitter::FixedMEdge=tempMedge;
798     EdgeFitter::DoFit(is_data, jzb_cut);
799 buchmann 1.21 if(EdgeFitter::FixedMEdgeChi2_H1<0) continue;
800     gr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1);
801     Rgr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1/EdgeFitter::FixedMEdgeChi2_H0);
802 buchmann 1.12 i++;
803 buchmann 1.11 }
804    
805     TCanvas *ScanCan = new TCanvas("ScanCan","ScanCan",500,500);
806     gr->GetXaxis()->SetTitle("m_{edge}");
807     gr->GetXaxis()->CenterTitle();
808 buchmann 1.12 gr->GetYaxis()->SetTitle("#Chi^{2} / NDF");
809 buchmann 1.11 gr->GetYaxis()->CenterTitle();
810 buchmann 1.13 gr->GetYaxis()->SetTitleOffset(0.95);
811     gr->GetXaxis()->SetTitleOffset(0.9);
812     gr->SetLineColor(kBlue);
813     gr->SetTitle("");
814     gr->Draw("AL");
815 buchmann 1.11 stringstream ScanCanSave;
816 buchmann 1.14 ScanCanSave << "Edge/MEdgeScan_"+EdgeFitter::Mode+"_" << jzb_cut;
817 buchmann 1.12 if(is_data) DrawPrelim();
818     else DrawMCPrelim();
819 buchmann 1.13 CompleteSave(ScanCan,ScanCanSave.str());
820 buchmann 1.21
821     Rgr->GetXaxis()->SetTitle("m_{edge}");
822     Rgr->GetXaxis()->CenterTitle();
823     Rgr->GetYaxis()->SetTitle("#Chi^{2} / NDF");
824     Rgr->GetYaxis()->CenterTitle();
825     Rgr->GetYaxis()->SetTitleOffset(0.95);
826     Rgr->GetXaxis()->SetTitleOffset(0.9);
827     Rgr->SetLineColor(kBlue);
828     Rgr->SetTitle("");
829     Rgr->Draw("AL");
830     ScanCanSave.str("");
831     ScanCanSave << "Edge/MEdgeScan_Ratio_"+EdgeFitter::Mode+"_" << jzb_cut;
832     if(is_data) DrawPrelim();
833     else DrawMCPrelim();
834     CompleteSave(ScanCan,ScanCanSave.str());
835 buchmann 1.13 fscan->cd();
836     gr->Write();
837 buchmann 1.11 delete ScanCan;
838 buchmann 1.12 fscan->Close();
839 buchmann 1.11 } else {
840     EdgeFitter::DoFit(is_data, jzb_cut);
841 buchmann 1.21 dout << "Chi^2 (H0) = " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
842     dout << "Chi^2 (H1) = " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
843 buchmann 1.11 }
844    
845    
846 buchmann 1.5 RooMsgService::instance().setGlobalKillBelow(msglevel);
847    
848 buchmann 1.2 f->Close();
849    
850     }
851    
852     void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
853 buchmann 1.14
854     EdgeFitter::Mode="JZB";
855     if(mcjzb=="met[4]") EdgeFitter::Mode="MET";
856    
857 buchmann 1.7 for(int icut=0;icut<(int)jzb_cut.size();icut++) {
858 buchmann 1.2 stringstream addcut;
859     if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
860     if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";
861     TCut jcut(addcut.str().c_str());
862    
863 buchmann 1.5
864 buchmann 1.2 EdgeFitter::DoEdgeFit(mcjzb, datajzb, DataPeakError, MCPeakError, jzb_cut[icut], icut, is_data, jcut&&cut, signalevents);
865    
866     }
867     }