Selection/src/applySelection.cc

//
// System headers
//
#include <vector>                   // STL vector class
#include <iostream>                 // standard I/O
#include <iomanip>                  // functions to format standard I/O
#include <bitset>  
#include <map>  
using namespace std;

//
// root headers
//
#include <TFile.h>
#include <TRegexp.h>
#include <TTree.h>                  
#include <TChain.h>                 
#include <TBranch.h>                
#include <TClonesArray.h>           


//
// TMVA
//
#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/Config.h"
#include "TMVA/MethodBDT.h"

#include "FactorizedJetCorrector.h"
#include "JetCorrectorParameters.h"

//
// ntuple format headers
//
#include "EventHeader.h"
#include "Electron.h"
#include "Muon.h"
#include "Vertex.h"
#include "PFMet.h"
#include "PFJetCol.h"
#include "PFCandidate.h"
#include "PFCandidateCol.h"
#include "PileupInfoCol.h"
#include "PileupEnergyDensity.h"
#include "MCParticle.h"
#include "TriggerMask.h"
#include "TriggerTable.h"
#include "TriggerObject.h"
#include "TriggerObjectRel.h"
#include "Names.h"
#include "BaseMod.h"
#include "TriggerObjectsTable.h"
#include "JetIDMVA.h"

//
// our headers
//
#include "ParseArgs.h"
#include "MuonSelection.h"
#include "ElectronSelection.h"
#include "ElectronMomentumCorrection.h"
#include "JetSelection.h"
#include "IsolationSelection.h"
#include "ReferenceSelection.h"
#include "MuonMomentumCorrection.h"

#include "TriggerUtils.h"
#include "PassHLT.h"
#include "Angles.h"
#include "KinematicsStruct.h"
#include "InfoStruct.h"
#include "GenInfoStruct.h"
#include "WeightStruct.h"
#include "AngleTuple.h"
#include "FOTuple.h"

#include "RunLumiRangeMap.h"
#include "SampleWeight.h"
#include "EfficiencyWeightsInterface.h"
#include "SelectionFuncs.h"

#include "SimpleLepton.h"

#ifndef CMSSW_BASE
#define CMSSW_BASE "../"
#endif

using namespace mithep;

//
// globals
//----------------------------------------------------------------------------
TH1D *hpu_2011, *hpu_2012;
RunLumiRangeMap rlrm;
vector<SimpleLepton> failingLeptons ; // for fake estimation
vector<SimpleLepton> passingLeptons;      // for fake estimation
vector<unsigned> cutvec;
vector<vector<unsigned> > zcutvec;
vector<vector<unsigned> > zzcutvec;
map<unsigned,float>       evtrhoMap;
vector<string> cutstrs;
bool passes_HLT;
vector<bool>   PFnoPUflag;;
map<TString, map<TString,int>* > counts;
ElectronMomentumCorrection electron_momentum_correction;
MuCorr *muCorr;
//----------------------------------------------------------------------------
void initRunLumiRangeMap(ControlFlags &ctrl) 
//----------------------------------------------------------------------------
{
  if(ctrl.jsonFile==""){
    cout << "WARNING: running with empty jsonFile, so setting noJSON to true" << endl;
    ctrl.noJSON = true;
  }
  rlrm.AddJSONFile(string(ctrl.jsonFile));  
}

//
// prototypes
//----------------------------------------------------------------------------

