CMS 3D CMS Logo

CastorClusterProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Castor
4 // Class: CastorClusterProducer
5 //
12 //
13 // Original Author: Hans Van Haevermaet, Benoit Roland
14 // Created: Wed Jul 9 14:00:40 CEST 2008
15 //
16 //
17 
18 // system include
19 #include <memory>
20 #include <vector>
21 #include <iostream>
22 #include <TMath.h>
23 #include <TRandom3.h>
24 
25 // user include
32 
33 // Castor Object include
43 
44 //
45 // class decleration
46 //
47 
49 public:
51  ~CastorClusterProducer() override;
52 
53 private:
54  void beginJob() override;
55  void produce(edm::Event&, const edm::EventSetup&) override;
56  double phiangle(double testphi);
57  void endJob() override;
58 
59  // ----------member data ---------------------------
61  typedef ROOT::Math::RhoEtaPhiPoint TowerPoint;
62  typedef ROOT::Math::RhoZPhiPoint CellPoint;
63  typedef std::vector<reco::CastorTower> CastorTowerCollection;
64  typedef std::vector<reco::CastorCluster> CastorClusterCollection;
70 };
71 
72 //
73 // constants, enums and typedefs
74 //
75 
76 //
77 // static data member definitions
78 //
79 
80 //
81 // constructor and destructor
82 //
83 
85  : input_(iConfig.getUntrackedParameter<std::string>("inputtowers", "")),
86  basicjets_(iConfig.getUntrackedParameter<std::string>("basicjets", "")),
87  clusteralgo_(iConfig.getUntrackedParameter<bool>("ClusterAlgo", false)) {
88  // register for data access
89  tok_input_ = consumes<CastorTowerCollection>(edm::InputTag(input_));
90  tok_jets_ = consumes<reco::BasicJetCollection>(edm::InputTag(basicjets_));
91  tok_tower_ = consumes<CastorTowerCollection>(edm::InputTag("CastorTowerReco"));
92  // register your products
93  produces<CastorClusterCollection>();
94  // now do what ever other initialization is needed
95 }
96 
98  // do anything here that needs to be done at desctruction time
99  // (e.g. close files, deallocate resources etc.)
100 }
101 
102 //
103 // member functions
104 //
105 
106 // ------------ method called to produce the data ------------
108  using namespace edm;
109  using namespace reco;
110  using namespace TMath;
111 
112  LogDebug("CastorClusterProducer") << "3. entering CastorClusterProducer";
113 
114  if (!input_.empty()) {
115  // Produce CastorClusters from CastorTowers
116 
118  iEvent.getByToken(tok_input_, InputTowers);
119 
120  auto OutputClustersfromClusterAlgo = std::make_unique<CastorClusterCollection>();
121 
122  // get and check input size
123  int nTowers = InputTowers->size();
124 
125  if (nTowers == 0)
126  LogDebug("CastorClusterProducer") << "Warning: You are trying to run the Cluster algorithm with 0 input towers.";
127 
128  CastorTowerRefVector posInputTowers, negInputTowers;
129 
130  for (size_t i = 0; i < InputTowers->size(); ++i) {
131  reco::CastorTowerRef tower_p = reco::CastorTowerRef(InputTowers, i);
132  if (tower_p->eta() > 0.)
133  posInputTowers.push_back(tower_p);
134  if (tower_p->eta() < 0.)
135  negInputTowers.push_back(tower_p);
136  }
137 
138  // build cluster from ClusterAlgo
139  if (clusteralgo_ == true) {
140  // code
141  iEvent.put(std::move(OutputClustersfromClusterAlgo));
142  }
143  }
144 
145  if (!basicjets_.empty()) {
146  Handle<BasicJetCollection> bjCollection;
147  iEvent.getByToken(tok_jets_, bjCollection);
148 
149  Handle<CastorTowerCollection> ctCollection;
150  iEvent.getByToken(tok_tower_, ctCollection);
151 
152  auto OutputClustersfromBasicJets = std::make_unique<CastorClusterCollection>();
153 
154  if (bjCollection->empty())
155  LogDebug("CastorClusterProducer")
156  << "Warning: You are trying to run the Cluster algorithm with 0 input basicjets.";
157 
158  for (unsigned i = 0; i < bjCollection->size(); i++) {
159  const BasicJet* bj = &(*bjCollection)[i];
160 
161  double energy = bj->energy();
162  TowerPoint temp(88.5, bj->eta(), bj->phi());
164  double emEnergy = 0.;
165  double hadEnergy = 0.;
166  double width = 0.;
167  double depth = 0.;
168  double fhot = 0.;
169  double sigmaz = 0.;
170  CastorTowerRefVector usedTowers;
171  double zmean = 0.;
172  double z2mean = 0.;
173 
174  std::vector<CandidatePtr> ccp = bj->getJetConstituents();
175  std::vector<CandidatePtr>::const_iterator itParticle;
176  for (itParticle = ccp.begin(); itParticle != ccp.end(); ++itParticle) {
177  const CastorTower* castorcand = dynamic_cast<const CastorTower*>(itParticle->get());
178  //cout << " castortowercandidate reference energy = " << castorcand->castorTower()->energy() << endl;
179  //cout << " castortowercandidate reference eta = " << castorcand->castorTower()->eta() << endl;
180  //cout << " castortowercandidate reference phi = " << castorcand->castorTower()->phi() << endl;
181  //cout << " castortowercandidate reference depth = " << castorcand->castorTower()->depth() << endl;
182 
183  //CastorTowerCollection *ctc = new CastorTowerCollection();
184  //ctc->push_back(*castorcand);
185  //CastorTowerRef towerref = CastorTowerRef(ctc,0);
186 
187  size_t thisone = 0;
188  for (size_t l = 0; l < ctCollection->size(); l++) {
189  const CastorTower ct = (*ctCollection)[l];
190  if (std::abs(ct.phi() - castorcand->phi()) < 0.0001) {
191  thisone = l;
192  }
193  }
194 
195  CastorTowerRef towerref(ctCollection, thisone);
196  usedTowers.push_back(towerref);
197  emEnergy += castorcand->emEnergy();
198  hadEnergy += castorcand->hadEnergy();
199  depth += castorcand->depth() * castorcand->energy();
200  width += pow(phiangle(castorcand->phi() - bj->phi()), 2) * castorcand->energy();
201  fhot += castorcand->fhot() * castorcand->energy();
202 
203  // loop over rechits
204  for (edm::RefVector<edm::SortedCollection<CastorRecHit> >::iterator it = castorcand->rechitsBegin();
205  it != castorcand->rechitsEnd();
206  it++) {
208  double Erechit = rechit_p->energy();
209  HcalCastorDetId id = rechit_p->id();
210  int module = id.module();
211  double zrechit = 0;
212  if (module < 3)
213  zrechit = -14390 - 24.75 - 49.5 * (module - 1);
214  if (module > 2)
215  zrechit = -14390 - 99 - 49.5 - 99 * (module - 3);
216  zmean += Erechit * zrechit;
217  z2mean += Erechit * zrechit * zrechit;
218  } // end loop over rechits
219  }
220  //cout << "" << endl;
221 
222  depth = depth / energy;
223  width = sqrt(width / energy);
224  fhot = fhot / energy;
225 
226  zmean = zmean / energy;
227  z2mean = z2mean / energy;
228  double sigmaz2 = z2mean - zmean * zmean;
229  if (sigmaz2 > 0)
230  sigmaz = sqrt(sigmaz2);
231 
232  CastorCluster cc(
233  energy, position, emEnergy, hadEnergy, emEnergy / energy, width, depth, fhot, sigmaz, usedTowers);
234  OutputClustersfromBasicJets->push_back(cc);
235  }
236 
237  iEvent.put(std::move(OutputClustersfromBasicJets));
238  }
239 }
240 
241 // help function to calculate phi within [-pi,+pi]
242 double CastorClusterProducer::phiangle(double testphi) {
243  double phi = testphi;
244  while (phi > M_PI)
245  phi -= (2 * M_PI);
246  while (phi < -M_PI)
247  phi += (2 * M_PI);
248  return phi;
249 }
250 
251 // ------------ method called once each job just before starting event loop ------------
252 void CastorClusterProducer::beginJob() { LogDebug("CastorClusterProducer") << "Starting CastorClusterProducer"; }
253 
254 // ------------ method called once each job just after ending the event loop ------------
255 void CastorClusterProducer::endJob() { LogDebug("CastorClusterProducer") << "Ending CastorClusterProducer"; }
256 
257 //define this as a plug-in
edm::EDGetTokenT< CastorTowerCollection > tok_input_
ProductID id() const
Accessor for product ID.
Definition: Ref.h:244
edm::EDGetTokenT< CastorTowerCollection > tok_tower_
CastorRecHitRefs::iterator rechitsBegin() const
fist iterator over CastorRecHit constituents
Definition: CastorTower.h:78
double depth() const
tower depth in z
Definition: CastorTower.h:69
double phiangle(double testphi)
edm::Ref< CastorTowerCollection > CastorTowerRef
Definition: CastorTower.h:140
Jets made from CaloTowers.
Definition: BasicJet.h:19
CastorClusterProducer(const edm::ParameterSet &)
virtual Constituents getJetConstituents() const
list of constituents
double emEnergy() const
tower em energy
Definition: CastorTower.h:60
double hadEnergy() const
tower had energy
Definition: CastorTower.h:63
int iEvent
Definition: GenABIO.cc:224
std::vector< reco::CastorCluster > CastorClusterCollection
T sqrt(T t)
Definition: SSEVec.h:19
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< reco::CastorTower > CastorTowerCollection
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::EDGetTokenT< reco::BasicJetCollection > tok_jets_
ROOT::Math::RhoEtaPhiPoint TowerPoint
double fhot() const
tower hotcell/tot ratio
Definition: CastorTower.h:72
void produce(edm::Event &, const edm::EventSetup &) override
#define M_PI
ROOT::Math::RhoZPhiPoint CellPoint
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double > > XYZPointD
point in space with cartesian internal representation
Definition: Point3D.h:8
fixed size matrix
HLT enums.
CastorRecHitRefs::iterator rechitsEnd() const
last iterator over CastorRecHit constituents
Definition: CastorTower.h:81
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
static int position[264][3]
Definition: ReadPGInfo.cc:289
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:67
double phi() const final
momentum azimuthal angle
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
def move(src, dest)
Definition: eostools.py:511
#define LogDebug(id)
double energy() const final
energy
double eta() const final
momentum pseudorapidity