ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/LJMet/MultivariateAnalysis/root/toyMC.C
(Generate patch)

Comparing UserCode/LJMet/MultivariateAnalysis/root/toyMC.C (file contents):
Revision 1.8 by kukartse, Wed May 20 01:20:04 2009 UTC vs.
Revision 1.9 by kukartse, Fri May 22 16:39:00 2009 UTC

# Line 1 | Line 1
1   #include <iostream>
2 < #include <pair.h>
2 > #include <utility>
3   #include "TH1F.h"
4   #include "TObjArray.h"
5   #include "TFractionFitter.h"
# Line 26 | Line 26 | public:
26    
27    void load_templates();
28    void load_ejets_templates();
29 <  int do_fit(TH1F * _data);
29 >
30 >  //
31 >  //_____ Fix bkg2 fraction to qcd_frac _________________________________
32 >  //      Do not fix if qcd_frac is negative
33 >  //
34 >  int do_fit(TH1F * _data, double qcd_frac=-1.0);
35    int generate_toy(int n_sig, int n_phys, int n_qcd);
36    void set_overflow_bins(TH1F * h);
37    void debug_hist(TH1F * h);
# Line 40 | Line 45 | public:
45    TH1F * wzjets;
46    TH1F * qcd;
47  
48 +  // templates for systematics
49 +  TH1F * ttbar_s;
50 +  TH1F * wzjets_s;
51 +  TH1F * qcd_s;
52 +
53    int n_bins;
54  
55    TObjArray * mc;
# Line 222 | Line 232 | void toyMC::load_ejets_templates(){
232   }
233  
234  
235 < int toyMC::do_fit(TH1F * _data){
235 > int toyMC::do_fit(TH1F * _data, double qcd_frac){
236    delete fit;
237    fit = new TFractionFitter(_data, mc); // initialise
238  
239    fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
240    fit->Constrain(2,0.0,1.0);
241 <  //fit->Constrain(3,0.0,1.0);
242 <
241 >  if (qcd_frac<0.0){
242 >    fit->Constrain(3,0.0,1.0);
243 >  }
244 >  else{
245 >    fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
246 >  }
247    //fit->Constrain(2,0.1651865,0.1651866);
248    //fit->Constrain(3,0.024999,0.025001);
249 <  fit->Constrain(3,0.176895306858,0.176895306860);
249 >  //fit->Constrain(3,0.176895306858,0.176895306860);
250    //fit->Constrain(3,0.37259,0.37261);
251    //fit->Constrain(3,0.036759,0.036761); // mu+jets true
252    //fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
253    //fit->Constrain(3,0.4422735,0.4422736);
254    //
241  //double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
242  //fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
243  //
255    fit->SetRangeX(1,n_bins);                    // use only some bins to fit
256    Int_t status = fit->Fit();               // perform the fit
257    if (status!=0) status = fit->Fit();
# Line 264 | Line 275 | int toyMC::fit1(TH1F * _data, double bkg
275   }
276   */
277  
278 + //
279 + //_____ generates a pseudoexperiment, given n_events, smears by Poisson _
280 + //
281   int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
282    int n_events=0;
283    delete toy_data;
# Line 286 | Line 300 | int toyMC::generate_toy(int n_sig, int n
300    for (int i=0;i<_nbg;i++){
301      toy_data->Fill(wzjets->GetRandom());
302    }
303 +
304    for (int i=0;i<_nqcd;i++){
305      toy_data->Fill(qcd->GetRandom());
306    }
# Line 365 | Line 380 | std::pair<double,double> toyMC::do_toys(
380    toy_error->SetFillColor(44);
381    toy_pull->SetFillColor(44);
382    for (int i=0;i<n_exp;i++){
383 <    generate_toy(n_sig,n_bkg1,n_bkg2);
383 >    double n_bkg2_smeared = n_bkg2;
384 >    /************
385 >    //_____ smear qcd rate due to xsec uncertainty
386 >    //
387 >    n_bkg2_smeared = -1.0;
388 >    while (n_bkg2_smeared < 0.0){
389 >      n_bkg2_smeared = _random.Gaus(n_bkg2,n_bkg2);
390 >    }
391 >    *************/
392 >    generate_toy(n_sig,n_bkg1,n_bkg2_smeared);
393      int _ttbar = n_gen_sig;
394 <    int status = do_fit(toy_data);
394 >    //double qcd_rate = 9.5;
395 >    //double qcd_rate = 16.6;
396 >    double qcd_rate = 2.4;
397 >    double qcd_rate_err = 7.1;
398 >    double qcd_rate_smeared = qcd_rate;
399 >    /*************
400 >    //_____ smearing of the ABCD prediction
401 >    //
402 >    double qcd_rate_smeared = -1.0;
403 >    while (qcd_rate_smeared < 0.0){
404 >      qcd_rate_smeared = _random.Gaus(qcd_rate,qcd_rate_err);
405 >    }    
406 >    **************/
407 >    double qcd_frac = qcd_rate_smeared/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
408 >    int status = do_fit(toy_data, qcd_frac);
409      if (status==0){
410        Double_t value,error;
411        fit->GetResult(0,value,error);
# Line 383 | Line 421 | std::pair<double,double> toyMC::do_toys(
421    toy_mean->Fit("gaus");
422    TF1 * _fit = toy_mean->GetFunction("gaus");
423    result.first = _fit->GetParameter(1);
424 <  result.second = _fit->GetParameter(2);
424 >  result.second = _fit->GetParError(1);
425 >  //result.second = _fit->GetParameter(2);
426    pad->SaveAs("mean.eps");
427    pad = c->cd(2);
428    toy_error->Draw();

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines