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.2 by buchmann, Thu Apr 12 14:58:25 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 <
59 > write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
60   datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
61   datacard << "Writing this to a file.\n";
62   datacard << "imax 1\n"; // number of channels
# 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 << "       1.0       1.0    luminosity uncertainty\n"; // only affects MC -> signal!
77 < // datacard << "bgnorm   lnN    1.00       1.4  uncertainty on our prediction (40%)\n";
78 < datacard << "Stat   shape    1          1          1    statistical uncertainty\n";
79 < datacard << "Sys    shape    -          1          1    systematic uncertainty\n";
76 > datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
77 > datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
78 > datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
79 > datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
80 > datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
81 > datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
82   datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
83   datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
84   if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
85  
86 <
90 < // datacard << "peak   shape    1          1    uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, \n";
91 < // datacard << "#                                so divide the unit gaussian by 2 before doing the interpolation\n";
86 > // datacard << "peak   shape    1          1    uncertainty on signal resolution.
87   datacard.close();
88   }
89  
# Line 109 | Line 104 | TH1F* QuickDraw(TTree *events, string hn
104        stringstream mSUGRAweight;
105        float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
106        totalxs+=xs;
107 +      mSUGRAweight << "(0.5*(2*(process==" << i << "))*" << xs << ")";
108        events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
109      }
110 <    histo->Scale(1.0/totalxs);
110 >    histo->Scale(1.0/totalxs); // to get back to 1 /pb normalization
111    }
112    events->Draw("eventNum",addcut.c_str(),"goff");
113    float nevents = events->GetSelectedRows();
# Line 147 | Line 143 | TH1F* Multiply(TH1F *h1, TH1F *h2) {
143   }
144  
145   void generate_shapes_for_systematic(bool signalonly, bool dataonly,TFile *limfile, TTree *signalevents, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan, string addcut, map <  pair<float, float>, map<string, float>  > &xsec) {
146 +  
147    binning.push_back(8000);
148    limfile->cd();
149 <  dout << "Creatig shape templates ";
149 >  dout << "Creating shape templates ";
150    if(identifier!="") dout << "for "<<identifier;
151    dout << " ... " << endl;
152    
153    TCut limitnJetcut;
154 +  
155    if(JES==noJES) limitnJetcut=cutnJets;
156    else {
157      if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
158      if(JES==JESup) limitnJetcut=cutnJetsJESup;
159    }
160    
163 –  stringstream mSUGRAweight;
164 –  if(SUSYScanSpace::SUSYscantype==mSUGRA) {
165 –    //for(int i<
166 –    mSUGRAweight << "bla";
167 –  } else mSUGRAweight << "(1)";
168 –  
169 –  write_warning(__FUNCTION__,"Not done yet for mSUGRA");
170 –  
161    float zjetsestimateuncert=zjetsestimateuncertONPEAK;
162    float emuncert=emuncertONPEAK;
163    float emsidebanduncert=emsidebanduncertONPEAK;
# Line 181 | 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 213 | 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 <    
213 >
214      
215      flippedLobs->Add(ZOSSFN);
216      flippedLpred->Add(ZOSSFP);
# Line 272 | 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 <        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 293 | 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 325 | 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 459 | 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 478 | 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 493 | 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 508 | 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 524 | 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 539 | 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 561 | 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 579 | 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 594 | 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 608 | 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 625 | 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 640 | 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 667 | Line 686 | void generate_shapes_for_systematic(bool
686    
687   }
688  
689 < ShapeDroplet LimitsFromShapes(TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
689 > ShapeDroplet LimitsFromShapes(bool asymptotic, float firstGuess, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
690    map <  pair<float, float>, map<string, float>  >  xsec;
691    if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
692    
674 –  write_info(__FUNCTION__,"Implementing the shape function");
693    time_t timestamp;
694    tm *now;
695    timestamp = time(0);
# Line 682 | Line 700 | ShapeDroplet LimitsFromShapes(TTree *eve
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 720 | Line 738 | ShapeDroplet LimitsFromShapes(TTree *eve
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;
746 >  float PDFuncert=0;
747    int NPdfs = get_npdfs(events);
748    if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
749  
729 –
750    float JZBscale=scaledown;
751    if(scaleup>scaledown) JZBscale=scaleup;
752    dout << endl;
# Line 735 | Line 755 | ShapeDroplet LimitsFromShapes(TTree *eve
755    dout << "       JZB Scale (JSU) : " << JZBscale << endl;
756    dout << "       PDF             : " << PDFuncert << endl;
757    
758 <  prepare_limit_datacard(RunDirectory,limfile,JZBscale,PDFuncert,flipped);
758 >  float SignalIntegral;
759 >  if(flipped) {
760 >    TH1F *flisi = (TH1F*)(limfile->Get("flipped_signal"));
761 >    SignalIntegral= flisi ->Integral();
762 >    delete flisi;
763 >  }
764 >  else {
765 >    TH1F *flisi = (TH1F*)(limfile->Get("signal"));
766 >    SignalIntegral= flisi ->Integral();
767 >    delete flisi;
768 >  }
769 >    
770    TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
771    
772    TIter nnextkey(limfile->GetListOfKeys());
773    TKey *nkey;
774 +  
775 +  if(SignalIntegral==0) SignalIntegral=1; //case when the integral is zero - don't try to normalize, just leave it as it is
776 +  
777    while ((nkey = (TKey*)nnextkey()))
778    {
779        TH1F *obj =(TH1F*) nkey->ReadObj();
# Line 747 | Line 781 | ShapeDroplet LimitsFromShapes(TTree *eve
781        if(!flipped&&Contains(obj->GetName(),"flipped")) continue;
782        if(flipped) obj->SetName(strip_flip_away(obj->GetName()));
783        final_limfile->cd();
784 +      if(Contains(obj->GetName(),"signal")) {
785 +        obj->Scale(1.0/SignalIntegral);
786 +      }
787        obj->Write();
788    }
789  
790 +  prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
791 +  
792    final_limfile->Close();
793    limfile->Close();
794 <  write_info(__FUNCTION__,"Shape root file and datacard have been generated in "+RunDirectory);
795 <  
796 <  int CreatedModelFileExitCode = gSystem->Exec(("bash CreateModel.sh "+RunDirectory+" susydatacard.txt").c_str());
797 <  cout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
798 <  assert(CreatedModelFileExitCode==0);
794 >  dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
795 >  stringstream command;
796 >  if(asymptotic) {
797 >    if(firstGuess>0) command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
798 >    else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
799 >  }
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 >  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;
814    alpha.JSU=JZBscale;
815 +  alpha.SignalIntegral=SignalIntegral;
816    dout << "Done reading limit information from droplet" << endl;
817    dout << alpha << endl;
818  
767 –  write_warning(__FUNCTION__,"Once you're done implementing the limit computation, please uncomment the following line to enable clean up.");
819    dout << "Everything is saved in " << RunDirectory << endl;
820 < //  gSystem->Exec(("rm -r "+RunDirectory).c_str());
821 <
822 <    return alpha;
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;
825 >  stringstream delcommand;
826 >  delcommand << "rm models/model_" << name << ".root";
827 >  gSystem->Exec(delcommand.str().c_str());
828 >  stringstream delcommand2;
829 >  delcommand2 << "rm models/model_" << name << ".txt";
830 >  gSystem->Exec(delcommand2.str().c_str());
831 >  dout << " ok!" << endl;
832 >  dout << "   3) Transfering the new model files ... " << std::flush;
833 >  stringstream copycommand;
834 >  copycommand << "cp " << RunDirectory << "/model.root models/model_" << name << ".root";
835 >  gSystem->Exec(copycommand.str().c_str());
836 >  stringstream copycommand2;
837 >  copycommand2 << "cp " << RunDirectory << "/susydatacard.txt models/model_" << name << ".txt";
838 >  gSystem->Exec(copycommand2.str().c_str());
839 >  stringstream copycommand3;
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 >  gSystem->Exec(("rm -r "+RunDirectory).c_str());
845 >  dout << " ok!" << endl;
846 >  delete limcan;
847 >  return alpha;
848   }
849  
850   void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {

Diff Legend

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