ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C
Revision: 1.11
Committed: Wed Jun 12 05:43:20 2013 UTC (11 years, 10 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.10: +52 -10 lines
Log Message: 
Added possibility to scan whole mass range and extract Chi2 for each mass point

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