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

float deltaPhi( float phi1, float phi2) {
        float result = phi1 - phi2;
        while (result > M_PI) result -= 2*M_PI;
        while (result <= -M_PI) result += 2*M_PI;
        return result;
}

// useful functions
float deltaR( const TLorentzVector& v1, const TLorentzVector& v2 ) {
        // deltaR  = sqrt ( deltaEta^2 + deltaPhi^2 )
        return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(deltaPhi(v1.Phi(),v2.Phi()), 2) );
}

float effectiveAreaElectron( float eta ) {
        // needed by calculating the isolation for electrons
        // see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
        // only for Delta R = 0.3 on 2012 Data
        eta = fabs( eta );
        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("weight", &weight, "weight/D");
        tree->Branch("genElectron", &genElectron);
        tree->Branch("genPhoton", &genPhoton);

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

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

                // weights
                if (pileupHisto == 0) {
                        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;
                        }
                        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->isEE() || it->isEB()) && it->momentum.Pt()<20 && it->isEBEtaGap() && it->isEBPhiGap() && it->isEERingGap() && it->isEEDeeGap() && it->isEBEEGap() && skim )
                                continue;
                        tree::Photon thisphoton;

                        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 = it->isEB()
                                && 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 = it->isEE()
                                && it->hadTowOverEm<0.05
                                && it->sigmaIetaIeta<0.034
                                && thisphoton.chargedIso<2.3
                                && thisphoton.neutralIso<2.9+0.04*thisphoton.pt;

                        thisphoton.ptJet = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
                        if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.ptJet > 75 ) && skim )
                                continue;
                        thisphoton.pt = it->momentum.Pt();
                        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;
                        thisphoton.genInformation = 0;
                        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;

                if( ht < 450 && skim)
                        continue;


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

                tree::Particle thisGenParticle;
                for( std::vector<susy::Particle>::iterator it = event->genParticles.begin(); it != event->genParticles.end(); ++it ) {
                        if( it->momentum.Pt() < 20 ) continue;
                        thisGenParticle.pt = it->momentum.Pt();
                        thisGenParticle.eta = it->momentum.Eta();
                        thisGenParticle.phi = it->momentum.Phi();
                        switch( std::abs(it->pdgId) ) {
                                case 22: // photon
                                        genPhoton.push_back( thisGenParticle );
                                        break;
                                case 11: // electron
                                        // Demand a W boson as mother particle of electron
                                        if( abs(event->genParticles[it->motherIndex].pdgId) == 24 )
                                                genElectron.push_back( thisGenParticle );
                                        break;
                        }
                }

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


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

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