NonMCBackground/src/applyEMU.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 <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"

//
// ntuple format headers
//
#include "EventHeader.h"
#include "PFMet.h"
#include "Electron.h"
#include "Muon.h"
#include "Vertex.h"
#include "PFJet.h"
#include "PFCandidate.h"
#include "PFCandidateCol.h"
#include "PileupInfoCol.h"
#include "PileupEnergyDensity.h"
#include "MCParticle.h"
#include "TriggerMask.h"
#include "TriggerTable.h"
#include "Names.h"
#include "BaseMod.h"

//
// our headers
//
#include "ParseArgs.h"
#include "MuonSelection.h"
#include "ElectronSelection.h"
#include "IsolationSelection.h"
#include "SelectionEMU.h"

#include "TriggerUtils.h"
#include "PassHLT.h"
#include "Angles.h"
#include "KinematicsStruct.h"
#include "InfoStruct.h"
#include "GenInfoStruct.h"
#include "WeightStruct.h"
//#include "GlobalDataAndFuncs.h"
#include "RunLumiRangeMap.h"
#include "SelectionFuncs.h"

#include "AngleTuple.h"
#include "FOTuple.h"

#include "SampleWeight.h"
#include "EfficiencyWeightsInterface.h"

#include "SimpleLepton.h"

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

using namespace mithep;

//
// globals
//----------------------------------------------------------------------------
TH1D *hpu_2011, *hpu_2012; 
mithep::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_MC;
vector<bool>   PFnoPUflag;;
vector<int> muTrigObjs,eleTrigObjs,muTriggers,eleTriggers;
std::bitset<TRIGGER_BIG_NUMBER> triggerBits;

//
// prototypes
//----------------------------------------------------------------------------
void initRunLumiRangeMap(ControlFlags &ctrl); 

