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.6 by yilmaz, Mon Nov 1 21:48:31 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 60 | Line 61
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  
# Line 68 | Line 72 | using namespace std;
72   #define MAXHITS 1000000
73  
74   struct MyRecHit{
75 <
75 >  int depth[MAXHITS];
76    int n;
77  
78    float e[MAXHITS];
# Line 76 | Line 80 | struct MyRecHit{
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;
# Line 84 | Line 92 | struct MyRecHit{
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 96 | 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 117 | Line 133 | class RecHitTreeProducer : public edm::E
133  
134     edm::Handle<reco::BasicClusterCollection> bClusters;
135     edm::Handle<CaloTowerCollection> towers;
136 <
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;
# Line 129 | Line 148 | class RecHitTreeProducer : public edm::E
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     TTree* bcTree;
159     TTree* towerTree;
160 +   TTree* bkgTree;
161  
162    double cone;
163 <
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 148 | Line 182 | class RecHitTreeProducer : public edm::E
182    edm::InputTag EBSrc_;
183    edm::InputTag EESrc_;
184     edm::InputTag BCSrc_;
151 –
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 181 | Line 223 | RecHitTreeProducer::RecHitTreeProducer(c
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 <  JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
227 <  useJets_ = iConfig.getUntrackedParameter<bool>("useJets",false);
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 <
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 210 | 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 <
283 >  }
284    if(useJets_) {
285      ev.getByLabel(JetSrc_,jets);
286    }
# Line 225 | Line 289 | RecHitTreeProducer::analyze(const edm::E
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  
# Line 234 | Line 310 | RecHitTreeProducer::analyze(const edm::E
310        iSetup.get<CaloGeometryRecord>().get(pGeo);
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();
# Line 242 | Line 328 | RecHitTreeProducer::analyze(const edm::E
328       hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
329       hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
330       hfRecHit.isjet[hfRecHit.n] = false;
331 <     if(!useJets_){
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 <     hfRecHit.n++;
348 >     if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
349     }
350 <  
350 >   if(!doEbyEonly_){
351     for(unsigned int i = 0; i < hbheHits->size(); ++i){
352       const HBHERecHit & hit= (*hbheHits)[i];
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_){
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());
# Line 268 | Line 365 | RecHitTreeProducer::analyze(const edm::E
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_){
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());
# Line 288 | Line 388 | RecHitTreeProducer::analyze(const edm::E
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_){
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());
# Line 302 | Line 403 | RecHitTreeProducer::analyze(const edm::E
403       }
404       eeRecHit.n++;
405     }
406 +   }
407 +
408 +   if(doTowers_){
409 +
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 +   }
437  
438 <   if(doBasicClusters_){
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;
# Line 326 | Line 457 | RecHitTreeProducer::analyze(const edm::E
457           bcTree->Fill();
458        }
459     }
460 <  
461 <   eeTree->Fill();
462 <   ebTree->Fill();
463 <  
464 <   hbheTree->Fill();
465 <   hfTree->Fill();
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  
# Line 355 | Line 496 | RecHitTreeProducer::beginJob()
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","");
# Line 373 | Line 515 | RecHitTreeProducer::beginJob()
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");
# Line 386 | Line 541 | RecHitTreeProducer::beginJob()
541       //     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
542    }
543  
544 <
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 413 | 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