|
|
|
#include <FastSimulation/ForwardDetectors/plugins/CastorFastTowerProducer.cc>
Public Member Functions | |
| CastorFastTowerProducer (const edm::ParameterSet &) | |
| ~CastorFastTowerProducer () | |
Private Types | |
| typedef std::vector < reco::CastorTower > | CastorTowerCollection |
| typedef math::XYZPointD | Point |
| typedef ROOT::Math::RhoEtaPhiPoint | TowerPoint |
Private Member Functions | |
| virtual void | beginRun (edm::Run &, edm::EventSetup const &) |
| virtual void | endRun () |
| double | make_noise () |
| virtual void | produce (edm::Event &, const edm::EventSetup &) |
Description: <one line="" class="" summary>="">
Implementation: <Notes on="" implementation>="">
Definition at line 17 of file CastorFastTowerProducer.h.
typedef std::vector<reco::CastorTower> CastorFastTowerProducer::CastorTowerCollection [private] |
Definition at line 31 of file CastorFastTowerProducer.h.
typedef math::XYZPointD CastorFastTowerProducer::Point [private] |
Definition at line 29 of file CastorFastTowerProducer.h.
typedef ROOT::Math::RhoEtaPhiPoint CastorFastTowerProducer::TowerPoint [private] |
Definition at line 30 of file CastorFastTowerProducer.h.
| CastorFastTowerProducer::CastorFastTowerProducer | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 55 of file CastorFastTowerProducer.cc.
{
//register your products
produces<CastorTowerCollection>();
//now do what ever other initialization is needed
}
| CastorFastTowerProducer::~CastorFastTowerProducer | ( | ) |
Definition at line 65 of file CastorFastTowerProducer.cc.
{
// do anything here that needs to be done at desctruction time
// (e.g. close files, deallocate resources etc.)
}
| void CastorFastTowerProducer::beginRun | ( | edm::Run & | run, |
| edm::EventSetup const & | es | ||
| ) | [private, virtual] |
| void CastorFastTowerProducer::endRun | ( | void | ) | [private, virtual] |
Definition at line 390 of file CastorFastTowerProducer.cc.
{
}
| double CastorFastTowerProducer::make_noise | ( | ) | [private] |
Definition at line 368 of file CastorFastTowerProducer.cc.
References query::result.
Referenced by produce().
| void CastorFastTowerProducer::produce | ( | edm::Event & | iEvent, |
| const edm::EventSetup & | iSetup | ||
| ) | [private, virtual] |
Implements edm::EDProducer.
Definition at line 80 of file CastorFastTowerProducer.cc.
References a, abs, b, debug_cff::d0, debug_cff::d1, reco::Candidate::energy(), relval_parameters_module::energy, reco::Candidate::eta(), Gamma, genParticleCandidates2GenParticles_cfi::genParticles, edm::Event::getByLabel(), i, j, funct::log(), M_PI, make_noise(), plotscripts::mean(), L1TEmulatorMonitor_cff::p, reco::Candidate::pdgId(), reco::Candidate::phi(), funct::pow(), edm::RefVector< C, T, F >::push_back(), edm::Event::put(), csvReporter::r, dt_offlineAnalysis_common_cff::reco, ExpressReco_HICollisions_FallBack::sigma, tmax, and ExpressReco_HICollisions_FallBack::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
auto_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;
CastorCellRefVector 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;
CastorCellRefVector refvector;
CastorTowers->push_back(reco::CastorTower(castormin[0][j],pt2,castormin[1][j],castormin[2][j],fem,depth_mean,fhot_mean,refvector));
}
}
iEvent.put(CastorTowers);
}
1.7.3