Back to home page

Project CMSSW displayed by LXR
 
 

    

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 &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   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   // construct beam line for FP420: .
0096   if (m_FP420Transport && lengthfp420 > 0.) {
0097     m_beamlineFP4201 = new H_BeamLine(1, lengthfp420 + 0.1);   // (direction, length)
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     // construct beam line for ZDC: .
0112     m_beamlineZDC1 = new H_BeamLine(1, lengthzdc + 0.1);   // (direction, length)
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     // construct beam line for D1: .
0122     m_beamlineD11 = new H_BeamLine(1, lengthd1 + 0.1);   // (direction, length)
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;  // neutrals
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;  // charged
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           // from mm to um
0200           double XforPosition = (*eventParticle)->production_vertex()->position().x() / micrometer;  // um
0201           double YforPosition = (*eventParticle)->production_vertex()->position().y() / micrometer;  // um
0202           double ZforPosition = (*eventParticle)->production_vertex()->position().z() / meter;       // m
0203           // crossing angle (beam tilt) is not known a priory; keep now 0.0 but,
0204           // in principle, can be entered as parameters
0205           double TXforPosition = 0.0, TYforPosition = 0.0;  // urad
0206 
0207           // Clears H_BeamParticle::positions and sets the initial one
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         }  // if find line
0225       }    // if eta > 8.2
0226     }      // if status
0227   }        // for loop
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           // the beam transverse direction is centered on (TXforPosition,
0253           // TYforPosition) at IP
0254           if (m_sigmaSTX > 0. && m_sigmaSTY > 0.) {
0255             part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0256           } else {
0257             // for smearAng() in urad, default are (STX=30.23, STY=30.23)
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);  // in GeV, default is SBE=0.79
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         // propagating
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       }  // if isCharged
0311       else {
0312         m_isStoppedfp420[line] = true;  // imply that neutral particles stopped to reach 420m
0313         if (m_verbosity)
0314           LogDebug("HectorEventProcessing")
0315               << "Hector:filterFP420: barcode = " << line << " isStopped=" << (*m_isStoppedfp420.find(line)).second;
0316       }
0317 
0318     }  // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
0319   }    // if ( m_beamPart.size() )
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             // the beam transverse direction is centered on (TXforPosition,
0350             // TYforPosition) at IP
0351             part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0352           } else {
0353             // for smearAng() in urad, default are (STX=30.23, STY=30.23)
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);  // in GeV, default is SBE=0.79
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       // if stopfp420 charged particles
0386       /*
0387         else if ( ((*m_isCharged.find(line)).second) ){
0388         m_isStoppedzdc[line] = false;// not stopped in propagating to FP420 and
0389         therefore in  propagation to ZDC too. if(m_verbosity)
0390         LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode = " <<
0391         line << " isStopped=" << (*m_isStoppedzdc.find(line)).second;
0392         } */
0393       else {
0394         m_isStoppedzdc[line] = true;  // neutrals particles considered as stopped
0395                                       // in propagating to ZDC
0396         if (m_verbosity)
0397           LogDebug("HectorEventProcessing")
0398               << "Hector:filterZDC: barcode = " << line << " isStopped=" << (*m_isStoppedzdc.find(line)).second;
0399       }
0400 
0401     }  // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
0402   }    // if ( m_beamPart.size() )
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             // the beam transverse direction is centered on (TXforPosition,
0433             // TYforPosition) at IP
0434             part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
0435           } else {
0436             // for smearAng() in urad, default are (STX=30.23, STY=30.23)
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);  // in GeV, default is SBE=0.79
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         // propagating
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       }  // if stopzdc
0480       else {
0481         m_isStoppedd1[line] = false;  // not stopped in propagating to ZDC and
0482                                       // therefore in  propagation to D1 too.
0483         if (m_verbosity)
0484           LogDebug("HectorEventProcessing")
0485               << "Hector:filterD1: barcode = " << line << " isStopped=" << (*m_isStoppedd1.find(line)).second;
0486       }
0487 
0488     }  // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
0489   }    // if ( m_beamPart.size() )
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);  // tx, ty never == 0?
0553         energy = (*m_eAtTrPoint.find(line)).second;
0554 
0555         time = (ddd * meter - gpart->production_vertex()->position().z() * mm);  // 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         }  // ddd
0593       }    // if gpart
0594     }      // if !isStopped
0595 
0596     else {
0597       gpart = evt->barcode_to_particle(line);
0598       if (gpart) {
0599         //        HepMC::GenVertex * vert= new HepMC::GenVertex();
0600         gpart->set_status(2);
0601         //        vert->add_particle_in( gpart );
0602         //        vert->add_particle_out( gpart );
0603         //        evt->add_vertex( vert );
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   }  // for
0611 
0612   return evt;
0613 }
[ 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!