CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
SimAnalyzerMinbias.cc
Go to the documentation of this file.
1 // system include files
2 #include <memory>
3 #include <string>
4 #include <iostream>
5 #include <fstream>
6 #include <sstream>
7 #include <vector>
8 #include <map>
9 
10 // user include files
21 
26 
32 
34 
35 #include "TFile.h"
36 #include "TTree.h"
37 
38 // class declaration
39 class SimAnalyzerMinbias : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
40 public:
41  explicit SimAnalyzerMinbias(const edm::ParameterSet&);
42  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
43  void analyze(const edm::Event&, const edm::EventSetup&) override;
44  void beginJob() override;
45  void endJob() override;
46  void beginRun(const edm::Run& r, const edm::EventSetup& iSetup) override{};
47  void endRun(const edm::Run& r, const edm::EventSetup& iSetup) override{};
48 
49 private:
50  // ----------member data ---------------------------
51  const double timeCut_;
52  TTree* myTree_;
53 
54  // Root tree members
57  struct myInfo {
59  myInfo() { theMB0 = theMB1 = theMB2 = theMB3 = theMB4 = 0; }
60  };
61  std::map<HcalDetId, myInfo> myMap_;
64 };
65 
66 // constructors and destructor
67 
69  : timeCut_(iConfig.getUntrackedParameter<double>("TimeCut", 500)),
70  tok_evt_(consumes<edm::HepMCProduct>(edm::InputTag("generator"))),
71  tok_hcal_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"))) {
72  usesResource(TFileService::kSharedResource);
73  edm::LogVerbatim("AnalyzerMB") << "Use Time cut of " << timeCut_ << " ns";
74 }
75 
78  desc.addUntracked<double>("TimeCut", 500);
79  descriptions.add("simAnalyzerMinbias", desc);
80 }
81 
84  myTree_ = fs->make<TTree>("SimJet", "SimJet Tree");
85  myTree_->Branch("mydet", &mydet, "mydet/I");
86  myTree_->Branch("mysubd", &mysubd, "mysubd/I");
87  myTree_->Branch("cells", &cells, "cells");
88  myTree_->Branch("depth", &depth, "depth/I");
89  myTree_->Branch("ieta", &ieta, "ieta/I");
90  myTree_->Branch("iphi", &iphi, "iphi/I");
91  myTree_->Branch("mom0_MB", &mom0_MB, "mom0_MB/F");
92  myTree_->Branch("mom1_MB", &mom1_MB, "mom1_MB/F");
93  myTree_->Branch("mom2_MB", &mom2_MB, "mom2_MB/F");
94  myTree_->Branch("mom3_MB", &mom3_MB, "mom3_MB/F");
95  myTree_->Branch("mom4_MB", &mom4_MB, "mom4_MB/F");
96 
97  myMap_.clear();
98 }
99 
100 // EndJob
101 //
103  cells = 0;
104  for (std::map<HcalDetId, myInfo>::const_iterator itr = myMap_.begin(); itr != myMap_.end(); ++itr) {
105  mysubd = itr->first.subdet();
106  depth = itr->first.depth();
107  iphi = itr->first.iphi();
108  ieta = itr->first.ieta();
109  myInfo info = itr->second;
110  if (info.theMB0 > 0) {
111  mom0_MB = info.theMB0;
112  mom1_MB = info.theMB1;
113  mom2_MB = info.theMB2;
114  mom3_MB = info.theMB3;
115  mom4_MB = info.theMB4;
116  cells++;
117 
118  edm::LogVerbatim("SimAnalyzerMinbias")
119  << " Result= " << mysubd << " " << ieta << " " << iphi << " mom0 " << mom0_MB << " mom1 " << mom1_MB
120  << " mom2 " << mom2_MB << " mom3 " << mom3_MB << " mom4 " << mom4_MB;
121  myTree_->Fill();
122  }
123  }
124  edm::LogVerbatim("SimAnalyzerMinbias") << "cells " << cells;
125 }
126 
127 //
128 // member functions
129 //
130 
131 // ------------ method called to produce the data ------------
132 
134  edm::LogVerbatim("SimAnalyzerMinbias") << " Start SimAnalyzerMinbias::analyze " << iEvent.id().run() << ":"
135  << iEvent.id().event();
136 
137  const edm::Handle<edm::HepMCProduct>& evtMC = iEvent.getHandle(tok_evt_);
138  if (!evtMC.isValid()) {
139  edm::LogWarning("SimAnalyzerMinbias") << "no HepMCProduct found";
140  } else {
141  const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
142  edm::LogVerbatim("SimAnalyzerMinbias") << "Event with " << myGenEvent->particles_size() << " particles + "
143  << myGenEvent->vertices_size() << " vertices";
144  }
145 
146  const edm::Handle<edm::PCaloHitContainer>& hcalHits = iEvent.getHandle(tok_hcal_);
147  if (!hcalHits.isValid()) {
148  edm::LogWarning("SimAnalyzerMinbias") << "Error! can't get HcalHits product!";
149  return;
150  }
151 
152  const edm::PCaloHitContainer* HitHcal = hcalHits.product();
153  std::map<HcalDetId, double> hitMap;
154  for (std::vector<PCaloHit>::const_iterator hcalItr = HitHcal->begin(); hcalItr != HitHcal->end(); ++hcalItr) {
155  double time = hcalItr->time();
156  if (time < timeCut_) {
157  double energyhit = hcalItr->energy();
158  HcalDetId hid = HcalDetId(hcalItr->id());
159  std::map<HcalDetId, double>::iterator itr1 = hitMap.find(hid);
160  if (itr1 == hitMap.end()) {
161  hitMap[hid] = 0;
162  itr1 = hitMap.find(hid);
163  }
164  itr1->second += energyhit;
165  }
166  }
167  edm::LogVerbatim("SimAnalyzerMinbias") << "extract information of " << hitMap.size() << " towers from "
168  << HitHcal->size() << " hits";
169 
170  for (std::map<HcalDetId, double>::const_iterator hcalItr = hitMap.begin(); hcalItr != hitMap.end(); ++hcalItr) {
171  HcalDetId hid = hcalItr->first;
172  double energyhit = hcalItr->second;
173  std::map<HcalDetId, myInfo>::iterator itr1 = myMap_.find(hid);
174  if (itr1 == myMap_.end()) {
175  myInfo info;
176  myMap_[hid] = info;
177  itr1 = myMap_.find(hid);
178  }
179  itr1->second.theMB0++;
180  itr1->second.theMB1 += energyhit;
181  itr1->second.theMB2 += (energyhit * energyhit);
182  itr1->second.theMB3 += (energyhit * energyhit * energyhit);
183  itr1->second.theMB4 += (energyhit * energyhit * energyhit * energyhit);
184  edm::LogVerbatim("SimAnalyzerMinbias") << "ID " << hid << " with energy " << energyhit;
185  }
186 }
187 
188 //define this as a plug-in
static const std::string kSharedResource
Definition: TFileService.h:76
Log< level::Info, true > LogVerbatim
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void beginJob() override
std::vector< PCaloHit > PCaloHitContainer
static const TGPicture * info(bool iBackgroundIsBlack)
T const * product() const
Definition: Handle.h:70
void endRun(const edm::Run &r, const edm::EventSetup &iSetup) override
void beginRun(const edm::Run &r, const edm::EventSetup &iSetup) override
int iEvent
Definition: GenABIO.cc:224
void endJob() override
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
const edm::EDGetTokenT< edm::HepMCProduct > tok_evt_
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:37
void analyze(const edm::Event &, const edm::EventSetup &) override
void add(std::string const &label, ParameterSetDescription const &psetDescription)
const edm::EDGetTokenT< edm::PCaloHitContainer > tok_hcal_
bool isValid() const
Definition: HandleBase.h:70
HLT enums.
std::map< HcalDetId, myInfo > myMap_
SimAnalyzerMinbias(const edm::ParameterSet &)
Log< level::Warning, false > LogWarning
Definition: Run.h:45