//
// MAIN
//----------------------------------------------------------------------------
int main(int argc, char** argv) 
//----------------------------------------------------------------------------
{
  vector<TString> cuts;
  cuts.push_back("start"); // NOTE: is incremented in ReferenceSelection, so if you turn on JSONS be careful
  cuts.push_back("trigger");
  cuts.push_back("sfOsHiPt");
  cuts.push_back("4<mZ2<120");
  cuts.push_back("pt>20,10");
  cuts.push_back("mll>4");
  cuts.push_back("m4l>70");
  cuts.push_back("mZ2>12");
  cuts.push_back("m4l>100");

  for(unsigned icut=0; icut<cuts.size(); icut++)
    counts[cuts[icut]] = new map<TString,int>;

  map<TString,map<TString,int>* >::iterator cit;
  for(cit=counts.begin(); cit!=counts.end(); cit++) {
    map<TString,int> *cts = (*cit).second;
    (*cts)["all"]  = int(0);
    (*cts)["4e"]   = int(0);
    (*cts)["4m"]   = int(0);
    (*cts)["2e2m"] = int(0);
  }
  
  string cmsswpath(CMSSW_BASE);
  cmsswpath.append("/src");
  bitset<1024> triggerBits;
  vector<vector<string> > inputFiles;
  inputFiles.push_back(vector<string>());
  map<string,unsigned> entrymap; // number of unskimmed entries in each file

  //
  // args
  //--------------------------------------------------------------------------------------------------------------
  ControlFlags ctrl;
  parse_args( argc, argv, ctrl );
  if( ctrl.inputfiles.empty() &&ctrl.inputfile.empty() ) 
    {
      cerr << "usage: applySelection.exe <flags> " << endl;
      cerr << "\tmandoatory flags : " << endl;
      cerr << "\t--inputfiles |  file containing a list of ntuples to run over" << endl;
      cerr << "\t--inputfile |  a file to run over" << endl;
      cerr << "\t--outputfile | file to store  selected evet" << endl;
      return 1;
    }
  ctrl.dump();
  cout << "\n\nNOTE: now setting era with arguments... double check this doesnt mess something else up\n\n" << endl;
  bool doJets = false;//TString(ctrl.mcfmfname).Contains("doJets");
  assert(ctrl.doFSR);

  //
  // File I/O
  //--------------------------------------------------------------------------------------------------------------
  TChain * chain = new TChain("Events");
  TChain * hltchain = new TChain("HLT");
  TChain * runchain = new TChain(Names::gkRunTreeName);

  string fname;
  unsigned total_unskimmed=0;
  if( !ctrl.inputfiles.empty() ) { 
    ifstream f(ctrl.inputfiles.c_str());
    while (f >> fname) {
      if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
      if(fname.find("/store/") != string::npos) fname = "root://eoscms//"+fname;
      cout << "adding inputfile : " << fname.c_str() << endl;
      entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
      cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
      total_unskimmed += entrymap[string(fname.c_str())];
      chain->AddFile(fname.c_str());
      hltchain->AddFile(fname.c_str());
      runchain->AddFile(fname.c_str());
    }
  } else { 
    if(ctrl.inputfile.find("/store/") != string::npos) ctrl.inputfile = "root://eoscms//"+ctrl.inputfile;
    cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
    unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
    entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
    cout << "unskimmed entries: " << unsk_ents << endl;
    total_unskimmed += unsk_ents;
    chain->AddFile(ctrl.inputfile.c_str());
    hltchain->AddFile(ctrl.inputfile.c_str());
    runchain->AddFile(ctrl.inputfile.c_str());
  }
  // write the total number of unskimmed events that went into making this output file to a text file
  writeEntries(ctrl,total_unskimmed);
  ctrl.totalMC = total_unskimmed;

  const char * ofname;
  if( strcmp( ctrl.outputfile.c_str(), "") ) {
    ofname = ctrl.outputfile.c_str();
  } else {
    ofname = "tmp.root";
  }
  // table of cross section values   
  string xspath = (cmsswpath+string("/MitPhysics/data/xs.dat"));
  cout << "xspath: " << xspath.c_str() << endl;
  SimpleTable xstab(xspath.c_str());

  //
  // Setup
  //--------------------------------------------------------------------------------------------------------------
  Angles angles;
  KinematicsStruct kine;
  InfoStruct evtinfo;
  WeightStruct weights;
  GenInfoStruct geninfo;
  JetInfoStruct ji;
  vector<SimpleLepton> simpleLeptonJets;

  AngleTuple nt( (const char*)ofname, (const char*)"zznt");
  nt.makeAngleBranch(angles);
  nt.makeKinematicsBranch(kine);
  nt.makeInfoBranch(evtinfo);
  if(doJets)
    nt.makeJetInfoBranch(ji);

  int genSampleType(0);
  if(ctrl.mc) { 
    nt.makeWeightBranch(weights);
    if(ctrl.fillGen) {
      genSampleType = getGenSampleType(ofname);
      if(genSampleType == -1) {
        cout << "\ngenSampleType returned as -1, setting fillGen to false" << endl;
        ctrl.fillGen = false;
      } else
        nt.makeGenInfoBranch(geninfo);
    }
    initEfficiencyWeights();
    // initPUWeights();
  } else { 
    if(!(ctrl.noJSON) )
      initRunLumiRangeMap(ctrl);
  }

  initElectronIDMVAV1();
  initTrigger();
  TrigInfo ti;
  initAnalysisTriggers(ti);
  if(ctrl.do_escale) {
    electron_momentum_correction.initialize_regression("../EGamma/EGammaAnalysisTools/data/eleEnergyRegWeights_V1.root");
    if(ctrl.mc) {
      if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Fall11";
      else                 ctrl.smearing_scaling_dataset = "Summer12_DR53X_HCP2012";
    } else {
      if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Jan16ReReco";
      else                 ctrl.smearing_scaling_dataset = "2012Jul13ReReco";
    }
  }
  if(ctrl.correct_muon_momentum) {
    muCorr = new MuCorr;
    muCorr->corr2011 = new rochcor;
    muCorr->corr2012 = new rochcor2012;
  }

  vector<TString> fJetCorrParsv;
  vector<JetCorrectorParameters> correctionParameters;
  FactorizedJetCorrector *fJetCorrector;
  JetIDMVA *fJetIDMVA;
  if(doJets) {
    initJets(ctrl, cmsswpath, fJetCorrParsv, correctionParameters, fJetCorrector, fJetIDMVA);
  }

  //
  // Setup tree I/O
  //--------------------------------------------------------------------------------------------------------------
  TFile *inputFile=0;
  TTree *eventTree=0;  
  string currentFile("");

  UInt_t fNMaxTriggers = TRIGGER_BIG_NUMBER;
  EventHeader                   *info                    = new EventHeader();
  Array<PFMet>                  *metArr                  = new Array<PFMet>();
  Array<Electron>               *electronArr             = new Array<Electron>();
  Array<Muon>                   *muonArr                 = new Array<Muon>();
  Array<Vertex>                 *vtxArr                  = new Array<Vertex>();
  Array<PFCandidate>            *pfArr                   = new Array<PFCandidate>();
  Array<PFJet>                  *jetArr                  = new Array<PFJet>();
  Array<PileupInfo>             *puArr                   = new Array<PileupInfo>();
  Array<PileupEnergyDensity>    *puDArr                  = new Array<PileupEnergyDensity>();
  Array<Track>                  *trkArr                  = new Array<Track>();
  Array<MCParticle>             *mcArr                   = new Array<MCParticle>();
  MCEventInfo                   *mcEvtInfo               = new MCEventInfo();
  TriggerMask                   *trigMask                = new TriggerMask();
  TriggerTable                  *hltTable                = new TriggerTable(fNMaxTriggers);
  vector<string>                *hltTableStrings         = new vector<string>();
  vector<string>                *hltLabelStrings         = new vector<string>();
  Array<TriggerObject>          *hltObjArr               = new Array<TriggerObject>();
  Array<TriggerObjectRel>       *hltRelsArr              = new Array<TriggerObjectRel>();
  TriggerObjectsTable           *fTrigObjs               = new TriggerObjectsTable(hltTable,fNMaxTriggers);
  RunInfo                       *runInfo                 = new RunInfo();
  
  chain->SetBranchAddress(Names::gkEvtHeaderBrn,                &info);
  chain->SetBranchAddress("PFMet",                              &metArr);
  chain->SetBranchAddress(Names::gkElectronBrn,                 &electronArr);
  chain->SetBranchAddress(Names::gkMuonBrn,                     &muonArr);
  chain->SetBranchAddress(Names::gkPVBrn,                       &vtxArr);
  chain->SetBranchAddress(Names::gkPFCandidatesBrn,             &pfArr);
  chain->SetBranchAddress(Names::gkPFJetBrn,                    &jetArr);
  if( ctrl.mc ) { 
    chain->SetBranchAddress(Names::gkPileupInfoBrn,             &puArr);
    chain->SetBranchAddress(Names::gkMCPartBrn,                 &mcArr);
    chain->SetBranchAddress(Names::gkMCEvtInfoBrn,              &mcEvtInfo);
  }
  chain->SetBranchAddress(Names::gkPileupEnergyDensityBrn,      &puDArr);
  chain->SetBranchAddress(Names::gkTrackBrn,                    &trkArr);
  chain->SetBranchAddress(Names::gkHltBitBrn,                   &trigMask);
  chain->SetBranchAddress(Names::gkHltObjBrn,                   &hltObjArr);
  chain->SetBranchAddress("HLTObjectsRelation",                 &hltRelsArr);
  hltchain->SetBranchAddress(Names::gkHltTableBrn,              &hltTableStrings);
  hltchain->SetBranchAddress(Names::gkHltLabelBrn,              &hltLabelStrings);

  runchain->SetBranchAddress(Names::gkRunInfoBrn, &runInfo);

  Vertex              vtx;          // best primary vertex in the event

  hltchain->GetEntry(0);

  int imax = (ctrl.n_events_to_process < 0) ? chain->GetEntries() : ctrl.n_events_to_process;
  cerr << "processing: " << imax << " ( / " << chain->GetEntries() << " )" << endl;

  if(ctrl.mc && ctrl.debugOffMc) // necessary for batch mode -- we need debug output for, data but size is too large to store MC output
    ctrl.debug = false;

  int lastHltEntry=-2;
  for(UInt_t ientry=0; ientry<imax; ientry++) 
    {
      chain->GetEntry(ientry);
      if(!(ientry%1000)) cerr << "entry: " << ientry << " ( / " << imax << " )" << endl;

      if( ctrl.debug ) {
        cout << "-----------------------------------------------------------------" << endl;
        cout << "-----------------------------------------------------------------" << endl;
        cout << "Run: " << info->RunNum() 
             << "\tEvt: " << info->EvtNum()
             << "\tLumi: " << info->LumiSec()
             << endl;
        cerr << "Run: " << info->RunNum() 
             << "\tEvt: " << info->EvtNum()
             << "\tLumi: " << info->LumiSec()
             << endl;
        cout << "-----------------------------------------------------------------" << endl;
      }

      runchain->GetEvent(info->RunEntry());
      hltchain->GetEntry(runInfo->HltEntry());
      hltTable->Clear();
      fillTriggerNames(hltTable, hltTableStrings );
      if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print();
      lastHltEntry = runInfo->HltEntry();
      fillTriggerBits( hltTable, trigMask, triggerBits );
      setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
      // printTriggerObjs( hltTable, fTrigObjs);
      
      string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
      if(ctrl.debug) cout << "era is " << ctrl.era << endl;
      if( ctrl.era == 0  ) {
        setEra( fname, ctrl );
        if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
      }

      // pfNoPU 
      PFnoPUflag.clear();
      int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
      assert(pfnopu_size == pfArr->GetEntries());

      // trigger
      if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
        currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
      }
      if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
      passes_HLT = passHLT(ctrl, triggerBits, ti, info->RunNum(), hltTable, fTrigObjs); // NOTE: run value is not used if it's MC

      
      // data/MC
      if(ctrl.mc) {
        if(ctrl.fillGen) fillGenInfo( mcArr, mcEvtInfo, geninfo, genSampleType, ctrl);
        // NOTE: xs weight and npuw are reset in merge.cc
        weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
        if(ctrl.debug) cout ;
        weights.npu = getNPU(puArr, 0);
        weights.npuw = getPUWeight(ctrl.era, fname, weights.npu);
      } else { 
        // JSON
        if(!(ctrl.noJSON) ) {
          RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
          if( !(rlrm.HasRunLumi(rl)) )  {
            if( ctrl.debug ) cout << "fails JSON" << endl;
            continue;
          }
        }
      }

      //
      // lepton/kinematic selection ...
      //
      failingLeptons.clear();
      passingLeptons.clear();
      EventData ret4l = apply_HZZ4L_reference_selection(ctrl, info, 
                                                        vtxArr,
                                                        pfArr,
                                                        puDArr,
                                                        electronArr, 
                                                        &electronReferencePreSelection,
                                                        &electronReferenceIDMVASelectionV1, // not used anymore! (Blame the damn regression implementation, not me!)
                                                        &electronReferenceIsoSelection,
                                                        muonArr, 
                                                        &muonReferencePreSelection,
                                                        &muonIDPFSelection,
                                                        &muonReferenceIsoSelection);
      
      if( ctrl.debug ) cout << endl;

      if( ret4l.status.pass() ) { 
        fillAngles(ret4l,angles);

        fillKinematics(ret4l,kine);
        // fillMassErrors(ret4l,muonArr,electronArr,kine);

        if(doJets) {
          fillJets(ctrl, jetArr, simpleLeptonJets, fJetCorrector, puDArr, fJetIDMVA, vtxArr, ret4l);
          fillJetInfo(simpleLeptonJets,ji,ctrl);
        }

        TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
        fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries());
        if( ctrl.mc) { 
          setEfficiencyWeights(ctrl, ctrl.era, ret4l, triggerBits, hltTable, hltObjArr, fTrigObjs, weights); 
          if(ctrl.debug) 
            cout << "w: " << weights.w << "\t" << "npuw: " << weights.npuw << "\t" << "won: " << weights.won << "\t" << "woff: " << weights.woff << endl;
        }
        
        nt.Fill();
        
        cerr  << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
              << kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;
        cout  << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
              << kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;

      } else { 
        if(ctrl.debug)
          cout << "failing with code : " << ret4l.status.selectionBits << endl;
      }
    }

  nt.WriteClose();

  for(unsigned icut=0; icut<cuts.size(); icut++) {
    map<TString,int> thisCount(*(counts[cuts[icut]]));
    cout << setw(20) << cuts[icut] << "\t" << setw(9) << thisCount["all"] << "  ("
         << setw(9) << thisCount["4e"] << setw(9) << thisCount["4m"] << setw(9) << thisCount["2e2m"] << ")" << endl;
  }
} 
Revision:	1.61
Committed:	Wed Oct 31 14:28:14 2012 UTC (12 years, 6 months ago) by anlevin
Content type:	text/plain
Branch:	MAIN
Changes since 1.60:	+6 -5 lines
Log Message:	added trigger object matching for 2011
#	User	Rev	Content
1	khahn	1.1	//
2			// System headers
3			//
4			#include <vector> // STL vector class
5			#include <iostream> // standard I/O
6			#include <iomanip> // functions to format standard I/O
7	khahn	1.13	#include <bitset>
8	khahn	1.19	#include <map>
9	khahn	1.1	using namespace std;
10
11			//
12	khahn	1.13	// root headers
13	khahn	1.1	//
14	dkralph	1.60	#include <TFile.h>
15			#include <TRegexp.h>
16	khahn	1.13	#include <TTree.h>
17			#include <TChain.h>
18			#include <TBranch.h>
19			#include <TClonesArray.h>
20	khahn	1.1
21	khahn	1.11
22	khahn	1.1	//
23	khahn	1.13	// TMVA
24			//
25			#include "TMVA/Reader.h"
26			#include "TMVA/Tools.h"
27			#include "TMVA/Config.h"
28			#include "TMVA/MethodBDT.h"
29
30	dkralph	1.60	#include "FactorizedJetCorrector.h"
31			#include "JetCorrectorParameters.h"
32
33	khahn	1.13	//
34	khahn	1.1	// ntuple format headers
35			//
36	khahn	1.13	#include "EventHeader.h"
37			#include "Electron.h"
38			#include "Muon.h"
39			#include "Vertex.h"
40	dkralph	1.48	#include "PFMet.h"
41	dkralph	1.60	#include "PFJetCol.h"
42	khahn	1.13	#include "PFCandidate.h"
43			#include "PFCandidateCol.h"
44	khahn	1.16	#include "PileupInfoCol.h"
45			#include "PileupEnergyDensity.h"
46	khahn	1.23	#include "MCParticle.h"
47	khahn	1.13	#include "TriggerMask.h"
48			#include "TriggerTable.h"
49	khahn	1.41	#include "TriggerObject.h"
50			#include "TriggerObjectRel.h"
51	khahn	1.13	#include "Names.h"
52			#include "BaseMod.h"
53	dkralph	1.44	#include "TriggerObjectsTable.h"
54	dkralph	1.60	#include "JetIDMVA.h"
55	khahn	1.43
56	khahn	1.1	//
57	khahn	1.13	// our headers
58	khahn	1.1	//
59			#include "ParseArgs.h"
60	khahn	1.13	#include "MuonSelection.h"
61			#include "ElectronSelection.h"
62	anlevin	1.56	#include "ElectronMomentumCorrection.h"
63	dkralph	1.60	#include "JetSelection.h"
64	khahn	1.13	#include "IsolationSelection.h"
65	khahn	1.29	#include "ReferenceSelection.h"
66	dkralph	1.60	#include "MuonMomentumCorrection.h"
67	khahn	1.14
68	khahn	1.13	#include "TriggerUtils.h"
69			#include "PassHLT.h"
70			#include "Angles.h"
71			#include "KinematicsStruct.h"
72			#include "InfoStruct.h"
73	khahn	1.28	#include "GenInfoStruct.h"
74	khahn	1.13	#include "WeightStruct.h"
75			#include "AngleTuple.h"
76	khahn	1.15	#include "FOTuple.h"
77	khahn	1.13
78	khahn	1.33	#include "RunLumiRangeMap.h"
79	dkralph	1.2	#include "SampleWeight.h"
80	khahn	1.13	#include "EfficiencyWeightsInterface.h"
81	dkralph	1.38	#include "SelectionFuncs.h"
82	khahn	1.1
83	khahn	1.16	#include "SimpleLepton.h"
84
85	khahn	1.13	#ifndef CMSSW_BASE
86			#define CMSSW_BASE "../"
87			#endif
88	khahn	1.1
89	khahn	1.13	using namespace mithep;
90	khahn	1.11
91	khahn	1.16	//
92			// globals
93			//----------------------------------------------------------------------------
94	khahn	1.47	TH1D hpu_2011, hpu_2012;
95	khahn	1.16	RunLumiRangeMap rlrm;
96			vector<SimpleLepton> failingLeptons ; // for fake estimation
97			vector<SimpleLepton> passingLeptons; // for fake estimation
98	khahn	1.19	vector<unsigned> cutvec;
99			vector<vector<unsigned> > zcutvec;
100			vector<vector<unsigned> > zzcutvec;
101			map<unsigned,float> evtrhoMap;
102			vector<string> cutstrs;
103	dkralph	1.60	bool passes_HLT;
104	khahn	1.22	vector<bool> PFnoPUflag;;
105	anlevin	1.57	map<TString, map<TString,int>* > counts;
106	dkralph	1.60	ElectronMomentumCorrection electron_momentum_correction;
107			MuCorr *muCorr;
108			//----------------------------------------------------------------------------
109			void initRunLumiRangeMap(ControlFlags &ctrl)
110			//----------------------------------------------------------------------------
111			{
112			if(ctrl.jsonFile==""){
113			cout << "WARNING: running with empty jsonFile, so setting noJSON to true" << endl;
114			ctrl.noJSON = true;
115			}
116			rlrm.AddJSONFile(string(ctrl.jsonFile));
117			}
118	khahn	1.11
119	khahn	1.16	//
120			// prototypes
121	khahn	1.13	//----------------------------------------------------------------------------
122	khahn	1.11
123	khahn	1.16	//
124			// MAIN
125			//----------------------------------------------------------------------------
126	khahn	1.1	int main(int argc, char** argv)
127	khahn	1.16	//----------------------------------------------------------------------------
128	khahn	1.1	{
129	anlevin	1.57	vector<TString> cuts;
130	dkralph	1.60	cuts.push_back("start"); // NOTE: is incremented in ReferenceSelection, so if you turn on JSONS be careful
131	anlevin	1.57	cuts.push_back("trigger");
132			cuts.push_back("sfOsHiPt");
133			cuts.push_back("4<mZ2<120");
134			cuts.push_back("pt>20,10");
135			cuts.push_back("mll>4");
136			cuts.push_back("m4l>70");
137	dkralph	1.60	cuts.push_back("mZ2>12");
138			cuts.push_back("m4l>100");
139	anlevin	1.57
140			for(unsigned icut=0; icut<cuts.size(); icut++)
141			counts[cuts[icut]] = new map<TString,int>;
142
143			map<TString,map<TString,int>* >::iterator cit;
144			for(cit=counts.begin(); cit!=counts.end(); cit++) {
145			map<TString,int> cts = (cit).second;
146			(*cts)["all"] = int(0);
147			(*cts)["4e"] = int(0);
148			(*cts)["4m"] = int(0);
149			(*cts)["2e2m"] = int(0);
150			}
151	dkralph	1.60
152	khahn	1.13	string cmsswpath(CMSSW_BASE);
153			cmsswpath.append("/src");
154	dkralph	1.60	bitset<1024> triggerBits;
155	khahn	1.1	vector<vector<string> > inputFiles;
156			inputFiles.push_back(vector<string>());
157	khahn	1.16	map<string,unsigned> entrymap; // number of unskimmed entries in each file
158	khahn	1.1
159			//
160			// args
161			//--------------------------------------------------------------------------------------------------------------
162			ControlFlags ctrl;
163			parse_args( argc, argv, ctrl );
164	khahn	1.13	if( ctrl.inputfiles.empty() &&ctrl.inputfile.empty() )
165			{
166			cerr << "usage: applySelection.exe <flags> " << endl;
167			cerr << "\tmandoatory flags : " << endl;
168			cerr << "\t--inputfiles \| file containing a list of ntuples to run over" << endl;
169			cerr << "\t--inputfile \| a file to run over" << endl;
170			cerr << "\t--outputfile \| file to store selected evet" << endl;
171			return 1;
172			}
173	khahn	1.1	ctrl.dump();
174	dkralph	1.60	cout << "\n\nNOTE: now setting era with arguments... double check this doesnt mess something else up\n\n" << endl;
175	anlevin	1.61	bool doJets = false;//TString(ctrl.mcfmfname).Contains("doJets");
176	dkralph	1.60	assert(ctrl.doFSR);
177	dkralph	1.4
178	khahn	1.1	//
179			// File I/O
180			//--------------------------------------------------------------------------------------------------------------
181			TChain * chain = new TChain("Events");
182	khahn	1.13	TChain * hltchain = new TChain("HLT");
183	dkralph	1.53	TChain * runchain = new TChain(Names::gkRunTreeName);
184	khahn	1.13
185	khahn	1.1	string fname;
186	dkralph	1.44	unsigned total_unskimmed=0;
187	khahn	1.13	if( !ctrl.inputfiles.empty() ) {
188			ifstream f(ctrl.inputfiles.c_str());
189	dkralph	1.60	while (f >> fname) {
190	khahn	1.13	if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
191	dkralph	1.60	if(fname.find("/store/") != string::npos) fname = "root://eoscms//"+fname;
192	khahn	1.13	cout << "adding inputfile : " << fname.c_str() << endl;
193			entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
194	khahn	1.16	cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
195	dkralph	1.44	total_unskimmed += entrymap[string(fname.c_str())];
196	khahn	1.13	chain->AddFile(fname.c_str());
197			hltchain->AddFile(fname.c_str());
198	dkralph	1.53	runchain->AddFile(fname.c_str());
199	khahn	1.13	}
200			} else {
201	dkralph	1.60	if(ctrl.inputfile.find("/store/") != string::npos) ctrl.inputfile = "root://eoscms//"+ctrl.inputfile;
202	khahn	1.13	cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
203	dkralph	1.44	unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
204			entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
205			cout << "unskimmed entries: " << unsk_ents << endl;
206			total_unskimmed += unsk_ents;
207	khahn	1.13	chain->AddFile(ctrl.inputfile.c_str());
208			hltchain->AddFile(ctrl.inputfile.c_str());
209	dkralph	1.53	runchain->AddFile(ctrl.inputfile.c_str());
210	khahn	1.1	}
211	dkralph	1.60	// write the total number of unskimmed events that went into making this output file to a text file
212	dkralph	1.53	writeEntries(ctrl,total_unskimmed);
213	dkralph	1.60	ctrl.totalMC = total_unskimmed;
214	khahn	1.1
215	dkralph	1.10	const char * ofname;
216			if( strcmp( ctrl.outputfile.c_str(), "") ) {
217			ofname = ctrl.outputfile.c_str();
218			} else {
219			ofname = "tmp.root";
220			}
221	khahn	1.13	// table of cross section values
222			string xspath = (cmsswpath+string("/MitPhysics/data/xs.dat"));
223			cout << "xspath: " << xspath.c_str() << endl;
224			SimpleTable xstab(xspath.c_str());
225	khahn	1.11
226	khahn	1.1	//
227	khahn	1.13	// Setup
228	khahn	1.1	//--------------------------------------------------------------------------------------------------------------
229	khahn	1.13	Angles angles;
230	anlevin	1.57	KinematicsStruct kine;
231	khahn	1.13	InfoStruct evtinfo;
232			WeightStruct weights;
233	khahn	1.28	GenInfoStruct geninfo;
234	dkralph	1.60	JetInfoStruct ji;
235			vector<SimpleLepton> simpleLeptonJets;
236	khahn	1.13
237			AngleTuple nt( (const char)ofname, (const char)"zznt");
238			nt.makeAngleBranch(angles);
239	anlevin	1.57	nt.makeKinematicsBranch(kine);
240	khahn	1.13	nt.makeInfoBranch(evtinfo);
241	dkralph	1.60	if(doJets)
242			nt.makeJetInfoBranch(ji);
243	khahn	1.17
244	dkralph	1.60	int genSampleType(0);
245	khahn	1.13	if(ctrl.mc) {
246			nt.makeWeightBranch(weights);
247	dkralph	1.60	if(ctrl.fillGen) {
248			genSampleType = getGenSampleType(ofname);
249			if(genSampleType == -1) {
250			cout << "\ngenSampleType returned as -1, setting fillGen to false" << endl;
251			ctrl.fillGen = false;
252			} else
253			nt.makeGenInfoBranch(geninfo);
254	khahn	1.13	}
255	dkralph	1.60	initEfficiencyWeights();
256			// initPUWeights();
257	khahn	1.13	} else {
258	dkralph	1.60	if(!(ctrl.noJSON) )
259			initRunLumiRangeMap(ctrl);
260	khahn	1.13	}
261
262	khahn	1.33	initElectronIDMVAV1();
263	khahn	1.13	initTrigger();
264	dkralph	1.53	TrigInfo ti;
265			initAnalysisTriggers(ti);
266	dkralph	1.60	if(ctrl.do_escale) {
267	anlevin	1.61	electron_momentum_correction.initialize_regression("../EGamma/EGammaAnalysisTools/data/eleEnergyRegWeights_V1.root");
268	dkralph	1.60	if(ctrl.mc) {
269	anlevin	1.61	if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Fall11";
270	dkralph	1.60	else ctrl.smearing_scaling_dataset = "Summer12_DR53X_HCP2012";
271			} else {
272	anlevin	1.61	if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Jan16ReReco";
273	dkralph	1.60	else ctrl.smearing_scaling_dataset = "2012Jul13ReReco";
274			}
275			}
276			if(ctrl.correct_muon_momentum) {
277			muCorr = new MuCorr;
278			muCorr->corr2011 = new rochcor;
279			muCorr->corr2012 = new rochcor2012;
280			}
281
282			vector<TString> fJetCorrParsv;
283			vector<JetCorrectorParameters> correctionParameters;
284			FactorizedJetCorrector *fJetCorrector;
285			JetIDMVA *fJetIDMVA;
286			if(doJets) {
287			initJets(ctrl, cmsswpath, fJetCorrParsv, correctionParameters, fJetCorrector, fJetIDMVA);
288			}
289	khahn	1.19
290	khahn	1.16	//
291	khahn	1.13	// Setup tree I/O
292	khahn	1.16	//--------------------------------------------------------------------------------------------------------------
293	khahn	1.1	TFile *inputFile=0;
294			TTree *eventTree=0;
295	khahn	1.13	string currentFile("");
296
297	dkralph	1.60	UInt_t fNMaxTriggers = TRIGGER_BIG_NUMBER;
298			EventHeader *info = new EventHeader();
299			Array<PFMet> *metArr = new Array<PFMet>();
300			Array<Electron> *electronArr = new Array<Electron>();
301			Array<Muon> *muonArr = new Array<Muon>();
302			Array<Vertex> *vtxArr = new Array<Vertex>();
303			Array<PFCandidate> *pfArr = new Array<PFCandidate>();
304			Array<PFJet> *jetArr = new Array<PFJet>();
305			Array<PileupInfo> *puArr = new Array<PileupInfo>();
306			Array<PileupEnergyDensity> *puDArr = new Array<PileupEnergyDensity>();
307			Array<Track> *trkArr = new Array<Track>();
308			Array<MCParticle> *mcArr = new Array<MCParticle>();
309			MCEventInfo *mcEvtInfo = new MCEventInfo();
310			TriggerMask *trigMask = new TriggerMask();
311			TriggerTable *hltTable = new TriggerTable(fNMaxTriggers);
312			vector<string> *hltTableStrings = new vector<string>();
313			vector<string> *hltLabelStrings = new vector<string>();
314			Array<TriggerObject> *hltObjArr = new Array<TriggerObject>();
315			Array<TriggerObjectRel> *hltRelsArr = new Array<TriggerObjectRel>();
316			TriggerObjectsTable *fTrigObjs = new TriggerObjectsTable(hltTable,fNMaxTriggers);
317			RunInfo *runInfo = new RunInfo();
318	dkralph	1.44
319	dkralph	1.60	chain->SetBranchAddress(Names::gkEvtHeaderBrn, &info);
320			chain->SetBranchAddress("PFMet", &metArr);
321			chain->SetBranchAddress(Names::gkElectronBrn, &electronArr);
322			chain->SetBranchAddress(Names::gkMuonBrn, &muonArr);
323			chain->SetBranchAddress(Names::gkPVBrn, &vtxArr);
324			chain->SetBranchAddress(Names::gkPFCandidatesBrn, &pfArr);
325			chain->SetBranchAddress(Names::gkPFJetBrn, &jetArr);
326	khahn	1.23	if( ctrl.mc ) {
327	dkralph	1.60	chain->SetBranchAddress(Names::gkPileupInfoBrn, &puArr);
328			chain->SetBranchAddress(Names::gkMCPartBrn, &mcArr);
329			chain->SetBranchAddress(Names::gkMCEvtInfoBrn, &mcEvtInfo);
330	khahn	1.23	}
331	dkralph	1.60	chain->SetBranchAddress(Names::gkPileupEnergyDensityBrn, &puDArr);
332			chain->SetBranchAddress(Names::gkTrackBrn, &trkArr);
333			chain->SetBranchAddress(Names::gkHltBitBrn, &trigMask);
334			chain->SetBranchAddress(Names::gkHltObjBrn, &hltObjArr);
335			chain->SetBranchAddress("HLTObjectsRelation", &hltRelsArr);
336			hltchain->SetBranchAddress(Names::gkHltTableBrn, &hltTableStrings);
337			hltchain->SetBranchAddress(Names::gkHltLabelBrn, &hltLabelStrings);
338	khahn	1.1
339	dkralph	1.53	runchain->SetBranchAddress(Names::gkRunInfoBrn, &runInfo);
340
341	dkralph	1.60	Vertex vtx; // best primary vertex in the event
342	khahn	1.11
343	khahn	1.13	hltchain->GetEntry(0);
344
345	dkralph	1.60	int imax = (ctrl.n_events_to_process < 0) ? chain->GetEntries() : ctrl.n_events_to_process;
346			cerr << "processing: " << imax << " ( / " << chain->GetEntries() << " )" << endl;
347	anlevin	1.54
348	dkralph	1.60	if(ctrl.mc && ctrl.debugOffMc) // necessary for batch mode -- we need debug output for, data but size is too large to store MC output
349			ctrl.debug = false;
350	khahn	1.11
351	dkralph	1.53	int lastHltEntry=-2;
352	khahn	1.13	for(UInt_t ientry=0; ientry<imax; ientry++)
353			{
354			chain->GetEntry(ientry);
355	dkralph	1.60	if(!(ientry%1000)) cerr << "entry: " << ientry << " ( / " << imax << " )" << endl;
356	dkralph	1.44
357	khahn	1.41	if( ctrl.debug ) {
358			cout << "-----------------------------------------------------------------" << endl;
359			cout << "-----------------------------------------------------------------" << endl;
360			cout << "Run: " << info->RunNum()
361			<< "\tEvt: " << info->EvtNum()
362			<< "\tLumi: " << info->LumiSec()
363			<< endl;
364	dkralph	1.53	cerr << "Run: " << info->RunNum()
365			<< "\tEvt: " << info->EvtNum()
366			<< "\tLumi: " << info->LumiSec()
367			<< endl;
368	khahn	1.41	cout << "-----------------------------------------------------------------" << endl;
369			}
370
371	dkralph	1.53	runchain->GetEvent(info->RunEntry());
372			hltchain->GetEntry(runInfo->HltEntry());
373			hltTable->Clear();
374			fillTriggerNames(hltTable, hltTableStrings );
375			if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print();
376			lastHltEntry = runInfo->HltEntry();
377			fillTriggerBits( hltTable, trigMask, triggerBits );
378			setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
379			// printTriggerObjs( hltTable, fTrigObjs);
380
381	khahn	1.28	string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
382	khahn	1.49	if(ctrl.debug) cout << "era is " << ctrl.era << endl;
383			if( ctrl.era == 0 ) {
384			setEra( fname, ctrl );
385			if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
386			}
387	khahn	1.41
388	khahn	1.22	// pfNoPU
389			PFnoPUflag.clear();
390			int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
391			assert(pfnopu_size == pfArr->GetEntries());
392	khahn	1.13
393	dkralph	1.60	// trigger
394	dkralph	1.53	if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
395			currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
396			}
397			if( ctrl.debug ) cout << "file is : " << currentFile << endl;
398	anlevin	1.61	passes_HLT = passHLT(ctrl, triggerBits, ti, info->RunNum(), hltTable, fTrigObjs); // NOTE: run value is not used if it's MC
399	anlevin	1.58
400	anlevin	1.61
401	khahn	1.16	// data/MC
402	khahn	1.13	if(ctrl.mc) {
403	dkralph	1.60	if(ctrl.fillGen) fillGenInfo( mcArr, mcEvtInfo, geninfo, genSampleType, ctrl);
404			// NOTE: xs weight and npuw are reset in merge.cc
405	khahn	1.31	weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
406	dkralph	1.60	if(ctrl.debug) cout ;
407			weights.npu = getNPU(puArr, 0);
408			weights.npuw = getPUWeight(ctrl.era, fname, weights.npu);
409	khahn	1.13	} else {
410	khahn	1.16	// JSON
411	khahn	1.40	if(!(ctrl.noJSON) ) {
412	khahn	1.36	RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
413			if( !(rlrm.HasRunLumi(rl)) ) {
414	khahn	1.41	if( ctrl.debug ) cout << "fails JSON" << endl;
415	khahn	1.36	continue;
416	khahn	1.40	}
417	khahn	1.16	}
418	khahn	1.13	}
419	khahn	1.16
420	khahn	1.13	//
421			// lepton/kinematic selection ...
422			//
423	khahn	1.16	failingLeptons.clear();
424			passingLeptons.clear();
425	dkralph	1.60	EventData ret4l = apply_HZZ4L_reference_selection(ctrl, info,
426			vtxArr,
427			pfArr,
428			puDArr,
429			electronArr,
430			&electronReferencePreSelection,
431			&electronReferenceIDMVASelectionV1, // not used anymore! (Blame the damn regression implementation, not me!)
432			&electronReferenceIsoSelection,
433			muonArr,
434			&muonReferencePreSelection,
435			&muonIDPFSelection,
436			&muonReferenceIsoSelection);
437	khahn	1.33
438	khahn	1.13	if( ctrl.debug ) cout << endl;
439	khahn	1.14
440	khahn	1.13	if( ret4l.status.pass() ) {
441			fillAngles(ret4l,angles);
442	dkralph	1.60
443	anlevin	1.57	fillKinematics(ret4l,kine);
444	dkralph	1.60	// fillMassErrors(ret4l,muonArr,electronArr,kine);
445
446			if(doJets) {
447			fillJets(ctrl, jetArr, simpleLeptonJets, fJetCorrector, puDArr, fJetIDMVA, vtxArr, ret4l);
448			fillJetInfo(simpleLeptonJets,ji,ctrl);
449			}
450
451	dkralph	1.48	TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
452	dkralph	1.60	fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries());
453	khahn	1.16	if( ctrl.mc) {
454	dkralph	1.60	setEfficiencyWeights(ctrl, ctrl.era, ret4l, triggerBits, hltTable, hltObjArr, fTrigObjs, weights);
455	dkralph	1.44	if(ctrl.debug)
456	dkralph	1.60	cout << "w: " << weights.w << "\t" << "npuw: " << weights.npuw << "\t" << "won: " << weights.won << "\t" << "woff: " << weights.woff << endl;
457	khahn	1.16	}
458	khahn	1.13
459			nt.Fill();
460
461	dkralph	1.60	cerr << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
462			<< kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;
463			cout << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
464			<< kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;
465
466	khahn	1.36	} else {
467	dkralph	1.48	if(ctrl.debug)
468	dkralph	1.60	cout << "failing with code : " << ret4l.status.selectionBits << endl;
469	khahn	1.15	}
470	dkralph	1.53	}
471	khahn	1.1
472	khahn	1.13	nt.WriteClose();
473	khahn	1.19
474	anlevin	1.57	for(unsigned icut=0; icut<cuts.size(); icut++) {
475	dkralph	1.60	map<TString,int> thisCount(*(counts[cuts[icut]]));
476			cout << setw(20) << cuts[icut] << "\t" << setw(9) << thisCount["all"] << " ("
477			<< setw(9) << thisCount["4e"] << setw(9) << thisCount["4m"] << setw(9) << thisCount["2e2m"] << ")" << endl;
478	anlevin	1.57	}
479	khahn	1.13	}