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

template<class Type>
std::vector<Type>* getVectorFromMap( map<TString, vector<Type> > myMap, TString search ){
        // if nothing is found, return a empty vector (which has to be deleted never the less)
        std::vector<Type>* vec;
        typename map<TString, vector<Type> >::iterator mapIt = myMap.find( search );
        if( mapIt == myMap.end() ) {
                cout << "ERROR: Collection \"" << search << "\" not found!" << endl;
                vec = new std::vector<Type>;
        } else
                vec = &mapIt->second;
        cout << vec->size() << endl;
        return vec;
}

// useful functions
float deltaR( const TLorentzVector& v1, const TLorentzVector& v2 ) {
        // deltaR  = sqrt ( deltaEta^2 + deltaPhi^2 )
        return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(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 );
        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(const 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(const 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(const 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 ) {
        // 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());
        return d0;
}

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

        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( const susy::Photon& myPhoton, const 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>::const_iterator it = jetColl.begin();
                        it != jetColl.end(); ++it) {
                float scale = 1.;
                std::map<TString,Float_t>::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;
                }
                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() {
        /**
         * \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;

        // 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 (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();
                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
                std::vector<susy::PFJet> jetVector = event->pfJets["ak5"];

                // photons
                std::vector<susy::Photon> photonVector = event->photons["photons"];

                for(std::vector<susy::Photon>::iterator it = photonVector.begin();
                                it != photonVector.end(); ++it ) {
                        if( !it->isEB() && 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_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


                for(std::vector<susy::PFJet>::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(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;

                        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;


                // 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 = event->electrons["gsfElectrons"];
                for(std::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)0. )
                                                ) / it->momentum.Pt();
                        float d0 = d0correction( *it, *event );
                        float dZ = std::abs( dZcorrection( *it, *event ) );
                        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;

                        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();
                        electron.push_back( thiselectron );
                }
                if( loggingVerbosity > 1 )
                        std::cout << "Found " << electron.size() << " electrons" << std::endl;

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