58 produces<CastorTowerCollection>();
85 using namespace TMath;
100 double castorplus [4][16];
101 double castormin [4][16];
103 for (
int j = 0;
j < 16;
j++) {
104 castorplus[3][
j] = -2.94524 +
j*0.3927;
105 castormin[3][
j] = -2.94524 +
j*0.3927;
106 castorplus[0][
j] = 0.;
107 castormin[0][
j] = 0.;
108 castorplus[1][
j] = 0.;
109 castormin[1][
j] = 0.;
110 castorplus[2][
j] = 0.;
111 castormin[2][
j] = 0.;
117 vector<double> depthplus[16];
118 vector<double> depthmin[16];
119 vector<double> fhotplus [16];
120 vector<double> fhotmin [16];
121 vector<double> energyplus [16];
122 vector<double> energymin [16];
126 for (
int i=0;
i<16;
i++) {
127 depthplus[
i].clear();
131 energyplus[
i].clear();
132 energymin[
i].clear();
144 if ( fabs(p.
eta()) > 5.2 && fabs(p.
eta()) < 6.6) {
151 double emEnergy = 0.;
152 double hadEnergy = 0.;
161 while ( energy < 0.) {
167 energy = r.Gaus(mean,sigma);
174 if ( p.
pdgId() == 11) { cte = -0.5; }
else { cte = 0.5; }
175 tmax = 1.0*(
log(energy/0.0015)+cte);
178 double x = 0.5*19.38;
179 leakage = energy - energy*
Gamma(a,x);
182 emEnergy = energy - leakage;
188 double d1 = 5.4336 *
pow(p.
energy(),0.2410) + 14408.1025;
189 double d2 = 1.4692 *
pow(p.
energy(),0.1307) - 0.5216;
190 if (d0 < 0.) d0 =
abs(d0);
192 TF1 *fdepth =
new TF1(
"fdepth",
"[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2])))",14400.,14460.);
193 fdepth->SetParameters(d0,d1,d2);
194 depth = fdepth->GetRandom();
196 if (p.
eta() < 0.) depth = -1*depth;
201 while (energy < 0.) {
207 energy = r.Gaus(mean,sigma);
216 double d0 = -0.000012 * p.
energy() + 0.0661;
217 double d1 = 785.7524 *
pow(p.
energy(),0.0262) + 13663.4262;
218 double d2 = 9.8748 *
pow(p.
energy(),0.1720) + 37.0187;
219 if (d0 < 0.) d0 =
abs(d0);
221 TF1 *fdepth =
new TF1(
"fdepth",
"[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2]) ))",14400.,15500.);
222 fdepth->SetParameters(d0,d1,d2);
223 depth = fdepth->GetRandom();
225 if (p.
eta() < 0.) depth = -1*depth;
234 for (
int j = 0;
j < 16;
j++) {
235 double a = -
M_PI +
j*0.3927;
236 double b = -
M_PI + (
j+1)*0.3927;
237 if ( (p.
phi() >
a) && (p.
phi() <
b)) {
244 castorplus[0][sector] = castorplus[0][sector] +
energy;
245 castorplus[1][sector] = castorplus[1][sector] + emEnergy;
246 castorplus[2][sector] = castorplus[2][sector] + hadEnergy;
248 depthplus[sector].push_back(depth);
249 fhotplus[sector].push_back(fhot);
250 energyplus[sector].push_back(energy);
257 castormin[0][sector] = castormin[0][sector] +
energy;
258 castormin[1][sector] = castormin[1][sector] + emEnergy;
259 castormin[2][sector] = castormin[2][sector] + hadEnergy;
262 depthmin[sector].push_back(depth);
263 fhotmin[sector].push_back(fhot);
264 energymin[sector].push_back(energy);
275 for (
int j = 0;
j < 16;
j++) {
276 double hadnoise = 0.;
278 for (
int i=0;
i<12;
i++) {
283 hadnoise = hadnoise - 12*0.053;
284 emnoise = emnoise - 2*0.053;
285 if ( hadnoise < 0.) hadnoise = 0.;
286 if ( emnoise < 0.) emnoise = 0.;
287 double totnoise = hadnoise + emnoise;
290 castorplus[0][
j] = castorplus[0][
j] + totnoise;
291 castormin[0][
j] = castormin[0][
j] + totnoise;
292 castorplus[1][
j] = castorplus[1][
j] + emnoise;
293 castormin[1][
j] = castormin[1][
j] + emnoise;
294 castorplus[2][
j] = castorplus[2][
j] + hadnoise;
295 castormin[2][
j] = castormin[2][
j] + hadnoise;
306 for (
int j=0;
j<16;
j++) {
307 if (castorplus[0][
j] != 0.) {
310 fem = castorplus[1][
j]/castorplus[0][
j];
316 double depth_mean = 0.;
317 double fhot_mean = 0.;
318 double sum_energy = 0.;
321 for (
size_t p = 0;
p<energyplus[
j].size();
p++) {
322 depth_mean = depth_mean + depthplus[
j][
p]*energyplus[
j][
p];
323 fhot_mean = fhot_mean + fhotplus[
j][
p]*energyplus[
j][
p];
324 sum_energy = sum_energy + energyplus[
j][
p];
326 depth_mean = depth_mean/sum_energy;
327 fhot_mean = fhot_mean/sum_energy;
333 CastorTowers->
push_back(
reco::CastorTower(castorplus[0][j],pt2,castorplus[1][j],castorplus[2][j],fem,depth_mean,fhot_mean,refvector));
337 for (
int j=0;
j<16;
j++) {
338 if (castormin[0][
j] != 0.) {
340 fem = castormin[1][
j]/castormin[0][
j];
344 double depth_mean = 0.;
345 double fhot_mean = 0.;
346 double sum_energy = 0.;
349 for (
size_t p = 0;
p<energymin[
j].size();
p++) {
350 depth_mean = depth_mean + depthmin[
j][
p]*energymin[
j][
p];
351 fhot_mean = fhot_mean + fhotmin[
j][
p]*energymin[
j][
p];
352 sum_energy = sum_energy + energymin[
j][
p];
354 depth_mean = depth_mean/sum_energy;
355 fhot_mean = fhot_mean/sum_energy;
360 CastorTowers->
push_back(
reco::CastorTower(castormin[0][j],pt2,castormin[1][j],castormin[2][j],fem,depth_mean,fhot_mean,refvector));
364 iEvent.
put(CastorTowers);
371 double mu_noise = 0.053;
372 double sigma_noise = 0.027;
374 while (result < 0.) {
375 result = r2.Gaus(mu_noise,sigma_noise);
virtual double energy() const =0
energy
virtual void beginRun(edm::Run &, edm::EventSetup const &)
CastorFastTowerProducer(const edm::ParameterSet &)
std::vector< reco::CastorTower > CastorTowerCollection
~CastorFastTowerProducer()
#define DEFINE_FWK_MODULE(type)
ROOT::Math::RhoEtaPhiPoint TowerPoint
virtual void produce(edm::Event &, const edm::EventSetup &)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
virtual int pdgId() const =0
PDG identifier.
static const double tmax[3]
Log< T >::type log(const T &t)
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Power< A, B >::type pow(const A &a, const B &b)
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity