kiesel/TreeWriter/treeWriter.cc

#include "treeWriter.h"

using namespace std;

TreeWriter::TreeWriter( std::string inputName, std::string outputName, int loggingVerbosity_ ) {
        // read the input file
        inputTree = new TChain("susyTree");
        if (loggingVerbosity_ > 0)
                std::cout << "Add files to chain" << std::endl;
        inputTree->Add( inputName.c_str() );
        Init( outputName, loggingVerbosity_ );
}

TreeWriter::TreeWriter( TChain* inputTree_, std::string outputName, int loggingVerbosity_ ) {
        inputTree = inputTree_;
        Init( outputName, loggingVerbosity_ );
}


void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) {

        if (loggingVerbosity_ > 0)
                std::cout << "Set Branch Address of susy::Event" << std::endl;
        event = new susy::Event;
        inputTree->SetBranchAddress("susyEvent", &event);

        // Here the number of proceeded events will be stored. For plotting, simply use L*sigma/eventNumber
        eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated events", 1, 0, 1);
        eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated events");

        // open the output file
        if (loggingVerbosity_>0)
                std::cout << "Open file " << outputName << " for writing." << std::endl;
        outFile = new TFile( outputName.c_str(), "recreate" );
        tree = new TTree("susyTree","Tree for single photon analysis");

        // set default parameter
        processNEvents = -1;
        reportEvery = 1000;
        loggingVerbosity = loggingVerbosity_;
        skim = true;
        pileupHisto = 0;
}

void TreeWriter::PileUpWeightFile( string pileupFileName ) {
        TFile *puFile = new TFile( pileupFileName.c_str() );
        pileupHisto = (TH1F*) puFile->Get("pileup");
}

TreeWriter::~TreeWriter() {
        if (pileupHisto != 0 )
                delete pileupHisto;
        inputTree->GetCurrentFile()->Close();
        delete inputTree;
        delete event;
        delete outFile;
        delete tree;
}

// useful functions
float deltaR( TLorentzVector v1, TLorentzVector v2 ) {
        return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
}

float effectiveAreaElectron( float eta ) {
        // see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
        // only for Delta R = 0.3 on 2012 Data
        eta = fabs( eta );
        float ea;
        if( eta < 1.0 ) ea = 0.13;
        else if( eta < 1.479 ) ea = 0.14;
        else if( eta < 2.0 ) ea = 0.07;
        else if( eta < 2.2 ) ea = 0.09;
        else if( eta < 2.3 ) ea = 0.11;
        else if( eta < 2.4 ) ea = 0.11;
        else ea = 0.14;
        return ea;
}

// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
        float eta = fabs(gamma.caloPosition.Eta());
        float ea;

        if(eta < 1.0) ea = 0.012;
        else if(eta < 1.479) ea = 0.010;
        else if(eta < 2.0) ea = 0.014;
        else if(eta < 2.2) ea = 0.012;
        else if(eta < 2.3) ea = 0.016;
        else if(eta < 2.4) ea = 0.020;
        else ea = 0.012;

        float iso = gamma.chargedHadronIso;
        iso = max(iso - rho*ea, (float)0.);

        return iso;
}

float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
        float eta = fabs(gamma.caloPosition.Eta());
        float ea;

        if(eta < 1.0) ea = 0.030;
        else if(eta < 1.479) ea = 0.057;
        else if(eta < 2.0) ea = 0.039;
        else if(eta < 2.2) ea = 0.015;
        else if(eta < 2.3) ea = 0.024;
        else if(eta < 2.4) ea = 0.039;
        else ea = 0.072;

        float iso = gamma.neutralHadronIso;
        iso = max(iso - rho*ea, (float)0.);

        return iso;
}

float photonIso_corrected(susy::Photon gamma, float rho) {
        float eta = fabs(gamma.caloPosition.Eta());
        float ea;

        if(eta < 1.0) ea = 0.148;
        else if(eta < 1.479) ea = 0.130;
        else if(eta < 2.0) ea = 0.112;
        else if(eta < 2.2) ea = 0.216;
        else if(eta < 2.3) ea = 0.262;
        else if(eta < 2.4) ea = 0.260;
        else ea = 0.266;

        float iso = gamma.photonIso;
        iso = max(iso - rho*ea, (float)0.);

        return iso;
}