//
// MAIN
//----------------------------------------------------------------------------
int main(int argc, char** argv) 
//----------------------------------------------------------------------------
{

  string cmsswpath(CMSSW_BASE);
  cmsswpath.append("/src");
  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();


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

  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
      cout << "adding inputfile : " << fname.c_str() << endl;
      entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
      total_unskimmed += entrymap[string(fname.c_str())];
      chain->AddFile(fname.c_str());
      hltchain->AddFile(fname.c_str());
    }
  } else { 
    cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
    unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
    entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
    total_unskimmed += unsk_ents;
    chain->AddFile(ctrl.inputfile.c_str());
    hltchain->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);

  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 kinematics;
  InfoStruct evtinfo;
  WeightStruct weights;
  GenInfoStruct geninfo;

  AngleTuple nt( (const char*)ofname, (const char*)"zznt");
  nt.makeAngleBranch(angles);
  nt.makeKinematicsBranch(kinematics);
  nt.makeInfoBranch(evtinfo);


  FOTuple foTree( nt.getFile(), (const char*)"FOtree");
  foTree.makeFailedLeptonBranch(failingLeptons);
  foTree.makePassedLeptonBranch(passingLeptons);

  TH1F * h_wf11_hpt;
  if(ctrl.mc) { 
    nt.makeWeightBranch(weights);
    if(ctrl.fillGen ) nt.makeGenInfoBranch(geninfo);
    // initEfficiencyWeights();
    cout << "\n\nWARNING: not initializing eff weights\n\n" << endl;
    initPUWeights();

    // Higgs only, pt reweighting
    if( ctrl.mH <= 140 ) { 
      char wptname[256];
      sprintf( wptname, "/scratch/dkralph/pharris/Flat/ntupler/root/weight_ptH_%i.root", ctrl.mH );
      TFile * f_hpt = new TFile(wptname);
      f_hpt->Print();
      sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
      h_wf11_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
    }

  } else { 
    initRunLumiRangeMap(ctrl);
  }

  //  initMuonIDMVA();
  initElectronIDMVAV1();
  initMuonIsoMVA();
  initElectronIDMVA();
  initElectronIsoMVA();
  initTrigger();

  muTriggers.push_back(kHLT_IsoMu24_eta2p1);    muTrigObjs.push_back(kHLT_IsoMu24_eta2p1_MuObj);
  eleTriggers.push_back(kHLT_Ele27_WP80);       eleTrigObjs.push_back(kHLT_Ele27_WP80_EleObj);

  // hack 
  initEvtRhoMap(evtrhoMap);


  //
  // Setup tree I/O
  //--------------------------------------------------------------------------------------------------------------
  string currentFile("");

  UInt_t fNMaxTriggers = TRIGGER_BIG_NUMBER;
  mithep::EventHeader *info    = new mithep::EventHeader();
  mithep::Array<mithep::PFMet>                *metArr        = new mithep::Array<mithep::PFMet>();
  mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
  mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
  mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
  mithep::Array<mithep::PFCandidate>          *pfArr         = new mithep::Array<mithep::PFCandidate>();
  mithep::Array<mithep::PFJet>                *jetArr        = new mithep::Array<mithep::PFJet>();
  mithep::Array<mithep::PileupInfo>           *puArr         = new mithep::Array<mithep::PileupInfo>();
  mithep::Array<mithep::PileupEnergyDensity>  *puDArr        = new mithep::Array<mithep::PileupEnergyDensity>();
  mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
  mithep::Array<mithep::MCParticle>           *mcArr         = new mithep::Array<mithep::MCParticle>();
  mithep::MCEventInfo                         *mcEvtInfo     = new mithep::MCEventInfo();
  mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
  mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
  vector<string>                              *hltTableStrings  = new vector<string>();
  vector<string>                              *hltLabelStrings  = new vector<string>();
  mithep::Array<mithep::TriggerObject>        *hltObjArr     = new mithep::Array<mithep::TriggerObject>();
  mithep::Array<mithep::TriggerObjectRel>     *hltRelsArr    = new mithep::Array<mithep::TriggerObjectRel>();
  std::vector<std::string>                    *hltTab        = new vector<string>();
  std::vector<std::string>                    *hltLab        = new vector<string>();
  mithep::TriggerObjectsTable                 *fTrigObjs     = new TriggerObjectsTable(hltTable,fNMaxTriggers);
  mithep::Array<mithep::DecayParticle>        *fConversions  = new mithep::Array<mithep::DecayParticle>();

  TString fTriggerTableName(Names::gkHltTableBrn);
  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);
  chain->SetBranchAddress(Names::gkMvfConversionBrn,            &fConversions);
  hltchain->SetBranchAddress(fTriggerTableName,                 &hltTableStrings);
  hltchain->SetBranchAddress(Names::gkHltLabelBrn,              &hltLabelStrings);

  mithep::Vertex              vtx;          // best primary vertex in the event

  //  ginfo = NULL;
  // if( ctrl.mc ) { chain->SetBranchAddress("Gen",        &ginfo);}

  // only 1 HLT table / file ???
  hltchain->GetEntry(0);

  int imax = chain->GetEntries();
  cout << "nEntries: " << imax << endl;


  //
  // Loop !!!!!!!!! 
  //--------------------------------------------------------------------------------------------------------------
  for(UInt_t ientry=0; ientry<imax; ientry++) {
    chain->GetEntry(ientry);
    if(!(ientry%1000)) cerr << "entry: " << ientry << "\t" << float(ientry)/imax << endl;

    string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
    if( ctrl.era == 0  )
      setEra( fname, ctrl );

    // 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());
      hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
      hltchain->GetEntry(0);
      hltTable->Clear();
      fillTriggerNames(hltTable, hltTableStrings );
    }
    fillTriggerBits( hltTable, trigMask, triggerBits );
    setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );

    //
    // data/MC
    //
    if(ctrl.mc) {
      if( ctrl.fillGen )
        fillGenInfo( mcArr, mcEvtInfo, geninfo, ESampleType::kHZZ, ctrl);
    } else { 
      if(!(ctrl.noJSON) ) {
        RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
        if( !(rlrm.HasRunLumi(rl)) ) continue;
      }
      // if( !passHLTEMU(ctrl,triggerBits) ) {
      //   continue;
      // }
    }

    // lepton/kinematic selection ...
    failingLeptons.clear();
    passingLeptons.clear();
    EventData ret4l;
        
    // if(ctrl.fakeScheme.Contains("oneZ")) 
    //   ret4l = apply_HZZ4L_Z_selection(ctrl, info,
    //                                hltTable,
    //                                hltObjArr,
    //                                fTrigObjs,
    //                                vtxArr,
    //                                pfArr,
    //                                puDArr,
    //                                electronArr, 
    //                                &electronReferencePreSelection,
    //                                &electronReferenceIDMVASelectionV1,
    //                                &electronReferenceIsoSelection,
    //                                muonArr, 
    //                                &muonReferencePreSelection,
    //                                &muonIDPFSelection,
    //                                &muonReferenceIsoSelection);
    // else
      ret4l = apply_HZZ4L_EMU_selection(ctrl, info,
                                        hltTable,
                                        hltObjArr,
                                        fTrigObjs,
                                        vtxArr,
                                        pfArr,
                                        puDArr,
                                        fConversions,
                                        electronArr, 
                                        &electronReferencePreSelection,
                                        &electronReferenceIDMVASelectionV1,
                                        &electronReferenceIsoSelection,
                                        muonArr, 
                                        &muonReferencePreSelection,
                                        &muonIDPFSelection,
                                        &muonReferenceIsoSelection);
      
    if( ret4l.status.pass() ) { 
      fillKinematics(ret4l,kinematics);
      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());
      foTree.Fill();
      nt.Fill();
    }
    
  }
  foTree.getFile()->cd();
  foTree.getTree()->Write();
  nt.WriteClose();
}
//----------------------------------------------------------------------------
void initRunLumiRangeMap(ControlFlags &ctrl) 
//----------------------------------------------------------------------------
{
  /*
    rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
    //  rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
    rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
    rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));  
    rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));  
    // r11b
    rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));  
  */

  // 2012 only for now ...
  rlrm.AddJSONFile(string(ctrl.jsonFile));

};
Revision:	1.7
Committed:	Mon Jul 9 09:19:11 2012 UTC (12 years, 10 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+17 -18 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	dkralph	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			#include <bitset>
8			#include <map>
9			using namespace std;
10
11			//
12			// root headers
13			//
14			#include <TFile.h>
15			#include <TTree.h>
16			#include <TChain.h>
17			#include <TBranch.h>
18			#include <TClonesArray.h>
19
20			//
21			// TMVA
22			//
23			#include "TMVA/Reader.h"
24			#include "TMVA/Tools.h"
25			#include "TMVA/Config.h"
26			#include "TMVA/MethodBDT.h"
27
28			//
29			// ntuple format headers
30			//
31			#include "EventHeader.h"
32			#include "PFMet.h"
33			#include "Electron.h"
34			#include "Muon.h"
35			#include "Vertex.h"
36	dkralph	1.5	#include "PFJet.h"
37	dkralph	1.1	#include "PFCandidate.h"
38			#include "PFCandidateCol.h"
39			#include "PileupInfoCol.h"
40			#include "PileupEnergyDensity.h"
41			#include "MCParticle.h"
42			#include "TriggerMask.h"
43			#include "TriggerTable.h"
44			#include "Names.h"
45			#include "BaseMod.h"
46
47			//
48			// our headers
49			//
50			#include "ParseArgs.h"
51			#include "MuonSelection.h"
52			#include "ElectronSelection.h"
53			#include "IsolationSelection.h"
54			#include "SelectionEMU.h"
55
56			#include "TriggerUtils.h"
57			#include "PassHLT.h"
58			#include "Angles.h"
59			#include "KinematicsStruct.h"
60			#include "InfoStruct.h"
61			#include "GenInfoStruct.h"
62			#include "WeightStruct.h"
63			//#include "GlobalDataAndFuncs.h"
64			#include "RunLumiRangeMap.h"
65			#include "SelectionFuncs.h"
66
67			#include "AngleTuple.h"
68			#include "FOTuple.h"
69
70			#include "SampleWeight.h"
71			#include "EfficiencyWeightsInterface.h"
72
73			#include "SimpleLepton.h"
74
75			#ifndef CMSSW_BASE
76			#define CMSSW_BASE "../"
77			#endif
78
79			using namespace mithep;
80
81			//
82			// globals
83			//----------------------------------------------------------------------------
84	khahn	1.3	TH1D hpu_2011, hpu_2012;
85	dkralph	1.1	mithep::RunLumiRangeMap rlrm;
86			vector<SimpleLepton> failingLeptons ; // for fake estimation
87			vector<SimpleLepton> passingLeptons; // for fake estimation
88			vector<unsigned> cutvec;
89			vector<vector<unsigned> > zcutvec;
90			vector<vector<unsigned> > zzcutvec;
91			map<unsigned,float> evtrhoMap;
92			vector<string> cutstrs;
93			bool passes_HLT_MC;
94			vector<bool> PFnoPUflag;;
95	dkralph	1.5	vector<int> muTrigObjs,eleTrigObjs,muTriggers,eleTriggers;
96			std::bitset<TRIGGER_BIG_NUMBER> triggerBits;
97	dkralph	1.1
98			//
99			// prototypes
100			//----------------------------------------------------------------------------
101	dkralph	1.4	void initRunLumiRangeMap(ControlFlags &ctrl);
102	dkralph	1.1
103			//
104			// MAIN
105			//----------------------------------------------------------------------------
106			int main(int argc, char** argv)
107			//----------------------------------------------------------------------------
108			{
109
110			string cmsswpath(CMSSW_BASE);
111			cmsswpath.append("/src");
112			map<string,unsigned> entrymap; // number of unskimmed entries in each file
113
114			//
115			// args
116			//--------------------------------------------------------------------------------------------------------------
117			ControlFlags ctrl;
118			parse_args( argc, argv, ctrl );
119			if( ctrl.inputfiles.empty() &&ctrl.inputfile.empty() )
120			{
121			cerr << "usage: applySelection.exe <flags> " << endl;
122			cerr << "\tmandoatory flags : " << endl;
123			cerr << "\t--inputfiles \| file containing a list of ntuples to run over" << endl;
124			cerr << "\t--inputfile \| a file to run over" << endl;
125			cerr << "\t--outputfile \| file to store selected evet" << endl;
126			return 1;
127			}
128			ctrl.dump();
129
130
131
132			//
133			// File I/O
134			//--------------------------------------------------------------------------------------------------------------
135			TChain * chain = new TChain("Events");
136			TChain * hltchain = new TChain("HLT");
137
138			string fname;
139			unsigned total_unskimmed=0;
140			if( !ctrl.inputfiles.empty() ) {
141			ifstream f(ctrl.inputfiles.c_str());
142			while (f >> fname) {
143			if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
144			cout << "adding inputfile : " << fname.c_str() << endl;
145			entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
146			total_unskimmed += entrymap[string(fname.c_str())];
147			chain->AddFile(fname.c_str());
148			hltchain->AddFile(fname.c_str());
149			}
150			} else {
151			cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
152			unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
153			entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
154			total_unskimmed += unsk_ents;
155			chain->AddFile(ctrl.inputfile.c_str());
156			hltchain->AddFile(ctrl.inputfile.c_str());
157			}
158			// write the total number of unskimmed events that went into making this output file to a text file
159			writeEntries(ctrl,total_unskimmed);
160
161			const char * ofname;
162			if( strcmp( ctrl.outputfile.c_str(), "") ) {
163			ofname = ctrl.outputfile.c_str();
164			} else {
165			ofname = "tmp.root";
166			}
167			// table of cross section values
168			string xspath = (cmsswpath+string("/MitPhysics/data/xs.dat"));
169			cout << "xspath: " << xspath.c_str() << endl;
170			SimpleTable xstab(xspath.c_str());
171
172			//
173			// Setup
174			//--------------------------------------------------------------------------------------------------------------
175			Angles angles;
176			KinematicsStruct kinematics;
177			InfoStruct evtinfo;
178			WeightStruct weights;
179			GenInfoStruct geninfo;
180
181			AngleTuple nt( (const char)ofname, (const char)"zznt");
182			nt.makeAngleBranch(angles);
183			nt.makeKinematicsBranch(kinematics);
184			nt.makeInfoBranch(evtinfo);
185
186
187			FOTuple foTree( nt.getFile(), (const char*)"FOtree");
188			foTree.makeFailedLeptonBranch(failingLeptons);
189			foTree.makePassedLeptonBranch(passingLeptons);
190
191			TH1F * h_wf11_hpt;
192			if(ctrl.mc) {
193			nt.makeWeightBranch(weights);
194			if(ctrl.fillGen ) nt.makeGenInfoBranch(geninfo);
195	dkralph	1.5	// initEfficiencyWeights();
196			cout << "\n\nWARNING: not initializing eff weights\n\n" << endl;
197	dkralph	1.1	initPUWeights();
198
199			// Higgs only, pt reweighting
200			if( ctrl.mH <= 140 ) {
201			char wptname[256];
202			sprintf( wptname, "/scratch/dkralph/pharris/Flat/ntupler/root/weight_ptH_%i.root", ctrl.mH );
203			TFile * f_hpt = new TFile(wptname);
204			f_hpt->Print();
205			sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
206			h_wf11_hpt = (TH1F*)(f_hpt->FindObjectAny(wptname));
207			}
208
209			} else {
210	dkralph	1.4	initRunLumiRangeMap(ctrl);
211	dkralph	1.1	}
212
213			// initMuonIDMVA();
214			initElectronIDMVAV1();
215			initMuonIsoMVA();
216			initElectronIDMVA();
217			initElectronIsoMVA();
218			initTrigger();
219
220	dkralph	1.5	muTriggers.push_back(kHLT_IsoMu24_eta2p1); muTrigObjs.push_back(kHLT_IsoMu24_eta2p1_MuObj);
221			eleTriggers.push_back(kHLT_Ele27_WP80); eleTrigObjs.push_back(kHLT_Ele27_WP80_EleObj);
222
223	dkralph	1.1	// hack
224			initEvtRhoMap(evtrhoMap);
225
226
227			//
228			// Setup tree I/O
229			//--------------------------------------------------------------------------------------------------------------
230			string currentFile("");
231
232	dkralph	1.5	UInt_t fNMaxTriggers = TRIGGER_BIG_NUMBER;
233	dkralph	1.1	mithep::EventHeader *info = new mithep::EventHeader();
234			mithep::Array<mithep::PFMet> *metArr = new mithep::Array<mithep::PFMet>();
235			mithep::Array<mithep::Electron> *electronArr = new mithep::Array<mithep::Electron>();
236			mithep::Array<mithep::Muon> *muonArr = new mithep::Array<mithep::Muon>();
237			mithep::Array<mithep::Vertex> *vtxArr = new mithep::Array<mithep::Vertex>();
238			mithep::Array<mithep::PFCandidate> *pfArr = new mithep::Array<mithep::PFCandidate>();
239	dkralph	1.5	mithep::Array<mithep::PFJet> *jetArr = new mithep::Array<mithep::PFJet>();
240	dkralph	1.1	mithep::Array<mithep::PileupInfo> *puArr = new mithep::Array<mithep::PileupInfo>();
241			mithep::Array<mithep::PileupEnergyDensity> *puDArr = new mithep::Array<mithep::PileupEnergyDensity>();
242			mithep::Array<mithep::Track> *trkArr = new mithep::Array<mithep::Track>();
243			mithep::Array<mithep::MCParticle> *mcArr = new mithep::Array<mithep::MCParticle>();
244			mithep::MCEventInfo *mcEvtInfo = new mithep::MCEventInfo();
245			mithep::TriggerMask *trigMask = new mithep::TriggerMask();
246			mithep::TriggerTable *hltTable = new mithep::TriggerTable();
247			vector<string> *hltTableStrings = new vector<string>();
248	dkralph	1.5	vector<string> *hltLabelStrings = new vector<string>();
249			mithep::Array<mithep::TriggerObject> *hltObjArr = new mithep::Array<mithep::TriggerObject>();
250			mithep::Array<mithep::TriggerObjectRel> *hltRelsArr = new mithep::Array<mithep::TriggerObjectRel>();
251			std::vector<std::string> *hltTab = new vector<string>();
252			std::vector<std::string> *hltLab = new vector<string>();
253			mithep::TriggerObjectsTable *fTrigObjs = new TriggerObjectsTable(hltTable,fNMaxTriggers);
254			mithep::Array<mithep::DecayParticle> *fConversions = new mithep::Array<mithep::DecayParticle>();
255	dkralph	1.1
256			TString fTriggerTableName(Names::gkHltTableBrn);
257	dkralph	1.5	chain->SetBranchAddress(Names::gkEvtHeaderBrn, &info);
258			chain->SetBranchAddress("PFMet", &metArr);
259			chain->SetBranchAddress(Names::gkElectronBrn, &electronArr);
260			chain->SetBranchAddress(Names::gkMuonBrn, &muonArr);
261			chain->SetBranchAddress(Names::gkPVBrn, &vtxArr);
262			chain->SetBranchAddress(Names::gkPFCandidatesBrn, &pfArr);
263			chain->SetBranchAddress(Names::gkPFJetBrn, &jetArr);
264	dkralph	1.1	if( ctrl.mc ) {
265	dkralph	1.5	chain->SetBranchAddress(Names::gkPileupInfoBrn, &puArr);
266			chain->SetBranchAddress(Names::gkMCPartBrn, &mcArr);
267			chain->SetBranchAddress(Names::gkMCEvtInfoBrn, &mcEvtInfo);
268	dkralph	1.1	}
269	dkralph	1.5	chain->SetBranchAddress(Names::gkPileupEnergyDensityBrn, &puDArr);
270			chain->SetBranchAddress(Names::gkTrackBrn, &trkArr);
271			chain->SetBranchAddress(Names::gkHltBitBrn, &trigMask);
272			chain->SetBranchAddress(Names::gkHltObjBrn, &hltObjArr);
273			chain->SetBranchAddress("HLTObjectsRelation", &hltRelsArr);
274			chain->SetBranchAddress(Names::gkMvfConversionBrn, &fConversions);
275			hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
276			hltchain->SetBranchAddress(Names::gkHltLabelBrn, &hltLabelStrings);
277	dkralph	1.1
278			mithep::Vertex vtx; // best primary vertex in the event
279
280			// ginfo = NULL;
281			// if( ctrl.mc ) { chain->SetBranchAddress("Gen", &ginfo);}
282
283			// only 1 HLT table / file ???
284			hltchain->GetEntry(0);
285
286			int imax = chain->GetEntries();
287			cout << "nEntries: " << imax << endl;
288
289
290			//
291			// Loop !!!!!!!!!
292			//--------------------------------------------------------------------------------------------------------------
293			for(UInt_t ientry=0; ientry<imax; ientry++) {
294			chain->GetEntry(ientry);
295			if(!(ientry%1000)) cerr << "entry: " << ientry << "\t" << float(ientry)/imax << endl;
296
297			string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
298	dkralph	1.6	if( ctrl.era == 0 )
299			setEra( fname, ctrl );
300	dkralph	1.1
301			// pfNoPU
302			PFnoPUflag.clear();
303			int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
304			assert(pfnopu_size == pfArr->GetEntries());
305
306			// trigger
307			if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
308			currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
309			hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
310			hltchain->GetEntry(0);
311			hltTable->Clear();
312			fillTriggerNames(hltTable, hltTableStrings );
313			}
314			fillTriggerBits( hltTable, trigMask, triggerBits );
315	dkralph	1.5	setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
316	dkralph	1.1
317			//
318			// data/MC
319			//
320			if(ctrl.mc) {
321			if( ctrl.fillGen )
322			fillGenInfo( mcArr, mcEvtInfo, geninfo, ESampleType::kHZZ, ctrl);
323			} else {
324			if(!(ctrl.noJSON) ) {
325	dkralph	1.4	RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
326			if( !(rlrm.HasRunLumi(rl)) ) continue;
327	dkralph	1.1	}
328			// if( !passHLTEMU(ctrl,triggerBits) ) {
329			// continue;
330			// }
331			}
332
333			// lepton/kinematic selection ...
334			failingLeptons.clear();
335			passingLeptons.clear();
336	dkralph	1.6	EventData ret4l;
337	dkralph	1.1
338	dkralph	1.7	// if(ctrl.fakeScheme.Contains("oneZ"))
339			// ret4l = apply_HZZ4L_Z_selection(ctrl, info,
340			// hltTable,
341			// hltObjArr,
342			// fTrigObjs,
343			// vtxArr,
344			// pfArr,
345			// puDArr,
346			// electronArr,
347			// &electronReferencePreSelection,
348			// &electronReferenceIDMVASelectionV1,
349			// &electronReferenceIsoSelection,
350			// muonArr,
351			// &muonReferencePreSelection,
352			// &muonIDPFSelection,
353			// &muonReferenceIsoSelection);
354			// else
355	dkralph	1.6	ret4l = apply_HZZ4L_EMU_selection(ctrl, info,
356			hltTable,
357			hltObjArr,
358			fTrigObjs,
359			vtxArr,
360			pfArr,
361			puDArr,
362			fConversions,
363			electronArr,
364			&electronReferencePreSelection,
365			&electronReferenceIDMVASelectionV1,
366			&electronReferenceIsoSelection,
367			muonArr,
368			&muonReferencePreSelection,
369			&muonIDPFSelection,
370			&muonReferenceIsoSelection);
371	dkralph	1.1
372			if( ret4l.status.pass() ) {
373	dkralph	1.6	fillKinematics(ret4l,kinematics);
374	dkralph	1.1	TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
375	dkralph	1.5	fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries());
376	dkralph	1.1	foTree.Fill();
377			nt.Fill();
378			}
379	dkralph	1.6
380	dkralph	1.1	}
381			foTree.getFile()->cd();
382			foTree.getTree()->Write();
383			nt.WriteClose();
384			}
385			//----------------------------------------------------------------------------
386	dkralph	1.4	void initRunLumiRangeMap(ControlFlags &ctrl)
387	dkralph	1.1	//----------------------------------------------------------------------------
388			{
389			/*
390			rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
391			// rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
392			rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
393			rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));
394			rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));
395			// r11b
396			rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));
397			*/
398
399			// 2012 only for now ...
400	dkralph	1.4	rlrm.AddJSONFile(string(ctrl.jsonFile));
401	dkralph	1.1
402			};