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

Comparing UserCode/cbrown/Development/Plotting/Modules/StudyModule.C (file contents):
Revision 1.1 by buchmann, Mon Jan 30 14:46:26 2012 UTC vs.
Revision 1.10 by buchmann, Wed Aug 15 13:52:13 2012 UTC

# Line 16 | Line 16
16   #include <TProfile.h>
17   #include <TPaveStats.h>
18  
19 + #include "GeneratorLevelStudyModule.C"
20 +
21   //#include "TTbar_stuff.C"
22   using namespace std;
23  
# Line 36 | Line 38 | void do_experimental_pred_obs_calculatio
38    TH1F *SBOSOFN;
39    
40    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
41 <  if(PlottingSetup::RestrictToMassPeak) {
41 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
42      SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
43      SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
44      SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
# Line 45 | Line 47 | void do_experimental_pred_obs_calculatio
47      
48      
49    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
50 <  if(PlottingSetup::RestrictToMassPeak) {
50 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
51      dout << "   Observed : " << ZOSSFP->Integral() << endl;
52      dout << "   Predicted: " << ZOSSFN->Integral() << " + (1/3)*(" << ZOSOFP->Integral() << "-" << ZOSOFN->Integral()<<") + (1/3)*(" << SBOSSFP->Integral() << "-" << SBOSSFN->Integral()<<") + (1/3)*(" << SBOSOFP->Integral() << "-" << SBOSOFN->Integral()<<")" << endl;
53      dout << "        P(ZJets ) \t " << ZOSSFN->Integral() << endl;
# Line 74 | Line 76 | void do_experimental_pred_obs_calculatio
76   void look_at_sidebands(string mcjzb, string datajzb, bool includejetcut, float cutat=0) {
77    
78    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak -- this funciton is meaningless for the offpeak case
79 <  if(!PlottingSetup::RestrictToMassPeak) return;
79 >  if(!PlottingSetup::RestrictToMassPeak||!PlottingSetup::UseSidebandsForcJZB) return;
80    dout << "Looking at sidebands ... " << endl;
81    int mcordata=data;//data     // you can perform this study for mc or data ...
82    
# Line 244 | Line 246 | int i=0;
246    
247   void find_sideband_definition() {
248    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
249 <  if(!PlottingSetup::RestrictToMassPeak) return; // this function is meaningless for the offpeak analysis
249 >  if(!PlottingSetup::RestrictToMassPeak||!PlottingSetup::UseSidebandsForcJZB) return; // this function is meaningless for the offpeak analysis
250  
251    TH1F *mllttbar  = allsamples.Draw("mllttbar","mll",145,55,200, "mll [GeV]", "events",cutOSSF&&cutnJets&&!cutmass,mc,luminosity,allsamples.FindSample("TTJets"));
252    TH1F *mllttbarz  = allsamples.Draw("mllttbarz","mll",1,50,200, "mll [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("TTJets"));
# Line 318 | Line 320 | void run_check() {
320    TCanvas *c1 = new TCanvas("runnum","runnum",800,1000);
321    c1->Divide(2,4);
322    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
323 <  if(PlottingSetup::RestrictToMassPeak) c1->Divide(2,2); // there are only four regions ...
323 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) c1->Divide(2,2); // there are only four regions ...
324  
325    c1->cd(1);
326    TH1F *ossfp = runcheckhisto((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
# Line 340 | Line 342 | void run_check() {
342    TText *t4 = write_title("OSOF,N");t4->Draw();
343    
344    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
345 <  if(PlottingSetup::RestrictToMassPeak) {
345 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
346      c1->cd(5);
347      TH1F *sbofp = runcheckhisto((const char*)(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
348      sbofp->Draw();
# Line 376 | Line 378 | TH1F *give_boson_pred(TCut bcut,string m
378    TH1F *jzbpss;
379    
380    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
381 <  if(PlottingSetup::RestrictToMassPeak) {
381 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
382      jzbnos  = allsamples.Draw("jzbnos","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",bcut&&cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
383      jzbpos  = allsamples.Draw("jzbpos",mcjzb,nbins,0,350, "JZB [GeV]", "events",    bcut&&cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
384      jzbnss  = allsamples.Draw("jzbnss","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",bcut&&cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
# Line 386 | Line 388 | TH1F *give_boson_pred(TCut bcut,string m
388      
389    TH1F *pred = (TH1F*)jzbn->Clone("pred");
390    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
391 <  if(PlottingSetup::RestrictToMassPeak) {
391 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
392      pred->Add(jzbno,-1.0/3);
393      pred->Add(jzbpo,1.0/3);
394      pred->Add(jzbnos,-1.0/3);
# Line 610 | Line 612 | TH1F *give_lm0_pred(TCut bcut,string mcj
612    TH1F *jzbnss;
613    TH1F *jzbpss;
614    
615 <  if(PlottingSetup::RestrictToMassPeak) {
615 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
616      jzbnos  = signalsamples.Draw("jzbnos","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
617      jzbpos  = signalsamples.Draw("jzbpos",mcjzb,nbins,0,350, "JZB [GeV]", "events",    cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
618      jzbnss  = signalsamples.Draw("jzbnss","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
# Line 619 | Line 621 | TH1F *give_lm0_pred(TCut bcut,string mcj
621    
622    TH1F *pred = (TH1F*)jzbn->Clone("pred");
623    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
624 <  if(PlottingSetup::RestrictToMassPeak) {
624 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
625      pred->Add(jzbno,-1.0/3);
626      pred->Add(jzbpo,1.0/3);
627      pred->Add(jzbnos,-1.0/3);
# Line 757 | Line 759 | void kinematic_dist_of_pred_and_obs() {/
759    make_double_plot("pfJetGoodNum",8,0.5,8.5,20,nolog,"nJets","nJets_MC",observed,predictedMC,mc);
760  
761   }
762 +
763 +
764 + void generator_study_plots() {
765 +
766 + //uncomment whichever one you want to see :-)
767 +  
768 + /*
769 + GenLevelStudy::X_vs_generation_lm4();
770 + GenLevelStudy::compare_sms_lm4();
771 + */
772 +
773 + ///GenLevelStudy::MomentumFraction();
774 + ///GenLevelStudy::AngleLSPLSP();
775 + ///GenLevelStudy::DeltaLSPmomentum();
776 + ///GenLevelStudy::ZDecayIllustration();
777 + /*
778 + GenLevelStudy::AngleMETsumLSP();
779 + GenLevelStudy::RatioMETsumLSP();
780 + GenLevelStudy::AngleLSPZ();
781 + GenLevelStudy::WidthIllustration();
782 +
783 +
784 + GenLevelStudy::pStarIllustration(0.5);
785 + //GenLevelStudy::pStarIllustration(0.25);
786 + //GenLevelStudy::pStarIllustration(0.75);
787 + GenLevelStudy::DrawJetBand();
788 + GenLevelStudy::DrawOnly100to150inJetBand();*/
789 + GenLevelStudy::ImpactOfGluinoChi2MassDifference();
790 +
791 + }
792 +
793 +
794 + void jzb_negative_generator_study() {
795 +  write_warning(__FUNCTION__,"We are going to use a t5zz scan file, and \033[1;31m WON'T \033[1;35m cut on MassGlu/MassLSP in order to improve statistics. This is ok for small studies, for a real study you'll need to look at points individually ...");
796 +
797 +
798 +  scansample.AddSample("/scratch/buchmann/ntuples/GeneratorInformationInJZB___JZBplusSamples_TestingSMS_v7/SMS-T5zz_x-05_Mgluino-150to1200_mLSP-50to1150_7TeV-Pythia6Z__Summer11-PU_START42_V11_FastSim-v2__AODSIM___inclindex_v2.root","SMST5zz",1,1,false,false,0,kRed);
799 +  TCanvas *jcan = new TCanvas("jcan","jcan");
800 +  scansample.collection[scansample.collection.size()-1].events->Draw("(LSP1pt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF");
801 +  TH1F *h = new TH1F("h","h",100,-500,500);
802 +  h->SetLineColor(kBlack);
803 +  TProfile *p = (TProfile*)jcan->GetPrimitive("htemp");
804 +  p->GetXaxis()->SetTitle("Generator JZB");
805 +  p->GetXaxis()->CenterTitle();
806 +  p->GetYaxis()->SetTitle("( LSP p_{T} ) / ( LSP mother p_{T} )");
807 +  p->GetYaxis()->CenterTitle();
808 +  p->SetLineColor(kBlue);
809 +  TLegend* leg = make_legend("", 0.40, 0.75, false);
810 +  leg->AddEntry(p,"LSP1 pt / LSP1 Mother pt","l");
811 +  leg->AddEntry(h,"Z pt / Z Mother pt","l");
812 +  leg->Draw();
813 +  TText *title = write_title("JZB as a function of the first LSP's momentum transfer");
814 +  title->Draw();
815 +  scansample.collection[scansample.collection.size()-1].events->Draw("(genZPt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF,same");
816 +
817 +  CompleteSave(jcan,"NegativeJZBStudy/LSPpt_LSPMopt");
818 +
819 +  scansample.collection[scansample.collection.size()-1].events->Draw("angleLSPLSP:pureGeneratorJZB","genNjets>2","PROF");
820 +  TProfile *p1 = (TProfile*)jcan->GetPrimitive("htemp");
821 +  p1->GetXaxis()->SetTitle("Generator JZB");
822 +  p1->GetXaxis()->CenterTitle();
823 +  p1->GetYaxis()->SetTitle("#angle(LSP1,LSP2)");
824 +  p1->GetYaxis()->CenterTitle();
825 +  TText *title1 = write_title("JZB as a function of the angle between the two LSPs");
826 +  title1->Draw();
827 +  CompleteSave(jcan,"NegativeJZBStudy/AngleLSPLSP");
828 +
829 +  TH1F *jzbdistributionvsz[5];
830 +  THStack zstack;
831 +  jcan->SetLogy(1);
832 +  TLegend* leg2 = make_legend("", 0.15, 0.75, false);
833 +  leg2->SetX2(0.4);
834 +  TLegend* leg3 = make_legend("", 0.15, 0.75, false);
835 +  leg3->SetX2(0.4);
836 +  for(int z=0;z<5;z++) {
837 +   stringstream specialcut;
838 +   if(z<4) specialcut << "genNjets>2&&(LSP1pt/LSP1Mopt)>" << z*0.2 << "&&(LSP1pt/LSP1Mopt)<" << (z+1)*0.2;
839 +   else specialcut << "genNjets>2&&(LSP1pt/LSP1Mopt)>" << z*0.2;
840 +   stringstream histtitle;
841 +   histtitle << "splitup_" << z;
842 +   stringstream ntitle;
843 +   if(z<4) ntitle << z*0.2 << " < z < " << (z+1)*0.2;
844 +   else ntitle << z*0.2 << " < z";
845 +   jzbdistributionvsz[z] = scansample.Draw(histtitle.str().c_str(),"pureGeneratorJZB",100,-500,500, "generator JZB (GeV)", "events",specialcut.str().c_str(),mc,1.0);
846 +   jzbdistributionvsz[z]->SetLineColor(z+1);
847 +   jzbdistributionvsz[z]->SetMarkerSize(0);
848 +   zstack.Add(jzbdistributionvsz[z]);
849 +   leg2->AddEntry(jzbdistributionvsz[z],ntitle.str().c_str(),"f");
850 +   leg3->AddEntry(jzbdistributionvsz[z],ntitle.str().c_str(),"l");
851 +  }
852 +
853 +  jzbdistributionvsz[0]->GetYaxis()->SetRangeUser(2,800);
854 +  jzbdistributionvsz[0]->DrawNormalized("histo");
855 +  for(int z=0;z<5;z++) {
856 +   jzbdistributionvsz[z]->DrawNormalized("histo,same");
857 +  }
858 +
859 + //  zstack.Draw("nostack,histo");
860 +  leg3->Draw("same");
861 +  CompleteSave(jcan,"NegativeJZBStudy/StackedAccordingToMomentumFractionIndividual");
862 +
863 +  for(int z=0;z<5;z++) {
864 +   jzbdistributionvsz[z]->SetFillColor(z+1);
865 +  }
866 +
867 +  zstack.Draw("histo");
868 +  leg2->Draw("same");
869 +  CompleteSave(jcan,"NegativeJZBStudy/StackedAccordingToMomentumFractionStacked");
870 +
871 + //varangle vasysyn
872 +
873 + //  scansample.collection[scansample.collection.size()-1].events->Draw("(LSPPromptnessLevel[0]/4.0)*angleLSPLSP/(LSP1pt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF");
874 +
875 +
876 +
877 + }
878 +
879 + void compare_lm4_sms_variable(TTree *eventsLM4, TTree *eventsSMS, string variable, int nbins, float xmin, float xmax, TCut cut, string saveas, bool dology=false) {
880 +  TCanvas *can = new TCanvas("can","can");
881 +  can->SetLogy(dology);
882 +  TH1F *hlm4 = new TH1F("hlm4","hlm4",nbins,xmin,xmax);
883 +  TH1F *hsms = new TH1F("hsms","hsms",nbins,xmin,xmax);
884 +  eventsLM4->Draw((variable+">>hlm4").c_str(),cut,"goff");
885 +  eventsSMS->Draw((variable+">>hsms").c_str(),cut,"goff");
886 +  hlm4->SetLineColor(kBlue);
887 +  hsms->SetLineColor(kRed);
888 +
889 +  if(hlm4->Integral()>0) hlm4->Scale(1.0/hlm4->Integral());
890 +  else write_warning(__FUNCTION__,"Watch out, LM4 histo is empty!");
891 +  if(hsms->Integral()>0) hsms->Scale(1.0/hsms->Integral());
892 +  else write_warning(__FUNCTION__,"Watch out, SMS histo is empty!");
893 +
894 +  float min=get_nonzero_minimum(hlm4);
895 +  float max=hlm4->GetMaximum();
896 +  if(get_nonzero_minimum(hsms)<min) min=get_nonzero_minimum(hsms);
897 +  if(hsms->GetMaximum()>max) max=hsms->GetMaximum();
898 +  if(dology) max*=5;
899 +  else max*=2;
900 +
901 +  hlm4->GetYaxis()->SetRangeUser(min,max);
902 +  hlm4->GetXaxis()->SetTitle(variable.c_str());
903 +  hlm4->GetXaxis()->CenterTitle();
904 +  hlm4->Draw("histo");
905 +  hsms->Draw("histo,same");
906 +  TLegend *leg = make_legend("",0.2,0.98,false);
907 +  leg->SetY2(1.0);
908 +  leg->SetNColumns(2);
909 +  leg->AddEntry(hlm4,"LM4","l");
910 +  leg->AddEntry(hsms,"\"LM4\" SMS","l");
911 +  leg->Draw();
912 +  stringstream saveas2;
913 +  saveas2 << "ComparingLM4_SMS/" << removefunnystring(saveas);
914 +  CompleteSave(can,saveas2.str());
915 +  delete can;
916 +  delete hlm4;
917 +  delete hsms;
918 + }
919 +
920 +
921 + void compare_LM4_and_SMS() {
922 +  TFile *f1 = new TFile("/shome/lbaeni/jzb/LM4_SMS/SMS_LM4_JZB.root");
923 +  TTree *LM4events = (TTree*)f1->Get("events");
924 +  TFile *f2 = new TFile("/scratch/buchmann/ntuples/GeneratorInformationInJZB___JZBplusSamples_TestingSMS_v5/LM4_SUSY_sftsht_7TeV-pythia6__Summer11-PU_S4_START42_V11-v2__withIndex.root");
925 +  TTree *SMSevents = (TTree*)f2->Get("events");
926 +  
927 +  compare_lm4_sms_variable(LM4events, SMSevents, "mll",100,50,150,cutOSSF&&cutnJets&&basiccut,"mll",true);
928 +  compare_lm4_sms_variable(LM4events, SMSevents, "jzb[1]+0.04*pt-2.32212",100,-300,700,cutmass&&cutOSSF&&cutnJets&&basiccut,"jzb",true);
929 +  compare_lm4_sms_variable(LM4events, SMSevents, "pureGeneratorJZB",100,-300.0,700.0,cutOSSF&&basiccut,"pureGeneratorJZB",true);
930 +  compare_lm4_sms_variable(LM4events, SMSevents, "pfJetGoodNum",10,-0.5,9.5,cutmass&&cutOSSF&&basiccut,"pfJetGoodNum",true);
931 +  compare_lm4_sms_variable(LM4events, SMSevents, "pt",100,15.0,200.0,cutOSSF&&basiccut,"pt",false);
932 +  compare_lm4_sms_variable(LM4events, SMSevents, "pt1",100,15.0,100.0,cutOSSF&&basiccut,"pt1",false);
933 +  compare_lm4_sms_variable(LM4events, SMSevents, "pt2",100,15.0,100.0,cutOSSF&&basiccut,"pt2",false);
934 +  compare_lm4_sms_variable(LM4events, SMSevents, "met[4]",100,0.0,600.0,cutOSSF&&basiccut,"met",false);
935 +  compare_lm4_sms_variable(LM4events, SMSevents, "genMET",100,0.0,600.0,cutOSSF&&basiccut,"genMET",false);
936 +  compare_lm4_sms_variable(LM4events, SMSevents, "genNjets",10,-0.5,9.5,cutOSSF&&basiccut,"genNjets",true);
937 + }
938 +
939 +
940 + void compare_onpeak_offpeak_signal_distributions(string mcjzb, string datajzb, int sampleindex) {
941 +  cout << (signalsamples.collection)[sampleindex].filename << endl;
942 + }
943 +
944 + void compare_onpeak_offpeak_signal_distributions(string mcjzb, string datajzb) {
945 +  for(int i=0;i<signalsamples.collection.size();i++) compare_onpeak_offpeak_signal_distributions(mcjzb,datajzb,i);
946 + }
947 +
948 + bool load_adequate_npred_nobs(string code){
949 +
950 + if(code=="201540") {
951 +  Nobs[0]=1108.67;
952 +  Npred[0]=1106.69;
953 +  Nprederr[0]=275.424852964089;
954 +  Nobs[1]=291.679;
955 +  Npred[1]=288.979;
956 +  Nprederr[1]=74.0248355630055;
957 +  Nobs[2]=76.0776;
958 +  Npred[2]=76.8332;
959 +  Nprederr[2]=20.9768252604726;
960 +  Nobs[3]=22.4581;
961 +  Npred[3]=22.5385;
962 +  Nprederr[3]=7.46758781734772;
963 +  Nobs[4]=5.41088;
964 +  Npred[4]=9.51903;
965 +  Nprederr[4]=3.86595196996807;
966 +  return true;
967 + }
968 +
969 + if(code=="201040") {
970 +  Nobs[0]=1210.02;
971 +  Npred[0]=1228.62;
972 +  Nprederr[0]=305.642034977193;
973 +  Nobs[1]=324.008;
974 +  Npred[1]=323.042;
975 +  Nprederr[1]=82.6094009556418;
976 +  Nobs[2]=84.8866;
977 +  Npred[2]=86.4243;
978 +  Nprederr[2]=23.4841694480367;
979 +  Nobs[3]=25.4339;
980 +  Npred[3]=25.6437;
981 +  Nprederr[3]=8.32541119033168;
982 +  Nobs[4]=6.24946;
983 +  Npred[4]=9.97731;
984 +  Nprederr[4]=4.07510841593202;
985 +  return true;
986 + }
987 +
988 + if(code=="151040") {
989 +  Nobs[0]=1218.67;
990 +  Npred[0]=1239.47;
991 +  Nprederr[0]=308.268321394268;
992 +  Nobs[1]=326.245;
993 +  Npred[1]=325.731;
994 +  Nprederr[1]=83.2795749634327;
995 +  Nobs[2]=85.3292;
996 +  Npred[2]=87.239;
997 +  Nprederr[2]=23.6875859996856;
998 +  Nobs[3]=25.4339;
999 +  Npred[3]=25.9456;
1000 +  Nprederr[3]=8.40219371514963;
1001 +  Nobs[4]=6.24946;
1002 +  Npred[4]=9.97731;
1003 +  Nprederr[4]=4.07510841593202;
1004 +  return true;
1005 + }
1006 +
1007 + if(code=="151550") {
1008 +  Nobs[0]=1009.57;
1009 +  Npred[0]=1015.24;
1010 +  Nprederr[0]=252.027368634063;
1011 +  Nobs[1]=266.851;
1012 +  Npred[1]=264.99;
1013 +  Nprederr[1]=67.9766170833765;
1014 +  Nobs[2]=68.2042;
1015 +  Npred[2]=69.3185;
1016 +  Nprederr[2]=19.0998669414868;
1017 +  Nobs[3]=20.873;
1018 +  Npred[3]=20.1324;
1019 +  Nprederr[3]=6.85046567121535;
1020 +  Nobs[4]=4.67275;
1021 +  Npred[4]=8.42002;
1022 +  Nprederr[4]=3.55560434037591;
1023 +  return true;
1024 + }
1025 +
1026 + if(code=="151540") {
1027 +  Nobs[0]=1112.51;
1028 +  Npred[0]=1111.01;
1029 +  Nprederr[0]=276.447968514945;
1030 +  Nobs[1]=292.622;
1031 +  Npred[1]=290.731;
1032 +  Nprederr[1]=74.461488188526;
1033 +  Nobs[2]=76.1949;
1034 +  Npred[2]=77.1915;
1035 +  Nprederr[2]=21.0662269219811;
1036 +  Nobs[3]=22.4581;
1037 +  Npred[3]=22.8404;
1038 +  Nprederr[3]=7.54484166514447;
1039 +  Nobs[4]=5.41088;
1040 +  Npred[4]=9.51903;
1041 +  Nprederr[4]=3.86595196996807;
1042 +  return true;
1043 + }
1044 +
1045 + if(code=="151530") {
1046 +  Nobs[0]=1182.29;
1047 +  Npred[0]=1177.48;
1048 +  Nprederr[0]=293.397396547567;
1049 +  Nobs[1]=310.125;
1050 +  Npred[1]=308.279;
1051 +  Nprederr[1]=78.8363941324691;
1052 +  Nobs[2]=83.376;
1053 +  Npred[2]=80.804;
1054 +  Nprederr[2]=21.9686562734479;
1055 +  Nobs[3]=24.9171;
1056 +  Npred[3]=23.4465;
1057 +  Nprederr[3]=7.69983840500565;
1058 +  Nobs[4]=6.151;
1059 +  Npred[4]=9.7932;
1060 +  Nprederr[4]=3.94255080734542;
1061 +  return true;
1062 + }
1063 +
1064 + if(code=="101040") {
1065 +  Nobs[0]=1219.78;
1066 +  Npred[0]=1241.65;
1067 +  Nprederr[0]=308.795385452892;
1068 +  Nobs[1]=326.577;
1069 +  Npred[1]=327.172;
1070 +  Nprederr[1]=83.6388357552878;
1071 +  Nobs[2]=85.3292;
1072 +  Npred[2]=87.5409;
1073 +  Nprederr[2]=23.7629564995099;
1074 +  Nobs[3]=25.4339;
1075 +  Npred[3]=26.2475;
1076 +  Nprederr[3]=8.47896339395919;
1077 +  Nobs[4]=6.24946;
1078 +  Npred[4]=9.97731;
1079 +  Nprederr[4]=4.07510841593202;
1080 +  return true;
1081 + }
1082 +
1083 + if(code=="101030") {
1084 +  Nobs[0]=1335.36;
1085 +  Npred[0]=1342.37;
1086 +  Nprederr[0]=334.276222135228;
1087 +  Nobs[1]=355.554;
1088 +  Npred[1]=355.065;
1089 +  Nprederr[1]=90.5352675557984;
1090 +  Nobs[2]=94.9863;
1091 +  Npred[2]=93.7287;
1092 +  Nprederr[2]=25.3082516747483;
1093 +  Nobs[3]=28.3785;
1094 +  Npred[3]=28.3783;
1095 +  Nprederr[3]=9.02022017545026;
1096 +  Nobs[4]=6.98957;
1097 +  Npred[4]=10.8503;
1098 +  Nprederr[4]=4.31460183234792;
1099 +  return true;
1100 + }
1101 +
1102 + if(code=="101020") {
1103 +  Nobs[0]=1418.83;
1104 +  Npred[0]=1421.23;
1105 +  Nprederr[0]=353.214726699284;
1106 +  Nobs[1]=387.594;
1107 +  Npred[1]=380.483;
1108 +  Nprederr[1]=96.9195932538411;
1109 +  Nobs[2]=106.565;
1110 +  Npred[2]=99.4694;
1111 +  Nprederr[2]=26.7421279631969;
1112 +  Nobs[3]=32.9309;
1113 +  Npred[3]=31.6286;
1114 +  Nprederr[3]=9.84399940108186;
1115 +  Nobs[4]=9.25395;
1116 +  Npred[4]=12.5158;
1117 +  Nprederr[4]=4.76585512033255;
1118 +  return true;
1119 + }
1120 +
1121 + if(code=="202040") {
1122 +  Nobs[0]=960.506;
1123 +  Npred[0]=954.77;
1124 +  Nprederr[0]=238.5811858172;
1125 +  Nobs[1]=251.453;
1126 +  Npred[1]=250.838;
1127 +  Nprederr[1]=64.2967104354;
1128 +  Nobs[2]=64.4903;
1129 +  Npred[2]=66.6518;
1130 +  Nprederr[2]=18.4320150537;
1131 +  Nobs[3]=19.1953;
1132 +  Npred[3]=20.1075;
1133 +  Nprederr[3]=6.844063617;
1134 +  Nobs[4]=5.09139;
1135 +  Npred[4]=8.67175;
1136 +  Nprederr[4]=3.6271903178;
1137 +  return true;
1138 + }
1139 +
1140 + if(code=="onpeak") {
1141 +  Nobs[0]=387.268;
1142 +  Npred[0]=377.146;
1143 +  Nprederr[0]=55.4000603791;
1144 +  Nobs[1]=93.6071;
1145 +  Npred[1]=89.2911;
1146 +  Nprederr[1]=13.7946793213;
1147 +  Nobs[2]=22.4222;
1148 +  Npred[2]=22.2462;
1149 +  Nprederr[2]=4.252924039;
1150 +  Nobs[3]=7.49126;
1151 +  Npred[3]=6.98267;
1152 +  Nprederr[3]=1.913797227;
1153 +  Nobs[4]=1.80426;
1154 +  Npred[4]=3.00437;
1155 +  Nprederr[4]=1.2381643019;
1156 +  return true;
1157 + }
1158 +
1159 + if(code=="151550") {
1160 +  Nobs[0]=1009.57;
1161 +  Npred[0]=1015.24;
1162 +  Nprederr[0]=252.0273686341;
1163 +  Nobs[1]=266.851;
1164 +  Npred[1]=264.99;
1165 +  Nprederr[1]=67.9766170834;
1166 +  Nobs[2]=68.2042;
1167 +  Npred[2]=69.3185;
1168 +  Nprederr[2]=19.0998669415;
1169 +  Nobs[3]=20.873;
1170 +  Npred[3]=20.1324;
1171 +  Nprederr[3]=6.8504656712;
1172 +  Nobs[4]=4.67275;
1173 +  Npred[4]=8.42002;
1174 +  Nprederr[4]=3.5556043404;
1175 +  return true;
1176 + }
1177 +
1178 + if(code=="101040") {
1179 +  Nobs[0]=1108.67;
1180 +  Npred[0]=1106.69;
1181 +  Nprederr[0]=275.4248529641;
1182 +  Nobs[1]=291.679;
1183 +  Npred[1]=288.979;
1184 +  Nprederr[1]=74.024835563;
1185 +  Nobs[2]=76.0776;
1186 +  Npred[2]=76.8332;
1187 +  Nprederr[2]=20.9768252605;
1188 +  Nobs[3]=22.4581;
1189 +  Npred[3]=22.5385;
1190 +  Nprederr[3]=7.4675878173;
1191 +  Nobs[4]=5.41088;
1192 +  Npred[4]=9.51903;
1193 +  Nprederr[4]=3.86595197;
1194 +  return true;
1195 + }
1196 +
1197 + if(code=="151040") {
1198 +  Nobs[0]=1218.67;
1199 +  Npred[0]=1239.47;
1200 +  Nprederr[0]=308.2683213943;
1201 +  Nobs[1]=326.245;
1202 +  Npred[1]=325.731;
1203 +  Nprederr[1]=83.2795749634;
1204 +  Nobs[2]=85.3292;
1205 +  Npred[2]=87.239;
1206 +  Nprederr[2]=23.6875859997;
1207 +  Nobs[3]=25.4339;
1208 +  Npred[3]=25.9456;
1209 +  Nprederr[3]=8.4021937151;
1210 +  Nobs[4]=6.24946;
1211 +  Npred[4]=9.97731;
1212 +  Nprederr[4]=4.0751084159;
1213 +  return true;
1214 + }
1215 +
1216 + if(code=="201040") {
1217 +  Nobs[0]=1210.02;
1218 +  Npred[0]=1228.62;
1219 +  Nprederr[0]=305.6420349772;
1220 +  Nobs[1]=324.008;
1221 +  Npred[1]=323.042;
1222 +  Nprederr[1]=82.6094009556;
1223 +  Nobs[2]=84.8866;
1224 +  Npred[2]=86.4243;
1225 +  Nprederr[2]=23.484169448;
1226 +  Nobs[3]=25.4339;
1227 +  Npred[3]=25.6437;
1228 +  Nprederr[3]=8.3254111903;
1229 +  Nobs[4]=6.24946;
1230 +  Npred[4]=9.97731;
1231 +  Nprederr[4]=4.0751084159;
1232 +  return true;
1233 + }
1234 +
1235 + if(code=="151540") {
1236 +  Nobs[0]=1112.51;
1237 +  Npred[0]=1111.01;
1238 +  Nprederr[0]=276.4479685149;
1239 +  Nobs[1]=292.622;
1240 +  Npred[1]=290.731;
1241 +  Nprederr[1]=74.4614881885;
1242 +  Nobs[2]=76.1949;
1243 +  Npred[2]=77.1915;
1244 +  Nprederr[2]=21.066226922;
1245 +  Nobs[3]=22.4581;
1246 +  Npred[3]=22.8404;
1247 +  Nprederr[3]=7.5448416651;
1248 +  Nobs[4]=5.41088;
1249 +  Npred[4]=9.51903;
1250 +  Nprederr[4]=3.86595197;
1251 +  return true;
1252 + }
1253 +
1254 +
1255 + return false;
1256 + }
1257 +
1258 + vector<float> do_simulate_upper_limits(string mcjzb, string datajzb, float MCPeakError, vector<float> jzb_cut, string code) {
1259 +  vector<vector<float> > all_systematics;
1260 +  vector<float> bestUL;
1261 +  //we need to set the correct npred, nprederr, and nobs.
1262 +  bool loaded_info = load_adequate_npred_nobs(code);
1263 +  if(!loaded_info) {
1264 +    write_warning(__FUNCTION__,"obs/pred/prederr could not be loaded for this configuration.");
1265 +    cout << "Configuration " << code << " (pt1/pt2/mll) caused problems. will be skipped." << endl;
1266 +  } else {
1267 +    cout << "Loaded configuration for code " << code << " successfully, e.g. Npred[0] is " << Npred[0] << endl;
1268 +    bool doobserved=false;
1269 +    all_systematics=compute_systematics(mcjzb,MCPeakError,alwaysflip,datajzb,signalsamples,jzb_cut,requireZ);
1270 +    bestUL = compute_upper_limits_from_counting_experiment(all_systematics,jzb_cut,mcjzb,doobserved,alwaysflip);
1271 +    
1272 +  }
1273 +  return bestUL;
1274 + }  
1275 +
1276 + void decipherM0M12sample(string samplename, float &M0, float &M12) {
1277 +  int position = samplename.find("M0_");
1278 +  string interestingpart = samplename.substr(position+3,4);
1279 +  position = interestingpart.find("_");
1280 +  if(position>0&&position<5) interestingpart=interestingpart.substr(0,position);
1281 +  stringstream M0a;
1282 +  M0a << interestingpart;
1283 +  M0a>>M0;
1284 +  position = samplename.find("M12_");
1285 +  interestingpart = samplename.substr(position+4,4);
1286 +  position = interestingpart.find("_");
1287 +  if(position>0&&position<5) interestingpart=interestingpart.substr(0,position);
1288 +  stringstream M12a;
1289 +  M12a << interestingpart;
1290 +  M12a>>M12;
1291 + }
1292 +
1293 + vector<TGraph*> create_limit_graphs(vector<float> limits) {
1294 +  float x[limits.size()];
1295 +  float y[limits.size()];
1296 +  int nx=0;
1297 +  float x2[limits.size()];
1298 +  float y2[limits.size()];
1299 +  int nx2=0;
1300 +  for(int i=0;i<limits.size();i++) {
1301 +    string samplename=signalsamples.collection[i].samplename;
1302 +    float m0,m12;
1303 +    decipherM0M12sample(samplename,m0,m12);
1304 +    if(m12==300) {
1305 +      x[nx]=m0;
1306 +      y[nx]=limits[i];
1307 +      nx++;
1308 +    }
1309 +    if(m12==400) {
1310 +      x2[nx2]=m0;
1311 +      y2[nx2]=limits[i];
1312 +      nx2++;
1313 +    }
1314 +  }
1315 +  
1316 +  TGraph *gra = new TGraph(nx,x,y);
1317 +  TGraph *grb = new TGraph(nx2,x2,y2);
1318 +  vector<TGraph*> graphs;
1319 +  graphs.push_back(gra);
1320 +  graphs.push_back(grb);
1321 +  return graphs;
1322 + }
1323 +  
1324 +  
1325 +
1326 + void interpret_bestULs_for_cuts(vector<pair<vector<float>,string> > bestUL) {
1327 +  TCanvas *can = new TCanvas("can","can",1800,900);
1328 +  can->Divide(2,1);
1329 +  can->cd(1);
1330 +  TGraph *graphs[2][bestUL.size()];
1331 +  TLegend *leg = make_legend();
1332 +  for(int i=0;i<bestUL.size();i++) {
1333 +    vector<TGraph*> grs=create_limit_graphs(bestUL[i].first);
1334 +    graphs[0][i]=grs[0];
1335 +    graphs[1][i]=grs[1];
1336 +    graphs[0][i]->SetLineColor(GenLevelStudy::diversehistocolor(i));
1337 +    graphs[1][i]->SetLineColor(GenLevelStudy::diversehistocolor(i));
1338 +    graphs[0][i]->SetMarkerColor(GenLevelStudy::diversehistocolor(i));
1339 +    graphs[1][i]->SetMarkerColor(GenLevelStudy::diversehistocolor(i));
1340 +    leg->AddEntry(graphs[0][i],bestUL[i].second.c_str(),"lp");
1341 +  }
1342 +  
1343 +  
1344 +  for(int i=0;i<bestUL.size();i++) {
1345 +    graphs[0][i]->GetXaxis()->SetTitle("M_{0}");
1346 +    graphs[0][i]->GetYaxis()->SetTitle("Upper Limit (M_{1/2}=300)");
1347 +    graphs[0][i]->GetXaxis()->CenterTitle();
1348 +    graphs[0][i]->GetYaxis()->CenterTitle();
1349 +    
1350 +    graphs[1][i]->GetXaxis()->SetTitle("M_{0}");
1351 +    graphs[1][i]->GetYaxis()->SetTitle("Upper Limit (M_{1/2}=400)");
1352 +    graphs[1][i]->GetXaxis()->CenterTitle();
1353 +    graphs[1][i]->GetYaxis()->CenterTitle();
1354 +    
1355 +    can->cd(1);
1356 +    if(i==0) graphs[0][i]->Draw("ALP");
1357 +    else graphs[0][i]->Draw("LP");
1358 +    can->cd(2);
1359 +    if(i==0) graphs[1][i]->Draw("ALP");
1360 +    else graphs[1][i]->Draw("LP");
1361 +  }
1362 +  can->cd(1);
1363 +  leg->Draw();
1364 +  can->cd(2);
1365 +  leg->Draw();
1366 +  CompleteSave(can,"mSUGRAlimits");
1367 +  delete can;
1368 + }
1369 +
1370 + void do_simulate_upper_limits(string mcjzb, string datajzb, float MCPeakError, vector<float> jzb_cut) {
1371 +  
1372 +  cout << "Going to simulate upper limits (i.e. expected limits) USING MONTE CARLO!!!!" << endl;
1373 +  while(jzb_cut.size()>Nobs.size()||jzb_cut.size()>Npred.size()||jzb_cut.size()>Nprederr.size()) {
1374 +    Nobs.push_back(0);
1375 +    Npred.push_back(0);
1376 +    Nprederr.push_back(0);
1377 +  }
1378 +  
1379 +  vector<pair<vector<float>,string> > bestUL;
1380 +  
1381 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"201540"),"201540"));
1382 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"201040"),"201040"));/*
1383 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151040"),"151040"));
1384 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151550"),"151550"));
1385 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151540"),"151540"));
1386 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151530"),"151530"));
1387 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"101040"),"101040"));
1388 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"101030"),"101030"));
1389 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"101020"),"101020"));
1390 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151540"),"151540"));
1391 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151550"),"151550"));
1392 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"201040"),"201040"));
1393 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151040"),"151040"));
1394 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"101040"),"101040"));
1395 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"onpeak"),"onpeak"));
1396 +  bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"202040"),"202040"));
1397 + */
1398 + /*
1399 +  //-------------------
1400 +  write_warning(__FUNCTION__,"Just testing the components now ...");
1401 +  vector<float> limits;
1402 +  for(int i=0;i<signalsamples.collection.size();i++) limits.push_back(1/(0.01*(i+1)));
1403 +  vector<float> limits2;
1404 +  for(int i=0;i<signalsamples.collection.size();i++) limits2.push_back(1/(0.01*(2*i+2)));
1405 +  bestUL.push_back(make_pair(limits,"abcd"));
1406 +  bestUL.push_back(make_pair(limits2,"efgh"));
1407 +  
1408 +  //-------------------
1409 +  */
1410 +  
1411 +  
1412 +  interpret_bestULs_for_cuts(bestUL);
1413 + }

Diff Legend

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