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.9 by yilmaz, Tue Nov 9 22:54:45 2010 UTC

# Line 68 | Line 68 | using namespace std;
68   #define MAXHITS 1000000
69  
70   struct MyRecHit{
71 <
71 >  int depth[MAXHITS];
72    int n;
73  
74    float e[MAXHITS];
# Line 84 | Line 84 | struct MyRecHit{
84   };
85  
86  
87 + struct MyBkg{
88 +   int n;
89 +   float rho[50];
90 +   float sigma[50];
91 + };
92 +
93  
94   //
95   // class declaration
# Line 122 | Line 128 | class RecHitTreeProducer : public edm::E
128    typedef vector<EcalRecHit>::const_iterator EcalIterator;
129    
130    edm::Handle<reco::CaloJetCollection> jets;
131 +
132 +   edm::Handle<std::vector<double> > rhos;
133 +   edm::Handle<std::vector<double> > sigmas;
134    
135    MyRecHit hbheRecHit;
136    MyRecHit hfRecHit;
# Line 129 | Line 138 | class RecHitTreeProducer : public edm::E
138    MyRecHit eeRecHit;
139     MyRecHit myBC;
140     MyRecHit myTowers;
141 +   MyBkg bkg;
142  
143 +  TNtuple* nt;
144    TTree* hbheTree;
145    TTree* hfTree;
146    TTree* ebTree;
147    TTree* eeTree;
148     TTree* bcTree;
149     TTree* towerTree;
150 +   TTree* bkgTree;
151  
152    double cone;
153 +  double hfTowerThreshold_;
154 +  double hfLongThreshold_;
155 +  double hfShortThreshold_;
156 +  double hbheThreshold_;
157 +  double ebThreshold_;
158 +  double eeThreshold_;
159  
160     edm::Service<TFileService> fs;
161     const CentralityBins * cbins_;
# Line 148 | Line 166 | class RecHitTreeProducer : public edm::E
166    edm::InputTag EBSrc_;
167    edm::InputTag EESrc_;
168     edm::InputTag BCSrc_;
151 –
169     edm::InputTag TowerSrc_;
170  
171    edm::InputTag JetSrc_;
172  
173 +   edm::InputTag FastJetTag_;
174 +
175    bool useJets_;
176     bool doBasicClusters_;
177 +   bool doTowers_;
178 +   bool doEcal_;
179 +   bool doHcal_;
180 +
181 +   bool doFastJet_;
182  
183 +  bool doEbyEonly_;
184   };
185  
186   //
# Line 182 | Line 207 | RecHitTreeProducer::RecHitTreeProducer(c
207    BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
208    TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
209    JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
210 <  useJets_ = iConfig.getUntrackedParameter<bool>("useJets",false);
210 >  useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
211    doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
212 <
212 >  doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
213 >  doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
214 >  doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
215 >  doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
216 >  FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
217 >  doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
218 >  hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
219 >  hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",3.);
220 >  hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",3.);
221   }
222  
223  
# Line 210 | Line 243 | RecHitTreeProducer::analyze(const edm::E
243    hbheRecHit.n = 0;
244    ebRecHit.n = 0;
245    eeRecHit.n = 0;
246 +  myBC.n = 0;
247 +   myTowers.n = 0;
248 +   bkg.n = 0;
249  
250 +   if(doEcal_){
251    ev.getByLabel(EBSrc_,ebHits);
252    ev.getByLabel(EESrc_,eeHits);
253 <
253 >   }
254 >  if(doHcal_){
255    ev.getByLabel(HcalRecHitHFSrc_,hfHits);
256    ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
257 <
257 >  }
258    if(useJets_) {
259      ev.getByLabel(JetSrc_,jets);
260    }
# Line 225 | Line 263 | RecHitTreeProducer::analyze(const edm::E
263       ev.getByLabel(BCSrc_,bClusters);
264    }
265  
266 +  if(doTowers_){
267 +     ev.getByLabel(TowerSrc_,towers);
268 +  }
269  
270 +  if(doFastJet_){
271 +     ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
272 +     ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
273 +     bkg.n = rhos->size();
274 +     for(unsigned int i = 0; i < rhos->size(); ++i){
275 +        bkg.rho[i] = (*rhos)[i];
276 +        bkg.sigma[i] = (*sigmas)[i];
277 +     }
278 +  }
279    
280    if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
281  
# Line 234 | Line 284 | RecHitTreeProducer::analyze(const edm::E
284        iSetup.get<CaloGeometryRecord>().get(pGeo);
285        geo = pGeo.product();
286     }
287 +
288 +   int nHFlong = 0;
289 +   int nHFshort = 0;
290 +   int nHFtower = 0;
291    
292 +   if(doHcal_){
293     for(unsigned int i = 0; i < hfHits->size(); ++i){
294       const HFRecHit & hit= (*hfHits)[i];
295       hfRecHit.e[hfRecHit.n] = hit.energy();
# Line 242 | Line 297 | RecHitTreeProducer::analyze(const edm::E
297       hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
298       hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
299       hfRecHit.isjet[hfRecHit.n] = false;
300 <     if(!useJets_){
300 >     hfRecHit.depth[hfRecHit.n] = hit.id().depth();
301 >     if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshort++;
302 >     if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlong++;
303 >
304 >     if(useJets_){
305         for(unsigned int j = 0 ; j < jets->size(); ++j){
306           const reco::Jet& jet = (*jets)[j];
307           double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
# Line 251 | Line 310 | RecHitTreeProducer::analyze(const edm::E
310       }
311       hfRecHit.n++;
312     }
313 <  
313 >   if(!doEbyEonly_){
314     for(unsigned int i = 0; i < hbheHits->size(); ++i){
315       const HBHERecHit & hit= (*hbheHits)[i];
316       hbheRecHit.e[hbheRecHit.n] = hit.energy();
# Line 259 | Line 318 | RecHitTreeProducer::analyze(const edm::E
318       hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
319       hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
320       hbheRecHit.isjet[hbheRecHit.n] = false;
321 <     if(!useJets_){
321 >     if(useJets_){
322         for(unsigned int j = 0 ; j < jets->size(); ++j){
323           const reco::Jet& jet = (*jets)[j];
324           double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
# Line 268 | Line 327 | RecHitTreeProducer::analyze(const edm::E
327       }
328       hbheRecHit.n++;
329     }
330 <  
330 >   }
331 >   }
332 >   if(doEcal_ && !doEbyEonly_){
333     for(unsigned int i = 0; i < ebHits->size(); ++i){
334       const EcalRecHit & hit= (*ebHits)[i];
335       ebRecHit.e[ebRecHit.n] = hit.energy();
# Line 276 | Line 337 | RecHitTreeProducer::analyze(const edm::E
337       ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
338       ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
339       ebRecHit.isjet[ebRecHit.n] = false;
340 <     if(!useJets_){
340 >     if(useJets_){
341         for(unsigned int j = 0 ; j < jets->size(); ++j){
342           const reco::Jet& jet = (*jets)[j];
343           double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
# Line 293 | Line 354 | RecHitTreeProducer::analyze(const edm::E
354       eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
355       eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
356       eeRecHit.isjet[eeRecHit.n] = false;
357 <     if(!useJets_){
357 >     if(useJets_){
358         for(unsigned int j = 0 ; j < jets->size(); ++j){
359           const reco::Jet& jet = (*jets)[j];
360           double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
# Line 302 | Line 363 | RecHitTreeProducer::analyze(const edm::E
363       }
364       eeRecHit.n++;
365     }
366 +   }
367 +
368 +   if(doTowers_){
369 +
370 +      for(unsigned int i = 0; i < towers->size(); ++i){
371 +      const CaloTower & hit= (*towers)[i];
372 +      myTowers.e[myTowers.n] = hit.energy();
373 +      myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
374 +      myTowers.eta[myTowers.n] = getEta(hit.id());
375 +      myTowers.phi[myTowers.n] = getPhi(hit.id());
376 +      myTowers.isjet[myTowers.n] = false;
377 +
378 +      if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtower++;
379 +
380 +      if(useJets_){
381 +         for(unsigned int j = 0 ; j < jets->size(); ++j){
382 +            const reco::Jet& jet = (*jets)[j];
383 +            double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
384 +            if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
385 +         }
386 +      }
387 +      myTowers.n++;
388 +      }
389  
390 +   }
391  
392 <   if(doBasicClusters_){
392 >   if(doBasicClusters_ && !doEbyEonly_){
393        for(unsigned int j = 0 ; j < jets->size(); ++j){
394           const reco::Jet& jet = (*jets)[j];
395           myBC.n = 0;
# Line 326 | Line 411 | RecHitTreeProducer::analyze(const edm::E
411           bcTree->Fill();
412        }
413     }
414 <  
415 <   eeTree->Fill();
416 <   ebTree->Fill();
417 <  
418 <   hbheTree->Fill();
419 <   hfTree->Fill();
414 >
415 >   if(!doEbyEonly_){
416 >     towerTree->Fill();
417 >    
418 >     eeTree->Fill();
419 >     ebTree->Fill();
420 >    
421 >     hbheTree->Fill();
422 >     hfTree->Fill();
423 >      
424 >     if (doFastJet_) {
425 >       bkgTree->Fill();
426 >     }
427 >   }
428 >
429 >   nt->Fill(nHFtower,nHFlong,nHFshort);
430    
431   }
432  
# Line 355 | Line 450 | RecHitTreeProducer::beginJob()
450    hfTree->Branch("et",hfRecHit.et,"et[n]/F");
451    hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
452    hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
453 +  hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
454    hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
455  
456    eeTree = fs->make<TTree>("ee","");
# Line 373 | Line 469 | RecHitTreeProducer::beginJob()
469    ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
470    ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
471  
472 +  towerTree = fs->make<TTree>("tower","");
473 +  towerTree->Branch("n",&myTowers.n,"n/I");
474 +  towerTree->Branch("e",myTowers.e,"e[n]/F");
475 +  towerTree->Branch("et",myTowers.et,"et[n]/F");
476 +  towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
477 +  towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
478 +  towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
479 +
480 +
481    if(doBasicClusters_){
482       bcTree = fs->make<TTree>("bc","");
483       bcTree->Branch("n",&myBC.n,"n/I");
# Line 386 | Line 491 | RecHitTreeProducer::beginJob()
491       //     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
492    }
493  
494 <
494 >  if(doFastJet_){
495 >     bkgTree = fs->make<TTree>("bkg","");
496 >     bkgTree->Branch("n",&bkg.n,"n/I");
497 >     bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
498 >     bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
499 >  }
500 >  
501 >  nt = fs->make<TNtuple>("ntEvent","","nHF:nHFlong:nHFshort");
502  
503   }
504  

Diff Legend

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