CMS 3D CMS Logo
Functions
HadronDecayer.h File Reference
#include "DatabasePDG.h"
#include "Particle.h"

Go to the source code of this file.

Functions

void Decay (List_t &output, Particle &p, ParticleAllocator &allocator, DatabasePDG *database)
 
double GetDecayTime (const Particle &p, double weakDecayLimit)
 

Function Documentation

◆ Decay()

void Decay ( List_t &  output,
Particle &  p,
ParticleAllocator &  allocator,
DatabasePDG *  database 
)

Definition at line 56 of file HadronDecayer.cc.

56  {
57  // Get the PDG properties of the particle
58  ParticlePDG *pDef = parent.Def();
59 
60  // Get the number of posible decay channels
61  int nDecayChannel = pDef->GetNDecayChannels();
62  if (pDef->GetWidth() > 0 && nDecayChannel == 0) {
63  LogDebug("HadronDecayer") << "Decay(): WARNING!! Particle " << pDef->GetPDG() << " has finite width ("
64  << pDef->GetWidth()
65  << ") but NO decay channels specified in the database. Check it out!!";
66  return;
67  }
68 
69  // check the full branching of this specie
70  double fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays
71 
72  // return if particle has no branching
73  if (fullBranching < 0.00001)
74  return;
75 
76  // get the PDG mass of the specie
77  double PDGmass = pDef->GetMass();
78  int ComprCodePyth = 0;
79  float Delta = 0;
80 
81  bool success = kFALSE;
82  int iterations = 0;
83  // Try to decay the particle
84  while (!success) {
85  if (iterations > 1000) { //???
86  LogDebug("HadronDecayer") << "Decay(): WARNING!!! iterations to decay " << pDef->GetPDG() << " : " << iterations
87  << " Check it out!!!";
88  }
89 
90  // get a random mass using the Breit-Wigner distribution
91  double BWmass = CLHEP::RandBreitWigner::shoot(hjRandomEngine, PDGmass, pDef->GetWidth());
92 
93  // Try to cut the Breit Wigner tail of the particle using the cuts from pythia
94  // The Delta variable is obtained from pythia based on the specie
95  int encoding = pDef->GetPDG();
96  SERVICEEV.ipdg = encoding;
97  mydelta_();
98  ComprCodePyth = SERVICEEV.KC;
99  Delta = SERVICEEV.delta; // PYDAT2.PMAS[KC][3];
100 
101  //if there are no such particle in PYTHIA particle table, we take Delta=0.4
102  if (ComprCodePyth == 0) {
103  BWmass = PDGmass;
104  Delta = 0.0;
105  }
106 
107  //bad delta - an exception
108  if (ComprCodePyth == 254) {
109  BWmass = PDGmass;
110  Delta = 0.0;
111  }
112 
113  // K0 decay into K0s or K0l
114  if (TMath::Abs(encoding) == 311) {
115  BWmass = PDGmass;
116  Delta = 0.0;
117  }
118 
119  //for particles from PYTHIA table only, if the BW mass is outside the cut range then quit this iteration and generate another BW mass
120  if (ComprCodePyth != 0 && Delta > 0 && (BWmass < PDGmass - Delta || BWmass > PDGmass + Delta)) {
121  iterations++;
122  continue;
123  }
124 
125  if (BWmass > 5)
126  LogDebug("HadronDecayer") << "Decay(): Breit-Wigner mass > 5GeV for encoding: " << encoding
127  << "; PDG mass: " << PDGmass << "; delta: " << Delta << "; width: " << pDef->GetWidth()
128  << "; mass: " << BWmass << "; ComprCodePyth: " << ComprCodePyth;
129 
130  // check how many decay channels are allowed with the generated mass
131  int nAllowedChannels = database->GetNAllowedChannels(pDef, BWmass);
132  // if no decay channels are posible with this mass, then generate another BW mass
133  if (nAllowedChannels == 0) {
134  iterations++;
135  continue;
136  }
137 
138  std::vector<Particle> apDaughter;
139  std::vector<double> dMass; //daughters'mass
140  std::vector<double> dMom;
141  std::vector<double> sm;
142  std::vector<double> rd;
143 
144  // we need to choose an allowed decay channel
145  double randValue = CLHEP::RandFlat::shoot(hjRandomEngine) * fullBranching;
146 
147  int chosenChannel = 1000;
148  bool found = kFALSE;
149  int channelIterations = 0;
150  while (!found) {
151  if (channelIterations > 1000) {
152  LogDebug("HadronDecayer") << "Decay(): More than 1000 iterations to choose a decay channel. Check it out !!";
153  }
154  for (int nChannel = 0; nChannel < nDecayChannel; ++nChannel) {
155  randValue -= pDef->GetDecayChannel(nChannel)->GetBranching();
156  if (randValue <= 0. && database->IsChannelAllowed(pDef->GetDecayChannel(nChannel), BWmass)) {
157  chosenChannel = nChannel;
158  found = kTRUE;
159  break;
160  }
161  }
162  channelIterations++;
163  }
164 
165  // get the PDG information for the chosen decay channel
166  DecayChannel *dc = pDef->GetDecayChannel(chosenChannel);
167  int nSec = dc->GetNDaughters();
168 
169  // Adjust the parent momentum four-vector for the MC generated Breit-Wigner mass
170  Particle parentBW(database->GetPDGParticle(parent.Encoding()));
171  parentBW.Pos(parent.Pos());
172  double BWenergy = TMath::Sqrt(parent.Mom().X() * parent.Mom().X() + parent.Mom().Y() * parent.Mom().Y() +
173  parent.Mom().Z() * parent.Mom().Z() + BWmass * BWmass);
174 
175  int NB = (int)parent.GetType(); //particle from jets
176 
177  TLorentzVector MomparentBW(parent.Mom().X(), parent.Mom().Y(), parent.Mom().Z(), BWenergy);
178  parentBW.Mom(MomparentBW);
179  // take into account BW when calculating boost velocity (for wide resonances it matters)
180  TVector3 velocityBW(parentBW.Mom().BoostVector());
181 
182  // now we have an allowed decay
183  // first case: one daughter particle
184  if (nSec == 1) {
185  // initialize the daughter particle
186  Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
187  p1.Pos(parentBW.Pos());
188  p1.Mom(parent.Mom());
189  p1.SetLastMotherPdg(parentBW.Encoding());
190  p1.SetLastMotherDecayCoor(parentBW.Pos());
191  p1.SetLastMotherDecayMom(parentBW.Mom());
192  p1.SetType(NB);
193  p1.SetPythiaStatusCode(parent.GetPythiaStatusCode());
194 
195  // add the daughter particle to the list of secondaries
196  int parentIndex = parent.GetIndex();
197  int p1Index = p1.SetIndex();
198  p1.SetMother(parentIndex);
199  parent.SetFirstDaughterIndex(p1Index);
200  parent.SetLastDaughterIndex(p1Index);
201  allocator.AddParticle(p1, output);
202  success = kTRUE;
203  }
204  // second case: two daughter particles
205  else if (nSec == 2) {
206  // initialize the daughter particles
207  Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
208  p1.Pos(parentBW.Pos());
209  Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
210  p2.Pos(parentBW.Pos());
211 
212  // calculate the momenta in rest frame of mother for the two particles (theta and phi are isotropic)
213  MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), BWmass);
214 
215  // boost to the laboratory system (to the mother velocity)
216  p1.Mom().Boost(velocityBW);
217  p2.Mom().Boost(velocityBW);
218 
219  //store information about mother
220  p1.SetLastMotherPdg(parentBW.Encoding());
221  p1.SetLastMotherDecayCoor(parentBW.Pos());
222  p1.SetLastMotherDecayMom(parentBW.Mom());
223  p2.SetLastMotherPdg(parentBW.Encoding());
224  p2.SetLastMotherDecayCoor(parentBW.Pos());
225  p2.SetLastMotherDecayMom(parentBW.Mom());
226  //set to daughters the same type as has mother
227  p1.SetType(NB);
228  p2.SetType(NB);
229  p1.SetPythiaStatusCode(parent.GetPythiaStatusCode());
230  p2.SetPythiaStatusCode(parent.GetPythiaStatusCode());
231 
232  // check the kinematics in the lab system
233  double deltaS = TMath::Sqrt(
234  (parentBW.Mom().X() - p1.Mom().X() - p2.Mom().X()) * (parentBW.Mom().X() - p1.Mom().X() - p2.Mom().X()) +
235  (parentBW.Mom().Y() - p1.Mom().Y() - p2.Mom().Y()) * (parentBW.Mom().Y() - p1.Mom().Y() - p2.Mom().Y()) +
236  (parentBW.Mom().Z() - p1.Mom().Z() - p2.Mom().Z()) * (parentBW.Mom().Z() - p1.Mom().Z() - p2.Mom().Z()) +
237  (parentBW.Mom().E() - p1.Mom().E() - p2.Mom().E()) * (parentBW.Mom().E() - p1.Mom().E() - p2.Mom().E()));
238  // if deltaS is too big then repeat the kinematic procedure
239 
240  if (deltaS > 0.001) {
241  LogDebug("HadronDecayer|deltaS") << "2-body decay kinematic check in lab system: " << pDef->GetPDG() << " >>> "
242  << p1.Encoding() << " + " << p2.Encoding() << endl
243  << " Mother (e,px,py,pz): " << parentBW.Mom().E() << " : "
244  << parentBW.Mom().X() << " : " << parentBW.Mom().Y() << " : "
245  << parentBW.Mom().Z() << endl
246  << " Mother (x,y,z,t): " << parentBW.Pos().X() << " : "
247  << parentBW.Pos().Y() << " : " << parentBW.Pos().Z() << " : "
248  << parentBW.Pos().T() << endl
249  << " Daughter1 (e,px,py,pz): " << p1.Mom().E() << " : " << p1.Mom().X()
250  << " : " << p1.Mom().Y() << " : " << p1.Mom().Z() << endl
251  << " Daughter2 (e,px,py,pz): " << p2.Mom().E() << " : " << p2.Mom().X()
252  << " : " << p2.Mom().Y() << " : " << p2.Mom().Z() << endl
253  << " 2-body decay delta(sqrtS): " << deltaS << endl
254  << " Repeating the decay algorithm";
255 
256  iterations++;
257  continue;
258  }
259  // push particles to the list of secondaries
260  int parentIndex = parent.GetIndex();
261  p1.SetIndex();
262  p2.SetIndex();
263  p1.SetMother(parentIndex);
264  p2.SetMother(parentIndex);
265  parent.SetFirstDaughterIndex(p1.GetIndex());
266  parent.SetLastDaughterIndex(p2.GetIndex());
267  allocator.AddParticle(p1, output);
268  allocator.AddParticle(p2, output);
269  success = kTRUE;
270  }
271 
272  // third case: three daughter particle
273  else if (nSec == 3) {
274  // initialize the daughter particle
275  Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
276  p1.Pos(parentBW.Pos());
277  Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
278  p2.Pos(parentBW.Pos());
279  Particle p3(database->GetPDGParticle(dc->GetDaughterPDG(2)));
280  p3.Pos(parentBW.Pos());
281  // calculate the momenta in the rest frame of the mother particle
282  double pAbs1 = 0., pAbs2 = 0., pAbs3 = 0., sumPabs = 0., maxPabs = 0.;
283  double mass1 = p1.TableMass(), mass2 = p2.TableMass(), mass3 = p3.TableMass();
284  TLorentzVector &mom1 = p1.Mom(), &mom2 = p2.Mom(), &mom3 = p3.Mom();
285  double deltaMass = BWmass - mass1 - mass2 - mass3;
286 
287  do {
288  double rd1 = CLHEP::RandFlat::shoot(hjRandomEngine);
289  double rd2 = CLHEP::RandFlat::shoot(hjRandomEngine);
290  if (rd2 > rd1)
291  std::swap(rd1, rd2);
292  // 1
293  double e = rd2 * deltaMass;
294  pAbs1 = TMath::Sqrt(e * e + 2 * e * mass1);
295  sumPabs = pAbs1;
296  maxPabs = sumPabs;
297  // 2
298  e = (1 - rd1) * deltaMass;
299  pAbs2 = TMath::Sqrt(e * e + 2 * e * mass2);
300 
301  if (pAbs2 > maxPabs)
302  maxPabs = pAbs2;
303 
304  sumPabs += pAbs2;
305  // 3
306  e = (rd1 - rd2) * deltaMass;
307  pAbs3 = TMath::Sqrt(e * e + 2 * e * mass3);
308 
309  if (pAbs3 > maxPabs)
310  maxPabs = pAbs3;
311  sumPabs += pAbs3;
312  } while (maxPabs > sumPabs - maxPabs);
313 
314  // isotropic sample first particle 3-momentum
315  double cosTheta = 2 * (CLHEP::RandFlat::shoot(hjRandomEngine)) - 1;
316  double sinTheta = TMath::Sqrt(1 - cosTheta * cosTheta);
317  double phi = TMath::TwoPi() * (CLHEP::RandFlat::shoot(hjRandomEngine));
318  double sinPhi = TMath::Sin(phi);
319  double cosPhi = TMath::Cos(phi);
320 
321  mom1.SetPxPyPzE(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta, 0);
322  mom1 *= pAbs1;
323  // sample rest particle 3-momentum
324  double cosThetaN = (pAbs2 * pAbs2 - pAbs3 * pAbs3 - pAbs1 * pAbs1) / (2 * pAbs1 * pAbs3);
325  double sinThetaN = TMath::Sqrt(1 - cosThetaN * cosThetaN);
326  double phiN = TMath::TwoPi() * (CLHEP::RandFlat::shoot(hjRandomEngine));
327  double sinPhiN = TMath::Sin(phiN);
328  double cosPhiN = TMath::Cos(phiN);
329 
330  mom3.SetPxPyPzE(
331  sinThetaN * cosPhiN * cosTheta * cosPhi - sinThetaN * sinPhiN * sinPhi + cosThetaN * sinTheta * cosPhi,
332  sinThetaN * cosPhiN * cosTheta * sinPhi + sinThetaN * sinPhiN * cosPhi + cosThetaN * sinTheta * sinPhi,
333  -sinThetaN * cosPhiN * sinTheta + cosThetaN * cosTheta,
334  0.);
335 
336  mom3 *= pAbs3 * mom3.P();
337  mom2 = mom1;
338  mom2 += mom3;
339  mom2 *= -1.;
340  // calculate energy
341  mom1.SetE(TMath::Sqrt(mom1.P() * mom1.P() + mass1 * mass1));
342  mom2.SetE(TMath::Sqrt(mom2.P() * mom2.P() + mass2 * mass2));
343  mom3.SetE(TMath::Sqrt(mom3.P() * mom3.P() + mass3 * mass3));
344 
345  // boost to Lab system
346  mom1.Boost(velocityBW);
347  mom2.Boost(velocityBW);
348  mom3.Boost(velocityBW);
349 
350  p1.SetLastMotherPdg(parentBW.Encoding());
351  p1.SetLastMotherDecayCoor(parentBW.Pos());
352  p1.SetLastMotherDecayMom(parentBW.Mom());
353  p2.SetLastMotherPdg(parentBW.Encoding());
354  p2.SetLastMotherDecayCoor(parentBW.Pos());
355  p2.SetLastMotherDecayMom(parentBW.Mom());
356  p3.SetLastMotherPdg(parentBW.Encoding());
357  p3.SetLastMotherDecayCoor(parentBW.Pos());
358  p3.SetLastMotherDecayMom(parentBW.Mom());
359 
360  //set to daughters the same type as has mother
361  p1.SetType(NB);
362  p2.SetType(NB);
363  p3.SetType(NB);
364  p1.SetPythiaStatusCode(parent.GetPythiaStatusCode());
365  p2.SetPythiaStatusCode(parent.GetPythiaStatusCode());
366  p3.SetPythiaStatusCode(parent.GetPythiaStatusCode());
367 
368  // energy conservation check in the lab system
369  double deltaS = TMath::Sqrt((parentBW.Mom().X() - p1.Mom().X() - p2.Mom().X() - p3.Mom().X()) *
370  (parentBW.Mom().X() - p1.Mom().X() - p2.Mom().X() - p3.Mom().X()) +
371  (parentBW.Mom().Y() - p1.Mom().Y() - p2.Mom().Y() - p3.Mom().Y()) *
372  (parentBW.Mom().Y() - p1.Mom().Y() - p2.Mom().Y() - p3.Mom().Y()) +
373  (parentBW.Mom().Z() - p1.Mom().Z() - p2.Mom().Z() - p3.Mom().Z()) *
374  (parentBW.Mom().Z() - p1.Mom().Z() - p2.Mom().Z() - p3.Mom().Z()) +
375  (parentBW.Mom().E() - p1.Mom().E() - p2.Mom().E() - p3.Mom().E()) *
376  (parentBW.Mom().E() - p1.Mom().E() - p2.Mom().E() - p3.Mom().E()));
377  // if deltaS is too big then repeat the kinematic procedure
378  if (deltaS > 0.001) {
379  LogDebug("HadronDecayer|deltaS") << "3-body decay kinematic check in lab system: " << pDef->GetPDG() << " >>> "
380  << p1.Encoding() << " + " << p2.Encoding() << " + " << p3.Encoding() << endl
381  << "Mother (e,px,py,pz): " << parentBW.Mom().E() << " : "
382  << parentBW.Mom().X() << " : " << parentBW.Mom().Y() << " : "
383  << parentBW.Mom().Z() << endl
384  << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << " : " << p1.Mom().X() << " : "
385  << p1.Mom().Y() << " : " << p1.Mom().Z() << endl
386  << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << " : " << p2.Mom().X() << " : "
387  << p2.Mom().Y() << " : " << p2.Mom().Z() << endl
388  << "Daughter3 (e,px,py,pz): " << p3.Mom().E() << " : " << p3.Mom().X() << " : "
389  << p3.Mom().Y() << " : " << p3.Mom().Z() << endl
390  << "3-body decay delta(sqrtS): " << deltaS << endl
391  << "Repeating the decay algorithm";
392 
393  iterations++;
394  continue;
395  }
396 
397  // put particles in the list
398  int parentIndex = parent.GetIndex();
399  p1.SetIndex();
400  p2.SetIndex();
401  p3.SetIndex();
402  p1.SetMother(parentIndex);
403  p2.SetMother(parentIndex);
404  p3.SetMother(parentIndex);
405  parent.SetFirstDaughterIndex(p1.GetIndex());
406  parent.SetLastDaughterIndex(p3.GetIndex());
407  allocator.AddParticle(p1, output);
408  allocator.AddParticle(p2, output);
409  allocator.AddParticle(p3, output);
410  success = kTRUE;
411  }
412  }
413 
414  return;
415 }

