ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/RecHitTreeProducer.cc
(Generate patch)

Comparing UserCode/CmsHi/JetAnalysis/src/RecHitTreeProducer.cc (file contents):
Revision 1.2 by yilmaz, Mon Oct 18 16:13:37 2010 UTC vs.
Revision 1.14 by frankma, Thu Oct 13 10:40:53 2011 UTC

# Line 12 | Line 12
12   */
13   //
14   // Original Author:  Yetkin Yilmaz
15 + // Modified: Frank Ma
16   //         Created:  Tue Sep  7 11:38:19 EDT 2010
17   // $Id$
18   //
# Line 57 | Line 58
58   #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
59   #include "DataFormats/HcalDetId/interface/HcalDetId.h"
60   #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
61 + #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
62   #include "DataFormats/DetId/interface/DetId.h"
63  
64 + #include "DataFormats/VertexReco/interface/Vertex.h"
65 + #include "DataFormats/VertexReco/interface/VertexFwd.h"
66 +
67  
68   #include "TNtuple.h"
69  
70   using namespace std;
71  
72 < #define MAXHITS 100000
72 > #define MAXHITS 1000000
73  
74   struct MyRecHit{
75 <
75 >  int depth[MAXHITS];
76    int n;
77  
78    float e[MAXHITS];
79    float et[MAXHITS];
80    float eta[MAXHITS];
81    float phi[MAXHITS];
82 +  bool isjet[MAXHITS];
83 +  float etvtx[MAXHITS];
84 +  float etavtx[MAXHITS];
85 +  float emEtVtx[MAXHITS];
86 +  float hadEtVtx[MAXHITS];
87 +
88 +   float jtpt;
89 +   float jteta;
90 +   float jtphi;
91  
92   };
93  
94  
95 + struct MyBkg{
96 +   int n;
97 +   float rho[50];
98 +   float sigma[50];
99 + };
100 +
101  
102   //
103   // class declaration
# Line 90 | Line 110 | class RecHitTreeProducer : public edm::E
110    double       getEt(const DetId &id, double energy);
111    double       getEta(const DetId &id);
112    double       getPhi(const DetId &id);
113 +  reco::Vertex::Point getVtx(const edm::Event& ev);
114 +
115  
116  
117     private:
# Line 99 | Line 121 | class RecHitTreeProducer : public edm::E
121  
122        // ----------member data ---------------------------
123  
102
124     edm::Handle<reco::Centrality> cent;
125     edm::Handle<vector<double> > ktRhos;
126     edm::Handle<vector<double> > akRhos;
# Line 110 | Line 131 | class RecHitTreeProducer : public edm::E
131    edm::Handle<HFRecHitCollection> hfHits;
132    edm::Handle<HBHERecHitCollection> hbheHits;
133  
134 <   typedef vector<EcalRecHit>::const_iterator EcalIterator;
135 <
136 <   edm::Handle<reco::CaloJetCollection> signalJets;
137 <
134 >   edm::Handle<reco::BasicClusterCollection> bClusters;
135 >   edm::Handle<CaloTowerCollection> towers;
136 >  edm::Handle<reco::VertexCollection> vtxs;
137 >
138 >  typedef vector<EcalRecHit>::const_iterator EcalIterator;
139 >  
140 >  edm::Handle<reco::CaloJetCollection> jets;
141 >
142 >   edm::Handle<std::vector<double> > rhos;
143 >   edm::Handle<std::vector<double> > sigmas;
144 >  
145    MyRecHit hbheRecHit;
146    MyRecHit hfRecHit;
147    MyRecHit ebRecHit;
148    MyRecHit eeRecHit;
149 +   MyRecHit myBC;
150 +   MyRecHit myTowers;
151 +   MyBkg bkg;
152  
153 +  TNtuple* nt;
154    TTree* hbheTree;
155    TTree* hfTree;
156    TTree* ebTree;
157    TTree* eeTree;
158 <  double cone;
158 >   TTree* bcTree;
159 >   TTree* towerTree;
160 >   TTree* bkgTree;
161  
162 +  double cone;
163 +  double hfTowerThreshold_;
164 +  double hfLongThreshold_;
165 +  double hfShortThreshold_;
166 +  double hbheThreshold_;
167 +  double ebThreshold_;
168 +  double eeThreshold_;
169 +  
170 +  double hbhePtMin_;
171 +  double hfPtMin_;
172 +  double ebPtMin_;
173 +  double eePtMin_;
174 +  double towerPtMin_;
175 +  
176     edm::Service<TFileService> fs;
177     const CentralityBins * cbins_;
178     const CaloGeometry *geo;
# Line 133 | Line 181 | class RecHitTreeProducer : public edm::E
181    edm::InputTag HcalRecHitHBHESrc_;
182    edm::InputTag EBSrc_;
183    edm::InputTag EESrc_;
184 +   edm::InputTag BCSrc_;
185 +   edm::InputTag TowerSrc_;
186 +  edm::InputTag VtxSrc_;
187 +
188 +  edm::InputTag JetSrc_;
189  
190 +   edm::InputTag FastJetTag_;
191  
192 +  bool useJets_;
193 +   bool doBasicClusters_;
194 +   bool doTowers_;
195 +   bool doEcal_;
196 +   bool doHcal_;
197 +
198 +   bool doFastJet_;
199 +
200 +  bool doEbyEonly_;
201   };
202  
203   //
# Line 154 | Line 217 | RecHitTreeProducer::RecHitTreeProducer(c
217     cone(0.5)
218   {
219     //now do what ever initialization is needed
157
220    EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
221    EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
222    HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
223    HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
224 +  BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
225 +  TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
226 +  VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex"));
227 +  JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
228 +  useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
229 +  doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
230 +  doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
231 +  doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
232 +  doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
233 +  doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
234 +  FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
235 +  doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
236 +  hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
237 +  hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
238 +  hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
239 +  hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
240 +  hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
241 +  ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
242 +  eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
243 +  towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
244   }
245  
246  
# Line 184 | Line 266 | RecHitTreeProducer::analyze(const edm::E
266    hbheRecHit.n = 0;
267    ebRecHit.n = 0;
268    eeRecHit.n = 0;
269 +  myBC.n = 0;
270 +   myTowers.n = 0;
271 +   bkg.n = 0;
272  
273 +  // get vertex
274 +  reco::Vertex::Point vtx = getVtx(ev);
275 +
276 +   if(doEcal_){
277    ev.getByLabel(EBSrc_,ebHits);
278    ev.getByLabel(EESrc_,eeHits);
279 <
279 >   }
280 >  if(doHcal_){
281    ev.getByLabel(HcalRecHitHFSrc_,hfHits);
282    ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
283 <
284 <   if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
283 >  }
284 >  if(useJets_) {
285 >    ev.getByLabel(JetSrc_,jets);
286 >  }
287 >
288 >  if(doBasicClusters_){
289 >     ev.getByLabel(BCSrc_,bClusters);
290 >  }
291 >
292 >  if(doTowers_){
293 >     ev.getByLabel(TowerSrc_,towers);
294 >  }
295 >
296 >  if(doFastJet_){
297 >     ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
298 >     ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
299 >     bkg.n = rhos->size();
300 >     for(unsigned int i = 0; i < rhos->size(); ++i){
301 >        bkg.rho[i] = (*rhos)[i];
302 >        bkg.sigma[i] = (*sigmas)[i];
303 >     }
304 >  }
305 >  
306 >  if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
307  
308     if(!geo){
309        edm::ESHandle<CaloGeometry> pGeo;
# Line 199 | Line 311 | RecHitTreeProducer::analyze(const edm::E
311        geo = pGeo.product();
312     }
313  
314 +   int nHFlongPlus = 0;
315 +   int nHFshortPlus = 0;
316 +   int nHFtowerPlus = 0;
317 +   int nHFlongMinus = 0;
318 +   int nHFshortMinus = 0;
319 +   int nHFtowerMinus = 0;
320 +
321 +
322 +  
323 +   if(doHcal_){
324     for(unsigned int i = 0; i < hfHits->size(); ++i){
325       const HFRecHit & hit= (*hfHits)[i];
326       hfRecHit.e[hfRecHit.n] = hit.energy();
327       hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
328       hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
329       hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
330 <     hfRecHit.n++;
331 <   }
330 >     hfRecHit.isjet[hfRecHit.n] = false;
331 >     hfRecHit.depth[hfRecHit.n] = hit.id().depth();
332  
333 +     if(hit.id().ieta() > 0){
334 +     if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
335 +     if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
336 +     }else{
337 +       if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
338 +       if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
339 +     }
340 +
341 +     if(useJets_){
342 +       for(unsigned int j = 0 ; j < jets->size(); ++j){
343 +         const reco::Jet& jet = (*jets)[j];
344 +         double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
345 +         if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
346 +       }
347 +     }
348 +     if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
349 +   }
350 +   if(!doEbyEonly_){
351     for(unsigned int i = 0; i < hbheHits->size(); ++i){
352       const HBHERecHit & hit= (*hbheHits)[i];
353 <     hbheRecHit.e[hfRecHit.n] = hit.energy();
354 <     hbheRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
355 <     hbheRecHit.eta[hfRecHit.n] = getEta(hit.id());
356 <     hbheRecHit.phi[hfRecHit.n] = getPhi(hit.id());
353 >     if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
354 >     hbheRecHit.e[hbheRecHit.n] = hit.energy();
355 >     hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
356 >     hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
357 >     hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
358 >     hbheRecHit.isjet[hbheRecHit.n] = false;
359 >     if(useJets_){
360 >       for(unsigned int j = 0 ; j < jets->size(); ++j){
361 >         const reco::Jet& jet = (*jets)[j];
362 >         double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
363 >         if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
364 >       }
365 >     }
366       hbheRecHit.n++;
367     }
368 <
368 >   }
369 >   }
370 >   if(doEcal_ && !doEbyEonly_){
371     for(unsigned int i = 0; i < ebHits->size(); ++i){
372       const EcalRecHit & hit= (*ebHits)[i];
373 +     if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
374       ebRecHit.e[ebRecHit.n] = hit.energy();
375       ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
376       ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
377       ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
378 +     ebRecHit.isjet[ebRecHit.n] = false;
379 +     if(useJets_){
380 +       for(unsigned int j = 0 ; j < jets->size(); ++j){
381 +         const reco::Jet& jet = (*jets)[j];
382 +         double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
383 +         if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
384 +       }
385 +     }
386       ebRecHit.n++;
387     }
388 <
388 >  
389     for(unsigned int i = 0; i < eeHits->size(); ++i){
390       const EcalRecHit & hit= (*eeHits)[i];
391 +     if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
392       eeRecHit.e[eeRecHit.n] = hit.energy();
393       eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
394       eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
395       eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
396 +     eeRecHit.isjet[eeRecHit.n] = false;
397 +     if(useJets_){
398 +       for(unsigned int j = 0 ; j < jets->size(); ++j){
399 +         const reco::Jet& jet = (*jets)[j];
400 +         double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
401 +         if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
402 +       }
403 +     }
404       eeRecHit.n++;
405     }
406 +   }
407 +
408 +   if(doTowers_){
409  
410 <   eeTree->Fill();
411 <   ebTree->Fill();
410 >      for(unsigned int i = 0; i < towers->size(); ++i){
411 >      const CaloTower & hit= (*towers)[i];
412 >      if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
413 >      myTowers.e[myTowers.n] = hit.energy();
414 >      myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
415 >      myTowers.eta[myTowers.n] = getEta(hit.id());
416 >      myTowers.phi[myTowers.n] = getPhi(hit.id());
417 >      myTowers.isjet[myTowers.n] = false;
418 >      myTowers.etvtx[myTowers.n] = hit.p4(vtx).Et();
419 >      myTowers.etavtx[myTowers.n] = hit.p4(vtx).Eta();
420 >      myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
421 >      myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
422 >
423 >      if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
424 >      if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
425 >
426 >      if(useJets_){
427 >         for(unsigned int j = 0 ; j < jets->size(); ++j){
428 >            const reco::Jet& jet = (*jets)[j];
429 >            double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
430 >            if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
431 >         }
432 >      }
433 >      myTowers.n++;
434 >      }
435  
436 <   hbheTree->Fill();
437 <   hfTree->Fill();
436 >   }
437 >
438 >   if(doBasicClusters_ && !doEbyEonly_){
439 >      for(unsigned int j = 0 ; j < jets->size(); ++j){
440 >         const reco::Jet& jet = (*jets)[j];
441 >         myBC.n = 0;
442 >         myBC.jtpt = jet.pt();
443 >         myBC.jteta = jet.eta();
444 >         myBC.jtphi = jet.phi();
445 >
446 >         for(unsigned int i = 0; i < bClusters->size(); ++i){
447 >            const reco::BasicCluster & bc= (*bClusters)[i];
448 >            double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
449 >            if(dr < cone){
450 >               myBC.e[myBC.n] = bc.energy();
451 >               myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
452 >               myBC.eta[myBC.n] = bc.eta();
453 >               myBC.phi[myBC.n] = bc.phi();
454 >               myBC.n++;
455 >            }
456 >         }
457 >         bcTree->Fill();
458 >      }
459 >   }
460 >
461 >   if(!doEbyEonly_){
462 >     towerTree->Fill();
463 >    
464 >     eeTree->Fill();
465 >     ebTree->Fill();
466 >    
467 >     hbheTree->Fill();
468 >     hfTree->Fill();
469 >      
470 >     if (doFastJet_) {
471 >       bkgTree->Fill();
472 >     }
473 >   }
474  
475 +   nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
476 +  
477   }
478  
479  
# Line 248 | Line 481 | RecHitTreeProducer::analyze(const edm::E
481   void
482   RecHitTreeProducer::beginJob()
483   {
484 <
484 >  
485    hbheTree = fs->make<TTree>("hbhe","");
486    hbheTree->Branch("n",&hbheRecHit.n,"n/I");
487    hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
488    hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
489    hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
490    hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
491 <
491 >  hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
492 >  
493    hfTree = fs->make<TTree>("hf","");
494    hfTree->Branch("n",&hfRecHit.n,"n/I");
495    hfTree->Branch("e",hfRecHit.e,"e[n]/F");
496    hfTree->Branch("et",hfRecHit.et,"et[n]/F");
497    hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
498    hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
499 +  hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
500 +  hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
501  
502    eeTree = fs->make<TTree>("ee","");
503    eeTree->Branch("n",&eeRecHit.n,"n/I");
# Line 269 | Line 505 | RecHitTreeProducer::beginJob()
505    eeTree->Branch("et",eeRecHit.et,"et[n]/F");
506    eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
507    eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
508 <
508 >  eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
509 >
510    ebTree = fs->make<TTree>("eb","");
511    ebTree->Branch("n",&ebRecHit.n,"n/I");
512    ebTree->Branch("e",ebRecHit.e,"e[n]/F");
513    ebTree->Branch("et",ebRecHit.et,"et[n]/F");
514    ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
515    ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
516 +  ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
517 +
518 +  towerTree = fs->make<TTree>("tower","");
519 +  towerTree->Branch("n",&myTowers.n,"n/I");
520 +  towerTree->Branch("e",myTowers.e,"e[n]/F");
521 +  towerTree->Branch("et",myTowers.et,"et[n]/F");
522 +  towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
523 +  towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
524 +  towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
525 +  towerTree->Branch("etvtx",myTowers.etvtx,"etvtx[n]/F");
526 +  towerTree->Branch("etavtx",myTowers.etavtx,"etavtx[n]/F");
527 +  towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
528 +  towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
529 +
530 +
531 +  if(doBasicClusters_){
532 +     bcTree = fs->make<TTree>("bc","");
533 +     bcTree->Branch("n",&myBC.n,"n/I");
534 +     bcTree->Branch("e",myBC.e,"e[n]/F");
535 +     bcTree->Branch("et",myBC.et,"et[n]/F");
536 +     bcTree->Branch("eta",myBC.eta,"eta[n]/F");
537 +     bcTree->Branch("phi",myBC.phi,"phi[n]/F");
538 +     bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
539 +     bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
540 +     bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
541 +     //     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
542 +  }
543 +
544 +  if(doFastJet_){
545 +     bkgTree = fs->make<TTree>("bkg","");
546 +     bkgTree->Branch("n",&bkg.n,"n/I");
547 +     bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
548 +     bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
549 +  }
550 +  
551 +  nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
552  
553   }
554  
# Line 289 | Line 562 | double RecHitTreeProducer::getEt(const D
562    double et = energy*sin(pos.theta());
563    return et;
564   }
565 <
565 >
566   double RecHitTreeProducer::getEta(const DetId &id){
567    const GlobalPoint& pos=geo->getPosition(id);
568    double et = pos.eta();
# Line 302 | Line 575 | double RecHitTreeProducer::getPhi(const
575    return et;
576   }
577  
578 <
578 > reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
579 > {
580 >  ev.getByLabel(VtxSrc_,vtxs);
581 >  int greatestvtx = 0;
582 >  int nVertex = vtxs->size();
583 >  
584 >  for (unsigned int i = 0 ; i< vtxs->size(); ++i){
585 >    unsigned int daughter = (*vtxs)[i].tracksSize();
586 >    if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i;
587 >    //cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl;
588 >  }
589 >  
590 >  if(nVertex<=0){
591 >    return reco::Vertex::Point(-999,-999,-999);
592 >  }
593 >  return (*vtxs)[greatestvtx].position();
594 > }
595  
596   //define this as a plug-in
597   DEFINE_FWK_MODULE(RecHitTreeProducer);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines