ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
(Generate patch)

Comparing UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C (file contents):
Revision 1.4 by buchmann, Wed Apr 25 12:59:53 2012 UTC vs.
Revision 1.7 by buchmann, Fri Apr 27 14:07:39 2012 UTC

# Line 38 | Line 38 | const char* strip_flip_away(string flipp
38    return flipped_name.c_str();
39   }
40  
41 < void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF, int flipped) {
42 < TH1F *dataob;
43 < if(flipped) dataob = (TH1F*)f->Get("flipped_data_obs");
44 < else dataob = (TH1F*)f->Get("data_obs");
45 < TH1F *signal;
46 <
47 < if(flipped) signal = (TH1F*)f->Get("flipped_signal");
48 < else signal = (TH1F*)f->Get("signal");
49 <
50 < TH1F *background1;
51 < if(flipped) background1 = (TH1F*)f->Get("flipped_TTbarBackground");
52 < else background1 = (TH1F*)f->Get("TTbarBackground");
53 <
54 < TH1F *background2;
55 < if(flipped) background2 = (TH1F*)f->Get("flipped_ZJetsBackground");
56 < else background2 = (TH1F*)f->Get("ZJetsBackground");
41 > void EliminateNegativeEntries(TH1F *histo) {
42 >  for(int i=1;i<=histo->GetNbinsX();i++) {
43 >    if(histo->GetBinContent(i)<0) histo->SetBinContent(i,0);
44 >  }
45 > }
46 >
47 > void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
48 > TH1F *dataob = (TH1F*)f->Get("data_obs");
49 > TH1F *signal = (TH1F*)f->Get("signal");
50 > TH1F *Tbackground = (TH1F*)f->Get("TTbarBackground");
51 > TH1F *Zbackground  = (TH1F*)f->Get("ZJetsBackground");
52  
53   assert(dataob);
54   assert(signal);
55 < assert(background1);
56 < assert(background2);
55 > assert(Tbackground);
56 > assert(Zbackground);
57  
58   ofstream datacard;
59   write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
# Line 76 | Line 71 | void prepare_limit_datacard(string RunDi
71   datacard << "bin             1          1          1\n";
72   datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
73   datacard << "process         0          1          2\n";
74 < datacard << "rate            "<<signal->Integral()<<"         "<<background1->Integral()<<"         "<<background2->Integral()<<"\n";
74 > datacard << "rate            "<<signal->Integral()<<"         "<<Tbackground->Integral()<<"         "<<Zbackground->Integral()<<"\n";
75   datacard << "--------------------------------\n";
76   datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
77   datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
# Line 176 | Line 171 | void generate_shapes_for_systematic(bool
171    }
172      
173    if(signalonly) {
174 <    cout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: " << identifier << ")" << endl;
174 >    dout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: '" << identifier << "')" << endl;
175      TH1F *ZOSSFP = QuickDraw(signalevents,"ZOSSFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
176      TH1F *ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
177      TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
# Line 208 | Line 203 | void generate_shapes_for_systematic(bool
203      Lobs->Add(ZOSSFP);
204      Lpred->Add(ZOSSFN);
205      
206 <    cout << "SITUATION FOR SIGNAL: " << endl;
207 <    cout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
208 <    cout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
206 >    dout << "SITUATION FOR SIGNAL: " << endl;
207 >    dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
208 >    dout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
209      if(PlottingSetup::RestrictToMassPeak) {
210 <      cout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
211 <      cout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
210 >      dout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
211 >      dout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
212      }
213  
214      
# Line 267 | Line 262 | void generate_shapes_for_systematic(bool
262        
263        for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
264          float staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
265 <        cout << "Stat err in bin " << i << " : " << staterr << endl;
266 <        cout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
267 <        cout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
268 <        signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
265 >        dout << "Stat err in bin " << i << " : " << staterr << endl;
266 >        dout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
267 >        dout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
268 >        if(signal->GetBinContent(i)-staterr>0) signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
269 >        else signalStatDn->SetBinContent(i,0);
270          signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
271          signal->SetBinError(i,staterr);
272        }
# Line 291 | Line 287 | void generate_shapes_for_systematic(bool
287        
288        for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
289          float staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
290 <        flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
290 >        if(flippedsignal->GetBinContent(i)-staterr>0) flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
291 >        else flippedsignalStatDn->SetBinContent(i,0);
292          flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
293          flippedsignal->SetBinError(i,staterr);
294        }
# Line 323 | Line 320 | void generate_shapes_for_systematic(bool
320    
321      
322    if(dataonly) {
323 <    cout << "Processing data with datajzb: " << datajzb << endl;
323 >    dout << "Processing data with datajzb: " << datajzb << endl;
324      TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
325      TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
326      TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
# Line 457 | Line 454 | void generate_shapes_for_systematic(bool
454        SQRT(predstaterr);
455        TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
456        bgStatUp->Add(predstaterr);
457 +      EliminateNegativeEntries(bgStatUp);
458        bgStatUp->Write();
459        TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
460        bgStatDn->Add(predstaterr,-1);
461 +      EliminateNegativeEntries(bgStatDn);
462        bgStatDn->Write();
463   //      delete bgStatDn;
464   //      delete bgStatUp;
# Line 476 | Line 475 | void generate_shapes_for_systematic(bool
475        SQRT(flippedpredstaterr);
476        TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
477        fbgStatUp->Add(predstaterr);
478 +      EliminateNegativeEntries(fbgStatUp);
479        fbgStatUp->Write();
480        TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
481        fbgStatDn->Add(predstaterr,-1);
482 +      EliminateNegativeEntries(fbgStatDn);
483        fbgStatDn->Write();
484   //      delete fbgStatDn;
485   //      delete fbgStatUp;
# Line 491 | Line 492 | void generate_shapes_for_systematic(bool
492        SQRT(Tpredstaterr);
493        TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
494        TpredStatUp->Add(Tpredstaterr);
495 +      EliminateNegativeEntries(TpredStatUp);
496        TpredStatUp->Write();
497        TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
498        TpredStatDn->Add(Tpredstaterr,-1);
499 +      EliminateNegativeEntries(TpredStatDn);
500        TpredStatDn->Write();
501   //      delete TpredStatDn;
502   //      delete TpredStatUp;
# Line 506 | Line 509 | void generate_shapes_for_systematic(bool
509        SQRT(fTpredstaterr);
510        TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
511        fTpredStatUp->Add(fTpredstaterr);
512 +      EliminateNegativeEntries(fTpredStatUp);
513        fTpredStatUp->Write();
514        TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
515        fTpredStatDn->Add(fTpredstaterr,-1);
516 +      EliminateNegativeEntries(fTpredStatDn);
517        fTpredStatDn->Write();
518   //      delete fTpredStatDn;
519   //      delete fTpredStatUp;
# Line 522 | Line 527 | void generate_shapes_for_systematic(bool
527        SQRT(Zpredstaterr);
528        TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
529        ZpredStatUp->Add(Zpredstaterr);
530 +      EliminateNegativeEntries(ZpredStatUp);
531        ZpredStatUp->Write();
532        TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
533        ZpredStatDn->Add(Zpredstaterr,-1);
534 +      EliminateNegativeEntries(ZpredStatDn);
535        ZpredStatDn->Write();
536   //      delete ZpredStatDn;
537   //      delete ZpredStatUp;
# Line 537 | Line 544 | void generate_shapes_for_systematic(bool
544        SQRT(fTpredstaterr);
545        TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
546        fZpredStatUp->Add(fZpredstaterr);
547 +      EliminateNegativeEntries(fZpredStatUp);
548        fZpredStatUp->Write();
549        TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
550        fZpredStatDn->Add(fZpredstaterr,-1);
551 +      EliminateNegativeEntries(fZpredStatDn);
552        fZpredStatDn->Write();
553   //      delete fZpredStatDn;
554   //      delete fZpredStatUp;
# Line 559 | Line 568 | void generate_shapes_for_systematic(bool
568        SQRT(predsyserr);
569        TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
570        bgSysUp->Add(predsyserr);
571 +      EliminateNegativeEntries(bgSysUp);
572        bgSysUp->Write();
573        TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
574        bgSysDn->Add(predsyserr,-1);
575 +      EliminateNegativeEntries(bgSysDn);
576        bgSysDn->Write();
577        delete predsyserr;
578        
# Line 577 | Line 588 | void generate_shapes_for_systematic(bool
588        SQRT(fpredsyserr);
589        TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
590        fbgSysUp->Add(fpredsyserr);
591 +      EliminateNegativeEntries(fbgSysUp);
592        fbgSysUp->Write();
593        TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
594        fbgSysDn->Add(fpredsyserr,-1);
595 +      EliminateNegativeEntries(fbgSysDn);
596        fbgSysDn->Write();
597        delete fpredsyserr;
598        
# Line 592 | Line 605 | void generate_shapes_for_systematic(bool
605        SQRT(Tpredsyserr);
606        TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
607        TpredSysUp->Add(Tpredsyserr);
608 +      EliminateNegativeEntries(TpredSysUp);
609        TpredSysUp->Write();
610        TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
611        TpredSysDn->Add(Tpredsyserr,-1);
612 +      EliminateNegativeEntries(TpredSysDn);
613        TpredSysDn->Write();
614        delete Tpredsyserr;
615        
# Line 606 | Line 621 | void generate_shapes_for_systematic(bool
621        SQRT(fTpredsyserr);
622        TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
623        fTpredSysUp->Add(fTpredsyserr);
624 +      EliminateNegativeEntries(fTpredSysUp);
625        fTpredSysUp->Write();
626        TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
627        fTpredSysDn->Add(fTpredsyserr,-1);
628 +      EliminateNegativeEntries(fTpredSysDn);
629        fTpredSysDn->Write();
630        delete fTpredsyserr;
631        
# Line 623 | Line 640 | void generate_shapes_for_systematic(bool
640        SQRT(Zpredsyserr);
641        TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
642        ZpredSysUp->Add(Zpredsyserr);
643 +      EliminateNegativeEntries(ZpredSysUp);
644        ZpredSysUp->Write();
645        TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
646        ZpredSysDn->Add(Zpredsyserr,-1);
647 +      EliminateNegativeEntries(ZpredSysDn);
648        ZpredSysDn->Write();
649        delete Zpredsyserr;
650        
# Line 638 | Line 657 | void generate_shapes_for_systematic(bool
657        SQRT(fZpredsyserr);
658        TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
659        fZpredSysUp->Add(fZpredsyserr);
660 +      EliminateNegativeEntries(fZpredSysUp);
661        fZpredSysUp->Write();
662        TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
663        fZpredSysDn->Add(fZpredsyserr,-1);
664 +      EliminateNegativeEntries(fZpredSysDn);
665        fZpredSysDn->Write();
666        delete fZpredsyserr;
667      }
668      
669   /*if(identifier=="") {
670    for(int i=0;i<binning.size()-1;i++) {
671 <    cout << "[ " << binning[i] << " , " << binning[i+1] << "] : O " << obs->GetBinContent(i+1) << " P " << pred->GetBinContent(i+1) << " (Z: " << Zpred->GetBinContent(i+1) << " , T: " << Tpred->GetBinContent(i+1) << ")" << endl;
671 >    dout << "[ " << binning[i] << " , " << binning[i+1] << "] : O " << obs->GetBinContent(i+1) << " P " << pred->GetBinContent(i+1) << " (Z: " << Zpred->GetBinContent(i+1) << " , T: " << Tpred->GetBinContent(i+1) << ")" << endl;
672    }
673   }*/
674      delete ZOSSFP;
# Line 679 | Line 700 | ShapeDroplet LimitsFromShapes(bool asymp
700    
701    ensure_directory_exists(RunDirectory);
702    
703 <  TFile *datafile = new TFile("../StoredShapes.root","READ");
703 >  TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str(),"READ");
704    if(datafile->IsZombie()) {
705      write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
706      assert(!datafile->IsZombie());
707    }
708 <  cout << "Run Directory: " << RunDirectory << endl;
708 >  dout << "Run Directory: " << RunDirectory << endl;
709    TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
710    
711    TIter nextkey(datafile->GetListOfKeys());
# Line 717 | Line 738 | ShapeDroplet LimitsFromShapes(bool asymp
738    bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
739    
740    MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
741 <  if(mceff<0) flipped=1;
741 >  if(mceff<0) {
742 >    flipped=1;
743 >    write_info(__FUNCTION__,"Doing flipping!");
744 >  }
745    doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
746    float PDFuncert=0;
747    int NPdfs = get_npdfs(events);
# Line 763 | Line 787 | ShapeDroplet LimitsFromShapes(bool asymp
787        obj->Write();
788    }
789  
790 <  prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert,flipped);
790 >  prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
791    
792    final_limfile->Close();
793    limfile->Close();
# Line 776 | Line 800 | ShapeDroplet LimitsFromShapes(bool asymp
800    else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.5 * firstGuess) << " " << int(2*firstGuess); // ASYMPTOTIC LIMITS
801    dout <<"Going to run : " << command.str() << endl;
802    int CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
803 <  cout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
804 <  assert(CreatedModelFileExitCode==0);
803 >  dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
804 >  if(!(CreatedModelFileExitCode==0)) {
805 >    write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
806 >    ShapeDroplet alpha;
807 >    alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
808 >    alpha.SignalIntegral=1;
809 >    return alpha;
810 >  }
811    ShapeDroplet alpha;
812    alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
813    alpha.PDF=PDFuncert;
# Line 787 | Line 817 | ShapeDroplet LimitsFromShapes(bool asymp
817    dout << alpha << endl;
818  
819    dout << "Everything is saved in " << RunDirectory << endl;
820 <  dout << "Will transfer model and datacard over for possible post-processing" << endl;
821 <  dout << "   1) Make sure models directory exists ... " << std::flush;
820 >  dout << "Cleaning up ... " << std::flush;
821 > /*  dout << "   1) Make sure models directory exists ... " << std::flush;
822    gSystem->Exec("mkdir -p models/");
823    dout << " ok!" << endl;
824    dout << "   2) Deleting any previous model files with the same name ... " << std::flush;
# Line 810 | Line 840 | ShapeDroplet LimitsFromShapes(bool asymp
840    copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root";
841    gSystem->Exec(copycommand3.str().c_str());
842    dout << " ok!" << endl;
843 <  dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;
844 < write_warning(__FUNCTION__,"Watch out : need to uncomment the line below to remove the original working directory again");
815 < //  gSystem->Exec(("rm -r "+RunDirectory).c_str());
843 >  dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;*/
844 >  gSystem->Exec(("rm -r "+RunDirectory).c_str());
845    dout << " ok!" << endl;
846    delete limcan;
847    return alpha;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines