CastorFastTowerProducer Class Reference

#include <FastSimulation/ForwardDetectors/plugins/CastorFastTowerProducer.cc>

Inheritance diagram for CastorFastTowerProducer:

Public Member Functions
	CastorFastTowerProducer (const edm::ParameterSet &)
	~CastorFastTowerProducer () override
Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default
	EDProducer (const EDProducer &)=delete
bool	hasAbilityToProduceInBeginLumis () const final
bool	hasAbilityToProduceInBeginProcessBlocks () const final
bool	hasAbilityToProduceInBeginRuns () const final
bool	hasAbilityToProduceInEndLumis () const final
bool	hasAbilityToProduceInEndProcessBlocks () const final
bool	hasAbilityToProduceInEndRuns () const final
const EDProducer &	operator= (const EDProducer &)=delete

Private Types
typedef std::vector < reco::CastorTower >	CastorTowerCollection
typedef math::XYZPointD	Point
typedef ROOT::Math::RhoEtaPhiPoint	TowerPoint

Private Member Functions
double	make_noise ()
void	produce (edm::Event &, const edm::EventSetup &) override

Additional Inherited Members
Public Types inherited from edm::stream::EDProducer<>
using	CacheTypes = CacheContexts< T...>
using	GlobalCache = typename CacheTypes::GlobalCache
using	HasAbility = AbilityChecker< T...>
using	InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache
using	LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache
using	LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >
using	LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache
using	RunCache = typename CacheTypes::RunCache
using	RunContext = RunContextT< RunCache, GlobalCache >
using	RunSummaryCache = typename CacheTypes::RunSummaryCache

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Definition at line 17 of file CastorFastTowerProducer.h.

Member Typedef Documentation

typedef std::vector<reco::CastorTower> CastorFastTowerProducer::CastorTowerCollection

private

Definition at line 29 of file CastorFastTowerProducer.h.

typedef math::XYZPointD CastorFastTowerProducer::Point

private

Definition at line 27 of file CastorFastTowerProducer.h.

typedef ROOT::Math::RhoEtaPhiPoint CastorFastTowerProducer::TowerPoint

private

Definition at line 28 of file CastorFastTowerProducer.h.

Constructor & Destructor Documentation

CastorFastTowerProducer::CastorFastTowerProducer ( const edm::ParameterSet &  iConfig )

explicit

Definition at line 54 of file CastorFastTowerProducer.cc.

                                                                               {
  //register your products
  produces<CastorTowerCollection>();

  //now do what ever other initialization is needed
}

CastorFastTowerProducer::~CastorFastTowerProducer ( )

override

Definition at line 61 of file CastorFastTowerProducer.cc.

                                                  {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}

Member Function Documentation

double CastorFastTowerProducer::make_noise ( )

private

Definition at line 356 of file CastorFastTowerProducer.cc.

References diffTwoXMLs::r2, and mps_fire::result.

Referenced by produce().

                                           {
  double result = -1.;
  TRandom3 r2(0);
  double mu_noise = 0.053;     // GeV (from 1.214 ADC) per channel
  double sigma_noise = 0.027;  // GeV (from 0.6168 ADC) per channel

  while (result < 0.) {
    result = r2.Gaus(mu_noise, sigma_noise);
  }

  return result;
}

void CastorFastTowerProducer::produce ( edm::Event &  iEvent, const edm::EventSetup &  iSetup  )

overrideprivate

Definition at line 71 of file CastorFastTowerProducer.cc.

References a, funct::abs(), b, d0, d1, HLT_FULL_cff::depth, relval_parameters_module::energy, reco::Candidate::energy(), reco::Candidate::eta(), hcaldqm::quantity::fdepth, genParticleCandidates2GenParticles_cfi::genParticles, edm::Event::getByLabel(), mps_fire::i, dqmiolumiharvest::j, log, M_PI, make_noise(), SiStripPI::mean, eostools::move(), AlCaHLTBitMon_ParallelJobs::p, reco::Candidate::pdgId(), reco::Candidate::phi(), funct::pow(), edm::RefVector< C, T, F >::push_back(), edm::Event::put(), alignCSCRings::r, dt_dqm_sourceclient_common_cff::reco, tmax, and x.

                                                                                   {
  using namespace edm;
  using namespace reco;
  using namespace std;
  using namespace TMath;

  //
  // Make CastorTower objects
  //

  //cout << "entering event" << endl;

  Handle<GenParticleCollection> genParticles;
  iEvent.getByLabel("genParticles", genParticles);

  // make pointer to towers that will be made
  unique_ptr<CastorTowerCollection> CastorTowers(new CastorTowerCollection);

  // declare castor array
  double castorplus[4][16];  // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = 5.9
  double castormin[4][16];  // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = -5.9
  // set phi values of array sectors and everything else to zero
  for (int j = 0; j < 16; j++) {
    castorplus[3][j] = -2.94524 + j * 0.3927;
    castormin[3][j] = -2.94524 + j * 0.3927;
    castorplus[0][j] = 0.;
    castormin[0][j] = 0.;
    castorplus[1][j] = 0.;
    castormin[1][j] = 0.;
    castorplus[2][j] = 0.;
    castormin[2][j] = 0.;
    //castorplus[4][j] = 0.;
    //castormin[4][j] = 0.;
  }

  // declare properties vectors
  vector<double> depthplus[16];
  vector<double> depthmin[16];
  vector<double> fhotplus[16];
  vector<double> fhotmin[16];
  vector<double> energyplus[16];
  vector<double> energymin[16];
  //vector<double> femplus [16];
  //vector<double> femmin [16];

  for (int i = 0; i < 16; i++) {
    depthplus[i].clear();
    depthmin[i].clear();
    fhotplus[i].clear();
    fhotmin[i].clear();
    energyplus[i].clear();
    energymin[i].clear();
    //femplus[i].clear();
    //femmin[i].clear();
  }

  //cout << "declared everything" << endl;

  // start particle loop
  for (size_t i = 0; i < genParticles->size(); ++i) {
    const Candidate& p = (*genParticles)[i];

    // select particles in castor
    if (fabs(p.eta()) > 5.2 && fabs(p.eta()) < 6.6) {
      //cout << "found particle in castor, start calculating" << endl;

      // declare energies
      //double gaus_E = -1.;
      double energy = -1.;
      double emEnergy = 0.;
      double hadEnergy = 0.;
      double fhot = 0.;
      double depth = 0.;
      //double fem = 0.;

      // add energies - em: if particle is e- or gamma
      if (p.pdgId() == 11 || p.pdgId() == 22) {
        // calculate primary tower energy for electrons
        while (energy < 0.) {
          // apply energy smearing with gaussian random generator
          TRandom3 r(0);
          // use sigma/E parametrization from the full simulation
          double mean = 1.0024 * p.energy() - 0.3859;
          double sigma = 0.0228 * p.energy() + 2.1061;
          energy = r.Gaus(mean, sigma);
        }

        // calculate electromagnetic electron/photon energy leakage
        double tmax;
        double a;
        double cte;
        if (p.pdgId() == 11) {
          cte = -0.5;
        } else {
          cte = 0.5;
        }
        tmax = 1.0 * (log(energy / 0.0015) + cte);
        a = tmax * 0.5 + 1;
        double leakage;
        double x = 0.5 * 19.38;
        leakage = energy - energy * Gamma(a, x);

        // add emEnergy
        emEnergy = energy - leakage;
        // add hadEnergy leakage
        hadEnergy = leakage;

        // calculate EM depth from parametrization
        double d0 = 0.2338 * pow(p.energy(), -0.1634);
        double d1 = 5.4336 * pow(p.energy(), 0.2410) + 14408.1025;
        double d2 = 1.4692 * pow(p.energy(), 0.1307) - 0.5216;
        if (d0 < 0.)
          d0 = abs(d0);

        TF1* fdepth =
            new TF1("fdepth", "[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2])))", 14400., 14460.);
        fdepth->SetParameters(d0, d1, d2);
        depth = fdepth->GetRandom();
        fdepth->Delete();
        if (p.eta() < 0.)
          depth = -1 * depth;

      } else {
        // calculate primary tower energy for hadrons
        while (energy < 0.) {
          // apply energy smearing with gaussian random generator
          TRandom3 r(0);
          // use sigma/E parametrization from the full simulation
          double mean = 0.8340 * p.energy() - 8.5054;
          double sigma = 0.1595 * p.energy() + 3.1183;
          energy = r.Gaus(mean, sigma);
        }

        // add hadEnergy
        hadEnergy = energy;

        // in the near future add fem parametrization

        // calculate depth for HAD particle from parametrization
        double d0 = -0.000012 * p.energy() + 0.0661;
        double d1 = 785.7524 * pow(p.energy(), 0.0262) + 13663.4262;
        double d2 = 9.8748 * pow(p.energy(), 0.1720) + 37.0187;
        if (d0 < 0.)
          d0 = abs(d0);

        TF1* fdepth =
            new TF1("fdepth", "[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2]) ))", 14400., 15500.);
        fdepth->SetParameters(d0, d1, d2);
        depth = fdepth->GetRandom();
        fdepth->Delete();
        if (p.eta() < 0.)
          depth = -1 * depth;
      }

      // make tower

      // set sector
      int sector = -1;
      for (int j = 0; j < 16; j++) {
        double a = -M_PI + j * 0.3927;
        double b = -M_PI + (j + 1) * 0.3927;
        if ((p.phi() > a) && (p.phi() < b)) {
          sector = j;
        }
      }

      // set eta
      if (p.eta() > 0) {
        castorplus[0][sector] = castorplus[0][sector] + energy;
        castorplus[1][sector] = castorplus[1][sector] + emEnergy;
        castorplus[2][sector] = castorplus[2][sector] + hadEnergy;

        depthplus[sector].push_back(depth);
        fhotplus[sector].push_back(fhot);
        energyplus[sector].push_back(energy);
        //cout << "filled vectors" << endl;
        //cout << "energyplus size = " << energyplus[sector].size() << endl;
        //cout << "depthplus size = " << depthplus[sector].size() << endl;
        //cout << "fhotplus size = " << fhotplus[sector].size() << endl;

      } else {
        castormin[0][sector] = castormin[0][sector] + energy;
        castormin[1][sector] = castormin[1][sector] + emEnergy;
        castormin[2][sector] = castormin[2][sector] + hadEnergy;

        depthmin[sector].push_back(depth);
        fhotmin[sector].push_back(fhot);
        energymin[sector].push_back(energy);
        //cout << "filled vectors" << endl;
      }
    }
  }

  // add and substract pedestals/noise
  for (int j = 0; j < 16; j++) {
    double hadnoise = 0.;
    double emnoise = 0.;
    for (int i = 0; i < 12; i++) {
      hadnoise = hadnoise + make_noise();
      if (i < 2)
        emnoise = emnoise + make_noise();
    }

    hadnoise = hadnoise - 12 * 0.053;
    emnoise = emnoise - 2 * 0.053;
    if (hadnoise < 0.)
      hadnoise = 0.;
    if (emnoise < 0.)
      emnoise = 0.;
    double totnoise = hadnoise + emnoise;

    // add random noise
    castorplus[0][j] = castorplus[0][j] + totnoise;
    castormin[0][j] = castormin[0][j] + totnoise;
    castorplus[1][j] = castorplus[1][j] + emnoise;
    castormin[1][j] = castormin[1][j] + emnoise;
    castorplus[2][j] = castorplus[2][j] + hadnoise;
    castormin[2][j] = castormin[2][j] + hadnoise;

    //cout << "after constant substraction" << endl;
    //cout << "total noise = " << castorplus[0][j] << " em noise = " << castorplus[1][j] << " had noise = " << castorplus[2][j] << endl;
    //cout << "fem should be = " << castorplus[1][j]/castorplus[0][j] << endl;
  }

  // store towers from castor arrays
  // eta = 5.9
  for (int j = 0; j < 16; j++) {
    if (castorplus[0][j] != 0.) {
      double fem = 0.;
      fem = castorplus[1][j] / castorplus[0][j];
      TowerPoint pt1(88.5, 5.9, castorplus[3][j]);
      Point pt2(pt1);

      //cout << "fem became = " << castorplus[1][j]/castorplus[0][j] << endl;

      double depth_mean = 0.;
      double fhot_mean = 0.;
      double sum_energy = 0.;

      //cout << "energyplus size = " << energyplus[j].size()<< endl;
      for (size_t p = 0; p < energyplus[j].size(); p++) {
        depth_mean = depth_mean + depthplus[j][p] * energyplus[j][p];
        fhot_mean = fhot_mean + fhotplus[j][p] * energyplus[j][p];
        sum_energy = sum_energy + energyplus[j][p];
      }
      depth_mean = depth_mean / sum_energy;
      fhot_mean = fhot_mean / sum_energy;
      //cout << "computed depth/fhot" << endl;

      //std::vector<CastorCell> usedCells;
      edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
      CastorTowers->push_back(reco::CastorTower(
          castorplus[0][j], pt2, castorplus[1][j], castorplus[2][j], fem, depth_mean, fhot_mean, refvector));
    }
  }
  // eta = -5.9
  for (int j = 0; j < 16; j++) {
    if (castormin[0][j] != 0.) {
      double fem = 0.;
      fem = castormin[1][j] / castormin[0][j];
      TowerPoint pt1(88.5, -5.9, castormin[3][j]);
      Point pt2(pt1);

      double depth_mean = 0.;
      double fhot_mean = 0.;
      double sum_energy = 0.;

      for (size_t p = 0; p < energymin[j].size(); p++) {
        depth_mean = depth_mean + depthmin[j][p] * energymin[j][p];
        fhot_mean = fhot_mean + fhotmin[j][p] * energymin[j][p];
        sum_energy = sum_energy + energymin[j][p];
      }
      depth_mean = depth_mean / sum_energy;
      fhot_mean = fhot_mean / sum_energy;

      //std::vector<CastorCell> usedCells;
      edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
      CastorTowers->push_back(reco::CastorTower(
          castormin[0][j], pt2, castormin[1][j], castormin[2][j], fem, depth_mean, fhot_mean, refvector));
    }
  }

  iEvent.put(std::move(CastorTowers));
}