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.6 by kiesel, Tue Apr 9 08:24:40 2013 UTC vs.
Revision 1.9 by kiesel, Sat Apr 13 09:53:02 2013 UTC

# Line 1 | Line 1
1   #include<iostream>
2   #include<math.h>
3 + #include<string>
4  
5   #include "TSystem.h"
6  
# Line 7 | Line 8
8  
9   using namespace std;
10  
11 < TreeWriter::TreeWriter(TString inputName, TString outputName, int loggingVerbosity_ ) {
11 > TreeWriter::TreeWriter( TString inputName, TString outputName, int loggingVerbosity_ ) {
12          // read the input file
13          inputTree = new TChain("susyTree");
14          if (loggingVerbosity_ > 0)
# Line 20 | Line 21 | TreeWriter::TreeWriter(TString inputName
21          inputTree->SetBranchAddress("susyEvent", &event);
22  
23          // open the output file
24 +        if (loggingVerbosity_>0)
25 +                std::cout << "Open file " << outputName << " for writing." << std::endl;
26          outFile = new TFile( outputName, "recreate" );
27          tree = new TTree("susyTree","Tree for single photon analysis");
28  
# Line 28 | Line 31 | TreeWriter::TreeWriter(TString inputName
31          reportEvery = 1000;
32          loggingVerbosity = loggingVerbosity_;
33          skim = true;
34 +        pileupHisto = 0;
35 + }
36  
37 + void TreeWriter::PileUpWeightFile( string pileupFileName ) {
38 +        TFile *puFile = new TFile( pileupFileName.c_str() );
39 +        pileupHisto = (TH1F*) puFile->Get("pileup");
40   }
41  
42 +
43 +
44   TreeWriter::~TreeWriter() {
45 <        if (!inputTree) return;
46 <        delete inputTree->GetCurrentFile();
45 >        if (pileupHisto != 0 )
46 >                delete pileupHisto;
47 >        inputTree->GetCurrentFile()->Close();
48   }
49  
50 + // useful functions
51   float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
52          return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
53   }
54  
55 + // correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
56 + float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
57 +  float eta = fabs(gamma.caloPosition.Eta());
58 +  float ea;
59 +
60 +  if(eta < 1.0) ea = 0.012;
61 +  else if(eta < 1.479) ea = 0.010;
62 +  else if(eta < 2.0) ea = 0.014;
63 +  else if(eta < 2.2) ea = 0.012;
64 +  else if(eta < 2.3) ea = 0.016;
65 +  else if(eta < 2.4) ea = 0.020;
66 +  else ea = 0.012;
67 +
68 +  float iso = gamma.chargedHadronIso;
69 +  iso = max(iso - rho*ea, (float)0.);
70 +
71 +  return iso;
72 + }
73 +
74 + float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
75 +  float eta = fabs(gamma.caloPosition.Eta());
76 +  float ea;
77 +
78 +  if(eta < 1.0) ea = 0.030;
79 +  else if(eta < 1.479) ea = 0.057;
80 +  else if(eta < 2.0) ea = 0.039;
81 +  else if(eta < 2.2) ea = 0.015;
82 +  else if(eta < 2.3) ea = 0.024;
83 +  else if(eta < 2.4) ea = 0.039;
84 +  else ea = 0.072;
85 +
86 +  float iso = gamma.neutralHadronIso;
87 +  iso = max(iso - rho*ea, (float)0.);
88 +
89 +  return iso;
90 + }
91 +
92 + float photonIso_corrected(susy::Photon gamma, float rho) {
93 +  float eta = fabs(gamma.caloPosition.Eta());
94 +  float ea;
95 +
96 +  if(eta < 1.0) ea = 0.148;
97 +  else if(eta < 1.479) ea = 0.130;
98 +  else if(eta < 2.0) ea = 0.112;
99 +  else if(eta < 2.2) ea = 0.216;
100 +  else if(eta < 2.3) ea = 0.262;
101 +  else if(eta < 2.4) ea = 0.260;
102 +  else ea = 0.266;
103 +
104 +  float iso = gamma.photonIso;
105 +  iso = max(iso - rho*ea, (float)0.);
106 +
107 +  return iso;
108 + }
109 +
110   float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
111          /**
112           * \brief Takes jet p_T as photon p_T
# Line 58 | Line 125 | float TreeWriter::getPtFromMatchedJet( s
125          } else {
126                  susy::PFJetCollection& jetColl = pfJets_it->second;
127                  for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
128 <                                it != jetColl.end(); it++) {
128 >                                it != jetColl.end(); ++it) {
129                          std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
130                          if (s_it == it->jecScaleFactors.end()) {
131                                  std::cout << "JEC is not available for this jet!!!" << std::endl;
# Line 82 | Line 149 | float TreeWriter::getPtFromMatchedJet( s
149          float pt = 0;
150          float minPtDifferenz = 1E20; // should be very high
151          for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
152 <                        it != jetEnd; it++ ) {
152 >                        it != jetEnd; ++it ) {
153                  float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
154                  if (  ptDiff < minPtDifferenz ) {
155                          minPtDifferenz = ptDiff;
# Line 116 | Line 183 | void TreeWriter::Loop() {
183  
184          tree->Branch("photon", &photon);
185          tree->Branch("jet", &jet);
186 +        tree->Branch("electron", &electron);
187 +        tree->Branch("muon", &muon);
188          tree->Branch("met", &met, "met/F");
189          tree->Branch("ht", &ht, "ht/F");
190          tree->Branch("nVertex", &nVertex, "nVertex/I");
191 <        tree->Branch("nElectron", &nElectron, "nElectron/I");
191 >        tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
192  
193  
194 <        for (long jentry=0; jentry < processNEvents; jentry++) {
194 >        for (long jentry=0; jentry < processNEvents; ++jentry) {
195                  if ( loggingVerbosity>0 && jentry%reportEvery==0 )
196                          std::cout << jentry << " / " << processNEvents << std :: endl;
197                  inputTree->GetEntry(jentry);
198  
199                  photon.clear();
200                  jet.clear();
201 +                electron.clear();
202 +                muon.clear();
203                  ht = 0;
204  
205 +                // weights
206 +                if (pileupHisto == 0) {
207 +                        pu_weight = 1.;
208 +                } else {
209 +                        float trueNumInteractions = -1;
210 +                        for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
211 +                                        iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
212 +                                if (iBX->BX == 0)
213 +                                        trueNumInteractions = iBX->trueNumInteractions;
214 +                        }
215 +                        pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
216 +                }
217 +
218 +
219                  // photons
220                  if( loggingVerbosity > 1 )
221                          std::cout << "Process photons" << std::endl;
222                  std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
223                  for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
224 <                                it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
224 >                                it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
225                          if( ! it->isEB() && skim )
226                                  continue;
227                          thisphoton->pt = getPtFromMatchedJet( *it, *event );
# Line 144 | Line 229 | void TreeWriter::Loop() {
229                                  continue;
230                          thisphoton->eta = it->momentum.Eta();
231                          thisphoton->phi = it->momentum.Phi();
232 <                        thisphoton->chargedIso = it->chargedHadronIso;
233 <                        thisphoton->neutralIso = it->neutralHadronIso;
234 <                        thisphoton->photonIso = it->photonIso;
235 <                        if ( it->r9 > 1 && skim ) // if == 1 ?
232 >                        thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
233 >                        thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
234 >                        thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
235 >                        if ( it->r9 >= 1 && skim )
236                                  continue;
237                          thisphoton->r9 = it->r9;
238                          thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
239                          thisphoton->hadTowOverEm = it->hadTowOverEm;
240                          thisphoton->pixelseed = it->nPixelSeeds;
241                          photon.push_back( *thisphoton );
157 –                        ht += thisphoton->pt;
242                          if( loggingVerbosity > 2 )
243                                  std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
244                  }
# Line 174 | Line 258 | void TreeWriter::Loop() {
258                          susy::PFJetCollection& jetColl = pfJets_it->second;
259  
260                          for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
261 <                                        it != jetColl.end(); it++) {
261 >                                        it != jetColl.end(); ++it) {
262                                  std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
263                                  if (s_it == it->jecScaleFactors.end()) {
264                                          std::cout << "JEC is not available for this jet!!!" << std::endl;
# Line 184 | Line 268 | void TreeWriter::Loop() {
268                                  TLorentzVector corrP4 = scale * it->momentum;
269  
270                                  if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
271 +                                if(corrP4.Et() < 30 && skim ) continue;
272                                  thisjet->pt = corrP4.Et();
273                                  thisjet->eta = corrP4.Eta();
274                                  thisjet->phi = corrP4.Phi();
275 +                                thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
276 +                                if( loggingVerbosity > 2 )
277 +                                        std::cout << " p_T, jet = " << thisjet->pt << std::endl;
278 +
279                                  jet.push_back( *thisjet );
280                                  ht += thisjet->pt;
281                          }// for jet
# Line 194 | Line 283 | void TreeWriter::Loop() {
283                  if( jet.size() < 2 && skim )
284                          continue;
285                  std::sort( jet.begin(), jet.end(), tree::EtGreater);
286 +                if( loggingVerbosity > 1 )
287 +                        std::cout << "Found " << jet.size() << " jets" << std::endl;
288 +
289  
290                  // met
291                  std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
292                  susy::MET* metobj = &(met_it->second);
293                  met = metobj->met();
294 +                met_phi = metobj->mEt.Phi();
295 +                if( loggingVerbosity > 2 )
296 +                        std::cout << " met = " << met << std::endl;
297 +
298 +                std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
299 +                susy::MET* type1metobj = &(type1met_it->second);
300 +                type1met = type1metobj->met();
301 +                type1met_phi = type1metobj->mEt.Phi();
302 +                if( loggingVerbosity > 2 )
303 +                        std::cout << " type1met = " << type1met << std::endl;
304  
305                  // electrons
306 <                //std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
205 <                //nElectron = eVector.size();
306 >                tree::Particle* thiselectron = new tree::Particle();
307  
308 <                // vertices
308 >                map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
309 >                if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
310 >                        cout << "gsfElectrons not found!" << endl;
311 >                } else {
312 >                        for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
313 >                                thiselectron->pt = it->momentum.Et();
314 >                                if( thiselectron->pt < 20 )
315 >                                        continue;
316 >                                if( loggingVerbosity > 2 )
317 >                                        std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
318 >                                thiselectron->eta = it->momentum.Eta();
319 >                                thiselectron->phi = it->momentum.Phi();
320 >                                electron.push_back( *thiselectron );
321 >                        }
322 >                }
323 >                if( loggingVerbosity > 1 )
324 >                        std::cout << "Found " << electron.size() << " electrons" << std::endl;
325  
326 +                // this seems not to work yet, where is the bug?
327 +                /*
328 +                std::vector<susy::Electron> eVector = event->electrons["gsfElectronsx"];
329 +                for( std::vector<susy::Electron>::iterator it = eVector.begin(); it != eVector.end(); ++it) {
330 +                        thiselectron->pt = it->momentum.Et();
331 +                        if( thiselectron->pt < 20 )
332 +                                continue;
333 +                        if( loggingVerbosity > 2 )
334 +                                std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
335 +                        thiselectron->eta = it->momentum.Eta();
336 +                        thiselectron->phi = it->momentum.Phi();
337 +                        electron.push_back( *thiselectron );
338 +                }
339 +                */
340 +
341 +                // muons
342 +                std::vector<susy::Muon> mVector = event->muons["muons"];
343 +                tree::Particle* thismuon = new tree::Particle();
344 +                for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
345 +                        thismuon->pt = it->momentum.Et();
346 +                        if( thismuon->pt < 20 )
347 +                                continue;
348 +                        thismuon->eta = it->momentum.Eta();
349 +                        thismuon->phi = it->momentum.Phi();
350 +                        muon.push_back( *thismuon );
351 +                }
352 +                if( loggingVerbosity > 1 )
353 +                        std::cout << "Found " << muon.size() << " muons" << std::endl;
354 +
355 +
356 +                // vertices
357                  nVertex = event->vertices.size();
358  
359                  if( ht < 450 && skim)
# Line 215 | Line 363 | void TreeWriter::Loop() {
363          } // for jentry
364  
365  
218 –
219 –        tree->Write();
366          outFile->cd();
367 +        tree->Write();
368          outFile->Write();
369          outFile->Close();
370   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines