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