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.15 by kiesel, Thu Apr 18 13:38:55 2013 UTC vs.
Revision 1.25 by kiesel, Wed Apr 24 15:25:35 2013 UTC

# Line 16 | Line 16 | TreeWriter::TreeWriter( TChain* inputTre
16          Init( outputName, loggingVerbosity_ );
17   }
18  
19
19   void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) {
20  
21          if (loggingVerbosity_ > 0)
# Line 57 | Line 56 | TreeWriter::~TreeWriter() {
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 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 78 | 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 96 | 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 114 | 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 133 | Line 148 | float photonIso_corrected(susy::Photon g
148   }
149  
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());
# Line 140 | Line 156 | float d0correction( const susy::Electron
156   }
157  
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 148 | Line 165 | float dZcorrection( const susy::Electron
165          return dz;
166   }
167  
168 < float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) {
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 160 | Line 177 | float getPtFromMatchedJet( susy::Photon
177          std::vector<susy::PFJet> nearJets;
178          nearJets.clear();
179  
180 <        for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
180 >        for(std::vector<susy::PFJet>::const_iterator it = jetColl.begin();
181                          it != jetColl.end(); ++it) {
182 <                std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
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                  }
170                float scale = s_it->second;
190                  TLorentzVector corrP4 = scale * it->momentum;
191                  float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
192                  if (deltaR_ > 0.3) continue;
# Line 202 | Line 221 | float getPtFromMatchedJet( susy::Photon
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 209 | 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 )
# Line 227 | Line 253 | void TreeWriter::Loop() {
253          tree->Branch("ht", &ht, "ht/F");
254          tree->Branch("nVertex", &nVertex, "nVertex/I");
255          tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
256 +        tree->Branch("genElectron", &genElectron);
257 +        tree->Branch("genPhoton", &genPhoton);
258  
259 <
260 <        for (long jentry=0; jentry < processNEvents; ++jentry) {
233 <                if ( loggingVerbosity>2 || jentry%reportEvery==0 )
259 >        for (unsigned long jentry=0; jentry < processNEvents; ++jentry) {
260 >                if ( loggingVerbosity>1 || jentry%reportEvery==0 )
261                          std::cout << jentry << " / " << processNEvents << std :: endl;
262 +                inputTree->LoadTree( jentry );
263 +                inputTree->GetEntry( jentry );
264  
265                  photon.clear();
266                  jet.clear();
267                  electron.clear();
268                  muon.clear();
269 +                genElectron.clear();
270 +                genPhoton.clear();
271                  ht = 0;
272  
273                  // weights
# Line 253 | Line 284 | void TreeWriter::Loop() {
284                  }
285  
286                  // get ak5 jets
287 <                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;
287 >                std::vector<susy::PFJet> jetVector = event->pfJets["ak5"];
288  
289                  // photons
290 <                if( loggingVerbosity > 1 )
264 <                        std::cout << "Process photons" << std::endl;
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;
290 >                std::vector<susy::Photon> photonVector = event->photons["photons"];
291  
292 <                for(susy::PhotonCollection :: iterator it = photonVector.begin();
292 >                for(std::vector<susy::Photon>::iterator it = photonVector.begin();
293                                  it != photonVector.end(); ++it ) {
294 <                        if( !(it->isEB() || it->isEE()) && skim )
294 >                        if( !(it->isEE() || it->isEB()) && it->isEBEtaGap() && it->isEBPhiGap() && it->isEERingGap() && it->isEEDeeGap() && it->isEBEEGap() && skim )
295                                  continue;
296                          tree::Photon thisphoton;
297                          thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
# Line 294 | Line 315 | void TreeWriter::Loop() {
315                                  && it->sigmaIetaIeta<0.034
316                                  && thisphoton.chargedIso<2.3
317                                  && thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
318 +
319                          if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
320                                  continue;
321                          thisphoton.eta = it->momentum.Eta();
# Line 315 | Line 337 | void TreeWriter::Loop() {
337                          std::cout << "Found " << photon.size() << " photons" << std::endl;
338  
339                  // jets
318                std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
319                if(pfJets_it == event->pfJets.end()){
320                        if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
321                } else {
340  
323                        susy::PFJetCollection& jetColl = pfJets_it->second;
341  
342 <                        for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
343 <                                        it != jetColl.end(); ++it) {
344 <                                tree::Jet thisjet;
345 <                                /*std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
346 <                                if (s_it == it->jecScaleFactors.end()) {
347 <                                        std::cout << "JEC is not available for this jet!!!" << std::endl;
348 <                                        continue;
349 <                                }
350 <                                float scale = s_it->second;
351 <                                */
352 <                                float scale = 1.;
353 <                                if(it->jecScaleFactors.count("L2L3") == 0)
354 <                                        std::cout << "ERROR: JEC is not available for this jet" << std::endl;
355 <                                else
356 <                                        scale = it->jecScaleFactors.find("L2L3")->second;
357 <                                TLorentzVector corrP4 = scale * it->momentum;
341 <
342 <                                if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
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;
342 >                for(std::vector<susy::PFJet>::iterator it = jetVector.begin();
343 >                                it != jetVector.end(); ++it) {
344 >                        tree::Jet thisjet;
345 >
346 >                        // scale with JEC
347 >                        float scale = 1.;
348 >                        if(it->jecScaleFactors.count("L2L3") == 0)
349 >                                std::cout << "ERROR: JEC is not available for this jet" << std::endl;
350 >                        else
351 >                                scale = it->jecScaleFactors.find("L2L3")->second;
352 >                        TLorentzVector corrP4 = scale * it->momentum;
353 >
354 >                        // Calculate HT.
355 >                        // The definiton differs from the saved jet, since trigger is described better
356 >                        if( std::abs( corrP4.Eta() ) < 3 && corrP4.Pt() > 40 )
357 >                                ht += thisjet.pt;
358  
359 <                                if( loggingVerbosity > 2 )
360 <                                        std::cout << " p_T, jet = " << thisjet.pt << std::endl;
359 >                        if(std::abs(corrP4.Eta()) > 2.6 && skim ) continue;
360 >                        if(corrP4.Pt() < 30 && skim ) continue;
361 >                        thisjet.pt = corrP4.Pt();
362 >                        thisjet.eta = corrP4.Eta();
363 >                        thisjet.phi = corrP4.Phi();
364 >                        thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
365 >                        // jet composition
366 >                        thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
367 >                        thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
368 >                        thisjet.photonEnergy = it->photonEnergy;
369 >                        thisjet.electronEnergy = it->electronEnergy;
370 >                        thisjet.muonEnergy = it->muonEnergy;
371 >                        thisjet.HFHadronEnergy = it->HFHadronEnergy;
372 >                        thisjet.HFEMEnergy = it->HFEMEnergy;
373 >                        thisjet.chargedEmEnergy = it->chargedEmEnergy;
374 >                        thisjet.chargedMuEnergy = it->chargedMuEnergy;
375 >                        thisjet.neutralEmEnergy = it->neutralEmEnergy;
376  
377 <                                jet.push_back( thisjet );
378 <                                ht += thisjet.pt;
379 <                        }// for jet
380 <                }// if, else
377 >                        if( loggingVerbosity > 2 )
378 >                                std::cout << " p_T, jet = " << thisjet.pt << std::endl;
379 >
380 >                        jet.push_back( thisjet );
381 >                }// for jet
382                  if( jet.size() < 2 && skim )
383                          continue;
384                  std::sort( jet.begin(), jet.end(), tree::EtGreater);
385                  if( loggingVerbosity > 1 )
386                          std::cout << "Found " << jet.size() << " jets" << std::endl;
387  
388 +                if( ht < 450 && skim)
389 +                        continue;
390 +
391 +
392  
393                  // met
394                  std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
# Line 387 | Line 406 | void TreeWriter::Loop() {
406                          std::cout << " type1met = " << type1met << std::endl;
407  
408                  // electrons
409 <                std::vector<susy::Electron>* eVector;
410 <                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 <                for(vector<susy::Electron>::iterator it = eVector->begin(); it < eVector->end(); ++it) {
409 >                std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
410 >                for(std::vector<susy::Electron>::iterator it = eVector.begin(); it < eVector.end(); ++it) {
411                          tree::Particle thiselectron;
412                          if( loggingVerbosity > 2 )
413                                  cout << " electron pt = " << it->momentum.Pt() << endl;
# Line 402 | Line 415 | void TreeWriter::Loop() {
415                          // use veto electrons
416                          if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
417                                  continue; // spike rejection
418 <                        float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
406 <                                                                                                                - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
418 >                        float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (float)0. )
419                                                  ) / it->momentum.Pt();
408                        cout << iso << endl;
420                          float d0 = d0correction( *it, *event );
410                        cout << d0 << endl;
421                          float dZ = std::abs( dZcorrection( *it, *event ) );
412                        cout << dZ << endl;
422                          if ( it->isEB() ){
414                                cout << "electron in barrel" << endl;
423                                  if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
424                                                  || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
425                                                  || it->sigmaIetaIeta > 0.01
426                                                  || it->hcalOverEcalBc > 0.15
427                                                  || d0 > 0.04
428                                                  || dZ > 0.2
429 <                                                || iso > 0.15 ){
422 <                                        cout << " no real electron" << endl;
429 >                                                || iso > 0.15 )
430                                          continue;
431                                  }
425                                }
432                          else if( it->isEE() ) {
427                                cout << "electron in endcap" << endl;
433                                  if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
434                                                  || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
435                                                  || it->sigmaIetaIeta > 0.03
436                                                  || d0 > 0.04
437                                                  || dZ > 0.2
438 <                                                || iso > 0.15 ){
434 <                                        cout << "no real electron" << endl;
438 >                                                || iso > 0.15 )
439                                          continue;
440                                  }
441 <                                }
438 <                        else{ // not in barrel nor in endcap
439 <                                cout << " electron in gap" << endl;
441 >                        else // not in barrel nor in endcap
442                                  continue;
443 <                        // TODO: conversion rejection information not implemented yet, see twiki for more details
442 <                        }
443 >
444                          thiselectron.pt = it->momentum.Pt();
445                          if( loggingVerbosity > 2 )
446                                  std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
447                          thiselectron.eta = it->momentum.Eta();
448                          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"];
455                  tree::Particle thismuon;
456 +                std::vector<susy::Muon> mVector = event->muons["muons"];
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
# Line 467 | Line 467 | void TreeWriter::Loop() {
467                  if( loggingVerbosity > 1 )
468                          std::cout << "Found " << muon.size() << " muons" << std::endl;
469  
470
470                  // vertices
471                  nVertex = event->vertices.size();
472  
473 <                if( ht < 450 && skim)
474 <                        continue;
475 <
473 >                tree::Particle thisGenParticle;
474 >                for( std::vector<susy::Particle>::iterator it = event->genParticles.begin(); it != event->genParticles.end(); ++it ) {
475 >                        if( it->momentum.Pt() < 20 ) continue;
476 >                        thisGenParticle.pt = it->momentum.Pt();
477 >                        thisGenParticle.eta = it->momentum.Eta();
478 >                        thisGenParticle.phi = it->momentum.Phi();
479 >                        switch( std::abs(it->pdgId) ) {
480 >                                case 22: // photon
481 >                                        genPhoton.push_back( thisGenParticle );
482 >                                        break;
483 >                                case 11: // electron
484 >                                        // Demand a W boson as mother particle of electron
485 >                                        if( abs(event->genParticles[it->motherIndex].pdgId) == 24 )
486 >                                                genElectron.push_back( thisGenParticle );
487 >                                        break;
488 >                        }
489 >                }
490  
491                  tree->Fill();
492          } // for jentry

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines