Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:21:57 

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