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