References Abs(), ParticleAllocator::AddParticle(), MillePedeFileConverter_cfg::e, newFWLiteAna::found, DecayChannel::GetBranching(), DecayChannel::GetDaughterPDG(), ParticlePDG::GetDecayChannel(), ParticlePDG::GetFullBranching(), ParticlePDG::GetMass(), DatabasePDG::GetNAllowedChannels(), DecayChannel::GetNDaughters(), ParticlePDG::GetNDecayChannels(), ParticlePDG::GetPDG(), DatabasePDG::GetPDGParticle(), ParticlePDG::GetWidth(), hjRandomEngine, createfilelist::int, LogDebug, MomAntiMom(), mydelta_(), convertSQLitetoXML_cfg::output, p1, p2, p3, class-composition::parent, SERVICEEV, summarizeEdmComparisonLogfiles::success, std::swap(), and TwoPi.

Referenced by InitialState::Evolve().

◆ GetDecayTime()

double GetDecayTime ( const Particle &  p,
double  weakDecayLimit 
)

Definition at line 34 of file HadronDecayer.cc.

34  {
35  ParticlePDG *pDef = parent.Def();
36  double width = pDef->GetWidth(); //GeV
37 
38  double fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays
39  // check if full branching is 100%
40  if (pDef->GetNDecayChannels() > 0 && fullBranching < 0.9999) {
41  LogDebug("HadronDecayer") << "GetDecayTime(): WARNING!!! The full branching for specie " << pDef->GetPDG()
42  << " is less than 100% (" << fullBranching * 100 << "%) and decay channels exist!";
43  }
44 
45  if (width > weakDecayLimit) {
46  double slope = parent.E() * 0.1973 / (pDef->GetMass() * width);
47  return -slope * TMath::Log(CLHEP::RandFlat::shoot(hjRandomEngine)); //in fm/c
48  }
49 
50  return 0.;
51 }

