Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:19:20 

0001 #include "HepMCNodeReader.h"
0002 #include "PHG4InEvent.h"
0003 #include "PHG4Particle.h"
0004 #include "PHG4Particlev1.h"
0005 
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 
0008 #include <phhepmc/PHHepMCDefs.h>
0009 #include <phhepmc/PHHepMCGenEvent.h>
0010 #include <phhepmc/PHHepMCGenEventMap.h>
0011 
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHDataNode.h>
0014 #include <phool/PHNode.h>  // for PHNode
0015 #include <phool/PHNodeIterator.h>
0016 #include <phool/PHObject.h>
0017 #include <phool/PHRandomSeed.h>
0018 #include <phool/getClass.h>
0019 #include <phool/phool.h>
0020 #include <phool/recoConsts.h>
0021 
0022 #pragma GCC diagnostic push
0023 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0024 #include <HepMC/GenEvent.h>
0025 #pragma GCC diagnostic pop
0026 #include <HepMC/GenParticle.h>
0027 #include <HepMC/GenVertex.h>
0028 #include <HepMC/IteratorRange.h>
0029 #include <HepMC/SimpleVector.h>
0030 #include <HepMC/Units.h>
0031 
0032 #include <CLHEP/Vector/LorentzRotation.h>
0033 #include <CLHEP/Vector/LorentzVector.h>
0034 
0035 #include <gsl/gsl_const.h>
0036 #include <gsl/gsl_randist.h>
0037 #include <gsl/gsl_rng.h>
0038 
0039 #include <TDatabasePDG.h>
0040 #include <TLorentzVector.h>
0041 
0042 #include <cassert>
0043 #include <cmath>
0044 #include <cstdlib>
0045 #include <iostream>
0046 #include <iterator>
0047 #include <list>
0048 #include <utility>
0049 
0050 //// All length Units are in cm, no conversion to G4 internal units since
0051 //// this is filled into our objects (PHG4VtxPoint and PHG4Particle)
0052 //
0053 //// pythia vtx time seems to be in mm/c
0054 // const double mm_over_c_to_sec = 0.1 / GSL_CONST_CGS_SPEED_OF_LIGHT;
0055 //// pythia vtx time seems to be in mm/c
0056 // const double mm_over_c_to_nanosecond = mm_over_c_to_sec * 1e9;
0057 /// \class  IsStateFinal
0058 
0059 /// this predicate returns true if the input has no decay vertex
0060 class IsStateFinal
0061 {
0062  public:
0063   /// returns true if the GenParticle does not decay
0064   bool operator()(const HepMC::GenParticle *p)
0065   {
0066     if (!p->end_vertex() && p->status() == 1)
0067     {
0068       return true;
0069     }
0070     return false;
0071   }
0072 };
0073 
0074 static IsStateFinal isfinal;
0075 
0076 HepMCNodeReader::HepMCNodeReader(const std::string &name)
0077   : SubsysReco(name)
0078 {
0079   RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
0080   return;
0081 }
0082 
0083 HepMCNodeReader::~HepMCNodeReader()
0084 {
0085   gsl_rng_free(RandomGenerator);
0086 }
0087 
0088 int HepMCNodeReader::Init(PHCompositeNode *topNode)
0089 {
0090   PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0091   if (!ineve)
0092   {
0093     PHNodeIterator iter(topNode);
0094     PHCompositeNode *dstNode;
0095     dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0096 
0097     ineve = new PHG4InEvent();
0098     PHDataNode<PHObject> *newNode =
0099         new PHDataNode<PHObject>(ineve, "PHG4INEVENT", "PHObject");
0100     dstNode->addNode(newNode);
0101   }
0102   unsigned int phseed = PHRandomSeed();  // fixed seed is handled in this funtcion, need to call it to preserve random numbder order even if we override it
0103   if (use_seed)
0104   {
0105     phseed = seed;
0106     std::cout << Name() << " override random seed: " << phseed << std::endl;
0107   }
0108   gsl_rng_set(RandomGenerator, phseed);
0109 
0110   if (addfraction < 0.0)
0111   {
0112     std::cout << "[WARNING] addfraction is negative which is not allowed. Setting addfraction to 0." << std::endl;
0113     addfraction = 0.0;
0114   }
0115 
0116   if (addfraction != 0.0)
0117   {
0118     fpt = new TF1("fpt", EMGFunction, -10, 10, 4);
0119     fpt->SetParameter(0, 1);            // normalization
0120     fpt->SetParameter(1, 4.71648e-01);  // mu
0121     fpt->SetParameter(2, 1.89602e-01);  // sigma
0122     fpt->SetParameter(3, 2.26981e+00);  // lambda
0123 
0124     std::cout << "[INFO] eta is sampled between [" << -1 * sel_eta << "," << sel_eta << "]" << std::endl;
0125     feta = new TF1("feta", DBGFunction, -1*sel_eta, sel_eta, 4);
0126     feta->SetParameter(0, 1);             // normalization
0127     feta->SetParameter(1, -4.08301e-01);  // mu1
0128     feta->SetParameter(2, 4.11930e-01);   // mu2
0129     feta->SetParameter(3, 3.59063e-01);   // sigma
0130 
0131     std::vector<std::pair<int, double>> l_PIDProb;
0132     for (size_t i = 0; i < list_strangePID.size(); i++)
0133     {
0134       l_PIDProb.emplace_back(list_strangePID[i], list_strangePIDprob[i]);
0135     }
0136 
0137     std::sort(l_PIDProb.begin(), l_PIDProb.end(),
0138               [](const std::pair<int, double> &a, const std::pair<int, double> &b) -> bool
0139               {
0140                 return a.second < b.second;
0141               });
0142 
0143     double sum = 0.0;
0144     for (const auto &it : l_PIDProb)
0145     {
0146       sum += it.second;
0147       list_strangePID_probrange.emplace_back(it.first, std::make_pair(sum - it.second, sum));
0148     }
0149 
0150     if (Verbosity() > 0)
0151     {
0152       std::cout << "[INFO] Sorted list of strange particles and their probabilities: " << std::endl;
0153 
0154       for (const auto &it : l_PIDProb)
0155       {
0156         std::cout << "PID: " << it.first << " Probability: " << it.second << std::endl;
0157       }
0158 
0159       std::cout << "[INFO] List of strange particles and their probability ranges: " << std::endl;
0160       for (const auto &it : list_strangePID_probrange)
0161       {
0162         std::cout << "PID: " << it.first << " Probability range: [" << it.second.first << "," << it.second.second << "]" << std::endl;
0163       }
0164     }
0165   }
0166 
0167   return 0;
0168 }
0169 
0170 int HepMCNodeReader::process_event(PHCompositeNode *topNode)
0171 {
0172   // For pile-up simulation: define GenEventMap
0173   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0174 
0175   if (!genevtmap)
0176   {
0177     static bool once = true;
0178 
0179     if (once and Verbosity())
0180     {
0181       once = false;
0182 
0183       std::cout << "HepMCNodeReader::process_event - No PHHepMCGenEventMap node. Do not perform HepMC->Geant4 input" << std::endl;
0184     }
0185 
0186     return Fun4AllReturnCodes::DISCARDEVENT;
0187   }
0188 
0189   PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0190   if (!ineve)
0191   {
0192     std::cout << PHWHERE << "no PHG4INEVENT node" << std::endl;
0193     return Fun4AllReturnCodes::ABORTEVENT;
0194   }
0195 
0196   recoConsts *rc = recoConsts::instance();
0197 
0198   float worldsizex = rc->get_FloatFlag("WorldSizex");
0199   float worldsizey = rc->get_FloatFlag("WorldSizey");
0200   float worldsizez = rc->get_FloatFlag("WorldSizez");
0201   std::string worldshape = rc->get_StringFlag("WorldShape");
0202 
0203   enum
0204   {
0205     ShapeG4Tubs = 0,
0206     ShapeG4Box = 1
0207   };
0208 
0209   int ishape;
0210   if (worldshape == "G4Tubs")
0211   {
0212     ishape = ShapeG4Tubs;
0213   }
0214   else if (worldshape == "G4Box")
0215   {
0216     ishape = ShapeG4Box;
0217   }
0218   else
0219   {
0220     std::cout << PHWHERE << " unknown world shape " << worldshape << std::endl;
0221     exit(1);
0222   }
0223 
0224   // For pile-up simulation: loop over PHHepMC event map
0225   // insert highest embedding ID event first, whose vertex maybe resued in  PHG4ParticleGeneratorBase::ReuseExistingVertex()
0226   int vtxindex = -1;
0227   bool use_embedding_vertex = false;
0228   HepMC::FourVector collisionVertex;
0229   PHHepMCGenEvent *evtvertex = genevtmap->get(PHHepMCDefs::DataVertexIndex);
0230   if (evtvertex)
0231   {
0232     collisionVertex = evtvertex->get_collision_vertex();
0233     genevtmap->erase(PHHepMCDefs::DataVertexIndex);
0234     use_embedding_vertex = true;
0235   }
0236   for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
0237   {
0238     PHHepMCGenEvent *genevt = iter->second;
0239     assert(genevt);
0240 
0241     if (genevt->is_simulated())
0242     {
0243       if (Verbosity())
0244       {
0245         std::cout << "HepMCNodeReader::process_event - this event is already simulated. Move on: ";
0246         genevt->identify();
0247       }
0248 
0249       continue;
0250     }
0251 
0252     if (Verbosity())
0253     {
0254       std::cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found PHHepMCGenEvent:" << std::endl;
0255       genevt->identify();
0256     }
0257 
0258     if (!use_embedding_vertex)
0259     {
0260       collisionVertex = genevt->get_collision_vertex();
0261     }
0262     else
0263     {
0264       genevt->set_collision_vertex(collisionVertex);  // save used vertex in HepMC
0265     }
0266     const int embed_flag = genevt->get_embedding_id();
0267     HepMC::GenEvent *evt = genevt->getEvent();
0268     if (!evt)
0269     {
0270       std::cout << PHWHERE << " no evt pointer under HEPMC Node found";
0271       genevt->identify();
0272       return Fun4AllReturnCodes::ABORTEVENT;
0273     }
0274 
0275     if (Verbosity())
0276     {
0277       std::cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found HepMC::GenEvent:" << std::endl;
0278       evt->print();
0279     }
0280 
0281     genevt->is_simulated(true);
0282     double xshift = vertex_pos_x + genevt->get_collision_vertex().x();
0283     double yshift = vertex_pos_y + genevt->get_collision_vertex().y();
0284     double zshift = vertex_pos_z + genevt->get_collision_vertex().z();
0285     double tshift = vertex_t0 + genevt->get_collision_vertex().t();
0286     const CLHEP::HepLorentzRotation lortentz_rotation(genevt->get_LorentzRotation_EvtGen2Lab());
0287 
0288     if (width_vx > 0.0)
0289     {
0290       xshift += smeargauss(width_vx);
0291     }
0292     else if (width_vx < 0.0)
0293     {
0294       xshift += smearflat(width_vx);
0295     }
0296 
0297     if (width_vy > 0.0)
0298     {
0299       yshift += smeargauss(width_vy);
0300     }
0301     else if (width_vy < 0.0)
0302     {
0303       yshift += smearflat(width_vy);
0304     }
0305 
0306     if (width_vz > 0.0)
0307     {
0308       zshift += smeargauss(width_vz);
0309     }
0310     else if (width_vz < 0.0)
0311     {
0312       zshift += smearflat(width_vz);
0313     }
0314 
0315     std::list<HepMC::GenParticle *> finalstateparticles;
0316     std::list<HepMC::GenParticle *>::const_iterator fiter;
0317 
0318     int Nstrange = 0;  // count the number of strange particles with PID in list_strangePID per event
0319 
0320     // units in G4 interface are GeV and CM as in PHENIX convention
0321     const double mom_factor = HepMC::Units::conversion_factor(evt->momentum_unit(), HepMC::Units::GEV);
0322     const double length_factor = HepMC::Units::conversion_factor(evt->length_unit(), HepMC::Units::CM);
0323     const double time_factor = HepMC::Units::conversion_factor(evt->length_unit(), HepMC::Units::CM) / GSL_CONST_CGS_SPEED_OF_LIGHT * 1e9;  // from length_unit()/c to ns
0324 
0325     for (HepMC::GenEvent::vertex_iterator v = evt->vertices_begin();
0326          v != evt->vertices_end();
0327          ++v)
0328     {
0329       if (Verbosity() > 1)
0330       {
0331         std::cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found vertex:" << std::endl;
0332         (*v)->print();
0333       }
0334 
0335       finalstateparticles.clear();
0336       for (HepMC::GenVertex::particle_iterator p =
0337                (*v)->particles_begin(HepMC::children);
0338            p != (*v)->particles_end(HepMC::children); ++p)
0339       {
0340         if (Verbosity() > 1)
0341         {
0342           std::cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found particle:" << std::endl;
0343           (*p)->print();
0344           std::cout << "end vertex " << (*p)->end_vertex() << std::endl;
0345         }
0346         if (isfinal(*p))
0347         {
0348           if (Verbosity() > 1)
0349           {
0350             std::cout << __PRETTY_FUNCTION__ << " " << __LINE__ << std::endl;
0351             std::cout << "\tparticle passed " << std::endl;
0352           }
0353           finalstateparticles.push_back(*p);
0354 
0355           if (std::find(list_strangePID.begin(), list_strangePID.end(), std::abs((*p)->pdg_id())) != list_strangePID.end())
0356           {
0357             // count number of strange particles with PID in list_strangePID within the kinematic selection (sPHNEIX acceptance)
0358             // |eta|<=1 and momentum within sel_ptmin and sel_ptmax
0359             if ((fabs((*p)->momentum().eta()) <= sel_eta) && //
0360                 ((*p)->momentum().perp() >= sel_ptmin && (*p)->momentum().perp() <= sel_ptmax) //
0361             )
0362             {
0363               Nstrange++;
0364             }
0365           }
0366         }
0367         else
0368         {
0369           if (Verbosity() > 1)
0370           {
0371             std::cout << __PRETTY_FUNCTION__ << " " << __LINE__ << std::endl;
0372             std::cout << "\tparticle failed " << std::endl;
0373           }
0374         }
0375       }  // for (HepMC::GenVertex::particle_iterator p = (*v)->particles_begin(HepMC::children); p != (*v)->particles_end(HepMC::children); ++p)
0376 
0377       // Add additional strange particles if addfraction is not zero; round up to the ceiling integer
0378       Nstrange_add = static_cast<int>(std::ceil(Nstrange * addfraction * 0.01));
0379 
0380       if (Verbosity() > 1)
0381       {
0382         std::cout << "[DEBUG] Number of original strange particles with kinematic selection (abs(eta) cut: " << sel_eta << ", pT cut: " << sel_ptmin << " - " << sel_ptmax << "): " << Nstrange << std::endl;
0383         std::cout << "[DEBUG] addfraction: " << addfraction << "%; Number of strange particles to be added: " << Nstrange_add << std::endl;
0384       }
0385 
0386       if (!finalstateparticles.empty())
0387       {
0388         CLHEP::HepLorentzVector lv_vertex((*v)->position().x(),
0389                                           (*v)->position().y(),
0390                                           (*v)->position().z(),
0391                                           (*v)->position().t());
0392         if (is_pythia)
0393         {
0394           lv_vertex.setX(collisionVertex.x());
0395           lv_vertex.setY(collisionVertex.y());
0396           lv_vertex.setZ(collisionVertex.z());
0397           lv_vertex.setT(collisionVertex.t());
0398           if (Verbosity() > 1)
0399           {
0400             std::cout << __PRETTY_FUNCTION__ << " " << __LINE__
0401                       << std::endl;
0402             std::cout << "\t vertex reset to collision vertex: "
0403                       << lv_vertex << std::endl;
0404           }
0405         }
0406 
0407         // event gen frame to lab frame
0408         lv_vertex = lortentz_rotation(lv_vertex);
0409 
0410         double xpos = lv_vertex.x() * length_factor + xshift;
0411         double ypos = lv_vertex.y() * length_factor + yshift;
0412         double zpos = lv_vertex.z() * length_factor + zshift;
0413         double time = lv_vertex.t() * time_factor + tshift;
0414 
0415         if (Verbosity() > 1)
0416         {
0417           std::cout << __PRETTY_FUNCTION__ << " " << __LINE__ << std::endl;
0418           std::cout << "Vertex : " << std::endl;
0419           (*v)->print();
0420           std::cout << "id: " << (*v)->barcode() << std::endl;
0421           std::cout << "x: " << xpos << std::endl;
0422           std::cout << "y: " << ypos << std::endl;
0423           std::cout << "z: " << zpos << std::endl;
0424           std::cout << "t: " << time << std::endl;
0425           std::cout << "Particles" << std::endl;
0426         }
0427 
0428         if (ishape == ShapeG4Tubs)
0429         {
0430           if (sqrt(xpos * xpos + ypos * ypos) > worldsizey / 2 ||
0431               fabs(zpos) > worldsizez / 2)
0432           {
0433             std::cout << "vertex x/y/z " << xpos << "/" << ypos << "/" << zpos
0434                       << " id: " << (*v)->barcode()
0435                       << " outside world volume radius/z (+-) " << worldsizex / 2
0436                       << "/" << worldsizez / 2 << ", dropping it and its particles"
0437                       << std::endl;
0438             continue;
0439           }
0440         }
0441         else if (ishape == ShapeG4Box)
0442         {
0443           if (fabs(xpos) > worldsizex / 2 || fabs(ypos) > worldsizey / 2 ||
0444               fabs(zpos) > worldsizez / 2)
0445           {
0446             std::cout << "Vertex x/y/z " << xpos << "/" << ypos << "/" << zpos
0447                       << " outside world volume x/y/z (+-) " << worldsizex / 2 << "/"
0448                       << worldsizey / 2 << "/" << worldsizez / 2
0449                       << ", dropping it and its particles" << std::endl;
0450             continue;
0451           }
0452         }
0453         else
0454         {
0455           std::cout << PHWHERE << " shape " << ishape << " not implemented. exiting"
0456                     << std::endl;
0457           exit(1);
0458         }
0459 
0460         // For pile-up simulation: vertex position
0461         vtxindex = ineve->AddVtx(xpos, ypos, zpos, time);
0462         for (fiter = finalstateparticles.begin();
0463              fiter != finalstateparticles.end();
0464              ++fiter)
0465         {
0466           if (Verbosity() > 1)
0467           {
0468             std::cout << __PRETTY_FUNCTION__ << " " << __LINE__ << std::endl;
0469             (*fiter)->print();
0470           }
0471 
0472           CLHEP::HepLorentzVector lv_momentum((*fiter)->momentum().px(),
0473                                               (*fiter)->momentum().py(),
0474                                               (*fiter)->momentum().pz(),
0475                                               (*fiter)->momentum().e());
0476 
0477           // event gen frame to lab frame
0478           lv_momentum = lortentz_rotation(lv_momentum);
0479 
0480           PHG4Particle *particle = new PHG4Particlev1();
0481           particle->set_pid((*fiter)->pdg_id());
0482           particle->set_px(lv_momentum.x() * mom_factor);
0483           particle->set_py(lv_momentum.y() * mom_factor);
0484           particle->set_pz(lv_momentum.z() * mom_factor);
0485           particle->set_barcode((*fiter)->barcode());
0486 
0487           ineve->AddParticle(vtxindex, particle);
0488 
0489           if (embed_flag != 0)
0490           {
0491             ineve->AddEmbeddedParticle(particle, embed_flag);
0492           }
0493         }  // for (fiter = finalstateparticles.begin(); fiter != finalstateparticles.end(); ++fiter)
0494 
0495         // add strange particles given Nstrange_add
0496         if (addfraction > 0)
0497         {
0498           if (Verbosity() > 1)
0499           {
0500             std::cout << "[INFO] Add strange particles. Number of strange particles to be added: " << Nstrange_add << std::endl;
0501           }
0502 
0503           // Add strange particles given Nstrange_add
0504           for (int i = 0; i < Nstrange_add; i++)
0505           {
0506             int pid = list_strangePID[0];  // default to the first PID in the list, i.e K_s0
0507             double prob = gsl_rng_uniform_pos(RandomGenerator);
0508             for (const auto &it : list_strangePID_probrange)
0509             {
0510               if (prob >= it.second.first && prob < it.second.second)
0511               {
0512                 pid = it.first;
0513                 break;
0514               }
0515             }
0516 
0517             // sample pt and eta from the EMG and DBG functions; phi between -pi and pi
0518             double pt = fpt->GetRandom();
0519             double eta = feta->GetRandom();
0520             double phi = (gsl_rng_uniform_pos(RandomGenerator) * 2 * M_PI) - M_PI;
0521             double mass = TDatabasePDG::Instance()->GetParticle(pid)->Mass();
0522 
0523             TLorentzVector lv;
0524             lv.SetPtEtaPhiM(pt, eta, phi, mass);
0525 
0526             // create a new particle
0527             PHG4Particle *particle = new PHG4Particlev1();
0528             particle->set_pid(pid);
0529             particle->set_px(lv.Px());
0530             particle->set_py(lv.Py());
0531             particle->set_pz(lv.Pz());
0532             particle->set_barcode(std::numeric_limits<int>::max() - i);  // set the barcode to be distinct from the existing particles; backward counting from the maximum integer value
0533 
0534             ineve->AddParticle(vtxindex, particle);
0535           }
0536         }
0537 
0538       }  // if (!finalstateparticles.empty())
0539     }    // for (HepMC::GenEvent::vertex_iterator v = evt->vertices_begin();
0540   }      // For pile-up simulation: loop end for PHHepMC event map
0541   if (Verbosity() > 0)
0542   {
0543     ineve->identify();
0544   }
0545 
0546   return Fun4AllReturnCodes::EVENT_OK;
0547 }
0548 
0549 double HepMCNodeReader::smeargauss(const double width)
0550 {
0551   if (width == 0)
0552   {
0553     return 0;
0554   }
0555   return gsl_ran_gaussian(RandomGenerator, width);
0556 }
0557 
0558 double HepMCNodeReader::smearflat(const double width)
0559 {
0560   if (width == 0)
0561   {
0562     return 0;
0563   }
0564   return 2.0 * width * (gsl_rng_uniform_pos(RandomGenerator) - 0.5);
0565 }
0566 
0567 void HepMCNodeReader::VertexPosition(const double v_x, const double v_y,
0568                                      const double v_z)
0569 {
0570   std::cout << "HepMCNodeReader::VertexPosition - WARNING - this function is depreciated. "
0571             << "HepMCNodeReader::VertexPosition() move all HEPMC subevents to a new vertex location. "
0572             << "This also leads to a different vertex is used for HepMC subevent in Geant4 than that recorded in the HepMCEvent Node."
0573             << "Recommendation: the vertex shifts are better controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
0574             << std::endl;
0575 
0576   vertex_pos_x = v_x;
0577   vertex_pos_y = v_y;
0578   vertex_pos_z = v_z;
0579   return;
0580 }
0581 
0582 void HepMCNodeReader::SmearVertex(const double s_x, const double s_y,
0583                                   const double s_z)
0584 {
0585   std::cout << "HepMCNodeReader::SmearVertex - WARNING - this function is depreciated. "
0586             << "HepMCNodeReader::SmearVertex() smear each HEPMC subevents to a new vertex location. "
0587             << "This also leads to a different vertex is used for HepMC subevent in Geant4 than that recorded in the HepMCEvent Node."
0588             << "Recommendation: the vertex smears are better controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
0589             << std::endl;
0590 
0591   width_vx = s_x;
0592   width_vy = s_y;
0593   width_vz = s_z;
0594   return;
0595 }
0596 
0597 double HepMCNodeReader::EMGFunction(double *x, double *par)
0598 {
0599   // parameterization: https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution
0600 
0601   double N = par[0];       // Normalization
0602   double mu = par[1];      // Mean of the Gaussian
0603   double sigma = par[2];   // Width of the Gaussian
0604   double lambda = par[3];  // Decay constant of the Exponential
0605 
0606   double t = x[0];
0607   double z = (mu + lambda * sigma * sigma - t) / (sqrt(2) * sigma);
0608 
0609   double prefactor = lambda / 2.0;
0610   double exp_part = exp((lambda / 2.0) * (2.0 * mu + lambda * sigma * sigma - 2.0 * t));
0611   double erfc_part = TMath::Erfc(z);
0612 
0613   return N * prefactor * exp_part * erfc_part;
0614 }
0615 
0616 double HepMCNodeReader::DBGFunction(double *x, double *par)
0617 {
0618   double N = par[0];      // Normalization
0619   double mu1 = par[1];    // Mean of the first Gaussian
0620   double mu2 = par[2];    // Mean of the second Gaussian
0621   double sigma = par[3];  // Width of the Gaussian
0622 
0623   return N * (TMath::Gaus(x[0], mu1, sigma) + TMath::Gaus(x[0], mu2, sigma));
0624 }
0625 
0626 void HepMCNodeReader::Embed(const int /*unused*/)
0627 {
0628   std::cout << "HepMCNodeReader::Embed - WARNING - this function is depreciated. "
0629             << "Embedding IDs are controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
0630             << std::endl;
0631 
0632   return;
0633 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!