ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C
Revision: 1.10
Committed: Tue Jun 11 19:27:47 2013 UTC (11 years, 10 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.9: +174 -104 lines
Log Message: 
Avoided memory problems; fixed notation (there was ee, mm, and em even though ee contained ee&mm, mm was empty); fixed several minor things to bring plots more in line with our standard ones (still working on that, though); optimized several parameters in the fit and named them correctly

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2    
3 buchmann 1.10 #include <TVirtualIndex.h>
4    
5 buchmann 1.2 #include <RooRealVar.h>
6     #include <RooArgSet.h>
7     #include <RooDataSet.h>
8 buchmann 1.4 #include <RooMCStudy.h>
9     #include <RooCategory.h>
10    
11 buchmann 1.5 #include <RooPlot.h>
12     #include <RooSimultaneous.h>
13     #include <RooAddPdf.h>
14     #include <RooFitResult.h>
15     #include <RooVoigtian.h>
16     #include <RooMsgService.h>
17    
18 buchmann 1.3 #include <RooStats/ModelConfig.h>
19     #include "RooStats/ProfileLikelihoodCalculator.h"
20     #include "RooStats/LikelihoodInterval.h"
21     #include "RooStats/HypoTestResult.h"
22     #include "RooStats/SimpleLikelihoodRatioTestStat.h"
23     #include "RooStats/ProfileLikelihoodTestStat.h"
24     #include "RooStats/HybridCalculatorOriginal.h"
25     #include "RooStats/HypoTestInverterOriginal.h"
26    
27 buchmann 1.5 //#include "ParametrizedEdge.C"
28 buchmann 1.8 #include "EdgeModules/RooSUSYTPdf.cxx"
29     #include "EdgeModules/RooSUSYBkgPdf.cxx"
30 buchmann 1.5
31 buchmann 1.2
32 buchmann 1.1 using namespace std;
33     using namespace PlottingSetup;
34    
35    
36 buchmann 1.2
37    
38 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) {
39     write_error(__FUNCTION__,"Not implemented edge limits yet (returning crap)");
40     ShapeDroplet beta;beta.observed=-12345;beta.SignalIntegral=1;return beta;
41     }
42    
43    
44     void PrepareEdgeShapes(string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
45     write_error(__FUNCTION__,"Not implemented edge shape storage yet");
46     }
47    
48    
49 buchmann 1.2 ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
50    
51    
52     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 buchmann 1.7 void getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
57 buchmann 1.2 void InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut);
58     void PrepareDatasets(int);
59 buchmann 1.5 void DoFit(int is_data, float jzb_cut);
60 buchmann 1.2 string RandomStorageFile();
61     Yield Get_Z_estimate(float,int);
62     Yield Get_T_estimate(float,int);
63 buchmann 1.7 float calcExclusion(RooWorkspace *ws, RooDataSet data, bool calcExclusion);
64     vector<RooDataSet> generateToys(RooWorkspace *ws, int nToys);
65 buchmann 1.4 void prepareLimits(RooWorkspace *ws, bool ComputeBands);
66     TGraph* prepareLM(float mass, float nEv);
67 buchmann 1.2
68     float jzbmax;
69     float mllmin;
70     float mllmax;
71     TCut cut;
72    
73     RooDataSet* AllData;
74 buchmann 1.10 RooDataSet* SFSample;
75     RooDataSet* OFSample;
76 buchmann 1.2
77 buchmann 1.6 bool MarcoDebug=true;
78 buchmann 1.2 }
79    
80 buchmann 1.4 TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
81     float massLM[1];
82     massLM[0]=mass;
83     float accEffLM[1];
84     accEffLM[0]=nEv/PlottingSetup::luminosity;
85     TGraph *lm = new TGraph(1, massLM, accEffLM);
86     lm->GetXaxis()->SetNoExponent(true);
87     lm->GetXaxis()->SetTitle("m_{cut} [GeV]");
88     lm->GetYaxis()->SetTitle("#sigma #times A [pb] 95% CL UL");
89     lm->GetXaxis()->SetLimits(0.,300.);
90     lm->GetYaxis()->SetRangeUser(0.,0.08);
91     lm->SetMarkerStyle(34);
92     lm->SetMarkerColor(kRed);
93     return lm;
94     }
95    
96 buchmann 1.7 vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
97 buchmann 1.10 ws->ls();
98 buchmann 1.4 ws->var("nSig")->setVal(0.);
99     ws->var("nSig")->setConstant(true);
100     RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
101 buchmann 1.7 vector<RooDataSet> theToys;
102 buchmann 1.4
103     RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
104     mcEE.generate(nToys,44,true);
105     RooMCStudy mcMM(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"MM"));
106     mcMM.generate(nToys,44,true);
107     RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
108     mcOSOF.generate(nToys,44,true);
109    
110 buchmann 1.10 RooRealVar mll("m_{ll}","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
111 buchmann 1.4 RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
112     RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
113     RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
114     RooRealVar weight("weight","weight",0,1000,"");
115     RooArgSet observables(mll,jzb,id1,id2,weight);
116    
117     for(int i=0;i<nToys;i++) {
118     RooDataSet* toyEE = (RooDataSet*)mcEE.genData(i);
119     RooDataSet* toyMM = (RooDataSet*)mcMM.genData(i);
120     RooDataSet* toyOSOF = (RooDataSet*)mcOSOF.genData(i);
121     stringstream toyname;
122     toyname << "theToy_" << i;
123     write_warning(__FUNCTION__,"Problem while adding toys");
124 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));
125     theToys.push_back(toyData);
126 buchmann 1.4 }
127     ws->var("nSig")->setVal(17.0);
128     ws->var("nSig")->setConstant(false);
129     return theToys;
130     }
131    
132 buchmann 1.7 void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
133 buchmann 1.4 TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
134     vector<float> theLimits;
135 buchmann 1.9 for(int itoy=0;itoy<(int)theToys.size();itoy++) {
136 buchmann 1.7 float theLimit = calcExclusion(ws,theToys[itoy],false);
137 buchmann 1.4 if(theLimit > 0 ) gauLimit->Fill(theLimit);
138     }
139     const Int_t nQ = 4;
140     Double_t yQ[nQ] = {0.,0.,0.,0.};
141     Double_t xQ[nQ] = {0.022750132,0.1586552539,0.84134474609999998,0.977249868};
142     gauLimit->GetQuantiles(nQ,yQ,xQ);
143     median = gauLimit->GetMean();
144     // median = median1(gauLimit);
145     twoSigmaDown = abs(yQ[0]-median);
146     sigmaDown = abs(yQ[1]-median);
147     sigmaUp = abs(yQ[2]-median);
148     twoSigmaUp = abs(yQ[3]-median);
149     cout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
150     }
151    
152     void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
153     if(ComputeBands) {
154 buchmann 1.7 vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
155 buchmann 1.4 vector<float> medVals;
156     vector<float> medLimits;
157     vector<float> sigmaLimitsDown;
158     vector<float> sigmaLimitsUp;
159     vector<float> twoSigmaLimitsDown;
160     vector<float> twoSigmaLimitsUp;
161     for(int i=20;i<=320;i+=40) {
162     ws->var("nSig")->setVal(10.0);
163     medVals.push_back((float)i);
164     ws->var("m0")->setVal((float)i);
165     ws->var("m0")->setConstant(true);
166     float Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp;
167     EdgeFitter::getMedianLimit(ws,theToys,Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp);
168     medLimits.push_back(Smedian);
169     sigmaLimitsDown.push_back(SsigmaDown);
170     sigmaLimitsUp.push_back(SsigmaUp);
171     twoSigmaLimitsDown.push_back(StwoSigmaDown);
172     twoSigmaLimitsUp.push_back(StwoSigmaUp);
173     }
174     write_warning(__FUNCTION__,"Still need to store limits");
175     } else {
176     vector<float> theVals;
177     vector<float> theLimits;
178     for(int i=20;i<300;i+=5) {
179     ws->var("nSig")->setVal(0.0);
180     theVals.push_back((float)i);
181     ws->var("m0")->setVal((float)i);
182     ws->var("m0")->setConstant(true);
183 buchmann 1.7 // theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
184     write_error(__FUNCTION__,"Error while casting roo data set");
185 buchmann 1.4 }
186    
187 buchmann 1.7 for(int i=0;i<(int)theLimits.size();i++) {
188 buchmann 1.4 if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
189     cout << i << " : " << theVals[i] << endl;
190     theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
191     write_warning(__FUNCTION__,"Need to store limits");
192     }
193     write_warning(__FUNCTION__,"Need to store limits");
194     }
195     }
196     }
197    
198    
199 buchmann 1.7 float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
200     int numberOfToys=50;
201 buchmann 1.3 RooRealVar mu("mu","nSig",0,10000,"");
202     RooArgSet poi = RooArgSet(mu);
203     RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
204     nullParams->setRealValue("nSig",0);
205     RooStats::ModelConfig *model = new RooStats::ModelConfig();
206     model->SetWorkspace(*ws);
207     model->SetPdf("combModel");
208     model->SetParametersOfInterest(poi);
209 buchmann 1.7 // if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
210 buchmann 1.3
211 buchmann 1.7 RooStats::ProfileLikelihoodCalculator plc(data, *model);
212 buchmann 1.3 plc.SetNullParameters(*nullParams);
213     plc.SetTestSize(0.05);
214 buchmann 1.7
215 buchmann 1.3 RooStats::LikelihoodInterval* interval = plc.GetInterval();
216     RooStats::HypoTestResult *htr = plc.GetHypoTest();
217     double theLimit = interval->UpperLimit( mu );
218 buchmann 1.7 // double significance = htr->Significance();
219 buchmann 1.3
220     ws->defineSet("obs","nB");
221     ws->defineSet("poi","nSig");
222    
223     RooStats::ModelConfig b_model = RooStats::ModelConfig("B_model", ws);
224     b_model.SetPdf(*ws->pdf("combModel"));
225     b_model.SetObservables(*ws->set("obs"));
226     b_model.SetParametersOfInterest(*ws->set("poi"));
227     ws->var("nSig")->setVal(0.0); //# important!
228     b_model.SetSnapshot(*ws->set("poi"));
229    
230     RooStats::ModelConfig sb_model = RooStats::ModelConfig("S+B_model", ws);
231     sb_model.SetPdf(*ws->pdf("combModel"));
232     sb_model.SetObservables(*ws->set("obs"));
233     sb_model.SetParametersOfInterest(*ws->set("poi"));
234     ws->var("nSig")->setVal(64.0); //# important!
235     sb_model.SetSnapshot(*ws->set("poi"));
236    
237     RooStats::SimpleLikelihoodRatioTestStat slrts = RooStats::SimpleLikelihoodRatioTestStat(*b_model.GetPdf(),*sb_model.GetPdf());
238     slrts.SetNullParameters(*b_model.GetSnapshot());
239     slrts.SetAltParameters(*sb_model.GetSnapshot());
240     RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
241    
242 buchmann 1.7 RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
243 buchmann 1.3 hc.SetTestStatistic(2);
244 buchmann 1.7 hc.SetNumberOfToys(numberOfToys);
245 buchmann 1.3
246     RooStats::HypoTestInverterOriginal hcInv = RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
247     hcInv.SetTestSize(0.05);
248     hcInv.UseCLs(true);
249     hcInv.RunFixedScan(5,theLimit-0.5,theLimit+0.5);;
250     RooStats::HypoTestInverterResult* hcInt = hcInv.GetInterval();
251     float ulError = hcInt->UpperLimitEstimatedError();
252     theLimit = hcInt->UpperLimit();
253     cout << "Found upper limit : " << theLimit << "+/-" << ulError << endl;
254    
255     return theLimit/PlottingSetup::luminosity;
256    
257     }
258    
259 buchmann 1.2 TTree* SkimTree(int isample) {
260     TTree* newTree = allsamples.collection[isample].events->CloneTree(0);
261     float xsweight=1.0;
262     if(allsamples.collection[isample].is_data==false) xsweight=luminosity*allsamples.collection[isample].weight;
263     if(EdgeFitter::MarcoDebug) {
264     cout << " Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
265     cout << " Going to reduce it with cut " << EdgeFitter::cut << endl;
266     }
267 buchmann 1.6 float edgeWeight;
268     newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
269     float tmll;
270     allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
271 buchmann 1.7 // int id1,id2;
272 buchmann 1.6
273 buchmann 1.2 TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
274 buchmann 1.10 TTreeFormula *Weight = new TTreeFormula("Weight", cutWeight, allsamples.collection[isample].events);
275 buchmann 1.2 float wgt=1.0;
276 buchmann 1.10 // allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
277 buchmann 1.2 for (Int_t entry = 0 ; entry < allsamples.collection[isample].events->GetEntries() ; entry++) {
278     allsamples.collection[isample].events->LoadTree(entry);
279     if (select->EvalInstance()) {
280     allsamples.collection[isample].events->GetEntry(entry);
281 buchmann 1.10 wgt=Weight->EvalInstance();
282 buchmann 1.6 edgeWeight=wgt*xsweight;
283 buchmann 1.2 newTree->Fill();
284     }
285     }
286    
287     if(EdgeFitter::MarcoDebug) cout << " Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
288     return newTree;
289     }
290    
291     void EdgeFitter::InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut) {
292     mllmin=_mllmin;
293     mllmax=_mllmax;
294     jzbmax=_jzbmax;
295     cut=_cut;
296     }
297 buchmann 1.10
298     TTree* MergeTrees(vector<TTree*> trees) {
299     TTree * newtree = (TTree*)trees[0]->CloneTree();
300     trees[0]->GetListOfClones()->Remove(newtree);
301     trees[0]->ResetBranchAddresses();
302     newtree->ResetBranchAddresses();
303    
304     for(int itree=1;itree<trees.size();itree++) {
305     newtree->CopyAddresses(trees[itree]);
306     Long64_t nentries = trees[itree]->GetEntries();
307     for (Long64_t iEntry=0;iEntry<nentries;iEntry++) {
308     trees[itree]->GetEntry(iEntry);
309     newtree->Fill();
310     }
311     trees[itree]->ResetBranchAddresses(); // Disconnect from new tree
312     if (newtree->GetTreeIndex()) {
313     newtree->GetTreeIndex()->Append(trees[itree]->GetTreeIndex(),kTRUE);
314     }
315     if (newtree && newtree->GetTreeIndex()) {
316     newtree->GetTreeIndex()->Append(0,kFALSE); // Force the sorting
317     }
318     }
319     return newtree;
320     }
321    
322    
323    
324 buchmann 1.2 void EdgeFitter::PrepareDatasets(int is_data) {
325 buchmann 1.5 write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
326 buchmann 1.10 // TFile *tempout = new TFile("tempout.root","RECREATE");
327     vector<TTree*> SkimmedTrees;
328     TTree *SkimmedTree[(int)allsamples.collection.size()];
329 buchmann 1.7 for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
330 buchmann 1.2 if(!allsamples.collection[isample].is_active) continue;
331     if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
332     if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
333     if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
334     if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
335     if(EdgeFitter::MarcoDebug) cout << "Considering : " << allsamples.collection[isample].samplename << endl;
336 buchmann 1.10 SkimmedTrees.push_back(SkimTree(isample));
337     // SkimmedTree[isample] = SkimTree(isample);
338     // tempout->cd();
339     // SkimmedTree[isample]->Write();
340     // treelist->Add(SkimmedTree[isample]);
341     //treelist->Add(SkimTree(isample));
342     // allsamples.collection[isample].tfile->Close();
343 buchmann 1.2 }
344 buchmann 1.10
345     TTree *completetree = MergeTrees(SkimmedTrees);
346    
347     // for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
348     // if(SkimmedTree[isample]) SkimmedTree[isample]->Delete();
349     // }
350    
351 buchmann 1.6 if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
352 buchmann 1.2
353 buchmann 1.10 RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
354 buchmann 1.2 RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
355     RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
356 buchmann 1.6 //RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
357     RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
358     RooArgSet observables(mll,id1,id2,edgeWeight);
359 buchmann 1.2
360     string title="CMS Data";
361     if(is_data!=1) title="CMS SIMULATION";
362 buchmann 1.6 RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
363 buchmann 1.2 completetree->Write();
364 buchmann 1.6 delete completetree;
365 buchmann 1.10 // tempout->Close();
366 buchmann 1.2
367 buchmann 1.10 EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
368     EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
369 buchmann 1.2 EdgeFitter::AllData = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
370    
371 buchmann 1.10 SFSample->SetName("SFSample");
372     OFSample->SetName("OFSample");
373 buchmann 1.2 AllData->SetName("AllData");
374    
375     if(EdgeFitter::MarcoDebug) {
376     cout << "Number of events in data sample = " << AllData->numEntries() << endl;
377 buchmann 1.10 cout << "Number of events in eemm sample = " << SFSample->numEntries() << endl;
378     cout << "Number of events in em sample = " << OFSample->numEntries() << endl;
379 buchmann 1.2 }
380 buchmann 1.6
381 buchmann 1.2 }
382    
383 buchmann 1.10 string WriteWithError(float central, float error, int digits) {
384     float ref=central;
385     if(ref<0) ref=-central;
386     int HighestSigDigit = 0;
387     if(ref>1) HighestSigDigit = int(log(ref)/log(10))+1;
388     else HighestSigDigit = int(log(ref)/(log(10)));
389    
390     float divider=pow(10.0,(double(HighestSigDigit-digits)));
391    
392     stringstream result;
393     result << divider*int(central/divider+0.5) << " #pm " << divider*int(error/divider+0.5);
394     return result.str();
395     }
396    
397    
398 buchmann 1.2 string EdgeFitter::RandomStorageFile() {
399     TRandom3 *r = new TRandom3(0);
400     int rho = (int)r->Uniform(1,10000000);
401     stringstream RandomFile;
402     RandomFile << PlottingSetup::cbafbasedir << "/exchange/TempEdgeFile_" << rho << ".root";
403     delete r;
404     return RandomFile.str();
405     }
406    
407     Yield EdgeFitter::Get_Z_estimate(float jzb_cut, int icut) {
408     if(MarcoDebug) write_error(__FUNCTION__,"Returning random Z yield");
409     Yield a(123,45,67); return a;
410     return PlottingSetup::allresults.predictions[icut].Zbkg;
411     }
412    
413     Yield EdgeFitter::Get_T_estimate(float jzb_cut, int icut) {
414     if(MarcoDebug) write_error(__FUNCTION__,"Returning random TTbar yield");
415     Yield a(1234,56,78); return a;
416     return PlottingSetup::allresults.predictions[icut].Flavorsym;
417     }
418    
419 buchmann 1.5 void EdgeFitter::DoFit(int is_data, float jzb_cut) {
420 buchmann 1.10 RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
421 buchmann 1.6 RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
422 buchmann 1.5 RooCategory sample("sample","sample") ;
423 buchmann 1.10 sample.defineType("SF");
424 buchmann 1.5 //sample.defineType("mm");
425 buchmann 1.10 sample.defineType("OF");
426 buchmann 1.6
427 buchmann 1.10 //RooDataSet combData("combData","combined data",mll,Index(sample),Import("SF",*SFSample),Import("mm",*mmSample),Import("OF",*OFSample));
428     RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("SF",*SFSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
429 buchmann 1.6
430 buchmann 1.5
431     //First we make a fit to opposite flavor
432 buchmann 1.10 RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 10000);
433     RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
434     RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
435     RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
436     RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
437     RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
438     RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
439 buchmann 1.5 RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
440 buchmann 1.10 simPdfOF.addPdf(model_OF,"OF");
441 buchmann 1.5 RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save());
442     resultOF->Print();
443    
444 buchmann 1.10 RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarOF");
445 buchmann 1.5 float resultOFpar1 = resultOFpar1_->getVal();
446 buchmann 1.10 //RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarOF");
447 buchmann 1.5 //float resultOFpar2 = resultOFpar2_->getVal();
448     //cout << "caca2.txt" << endl;
449    
450 buchmann 1.10 RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarOF");
451 buchmann 1.5 float resultOFpar3 = resultOFpar3_->getVal();
452    
453 buchmann 1.10 //RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarOF");
454 buchmann 1.5 //float resultOFpar4 = resultOFpar4_->getVal();
455     //cout << "caca4.txt" << endl;
456    
457    
458     // Now same flavor
459 buchmann 1.10 RooRealVar fzSF("fzSF", "fzSF", 5, 0, 100000);
460     RooRealVar meanzSF("meanzSF", "meanzSF", 91.1876, 89, 95);
461     //RooRealVar sigmazSF("sigmazSF", "sigmazSF", 0.5);
462     RooRealVar sigmazSF("sigmazSF", "sigmazSF", 5, 0.5, 5);
463     RooRealVar widthzSF("widthzSF", "widthzSF", 2.94);
464    
465     RooRealVar fttbarSF("fttbarSF", "fttbarSF", 100, 0, 100000);
466     RooRealVar par1ttbarSF("par1ttbarSF", "par1ttbarSF", resultOFpar1, 0, 100);
467     RooRealVar par2ttbarSF("par2ttbarSF", "par2ttbarSF", 1.0);
468     RooRealVar par3ttbarSF("par3ttbarSF", "par3ttbarSF", resultOFpar3, 0, 100);
469     RooRealVar par4ttbarSF("par4ttbarSF", "par4ttbarSF", 2.0);
470 buchmann 1.5
471 buchmann 1.10 RooRealVar fsignalSF("fsignalSF", "fsignalSF", 10, 0, 300);
472     RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
473     RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
474     RooRealVar par3signalSF("par3signalSF", "par3signalSF", 70, 0, 300);
475 buchmann 1.5
476 buchmann 1.10 RooVoigtian zSF("zSF", "zSF", mll, meanzSF, widthzSF, sigmazSF);
477 buchmann 1.5
478    
479 buchmann 1.10 RooSUSYBkgPdf ttbarSF("ttbarSF","ttbarSF", mll , par1ttbarSF, par2ttbarSF, par3ttbarSF, par4ttbarSF);
480     //RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, par2signalSF, par3signalSF);
481     RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, sigmazSF, par3signalSF);
482    
483     /* par1ttbarSF.setConstant(true);
484     par2ttbarSF.setConstant(true);
485     par3ttbarSF.setConstant(true);
486     par4ttbarSF.setConstant(true);*/
487    
488 buchmann 1.5
489 buchmann 1.10 //RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
490     RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
491     RooAddPdf model_em("model_em","model_em", RooArgList(ttbarSF), RooArgList(fttbarSF));
492 buchmann 1.5
493    
494     RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
495 buchmann 1.10 simPdf.addPdf(model_SF,"SF") ;
496     simPdf.addPdf(model_em,"em") ;
497 buchmann 1.5
498     RooFitResult *result = simPdf.fitTo(combData, RooFit::Save());
499     result->Print();
500    
501 buchmann 1.10 RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
502     frame1->GetXaxis()->CenterTitle(1);
503     combData.plotOn(frame1,RooFit::Cut("sample==sample::SF")) ;
504     simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
505     simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Components("ttbarSF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
506     simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Components("zSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
507     simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Components("signalSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
508    
509 buchmann 1.5
510    
511     cout << "Result : " << endl;
512 buchmann 1.10 cout << "f signal : " << fsignalSF.getVal() << " +/- " << fsignalSF.getError() << endl;
513     cout << "f ttbar : " << fttbarSF.getVal() << " +/- " << fttbarSF.getError() << endl;
514     cout << "f tt OF : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
515     cout << "f z SF : " << fzSF.getVal() << " +/- " << fzSF.getError() << endl;
516 buchmann 1.5
517     // The same plot for the cointrol sample slice
518 buchmann 1.10 RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
519     frame3->GetXaxis()->CenterTitle(1);
520     frame3->SetMaximum(frame1->GetMaximum());
521     combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
522     simPdfOF.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
523     simPdfOF.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
524    
525 buchmann 1.5
526     stringstream prefix;
527     if(is_data==data) prefix << "data_";
528     if(is_data==mc) prefix << "mc_";
529     if(is_data==mcwithsignal) prefix << "mcwithS_";
530    
531     prefix << "JZB_" << jzb_cut;
532    
533    
534    
535     TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
536     c->cd() ;
537     gPad->SetLeftMargin(0.15);
538     frame1->GetYaxis()->SetTitleOffset(1.4);
539     frame1->Draw();
540     if(is_data==data) DrawPrelim();
541     else DrawPrelim(PlottingSetup::luminosity,true);
542 buchmann 1.10 stringstream infotext;
543     infotext << "#splitline{Fit results (JZB>" << jzb_cut << "): }{#splitline{";
544     infotext << "N(Data) = " << combData.sumEntries() << "}{#splitline{";
545     infotext << "N(Z+Jets) = " << WriteWithError(fzSF.getVal(),fzSF.getError(),3) << "}{#splitline{";
546     infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
547     infotext << "N(signal) = " << WriteWithError(fsignalSF.getVal(),fsignalSF.getError(),3) << "}{";
548     infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";
549    
550     TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
551     infobox->SetNDC();
552     infobox->SetTextSize(0.03);
553     infobox->Draw();
554     CompleteSave(c,"Edge/"+prefix.str()+"_SF");
555 buchmann 1.5 delete c;
556    
557     TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
558     e->cd();
559     gPad->SetLeftMargin(0.15);
560     frame3->GetYaxis()->SetTitleOffset(1.4);
561     frame3->Draw();
562     if(is_data==data) DrawPrelim();
563     else DrawPrelim(PlottingSetup::luminosity,true);
564 buchmann 1.10 CompleteSave(e,"Edge/"+prefix.str()+"_OF");
565 buchmann 1.5 delete e;
566    
567 buchmann 1.6
568    
569    
570 buchmann 1.5 /* TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
571     f->cd();
572     gPad->SetLeftMargin(0.15);
573     frame4->GetYaxis()->SetTitleOffset(1.4);
574     frame4->Draw();
575     if(is_data==data) DrawPrelim();
576     else DrawPrelim(PlottingSetup::luminosity,true);
577 buchmann 1.10 CompleteSave(f,"Edge/"+prefix.str()+"_SF");
578 buchmann 1.5 delete f;*/
579 buchmann 1.7
580 buchmann 1.10
581     /*
582     float maxZ=200;
583     RooWorkspace* wspace = new RooWorkspace();
584     stringstream mllvar;
585     mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
586     wspace->factory(mllvar.str().c_str());
587     wspace->var("mll")->setBins(30);
588     wspace->factory("nSig[1.,0.,100.]");
589     wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
590     wspace->factory("rME[1.12,1.05,1.19]");
591     wspace->factory("effUncert[1.]");
592     EdgeFitter::prepareLimits(wspace, true);
593     */
594    
595     write_warning(__FUNCTION__," A lot missing here to calculate limits");
596    
597 buchmann 1.5 }
598    
599 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) {
600    
601 buchmann 1.5 TCut _cut(cut&&PlottingSetup::basiccut&&PlottingSetup::passtrig);
602    
603     TFile *f = new TFile("workingfile.root","RECREATE");
604    
605     EdgeFitter::InitializeVariables(PlottingSetup::iMllLow,PlottingSetup::iMllHigh,PlottingSetup::jzbHigh,_cut);
606 buchmann 1.2
607     EdgeFitter::PrepareDatasets(is_data);
608 buchmann 1.5
609     RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
610     RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
611 buchmann 1.6 write_warning(__FUNCTION__,"Deactivated actual fitting procedure ATM");
612 buchmann 1.5 EdgeFitter::DoFit(is_data, jzb_cut);
613     RooMsgService::instance().setGlobalKillBelow(msglevel);
614    
615 buchmann 1.2
616     f->Close();
617    
618     }
619    
620     void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
621 buchmann 1.7 for(int icut=0;icut<(int)jzb_cut.size();icut++) {
622 buchmann 1.2 stringstream addcut;
623     if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
624     if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";
625     TCut jcut(addcut.str().c_str());
626    
627 buchmann 1.5
628 buchmann 1.2 EdgeFitter::DoEdgeFit(mcjzb, datajzb, DataPeakError, MCPeakError, jzb_cut[icut], icut, is_data, jcut&&cut, signalevents);
629    
630     }
631     }