![]() |
|
|||
File indexing completed on 2024-04-06 12:11:16
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);
[ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
This page was automatically generated by the 2.2.1 LXR engine. The LXR team |
![]() ![]() |