ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/kiesel/TreeWriter/treeWriter.cc
(Generate patch)

Comparing UserCode/kiesel/TreeWriter/treeWriter.cc (file contents):
Revision 1.14 by kiesel, Wed Apr 17 09:57:52 2013 UTC vs.
Revision 1.17 by kiesel, Fri Apr 19 10:13:02 2013 UTC

# Line 1 | Line 1
1 #include<iostream>
2 #include<math.h>
3 #include<string>
4
5 #include "TSystem.h"
6
1   #include "treeWriter.h"
2  
3   using namespace std;
4  
5 < TreeWriter::TreeWriter( TString inputName, TString outputName, int loggingVerbosity_ ) {
5 > TreeWriter::TreeWriter( std::string inputName, std::string outputName, int loggingVerbosity_ ) {
6          // read the input file
7          inputTree = new TChain("susyTree");
8          if (loggingVerbosity_ > 0)
9                  std::cout << "Add files to chain" << std::endl;
10 <        inputTree->Add( inputName );
10 >        inputTree->Add( inputName.c_str() );
11          Init( outputName, loggingVerbosity_ );
12   }
13  
14 < TreeWriter::TreeWriter( TChain* inputTree_, TString outputName, int loggingVerbosity_ ) {
14 > TreeWriter::TreeWriter( TChain* inputTree_, std::string outputName, int loggingVerbosity_ ) {
15          inputTree = inputTree_;
16          Init( outputName, loggingVerbosity_ );
17   }
18  
19 <
26 < void TreeWriter::Init( TString outputName, int loggingVerbosity_ ) {
19 > void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) {
20  
21          if (loggingVerbosity_ > 0)
22                  std::cout << "Set Branch Address of susy::Event" << std::endl;
# Line 31 | Line 24 | void TreeWriter::Init( TString outputNam
24          inputTree->SetBranchAddress("susyEvent", &event);
25  
26          // Here the number of proceeded events will be stored. For plotting, simply use L*sigma/eventNumber
27 <        eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated Events", 1, 0, 1);
28 <        eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated Events");
27 >        eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated events", 1, 0, 1);
28 >        eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated events");
29  
30          // open the output file
31          if (loggingVerbosity_>0)
32                  std::cout << "Open file " << outputName << " for writing." << std::endl;
33 <        outFile = new TFile( outputName, "recreate" );
33 >        outFile = new TFile( outputName.c_str(), "recreate" );
34          tree = new TTree("susyTree","Tree for single photon analysis");
35  
36          // set default parameter
# Line 57 | Line 50 | TreeWriter::~TreeWriter() {
50          if (pileupHisto != 0 )
51                  delete pileupHisto;
52          inputTree->GetCurrentFile()->Close();
53 +        delete inputTree;
54 +        delete event;
55 +        delete outFile;
56 +        delete tree;
57 + }
58 +
59 + 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   // useful functions
74 < float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
74 > float deltaR( const TLorentzVector& v1, const TLorentzVector& v2 ) {
75 >        // deltaR  = sqrt ( deltaEta^2 + deltaPhi^2 )
76          return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
77   }
78  
79   float effectiveAreaElectron( float eta ) {
80 +        // needed by calculating the isolation for electrons
81          // see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
82          // only for Delta R = 0.3 on 2012 Data
83          eta = fabs( eta );
# Line 80 | Line 93 | float effectiveAreaElectron( float eta )
93   }
94  
95   // correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
96 < float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
96 > float chargedHadronIso_corrected(const susy::Photon& gamma, float rho) {
97          float eta = fabs(gamma.caloPosition.Eta());
98          float ea;
99  
# Line 98 | Line 111 | float chargedHadronIso_corrected(susy::P
111          return iso;
112   }
113  
114 < float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
114 > float neutralHadronIso_corrected(const susy::Photon& gamma, float rho) {
115          float eta = fabs(gamma.caloPosition.Eta());
116          float ea;
117  
# Line 116 | Line 129 | float neutralHadronIso_corrected(susy::P
129          return iso;
130   }
131  
132 < float photonIso_corrected(susy::Photon gamma, float rho) {
132 > float photonIso_corrected(const susy::Photon& gamma, float rho) {
133          float eta = fabs(gamma.caloPosition.Eta());
134          float ea;
135  
# Line 134 | Line 147 | float photonIso_corrected(susy::Photon g
147          return iso;
148   }
149  
150 < float d0correction( susy::Electron electron, susy::Event event ) {
150 > float d0correction( const susy::Electron& electron, const susy::Event& event ) {
151 >        // copied from Brian Francis
152          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());
141        cout << "return " << d0 << endl;
155          return d0;
156   }
157  
158 < float dZcorrection( susy::Electron electron, susy::Event event ) {
158 > float dZcorrection( const susy::Electron& electron, const susy::Event& event ) {
159 >        // copied from Brian Francis
160          TVector3 beamspot = event.vertices[0].position;
161          susy::Track track = event.tracks[electron.gsfTrackIndex];
162  
# Line 151 | Line 165 | float dZcorrection( susy::Electron elect
165          return dz;
166   }
167  
168 < float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
168 > float getPtFromMatchedJet( const susy::Photon& myPhoton, const susy::PFJetCollection& jetColl, int loggingVerbosity = 0 ) {
169          /**
170           * \brief Takes jet p_T as photon p_T
171           *
# Line 163 | Line 177 | float TreeWriter::getPtFromMatchedJet( s
177          std::vector<susy::PFJet> nearJets;
178          nearJets.clear();
179  
180 <        std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
181 <        if(pfJets_it == myEvent.pfJets.end()){
182 <                if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
183 <        } else {
184 <                susy::PFJetCollection& jetColl = pfJets_it->second;
185 <                for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
186 <                                it != jetColl.end(); ++it) {
187 <                        std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
188 <                        if (s_it == it->jecScaleFactors.end()) {
189 <                                std::cout << "JEC is not available for this jet!!!" << std::endl;
190 <                                continue;
191 <                        }
192 <                        float scale = s_it->second;
193 <                        TLorentzVector corrP4 = scale * it->momentum;
194 <                        float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
195 <                        if (deltaR_ > 0.3) continue;
196 <                        if( loggingVerbosity > 0 )
183 <                                std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
184 <                        nearJets.push_back( *it );
185 <                }// for jet
186 <        }// if, else
180 >        for(std::vector<susy::PFJet>::const_iterator it = jetColl.begin();
181 >                        it != jetColl.end(); ++it) {
182 >                float scale = 1.;
183 >                std::map<TString,Float_t>::const_iterator s_it = it->jecScaleFactors.find("L2L3");
184 >                if (s_it == it->jecScaleFactors.end()) {
185 >                        std::cout << "JEC is not available for this jet!!!" << std::endl;
186 >                        continue;
187 >                } else {
188 >                        scale = s_it->second;
189 >                }
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  
198          if ( nearJets.size() == 0 ) {
199                  if( loggingVerbosity > 1 )
# Line 211 | Line 221 | float TreeWriter::getPtFromMatchedJet( s
221  
222  
223   void TreeWriter::Loop() {
224 +        /**
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  
232          // here the event loop is implemented and the tree is filled
233          if (inputTree == 0) return;
# Line 218 | Line 235 | void TreeWriter::Loop() {
235          // get number of events to be proceeded
236          Long64_t nentries = inputTree->GetEntries();
237          // store them in histo
238 <        eventNumbers->Fill( "Number of generated Events", nentries );
238 >        eventNumbers->Fill( "Number of generated events", nentries );
239          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  
228        tree::Photon *thisphoton = new tree::Photon();
229        tree::Jet *thisjet = new tree::Jet();
230
245          tree->Branch("photon", &photon);
246          tree->Branch("jet", &jet);
247          tree->Branch("electron", &electron);
# Line 241 | Line 255 | void TreeWriter::Loop() {
255          tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
256  
257  
258 <        for (long jentry=0; jentry < processNEvents; ++jentry) {
259 <                if ( loggingVerbosity>0 && jentry%reportEvery==0 )
258 >        for (unsigned long jentry=0; jentry < processNEvents; ++jentry) {
259 >                if ( loggingVerbosity>1 || jentry%reportEvery==0 )
260                          std::cout << jentry << " / " << processNEvents << std :: endl;
261 <                inputTree->GetEntry(jentry);
261 >                inputTree->LoadTree( jentry );
262 >                inputTree->GetEntry( jentry );
263  
264                  photon.clear();
265                  jet.clear();
# Line 265 | Line 280 | void TreeWriter::Loop() {
280                          pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
281                  }
282  
283 +                // get ak5 jets
284 +                std::vector<susy::PFJet> jetVector = event->pfJets["ak5"];
285  
286                  // photons
287 <                if( loggingVerbosity > 1 )
288 <                        std::cout << "Process photons" << std::endl;
289 <                std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
290 <                for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
274 <                                it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
287 >                std::vector<susy::Photon> photonVector = event->photons["photons"];
288 >
289 >                for(std::vector<susy::Photon>::iterator it = photonVector.begin();
290 >                                it != photonVector.end(); ++it ) {
291                          if( !(it->isEB() || it->isEE()) && skim )
292                                  continue;
293 <                        thisphoton->pt = getPtFromMatchedJet( *it, *event );
293 >                        tree::Photon thisphoton;
294 >                        thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
295  
296 <                        thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
297 <                        thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
298 <                        thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
296 >                        thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
297 >                        thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
298 >                        thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
299  
300 <                        bool loose_photon_barrel = thisphoton->pt>20
300 >                        bool loose_photon_barrel = thisphoton.pt>20
301                                  && it->isEB()
302                                  && it->passelectronveto
303                                  && it->hadTowOverEm<0.05
304                                  && it->sigmaIetaIeta<0.012
305 <                                && thisphoton->chargedIso<2.6
306 <                                && thisphoton->neutralIso<3.5+0.04*thisphoton->pt
307 <                                && thisphoton->photonIso<1.3+0.005*thisphoton->pt;
308 <                        bool loose_photon_endcap = thisphoton->pt > 20
305 >                                && thisphoton.chargedIso<2.6
306 >                                && thisphoton.neutralIso<3.5+0.04*thisphoton.pt
307 >                                && thisphoton.photonIso<1.3+0.005*thisphoton.pt;
308 >                        bool loose_photon_endcap = thisphoton.pt > 20
309                                  && it->isEE()
310                                  && it->passelectronveto
311                                  && it->hadTowOverEm<0.05
312                                  && it->sigmaIetaIeta<0.034
313 <                                && thisphoton->chargedIso<2.3
314 <                                && thisphoton->neutralIso<2.9+0.04*thisphoton->pt;
315 <                        if(!(loose_photon_endcap || loose_photon_barrel || thisphoton->pt > 75 ) && skim )
313 >                                && thisphoton.chargedIso<2.3
314 >                                && thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
315 >                        if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
316                                  continue;
317 <                        thisphoton->eta = it->momentum.Eta();
318 <                        thisphoton->phi = it->momentum.Phi();
319 <                        thisphoton->r9 = it->r9;
320 <                        thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
321 <                        thisphoton->hadTowOverEm = it->hadTowOverEm;
322 <                        thisphoton->pixelseed = it->nPixelSeeds;
323 <                        thisphoton->conversionSafeVeto = it->passelectronveto;
324 <                        photon.push_back( *thisphoton );
317 >                        thisphoton.eta = it->momentum.Eta();
318 >                        thisphoton.phi = it->momentum.Phi();
319 >                        thisphoton.r9 = it->r9;
320 >                        thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
321 >                        thisphoton.hadTowOverEm = it->hadTowOverEm;
322 >                        thisphoton.pixelseed = it->nPixelSeeds;
323 >                        thisphoton.conversionSafeVeto = it->passelectronveto;
324 >                        photon.push_back( thisphoton );
325                          if( loggingVerbosity > 2 )
326 <                                std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
326 >                                std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
327                  }
328  
329                  if( photon.size() == 0 && skim )
# Line 316 | Line 333 | void TreeWriter::Loop() {
333                          std::cout << "Found " << photon.size() << " photons" << std::endl;
334  
335                  // jets
319                std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
320                if(pfJets_it == event->pfJets.end()){
321                        if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
322                } else {
336  
324                        susy::PFJetCollection& jetColl = pfJets_it->second;
337  
338 <                        for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
339 <                                        it != jetColl.end(); ++it) {
340 <                                std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
341 <                                if (s_it == it->jecScaleFactors.end()) {
342 <                                        std::cout << "JEC is not available for this jet!!!" << std::endl;
343 <                                        continue;
344 <                                }
345 <                                float scale = s_it->second;
346 <                                TLorentzVector corrP4 = scale * it->momentum;
338 >                for(std::vector<susy::PFJet>::iterator it = jetVector.begin();
339 >                                it != jetVector.end(); ++it) {
340 >                        tree::Jet thisjet;
341 >
342 >                        // scale with JEC
343 >                        float scale = 1.;
344 >                        if(it->jecScaleFactors.count("L2L3") == 0)
345 >                                std::cout << "ERROR: JEC is not available for this jet" << std::endl;
346 >                        else
347 >                                scale = it->jecScaleFactors.find("L2L3")->second;
348 >                        TLorentzVector corrP4 = scale * it->momentum;
349  
350 <                                if(std::fabs(corrP4.Eta()) > 3.0 && skim ) continue;
351 <                                if(corrP4.Et() < 30 && skim ) continue;
352 <                                thisjet->pt = corrP4.Et();
353 <                                thisjet->eta = corrP4.Eta();
354 <                                thisjet->phi = corrP4.Phi();
355 <                                thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
356 <                                // jet composition
357 <                                thisjet->chargedHadronEnergy = it->chargedHadronEnergy;
358 <                                thisjet->neutralHadronEnergy = it->neutralHadronEnergy;
359 <                                thisjet->photonEnergy = it->photonEnergy;
360 <                                thisjet->electronEnergy = it->electronEnergy;
361 <                                thisjet->muonEnergy = it->muonEnergy;
362 <                                thisjet->HFHadronEnergy = it->HFHadronEnergy;
363 <                                thisjet->HFEMEnergy = it->HFEMEnergy;
364 <                                thisjet->chargedEmEnergy = it->chargedEmEnergy;
365 <                                thisjet->chargedMuEnergy = it->chargedMuEnergy;
366 <                                thisjet->neutralEmEnergy = it->neutralEmEnergy;
367 <
368 <                                if( loggingVerbosity > 2 )
369 <                                        std::cout << " p_T, jet = " << thisjet->pt << std::endl;
370 <
371 <                                jet.push_back( *thisjet );
372 <                                ht += thisjet->pt;
373 <                        }// for jet
360 <                }// if, else
350 >                        if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
351 >                        if(corrP4.Pt() < 30 && skim ) continue;
352 >                        thisjet.pt = corrP4.Pt();
353 >                        thisjet.eta = corrP4.Eta();
354 >                        thisjet.phi = corrP4.Phi();
355 >                        thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
356 >                        // jet composition
357 >                        thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
358 >                        thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
359 >                        thisjet.photonEnergy = it->photonEnergy;
360 >                        thisjet.electronEnergy = it->electronEnergy;
361 >                        thisjet.muonEnergy = it->muonEnergy;
362 >                        thisjet.HFHadronEnergy = it->HFHadronEnergy;
363 >                        thisjet.HFEMEnergy = it->HFEMEnergy;
364 >                        thisjet.chargedEmEnergy = it->chargedEmEnergy;
365 >                        thisjet.chargedMuEnergy = it->chargedMuEnergy;
366 >                        thisjet.neutralEmEnergy = it->neutralEmEnergy;
367 >
368 >                        if( loggingVerbosity > 2 )
369 >                                std::cout << " p_T, jet = " << thisjet.pt << std::endl;
370 >
371 >                        jet.push_back( thisjet );
372 >                        ht += thisjet.pt;
373 >                }// for jet
374                  if( jet.size() < 2 && skim )
375                          continue;
376                  std::sort( jet.begin(), jet.end(), tree::EtGreater);
# Line 381 | Line 394 | void TreeWriter::Loop() {
394                          std::cout << " type1met = " << type1met << std::endl;
395  
396                  // electrons
397 <                tree::Particle* thiselectron = new tree::Particle();
398 <                map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
399 <                if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
400 <                        cout << "gsfElectrons not found!" << endl;
401 <                } else {
402 <                        cout << "now begin electrons " << endl;
403 <                        for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
404 <                                // for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification for veto electrons
405 <                                if( it->momentum.Pt() < 1  || it->momentum.Pt() > 1e6 )
406 <                                        continue; // spike rejection
407 <                                float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
408 <                                                                                                                        - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
409 <                                                        ) / it->momentum.Pt();
410 <                                cout << iso << endl;
411 <                                cout << d0correction( *it, *event ) << endl;
412 <                                float test = d0correction( *it, *event );
413 <                                cout << test << endl;
414 <                                float dZ = std::fabs( dZcorrection( *it, *event ) );
415 <                                cout <<" find e" << endl;
416 <                                float d0 = 0.;
417 <                                if ( it->isEB() ){
418 <                                        if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
406 <                                                        || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
407 <                                                        || it->sigmaIetaIeta > 0.01
408 <                                                        || it->hcalOverEcalBc > 0.15
409 <                                                        || d0 > 0.04
410 <                                                        || dZ > 0.2
411 <                                                        || iso > 0.15 )
412 <                                                continue;
413 <                                        }
414 <                                else if( it->isEE() ) {
415 <                                        if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
416 <                                                        || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
417 <                                                        || it->sigmaIetaIeta > 0.03
418 <                                                        || d0 > 0.04
419 <                                                        || dZ > 0.2
420 <                                                        || iso > 0.15 )
421 <                                                continue;
422 <                                        }
423 <                                else // not in barrel nor in endcap
397 >                std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
398 >                for(std::vector<susy::Electron>::iterator it = eVector.begin(); it < eVector.end(); ++it) {
399 >                        tree::Particle thiselectron;
400 >                        if( loggingVerbosity > 2 )
401 >                                cout << " electron pt = " << it->momentum.Pt() << endl;
402 >                        // for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
403 >                        // use veto electrons
404 >                        if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
405 >                                continue; // spike rejection
406 >                        float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
407 >                                                                                                                - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
408 >                                                ) / it->momentum.Pt();
409 >                        float d0 = d0correction( *it, *event );
410 >                        float dZ = std::abs( dZcorrection( *it, *event ) );
411 >                        if ( it->isEB() ){
412 >                                if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
413 >                                                || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
414 >                                                || it->sigmaIetaIeta > 0.01
415 >                                                || it->hcalOverEcalBc > 0.15
416 >                                                || d0 > 0.04
417 >                                                || dZ > 0.2
418 >                                                || iso > 0.15 )
419                                          continue;
420 <                                // TODO: conversion rejection information not implemented yet, see twiki for more details
421 <
422 <                                thiselectron->pt = it->momentum.Pt();
423 <                                if( thiselectron->pt < 20 )
420 >                                }
421 >                        else if( it->isEE() ) {
422 >                                if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
423 >                                                || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
424 >                                                || it->sigmaIetaIeta > 0.03
425 >                                                || d0 > 0.04
426 >                                                || dZ > 0.2
427 >                                                || iso > 0.15 )
428                                          continue;
429 <                                if( loggingVerbosity > 2 )
430 <                                        std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
431 <                                thiselectron->eta = it->momentum.Eta();
432 <                                thiselectron->phi = it->momentum.Phi();
433 <                                electron.push_back( *thiselectron );
434 <                        }
429 >                                }
430 >                        else // not in barrel nor in endcap
431 >                                continue;
432 >
433 >                        thiselectron.pt = it->momentum.Pt();
434 >                        if( loggingVerbosity > 2 )
435 >                                std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
436 >                        thiselectron.eta = it->momentum.Eta();
437 >                        thiselectron.phi = it->momentum.Phi();
438 >                        electron.push_back( thiselectron );
439                  }
440                  if( loggingVerbosity > 1 )
441                          std::cout << "Found " << electron.size() << " electrons" << std::endl;
442  
443                  // muons
444 +                tree::Particle thismuon;
445                  std::vector<susy::Muon> mVector = event->muons["muons"];
442                tree::Particle* thismuon = new tree::Particle();
446                  for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
447                          if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
448                                  continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
449 <                        thismuon->pt = it->momentum.Et();
450 <                        if( thismuon->pt < 20 )
449 >                        thismuon.pt = it->momentum.Et();
450 >                        if( thismuon.pt < 20 )
451                                  continue;
452 <                        thismuon->eta = it->momentum.Eta();
453 <                        thismuon->phi = it->momentum.Phi();
454 <                        muon.push_back( *thismuon );
452 >                        thismuon.eta = it->momentum.Eta();
453 >                        thismuon.phi = it->momentum.Phi();
454 >                        muon.push_back( thismuon );
455                  }
456                  if( loggingVerbosity > 1 )
457                          std::cout << "Found " << muon.size() << " muons" << std::endl;
458  
456
459                  // vertices
460                  nVertex = event->vertices.size();
461  
462                  if( ht < 450 && skim)
463                          continue;
464  
465 +
466                  tree->Fill();
467          } // for jentry
468  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines