CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CastorTowerProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Castor
4 // Class: CastorTowerProducer
5 //
6 
13 //
14 // Original Author: Hans Van Haevermaet, Benoit Roland
15 // Created: Wed Jul 9 14:00:40 CEST 2008
16 // $Id: CastorTowerProducer.cc,v 1.11 2012/07/18 00:31:00 slava77 Exp $
17 //
18 //
19 
20 
21 // system include
22 #include <memory>
23 #include <vector>
24 #include <iostream>
25 #include <TMath.h>
26 #include <TRandom3.h>
27 
28 // user include
37 
38 // Castor Object include
42 
43 // Channel quality
47 
48 //
49 // class declaration
50 //
51 
53  public:
54  explicit CastorTowerProducer(const edm::ParameterSet&);
56 
57  private:
58  virtual void beginJob() ;
59  virtual void produce(edm::Event&, const edm::EventSetup&);
60  virtual void endJob() ;
61  virtual void ComputeTowerVariable(const edm::RefVector<edm::SortedCollection<CastorRecHit> >& usedRecHits, double& Ehot, double& depth);
62 
63  // member data
65  typedef ROOT::Math::RhoEtaPhiPoint TowerPoint;
66  typedef ROOT::Math::RhoZPhiPoint CellPoint;
68  typedef std::vector<reco::CastorTower> CastorTowerCollection;
70  std::string input_;
71  double towercut_;
72  double mintime_;
73  double maxtime_;
74 };
75 
76 //
77 // constants, enums and typedefs
78 //
79 
80 const double MYR2D = 180/M_PI;
81 
82 //
83 // static data member definitions
84 //
85 
86 //
87 // constructor and destructor
88 //
89 
91  input_(iConfig.getUntrackedParameter<std::string>("inputprocess","castorreco")),
92  towercut_(iConfig.getUntrackedParameter<double>("towercut",0.65)),
93  mintime_(iConfig.getUntrackedParameter<double>("mintime",-999)),
94  maxtime_(iConfig.getUntrackedParameter<double>("maxtime",999))
95 {
96  //register your products
97  produces<CastorTowerCollection>();
98  //now do what ever other initialization is needed
99 }
100 
101 
103 {
104  // do anything here that needs to be done at desctruction time
105  // (e.g. close files, deallocate resources etc.)
106 }
107 
108 
109 //
110 // member functions
111 //
112 
113 // ------------ method called to produce the data ------------
115 
116  using namespace edm;
117  using namespace reco;
118  using namespace TMath;
119 
120  // Produce CastorTowers from CastorCells
121 
123  iEvent.getByLabel(input_,InputRecHits);
124 
125  std::auto_ptr<CastorTowerCollection> OutputTowers (new CastorTowerCollection);
126 
127  // get and check input size
128  int nRecHits = InputRecHits->size();
129 
130  LogDebug("CastorTowerProducer")
131  <<"2. entering CastorTowerProducer"<<std::endl;
132 
133  if (nRecHits==0)
134  LogDebug("CastorTowerProducer") <<"Warning: You are trying to run the Tower algorithm with 0 input rechits.";
135 
136  // declare castor array
137  // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position
138 
139  double poscastortowerarray[4][16];
140  double negcastortowerarray[4][16];
141 
142  CastorRecHitRefVector poscastorusedrechits[16];
143  CastorRecHitRefVector negcastorusedrechits[16];
144 
145  // set phi values and everything else to zero
146  for (int j = 0; j < 16; j++) {
147  poscastortowerarray[3][j] = -2.94524 + j*0.3927;
148  poscastortowerarray[0][j] = 0.;
149  poscastortowerarray[1][j] = 0.;
150  poscastortowerarray[2][j] = 0.;
151 
152  negcastortowerarray[3][j] = -2.94524 + j*0.3927;
153  negcastortowerarray[0][j] = 0.;
154  negcastortowerarray[1][j] = 0.;
155  negcastortowerarray[2][j] = 0.;
156  }
157 
158  // retrieve the channel quality lists from database
160  iSetup.get<CastorChannelQualityRcd>().get(p);
161  std::vector<DetId> channels = p->getAllChannels();
162 
163  // loop over rechits to build castortowerarray[4][16] and castorusedrechits[16]
164  for (unsigned int i = 0; i < InputRecHits->size(); i++) {
165 
167 
168  HcalCastorDetId id = rechit_p->id();
169  DetId genericID=(DetId)id;
170 
171  // first check if the rechit is in the BAD channel list
172  bool bad = false;
173  for (std::vector<DetId>::iterator channel = channels.begin();channel != channels.end();channel++) {
174  if (channel->rawId() == genericID.rawId()) {
175  // if the rechit is found in the list, set it bad
176  bad = true; break;
177  }
178  }
179  // if bad, continue the loop to the next rechit
180  if (bad) continue;
181 
182  double Erechit = rechit_p->energy();
183  int module = id.module();
184  int sector = id.sector();
185  double zrechit = 0;
186  if (module < 3) zrechit = -14390 - 24.75 - 49.5*(module-1);
187  if (module > 2) zrechit = -14390 - 99 - 49.5 - 99*(module-3);
188  double phirechit = -100;
189  if (sector < 9) phirechit = 0.19635 + (sector-1)*0.3927;
190  if (sector > 8) phirechit = -2.94524 + (sector - 9)*0.3927;
191 
192  // add time conditions for the rechit
193  if (rechit_p->time() > mintime_ && rechit_p->time() < maxtime_) {
194 
195  // loop over the 16 towers possibilities
196  for ( int j=0;j<16;j++) {
197 
198  // phi matching condition
199  if (TMath::Abs(phirechit - poscastortowerarray[3][j]) < 0.0001) {
200 
201  // condition over rechit z value
202  if (zrechit > 0.) {
203  poscastortowerarray[0][j]+=Erechit;
204  if (module < 3) {poscastortowerarray[1][j]+=Erechit;} else {poscastortowerarray[2][j]+=Erechit;}
205  poscastorusedrechits[j].push_back(rechit_p);
206  } else {
207  negcastortowerarray[0][j]+=Erechit;
208  if (module < 3) {negcastortowerarray[1][j]+=Erechit;} else {negcastortowerarray[2][j]+=Erechit;}
209  negcastorusedrechits[j].push_back(rechit_p);
210  } // end condition over rechit z value
211  } // end phi matching condition
212  } // end loop over the 16 towers possibilities
213  } // end time conditions
214 
215  } // end loop over rechits to build castortowerarray[4][16] and castorusedrechits[16]
216 
217  // make towers of the arrays
218 
219  double fem, Ehot, depth;
220  double rhoTower = 88.5;
221 
222  // loop over the 16 towers possibilities
223  for (int k=0;k<16;k++) {
224 
225  fem = 0;
226  Ehot = 0;
227  depth = 0;
228 
229  // select the positive towers with E > sqrt(Nusedrechits)*Ecut
230  if (poscastortowerarray[0][k] > sqrt(poscastorusedrechits[k].size())*towercut_) {
231 
232  fem = poscastortowerarray[1][k]/poscastortowerarray[0][k];
233  CastorRecHitRefVector usedRecHits = poscastorusedrechits[k];
234  ComputeTowerVariable(usedRecHits,Ehot,depth);
235 
236  LogDebug("CastorTowerProducer")
237  <<"tower "<<k+1<<": fem = "<<fem<<" ,depth = "<<depth<<" ,Ehot = "<<Ehot<<std::endl;
238 
239  TowerPoint temptowerposition(rhoTower,5.9,poscastortowerarray[3][k]);
240  Point towerposition(temptowerposition);
241 
242  CastorTower newtower(poscastortowerarray[0][k],towerposition,poscastortowerarray[1][k],poscastortowerarray[2][k],fem,depth,Ehot,
243  poscastorusedrechits[k]);
244  OutputTowers->push_back(newtower);
245  } // end select the positive towers with E > Ecut
246 
247  // select the negative towers with E > sqrt(Nusedrechits)*Ecut
248  if (negcastortowerarray[0][k] > sqrt(negcastorusedrechits[k].size())*towercut_) {
249 
250  fem = negcastortowerarray[1][k]/negcastortowerarray[0][k];
251  CastorRecHitRefVector usedRecHits = negcastorusedrechits[k];
252  ComputeTowerVariable(usedRecHits,Ehot,depth);
253 
254  LogDebug("CastorTowerProducer")
255  <<"tower "<<k+1 << " energy = " << negcastortowerarray[0][k] << "EM = " << negcastortowerarray[1][k] << "HAD = " << negcastortowerarray[2][k] << "phi = " << negcastortowerarray[3][k] << ": fem = "<<fem<<" ,depth = "<<depth<<" ,Ehot = "<<Ehot<<std::endl;
256 
257  TowerPoint temptowerposition(rhoTower,-5.9,negcastortowerarray[3][k]);
258  Point towerposition(temptowerposition);
259 
260  CastorTower newtower(negcastortowerarray[0][k],towerposition,negcastortowerarray[1][k],negcastortowerarray[2][k],fem,depth,Ehot,
261  negcastorusedrechits[k]);
262  OutputTowers->push_back(newtower);
263  } // end select the negative towers with E > Ecut
264 
265  } // end loop over the 16 towers possibilities
266 
267  iEvent.put(OutputTowers);
268 }
269 
270 
271 // ------------ method called once each job just before starting event loop ------------
273  LogDebug("CastorTowerProducer")
274  <<"Starting CastorTowerProducer";
275 }
276 
277 // ------------ method called once each job just after ending the event loop ------------
279  LogDebug("CastorTowerProducer")
280  <<"Ending CastorTowerProducer";
281 }
282 
284 
285  using namespace reco;
286 
287  double Etot = 0;
288 
289  // loop over the cells used in the tower k
290  for (CastorRecHitRefVector::iterator it = usedRecHits.begin(); it != usedRecHits.end(); it++) {
291  edm::Ref<CastorRecHitCollection> rechit_p = *it;
292 
293  double Erechit = rechit_p->energy();
294  HcalCastorDetId id = rechit_p->id();
295  int module = id.module();
296  double zrechit = 0;
297  if (module < 3) zrechit = -14390 - 24.75 - 49.5*(module-1);
298  if (module > 2) zrechit = -14390 - 99 - 49.5 - 99*(module-3);
299 
300  if(Erechit > Ehot) Ehot = Erechit;
301  depth+=Erechit*zrechit;
302  Etot+=Erechit;
303  }
304 
305  depth/=Etot;
306  Ehot/=Etot;
307 }
308 
309 //define this as a plug-in
#define LogDebug(id)
int i
Definition: DBlmapReader.cc:9
edm::SortedCollection< CastorRecHit > CastorRecHitCollection
ROOT::Math::RhoZPhiPoint CellPoint
module()
Definition: vlib.cc:994
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< reco::CastorTower > CastorTowerCollection
ROOT::Math::RhoEtaPhiPoint TowerPoint
int bad(Items const &cont)
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
int iEvent
Definition: GenABIO.cc:243
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:85
T sqrt(T t)
Definition: SSEVec.h:46
virtual void produce(edm::Event &, const edm::EventSetup &)
edm::RefVector< CastorRecHitCollection > CastorRecHitRefVector
int j
Definition: DBlmapReader.cc:9
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
int k[5][pyjets_maxn]
Definition: DetId.h:20
#define M_PI
Definition: BFit3D.cc:3
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double > > XYZPointD
point in space with cartesian internal representation
Definition: Point3D.h:9
const double MYR2D
CastorTowerProducer(const edm::ParameterSet &)
virtual void ComputeTowerVariable(const edm::RefVector< edm::SortedCollection< CastorRecHit > > &usedRecHits, double &Ehot, double &depth)
const T & get() const
Definition: EventSetup.h:55
ProductID id() const
Accessor for product ID.
Definition: Ref.h:256
Definition: vlib.h:209
tuple size
Write out results.