float d0correction( const susy::Electron& electron, const susy::Event& event ) {
        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());
        return d0;
}

float dZcorrection( const susy::Electron& electron, const susy::Event& event ) {
        TVector3 beamspot = event.vertices[0].position;
        susy::Track track = event.tracks[electron.gsfTrackIndex];

        if(track.momentum.Pt() == 0.) return 1.e6;
        float dz = (track.vertex.Z() - beamspot.Z()) - ((track.vertex.X() - beamspot.X())*track.momentum.Px() + (track.vertex.Y() - beamspot.Y())*track.momentum.Py()) / track.momentum.Pt() * (track.momentum.Pz() / track.momentum.Pt());
        return dz;
}

float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) {
        /**
         * \brief Takes jet p_T as photon p_T
         *
         * At first all jets with DeltaR < 0.3 (isolation cone) are searched.
         * If several jets are found, take the one with the minimal pt difference
         * compared to the photon. If no such jets are found, keep the photon_pt
         * TODO: remove photon matched jet from jet-selection?
         */
        std::vector<susy::PFJet> nearJets;
        nearJets.clear();

        for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
                        it != jetColl.end(); ++it) {
                std::map<TString,Float_t>::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;
                TLorentzVector corrP4 = scale * it->momentum;
                float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
                if (deltaR_ > 0.3) continue;
                if( loggingVerbosity > 2 )
                        std::cout << " pT_jet / pT_gamma = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
                nearJets.push_back( *it );
        }// for jet

        if ( nearJets.size() == 0 ) {
                if( loggingVerbosity > 1 )
                        std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
                return myPhoton.momentum.Et();
        }

        float pt = 0;
        float minPtDifferenz = 1E20; // should be very high
        for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
                        it != jetEnd; ++it ) {
                float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
                if (  ptDiff < minPtDifferenz ) {
                        minPtDifferenz = ptDiff;
                        pt = it->momentum.Et();
                }
        }

        // testing
        if( nearJets.size() > 1 && loggingVerbosity > 0 )
                std::cout << "There are several jets matching to this photon. "
                                        << "Please check if jet-matching is correct." << std::endl;
        return pt;
}


