CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MinBias.cc
Go to the documentation of this file.
1 // user include files
5 
6 using namespace std;
7 namespace cms
8 {
9 MinBias::MinBias(const edm::ParameterSet& iConfig)
10 {
11  // get names of modules, producing object collections
12  hbheLabel_= iConfig.getParameter<std::string>("hbheInput");
13  hoLabel_=iConfig.getParameter<std::string>("hoInput");
14  hfLabel_=iConfig.getParameter<std::string>("hfInput");
15  allowMissingInputs_=iConfig.getUntrackedParameter<bool>("AllowMissingInputs",false);
16  // get name of output file with histogramms
17  fOutputFileName = iConfig.getUntrackedParameter<std::string>("HistOutFile");
18 
19 }
20 
21 
22 MinBias::~MinBias()
23 {
24 
25  // do anything here that needs to be done at desctruction time
26  // (e.g. close files, deallocate resources etc.)
27 
28 }
29 
31 {
32  hOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
33  myTree = new TTree("RecJet","RecJet Tree");
34  myTree->Branch("mydet", &mydet, "mydet/I");
35  myTree->Branch("mysubd", &mysubd, "mysubd/I");
36  myTree->Branch("depth", &depth, "depth/I");
37  myTree->Branch("ieta", &ieta, "ieta/I");
38  myTree->Branch("iphi", &iphi, "iphi/I");
39  myTree->Branch("eta", &eta, "eta/F");
40  myTree->Branch("phi", &phi, "phi/F");
41  myTree->Branch("mom1", &mom1, "mom1/F");
42  myTree->Branch("mom2", &mom2, "mom2/F");
43  myTree->Branch("mom3", &mom3, "mom3/F");
44  myTree->Branch("mom4", &mom4, "mom4/F");
45 
46  ievent = 0;
47 
48 }
49 
50 void MinBias::endJob()
51 {
52 
53  const std::vector<DetId>& did = geo->getSubdetectorGeometry( DetId::Hcal, 1 )->getValidDetIds() ;
54  int i=0;
55  for(std::vector<DetId>::const_iterator id=did.begin(); id != did.end(); id++)
56  {
57 // if( (*id).det() == DetId::Hcal ) {
58  GlobalPoint pos = geo->getPosition(*id);
59  mydet = ((*id).rawId()>>28)&0xF;
60  mysubd = ((*id).rawId()>>25)&0x7;
61  depth = HcalDetId(*id).depth();
62  ieta = HcalDetId(*id).ieta();
63  iphi = HcalDetId(*id).iphi();
64  phi = pos.phi();
65  eta = pos.eta();
66  if ( theFillDetMap0[*id] > 0. )
67  {
68  mom1 = theFillDetMap1[*id]/theFillDetMap0[*id];
69  mom2 = theFillDetMap2[*id]/theFillDetMap0[*id]-(mom1*mom1);
70  mom3 = theFillDetMap3[*id]/theFillDetMap0[*id]-3.*mom1*theFillDetMap2[*id]/theFillDetMap0[*id]+
71  2.*pow(mom2,3);
72  mom4 = (theFillDetMap4[*id]-4.*mom1*theFillDetMap3[*id]+6.*pow(mom1,2)*theFillDetMap2[*id])/theFillDetMap0[*id]-3.*pow(mom1,4);
73 
74  } else
75  {
76  mom1 = 0.; mom2 = 0.; mom3 = 0.; mom4 = 0.;
77  }
78  std::cout<<" Detector "<<(*id).rawId()<<" mydet "<<mydet<<" "<<mysubd<<" "<<depth<<" "<<
79  HcalDetId(*id).subdet()<<" "<<ieta<<" "<<iphi<<" "<<pos.eta()<<" "<<pos.phi()<<std::endl;
80  std::cout<<" Energy "<<mom1<<" "<<mom2<<std::endl;
81  myTree->Fill();
82  i++;
83 // }
84  }
85  std::cout<<" The number of CaloDet records "<<did.size()<<std::endl;
86  std::cout<<" The number of Hcal records "<<i<<std::endl;
87 
88 
89  std::cout << "===== Start writing user histograms =====" << std::endl;
90  hOutputFile->SetCompressionLevel(2);
91  hOutputFile->cd();
92  myTree->Write();
93  hOutputFile->Close() ;
94  std::cout << "===== End writing user histograms =======" << std::endl;
95 }
96 
97 
98 //
99 // member functions
100 //
101 
102 // ------------ method called to produce the data ------------
103 void
105 {
106 
107  using namespace edm;
108  if(ievent == 0 ){
110  iSetup.get<CaloGeometryRecord>().get(pG);
111  geo = pG.product();
112  std::vector<DetId> did = geo->getValidDetIds();
113 
114  for(std::vector<DetId>::const_iterator id=did.begin(); id != did.end(); id++)
115  {
116  if( (*id).det() == DetId::Hcal ) {
117  theFillDetMap0[*id] = 0.;
118  theFillDetMap1[*id] = 0.;
119  theFillDetMap2[*id] = 0.;
120  theFillDetMap3[*id] = 0.;
121  theFillDetMap4[*id] = 0.;
122  }
123  }
124  }
125 
126 
127 
128  if (!hbheLabel_.empty()) {
130  iEvent.getByLabel(hbheLabel_,hbhe);
131  if (!hbhe.isValid()) {
132  // can't find it!
133  if (!allowMissingInputs_) {
134  *hbhe; // will throw the proper exception
135  }
136  } else {
137  for(HBHERecHitCollection::const_iterator hbheItr = (*hbhe).begin();
138  hbheItr != (*hbhe).end(); ++hbheItr)
139  {
140  DetId id = (hbheItr)->detid();
141  if( (*hbheItr).energy() > 0. ) std::cout<<" Energy = "<<(*hbheItr).energy()<<std::endl;
142  theFillDetMap0[id] = theFillDetMap0[id]+ 1.;
143  theFillDetMap1[id] = theFillDetMap1[id]+(*hbheItr).energy();
144  theFillDetMap2[id] = theFillDetMap2[id]+pow((*hbheItr).energy(),2);
145  theFillDetMap3[id] = theFillDetMap3[id]+pow((*hbheItr).energy(),3);
146  theFillDetMap4[id] = theFillDetMap4[id]+pow((*hbheItr).energy(),4);
147  }
148  }
149  }
150 
151  if (!hoLabel_.empty()) {
153  iEvent.getByLabel(hoLabel_,ho);
154  if (!ho.isValid()) {
155  // can't find it!
156  if (!allowMissingInputs_) {
157  *ho; // will throw the proper exception
158  }
159  } else {
160  for(HORecHitCollection::const_iterator hoItr = (*ho).begin();
161  hoItr != (*ho).end(); ++hoItr)
162  {
163  DetId id = (hoItr)->detid();
164  theFillDetMap0[id] = theFillDetMap0[id]+ 1.;
165  theFillDetMap1[id] = theFillDetMap1[id]+(*hoItr).energy();
166  theFillDetMap2[id] = theFillDetMap2[id]+pow((*hoItr).energy(),2);
167  theFillDetMap3[id] = theFillDetMap3[id]+pow((*hoItr).energy(),3);
168  theFillDetMap4[id] = theFillDetMap4[id]+pow((*hoItr).energy(),4);
169  }
170  }
171  }
172 
173  if (!hfLabel_.empty()) {
175  iEvent.getByLabel(hfLabel_,hf);
176  if (!hf.isValid()) {
177  // can't find it!
178  if (!allowMissingInputs_) {
179  *hf; // will throw the proper exception
180  }
181  } else {
182  for(HFRecHitCollection::const_iterator hfItr = (*hf).begin();
183  hfItr != (*hf).end(); ++hfItr)
184  {
185  DetId id = (hfItr)->detid();
186  theFillDetMap0[id] = theFillDetMap0[id]+ 1.;
187  theFillDetMap1[id] = theFillDetMap1[id]+(*hfItr).energy();
188  theFillDetMap2[id] = theFillDetMap2[id]+pow((*hfItr).energy(),2);
189  theFillDetMap3[id] = theFillDetMap3[id]+pow((*hfItr).energy(),3);
190  theFillDetMap4[id] = theFillDetMap4[id]+pow((*hfItr).energy(),4);
191  }
192  }
193  }
194 
195 }
196 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:32
Geom::Phi< T > phi() const
Definition: PV3DBase.h:68
std::vector< T >::const_iterator const_iterator
T eta() const
void beginJob()
Definition: Breakpoints.cc:15
int depth() const
get the tower depth
Definition: HcalDetId.h:42
int iEvent
Definition: GenABIO.cc:243
int ieta() const
get the cell ieta
Definition: HcalDetId.h:38
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
int iphi() const
get the cell iphi
Definition: HcalDetId.h:40
Definition: DetId.h:20
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
T eta() const
Definition: PV3DBase.h:75
tuple cout
Definition: gather_cfg.py:121
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
Definition: DDAxes.h:10