References ParticlePDG::GetFullBranching(), ParticlePDG::GetMass(), ParticlePDG::GetNDecayChannels(), ParticlePDG::GetPDG(), ParticlePDG::GetWidth(), hjRandomEngine, LogDebug, class-composition::parent, slope, and ApeEstimator_cff::width.

Referenced by InitialState::Evolve().

ApeEstimator_cff.width
width
Definition: ApeEstimator_cff.py:24
TwoPi
const double TwoPi
Definition: CosmicMuonParameters.h:19
ParticlePDG::GetFullBranching
double GetFullBranching()
Definition: ParticlePDG.cc:49
ParticlePDG::GetDecayChannel
DecayChannel * GetDecayChannel(int i)
Definition: ParticlePDG.h:92
convertSQLitetoXML_cfg.output
output
Definition: convertSQLitetoXML_cfg.py:72
DecayChannel::GetDaughterPDG
int GetDaughterPDG(int i)
Definition: DecayChannel.cc:69
hjRandomEngine
CLHEP::HepRandomEngine * hjRandomEngine
Definition: Hydjet2Hadronizer.h:55
newFWLiteAna.found
found
Definition: newFWLiteAna.py:118
ParticlePDG::GetNDecayChannels
int GetNDecayChannels()
Definition: ParticlePDG.h:72
Abs
T Abs(T a)
Definition: MathUtil.h:49
std::swap
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
Definition: DataFrameContainer.h:209
DecayChannel
Definition: DecayChannel.h:22
p2
double p2[4]
Definition: TauolaWrapper.h:90
summarizeEdmComparisonLogfiles.success
success
Definition: summarizeEdmComparisonLogfiles.py:115
ParticlePDG
Definition: ParticlePDG.h:24
mydelta_
void mydelta_()
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
PVValHelper::phi
Definition: PVValidationHelpers.h:69
createfilelist.int
int
Definition: createfilelist.py:10
p1
double p1[4]
Definition: TauolaWrapper.h:89
ParticlePDG::GetMass
double GetMass()
Definition: ParticlePDG.h:70
ParticlePDG::GetWidth
double GetWidth()
Definition: ParticlePDG.h:71
SERVICEEV
#define SERVICEEV
Definition: HYJET_COMMONS.h:52
MomAntiMom
void MomAntiMom(TLorentzVector &mom, double mass, TLorentzVector &antiMom, double antiMass, double initialMass)
Definition: UKUtility.cc:45
DecayChannel::GetBranching
double GetBranching()
Definition: DecayChannel.h:40
Particle
Definition: Particle.py:1
p3
double p3[4]
Definition: TauolaWrapper.h:91
DecayChannel::GetNDaughters
int GetNDaughters()
Definition: DecayChannel.h:41
ParticlePDG::GetPDG
int GetPDG()
Definition: ParticlePDG.h:69
DatabasePDG::GetNAllowedChannels
int GetNAllowedChannels(ParticlePDG *particle, double motherMass)
Definition: DatabasePDG.cc:657
DatabasePDG::GetPDGParticle
ParticlePDG * GetPDGParticle(int pdg)
Definition: DatabasePDG.cc:218
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
class-composition.parent
parent
Definition: class-composition.py:88
ParticleAllocator::AddParticle
void AddParticle(const Particle &particle, List_t &list)
Definition: Particle.cc:107
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16