|
|
|
00001 // -*- C++ -*- 00002 // 00003 // Package: CastorFastClusterProducer 00004 // Class: CastorFastClusterProducer 00005 // 00013 // 00014 // Original Author: Hans Van Haevermaet 00015 // Created: Thu Mar 13 12:00:56 CET 2008 00016 // $Id: CastorFastClusterProducer.cc,v 1.2 2011/03/04 11:00:55 hvanhaev Exp $ 00017 // 00018 // 00019 00020 00021 // system include files 00022 #include <memory> 00023 #include <vector> 00024 #include <iostream> 00025 #include <sstream> 00026 #include <TMath.h> 00027 #include <TRandom3.h> 00028 #include <TF1.h> 00029 00030 // user include files 00031 #include "FWCore/Framework/interface/Frameworkfwd.h" 00032 #include "FWCore/Framework/interface/EDProducer.h" 00033 00034 #include "FWCore/Framework/interface/Event.h" 00035 #include "FWCore/Framework/interface/MakerMacros.h" 00036 00037 #include "FWCore/ParameterSet/interface/ParameterSet.h" 00038 00039 #include "DataFormats/Math/interface/Point3D.h" 00040 00041 // Castorobject includes 00042 #include "DataFormats/CastorReco/interface/CastorCluster.h" 00043 #include "DataFormats/HcalRecHit/interface/CastorRecHit.h" 00044 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h" 00045 00046 // genCandidate particle includes 00047 #include "DataFormats/Candidate/interface/Candidate.h" 00048 #include "DataFormats/HepMCCandidate/interface/GenParticle.h" 00049 00050 #include "FastSimulation/ForwardDetectors/plugins/CastorFastClusterProducer.h" 00051 00052 00053 // 00054 // constructors and destructor 00055 // 00056 CastorFastClusterProducer::CastorFastClusterProducer(const edm::ParameterSet& iConfig) 00057 { 00058 //register your products 00059 produces<CastorClusterCollection>(); 00060 00061 //now do what ever other initialization is needed 00062 00063 } 00064 00065 00066 CastorFastClusterProducer::~CastorFastClusterProducer() 00067 { 00068 00069 // do anything here that needs to be done at desctruction time 00070 // (e.g. close files, deallocate resources etc.) 00071 00072 } 00073 00074 00075 // 00076 // member functions 00077 // 00078 00079 // ------------ method called to produce the data ------------ 00080 void 00081 CastorFastClusterProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) 00082 { 00083 using namespace edm; 00084 using namespace reco; 00085 using namespace std; 00086 using namespace TMath; 00087 00088 // 00089 // Make CastorCluster objects 00090 // 00091 00092 //cout << "entering event" << endl; 00093 00094 Handle<GenParticleCollection> genParticles; 00095 iEvent.getByLabel("genParticles", genParticles); 00096 00097 // make pointer to towers that will be made 00098 auto_ptr<CastorClusterCollection> CastorClusters (new CastorClusterCollection); 00099 00100 /* 00101 // declare castor array 00102 double castorplus [4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = 5.9 00103 double castormin [4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = -5.9 00104 // set phi values of array sectors and everything else to zero 00105 for (int j = 0; j < 16; j++) { 00106 castorplus[3][j] = -2.94524 + j*0.3927; 00107 castormin[3][j] = -2.94524 + j*0.3927; 00108 castorplus[0][j] = 0.; 00109 castormin[0][j] = 0.; 00110 castorplus[1][j] = 0.; 00111 castormin[1][j] = 0.; 00112 castorplus[2][j] = 0.; 00113 castormin[2][j] = 0.; 00114 //castorplus[4][j] = 0.; 00115 //castormin[4][j] = 0.; 00116 } 00117 00118 // declare properties vectors 00119 vector<double> depthplus[16]; 00120 vector<double> depthmin[16]; 00121 vector<double> fhotplus [16]; 00122 vector<double> fhotmin [16]; 00123 vector<double> energyplus [16]; 00124 vector<double> energymin [16]; 00125 00126 for (int i=0;i<16;i++) { 00127 depthplus[i].clear(); 00128 depthmin[i].clear(); 00129 fhotplus[i].clear(); 00130 fhotmin[i].clear(); 00131 energyplus[i].clear(); 00132 energymin[i].clear(); 00133 } 00134 */ 00135 00136 //cout << "declared everything" << endl; 00137 00138 // start particle loop 00139 for (size_t i = 0; i < genParticles->size(); ++i) { 00140 const Candidate & p = (*genParticles)[i]; 00141 00142 // select particles in castor 00143 if ( fabs(p.eta()) > 5.2 && fabs(p.eta()) < 6.6) { 00144 00145 //cout << "found particle in castor, start calculating" << endl; 00146 00147 // declare energies 00148 double gaus_E = -1.; 00149 double emEnergy = 0.; 00150 double hadEnergy = 0.; 00151 //double fhot = 0.; 00152 //double depth = 0.; 00153 00154 // add energies - em: if particle is e- or gamma 00155 if (p.pdgId() == 11 || p.pdgId() == 22) { 00156 00157 while ( gaus_E < 0.) { 00158 // apply energy smearing with gaussian random generator 00159 TRandom3 r(0); 00160 // use sigma/E parametrization for the EM sections of CASTOR TB 2007 results 00161 double sigma = p.energy()*(sqrt(pow(0.044,2) + pow(0.513/sqrt(p.energy()),2))); 00162 gaus_E = r.Gaus(p.energy(),sigma); 00163 } 00164 00165 // calculate electromagnetic electron/photon energy leakage 00166 double tmax; 00167 double a; 00168 double cte; 00169 if ( p.pdgId() == 11) { cte = -0.5; } else { cte = 0.5; } 00170 tmax = 1.0*(log(gaus_E/0.0015)+cte); 00171 a = tmax*0.5 + 1; 00172 double leakage; 00173 double x = 0.5*19.38; 00174 leakage = gaus_E - gaus_E*Gamma(a,x); 00175 00176 // add emEnergy 00177 emEnergy = gaus_E - leakage; 00178 // add hadEnergy leakage 00179 hadEnergy = leakage; 00180 00181 // make cluster 00182 ClusterPoint pt1; 00183 if (p.eta() > 0.) {ClusterPoint temp(88.5,5.9,p.phi()); pt1 = temp;} 00184 if (p.eta() < 0.) {ClusterPoint temp(88.5,-5.9,p.phi()); pt1 = temp;} 00185 Point pt2(pt1); 00186 CastorTowerRefVector refvector; 00187 CastorClusters->push_back(reco::CastorCluster(gaus_E,pt2,emEnergy,hadEnergy,emEnergy/gaus_E,0.,0.,0.,0.,refvector)); 00188 00189 } else { 00190 00191 while (gaus_E < 0.) { 00192 // apply energy smearing with gaussian random generator 00193 TRandom3 r(0); 00194 // use sigma/E parametrization for the HAD sections of CASTOR TB 2007 results 00195 double sigma = p.energy()*(sqrt(pow(0.121,2) + pow(1.684/sqrt(p.energy()),2))); 00196 gaus_E = r.Gaus(p.energy(),sigma); 00197 } 00198 00199 // add hadEnergy 00200 hadEnergy = gaus_E; 00201 00202 // make cluster 00203 ClusterPoint pt1; 00204 if (p.eta() > 0.) {ClusterPoint temp(88.5,5.9,p.phi()); pt1 = temp;} 00205 if (p.eta() < 0.) {ClusterPoint temp(88.5,-5.9,p.phi()); pt1 = temp;} 00206 Point pt2(pt1); 00207 CastorTowerRefVector refvector; 00208 CastorClusters->push_back(reco::CastorCluster(gaus_E,pt2,0.,hadEnergy,0.,0.,0.,0.,0.,refvector)); 00209 } 00210 00211 /* 00212 // make tower 00213 00214 // set sector 00215 int sector = -1; 00216 for (int j = 0; j < 16; j++) { 00217 double a = -M_PI + j*0.3927; 00218 double b = -M_PI + (j+1)*0.3927; 00219 if ( (p.phi() > a) && (p.phi() < b)) { 00220 sector = j; 00221 } 00222 } 00223 00224 // set eta 00225 if (p.eta() > 0) { 00226 castorplus[0][sector] = castorplus[0][sector] + gaus_E; 00227 castorplus[1][sector] = castorplus[1][sector] + emEnergy; 00228 castorplus[2][sector] = castorplus[2][sector] + hadEnergy; 00229 00230 depthplus[sector].push_back(depth); 00231 fhotplus[sector].push_back(fhot); 00232 energyplus[sector].push_back(gaus_E); 00233 //cout << "filled vectors" << endl; 00234 //cout << "energyplus size = " << energyplus[sector].size() << endl; 00235 //cout << "depthplus size = " << depthplus[sector].size() << endl; 00236 //cout << "fhotplus size = " << fhotplus[sector].size() << endl; 00237 00238 } else { 00239 castormin[0][sector] = castormin[0][sector] + gaus_E; 00240 castormin[1][sector] = castormin[1][sector] + emEnergy; 00241 castormin[2][sector] = castormin[2][sector] + hadEnergy; 00242 00243 00244 depthmin[sector].push_back(depth); 00245 fhotmin[sector].push_back(fhot); 00246 energymin[sector].push_back(gaus_E); 00247 //cout << "filled vectors" << endl; 00248 00249 } 00250 */ 00251 00252 } 00253 00254 } 00255 00256 /* 00257 // substract pedestals/noise 00258 for (int j = 0; j < 16; j++) { 00259 double hadnoise = 0.; 00260 for (int i=0;i<12;i++) { 00261 hadnoise = hadnoise + make_noise(); 00262 } 00263 castorplus[0][j] = castorplus[0][j] - hadnoise - make_noise() - make_noise(); 00264 castormin[0][j] = castormin[0][j] - hadnoise - make_noise() - make_noise(); 00265 castorplus[1][j] = castorplus[1][j] - make_noise() - make_noise(); 00266 castormin[1][j] = castormin[1][j] - make_noise() - make_noise(); 00267 castorplus[2][j] = castorplus[2][j] - hadnoise; 00268 castormin[2][j] = castormin[2][j] - hadnoise; 00269 00270 // set possible negative values to zero 00271 if (castorplus[0][j] < 0.) castorplus[0][j] = 0.; 00272 if (castormin[0][j] < 0.) castormin[0][j] = 0.; 00273 if (castorplus[1][j] < 0.) castorplus[1][j] = 0.; 00274 if (castormin[1][j] < 0.) castormin[1][j] = 0.; 00275 if (castorplus[2][j] < 0.) castorplus[2][j] = 0.; 00276 if (castormin[2][j] < 0.) castormin[2][j] = 0.; 00277 } 00278 */ 00279 00280 /* 00281 // store towers from castor arrays 00282 // eta = 5.9 00283 for (int j=0;j<16;j++) { 00284 if (castorplus[0][j] > 0.) { 00285 00286 double fem = 0.; 00287 fem = castorplus[1][j]/castorplus[0][j]; 00288 ClusterPoint pt1(88.5,5.9,castorplus[3][j]); 00289 Point pt2(pt1); 00290 00291 // parametrize depth and fhot from full sim 00292 // get fit parameters from energy 00293 // get random number according to distribution with fit parameters 00294 double depth_mean = 0.; 00295 double fhot_mean = 0.; 00296 double sum_energy = 0.; 00297 00298 //cout << "energyplus size = " << energyplus[j].size()<< endl; 00299 for (size_t p = 0; p<energyplus[j].size();p++) { 00300 depth_mean = depth_mean + depthplus[j][p]*energyplus[j][p]; 00301 fhot_mean = fhot_mean + fhotplus[j][p]*energyplus[j][p]; 00302 sum_energy = sum_energy + energyplus[j][p]; 00303 } 00304 depth_mean = depth_mean/sum_energy; 00305 fhot_mean = fhot_mean/sum_energy; 00306 cout << "computed depth/fhot" << endl; 00307 00308 00309 edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector; 00310 CastorClusters->push_back(reco::CastorCluster(castorplus[0][j],pt2,castorplus[1][j],castorplus[2][j],fem,depth_mean,fhot_mean,refvector)); 00311 } 00312 } 00313 // eta = -5.9 00314 for (int j=0;j<16;j++) { 00315 if (castormin[0][j] > 0.) { 00316 double fem = 0.; 00317 fem = castormin[1][j]/castormin[0][j]; 00318 ClusterPoint pt1(88.5,-5.9,castormin[3][j]); 00319 Point pt2(pt1); 00320 00321 // parametrize depth and fhot from full sim 00322 // get fit parameters from energy 00323 // get random number according to distribution with fit parameters 00324 double depth_mean = 0.; 00325 double fhot_mean = 0.; 00326 double sum_energy = 0.; 00327 00328 00329 for (size_t p = 0; p<energymin[j].size();p++) { 00330 depth_mean = depth_mean + depthmin[j][p]*energymin[j][p]; 00331 fhot_mean = fhot_mean + fhotmin[j][p]*energymin[j][p]; 00332 sum_energy = sum_energy + energymin[j][p]; 00333 } 00334 depth_mean = depth_mean/sum_energy; 00335 fhot_mean = fhot_mean/sum_energy; 00336 00337 00338 edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector; 00339 CastorClusters->push_back(reco::CastorCluster(castormin[0][j],pt2,castormin[1][j],castormin[2][j],fem,depth_mean,fhot_mean,refvector)); 00340 } 00341 } 00342 */ 00343 00344 iEvent.put(CastorClusters); 00345 00346 } 00347 00348 double CastorFastClusterProducer::make_noise() { 00349 double result = -1.; 00350 TRandom3 r2(0); 00351 double mu_noise = 0.053; // GeV (from 1.214 ADC) per channel 00352 double sigma_noise = 0.027; // GeV (from 0.6168 ADC) per channel 00353 00354 while (result < 0.) { 00355 result = r2.Gaus(mu_noise,sigma_noise); 00356 } 00357 00358 return result; 00359 } 00360 00361 00362 // ------------ method called once each job just before starting event loop ------------ 00363 void 00364 CastorFastClusterProducer::beginRun(edm::Run& run, edm::EventSetup const& es) 00365 { 00366 } 00367 00368 // ------------ method called once each job just after ending the event loop ------------ 00369 void 00370 CastorFastClusterProducer::endRun() { 00371 } 00372 00373 //define this as a plug-in 00374 DEFINE_FWK_MODULE(CastorFastClusterProducer);
1.7.3