File indexing completed on 2024-09-10 02:59:14
0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003
0004 #include "SimTransport/HectorProducer/interface/Hector.h"
0005
0006 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0007 #include <CLHEP/Units/SystemOfUnits.h>
0008 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0009 #include "HepMC/SimpleVector.h"
0010
0011 #include "TRandom3.h"
0012
0013 #include "H_Parameters.h"
0014
0015 #include <cmath>
0016
0017 Hector::Hector(const edm::ParameterSet& param,
0018 const edm::ESGetToken<HepPDT::ParticleDataTable, PDTRecord>& token,
0019 bool verbosity,
0020 bool FP420Transport,
0021 bool ZDCTransport)
0022 : tok_pdt_(token), m_verbosity(verbosity), m_FP420Transport(FP420Transport), m_ZDCTransport(ZDCTransport) {
0023
0024 edm::ParameterSet hector_par = param.getParameter<edm::ParameterSet>("Hector");
0025
0026
0027 lengthfp420 = hector_par.getParameter<double>("BeamLineLengthFP420");
0028 m_rpp420_f = (float)hector_par.getParameter<double>("RP420f");
0029 m_rpp420_b = (float)hector_par.getParameter<double>("RP420b");
0030 lengthzdc = hector_par.getParameter<double>("BeamLineLengthZDC");
0031 lengthd1 = hector_par.getParameter<double>("BeamLineLengthD1");
0032 beam1filename = hector_par.getParameter<string>("Beam1");
0033 beam2filename = hector_par.getParameter<string>("Beam2");
0034 m_rppzdc = (float)lengthzdc;
0035 m_rppd1 = (float)lengthd1;
0036 m_smearAng = hector_par.getParameter<bool>("smearAng");
0037 m_sigmaSTX = hector_par.getParameter<double>("sigmaSTX");
0038 m_sigmaSTY = hector_par.getParameter<double>("sigmaSTY");
0039 m_smearE = hector_par.getParameter<bool>("smearEnergy");
0040 m_sig_e = hector_par.getParameter<double>("sigmaEnergy");
0041 etacut = hector_par.getParameter<double>("EtaCutForHector");
0042
0043 theCorrespondenceMap.clear();
0044
0045 if (m_verbosity) {
0046 edm::LogInfo("HectorSetup") << "===================================================================\n"
0047 << " * * * * * * * * * * * * * * * * * * * * * * * * * * * * \n"
0048 << " * * \n"
0049 << " * --<--<-- A fast simulator --<--<-- * \n"
0050 << " * | --<--<-- of particle --<--<-- * \n"
0051 << " * ----HECTOR----< * \n"
0052 << " * | -->-->-- transport through-->-->-- * \n"
0053 << " * -->-->-- generic beamlines -->-->-- * \n"
0054 << " * * \n"
0055 << " * JINST 2:P09005 (2007) * \n"
0056 << " * X Rouby, J de Favereau, K Piotrzkowski (CP3) * \n"
0057 << " * http://www.fynu.ucl.ac.be/hector.html * \n"
0058 << " * * \n"
0059 << " * Center for Cosmology, Particle Physics and Phenomenology * \n"
0060 << " * Universite catholique de Louvain * \n"
0061 << " * Louvain-la-Neuve, Belgium * \n"
0062 << " * * \n"
0063 << " * * * * * * * * * * * * * * * * * * * * * * * * * * * * \n"
0064 << " Hector configuration: \n"
0065 << " m_FP420Transport = " << m_FP420Transport << "\n"
0066 << " m_ZDCTransport = " << m_ZDCTransport << "\n"
0067 << " lengthfp420 = " << lengthfp420 << "\n"
0068 << " m_rpp420_f = " << m_rpp420_f << "\n"
0069 << " m_rpp420_b = " << m_rpp420_b << "\n"
0070 << " lengthzdc = " << lengthzdc << "\n"
0071 << " lengthd1 = " << lengthd1 << "\n"
0072 << "===================================================================\n";
0073 }
0074 edm::FileInPath b1(beam1filename.c_str());
0075 edm::FileInPath b2(beam2filename.c_str());
0076
0077
0078 if (m_FP420Transport && lengthfp420 > 0.) {
0079 m_beamlineFP4201 = new H_BeamLine(1, lengthfp420 + 0.1);
0080 m_beamlineFP4202 = new H_BeamLine(-1, lengthfp420 + 0.1);
0081 m_beamlineFP4201->fill(b1.fullPath(), 1, "IP5");
0082 m_beamlineFP4202->fill(b2.fullPath(), -1, "IP5");
0083 m_beamlineFP4201->offsetElements(120, -0.097);
0084 m_beamlineFP4202->offsetElements(120, +0.097);
0085 m_beamlineFP4201->calcMatrix();
0086 m_beamlineFP4202->calcMatrix();
0087 } else {
0088 if (m_verbosity)
0089 edm::LogVerbatim("HectorSetup") << "Hector: WARNING: lengthfp420= " << lengthfp420;
0090 }
0091
0092 if (m_ZDCTransport && lengthzdc > 0. && lengthd1 > 0.) {
0093
0094 m_beamlineZDC1 = new H_BeamLine(1, lengthzdc + 0.1);
0095 m_beamlineZDC2 = new H_BeamLine(-1, lengthzdc + 0.1);
0096 m_beamlineZDC1->fill(b1.fullPath(), 1, "IP5");
0097 m_beamlineZDC2->fill(b2.fullPath(), -1, "IP5");
0098 m_beamlineZDC1->offsetElements(120, -0.097);
0099 m_beamlineZDC2->offsetElements(120, +0.097);
0100 m_beamlineZDC1->calcMatrix();
0101 m_beamlineZDC2->calcMatrix();
0102
0103
0104 m_beamlineD11 = new H_BeamLine(1, lengthd1 + 0.1);
0105 m_beamlineD12 = new H_BeamLine(-1, lengthd1 + 0.1);
0106 m_beamlineD11->fill(b1.fullPath(), 1, "IP5");
0107 m_beamlineD12->fill(b2.fullPath(), -1, "IP5");
0108 m_beamlineD11->offsetElements(120, -0.097);
0109 m_beamlineD12->offsetElements(120, +0.097);
0110 m_beamlineD11->calcMatrix();
0111 m_beamlineD12->calcMatrix();
0112 } else {
0113 if (m_verbosity)
0114 edm::LogVerbatim("HectorSetup") << "Hector: WARNING: lengthzdc= " << lengthzdc << "lengthd1= " << lengthd1;
0115 }
0116 }
0117
0118 Hector::~Hector() {
0119 for (std::map<unsigned int, H_BeamParticle*>::iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
0120 delete (*it).second;
0121 }
0122
0123 delete m_beamlineFP4201;
0124 delete m_beamlineFP4202;
0125 delete m_beamlineZDC1;
0126 delete m_beamlineZDC2;
0127 delete m_beamlineD11;
0128 delete m_beamlineD12;
0129 }
0130
0131 void Hector::clearApertureFlags() {
0132 m_isStoppedfp420.clear();
0133 m_isStoppedzdc.clear();
0134 m_isStoppedd1.clear();
0135 }
0136
0137 void Hector::clear() {
0138 for (std::map<unsigned int, H_BeamParticle*>::iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
0139 delete (*it).second;
0140 };
0141 m_beamPart.clear();
0142 m_direct.clear();
0143 m_eta.clear();
0144 m_pdg.clear();
0145 m_pz.clear();
0146 m_isCharged.clear();
0147 }
0148
0149 void Hector::add(const HepMC::GenEvent* evt, const edm::EventSetup& iSetup) {
0150 H_BeamParticle* h_p;
0151 unsigned int line;
0152
0153 for (HepMC::GenEvent::particle_const_iterator eventParticle = evt->particles_begin();
0154 eventParticle != evt->particles_end();
0155 ++eventParticle) {
0156 if ((*eventParticle)->status() == 1) {
0157 if (abs((*eventParticle)->momentum().eta()) > etacut) {
0158 line = (*eventParticle)->barcode();
0159 if (m_beamPart.find(line) == m_beamPart.end()) {
0160 double charge = 1.;
0161 m_isCharged[line] = false;
0162 HepMC::GenParticle* g = (*eventParticle);
0163 pdt = iSetup.getHandle(tok_pdt_);
0164 const ParticleData* part = pdt->particle(g->pdg_id());
0165 if (part) {
0166 charge = part->charge();
0167 }
0168 if (charge != 0)
0169 m_isCharged[line] = true;
0170 double mass = (*eventParticle)->generatedMass();
0171
0172 h_p = new H_BeamParticle(mass, charge);
0173
0174 double px, py, pz;
0175 px = (*eventParticle)->momentum().px();
0176 py = (*eventParticle)->momentum().py();
0177 pz = (*eventParticle)->momentum().pz();
0178
0179 h_p->set4Momentum(px, py, pz, (*eventParticle)->momentum().e());
0180
0181
0182 double XforPosition = (*eventParticle)->production_vertex()->position().x() / CLHEP::micrometer;
0183 double YforPosition = (*eventParticle)->production_vertex()->position().y() / CLHEP::micrometer;
0184 double ZforPosition = (*eventParticle)->production_vertex()->position().z() / CLHEP::meter;
0185
0186 double TXforPosition = 0.0, TYforPosition = 0.0;
0187
0188
0189 h_p->setPosition(XforPosition, YforPosition, TXforPosition, TYforPosition, ZforPosition);
0190
0191 m_beamPart[line] = h_p;
0192 m_direct[line] = 0;
0193 m_direct[line] = (pz > 0) ? 1 : -1;
0194
0195 m_eta[line] = (*eventParticle)->momentum().eta();
0196 m_pdg[line] = (*eventParticle)->pdg_id();
0197 m_pz[line] = (*eventParticle)->momentum().pz();
0198
0199 if (m_verbosity) {
0200 edm::LogVerbatim("HectorEventProcessing")
0201 << "Hector:add: barcode = " << line << " status = " << g->status() << " PDG Id = " << g->pdg_id()
0202 << " mass = " << mass << " pz = " << pz << " charge = " << charge
0203 << " m_isCharged[line] = " << m_isCharged[line];
0204 }
0205 }
0206 }
0207 }
0208 }
0209 }
0210
0211 void Hector::filterFP420(TRandom3* rootEngine) {
0212 unsigned int line;
0213 H_BeamParticle* part;
0214 std::map<unsigned int, H_BeamParticle*>::iterator it;
0215
0216 bool is_stop;
0217 int direction;
0218
0219 float x1_420;
0220 float y1_420;
0221
0222 if (!m_beamPart.empty() && lengthfp420 > 0.) {
0223 for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
0224 line = (*it).first;
0225 part = (*it).second;
0226
0227 if (m_verbosity) {
0228 edm::LogVerbatim("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line;
0229 }
0230 if ((*m_isCharged.find(line)).second) {
0231 direction = (*m_direct.find(line)).second;
0232 if (m_smearAng) {
0233
0234 if (m_sigmaSTX > 0. && m_sigmaSTY > 0.) {
0235 part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0236 } else {
0237
0238 part->smearAng(STX, STY, rootEngine);
0239 }
0240 }
0241 if (m_smearE) {
0242 if (m_sig_e) {
0243 part->smearE(m_sig_e, rootEngine);
0244 } else {
0245 part->smearE(SBE, rootEngine);
0246 }
0247 }
0248 if (direction == 1 && m_beamlineFP4201 != nullptr) {
0249 part->computePath(m_beamlineFP4201);
0250 is_stop = part->stopped(m_beamlineFP4201);
0251 if (m_verbosity)
0252 edm::LogVerbatim("HectorEventProcessing")
0253 << "Hector:filterFP420: barcode = " << line << " positive is_stop= " << is_stop;
0254 } else if (direction == -1 && m_beamlineFP4202 != nullptr) {
0255 part->computePath(m_beamlineFP4202);
0256 is_stop = part->stopped(m_beamlineFP4202);
0257 if (m_verbosity)
0258 edm::LogVerbatim("HectorEventProcessing")
0259 << "Hector:filterFP420: barcode = " << line << " negative is_stop= " << is_stop;
0260 } else {
0261 is_stop = true;
0262 if (m_verbosity)
0263 edm::LogVerbatim("HectorEventProcessing")
0264 << "Hector:filterFP420: barcode = " << line << " 0 is_stop= " << is_stop;
0265 }
0266
0267
0268 m_isStoppedfp420[line] = is_stop;
0269 if (m_verbosity)
0270 edm::LogVerbatim("HectorEventProcessing")
0271 << "Hector:filterFP420: barcode = " << line << " isStopped=" << (*m_isStoppedfp420.find(line)).second;
0272
0273 if (!is_stop) {
0274 if (direction == 1)
0275 part->propagate(m_rpp420_f);
0276 if (direction == -1)
0277 part->propagate(m_rpp420_b);
0278 x1_420 = part->getX();
0279 y1_420 = part->getY();
0280 if (m_verbosity)
0281 edm::LogVerbatim("HectorEventProcessing")
0282 << "Hector:filterFP420: barcode = " << line << " x= " << x1_420 << " y= " << y1_420;
0283
0284 m_xAtTrPoint[line] = x1_420;
0285 m_yAtTrPoint[line] = y1_420;
0286 m_TxAtTrPoint[line] = part->getTX();
0287 m_TyAtTrPoint[line] = part->getTY();
0288 m_eAtTrPoint[line] = part->getE();
0289 }
0290 }
0291 else {
0292 m_isStoppedfp420[line] = true;
0293 if (m_verbosity)
0294 edm::LogVerbatim("HectorEventProcessing")
0295 << "Hector:filterFP420: barcode = " << line << " isStopped=" << (*m_isStoppedfp420.find(line)).second;
0296 }
0297
0298 }
0299 }
0300 }
0301
0302 void Hector::filterZDC(TRandom3* rootEngine) {
0303 unsigned int line;
0304 H_BeamParticle* part;
0305 std::map<unsigned int, H_BeamParticle*>::iterator it;
0306
0307 bool is_stop_zdc = false;
0308 int direction;
0309
0310 if (!m_beamPart.empty() && lengthzdc > 0.) {
0311 for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
0312 line = (*it).first;
0313 part = (*it).second;
0314 if (m_verbosity) {
0315 edm::LogVerbatim("HectorEventProcessing")
0316 << "Hector:filterZDC: barcode = " << line << " charge = " << (*m_isCharged.find(line)).second;
0317 if (m_FP420Transport)
0318 edm::LogVerbatim("HectorEventProcessing") << " isStoppedFP420 =" << (*m_isStoppedfp420.find(line)).second;
0319 }
0320
0321 if (((*m_isCharged.find(line)).second)) {
0322 if (m_verbosity)
0323 edm::LogVerbatim("HectorEventProcessing") << "Hector:filterZDC: barcode = " << line << " propagated ";
0324
0325 direction = (*m_direct.find(line)).second;
0326 if (m_verbosity)
0327 edm::LogVerbatim("HectorEventProcessing")
0328 << "Hector:filterZDC: barcode = " << line << " direction = " << direction;
0329 if (m_smearAng) {
0330 if (m_sigmaSTX > 0. && m_sigmaSTY > 0.) {
0331
0332 part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0333 } else {
0334
0335 part->smearAng(STX, STY, rootEngine);
0336 }
0337 }
0338 if (m_smearE) {
0339 if (m_sig_e) {
0340 part->smearE(m_sig_e, rootEngine);
0341 } else {
0342 part->smearE(SBE, rootEngine);
0343 }
0344 }
0345 if (direction == 1 && m_beamlineZDC1 != nullptr) {
0346 part->computePath(m_beamlineZDC1);
0347 is_stop_zdc = part->stopped(m_beamlineZDC1);
0348 m_isStoppedzdc[line] = is_stop_zdc;
0349 if (m_verbosity)
0350 edm::LogVerbatim("HectorEventProcessing")
0351 << "Hector:filterZDC: barcode " << line << " positive is_stop_zdc= " << is_stop_zdc;
0352 } else if (direction == -1 && m_beamlineZDC2 != nullptr) {
0353 part->computePath(m_beamlineZDC2);
0354 is_stop_zdc = part->stopped(m_beamlineZDC2);
0355 m_isStoppedzdc[line] = is_stop_zdc;
0356 if (m_verbosity)
0357 edm::LogVerbatim("HectorEventProcessing")
0358 << "Hector:filterZDC: barcode " << line << " negative is_stop_zdc= " << is_stop_zdc;
0359 } else {
0360 m_isStoppedzdc[line] = true;
0361 if (m_verbosity)
0362 edm::LogVerbatim("HectorEventProcessing")
0363 << "Hector:filterZDC: barcode " << line << " 0 is_stop_zdc= " << is_stop_zdc;
0364 }
0365 }
0366
0367 else {
0368 m_isStoppedzdc[line] = true;
0369 if (m_verbosity)
0370 edm::LogVerbatim("HectorEventProcessing")
0371 << "Hector:filterZDC: barcode = " << line << " isStopped=" << (*m_isStoppedzdc.find(line)).second;
0372 }
0373
0374 }
0375 }
0376 }
0377
0378 void Hector::filterD1(TRandom3* rootEngine) {
0379 unsigned int line;
0380 H_BeamParticle* part;
0381 std::map<unsigned int, H_BeamParticle*>::iterator it;
0382
0383 bool is_stop_d1;
0384 int direction;
0385
0386 float x1_d1;
0387 float y1_d1;
0388
0389 if (!m_beamPart.empty() && lengthd1 > 0.) {
0390 for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
0391 line = (*it).first;
0392 part = (*it).second;
0393 if (m_verbosity)
0394 edm::LogVerbatim("HectorEventProcessing")
0395 << "Hector:filterD1: barcode = " << line << " isStoppedZDC =" << (*m_isStoppedzdc.find(line)).second;
0396 if (((*m_isStoppedzdc.find(line)).second) || !((*m_isCharged.find(line)).second)) {
0397 if (m_verbosity)
0398 edm::LogVerbatim("HectorEventProcessing") << "Hector:filterD1: barcode = " << line << " propagated ";
0399
0400 direction = (*m_direct.find(line)).second;
0401 if (m_verbosity)
0402 edm::LogVerbatim("HectorEventProcessing") << "Hector:filterD1:direction=" << direction;
0403 if (m_smearAng) {
0404 if (m_sigmaSTX > 0. && m_sigmaSTY > 0.) {
0405
0406 part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0407 } else {
0408
0409 part->smearAng(STX, STY, rootEngine);
0410 }
0411 }
0412 if (m_smearE) {
0413 if (m_sig_e) {
0414 part->smearE(m_sig_e, rootEngine);
0415 } else {
0416 part->smearE(SBE, rootEngine);
0417 }
0418 }
0419 if (direction == 1 && m_beamlineD11 != nullptr) {
0420 part->computePath(m_beamlineD11);
0421 is_stop_d1 = part->stopped(m_beamlineD11);
0422 m_isStoppedd1[line] = is_stop_d1;
0423 if (m_verbosity)
0424 edm::LogVerbatim("HectorEventProcessing")
0425 << "Hector:filterD1 barcode " << line << " positive is_stop_d1 = " << is_stop_d1;
0426 } else if (direction == -1 && m_beamlineD12 != nullptr) {
0427 part->computePath(m_beamlineD12);
0428 is_stop_d1 = part->stopped(m_beamlineD12);
0429 m_isStoppedd1[line] = is_stop_d1;
0430 if (m_verbosity)
0431 edm::LogVerbatim("HectorEventProcessing")
0432 << "Hector:filterD1 barcode " << line << " negative is_stop_d1 = " << is_stop_d1;
0433 } else {
0434 is_stop_d1 = true;
0435 m_isStoppedd1[line] = is_stop_d1;
0436 if (m_verbosity)
0437 edm::LogVerbatim("HectorEventProcessing")
0438 << "Hector:filterD1 barcode " << line << " 0 is_stop_d1 = " << is_stop_d1;
0439 }
0440
0441 if (!is_stop_d1) {
0442 part->propagate(lengthd1);
0443 x1_d1 = part->getX();
0444 y1_d1 = part->getY();
0445 m_xAtTrPoint[line] = x1_d1;
0446 m_yAtTrPoint[line] = y1_d1;
0447 m_TxAtTrPoint[line] = part->getTX();
0448 m_TyAtTrPoint[line] = part->getTY();
0449 m_eAtTrPoint[line] = part->getE();
0450 }
0451 }
0452 else {
0453 m_isStoppedd1[line] = false;
0454 if (m_verbosity)
0455 edm::LogVerbatim("HectorEventProcessing")
0456 << "Hector:filterD1: barcode = " << line << " isStopped=" << (*m_isStoppedd1.find(line)).second;
0457 }
0458
0459 }
0460 }
0461 }
0462
0463 int Hector::getDirect(unsigned int part_n) const {
0464 std::map<unsigned int, int>::const_iterator it = m_direct.find(part_n);
0465 if (it != m_direct.end()) {
0466 return (*it).second;
0467 }
0468 return 0;
0469 }
0470
0471 void Hector::print() const {
0472 for (std::map<unsigned int, H_BeamParticle*>::const_iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
0473 (*it).second->printProperties();
0474 };
0475 }
0476
0477 HepMC::GenEvent* Hector::addPartToHepMC(HepMC::GenEvent* evt) {
0478 theCorrespondenceMap.clear();
0479
0480 unsigned int line;
0481
0482 HepMC::GenParticle* gpart;
0483 long double tx, ty, theta, fi, energy, time = 0;
0484 std::map<unsigned int, H_BeamParticle*>::iterator it;
0485
0486 for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
0487 line = (*it).first;
0488 if (!m_FP420Transport)
0489 m_isStoppedfp420[line] = true;
0490 if (!m_ZDCTransport) {
0491 m_isStoppedzdc[line] = false;
0492 m_isStoppedd1[line] = true;
0493 }
0494 if (m_verbosity) {
0495 edm::LogVerbatim("HectorEventProcessing")
0496 << "Hector:addPartToHepMC: barcode = " << line << "\n"
0497 << "Hector:addPartToHepMC: isStoppedfp420=" << (*m_isStoppedfp420.find(line)).second << "\n"
0498 << "Hector:addPartToHepMC: isStoppedzdc=" << (*m_isStoppedzdc.find(line)).second << "\n"
0499 << "Hector:addPartToHepMC: isStoppedd1=" << (*m_isStoppedd1.find(line)).second;
0500 }
0501 if (!((*m_isStoppedfp420.find(line)).second) ||
0502 (!((*m_isStoppedd1.find(line)).second) && ((*m_isStoppedzdc.find(line)).second))) {
0503 gpart = evt->barcode_to_particle(line);
0504 if (gpart) {
0505 tx = (*m_TxAtTrPoint.find(line)).second / 1000000.;
0506 ty = (*m_TyAtTrPoint.find(line)).second / 1000000.;
0507 theta = sqrt((tx * tx) + (ty * ty));
0508 double ddd = 0.;
0509 if (!((*m_isStoppedfp420.find(line)).second)) {
0510 if ((*m_direct.find(line)).second > 0) {
0511 ddd = m_rpp420_f;
0512 } else if ((*m_direct.find(line)).second < 0) {
0513 ddd = m_rpp420_b;
0514 theta = pi - theta;
0515 }
0516 } else {
0517 ddd = lengthd1;
0518 if ((*m_direct.find(line)).second < 0)
0519 theta = pi - theta;
0520 }
0521
0522 fi = std::atan2(tx, ty);
0523 energy = (*m_eAtTrPoint.find(line)).second;
0524
0525 time = (ddd * CLHEP::meter - gpart->production_vertex()->position().z() * CLHEP::mm);
0526
0527 if (ddd != 0.) {
0528 if (m_verbosity) {
0529 edm::LogVerbatim("HectorEventProcessing")
0530 << "Hector:: x= " << (*(m_xAtTrPoint.find(line))).second * 0.001 << "\n"
0531 << "Hector:: y= " << (*(m_yAtTrPoint.find(line))).second * 0.001 << "\n"
0532 << "Hector:: z= " << ddd * (*(m_direct.find(line))).second * 1000. << "\n"
0533 << "Hector:: t= " << time;
0534 }
0535
0536 HepMC::GenVertex* vert = new HepMC::GenVertex(HepMC::FourVector((*(m_xAtTrPoint.find(line))).second * 0.001,
0537 (*(m_yAtTrPoint.find(line))).second * 0.001,
0538 ddd * (*(m_direct.find(line))).second * 1000.,
0539 time + .001 * time));
0540
0541 gpart->set_status(2);
0542 vert->add_particle_in(gpart);
0543 vert->add_particle_out(new HepMC::GenParticle(HepMC::FourVector(energy * std::sin(theta) * std::sin(fi),
0544 energy * std::sin(theta) * std::cos(fi),
0545 energy * std::cos(theta),
0546 energy),
0547 gpart->pdg_id(),
0548 1,
0549 gpart->flow()));
0550 evt->add_vertex(vert);
0551
0552 int ingoing = (*vert->particles_in_const_begin())->barcode();
0553 int outgoing = (*vert->particles_out_const_begin())->barcode();
0554 LHCTransportLink theLink(ingoing, outgoing);
0555 if (m_verbosity)
0556 edm::LogVerbatim("HectorEventProcessing") << "Hector:addPartToHepMC: LHCTransportLink " << theLink;
0557 theCorrespondenceMap.push_back(theLink);
0558
0559 if (m_verbosity)
0560 edm::LogVerbatim("HectorEventProcessing")
0561 << "Hector::TRANSPORTED pz= " << gpart->momentum().pz() << " eta= " << gpart->momentum().eta()
0562 << " status= " << gpart->status();
0563
0564 }
0565 }
0566 }
0567
0568 else {
0569 gpart = evt->barcode_to_particle(line);
0570 if (gpart) {
0571
0572 gpart->set_status(2);
0573
0574
0575
0576 if (m_verbosity)
0577 edm::LogVerbatim("HectorEventProcessing")
0578 << "Hector::NON-transp. pz= " << gpart->momentum().pz() << " eta= " << gpart->momentum().eta()
0579 << " status= " << gpart->status();
0580 }
0581 }
0582
0583 }
0584
0585 return evt;
0586 }