ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/paUpcCode/src/UpcAna.cc
Revision: 1.1
Committed: Wed Jan 9 16:04:32 2013 UTC (12 years, 3 months ago) by mojoe137
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message: 
First check in of upc pa code.

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: UpcAna
4 // Class: UpcAna
5 //
6 /**\class UpcAna UpcAna.cc ana/UpcAna/src/UpcAna.cc
7
8 Description: [one line class summary]
9
10 Implementation:
11 [Notes on implementation]
12 */
13 //
14 // Original Author: Magdalena Malek
15 // Created: Tue Mar 20 09:53:57 EDT 2012
16 // $Id$
17 //
18 //
19
20
21 // system include files
22 #include <memory>
23 #include <iostream>
24 #include <string>
25 #include <sstream>
26 #include <vector>
27 #include <utility>
28
29 #include <TTree.h>
30 #include <TLorentzVector.h>
31 #include <TClonesArray.h>
32 #include "TFile.h"
33
34
35 // user include files
36 #include "FWCore/Framework/interface/Frameworkfwd.h"
37 #include "FWCore/Framework/interface/EDAnalyzer.h"
38 #include "FWCore/Framework/interface/Event.h"
39 #include "FWCore/Framework/interface/MakerMacros.h"
40 #include "FWCore/ParameterSet/interface/ParameterSet.h"
41
42 #include "DataFormats/PatCandidates/interface/Muon.h"
43 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
44 #include "DataFormats/VertexReco/interface/Vertex.h"
45 #include "DataFormats/VertexReco/interface/VertexFwd.h"
46 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
47 #include "DataFormats/TrackReco/interface/Track.h"
48 #include "DataFormats/TrackReco/interface/TrackFwd.h"
49
50 #include "DataFormats/HeavyIonEvent/interface/CentralityProvider.h"
51
52 #include "FWCore/ServiceRegistry/interface/Service.h"
53 #include "CommonTools/UtilAlgos/interface/TFileService.h"
54
55 #include "Geometry/Records/interface/IdealGeometryRecord.h"
56 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
57 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
58 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
59 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
60 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
61 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
62 #include "Geometry/Records/interface/CaloGeometryRecord.h"
63 #include "Geometry/CaloTopology/interface/HcalTopology.h"
64
65
66
67
68
69 #include "FWCore/Framework/interface/ESHandle.h"
70 #include "FWCore/Framework/interface/EventSetup.h"
71
72 #include "DataFormats/HcalRecHit/interface/CastorRecHit.h"
73 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
74 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
75 #include "DataFormats/Candidate/interface/Candidate.h"
76 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
77
78 #include "CommonTools/UtilAlgos/interface/TFileService.h"
79 #include "DataFormats/HcalDigi/interface/CastorDataFrame.h"
80 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
81
82 #include "Geometry/Records/interface/IdealGeometryRecord.h"
83 #include "Geometry/Records/interface/CastorGeometryRecord.h"
84 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
85
86 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
87 #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
88 #include "DataFormats/DetId/interface/DetId.h"
89 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
90 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
91 #include "CondFormats/HcalObjects/interface/HcalQIECoder.h"
92 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
93 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
94 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
95
96 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
97 #include "DataFormats/CaloTowers/interface/CaloTower.h"
98 #include "DataFormats/Candidate/interface/Candidate.h"
99
100
101 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
102 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
103 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
104 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
105 #include "Geometry/Records/interface/CaloGeometryRecord.h"
106 #include "Geometry/CaloTopology/interface/HcalTopology.h"
107
108
109 ///
110 #include "FWCore/Common/interface/TriggerNames.h"
111 #include "HLTrigger/HLTanalyzers/interface/HLTrigReport.h"
112 #include "DataFormats/Common/interface/Handle.h"
113 #include "DataFormats/Common/interface/TriggerResults.h"
114 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
115 #include "L1Trigger/L1ExtraFromDigis/interface/L1ExtraParticleMapProd.h"
116 #include "DataFormats/L1Trigger/interface/L1ParticleMap.h"
117 #include "DataFormats/L1Trigger/interface/L1ParticleMapFwd.h"
118 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
119 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
120 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
121 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
122 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
123 #include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h"
124 ///
125
126 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
127 #include "DataFormats/Math/interface/deltaPhi.h"
128
129
130 #include "DataFormats/CaloTowers/interface/CaloTower.h"
131 #include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
132 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
133 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
134 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
135 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
136
137 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
138
139
140 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
141 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
142 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
143 #include "DataFormats/EcalDetId/interface/EEDetId.h"
144 #include "DataFormats/EcalDetId/interface/EBDetId.h"
145
146 using namespace edm;
147 using namespace std;
148 using namespace reco;
149
150
151 //
152 // class declaration
153 //
154
155 class UpcAna : public edm::EDAnalyzer {
156 public:
157 explicit UpcAna(const edm::ParameterSet&);
158 ~UpcAna();
159
160 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
161
162
163 private:
164 virtual void beginJob() ;
165 virtual void analyze(const edm::Event&, const edm::EventSetup&);
166 virtual void endJob() ;
167
168 virtual void beginRun(edm::Run const&, edm::EventSetup const&);
169 virtual void endRun(edm::Run const&, edm::EventSetup const&);
170 virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
171 virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
172
173 const CaloGeometry* fGeo;
174 edm::ESHandle<HcalDbService> conditions;
175
176 TLorentzVector lorentzMomentum(const reco::Candidate::LorentzVector& p);
177 Float_t CalcDeltaPhi(Float_t phi1, Float_t phi2);
178
179 // ----------member data ---------------------------
180
181 edm::InputTag _thePVs;
182 string _trackCollectionTag;
183
184 TTree* upcTree;
185
186 edm::Service<TFileService> mFileServer;
187
188 //general event info
189 unsigned int nEvents;
190 unsigned int runNb;
191 unsigned int eventNb;
192 unsigned int lumiSection;
193
194 //centrality of the event
195 CentralityProvider * centrality_;
196 int centBin;
197 int theCentralityBin;
198
199 //vertex info
200 math::XYZPoint RefVtx;
201 float xVtx;
202 float yVtx;
203 float zVtx;
204 float nPV;
205
206 //tracks
207 int nTracks;
208 double track_chi2[100000];
209 double track_ndof[100000];
210 double track_charge[100000];
211 double track_x[100000];
212 double track_y[100000];
213 double track_z[100000];
214 double track_p[100000];
215 double track_pt[100000];
216 double track_px[100000];
217 double track_py[100000];
218 double track_pz[100000];
219 double track_qoverp[100000];
220 double track_lambda[100000];
221 double track_phi[100000];
222 double track_eta[100000];
223 double track_theta[100000];
224 int fnbOfTrax;
225
226 //RecHit
227 //HF detector
228 int fNumRecHits;
229 int fNumHBHERecHits;
230 int fNumECALRecHits;
231 int fNumECALEERecHits;
232
233 float fEnergy[1728];
234 float fTransverseEnergy[1728];
235 float fEta[1728];
236
237 float fHBHERecHitEta[10000];
238 float fHBHERecHitEnergy[10000];
239
240 float fECALRecHitEta[100000];
241 float fECALRecHitEnergy[100000];
242
243 float fECALEERecHitEta[100000];
244 float fECALEERecHitEnergy[100000];
245
246
247 //add calotowers
248 int fNumCaloTowers;
249 float fCaloTowerEnergy[100000];
250 float fCaloTowerEta[100000];
251 float fCaloTowerTransverseEnergy[100000];
252 int fieta[100000];
253
254
255 float fCaloTowerEMEnergy[100000];
256 float fCaloTowerHADEnergy[100000];
257
258
259 int fNumCaloTowersHF;
260 float fHFCaloTowerEnergy[100000];
261 float fHFCaloTowerEta[100000];
262
263
264 int fNumCaloTowersHB;
265 float fHBCaloTowerEnergy[100000];
266 float fHBCaloTowerEta[100000];
267
268 int fNumCaloTowersHE;
269 float fHECaloTowerEnergy[100000];
270 float fHECaloTowerEta[100000];
271
272 int fNumCaloTowersEB;
273 float fEBCaloTowerEnergy[100000];
274 float fEBCaloTowerEta[100000];
275
276
277 int fNumCaloTowersEE;
278 float fEECaloTowerEnergy[100000];
279 float fEECaloTowerEta[100000];
280
281
282 int l1AcceptBit;
283
284
285 //tracker muons
286 float Reco_mu_pt[100000];
287 float Reco_mu_eta[100000];
288 float Reco_mu_phi[100000];
289 int Reco_mu_charge[100000];
290 float Reco_mu_rapidity[100000];
291 float Reco_mu_p[100000];
292 float Reco_mu_energy[100000];
293 float Reco_mu_et[100000];
294 float Reco_mu_px[100000];
295 float Reco_mu_py[100000];
296 float Reco_mu_pz[100000];
297 float Reco_mu_theta[100000];
298 int fNBOfMuons;
299
300
301
302 float ResMass;
303 float ResPt;
304 float ResY;
305 float ResPhi;
306
307
308
309
310
311
312 };
313
314 //
315 // constants, enums and typedefs
316 //
317
318 //
319 // static data member definitions
320 //
321
322 //
323 // constructors and destructor
324 //
325 UpcAna::UpcAna(const edm::ParameterSet& iConfig)
326 {
327 //now do what ever initialization is needed
328 _thePVs = iConfig.getParameter<edm::InputTag>("primaryVertexTag");
329 _trackCollectionTag = iConfig.getParameter<std::string>("trackCollection");
330 }
331
332
333 UpcAna::~UpcAna()
334 {
335
336 // do anything here that needs to be done at desctruction time
337 // (e.g. close files, deallocate resources etc.)
338
339 }
340
341
342 //
343 // member functions
344 //
345
346 // ------------ method called for each event ------------
347 void
348 UpcAna::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
349 {
350
351
352 for (int i = 0; i < 100000; ++i)
353 {
354 track_chi2[i] = -999.;
355 track_ndof[i] = 999.;
356 track_charge[i] = 999.;
357 track_x[i] = -999.;
358 track_y[i] = -999.;
359 track_z[i] = -999.;
360 track_p[i] = -999.;
361 track_pt[i] = -999.;
362 track_px[i] = -999.;
363 track_py[i] = -999.;
364 track_pz[i] = -999.;
365 track_qoverp[i] = -999.;
366 track_lambda[i] = -999.;
367 track_phi[i] = -999.;
368 track_eta[i] = -999.;
369 track_theta[i] = 999.;
370
371
372 Reco_mu_pt[i] = 999.;
373 Reco_mu_eta[i] = 999.;
374 Reco_mu_phi[i] = 999.;
375 Reco_mu_charge[i] = 999.;
376 Reco_mu_rapidity[i] = 999.;
377 Reco_mu_p[i] = 999.;
378 Reco_mu_energy[i] = 999.;
379 Reco_mu_et[i] = 999.;
380 Reco_mu_px[i] = 999.;
381 Reco_mu_py[i] = 999.;
382 Reco_mu_pz[i] = 999.;
383 Reco_mu_theta[i] = 999.;
384
385
386
387 }
388
389
390
391
392 ResMass=-999.;
393 ResPt=-999.;
394 ResY=-999.;
395 ResPhi=-999.;
396
397 nEvents++;
398
399 runNb = iEvent.id().run();
400 eventNb = iEvent.id().event();
401 lumiSection = iEvent.luminosityBlock();
402
403 if(!centrality_) centrality_ = new CentralityProvider(iSetup);
404 centrality_->newEvent(iEvent,iSetup);
405 centBin = centrality_->getBin();
406
407
408 //trigger stuff
409
410
411 edm::Handle<L1GlobalTriggerReadoutRecord> l1GtRR;
412 iEvent.getByLabel("gtDigis",l1GtRR);
413 edm::Handle<L1GlobalTriggerObjectMapRecord> l1GtOMRec;
414 iEvent.getByLabel("hltL1GtObjectMap",l1GtOMRec);
415
416 L1GlobalTriggerReadoutRecord L1GTRR = *l1GtRR.product();
417 L1GlobalTriggerObjectMapRecord L1GTOMRec = *l1GtOMRec.product();
418
419 DecisionWord gtDecisionWord = L1GTRR.decisionWord();
420 string l1BitName;
421 int l1Accept=0;
422 // get ObjectMaps from ObjectMapRecord
423 const vector<L1GlobalTriggerObjectMap>& objMapVec = L1GTOMRec.gtObjectMap();
424 for (vector<L1GlobalTriggerObjectMap>::const_iterator itMap = objMapVec.begin();
425 itMap != objMapVec.end(); ++itMap) {
426 int iBit = (*itMap).algoBitNumber();
427 l1BitName = string( (*itMap).algoName() );
428 l1Accept = gtDecisionWord[iBit];
429
430 if (l1BitName.compare("L1_SingleMuOpen_BptxAND")==0) {
431 //cout<<"l1Decision "<<l1BitName<<" "<<l1Accept<<endl;
432 l1AcceptBit=l1Accept;
433
434 }
435 }
436
437
438
439
440
441 //VERTEX
442 edm::Handle<reco::VertexCollection> privtxs;
443 iEvent.getByLabel(_thePVs, privtxs);
444 reco::VertexCollection::const_iterator privtx;
445 nPV = privtxs->size();
446
447 if ( privtxs->begin() != privtxs->end() ) {
448 privtx=privtxs->begin();
449 RefVtx = privtx->position();
450 } else {
451 RefVtx.SetXYZ(0.,0.,0.);
452 }
453
454 xVtx = RefVtx.X();
455 yVtx = RefVtx.Y();
456 zVtx = RefVtx.Z();
457
458
459
460 //Tracks
461 int nbOfTrax=0;
462 fnbOfTrax=0;
463 Handle<TrackCollection> hiTrax;
464 iEvent.getByLabel(_trackCollectionTag,hiTrax);
465 //nTracks=hiTrax->size();
466
467 for(TrackCollection::const_iterator trax = hiTrax->begin();
468 trax!=hiTrax->end();trax++){
469
470 track_chi2[nbOfTrax] = trax->chi2();
471 track_ndof[nbOfTrax] = trax->ndof();
472 track_charge[nbOfTrax] = trax->charge();
473 track_x[nbOfTrax] = trax->vx();
474 track_y[nbOfTrax] = trax->vy();
475 track_z[nbOfTrax] = trax->vz();
476 track_p[nbOfTrax] = trax->p();
477 track_pt[nbOfTrax] = trax->pt();
478 track_px[nbOfTrax] = trax->px();
479 track_py[nbOfTrax] = trax->py();
480 track_pz[nbOfTrax] = trax->pz();
481 track_qoverp[nbOfTrax] = trax->qoverp();
482 track_lambda[nbOfTrax] = trax->lambda();
483 track_phi[nbOfTrax] = trax->phi();
484 track_eta[nbOfTrax] = trax->eta();
485 track_theta[nbOfTrax] = trax->theta();
486
487 nbOfTrax++;
488 }
489 fnbOfTrax = nbOfTrax;
490
491
492
493 //HF detector
494
495
496 iSetup.get < HcalDbRecord > ().get(conditions);
497 edm::ESHandle < CaloGeometry > caloGeom;
498 iSetup.get < CaloGeometryRecord > ().get(caloGeom);
499 fGeo = caloGeom.product();
500
501
502
503
504
505 //The caloGeom needed to get eta of the rechit
506
507 //HF RecHits
508 for (int i = 0; i < 1728; ++i)
509 {
510 fEnergy[i] = -999.;
511 fTransverseEnergy[i] = -999.;
512 fEta[i] = -999.;
513 }
514
515 int nHits = 0;
516 fNumRecHits = 0;
517
518 Handle < HFRecHitCollection > hf_recHits;
519 iEvent.getByLabel("hfreco", hf_recHits);
520 //cout<<"rechit coll size "<<hf_recHits->size()<<endl;
521
522 for (HFRecHitCollection::const_iterator hhit = hf_recHits->begin(); hhit != hf_recHits->end(); hhit++) {
523 GlobalPoint fPos = fGeo->getPosition(hhit->detid());
524 fEnergy[nHits] = hhit->energy();
525 fEta[nHits] = fPos.eta();
526 fTransverseEnergy[nHits] = hhit->energy()/TMath::CosH(fPos.eta());
527
528 //cout<<nHits<<"\t"<<fEta[nHits]<<"\t"<<fEnergy[nHits]<<endl;
529
530 ++nHits;
531 }
532 fNumRecHits = nHits;
533 // cout<<"fNumRecHits "<<fNumRecHits<<endl;
534
535 for (int iii = 0; iii < 10000; ++iii)
536 {
537 fHBHERecHitEta[iii]=-999.;
538 fHBHERecHitEnergy[iii]=-999.;
539
540 fECALRecHitEta[iii]=-999.;
541 fECALRecHitEnergy[iii]=-999.;
542
543 fECALEERecHitEta[iii]=-999.;
544 fECALEERecHitEnergy[iii]=-999.;
545 }
546
547
548 //HBHE rechits
549 edm::Handle<HBHERecHitCollection> hbhe;
550 iEvent.getByLabel("hbhereco",hbhe);
551 const HBHERecHitCollection *hbhe_recHits = hbhe.failedToGet()? 0 : &*hbhe;
552 int NumHBHERecHits = 0;
553 fNumHBHERecHits = 0;
554 for(HBHERecHitCollection::const_iterator hbheItr = hbhe_recHits->begin(); hbheItr != hbhe_recHits->end(); hbheItr++)
555 {
556 GlobalPoint fPos = fGeo->getPosition(hbheItr->detid());
557 HcalDetId idHBHE = hbheItr->id();
558 fHBHERecHitEta[NumHBHERecHits]=fPos.eta();
559 fHBHERecHitEnergy[NumHBHERecHits]=hbheItr->energy();
560 ++NumHBHERecHits;
561 }
562 fNumHBHERecHits=NumHBHERecHits;
563
564
565 //ECAL rechits
566 edm::Handle<EcalRecHitCollection> EBRecHits;
567 iEvent.getByLabel("ecalRecHit","EcalRecHitsEB", EBRecHits);
568 const EcalRecHitCollection *rechits = EBRecHits.product();
569 int NumECALRecHits = 0;
570 fNumECALRecHits = 0;
571
572 for(EcalRecHitCollection::const_iterator ecalit=rechits->begin(); ecalit!=rechits->end(); ecalit++) {
573 GlobalPoint fPos = fGeo->getPosition(ecalit->detid());
574 DetId idECAL = ecalit->id();
575 fECALRecHitEta[NumECALRecHits]=fPos.eta();
576 fECALRecHitEnergy[NumECALRecHits]=ecalit->energy();
577 ++NumECALRecHits;
578 }
579 fNumECALRecHits=NumECALRecHits;
580
581 //ECAL EE rechits
582 edm::Handle<EcalRecHitCollection> EERecHits;
583 iEvent.getByLabel("ecalRecHit","EcalRecHitsEE", EERecHits);
584 const EcalRecHitCollection *rechitsEE = EERecHits.product();
585 int NumECALEERecHits = 0;
586 fNumECALEERecHits = 0;
587
588 for(EcalRecHitCollection::const_iterator ecaleeit=rechitsEE->begin(); ecaleeit!=rechitsEE->end(); ecaleeit++) {
589 GlobalPoint fPos = fGeo->getPosition(ecaleeit->detid());
590 DetId idECALEE = ecaleeit->id();
591 fECALEERecHitEta[NumECALEERecHits]=fPos.eta();
592 fECALEERecHitEnergy[NumECALEERecHits]=ecaleeit->energy();
593 ++NumECALEERecHits;
594 }
595 fNumECALEERecHits=NumECALEERecHits;
596
597
598
599
600 //calo towers
601
602 for (int i = 0; i < 10000; ++i) {
603 fCaloTowerEnergy[i] = -999.;
604 fCaloTowerEta[i] = -999;
605 fCaloTowerTransverseEnergy[i] = -999;
606 fieta[i] =-999;
607
608 fHFCaloTowerEnergy[i]=-999.;
609 fHFCaloTowerEta[i]=-999;
610
611 fCaloTowerEMEnergy[i]=-999.;
612 fCaloTowerHADEnergy[i]=-999.;
613
614
615
616 }
617
618 Handle < CaloTowerCollection > calo_h;
619 iEvent.getByLabel("towerMaker", calo_h);
620
621
622 const CaloTowerCollection *towers = calo_h.failedToGet()? 0 : &*calo_h;
623 // cout<< calo_h->size()<<endl;
624
625 int NumCaloTowers = 0;
626 fNumCaloTowers = 0;
627
628
629
630 int nTowersHF = 0;
631 int nTowersHE = 0;
632 int nTowersHB = 0;
633 int nTowersHO = 0;
634
635 int nTowersEE = 0;
636 int nTowersEB = 0;
637
638
639
640
641 for (CaloTowerCollection::const_iterator cal = towers->begin(); cal != towers->end(); cal++) {
642 // cout<<"NumCaloTowers "<<NumCaloTowers<<" consize "<<cal->constituentsSize()<<endl;
643 bool hasHCAL = false;
644 bool hasHF = false;
645 bool hasHE = false;
646 bool hasHB = false;
647 bool hasHO = false;
648 bool hasECAL = false;
649 bool hasEE = false;
650 bool hasEB = false;
651
652
653 for(size_t iconst = 0; iconst < cal->constituentsSize(); iconst++){
654 DetId detId = cal->constituent(iconst);
655 //cout<<detId.det()<<endl;
656 if(detId.det()==DetId::Hcal)
657 {
658 //cout<<"hcal "<<detId.det()<<endl;
659 HcalDetId hcalDetId(detId);
660 if(hcalDetId.subdet()==HcalForward)
661 {
662 hasHF = true;
663 //++nTowersHF;
664 //cout<<"HcalForward "<<hcalDetId.subdet()<<endl;
665 }
666 if(hcalDetId.subdet()==HcalEndcap)
667 {
668 hasHE = true;
669 //++nTowersHE;
670 //cout<<"HcalEndcap "<<hcalDetId.subdet()<<endl;
671 }
672 if(hcalDetId.subdet()==HcalBarrel)
673 {
674 hasHB = true;
675 // ++nTowersHB;
676 //cout<<"HcalBarrel "<<hcalDetId.subdet()<<" energy "<<cal->energy()<<" eta "<<cal->eta()<<endl;
677 }
678 if(hcalDetId.subdet()==HcalOuter)
679 {
680 hasHO = true;
681 //++nTowersHO;
682 //cout<<"HcalOuter "<<hcalDetId.subdet()<<" energy "<<cal->energy()<<" eta "<<cal->eta()<<endl;
683 }
684 }
685
686 if(detId.det()==DetId::Ecal)
687 {
688 //cout<<"ecal "<<detId.det()<<endl;
689 EcalSubdetector ecalSubDet = (EcalSubdetector)detId.subdetId();
690 if(ecalSubDet == EcalEndcap)
691 {
692 hasEE = true;
693 // cout<<"EcalEndcap "<<ecalSubDet<<endl;
694 }
695 if(ecalSubDet == EcalBarrel)
696 {
697 hasEB = true;
698 //cout<<"EcalBarrel "<<ecalSubDet<<endl;
699 }
700 }
701 }//costituants
702
703
704 if( hasHF && !hasHE )
705 {
706 //cout<<" hasHF "<<hasHF<<" hasHE "<<hasHE<<endl;
707 //cout<<cal->energy()<<"\t"<<cal->eta()<<endl;
708 fHFCaloTowerEnergy[nTowersHF]=cal->energy();
709 fHFCaloTowerEta[nTowersHF]=cal->eta();
710 ++nTowersHF;
711 }
712
713 if( hasHE && !hasHF && !hasHB )
714 {
715 //cout<<" hasHE "<<hasHE<<" hasHF "<<hasHF<<" hasHB "<<hasHB<<endl;
716 //cout<<cal->energy()<<"\t"<<cal->eta()<<endl;
717 fHECaloTowerEnergy[nTowersHE]=cal->energy();
718 fHECaloTowerEta[nTowersHE]=cal->eta();
719 ++nTowersHE;
720 }
721
722
723 if( hasHB && !hasHE)
724 {
725 //cout<<" hasHB "<<hasHB<<" hasHE "<<hasHE<<endl;
726 //cout<<cal->energy()<<"\t"<<cal->eta()<<endl;
727 fHBCaloTowerEnergy[nTowersHB]=cal->energy();
728 fHBCaloTowerEta[nTowersHB]=cal->eta();
729 ++nTowersHB;
730 }
731
732
733 if( hasEE && !hasEB)
734 {
735 //cout<<" hasEE "<<hasEE<<" hasEB "<<hasEB<<endl;
736 //cout<<cal->energy()<<"\t"<<cal->eta()<<endl;
737 fEECaloTowerEnergy[nTowersEE]=cal->energy();
738 fEECaloTowerEta[nTowersEE]=cal->eta();
739 ++nTowersEE;
740 }
741
742 if( hasEB && !hasEE)
743 {
744 //cout<<" hasEB "<<hasEB<<" hasEE "<<hasEE<<endl;
745 //cout<<cal->energy()<<"\t"<<cal->eta()<<endl;
746 fEBCaloTowerEnergy[nTowersEB]=cal->energy();
747 fEBCaloTowerEta[nTowersEB]=cal->eta();
748 ++nTowersEB;
749 }
750
751
752 fCaloTowerEnergy[NumCaloTowers] = cal->energy();
753 fCaloTowerTransverseEnergy[NumCaloTowers] = cal->et();
754 fCaloTowerEta[NumCaloTowers] = cal->eta();
755 fieta[NumCaloTowers] = cal->ieta();
756 fCaloTowerEMEnergy[NumCaloTowers]=cal->emEnergy();
757 fCaloTowerHADEnergy[NumCaloTowers]=cal->hadEnergy();
758
759 ++NumCaloTowers;
760 }
761 fNumCaloTowers = NumCaloTowers;
762
763 fNumCaloTowersHF = nTowersHF;
764 //cout<<"fNumCaloTowersHF "<<fNumCaloTowersHF<<endl;
765 fNumCaloTowersHE = nTowersHE;
766 //cout<<"fNumCaloTowersHE "<<fNumCaloTowersHE<<endl;
767 fNumCaloTowersHB=nTowersHB;
768 //cout<<"fNumCaloTowersHB "<<fNumCaloTowersHB<<endl;
769 fNumCaloTowersEE=nTowersEE;
770 //cout<<"fNumCaloTowersEE "<<fNumCaloTowersEE<<endl;
771 fNumCaloTowersEB=nTowersEB;
772 //cout<<"fNumCaloTowersEB "<<fNumCaloTowersEB<<endl;
773
774
775 // &&&&&&&&&& ZDC
776
777
778
779
780 //cout<<"nb of tracks : "<<fnbOfTrax<<endl;
781
782
783 TLorentzVector vtmuon;
784 int NBOfMuons=0;
785 Handle<reco::MuonCollection> gMuons;
786 iEvent.getByLabel("muons",gMuons);
787 for (reco::MuonCollection::const_iterator global = gMuons->begin(); global != gMuons->end(); ++global ) {
788 if(global->isTrackerMuon())
789 {
790 //cout<<"tmuon "<<global->pt()<<"\t"<<global->eta()<<"\t"<<global->phi()<<"\t"<<global->energy()<<endl;
791 vtmuon.SetPtEtaPhiE(global->pt(),global->eta(),global->phi(),global->energy());
792 //cout<<"vtmuon[0] "<<vtmuon[0]<<endl;
793
794
795
796
797 Reco_mu_pt[NBOfMuons] = global->pt();
798 Reco_mu_eta[NBOfMuons] = global->eta();
799 Reco_mu_phi[NBOfMuons] = global->phi();
800 Reco_mu_charge[NBOfMuons] = global->charge();
801 Reco_mu_rapidity[NBOfMuons] = global->rapidity();
802 Reco_mu_p[NBOfMuons] = global->p();
803 Reco_mu_energy[NBOfMuons] = global->energy();
804 Reco_mu_et[NBOfMuons] = global->et();
805 Reco_mu_px[NBOfMuons] = global->px();
806 Reco_mu_py[NBOfMuons] = global->py();
807 Reco_mu_pz[NBOfMuons] = global->pz();
808 Reco_mu_theta[NBOfMuons] = global->theta();
809
810 NBOfMuons++;
811 }
812
813 }
814
815 fNBOfMuons=NBOfMuons;
816
817 if( (fnbOfTrax==2) && (fNBOfMuons==2) )
818 {
819 TLorentzVector res(Reco_mu_px[0]+Reco_mu_py[1],Reco_mu_py[0]+Reco_mu_py[1],Reco_mu_pz[0]+Reco_mu_pz[1],Reco_mu_energy[0]+Reco_mu_energy[1]);
820 ResMass=res.M();
821 ResPt=res.Pt();
822 ResY=res.Rapidity();
823 ResPhi=res.Phi();
824 }
825 upcTree->Fill();
826
827 }
828
829
830
831 TLorentzVector
832 UpcAna::lorentzMomentum(const reco::Candidate::LorentzVector& p) {
833 TLorentzVector res;
834 res.SetPtEtaPhiM(p.pt(), p.eta(), p.phi(), p.mass());
835
836 return res;
837 }
838
839
840 Float_t
841 UpcAna::CalcDeltaPhi(Float_t phi1, Float_t phi2) {
842 Float_t result = phi1 - phi2;
843 if (result > TMath::Pi()) result -= 2*TMath::Pi();
844 while (result <= -TMath::Pi()) result += 2*TMath::Pi();
845 return result;
846 }
847
848 // ------------ method called once each job just before starting event loop ------------
849 void
850 UpcAna::beginJob()
851 {
852
853
854
855 nEvents=0;
856 centrality_ = 0;
857
858 edm::Service<TFileService> mFileServer;
859 mFileServer->file().cd();
860
861
862
863
864
865
866
867
868 //ttree per event...filled only when there is at least 1 composite candidates
869 upcTree = new TTree("AnaTreeUPC","AnaTreeUPC");
870
871 upcTree->Branch("eventNb", &eventNb, "eventNb/i");
872 upcTree->Branch("runNb", &runNb, "runNb/i");
873 upcTree->Branch("LS", &lumiSection, "LS/i");
874 upcTree->Branch("Centrality", &centBin, "Centrality/I");
875 upcTree->Branch("nEvents", &nEvents, "nEvents/I");
876
877
878
879
880
881
882 upcTree->Branch("xVtx", &xVtx, "xVtx/F");
883 upcTree->Branch("yVtx", &yVtx, "yVtx/F");
884 upcTree->Branch("zVtx", &zVtx, "zVtx/F");
885
886 upcTree->Branch("fnbOfTrax",&fnbOfTrax,"fnbOfTrax/I");
887 upcTree->Branch("track_ndof",track_ndof,"track_ndof[fnbOfTrax]/D");
888 upcTree->Branch("track_chi2",track_chi2,"track_chi2[fnbOfTrax]/D");
889 upcTree->Branch("track_charge",track_charge,"track_charge[fnbOfTrax]/D");
890 upcTree->Branch("track_x",track_x,"track_x[fnbOfTrax]/D");
891 upcTree->Branch("track_y",track_y,"track_y[fnbOfTrax]/D");
892 upcTree->Branch("track_z",track_z,"track_z[fnbOfTrax]/D");
893 upcTree->Branch("track_p",track_p,"track_p[fnbOfTrax]/D");
894 upcTree->Branch("track_pt",track_pt,"track_pt[fnbOfTrax]/D");
895 upcTree->Branch("track_px",track_px,"track_px[fnbOfTrax]/D");
896 upcTree->Branch("track_py",track_py,"track_py[fnbOfTrax]/D");
897 upcTree->Branch("track_pz",track_pz,"track_pz[fnbOfTrax]/D");
898 upcTree->Branch("track_qoverp",track_qoverp,"track_qoverp[fnbOfTrax]/D");
899 upcTree->Branch("track_lambda",track_lambda,"track_lambda[fnbOfTrax]/D");
900 upcTree->Branch("track_phi",track_phi,"track_phi[fnbOfTrax]/D");
901 upcTree->Branch("track_eta",track_eta,"track_eta[fnbOfTrax]/D");
902 upcTree->Branch("track_theta",track_theta,"track_theta[fnbOfTrax]/D");
903
904
905
906
907
908
909 upcTree->Branch("fNumRecHits",&fNumRecHits,"fNumRecHits/I");
910 upcTree->Branch("fEnergy",fEnergy,"fEnergy[fNumRecHits]/F");
911 upcTree->Branch("fTransverseEnergy",fTransverseEnergy,"fTransverseEnergy[fNumRecHits]/F");
912 upcTree->Branch("fEta",fEta,"fEta[fNumRecHits]/F");
913
914 upcTree->Branch("NumCaloTowers", &fNumCaloTowers, "NumCaloTowers/I");
915 upcTree->Branch("CaloTowerTransverseEnergy", fCaloTowerTransverseEnergy, "fCaloTowerTransverseEnergy[NumCaloTowers]/F");
916 upcTree->Branch("CaloTowerEnergy", fCaloTowerEnergy, "CaloTowerEnergy[NumCaloTowers]/F");
917 upcTree->Branch("CaloTowerEta", fCaloTowerEta, "CaloTowerEta[NumCaloTowers]/F");
918 upcTree->Branch("CaloToweriEta", fieta, "fieta[NumCaloTowers]/I");
919 upcTree->Branch("fCaloTowerEMEnergy", fCaloTowerEMEnergy, "fCaloTowerEMEnergy[NumCaloTowers]/F");
920 upcTree->Branch("fCaloTowerHADEnergy", fCaloTowerHADEnergy, "fCaloTowerHADEnergy[NumCaloTowers]/F");
921
922
923 upcTree->Branch("fNumCaloTowersHF", &fNumCaloTowersHF, "fNumCaloTowersHF/I");
924 upcTree->Branch("HFCaloTowerEnergy", fHFCaloTowerEnergy, "fHFCaloTowerEnergy[fNumCaloTowersHF]/F");
925 upcTree->Branch("HFCaloTowerEta", fHFCaloTowerEta, "fHFCaloTowerEta[fNumCaloTowersHF]/F");
926
927 upcTree->Branch("fNumCaloTowersHE", &fNumCaloTowersHE, "fNumCaloTowersHE/I");
928 upcTree->Branch("HECaloTowerEnergy", fHECaloTowerEnergy, "fHECaloTowerEnergy[fNumCaloTowersHE]/F");
929 upcTree->Branch("HECaloTowerEta", fHECaloTowerEta, "fHECaloTowerEta[fNumCaloTowersHE]/F");
930
931 upcTree->Branch("fNumCaloTowersHB", &fNumCaloTowersHB, "fNumCaloTowersHB/I");
932 upcTree->Branch("HBCaloTowerEnergy", fHBCaloTowerEnergy, "fHBCaloTowerEnergy[fNumCaloTowersHB]/F");
933 upcTree->Branch("HBCaloTowerEta", fHBCaloTowerEta, "fHBCaloTowerEta[fNumCaloTowersHB]/F");
934
935 upcTree->Branch("fNumCaloTowersEE", &fNumCaloTowersEE, "fNumCaloTowersEE/I");
936 upcTree->Branch("EECaloTowerEnergy", fEECaloTowerEnergy, "fEECaloTowerEnergy[fNumCaloTowersEE]/F");
937 upcTree->Branch("EECaloTowerEta", fEECaloTowerEta, "fEECaloTowerEta[fNumCaloTowersEE]/F");
938
939 upcTree->Branch("fNumCaloTowersEB", &fNumCaloTowersEB, "fNumCaloTowersEB/I");
940 upcTree->Branch("EBCaloTowerEnergy", fEBCaloTowerEnergy, "fEBCaloTowerEnergy[fNumCaloTowersEB]/F");
941 upcTree->Branch("EBCaloTowerEta", fEBCaloTowerEta, "fEBCaloTowerEta[fNumCaloTowersEB]/F");
942
943 upcTree->Branch("fNBOfMuons", &fNBOfMuons, "fNBOfMuons/I");
944 upcTree->Branch("Reco_mu_pt", Reco_mu_pt, "Reco_mu_pt[fNBOfMuons]/F");
945 upcTree->Branch("Reco_mu_eta", Reco_mu_eta, "Reco_mu_eta[fNBOfMuons]/F");
946 upcTree->Branch("Reco_mu_phi", Reco_mu_phi, "Reco_mu_phi[fNBOfMuons]/F");
947 upcTree->Branch("Reco_mu_charge", Reco_mu_charge, "Reco_mu_charge[fNBOfMuons]/I");
948 upcTree->Branch("Reco_mu_rapidity", Reco_mu_rapidity, "Reco_mu_rapidity[fNBOfMuons]/F");
949 upcTree->Branch("Reco_mu_p", Reco_mu_p, "Reco_mu_p[fNBOfMuons]/F");
950 upcTree->Branch("Reco_mu_energy", Reco_mu_energy, "Reco_mu_energy[fNBOfMuons]/F");
951 upcTree->Branch("Reco_mu_et", Reco_mu_et, "Reco_mu_et[fNBOfMuons]/F");
952 upcTree->Branch("Reco_mu_px", Reco_mu_px, "Reco_mu_px[fNBOfMuons]/F");
953 upcTree->Branch("Reco_mu_py", Reco_mu_py, "Reco_mu_py[fNBOfMuons]/F");
954 upcTree->Branch("Reco_mu_pz", Reco_mu_pz, "Reco_mu_pz[fNBOfMuons]/F");
955 upcTree->Branch("Reco_mu_theta", Reco_mu_theta, "Reco_mu_theta[fNBOfMuons]/F");
956
957
958 upcTree->Branch("ResMass", &ResMass, "ResMass/F");
959 upcTree->Branch("ResPt", &ResPt, "ResPt/F");
960 upcTree->Branch("ResY", &ResY, "ResY/F");
961 upcTree->Branch("ResPhi", &ResPhi, "ResPhi/F");
962
963
964
965
966
967
968
969
970
971 }
972
973 // ------------ method called once each job just after ending the event loop ------------
974 void
975 UpcAna::endJob()
976 {
977
978 }
979
980 // ------------ method called when starting to processes a run ------------
981 void
982 UpcAna::beginRun(edm::Run const&, edm::EventSetup const&)
983 {
984 }
985
986 // ------------ method called when ending the processing of a run ------------
987 void
988 UpcAna::endRun(edm::Run const&, edm::EventSetup const&)
989 {
990 }
991
992 // ------------ method called when starting to processes a luminosity block ------------
993 void
994 UpcAna::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
995 {
996 }
997
998 // ------------ method called when ending the processing of a luminosity block ------------
999 void
1000 UpcAna::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
1001 {
1002 }
1003
1004 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
1005 void
1006 UpcAna::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1007 //The following says we do not know what parameters are allowed so do no validation
1008 // Please change this to state exactly what you do use, even if it is no parameters
1009 edm::ParameterSetDescription desc;
1010 desc.setUnknown();
1011 descriptions.addDefault(desc);
1012 }
1013
1014 //define this as a plug-in
1015 DEFINE_FWK_MODULE(UpcAna);
1016
1017