void TreeWriter::Loop() {

        // here the event loop is implemented and the tree is filled
        if (inputTree == 0) return;

        // get number of events to be proceeded
        Long64_t nentries = inputTree->GetEntries();
        // store them in histo
        eventNumbers->Fill( "Number of generated Events", nentries );
        if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;

        if( loggingVerbosity > 0 )
                std::cout << "Processing " << processNEvents << " ouf of "
                        << nentries << " events. " << std::endl;

        tree->Branch("photon", &photon);
        tree->Branch("jet", &jet);
        tree->Branch("electron", &electron);
        tree->Branch("muon", &muon);
        tree->Branch("met", &met, "met/F");
        tree->Branch("metPhi", &met_phi, "metPhi/F");
        tree->Branch("type1met", &type1met, "type1met/F");
        tree->Branch("type1metPhi", &type1met_phi, "type1metPhi/F");
        tree->Branch("ht", &ht, "ht/F");
        tree->Branch("nVertex", &nVertex, "nVertex/I");
        tree->Branch("pu_weight", &pu_weight, "pu_weight/F");


        for (long jentry=0; jentry < processNEvents; ++jentry) {
                if ( loggingVerbosity>2 || jentry%reportEvery==0 )
                        std::cout << jentry << " / " << processNEvents << std :: endl;

                photon.clear();
                jet.clear();
                electron.clear();
                muon.clear();
                ht = 0;

                // weights
                if (pileupHisto == 0) {
                        pu_weight = 1.;
                } else {
                        float trueNumInteractions = -1;
                        for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
                                        iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
                                if (iBX->BX == 0)
                                        trueNumInteractions = iBX->trueNumInteractions;
                        }
                        pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
                }

                // 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;

                // 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;

                for(susy::PhotonCollection :: iterator it = photonVector.begin();
                                it != photonVector.end(); ++it ) {
                        if( !(it->isEB() || it->isEE()) && skim )
                                continue;
                        tree::Photon thisphoton;
                        thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );

                        thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
                        thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
                        thisphoton.photonIso = photonIso_corrected(*it, event->rho25);

                        bool loose_photon_barrel = thisphoton.pt>20
                                && it->isEB()
                                && it->passelectronveto
                                && it->hadTowOverEm<0.05
                                && it->sigmaIetaIeta<0.012
                                && thisphoton.chargedIso<2.6
                                && thisphoton.neutralIso<3.5+0.04*thisphoton.pt
                                && thisphoton.photonIso<1.3+0.005*thisphoton.pt;
                        bool loose_photon_endcap = thisphoton.pt > 20
                                && it->isEE()
                                && it->passelectronveto
                                && it->hadTowOverEm<0.05
                                && 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();
                        thisphoton.phi = it->momentum.Phi();
                        thisphoton.r9 = it->r9;
                        thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
                        thisphoton.hadTowOverEm = it->hadTowOverEm;
                        thisphoton.pixelseed = it->nPixelSeeds;
                        thisphoton.conversionSafeVeto = it->passelectronveto;
                        photon.push_back( thisphoton );
                        if( loggingVerbosity > 2 )
                                std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
                }

                if( photon.size() == 0 && skim )
                        continue;
                std::sort( photon.begin(), photon.end(), tree::EtGreater);
                if( loggingVerbosity > 1 )
                        std::cout << "Found " << photon.size() << " photons" << std::endl;

                // jets
                std::map<TString,susy::PFJetCollection>::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<susy::PFJet>::iterator it = jetColl.begin();
                                        it != jetColl.end(); ++it) {
                                tree::Jet thisjet;
                                /*std::map<TString,Float_t>::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;

                                if( loggingVerbosity > 2 )
                                        std::cout << " p_T, jet = " << thisjet.pt << std::endl;

                                jet.push_back( thisjet );
                                ht += thisjet.pt;
                        }// for jet
                }// if, else
                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;


                // met
                std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
                susy::MET* metobj = &(met_it->second);
                met = metobj->met();
                met_phi = metobj->mEt.Phi();
                if( loggingVerbosity > 2 )
                        std::cout << " met = " << met << std::endl;

                std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
                susy::MET* type1metobj = &(type1met_it->second);
                type1met = type1metobj->met();
                type1met_phi = type1metobj->mEt.Phi();
                if( loggingVerbosity > 2 )
                        std::cout << " type1met = " << type1met << std::endl;

                // electrons
                std::vector<susy::Electron>* eVector;
                if( event->electrons.count("gsfElectrons") == 0) {
                        cout << "EROR: no electron collection found" << endl;
                        continue;
                } else
                        eVector = &event->electrons.find("gsfElectrons")->second;

                for(vector<susy::Electron>::iterator it = eVector->begin(); it < eVector->end(); ++it) {
                        tree::Particle thiselectron;
                        if( loggingVerbosity > 2 )
                                cout << " electron pt = " << it->momentum.Pt() << endl;
                        // for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
                        // 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. )
                                                ) / 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;
                                        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;
                                        continue;
                                }
                                }
                        else{ // not in barrel nor in endcap
                                cout << " electron in gap" << endl;
                                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<susy::Muon> mVector = event->muons["muons"];
                tree::Particle thismuon;
                for( std::vector<susy::Muon>::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
                        thismuon.pt = it->momentum.Et();
                        if( thismuon.pt < 20 )
                                continue;
                        thismuon.eta = it->momentum.Eta();
                        thismuon.phi = it->momentum.Phi();
                        muon.push_back( thismuon );
                }
                if( loggingVerbosity > 1 )
                        std::cout << "Found " << muon.size() << " muons" << std::endl;


                // vertices
                nVertex = event->vertices.size();

                if( ht < 450 && skim)
                        continue;


                tree->Fill();
        } // for jentry


        outFile->cd();
        eventNumbers->Write();
        tree->Write();
        outFile->Write();
        outFile->Close();
}

Revision:	1.15
Committed:	Thu Apr 18 13:38:55 2013 UTC (12 years ago) by kiesel
Content type:	text/plain
Branch:	MAIN
Changes since 1.14:	+158 -143 lines
Log Message:	treeWriter can now be compiled with g++
#	User	Rev	Content
1	kiesel	1.5	#include "treeWriter.h"
2
3			using namespace std;
4	kiesel	1.1
5	kiesel	1.15	TreeWriter::TreeWriter( std::string inputName, std::string outputName, int loggingVerbosity_ ) {
6	kiesel	1.1	// read the input file
7	kiesel	1.5	inputTree = new TChain("susyTree");
8	kiesel	1.6	if (loggingVerbosity_ > 0)
9			std::cout << "Add files to chain" << std::endl;
10	kiesel	1.15	inputTree->Add( inputName.c_str() );
11	kiesel	1.10	Init( outputName, loggingVerbosity_ );
12			}
13
14	kiesel	1.15	TreeWriter::TreeWriter( TChain* inputTree_, std::string outputName, int loggingVerbosity_ ) {
15	kiesel	1.10	inputTree = inputTree_;
16			Init( outputName, loggingVerbosity_ );
17			}
18
19
20	kiesel	1.15	void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) {
21	kiesel	1.5
22	kiesel	1.6	if (loggingVerbosity_ > 0)
23			std::cout << "Set Branch Address of susy::Event" << std::endl;
24	kiesel	1.1	event = new susy::Event;
25			inputTree->SetBranchAddress("susyEvent", &event);
26
27	kiesel	1.14	// Here the number of proceeded events will be stored. For plotting, simply use L*sigma/eventNumber
28	kiesel	1.15	eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated events", 1, 0, 1);
29			eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated events");
30	kiesel	1.14
31	kiesel	1.9	// open the output file
32	kiesel	1.7	if (loggingVerbosity_>0)
33			std::cout << "Open file " << outputName << " for writing." << std::endl;
34	kiesel	1.15	outFile = new TFile( outputName.c_str(), "recreate" );
35	kiesel	1.1	tree = new TTree("susyTree","Tree for single photon analysis");
36
37			// set default parameter
38			processNEvents = -1;
39			reportEvery = 1000;
40	kiesel	1.6	loggingVerbosity = loggingVerbosity_;
41	kiesel	1.5	skim = true;
42	kiesel	1.9	pileupHisto = 0;
43			}
44	kiesel	1.1
45	kiesel	1.9	void TreeWriter::PileUpWeightFile( string pileupFileName ) {
46			TFile *puFile = new TFile( pileupFileName.c_str() );
47			pileupHisto = (TH1F*) puFile->Get("pileup");
48	kiesel	1.1	}
49
50			TreeWriter::~TreeWriter() {
51	kiesel	1.9	if (pileupHisto != 0 )
52			delete pileupHisto;
53			inputTree->GetCurrentFile()->Close();
54	kiesel	1.15	delete inputTree;
55			delete event;
56			delete outFile;
57			delete tree;
58	kiesel	1.1	}
59
60	kiesel	1.8	// useful functions
61	kiesel	1.15	float deltaR( TLorentzVector v1, TLorentzVector v2 ) {
62	kiesel	1.4	return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
63			}
64
65	kiesel	1.12	float effectiveAreaElectron( float eta ) {
66			// see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
67			// only for Delta R = 0.3 on 2012 Data
68			eta = fabs( eta );
69			float ea;
70			if( eta < 1.0 ) ea = 0.13;
71			else if( eta < 1.479 ) ea = 0.14;
72			else if( eta < 2.0 ) ea = 0.07;
73			else if( eta < 2.2 ) ea = 0.09;
74			else if( eta < 2.3 ) ea = 0.11;
75			else if( eta < 2.4 ) ea = 0.11;
76			else ea = 0.14;
77			return ea;
78			}
79
80	kiesel	1.8	// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
81			float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
82	kiesel	1.12	float eta = fabs(gamma.caloPosition.Eta());
83			float ea;
84	kiesel	1.8
85	kiesel	1.12	if(eta < 1.0) ea = 0.012;
86			else if(eta < 1.479) ea = 0.010;
87			else if(eta < 2.0) ea = 0.014;
88			else if(eta < 2.2) ea = 0.012;
89			else if(eta < 2.3) ea = 0.016;
90			else if(eta < 2.4) ea = 0.020;
91			else ea = 0.012;
92	kiesel	1.8
93	kiesel	1.12	float iso = gamma.chargedHadronIso;
94			iso = max(iso - rho*ea, (float)0.);
95	kiesel	1.8
96	kiesel	1.12	return iso;
97	kiesel	1.8	}
98
99			float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
100	kiesel	1.12	float eta = fabs(gamma.caloPosition.Eta());
101			float ea;
102	kiesel	1.8
103	kiesel	1.12	if(eta < 1.0) ea = 0.030;
104			else if(eta < 1.479) ea = 0.057;
105			else if(eta < 2.0) ea = 0.039;
106			else if(eta < 2.2) ea = 0.015;
107			else if(eta < 2.3) ea = 0.024;
108			else if(eta < 2.4) ea = 0.039;
109			else ea = 0.072;
110	kiesel	1.8
111	kiesel	1.12	float iso = gamma.neutralHadronIso;
112			iso = max(iso - rho*ea, (float)0.);
113	kiesel	1.8
114	kiesel	1.12	return iso;
115	kiesel	1.8	}
116
117			float photonIso_corrected(susy::Photon gamma, float rho) {
118	kiesel	1.12	float eta = fabs(gamma.caloPosition.Eta());
119			float ea;
120	kiesel	1.8
121	kiesel	1.12	if(eta < 1.0) ea = 0.148;
122			else if(eta < 1.479) ea = 0.130;
123			else if(eta < 2.0) ea = 0.112;
124			else if(eta < 2.2) ea = 0.216;
125			else if(eta < 2.3) ea = 0.262;
126			else if(eta < 2.4) ea = 0.260;
127			else ea = 0.266;
128	kiesel	1.8
129	kiesel	1.12	float iso = gamma.photonIso;
130			iso = max(iso - rho*ea, (float)0.);
131	kiesel	1.8
132	kiesel	1.12	return iso;
133	kiesel	1.8	}
134
135	kiesel	1.15	float d0correction( const susy::Electron& electron, const susy::Event& event ) {
136	kiesel	1.14	TVector3 beamspot = event.vertices[0].position;
137			susy::Track track = event.tracks[electron.gsfTrackIndex];
138			float d0 = track.d0() - beamspot.X()sin(track.phi()) + beamspot.Y()cos(track.phi());
139			return d0;
140			}
141
142	kiesel	1.15	float dZcorrection( const susy::Electron& electron, const susy::Event& event ) {
143	kiesel	1.14	TVector3 beamspot = event.vertices[0].position;
144			susy::Track track = event.tracks[electron.gsfTrackIndex];
145
146			if(track.momentum.Pt() == 0.) return 1.e6;
147			float dz = (track.vertex.Z() - beamspot.Z()) - ((track.vertex.X() - beamspot.X())track.momentum.Px() + (track.vertex.Y() - beamspot.Y())track.momentum.Py()) / track.momentum.Pt() * (track.momentum.Pz() / track.momentum.Pt());
148			return dz;
149			}
150
151	kiesel	1.15	float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) {
152	kiesel	1.4	/**
153			* \brief Takes jet p_T as photon p_T
154			*
155			* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
156			* If several jets are found, take the one with the minimal pt difference
157	kiesel	1.5	* compared to the photon. If no such jets are found, keep the photon_pt
158	kiesel	1.4	* TODO: remove photon matched jet from jet-selection?
159			*/
160			std::vector<susy::PFJet> nearJets;
161			nearJets.clear();
162
163	kiesel	1.15	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
164			it != jetColl.end(); ++it) {
165			std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
166			if (s_it == it->jecScaleFactors.end()) {
167			std::cout << "JEC is not available for this jet!!!" << std::endl;
168			continue;
169			}
170			float scale = s_it->second;
171			TLorentzVector corrP4 = scale * it->momentum;
172			float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
173			if (deltaR_ > 0.3) continue;
174			if( loggingVerbosity > 2 )
175			std::cout << " pT_jet / pT_gamma = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
176			nearJets.push_back( *it );
177			}// for jet
178	kiesel	1.4
179			if ( nearJets.size() == 0 ) {
180	kiesel	1.12	if( loggingVerbosity > 1 )
181			std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
182	kiesel	1.5	return myPhoton.momentum.Et();
183	kiesel	1.4	}
184
185			float pt = 0;
186			float minPtDifferenz = 1E20; // should be very high
187			for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
188	kiesel	1.9	it != jetEnd; ++it ) {
189	kiesel	1.5	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
190	kiesel	1.4	if ( ptDiff < minPtDifferenz ) {
191			minPtDifferenz = ptDiff;
192			pt = it->momentum.Et();
193			}
194			}
195
196			// testing
197	kiesel	1.6	if( nearJets.size() > 1 && loggingVerbosity > 0 )
198	kiesel	1.4	std::cout << "There are several jets matching to this photon. "
199	kiesel	1.6	<< "Please check if jet-matching is correct." << std::endl;
200	kiesel	1.4	return pt;
201			}
202
203
204	kiesel	1.1	void TreeWriter::Loop() {
205	kiesel	1.4
206	kiesel	1.1	// here the event loop is implemented and the tree is filled
207			if (inputTree == 0) return;
208
209			// get number of events to be proceeded
210			Long64_t nentries = inputTree->GetEntries();
211	kiesel	1.14	// store them in histo
212			eventNumbers->Fill( "Number of generated Events", nentries );
213	kiesel	1.1	if(processNEvents <= 0 \|\| processNEvents > nentries) processNEvents = nentries;
214
215			if( loggingVerbosity > 0 )
216			std::cout << "Processing " << processNEvents << " ouf of "
217			<< nentries << " events. " << std::endl;
218
219			tree->Branch("photon", &photon);
220			tree->Branch("jet", &jet);
221	kiesel	1.12	tree->Branch("electron", &electron);
222	kiesel	1.9	tree->Branch("muon", &muon);
223	kiesel	1.1	tree->Branch("met", &met, "met/F");
224	kiesel	1.11	tree->Branch("metPhi", &met_phi, "metPhi/F");
225			tree->Branch("type1met", &type1met, "type1met/F");
226			tree->Branch("type1metPhi", &type1met_phi, "type1metPhi/F");
227	kiesel	1.5	tree->Branch("ht", &ht, "ht/F");
228	kiesel	1.1	tree->Branch("nVertex", &nVertex, "nVertex/I");
229	kiesel	1.9	tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
230	kiesel	1.1
231
232	kiesel	1.9	for (long jentry=0; jentry < processNEvents; ++jentry) {
233	kiesel	1.15	if ( loggingVerbosity>2 \|\| jentry%reportEvery==0 )
234	kiesel	1.1	std::cout << jentry << " / " << processNEvents << std :: endl;
235
236			photon.clear();
237			jet.clear();
238	kiesel	1.9	electron.clear();
239			muon.clear();
240	kiesel	1.5	ht = 0;
241	kiesel	1.1
242	kiesel	1.9	// weights
243			if (pileupHisto == 0) {
244			pu_weight = 1.;
245			} else {
246			float trueNumInteractions = -1;
247			for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
248			iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
249			if (iBX->BX == 0)
250			trueNumInteractions = iBX->trueNumInteractions;
251			}
252			pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
253			}
254
255	kiesel	1.15	// get ak5 jets
256			susy::PFJetCollection jetVector;
257			if(event->pfJets.count("ak5") == 0)
258			cout << "ERROR: Jet collection 'ak5' not found" << endl;
259			else
260			jetVector = event->pfJets.find("ak5")->second;
261	kiesel	1.9
262	kiesel	1.1	// photons
263	kiesel	1.5	if( loggingVerbosity > 1 )
264	kiesel	1.6	std::cout << "Process photons" << std::endl;
265	kiesel	1.15	susy::PhotonCollection photonVector;
266			if(event->photons.count("photons") == 0)
267			cout << "ERROR: Photon collection 'photons' not found" << endl;
268			else
269			photonVector = event->photons.find("photons")->second;
270
271			for(susy::PhotonCollection :: iterator it = photonVector.begin();
272			it != photonVector.end(); ++it ) {
273	kiesel	1.12	if( !(it->isEB() \|\| it->isEE()) && skim )
274	kiesel	1.1	continue;
275	kiesel	1.15	tree::Photon thisphoton;
276			thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
277	kiesel	1.12
278	kiesel	1.15	thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
279			thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
280			thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
281	kiesel	1.12
282	kiesel	1.15	bool loose_photon_barrel = thisphoton.pt>20
283	kiesel	1.12	&& it->isEB()
284			&& it->passelectronveto
285			&& it->hadTowOverEm<0.05
286			&& it->sigmaIetaIeta<0.012
287	kiesel	1.15	&& thisphoton.chargedIso<2.6
288			&& thisphoton.neutralIso<3.5+0.04*thisphoton.pt
289			&& thisphoton.photonIso<1.3+0.005*thisphoton.pt;
290			bool loose_photon_endcap = thisphoton.pt > 20
291	kiesel	1.12	&& it->isEE()
292			&& it->passelectronveto
293			&& it->hadTowOverEm<0.05
294			&& it->sigmaIetaIeta<0.034
295	kiesel	1.15	&& thisphoton.chargedIso<2.3
296			&& thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
297			if(!(loose_photon_endcap \|\| loose_photon_barrel \|\| thisphoton.pt > 75 ) && skim )
298	kiesel	1.1	continue;
299	kiesel	1.15	thisphoton.eta = it->momentum.Eta();
300			thisphoton.phi = it->momentum.Phi();
301			thisphoton.r9 = it->r9;
302			thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
303			thisphoton.hadTowOverEm = it->hadTowOverEm;
304			thisphoton.pixelseed = it->nPixelSeeds;
305			thisphoton.conversionSafeVeto = it->passelectronveto;
306			photon.push_back( thisphoton );
307	kiesel	1.5	if( loggingVerbosity > 2 )
308	kiesel	1.15	std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
309	kiesel	1.1	}
310	kiesel	1.5
311	kiesel	1.4	if( photon.size() == 0 && skim )
312	kiesel	1.1	continue;
313			std::sort( photon.begin(), photon.end(), tree::EtGreater);
314	kiesel	1.5	if( loggingVerbosity > 1 )
315			std::cout << "Found " << photon.size() << " photons" << std::endl;
316	kiesel	1.1
317			// jets
318			std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
319			if(pfJets_it == event->pfJets.end()){
320			if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
321			} else {
322
323			susy::PFJetCollection& jetColl = pfJets_it->second;
324
325			for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
326	kiesel	1.9	it != jetColl.end(); ++it) {
327	kiesel	1.15	tree::Jet thisjet;
328			/*std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
329	kiesel	1.1	if (s_it == it->jecScaleFactors.end()) {
330			std::cout << "JEC is not available for this jet!!!" << std::endl;
331			continue;
332			}
333			float scale = s_it->second;
334	kiesel	1.15	*/
335			float scale = 1.;
336			if(it->jecScaleFactors.count("L2L3") == 0)
337			std::cout << "ERROR: JEC is not available for this jet" << std::endl;
338			else
339			scale = it->jecScaleFactors.find("L2L3")->second;
340	kiesel	1.1	TLorentzVector corrP4 = scale * it->momentum;
341
342	kiesel	1.15	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
343			if(corrP4.Pt() < 30 && skim ) continue;
344			thisjet.pt = corrP4.Pt();
345			thisjet.eta = corrP4.Eta();
346			thisjet.phi = corrP4.Phi();
347			thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
348	kiesel	1.12	// jet composition
349	kiesel	1.15	thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
350			thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
351			thisjet.photonEnergy = it->photonEnergy;
352			thisjet.electronEnergy = it->electronEnergy;
353			thisjet.muonEnergy = it->muonEnergy;
354			thisjet.HFHadronEnergy = it->HFHadronEnergy;
355			thisjet.HFEMEnergy = it->HFEMEnergy;
356			thisjet.chargedEmEnergy = it->chargedEmEnergy;
357			thisjet.chargedMuEnergy = it->chargedMuEnergy;
358			thisjet.neutralEmEnergy = it->neutralEmEnergy;
359	kiesel	1.12
360	kiesel	1.9	if( loggingVerbosity > 2 )
361	kiesel	1.15	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
362	kiesel	1.9
363	kiesel	1.15	jet.push_back( thisjet );
364			ht += thisjet.pt;
365	kiesel	1.1	}// for jet
366			}// if, else
367	kiesel	1.4	if( jet.size() < 2 && skim )
368			continue;
369			std::sort( jet.begin(), jet.end(), tree::EtGreater);
370	kiesel	1.9	if( loggingVerbosity > 1 )
371			std::cout << "Found " << jet.size() << " jets" << std::endl;
372
373	kiesel	1.1
374			// met
375			std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
376			susy::MET* metobj = &(met_it->second);
377			met = metobj->met();
378	kiesel	1.9	met_phi = metobj->mEt.Phi();
379			if( loggingVerbosity > 2 )
380			std::cout << " met = " << met << std::endl;
381
382			std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
383			susy::MET* type1metobj = &(type1met_it->second);
384			type1met = type1metobj->met();
385			type1met_phi = type1metobj->mEt.Phi();
386			if( loggingVerbosity > 2 )
387			std::cout << " type1met = " << type1met << std::endl;
388	kiesel	1.1
389	kiesel	1.4	// electrons
390	kiesel	1.15	std::vector<susy::Electron>* eVector;
391			if( event->electrons.count("gsfElectrons") == 0) {
392			cout << "EROR: no electron collection found" << endl;
393			continue;
394			} else
395			eVector = &event->electrons.find("gsfElectrons")->second;
396
397			for(vector<susy::Electron>::iterator it = eVector->begin(); it < eVector->end(); ++it) {
398			tree::Particle thiselectron;
399			if( loggingVerbosity > 2 )
400			cout << " electron pt = " << it->momentum.Pt() << endl;
401			// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
402			// use veto electrons
403			if( it->momentum.Pt() < 20 \|\| it->momentum.Pt() > 1e6 )
404			continue; // spike rejection
405			float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
406			- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
407			) / it->momentum.Pt();
408			cout << iso << endl;
409			float d0 = d0correction( it, event );
410			cout << d0 << endl;
411			float dZ = std::abs( dZcorrection( it, event ) );
412			cout << dZ << endl;
413			if ( it->isEB() ){
414			cout << "electron in barrel" << endl;
415			if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
416			\|\| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
417			\|\| it->sigmaIetaIeta > 0.01
418			\|\| it->hcalOverEcalBc > 0.15
419			\|\| d0 > 0.04
420			\|\| dZ > 0.2
421			\|\| iso > 0.15 ){
422			cout << " no real electron" << endl;
423	kiesel	1.12	continue;
424	kiesel	1.15	}
425			}
426			else if( it->isEE() ) {
427			cout << "electron in endcap" << endl;
428			if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
429			\|\| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
430			\|\| it->sigmaIetaIeta > 0.03
431			\|\| d0 > 0.04
432			\|\| dZ > 0.2
433			\|\| iso > 0.15 ){
434			cout << "no real electron" << endl;
435	kiesel	1.9	continue;
436	kiesel	1.15	}
437			}
438			else{ // not in barrel nor in endcap
439			cout << " electron in gap" << endl;
440			continue;
441			// TODO: conversion rejection information not implemented yet, see twiki for more details
442	kiesel	1.9	}
443	kiesel	1.15	thiselectron.pt = it->momentum.Pt();
444			if( loggingVerbosity > 2 )
445			std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
446			thiselectron.eta = it->momentum.Eta();
447			thiselectron.phi = it->momentum.Phi();
448			cout << " adde electron" << endl;
449			electron.push_back( thiselectron );
450	kiesel	1.9	}
451			if( loggingVerbosity > 1 )
452			std::cout << "Found " << electron.size() << " electrons" << std::endl;
453
454			// muons
455			std::vector<susy::Muon> mVector = event->muons["muons"];
456	kiesel	1.15	tree::Particle thismuon;
457	kiesel	1.9	for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
458	kiesel	1.10	if( !( it->isPFMuon() && ( it->isGlobalMuon() \|\| it->isTrackerMuon() ) ) )
459			continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
460	kiesel	1.15	thismuon.pt = it->momentum.Et();
461			if( thismuon.pt < 20 )
462	kiesel	1.9	continue;
463	kiesel	1.15	thismuon.eta = it->momentum.Eta();
464			thismuon.phi = it->momentum.Phi();
465			muon.push_back( thismuon );
466	kiesel	1.9	}
467			if( loggingVerbosity > 1 )
468			std::cout << "Found " << muon.size() << " muons" << std::endl;
469
470	kiesel	1.4
471	kiesel	1.1	// vertices
472			nVertex = event->vertices.size();
473
474	kiesel	1.5	if( ht < 450 && skim)
475			continue;
476	kiesel	1.4
477	kiesel	1.15
478	kiesel	1.1	tree->Fill();
479			} // for jentry
480
481
482	kiesel	1.9	outFile->cd();
483	kiesel	1.14	eventNumbers->Write();
484	kiesel	1.1	tree->Write();
485			outFile->Write();
486			outFile->Close();
487			}
488	kiesel	1.3