CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HadronDecayer.cc
Go to the documentation of this file.
1 /*
2 
3 July 2008 BW mass is limited by "PYTHIA method", by I. Lokhtin and L. Malinina
4 
5 Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
6 amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
7 November. 2, 2005
8 
9 */
10 
11 #include <functional>
12 #include <algorithm>
13 #include <vector>
14 #include <iostream>
15 #include <TMath.h>
16 
17 #include "GeneratorInterface/Hydjet2Interface/interface/DatabasePDG.h"
18 #include "GeneratorInterface/Hydjet2Interface/interface/ParticlePDG.h"
19 #include "GeneratorInterface/Hydjet2Interface/interface/DecayChannel.h"
20 #include "GeneratorInterface/Hydjet2Interface/interface/HadronDecayer.h"
21 #include "GeneratorInterface/Hydjet2Interface/interface/UKUtility.h"
22 #include "GeneratorInterface/Hydjet2Interface/interface/Particle.h"
23 #include "GeneratorInterface/Hydjet2Interface/interface/HYJET_COMMONS.h"
24 
25 //calculates decay time in fm/c
26 //calculates 1,2 and 3 body decays
27 
28 using namespace std;
29 
30 #include "CLHEP/Random/RandomEngine.h"
31 #include "CLHEP/Random/RandBreitWigner.h"
32 extern CLHEP::HepRandomEngine* hjRandomEngine;
33 
34 
35 double GetDecayTime(const Particle &parent, double weakDecayLimit) {
36  ParticlePDG *pDef = parent.Def();
37  double width = pDef->GetWidth(); //GeV
38 
39  double fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays
40  // check if full branching is 100%
41  if(pDef->GetNDecayChannels()>0 && fullBranching<0.9999) {
42  LogDebug("HadronDecayer") << "GetDecayTime(): WARNING!!! The full branching for specie " << pDef->GetPDG() << " is less than 100% ("
43  << fullBranching*100 << "%) and decay channels exist!";
44  }
45 
46  if(width > weakDecayLimit) {
47  double slope = parent.E() * 0.1973 / (pDef->GetMass() * width);
48  return -slope * TMath::Log(CLHEP::RandFlat::shoot(hjRandomEngine));//in fm/c
49  }
50 
51  return 0.;
52 }
53 
54 extern "C" void mydelta_();
55 extern SERVICEEVCommon SERVICEEV;
56 
57 void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, DatabasePDG* database) {
58  // Get the PDG properties of the particle
59  ParticlePDG *pDef = parent.Def();
60 
61  // Get the number of posible decay channels
62  int nDecayChannel = pDef->GetNDecayChannels();
63  if(pDef->GetWidth()>0 && nDecayChannel==0) {
64  LogDebug("HadronDecayer") << "Decay(): WARNING!! Particle " << pDef->GetPDG() << " has finite width (" << 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 "
87  << pDef->GetPDG() << " : " << iterations << " 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: "
128  << pDef->GetWidth() << "; 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() +
173  parent.Mom().Y()*parent.Mom().Y() +
174  parent.Mom().Z()*parent.Mom().Z() +
175  BWmass*BWmass);
176 
177  int NB = (int)parent.GetType(); //particle from jets
178 
179  TLorentzVector MomparentBW(parent.Mom().X(), parent.Mom().Y(), parent.Mom().Z(), BWenergy);
180  parentBW.Mom(MomparentBW);
181  // take into account BW when calculating boost velocity (for wide resonances it matters)
182  TVector3 velocityBW(parentBW.Mom().BoostVector());
183 
184  // now we have an allowed decay
185  // first case: one daughter particle
186  if(nSec == 1) {
187  // initialize the daughter particle
188  Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
189  p1.Pos(parentBW.Pos());
190  p1.Mom(parent.Mom());
191  p1.SetLastMotherPdg(parentBW.Encoding());
192  p1.SetLastMotherDecayCoor(parentBW.Pos());
193  p1.SetLastMotherDecayMom(parentBW.Mom());
194  p1.SetType(NB);
195  p1.SetPythiaStatusCode(parent.GetPythiaStatusCode());
196 
197  // add the daughter particle to the list of secondaries
198  int parentIndex = parent.GetIndex();
199  int p1Index = p1.SetIndex();
200  p1.SetMother(parentIndex);
201  parent.SetFirstDaughterIndex(p1Index);
202  parent.SetLastDaughterIndex(p1Index);
203  allocator.AddParticle(p1, output);
204  success = kTRUE;
205  }
206  // second case: two daughter particles
207  else if(nSec == 2) {
208  // initialize the daughter particles
209  Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
210  p1.Pos(parentBW.Pos());
211  Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
212  p2.Pos(parentBW.Pos());
213 
214  // calculate the momenta in rest frame of mother for the two particles (theta and phi are isotropic)
215  MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), BWmass);
216 
217  // boost to the laboratory system (to the mother velocity)
218  p1.Mom().Boost(velocityBW);
219  p2.Mom().Boost(velocityBW);
220 
221  //store information about mother
222  p1.SetLastMotherPdg(parentBW.Encoding());
223  p1.SetLastMotherDecayCoor(parentBW.Pos());
224  p1.SetLastMotherDecayMom(parentBW.Mom());
225  p2.SetLastMotherPdg(parentBW.Encoding());
226  p2.SetLastMotherDecayCoor(parentBW.Pos());
227  p2.SetLastMotherDecayMom(parentBW.Mom());
228  //set to daughters the same type as has mother
229  p1.SetType(NB);
230  p2.SetType(NB);
231  p1.SetPythiaStatusCode(parent.GetPythiaStatusCode());
232  p2.SetPythiaStatusCode(parent.GetPythiaStatusCode());
233 
234  // check the kinematics in the lab system
235  double deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())+
236  (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())+
237  (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())+
238  (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()));
239  // if deltaS is too big then repeat the kinematic procedure
240 
241 
242  if(deltaS>0.001) {
243 
244  LogDebug("HadronDecayer|deltaS") << "2-body decay kinematic check in lab system: " << pDef->GetPDG() << " >>> " << p1.Encoding() << " + " << p2.Encoding() << endl
245  << " Mother (e,px,py,pz): " << parentBW.Mom().E() << " : " << parentBW.Mom().X() << " : " << parentBW.Mom().Y() << " : " << parentBW.Mom().Z() << endl
246  << " Mother (x,y,z,t): " << parentBW.Pos().X() << " : " << parentBW.Pos().Y() << " : " << parentBW.Pos().Z() << " : " << parentBW.Pos().T() << endl
247  << " Daughter1 (e,px,py,pz): " << p1.Mom().E() << " : " << p1.Mom().X() << " : " << p1.Mom().Y() << " : " << p1.Mom().Z() << endl
248  << " Daughter2 (e,px,py,pz): " << p2.Mom().E() << " : " << p2.Mom().X() << " : " << p2.Mom().Y() << " : " << p2.Mom().Z() << endl
249  << " 2-body decay delta(sqrtS): " << deltaS << endl
250  << " Repeating the decay algorithm";
251 
252  iterations++;
253  continue;
254  }
255  // push particles to the list of secondaries
256  int parentIndex = parent.GetIndex();
257  p1.SetIndex();
258  p2.SetIndex();
259  p1.SetMother(parentIndex);
260  p2.SetMother(parentIndex);
261  parent.SetFirstDaughterIndex(p1.GetIndex());
262  parent.SetLastDaughterIndex(p2.GetIndex());
263  allocator.AddParticle(p1, output);
264  allocator.AddParticle(p2, output);
265  success = kTRUE;
266  }
267 
268  // third case: three daughter particle
269  else if(nSec == 3) {
270  // initialize the daughter particle
271  Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
272  p1.Pos(parentBW.Pos());
273  Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
274  p2.Pos(parentBW.Pos());
275  Particle p3(database->GetPDGParticle(dc->GetDaughterPDG(2)));
276  p3.Pos(parentBW.Pos());
277  // calculate the momenta in the rest frame of the mother particle
278  double pAbs1 = 0., pAbs2 = 0., pAbs3 = 0., sumPabs = 0., maxPabs = 0.;
279  double mass1 = p1.TableMass(), mass2 = p2.TableMass(), mass3 = p3.TableMass();
280  TLorentzVector &mom1 = p1.Mom(), &mom2 = p2.Mom(), &mom3 = p3.Mom();
281  double deltaMass = BWmass - mass1 - mass2 - mass3;
282 
283  do {
284  double rd1 = CLHEP::RandFlat::shoot(hjRandomEngine);
285  double rd2 = CLHEP::RandFlat::shoot(hjRandomEngine);
286  if (rd2 > rd1)
287  std::swap(rd1, rd2);
288  // 1
289  double e = rd2*deltaMass;
290  pAbs1 = TMath::Sqrt(e*e + 2*e*mass1);
291  sumPabs = pAbs1;
292  maxPabs = sumPabs;
293  // 2
294  e = (1-rd1)*deltaMass;
295  pAbs2 = TMath::Sqrt(e*e + 2*e*mass2);
296 
297  if(pAbs2 > maxPabs)
298  maxPabs = pAbs2;
299 
300  sumPabs += pAbs2;
301  // 3
302  e = (rd1-rd2)*deltaMass;
303  pAbs3 = TMath::Sqrt(e*e + 2*e*mass3);
304 
305  if (pAbs3 > maxPabs)
306  maxPabs = pAbs3;
307  sumPabs += pAbs3;
308  } while(maxPabs > sumPabs - maxPabs);
309 
310  // isotropic sample first particle 3-momentum
311  double cosTheta = 2*(CLHEP::RandFlat::shoot(hjRandomEngine)) - 1;
312  double sinTheta = TMath::Sqrt(1 - cosTheta*cosTheta);
313  double phi = TMath::TwoPi()*(CLHEP::RandFlat::shoot(hjRandomEngine));
314  double sinPhi = TMath::Sin(phi);
315  double cosPhi = TMath::Cos(phi);
316 
317  mom1.SetPxPyPzE(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta, 0);
318  mom1 *= pAbs1;
319  // sample rest particle 3-momentum
320  double cosThetaN = (pAbs2*pAbs2 - pAbs3*pAbs3 - pAbs1*pAbs1)/(2*pAbs1*pAbs3);
321  double sinThetaN = TMath::Sqrt(1 - cosThetaN*cosThetaN);
322  double phiN = TMath::TwoPi()*(CLHEP::RandFlat::shoot(hjRandomEngine));
323  double sinPhiN = TMath::Sin(phiN);
324  double cosPhiN = TMath::Cos(phiN);
325 
326  mom3.SetPxPyPzE(sinThetaN*cosPhiN*cosTheta*cosPhi - sinThetaN*sinPhiN*sinPhi + cosThetaN*sinTheta*cosPhi,
327  sinThetaN*cosPhiN*cosTheta*sinPhi + sinThetaN*sinPhiN*cosPhi + cosThetaN*sinTheta*sinPhi,
328  -sinThetaN*cosPhiN*sinTheta + cosThetaN*cosTheta,
329  0.);
330 
331  mom3 *= pAbs3*mom3.P();
332  mom2 = mom1;
333  mom2 += mom3;
334  mom2 *= -1.;
335  // calculate energy
336  mom1.SetE(TMath::Sqrt(mom1.P()*mom1.P() + mass1*mass1));
337  mom2.SetE(TMath::Sqrt(mom2.P()*mom2.P() + mass2*mass2));
338  mom3.SetE(TMath::Sqrt(mom3.P()*mom3.P() + mass3*mass3));
339 
340  // boost to Lab system
341  mom1.Boost(velocityBW);
342  mom2.Boost(velocityBW);
343  mom3.Boost(velocityBW);
344 
345  p1.SetLastMotherPdg(parentBW.Encoding());
346  p1.SetLastMotherDecayCoor(parentBW.Pos());
347  p1.SetLastMotherDecayMom(parentBW.Mom());
348  p2.SetLastMotherPdg(parentBW.Encoding());
349  p2.SetLastMotherDecayCoor(parentBW.Pos());
350  p2.SetLastMotherDecayMom(parentBW.Mom());
351  p3.SetLastMotherPdg(parentBW.Encoding());
352  p3.SetLastMotherDecayCoor(parentBW.Pos());
353  p3.SetLastMotherDecayMom(parentBW.Mom());
354 
355  //set to daughters the same type as has mother
356  p1.SetType(NB);
357  p2.SetType(NB);
358  p3.SetType(NB);
359  p1.SetPythiaStatusCode(parent.GetPythiaStatusCode());
360  p2.SetPythiaStatusCode(parent.GetPythiaStatusCode());
361  p3.SetPythiaStatusCode(parent.GetPythiaStatusCode());
362 
363  // energy conservation check in the lab system
364  double deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X()-p3.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X()-p3.Mom().X()) +
365  (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y()) +
366  (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z()) +
367  (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E()));
368  // if deltaS is too big then repeat the kinematic procedure
369  if(deltaS>0.001) {
370 
371  LogDebug("HadronDecayer|deltaS") << "3-body decay kinematic check in lab system: " << pDef->GetPDG() << " >>> " << p1.Encoding() << " + " << p2.Encoding() << " + " << p3.Encoding() << endl
372  << "Mother (e,px,py,pz): " << parentBW.Mom().E() << " : " << parentBW.Mom().X() << " : " << parentBW.Mom().Y() << " : " << parentBW.Mom().Z() << endl
373  << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << " : " << p1.Mom().X() << " : " << p1.Mom().Y() << " : " << p1.Mom().Z() << endl
374  << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << " : " << p2.Mom().X() << " : " << p2.Mom().Y() << " : " << p2.Mom().Z() << endl
375  << "Daughter3 (e,px,py,pz): " << p3.Mom().E() << " : " << p3.Mom().X() << " : " << p3.Mom().Y() << " : " << p3.Mom().Z() << endl
376  << "3-body decay delta(sqrtS): " << deltaS << endl
377  << "Repeating the decay algorithm";
378 
379  iterations++;
380  continue;
381  }
382 
383  // put particles in the list
384  int parentIndex = parent.GetIndex();
385  p1.SetIndex();
386  p2.SetIndex();
387  p3.SetIndex();
388  p1.SetMother(parentIndex);
389  p2.SetMother(parentIndex);
390  p3.SetMother(parentIndex);
391  parent.SetFirstDaughterIndex(p1.GetIndex());
392  parent.SetLastDaughterIndex(p3.GetIndex());
393  allocator.AddParticle(p1, output);
394  allocator.AddParticle(p2, output);
395  allocator.AddParticle(p3, output);
396  success = kTRUE;
397  }
398  }
399 
400  return;
401 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
TwoPi
const double TwoPi
Definition: CosmicMuonParameters.h:19
DatabasePDG::GetPDGParticle
ParticlePDG * GetPDGParticle(int pdg)
Definition: DatabasePDG.cc:222
create_public_lumi_plots.width
int width
Definition: create_public_lumi_plots.py:1009
DecayChannel::GetNDaughters
int GetNDaughters()
Definition: DecayChannel.h:41
ParticlePDG::GetPDG
int GetPDG()
Definition: ParticlePDG.h:67
DecayChannel::GetBranching
double GetBranching()
Definition: DecayChannel.h:40
confdb.parent
parent
Definition: confdb.py:1052
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
HYJET_COMMONS.h
ParticlePDG::GetDecayChannel
DecayChannel * GetDecayChannel(int i)
Definition: ParticlePDG.h:87
SERVICEEV
#define SERVICEEV
Definition: HYJET_COMMONS.h:53
List_t
std::list< Particle > List_t
Definition: Particle.h:141
mydelta_
void mydelta_()
ParticlePDG::GetWidth
double GetWidth()
Definition: ParticlePDG.h:69
std::swap
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
Definition: DataFrameContainer.h:234
Decay
void Decay(List_t &output, Particle &p, ParticleAllocator &allocator, DatabasePDG *database)
Definition: HadronDecayer.cc:57
ParticlePDG::GetNDecayChannels
int GetNDecayChannels()
Definition: ParticlePDG.h:70
Abs
T Abs(T a)
Definition: MathUtil.h:49
DecayChannel::GetDaughterPDG
int GetDaughterPDG(int i)
Definition: DecayChannel.cc:64
p2
double p2[4]
Definition: TauolaWrapper.h:90
hjRandomEngine
CLHEP::HepRandomEngine * hjRandomEngine
Definition: Hydjet2Hadronizer.h:55
convertSQLitetoXML_cfg.output
tuple output
Definition: convertSQLitetoXML_cfg.py:32
phi
Geom::Phi< T > phi() const
Definition: Basic3DVectorLD.h:163
MomAntiMom
void MomAntiMom(TLorentzVector &mom, double mass, TLorentzVector &antiMom, double antiMass, double initialMass)
Definition: UKUtility.cc:45
GetDecayTime
double GetDecayTime(const Particle &p, double weakDecayLimit)
Definition: HadronDecayer.cc:35
p1
double p1[4]
Definition: TauolaWrapper.h:89
ParticlePDG::GetMass
double GetMass()
Definition: ParticlePDG.h:68
ParticleAllocator::AddParticle
void AddParticle(const Particle &particle, List_t &list)
Definition: Particle.cc:109
ParticlePDG::GetFullBranching
double GetFullBranching()
Definition: ParticlePDG.cc:47
DatabasePDG::GetNAllowedChannels
int GetNAllowedChannels(ParticlePDG *particle, double motherMass)
Definition: DatabasePDG.cc:650
p3
double p3[4]
Definition: TauolaWrapper.h:91
Particle::Particle
Particle(ParticlePDG *pdg=0)
Definition: Particle.cc:13
Generated for CMSSW Reference Manual by doxygen 1.8.5