ForwardDetectors/plugins/CastorFastTowerProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CastorFastTowerProducer
0004 // Class:      CastorFastTowerProducer
0005 //
0006 /**\class CastorFastTowerProducer CastorFastTowerProducer.cc FastSimulation/ForwardDetectors/plugins/CastorFastTowerProducer.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Hans Van Haevermaet
0015 //         Created:  Thu Mar 13 12:00:56 CET 2008
0016 // $Id: CastorFastTowerProducer.cc,v 1.3 2011/03/04 10:52:06 hvanhaev Exp $
0017 //
0018 //
0019
0020 // system include files
0021 #include <memory>
0022 #include <vector>
0023 #include <iostream>
0024 #include <sstream>
0025 #include <TMath.h>
0026 #include <TRandom3.h>
0027 #include <TF1.h>
0028
0029 // user include files
0030 #include "FWCore/Framework/interface/Frameworkfwd.h"
0031
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036
0037 #include "DataFormats/Math/interface/Point3D.h"
0038
0039 // Castorobject includes
0040 #include "DataFormats/CastorReco/interface/CastorTower.h"
0041 #include "DataFormats/HcalRecHit/interface/CastorRecHit.h"
0042 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
0043
0044 // genCandidate particle includes
0045 #include "DataFormats/Candidate/interface/Candidate.h"
0046
0047 #include "FastSimulation/ForwardDetectors/plugins/CastorFastTowerProducer.h"
0048
0049 //
0050 // constructors and destructor
0051 //
0052 CastorFastTowerProducer::CastorFastTowerProducer(const edm::ParameterSet& iConfig)
0053     : tokGenPart_(consumes<reco::GenParticleCollection>(edm::InputTag{"genParticles"})) {
0054   //register your products
0055   produces<CastorTowerCollection>();
0056
0057   //now do what ever other initialization is needed
0058 }
0059
0060 CastorFastTowerProducer::~CastorFastTowerProducer() {
0061   // do anything here that needs to be done at desctruction time
0062   // (e.g. close files, deallocate resources etc.)
0063 }
0064
0065 //
0066 // member functions
0067 //
0068
0069 // ------------ method called to produce the data  ------------
0070 void CastorFastTowerProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0071   using namespace edm;
0072   using namespace reco;
0073   using namespace std;
0074   using namespace TMath;
0075
0076   //
0077   // Make CastorTower objects
0078   //
0079
0080   //cout << "entering event" << endl;
0081
0082   const edm::Handle<GenParticleCollection>& genParticles = iEvent.getHandle(tokGenPart_);
0083
0084   // make pointer to towers that will be made
0085   unique_ptr<CastorTowerCollection> CastorTowers(new CastorTowerCollection);
0086
0087   // declare castor array
0088   double castorplus[4][16];  // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = 5.9
0089   double castormin[4][16];  // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = -5.9
0090   // set phi values of array sectors and everything else to zero
0091   for (int j = 0; j < 16; j++) {
0092     castorplus[3][j] = -2.94524 + j * 0.3927;
0093     castormin[3][j] = -2.94524 + j * 0.3927;
0094     castorplus[0][j] = 0.;
0095     castormin[0][j] = 0.;
0096     castorplus[1][j] = 0.;
0097     castormin[1][j] = 0.;
0098     castorplus[2][j] = 0.;
0099     castormin[2][j] = 0.;
0100     //castorplus[4][j] = 0.;
0101     //castormin[4][j] = 0.;
0102   }
0103
0104   // declare properties vectors
0105   vector<double> depthplus[16];
0106   vector<double> depthmin[16];
0107   vector<double> fhotplus[16];
0108   vector<double> fhotmin[16];
0109   vector<double> energyplus[16];
0110   vector<double> energymin[16];
0111   //vector<double> femplus [16];
0112   //vector<double> femmin [16];
0113
0114   for (int i = 0; i < 16; i++) {
0115     depthplus[i].clear();
0116     depthmin[i].clear();
0117     fhotplus[i].clear();
0118     fhotmin[i].clear();
0119     energyplus[i].clear();
0120     energymin[i].clear();
0121     //femplus[i].clear();
0122     //femmin[i].clear();
0123   }
0124
0125   //cout << "declared everything" << endl;
0126
0127   // start particle loop
0128   for (size_t i = 0; i < genParticles->size(); ++i) {
0129     const Candidate& p = (*genParticles)[i];
0130
0131     // select particles in castor
0132     if (fabs(p.eta()) > 5.2 && fabs(p.eta()) < 6.6) {
0133       //cout << "found particle in castor, start calculating" << endl;
0134
0135       // declare energies
0136       //double gaus_E = -1.;
0137       double energy = -1.;
0138       double emEnergy = 0.;
0139       double hadEnergy = 0.;
0140       double fhot = 0.;
0141       double depth = 0.;
0142       //double fem = 0.;
0143
0144       // add energies - em: if particle is e- or gamma
0145       if (p.pdgId() == 11 || p.pdgId() == 22) {
0146         // calculate primary tower energy for electrons
0147         while (energy < 0.) {
0148           // apply energy smearing with gaussian random generator
0149           TRandom3 r(0);
0150           // use sigma/E parametrization from the full simulation
0151           double mean = 1.0024 * p.energy() - 0.3859;
0152           double sigma = 0.0228 * p.energy() + 2.1061;
0153           energy = r.Gaus(mean, sigma);
0154         }
0155
0156         // calculate electromagnetic electron/photon energy leakage
0157         double tmax;
0158         double a;
0159         double cte;
0160         if (p.pdgId() == 11) {
0161           cte = -0.5;
0162         } else {
0163           cte = 0.5;
0164         }
0165         tmax = 1.0 * (log(energy / 0.0015) + cte);
0166         a = tmax * 0.5 + 1;
0167         double leakage;
0168         double x = 0.5 * 19.38;
0169         leakage = energy - energy * Gamma(a, x);
0170
0171         // add emEnergy
0172         emEnergy = energy - leakage;
0173         // add hadEnergy leakage
0174         hadEnergy = leakage;
0175
0176         // calculate EM depth from parametrization
0177         double d0 = 0.2338 * pow(p.energy(), -0.1634);
0178         double d1 = 5.4336 * pow(p.energy(), 0.2410) + 14408.1025;
0179         double d2 = 1.4692 * pow(p.energy(), 0.1307) - 0.5216;
0180         if (d0 < 0.)
0181           d0 = abs(d0);
0182
0183         TF1* fdepth =
0184             new TF1("fdepth", "[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2])))", 14400., 14460.);
0185         fdepth->SetParameters(d0, d1, d2);
0186         depth = fdepth->GetRandom();
0187         fdepth->Delete();
0188         if (p.eta() < 0.)
0189           depth = -1 * depth;
0190
0191       } else {
0192         // calculate primary tower energy for hadrons
0193         while (energy < 0.) {
0194           // apply energy smearing with gaussian random generator
0195           TRandom3 r(0);
0196           // use sigma/E parametrization from the full simulation
0197           double mean = 0.8340 * p.energy() - 8.5054;
0198           double sigma = 0.1595 * p.energy() + 3.1183;
0199           energy = r.Gaus(mean, sigma);
0200         }
0201
0202         // add hadEnergy
0203         hadEnergy = energy;
0204
0205         // in the near future add fem parametrization
0206
0207         // calculate depth for HAD particle from parametrization
0208         double d0 = -0.000012 * p.energy() + 0.0661;
0209         double d1 = 785.7524 * pow(p.energy(), 0.0262) + 13663.4262;
0210         double d2 = 9.8748 * pow(p.energy(), 0.1720) + 37.0187;
0211         if (d0 < 0.)
0212           d0 = abs(d0);
0213
0214         TF1* fdepth =
0215             new TF1("fdepth", "[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2]) ))", 14400., 15500.);
0216         fdepth->SetParameters(d0, d1, d2);
0217         depth = fdepth->GetRandom();
0218         fdepth->Delete();
0219         if (p.eta() < 0.)
0220           depth = -1 * depth;
0221       }
0222
0223       // make tower
0224
0225       // set sector
0226       int sector = -1;
0227       for (int j = 0; j < 16; j++) {
0228         double a = -M_PI + j * 0.3927;
0229         double b = -M_PI + (j + 1) * 0.3927;
0230         if ((p.phi() > a) && (p.phi() < b)) {
0231           sector = j;
0232         }
0233       }
0234
0235       // set eta
0236       if (p.eta() > 0) {
0237         castorplus[0][sector] = castorplus[0][sector] + energy;
0238         castorplus[1][sector] = castorplus[1][sector] + emEnergy;
0239         castorplus[2][sector] = castorplus[2][sector] + hadEnergy;
0240
0241         depthplus[sector].push_back(depth);
0242         fhotplus[sector].push_back(fhot);
0243         energyplus[sector].push_back(energy);
0244         //cout << "filled vectors" << endl;
0245         //cout << "energyplus size = " << energyplus[sector].size() << endl;
0246         //cout << "depthplus size = " << depthplus[sector].size() << endl;
0247         //cout << "fhotplus size = " << fhotplus[sector].size() << endl;
0248
0249       } else {
0250         castormin[0][sector] = castormin[0][sector] + energy;
0251         castormin[1][sector] = castormin[1][sector] + emEnergy;
0252         castormin[2][sector] = castormin[2][sector] + hadEnergy;
0253
0254         depthmin[sector].push_back(depth);
0255         fhotmin[sector].push_back(fhot);
0256         energymin[sector].push_back(energy);
0257         //cout << "filled vectors" << endl;
0258       }
0259     }
0260   }
0261
0262   // add and substract pedestals/noise
0263   for (int j = 0; j < 16; j++) {
0264     double hadnoise = 0.;
0265     double emnoise = 0.;
0266     for (int i = 0; i < 12; i++) {
0267       hadnoise = hadnoise + make_noise();
0268       if (i < 2)
0269         emnoise = emnoise + make_noise();
0270     }
0271
0272     hadnoise = hadnoise - 12 * 0.053;
0273     emnoise = emnoise - 2 * 0.053;
0274     if (hadnoise < 0.)
0275       hadnoise = 0.;
0276     if (emnoise < 0.)
0277       emnoise = 0.;
0278     double totnoise = hadnoise + emnoise;
0279
0280     // add random noise
0281     castorplus[0][j] = castorplus[0][j] + totnoise;
0282     castormin[0][j] = castormin[0][j] + totnoise;
0283     castorplus[1][j] = castorplus[1][j] + emnoise;
0284     castormin[1][j] = castormin[1][j] + emnoise;
0285     castorplus[2][j] = castorplus[2][j] + hadnoise;
0286     castormin[2][j] = castormin[2][j] + hadnoise;
0287
0288     //cout << "after constant substraction" << endl;
0289     //cout << "total noise = " << castorplus[0][j] << " em noise = " << castorplus[1][j] << " had noise = " << castorplus[2][j] << endl;
0290     //cout << "fem should be = " << castorplus[1][j]/castorplus[0][j] << endl;
0291   }
0292
0293   // store towers from castor arrays
0294   // eta = 5.9
0295   for (int j = 0; j < 16; j++) {
0296     if (castorplus[0][j] != 0.) {
0297       double fem = 0.;
0298       fem = castorplus[1][j] / castorplus[0][j];
0299       TowerPoint pt1(88.5, 5.9, castorplus[3][j]);
0300       Point pt2(pt1);
0301
0302       //cout << "fem became = " << castorplus[1][j]/castorplus[0][j] << endl;
0303
0304       double depth_mean = 0.;
0305       double fhot_mean = 0.;
0306       double sum_energy = 0.;
0307
0308       //cout << "energyplus size = " << energyplus[j].size()<< endl;
0309       for (size_t p = 0; p < energyplus[j].size(); p++) {
0310         depth_mean = depth_mean + depthplus[j][p] * energyplus[j][p];
0311         fhot_mean = fhot_mean + fhotplus[j][p] * energyplus[j][p];
0312         sum_energy = sum_energy + energyplus[j][p];
0313       }
0314       depth_mean = depth_mean / sum_energy;
0315       fhot_mean = fhot_mean / sum_energy;
0316       //cout << "computed depth/fhot" << endl;
0317
0318       //std::vector<CastorCell> usedCells;
0319       edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
0320       CastorTowers->push_back(reco::CastorTower(
0321           castorplus[0][j], pt2, castorplus[1][j], castorplus[2][j], fem, depth_mean, fhot_mean, refvector));
0322     }
0323   }
0324   // eta = -5.9
0325   for (int j = 0; j < 16; j++) {
0326     if (castormin[0][j] != 0.) {
0327       double fem = 0.;
0328       fem = castormin[1][j] / castormin[0][j];
0329       TowerPoint pt1(88.5, -5.9, castormin[3][j]);
0330       Point pt2(pt1);
0331
0332       double depth_mean = 0.;
0333       double fhot_mean = 0.;
0334       double sum_energy = 0.;
0335
0336       for (size_t p = 0; p < energymin[j].size(); p++) {
0337         depth_mean = depth_mean + depthmin[j][p] * energymin[j][p];
0338         fhot_mean = fhot_mean + fhotmin[j][p] * energymin[j][p];
0339         sum_energy = sum_energy + energymin[j][p];
0340       }
0341       depth_mean = depth_mean / sum_energy;
0342       fhot_mean = fhot_mean / sum_energy;
0343
0344       //std::vector<CastorCell> usedCells;
0345       edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
0346       CastorTowers->push_back(reco::CastorTower(
0347           castormin[0][j], pt2, castormin[1][j], castormin[2][j], fem, depth_mean, fhot_mean, refvector));
0348     }
0349   }
0350
0351   iEvent.put(std::move(CastorTowers));
0352 }
0353
0354 double CastorFastTowerProducer::make_noise() {
0355   double result = -1.;
0356   TRandom3 r2(0);
0357   double mu_noise = 0.053;     // GeV (from 1.214 ADC) per channel
0358   double sigma_noise = 0.027;  // GeV (from 0.6168 ADC) per channel
0359
0360   while (result < 0.) {
0361     result = r2.Gaus(mu_noise, sigma_noise);
0362   }
0363
0364   return result;
0365 }
0366
0367 //define this as a plug-in
0368 DEFINE_FWK_MODULE(CastorFastTowerProducer);