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.13 by buchmann, Thu May 17 19:52:55 2012 UTC vs.
Revision 1.15 by buchmann, Tue Jul 3 06:33:52 2012 UTC

# Line 55 | Line 55 | vector<float> get_expected_points(float
55      if(gaus->Integral(0,currentpoint)>((points.size()+1.0)/(npoints+2.0))) {
56        if(Verbosity>0) dout << "Added points for calculation at " << currentpoint << " (integral is " << gaus->Integral(0,currentpoint) << ")" << endl;
57        points.push_back(currentpoint);
58 <      if(points.size()==npoints) break;
58 >      if((int)points.size()==npoints) break;
59      }
60      currentpoint+=stepsize;
61    }
# Line 296 | Line 296 | TH1F* QuickDraw(TTree *events, string hn
296    return histo;
297   }
298  
299 + TH2F* empty2d(string hname) {
300 +    TH2F *h = new TH2F(hname.c_str(),hname.c_str(),1000,0,1000,1000,0,1000);
301 +    return h;
302 + }
303  
304 < void SQRT(TH1F *h) {
305 <  for (int i=1;i<=h->GetNbinsX();i++) {
306 <    h->SetBinContent(i,TMath::Sqrt(h->GetBinContent(i)));
307 <  }
304 > TH2F* QuickDraw2d(TTree *events, string hname, string variable, string xlabel, string ylabel, TCut cut, string addcut, bool dataormc, float luminosity, map <  pair<float, float>, map<string, float>  > &xsec) {
305 >    TH2F *histo = empty2d(hname);
306 >    histo->Sumw2();
307 >    stringstream drawcommand;
308 >    drawcommand << variable << ":mll>>" << hname;
309 >    if(SUSYScanSpace::SUSYscantype!=mSUGRA) {
310 >        events->Draw(drawcommand.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight);
311 >    } else {
312 >        float totalxs=0;
313 >        for(int i=0;i<12;i++) {
314 >            stringstream drawcommand2;
315 >            drawcommand2 << variable << ">>+" << hname;
316 >            stringstream mSUGRAweight;
317 >            float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
318 >            totalxs+=xs;
319 >            mSUGRAweight << "(0.5*(2*(process==" << i << "))*" << xs << ")";
320 >            events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
321 >        }
322 >        histo->Scale(1.0/totalxs); // to get back to 1 /pb normalization
323 >    }
324 >    events->Draw("eventNum",addcut.c_str(),"goff");
325 >    float nevents = events->GetSelectedRows();
326 >    QuickDrawNevents=nevents;
327 >    histo->Scale(luminosity/nevents); // this will give a histogram which is normalized to 1 pb so that any limit will be with respect to 1 pb
328 >    
329 >    return histo;
330   }
331  
332   TH1F* Multiply(TH1F *h1, TH1F *h2) {
# Line 326 | Line 352 | void generate_mc_shapes(bool dotree, TTr
352    mbinning.push_back(-8000);
353    for(int i=binning.size()-1;i>=0;i--) mbinning.push_back(-binning[i]);
354    mbinning.push_back(0);
355 <  for(int i=0;i<binning.size();i++) mbinning.push_back(binning[i]);
355 >  for(int i=0;i<(int)binning.size();i++) mbinning.push_back(binning[i]);
356    mbinning.push_back(8000);
357  
358    
# Line 417 | Line 443 | void generate_shapes_for_systematic(bool
443      TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
444      TH1F *ZOSOFP;
445      TH1F *ZOSOFN;
446 <
446 >    
447 >    TH2F *ZOSSFP2d = QuickDraw2d(signalevents,"ZOSSFP2d",mcjzb,    "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
448 >    TH2F *ZOSSFN2d = QuickDraw2d(signalevents,"ZOSSFN2d","-"+mcjzb,"JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
449 >    TH2F *ZOSOFP2d;
450 >    TH2F *ZOSOFN2d;
451 >    
452      if(!PlottingSetup::FullMCAnalysis) {
453        ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
454        ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
455 +      
456 +      ZOSOFP2d = QuickDraw2d(signalevents,"ZOSOFP2d",mcjzb,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
457 +      ZOSOFN2d = QuickDraw2d(signalevents,"ZOSOFN2d","-"+mcjzb, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
458      }
459      
460      TH1F *SBOSSFP;
# Line 428 | Line 462 | void generate_shapes_for_systematic(bool
462      TH1F *SBOSSFN;
463      TH1F *SBOSOFN;
464      
465 +    TH2F *SBOSSFP2d;
466 +    TH2F *SBOSOFP2d;
467 +    TH2F *SBOSSFN2d;
468 +    TH2F *SBOSOFN2d;
469 +      
470      if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
471        SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
472        SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
473        SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
474        SBOSOFN = QuickDraw(signalevents,"SBOSOFN","-"+mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
475 +      
476 +      SBOSSFP2d = QuickDraw2d(signalevents,"SBOSSFP2d",mcjzb,"JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
477 +      SBOSOFP2d = QuickDraw2d(signalevents,"SBOSOFP2d",mcjzb,"JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
478 +      SBOSSFN2d = QuickDraw2d(signalevents,"SBOSSFN2d","-"+mcjzb,"JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
479 +      SBOSOFN2d = QuickDraw2d(signalevents,"SBOSOFN2d","-"+mcjzb,"JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
480      }
481      
482      string signalname="signal";
483 <    if(identifier!="") {
440 <      signalname="signal_"+identifier;
441 <    }
483 >    if(identifier!="") signalname="signal_"+identifier;
484      
485      TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
486      TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
487      
488 +    TH2F *Lobs2d = empty2d("Lobs2d");
489 +    TH2F *Lpred2d = empty2d("Lobs2d");
490 +    
491      TH1F *flippedLobs = new TH1F("flippedLobs","flippedLobs",binning.size()-1,&binning[0]);
492      TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
493      
494 +    TH2F *flippedLobs2d = empty2d("flippedLobs2d");
495 +    TH2F *flippedLpred2d = empty2d("flippedLpred2d");
496 +    
497      Lobs->Add(ZOSSFP);
498 <    if(!PlottingSetup::FullMCAnalysis) Lpred->Add(ZOSSFN);
498 >    Lobs2d->Add(ZOSSFP2d);
499 >    if(!PlottingSetup::FullMCAnalysis) {
500 >      Lpred->Add(ZOSSFN);
501 >      Lpred2d->Add(ZOSSFN2d);
502 >    }
503      
504      dout << "SITUATION FOR SIGNAL: " << endl;
453 –    
454 –    
505      if(PlottingSetup::FullMCAnalysis) {
506 <      dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << endl;
506 >        dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << endl;
507      } else {
508 <      dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
509 <      dout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
510 <      if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
511 <        dout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
512 <        dout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
513 <      }
508 >        dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
509 >        dout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
510 >      
511 >        if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
512 >            dout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
513 >            dout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
514 >        }
515      }
465 –
516      
517      flippedLobs->Add(ZOSSFN);
518 <    if(!PlottingSetup::FullMCAnalysis) flippedLpred->Add(ZOSSFP);
518 >    flippedLobs2d->Add(ZOSSFN2d);
519 >    
520 >    if(!PlottingSetup::FullMCAnalysis) {
521 >      flippedLpred->Add(ZOSSFP);
522 >      flippedLpred2d->Add(ZOSSFP2d);
523 >    }
524      
525      if(!PlottingSetup::FullMCAnalysis) {
526        if(PlottingSetup::RestrictToMassPeak) {
527 <        Lpred->Add(ZOSOFP,1.0/3);
528 <        Lpred->Add(ZOSOFN,-1.0/3);
529 <        Lpred->Add(SBOSSFP,1.0/3);
530 <        Lpred->Add(SBOSSFN,-1.0/3);
531 <        Lpred->Add(SBOSOFP,1.0/3);
532 <        Lpred->Add(SBOSOFN,-1.0/3);
533 <        
534 <        //flipped prediction
535 <        flippedLpred->Add(ZOSOFP,-1.0/3);
536 <        flippedLpred->Add(ZOSOFN,1.0/3);
537 <        flippedLpred->Add(SBOSSFP,-1.0/3);
538 <        flippedLpred->Add(SBOSSFN,1.0/3);
539 <        flippedLpred->Add(SBOSOFP,-1.0/3);
540 <        flippedLpred->Add(SBOSOFN,1.0/3);
527 >          Lpred->Add(ZOSOFP,1.0/3);
528 >          Lpred->Add(ZOSOFN,-1.0/3);
529 >          Lpred->Add(SBOSSFP,1.0/3);
530 >          Lpred->Add(SBOSSFN,-1.0/3);
531 >          Lpred->Add(SBOSOFP,1.0/3);
532 >          Lpred->Add(SBOSOFN,-1.0/3);
533 >          
534 >          //flipped prediction
535 >          flippedLpred->Add(ZOSOFP,-1.0/3);
536 >          flippedLpred->Add(ZOSOFN,1.0/3);
537 >          flippedLpred->Add(SBOSSFP,-1.0/3);
538 >          flippedLpred->Add(SBOSSFN,1.0/3);
539 >          flippedLpred->Add(SBOSOFP,-1.0/3);
540 >          flippedLpred->Add(SBOSOFN,1.0/3);
541 >          
542 >          //2d stuff: regular
543 >          Lpred2d->Add(ZOSOFP2d,1.0/3);
544 >          Lpred2d->Add(ZOSOFN2d,-1.0/3);
545 >          Lpred2d->Add(SBOSSFP2d,1.0/3);
546 >          Lpred2d->Add(SBOSSFN2d,-1.0/3);
547 >          Lpred2d->Add(SBOSOFP2d,1.0/3);
548 >          Lpred2d->Add(SBOSOFN2d,-1.0/3);
549 >          
550 >          //3d stuff: flipped prediction
551 >          flippedLpred2d->Add(ZOSOFP2d,-1.0/3);
552 >          flippedLpred2d->Add(ZOSOFN2d,1.0/3);
553 >          flippedLpred2d->Add(SBOSSFP2d,-1.0/3);
554 >          flippedLpred2d->Add(SBOSSFN2d,1.0/3);
555 >          flippedLpred2d->Add(SBOSOFP2d,-1.0/3);
556 >          flippedLpred2d->Add(SBOSOFN2d,1.0/3);
557        } else {
558 <        Lpred->Add(ZOSOFP,1.0);
559 <        Lpred->Add(ZOSOFN,-1.0);
560 <        
561 <        //flipped prediction
562 <        flippedLpred->Add(ZOSOFP,-1.0);
563 <        flippedLpred->Add(ZOSOFN,1.0);
558 >          Lpred->Add(ZOSOFP,1.0);
559 >          Lpred->Add(ZOSOFN,-1.0);
560 >          
561 >          //flipped prediction
562 >          flippedLpred->Add(ZOSOFP,-1.0);
563 >          flippedLpred->Add(ZOSOFN,1.0);
564 >          
565 >          //2d stuff
566 >          Lpred2d->Add(ZOSOFP2d,1.0);
567 >          Lpred2d->Add(ZOSOFN2d,-1.0);
568 >          
569 >          //flipped prediction
570 >          flippedLpred2d->Add(ZOSOFP2d,-1.0);
571 >          flippedLpred2d->Add(ZOSOFN2d,1.0);
572        }
573      }
574      
575      TH1F *signal = (TH1F*)Lobs->Clone("signal");
576 +    TH1F *signal2d = (TH1F*)Lobs2d->Clone("signal2d");
577 +    
578      if(!PlottingSetup::FullMCAnalysis) signal->Add(Lpred,-1);
579      signal->SetName(signalname.c_str());
580      signal->SetTitle(signalname.c_str());
581      signal->Write();
582      
583 +    if(!PlottingSetup::FullMCAnalysis) signal2d->Add(Lpred2d,-1);
584 +    signal2d->SetName((signalname+"2d").c_str());
585 +    signal2d->SetTitle((signalname+"2d").c_str());
586 +    signal2d->Write();
587 +      
588      TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
589      if(!PlottingSetup::FullMCAnalysis) flippedsignal->Add(flippedLpred,-1);
590      flippedsignal->SetName(("flipped_"+signalname).c_str());
591      flippedsignal->Write();
592      
593 +    TH2F *flippedsignal2d = (TH2F*)flippedLobs2d->Clone();
594 +    if(!PlottingSetup::FullMCAnalysis) flippedsignal2d->Add(flippedLpred,-1);
595 +    flippedsignal2d->SetName(("flipped2d_"+signalname).c_str());
596 +    flippedsignal2d->Write();
597 +      
598      if(identifier=="") {
599        TH1F *signalStatUp = (TH1F*)signal->Clone();
600        signalStatUp->SetName(((string)signal->GetName()+"_StatUp").c_str());
# Line 513 | Line 604 | void generate_shapes_for_systematic(bool
604        signalStatDn->SetName(((string)signal->GetName()+"_StatDown").c_str());
605        signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
606        
607 +      TH2F *signalStatUp2d = (TH2F*)signal->Clone();
608 +      signalStatUp2d->SetName(((string)signal2d->GetName()+"_StatUp").c_str());
609 +      signalStatUp2d->SetTitle(((string)signal2d->GetTitle()+"_StatUp").c_str());
610 +      
611 +      TH2F *signalStatDn2d = (TH2F*)signal->Clone();
612 +      signalStatDn2d->SetName(((string)signal2d->GetName()+"_StatDown").c_str());
613 +      signalStatDn2d->SetTitle(((string)signal2d->GetTitle()+"_StatDown").c_str());
614 +      
615 +        
616        for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
617          float staterr;
618          if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
# Line 529 | Line 629 | void generate_shapes_for_systematic(bool
629          signal->SetBinError(i,staterr);
630        }
631        
632 +      for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
633 +        for(int j=1;j<=signalStatDn->GetNbinsY();j++) {
634 +          float staterr;
635 +          if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred2d->GetBinContent(i,j) + Lobs2d->GetBinContent(i,j));
636 +          else staterr = TMath::Sqrt(Lobs2d->GetBinContent(i,j));
637 +          
638 +          if(!PlottingSetup::FullMCAnalysis) {
639 +            dout << "Stat err in bin " << i << ", " << j << " : " << staterr << endl;
640 +            dout << "    prediction: " << Lpred2d->GetBinContent(i,j) << " , observation: " << Lobs2d->GetBinContent(i,j) << " --> signal: " << signal2d->GetBinContent(i,j) << endl;
641 +            dout << "    we obtain : " << signal2d->GetBinContent(i,j)-staterr << " , " << signal2d->GetBinContent(i,j)+staterr << endl;
642 +          }
643 +          if(signal2d->GetBinContent(i,j)-staterr>0) signalStatDn2d->SetBinContent(i,j,signal2d->GetBinContent(i,j)-staterr);
644 +          else signalStatDn2d->SetBinContent(i,j,0);
645 +          signalStatUp2d->SetBinContent(i,signal2d->GetBinContent(i,j)+staterr);
646 +          signal2d->SetBinError(i,j,staterr);
647 +        }
648 +      }
649 +      
650 +    //----------------------------------- ******** GOTTEN HERE *******---------------------------------------
651 +      
652        signalStatDn->Write();
653        signalStatUp->Write();
654        
# Line 943 | Line 1063 | void generate_shapes_for_systematic(bool
1063        delete SBOSOFN;
1064      }
1065    }
1066 <  
1066 >    
1067 >    
1068 >    //----------------------------------- ******** GOTTEN HERE *******---------------------------------------
1069 >    
1070 >    
1071 >    
1072 >    
1073 >
1074   }
1075  
1076   void DoFullCls(string RunDirectory,float firstGuess) {
# Line 983 | Line 1110 | void DoFullCls(string RunDirectory,float
1110    RealScript << "time for i in {";
1111    RealScript << epoints[0]/4 << ","; // adding a VERY VERY low point just to be totally safe
1112    RealScript << epoints[0]/2 << ","; // adding a VERY low point just to be safe
1113 <  for(int i=0;i<epoints.size();i++) RealScript << epoints[i] << ",";
1113 >  for(int i=0;i<(int)epoints.size();i++) RealScript << epoints[i] << ",";
1114    RealScript << 2*epoints[epoints.size()-1] << ","; // adding a VERY high point just to be safe
1115    RealScript << 4*epoints[epoints.size()-1]; // adding a VERY VERY high point just to be totally safe
1116    RealScript << "}; do \n";
# Line 1052 | Line 1179 | ShapeDroplet LimitsFromShapes(float low_
1179    
1180    ensure_directory_exists(RunDirectory);
1181    
1182 <  TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str(),"READ");
1182 >  TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str());
1183    if(datafile->IsZombie()) {
1184      write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
1185      assert(!datafile->IsZombie());
# Line 1183 | Line 1310 | ShapeDroplet LimitsFromShapes(float low_
1310    dout << "Everything is saved in " << RunDirectory << endl;
1311    dout << "Cleaning up ... " << std::flush;
1312    write_warning(__FUNCTION__,"NOT CLEANING UP");
1313 <  gSystem->Exec(("rm -rf "+RunDirectory).c_str());
1313 > //  gSystem->Exec(("rm -rf "+RunDirectory).c_str());
1314    dout << " ok!" << endl;
1315    delete limcan;
1316    return alpha;
1317   }
1318  
1319   void PrepareDataShapes(string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
1320 <  
1320 >  dout << "Preparing data shapes ... " << endl;
1321    map <  pair<float, float>, map<string, float>  >  xsec;
1322    TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
1323    TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
# Line 1207 | Line 1334 | void PrepareDataShapes(string mcjzb, str
1334    generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffDown"),cutnJets,"","EffDown",xsec);
1335    generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatUp",xsec);
1336    generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatDown",xsec);
1337 <  /*
1337 >  
1338    datafile->Close();
1339 <  */
1339 >  
1340   }
1341    

Diff Legend

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