ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C
Revision: 1.8
Committed: Mon Jul 16 02:31:01 2012 UTC (12 years, 9 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.7: +2 -2 lines
Log Message: 
Re-established possibility to compile CBAF on private computers by uploading SUSY PDF functions

File Contents

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