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

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

// 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 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 ) {
                //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();
        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("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() && skim )
                                continue;
                        thisphoton->pt = getPtFromMatchedJet( *it, *event );
                        if( thisphoton->pt < 80 && skim )
                                continue;
                        thisphoton->eta = it->momentum.Eta();
                        thisphoton->phi = it->momentum.Phi();
                        thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
                        thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
                        thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
                        thisphoton->r9 = it->r9;
                        thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
                        thisphoton->hadTowOverEm = it->hadTowOverEm;
                        thisphoton->pixelseed = it->nPixelSeeds;
                        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::abs(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];
                                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 {
                        for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
                                thiselectron->pt = it->momentum.Et();
                                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;
                */

                // this seems not to work yet, where is the bug?
                /*
                std::vector<susy::Electron> eVector = event->electrons["gsfElectronsx"];
                for( std::vector<susy::Electron>::iterator it = eVector.begin(); it != eVector.end(); ++it) {
                        thiselectron->pt = it->momentum.Et();
                        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 );
                }
                */

                // 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();
        tree->Write();
        outFile->Write();
        outFile->Close();
}

Revision:	1.10
Committed:	Mon Apr 15 09:02:50 2013 UTC (12 years ago) by kiesel
Content type:	text/plain
Branch:	MAIN
Changes since 1.9:	+15 -6 lines
Log Message:	added second constructon function to handle input TChains
#	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.9	// open the output file
34	kiesel	1.7	if (loggingVerbosity_>0)
35			std::cout << "Open file " << outputName << " for writing." << std::endl;
36	kiesel	1.1	outFile = new TFile( outputName, "recreate" );
37			tree = new TTree("susyTree","Tree for single photon analysis");
38
39			// set default parameter
40			processNEvents = -1;
41			reportEvery = 1000;
42	kiesel	1.6	loggingVerbosity = loggingVerbosity_;
43	kiesel	1.5	skim = true;
44	kiesel	1.9	pileupHisto = 0;
45			}
46	kiesel	1.1
47	kiesel	1.9	void TreeWriter::PileUpWeightFile( string pileupFileName ) {
48			TFile *puFile = new TFile( pileupFileName.c_str() );
49			pileupHisto = (TH1F*) puFile->Get("pileup");
50	kiesel	1.1	}
51
52			TreeWriter::~TreeWriter() {
53	kiesel	1.9	if (pileupHisto != 0 )
54			delete pileupHisto;
55			inputTree->GetCurrentFile()->Close();
56	kiesel	1.1	}
57
58	kiesel	1.8	// useful functions
59	kiesel	1.4	float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
60			return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
61			}
62
63	kiesel	1.8	// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
64			float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
65			float eta = fabs(gamma.caloPosition.Eta());
66			float ea;
67
68			if(eta < 1.0) ea = 0.012;
69			else if(eta < 1.479) ea = 0.010;
70			else if(eta < 2.0) ea = 0.014;
71			else if(eta < 2.2) ea = 0.012;
72			else if(eta < 2.3) ea = 0.016;
73			else if(eta < 2.4) ea = 0.020;
74			else ea = 0.012;
75
76			float iso = gamma.chargedHadronIso;
77			iso = max(iso - rho*ea, (float)0.);
78
79			return iso;
80			}
81
82			float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
83			float eta = fabs(gamma.caloPosition.Eta());
84			float ea;
85
86			if(eta < 1.0) ea = 0.030;
87			else if(eta < 1.479) ea = 0.057;
88			else if(eta < 2.0) ea = 0.039;
89			else if(eta < 2.2) ea = 0.015;
90			else if(eta < 2.3) ea = 0.024;
91			else if(eta < 2.4) ea = 0.039;
92			else ea = 0.072;
93
94			float iso = gamma.neutralHadronIso;
95			iso = max(iso - rho*ea, (float)0.);
96
97			return iso;
98			}
99
100			float photonIso_corrected(susy::Photon gamma, float rho) {
101			float eta = fabs(gamma.caloPosition.Eta());
102			float ea;
103
104			if(eta < 1.0) ea = 0.148;
105			else if(eta < 1.479) ea = 0.130;
106			else if(eta < 2.0) ea = 0.112;
107			else if(eta < 2.2) ea = 0.216;
108			else if(eta < 2.3) ea = 0.262;
109			else if(eta < 2.4) ea = 0.260;
110			else ea = 0.266;
111
112			float iso = gamma.photonIso;
113			iso = max(iso - rho*ea, (float)0.);
114
115			return iso;
116			}
117
118	kiesel	1.5	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
119	kiesel	1.4	/**
120			* \brief Takes jet p_T as photon p_T
121			*
122			* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
123			* If several jets are found, take the one with the minimal pt difference
124	kiesel	1.5	* compared to the photon. If no such jets are found, keep the photon_pt
125	kiesel	1.4	* TODO: remove photon matched jet from jet-selection?
126			*/
127			std::vector<susy::PFJet> nearJets;
128			nearJets.clear();
129
130	kiesel	1.5	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
131			if(pfJets_it == myEvent.pfJets.end()){
132			if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
133	kiesel	1.4	} else {
134			susy::PFJetCollection& jetColl = pfJets_it->second;
135			for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
136	kiesel	1.9	it != jetColl.end(); ++it) {
137	kiesel	1.4	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
138			if (s_it == it->jecScaleFactors.end()) {
139			std::cout << "JEC is not available for this jet!!!" << std::endl;
140			continue;
141			}
142			float scale = s_it->second;
143			TLorentzVector corrP4 = scale * it->momentum;
144	kiesel	1.5	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
145	kiesel	1.4	if (deltaR_ > 0.3) continue;
146	kiesel	1.5	if( loggingVerbosity > 0 )
147			std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
148	kiesel	1.4	nearJets.push_back( *it );
149			}// for jet
150			}// if, else
151
152			if ( nearJets.size() == 0 ) {
153	kiesel	1.5	//std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
154			return myPhoton.momentum.Et();
155	kiesel	1.4	}
156
157			float pt = 0;
158			float minPtDifferenz = 1E20; // should be very high
159			for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
160	kiesel	1.9	it != jetEnd; ++it ) {
161	kiesel	1.5	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
162	kiesel	1.4	if ( ptDiff < minPtDifferenz ) {
163			minPtDifferenz = ptDiff;
164			pt = it->momentum.Et();
165			}
166			}
167
168			// testing
169	kiesel	1.6	if( nearJets.size() > 1 && loggingVerbosity > 0 )
170	kiesel	1.4	std::cout << "There are several jets matching to this photon. "
171	kiesel	1.6	<< "Please check if jet-matching is correct." << std::endl;
172	kiesel	1.4	return pt;
173			}
174
175
176	kiesel	1.1	void TreeWriter::Loop() {
177	kiesel	1.4
178	kiesel	1.1	// here the event loop is implemented and the tree is filled
179			if (inputTree == 0) return;
180
181			// get number of events to be proceeded
182			Long64_t nentries = inputTree->GetEntries();
183			if(processNEvents <= 0 \|\| processNEvents > nentries) processNEvents = nentries;
184
185			if( loggingVerbosity > 0 )
186			std::cout << "Processing " << processNEvents << " ouf of "
187			<< nentries << " events. " << std::endl;
188
189			tree::Photon *thisphoton = new tree::Photon();
190			tree::Jet *thisjet = new tree::Jet();
191
192			tree->Branch("photon", &photon);
193			tree->Branch("jet", &jet);
194	kiesel	1.10	//tree->Branch("electron", &electron);
195	kiesel	1.9	tree->Branch("muon", &muon);
196	kiesel	1.1	tree->Branch("met", &met, "met/F");
197	kiesel	1.5	tree->Branch("ht", &ht, "ht/F");
198	kiesel	1.1	tree->Branch("nVertex", &nVertex, "nVertex/I");
199	kiesel	1.9	tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
200	kiesel	1.1
201
202	kiesel	1.9	for (long jentry=0; jentry < processNEvents; ++jentry) {
203	kiesel	1.4	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
204	kiesel	1.1	std::cout << jentry << " / " << processNEvents << std :: endl;
205			inputTree->GetEntry(jentry);
206
207			photon.clear();
208			jet.clear();
209	kiesel	1.9	electron.clear();
210			muon.clear();
211	kiesel	1.5	ht = 0;
212	kiesel	1.1
213	kiesel	1.9	// weights
214			if (pileupHisto == 0) {
215			pu_weight = 1.;
216			} else {
217			float trueNumInteractions = -1;
218			for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
219			iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
220			if (iBX->BX == 0)
221			trueNumInteractions = iBX->trueNumInteractions;
222			}
223			pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
224			}
225
226
227	kiesel	1.1	// photons
228	kiesel	1.5	if( loggingVerbosity > 1 )
229	kiesel	1.6	std::cout << "Process photons" << std::endl;
230	kiesel	1.1	std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
231			for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
232	kiesel	1.9	it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
233	kiesel	1.4	if( ! it->isEB() && skim )
234	kiesel	1.1	continue;
235	kiesel	1.4	thisphoton->pt = getPtFromMatchedJet( it, event );
236			if( thisphoton->pt < 80 && skim )
237	kiesel	1.1	continue;
238			thisphoton->eta = it->momentum.Eta();
239	kiesel	1.4	thisphoton->phi = it->momentum.Phi();
240	kiesel	1.8	thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
241			thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
242			thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
243	kiesel	1.1	thisphoton->r9 = it->r9;
244			thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
245			thisphoton->hadTowOverEm = it->hadTowOverEm;
246			thisphoton->pixelseed = it->nPixelSeeds;
247			photon.push_back( *thisphoton );
248	kiesel	1.5	if( loggingVerbosity > 2 )
249	kiesel	1.6	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
250	kiesel	1.1	}
251	kiesel	1.5
252	kiesel	1.4	if( photon.size() == 0 && skim )
253	kiesel	1.1	continue;
254			std::sort( photon.begin(), photon.end(), tree::EtGreater);
255	kiesel	1.5	if( loggingVerbosity > 1 )
256			std::cout << "Found " << photon.size() << " photons" << std::endl;
257	kiesel	1.1
258			// jets
259			std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
260			if(pfJets_it == event->pfJets.end()){
261			if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
262			} else {
263
264			susy::PFJetCollection& jetColl = pfJets_it->second;
265
266			for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
267	kiesel	1.9	it != jetColl.end(); ++it) {
268	kiesel	1.1	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
269			if (s_it == it->jecScaleFactors.end()) {
270			std::cout << "JEC is not available for this jet!!!" << std::endl;
271			continue;
272			}
273			float scale = s_it->second;
274			TLorentzVector corrP4 = scale * it->momentum;
275
276	kiesel	1.4	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
277	kiesel	1.8	if(corrP4.Et() < 30 && skim ) continue;
278	kiesel	1.1	thisjet->pt = corrP4.Et();
279			thisjet->eta = corrP4.Eta();
280	kiesel	1.4	thisjet->phi = corrP4.Phi();
281	kiesel	1.9	thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
282			if( loggingVerbosity > 2 )
283			std::cout << " p_T, jet = " << thisjet->pt << std::endl;
284
285	kiesel	1.1	jet.push_back( *thisjet );
286	kiesel	1.5	ht += thisjet->pt;
287	kiesel	1.1	}// for jet
288			}// if, else
289	kiesel	1.4	if( jet.size() < 2 && skim )
290			continue;
291			std::sort( jet.begin(), jet.end(), tree::EtGreater);
292	kiesel	1.9	if( loggingVerbosity > 1 )
293			std::cout << "Found " << jet.size() << " jets" << std::endl;
294
295	kiesel	1.1
296			// met
297			std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
298			susy::MET* metobj = &(met_it->second);
299			met = metobj->met();
300	kiesel	1.9	met_phi = metobj->mEt.Phi();
301			if( loggingVerbosity > 2 )
302			std::cout << " met = " << met << std::endl;
303
304			std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
305			susy::MET* type1metobj = &(type1met_it->second);
306			type1met = type1metobj->met();
307			type1met_phi = type1metobj->mEt.Phi();
308			if( loggingVerbosity > 2 )
309			std::cout << " type1met = " << type1met << std::endl;
310	kiesel	1.1
311	kiesel	1.4	// electrons
312	kiesel	1.10	/*
313	kiesel	1.9	tree::Particle* thiselectron = new tree::Particle();
314			map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
315			if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
316			cout << "gsfElectrons not found!" << endl;
317			} else {
318			for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
319			thiselectron->pt = it->momentum.Et();
320			if( thiselectron->pt < 20 )
321			continue;
322			if( loggingVerbosity > 2 )
323			std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
324			thiselectron->eta = it->momentum.Eta();
325			thiselectron->phi = it->momentum.Phi();
326			electron.push_back( *thiselectron );
327			}
328			}
329			if( loggingVerbosity > 1 )
330			std::cout << "Found " << electron.size() << " electrons" << std::endl;
331	kiesel	1.10	*/
332	kiesel	1.9
333			// this seems not to work yet, where is the bug?
334			/*
335			std::vector<susy::Electron> eVector = event->electrons["gsfElectronsx"];
336			for( std::vector<susy::Electron>::iterator it = eVector.begin(); it != eVector.end(); ++it) {
337			thiselectron->pt = it->momentum.Et();
338			if( thiselectron->pt < 20 )
339			continue;
340			if( loggingVerbosity > 2 )
341			std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
342			thiselectron->eta = it->momentum.Eta();
343			thiselectron->phi = it->momentum.Phi();
344			electron.push_back( *thiselectron );
345			}
346			*/
347
348			// muons
349			std::vector<susy::Muon> mVector = event->muons["muons"];
350			tree::Particle* thismuon = new tree::Particle();
351			for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
352	kiesel	1.10	if( !( it->isPFMuon() && ( it->isGlobalMuon() \|\| it->isTrackerMuon() ) ) )
353			continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
354	kiesel	1.9	thismuon->pt = it->momentum.Et();
355			if( thismuon->pt < 20 )
356			continue;
357			thismuon->eta = it->momentum.Eta();
358			thismuon->phi = it->momentum.Phi();
359			muon.push_back( *thismuon );
360			}
361			if( loggingVerbosity > 1 )
362			std::cout << "Found " << muon.size() << " muons" << std::endl;
363
364	kiesel	1.4
365	kiesel	1.1	// vertices
366			nVertex = event->vertices.size();
367
368	kiesel	1.5	if( ht < 450 && skim)
369			continue;
370	kiesel	1.4
371	kiesel	1.1	tree->Fill();
372			} // for jentry
373
374
375	kiesel	1.9	outFile->cd();
376	kiesel	1.1	tree->Write();
377			outFile->Write();
378			outFile->Close();
379			}
380	kiesel	1.3