ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/applySelection.cc
(Generate patch)

Comparing UserCode/MitHzz4l/Selection/src/applySelection.cc (file contents):
Revision 1.23 by khahn, Sun May 20 19:08:59 2012 UTC vs.
Revision 1.37 by dkralph, Wed Jun 13 13:39:12 2012 UTC

# Line 50 | Line 50 | using namespace std;
50   #include "MuonSelection.h"
51   #include "ElectronSelection.h"
52   #include "IsolationSelection.h"
53 <
54 < //#include "Selection.h"
55 < //#include "ReferenceSelection.h"
56 < #include "ReferenceSelection2.h"
53 > #include "ReferenceSelection.h"
54  
55   #include "TriggerUtils.h"
56   #include "PassHLT.h"
57   #include "Angles.h"
58   #include "KinematicsStruct.h"
59   #include "InfoStruct.h"
60 < //#include "GenInfoStruct.h"
60 > #include "GenInfoStruct.h"
61   #include "WeightStruct.h"
65 //#include "GlobalDataAndFuncs.h"
66 #include "RunLumiRangeMap.h"
67
62   #include "AngleTuple.h"
63   #include "FOTuple.h"
64  
65 + #include "RunLumiRangeMap.h"
66   #include "SampleWeight.h"
67   #include "EfficiencyWeightsInterface.h"
68  
# Line 97 | Line 92 | vector<bool>   PFnoPUflag;;
92   //
93   // prototypes
94   //----------------------------------------------------------------------------
100 bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& );
95   void initPUWeights();
96   double getPUWeight(unsigned npu);
97 + void setEra(string, ControlFlags&);
98   void setEffiencyWeights(EventData&, WeightStruct& );
99   void initRunLumiRangeMap();
100   void initEvtRhoMap(map<unsigned,float> &);
101   unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
102                           vector<bool> &pfNoPileUpFlag,
103 <                         mithep::Array<Vertex>      * vtxArr );
103 >                         const mithep::Array<mithep::Vertex>     * vtxArr );
104   //----------------------------------------------------------------------------
105  
106  
# Line 117 | Line 112 | int main(int argc, char** argv)
112   {
113    cutstrs.push_back(string("skim0"));              //0
114    cutstrs.push_back(string("skim1"));              //1
115 +  cutstrs.push_back(string("skim2"));              //1
116    cutstrs.push_back(string("trigger"));            //2
117    // -------------------------------------------------
118    cutstrs.push_back(string("Z candidate"));        //3
# Line 128 | Line 124 | int main(int argc, char** argv)
124    cutstrs.push_back(string("ZZ 20/10"));           //8
125    cutstrs.push_back(string("resonance"));          //9
126    cutstrs.push_back(string("m4l>70"));             //10
127 <  cutstrs.push_back(string("m4l>100"));            //11
127 >  cutstrs.push_back(string("m4l>100, mZ2 > 12"));                  //11
128  
129  
130  
# Line 214 | Line 210 | int main(int argc, char** argv)
210    KinematicsStruct kinematics;
211    InfoStruct evtinfo;
212    WeightStruct weights;
213 <  //  GenInfoStruct geninfo;
213 >  GenInfoStruct geninfo;
214  
215    AngleTuple nt( (const char*)ofname, (const char*)"zznt");
216    nt.makeAngleBranch(angles);
217    nt.makeKinematicsBranch(kinematics);
218    nt.makeInfoBranch(evtinfo);
219  
220 +  
221    FOTuple foTree( nt.getFile(), (const char*)"FOtree");
222    foTree.makeFailedLeptonBranch(failingLeptons);
223 <  foTree.makeZLeptonBranch(passingLeptons);
223 >  foTree.makePassedLeptonBranch(passingLeptons);
224  
225 <  TH1F * h_w_hpt;
225 >  TH1F * h_wf11_hpt;
226    if(ctrl.mc) {
227      nt.makeWeightBranch(weights);
228 <    //    nt.makeGenInfoBranch(geninfo);
228 >    if(ctrl.fillGen ) nt.makeGenInfoBranch(geninfo);
229      initEfficiencyWeights();
230      initPUWeights();
231 <    /*
231 >
232      // Higgs only, pt reweighting
233      if( ctrl.mH <= 140 ) {
234        char wptname[256];
# Line 239 | Line 236 | int main(int argc, char** argv)
236        TFile * f_hpt = new TFile(wptname);
237        f_hpt->Print();
238        sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
239 <      h_w_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
243 <      h_w_hpt->Print();
239 >      h_wf11_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
240      }
241 <    */
241 >
242    } else {
243      initRunLumiRangeMap();
244    }
245  
246    //  initMuonIDMVA();
247 <  initMuonIsoMVA();
252 <  initElectronIDMVA();
253 <  initElectronIsoMVA();
247 >  initElectronIDMVAV1();
248    initTrigger();
249  
250    // hack
# Line 265 | Line 259 | int main(int argc, char** argv)
259    string currentFile("");
260  
261    mithep::EventHeader *info    = new mithep::EventHeader();
268  //  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
262    mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
263    mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
264    mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
# Line 274 | Line 267 | int main(int argc, char** argv)
267    mithep::Array<mithep::PileupEnergyDensity>  *puDArr        = new mithep::Array<mithep::PileupEnergyDensity>();
268    mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
269    mithep::Array<mithep::MCParticle>           *mcArr         = new mithep::Array<mithep::MCParticle>();
270 +  mithep::MCEventInfo                         *mcEvtInfo     = new mithep::MCEventInfo();
271    mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
272    mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
273    vector<string>                              *hltTableStrings  = new vector<string>();
# Line 287 | Line 281 | int main(int argc, char** argv)
281    TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
282    TString fTrackName(Names::gkTrackBrn);
283    TString fMCParticleName(Names::gkMCPartBrn);
284 +  TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
285    TString fTriggerMaskName(Names::gkHltBitBrn);
286    TString fTriggerTableName(Names::gkHltTableBrn);
287    
# Line 298 | Line 293 | int main(int argc, char** argv)
293    if( ctrl.mc ) {
294      chain->SetBranchAddress(fPileupInfoName,  &puArr);
295      chain->SetBranchAddress(fMCParticleName,  &mcArr);
296 +    chain->SetBranchAddress(fMCEvtInfoName,  &mcEvtInfo);
297    }
298    chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
299    chain->SetBranchAddress(fTrackName, &trkArr);
# Line 330 | Line 326 | int main(int argc, char** argv)
326    for(UInt_t ientry=0; ientry<imax; ientry++)
327      {
328        chain->GetEntry(ientry);
333
334      /*
335      bool found_tau=false;
336      if( ctrl.mc) {  // && ctrl.debug
337        cout << "nparticles: " << mcArr->GetEntries() << endl;
338        cout << "-----------------------------------" << endl;
339        for( int i=0; i<mcArr->GetEntries(); i++ ) {
340          const mithep::MCParticle *mc = mcArr->At(i);
341          cout << "i: " << i << "\t"
342               << "id: " << mc->AbsPdgId()
343               << endl;
344          if( mc->Is(MCParticle::kTau) ) {
345            found_tau=true;
346          }
347        }
348      }
349      cout << "-----------------------------------" << endl;
350      cout << "found tau: " << found_tau << endl;
351      */
352
353      // events that fail bc ele overlaps w/ reco mu
354      /*
355      if( info->EvtNum() != 67315 &&
356          info->EvtNum() != 288693 ) continue;
357      */
358
359      // events where they apparently don't veto electrons ...
360      /*
361      if( info->EvtNum() != 141713 &&
362          info->EvtNum() != 178019 &&
363          info->EvtNum() != 60823 &&
364          info->EvtNum() != 192795 ) continue;
365      */
366
329        if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
368      setPV( ctrl, vtxArr, vtx);
369      //if (!(setPV( ctrl, vtxArr, vtx)) ) continue;
330  
331 <      float rho = evtrhoMap[info->EvtNum()];
331 >      string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
332 >      setEra( fname, ctrl );
333 >
334  
335        //
336        // pfNoPU
# Line 383 | Line 345 | int main(int argc, char** argv)
345        //
346        if(ctrl.mc) {
347          //
348 +        // gen info
349 +        //
350 +        if( ctrl.fillGen )
351 +          fillGenInfo( mcArr, mcEvtInfo, geninfo, ESampleType::kHZZ, ctrl);
352 +
353 +        //
354          // xsec & PU weights
355          //
356          weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
# Line 395 | Line 363 | int main(int argc, char** argv)
363          assert(npu>=0);
364          evtinfo.npu;
365          weights.npuw = getPUWeight(evtinfo.npu);
366 <        //      cout << "weight: " << weights.w << endl;
366 >        
367 >        //
368 >        // pt reweighting for Higgs < 140, F11
369 >        //
370 >        if( ctrl.fillGen ) {
371 >          if( (fname.find("f11-h115zz4l") != string::npos) ||
372 >              (fname.find("f11-h120zz4l") != string::npos) ||
373 >              (fname.find("f11-h130zz4l") != string::npos)  ) {
374 >            weights.w *= h_wf11_hpt->GetBinContent(h_wf11_hpt->FindBin(geninfo.pt_zz));
375 >          }
376 >        }
377 >
378          //
379          // trigger
380          //
# Line 404 | Line 383 | int main(int argc, char** argv)
383            hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
384            hltchain->GetEntry(0);
385            hltTable->Clear();
386 <          fillTriggerBits(hltTable, hltTableStrings );
386 >          fillTriggerNames(hltTable, hltTableStrings );
387          }
388          if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
389          fillTriggerBits( hltTable, trigMask, triggerBits );
390          passes_HLT_MC = true; // passed to apply as extern global, just for sync
391 <        if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum() ) ) {
391 >        if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum(), k2012_MC ) ) {
392            passes_HLT_MC = false;
393          }
394        } else {
395          //
396          // JSON
397          //
398 <        /*
399 <        RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
400 <        if( !(rlrm.HasRunLumi(rl)) )  {
401 <          if( ctrl.debug ) cout << "\tfails JSON" << endl;
402 <          continue;
398 >        // if(!(ctrl.noJSON) ) {
399 >        cout << "noJSON flag not in Control Flags" << endl;
400 >        assert(0);
401 >          RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
402 >          if( !(rlrm.HasRunLumi(rl)) )  {
403 >            if( ctrl.debug ) cout << "\tfails JSON" << endl;
404 >            continue;
405 >          // }
406          }
407 <        */
407 >
408          //
409          // trigger
410          //
# Line 431 | Line 413 | int main(int argc, char** argv)
413            hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
414            hltchain->GetEntry(0);
415            hltTable->Clear();
416 <          fillTriggerBits(hltTable, hltTableStrings );
416 >          fillTriggerNames(hltTable, hltTableStrings );
417          }
418          if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
419          /*
# Line 452 | Line 434 | int main(int argc, char** argv)
434  
435  
436          // reference
437 <        /*
438 <        apply_HZZ4L_reference2_selection(ctrl, info,
439 <                                         vtx,
440 <                                         pfArr,
441 <                                         //rho,
442 <                                         puDArr,
443 <                                         electronArr,
444 <                                         &electronReferencePreSelection,
445 <                                         &electronReferenceIDMVASelection,
446 <                                         &electronReferenceIsoSelection,
447 <                                         muonArr,
448 <                                         &muonReferencePreSelection,
449 <                                         &muonIDPFSelection,
468 <                                         &muonReferenceIsoSelection);
469 <        */
470 <
471 <        // evolved
472 <        apply_HZZ4L_reference2_selection(ctrl, info,
473 <                                         vtx,
474 <                                         pfArr,
475 <                                         //rho,
476 <                                         puDArr,
477 <                                         electronArr,
478 <                                         &electronReferencePreSelection,
479 <                                         //&electronPreSelection,
480 <                                         &electronReferenceIDMVASelection,
481 <                                         &electronIsoMVASelection,
482 <                                         muonArr,
483 <                                         &muonReferencePreSelection,
484 <                                         //&muonPreSelection,
485 <                                         &muonIDPFSelection,
486 <                                         &muonIsoMVASelection);
487 <
488 <
437 >        apply_HZZ4L_reference_selection(ctrl, info,
438 >                                        vtxArr,
439 >                                        pfArr,
440 >                                        puDArr,
441 >                                        electronArr,
442 >                                        &electronReferencePreSelection,
443 >                                        &electronReferenceIDMVASelectionV1,
444 >                                        &electronReferenceIsoSelection,
445 >                                        muonArr,
446 >                                        &muonReferencePreSelection,
447 >                                        &muonIDPFSelection,
448 >                                        &muonReferenceIsoSelection);
449 >      
450        if( ctrl.debug ) cout << endl;
451  
452  
# Line 509 | Line 470 | int main(int argc, char** argv)
470        if( ret4l.status.pass() ) {
471          fillAngles(ret4l,angles);
472          fillKinematics(ret4l,kinematics);
473 +        fillMassErrors(ret4l,muonArr,electronArr,kinematics);
474          fillEventInfo(info,evtinfo);
475          if( ctrl.mc) {
514        // fillGenInfo(ginfo,geninfo);
476            setEffiencyWeights(ret4l, weights);
477             if(ctrl.debug)
478               cout << "w: " << weights.w << "\t"
# Line 519 | Line 480 | int main(int argc, char** argv)
480                    << "woff: " << weights.woff << endl;
481          }
482          
483 <        /*
523 <        // only for Higgs < 140
524 <        geninfo.weight *= h_w_hpt->GetBinContent(h_w_hpt->FindBin(geninfo.pt_zz));
525 <        angles.bdt = reader->EvaluateMVA("BDTG");
526 <        */
483 >
484          nt.Fill();
485          
486          cerr  << "PASS:: "
# Line 538 | Line 495 | int main(int argc, char** argv)
495          pass++;
496          //      if( pass > 3 ) break;
497  
498 <      } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) ) { // save for fakes ...
499 <        //      if(ctrl.debug) {
498 >      } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) &&
499 >                 ret4l.status.selectionBits.to_ulong() < (1UL<<PASS_4L)) { // save for fakes ...
500 >        if(ctrl.debug) {
501            cout << "failingLeptons : " << failingLeptons.size() << endl;
502            for( int f=0; f<failingLeptons.size(); f++ ) {
503              SimpleLepton  l = failingLeptons[f];
# Line 548 | Line 506 | int main(int argc, char** argv)
506                   << "pT: " << l.vec.Pt()  << "\t"
507                   << "eta: " << l.vec.Eta()
508                   << endl;
509 <            //    }
509 >          }
510          }
511          foTree.Fill();
512 +      } else {
513 +        cout << "failing with some other code : " << ret4l.status.selectionBits << endl;
514        }
515      }  
516  
# Line 581 | Line 541 | int main(int argc, char** argv)
541  
542   }
543  
584 //----------------------------------------------------------------------------
585 //
586 // Get primary vertices
587 // Assumes primary vertices are ordered by sum-pT^2 (as should be in CMSSW)
588 // NOTE: if no PV is found from fitting tracks, the beamspot is used
589 //
590 //----------------------------------------------------------------------------
591 bool setPV(ControlFlags ctrl,
592           const mithep::Array<mithep::Vertex> * vtxArr,
593           mithep::Vertex & vtx)
594 //----------------------------------------------------------------------------
595 {
596  const Vertex *bestPV = 0;
597  float best_sumpt=-1;
598
599  // good PV requirements
600  const UInt_t   fMinNTracksFit = 0;
601  const Double_t fMinNdof       = 4;
602  const Double_t fMaxAbsZ       = 24;
603  const Double_t fMaxRho        = 2;
604  
605  for(int i=0; i<vtxArr->GetEntries(); ++i) {
606    const Vertex *pv = (Vertex*)(vtxArr->At(i));
607    if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
608    
609    // Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
610    if(!pv->IsValid())                                continue;
611    if(pv->NTracksFit()       < fMinNTracksFit)       continue;
612    if(pv->Ndof()                 < fMinNdof)         continue;
613    if(fabs(pv->Z()) > fMaxAbsZ)                      continue;
614    if(pv->Position().Rho()   > fMaxRho)              continue;    
615    
616    // take the first ...
617    bestPV = pv;
618    break;
619
620    // this never reached ...    
621    float tmp_sumpt=0;
622    for( int t=0; t<pv->NTracks(); t++ )
623      tmp_sumpt += pv->Trk(t)->Pt();
624    
625    if( tmp_sumpt > best_sumpt  ) {
626      bestPV = pv;
627      best_sumpt = tmp_sumpt;
628      if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
629    }
630  }
631  // sync
632  if(!bestPV) bestPV = vtxArr->At(0);
633  //if(!bestPV) return false;
634  vtx.SetPosition(bestPV->X(),bestPV->Y(),bestPV->Z());
635  vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr());
636  return true;
637 };
638
544  
545  
546   //----------------------------------------------------------------------------
# Line 700 | Line 605 | void initRunLumiRangeMap()
605    // r11b
606    rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));  
607    */
608 +
609 +  // 2012 only for now ...
610 +  rlrm.AddJSONFile(std::string("./data/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt"));  
611 +
612   };
613  
614  
# Line 728 | Line 637 | void initEvtRhoMap( map<unsigned, float>
637    cout << "initialzing EvtRhoMap ";
638    //FILE * f = fopen("evtrho.dat", "r");
639    //  FILE * f = fopen("allrho.cali.uniq.dat", "r");
640 <  FILE * f = fopen("allrhoAA.cali.uniq.dat", "r");
640 >  FILE * f = fopen("/data/blue/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
641    int lcount=0;
642    while( fscanf( f, "%u:%f", &evt, &rho ) ) {
643      if( feof(f) ) break;
# Line 741 | Line 650 | void initEvtRhoMap( map<unsigned, float>
650  
651   //----------------------------------------------------------------------------
652   unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
653 <                     vector<bool> &pfNoPileUpflag,
654 <                     mithep::Array<Vertex>      * vtxArr )
653 >                         vector<bool> &pfNoPileUpflag,
654 >                         const mithep::Array<mithep::Vertex>      * vtxArr )
655   //----------------------------------------------------------------------------
656   {
657    for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
# Line 810 | Line 719 | unsigned makePFnoPUArray(mithep::Array<P
719    return pfNoPileUpflag.size();
720   }
721  
722 + void setEra(string fname, ControlFlags &ctrl)
723 + {
724 +  if( ctrl.mc && (fname.find("f11-") != string::npos)  ) ctrl.era=2011 ;
725 +  if( ctrl.mc && (fname.find("s12-") != string::npos)  ) ctrl.era=2012 ;
726 +  if( !ctrl.mc && (fname.find("r11-") != string::npos)  ) ctrl.era=2011 ;
727 +  if( !ctrl.mc && (fname.find("r12-") != string::npos)  ) ctrl.era=2012 ;
728 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines