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
bool	hasAbilityToProduceInBeginLumis () const final
bool	hasAbilityToProduceInBeginRuns () const final
bool	hasAbilityToProduceInEndLumis () const final
bool	hasAbilityToProduceInEndRuns () const final

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<>
typedef CacheContexts< T... >	CacheTypes
typedef CacheTypes::GlobalCache	GlobalCache
typedef AbilityChecker< T... >	HasAbility
typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache
typedef LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext
typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache
typedef CacheTypes::RunCache	RunCache
typedef RunContextT< RunCache, GlobalCache >	RunContext
typedef CacheTypes::RunSummaryCache	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 DEFINE_FWK_MODULE, 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, HLTMuonOfflineAnalyzer_cfi::d0, LEDCalibrationChannels::depth, HCALHighEnergyHPDFilter_cfi::energy, reco::Candidate::energy(), reco::Candidate::eta(), hcaldqm::quantity::fdepth, genParticles2HepMC_cfi::genParticles, edm::Event::getByLabel(), mps_fire::i, dqmiolumiharvest::j, dqm-mbProfile::log, M_PI, make_noise(), SiStripPI::mean, eostools::move(), AlCaHLTBitMon_ParallelJobs::p, reco::Candidate::pdgId(), reco::Candidate::phi(), funct::pow(), HLT_2018_cff::pt1, HLT_2018_cff::pt2, edm::RefVector< C, T, F >::push_back(), edm::Event::put(), alignCSCRings::r, tmax, and x.

Referenced by JSONExport.JsonExport::export(), HTMLExport.HTMLExport::export(), and HTMLExport.HTMLExportStatic::export().

                                                                                   {
  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));
}