kiesel/TreeWriter/treeWriter.cc

#include<iostream>
#include<math.h>
#include<string>

#include "TSystem.h"

#include "treeWriter.h"

using namespace std;

TreeWriter::TreeWriter( TString inputName, TString 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 );
        Init( outputName, loggingVerbosity_ );
}

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


void TreeWriter::Init( TString 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, "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();
}

// useful functions
float TreeWriter::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( susy::Electron electron, 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());
        cout << "return " << d0 << endl;
        return d0;
}

float dZcorrection( susy::Electron electron, 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 TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
        /**
         * \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();

        std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
        if(pfJets_it == myEvent.pfJets.end()){
                if(myEvent.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) {
                        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 > 0 )
                                std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
                        nearJets.push_back( *it );
                }// for jet
        }// if, else

        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::Photon *thisphoton = new tree::Photon();
        tree::Jet *thisjet = new tree::Jet();

        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>0 && jentry%reportEvery==0 )
                        std::cout << jentry << " / " << processNEvents << std :: endl;
                inputTree->GetEntry(jentry);

                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 ) );
                }


                // photons
                if( loggingVerbosity > 1 )
                        std::cout << "Process photons" << std::endl;
                std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
                for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
                                it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
                        if( !(it->isEB() || it->isEE()) && skim )
                                continue;
                        thisphoton->pt = getPtFromMatchedJet( *it, *event );

                        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) {
                                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;

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