From d2319ddc5b0137a41126f738d0a842ec400de648 Mon Sep 17 00:00:00 2001 From: Davide Lancierini Date: Wed, 24 Jun 2026 10:57:29 +0100 Subject: [PATCH] add possibility of B decay time saving to enable mixing effects in RapidSim + various other --- src/RapidConfig.cc | 3 + src/RapidConfig.h | 2 + src/RapidExternalEvtGen.cc | 211 +++++++++++++++++++++++++++++-------- src/RapidExternalEvtGen.h | 4 +- src/RapidParam.cc | 54 +++++++++- src/RapidParam.h | 7 +- src/RapidParticle.h | 14 ++- src/RapidParticleData.cc | 6 +- 8 files changed, 252 insertions(+), 49 deletions(-) diff --git a/src/RapidConfig.cc b/src/RapidConfig.cc index 4af94ce..482937a 100644 --- a/src/RapidConfig.cc +++ b/src/RapidConfig.cc @@ -439,6 +439,9 @@ bool RapidConfig::configParticle(unsigned int part, TString command, TString val parts_[part]->setEvtGenDecayModel(value); std::cout << "INFO in RapidConfig::configParticle : set EvtGen decay model for particle " << parts_[part]->name() << std::endl << " : " << value << std::endl; + } else if(command=="Define"){ + std::cout<<"Define: "<setEvtGenDefinitions(value); } return true; diff --git a/src/RapidConfig.h b/src/RapidConfig.h index 404c22c..4c6d466 100644 --- a/src/RapidConfig.h +++ b/src/RapidConfig.h @@ -99,6 +99,8 @@ class RapidConfig { //Vtx smearing lookup for each smearing category, placeholder for now //std::map vtxSmearCategories_; + // std::vector> definitions_; + std::vector definitions_; //accept reject hist to sculpt kinematics TH1* accRejHisto_; RapidParam* accRejParameterX_; diff --git a/src/RapidExternalEvtGen.cc b/src/RapidExternalEvtGen.cc index 63405c1..e3deb16 100644 --- a/src/RapidExternalEvtGen.cc +++ b/src/RapidExternalEvtGen.cc @@ -6,10 +6,14 @@ #include "TRandom.h" #include "TSystem.h" +#include "TObjString.h" #ifdef RAPID_EVTGEN #include "EvtGen/EvtGen.hh" +#include "EvtGenBase/EvtConst.hh" #include "EvtGenBase/EvtRandomEngine.hh" +#include "EvtGenBase/EvtCPUtil.hh" +#include "EvtGenBase/EvtHepMCEvent.hh" #include "EvtGenBase/EvtAbsRadCorr.hh" #include "EvtGenBase/EvtDecayBase.hh" #include "EvtGenBase/EvtParticle.hh" @@ -21,78 +25,156 @@ #include "EvtGenExternal/EvtExternalGenList.hh" #endif +int B0Id( 511 ), B0barId( -511 ); +int BsId( 531 ), BsbarId( -531 ); + + bool RapidExternalEvtGen::decay(std::vector& parts) { #ifdef RAPID_EVTGEN - EvtParticle* theParent(0); if(parts.size() < 1) { std::cout << "WARNING in RapidExternalEvtGen::decay : There are no particles to decay." << std::endl; return false; } + + EvtSpinDensity* spinDensity = 0; + TString parentName = getEvtGenName(parts[0]->id()); EvtId theId = EvtPDL::evtIdFromLundKC(parts[0]->id()); + int PDGId = EvtPDL::getStdHep( theId ); + // std::cout<<"PDGId here: "<setDiag( EvtSpinType::getSpinStates( EvtSpinType::VECTOR ) ); + spinDensity->set( 1, 1, EvtComplex( 0.0, 0.0 ) ); + } + + // RapidVertex* vtx(0); + // vtx = parts[0]->getOriginVertex(); + + // Parent particle XYZ-position + // ROOT::Math::XYZPoint point = vtx->getVertex(true); + // EvtVector4R origin(0.0, point.X(), point.Y(), point.Z()); + EvtVector4R origin(0.0, 0.0, 0.0, 0.0); + // Parent particle 4-momentum TLorentzVector pIn = parts[0]->getP(); EvtVector4R pInit(pIn.E(),pIn.Px(),pIn.Py(),pIn.Pz()); + // double mass = EvtPDL::getMeanMass( theId ); + // EvtVector4R pInit(mass,0.0, 0.0, 0.0); - theParent = EvtParticleFactory::particleFactory(theId, pInit); - if (theParent->getSpinStates() == 3) {theParent->setVectorSpinDensity();} + EvtHepMCEvent* theEvent = + evtGen_->generateDecay( PDGId, pInit, origin, spinDensity ); + // Retrieve the HepMC event information + GenEvent* hepMCEvent = theEvent->getEvent(); - // Generate the event - evtGen_->generateDecay(theParent); + std::list allVertices; + + double flightTime; - // Store particles to read in the order RapidSim stores them - std::queue evtParts; - // Also store the number of children expected for each of these particles so we can remove PHOTOS photons - std::queue nExpectedChildren; - evtParts.push(theParent); - nExpectedChildren.push(parts[0]->nDaughters()); + int iVtx(0); + int iPart(1); + bool hasOsc(false); - EvtVector4R x4Evt; - EvtVector4R p4Evt; TLorentzVector p4TLV; + FourVector FourMom; + FourVector DecayVtx; + FourVector OrigVtx; - int iPart=1; // The momentum and origin vertex of the first particle are already set + std::queue nExpectedChildren; + nExpectedChildren.push(parts[0]->nDaughters()); - while(!evtParts.empty()) { - EvtParticle* theParticle = evtParts.front(); + for ( auto theVertex : hepMCEvent->vertices() ) { + if ( theVertex == 0 ) { + continue; + } + auto nin = theVertex->particles_in_size(); + auto nout = theVertex->particles_out_size(); - // B0 and Bs may mix in EvtGen - RapidSim ignores this step and only records the second state - while(theParticle->getNDaug()==1) { - theParticle = theParticle->getDaug(0); + while (nin == 1 && nout ==1){ + //B meson oscillation + continue; } uint nChildren = nExpectedChildren.front(); - // Loop over the daughter tracks - for (uint iChild = 0; iChild < nChildren; ++iChild) { - EvtParticle* child = theParticle->getDaug(iChild); - - if (child != 0) { - p4Evt = child->getP4Lab(); - x4Evt = child->get4Pos(); - p4TLV.SetPxPyPzE(p4Evt.get(1),p4Evt.get(2),p4Evt.get(3),p4Evt.get(0)); - if(parts.size() < iPart+1u) { - std::cout << "WARNING in RapidExternalEvtGen::decay : EvtGen has produced too many particles." << std::endl; - return false; - } + // For these, get the mother decay vertex position and the 4-momentum to calculate + // the flight time. + + for ( auto inParticle : theVertex->particles_in() ) { + + if ( inParticle == 0 ) { + continue; + } + + int inPDGId = inParticle->pdg_id(); + FourMom = inParticle->momentum(); + DecayVtx = theVertex->position(); + + if ( inPDGId == B0Id || inPDGId == BsId || inPDGId == B0barId || inPDGId == BsbarId ) { + + if (PDGId!=inPDGId){ + hasOsc=true; + } + else{ + hasOsc=false; + + } + + parts[iVtx]->setHasOsc(hasOsc); + + } + + flightTime = calcFlightTime( DecayVtx, FourMom ); + parts[iVtx]->setId(inPDGId); + parts[iVtx]->setDecaytime(flightTime); + parts[iVtx]->getDecayVertex()->setXYZ(DecayVtx.x(), DecayVtx.y(), DecayVtx.x()); + + } + + uint nRecorded(0); + for ( auto outParticle : theVertex->particles_out() ) { + + if (nRecorded < nChildren) { + + FourMom = outParticle->momentum(); + OrigVtx = theVertex->position(); + + p4TLV.SetPxPyPzE(FourMom.px(),FourMom.py(),FourMom.pz(),FourMom.e()); + + // std::cout<<"ID:"<pdg_id()<setP(p4TLV); - parts[iPart]->getOriginVertex()->setXYZ(x4Evt.get(1),x4Evt.get(2),x4Evt.get(3)); - evtParts.push(child); + + parts[iPart]->getOriginVertex()->setXYZ(OrigVtx.x(),OrigVtx.y(),OrigVtx.z()); + parts[iPart]->setId(outParticle->pdg_id()); + nExpectedChildren.push(parts[iPart]->nDaughters()); + ++iPart; - } - } - // Clean up any PHOTOS photons - for (uint iChild = nChildren; iChild < theParticle->getNDaug(); ++iChild) { - delete theParticle->getDaug(iChild); + } + ++nRecorded; + } - delete theParticle; - evtParts.pop(); nExpectedChildren.pop(); + ++iVtx; } + delete hepMCEvent; + return true; #else if(!suppressWarning_) { @@ -152,8 +234,11 @@ bool RapidExternalEvtGen::setupGenerator() { decPath += "/config/evtgen/DECAY.DEC"; } + int mixingType = EvtCPUtil::Incoherent; + // int mixingType = EvtCPUtil::Coherent; + evtGen_ = new EvtGen(decPath.Data(), evtPDLPath.Data(), randomEngine, - radCorrEngine, &extraModels); + radCorrEngine, &extraModels, mixingType); return true; #else @@ -173,11 +258,28 @@ void RapidExternalEvtGen::writeDecFile(TString fname, std::vectorevtGenDefinitions()){ + defString = parts[iPart]->evtGenDefinitions(); + TObjArray* tokens = defString.Tokenize(","); + for (int i = 0; i < tokens->GetEntries(); i++) { + TString token = tokens->At(i)->GetName(); + token = token.Strip(TString::kBoth); // removes leading and trailing whitespace + std::cout << "INFO: Defining in DecFile: " << token << std::endl; + fout << "Define " << token << std::endl; + } + delete tokens; + } + } + fout << "\n" ; + // Loop over all particles and write out Decay rule for each for(unsigned int iPart=0; iPartnDaughters(); @@ -226,3 +328,26 @@ TString RapidExternalEvtGen::getEvtGenConjName(int id) { return ""; #endif } + +double RapidExternalEvtGen::calcFlightTime( FourVector& DecayVtx, FourVector& P4mtm ) +{ + double flightTime( 0.0 ); + +#ifdef EVTGEN_HEPMC3 + double distance = DecayVtx.length(); // mm + double momentum = P4mtm.length(); // GeV/c + double mass = P4mtm.m(); +#else + double distance = DecayVtx.rho(); // mm + double momentum = P4mtm.rho(); + double mass = P4mtm.m(); +#endif + + double c0 = EvtConst::c ;//mm/s + if ( momentum > 0.0 ) { + flightTime = 1.0e12 * distance * mass / + ( momentum * c0 ); // picoseconds + } + + return flightTime; +} \ No newline at end of file diff --git a/src/RapidExternalEvtGen.h b/src/RapidExternalEvtGen.h index 15c8914..1ee3a69 100644 --- a/src/RapidExternalEvtGen.h +++ b/src/RapidExternalEvtGen.h @@ -4,6 +4,8 @@ #include "TString.h" #include "RapidExternalGenerator.h" +#include "EvtGenBase/EvtParticle.hh" +#include "EvtGenBase/EvtHepMCEvent.hh" #ifdef RAPID_EVTGEN class EvtGen; @@ -28,6 +30,7 @@ class RapidExternalEvtGen : public RapidExternalGenerator { static TString getEvtGenConjName(int id); private: + double calcFlightTime( FourVector& BDecayVtx, FourVector& B4mtm ); #ifdef RAPID_EVTGEN EvtGen* evtGen_; #else @@ -37,4 +40,3 @@ class RapidExternalEvtGen : public RapidExternalGenerator { }; #endif - diff --git a/src/RapidParam.cc b/src/RapidParam.cc index fe0787b..9ac2c92 100644 --- a/src/RapidParam.cc +++ b/src/RapidParam.cc @@ -19,6 +19,11 @@ double RapidParam::eval() { double minip = particles_[0]->getMinIP(); double ipSigma = particles_[0]->getSigmaIP(); double minipSigma = particles_[0]->getSigmaMinIP(); + + double decaytime = particles_[0]->t(); + int id = particles_[0]->id(); + bool hasOsc = particles_[0]->hasOsc(); + mom_.SetPxPyPzE(0.,0.,0.,0.); if(truth_) { for(unsigned int i=0; imassHypothesisName() == "") id = particles_[0]->id(); +// else id = particleData->pdgCode(particles_[0]->massHypothesisName()); + +// return id; +// } +// double RapidParam::evalDecayTime() { + +// RapidParticleData * particleData = RapidParticleData::getInstance(); +// int decaytime(0); +// decaytime = particles_[0]->t(); +// return decaytime; +// } + double RapidParam::evalCorrectedMass() { TLorentzVector momS, momT; @@ -458,6 +498,12 @@ TString RapidParam::typeName() { return "origY"; case RapidParam::OrigZ: return "origZ"; + case RapidParam::ID: + return "ID"; + case RapidParam::decaytime: + return "decaytime"; + case RapidParam::hasOsc: + return "hasOsc"; default: std::cout << "WARNING in RapidParam::typeName : unknown type " << type_ << "." << std::endl << " returning empty string." << std::endl; @@ -534,6 +580,12 @@ RapidParam::ParamType RapidParam::typeFromString(TString str) { return RapidParam::OrigY; } else if(str=="origZ") { return RapidParam::OrigZ; + } else if(str=="ID") { + return RapidParam::ID; + } else if(str=="decaytime") { + return RapidParam::decaytime; + } else if(str=="hasOsc") { + return RapidParam::hasOsc; } else { std::cout << "WARNING in RapidParam::typeFromString : unknown type name " << str << "." << std::endl << " returning mass parameter type." << std::endl; @@ -754,4 +806,4 @@ void RapidParam::setMassMinMax(const std::vector& parts, double& if(min<0.) min=0.; } -} +} \ No newline at end of file diff --git a/src/RapidParam.h b/src/RapidParam.h index 6d07336..b122da2 100644 --- a/src/RapidParam.h +++ b/src/RapidParam.h @@ -48,6 +48,9 @@ class RapidParam { OrigX, // x coordinate of origin vertex OrigY, // y coordinate of origin vertex OrigZ, // z coordinate of origin vertex + ID, // TRUE ID of the particle + hasOsc, // whether the particle has oscillated before decay + decaytime, // decaytime of the particle UNKNOWN //unused }; @@ -101,6 +104,8 @@ class RapidParam { double evalCorrectedMass(); double evalTheta(); double evalPID(); + // double evalID(); + // double evalDecayTime(); double evalPos(); void setDefaultMinMax() {setDefaultMinMax(particles_,minVal_,maxVal_);}//TODO remove/make virtual and give a warning in the base class @@ -122,4 +127,4 @@ class RapidParam { TLorentzVector mom_; }; -#endif +#endif \ No newline at end of file diff --git a/src/RapidParticle.h b/src/RapidParticle.h index 577d834..cc5013d 100644 --- a/src/RapidParticle.h +++ b/src/RapidParticle.h @@ -55,6 +55,8 @@ class RapidParticle { double charge() { return charge_; } double minMass() { return minMass_; } double maxMass() { return maxMass_; } + double t() {return t_; } + bool hasOsc() {return hasOsc_;} RapidVertex * getOriginVertex() {return originVertex_;} RapidVertex * getDecayVertex() {return decayVertex_;} @@ -92,7 +94,9 @@ class RapidParticle { void setMinIPSigma(double sigma) { sigmaminip_ = sigma; } // void setPtEtaPhi(double pt, double eta, double phi) { p_.SetPtEtaPhiM(pt,eta,phi,mass_); } - + void setDecaytime(double t) { t_ = t; } + void setHasOsc(bool hasOsc) { hasOsc_ = hasOsc; } + void setId(int id) { id_ = id; } void print(int index); void setMassShape(RooDataSet* ds, double minMass, double maxMass, TString varName); @@ -101,6 +105,9 @@ class RapidParticle { TString evtGenDecayModel() { return evtGenModel_; } void setEvtGenDecayModel(TString value) { evtGenModel_ = value; } + TString evtGenDefinitions() { return evtGenDefs_; } + void setEvtGenDefinitions(TString value) { evtGenDefs_ = value; } + private: bool hasFlavour(int flavour); @@ -131,6 +138,8 @@ class RapidParticle { double sigmaminip_; double ipSmeared_; double minipSmeared_; + double t_; + bool hasOsc_; TLorentzVector p_; TLorentzVector pSmeared_; @@ -145,6 +154,7 @@ class RapidParticle { TString varName_; TString evtGenModel_; + TString evtGenDefs_; //store alternative mass hypotheses std::vector massHypothesisNames_; @@ -155,4 +165,4 @@ class RapidParticle { RapidVertex * originVertex_; RapidVertex * decayVertex_; }; -#endif +#endif \ No newline at end of file diff --git a/src/RapidParticleData.cc b/src/RapidParticleData.cc index 1b92544..cd1c5ef 100644 --- a/src/RapidParticleData.cc +++ b/src/RapidParticleData.cc @@ -15,6 +15,9 @@ RapidParticleData* RapidParticleData::getInstance() { if(!instance_) { instance_ = new RapidParticleData(); TString path; + path=getenv("RAPIDSIM_PARTICLE_TABLE"); + if(path!="") instance_->loadData(path); + path=getenv("RAPIDSIM_CONFIG"); if(path!="") instance_->loadData(path+"/config/particles.dat"); @@ -215,6 +218,7 @@ void RapidParticleData::addEntry(int id, TString name, double mass, double width if(idToName_.find(id) != idToName_.end()) { std::cout << "INFO in RapidParticleData::addEntry : particle with ID " << id << " already defined with name " << idToName_[id] << std::endl; std::cout << " second definition will be ignored." << std::endl; + return; } idToCT_[id] = ctau; @@ -465,4 +469,4 @@ void RapidParticleData::getMaxAltHypothesisMassShifts(RapidParticle* part, doubl getMaxAltHypothesisMassShifts(daug,deltaDown,deltaUp); daug = daug->next(); } -} +} \ No newline at end of file