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

# User Rev Content
1 mojoe137 1.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