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.4 by buchmann, Mon Feb 20 17:06:18 2012 UTC vs.
Revision 1.10 by buchmann, Wed Aug 15 13:52:13 2012 UTC

# Line 38 | 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 47 | 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 76 | 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 246 | 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 320 | 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 342 | 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 378 | 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 388 | 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 612 | 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 621 | 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 768 | Line 768 | void generator_study_plots() {
768   /*
769   GenLevelStudy::X_vs_generation_lm4();
770   GenLevelStudy::compare_sms_lm4();
771 < GenLevelStudy::MomentumFraction();
771 > */
772 >
773 > ///GenLevelStudy::MomentumFraction();
774 > ///GenLevelStudy::AngleLSPLSP();
775 > ///GenLevelStudy::DeltaLSPmomentum();
776 > ///GenLevelStudy::ZDecayIllustration();
777 > /*
778   GenLevelStudy::AngleMETsumLSP();
779   GenLevelStudy::RatioMETsumLSP();
774 GenLevelStudy::AngleLSPLSP();
780   GenLevelStudy::AngleLSPZ();
781   GenLevelStudy::WidthIllustration();
782 < GenLevelStudy::ZDecayIllustration();
783 < GenLevelStudy::DeltaLSPmomentum();
782 >
783 >
784   GenLevelStudy::pStarIllustration(0.5);
785 < GenLevelStudy::pStarIllustration(0.25);
786 < GenLevelStudy::pStarIllustration(0.75);
782 < GenLevelStudy::pStarDistributions();*/
785 > //GenLevelStudy::pStarIllustration(0.25);
786 > //GenLevelStudy::pStarIllustration(0.75);
787   GenLevelStudy::DrawJetBand();
788 < //GenLevelStudy::DrawOnly100to150inJetBand();
788 > GenLevelStudy::DrawOnly100to150inJetBand();*/
789 > GenLevelStudy::ImpactOfGluinoChi2MassDifference();
790  
791   }
792  
# Line 871 | Line 876 | void jzb_negative_generator_study() {
876  
877   }
878  
874 string ReplaceCharacter(string originalstring,string replacethis,string replacewiththis)
875 {
876        int pos = originalstring.find(replacethis);
877        if(pos == -1) return originalstring;
878        originalstring.replace(pos, replacewiththis.length(), replacewiththis);
879        return originalstring;
880 }
881 string removefunnystring(string name) {
882  name=ReplaceCharacter(name,"[","_");
883  name=ReplaceCharacter(name,"]","_");
884  name=ReplaceCharacter(name,"{","_");
885  name=ReplaceCharacter(name,"}","_");
886  name=ReplaceCharacter(name,".","_");
887  name=ReplaceCharacter(name,",","_");
888  name=ReplaceCharacter(name,";","_");
889  name=ReplaceCharacter(name,":","_");
890  name=ReplaceCharacter(name,"'","_");
891  name=ReplaceCharacter(name,"$","_");
892  name=ReplaceCharacter(name,"@","_");
893  return name;
894 }
895
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);
# Line 952 | Line 935 | void compare_LM4_and_SMS() {
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