113 using namespace reco;
114 using namespace TMath;
121 auto OutputTowers = std::make_unique<CastorTowerCollection>();
124 int nRecHits = InputRecHits->
size();
127 <<
"2. entering CastorTowerProducer"<<std::endl;
130 LogDebug(
"CastorTowerProducer") <<
"Warning: You are trying to run the Tower algorithm with 0 input rechits.";
135 double poscastortowerarray[4][16];
136 double negcastortowerarray[4][16];
142 for (
int j = 0; j < 16; j++) {
143 poscastortowerarray[3][j] = -2.94524 + j*0.3927;
144 poscastortowerarray[0][j] = 0.;
145 poscastortowerarray[1][j] = 0.;
146 poscastortowerarray[2][j] = 0.;
148 negcastortowerarray[3][j] = -2.94524 + j*0.3927;
149 negcastortowerarray[0][j] = 0.;
150 negcastortowerarray[1][j] = 0.;
151 negcastortowerarray[2][j] = 0.;
160 for (
unsigned int i = 0;
i < InputRecHits->
size();
i++) {
169 for (std::vector<DetId>::iterator channel = channels.begin();channel != channels.end();channel++) {
170 if (channel->rawId() == genericID.
rawId()) {
178 double Erechit = rechit_p->energy();
180 int sector =
id.sector();
182 if (module < 3) zrechit = -14390 - 24.75 - 49.5*(module-1);
183 if (module > 2) zrechit = -14390 - 99 - 49.5 - 99*(module-3);
184 double phirechit = -100;
185 if (sector < 9) phirechit = 0.19635 + (sector-1)*0.3927;
186 if (sector > 8) phirechit = -2.94524 + (sector - 9)*0.3927;
192 for (
int j=0;j<16;j++) {
195 if (
TMath::Abs(phirechit - poscastortowerarray[3][j]) < 0.0001) {
199 poscastortowerarray[0][j]+=Erechit;
200 if (module < 3) {poscastortowerarray[1][j]+=Erechit;}
else {poscastortowerarray[2][j]+=Erechit;}
201 poscastorusedrechits[j].push_back(rechit_p);
203 negcastortowerarray[0][j]+=Erechit;
204 if (module < 3) {negcastortowerarray[1][j]+=Erechit;}
else {negcastortowerarray[2][j]+=Erechit;}
205 negcastorusedrechits[j].push_back(rechit_p);
215 double fem, Ehot,
depth;
216 double rhoTower = 88.5;
219 for (
int k=0;
k<16;
k++) {
228 fem = poscastortowerarray[1][
k]/poscastortowerarray[0][
k];
233 <<
"tower "<<
k+1<<
": fem = "<<fem<<
" ,depth = "<<depth<<
" ,Ehot = "<<Ehot<<std::endl;
235 TowerPoint temptowerposition(rhoTower,5.9,poscastortowerarray[3][
k]);
236 Point towerposition(temptowerposition);
238 CastorTower newtower(poscastortowerarray[0][
k],towerposition,poscastortowerarray[1][k],poscastortowerarray[2][k],fem,depth,Ehot,
239 poscastorusedrechits[k]);
240 OutputTowers->push_back(newtower);
246 fem = negcastortowerarray[1][
k]/negcastortowerarray[0][
k];
251 <<
"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;
253 TowerPoint temptowerposition(rhoTower,-5.9,negcastortowerarray[3][
k]);
254 Point towerposition(temptowerposition);
256 CastorTower newtower(negcastortowerarray[0][
k],towerposition,negcastortowerarray[1][k],negcastortowerarray[2][k],fem,depth,Ehot,
257 negcastorusedrechits[k]);
258 OutputTowers->push_back(newtower);
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
constexpr uint32_t rawId() const
get the raw id
std::vector< DetId > getAllChannels() const
ROOT::Math::RhoEtaPhiPoint TowerPoint
ProductID id() const
Accessor for product ID.
edm::RefVector< CastorRecHitCollection > CastorRecHitRefVector
edm::EDGetTokenT< CastorRecHitCollection > tok_input_
virtual void ComputeTowerVariable(const edm::RefVector< edm::SortedCollection< CastorRecHit > > &usedRecHits, double &Ehot, double &depth)