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.13 by kiesel, Tue Apr 16 18:53:55 2013 UTC vs.
Revision 1.15 by kiesel, Thu Apr 18 13:38:55 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  
20 < void TreeWriter::Init( TString outputName, int loggingVerbosity_ ) {
20 > void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) {
21  
22          if (loggingVerbosity_ > 0)
23                  std::cout << "Set Branch Address of susy::Event" << std::endl;
24          event = new susy::Event;
25          inputTree->SetBranchAddress("susyEvent", &event);
26  
27 +        // Here the number of proceeded events will be stored. For plotting, simply use L*sigma/eventNumber
28 +        eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated events", 1, 0, 1);
29 +        eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated events");
30 +
31          // open the output file
32          if (loggingVerbosity_>0)
33                  std::cout << "Open file " << outputName << " for writing." << std::endl;
34 <        outFile = new TFile( outputName, "recreate" );
34 >        outFile = new TFile( outputName.c_str(), "recreate" );
35          tree = new TTree("susyTree","Tree for single photon analysis");
36  
37          // set default parameter
# Line 53 | Line 51 | TreeWriter::~TreeWriter() {
51          if (pileupHisto != 0 )
52                  delete pileupHisto;
53          inputTree->GetCurrentFile()->Close();
54 +        delete inputTree;
55 +        delete event;
56 +        delete outFile;
57 +        delete tree;
58   }
59  
60   // useful functions
61 < float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
61 > float deltaR( TLorentzVector v1, TLorentzVector v2 ) {
62          return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
63   }
64  
# Line 130 | Line 132 | float photonIso_corrected(susy::Photon g
132          return iso;
133   }
134  
135 < float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
135 > float d0correction( const susy::Electron& electron, const susy::Event& event ) {
136 >        TVector3 beamspot = event.vertices[0].position;
137 >        susy::Track track = event.tracks[electron.gsfTrackIndex];
138 >        float d0 = track.d0() - beamspot.X()*sin(track.phi()) + beamspot.Y()*cos(track.phi());
139 >        return d0;
140 > }
141 >
142 > float dZcorrection( const susy::Electron& electron, const susy::Event& event ) {
143 >        TVector3 beamspot = event.vertices[0].position;
144 >        susy::Track track = event.tracks[electron.gsfTrackIndex];
145 >
146 >        if(track.momentum.Pt() == 0.) return 1.e6;
147 >        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());
148 >        return dz;
149 > }
150 >
151 > float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) {
152          /**
153           * \brief Takes jet p_T as photon p_T
154           *
# Line 142 | Line 160 | float TreeWriter::getPtFromMatchedJet( s
160          std::vector<susy::PFJet> nearJets;
161          nearJets.clear();
162  
163 <        std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
164 <        if(pfJets_it == myEvent.pfJets.end()){
165 <                if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
166 <        } else {
167 <                susy::PFJetCollection& jetColl = pfJets_it->second;
168 <                for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
169 <                                it != jetColl.end(); ++it) {
170 <                        std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
171 <                        if (s_it == it->jecScaleFactors.end()) {
172 <                                std::cout << "JEC is not available for this jet!!!" << std::endl;
173 <                                continue;
174 <                        }
175 <                        float scale = s_it->second;
176 <                        TLorentzVector corrP4 = scale * it->momentum;
177 <                        float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
160 <                        if (deltaR_ > 0.3) continue;
161 <                        if( loggingVerbosity > 0 )
162 <                                std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
163 <                        nearJets.push_back( *it );
164 <                }// for jet
165 <        }// if, else
163 >        for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
164 >                        it != jetColl.end(); ++it) {
165 >                std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
166 >                if (s_it == it->jecScaleFactors.end()) {
167 >                        std::cout << "JEC is not available for this jet!!!" << std::endl;
168 >                        continue;
169 >                }
170 >                float scale = s_it->second;
171 >                TLorentzVector corrP4 = scale * it->momentum;
172 >                float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
173 >                if (deltaR_ > 0.3) continue;
174 >                if( loggingVerbosity > 2 )
175 >                        std::cout << " pT_jet / pT_gamma = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
176 >                nearJets.push_back( *it );
177 >        }// for jet
178  
179          if ( nearJets.size() == 0 ) {
180                  if( loggingVerbosity > 1 )
# Line 196 | Line 208 | void TreeWriter::Loop() {
208  
209          // get number of events to be proceeded
210          Long64_t nentries = inputTree->GetEntries();
211 +        // store them in histo
212 +        eventNumbers->Fill( "Number of generated Events", nentries );
213          if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;
214  
215          if( loggingVerbosity > 0 )
216                  std::cout << "Processing " << processNEvents << " ouf of "
217                          << nentries << " events. " << std::endl;
218  
205        tree::Photon *thisphoton = new tree::Photon();
206        tree::Jet *thisjet = new tree::Jet();
207
219          tree->Branch("photon", &photon);
220          tree->Branch("jet", &jet);
221          tree->Branch("electron", &electron);
# Line 219 | Line 230 | void TreeWriter::Loop() {
230  
231  
232          for (long jentry=0; jentry < processNEvents; ++jentry) {
233 <                if ( loggingVerbosity>0 && jentry%reportEvery==0 )
233 >                if ( loggingVerbosity>2 || jentry%reportEvery==0 )
234                          std::cout << jentry << " / " << processNEvents << std :: endl;
224                inputTree->GetEntry(jentry);
235  
236                  photon.clear();
237                  jet.clear();
# Line 242 | Line 252 | void TreeWriter::Loop() {
252                          pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
253                  }
254  
255 +                // get ak5 jets
256 +                susy::PFJetCollection jetVector;
257 +                if(event->pfJets.count("ak5") == 0)
258 +                        cout << "ERROR: Jet collection 'ak5' not found" << endl;
259 +                else
260 +                        jetVector = event->pfJets.find("ak5")->second;
261  
262                  // photons
263                  if( loggingVerbosity > 1 )
264                          std::cout << "Process photons" << std::endl;
265 <                std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
266 <                for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
267 <                                it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
265 >                susy::PhotonCollection photonVector;
266 >                if(event->photons.count("photons") == 0)
267 >                        cout << "ERROR: Photon collection 'photons' not found" << endl;
268 >                else
269 >                        photonVector = event->photons.find("photons")->second;
270 >
271 >                for(susy::PhotonCollection :: iterator it = photonVector.begin();
272 >                                it != photonVector.end(); ++it ) {
273                          if( !(it->isEB() || it->isEE()) && skim )
274                                  continue;
275 <                        thisphoton->pt = getPtFromMatchedJet( *it, *event );
275 >                        tree::Photon thisphoton;
276 >                        thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
277  
278 <                        thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
279 <                        thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
280 <                        thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
278 >                        thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
279 >                        thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
280 >                        thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
281  
282 <                        bool loose_photon_barrel = thisphoton->pt>20
282 >                        bool loose_photon_barrel = thisphoton.pt>20
283                                  && it->isEB()
284                                  && it->passelectronveto
285                                  && it->hadTowOverEm<0.05
286                                  && it->sigmaIetaIeta<0.012
287 <                                && thisphoton->chargedIso<2.6
288 <                                && thisphoton->neutralIso<3.5+0.04*thisphoton->pt
289 <                                && thisphoton->photonIso<1.3+0.005*thisphoton->pt;
290 <                        bool loose_photon_endcap = thisphoton->pt > 20
287 >                                && thisphoton.chargedIso<2.6
288 >                                && thisphoton.neutralIso<3.5+0.04*thisphoton.pt
289 >                                && thisphoton.photonIso<1.3+0.005*thisphoton.pt;
290 >                        bool loose_photon_endcap = thisphoton.pt > 20
291                                  && it->isEE()
292                                  && it->passelectronveto
293                                  && it->hadTowOverEm<0.05
294                                  && it->sigmaIetaIeta<0.034
295 <                                && thisphoton->chargedIso<2.3
296 <                                && thisphoton->neutralIso<2.9+0.04*thisphoton->pt;
297 <                        if(!(loose_photon_endcap || loose_photon_barrel || thisphoton->pt > 75 ) && skim )
295 >                                && thisphoton.chargedIso<2.3
296 >                                && thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
297 >                        if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
298                                  continue;
299 <                        thisphoton->eta = it->momentum.Eta();
300 <                        thisphoton->phi = it->momentum.Phi();
301 <                        thisphoton->r9 = it->r9;
302 <                        thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
303 <                        thisphoton->hadTowOverEm = it->hadTowOverEm;
304 <                        thisphoton->pixelseed = it->nPixelSeeds;
305 <                        thisphoton->conversionSafeVeto = it->passelectronveto;
306 <                        photon.push_back( *thisphoton );
299 >                        thisphoton.eta = it->momentum.Eta();
300 >                        thisphoton.phi = it->momentum.Phi();
301 >                        thisphoton.r9 = it->r9;
302 >                        thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
303 >                        thisphoton.hadTowOverEm = it->hadTowOverEm;
304 >                        thisphoton.pixelseed = it->nPixelSeeds;
305 >                        thisphoton.conversionSafeVeto = it->passelectronveto;
306 >                        photon.push_back( thisphoton );
307                          if( loggingVerbosity > 2 )
308 <                                std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
308 >                                std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
309                  }
310  
311                  if( photon.size() == 0 && skim )
# Line 302 | Line 324 | void TreeWriter::Loop() {
324  
325                          for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
326                                          it != jetColl.end(); ++it) {
327 <                                std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
327 >                                tree::Jet thisjet;
328 >                                /*std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
329                                  if (s_it == it->jecScaleFactors.end()) {
330                                          std::cout << "JEC is not available for this jet!!!" << std::endl;
331                                          continue;
332                                  }
333                                  float scale = s_it->second;
334 +                                */
335 +                                float scale = 1.;
336 +                                if(it->jecScaleFactors.count("L2L3") == 0)
337 +                                        std::cout << "ERROR: JEC is not available for this jet" << std::endl;
338 +                                else
339 +                                        scale = it->jecScaleFactors.find("L2L3")->second;
340                                  TLorentzVector corrP4 = scale * it->momentum;
341  
342                                  if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
343 <                                if(corrP4.Et() < 30 && skim ) continue;
344 <                                thisjet->pt = corrP4.Et();
345 <                                thisjet->eta = corrP4.Eta();
346 <                                thisjet->phi = corrP4.Phi();
347 <                                thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
343 >                                if(corrP4.Pt() < 30 && skim ) continue;
344 >                                thisjet.pt = corrP4.Pt();
345 >                                thisjet.eta = corrP4.Eta();
346 >                                thisjet.phi = corrP4.Phi();
347 >                                thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
348                                  // jet composition
349 <                                thisjet->chargedHadronEnergy = it->chargedHadronEnergy;
350 <                                thisjet->neutralHadronEnergy = it->neutralHadronEnergy;
351 <                                thisjet->photonEnergy = it->photonEnergy;
352 <                                thisjet->electronEnergy = it->electronEnergy;
353 <                                thisjet->muonEnergy = it->muonEnergy;
354 <                                thisjet->HFHadronEnergy = it->HFHadronEnergy;
355 <                                thisjet->HFEMEnergy = it->HFEMEnergy;
356 <                                thisjet->chargedEmEnergy = it->chargedEmEnergy;
357 <                                thisjet->chargedMuEnergy = it->chargedMuEnergy;
358 <                                thisjet->neutralEmEnergy = it->neutralEmEnergy;
349 >                                thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
350 >                                thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
351 >                                thisjet.photonEnergy = it->photonEnergy;
352 >                                thisjet.electronEnergy = it->electronEnergy;
353 >                                thisjet.muonEnergy = it->muonEnergy;
354 >                                thisjet.HFHadronEnergy = it->HFHadronEnergy;
355 >                                thisjet.HFEMEnergy = it->HFEMEnergy;
356 >                                thisjet.chargedEmEnergy = it->chargedEmEnergy;
357 >                                thisjet.chargedMuEnergy = it->chargedMuEnergy;
358 >                                thisjet.neutralEmEnergy = it->neutralEmEnergy;
359  
360                                  if( loggingVerbosity > 2 )
361 <                                        std::cout << " p_T, jet = " << thisjet->pt << std::endl;
361 >                                        std::cout << " p_T, jet = " << thisjet.pt << std::endl;
362  
363 <                                jet.push_back( *thisjet );
364 <                                ht += thisjet->pt;
363 >                                jet.push_back( thisjet );
364 >                                ht += thisjet.pt;
365                          }// for jet
366                  }// if, else
367                  if( jet.size() < 2 && skim )
# Line 358 | Line 387 | void TreeWriter::Loop() {
387                          std::cout << " type1met = " << type1met << std::endl;
388  
389                  // electrons
390 <                tree::Particle* thiselectron = new tree::Particle();
391 <                map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
392 <                if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
393 <                        cout << "gsfElectrons not found!" << endl;
394 <                } else {
395 <                        for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
367 <                                // for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
368 <                                if( loggingVerbosity > 0)
369 <                                        std::cout << "TODO: piAtVtx is not calculated correctly" << std::endl;
370 <                                //float pAtVtx = it->trackMomentums.find("AtVtxWithConstraint")->second.P();
371 <                                float pAtVtx = it->ecalEnergy;
372 <                                if( it->momentum.Pt() < 1  || it->momentum.Pt() > 1e6 || pAtVtx == 0 )
373 <                                        continue;
374 <                                float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
375 <                                                                                                                        - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
376 <                                                        ) / it->momentum.Pt();
377 <                                if ( it->isEE() ){
378 <                                        if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
379 <                                                        || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.15
380 <                                                        || it->sigmaIetaIeta > 0.01
381 <                                                        || it->hcalOverEcalBc > 0.12
382 <                                                        || it->vertex.Perp() > 0.02
383 <                                                        || it->vertex.Z() > 0.2
384 <                                                        || fabs(1./it->ecalEnergy - 1./pAtVtx ) > 0.05
385 <                                                        || iso > 0.15 )
386 <                                                continue;
387 <                                        }
388 <                                else if( it->isEB() ) {
389 <                                        if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.009
390 <                                                        || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.10
391 <                                                        || it->sigmaIetaIeta > 0.03
392 <                                                        || it->hcalOverEcalBc > 0.10
393 <                                                        || it->vertex.Perp() > 0.02
394 <                                                        || it->vertex.Z() > 0.2
395 <                                                        || fabs(1./it->ecalEnergy - 1./pAtVtx) > 0.05
396 <                                                        || iso > 0.15 )
397 <                                                continue;
398 <                                        }
399 <                                else // not in barrel nor in endcap
400 <                                        continue;
401 <                                // TODO: conversion rejection information not implemented yet, see twiki for more details
390 >                std::vector<susy::Electron>* eVector;
391 >                if( event->electrons.count("gsfElectrons") == 0) {
392 >                        cout << "EROR: no electron collection found" << endl;
393 >                        continue;
394 >                } else
395 >                        eVector = &event->electrons.find("gsfElectrons")->second;
396  
397 <                                thiselectron->pt = it->momentum.Pt();
398 <                                if( thiselectron->pt < 20 )
397 >                for(vector<susy::Electron>::iterator it = eVector->begin(); it < eVector->end(); ++it) {
398 >                        tree::Particle thiselectron;
399 >                        if( loggingVerbosity > 2 )
400 >                                cout << " electron pt = " << it->momentum.Pt() << endl;
401 >                        // for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
402 >                        // use veto electrons
403 >                        if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
404 >                                continue; // spike rejection
405 >                        float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
406 >                                                                                                                - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
407 >                                                ) / it->momentum.Pt();
408 >                        cout << iso << endl;
409 >                        float d0 = d0correction( *it, *event );
410 >                        cout << d0 << endl;
411 >                        float dZ = std::abs( dZcorrection( *it, *event ) );
412 >                        cout << dZ << endl;
413 >                        if ( it->isEB() ){
414 >                                cout << "electron in barrel" << endl;
415 >                                if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
416 >                                                || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
417 >                                                || it->sigmaIetaIeta > 0.01
418 >                                                || it->hcalOverEcalBc > 0.15
419 >                                                || d0 > 0.04
420 >                                                || dZ > 0.2
421 >                                                || iso > 0.15 ){
422 >                                        cout << " no real electron" << endl;
423                                          continue;
424 <                                if( loggingVerbosity > 2 )
425 <                                        std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
426 <                                thiselectron->eta = it->momentum.Eta();
427 <                                thiselectron->phi = it->momentum.Phi();
428 <                                electron.push_back( *thiselectron );
424 >                                }
425 >                                }
426 >                        else if( it->isEE() ) {
427 >                                cout << "electron in endcap" << endl;
428 >                                if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
429 >                                                || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
430 >                                                || it->sigmaIetaIeta > 0.03
431 >                                                || d0 > 0.04
432 >                                                || dZ > 0.2
433 >                                                || iso > 0.15 ){
434 >                                        cout << "no real electron" << endl;
435 >                                        continue;
436 >                                }
437 >                                }
438 >                        else{ // not in barrel nor in endcap
439 >                                cout << " electron in gap" << endl;
440 >                                continue;
441 >                        // TODO: conversion rejection information not implemented yet, see twiki for more details
442                          }
443 +                        thiselectron.pt = it->momentum.Pt();
444 +                        if( loggingVerbosity > 2 )
445 +                                std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
446 +                        thiselectron.eta = it->momentum.Eta();
447 +                        thiselectron.phi = it->momentum.Phi();
448 +                        cout << " adde electron" << endl;
449 +                        electron.push_back( thiselectron );
450                  }
451                  if( loggingVerbosity > 1 )
452                          std::cout << "Found " << electron.size() << " electrons" << std::endl;
453  
454                  // muons
455                  std::vector<susy::Muon> mVector = event->muons["muons"];
456 <                tree::Particle* thismuon = new tree::Particle();
456 >                tree::Particle thismuon;
457                  for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
458                          if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
459                                  continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
460 <                        thismuon->pt = it->momentum.Et();
461 <                        if( thismuon->pt < 20 )
460 >                        thismuon.pt = it->momentum.Et();
461 >                        if( thismuon.pt < 20 )
462                                  continue;
463 <                        thismuon->eta = it->momentum.Eta();
464 <                        thismuon->phi = it->momentum.Phi();
465 <                        muon.push_back( *thismuon );
463 >                        thismuon.eta = it->momentum.Eta();
464 >                        thismuon.phi = it->momentum.Phi();
465 >                        muon.push_back( thismuon );
466                  }
467                  if( loggingVerbosity > 1 )
468                          std::cout << "Found " << muon.size() << " muons" << std::endl;
# Line 436 | Line 474 | void TreeWriter::Loop() {
474                  if( ht < 450 && skim)
475                          continue;
476  
477 +
478                  tree->Fill();
479          } // for jentry
480  
481  
482          outFile->cd();
483 +        eventNumbers->Write();
484          tree->Write();
485          outFile->Write();
486          outFile->Close();

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines