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.7 2010/07/03 19:24:25 hvanhaev 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
35 
36 // Castor Object include
40 
41 //
42 // class declaration
43 //
44 
46  public:
47  explicit CastorTowerProducer(const edm::ParameterSet&);
49 
50  private:
51  virtual void beginJob() ;
52  virtual void produce(edm::Event&, const edm::EventSetup&);
53  virtual void endJob() ;
54  virtual void ComputeTowerVariable(const edm::RefVector<edm::SortedCollection<CastorRecHit> >& usedRecHits, double& Ehot, double& depth);
55 
56  // member data
58  typedef ROOT::Math::RhoEtaPhiPoint TowerPoint;
59  typedef ROOT::Math::RhoZPhiPoint CellPoint;
61  typedef std::vector<reco::CastorTower> CastorTowerCollection;
63  std::string input_;
64  double towercut_;
65  double mintime_;
66  double maxtime_;
67 };
68 
69 //
70 // constants, enums and typedefs
71 //
72 
73 const double MYR2D = 180/M_PI;
74 
75 //
76 // static data member definitions
77 //
78 
79 //
80 // constructor and destructor
81 //
82 
84  input_(iConfig.getUntrackedParameter<std::string>("inputprocess","castorreco")),
85  towercut_(iConfig.getUntrackedParameter<double>("towercut",1.)),
86  mintime_(iConfig.getUntrackedParameter<double>("mintime",-999)),
87  maxtime_(iConfig.getUntrackedParameter<double>("maxtime",999))
88 {
89  //register your products
90  produces<CastorTowerCollection>();
91  //now do what ever other initialization is needed
92 }
93 
94 
96 {
97  // do anything here that needs to be done at desctruction time
98  // (e.g. close files, deallocate resources etc.)
99 }
100 
101 
102 //
103 // member functions
104 //
105 
106 // ------------ method called to produce the data ------------
108 
109  using namespace edm;
110  using namespace reco;
111  using namespace TMath;
112 
113  // Produce CastorTowers from CastorCells
114 
116  iEvent.getByLabel(input_,InputRecHits);
117 
118  std::auto_ptr<CastorTowerCollection> OutputTowers (new CastorTowerCollection);
119 
120  // get and check input size
121  int nRecHits = InputRecHits->size();
122 
123  LogDebug("CastorTowerProducer")
124  <<"2. entering CastorTowerProducer"<<std::endl;
125 
126  if (nRecHits==0)
127  LogDebug("CastorTowerProducer") <<"Warning: You are trying to run the Tower algorithm with 0 input rechits.";
128 
129  // declare castor array
130  // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position
131 
132  double poscastortowerarray[4][16];
133  double negcastortowerarray[4][16];
134 
135  CastorRecHitRefVector poscastorusedrechits[16];
136  CastorRecHitRefVector negcastorusedrechits[16];
137 
138  // set phi values and everything else to zero
139  for (int j = 0; j < 16; j++) {
140  poscastortowerarray[3][j] = -2.94524 + j*0.3927;
141  poscastortowerarray[0][j] = 0.;
142  poscastortowerarray[1][j] = 0.;
143  poscastortowerarray[2][j] = 0.;
144 
145  negcastortowerarray[3][j] = -2.94524 + j*0.3927;
146  negcastortowerarray[0][j] = 0.;
147  negcastortowerarray[1][j] = 0.;
148  negcastortowerarray[2][j] = 0.;
149  }
150 
151  // loop over rechits to build castortowerarray[4][16] and castorusedrechits[16]
152  for (unsigned int i = 0; i < InputRecHits->size(); i++) {
153 
155 
156  double Erechit = rechit_p->energy();
157  HcalCastorDetId id = rechit_p->id();
158  int module = id.module();
159  int sector = id.sector();
160  double zrechit = 0;
161  if (module < 3) zrechit = -14390 - 24.75 - 49.5*(module-1);
162  if (module > 2) zrechit = -14390 - 99 - 49.5 - 99*(module-3);
163  double phirechit = -100;
164  if (sector < 9) phirechit = 0.19635 + (sector-1)*0.3927;
165  if (sector > 8) phirechit = -2.94524 + (sector - 9)*0.3927;
166 
167  // add time conditions for the rechit
168  if (rechit_p->time() > mintime_ && rechit_p->time() < maxtime_) {
169 
170  // loop over the 16 towers possibilities
171  for ( int j=0;j<16;j++) {
172 
173  // phi matching condition
174  if (TMath::Abs(phirechit - poscastortowerarray[3][j]) < 0.0001) {
175 
176  // condition over rechit z value
177  if (zrechit > 0.) {
178  poscastortowerarray[0][j]+=Erechit;
179  if (module < 3) {poscastortowerarray[1][j]+=Erechit;} else {poscastortowerarray[2][j]+=Erechit;}
180  poscastorusedrechits[j].push_back(rechit_p);
181  } else {
182  negcastortowerarray[0][j]+=Erechit;
183  if (module < 3) {negcastortowerarray[1][j]+=Erechit;} else {negcastortowerarray[2][j]+=Erechit;}
184  negcastorusedrechits[j].push_back(rechit_p);
185  } // end condition over rechit z value
186  } // end phi matching condition
187  } // end loop over the 16 towers possibilities
188  } // end time conditions
189 
190  } // end loop over rechits to build castortowerarray[4][16] and castorusedrechits[16]
191 
192  // make towers of the arrays
193 
194  double fem, Ehot, depth;
195  double rhoTower = 88.5;
196 
197  // loop over the 16 towers possibilities
198  for (int k=0;k<16;k++) {
199 
200  fem = 0;
201  Ehot = 0;
202  depth = 0;
203 
204  // select the positive towers with E > Ecut
205  if (poscastortowerarray[0][k] > towercut_) {
206 
207  fem = poscastortowerarray[1][k]/poscastortowerarray[0][k];
208  CastorRecHitRefVector usedRecHits = poscastorusedrechits[k];
209  ComputeTowerVariable(usedRecHits,Ehot,depth);
210 
211  LogDebug("CastorTowerProducer")
212  <<"tower "<<k+1<<": fem = "<<fem<<" ,depth = "<<depth<<" ,Ehot = "<<Ehot<<std::endl;
213 
214  TowerPoint temptowerposition(rhoTower,5.9,poscastortowerarray[3][k]);
215  Point towerposition(temptowerposition);
216 
217  CastorTower newtower(poscastortowerarray[0][k],towerposition,poscastortowerarray[1][k],poscastortowerarray[2][k],fem,depth,Ehot,
218  poscastorusedrechits[k]);
219  OutputTowers->push_back(newtower);
220  } // end select the positive towers with E > Ecut
221 
222  // select the negative towers with E > Ecut
223  if (negcastortowerarray[0][k] > towercut_) {
224 
225  fem = negcastortowerarray[1][k]/negcastortowerarray[0][k];
226  CastorRecHitRefVector usedRecHits = negcastorusedrechits[k];
227  ComputeTowerVariable(usedRecHits,Ehot,depth);
228 
229  LogDebug("CastorTowerProducer")
230  <<"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;
231 
232  TowerPoint temptowerposition(rhoTower,-5.9,negcastortowerarray[3][k]);
233  Point towerposition(temptowerposition);
234 
235  CastorTower newtower(negcastortowerarray[0][k],towerposition,negcastortowerarray[1][k],negcastortowerarray[2][k],fem,depth,Ehot,
236  negcastorusedrechits[k]);
237  OutputTowers->push_back(newtower);
238  } // end select the negative towers with E > Ecut
239 
240  } // end loop over the 16 towers possibilities
241 
242  iEvent.put(OutputTowers);
243 }
244 
245 
246 // ------------ method called once each job just before starting event loop ------------
248  LogDebug("CastorTowerProducer")
249  <<"Starting CastorTowerProducer";
250 }
251 
252 // ------------ method called once each job just after ending the event loop ------------
254  LogDebug("CastorTowerProducer")
255  <<"Ending CastorTowerProducer";
256 }
257 
259 
260  using namespace reco;
261 
262  double Etot = 0;
263 
264  // loop over the cells used in the tower k
265  for (CastorRecHitRefVector::iterator it = usedRecHits.begin(); it != usedRecHits.end(); it++) {
266  edm::Ref<CastorRecHitCollection> rechit_p = *it;
267 
268  double Erechit = rechit_p->energy();
269  HcalCastorDetId id = rechit_p->id();
270  int module = id.module();
271  double zrechit = 0;
272  if (module < 3) zrechit = -14390 - 24.75 - 49.5*(module-1);
273  if (module > 2) zrechit = -14390 - 99 - 49.5 - 99*(module-3);
274 
275  if(Erechit > Ehot) Ehot = Erechit;
276  depth+=Erechit*zrechit;
277  Etot+=Erechit;
278  }
279 
280  depth/=Etot;
281  Ehot/=Etot;
282 }
283 
284 //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 iEvent
Definition: GenABIO.cc:243
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
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:355
int k[5][pyjets_maxn]
#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)
ProductID id() const
Accessor for product ID.
Definition: Ref.h:255
Definition: vlib.h:209