--- UserCode/kiesel/TreeWriter/treeWriter.cc 2013/04/18 13:38:55 1.15 +++ UserCode/kiesel/TreeWriter/treeWriter.cc 2013/04/25 17:37:32 1.26 @@ -16,7 +16,6 @@ TreeWriter::TreeWriter( TChain* inputTre Init( outputName, loggingVerbosity_ ); } - void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) { if (loggingVerbosity_ > 0) @@ -57,12 +56,35 @@ TreeWriter::~TreeWriter() { delete tree; } +template +std::vector* getVectorFromMap( map > myMap, TString search ){ + // if nothing is found, return a empty vector (which has to be deleted never the less) + std::vector* vec; + typename map >::iterator mapIt = myMap.find( search ); + if( mapIt == myMap.end() ) { + cout << "ERROR: Collection \"" << search << "\" not found!" << endl; + vec = new std::vector; + } else + vec = &mapIt->second; + cout << vec->size() << endl; + return vec; +} + +float deltaPhi( float phi1, float phi2) { + float result = phi1 - phi2; + while (result > M_PI) result -= 2*M_PI; + while (result <= -M_PI) result += 2*M_PI; + return result; +} + // useful functions -float deltaR( TLorentzVector v1, TLorentzVector v2 ) { - return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) ); +float deltaR( const TLorentzVector& v1, const TLorentzVector& v2 ) { + // deltaR = sqrt ( deltaEta^2 + deltaPhi^2 ) + return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(deltaPhi(v1.Phi(),v2.Phi()), 2) ); } float effectiveAreaElectron( float eta ) { + // needed by calculating the isolation for electrons // see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection // only for Delta R = 0.3 on 2012 Data eta = fabs( eta ); @@ -78,7 +100,7 @@ float effectiveAreaElectron( float eta ) } // correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012 -float chargedHadronIso_corrected(susy::Photon gamma, float rho) { +float chargedHadronIso_corrected(const susy::Photon& gamma, float rho) { float eta = fabs(gamma.caloPosition.Eta()); float ea; @@ -96,7 +118,7 @@ float chargedHadronIso_corrected(susy::P return iso; } -float neutralHadronIso_corrected(susy::Photon gamma, float rho) { +float neutralHadronIso_corrected(const susy::Photon& gamma, float rho) { float eta = fabs(gamma.caloPosition.Eta()); float ea; @@ -114,7 +136,7 @@ float neutralHadronIso_corrected(susy::P return iso; } -float photonIso_corrected(susy::Photon gamma, float rho) { +float photonIso_corrected(const susy::Photon& gamma, float rho) { float eta = fabs(gamma.caloPosition.Eta()); float ea; @@ -133,6 +155,7 @@ float photonIso_corrected(susy::Photon g } float d0correction( const susy::Electron& electron, const susy::Event& event ) { + // copied from Brian Francis TVector3 beamspot = event.vertices[0].position; susy::Track track = event.tracks[electron.gsfTrackIndex]; float d0 = track.d0() - beamspot.X()*sin(track.phi()) + beamspot.Y()*cos(track.phi()); @@ -140,6 +163,7 @@ float d0correction( const susy::Electron } float dZcorrection( const susy::Electron& electron, const susy::Event& event ) { + // copied from Brian Francis TVector3 beamspot = event.vertices[0].position; susy::Track track = event.tracks[electron.gsfTrackIndex]; @@ -148,7 +172,7 @@ float dZcorrection( const susy::Electron return dz; } -float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) { +float getPtFromMatchedJet( const susy::Photon& myPhoton, const susy::PFJetCollection& jetColl, int loggingVerbosity = 0 ) { /** * \brief Takes jet p_T as photon p_T * @@ -160,14 +184,16 @@ float getPtFromMatchedJet( susy::Photon std::vector nearJets; nearJets.clear(); - for(std::vector::iterator it = jetColl.begin(); + for(std::vector::const_iterator it = jetColl.begin(); it != jetColl.end(); ++it) { - std::map::iterator s_it = it->jecScaleFactors.find("L2L3"); + float scale = 1.; + std::map::const_iterator s_it = it->jecScaleFactors.find("L2L3"); if (s_it == it->jecScaleFactors.end()) { std::cout << "JEC is not available for this jet!!!" << std::endl; continue; + } else { + scale = s_it->second; } - float scale = s_it->second; TLorentzVector corrP4 = scale * it->momentum; float deltaR_ = deltaR(myPhoton.momentum, corrP4 ); if (deltaR_ > 0.3) continue; @@ -202,6 +228,13 @@ float getPtFromMatchedJet( susy::Photon void TreeWriter::Loop() { + /** + * \brief Loops over input chain and fills tree + * + * This is the major function of treeWriter, which initialize the output, loops + * over all events and fill the tree. In the end, the tree is saved to the + * output File + */ // here the event loop is implemented and the tree is filled if (inputTree == 0) return; @@ -209,7 +242,7 @@ void TreeWriter::Loop() { // get number of events to be proceeded Long64_t nentries = inputTree->GetEntries(); // store them in histo - eventNumbers->Fill( "Number of generated Events", nentries ); + eventNumbers->Fill( "Number of generated events", nentries ); if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries; if( loggingVerbosity > 0 ) @@ -227,16 +260,21 @@ void TreeWriter::Loop() { tree->Branch("ht", &ht, "ht/F"); tree->Branch("nVertex", &nVertex, "nVertex/I"); tree->Branch("pu_weight", &pu_weight, "pu_weight/F"); + tree->Branch("genElectron", &genElectron); + tree->Branch("genPhoton", &genPhoton); - - for (long jentry=0; jentry < processNEvents; ++jentry) { - if ( loggingVerbosity>2 || jentry%reportEvery==0 ) + for (unsigned long jentry=0; jentry < processNEvents; ++jentry) { + if ( loggingVerbosity>1 || jentry%reportEvery==0 ) std::cout << jentry << " / " << processNEvents << std :: endl; + inputTree->LoadTree( jentry ); + inputTree->GetEntry( jentry ); photon.clear(); jet.clear(); electron.clear(); muon.clear(); + genElectron.clear(); + genPhoton.clear(); ht = 0; // weights @@ -253,24 +291,14 @@ void TreeWriter::Loop() { } // get ak5 jets - susy::PFJetCollection jetVector; - if(event->pfJets.count("ak5") == 0) - cout << "ERROR: Jet collection 'ak5' not found" << endl; - else - jetVector = event->pfJets.find("ak5")->second; + std::vector jetVector = event->pfJets["ak5"]; // photons - if( loggingVerbosity > 1 ) - std::cout << "Process photons" << std::endl; - susy::PhotonCollection photonVector; - if(event->photons.count("photons") == 0) - cout << "ERROR: Photon collection 'photons' not found" << endl; - else - photonVector = event->photons.find("photons")->second; + std::vector photonVector = event->photons["photons"]; - for(susy::PhotonCollection :: iterator it = photonVector.begin(); + for(std::vector::iterator it = photonVector.begin(); it != photonVector.end(); ++it ) { - if( !(it->isEB() || it->isEE()) && skim ) + if( !(it->isEE() || it->isEB()) && it->isEBEtaGap() && it->isEBPhiGap() && it->isEERingGap() && it->isEEDeeGap() && it->isEBEEGap() && skim ) continue; tree::Photon thisphoton; thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity ); @@ -294,6 +322,7 @@ void TreeWriter::Loop() { && it->sigmaIetaIeta<0.034 && thisphoton.chargedIso<2.3 && thisphoton.neutralIso<2.9+0.04*thisphoton.pt; + if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim ) continue; thisphoton.eta = it->momentum.Eta(); @@ -315,61 +344,58 @@ void TreeWriter::Loop() { std::cout << "Found " << photon.size() << " photons" << std::endl; // jets - std::map::iterator pfJets_it = event->pfJets.find("ak5"); - if(pfJets_it == event->pfJets.end()){ - if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl; - } else { - susy::PFJetCollection& jetColl = pfJets_it->second; - for(std::vector::iterator it = jetColl.begin(); - it != jetColl.end(); ++it) { - tree::Jet thisjet; - /*std::map::iterator s_it = it->jecScaleFactors.find("L2L3"); - if (s_it == it->jecScaleFactors.end()) { - std::cout << "JEC is not available for this jet!!!" << std::endl; - continue; - } - float scale = s_it->second; - */ - float scale = 1.; - if(it->jecScaleFactors.count("L2L3") == 0) - std::cout << "ERROR: JEC is not available for this jet" << std::endl; - else - scale = it->jecScaleFactors.find("L2L3")->second; - TLorentzVector corrP4 = scale * it->momentum; - - if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue; - if(corrP4.Pt() < 30 && skim ) continue; - thisjet.pt = corrP4.Pt(); - thisjet.eta = corrP4.Eta(); - thisjet.phi = corrP4.Phi(); - thisjet.bCSV = it->bTagDiscriminators[susy::kCSV]; - // jet composition - thisjet.chargedHadronEnergy = it->chargedHadronEnergy; - thisjet.neutralHadronEnergy = it->neutralHadronEnergy; - thisjet.photonEnergy = it->photonEnergy; - thisjet.electronEnergy = it->electronEnergy; - thisjet.muonEnergy = it->muonEnergy; - thisjet.HFHadronEnergy = it->HFHadronEnergy; - thisjet.HFEMEnergy = it->HFEMEnergy; - thisjet.chargedEmEnergy = it->chargedEmEnergy; - thisjet.chargedMuEnergy = it->chargedMuEnergy; - thisjet.neutralEmEnergy = it->neutralEmEnergy; + for(std::vector::iterator it = jetVector.begin(); + it != jetVector.end(); ++it) { + tree::Jet thisjet; + + // scale with JEC + float scale = 1.; + if(it->jecScaleFactors.count("L2L3") == 0) + std::cout << "ERROR: JEC is not available for this jet" << std::endl; + else + scale = it->jecScaleFactors.find("L2L3")->second; + TLorentzVector corrP4 = scale * it->momentum; + + // Calculate HT. + // The definiton differs from the saved jet, since trigger is described better + if( std::abs( corrP4.Eta() ) < 3 && corrP4.Pt() > 40 ) + ht += thisjet.pt; - if( loggingVerbosity > 2 ) - std::cout << " p_T, jet = " << thisjet.pt << std::endl; + if(std::abs(corrP4.Eta()) > 2.6 && skim ) continue; + if(corrP4.Pt() < 30 && skim ) continue; + thisjet.pt = corrP4.Pt(); + thisjet.eta = corrP4.Eta(); + thisjet.phi = corrP4.Phi(); + thisjet.bCSV = it->bTagDiscriminators[susy::kCSV]; + // jet composition + thisjet.chargedHadronEnergy = it->chargedHadronEnergy; + thisjet.neutralHadronEnergy = it->neutralHadronEnergy; + thisjet.photonEnergy = it->photonEnergy; + thisjet.electronEnergy = it->electronEnergy; + thisjet.muonEnergy = it->muonEnergy; + thisjet.HFHadronEnergy = it->HFHadronEnergy; + thisjet.HFEMEnergy = it->HFEMEnergy; + thisjet.chargedEmEnergy = it->chargedEmEnergy; + thisjet.chargedMuEnergy = it->chargedMuEnergy; + thisjet.neutralEmEnergy = it->neutralEmEnergy; - jet.push_back( thisjet ); - ht += thisjet.pt; - }// for jet - }// if, else + if( loggingVerbosity > 2 ) + std::cout << " p_T, jet = " << thisjet.pt << std::endl; + + jet.push_back( thisjet ); + }// for jet if( jet.size() < 2 && skim ) continue; std::sort( jet.begin(), jet.end(), tree::EtGreater); if( loggingVerbosity > 1 ) std::cout << "Found " << jet.size() << " jets" << std::endl; + if( ht < 450 && skim) + continue; + + // met std::map::iterator met_it = event->metMap.find("pfMet"); @@ -387,14 +413,8 @@ void TreeWriter::Loop() { std::cout << " type1met = " << type1met << std::endl; // electrons - std::vector* eVector; - if( event->electrons.count("gsfElectrons") == 0) { - cout << "EROR: no electron collection found" << endl; - continue; - } else - eVector = &event->electrons.find("gsfElectrons")->second; - - for(vector::iterator it = eVector->begin(); it < eVector->end(); ++it) { + std::vector eVector = event->electrons["gsfElectrons"]; + for(std::vector::iterator it = eVector.begin(); it < eVector.end(); ++it) { tree::Particle thiselectron; if( loggingVerbosity > 2 ) cout << " electron pt = " << it->momentum.Pt() << endl; @@ -402,58 +422,45 @@ void TreeWriter::Loop() { // use veto electrons if( it->momentum.Pt() < 20 || it->momentum.Pt() > 1e6 ) continue; // spike rejection - float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso - - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. ) + float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (float)0. ) ) / it->momentum.Pt(); - cout << iso << endl; float d0 = d0correction( *it, *event ); - cout << d0 << endl; float dZ = std::abs( dZcorrection( *it, *event ) ); - cout << dZ << endl; if ( it->isEB() ){ - cout << "electron in barrel" << endl; if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007 || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8 || it->sigmaIetaIeta > 0.01 || it->hcalOverEcalBc > 0.15 || d0 > 0.04 || dZ > 0.2 - || iso > 0.15 ){ - cout << " no real electron" << endl; + || iso > 0.15 ) continue; } - } else if( it->isEE() ) { - cout << "electron in endcap" << endl; if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01 || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7 || it->sigmaIetaIeta > 0.03 || d0 > 0.04 || dZ > 0.2 - || iso > 0.15 ){ - cout << "no real electron" << endl; + || iso > 0.15 ) continue; } - } - else{ // not in barrel nor in endcap - cout << " electron in gap" << endl; + else // not in barrel nor in endcap continue; - // TODO: conversion rejection information not implemented yet, see twiki for more details - } + thiselectron.pt = it->momentum.Pt(); if( loggingVerbosity > 2 ) std::cout << " p_T, electron = " << it->momentum.Et() << std::endl; thiselectron.eta = it->momentum.Eta(); thiselectron.phi = it->momentum.Phi(); - cout << " adde electron" << endl; electron.push_back( thiselectron ); } if( loggingVerbosity > 1 ) std::cout << "Found " << electron.size() << " electrons" << std::endl; // muons - std::vector mVector = event->muons["muons"]; tree::Particle thismuon; + std::vector mVector = event->muons["muons"]; for( std::vector::iterator it = mVector.begin(); it != mVector.end(); ++it) { if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) ) continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon @@ -467,13 +474,26 @@ void TreeWriter::Loop() { if( loggingVerbosity > 1 ) std::cout << "Found " << muon.size() << " muons" << std::endl; - // vertices nVertex = event->vertices.size(); - if( ht < 450 && skim) - continue; - + tree::Particle thisGenParticle; + for( std::vector::iterator it = event->genParticles.begin(); it != event->genParticles.end(); ++it ) { + if( it->momentum.Pt() < 20 ) continue; + thisGenParticle.pt = it->momentum.Pt(); + thisGenParticle.eta = it->momentum.Eta(); + thisGenParticle.phi = it->momentum.Phi(); + switch( std::abs(it->pdgId) ) { + case 22: // photon + genPhoton.push_back( thisGenParticle ); + break; + case 11: // electron + // Demand a W boson as mother particle of electron + if( abs(event->genParticles[it->motherIndex].pdgId) == 24 ) + genElectron.push_back( thisGenParticle ); + break; + } + } tree->Fill(); } // for jentry