CMS 3D CMS Logo
EcalEBTrigPrimAnalyzer.cc
Go to the documentation of this file.
1 // system include files
2 #include <memory>
3 #include <utility>
4 
5 // user include files
6 #include "FWCore/Framework/interface/EDAnalyzer.h"
7 #include "FWCore/Framework/interface/Event.h"
8 #include "FWCore/Framework/interface/MakerMacros.h"
9 
10 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
11 #include "FWCore/Framework/interface/EventSetup.h"
12 #include "FWCore/Framework/interface/ESHandle.h"
13 #include "FWCore/ParameterSet/interface/ParameterSet.h"
14 
15 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
16 
17 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
18 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
19 
20 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
21 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
22 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
23 #include "Geometry/Records/interface/CaloGeometryRecord.h"
24 
25 //#include "CalibCalorimetry/EcalTPGTools/interface/EcalEBTPGScale.h"
26 #include "CalibCalorimetry/EcalTPGTools/interface/EcalTPGScale.h"
27 
28 #include "EcalEBTrigPrimAnalyzer.h"
29 
30 #include <TMath.h>
31 
32 using namespace edm;
33 class CaloSubdetectorGeometry;
34 
35 EcalEBTrigPrimAnalyzer::EcalEBTrigPrimAnalyzer(const edm::ParameterSet& iConfig)
36 
37 {
38  ecal_parts_.push_back("Barrel");
39 
40  histfile_ = new TFile("histos.root", "RECREATE");
41  tree_ = new TTree("TPGtree", "TPGtree");
42  tree_->Branch("tpIphi", &tpIphi_, "tpIphi/I");
43  tree_->Branch("tpIeta", &tpIeta_, "tpIeta/I");
44  tree_->Branch("rhIphi", &rhIphi_, "rhIphi/I");
45  tree_->Branch("rhIeta", &rhIeta_, "rhIeta/I");
46  tree_->Branch("eRec", &eRec_, "eRec/F");
47  tree_->Branch("tpgADC", &tpgADC_, "tpgADC/I");
48  tree_->Branch("tpgGeV", &tpgGeV_, "tpgGeV/F");
49  tree_->Branch("ttf", &ttf_, "ttf/I");
50  tree_->Branch("fg", &fg_, "fg/I");
51  for (unsigned int i = 0; i < ecal_parts_.size(); ++i) {
52  char title[30];
53  sprintf(title, "%s_Et", ecal_parts_[i].c_str());
54  ecal_et_[i] = new TH1I(title, "Et", 255, 0, 255);
55  sprintf(title, "%s_ttf", ecal_parts_[i].c_str());
56  ecal_tt_[i] = new TH1I(title, "TTF", 10, 0, 10);
57  sprintf(title, "%s_fgvb", ecal_parts_[i].c_str());
58  ecal_fgvb_[i] = new TH1I(title, "FGVB", 10, 0, 10);
59  }
60 
61  recHits_ = iConfig.getParameter<bool>("AnalyzeRecHits");
62  debug_ = iConfig.getParameter<bool>("Debug");
63  rechits_labelEB_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("inputRecHitsEB"));
64  primToken_ = consumes<EcalEBTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("inputTP"));
65  tokenEBdigi_ = consumes<EBDigiCollection>(iConfig.getParameter<edm::InputTag>("barrelEcalDigis"));
66  nEvents_ = 0;
67 
68  hTPvsTow_eta_ = new TH2F("TP_vs_Tow_eta", "TP vs Tow eta ; #eta(tow); #eta(tp)", 50, -2.5, 2.5, 50, -2.5, 2.5);
69  hAllTPperEvt_ = new TH1F("AllTPperEvt", "TP per Event; N_{TP}; ", 100, 0., 20000.);
70  hTPperEvt_ = new TH1F("TPperEvt", "N_{TP} per Event; N_{TP}; ", 100, 0., 500.);
71  hTP_iphiVsieta_ = new TH2F("TP_iphiVsieta", "TP i#phi vs i#eta ; i#eta(tp); i#phi(tp)", 10, 70, 80, 10, 340, 350);
72  hTP_iphiVsieta_fullrange_ =
73  new TH2F("TP_iphiVsieta_fullrange", "TP i#phi vs i#eta ; i#eta(tp); i#phi(tp)", 200, -100, 100, 350, 0, 350);
74 
75  if (recHits_) {
76  hTPvsTow_ieta_ =
77  new TH2F("TP_vs_Tow_ieta", "TP vs Tow ieta ; i#eta(tow); i#eta(tp)", 200, -100, 100, 200, -100, 100);
78 
79  hTPvsRechit_ = new TH2F("TP_vs_RecHit", "TP vs rechit Et;E_{T}(rh) (GeV);E_{T}(tp) (GeV)", 100, 0, 50, 100, 0, 50);
80  hDeltaEt_ = new TH1F("DeltaEt", "[Et(rh)-Et(TP)]/Et(rh); [E_{T}(rh)-E_{T}(tp)]/E_{T}(rh); Counts", 200, -1, 1);
81  hTPoverRechit_ = new TH1F("TP_over_RecHit", "Et(TP/rechit); E_{T}(tp)/E_{T}(rh); Counts", 200, 0, 2);
82  hRechitEt_ = new TH1F("RecHitEt", "E_{T};E_{T}(rh) (GeV);Counts", 100, 0, 50);
83  hTPEt_ = new TH1F("TPEt", "E_{T}{tp);E_{T}(rh) (GeV);Count", 100, 0, 50);
84  hRatioEt_ = new TH1F("RatioTPoverRH", "Et", 100, 0, 50);
85  hAllRechitEt_ = new TH1F("AllRecHit", "Et", 100, 0, 50);
86 
87  hRH_iphiVsieta_ = new TH2F("RH_iphiVsieta", "RH i#phi vs i#eta ; i#eta(rh); i#phi(rh)", 10, 70, 80, 10, 340, 350);
88  hRH_iphiVsieta_fullrange_ =
89  new TH2F("RH_iphiVsieta_fullrange", "RH i#phi vs i#eta ; i#eta(rh); i#phi(rh)", 200, -100, 100, 350, 0, 350);
90  }
91 }
92 
93 EcalEBTrigPrimAnalyzer::~EcalEBTrigPrimAnalyzer() {
94  // do anything here that needs to be done at desctruction time
95  // (e.g. close files, deallocate resources etc.)
96 
97  histfile_->Write();
98  histfile_->Close();
99 }
100 
101 void EcalEBTrigPrimAnalyzer::init(const edm::EventSetup& iSetup) { iSetup.get<IdealGeometryRecord>().get(eTTmap_); }
102 
103 //
104 // member functions
105 //
106 
107 // ------------ method called to analyze the data ------------
108 void EcalEBTrigPrimAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
109  using namespace edm;
110  using namespace std;
111  nEvents_++;
112 
113  if (nEvents_ == 1)
114  this->init(iSetup);
115 
116  // Get input
117  edm::Handle<EcalEBTrigPrimDigiCollection> tp;
118  iEvent.getByToken(primToken_, tp);
119  //
120  /*
121  edm::Handle<EBDigiCollection> barrelDigiHandle;
122  const EBDigiCollection *ebdigi=NULL;
123  iEvent.getByToken(tokenEBdigi_,barrelDigiHandle);
124  ebdigi=barrelDigiHandle.product();
125  */
126 
127  for (unsigned int i = 0; i < tp.product()->size(); i++) {
128  EcalEBTriggerPrimitiveDigi d = (*(tp.product()))[i];
129  int subdet = 0;
130  if (subdet == 0) {
131  ecal_et_[subdet]->Fill(d.encodedEt());
132  }
133  }
134 
135  // if (!recHits_) return;
136 
137  edm::Handle<EcalRecHitCollection> rechit_EB_col;
138  if (recHits_) {
139  // get the RecHits
140  iEvent.getByToken(rechits_labelEB_, rechit_EB_col);
141  }
142 
143  edm::ESHandle<CaloSubdetectorGeometry> theBarrelGeometry_handle;
144  iSetup.get<EcalBarrelGeometryRecord>().get("EcalBarrel", theBarrelGeometry_handle);
145  const CaloSubdetectorGeometry* theBarrelGeometry = theBarrelGeometry_handle.product();
146 
147  if (debug_) {
148  std::cout << " TP analyzer =================> Treating event " << iEvent.id() << " Number of TPs "
149  << tp.product()->size() << std::endl;
150  if (recHits_)
151  std::cout << " Number of EB rechits " << rechit_EB_col.product()->size() << std::endl;
152  }
153  hAllTPperEvt_->Fill(float(tp.product()->size()));
154 
155  //if ( iEvent.id().event() != 648) return;
156 
157  //EcalEBTPGScale ecalScale ;
158  EcalTPGScale ecalScale;
159  ecalScale.setEventSetup(iSetup);
160 
161  // for(unsigned int iDigi = 0; iDigi < ebdigi->size() ; ++iDigi) {
162  // EBDataFrame myFrame((*ebdigi)[iDigi]);
163  // const EBDetId & myId = myFrame.id();
164 
165  int nTP = 0;
166  for (unsigned int i = 0; i < tp.product()->size(); i++) {
167  EcalEBTriggerPrimitiveDigi d = (*(tp.product()))[i];
168  const EBDetId TPid = d.id();
169  // if ( myId != TPid ) continue;
170 
171  /*
172  int index=getIndex(ebdigi,coarser);
173  std::cout << " Same xTal " << myId << " " << TPid << " coarser " << coarser << " index " << index << std::endl;
174  double Et = ecalScale.getTPGInGeV(d.encodedEt(), coarser) ;
175  */
176  //this works if the energy is compressed into 8 bits float Et=d.compressedEt()/2.; // 2ADC counts/GeV
177  float Et = d.encodedEt() / 8.; // 8 ADCcounts/GeV
178  if (Et <= 5)
179  continue;
180  //if ( Et<= 0 ) continue;
181  nTP++;
182 
183  std::cout << " TP digi size " << d.size() << std::endl;
184  for (int iBx = 0; iBx < d.size(); iBx++) {
185  std::cout << " TP samples " << d.sample(iBx) << std::endl;
186  }
187 
188  // EcalTrigTowerDetId coarser=(*eTTmap_).towerOf(myId);
189  // does not work float etaTow = theBarrelGeometry->getGeometry(coarser)->getPosition().theta();
190  // float etaTP = theBarrelGeometry->getGeometry(TPid)->getPosition().eta();
191  // does not work hTPvsTow_eta_->Fill ( etaTow, etaTP );
192  // hTPvsTow_ieta_->Fill ( coarser.ieta(), TPid.ieta() );
193 
194  tpIphi_ = TPid.iphi();
195  tpIeta_ = TPid.ieta();
196  tpgADC_ = d.encodedEt();
197  tpgGeV_ = Et;
198 
199  hTP_iphiVsieta_->Fill(TPid.ieta(), TPid.iphi(), Et);
200  hTP_iphiVsieta_fullrange_->Fill(TPid.ieta(), TPid.iphi(), Et);
201 
202  if (recHits_) {
203  for (unsigned int j = 0; j < rechit_EB_col.product()->size(); j++) {
204  const EBDetId& myid1 = (*rechit_EB_col.product())[j].id();
205  float theta = theBarrelGeometry->getGeometry(myid1)->getPosition().theta();
206  float rhEt = ((*rechit_EB_col.product())[j].energy()) * sin(theta);
207  if (myid1 == TPid) {
208  if (debug_)
209  std::cout << " Analyzer same cristal " << myid1 << " " << TPid << std::endl;
210  // if ( rhEt < 1.5 && Et > 10 ) {
211  // std::cout << " TP analyzer =================> Treating event "<<iEvent.id()<< ", Number of EB rechits "<< rechit_EB_col.product()->size() << " Number of TPs " << tp.product()->size() << std::endl;
212  //std::cout << " TP compressed et " << d.encodedEt() << " Et in GeV " << Et << " RH Et " << rhEt << " Et/rhEt " << Et/rhEt << std::endl;
213  //}
214 
215  //std::cout << " TP out " << d << std::endl;
216 
217  // for (int isam=0;isam< d.size();++isam) {
218  // std::cout << " d[isam].raw() " << d[isam].raw() << std::endl;
219  //}
220 
221  rhIphi_ = myid1.iphi();
222  rhIeta_ = myid1.ieta();
223  hRH_iphiVsieta_->Fill(myid1.ieta(), myid1.iphi(), rhEt);
224  hRH_iphiVsieta_fullrange_->Fill(myid1.ieta(), myid1.iphi(), rhEt);
225 
226  hTPvsRechit_->Fill(rhEt, Et);
227  hTPoverRechit_->Fill(Et / rhEt);
228  hDeltaEt_->Fill((rhEt - Et) / rhEt);
229  if (debug_)
230  std::cout << " TP compressed et " << d.encodedEt() << " Et in GeV " << Et << " RH Et " << rhEt
231  << " Et/rhEt " << Et / rhEt << std::endl;
232  hRechitEt_->Fill(rhEt);
233  hTPEt_->Fill(Et);
234  if (rhEt < 1000000)
235  eRec_ = rhEt;
236  tree_->Fill();
237  }
238 
239  } // end loop of recHits
240  } // if recHits
241 
242  } // end loop over TP collection
243 
244  // } // end loop over digi collection
245 
246  hTPperEvt_->Fill(float(nTP));
247 
248  if (recHits_) {
249  hRatioEt_->Divide(hTPEt_, hRechitEt_);
250  for (unsigned int j = 0; j < rechit_EB_col.product()->size(); j++) {
251  const EBDetId& myid1 = (*rechit_EB_col.product())[j].id();
252  float theta = theBarrelGeometry->getGeometry(myid1)->getPosition().theta();
253  float rhEt = ((*rechit_EB_col.product())[j].energy()) * sin(theta);
254  if (rhEt > 0)
255  hAllRechitEt_->Fill(rhEt);
256  }
257  }
258 }
259 
260 void EcalEBTrigPrimAnalyzer::endJob() {
261  for (unsigned int i = 0; i < ecal_parts_.size(); ++i) {
262  ecal_et_[i]->Write();
263  ecal_tt_[i]->Write();
264  ecal_fgvb_[i]->Write();
265  }
266 
267  hAllTPperEvt_->Write();
268  hTPperEvt_->Write();
269 
270  if (recHits_) {
271  hTPvsTow_ieta_->Write();
272  hTPvsTow_eta_->Write();
273  hTPvsRechit_->Write();
274  hTPoverRechit_->Write();
275  hAllRechitEt_->Write();
276  hRechitEt_->Write();
277  hDeltaEt_->Write();
278  hTPEt_->Write();
279  hRatioEt_->Write();
280  hTP_iphiVsieta_->Write();
281  hRH_iphiVsieta_->Write();
282  hTP_iphiVsieta_fullrange_->Write();
283  hRH_iphiVsieta_fullrange_->Write();
284  }
285 }
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
init
int init
Definition: HydjetWrapper.h:64
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:49
mps_fire.i
i
Definition: mps_fire.py:355
EcalEBTrigPrimAnalyzer::EcalEBTrigPrimAnalyzer
EcalEBTrigPrimAnalyzer(const edm::ParameterSet &)
Definition: EcalEBTrigPrimAnalyzer.cc:35
EcalTPGScale
Definition: EcalTPGScale.h:8
edm::Handle::product
T const * product() const
Definition: Handle.h:70
ESHandle.h
EBDetId
Definition: EBDetId.h:17
edm
HLT enums.
Definition: AlignableModifier.h:19
gather_cfg.cout
cout
Definition: gather_cfg.py:144
edm::SortedCollection::size
size_type size() const
Definition: SortedCollection.h:215
EDAnalyzer.h
EcalTPGScale::setEventSetup
void setEventSetup(const edm::EventSetup &evtSetup)
Definition: EcalTPGScale.cc:16
EcalEBTrigPrimAnalyzer::init
void init(const edm::EventSetup &)
Definition: EcalEBTrigPrimAnalyzer.cc:101
edm::Handle
Definition: AssociativeIterator.h:50
EcalBarrelGeometryRecord
Definition: EcalBarrelGeometryRecord.h:23
EcalRecHitCollections.h
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MakerMacros.h
edm::EventSetup::get
T get() const
Definition: EventSetup.h:73
EcalDigiCollections.h
edm::ESHandle< CaloSubdetectorGeometry >
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
EcalEBTrigPrimAnalyzer::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: EcalEBTrigPrimAnalyzer.cc:108
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
cmsswSequenceInfo.tp
tp
Definition: cmsswSequenceInfo.py:17
CaloGeometryRecord.h
CaloSubdetectorGeometry.h
edm::ParameterSet
Definition: ParameterSet.h:36
Event.h
edm::get
T const & get(Event const &event, InputTag const &tag) noexcept(false)
Definition: Event.h:669
iEvent
int iEvent
Definition: GenABIO.cc:224
CaloSubdetectorGeometry::getGeometry
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
Definition: CaloSubdetectorGeometry.cc:36
edm::EventSetup
Definition: EventSetup.h:57
EcalRecHit.h
overlapproblemtsosanalyzer_cfi.title
title
Definition: overlapproblemtsosanalyzer_cfi.py:7
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
CaloCellGeometry.h
EcalTrigTowerConstituentsMap.h
EcalEBTriggerPrimitiveDigi
Definition: EcalEBTriggerPrimitiveDigi.h:15
std
Definition: JetResolutionObject.h:76
EcalEBTrigPrimAnalyzer.h
CaloGeometry.h
EventSetup.h
CaloSubdetectorGeometry
Definition: CaloSubdetectorGeometry.h:22
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:51
ztail.d
d
Definition: ztail.py:151
ParameterSet.h
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
edm::Event
Definition: Event.h:73
EcalEBTrigPrimAnalyzer::endJob
void endJob() override
Definition: EcalEBTrigPrimAnalyzer.cc:260
edm::InputTag
Definition: InputTag.h:15
EcalEBTrigPrimAnalyzer::~EcalEBTrigPrimAnalyzer
~EcalEBTrigPrimAnalyzer() override
Definition: EcalEBTrigPrimAnalyzer.cc:93
IdealGeometryRecord
Definition: IdealGeometryRecord.h:27
EcalTPGScale.h
Generated for CMSSW Reference Manual by doxygen 1.8.16