Back to home page

Project CMSSW displayed by LXR
 
 

    

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   //efftype = 0 sc->gsf, 1 gsf->wp95, 2 gsf->wp80, 3 wp95->hlt, 4 wp80->hlt
0007   //region = 0 all, 1 EB, 2 EE
0008   //bothElorPos = 0 all, 1 e-, 2 e+
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"); //note we only bother w/ template numbers for data for now at least... not MC
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      //use probe_passing for 95 and probe_passingId80 for 80
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     //if( !isData ){ entries = 100000; }
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       //this is for the non sc->gsf steps
0222       if( effType!= 0 && bothElorPos == 1 && probe_charge > 0 ){ continue;}
0223       else if(  effType!= 0 && bothElorPos == 2 && probe_charge < 0 ){ continue;}
0224 
0225       //this is for the sc->gsf steps
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       //note that for the MCTruth case we lose all charge info at this point...
0229 
0230       //float ebScaleShift = 1.0115 - 0.0;
0231       //float eeScaleShift = 1.0292 - 0.0;
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       //apply scale correction to iso variables
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; } //B+B
0269       else if( (isData || !isMCTruth) && fabs(probe_eta) > 1.566 &&  fabs(tag_eta) > 1.566 ){ mass_corr = eeScaleShift*mass; } //E+E
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       //verify that tag passes the tag criteria (WP80)
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 ){ //WP95 or HLT95
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 ){ //WP80 or HLT80
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       //DataPass_forTemplateMethod->Fill(probe_trkIso/probe_pt);
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       //DataFail_forTemplateMethod->Fill(probe_trkIso/probe_pt);
0331     }
0332     }
0333 
0334     //x.setBins(30);
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); //template method bkg hypothesis for pass pass case
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       //fracX.setVal(0.9);
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         //p1.setConstant(kFALSE);
0395         //p2.setConstant(kTRUE);
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     //RooFitResult *fitResult_pass = sum.fitTo( *data_pass, RooFit::Save(true),RooFit::Extended(false), RooFit::PrintLevel(1),RooFit::Minos(kTRUE));//binned
0416     RooFitResult *fitResult_pass = sum.fitTo(d_pass, RooFit::Save(true),RooFit::Extended(false), RooFit::PrintLevel(1),RooFit::Minos(kTRUE));//unbinned
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; // 1% bkg from fake rate method 
0437       nBkgPass = hGSF_pass->Integral() *(1.0 -  0.99 ); // 1% bkg from template method 
0438       nPassErr = hGSF_pass->Integral() * 0.014; //1.4% error on the bkg from template method
0439     }
0440 
0441     if( effType == 2 || effType == 4 || !isData ){ //no background passes
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); //include bkg for fail case
0455     p2.setConstant(kFALSE);
0456     gamma.setConstant(kFALSE);//TEST                                                                                                                               
0457     p1.setConstant(kTRUE);//TEST 
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       //gamma.setVal(6.3795e+00);                                                                                                                                                                 
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); //not enough data.. don't even bother                                                                                            
0506     }
0507     if( effType >=3 && region == 2){
0508       fracX.setVal(1);
0509       fracX.setConstant(kTRUE); //not enough data.. don't even bother                                                                                            
0510     }
0511 
0512 
0513     if( effType >= 3 && region == 1 && bothElorPos == 1 ){
0514       fracX.setVal(1.0);
0515       fracX.setConstant(kTRUE); //not enough data.. don't even bother
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     //except static for MC
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     //RooFitResult *fitResult_fail = sum.fitTo(*data_fail, RooFit::Save(true),RooFit::Extended(false), RooFit::PrintLevel(1), RooFit::Minos(kTRUE));//binned
0539     RooFitResult *fitResult_fail = sum.fitTo(d_fail, RooFit::Save(true),RooFit::Extended(false), RooFit::PrintLevel(1),RooFit::Minos(kTRUE));//unbinned
0540 
0541     RooPlot* frame2 = x.frame();
0542     frame2->SetMinimum(0);
0543     d_fail->plotOn(frame2,RooFit::(kBlack), RooFit::LineColor(kBlack));
0544     //data_fail->plotOn(frame2,RooFit::(kBlack), RooFit::LineColor(kBlack));
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     //std::cout<<"WARNING: NOTE THAT ANY ***FIT*** PRINTOUT BELOW IS FROM AN ENTIRELY DIFFERENT FIT, NOT RELATED TO THE LINESHAPE EFF FITS"<<std::endl;
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     //std::cout<<"nPassErr, nFailErr= "<< nPassErr <<", "<<nFailErr << std::endl;
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     //from data <q1q2> = -0.965 +/- 0.009
0623     // ---> qMisID = (1 - sqrt(fabs(<q1q2>))/2 = 0.00883 +/- 0.00229
0624     //...alright, but that is just statistical.  adding in both stat and syst errors we get:
0625     //<q1q2> = -0.965 +/- 0.015
0626 
0627 
0628     //float qMisID= 100.0;//
0629     float qMisID= 0.00759057; //for WP70 qmisid
0630     //float qMisID= 0.00883; //for WP85 qmisid
0631     if( !isData ){ qMisID= 0.0123; }
0632     //float qMisID= 0.00883; //for WP80 numbers
0633     float qMisID_fail= qMisID;//here we correct for the pass-->fail bias from MC
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     //float qMisID_biased= 0.00883 + 0.00229; //stat
0645     float qMisID_biased= qMisID + 0.00383; //stat+syst
0646     float qMisID_fail_biased= qMisID_fail + 1.0*1.97194399999999978e-02; //for WP95 numbers
0647     //float qMisID_fail_biased= qMisID + 1.97194399999999978e-02; //for WP80 numbers
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     //scratch the template method for now
0673     /*
0674     TH1F* DataPass_forTemplateMethod_test = ( (TH1F* ) templates_d->Get("RelTrkIso") ->Clone());
0675     DataPass_forTemplateMethod_test->Rebin(rebin);
0676 
0677     DataPass_forTemplateMethod_test->SetLineColor(kBlue);
0678     DataPass_forTemplateMethod_test->DrawNormalized();
0679     DataPass_forTemplateMethod->SetLineColor(kRed);
0680     DataPass_forTemplateMethod->DrawNormalized("same");
0681 
0682     TObjArray *Template = new TObjArray(2);
0683     Template->Add(SignalTemplate);
0684     Template->Add(BackgroundTemplate);
0685     TFractionFitter* fit = new TFractionFitter(DataPass_forTemplateMethod, Template);
0686 
0687     fit->Constrain(1, 0.0, 1.0);
0688     fit->Constrain(2, 0.0, 1.0);
0689 
0690     double theFrac[3], err[3];
0691     Int_t status = fit->Fit();
0692 
0693     cout << "PASS fit status: " << status << endl;
0694     if (status == 0 ) {                       // check on fit status                                                               
0695       TH2F* result = (TH2F*) fit->GetPlot();
0696       result->SetTitle("fit pass");
0697 
0698       for(int i=0; i<2; i++){
0699     fit->GetResult(i, theFrac[i], err[i]);
0700       }
0701 
0702       std::cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DONE FITTING PASSING !!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl;
0703       std::cout<<"PASS: fracS= "<< 1.0 - theFrac[1] <<" +/- "<< err[1] << std::endl;
0704     }
0705 
0706     TObjArray *TemplateF = new TObjArray(2);
0707     TemplateF->Add(SignalTemplate); //or could use SignalTemplate_fail for probes
0708     TemplateF->Add(BackgroundTemplate);//similarly...
0709     TFractionFitter* fitF = new TFractionFitter(DataFail_forTemplateMethod, TemplateF);
0710     fitF->Constrain(1, 0.0, 1.0);
0711     fitF->Constrain(2, 0.0, 1.0);
0712     double theFracF[3], errF[3];
0713     Int_t statusF = fitF->Fit();
0714     if (statusF == 0 ) {
0715       TH2F* resultF = (TH2F*) fitF->GetPlot();
0716       resultF->SetTitle("fit pass");
0717 
0718       for(int i=0; i<2; i++){
0719         fitF->GetResult(i, theFracF[i], errF[i]);
0720       }
0721       std::cout<<"FAIL: fracS= "<< 1.0 - theFracF[1] <<" +/- "<< errF[1] << std::endl;                                                                                                                                                 
0722 
0723       float nSigPass_temp =  DataPass_forTemplateMethod->GetEntries()*(1.0 - theFrac[1])/1.0;
0724       float nBkgPass_temp =  DataPass_forTemplateMethod->GetEntries()*theFrac[1]/1.0;
0725       float nSigFail_temp =  DataFail_forTemplateMethod->GetEntries()*(1.0 - theFracF[1])/1.0;
0726       float nBkgFail_temp =  DataFail_forTemplateMethod->GetEntries()*theFracF[1]/1.0;
0727 
0728       float nSigPass_temp_err =  DataPass_forTemplateMethod->GetEntries()*err[1]/1.0;
0729       float nBkgPass_temp_err =  DataPass_forTemplateMethod->GetEntries()*err[1]/1.0;
0730       float nSigFail_temp_err =  DataFail_forTemplateMethod->GetEntries()*errF[1]/1.0;
0731       float nBkgFail_temp_err =  DataFail_forTemplateMethod->GetEntries()*errF[1]/1.0;
0732 
0733       std::cout<<"template EFFICIENCY: "<<nSigPass_temp/(nSigPass_temp + nSigFail_temp)
0734            << " +/- " << (1.0 - nSigPass_temp/(nSigPass_temp + nSigFail_temp))*sqrt( (nSigPass_temp_err/nSigPass_temp)**2 + (nSigFail_temp_err/nSigFail_temp)**2 )<< std::endl;
0735 
0736       std::cout<<"template SIGNAL pass, fail= "<<nSigPass_temp<<" +/- " << nSigPass_temp_err  <<", "<< nSigFail_temp<<" +/- " << nSigFail_temp_err <<std::endl;
0737       std::cout<<"template BKG pass, fail= "<<nBkgPass_temp<<" +/- " << nBkgPass_temp_err  <<", "<< nBkgFail_temp<<" +/- " << nBkgFail_temp_err <<std::endl;
0738  
0739       }
0740     */
0741 
0742   fout->Write();
0743 }
0744
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!