#include "TDSuperAna.hh" //void XSroutine(lumtype, muon_data, Zmumu, HwwStn_zewk6m_180_0001.zlog, _zewk6m_runinfo) void TDSuperAna::XSroutine2() { char buffer[500]; //GoodRunv2 lumquery; //TDRunSet lumrunset=(*(_hwwlvlvGoodRun.GetRunSet(_dilType.Int2Type(0))));//-_runstoignore; TDRunSet lumrunset; if(_stn)lumrunset=(*(_hwwlvlvGoodRun13.GetRunSet(_dilType.Int2Type(0))));//-_runstoignore; if(!_stn)lumrunset=(*(_hwwlvlvGoodRun.GetRunSet(_dilType.Int2Type(0))));//-_runstoignore; lumrunset.SetMinMax(_minRun, _maxRun); //lumrunset.SetMax(186598); //lumrunset.SetMin(186599); //double lum_val=lumquery.Loop(&lumrunset)*1.019; TDRunInfoSet lumins("./v13goodrunquery/allreqruns.list_OffLumi_bhel.dat"); double lum_val=(lumins.GetSumTDRunInfo(lumrunset)._goodobsp)*1.019/1000; for(int dilnumber=0; dilnumber!=2; dilnumber++) { if(_stn)lumrunset=(*(_hwwlvlvGoodRun13.GetRunSet(_dilType.Int2Type(dilnumber))))-_runstoignore; if(!_stn)lumrunset=(*(_hwwlvlvGoodRun.GetRunSet(_dilType.Int2Type(dilnumber))))-_runstoignore; lumrunset.SetMinMax(_minRun, _maxRun); // lumrunset.SetMax(186598); // lumrunset.SetMin(186599); // double lum_val=(lumins.GetSumTDRunInfo(*(_hwwlvlvGoodRun.GetRunSet(_dilType.Int2Type(dilnumber))))._goodobsp)/1000; //double lum_val=lumquery.Loop(&lumrunset)*1.019; TDVnEEE Lum(lum_val, 0.0, 0.0, lum_val*0.059, "Lum", "Sample Luminosity"); TDVnEEE ElecmcEventsGenerated(GetSubAna("Zee")->GetGoodGenZ(lumrunset), 0, 0.0, 0.0, "ElecMcGen", "No. Z->ee Generated"); //TD the dilepton factor is taken care of in the weight for now. fitinfo<<(_dilType.Int2String(dilnumber)).c_str()<< " events: "<GetInputEvents(lumrunset)<< " Zs in mass window: "<GetGoodGenMassZ(lumrunset)<< " Zs in mass window and z0<60: "<GetGoodGenZ(lumrunset)<GetGoodGenMassZ(lumrunset)/GetSubAna("Zee")->GetInputEvents(lumrunset)<< " fz0&mass: "<GetGoodGenZ(lumrunset)/GetSubAna("Zee")->GetInputEvents(lumrunset)<< " fz0: "<GetGoodGenZ(lumrunset)/GetSubAna("Zee")->GetGoodGenMassZ(lumrunset)<GetRawDilNumber(dilnumber)<<" MC events failing good run: "<GetFailGoodRun(dilnumber)<< " MC events failing extra good run: "<GetFailExtraGoodRun(dilnumber)<GetGoodZDilepton(dilnumber), "Numdata", "No. Data events"); TDVnEEE dataSSInt(GetSubAna("data")->GetSSDilepton(dilnumber), "NumSSdata", "No. Same Sign Data events"); TDVnEEE ZeeSameSign(GetSubAna("Zee")->GetSSDilepton(dilnumber), "NumSSMC", "No. Same Sign MC events"); TDVnEEE ZeeInt(GetSubAna("Zee")->GetGoodZDilepton(dilnumber), "NumMC", "No. MC events"); // TDVnEEE ElecmcEventsGenerated(GetSubAna("Zee")->GetGoodGenZ(dilnumber), 0, 0.0, 0.0, "ElecMcGen", "No. Z->ee Generated"); sprintf(buffer,"electronZxs_%s", (_dilType.Int2String(dilnumber)).c_str()); TDVnEEE ElectronZxs=crossection(buffer, dataInt, ZeeInt, dataSSInt, Lum, ElecmcEventsGenerated, true, fitinfo); ElectronZxs.PrintValues(); _excelDataFile<<_dilType.Int2String(dilnumber)<<" "; ExcelDataOut(ElectronZxs, _excelDataFile); _excelDataFile<mumu Generated CMUP"); TDVnEEE LumCMX("LumCMX","Luminosity for CMX"); TDVnEEE dataIntCMX("NumdataCMX", "No. Data events CMX"); TDVnEEE dataSSIntCMX("NumSSdataCMX", "No. Same Sign Data events CMX"); TDVnEEE ZmumuSameSignCMX("NumSSMCCMX", "No. Same Sign MC events CMX"); TDVnEEE ZmumuIntCMX("NumMCCMX", "No. MC events CMX"); TDVnEEE MuonmcEventsGeneratedCMX("MuonMcGenCMX", "No. Z->mumu Generated CMX");*/ for(int dilnumber=17; dilnumber!=41; dilnumber++) { if(dilnumber>25 && dilnumber<30)continue; if(dilnumber==36 || dilnumber==37)continue;//no electron-SCMIO xs's. if(_stn)lumrunset=(*(_hwwlvlvGoodRun13.GetRunSet(_dilType.Int2Type(dilnumber))))-_runstoignore; if(!_stn)lumrunset=(*(_hwwlvlvGoodRun.GetRunSet(_dilType.Int2Type(dilnumber))))-_runstoignore; lumrunset.SetMinMax(_minRun, _maxRun); // lumrunset.SetMax(186598); //lumrunset.SetMin(186599); //double lum_val=lumquery.Loop(&lumrunset)*1.019; // double lum_val=(lumins.GetSumTDRunInfo(*(_hwwlvlvGoodRun.GetRunSet(_dilType.Int2Type(dilnumber))))._goodobsp)/1000; TDVnEEE Lum(lum_val, 0.0, 0.0, lum_val*0.059, "Lum", "Sample Luminosity"); TDVnEEE dataInt(GetSubAna("data")->GetGoodZDilepton(dilnumber), "Numdata", "No. Data events"); TDVnEEE dataSSInt(GetSubAna("data")->GetSSDilepton(dilnumber), "NumSSdata", "No. Same Sign Data events"); TDVnEEE ZmumuSameSign(GetSubAna("Zmumu")->GetSSDilepton(dilnumber), "NumSSMC", "No. Same Sign MC events"); TDVnEEE ZmumuInt(GetSubAna("Zmumu")->GetGoodZDilepton(dilnumber), "NumMC", "No. MC events"); TDVnEEE MuonmcEventsGenerated(GetSubAna("Zmumu")->GetGoodGenZ(lumrunset), 0, 0.0, 0.0, "MuonMcGen", "No. Z->mumu Generated"); //if(dilnumber==17){LumCMUP=Lum; MuonmcEventsGeneratedCMUP=MuonmcEventsGenerated;} //if(dilnumber==22){LumCMX=Lum; MuonmcEventsGeneratedCMX=MuonmcEventsGenerated;} fitinfo<<(_dilType.Int2String(dilnumber)).c_str()<< " events: "<GetInputEvents(lumrunset)<< " Zs in mass window: "<GetGoodGenMassZ(lumrunset)<< " Zs in mass window and z0<60: "<GetGoodGenZ(lumrunset)<GetGoodGenMassZ(lumrunset)/GetSubAna("Zmumu")->GetInputEvents(lumrunset)<< " fz0&mass: "<GetGoodGenZ(lumrunset)/GetSubAna("Zmumu")->GetInputEvents(lumrunset)<< " fz0: "<GetGoodGenZ(lumrunset)/GetSubAna("Zmumu")->GetGoodGenMassZ(lumrunset)<GetRawDilNumber(dilnumber)<<" MC events failing good run: "<GetFailGoodRun(dilnumber)<< " MC events failing extra good run: "<GetFailExtraGoodRun(dilnumber)<=17 && dilnumber<=21) { dataIntCMUP+=dataInt; dataSSIntCMUP+=dataSSInt; ZmumuSameSignCMUP+=ZmumuSameSign; ZmumuIntCMUP+=ZmumuInt; }else if(dilnumber>=22 && dilnumber<=25){ dataIntCMX+=dataInt; dataSSIntCMX+=dataSSInt; ZmumuSameSignCMX+=ZmumuSameSign; ZmumuIntCMX+=ZmumuInt; }*/ } /*TDVnEEE CMUPxs=crossection(muonZxs_AllCMUP, dataIntCMUP, ZmumuIntCMUP, dataSSIntCMUP, LumCMUP, MuonmcEventsGeneratedCMUP, true, fitinfo); _excelDataFile<<"ALLCMUP"<<" "; CMUPxs.PrintValues(); ExcelDataOut(CMUPxs, _excelDataFile); _excelDataFile<GetGoodGenZ(lumrunset), 0, 0.0, 0.0, "ElecMcGen", "No. Z->ee Generated"); fitinfo<<(_dilType.Int2String(dilnumber)).c_str()<< " events: "<GetInputEvents(lumrunset)<< " Zs in mass window: "<GetGoodGenMassZ(lumrunset)<< " Zs in mass window and z0<60: "<GetGoodGenZ(lumrunset)<GetGoodGenMassZ(lumrunset)/GetSubAna("Zee")->GetInputEvents(lumrunset)<< " fz0&mass: "<GetGoodGenZ(lumrunset)/GetSubAna("Zee")->GetInputEvents(lumrunset)<< " fz0: "<GetGoodGenZ(lumrunset)/GetSubAna("Zee")->GetGoodGenMassZ(lumrunset)<GetRawDilNumber(dilnumber)<<" MC events failing good run: "<GetFailGoodRun(dilnumber)<< " MC events failing extra good run: "<GetFailExtraGoodRun(dilnumber)<GetGoodZDilepton_mycut(dilnumber), "Numdata", "No. Data events"); TDVnEEE dataSSInt(GetSubAna("data")->GetSSDilepton_mycut(dilnumber), "NumSSdata", "No. Same Sign Data events"); TDVnEEE ZeeSameSign(GetSubAna("Zee")->GetSSDilepton_mycut(dilnumber), "NumSSMC", "No. Same Sign MC events"); TDVnEEE ZeeInt(GetSubAna("Zee")->GetGoodZDilepton_mycut(dilnumber), "NumMC", "No. MC events"); // TDVnEEE ElecmcEventsGenerated(GetSubAna("Zee")->GetGoodGenZ(dilnumber), 0, 0.0, 0.0, "ElecMcGen", "No. Z->ee Generated"); sprintf(buffer,"electronZxs_%s", (_dilType.Int2String(dilnumber)).c_str()); TDVnEEE ElectronZxs=crossection(buffer, dataInt, ZeeInt, dataSSInt, Lum, ElecmcEventsGenerated, true, fitinfo); ElectronZxs.PrintValues(); _excelDataFile<<_dilType.Int2String(dilnumber)<<" "; ExcelDataOut(ElectronZxs, _excelDataFile); _excelDataFile<GetGoodZDilepton_mycut(dilnumber), "Numdata", "No. Data events"); TDVnEEE dataSSInt(GetSubAna("data")->GetSSDilepton_mycut(dilnumber), "NumSSdata", "No. Same Sign Data events"); TDVnEEE ZmumuSameSign(GetSubAna("Zmumu")->GetSSDilepton_mycut(dilnumber), "NumSSMC", "No. Same Sign MC events"); TDVnEEE ZmumuInt(GetSubAna("Zmumu")->GetGoodZDilepton_mycut(dilnumber), "NumMC", "No. MC events"); TDVnEEE MuonmcEventsGenerated(GetSubAna("Zmumu")->GetGoodGenZ(lumrunset), 0, 0.0, 0.0, "MuonMcGen", "No. Z->mumu Generated"); fitinfo<<(_dilType.Int2String(dilnumber)).c_str()<< " events: "<GetInputEvents(lumrunset)<< " Zs in mass window: "<GetGoodGenMassZ(lumrunset)<< " Zs in mass window and z0<60: "<GetGoodGenZ(lumrunset)<GetGoodGenMassZ(lumrunset)/GetSubAna("Zmumu")->GetInputEvents(lumrunset)<< " fz0&mass: "<GetGoodGenZ(lumrunset)/GetSubAna("Zmumu")->GetInputEvents(lumrunset)<< " fz0: "<GetGoodGenZ(lumrunset)/GetSubAna("Zmumu")->GetGoodGenMassZ(lumrunset)<GetRawDilNumber(dilnumber)<<" MC events failing good run: "<GetFailGoodRun(dilnumber)<< " MC events failing extra good run: "<GetFailExtraGoodRun(dilnumber)<GetGoodGenZ(lumrunset), 0, 0.0, 0.0, "ElecMcGen", "No. Z->ee Generated"); fitinfo<<(_dilType.Int2String(dilnumber)).c_str()<< " events: "<GetInputEvents(lumrunset)<< " Zs in mass window: "<GetGoodGenMassZ(lumrunset)<< " Zs in mass window and z0<60: "<GetGoodGenZ(lumrunset)<GetGoodGenMassZ(lumrunset)/GetSubAna("Zee")->GetInputEvents(lumrunset)<< " fz0&mass: "<GetGoodGenZ(lumrunset)/GetSubAna("Zee")->GetInputEvents(lumrunset)<< " fz0: "<GetGoodGenZ(lumrunset)/GetSubAna("Zee")->GetGoodGenMassZ(lumrunset)<GetRawDilNumber(dilnumber)<<" MC events failing good run: "<GetFailGoodRun(dilnumber)<< " MC events failing extra good run: "<GetFailExtraGoodRun(dilnumber)<GetGoodZDilepton_1020cut(dilnumber), "Numdata", "No. Data events"); TDVnEEE dataSSInt(GetSubAna("data")->GetSSDilepton_1020cut(dilnumber), "NumSSdata", "No. Same Sign Data events"); TDVnEEE ZeeSameSign(GetSubAna("Zee")->GetSSDilepton_1020cut(dilnumber), "NumSSMC", "No. Same Sign MC events"); TDVnEEE ZeeInt(GetSubAna("Zee")->GetGoodZDilepton_1020cut(dilnumber), "NumMC", "No. MC events"); // TDVnEEE ElecmcEventsGenerated(GetSubAna("Zee")->GetGoodGenZ(dilnumber), 0, 0.0, 0.0, "ElecMcGen", "No. Z->ee Generated"); sprintf(buffer,"electronZxs_%s", (_dilType.Int2String(dilnumber)).c_str()); TDVnEEE ElectronZxs=crossection(buffer, dataInt, ZeeInt, dataSSInt, Lum, ElecmcEventsGenerated, true, fitinfo); ElectronZxs.PrintValues(); _excelDataFile<<_dilType.Int2String(dilnumber)<<" "; ExcelDataOut(ElectronZxs, _excelDataFile); _excelDataFile<GetGoodZDilepton_1020cut(dilnumber), "Numdata", "No. Data events"); TDVnEEE dataSSInt(GetSubAna("data")->GetSSDilepton_1020cut(dilnumber), "NumSSdata", "No. Same Sign Data events"); TDVnEEE ZmumuSameSign(GetSubAna("Zmumu")->GetSSDilepton_1020cut(dilnumber), "NumSSMC", "No. Same Sign MC events"); TDVnEEE ZmumuInt(GetSubAna("Zmumu")->GetGoodZDilepton_1020cut(dilnumber), "NumMC", "No. MC events"); TDVnEEE MuonmcEventsGenerated(GetSubAna("Zmumu")->GetGoodGenZ(lumrunset), 0, 0.0, 0.0, "MuonMcGen", "No. Z->mumu Generated"); fitinfo<<(_dilType.Int2String(dilnumber)).c_str()<< " events: "<GetInputEvents(lumrunset)<< " Zs in mass window: "<GetGoodGenMassZ(lumrunset)<< " Zs in mass window and z0<60: "<GetGoodGenZ(lumrunset)<GetGoodGenMassZ(lumrunset)/GetSubAna("Zmumu")->GetInputEvents(lumrunset)<< " fz0&mass: "<GetGoodGenZ(lumrunset)/GetSubAna("Zmumu")->GetInputEvents(lumrunset)<< " fz0: "<GetGoodGenZ(lumrunset)/GetSubAna("Zmumu")->GetGoodGenMassZ(lumrunset)<GetRawDilNumber(dilnumber)<<" MC events failing good run: "<GetFailGoodRun(dilnumber)<< " MC events failing extra good run: "<GetFailExtraGoodRun(dilnumber)<GetHistSet("data")->GetAnaHist("SameSignDilMass")->GetHist()->Integral(), "muondataSameSign", "same sign events in the muon data"); TDVnEEE electrondataSameSign(HistSets->GetHistSet("data")->GetAnaHist("SameSignDilMass")->GetHist()->Integral(), "electronsamesign", "same sign events in the electron data"); testmasspeak.RangeMassfitter(*((TH1F*)(HistSets->GetHistSet("data")->GetAnaHist("DileptonMass")->GetHist())), 87.0, 95.0); TDVnEEE muondataInt=testmasspeak.GetVnEIntegral(); double tempsyserror=HistSets->GetHistSet("Zmumu")->GetAnaHist("WeightError")->GetHist()->GetMean(); std::cout<<"mean"<GetHistSet("Zmumu")->GetAnaHist("WeightError")->GetHist()->GetMean()<GetHistSet("Zmumu")->GetAnaHist("Weight")->GetHist()->GetRMS(); //std::cout<<"tempsyserror"<GetHistSet("Zmumu")->GetAnaHist("Weight")->GetHist()->GetMean(); //muondataInt.SetsysterrorOnValue(tempsyserror*muondataInt.Getvalue()); muondataInt.Setname("muondataInt"); testmasspeak.RangeMassfitter(*((TH1F*)(HistSets->GetHistSet("Zmumu")->GetAnaHist("DileptonMass")->GetHist())), 87.0, 95.0); TDVnEEE ZmumuInt=testmasspeak.GetVnEIntegral(); ZmumuInt.Setname("ZmumuInt"); ZmumuInt.SetsysterrorOnValue(tempsyserror*ZmumuInt.Getvalue()); TDRunInfoSet _zewk6m_runinfo("../MCruns/zewk6m_runs_1.log"); _zewk6m_runinfo.AddFile("../MCruns/zewk6m_runs_2.log"); _zewk6m_runinfo.AddFile("../MCruns/zewk6m_runs_3.log"); _zewk6m_runinfo.AddFile("../MCruns/zewk6m_runs_4.log"); _zewk6m_runinfo.AddFile("../MCruns/zewk6m_runs_5.log"); _zewk6m_runinfo.AddFile("../MCruns/zewk6m_runs_6.log"); //#include "zewk6m/zewk6mzlog" /*TDRunInfoSet _zewk6m_runinfo("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0002.zlog"); //_zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0001.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0003.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0004.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0005.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0006.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0007.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0008.zlog"); //_zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0009.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0010.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0011.zlog"); //_zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0012.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0013.zlog"); //_zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_91/HwwStn_zewk6m_180_0014.zlog");*/ /*TDRunInfoSet _zewk6m_runinfo("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_0001.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_0002.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_0003.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_0004.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_0005.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_0006.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_0007.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_0008.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_00013.zlog"); _zewk6m_runinfo.AddFile("/data/nglas08/b/stdenis/zewk6m/v1_87/HwwStn_zewk6m_180_00014.zlog");*/ //std::cout<<"_zewk6m_runinfo.PrintValues();"<mumu Generated");//MCruns //TDVnEEE MuonmcEventsGenerated(1595781, 0, 0.0, 0.0, "MuonMcGen", "No. Z->mumu Generated");//MC char buffer[100]; //printf(buffer,"%s_%d", _dilset.GetDilString().c_str(), _weighting); //_excelMassFile<GetHistSet("data")->GetAnaHist("DileptonMass")->GetHist())), 87.0, 95.0); tempsyserror=HistSets->GetHistSet("Zee")->GetAnaHist("WeightError")->GetHist()->GetMean(); //std::cout<<"tempsyserror"<GetHistSet("Zee")->GetAnaHist("Weight")->GetHist()->GetRMS(); std::cout<<"RMS"<GetHistSet("Zee")->GetAnaHist("Weight")->GetHist()->GetRMS()<GetHistSet("Zee")->GetAnaHist("Weight")->GetHist()->GetMean(); //std::cout<<"tempsyserror"<GetHistSet("Zee")->GetAnaHist("DileptonMass")->GetHist())), 87.0, 95.0); TDVnEEE ZeeInt=testmasspeak.GetVnEIntegral(); ZeeInt.Setname("ZeeInt"); ZeeInt.SetsysterrorOnValue(tempsyserror*ZeeInt.Getvalue()); TDRunInfoSet _zewk6d_runinfo("../MCruns/zewk6d_runs_1.log"); _zewk6d_runinfo.AddFile("../MCruns/zewk6d_runs_2.log"); _zewk6d_runinfo.AddFile("../MCruns/zewk6d_runs_3.log"); _zewk6d_runinfo.AddFile("../MCruns/zewk6d_runs_4.log"); _zewk6d_runinfo.AddFile("../MCruns/zewk6d_runs_5.log"); _zewk6d_runinfo.AddFile("../MCruns/zewk6d_runs_6.log"); //#include "zewk6d/zewk6dzlog" // #include "zewkae/zewkaezlog" /*TDRunInfoSet _zewk6d_runinfo("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0001.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0002.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0003.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0004.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0005.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0006.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0007.zlog"); //_zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0008.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0009.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0010.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0011.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0012.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0013.zlog"); _zewk6d_runinfo.AddFile("/data/nglas08/b/stdenis/zewkae/v1_91/HwwStn_zewkae_180_0014.zlog");*/ /*TDRunInfoSet _zewk6d_runinfo("zewk6ae/HwwStn_zewkae_180_0001.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0002.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0003.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0004.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0005.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0006.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0007.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0008.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0009.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0011.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0012.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0013.zlog"); _zewk6d_runinfo.AddFile("zewk6ae/HwwStn_zewkae_180_0014.zlog");*/ goodobsvZ=(int)(_zewk6d_runinfo.GetSumTDRunInfo(MCGoodRuns)._goodobsp); TDVnEEE ElecmcEventsGenerated(goodobsvZ, 0, 0.0, 0.0, "ElecMcGen", "No. Z->ee Generated");//MCbadrun //TDVnEEE ElecmcEventsGenerated(1540477, 0, 0.0, 0.0, "ElecMcGen", "No. Z->ee Generated");//MC //char buffer[100]; sprintf(buffer,"%s_%d", _dilset.GetDilString().c_str(), _weighting); _excelDataFile<GetHistSet("Zmumu")->GetAnaHist("Weight")->GetHist()->GetMean(), subruns); thisxs.PrintValues(); } //excelplot.AddData(exceldata); //excelplot.PrintValues(); //excelplot.AddData(exceldata); //excelplot.AddData(exceldata); //excelplot.AddData(exceldata); //excelplot.AddData(exceldata); //excelplot.plot(CanvMan1); _excelDataFile<GetHistSet(dataname_s.c_str())->GetAnaHist("SameSignDilMass")->GetHist()->Integral(), (dataname_s+"dataSameSign").c_str(), ("same sign events in the "+dataname_s).c_str()); TDVnEEE dataInt(dataintegral); double tempsyserror=HistSets->GetHistSet(mc_name_s.c_str())->GetAnaHist("WeightError")->GetHist()->GetMean(); if(tempsyserror==0.0)tempsyserror=HistSets->GetHistSet(mc_name_s.c_str())->GetAnaHist("Weight")->GetHist()->GetRMS(); double meanweight=HistSets->GetHistSet(mc_name_s.c_str())->GetAnaHist("Weight")->GetHist()->GetMean(); dataInt.Setname((dataname_s+"Int").c_str()); testmasspeak.RangeMassfitter(*(HistSets->GetHistSet(mc_name_s.c_str())->GetAnaHist("DileptonMass")->GetHist()), 87.0, 95.0); TDVnEEE mcInt=testmasspeak.GetVnEIntegral(); mcInt.Setname((mc_name_s+"Int").c_str()); mcInt.SetsysterrorOnValue(tempsyserror*mcInt.Getvalue()); TDVnEEE Lum(_dilWeight.GetLum(_dilset.GetRepDilType()), 0.0, 0.0, _dilWeight.GetLumErr(_dilset.GetRepDilType()), "Lum", "Sample Luminosity"); TDVnEEE mcEventsGenerated(goodobsvZ, 0, 0.0, 0.0, (mc_name_s+"Gen").c_str(), ("No. Z Generated in"+mc_name_s).c_str()); TDVnEEE Zxs=crossection((dataname_s+"Zxs").c_str(), dataInt, mcInt, dataSameSign, Lum, mcEventsGenerated, true, fitinfo); Zxs.PrintValues(); if(Zxs.Getvalue()!=0.0) { Zxs.PrintValues(os1); ExcelXS(_dilset.GetDilString().c_str(), dataInt, mcInt, dataSameSign, Lum, mcEventsGenerated, meanweight, os2); } }*/ std::ostream& TDSuperAna::PrintAcceptances(std::ostream& os) { TDHistSetSet thisset(CanvMan1, false); thisset.AddHistSet(HistSets->GetHistSet("WW")); thisset.AddHistSet(HistSets->GetHistSet("WZ")); thisset.AddHistSet(HistSets->GetHistSet("ZZ")); thisset.AddHistSet(HistSets->GetHistSet("ttbar")); thisset.AddHistSet(HistSets->GetHistSet("wgamma")); TDHistSet DY(*(HistSets->GetHistSet("Zee"))); DY+=(*(HistSets->GetHistSet("Zmumu"))); DY+=(*(HistSets->GetHistSet("Ztautau"))); DY.SetDataSourcePrefix("DYll"); thisset.AddHistSet(&DY); TDHistSet BackGround(*(HistSets->GetHistSet("WW"))); BackGround+=(*(HistSets->GetHistSet("WZ"))); BackGround+=(*(HistSets->GetHistSet("ZZ"))); BackGround+=(*(HistSets->GetHistSet("ttbar"))); BackGround+=DY; BackGround+=(*(HistSets->GetHistSet("wgamma"))); BackGround.SetDataSourcePrefix("Background Total"); thisset.AddHistSet(&BackGround); thisset.AddHistSet(HistSets->GetHistSet("higgsww")); thisset.AddHistSet(HistSets->GetHistSet("data")); os<<" \\begin{table}[H]"<GetDataSource(),"Background Total"))) os<<"\\hline"<GetDataSource(),"data"))) os<<"\\hline"<GetDataSource(),"higgsww"))) os<<"\\hline"<GetDataSource()<<" & " <<(*((*(thisset.GetIt()))->GetAnaIt("dilDphiee")))->GetHist()->Integral(); if(strcmp((*(thisset.GetIt()))->GetDataSource(),"data")) os<<"+-"<GetAnaIt("dilDphiee")))->GetHist()->Integral()); os<<" & "; //<GetAnaIt("dilDphiemu")))->GetHist()->Integral(); if(strcmp((*(thisset.GetIt()))->GetDataSource(),"data")) os<<"+-"<GetAnaIt("dilDphiemu")))->GetHist()->Integral()); os<<" & "; os<<(*((*(thisset.GetIt()))->GetAnaIt("dilDphimumu")))->GetHist()->Integral(); if(strcmp((*(thisset.GetIt()))->GetDataSource(),"data")) os<<"+-"<GetAnaIt("dilDphimumu")))->GetHist()->Integral()); os<<" & "; os<<(*((*(thisset.GetIt()))->GetAnaIt("dilDphi")))->GetHist()->Integral(); if(strcmp((*(thisset.GetIt()))->GetDataSource(),"data")) os<<"+-"<GetAnaIt("dilDphi")))->GetHist()->Integral()); os<<" \\\\ "< TDSuperAna::getDatasets() { TDHistSetSet thisset(CanvMan1, false); thisset.AddHistSet(HistSets->GetHistSet("WW")); thisset.AddHistSet(HistSets->GetHistSet("WZ")); thisset.AddHistSet(HistSets->GetHistSet("ZZ")); thisset.AddHistSet(HistSets->GetHistSet("ttbar")); thisset.AddHistSet(HistSets->GetHistSet("wgamma")); TDHistSet DY(*(HistSets->GetHistSet("Zee"))); DY+=(*(HistSets->GetHistSet("Zmumu"))); DY+=(*(HistSets->GetHistSet("Ztautau"))); DY.SetDataSourcePrefix("DYll"); thisset.AddHistSet(&DY); TDHistSet BackGround(*(HistSets->GetHistSet("WW"))); BackGround+=(*(HistSets->GetHistSet("WZ"))); BackGround+=(*(HistSets->GetHistSet("ZZ"))); BackGround+=(*(HistSets->GetHistSet("ttbar"))); BackGround+=DY; BackGround+=(*(HistSets->GetHistSet("wgamma"))); BackGround.SetDataSourcePrefix("Background Total"); thisset.AddHistSet(&BackGround); thisset.AddHistSet(HistSets->GetHistSet("higgsww")); thisset.AddHistSet(HistSets->GetHistSet("data")); std::vector datsetvec; datsetvec.push_back("$M_H$"); if(!(thisset.InitiateIt()))do{ datsetvec.push_back((*(thisset.GetIt()))->GetDataSource()); }while(!(thisset.IncrementIt())); return datsetvec; } std::vector TDSuperAna::getAcceptanceVector() { TDHistSetSet thisset(CanvMan1, false); thisset.AddHistSet(HistSets->GetHistSet("WW")); thisset.AddHistSet(HistSets->GetHistSet("WZ")); thisset.AddHistSet(HistSets->GetHistSet("ZZ")); thisset.AddHistSet(HistSets->GetHistSet("ttbar")); thisset.AddHistSet(HistSets->GetHistSet("wgamma")); TDHistSet DY(*(HistSets->GetHistSet("Zee"))); DY+=(*(HistSets->GetHistSet("Zmumu"))); DY+=(*(HistSets->GetHistSet("Ztautau"))); DY.SetDataSourcePrefix("DYll"); thisset.AddHistSet(&DY); TDHistSet BackGround(*(HistSets->GetHistSet("WW"))); BackGround+=(*(HistSets->GetHistSet("WZ"))); BackGround+=(*(HistSets->GetHistSet("ZZ"))); BackGround+=(*(HistSets->GetHistSet("ttbar"))); BackGround+=DY; BackGround+=(*(HistSets->GetHistSet("wgamma"))); BackGround.SetDataSourcePrefix("Background Total"); thisset.AddHistSet(&BackGround); thisset.AddHistSet(HistSets->GetHistSet("higgsww")); thisset.AddHistSet(HistSets->GetHistSet("data")); char buffer[500]; sprintf(buffer,"%d",(int)_higgsMass); std::vector accvec; accvec.push_back(buffer); if(!(thisset.InitiateIt()))do{ //accvec.push_back((*((*(thisset.GetIt()))->GetAnaIt("dilDphi")))->GetHist()->Integral()); if(!strcmp((*(thisset.GetIt()))->GetDataSource(),"data")) { sprintf(buffer,"%d",(int)(*((*(thisset.GetIt()))->GetAnaIt("dilDphi")))->GetHist()->Integral()); }else{ sprintf(buffer,"%2.3f+-%2.3f",(*((*(thisset.GetIt()))->GetAnaIt("dilDphi")))->GetHist()->Integral(), sqrt((*((*(thisset.GetIt()))->GetAnaIt("dilDphi")))->GetHist()->Integral())); } accvec.push_back(buffer); }while(!(thisset.IncrementIt())); return accvec; } //void TDSuperAna::CreateSecondarythisset() //{ // _Secondarythisset=new (TDHistSetSet)(CanvMan1, false); // } void TDSuperAna::likelihood(bool sim) { /* const char *datset[10] = {"hww","wtop1w","wgamma","ttopel","ztopcz","wtop1z", "fakeBg","zexoll","dataSg","dataSg"};*/ //TDRunSet lumrunset; //if(_stn)lumrunset=(*(_hwwlvlvGoodRun13.GetRunSet(_dilType.Int2Type(0))));//-_runstoignore; //if(!_stn)lumrunset=(*(_hwwlvlvGoodRun.GetRunSet(_dilType.Int2Type(0)))); //lumrunset.SetMinMax(_minRun, _maxRun); //TDRunInfoSet lumins("./v13goodrunquery/allreqruns.list_OffLumi_bhel.dat"); // double Lum=(lumins.GetSumTDRunInfo(lumrunset)._goodobsp)*1.019/1000; double Lum=GetSubAna("data")->GetAvLum(); // Lum = 360.; double sLum = Lum*0.059; // const double sLum = 22.; // BR(WW->llvv) = 0.1068 in PYTHIA vs 0.11 in HERWIG /*const double nEvt[10] = {1069855.0, 1063469.0, 1127915.0, 1074251.0, 1089119.0, 1051020.0, 1033214.0, 1079780.0, 1092527.0, 1093662.0}; //int main() //void pseudolimits() { int mHig=170; char *hmas = gSystem->Getenv("HiggsMass"); if (hmas) for (int i=0; iGetenv("LimiLevel"); if (lvll) for (int i=0; i0) sprintf(filename,"results/%s_%d_phxCorr.root",datset[i],mHig); else sprintf(filename,"results/%s%d_%d_phxCorr.root",datset[i],mHig,mHig); f[i] = new TFile(filename,"READ"); } sprintf(filename,"dildPhi_all_%d",lLvl); */ //TH2F* g = (TH2F*)f[0]->Get("GC")->Clone(); //TD what is this? double nAcc=0.0; //for (int j=0; j<26; j++) nAcc += g->GetBinContent(lLvl+4,j+1); //TD what is this? nAcc=GetSubAna("higgsww")->GetAnaCutSet()->GetCutAnaCut()->Getpassed(); // double Hacc = 0.1068 * nAcc / nEvt[((int)_higgsMass-110)/10]; //TD where do the n Evt come from? int events=(int)(GetSubAna("higgsww")->GetInputEvents()); double xs=0.001; //TD nominal value if(events==0) { std::cout<<"using hard coded number of events"<Get(filename)->Clone(); h[i]->Rebin(20); if (i>2 && i<8) h[2]->Add(h[i]); }*/ // std::vector TDHistSet SigHS(*(HistSets->GetHistSet("higgsww"))); TDHistSet BkgHS(*(HistSets->GetHistSet("Zee"))); BkgHS+=(*(HistSets->GetHistSet("Zmumu"))); BkgHS+=(*(HistSets->GetHistSet("Ztautau"))); BkgHS+=(*(HistSets->GetHistSet("wgamma"))); BkgHS+=(*(HistSets->GetHistSet("ttbar"))); //BkgHS+=(*(HistSets->GetHistSet("WW"))); BkgHS+=(*(HistSets->GetHistSet("WZ"))); BkgHS+=(*(HistSets->GetHistSet("ZZ"))); TDHistSet WWHS(*(HistSets->GetHistSet("WW"))); TDHistSet DATAHS(*(HistSets->GetHistSet("data"))); SigHS.ScaleErrors(SigHS.GetAnaHist("RelWeightError")->GetHist()->GetMean()); BkgHS.ScaleErrors(BkgHS.GetAnaHist("RelWeightError")->GetHist()->GetMean()); WWHS.ScaleErrors(WWHS.GetAnaHist("RelWeightError")->GetHist()->GetMean()); DATAHS.ScaleErrors(DATAHS.GetAnaHist("RelWeightError")->GetHist()->GetMean()); /*TH1F* h[4]; h[0]=(TH1F*)HistSets->GetHistSet("higgsww")->GetAnaHist("dilDphilike")->GetHist()->Clone(); h[0]->Rebin(20); h[1]=(TH1F*)HistSets->GetHistSet("WW")->GetAnaHist("dilDphilike")->GetHist()->Clone(); h[1]->Rebin(20); h[2]=(TH1F*)HistSets->GetHistSet("Zmumu")->GetAnaHist("dilDphilike")->GetHist()->Clone(); h[2]->Add((TH1F*)HistSets->GetHistSet("Zee")->GetAnaHist("dilDphilike")->GetHist()->Clone()); h[2]->Add((TH1F*)HistSets->GetHistSet("Ztautau")->GetAnaHist("dilDphilike")->GetHist()->Clone()); h[2]->Add((TH1F*)HistSets->GetHistSet("WZ")->GetAnaHist("dilDphilike")->GetHist()->Clone()); h[2]->Add((TH1F*)HistSets->GetHistSet("ZZ")->GetAnaHist("dilDphilike")->GetHist()->Clone()); h[2]->Add((TH1F*)HistSets->GetHistSet("wgamma")->GetAnaHist("dilDphilike")->GetHist()->Clone()); h[2]->Add((TH1F*)HistSets->GetHistSet("ttbar")->GetAnaHist("dilDphilike")->GetHist()->Clone()); h[2]->Rebin(20); h[3]=(TH1F*)HistSets->GetHistSet("data")->GetAnaHist("dilDphilike")->GetHist()->Clone(); h[3]->Rebin(20);*/ TH1F* h[4]; h[0]=(TH1F*)SigHS.GetAnaHist("dilDphilike")->GetHist()->Clone(); h[0]->Rebin(20); h[1]=(TH1F*)WWHS.GetAnaHist("dilDphilike")->GetHist()->Clone(); h[1]->Rebin(20); h[2]=(TH1F*)BkgHS.GetAnaHist("dilDphilike")->GetHist()->Clone(); h[2]->Rebin(20); h[3]=(TH1F*)DATAHS.GetAnaHist("dilDphilike")->GetHist()->Clone(); h[3]->Rebin(20); /*if(!(HistSets->InitiateIt()))do{ if(h[0]==NULL && std::string((*(HistSets->GetIt()))->GetDataSource())==std::string("higgsww")) h[0]=(TH1F*)(*(HistSets->GetIt()))->GetAnaHist("dilDphilike")->GetHist()->Clone(); else if(h[1]==NULL && std::string((*(HistSets->GetIt()))->GetDataSource())==std::string("WW")) h[1]=(TH1F*)(*(HistSets->GetIt()))->GetAnaHist("dilDphilike")->GetHist()->Clone(); else if(h[2]==NULL) h[2]=(TH1F*)(*(HistSets->GetIt()))->GetAnaHist("dilDphilike")->GetHist()->Clone(); else h[2]->Add((TH1F*)(*(HistSets->GetIt()))->GetAnaHist("dilDphilike")->GetHist()->Clone()); }while(!HistSets->IncrementIt());*/ int nBins = h[0]->GetNbinsX(); //int nBins = HistSets->GetHistSet("higgsww")->GetAnaHist("dilDphilike")->GetHist()->GetNbinsX(); double nWW = h[1]->Integral(); double nBG = h[2]->Integral(); double sWW=0.0, sBG=0.0, norm_hww=0.0, norm_ww=0.0, norm_bg=0.0; double data_dil[200], hww_all[200], ww_sm[200], bg_sm[200]; for (int b=0; bGetBinContent(b+1); ww_sm[b] = h[1]->GetBinContent(b+1); bg_sm[b] = h[2]->GetBinContent(b+1); norm_hww += hww_all[b]; norm_ww += ww_sm[b]; norm_bg += bg_sm[b]; } sWW = 0.09 * nWW; sBG = 0.17 * nBG; cout << " Higgs Mass = " << (int)_higgsMass << " GeV; Cut Level = " << lLvl << endl; //TD someone has included a 'using std' some where cout << " Luminosity = " << Lum << "+-" << sLum << endl; cout << " HWW A*BR = " << Hacc << "+-" << sHacc << endl; cout << " WW BG = " << nWW << "+-" << sWW << endl; cout << " otherBGs = " << nBG << "+-" << sBG << endl << endl; fitinfo << " Higgs Mass = " << (int)_higgsMass << " GeV; Cut Level = " << lLvl << std::endl; fitinfo << " Luminosity = " << Lum << "+-" << sLum << std::endl; fitinfo << " HWW A*BR = " << Hacc << "+-" << sHacc << std::endl; fitinfo << " WW BG = " << nWW << "+-" << sWW << std::endl; fitinfo << " otherBGs = " << nBG << "+-" << sBG << std::endl << std::endl; norm_hww = 1./ norm_hww; norm_ww = 1./ norm_ww; norm_bg = 1./ norm_bg; char buffer[200]; char buffer2[200]; sprintf(filename,"results/psLm_phxCorr_%3d_%d.root",(int)_higgsMass,lLvl); TFile* v = new TFile(filename,"RECREATE"); if(sim)sprintf(buffer,"simcrosseclim"); else sprintf(buffer,"crosseclim"); if(sim)sprintf(buffer2,"simulated crossection limit at 95%s CL", "%"); else sprintf(buffer2,"crossection limit at 95%s CL","%"); TH1F* xsec = new TH1F(buffer,buffer2,1000,0.0,50.0); //TD memory leak! TH1F* xlik[1000]; TRandom* r = new TRandom(); int npseudo=1; if(sim)npseudo=1000; for (int mn=0; mnGaus(nWW, sWW); int nWW_poisson = r->Poisson(nWW_gauss); for (int j=0; jGetRandom(); data_dil[int(nBins*dphi/3.1416)]++; } double nBG_gauss = r->Gaus(nBG, sBG); int nBG_poisson = r->Poisson(nBG_gauss); for (int j=0; jGetRandom(); data_dil[int(nBins*dphi/3.1416)]++; } }else{ for (int b=0; bGetBinContent(b+1); } double par[5], lik[1000], intgratlik=0.0, intpartlik=0.0, xlim; for (int j1=0; j1<1000; j1++) { lik[j1]=0.0; for (int it=0; it<100; it++) { par[0] = 0.05*j1; par[1] = r->Gaus(Lum, sLum); par[2] = r->Gaus(Hacc,sHacc); par[3] = r->Gaus(nWW, sWW); par[4] = r->Gaus(nBG, sBG); double mu=0.0, loglik=0.0; for (int b=0; b0) loglik += (data_dil[b]*TMath::Log(mu) - mu - TMath::Log(TMath::Factorial((int)data_dil[b]))); } lik[j1] += TMath::Exp(loglik); } lik[j1] /= 100.0; intgratlik += lik[j1]; xlik[mn]->Fill(par[0], lik[j1]); } int j1=0; while ((intpartlik/intgratlik)<0.95) { intpartlik += lik[j1]; xlim = 0.05*j1; j1++; } xsec->Fill(xlim); //cout << " xlim[" << mn << "] = " << xlim << " pb" << endl; // fitinfo << " xlim[" << mn << "] = " << xlim << " pb" << endl; if (mn>0 && fmod((float)mn,(float)100)==0) { double xsec_median[1], x[1]; x[0]=0.5; xsec->GetQuantiles(1,xsec_median,x); // cout << " Crossec Lim Med = " << xsec_median[0] << " pb" << endl << endl; } } cout << " Cross-section Limit Mean = " << xsec->GetMean() << " pb" << endl; fitinfo << " Cross-section Limit Mean = " << xsec->GetMean() << " pb" << endl; double xsec_median[1], x[1]; x[0]=0.5; xsec->GetQuantiles(1,xsec_median,x); cout << " Crossection Limit Median = " << xsec_median[0] << " pb" << endl; cout << " Crossection Limit Median as ratio = " << xsec_median[0]/(GetHwwXS((int)_higgsMass, _stn)*0.1068) << endl << endl; if(sim) { fitinfo << "expected:"<Draw(); xsec->Write(); if(!sim)xlik[0]->Write(); //v->Write(); v->Close(); //for (int i=0; i<10; i++) f[i]->Close(); } void TDSuperAna::TomLikelihood7708() { std::cout<<"TomLikelihood7708"<GetHistSet("Zee"))); TomLikeBG+=(*(HistSets->GetHistSet("Zmumu"))); TomLikeBG+=(*(HistSets->GetHistSet("Ztautau"))); TomLikeBG+=(*(HistSets->GetHistSet("wgamma"))); TomLikeBG+=(*(HistSets->GetHistSet("ttbar"))); // TomLikeBG+=(*(HistSets->GetHistSet("WW"))); TomLikeBG+=(*(HistSets->GetHistSet("WZ"))); TomLikeBG+=(*(HistSets->GetHistSet("ZZ"))); TDHistSet TomLikeWW(*(HistSets->GetHistSet("WW"))); TH1F * WWbkgr=(TH1F*)(TomLikeWW.GetAnaHist("dilDphi7708")->GetHist()); TH1F * bkgr=(TH1F*)(TomLikeBG.GetAnaHist("dilDphi7708")->GetHist()); TDHistSet TomLikeSig(*(HistSets->GetHistSet("higgsww"))); TH1F * sig=(TH1F*)(TomLikeSig.GetAnaHist("dilDphi7708")->GetHist()); TH1F * data=(TH1F*)((HistSets->GetHistSet("data"))->GetAnaHist("dilDphi7708")->GetHist()); nps_low[1] = -0.08; nps_high[1] = 0.08; nullhyp->add_template(bkgr,1.0,2,ename,nps_low,nps_high, lowshape,lowsigma,highshape,highsigma,0,0,"Channel 1"); //TD one np for now. testhyp->add_template(bkgr,1.0,2,ename,nps_low,nps_high, lowshape,lowsigma,highshape,highsigma,0,0,"Channel 1"); nullhyp->add_template(WWbkgr,1.0,2,ename,WWnps_low,WWnps_high, lowshape,lowsigma,highshape,highsigma,0,0,"Channel 1"); testhyp->add_template(WWbkgr,1.0,2,ename,WWnps_low,WWnps_high, lowshape,lowsigma,highshape,highsigma,0,0,"Channel 1"); testhyp->add_template(sig,1.0,2,ename,signps_low,signps_high, lowshape,lowsigma,highshape,highsigma,0,1,"Channel 1"); csm_model* nullhyp_pe = nullhyp->Clone(); csm_model* testhyp_pe = testhyp->Clone(); TCanvas * mycanvas = (TCanvas *) new TCanvas("Canvas1","Canvas1"); mycanvas->Divide(2,1); mycanvas->cd(1); testhyp->plotwithdata("Channel 1",data ); csm* mycsm = new csm(); mycsm->set_htofit( data , "Channel 1"); mycsm->set_modeltofit(testhyp); double chisq = mycsm->chisquared(); std::cout<<"chisq: "<getbestmodel(); mycanvas->cd(2); bestsignalfit->plotwithdata("Channel 1", data); mycanvas->Print("results/tomlikeout.ps"); mclimit_csm* mymclimit = new mclimit_csm(); mymclimit->set_null_hypothesis(nullhyp); mymclimit->set_test_hypothesis(testhyp); mymclimit->set_null_hypothesis_pe(nullhyp_pe); mymclimit->set_test_hypothesis_pe(testhyp_pe); mymclimit->set_datahist(data, "Channel 1"); mymclimit->set_npe(1000); mymclimit->run_pseudoexperiments(); std::cout << "Getting results" << std::endl; std::cout << "ts: " << mymclimit->ts() << std::endl; std::cout << "tsbm2: " << mymclimit->tsbm2() << std::endl; std::cout << "tsbm1: " << mymclimit->tsbm1() << std::endl; std::cout << "tsbmed: " << mymclimit->tsbmed() << std::endl; std::cout << "tsbp1: " << mymclimit->tsbp1() << std::endl; std::cout << "tsbp2: " << mymclimit->tsbp2() << std::endl; std::cout << "tssm2: " << mymclimit->tssm2() << std::endl; std::cout << "tssm1: " << mymclimit->tssm1() << std::endl; std::cout << "tssmed: " << mymclimit->tssmed() << std::endl; std::cout << "tssp1: " << mymclimit->tssp1() << std::endl; std::cout << "tssp2: " << mymclimit->tssp2() << std::endl; std::cout << "CLs: " << mymclimit->cls() << std::endl; std::cout << "CLb: " << mymclimit->clb() << std::endl; std::cout << "CLsb: " << mymclimit->clsb() << std::endl; double conf95=mymclimit->s95med(); std::cout << "Median expected s95 in bg hypothesis: " << conf95 << std::endl; fitinfo << "Median expected s95 in bg hypothesis: " << conf95 << std::endl; std::cout<<"deleting Tomlike elements"<sysname.push_back(sname); ername[0] = "LUMI"; serr[0] = 0.04; sderr[0] = -0.04; wwerr[0] = 0.04; wwderr[0] = -0.04; bgerr[0] = 0.04; bgderr[0] = -0.04; // uncorrelated part of lumi error //sname = new char[strlen("CDFLUMI")+1]; //strcpy(sname,"CDFLUMI"); //chanresult->sysname.push_back(sname); ername[1] = "CDFLUMI"; serr[1] = 0.04; sderr[1] = -0.04; wwerr[1] = 0.04; wwderr[1] = -0.04; bgerr[1] = 0.04; bgderr[1] = -0.04; //sname = new char[strlen("CDFLEPTRIGGER")+1]; //strcpy(sname,"CDFLEPTRIGGER"); //chanresult->sysname.push_back(sname); ername[2] = "CDFLEPTRIGGER"; serr[2] = 0.01; sderr[2] = -0.01; wwerr[2] = 0.01; wwderr[2] = -0.01; bgerr[2] = 0.01; // only partly estimated in CDF 7708 bgderr[2] = -0.01; // only partly estimated in CDF 7708 //sname = new char[strlen("ISR")+1]; //strcpy(sname,"ISR"); //chanresult->sysname.push_back(sname); ername[3] = "ISR"; serr[3] = 0.03; sderr[3] = -0.03; wwerr[3] = 0.03; // not estimated in CDF 7708 wwderr[3] = -0.03; // not estimated in CDF 7708 bgerr[3] = 0.02; // partly estimated in CDF 7708 bgderr[3] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("PDF")+1]; //strcpy(sname,"PDF"); //chanresult->sysname.push_back(sname); ername[4] = "PDF"; serr[4] = 0.03; sderr[4] = -0.03; wwerr[4] = 0.06; wwderr[4] = -0.06; bgerr[4] = 0.04; bgderr[4] = -0.04; //sname = new char[strlen("CDFJES")+1]; //strcpy(sname,"CDFJES"); //chanresult->sysname.push_back(sname); ername[5] = "CDFJES"; serr[5] = 0.01; sderr[5] = -0.01; wwerr[5] = 0.01; wwderr[5] = -0.01; bgerr[5] = 0.02; // mostly from fakes bgderr[5] = -0.02; // mostly from fakes //sname = new char[strlen("CDFTRKISO")+1]; //strcpy(sname,"CDFTRKISO"); //chanresult->sysname.push_back(sname); ername[6] = "CDFTRKISO"; serr[6] = 0.02; sderr[6] = -0.02; wwerr[6] = 0.02; wwderr[6] = -0.02; bgerr[6] = 0.02; // partly estimated in CDF 7708 bgderr[6] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("ALPHAS")+1]; //strcpy(sname,"ALPHAS"); //chanresult->sysname.push_back(sname); ername[7] = "ALPHAS"; serr[7] = 0.03; sderr[7] = -0.03; wwerr[7] = 0.0; wwderr[7] = 0.0; bgerr[7] = 0.0; bgderr[7] = 0.0; //sname = new char[strlen("CDFLEPTID")+1]; //strcpy(sname,"CDFLEPTID"); //chanresult->sysname.push_back(sname); /*ername[8] = "CDFLEPTID"; serr[8] = 0.02; sderr[8] = -0.02; wwerr[8] = 0.02; wwderr[8] = -0.02; bgerr[8] = 0.02; bgderr[8] = -0.02;*/ //TD i think these are the the SF errors which I am including later, but I'm keeping these in as 0 in case I need to include them again ername[8] = "CDFLEPTID"; serr[8] = 0.00; sderr[8] = -0.00; wwerr[8] = 0.00; wwderr[8] = -0.00; bgerr[8] = 0.00; bgderr[8] = -0.00; //sname = new char[strlen("CDFFAKERATE")+1]; //strcpy(sname,"CDFFAKERATE"); //chanresult->sysname.push_back(sname); ername[9] = "CDFFAKERATE"; serr[9] = 0.0; sderr[9] = 0.0; wwerr[9] = 0.0; wwderr[9] = 0.0; bgerr[9] = 0.05; // 50% error on fakes, which are 10% of background bgderr[9] = -0.05; // 50% error on fakes, which are 10% of background //sname = new char[strlen("NNLO")+1]; //strcpy(sname,"NNLO"); //chanresult->sysname.push_back(sname); ername[10] = "NNLO"; serr[10] = 0.1; sderr[10] = -0.1; wwerr[10] = 0.0; wwderr[10] = 0.0; bgerr[10] = 0.0; bgderr[10] = 0.0; //sname = new char[strlen("CDFJETVETO")+1]; //strcpy(sname,"CDFJETVETO"); //chanresult->sysname.push_back(sname); ername[11] = "CDFJETVETO"; serr[11] = 0.0; sderr[11] = 0.0; wwerr[11] = 0.06; wwderr[11] = -0.06; bgerr[11] = 0.0; bgderr[11] = 0.0; //sname = new char[strlen("GENERATOR")+1]; //strcpy(sname,"GENERATOR"); //chanresult->sysname.push_back(sname); ername[12] = "GENERATOR"; serr[12] = 0.0; sderr[12] = 0.0; wwerr[12] = 0.04; wwderr[12] = -0.04; bgerr[12] = 0.06; bgderr[12] = -0.06; //sname = new char[strlen("CDFDYMET")+1]; //strcpy(sname,"CDFDYMET"); //chanresult->sysname.push_back(sname); ername[13] = "CDFDYMET"; serr[13] = 0.0; sderr[13] = 0.0; wwerr[13] = 0.0; wwderr[13] = 0.0; bgerr[13] = 0.05; bgderr[13] = -0.05; ername[14] = "SFERRORS"; TH1 * sigupshape[14]; TH1 * sigdownshape[14]; double sigupshapesigma[14]; double sigdownshapesigma[14]; TH1 * wwupshape[14]; TH1 * wwdownshape[14]; double wwupshapesigma[14]; double wwdownshapesigma[14]; TH1 * bgupshape[14]; TH1 * bgdownshape[14]; double bgupshapesigma[14]; double bgdownshapesigma[14]; for (int i=0;i<14;i++) { sigupshape[i] = 0; sigdownshape[i] = 0; sigupshapesigma[i] = 0; sigdownshapesigma[i] = 0; wwupshape[i] = 0; wwdownshape[i] = 0; wwupshapesigma[i] = 0; wwdownshapesigma[i] = 0; bgupshape[i] = 0; bgdownshape[i] = 0; bgupshapesigma[i] = 0; bgdownshapesigma[i] = 0; } csm_model* nullhyp = new csm_model(); csm_model* testhyp = new csm_model(); /* char *ename[5]; char en0[]="lumerror"; char en1[]="error1"; char en2[]="error2"; char en3[]="error3"; char en4[]="error4"; double nps_low[5]; double nps_high[5]; nps_low[0] = -0.059; nps_high[0] = 0.059; nps_low[1] = -0.08;//TD ? used later for weight errors. nps_high[1] = 0.08;//TD ? nps_low[2] = -0.05;//TD ? nps_high[2] = 0.05;//TD ? nps_low[3] = 0.08;//TD ? nps_high[3] = -0.08;//TD ? nps_low[4] = 0.05;//TD ? nps_high[4] = -0.05;//TD ?*/ /*double lowsigma[5]; double highsigma[5]; TH1 *lowshape[5]; //TD changed from TH1F for now. TH1 *highshape[5]; //TD changed from TH1F for now. for (int i=0;i<5;i++) { lowshape[i] = 0; highshape[i] = 0; lowsigma[i] =1; highsigma[i] =1; }*/ TDHistSet TomLikeBG(*(HistSets->GetHistSet("Zee"))); TomLikeBG+=(*(HistSets->GetHistSet("Zmumu"))); TomLikeBG+=(*(HistSets->GetHistSet("Ztautau"))); TomLikeBG+=(*(HistSets->GetHistSet("wgamma"))); TomLikeBG+=(*(HistSets->GetHistSet("ttbar"))); // TomLikeBG+=(*(HistSets->GetHistSet("WW"))); TomLikeBG+=(*(HistSets->GetHistSet("WZ"))); TomLikeBG+=(*(HistSets->GetHistSet("ZZ"))); TDHistSet TomLikeWW(*(HistSets->GetHistSet("WW"))); //TH1F * bkgr=(TH1F*)(TomLikeBG.GetAnaHist("dilDphi")->GetHist()); TDHistSet TomLikeSig(*(HistSets->GetHistSet("higgsww"))); //TH1F * sig=(TH1F*)(TomLikeSig.GetAnaHist("dilDphi")->GetHist()); //TH1F * data=(TH1F*)((HistSets->GetHistSet("data"))->GetAnaHist("dilDphi")->GetHist()); TDHistSet TomLikeDat(*(HistSets->GetHistSet("data"))); char buffer[200]; for(int type=DilType::TCE_TCE; type<=DilType::CMX_CMIO; ++type) //TD add each dilDphi from the diltypes separately { if(!(type==2 || type==3 || type==4 || type==15 || type==16)) //TD no pem, PHXPHX { sprintf(buffer,"Channel %s", (_dilType.Int2String(type)).c_str());//TD name the channels after the diltypes //Background //nps_high[1] bgerr[14] = ((TH1F*)(TomLikeBG.GetAnaHist(((_dilType.Int2String(type)) + " relweighterror").c_str())->GetHist()))->GetMean(); //nps_low[1] bgderr[14] = -bgerr[14]; //TD errors on the different dilDphi from the weightings is mean of the relative errors in their diltypes (TH1F*)(TomLikeBG.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist())->Rebin(nbin); //TD rebinned to a more healthy size nullhyp->add_template((TH1F*)(TomLikeBG.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist()), 1.0,15,ername,bgderr,bgerr, bgdownshape,bgdownshapesigma,bgupshape,bgupshapesigma,0,0,buffer); //TD add to null template for channel; 2 nps. testhyp->add_template((TH1F*)(TomLikeBG.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist()), 1.0,15,ername,bgderr,bgerr, bgdownshape,bgdownshapesigma,bgupshape,bgupshapesigma,0,0,buffer); //TD add to test template for channel; 2 nps. //WW Background wwerr[14] = ((TH1F*)(TomLikeWW.GetAnaHist(((_dilType.Int2String(type)) + " relweighterror").c_str())->GetHist()))->GetMean(); wwderr[14] = -wwerr[14]; (TH1F*)(TomLikeWW.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist())->Rebin(nbin); //TD rebinned to a more healthy size nullhyp->add_template((TH1F*)(TomLikeWW.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist()), 1.0,15,ername,wwderr,wwerr, wwdownshape,wwdownshapesigma,wwupshape,wwupshapesigma,0,0,buffer); //TD add to null template for channel; 2 nps. testhyp->add_template((TH1F*)(TomLikeWW.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist()), 1.0,15,ername,wwderr,wwerr, wwdownshape,wwdownshapesigma,wwupshape,wwupshapesigma,0,0,buffer); //TD add to test template for channel; 2 nps. //Signal //nps_high[1] i serr[14] = ((TH1F*)(TomLikeSig.GetAnaHist(((_dilType.Int2String(type)) + " relweighterror").c_str())->GetHist()))->GetMean(); //nps_low[1] = -nps_high[1]; //TD errors on the different dilDphi from the weightings is mean of the relative errors in their diltypes sderr[14] = -serr[14]; (TH1F*)(TomLikeSig.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist())->Rebin(nbin); testhyp->add_template((TH1F*)(TomLikeSig.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist()), 1.0,15,ername,sderr,serr, sigdownshape,sigdownshapesigma,sigupshape,sigupshapesigma,0,1,buffer); } } csm_model* nullhyp_pe = nullhyp->Clone(); csm_model* testhyp_pe = testhyp->Clone(); //TCanvas * mycanvas = (TCanvas *) new TCanvas("Canvas1","Canvas1"); //int dilcount=0; //for(int type=DilType::TCE_TCE; type<=DilType::CMX_CMIO; ++type)dilcount++; //mycanvas->Divide(2,1); //mycanvas->cd(1); /*for(int type=DilType::TCE_TCE; type<=DilType::CMX_CMIO; ++type) //TD add each dilDphi from the diltypes separately { sprintf(buffer,"Channel %s", _dilType.Int2String(type))//TD name the channels after the diltypes (TH1F*)(TomLikeDat.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist())->Rebin(nbin); testhyp->plotwithdata(buffer, (TH1F*)(TomLikeDat.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist())); }*/ //TD can't manage this as yet due to canvas problems csm* mycsm = new csm(); for(int type=DilType::TCE_TCE; type<=DilType::CMX_CMIO; ++type) //TD add each dilDphi from the diltypes separately { sprintf(buffer,"Channel %s", (_dilType.Int2String(type)).c_str()); (TH1F*)(TomLikeDat.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist())->Rebin(nbin); mycsm->set_htofit((TH1F*)(TomLikeDat.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist()), buffer); } mycsm->set_modeltofit(testhyp); double chisq = mycsm->chisquared(); std::cout<<"chisq: "<getbestmodel(); //TD not plotted, to not used.. //mycanvas->cd(2); //bestsignalfit->plotwithdata("Channel 1", data); //mycanvas->Print("results/tomlikeout.ps"); mclimit_csm* mymclimit = new mclimit_csm(); mymclimit->set_null_hypothesis(nullhyp); mymclimit->set_test_hypothesis(testhyp); mymclimit->set_null_hypothesis_pe(nullhyp_pe); mymclimit->set_test_hypothesis_pe(testhyp_pe); for(int type=DilType::TCE_TCE; type<=DilType::CMX_CMIO; ++type) //TD add each dilDphi from the diltypes separately { sprintf(buffer,"Channel %s", (_dilType.Int2String(type)).c_str()); mymclimit->set_datahist((TH1F*)(TomLikeDat.GetAnaHist(((_dilType.Int2String(type)) + " dilDphilike").c_str())->GetHist()), buffer); //mymclimit->set_datahist(data, "Channel 1"); } mymclimit->set_npe(1000); mymclimit->run_pseudoexperiments(); std::cout << "Getting results" << std::endl; std::cout << "ts: " << mymclimit->ts() << std::endl; std::cout << "tsbm2: " << mymclimit->tsbm2() << std::endl; std::cout << "tsbm1: " << mymclimit->tsbm1() << std::endl; std::cout << "tsbmed: " << mymclimit->tsbmed() << std::endl; std::cout << "tsbp1: " << mymclimit->tsbp1() << std::endl; std::cout << "tsbp2: " << mymclimit->tsbp2() << std::endl; std::cout << "tssm2: " << mymclimit->tssm2() << std::endl; std::cout << "tssm1: " << mymclimit->tssm1() << std::endl; std::cout << "tssmed: " << mymclimit->tssmed() << std::endl; std::cout << "tssp1: " << mymclimit->tssp1() << std::endl; std::cout << "tssp2: " << mymclimit->tssp2() << std::endl; std::cout << "CLs: " << mymclimit->cls() << std::endl; std::cout << "CLb: " << mymclimit->clb() << std::endl; std::cout << "CLsb: " << mymclimit->clsb() << std::endl; // std::cout << "Median expected s95 in bg hypothesis: " << mymclimit->s95med() << std::endl; // fitinfo << "Median expected s95 in bg hypothesis: " << mymclimit->s95med() << std::endl; // std::cout<<"deleting Tomlike elements"<bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); // limitcalc->bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, // npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); cout << "individual dilepton likelihood mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << "nbin:"<GetTotGoodWW()) xs=GetSubAna("data")->GetTotGoodWW()*/ std::cout<<"ratio between data and mc xs is"<< (HistSets->GetHistSet("data"))->GetAnaHist("WWmetmag")->GetHist()->Integral()/ (HistSets->GetHistSet("WW"))->GetAnaHist("WWmetmag")->GetHist()->Integral()<GetHistSet("data"))->GetAnaHist("WWmetmag")->GetHist()->Integral()/ (HistSets->GetHistSet("WW"))->GetAnaHist("WWmetmag")->GetHist()->Integral()<sysname.push_back(sname); ername[0] = "LUMI"; serr[0] = 0.04; sderr[0] = -0.04; wwerr[0] = 0.04; wwderr[0] = -0.04; bgerr[0] = 0.04; bgderr[0] = -0.04; // uncorrelated part of lumi error //sname = new char[strlen("CDFLUMI")+1]; //strcpy(sname,"CDFLUMI"); //chanresult->sysname.push_back(sname); ername[1] = "CDFLUMI"; serr[1] = 0.04; sderr[1] = -0.04; wwerr[1] = 0.04; wwderr[1] = -0.04; bgerr[1] = 0.04; bgderr[1] = -0.04; //sname = new char[strlen("CDFLEPTRIGGER")+1]; //strcpy(sname,"CDFLEPTRIGGER"); //chanresult->sysname.push_back(sname); ername[2] = "CDFLEPTRIGGER"; serr[2] = 0.01; sderr[2] = -0.01; wwerr[2] = 0.01; wwderr[2] = -0.01; bgerr[2] = 0.01; // only partly estimated in CDF 7708 bgderr[2] = -0.01; // only partly estimated in CDF 7708 //sname = new char[strlen("ISR")+1]; //strcpy(sname,"ISR"); //chanresult->sysname.push_back(sname); ername[3] = "ISR"; serr[3] = 0.03; sderr[3] = -0.03; wwerr[3] = 0.03; // not estimated in CDF 7708 wwderr[3] = -0.03; // not estimated in CDF 7708 bgerr[3] = 0.02; // partly estimated in CDF 7708 bgderr[3] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("PDF")+1]; //strcpy(sname,"PDF"); //chanresult->sysname.push_back(sname); ername[4] = "PDF"; serr[4] = 0.03; sderr[4] = -0.03; wwerr[4] = 0.06; wwderr[4] = -0.06; bgerr[4] = 0.04; bgderr[4] = -0.04; //sname = new char[strlen("CDFJES")+1]; //strcpy(sname,"CDFJES"); //chanresult->sysname.push_back(sname); ername[5] = "CDFJES"; serr[5] = 0.01; sderr[5] = -0.01; wwerr[5] = 0.01; wwderr[5] = -0.01; bgerr[5] = 0.02; // mostly from fakes bgderr[5] = -0.02; // mostly from fakes //sname = new char[strlen("CDFTRKISO")+1]; //strcpy(sname,"CDFTRKISO"); //chanresult->sysname.push_back(sname); ername[6] = "CDFTRKISO"; serr[6] = 0.02; sderr[6] = -0.02; wwerr[6] = 0.02; wwderr[6] = -0.02; bgerr[6] = 0.02; // partly estimated in CDF 7708 bgderr[6] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("ALPHAS")+1]; //strcpy(sname,"ALPHAS"); //chanresult->sysname.push_back(sname); ername[7] = "ALPHAS"; serr[7] = 0.03; sderr[7] = -0.03; wwerr[7] = 0.0; wwderr[7] = 0.0; bgerr[7] = 0.0; bgderr[7] = 0.0; //sname = new char[strlen("CDFLEPTID")+1]; //strcpy(sname,"CDFLEPTID"); //chanresult->sysname.push_back(sname); ername[8] = "CDFLEPTID"; serr[8] = 0.02; sderr[8] = -0.02; wwerr[8] = 0.02; wwderr[8] = -0.02; bgerr[8] = 0.02; bgderr[8] = -0.02; //sname = new char[strlen("CDFFAKERATE")+1]; //strcpy(sname,"CDFFAKERATE"); //chanresult->sysname.push_back(sname); ername[9] = "CDFFAKERATE"; serr[9] = 0.0; sderr[9] = 0.0; wwerr[9] = 0.0; wwderr[9] = 0.0; bgerr[9] = 0.05; // 50% error on fakes, which are 10% of background bgderr[9] = -0.05; // 50% error on fakes, which are 10% of background //sname = new char[strlen("NNLO")+1]; //strcpy(sname,"NNLO"); //chanresult->sysname.push_back(sname); ername[10] = "NNLO"; serr[10] = 0.1; sderr[10] = -0.1; wwerr[10] = 0.0; wwderr[10] = 0.0; bgerr[10] = 0.0; bgderr[10] = 0.0; //sname = new char[strlen("CDFJETVETO")+1]; //strcpy(sname,"CDFJETVETO"); //chanresult->sysname.push_back(sname); ername[11] = "CDFJETVETO"; serr[11] = 0.0; sderr[11] = 0.0; wwerr[11] = 0.06; wwderr[11] = -0.06; bgerr[11] = 0.0; bgderr[11] = 0.0; //sname = new char[strlen("GENERATOR")+1]; //strcpy(sname,"GENERATOR"); //chanresult->sysname.push_back(sname); ername[12] = "GENERATOR"; serr[12] = 0.0; sderr[12] = 0.0; wwerr[12] = 0.04; wwderr[12] = -0.04; bgerr[12] = 0.06; bgderr[12] = -0.06; //sname = new char[strlen("CDFDYMET")+1]; //strcpy(sname,"CDFDYMET"); //chanresult->sysname.push_back(sname); ername[13] = "CDFDYMET"; serr[13] = 0.0; sderr[13] = 0.0; wwerr[13] = 0.0; wwderr[13] = 0.0; bgerr[13] = 0.05; bgderr[13] = -0.05; TH1 * sigupshape[14]; TH1 * sigdownshape[14]; double sigupshapesigma[14]; double sigdownshapesigma[14]; TH1 * wwupshape[14]; TH1 * wwdownshape[14]; double wwupshapesigma[14]; double wwdownshapesigma[14]; TH1 * bgupshape[14]; TH1 * bgdownshape[14]; double bgupshapesigma[14]; double bgdownshapesigma[14]; for (int i=0;i<14;i++) { sigupshape[i] = 0; sigdownshape[i] = 0; sigupshapesigma[i] = 0; sigdownshapesigma[i] = 0; wwupshape[i] = 0; wwdownshape[i] = 0; wwupshapesigma[i] = 0; wwdownshapesigma[i] = 0; bgupshape[i] = 0; bgdownshape[i] = 0; bgupshapesigma[i] = 0; bgdownshapesigma[i] = 0; } TDHistSet TomLikeBG(*(HistSets->GetHistSet("Zee"))); TomLikeBG+=(*(HistSets->GetHistSet("Zmumu"))); TomLikeBG+=(*(HistSets->GetHistSet("Ztautau"))); TomLikeBG+=(*(HistSets->GetHistSet("wgamma"))); TomLikeBG+=(*(HistSets->GetHistSet("ttbar"))); //TomLikeBG+=(*(HistSets->GetHistSet("WW"))); TomLikeBG+=(*(HistSets->GetHistSet("WZ"))); TomLikeBG+=(*(HistSets->GetHistSet("ZZ"))); TDHistSet TomLikeWW(*(HistSets->GetHistSet("WW"))); TH1F * bkgr=(TH1F*)(TomLikeBG.GetAnaHist(thehist)->GetHist()); TH1F * WWbkgr=(TH1F*)(TomLikeWW.GetAnaHist(thehist)->GetHist()); TDHistSet TomLikeSig(*(HistSets->GetHistSet("higgsww"))); TH1F * sig=(TH1F*)(TomLikeSig.GetAnaHist(thehist)->GetHist()); TH1F * data=(TH1F*)((HistSets->GetHistSet("data"))->GetAnaHist(thehist)->GetHist()); nullhyp->add_template(bkgr,1.0,14,ername,bgderr,bgerr, bgdownshape,bgdownshapesigma,bgupshape,bgupshapesigma,0,0,"Channel 1"); testhyp->add_template(bkgr,1.0,14,ername,bgderr,bgerr, bgdownshape,bgdownshapesigma,bgupshape,bgupshapesigma,0,0,"Channel 1"); nullhyp->add_template(WWbkgr,1.0,14,ername,wwderr,wwerr, wwdownshape,wwdownshapesigma,wwupshape,wwupshapesigma,0,0,"Channel 1"); testhyp->add_template(WWbkgr,1.0,14,ername,wwderr,wwerr, wwdownshape,wwdownshapesigma,wwupshape,wwupshapesigma,0,0,"Channel 1"); testhyp->add_template(sig,1.0,14,ername,sderr,serr, sigdownshape,sigdownshapesigma,sigupshape,sigupshapesigma,0,1,"Channel 1"); csm_model* nullhyp_pe = nullhyp->Clone(); csm_model* testhyp_pe = testhyp->Clone(); TCanvas * mycanvas = (TCanvas *) new TCanvas("Canvas1","Canvas1"); mycanvas->Divide(2,1); mycanvas->cd(1); testhyp->plotwithdata("Channel 1",data ); csm* mycsm = new csm(); mycsm->set_htofit( data , "Channel 1"); mycsm->set_modeltofit(testhyp); double chisq = mycsm->chisquared(); std::cout<<"chisq: "<getbestmodel(); mycanvas->cd(2); bestsignalfit->plotwithdata("Channel 1", data); char buffer[200]; sprintf(buffer, "results/tomlikeouttop_%s.ps", thehist); mycanvas->Print(buffer); mclimit_csm* mymclimit = new mclimit_csm(); mymclimit->set_null_hypothesis(nullhyp); mymclimit->set_test_hypothesis(testhyp); mymclimit->set_null_hypothesis_pe(nullhyp_pe); mymclimit->set_test_hypothesis_pe(testhyp_pe); mymclimit->set_datahist(data, "Channel 1"); mymclimit->set_npe(1000); mymclimit->run_pseudoexperiments(); std::cout << "Getting results" << std::endl; std::cout << "ts: " << mymclimit->ts() << std::endl; std::cout << "tsbm2: " << mymclimit->tsbm2() << std::endl; std::cout << "tsbm1: " << mymclimit->tsbm1() << std::endl; std::cout << "tsbmed: " << mymclimit->tsbmed() << std::endl; std::cout << "tsbp1: " << mymclimit->tsbp1() << std::endl; std::cout << "tsbp2: " << mymclimit->tsbp2() << std::endl; std::cout << "tssm2: " << mymclimit->tssm2() << std::endl; std::cout << "tssm1: " << mymclimit->tssm1() << std::endl; std::cout << "tssmed: " << mymclimit->tssmed() << std::endl; std::cout << "tssp1: " << mymclimit->tssp1() << std::endl; std::cout << "tssp2: " << mymclimit->tssp2() << std::endl; std::cout << "CLs: " << mymclimit->cls() << std::endl; std::cout << "CLb: " << mymclimit->clb() << std::endl; std::cout << "CLsb: " << mymclimit->clsb() << std::endl; double fit95=mymclimit->s95med(); std::cout << "Median expected s95 in bg hypothesis: " << fit95 << std::endl; fitinfo << "Median expected s95 in bg hypothesis: " << fit95 << std::endl; double fit95act=mymclimit->s95(); //std::cout<<"deleting Tomlike elements"<bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); // limitcalc->bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, // npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); cout << "mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << endl; fitinfo << "mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << endl; delete mycanvas; /*if(nullhyp!=NULL)delete nullhyp; if(testhyp!=NULL)delete testhyp; if(mycsm!=NULL)delete mycsm; if(bestsignalfit!=NULL)delete bestsignalfit;*/ //TD tom's stuff does not like to be deleted... memory leak! std::cout<<"deleted Tomlike elements"<llvv) = 0.1068 in PYTHIA vs 0.11 in HERWIG const double nEvt[10] = {1069855.0, 1063469.0, 1127915.0, 1074251.0, 1089119.0, 1051020.0, 1033214.0, 1079780.0, 1092527.0, 1093662.0}; //int main() //void pseudolimits() { int mHig=170; //char *hmas = gSystem->Getenv("HiggsMass"); //if (hmas) for (int i=0; iGetenv("LimiLevel"); if (lvll) for (int i=0; i0) sprintf(filename,"mircresults/%s_%d_phxCorr.root",datset[i],mHig); else sprintf(filename,"mircresults/%s%d_%d_phxCorr.root",datset[i],mHig,mHig); f[i] = new TFile(filename,"READ"); } sprintf(filename,"dildPhi_all_%d",lLvl); TH2F* g = (TH2F*)f[0]->Get("GC")->Clone(); //TD what is this? double nAcc=0.0; for (int j=0; j<26; j++) nAcc += g->GetBinContent(lLvl+4,j+1); //TD what is this? double Hacc = 0.1068 * nAcc / nEvt[((int)_higgsMass-110)/10]; //TD where do the n Evt come from? //double Hacc = nAcc / nEvt[((int)_higgsMass-110)/10]; // double sHacc = 0.1068 * 0.06 * Hacc; double sHacc = 0.06 * Hacc; for (int i=0; i<9; i++) { h[i] = (TH1F*)f[i]->Get(filename)->Clone(); h[i]->Rebin(20); if (i>2 && i<8) h[2]->Add(h[i]); } int nBins = h[0]->GetNbinsX(); //int nBins = HistSets->GetHistSet("higgsww")->GetAnaHist("dilDphilike")->GetHist()->GetNbinsX(); double nWW = h[1]->Integral(); double nBG = h[2]->Integral(); double sWW=0.0, sBG=0.0, norm_hww=0.0, norm_ww=0.0, norm_bg=0.0; double data_dil[200], hww_all[200], ww_sm[200], bg_sm[200]; for (int b=0; bGetBinContent(b+1); ww_sm[b] = h[1]->GetBinContent(b+1); bg_sm[b] = h[2]->GetBinContent(b+1); norm_hww += hww_all[b]; norm_ww += ww_sm[b]; norm_bg += bg_sm[b]; } sWW = 0.09 * nWW; sBG = 0.17 * nBG; cout << " Higgs Mass = " << (int)_higgsMass << " GeV; Cut Level = " << lLvl << endl; //TD someone has included a 'using std' some where cout << " Luminosity = " << Lum << "+-" << sLum << endl; cout << " HWW A*BR = " << Hacc << "+-" << sHacc << endl; cout << " WW BG = " << nWW << "+-" << sWW << endl; cout << " otherBGs = " << nBG << "+-" << sBG << endl << endl; fitinfo << " Higgs Mass = " << (int)_higgsMass << " GeV; Cut Level = " << lLvl << std::endl; fitinfo << " Luminosity = " << Lum << "+-" << sLum << std::endl; fitinfo << " HWW A*BR = " << Hacc << "+-" << sHacc << std::endl; fitinfo << " WW BG = " << nWW << "+-" << sWW << std::endl; fitinfo << " otherBGs = " << nBG << "+-" << sBG << std::endl << std::endl; norm_hww = 1./ norm_hww; norm_ww = 1./ norm_ww; norm_bg = 1./ norm_bg; char buffer[200]; char buffer2[200]; sprintf(filename,"results/psLm_phxCorr_%3d_%d.root",(int)_higgsMass,lLvl); TFile* v = new TFile(filename,"RECREATE"); if(sim)sprintf(buffer,"simcrosseclim"); else sprintf(buffer,"crosseclim"); if(sim)sprintf(buffer2,"simulated crossection limit at 95%s CL", "%"); else sprintf(buffer2,"crossection limit at 95%s CL","%"); TH1F* xsec = new TH1F(buffer,buffer2,1000,0.0,50.0); //TD memory leak! TH1F* xlik[1000]; TRandom* r = new TRandom(); int npseudo=1; if(sim)npseudo=1000; for (int mn=0; mnGaus(nWW, sWW); int nWW_poisson = r->Poisson(nWW_gauss); for (int j=0; jGetRandom(); data_dil[int(nBins*dphi/3.1416)]++; } double nBG_gauss = r->Gaus(nBG, sBG); int nBG_poisson = r->Poisson(nBG_gauss); for (int j=0; jGetRandom(); data_dil[int(nBins*dphi/3.1416)]++; } }else{ for (int b=0; bGetBinContent(b+1); } double par[5], lik[1000], intgratlik=0.0, intpartlik=0.0, xlim; for (int j1=0; j1<1000; j1++) { lik[j1]=0.0; for (int it=0; it<100; it++) { par[0] = 0.05*j1; par[1] = r->Gaus(Lum, sLum); par[2] = r->Gaus(Hacc,sHacc); par[3] = r->Gaus(nWW, sWW); par[4] = r->Gaus(nBG, sBG); double mu=0.0, loglik=0.0; for (int b=0; b0) loglik += (data_dil[b]*TMath::Log(mu) - mu); } lik[j1] += TMath::Exp(loglik); } lik[j1] /= 100.0; intgratlik += lik[j1]; xlik[mn]->Fill(par[0], lik[j1]); } int j1=0; while ((intpartlik/intgratlik)<0.95) { intpartlik += lik[j1]; xlim = 0.05*j1; j1++; } xsec->Fill(xlim); //cout << " xlim[" << mn << "] = " << xlim << " pb" << endl; // fitinfo << " xlim[" << mn << "] = " << xlim << " pb" << endl; if (mn>0 && fmod((float)mn,(float)100)==0) { double xsec_median[1], x[1]; x[0]=0.5; xsec->GetQuantiles(1,xsec_median,x); // cout << " Crossec Lim Med = " << xsec_median[0] << " pb" << endl << endl; } } cout << " Cross-section Limit Mean = " << xsec->GetMean() << " pb" << endl; fitinfo << " Cross-section Limit Mean = " << xsec->GetMean() << " pb" << endl; double xsec_median[1], x[1]; x[0]=0.5; xsec->GetQuantiles(1,xsec_median,x); cout << " Crossection Limit Median = " << xsec_median[0] << " pb" << endl; cout << " Crossection Limit Median as ratio = " << xsec_median[0]/(GetHwwXS((int)_higgsMass, _stn)*0.1068) << endl << endl; if(sim) { fitinfo << "expected:"<Draw(); xsec->Write(); if(!sim)xlik[0]->Write(); //v->Write(); v->Close(); //for (int i=0; i<10; i++) f[i]->Close(); csm_model* nullhyp = new csm_model(); csm_model* testhyp = new csm_model(); char *ername[14]; double serr[14]; double sderr[14]; double wwerr[14]; double wwderr[14]; double bgerr[14]; double bgderr[14]; //n sname = new char[strlen("LUMI")+1]; // strcpy(sname,"LUMI"); // chanresult->sysname.push_back(sname); ername[0] = "LUMI"; serr[0] = 0.04; sderr[0] = -0.04; wwerr[0] = 0.04; wwderr[0] = -0.04; bgerr[0] = 0.04; bgderr[0] = -0.04; // uncorrelated part of lumi error //sname = new char[strlen("CDFLUMI")+1]; //strcpy(sname,"CDFLUMI"); //chanresult->sysname.push_back(sname); ername[1] = "CDFLUMI"; serr[1] = 0.04; sderr[1] = -0.04; wwerr[1] = 0.04; wwderr[1] = -0.04; bgerr[1] = 0.04; bgderr[1] = -0.04; //sname = new char[strlen("CDFLEPTRIGGER")+1]; //strcpy(sname,"CDFLEPTRIGGER"); //chanresult->sysname.push_back(sname); ername[2] = "CDFLEPTRIGGER"; serr[2] = 0.01; sderr[2] = -0.01; wwerr[2] = 0.01; wwderr[2] = -0.01; bgerr[2] = 0.01; // only partly estimated in CDF 7708 bgderr[2] = -0.01; // only partly estimated in CDF 7708 //sname = new char[strlen("ISR")+1]; //strcpy(sname,"ISR"); //chanresult->sysname.push_back(sname); ername[3] = "ISR"; serr[3] = 0.03; sderr[3] = -0.03; wwerr[3] = 0.03; // not estimated in CDF 7708 wwderr[3] = -0.03; // not estimated in CDF 7708 bgerr[3] = 0.02; // partly estimated in CDF 7708 bgderr[3] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("PDF")+1]; //strcpy(sname,"PDF"); //chanresult->sysname.push_back(sname); ername[4] = "PDF"; serr[4] = 0.03; sderr[4] = -0.03; wwerr[4] = 0.06; wwderr[4] = -0.06; bgerr[4] = 0.04; bgderr[4] = -0.04; //sname = new char[strlen("CDFJES")+1]; //strcpy(sname,"CDFJES"); //chanresult->sysname.push_back(sname); ername[5] = "CDFJES"; serr[5] = 0.01; sderr[5] = -0.01; wwerr[5] = 0.01; wwderr[5] = -0.01; bgerr[5] = 0.02; // mostly from fakes bgderr[5] = -0.02; // mostly from fakes //sname = new char[strlen("CDFTRKISO")+1]; //strcpy(sname,"CDFTRKISO"); //chanresult->sysname.push_back(sname); ername[6] = "CDFTRKISO"; serr[6] = 0.02; sderr[6] = -0.02; wwerr[6] = 0.02; wwderr[6] = -0.02; bgerr[6] = 0.02; // partly estimated in CDF 7708 bgderr[6] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("ALPHAS")+1]; //strcpy(sname,"ALPHAS"); //chanresult->sysname.push_back(sname); ername[7] = "ALPHAS"; serr[7] = 0.03; sderr[7] = -0.03; wwerr[7] = 0.0; wwderr[7] = 0.0; bgerr[7] = 0.0; bgderr[7] = 0.0; //sname = new char[strlen("CDFLEPTID")+1]; //strcpy(sname,"CDFLEPTID"); //chanresult->sysname.push_back(sname); ername[8] = "CDFLEPTID"; serr[8] = 0.02; sderr[8] = -0.02; wwerr[8] = 0.02; wwderr[8] = -0.02; bgerr[8] = 0.02; bgderr[8] = -0.02; //sname = new char[strlen("CDFFAKERATE")+1]; //strcpy(sname,"CDFFAKERATE"); //chanresult->sysname.push_back(sname); ername[9] = "CDFFAKERATE"; serr[9] = 0.0; sderr[9] = 0.0; wwerr[9] = 0.0; wwderr[9] = 0.0; bgerr[9] = 0.05; // 50% error on fakes, which are 10% of background bgderr[9] = -0.05; // 50% error on fakes, which are 10% of background //sname = new char[strlen("NNLO")+1]; //strcpy(sname,"NNLO"); //chanresult->sysname.push_back(sname); ername[10] = "NNLO"; serr[10] = 0.1; sderr[10] = -0.1; wwerr[10] = 0.0; wwderr[10] = 0.0; bgerr[10] = 0.0; bgderr[10] = 0.0; //sname = new char[strlen("CDFJETVETO")+1]; //strcpy(sname,"CDFJETVETO"); //chanresult->sysname.push_back(sname); ername[11] = "CDFJETVETO"; serr[11] = 0.0; sderr[11] = 0.0; wwerr[11] = 0.06; wwderr[11] = -0.06; bgerr[11] = 0.0; bgderr[11] = 0.0; //sname = new char[strlen("GENERATOR")+1]; //strcpy(sname,"GENERATOR"); //chanresult->sysname.push_back(sname); ername[12] = "GENERATOR"; serr[12] = 0.0; sderr[12] = 0.0; wwerr[12] = 0.04; wwderr[12] = -0.04; bgerr[12] = 0.06; bgderr[12] = -0.06; //sname = new char[strlen("CDFDYMET")+1]; //strcpy(sname,"CDFDYMET"); //chanresult->sysname.push_back(sname); ername[13] = "CDFDYMET"; serr[13] = 0.0; sderr[13] = 0.0; wwerr[13] = 0.0; wwderr[13] = 0.0; bgerr[13] = 0.05; bgderr[13] = -0.05; TH1 * sigupshape[14]; TH1 * sigdownshape[14]; double sigupshapesigma[14]; double sigdownshapesigma[14]; TH1 * wwupshape[14]; TH1 * wwdownshape[14]; double wwupshapesigma[14]; double wwdownshapesigma[14]; TH1 * bgupshape[14]; TH1 * bgdownshape[14]; double bgupshapesigma[14]; double bgdownshapesigma[14]; for (int i=0;i<14;i++) { sigupshape[i] = 0; sigdownshape[i] = 0; sigupshapesigma[i] = 0; sigdownshapesigma[i] = 0; wwupshape[i] = 0; wwdownshape[i] = 0; wwupshapesigma[i] = 0; wwdownshapesigma[i] = 0; bgupshape[i] = 0; bgdownshape[i] = 0; bgupshapesigma[i] = 0; bgdownshapesigma[i] = 0; } TH1F * bkgr=(TH1F*)h[2]; TH1F * WWbkgr=(TH1F*)h[1]; TH1F * sig=(TH1F*)h[0]; //sig->Scale(1.45); TH1F * data=(TH1F*)h[8]; nullhyp->add_template(bkgr,1.0,14,ername,bgderr,bgerr, bgdownshape,bgdownshapesigma,bgupshape,bgupshapesigma,0,0,"Channel 1"); testhyp->add_template(bkgr,1.0,14,ername,bgderr,bgerr, bgdownshape,bgdownshapesigma,bgupshape,bgupshapesigma,0,0,"Channel 1"); nullhyp->add_template(WWbkgr,1.0,14,ername,wwderr,wwerr, wwdownshape,wwdownshapesigma,wwupshape,wwupshapesigma,0,0,"Channel 1"); testhyp->add_template(WWbkgr,1.0,14,ername,wwderr,wwerr, wwdownshape,wwdownshapesigma,wwupshape,wwupshapesigma,0,0,"Channel 1"); testhyp->add_template(sig,1.0,14,ername,sderr,serr, sigdownshape,sigdownshapesigma,sigupshape,sigupshapesigma,0,1,"Channel 1"); csm_model* nullhyp_pe = nullhyp->Clone(); csm_model* testhyp_pe = testhyp->Clone(); TCanvas * mycanvas = (TCanvas *) new TCanvas("Canvas1","Canvas1"); mycanvas->Divide(2,1); mycanvas->cd(1); testhyp->plotwithdata("Channel 1",data ); csm* mycsm = new csm(); mycsm->set_htofit( data , "Channel 1"); mycsm->set_modeltofit(testhyp); double chisq = mycsm->chisquared(); std::cout<<"chisq: "<getbestmodel(); mycanvas->cd(2); bestsignalfit->plotwithdata("Channel 1", data); // mycanvas->Print("results/tomlikeouttop.ps"); mclimit_csm* mymclimit = new mclimit_csm(); mymclimit->set_null_hypothesis(nullhyp); mymclimit->set_test_hypothesis(testhyp); mymclimit->set_null_hypothesis_pe(nullhyp_pe); mymclimit->set_test_hypothesis_pe(testhyp_pe); mymclimit->set_datahist(data, "Channel 1"); mymclimit->set_npe(1000); mymclimit->run_pseudoexperiments(); std::cout << "Getting results" << std::endl; std::cout << "ts: " << mymclimit->ts() << std::endl; std::cout << "tsbm2: " << mymclimit->tsbm2() << std::endl; std::cout << "tsbm1: " << mymclimit->tsbm1() << std::endl; std::cout << "tsbmed: " << mymclimit->tsbmed() << std::endl; std::cout << "tsbp1: " << mymclimit->tsbp1() << std::endl; std::cout << "tsbp2: " << mymclimit->tsbp2() << std::endl; std::cout << "tssm2: " << mymclimit->tssm2() << std::endl; std::cout << "tssm1: " << mymclimit->tssm1() << std::endl; std::cout << "tssmed: " << mymclimit->tssmed() << std::endl; std::cout << "tssp1: " << mymclimit->tssp1() << std::endl; std::cout << "tssp2: " << mymclimit->tssp2() << std::endl; std::cout << "CLs: " << mymclimit->cls() << std::endl; std::cout << "CLb: " << mymclimit->clb() << std::endl; std::cout << "CLsb: " << mymclimit->clsb() << std::endl; double fit95=mymclimit->s95med(); std::cout << "Median expected s95 in bg hypothesis: " << fit95 << std::endl; fitinfo << "Median expected s95 in bg hypothesis: " << fit95 << std::endl; std::cout<<"deleting Tomlike elements"<bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); // limitcalc->bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, // npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); cout << "mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << endl; fitinfo << "mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << endl; delete mycanvas; /*if(nullhyp!=NULL)delete nullhyp; if(testhyp!=NULL)delete testhyp; if(mycsm!=NULL)delete mycsm; if(bestsignalfit!=NULL)delete bestsignalfit;*/ //TD tom's stuff does not like to be deleted... memory leak! // std::cout<<"deleted Tomlike elements"<sysname.push_back(sname); ername[0] = "LUMI"; serr[0] = 0.04; sderr[0] = -0.04; wwerr[0] = 0.04; wwderr[0] = -0.04; bgerr[0] = 0.04; bgderr[0] = -0.04; // uncorrelated part of lumi error //sname = new char[strlen("CDFLUMI")+1]; //strcpy(sname,"CDFLUMI"); //chanresult->sysname.push_back(sname); ername[1] = "CDFLUMI"; serr[1] = 0.04; sderr[1] = -0.04; wwerr[1] = 0.04; wwderr[1] = -0.04; bgerr[1] = 0.04; bgderr[1] = -0.04; //sname = new char[strlen("CDFLEPTRIGGER")+1]; //strcpy(sname,"CDFLEPTRIGGER"); //chanresult->sysname.push_back(sname); ername[2] = "CDFLEPTRIGGER"; serr[2] = 0.01; sderr[2] = -0.01; wwerr[2] = 0.01; wwderr[2] = -0.01; bgerr[2] = 0.01; // only partly estimated in CDF 7708 bgderr[2] = -0.01; // only partly estimated in CDF 7708 //sname = new char[strlen("ISR")+1]; //strcpy(sname,"ISR"); //chanresult->sysname.push_back(sname); ername[3] = "ISR"; serr[3] = 0.03; sderr[3] = -0.03; wwerr[3] = 0.03; // not estimated in CDF 7708 wwderr[3] = -0.03; // not estimated in CDF 7708 bgerr[3] = 0.02; // partly estimated in CDF 7708 bgderr[3] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("PDF")+1]; //strcpy(sname,"PDF"); //chanresult->sysname.push_back(sname); ername[4] = "PDF"; serr[4] = 0.03; sderr[4] = -0.03; wwerr[4] = 0.06; wwderr[4] = -0.06; bgerr[4] = 0.04; bgderr[4] = -0.04; //sname = new char[strlen("CDFJES")+1]; //strcpy(sname,"CDFJES"); //chanresult->sysname.push_back(sname); ername[5] = "CDFJES"; serr[5] = 0.01; sderr[5] = -0.01; wwerr[5] = 0.01; wwderr[5] = -0.01; bgerr[5] = 0.02; // mostly from fakes bgderr[5] = -0.02; // mostly from fakes //sname = new char[strlen("CDFTRKISO")+1]; //strcpy(sname,"CDFTRKISO"); //chanresult->sysname.push_back(sname); ername[6] = "CDFTRKISO"; serr[6] = 0.02; sderr[6] = -0.02; wwerr[6] = 0.02; wwderr[6] = -0.02; bgerr[6] = 0.02; // partly estimated in CDF 7708 bgderr[6] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("ALPHAS")+1]; //strcpy(sname,"ALPHAS"); //chanresult->sysname.push_back(sname); ername[7] = "ALPHAS"; serr[7] = 0.03; sderr[7] = -0.03; wwerr[7] = 0.0; wwderr[7] = 0.0; bgerr[7] = 0.0; bgderr[7] = 0.0; //sname = new char[strlen("CDFLEPTID")+1]; //strcpy(sname,"CDFLEPTID"); //chanresult->sysname.push_back(sname); ername[8] = "CDFLEPTID"; serr[8] = 0.02; sderr[8] = -0.02; wwerr[8] = 0.02; wwderr[8] = -0.02; bgerr[8] = 0.02; bgderr[8] = -0.02; //sname = new char[strlen("CDFFAKERATE")+1]; //strcpy(sname,"CDFFAKERATE"); //chanresult->sysname.push_back(sname); ername[9] = "CDFFAKERATE"; serr[9] = 0.0; sderr[9] = 0.0; wwerr[9] = 0.0; wwderr[9] = 0.0; bgerr[9] = 0.05; // 50% error on fakes, which are 10% of background bgderr[9] = -0.05; // 50% error on fakes, which are 10% of background //sname = new char[strlen("NNLO")+1]; //strcpy(sname,"NNLO"); //chanresult->sysname.push_back(sname); ername[10] = "NNLO"; serr[10] = 0.1; sderr[10] = -0.1; wwerr[10] = 0.0; wwderr[10] = 0.0; bgerr[10] = 0.0; bgderr[10] = 0.0; //sname = new char[strlen("CDFJETVETO")+1]; //strcpy(sname,"CDFJETVETO"); //chanresult->sysname.push_back(sname); ername[11] = "CDFJETVETO"; serr[11] = 0.0; sderr[11] = 0.0; wwerr[11] = 0.06; wwderr[11] = -0.06; bgerr[11] = 0.0; bgderr[11] = 0.0; //sname = new char[strlen("GENERATOR")+1]; //strcpy(sname,"GENERATOR"); //chanresult->sysname.push_back(sname); ername[12] = "GENERATOR"; serr[12] = 0.0; sderr[12] = 0.0; wwerr[12] = 0.04; wwderr[12] = -0.04; bgerr[12] = 0.06; bgderr[12] = -0.06; //sname = new char[strlen("CDFDYMET")+1]; //strcpy(sname,"CDFDYMET"); //chanresult->sysname.push_back(sname); ername[13] = "CDFDYMET"; serr[13] = 0.0; sderr[13] = 0.0; wwerr[13] = 0.0; wwderr[13] = 0.0; bgerr[13] = 0.05; bgderr[13] = -0.05; TH1 * sigupshape[14]; TH1 * sigdownshape[14]; double sigupshapesigma[14]; double sigdownshapesigma[14]; TH1 * wwupshape[14]; TH1 * wwdownshape[14]; double wwupshapesigma[14]; double wwdownshapesigma[14]; TH1 * bgupshape[14]; TH1 * bgdownshape[14]; double bgupshapesigma[14]; double bgdownshapesigma[14]; for (int i=0;i<14;i++) { sigupshape[i] = 0; sigdownshape[i] = 0; sigupshapesigma[i] = 0; sigdownshapesigma[i] = 0; wwupshape[i] = 0; wwdownshape[i] = 0; wwupshapesigma[i] = 0; wwdownshapesigma[i] = 0; bgupshape[i] = 0; bgdownshape[i] = 0; bgupshapesigma[i] = 0; bgdownshapesigma[i] = 0; } TFile f("/home/tdavies/tomlike/tomtar/cdfhww360.root"); TH1F * bkgr=(TH1F*)f.Get("bg_170;1"); TH1F * WWbkgr=(TH1F*)f.Get("ww_170;1"); TH1F * sig=(TH1F*)f.Get("hww_170;1"); sig->Scale(1.45); TH1F * data=(TH1F*)f.Get("data_170;1"); nullhyp->add_template(bkgr,1.0,14,ername,bgderr,bgerr, bgdownshape,bgdownshapesigma,bgupshape,bgupshapesigma,0,0,"Channel 1"); testhyp->add_template(bkgr,1.0,14,ername,bgderr,bgerr, bgdownshape,bgdownshapesigma,bgupshape,bgupshapesigma,0,0,"Channel 1"); nullhyp->add_template(WWbkgr,1.0,14,ername,wwderr,wwerr, wwdownshape,wwdownshapesigma,wwupshape,wwupshapesigma,0,0,"Channel 1"); testhyp->add_template(WWbkgr,1.0,14,ername,wwderr,wwerr, wwdownshape,wwdownshapesigma,wwupshape,wwupshapesigma,0,0,"Channel 1"); testhyp->add_template(sig,1.0,14,ername,sderr,serr, sigdownshape,sigdownshapesigma,sigupshape,sigupshapesigma,0,1,"Channel 1"); csm_model* nullhyp_pe = nullhyp->Clone(); csm_model* testhyp_pe = testhyp->Clone(); TCanvas * mycanvas = (TCanvas *) new TCanvas("Canvas1","Canvas1"); mycanvas->Divide(2,1); mycanvas->cd(1); testhyp->plotwithdata("Channel 1",data ); csm* mycsm = new csm(); mycsm->set_htofit( data , "Channel 1"); mycsm->set_modeltofit(testhyp); double chisq = mycsm->chisquared(); std::cout<<"chisq: "<getbestmodel(); mycanvas->cd(2); bestsignalfit->plotwithdata("Channel 1", data); // mycanvas->Print("results/tomlikeouttop.ps"); mclimit_csm* mymclimit = new mclimit_csm(); mymclimit->set_null_hypothesis(nullhyp); mymclimit->set_test_hypothesis(testhyp); mymclimit->set_null_hypothesis_pe(nullhyp_pe); mymclimit->set_test_hypothesis_pe(testhyp_pe); mymclimit->set_datahist(data, "Channel 1"); mymclimit->set_npe(1000); mymclimit->run_pseudoexperiments(); std::cout << "Getting results" << std::endl; std::cout << "ts: " << mymclimit->ts() << std::endl; std::cout << "tsbm2: " << mymclimit->tsbm2() << std::endl; std::cout << "tsbm1: " << mymclimit->tsbm1() << std::endl; std::cout << "tsbmed: " << mymclimit->tsbmed() << std::endl; std::cout << "tsbp1: " << mymclimit->tsbp1() << std::endl; std::cout << "tsbp2: " << mymclimit->tsbp2() << std::endl; std::cout << "tssm2: " << mymclimit->tssm2() << std::endl; std::cout << "tssm1: " << mymclimit->tssm1() << std::endl; std::cout << "tssmed: " << mymclimit->tssmed() << std::endl; std::cout << "tssp1: " << mymclimit->tssp1() << std::endl; std::cout << "tssp2: " << mymclimit->tssp2() << std::endl; std::cout << "CLs: " << mymclimit->cls() << std::endl; std::cout << "CLb: " << mymclimit->clb() << std::endl; std::cout << "CLsb: " << mymclimit->clsb() << std::endl; double fit95=mymclimit->s95med(); std::cout << "Median expected s95 in bg hypothesis: " << fit95 << std::endl; fitinfo << "Median expected s95 in bg hypothesis: " << fit95 << std::endl; std::cout<<"deleting Tomlike elements"<bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); // limitcalc->bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, // npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); cout << "mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << endl; fitinfo << "mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << endl; delete mycanvas; /*if(nullhyp!=NULL)delete nullhyp; if(testhyp!=NULL)delete testhyp; if(mycsm!=NULL)delete mycsm; if(bestsignalfit!=NULL)delete bestsignalfit;*/ //TD tom's stuff does not like to be deleted... memory leak! // std::cout<<"deleted Tomlike elements"<GetHistSet("Zee"))); TomLikeBG+=(*(HistSets->GetHistSet("Zmumu"))); TomLikeBG+=(*(HistSets->GetHistSet("Ztautau"))); TomLikeBG+=(*(HistSets->GetHistSet("wgamma"))); TomLikeBG+=(*(HistSets->GetHistSet("ttbar"))); //TomLikeBG+=(*(HistSets->GetHistSet("WW"))); TomLikeBG+=(*(HistSets->GetHistSet("WZ"))); TomLikeBG+=(*(HistSets->GetHistSet("ZZ"))); TDHistSet TomLikeWW(*(HistSets->GetHistSet("WW"))); TDHistSet TomLikeSig(*(HistSets->GetHistSet("higgsww"))); TH1F * bkgree=(TH1F*)(TomLikeBG.GetAnaHist("dilDphiee")->GetHist()); TH1F * WWbkgree=(TH1F*)(TomLikeWW.GetAnaHist("dilDphiee")->GetHist()); TH1F * sigee=(TH1F*)(TomLikeSig.GetAnaHist("dilDphiee")->GetHist()); TH1F * dataee=(TH1F*)((HistSets->GetHistSet("data"))->GetAnaHist("dilDphiee")->GetHist()); TH1F * bkgremu=(TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu")->GetHist()); TH1F * WWbkgremu=(TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu")->GetHist()); TH1F * sigemu=(TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu")->GetHist()); TH1F * dataemu=(TH1F*)((HistSets->GetHistSet("data"))->GetAnaHist("dilDphiemu")->GetHist()); TH1F * bkgrmumu=(TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu")->GetHist()); TH1F * WWbkgrmumu=(TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu")->GetHist()); TH1F * sigmumu=(TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu")->GetHist()); TH1F * datamumu=(TH1F*)((HistSets->GetHistSet("data"))->GetAnaHist("dilDphimumu")->GetHist()); std::cout<<"TomLikelihoodemu()"<sysname.push_back(sname); ername[0] = "LUMI"; serr_emu[0] = serr_mumu[0] = serr_ee[0] = 0.04; sderr_emu[0] = sderr_mumu[0] = sderr_ee[0] = -0.04; wwerr_emu[0] = wwerr_mumu[0] = wwerr_ee[0] = 0.04; wwderr_emu[0] = wwderr_mumu[0] = wwderr_ee[0] = -0.04; bgerr_emu[0] = bgerr_mumu[0] = bgerr_ee[0] = 0.04; bgderr_emu[0] = bgderr_mumu[0] = bgderr_ee[0] = -0.04; // uncorrelated part of lumi error //sname = new char[strlen("CDFLUMI")+1]; //strcpy(sname,"CDFLUMI"); //chanresult->sysname.push_back(sname); ername[1] = "CDFLUMI"; serr_emu[1] = serr_mumu[1] = serr_ee[1] = 0.04; sderr_emu[1] = sderr_mumu[1] = sderr_ee[1] = -0.04; wwerr_emu[1] = wwerr_mumu[1] = wwerr_ee[1] = 0.04; wwderr_emu[1] = wwderr_mumu[1] = wwderr_ee[1] = -0.04; bgerr_emu[1] = bgerr_mumu[1] = bgerr_ee[1] = 0.04; bgderr_emu[1] = bgderr_mumu[1] = bgderr_ee[1] = -0.04; //sname = new char[strlen("CDFLEPTRIGGER")+1]; //strcpy(sname,"CDFLEPTRIGGER"); //chanresult->sysname.push_back(sname); ername[2] = "CDFLEPTRIGGER"; serr_emu[2] = serr_mumu[2] = serr_ee[2] = 0.01; sderr_emu[2] = sderr_mumu[2] = sderr_ee[2] = -0.01; wwerr_emu[2] = wwerr_mumu[2] = wwerr_ee[2] = 0.01; wwderr_emu[2] = wwderr_mumu[2] = wwderr_ee[2] = -0.01; bgerr_emu[2] = bgerr_mumu[2] = bgerr_ee[2] = 0.01; // only partly estimated in CDF 7708 bgderr_emu[2] = bgderr_mumu[2] = bgderr_ee[2] = -0.01; // only partly estimated in CDF 7708 //sname = new char[strlen("ISR")+1]; //strcpy(sname,"ISR"); //chanresult->sysname.push_back(sname); ername[3] = "ISR"; serr_emu[3] = serr_mumu[3] = serr_ee[3] = 0.0;// 0.03; //TD now in shape error hists. sderr_emu[3] = sderr_mumu[3] = sderr_ee[3] = 0.0;//-0.03; wwerr_emu[3] = wwerr_mumu[3] = wwerr_ee[3] = 0.0;//0.03; // not estimated in CDF 7708 wwderr_emu[3] = wwderr_mumu[3] = wwderr_ee[3] = 0.0;//-0.03; // not estimated in CDF 7708 bgerr_emu[3] = bgerr_mumu[3] = bgerr_ee[3] = 0.0;//0.02; // partly estimated in CDF 7708 bgderr_emu[3] = bgderr_mumu[3] = bgderr_ee[3] = 0.0;//-0.02; // partly estimated in CDF 7708 //sname = new char[strlen("PDF")+1]; //strcpy(sname,"PDF"); //chanresult->sysname.push_back(sname); ername[4] = "PDF"; serr_emu[4] = serr_mumu[4] = serr_ee[4] = 0.03; sderr_emu[4] = sderr_mumu[4] = sderr_ee[4] = -0.03; wwerr_emu[4] = wwerr_mumu[4] = wwerr_ee[4] = 0.06; wwderr_emu[4] = wwderr_mumu[4] = wwderr_ee[4] = -0.06; bgerr_emu[4] = bgerr_mumu[4] = bgerr_ee[4] = 0.04; bgderr_emu[4] = bgderr_mumu[4] = bgderr_ee[4] = -0.04; //sname = new char[strlen("CDFJES")+1]; //strcpy(sname,"CDFJES"); //chanresult->sysname.push_back(sname); ername[5] = "CDFJES"; serr_emu[5] = serr_mumu[5] = serr_ee[5] = 0.0;//0.01; sderr_emu[5] = sderr_mumu[5] = sderr_ee[5] = 0.0;//-0.01; wwerr_emu[5] = wwerr_mumu[5] = wwerr_ee[5] = 0.0;//0.01; wwderr_emu[5] = wwderr_mumu[5] = wwderr_ee[5] = 0.0;//-0.01; bgerr_emu[5] = bgerr_mumu[5] = bgerr_ee[5] = 0.0;//0.02; // mostly from fakes bgderr_emu[5] = bgderr_mumu[5] = bgderr_ee[5] = 0.0;//-0.02; // mostly from fakes //sname = new char[strlen("CDFTRKISO")+1]; //strcpy(sname,"CDFTRKISO"); //chanresult->sysname.push_back(sname); ername[6] = "CDFTRKISO"; serr_emu[6] = serr_mumu[6] = serr_ee[6] = 0.02; sderr_emu[6] = sderr_mumu[6] = sderr_ee[6] = -0.02; wwerr_emu[6] = wwerr_mumu[6] = wwerr_ee[6] = 0.02; wwderr_emu[6] = wwderr_mumu[6] = wwderr_ee[6] = -0.02; bgerr_emu[6] = bgerr_mumu[6] = bgerr_ee[6] = 0.02; // partly estimated in CDF 7708 bgderr_emu[6] = bgderr_mumu[6] = bgderr_ee[6] = -0.02; // partly estimated in CDF 7708 //sname = new char[strlen("ALPHAS")+1]; //strcpy(sname,"ALPHAS"); //chanresult->sysname.push_back(sname); ername[7] = "ALPHAS"; serr_emu[7] = serr_mumu[7] = serr_ee[7] = 0.03; sderr_emu[7] = sderr_mumu[7] = sderr_ee[7] = -0.03; wwerr_emu[7] = wwerr_mumu[7] = wwerr_ee[7] = 0.0; wwderr_emu[7] = wwderr_mumu[7] = wwderr_ee[7] = 0.0; bgerr_emu[7] = bgerr_mumu[7] = bgerr_ee[7] = 0.0; bgderr_emu[7] = bgderr_mumu[7] = bgderr_ee[7] = 0.0; //sname = new char[strlen("CDFLEPTID")+1]; //strcpy(sname,"CDFLEPTID"); //chanresult->sysname.push_back(sname); ername[8] = "CDFLEPTID"; serr_emu[8] = ((TH1F*)(TomLikeSig.GetAnaHist("RelWeightError_emu")->GetHist()))->GetMean(); serr_mumu[8] = ((TH1F*)(TomLikeSig.GetAnaHist("RelWeightError_mumu")->GetHist()))->GetMean(); serr_ee[8] = ((TH1F*)(TomLikeSig.GetAnaHist("RelWeightError_ee")->GetHist()))->GetMean();//0.02; sderr_emu[8] = -serr_emu[8]; sderr_mumu[8] = -serr_mumu[8]; sderr_ee[8] = -serr_ee[8]; //-0.02; wwerr_emu[8] = ((TH1F*)(TomLikeWW.GetAnaHist("RelWeightError_emu")->GetHist()))->GetMean(); wwerr_mumu[8] = ((TH1F*)(TomLikeWW.GetAnaHist("RelWeightError_mumu")->GetHist()))->GetMean(); wwerr_ee[8] = ((TH1F*)(TomLikeWW.GetAnaHist("RelWeightError_ee")->GetHist()))->GetMean();//0.02; wwderr_emu[8] = -wwerr_emu[8]; wwderr_mumu[8] = -wwerr_mumu[8]; wwderr_ee[8] = -wwerr_ee[8];//-0.02; bgerr_emu[8] = ((TH1F*)(TomLikeBG.GetAnaHist("RelWeightError_emu")->GetHist()))->GetMean(); bgerr_mumu[8] = ((TH1F*)(TomLikeBG.GetAnaHist("RelWeightError_mumu")->GetHist()))->GetMean(); bgerr_ee[8] = ((TH1F*)(TomLikeBG.GetAnaHist("RelWeightError_ee")->GetHist()))->GetMean();//0.02; bgderr_emu[8] = -bgerr_emu[8]; bgderr_mumu[8] = -bgerr_mumu[8]; bgderr_ee[8] = -bgerr_ee[8];//-0.02; std::cout<sysname.push_back(sname); ername[9] = "CDFFAKERATE"; serr_emu[9] = serr_mumu[9] = serr_ee[9] = 0.0; sderr_emu[9] = sderr_mumu[9] = sderr_ee[9] = 0.0; wwerr_emu[9] = wwerr_mumu[9] = wwerr_ee[9] = 0.0; wwderr_emu[9] = wwderr_mumu[9] = wwderr_ee[9] = 0.0; bgerr_emu[9] = bgerr_mumu[9] = bgerr_ee[9] = 0.0; //TD fake in shape hists 0.05; // 50% error on fakes, which are 10% of background bgderr_emu[9] = bgderr_mumu[9] = bgderr_ee[9] = 0.0;//TD fake in shape hists -0.05; // 50% error on fakes, which are 10% of background //sname = new char[strlen("NNLO")+1]; //strcpy(sname,"NNLO"); //chanresult->sysname.push_back(sname); ername[10] = "NNLO"; serr_emu[10] = serr_mumu[10] = serr_ee[10] = 0.1; sderr_emu[10] = sderr_mumu[10] = sderr_ee[10] = -0.1; wwerr_emu[10] = wwerr_mumu[10] = wwerr_ee[10] = 0.0; wwderr_emu[10] = wwderr_mumu[10] = wwderr_ee[10] = 0.0; bgerr_emu[10] = bgerr_mumu[10] = bgerr_ee[10] = 0.0; bgderr_emu[10] = bgderr_mumu[10] = bgderr_ee[10] = 0.0; //sname = new char[strlen("CDFJETVETO")+1]; //strcpy(sname,"CDFJETVETO"); //chanresult->sysname.push_back(sname); ername[11] = "CDFJETVETO"; serr_emu[11] = serr_mumu[11] = serr_ee[11] = 0.0; sderr_emu[11] = sderr_mumu[11] = sderr_ee[11] = 0.0; wwerr_emu[11] = wwerr_mumu[11] = wwerr_ee[11] = 0.06; wwderr_emu[11] = wwderr_mumu[11] = wwderr_ee[11] = -0.06; bgerr_emu[11] = bgerr_mumu[11] = bgerr_ee[11] = 0.0; bgderr_emu[11] = bgderr_mumu[11] = bgderr_ee[11] = 0.0; //sname = new char[strlen("GENERATOR")+1]; //strcpy(sname,"GENERATOR"); //chanresult->sysname.push_back(sname); ername[12] = "GENERATOR"; serr_emu[12] = serr_mumu[12] = serr_ee[12] = 0.0; sderr_emu[12] = sderr_mumu[12] = sderr_ee[12] = 0.0; wwerr_emu[12] = wwerr_mumu[12] = wwerr_ee[12] = 0.04; wwderr_emu[12] = wwderr_mumu[12] = wwderr_ee[12] = -0.04; bgerr_emu[12] = bgerr_mumu[12] = bgerr_ee[12] = 0.06; bgderr_emu[12] = bgderr_mumu[12] = bgderr_ee[12] = -0.06; //sname = new char[strlen("CDFDYMET")+1]; //strcpy(sname,"CDFDYMET"); //chanresult->sysname.push_back(sname); ername[13] = "CDFDYMET"; serr_emu[13] = serr_mumu[13] = serr_ee[13] = 0.0; sderr_emu[13] = sderr_mumu[13] = sderr_ee[13] = 0.0; wwerr_emu[13] = wwerr_mumu[13] = wwerr_ee[13] = 0.0; wwderr_emu[13] = wwderr_mumu[13] = wwderr_ee[13] = 0.0; bgerr_emu[13] = bgerr_mumu[13] = bgerr_ee[13] = 0.05; bgderr_emu[13] = bgderr_mumu[13] = bgderr_ee[13] = -0.05; ername[14] = "LES";//TD only shape errors (below) serr_emu[14] = serr_mumu[14] = serr_ee[14] = 0.0; sderr_emu[14] = sderr_mumu[14] = sderr_ee[14] = 0.0; wwerr_emu[14] = wwerr_mumu[14] = wwerr_ee[14] = 0.0; wwderr_emu[14] = wwderr_mumu[14] = wwderr_ee[14] = 0.0; bgerr_emu[14] = bgerr_mumu[14] = bgerr_ee[14] = 0.00; bgderr_emu[14] = bgderr_mumu[14] = bgderr_ee[14] = 0.00; TH1 * sigupshape_ee[15]; TH1 * sigdownshape_ee[15]; double sigupshapesigma_ee[15]; double sigdownshapesigma_ee[15]; TH1 * wwupshape_ee[15]; TH1 * wwdownshape_ee[15]; double wwupshapesigma_ee[15]; double wwdownshapesigma_ee[15]; TH1 * bgupshape_ee[15]; TH1 * bgdownshape_ee[15]; double bgupshapesigma_ee[15]; double bgdownshapesigma_ee[15]; TH1 * sigupshape_emu[15]; TH1 * sigdownshape_emu[15]; double sigupshapesigma_emu[15]; double sigdownshapesigma_emu[15]; TH1 * wwupshape_emu[15]; TH1 * wwdownshape_emu[15]; double wwupshapesigma_emu[15]; double wwdownshapesigma_emu[15]; TH1 * bgupshape_emu[15]; TH1 * bgdownshape_emu[15]; double bgupshapesigma_emu[15]; double bgdownshapesigma_emu[15]; TH1 * sigupshape_mumu[15]; TH1 * sigdownshape_mumu[15]; double sigupshapesigma_mumu[15]; double sigdownshapesigma_mumu[15]; TH1 * wwupshape_mumu[15]; TH1 * wwdownshape_mumu[15]; double wwupshapesigma_mumu[15]; double wwdownshapesigma_mumu[15]; TH1 * bgupshape_mumu[15]; TH1 * bgdownshape_mumu[15]; double bgupshapesigma_mumu[15]; double bgdownshapesigma_mumu[15]; for (int i=0;i<15;i++) //TD omg so many empty histograms { sigupshape_ee[i] = 0; sigdownshape_ee[i] = 0; sigupshapesigma_ee[i] = 0; sigdownshapesigma_ee[i] = 0; wwupshape_ee[i] = 0; wwdownshape_ee[i] = 0; wwupshapesigma_ee[i] = 0; wwdownshapesigma_ee[i] = 0; bgupshape_ee[i] = 0; bgdownshape_ee[i] = 0; bgupshapesigma_ee[i] = 0; bgdownshapesigma_ee[i] = 0; sigupshape_emu[i] = 0; sigdownshape_emu[i] = 0; sigupshapesigma_emu[i] = 0; sigdownshapesigma_emu[i] = 0; wwupshape_emu[i] = 0; wwdownshape_emu[i] = 0; wwupshapesigma_emu[i] = 0; wwdownshapesigma_emu[i] = 0; bgupshape_emu[i] = 0; bgdownshape_emu[i] = 0; bgupshapesigma_emu[i] = 0; bgdownshapesigma_emu[i] = 0; sigupshape_mumu[i] = 0; sigdownshape_mumu[i] = 0; sigupshapesigma_mumu[i] = 0; sigdownshapesigma_mumu[i] = 0; wwupshape_mumu[i] = 0; wwdownshape_mumu[i] = 0; wwupshapesigma_mumu[i] = 0; wwdownshapesigma_mumu[i] = 0; bgupshape_mumu[i] = 0; bgdownshape_mumu[i] = 0; bgupshapesigma_mumu[i] = 0; bgdownshapesigma_mumu[i] = 0; } //JES shape errors sigupshape_ee[5] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_JESup")->GetHist()); sigdownshape_ee[5] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_JESdn")->GetHist()); sigupshapesigma_ee[5] = 1; sigdownshapesigma_ee[5] = -1; wwupshape_ee[5] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_JESup")->GetHist()); wwdownshape_ee[5] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_JESdn")->GetHist()); wwupshapesigma_ee[5] = 1; wwdownshapesigma_ee[5] = -1; bgupshape_ee[5] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_JESup")->GetHist()); bgdownshape_ee[5] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_JESdn")->GetHist()); bgupshapesigma_ee[5] = 1; bgdownshapesigma_ee[5] = -1; sigupshape_emu[5] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_JESup")->GetHist()); sigdownshape_emu[5] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_JESdn")->GetHist()); sigupshapesigma_emu[5] = 1; sigdownshapesigma_emu[5] = -1; wwupshape_emu[5] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_JESup")->GetHist()); wwdownshape_emu[5] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_JESdn")->GetHist()); wwupshapesigma_emu[5] = 1; wwdownshapesigma_emu[5] = -1; bgupshape_emu[5] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_JESup")->GetHist()); bgdownshape_emu[5] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_JESdn")->GetHist()); bgupshapesigma_emu[5] = 1; bgdownshapesigma_emu[5] = -1; sigupshape_mumu[5] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_JESup")->GetHist()); sigdownshape_mumu[5] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_JESdn")->GetHist()); sigupshapesigma_mumu[5] = 1; sigdownshapesigma_mumu[5] = -1; wwupshape_mumu[5] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_JESup")->GetHist()); wwdownshape_mumu[5] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_JESdn")->GetHist()); wwupshapesigma_mumu[5] = 1; wwdownshapesigma_mumu[5] = -1; bgupshape_mumu[5] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_JESup")->GetHist()); bgdownshape_mumu[5] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_JESdn")->GetHist()); bgupshapesigma_mumu[5] = 1; bgdownshapesigma_mumu[5] = -1; //TD Lepton Energy Scale Shape Errors sigupshape_ee[14] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_LESup")->GetHist()); sigdownshape_ee[14] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_LESdn")->GetHist()); sigupshapesigma_ee[14] = 1; sigdownshapesigma_ee[14] = -1; wwupshape_ee[14] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_LESup")->GetHist()); wwdownshape_ee[14] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_LESdn")->GetHist()); wwupshapesigma_ee[14] = 1; wwdownshapesigma_ee[14] = -1; bgupshape_ee[14] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_LESup")->GetHist()); bgdownshape_ee[14] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_LESdn")->GetHist()); bgupshapesigma_ee[14] = 1; bgdownshapesigma_ee[14] = -1; sigupshape_emu[14] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_LESup")->GetHist()); sigdownshape_emu[14] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_LESdn")->GetHist()); sigupshapesigma_emu[14] = 1; sigdownshapesigma_emu[14] = -1; wwupshape_emu[14] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_LESup")->GetHist()); wwdownshape_emu[14] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_LESdn")->GetHist()); wwupshapesigma_emu[14] = 1; wwdownshapesigma_emu[14] = -1; bgupshape_emu[14] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_LESup")->GetHist()); bgdownshape_emu[14] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_LESdn")->GetHist()); bgupshapesigma_emu[14] = 1; bgdownshapesigma_emu[14] = -1; sigupshape_mumu[14] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_LESup")->GetHist()); sigdownshape_mumu[14] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_LESdn")->GetHist()); sigupshapesigma_mumu[14] = 1; sigdownshapesigma_mumu[14] = -1; wwupshape_mumu[14] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_LESup")->GetHist()); wwdownshape_mumu[14] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_LESdn")->GetHist()); wwupshapesigma_mumu[14] = 1; wwdownshapesigma_mumu[14] = -1; bgupshape_mumu[14] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_LESup")->GetHist()); bgdownshape_mumu[14] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_LESdn")->GetHist()); bgupshapesigma_mumu[14] = 1; bgdownshapesigma_mumu[14] = -1; //TD PDF Shape errors /*sigupshape_ee[4] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_PDFup")->GetHist()); sigdownshape_ee[4] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_PDFdn")->GetHist()); sigupshapesigma_ee[4] = 1; sigdownshapesigma_ee[4] = -1; wwupshape_ee[4] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_PDFup")->GetHist()); wwdownshape_ee[4] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_PDFdn")->GetHist()); wwupshapesigma_ee[4] = 1; wwdownshapesigma_ee[4] = -1; bgupshape_ee[4] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_PDFup")->GetHist()); bgdownshape_ee[4] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_PDFdn")->GetHist()); bgupshapesigma_ee[4] = 1; bgdownshapesigma_ee[4] = -1; sigupshape_emu[4] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_PDFup")->GetHist()); sigdownshape_emu[4] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_PDFdn")->GetHist()); sigupshapesigma_emu[4] = 1; sigdownshapesigma_emu[4] = -1; wwupshape_emu[4] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_PDFup")->GetHist()); wwdownshape_emu[4] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_PDFdn")->GetHist()); wwupshapesigma_emu[4] = 1; wwdownshapesigma_emu[4] = -1; bgupshape_emu[4] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_PDFup")->GetHist()); bgdownshape_emu[4] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_PDFdn")->GetHist()); bgupshapesigma_emu[4] = 1; bgdownshapesigma_emu[4] = -1; sigupshape_mumu[4] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_PDFup")->GetHist()); sigdownshape_mumu[4] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_PDFdn")->GetHist()); sigupshapesigma_mumu[4] = 1; sigdownshapesigma_mumu[4] = -1; wwupshape_mumu[4] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_PDFup")->GetHist()); wwdownshape_mumu[4] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_PDFdn")->GetHist()); wwupshapesigma_mumu[4] = 1; wwdownshapesigma_mumu[4] = -1; bgupshape_mumu[4] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_PDFup")->GetHist()); bgdownshape_mumu[4] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_PDFdn")->GetHist()); bgupshapesigma_mumu[4] = 1; bgdownshapesigma_mumu[4] = -1;*/ //TD ISR Shape errors sigupshape_ee[3] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_ISRup")->GetHist()); sigdownshape_ee[3] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_ISRdn")->GetHist()); sigupshapesigma_ee[3] = 1; sigdownshapesigma_ee[3] = -1; wwupshape_ee[3] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_ISRup")->GetHist()); wwdownshape_ee[3] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_ISRdn")->GetHist()); wwupshapesigma_ee[3] = 1; wwdownshapesigma_ee[3] = -1; bgupshape_ee[3] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_ISRup")->GetHist()); bgdownshape_ee[3] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_ISRdn")->GetHist()); bgupshapesigma_ee[3] = 1; bgdownshapesigma_ee[3] = -1; sigupshape_emu[3] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_ISRup")->GetHist()); sigdownshape_emu[3] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_ISRdn")->GetHist()); sigupshapesigma_emu[3] = 1; sigdownshapesigma_emu[3] = -1; wwupshape_emu[3] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_ISRup")->GetHist()); wwdownshape_emu[3] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_ISRdn")->GetHist()); wwupshapesigma_emu[3] = 1; wwdownshapesigma_emu[3] = -1; bgupshape_emu[3] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_ISRup")->GetHist()); bgdownshape_emu[3] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_ISRdn")->GetHist()); bgupshapesigma_emu[3] = 1; bgdownshapesigma_emu[3] = -1; sigupshape_mumu[3] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_ISRup")->GetHist()); sigdownshape_mumu[3] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_ISRdn")->GetHist()); sigupshapesigma_mumu[3] = 1; sigdownshapesigma_mumu[3] = -1; wwupshape_mumu[3] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_ISRup")->GetHist()); wwdownshape_mumu[3] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_ISRdn")->GetHist()); wwupshapesigma_mumu[3] = 1; wwdownshapesigma_mumu[3] = -1; bgupshape_mumu[3] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_ISRup")->GetHist()); bgdownshape_mumu[3] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_ISRdn")->GetHist()); bgupshapesigma_mumu[3] = 1; bgdownshapesigma_mumu[3] = -1; //TD Fake shape errors (there mst be a way to handle these properly) sigupshape_ee[9] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_Fakeup")->GetHist()); sigdownshape_ee[9] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiee_Fakedn")->GetHist()); sigupshapesigma_ee[9] = 1; sigdownshapesigma_ee[9] = -1; wwupshape_ee[9] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_Fakeup")->GetHist()); wwdownshape_ee[9] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiee_Fakedn")->GetHist()); wwupshapesigma_ee[9] = 1; wwdownshapesigma_ee[9] = -1; bgupshape_ee[9] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_Fakeup")->GetHist()); bgdownshape_ee[9] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiee_Fakedn")->GetHist()); bgupshapesigma_ee[9] = 1; bgdownshapesigma_ee[9] = -1; sigupshape_emu[9] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_Fakeup")->GetHist()); sigdownshape_emu[9] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphiemu_Fakedn")->GetHist()); sigupshapesigma_emu[9] = 1; sigdownshapesigma_emu[9] = -1; wwupshape_emu[9] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_Fakeup")->GetHist()); wwdownshape_emu[9] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphiemu_Fakedn")->GetHist()); wwupshapesigma_emu[9] = 1; wwdownshapesigma_emu[9] = -1; bgupshape_emu[9] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_Fakeup")->GetHist()); bgdownshape_emu[9] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphiemu_Fakedn")->GetHist()); bgupshapesigma_emu[9] = 1; bgdownshapesigma_emu[9] = -1; sigupshape_mumu[9] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_Fakeup")->GetHist()); sigdownshape_mumu[9] = (TH1F*)(TomLikeSig.GetAnaHist("dilDphimumu_Fakedn")->GetHist()); sigupshapesigma_mumu[9] = 1; sigdownshapesigma_mumu[9] = -1; wwupshape_mumu[9] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_Fakeup")->GetHist()); wwdownshape_mumu[9] = (TH1F*)(TomLikeWW.GetAnaHist("dilDphimumu_Fakedn")->GetHist()); wwupshapesigma_mumu[9] = 1; wwdownshapesigma_mumu[9] = -1; bgupshape_mumu[9] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_Fakeup")->GetHist()); bgdownshape_mumu[9] = (TH1F*)(TomLikeBG.GetAnaHist("dilDphimumu_Fakedn")->GetHist()); bgupshapesigma_mumu[9] = 1; bgdownshapesigma_mumu[9] = -1; //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! no bnackground shapes yet !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! //ee channel nullhyp->add_template(bkgree,1.0,15,ername,bgderr_ee,bgerr_ee, bgdownshape_ee,bgdownshapesigma_ee,bgupshape_ee,bgupshapesigma_ee,0,0,"Channel ee"); testhyp->add_template(bkgree,1.0,15,ername,bgderr_ee,bgerr_ee, bgdownshape_ee,bgdownshapesigma_ee,bgupshape_ee,bgupshapesigma_ee,0,0,"Channel ee"); nullhyp->add_template(WWbkgree,1.0,15,ername,wwderr_ee,wwerr_ee, wwdownshape_ee,wwdownshapesigma_ee,wwupshape_ee,wwupshapesigma_ee,0,0,"Channel ee"); testhyp->add_template(WWbkgree,1.0,15,ername,wwderr_ee,wwerr_ee, wwdownshape_ee,wwdownshapesigma_ee,wwupshape_ee,wwupshapesigma_ee,0,0,"Channel ee"); testhyp->add_template(sigee,1.0,15,ername,sderr_ee,serr_ee, sigdownshape_ee,sigdownshapesigma_ee,sigupshape_ee,sigupshapesigma_ee,0,1,"Channel ee"); //emu channel nullhyp->add_template(bkgremu,1.0,15,ername,bgderr_emu,bgerr_emu, bgdownshape_emu,bgdownshapesigma_emu,bgupshape_emu,bgupshapesigma_emu,0,0,"Channel emu"); testhyp->add_template(bkgremu,1.0,15,ername,bgderr_emu,bgerr_emu, bgdownshape_emu,bgdownshapesigma_emu,bgupshape_emu,bgupshapesigma_emu,0,0,"Channel emu"); nullhyp->add_template(WWbkgremu,1.0,15,ername,wwderr_emu,wwerr_emu, wwdownshape_emu,wwdownshapesigma_emu,wwupshape_emu,wwupshapesigma_emu,0,0,"Channel emu"); testhyp->add_template(WWbkgremu,1.0,15,ername,wwderr_emu,wwerr_emu, wwdownshape_emu,wwdownshapesigma_emu,wwupshape_emu,wwupshapesigma_emu,0,0,"Channel emu"); testhyp->add_template(sigemu,1.0,15,ername,sderr_emu,serr_emu, sigdownshape_emu,sigdownshapesigma_emu,sigupshape_emu,sigupshapesigma_emu,0,1,"Channel emu"); //mumu channel nullhyp->add_template(bkgrmumu,1.0,15,ername,bgderr_mumu,bgerr_mumu, bgdownshape_mumu,bgdownshapesigma_mumu,bgupshape_mumu,bgupshapesigma_mumu,0,0,"Channel mumu"); testhyp->add_template(bkgrmumu,1.0,15,ername,bgderr_mumu,bgerr_mumu, bgdownshape_mumu,bgdownshapesigma_mumu,bgupshape_mumu,bgupshapesigma_mumu,0,0,"Channel mumu"); nullhyp->add_template(WWbkgrmumu,1.0,15,ername,wwderr_mumu,wwerr_mumu, wwdownshape_mumu,wwdownshapesigma_mumu,wwupshape_mumu,wwupshapesigma_mumu,0,0,"Channel mumu"); testhyp->add_template(WWbkgrmumu,1.0,15,ername,wwderr_mumu,wwerr_mumu, wwdownshape_mumu,wwdownshapesigma_mumu,wwupshape_mumu,wwupshapesigma_mumu,0,0,"Channel mumu"); testhyp->add_template(sigmumu,1.0,15,ername,sderr_mumu,serr_mumu, sigdownshape_mumu,sigdownshapesigma_mumu,sigupshape_mumu,sigupshapesigma_mumu,0,1,"Channel mumu"); csm_model* nullhyp_pe = nullhyp->Clone(); csm_model* testhyp_pe = testhyp->Clone(); TCanvas * mycanvas = (TCanvas *) new TCanvas("Canvasemu","Canvasemu"); mycanvas->Divide(3,2); mycanvas->cd(1); testhyp->plotwithdata("Channel ee",dataee ); mycanvas->cd(2); testhyp->plotwithdata("Channel emu",dataemu ); mycanvas->cd(3); testhyp->plotwithdata("Channel mumu",datamumu ); csm* mycsm = new csm(); mycsm->set_htofit( dataee , "Channel ee"); mycsm->set_htofit( dataemu , "Channel emu"); mycsm->set_htofit( datamumu , "Channel mumu"); mycsm->set_modeltofit(testhyp); double chisq = mycsm->chisquared(); std::cout<<"chisq: "<getbestmodel(); mycanvas->cd(4); bestsignalfit->plotwithdata("Channel ee", dataee); mycanvas->cd(5); bestsignalfit->plotwithdata("Channel emu", dataemu); mycanvas->cd(6); bestsignalfit->plotwithdata("Channel mumu", datamumu); mycanvas->Print("results/tomlikeoutemu.ps"); mclimit_csm* mymclimit = new mclimit_csm(); mymclimit->set_null_hypothesis(nullhyp); mymclimit->set_test_hypothesis(testhyp); mymclimit->set_null_hypothesis_pe(nullhyp_pe); mymclimit->set_test_hypothesis_pe(testhyp_pe); mymclimit->set_datahist(dataee, "Channel ee"); mymclimit->set_datahist(dataemu, "Channel emu"); mymclimit->set_datahist(datamumu, "Channel mumu"); mymclimit->set_npe(1000); mymclimit->run_pseudoexperiments(); std::cout << "Getting results" << std::endl; std::cout << "ts: " << mymclimit->ts() << std::endl; std::cout << "tsbm2: " << mymclimit->tsbm2() << std::endl; std::cout << "tsbm1: " << mymclimit->tsbm1() << std::endl; std::cout << "tsbmed: " << mymclimit->tsbmed() << std::endl; std::cout << "tsbp1: " << mymclimit->tsbp1() << std::endl; std::cout << "tsbp2: " << mymclimit->tsbp2() << std::endl; std::cout << "tssm2: " << mymclimit->tssm2() << std::endl; std::cout << "tssm1: " << mymclimit->tssm1() << std::endl; std::cout << "tssmed: " << mymclimit->tssmed() << std::endl; std::cout << "tssp1: " << mymclimit->tssp1() << std::endl; std::cout << "tssp2: " << mymclimit->tssp2() << std::endl; std::cout << "CLs: " << mymclimit->cls() << std::endl; std::cout << "CLb: " << mymclimit->clb() << std::endl; std::cout << "CLsb: " << mymclimit->clsb() << std::endl; //double fit95=mymclimit->s95med(); //std::cout << "Median expected s95 in bg hypothesis: " << fit95 << std::endl; //fitinfo << "Median expected s95 in bg hypothesis: " << fit95 << std::endl; //std::cout<<"deleting Tomlike elements"<bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); // limitcalc->bayes_heinrich_withexpect(0.95,&sflimit,&sfunc, // npxexpect,&sm2,&sm1,&smed,&sp1,&sp2); cout << "tom emu mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << endl; fitinfo << "tom emu mh: " << _higgsMass << " " << sflimit << " " << sm2 << " " << sm1 << " " << smed << " " << sp1 << " " << sp2 << endl; delete mycanvas; /*if(nullhyp!=NULL)delete nullhyp; if(testhyp!=NULL)delete testhyp; if(mycsm!=NULL)delete mycsm; if(bestsignalfit!=NULL)delete bestsignalfit;*/ //TD tom's stuff does not like to be deleted... memory leak! std::cout<<"deleted Tomlike elements"<