|
|
|
#include <FastSimulation/ForwardDetectors/plugins/CastorFastClusterProducer.cc>
Public Member Functions | |
| CastorFastClusterProducer (const edm::ParameterSet &) | |
| ~CastorFastClusterProducer () | |
Private Types | |
| typedef std::vector < reco::CastorCluster > | CastorClusterCollection |
| typedef ROOT::Math::RhoEtaPhiPoint | ClusterPoint |
| typedef math::XYZPointD | Point |
Private Member Functions | |
| double | make_noise () |
| virtual void | produce (edm::Event &, const edm::EventSetup &) override |
Description: <one line="" class="" summary>="">
Implementation: <Notes on="" implementation>="">
Definition at line 16 of file CastorFastClusterProducer.h.
typedef std::vector<reco::CastorCluster> CastorFastClusterProducer::CastorClusterCollection [private] |
Definition at line 28 of file CastorFastClusterProducer.h.
typedef ROOT::Math::RhoEtaPhiPoint CastorFastClusterProducer::ClusterPoint [private] |
Definition at line 27 of file CastorFastClusterProducer.h.
typedef math::XYZPointD CastorFastClusterProducer::Point [private] |
Definition at line 26 of file CastorFastClusterProducer.h.
| CastorFastClusterProducer::CastorFastClusterProducer | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 56 of file CastorFastClusterProducer.cc.
{
//register your products
produces<CastorClusterCollection>();
//now do what ever other initialization is needed
}
| CastorFastClusterProducer::~CastorFastClusterProducer | ( | ) |
Definition at line 66 of file CastorFastClusterProducer.cc.
{
// do anything here that needs to be done at desctruction time
// (e.g. close files, deallocate resources etc.)
}
| double CastorFastClusterProducer::make_noise | ( | ) | [private] |
Definition at line 348 of file CastorFastClusterProducer.cc.
References diffTwoXMLs::r2, and query::result.
| void CastorFastClusterProducer::produce | ( | edm::Event & | iEvent, |
| const edm::EventSetup & | iSetup | ||
| ) | [override, private, virtual] |
Implements edm::EDProducer.
Definition at line 81 of file CastorFastClusterProducer.cc.
References a, reco::Candidate::energy(), reco::Candidate::eta(), Gamma, genParticleCandidates2GenParticles_cfi::genParticles, edm::Event::getByLabel(), i, create_public_lumi_plots::log, AlCaHLTBitMon_ParallelJobs::p, reco::Candidate::pdgId(), reco::Candidate::phi(), funct::pow(), edm::RefVector< C, T, F >::push_back(), edm::Event::put(), alignCSCRings::r, L1Trigger_dataformats::reco, mathSSE::sqrt(), groupFilesInBlocks::temp, tmax, and x.
{
using namespace edm;
using namespace reco;
using namespace std;
using namespace TMath;
//
// Make CastorCluster objects
//
//cout << "entering event" << endl;
Handle<GenParticleCollection> genParticles;
iEvent.getByLabel("genParticles", genParticles);
// make pointer to towers that will be made
auto_ptr<CastorClusterCollection> CastorClusters (new CastorClusterCollection);
/*
// 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];
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();
}
*/
//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 emEnergy = 0.;
double hadEnergy = 0.;
//double fhot = 0.;
//double depth = 0.;
// add energies - em: if particle is e- or gamma
if (p.pdgId() == 11 || p.pdgId() == 22) {
while ( gaus_E < 0.) {
// apply energy smearing with gaussian random generator
TRandom3 r(0);
// use sigma/E parametrization for the EM sections of CASTOR TB 2007 results
double sigma = p.energy()*(sqrt(pow(0.044,2) + pow(0.513/sqrt(p.energy()),2)));
gaus_E = r.Gaus(p.energy(),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(gaus_E/0.0015)+cte);
a = tmax*0.5 + 1;
double leakage;
double x = 0.5*19.38;
leakage = gaus_E - gaus_E*Gamma(a,x);
// add emEnergy
emEnergy = gaus_E - leakage;
// add hadEnergy leakage
hadEnergy = leakage;
// make cluster
ClusterPoint pt1;
if (p.eta() > 0.) {ClusterPoint temp(88.5,5.9,p.phi()); pt1 = temp;}
if (p.eta() < 0.) {ClusterPoint temp(88.5,-5.9,p.phi()); pt1 = temp;}
Point pt2(pt1);
CastorTowerRefVector refvector;
CastorClusters->push_back(reco::CastorCluster(gaus_E,pt2,emEnergy,hadEnergy,emEnergy/gaus_E,0.,0.,0.,0.,refvector));
} else {
while (gaus_E < 0.) {
// apply energy smearing with gaussian random generator
TRandom3 r(0);
// use sigma/E parametrization for the HAD sections of CASTOR TB 2007 results
double sigma = p.energy()*(sqrt(pow(0.121,2) + pow(1.684/sqrt(p.energy()),2)));
gaus_E = r.Gaus(p.energy(),sigma);
}
// add hadEnergy
hadEnergy = gaus_E;
// make cluster
ClusterPoint pt1;
if (p.eta() > 0.) {ClusterPoint temp(88.5,5.9,p.phi()); pt1 = temp;}
if (p.eta() < 0.) {ClusterPoint temp(88.5,-5.9,p.phi()); pt1 = temp;}
Point pt2(pt1);
CastorTowerRefVector refvector;
CastorClusters->push_back(reco::CastorCluster(gaus_E,pt2,0.,hadEnergy,0.,0.,0.,0.,0.,refvector));
}
/*
// 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] + gaus_E;
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(gaus_E);
//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] + gaus_E;
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(gaus_E);
//cout << "filled vectors" << endl;
}
*/
}
}
/*
// substract pedestals/noise
for (int j = 0; j < 16; j++) {
double hadnoise = 0.;
for (int i=0;i<12;i++) {
hadnoise = hadnoise + make_noise();
}
castorplus[0][j] = castorplus[0][j] - hadnoise - make_noise() - make_noise();
castormin[0][j] = castormin[0][j] - hadnoise - make_noise() - make_noise();
castorplus[1][j] = castorplus[1][j] - make_noise() - make_noise();
castormin[1][j] = castormin[1][j] - make_noise() - make_noise();
castorplus[2][j] = castorplus[2][j] - hadnoise;
castormin[2][j] = castormin[2][j] - hadnoise;
// set possible negative values to zero
if (castorplus[0][j] < 0.) castorplus[0][j] = 0.;
if (castormin[0][j] < 0.) castormin[0][j] = 0.;
if (castorplus[1][j] < 0.) castorplus[1][j] = 0.;
if (castormin[1][j] < 0.) castormin[1][j] = 0.;
if (castorplus[2][j] < 0.) castorplus[2][j] = 0.;
if (castormin[2][j] < 0.) castormin[2][j] = 0.;
}
*/
/*
// 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];
ClusterPoint pt1(88.5,5.9,castorplus[3][j]);
Point pt2(pt1);
// parametrize depth and fhot from full sim
// get fit parameters from energy
// get random number according to distribution with fit parameters
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;
edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
CastorClusters->push_back(reco::CastorCluster(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];
ClusterPoint pt1(88.5,-5.9,castormin[3][j]);
Point pt2(pt1);
// parametrize depth and fhot from full sim
// get fit parameters from energy
// get random number according to distribution with fit parameters
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;
edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
CastorClusters->push_back(reco::CastorCluster(castormin[0][j],pt2,castormin[1][j],castormin[2][j],fem,depth_mean,fhot_mean,refvector));
}
}
*/
iEvent.put(CastorClusters);
}
1.7.3