ForwardDetectors/plugins/CastorFastClusterProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CastorFastClusterProducer
0004 // Class:      CastorFastClusterProducer
0005 //
0006 /**\class CastorFastClusterProducer CastorFastClusterProducer.cc FastSimulation/ForwardDetectors/plugins/CastorFastClusterProducer.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: CastorFastClusterProducer.cc,v 1.2 2011/03/04 11:00:55 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/CastorCluster.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/CastorFastClusterProducer.h"
0048
0049 //
0050 // constructors and destructor
0051 //
0052 CastorFastClusterProducer::CastorFastClusterProducer(const edm::ParameterSet& iConfig)
0053     : tokGenPart_(consumes<reco::GenParticleCollection>(edm::InputTag{"genParticles"})) {
0054   //register your products
0055   produces<CastorClusterCollection>();
0056
0057   //now do what ever other initialization is needed
0058 }
0059
0060 CastorFastClusterProducer::~CastorFastClusterProducer() {
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 CastorFastClusterProducer::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 CastorCluster objects
0078   //
0079
0080   //cout << "entering event" << endl;
0081
0082   const edm::Handle<reco::GenParticleCollection>& genParticles = iEvent.getHandle(tokGenPart_);
0083
0084   // make pointer to towers that will be made
0085   unique_ptr<CastorClusterCollection> CastorClusters(new CastorClusterCollection);
0086
0087   /*
0088    // declare castor array
0089    double castorplus [4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = 5.9
0090    double castormin [4][16];  // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = -5.9
0091    // set phi values of array sectors and everything else to zero
0092    for (int j = 0; j < 16; j++) {
0093     castorplus[3][j] = -2.94524 + j*0.3927;
0094     castormin[3][j] = -2.94524 + j*0.3927;
0095     castorplus[0][j] = 0.;
0096     castormin[0][j] = 0.;
0097     castorplus[1][j] = 0.;
0098     castormin[1][j] = 0.;
0099     castorplus[2][j] = 0.;
0100     castormin[2][j] = 0.;
0101     //castorplus[4][j] = 0.;
0102     //castormin[4][j] = 0.;
0103    }
0104
0105    // declare properties vectors
0106    vector<double> depthplus[16];
0107    vector<double> depthmin[16];
0108    vector<double> fhotplus [16];
0109    vector<double> fhotmin [16];
0110    vector<double> energyplus [16];
0111    vector<double> energymin [16];
0112
0113    for (int i=0;i<16;i++) {
0114     depthplus[i].clear();
0115     depthmin[i].clear();
0116     fhotplus[i].clear();
0117     fhotmin[i].clear();
0118     energyplus[i].clear();
0119     energymin[i].clear();
0120    }
0121    */
0122
0123   //cout << "declared everything" << endl;
0124
0125   // start particle loop
0126   for (size_t i = 0; i < genParticles->size(); ++i) {
0127     const Candidate& p = (*genParticles)[i];
0128
0129     // select particles in castor
0130     if (fabs(p.eta()) > 5.2 && fabs(p.eta()) < 6.6) {
0131       //cout << "found particle in castor, start calculating" << endl;
0132
0133       // declare energies
0134       double gaus_E = -1.;
0135       double emEnergy = 0.;
0136       double hadEnergy = 0.;
0137       //double fhot = 0.;
0138       //double depth = 0.;
0139
0140       // add energies - em: if particle is e- or gamma
0141       if (p.pdgId() == 11 || p.pdgId() == 22) {
0142         while (gaus_E < 0.) {
0143           // apply energy smearing with gaussian random generator
0144           TRandom3 r(0);
0145           // use sigma/E parametrization for the EM sections of CASTOR TB 2007 results
0146           double sigma = p.energy() * (sqrt(pow(0.044, 2) + pow(0.513 / sqrt(p.energy()), 2)));
0147           gaus_E = r.Gaus(p.energy(), sigma);
0148         }
0149
0150         // calculate electromagnetic electron/photon energy leakage
0151         double tmax;
0152         double a;
0153         double cte;
0154         if (p.pdgId() == 11) {
0155           cte = -0.5;
0156         } else {
0157           cte = 0.5;
0158         }
0159         tmax = 1.0 * (log(gaus_E / 0.0015) + cte);
0160         a = tmax * 0.5 + 1;
0161         double leakage;
0162         double x = 0.5 * 19.38;
0163         leakage = gaus_E - gaus_E * Gamma(a, x);
0164
0165         // add emEnergy
0166         emEnergy = gaus_E - leakage;
0167         // add hadEnergy leakage
0168         hadEnergy = leakage;
0169
0170         // make cluster
0171         ClusterPoint pt1;
0172         if (p.eta() > 0.) {
0173           ClusterPoint temp(88.5, 5.9, p.phi());
0174           pt1 = temp;
0175         }
0176         if (p.eta() < 0.) {
0177           ClusterPoint temp(88.5, -5.9, p.phi());
0178           pt1 = temp;
0179         }
0180         Point pt2(pt1);
0181         CastorTowerRefVector refvector;
0182         CastorClusters->push_back(
0183             reco::CastorCluster(gaus_E, pt2, emEnergy, hadEnergy, emEnergy / gaus_E, 0., 0., 0., 0., refvector));
0184
0185       } else {
0186         while (gaus_E < 0.) {
0187           // apply energy smearing with gaussian random generator
0188           TRandom3 r(0);
0189           // use sigma/E parametrization for the HAD sections of CASTOR TB 2007 results
0190           double sigma = p.energy() * (sqrt(pow(0.121, 2) + pow(1.684 / sqrt(p.energy()), 2)));
0191           gaus_E = r.Gaus(p.energy(), sigma);
0192         }
0193
0194         // add hadEnergy
0195         hadEnergy = gaus_E;
0196
0197         // make cluster
0198         ClusterPoint pt1;
0199         if (p.eta() > 0.) {
0200           ClusterPoint temp(88.5, 5.9, p.phi());
0201           pt1 = temp;
0202         }
0203         if (p.eta() < 0.) {
0204           ClusterPoint temp(88.5, -5.9, p.phi());
0205           pt1 = temp;
0206         }
0207         Point pt2(pt1);
0208         CastorTowerRefVector refvector;
0209         CastorClusters->push_back(reco::CastorCluster(gaus_E, pt2, 0., hadEnergy, 0., 0., 0., 0., 0., refvector));
0210       }
0211
0212       /*
0213         // make tower
0214
0215         // set sector
0216         int sector = -1;
0217         for (int j = 0; j < 16; j++) {
0218             double a = -M_PI + j*0.3927;
0219         double b = -M_PI + (j+1)*0.3927;
0220             if ( (p.phi() > a) && (p.phi() < b)) {
0221            sector = j;
0222         }
0223         }
0224
0225         // set eta
0226         if (p.eta() > 0) {
0227         castorplus[0][sector] = castorplus[0][sector] + gaus_E;
0228         castorplus[1][sector] = castorplus[1][sector] + emEnergy;
0229         castorplus[2][sector] = castorplus[2][sector] + hadEnergy;
0230
0231         depthplus[sector].push_back(depth);
0232         fhotplus[sector].push_back(fhot);
0233         energyplus[sector].push_back(gaus_E);
0234         //cout << "filled vectors" << endl;
0235         //cout << "energyplus size = " << energyplus[sector].size() << endl;
0236         //cout << "depthplus size = " << depthplus[sector].size() << endl;
0237         //cout << "fhotplus size = " << fhotplus[sector].size() << endl;
0238
0239         } else {
0240         castormin[0][sector] = castormin[0][sector] + gaus_E;
0241         castormin[1][sector] = castormin[1][sector] + emEnergy;
0242         castormin[2][sector] = castormin[2][sector] + hadEnergy;
0243
0244
0245         depthmin[sector].push_back(depth);
0246         fhotmin[sector].push_back(fhot);
0247         energymin[sector].push_back(gaus_E);
0248         //cout << "filled vectors" << endl;
0249
0250         }
0251         */
0252     }
0253   }
0254
0255   /*
0256    // substract pedestals/noise
0257    for (int j = 0; j < 16; j++) {
0258     double hadnoise = 0.;
0259     for (int i=0;i<12;i++) {
0260         hadnoise = hadnoise + make_noise();
0261     }
0262     castorplus[0][j] = castorplus[0][j] - hadnoise - make_noise() - make_noise();
0263     castormin[0][j] = castormin[0][j] - hadnoise - make_noise() - make_noise();
0264     castorplus[1][j] = castorplus[1][j] - make_noise() - make_noise();
0265     castormin[1][j] = castormin[1][j] - make_noise() - make_noise();
0266     castorplus[2][j] = castorplus[2][j] - hadnoise;
0267     castormin[2][j] = castormin[2][j] - hadnoise;
0268
0269     // set possible negative values to zero
0270     if (castorplus[0][j] < 0.) castorplus[0][j] = 0.;
0271     if (castormin[0][j] < 0.) castormin[0][j] = 0.;
0272     if (castorplus[1][j] < 0.) castorplus[1][j] = 0.;
0273     if (castormin[1][j] < 0.) castormin[1][j] = 0.;
0274     if (castorplus[2][j] < 0.) castorplus[2][j] = 0.;
0275     if (castormin[2][j] < 0.) castormin[2][j] = 0.;
0276    }
0277    */
0278
0279   /*
0280    // store towers from castor arrays
0281    // eta = 5.9
0282    for (int j=0;j<16;j++) {
0283     if (castorplus[0][j] > 0.) {
0284
0285         double fem = 0.;
0286         fem = castorplus[1][j]/castorplus[0][j];
0287         ClusterPoint pt1(88.5,5.9,castorplus[3][j]);
0288         Point pt2(pt1);
0289
0290         // parametrize depth and fhot from full sim
0291         // get fit parameters from energy
0292         // get random number according to distribution with fit parameters
0293         double depth_mean = 0.;
0294         double fhot_mean = 0.;
0295         double sum_energy = 0.;
0296
0297         //cout << "energyplus size = " << energyplus[j].size()<< endl;
0298         for (size_t p = 0; p<energyplus[j].size();p++) {
0299             depth_mean = depth_mean + depthplus[j][p]*energyplus[j][p];
0300             fhot_mean = fhot_mean + fhotplus[j][p]*energyplus[j][p];
0301         sum_energy = sum_energy + energyplus[j][p];
0302         }
0303         depth_mean = depth_mean/sum_energy;
0304         fhot_mean = fhot_mean/sum_energy;
0305         cout << "computed depth/fhot" << endl;
0306
0307
0308         edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
0309         CastorClusters->push_back(reco::CastorCluster(castorplus[0][j],pt2,castorplus[1][j],castorplus[2][j],fem,depth_mean,fhot_mean,refvector));
0310     }
0311    }
0312    // eta = -5.9
0313    for (int j=0;j<16;j++) {
0314     if (castormin[0][j] > 0.) {
0315         double fem = 0.;
0316         fem = castormin[1][j]/castormin[0][j];
0317         ClusterPoint pt1(88.5,-5.9,castormin[3][j]);
0318         Point pt2(pt1);
0319
0320         // parametrize depth and fhot from full sim
0321         // get fit parameters from energy
0322         // get random number according to distribution with fit parameters
0323         double depth_mean = 0.;
0324         double fhot_mean = 0.;
0325         double sum_energy = 0.;
0326
0327
0328         for (size_t p = 0; p<energymin[j].size();p++) {
0329             depth_mean = depth_mean + depthmin[j][p]*energymin[j][p];
0330             fhot_mean = fhot_mean + fhotmin[j][p]*energymin[j][p];
0331         sum_energy = sum_energy + energymin[j][p];
0332         }
0333         depth_mean = depth_mean/sum_energy;
0334         fhot_mean = fhot_mean/sum_energy;
0335
0336
0337         edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
0338         CastorClusters->push_back(reco::CastorCluster(castormin[0][j],pt2,castormin[1][j],castormin[2][j],fem,depth_mean,fhot_mean,refvector));
0339     }
0340    }
0341    */
0342
0343   iEvent.put(std::move(CastorClusters));
0344 }
0345
0346 double CastorFastClusterProducer::make_noise() {
0347   double result = -1.;
0348   TRandom3 r2(0);
0349   double mu_noise = 0.053;     // GeV (from 1.214 ADC) per channel
0350   double sigma_noise = 0.027;  // GeV (from 0.6168 ADC) per channel
0351
0352   while (result < 0.) {
0353     result = r2.Gaus(mu_noise, sigma_noise);
0354   }
0355
0356   return result;
0357 }
0358
0359 //define this as a plug-in
0360 DEFINE_FWK_MODULE(CastorFastClusterProducer);