Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // Create LHC beam line
0024   edm::ParameterSet hector_par = param.getParameter<edm::ParameterSet>("Hector");
0025 
0026   // User definitons
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   // construct beam line for FP420:                                                                                           .
0078   if (m_FP420Transport && lengthfp420 > 0.) {
0079     m_beamlineFP4201 = new H_BeamLine(1, lengthfp420 + 0.1);   // (direction, length)
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     // construct beam line for ZDC:
0094     m_beamlineZDC1 = new H_BeamLine(1, lengthzdc + 0.1);   // (direction, length)
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     // construct beam line for D1:
0104     m_beamlineD11 = new H_BeamLine(1, lengthd1 + 0.1);   // (direction, length)
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;  // neutrals
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;  //charged
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           // from mm to um
0182           double XforPosition = (*eventParticle)->production_vertex()->position().x() / CLHEP::micrometer;  //um
0183           double YforPosition = (*eventParticle)->production_vertex()->position().y() / CLHEP::micrometer;  //um
0184           double ZforPosition = (*eventParticle)->production_vertex()->position().z() / CLHEP::meter;       //m
0185           // crossing angle (beam tilt) is not known a priory; keep now 0.0 but, in principle, can be entered as parameters
0186           double TXforPosition = 0.0, TYforPosition = 0.0;  //urad
0187 
0188           // Clears H_BeamParticle::positions and sets the initial one
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         }  // if find line
0206       }  // if eta > 8.2
0207     }  // if status
0208   }  // for loop
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           // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
0234           if (m_sigmaSTX > 0. && m_sigmaSTY > 0.) {
0235             part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0236           } else {
0237             // for smearAng() in urad, default are (STX=30.23, STY=30.23)
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);  // in GeV, default is SBE=0.79
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         //propagating
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       }  // if isCharged
0291       else {
0292         m_isStoppedfp420[line] = true;  // imply that neutral particles stopped to reach 420m
0293         if (m_verbosity)
0294           edm::LogVerbatim("HectorEventProcessing")
0295               << "Hector:filterFP420: barcode = " << line << " isStopped=" << (*m_isStoppedfp420.find(line)).second;
0296       }
0297 
0298     }  // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
0299   }  // if ( m_beamPart.size() )
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       //      if ( ((*m_isStoppedfp420.find(line)).second) && ((*m_isCharged.find(line)).second) ) {
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             // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
0332             part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0333           } else {
0334             // for smearAng() in urad, default are (STX=30.23, STY=30.23)
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);  // in GeV, default is SBE=0.79
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       // if stopfp420 charged particles
0367       else {
0368         m_isStoppedzdc[line] = true;  // neutrals particles considered as stopped in propagating to ZDC
0369         if (m_verbosity)
0370           edm::LogVerbatim("HectorEventProcessing")
0371               << "Hector:filterZDC: barcode = " << line << " isStopped=" << (*m_isStoppedzdc.find(line)).second;
0372       }
0373 
0374     }  // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
0375   }  // if ( m_beamPart.size() )
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             // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
0406             part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0407           } else {
0408             // for smearAng() in urad, default are (STX=30.23, STY=30.23)
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);  // in GeV, default is SBE=0.79
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         //propagating
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       }  // if stopzdc
0452       else {
0453         m_isStoppedd1[line] = false;  // not stopped in propagating to ZDC and therefore in  propagation to D1 too.
0454         if (m_verbosity)
0455           edm::LogVerbatim("HectorEventProcessing")
0456               << "Hector:filterD1: barcode = " << line << " isStopped=" << (*m_isStoppedd1.find(line)).second;
0457       }
0458 
0459     }  // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
0460   }  // if ( m_beamPart.size() )
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);  // tx, ty never == 0?
0523         energy = (*m_eAtTrPoint.find(line)).second;
0524 
0525         time = (ddd * CLHEP::meter - gpart->production_vertex()->position().z() * CLHEP::mm);  // 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         }  // ddd
0565       }  // if gpart
0566     }  // if !isStopped
0567 
0568     else {
0569       gpart = evt->barcode_to_particle(line);
0570       if (gpart) {
0571         //        HepMC::GenVertex * vert= new HepMC::GenVertex();
0572         gpart->set_status(2);
0573         //        vert->add_particle_in( gpart );
0574         //        vert->add_particle_out( gpart );
0575         //        evt->add_vertex( vert );
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   }  //for
0584 
0585   return evt;
0586 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!