30 #include "CLHEP/Random/RandomEngine.h"
31 #include "CLHEP/Random/RandBreitWigner.h"
41 LogDebug(
"HadronDecayer") <<
"GetDecayTime(): WARNING!!! The full branching for specie " << pDef->
GetPDG()
42 <<
" is less than 100% (" << fullBranching * 100 <<
"%) and decay channels exist!";
45 if (
width > weakDecayLimit) {
62 if (pDef->
GetWidth() > 0 && nDecayChannel == 0) {
63 LogDebug(
"HadronDecayer") <<
"Decay(): WARNING!! Particle " << pDef->
GetPDG() <<
" has finite width ("
65 <<
") but NO decay channels specified in the database. Check it out!!";
73 if (fullBranching < 0.00001)
77 double PDGmass = pDef->
GetMass();
78 int ComprCodePyth = 0;
85 if (iterations > 1000) {
86 LogDebug(
"HadronDecayer") <<
"Decay(): WARNING!!! iterations to decay " << pDef->
GetPDG() <<
" : " << iterations
87 <<
" Check it out!!!";
95 int encoding = pDef->
GetPDG();
102 if (ComprCodePyth == 0) {
108 if (ComprCodePyth == 254) {
120 if (ComprCodePyth != 0 && Delta > 0 && (BWmass < PDGmass - Delta || BWmass > PDGmass + Delta)) {
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;
133 if (nAllowedChannels == 0) {
138 std::vector<Particle> apDaughter;
139 std::vector<double> dMass;
140 std::vector<double> dMom;
141 std::vector<double> sm;
142 std::vector<double> rd;
145 double randValue = CLHEP::RandFlat::shoot(
hjRandomEngine) * fullBranching;
147 int chosenChannel = 1000;
149 int channelIterations = 0;
151 if (channelIterations > 1000) {
152 LogDebug(
"HadronDecayer") <<
"Decay(): More than 1000 iterations to choose a decay channel. Check it out !!";
154 for (
int nChannel = 0; nChannel < nDecayChannel; ++nChannel) {
156 if (randValue <= 0. && database->IsChannelAllowed(pDef->
GetDecayChannel(nChannel), BWmass)) {
157 chosenChannel = nChannel;
171 parentBW.Pos(
parent.Pos());
177 TLorentzVector MomparentBW(
parent.Mom().X(),
parent.Mom().Y(),
parent.Mom().Z(), BWenergy);
178 parentBW.Mom(MomparentBW);
180 TVector3 velocityBW(parentBW.Mom().BoostVector());
187 p1.Pos(parentBW.Pos());
189 p1.SetLastMotherPdg(parentBW.Encoding());
190 p1.SetLastMotherDecayCoor(parentBW.Pos());
191 p1.SetLastMotherDecayMom(parentBW.Mom());
193 p1.SetPythiaStatusCode(
parent.GetPythiaStatusCode());
196 int parentIndex =
parent.GetIndex();
197 int p1Index =
p1.SetIndex();
198 p1.SetMother(parentIndex);
199 parent.SetFirstDaughterIndex(p1Index);
200 parent.SetLastDaughterIndex(p1Index);
205 else if (nSec == 2) {
208 p1.Pos(parentBW.Pos());
210 p2.Pos(parentBW.Pos());
216 p1.Mom().Boost(velocityBW);
217 p2.Mom().Boost(velocityBW);
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());
229 p1.SetPythiaStatusCode(
parent.GetPythiaStatusCode());
230 p2.SetPythiaStatusCode(
parent.GetPythiaStatusCode());
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()));
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";
260 int parentIndex =
parent.GetIndex();
263 p1.SetMother(parentIndex);
264 p2.SetMother(parentIndex);
265 parent.SetFirstDaughterIndex(
p1.GetIndex());
266 parent.SetLastDaughterIndex(
p2.GetIndex());
273 else if (nSec == 3) {
276 p1.Pos(parentBW.Pos());
278 p2.Pos(parentBW.Pos());
280 p3.Pos(parentBW.Pos());
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;
293 double e = rd2 * deltaMass;
294 pAbs1 = TMath::Sqrt(
e *
e + 2 *
e * mass1);
298 e = (1 - rd1) * deltaMass;
299 pAbs2 = TMath::Sqrt(
e *
e + 2 *
e * mass2);
306 e = (rd1 - rd2) * deltaMass;
307 pAbs3 = TMath::Sqrt(
e *
e + 2 *
e * mass3);
312 }
while (maxPabs > sumPabs - maxPabs);
315 double cosTheta = 2 * (CLHEP::RandFlat::shoot(
hjRandomEngine)) - 1;
316 double sinTheta = TMath::Sqrt(1 - cosTheta * cosTheta);
318 double sinPhi = TMath::Sin(phi);
319 double cosPhi = TMath::Cos(phi);
321 mom1.SetPxPyPzE(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta, 0);
324 double cosThetaN = (pAbs2 * pAbs2 - pAbs3 * pAbs3 - pAbs1 * pAbs1) / (2 * pAbs1 * pAbs3);
325 double sinThetaN = TMath::Sqrt(1 - cosThetaN * cosThetaN);
327 double sinPhiN = TMath::Sin(phiN);
328 double cosPhiN = TMath::Cos(phiN);
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,
336 mom3 *= pAbs3 * mom3.P();
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));
346 mom1.Boost(velocityBW);
347 mom2.Boost(velocityBW);
348 mom3.Boost(velocityBW);
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());
364 p1.SetPythiaStatusCode(
parent.GetPythiaStatusCode());
365 p2.SetPythiaStatusCode(
parent.GetPythiaStatusCode());
366 p3.SetPythiaStatusCode(
parent.GetPythiaStatusCode());
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()));
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";
398 int parentIndex =
parent.GetIndex();
402 p1.SetMother(parentIndex);
403 p2.SetMother(parentIndex);
404 p3.SetMother(parentIndex);
405 parent.SetFirstDaughterIndex(
p1.GetIndex());
406 parent.SetLastDaughterIndex(
p3.GetIndex());