File indexing completed on 2023-10-25 09:58:34
0001 #include <iostream>
0002 #include <iomanip>
0003
0004 effCalculator(const unsigned int effType=1, const unsigned int region=0, const unsigned int bothElorPos= 0, const bool isData = true, const bool isMCTruth = false){
0005
0006
0007
0008
0009
0010 std::cout<<"INPUT ARGS effType, region, bothElorPos, isData= "<<effType<< ", "<< region <<", "<< bothElorPos<<", "<< isData << std::endl;
0011
0012 gROOT->ProcessLine(".L ./tdrstyle.C");
0013 setTDRStyle();
0014 tdrStyle->SetErrorX(0.5);
0015 tdrStyle->SetPadLeftMargin(0.19);
0016 tdrStyle->SetPadRightMargin(0.10);
0017 tdrStyle->SetPadBottomMargin(0.15);
0018 tdrStyle->SetLegendBorderSize(0);
0019 tdrStyle->SetTitleYOffset(1.5);
0020 tdrStyle->SetOptTitle(0);
0021
0022 const unsigned int region = region;
0023 const unsigned int bothElorPos = bothElorPos;
0024
0025 using namespace RooFit;
0026
0027 TFile *templates_d = new TFile("templates.root");
0028
0029 const int rebin = 5;
0030
0031 TH1F* SignalTemplate = ( (TH1F* ) templates_d->Get("RelTrkIso_OS_tight") ->Clone());
0032 TH1F* BackgroundTemplate = ( (TH1F* ) templates_d->Get("RelTrkIso_SS_Sideband_loose") ->Clone());
0033
0034 TH1F* DataPass_forTemplateMethod = new TH1F("DataPass", "DataPass", 100/rebin,0.0,1.0);
0035 TH1F* DataFail_forTemplateMethod = new TH1F("DataFail", "DataFail", 100/rebin,0.0,1.0);
0036 SignalTemplate->Rebin(rebin);
0037 BackgroundTemplate->Rebin(rebin);
0038
0039
0040
0041 for(int k=0;k< 100/rebin + 1; k++){
0042 if(SignalTemplate->GetBinLowEdge(k)>= 0.15 ){
0043 SignalTemplate->SetBinContent(k,0);
0044 BackgroundTemplate->SetBinContent(k,0);
0045 }
0046
0047 }
0048
0049
0050 TString foutName = "tpHistos";
0051 if(effType == 1 ){ foutName += "_ID95"; }
0052 else if(effType == 2 ){ foutName += "_ID80"; }
0053 else if(effType == 0 ){ foutName += "_GSF"; }
0054 else if(effType == 3 ){ foutName += "_HLT95"; }
0055 else if(effType == 4 ){ foutName += "_HLT80"; }
0056 TString c1Name = "XXXXXXXXXXX";
0057 TString c2Name = "XXXXXXXXXXX";
0058 if( region == 1){ foutName += "_eb";}
0059 else if( region == 2){ foutName += "_ee";}
0060 if( bothElorPos == 1){ foutName += "_minus";}
0061 else if( bothElorPos == 2 ){ foutName += "_plus";}
0062 c1Name = foutName;
0063 c2Name = foutName;
0064 c1Name += "_pass";
0065 c2Name += "_fail";
0066 if( !isData ){ foutName += "_MONTECARLO"; }
0067 foutName += ".root";
0068
0069 TFile* fout = new TFile(foutName, "RECREATE");
0070 fout->cd();
0071 TH1F* hGSF_pass;
0072 TH1F* hGSF_fail;
0073 TChain* chain;
0074
0075 TString type = "scratch";
0076 TString chainName = "GsfToIso/fitter_tree";
0077
0078 if( !isData && isMCTruth ){
0079 chainName= "GsfToIso/mcUnbias_tree";
0080 }
0081
0082
0083 if( effType == 0 ){
0084 chainName= "PhotonToGsf/fitter_tree";
0085 if( !isData && isMCTruth ){
0086 chainName= "PhotonToGsf/mcUnbias_tree";
0087 }
0088 }
0089
0090 TString hname = type;
0091 std::cout<<hname<<std::endl;
0092 hname+= "_pass";
0093 hGSF_pass = new TH1F(hname, hname, 30, 60, 120);
0094 hname = type;
0095 hname+= "_fail";
0096 hGSF_fail = new TH1F(hname, hname, 12, 60, 120);
0097 std::cout<<"foo"<<std::endl;
0098
0099 chain = new TChain(chainName);
0100 if( isData ){
0101 chain->Add("allTPtrees.root");
0102 }else{ chain->Add("/uscmst1b_scratch/lpc1/old_scratch/lpctrig/jwerner/CMSSW_3_6_1_patch4/src/PhysicsTools/TagAndProbe/test/trees/TPtrees_Zee_All.root"); }
0103 float mass=0;
0104 float probe_et=0;
0105 float probe_eta=0;
0106 float probe_pt=0;
0107 float probe_charge=0;
0108 float probe_combIso=0;
0109 float probe_hcalIso=0;
0110 float probe_ecalIso=0;
0111 float probe_trkIso=0;
0112 int probe_ecalDriven=0;
0113 float probe_missHits=0;
0114 float probe_dist=0;
0115 float probe_dcot=0;
0116 float probe_sigIeta=0;
0117 float probe_deta=0;
0118 float probe_dphi=0;
0119 float probe_hoe=0;
0120
0121 float tag_et=0;
0122 float tag_eta=0;
0123 float tag_pt=0;
0124 float tag_charge=0;
0125 float tag_combIso=0;
0126 float tag_hcalIso=0;
0127 float tag_ecalIso=0;
0128 float tag_trkIso=0;
0129 float tag_dcot=0;
0130 float tag_dist=0;
0131 float tag_missHits=0;
0132 float tag_sigIeta=0;
0133 float tag_deta=0;
0134 float tag_dphi=0;
0135 float tag_hoe=0;
0136
0137 bool passing=false;
0138 bool passingALL=false;
0139 float mass = 0;
0140 float mass_corr = 0;
0141
0142
0143
0144 if( effType == 1 || effType == 3){
0145 chain->SetBranchAddress("probe_passingId", &passing);
0146 }else if( effType == 2 || effType == 4){
0147 chain->SetBranchAddress("probe_passingId80", &passing);
0148 }
0149
0150 if( effType != 0 ){
0151 chain->SetBranchAddress("probe_passingALL", &passingALL);
0152 chain->SetBranchAddress("probe_gsfEle_missingHits", &probe_missHits);
0153 chain->SetBranchAddress("probe_gsfEle_dcot", &probe_dcot);
0154 chain->SetBranchAddress("probe_gsfEle_dist", &probe_dist);
0155 chain->SetBranchAddress("probe_gsfEle_charge", &probe_charge);
0156 chain->SetBranchAddress("probe_gsfEle_hcaliso_dr03", &probe_hcalIso);
0157 chain->SetBranchAddress("probe_gsfEle_ecaliso_dr03", &probe_ecalIso);
0158 chain->SetBranchAddress("probe_gsfEle_trackiso_dr03", &probe_trkIso);
0159 chain->SetBranchAddress("probe_gsfEle_pt", &probe_pt);
0160 chain->SetBranchAddress("probe_sc_et", &probe_et);
0161 chain->SetBranchAddress("probe_sc_eta", &probe_eta);
0162 chain->SetBranchAddress("probe_gsfEle_ecalDrivenSeed", &probe_ecalDriven);
0163 chain->SetBranchAddress("probe_gsfEle_sigmaIetaIeta", &probe_sigIeta);
0164 chain->SetBranchAddress("probe_gsfEle_deltaEta", &probe_deta);
0165 chain->SetBranchAddress("probe_gsfEle_deltaPhi", &probe_dphi);
0166 chain->SetBranchAddress("probe_gsfEle_ecalDrivenSeed", &probe_ecalDriven);
0167 chain->SetBranchAddress("probe_gsfEle_sigmaIetaIeta", &probe_sigIeta);
0168 chain->SetBranchAddress("probe_gsfEle_deltaEta", &probe_deta);
0169 chain->SetBranchAddress("probe_gsfEle_deltaPhi", &probe_dphi);
0170 }
0171 else{
0172 chain->SetBranchAddress("probe_passing", &passing);
0173 chain->SetBranchAddress("probe_et", &probe_pt);
0174 chain->SetBranchAddress("probe_eta", &probe_eta);
0175 chain->SetBranchAddress("probe_hadronicOverEm", &probe_hoe);
0176 chain->SetBranchAddress("probe_hcalTowerSumEtConeDR03", &probe_hcalIso);
0177 chain->SetBranchAddress("probe_ecalRecHitSumEtConeDR03", &probe_ecalIso);
0178 chain->SetBranchAddress("probe_trkSumPtHollowConeDR03", &probe_trkIso);
0179 probe_ecalDriven = 1;
0180 }
0181 if( isData || !isMCTruth ){
0182 chain->SetBranchAddress("tag_sc_et", &tag_et);
0183 chain->SetBranchAddress("tag_sc_eta", &tag_eta);
0184 chain->SetBranchAddress("tag_gsfEle_charge", &tag_charge);
0185 chain->SetBranchAddress("tag_gsfEle_hcaliso_dr03", &tag_hcalIso);
0186 chain->SetBranchAddress("tag_gsfEle_ecaliso_dr03", &tag_ecalIso);
0187 chain->SetBranchAddress("tag_gsfEle_trackiso_dr03", &tag_trkIso);
0188 chain->SetBranchAddress("tag_gsfEle_pt", &tag_pt);
0189 chain->SetBranchAddress("tag_gsfEle_dcot", &tag_dcot);
0190 chain->SetBranchAddress("tag_gsfEle_dist", &tag_dist);
0191 chain->SetBranchAddress("tag_gsfEle_missingHits", &tag_missHits);
0192 chain->SetBranchAddress("tag_gsfEle_sigmaIetaIeta", &tag_sigIeta);
0193 chain->SetBranchAddress("tag_gsfEle_deltaEta", &tag_deta);
0194 chain->SetBranchAddress("tag_gsfEle_deltaPhi", &tag_dphi);
0195 chain->SetBranchAddress("tag_gsfEle_HoverE", &tag_hoe);
0196 chain->SetBranchAddress("mass", &mass);
0197 }
0198
0199
0200 int entries= chain->GetEntries();
0201
0202
0203 unsigned int nOSPass = 0;
0204 unsigned int nSSPass = 0;
0205 unsigned int nOSFail = 0;
0206 unsigned int nSSFail = 0;
0207
0208 RooRealVar x("x","M_{ee}",60,120, "GeV") ;
0209 RooDataSet d_pass("d_pass","d_pass",RooArgSet(x));
0210 RooDataSet d_fail("d_fail","d_fail",RooArgSet(x));
0211
0212 std::cout<<"entries= "<<entries<<std::endl;
0213 for(int i=0; i< entries; i++){
0214 chain->GetEntry(i);
0215
0216 if(effType == 0 ){ probe_et = probe_pt; }
0217
0218 if( region == 1 && fabs(probe_eta) > 1.4442 ){continue;}
0219 else if( region == 2 && fabs(probe_eta) < 1.566 ){continue;}
0220
0221
0222 if( effType!= 0 && bothElorPos == 1 && probe_charge > 0 ){ continue;}
0223 else if( effType!= 0 && bothElorPos == 2 && probe_charge < 0 ){ continue;}
0224
0225
0226 if( (isData || !isMCTruth) && effType== 0 && bothElorPos == 1 && tag_charge < 0 ){ continue;}
0227 else if( (isData || !isMCTruth) && effType== 0 && bothElorPos == 2 && tag_charge > 0 ){ continue;}
0228
0229
0230
0231
0232
0233 float ebScaleShift = 1.0145;
0234 float eeScaleShift = 1.0332;
0235
0236 if( !isData ){
0237 ebScaleShift = 1.0;
0238 eeScaleShift = 1.0;
0239 }
0240
0241
0242 if( fabs(probe_eta) < 1.4442 ){
0243 probe_trkIso = probe_trkIso/ebScaleShift;
0244 probe_ecalIso = probe_ecalIso/ebScaleShift;
0245 probe_hcalIso = probe_hcalIso/ebScaleShift;
0246 }else{
0247 probe_trkIso = probe_trkIso/eeScaleShift;
0248 probe_ecalIso = probe_ecalIso/eeScaleShift;
0249 probe_hcalIso = probe_hcalIso/eeScaleShift;
0250 }
0251 if( fabs(tag_eta) < 1.4442 ){
0252 tag_trkIso = tag_trkIso/ebScaleShift;
0253 tag_ecalIso = tag_ecalIso/ebScaleShift;
0254 tag_hcalIso = tag_hcalIso/ebScaleShift;
0255 }else{
0256 tag_trkIso = tag_trkIso/eeScaleShift;
0257 tag_ecalIso = tag_ecalIso/eeScaleShift;
0258 tag_hcalIso = tag_hcalIso/eeScaleShift;
0259 }
0260
0261
0262 if( probe_ecalDriven< 1){ continue; }
0263 if( fabs(probe_eta) < 1.4442 && probe_et*ebScaleShift < 20){ continue;}
0264 else if( (isData || !isMCTruth) && fabs(tag_eta) < 1.4442 && tag_et*ebScaleShift < 20){ continue;}
0265 else if( fabs(probe_eta) > 1.566 && probe_et*eeScaleShift < 20){ continue;}
0266 else if( (isData || !isMCTruth) && fabs(tag_eta) > 1.566 && tag_et*eeScaleShift < 20){ continue;}
0267
0268 if( (isData || !isMCTruth) && fabs(probe_eta) < 1.4442 && fabs(tag_eta) < 1.4442 ){ mass_corr = ebScaleShift*mass; }
0269 else if( (isData || !isMCTruth) && fabs(probe_eta) > 1.566 && fabs(tag_eta) > 1.566 ){ mass_corr = eeScaleShift*mass; }
0270 else if( (isData || !isMCTruth) ){ mass_corr = sqrt(ebScaleShift*eeScaleShift)*mass; }
0271 else{ mass = 90; mass_corr = mass; }
0272
0273 if( mass_corr > 120.0 || mass_corr < 60.0 ){ continue; }
0274
0275 tag_combIso = (tag_trkIso + max(0, tag_ecalIso - 1.0) + tag_hcalIso ) / tag_pt;
0276 if( fabs(tag_eta) > 1.566){
0277 tag_combIso = (tag_trkIso + tag_ecalIso + tag_hcalIso ) / tag_pt;
0278 }
0279
0280 probe_combIso = (probe_trkIso + max(0, probe_ecalIso - 1.0) + probe_hcalIso ) / probe_pt;
0281 if( fabs(probe_eta) > 1.566){
0282 probe_combIso = (probe_trkIso + probe_ecalIso + probe_hcalIso ) / probe_pt;
0283 }
0284
0285
0286
0287 if( (isData || !isMCTruth) && fabs(tag_eta) < 1.4442 && ( tag_missHits > 0 || ( fabs(tag_dist) < 0.02 && fabs(tag_dcot) < 0.02 ) ||
0288 tag_trkIso/tag_pt > 0.09 || tag_ecalIso/tag_pt > 0.07 || tag_hcalIso/tag_pt > 0.1 ||
0289 tag_sigIeta > 0.01 || fabs(tag_dphi) > 0.06 || fabs(tag_deta) > 0.004 || tag_hoe > 0.04 ) ){ continue;}
0290 if( (isData || !isMCTruth) && fabs(tag_eta) > 1.566 && ( tag_missHits > 0 || ( fabs(tag_dist) < 0.02 && fabs(tag_dcot) < 0.02 ) ||
0291 tag_trkIso/tag_pt > 0.04 || tag_ecalIso/tag_pt > 0.05 || tag_hcalIso/tag_pt > 0.025 ||
0292 tag_sigIeta > 0.03 || fabs(tag_dphi) > 0.03 || tag_hoe > 0.025 ) ){ continue;}
0293
0294
0295 if( effType == 1 || effType == 3 ){
0296
0297 if( fabs(probe_eta) < 1.4442 && ( probe_missHits > 1 ||
0298 probe_trkIso/probe_pt > 0.15 || probe_ecalIso/probe_pt > 2.0 || probe_hcalIso/probe_pt > 0.12 ||
0299 probe_sigIeta > 0.01 || fabs(probe_dphi) > 0.8 || fabs(probe_deta) > 0.007 || probe_hoe > 0.15 ) ){ passing= false;}
0300 if( fabs(probe_eta) > 1.566 && ( probe_missHits > 1 ||
0301 probe_trkIso/probe_pt > 0.08 || probe_ecalIso/probe_pt > 0.06 || probe_hcalIso/probe_pt > 0.05 ||
0302 probe_sigIeta > 0.03 || fabs(probe_dphi) > 0.7 || probe_hoe > 0.07 ) ){ passing= false;}
0303 }
0304
0305 if( effType == 2 || effType == 4 ){
0306
0307 if( fabs(probe_eta) < 1.4442 && ( probe_missHits > 0 || ( fabs(probe_dist) < 0.02 && fabs(probe_dcot) < 0.02 ) ||
0308 probe_trkIso/probe_pt > 0.09 || probe_ecalIso/probe_pt > 0.07 || probe_hcalIso/probe_pt > 0.1 ||
0309 probe_sigIeta > 0.01 || fabs(probe_dphi) > 0.06 || fabs(probe_deta) > 0.004 || probe_hoe > 0.04 ) ){ passing= false;}
0310 if( fabs(probe_eta) > 1.566 && ( probe_missHits > 0 || ( fabs(probe_dist) < 0.02 && fabs(probe_dcot) < 0.02 ) ||
0311 probe_trkIso/probe_pt > 0.04 || probe_ecalIso/probe_pt > 0.05 || probe_hcalIso/probe_pt > 0.025 ||
0312 probe_sigIeta > 0.03 || fabs(probe_dphi) > 0.03 || probe_hoe > 0.025 ) ){ passing= false;}
0313 }
0314
0315
0316 if( (effType >= 3 && passing && probe_hoe<0.15 && passingALL) || (effType < 3 && passing && probe_hoe<0.15) ){
0317 hGSF_pass->Fill(mass_corr);
0318 if( probe_charge*tag_charge > 0 ){ nSSPass++;}else{ nOSPass++;}
0319 DataPass_forTemplateMethod->Fill(tag_trkIso/tag_pt);
0320 x = mass_corr;
0321 d_pass.add( RooArgSet(x) );
0322
0323 }
0324 else if( (effType >= 3 && passing && probe_hoe<0.15 && !passingALL) || (effType < 3 && !passing && probe_hoe<0.15) ){
0325 hGSF_fail->Fill(mass_corr);
0326 if( probe_charge*tag_charge > 0 ){ nSSFail++;}else{ nOSFail++;}
0327 DataFail_forTemplateMethod->Fill(tag_trkIso/tag_pt);
0328 x = mass_corr;
0329 d_fail.add( RooArgSet(x) );
0330
0331 }
0332 }
0333
0334
0335 TFile* file_gen = new TFile("ZeeGenLevel.root", "READ");
0336 TH1* gen_hist = (TH1D*) file_gen->Get("Mass");
0337 gen_hist->Rebin(5);
0338 RooDataHist* rdh = new RooDataHist("rdh","", x, gen_hist);
0339 RooHistPdf* rdh_what = new RooHistPdf("rdh_what", "", RooArgSet(x), (const RooDataHist) rdh );
0340
0341 RooRealVar p1("p1","cb p1", 2.0946e-01 ) ;
0342 RooRealVar p2("p2","cb p2", 8.5695e-04 );
0343 RooRealVar p3("p3","cb p3", 3.8296e-04 );
0344 RooRealVar p4("p4","cp p4", 6.7489e+00 );
0345 RooRealVar sigma("sigma","width of gaussian", 2.5849e+00);
0346 RooRealVar frac("frac","fraction of gauss2", 6.5704e-01 , 0, 1);
0347 RooCBExGaussShape cbX("cbX", "crystal ball X", x, p1, p2, p3, p4, sigma, frac);
0348
0349 p1.setConstant(kFALSE);
0350 p2.setConstant(kFALSE);
0351 p3.setConstant(kTRUE);
0352 p4.setConstant(kTRUE);
0353 sigma.setConstant(kTRUE);
0354 frac.setConstant(kTRUE);
0355 RooFFTConvPdf* signalShapePdf = new RooFFTConvPdf("signalShapePdf","signalShapePdf", x, (const RooAbsPdf) rdh_what, cbX);
0356 float p2_init = p2.getVal();
0357 float p3_init = p3.getVal();
0358 RooRealVar alpha("alpha","alpha", 62.0,50,70 );
0359 RooRealVar beta("beta","beta", 0.001, 0.0, 0.1 );
0360 RooRealVar peak("peak","peak", 91.1876);
0361 RooRealVar gamma("gamma","gamma", 0.05, -0.1, 1.0 );
0362 CMSBkgLineShape bkgPdf("bkgPdf", "bkgPdf", x, alpha, beta, peak, gamma);
0363 alpha.setConstant(kTRUE);
0364 beta.setConstant(kTRUE);
0365 gamma.setConstant(kTRUE);
0366 peak.setConstant(kTRUE);
0367 float gamma_init = gamma.getVal();
0368 RooRealVar fracX("fracX","fraction of cbX", 1.0 , 0.0, 1.0);
0369 fracX.setConstant(kTRUE);
0370
0371 if(region == 2){
0372 gamma.setVal(1.4173e-02);
0373 p1.setVal(8.2633e-01);
0374 p1.setConstant(kTRUE);
0375 p2.setVal(8.5383e-04);
0376 p3.setVal(3.3236e-04);
0377 p4.setVal(5.8295e+01);
0378 }
0379
0380 if( effType== 0 ){
0381 gamma.setConstant(kFALSE);
0382
0383 fracX.setConstant(kFALSE);
0384 p2.setConstant(kTRUE);
0385 }
0386
0387 if((effType == 0 || effType == 1 || effType >= 3) && region != 0 && bothElorPos > 0 ){
0388 gamma.setRange(-10.0, 10.0);
0389 p1.setConstant(kTRUE);
0390 p2.setConstant(kFALSE);
0391 if( region == 2 ){
0392 fracX.setRange(-1.0,2.0);
0393 gamma.setConstant(kTRUE);
0394
0395
0396 }
0397 }
0398
0399 if( effType == 1 && region == 2 && bothElorPos == 1 ){
0400 p2.setConstant(kTRUE);
0401 p1.setConstant(kFALSE);
0402 }
0403
0404
0405 if( !isData ){
0406 p1.setConstant(kTRUE);
0407 gamma.setConstant(kTRUE);
0408 fracX.setConstant(kTRUE);
0409 }
0410
0411 RooAddPdf sum("sum", "sum", RooArgList(*signalShapePdf, bkgPdf), RooArgList(fracX) );
0412 RooDataHist* data_pass = new RooDataHist("data_pass","data_pass", RooArgList(x), hGSF_pass);
0413 RooDataHist* data_fail = new RooDataHist("data_fail","data_fail", RooArgList(x), hGSF_fail);
0414 x.setBins(30);
0415
0416 RooFitResult *fitResult_pass = sum.fitTo(d_pass, RooFit::Save(true),RooFit::Extended(false), RooFit::PrintLevel(1),RooFit::Minos(kTRUE));
0417
0418 RooPlot* frame1 = x.frame();
0419 frame1->SetMinimum(0);
0420 d_pass->plotOn(frame1,RooFit::(kBlack), RooFit::LineColor(kBlack));
0421 sum.plotOn(frame1, RooFit::(kBlue), RooFit::LineColor(kBlue));
0422 float chi2PerDof_pass = frame1->chiSquare(1);
0423 sum.plotOn(frame1,Components("bkgPdf"), RooFit::(kRed), RooFit::LineColor(kRed));
0424
0425 frame1->Draw("e0");
0426 std::cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! pass: chi2/dof= "<<chi2PerDof_pass<<std::endl;
0427
0428 std::cout<<"COUT PASSING FIT RESULT:"<<std::endl;
0429 fitResult_pass->Print("v");
0430 std::cout<<"DONE COUT'ING PASSING FIT RESULTS:"<<std::endl;
0431
0432 float nPass = hGSF_pass->Integral() * fracX->getVal();
0433 float nBkgPass = hGSF_pass->Integral() *(1.0 - fracX->getVal() );
0434 float nPassErr = hGSF_pass->Integral() * fracX->getError();
0435 if( effType == 1 || effType == 3 ){
0436 nPass = hGSF_pass->Integral() * 0.99;
0437 nBkgPass = hGSF_pass->Integral() *(1.0 - 0.99 );
0438 nPassErr = hGSF_pass->Integral() * 0.014;
0439 }
0440
0441 if( effType == 2 || effType == 4 || !isData ){
0442 nPass = hGSF_pass->Integral();
0443 nBkgPass = 0.0;
0444 nPassErr = 0;
0445 }
0446
0447 fracX.setVal(0.7);
0448 gamma.setVal(2.1027e-02);
0449 p1.setVal(2.8386e+00);
0450 p2.setVal(3.5282e-03);
0451 p3.setVal(4.0762e-04);
0452 p4.setVal(5.9024e+01);
0453
0454 fracX.setConstant(kFALSE);
0455 p2.setConstant(kFALSE);
0456 gamma.setConstant(kFALSE);
0457 p1.setConstant(kTRUE);
0458
0459 if(effType == 2 || effType == 4 ){
0460 p1.setConstant(kFALSE);
0461 }
0462
0463
0464 if ( (effType == 3 || effType == 4) && region!= 2 ){
0465 gamma.setVal(1.84739e-01);
0466 fracX.setRange(-1.0,2.0);
0467 p1.setVal(-4.4394e+00);
0468 p2.setVal(1.8673e-06);
0469 p3.setVal(5.4340e-04);
0470 p4.setVal(2.2096e+01);
0471 gamma.setConstant(kTRUE);
0472 p2.setConstant(kTRUE);
0473 p1.setConstant(kFALSE);
0474 }
0475
0476 if(region == 1){
0477 gamma.setVal(3.0834e-02);
0478 p1.setVal(1.1237e+00);
0479 p2.setVal(9.2403e-04);
0480 p3.setVal(3.6280e-04);
0481 p4.setVal(2.9168e+00);
0482 p1.setConstant(kTRUE);
0483 }
0484
0485
0486 if(region == 2){
0487 gamma.setRange(-1.0,100.0);
0488 gamma.setVal(1.4173e-02);
0489 fracX.setRange(-1.0,2.0);
0490 if( effType >= 3 ){
0491 fracX.setRange(0.0,2.0);
0492 gamma.setVal(50.0);
0493 }
0494
0495 p1.setVal(1.9879e+00);
0496 p2.setVal(2.6257e-03);
0497 p3.setVal(5.2198e-04);
0498 p4.setVal(5.9465e+01);
0499 p1.setConstant(kTRUE);
0500 gamma.setConstant(kTRUE);
0501 }
0502
0503 if(effType == 0 && region == 2 && bothElorPos >= 1 ){
0504 fracX.setVal(1);
0505 fracX.setConstant(kTRUE);
0506 }
0507 if( effType >=3 && region == 2){
0508 fracX.setVal(1);
0509 fracX.setConstant(kTRUE);
0510 }
0511
0512
0513 if( effType >= 3 && region == 1 && bothElorPos == 1 ){
0514 fracX.setVal(1.0);
0515 fracX.setConstant(kTRUE);
0516 p2.setConstant(kFALSE);
0517 }
0518
0519 if( effType == 2 && region == 2 && bothElorPos == 2 ){
0520 fracX.setRange(0,1.2);
0521 gamma.setVal(5.0);
0522 p2.setVal(2.4350e-03);
0523 p2.setConstant(kFALSE);
0524 }
0525
0526
0527 if( !isData ){
0528 gamma.setConstant(kTRUE);
0529 fracX.setConstant(kTRUE);
0530 p2.setConstant(kTRUE);
0531 p1.setConstant(kTRUE);
0532 }
0533
0534
0535
0536
0537 x.setBins(12);
0538
0539 RooFitResult *fitResult_fail = sum.fitTo(d_fail, RooFit::Save(true),RooFit::Extended(false), RooFit::PrintLevel(1),RooFit::Minos(kTRUE));
0540
0541 RooPlot* frame2 = x.frame();
0542 frame2->SetMinimum(0);
0543 d_fail->plotOn(frame2,RooFit::(kBlack), RooFit::LineColor(kBlack));
0544
0545 sum.plotOn(frame2, RooFit::(kBlue), RooFit::LineColor(kBlue));
0546 float chi2PerDof_fail = frame2->chiSquare(3);
0547 sum.plotOn(frame2,Components("bkgPdf"), RooFit::(kRed), RooFit::LineColor(kRed));
0548 std::cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! fail: chi2/dof= "<<chi2PerDof_fail<<std::endl;
0549
0550
0551 std::cout<<"COUT FAILING FIT RESULTS:"<<std::endl;
0552 fitResult_fail->Print("v");
0553 std::cout<<"DONE COUT'ING FAILING FIT RESULTS:"<<std::endl;
0554
0555 float nFail = hGSF_fail->Integral() * fracX->getVal();
0556 float nBkgFail = hGSF_fail->Integral() * (1 - fracX->getVal() );
0557 float nFailErr = hGSF_fail->Integral() * fracX->getError();
0558 if(!isData){
0559 nPass = hGSF_pass->Integral();
0560 nFail = hGSF_fail->Integral();
0561 }
0562 float eff = nPass / (nPass + nFail);
0563 float effErr = (1.0 - eff)*sqrt( (nPassErr/nPass)**2 + (nFailErr/nFail)**2 );
0564
0565 std::cout<< "LINESHAPE EFFICIENCY: "<< eff <<" +/- "<< effErr*sqrt(2.0)<<std::endl;
0566 std::cout<<"LINESHAPE RAW INTEGRAL pass, fail= "<< hGSF_pass->Integral() << std::setw(15) << hGSF_fail->Integral() << std::endl;
0567 std::cout<<"LINESHAPE SIGNAL pass, fail= "<< nPass <<" +/- "<< nPassErr <<std::setw(15)<<nFail <<" +/- "<< nFailErr << std::endl;
0568 std::cout<<"LINESHAPE BKG pass, fail= "<< nBkgPass <<" +/- "<<nPassErr <<", "<<nBkgFail<<" +/- "<< nFailErr << std::endl;
0569
0570
0571 TCanvas* c1 = new TCanvas(c1Name,c1Name,500,500);
0572 c1->cd();
0573 fout->cd();
0574 frame1->Draw("e0");
0575 char temp[100];
0576 sprintf(temp, "#chi^{2}/DOF = %.3f", chi2PerDof_pass);
0577 TPaveText *plotlabel8 = new TPaveText(0.65,0.80,0.85,0.77,"NDC");
0578 plotlabel8->AddText(temp);
0579 plotlabel8->SetTextColor(kBlack);
0580 plotlabel8->SetFillColor(kWhite);
0581 plotlabel8->SetBorderSize(0);
0582 plotlabel8->SetTextAlign(12);
0583 plotlabel8->SetTextSize(0.03);
0584 plotlabel8->Draw();
0585
0586 TCanvas* c2 = new TCanvas(c2Name,c2Name,500,500);
0587 c2->cd();
0588 fout->cd();
0589 frame2->Draw("e0");
0590 TPaveText *plotlabel7 = new TPaveText(0.65,0.82,0.85,0.87,"NDC");
0591 plotlabel7->SetTextColor(kBlack);
0592 plotlabel7->SetFillColor(kWhite);
0593 plotlabel7->SetBorderSize(0);
0594 plotlabel7->SetTextAlign(12);
0595 plotlabel7->SetTextSize(0.03);
0596 TPaveText *plotlabel9 = new TPaveText(0.65,0.80,0.85,0.77,"NDC");
0597 plotlabel9->SetTextColor(kBlack);
0598 plotlabel9->SetFillColor(kWhite);
0599 plotlabel9->SetBorderSize(0);
0600 plotlabel9->SetTextAlign(12);
0601 plotlabel9->SetTextSize(0.03);
0602 sprintf(temp, "#chi^{2}/DOF = %.3f", chi2PerDof_fail);
0603 plotlabel9->AddText(temp);
0604 sprintf(temp, "#epsilon = %.3f #pm %.3f", eff, effErr);
0605 plotlabel7->AddText(temp);
0606 plotlabel9->Draw();
0607 plotlabel7->Draw();
0608 c1->cd();
0609 plotlabel7->Draw();
0610 c2->cd();
0611 plotlabel9->Draw();
0612
0613 c1->Write();
0614 c2->Write();
0615
0616
0617
0618
0619 std::cout<<"nOSPass, nSSPass= "<<nOSPass<<", "<<nSSPass<<std::endl;
0620 std::cout<<"nOSFail, nSSFail= "<<nOSFail<<", "<<nSSFail<<std::endl;
0621
0622
0623
0624
0625
0626
0627
0628
0629 float qMisID= 0.00759057;
0630
0631 if( !isData ){ qMisID= 0.0123; }
0632
0633 float qMisID_fail= qMisID;
0634 float nSigPass_fromQ = (float(nOSPass) - float(nSSPass))/(1 - 2.0*qMisID)**2;
0635 float nBkgPass_fromQ = (float(nOSPass) + float(nSSPass)) - nSigPass_fromQ;
0636
0637
0638 std::cout<<"nOSFail, nSSFail= "<<nOSFail<<", "<<nSSFail<<std::endl;
0639
0640 float nSigFail_fromQ = (float(nOSFail) - float(nSSFail))/(1 - 2.0*qMisID_fail)**2;
0641 if(nSigFail_fromQ < 0 ){ nSigFail_fromQ = 0.0;}
0642 float nBkgFail_fromQ = (float(nOSFail) + float(nSSFail)) - nSigFail_fromQ;
0643
0644
0645 float qMisID_biased= qMisID + 0.00383;
0646 float qMisID_fail_biased= qMisID_fail + 1.0*1.97194399999999978e-02;
0647
0648 float nSigPass_fromQ_biased = (float(nOSPass) - float(nSSPass))/(1 - 2.0*qMisID_biased)**2;
0649 float nBkgPass_fromQ_biased = (float(nOSPass) + float(nSSPass)) - nSigPass_fromQ_biased;
0650
0651 float nSigFail_fromQ_biased = (float(nOSFail) - float(nSSFail))/(1 - 2.0*qMisID_fail_biased)**2;
0652 float nBkgFail_fromQ_biased = (float(nOSFail) + float(nSSFail)) - nSigFail_fromQ_biased;
0653
0654 float nSigPass_fromQ_err = fabs( nSigPass_fromQ - nSigPass_fromQ_biased );
0655 float nSigFail_fromQ_err = fabs( nSigFail_fromQ - nSigFail_fromQ_biased );
0656 float nBkgPass_fromQ_err = fabs( nBkgPass_fromQ - nBkgPass_fromQ_biased );
0657 float nBkgFail_fromQ_err = fabs( nBkgFail_fromQ - nBkgFail_fromQ_biased );
0658
0659
0660
0661
0662
0663 std::cout<<"qMisID EFFICIENCY: "<<nSigPass_fromQ/(nSigPass_fromQ + nSigFail_fromQ)
0664 << " +/- " <<
0665 (1.0/(nSigPass_fromQ + nSigFail_fromQ))*sqrt( nSigPass_fromQ*(1.0 - nSigPass_fromQ/(nSigPass_fromQ + nSigFail_fromQ)) )
0666 << " (stat)+/- "<<
0667 (1.0 - nSigPass_fromQ/(nSigPass_fromQ + nSigFail_fromQ))*sqrt( (nSigPass_fromQ_err/nSigPass_fromQ)**2 + (nSigFail_fromQ_err/((nSigFail_fromQ>0)?nSigFail_fromQ:nSigFail_fromQ_err) )**2 )<< " (syst)"<<std::endl;
0668
0669 std::cout<<"qMisID SIGNAL pass, fail= "<<nSigPass_fromQ<<" +/- " << nSigPass_fromQ_err <<", "<< nSigFail_fromQ<<" +/- " << nSigFail_fromQ_err <<std::endl;
0670 std::cout<<"qMisID BKG pass, fail= "<<nBkgPass_fromQ<<" +/- " << nBkgPass_fromQ_err <<", "<< nBkgFail_fromQ<<" +/- " << nBkgFail_fromQ_err <<std::endl;
0671
0672
0673
0674
0675
0676
0677
0678
0679
0680
0681
0682
0683
0684
0685
0686
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742 fout->Write();
0743 }
0744