1 /////////////////////////////////////////////////////////////////////////////////////
3 // July 2008 BW mass is limited by "PYTHIA method", by I. Lokhtin and L. Malinina //
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 //
9 /////////////////////////////////////////////////////////////////////////////////////
19 #include "DatabasePDG.h"
20 #include "ParticlePDG.h"
21 #include "DecayChannel.h"
22 #include "UKUtility.h"
24 #include "HYJET_COMMONS.h"
25 #include "HadronDecayer.h"
27 //calculates decay time in fm/c
28 //calculates 1,2 and 3 body decays
33 Double_t GetDecayTime(const Particle &parent, Double_t weakDecayLimit) {
35 // return a random decay time according to the particle's width
38 ParticlePDG *pDef = parent.Def();
39 Double_t width = pDef->GetWidth(); //GeV
41 // if particle is set to be stable then return 0
42 if(pDef->GetStableStatus())
45 if(width > weakDecayLimit && weakDecayLimit>=0.0) {
46 const Double_t slope = parent.E() * 0.1973 / (pDef->GetMass() * width);
47 return -slope * TMath::Log(gRandom->Rndm());//in fm/c
54 extern "C" void mydelta_();
55 extern SERVICEEVCommon SERVICEEV;
57 void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, const DatabasePDG* database) {
62 // Get the PDG properties of the particle
63 ParticlePDG *pDef = parent.Def();
65 // Get the number of posible decay channels
66 Int_t nDecayChannel = pDef->GetNDecayChannels();
68 // check the full branching of this specie
69 Double_t fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays
71 // return if particle has no branching
72 if(fullBranching < 0.00001)
75 // get the PDG mass of the specie
76 Double_t pdgMass = pDef->GetMass();
77 Int_t comprCodePyth=0;
80 Bool_t success = kFALSE;
82 // Try to decay the particle
84 // std::cout << "HadronDecayer::Decay() iteration #" << iterations << std::endl;
86 std::cout << "HadronDecayer::Decay() more than 1000 iterations to decay particle with code "
87 << pDef->GetPDG() << std::endl;
88 std::cout << " Will be left undecayed ... check it out!" << std::endl;
92 // get a random mass using the Breit Wigner distribution
93 Double_t bwMass = gRandom->BreitWigner(pdgMass, pDef->GetWidth());
95 // Try to cut the Breit Wigner tail of the particle using the cuts from pythia
96 // The delta variable is obtained from pythia based on the specie
97 int encoding =pDef->GetPDG();
98 SERVICEEV.ipdg = encoding;
100 comprCodePyth=SERVICEEV.KC;
101 delta = SERVICEEV.delta;// PYDAT2.PMAS[KC][3];
103 //if there are no such particle in PYTHIA particle table, we take delta=0.4
104 if(comprCodePyth==0){
109 //bad delta - an exception
110 if(comprCodePyth==254){
115 // K0 decay into K0s or K0l
116 if(TMath::Abs(encoding)==311) {
121 //for particles from PYTHIA table only, if the BW mass is outside the cut range then quit this iteration and generate another BW mass
122 if(comprCodePyth!=0 && delta>0 && (bwMass<pdgMass-delta || bwMass>pdgMass+delta)){
127 // check how many decay channels are allowed with the generated mass
128 Int_t nAllowedChannels = database->GetNAllowedChannels(pDef, bwMass);
129 // if no decay channels are posible with this mass, then generate another BW mass
130 if(nAllowedChannels==0) {
132 std::cout << "HadronDecayer::Decay() no decays allowed at this BW mass" << std::endl;
136 std::vector<Particle> apDaughter;
137 std::vector<Double_t> dMass; //daughters'mass
138 std::vector<Double_t> dMom;
139 std::vector<Double_t> sm;
140 std::vector<Double_t> rd;
142 // we need to choose an allowed decay channel
143 Double_t randValue = gRandom->Rndm() * fullBranching;
144 Int_t chosenChannel = 1000;
145 Bool_t found = kFALSE;
146 Int_t channelIterations = 0;
148 for(Int_t nChannel = 0; nChannel < nDecayChannel; nChannel++) {
149 randValue -= pDef->GetDecayChannel(nChannel)->GetBranching();
150 if(randValue <= 0. && database->IsChannelAllowed(pDef->GetDecayChannel(nChannel), bwMass)) {
151 chosenChannel = nChannel;
159 // get the PDG information for the chosen decay channel
160 DecayChannel *dc = pDef->GetDecayChannel(chosenChannel);
161 Int_t nSec = dc->GetNDaughters();
163 // Adjust the parent momentum four-vector for the MC generated Breit-Wigner mass
164 Particle parentBW(database->GetPDGParticle(parent.Encoding()));
165 parentBW.Pos(parent.Pos());
166 Double_t bwEnergy = TMath::Sqrt(parent.Mom().X()*parent.Mom().X() +
167 parent.Mom().Y()*parent.Mom().Y() +
168 parent.Mom().Z()*parent.Mom().Z() +
171 Int_t nb = (Int_t)parent.GetType(); //particle from jets
173 TLorentzVector MomparentBW(parent.Mom().X(), parent.Mom().Y(), parent.Mom().Z(), bwEnergy);
174 parentBW.Mom(MomparentBW);
175 // take into account BW when calculating boost velocity (for wide resonances it matters)
176 TVector3 velocityBW(parentBW.Mom().BoostVector());
178 // now we have an allowed decay
179 // first case: one daughter particle
181 // initialize the daughter particle
182 Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
183 p1.Pos(parentBW.Pos());
184 p1.Mom(parent.Mom());
185 p1.SetLastMotherPdg(parentBW.Encoding());
186 p1.SetLastMotherDecayCoor(parentBW.Pos());
187 p1.SetLastMotherDecayMom(parentBW.Mom());
190 // add the daughter to the list of secondaries
191 Int_t parentIndex = parent.GetIndex();
192 Int_t p1Index = p1.SetIndex(); // set the daughter index
193 p1.SetMother(parentIndex); // set the mother index for this daughter
194 parent.SetFirstDaughterIndex(p1Index);
195 parent.SetLastDaughterIndex(p1Index);
196 allocator.AddParticle(p1, output);
199 // second case: two daughter particles
201 // initialize the daughter particles
202 Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
203 p1.Pos(parentBW.Pos());
204 Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
205 p2.Pos(parentBW.Pos());
207 // calculate the momenta in rest frame of mother for the two particles (theta and phi are isotropic)
208 MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), bwMass);
210 // boost to the laboratory system (to the mother velocity)
211 p1.Mom().Boost(velocityBW);
212 p2.Mom().Boost(velocityBW);
214 //store information about mother
215 p1.SetLastMotherPdg(parentBW.Encoding());
216 p1.SetLastMotherDecayCoor(parentBW.Pos());
217 p1.SetLastMotherDecayMom(parentBW.Mom());
218 p2.SetLastMotherPdg(parentBW.Encoding());
219 p2.SetLastMotherDecayCoor(parentBW.Pos());
220 p2.SetLastMotherDecayMom(parentBW.Mom());
221 //set to daughters the same type as has mother
226 // check the kinematics in the lab system
227 Double_t deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())+
228 (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())+
229 (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())+
230 (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()));
231 // if deltaS is too big then repeat the kinematic procedure
237 // push particles to the list of secondaries
238 Int_t parentIndex = parent.GetIndex();
241 p1.SetMother(parentIndex);
242 p2.SetMother(parentIndex);
243 parent.SetFirstDaughterIndex(p1.GetIndex());
244 parent.SetLastDaughterIndex(p2.GetIndex());
245 allocator.AddParticle(p1, output);
246 allocator.AddParticle(p2, output);
251 // third case: three daughter particle
253 // initialize the daughter particle
254 Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
255 p1.Pos(parentBW.Pos());
256 Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
257 p2.Pos(parentBW.Pos());
258 Particle p3(database->GetPDGParticle(dc->GetDaughterPDG(2)));
259 p3.Pos(parentBW.Pos());
260 // calculate the momenta in the rest frame of the mother particle
261 Double_t pAbs1 = 0., pAbs2 = 0., pAbs3 = 0., sumPabs = 0., maxPabs = 0.;
262 Double_t mass1 = p1.TableMass(), mass2 = p2.TableMass(), mass3 = p3.TableMass();
263 TLorentzVector &mom1 = p1.Mom(), &mom2 = p2.Mom(), &mom3 = p3.Mom();
264 Double_t deltaMass = bwMass - mass1 - mass2 - mass3;
267 Double_t rd1 = gRandom->Rndm();
268 Double_t rd2 = gRandom->Rndm();
272 Double_t e = rd2*deltaMass;
273 pAbs1 = TMath::Sqrt(e*e + 2*e*mass1);
277 e = (1-rd1)*deltaMass;
278 pAbs2 = TMath::Sqrt(e*e + 2*e*mass2);
285 e = (rd1-rd2)*deltaMass;
286 pAbs3 = TMath::Sqrt(e*e + 2*e*mass3);
291 } while(maxPabs > sumPabs - maxPabs);
293 // isotropic sample first particle 3-momentum
294 Double_t cosTheta = 2*(gRandom->Rndm()) - 1;
295 Double_t sinTheta = TMath::Sqrt(1 - cosTheta*cosTheta);
296 Double_t phi = TMath::TwoPi()*(gRandom->Rndm());
297 Double_t sinPhi = TMath::Sin(phi);
298 Double_t cosPhi = TMath::Cos(phi);
300 mom1.SetPxPyPzE(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta, 0);
302 // sample rest particle 3-momentum
303 Double_t cosThetaN = (pAbs2*pAbs2 - pAbs3*pAbs3 - pAbs1*pAbs1)/(2*pAbs1*pAbs3);
304 Double_t sinThetaN = TMath::Sqrt(1 - cosThetaN*cosThetaN);
305 Double_t phiN = TMath::TwoPi()*(gRandom->Rndm());
306 Double_t sinPhiN = TMath::Sin(phiN);
307 Double_t cosPhiN = TMath::Cos(phiN);
309 mom3.SetPxPyPzE(sinThetaN*cosPhiN*cosTheta*cosPhi - sinThetaN*sinPhiN*sinPhi + cosThetaN*sinTheta*cosPhi,
310 sinThetaN*cosPhiN*cosTheta*sinPhi + sinThetaN*sinPhiN*cosPhi + cosThetaN*sinTheta*sinPhi,
311 -sinThetaN*cosPhiN*sinTheta + cosThetaN*cosTheta,
314 mom3 *= pAbs3*mom3.P();
319 mom1.SetE(TMath::Sqrt(mom1.P()*mom1.P() + mass1*mass1));
320 mom2.SetE(TMath::Sqrt(mom2.P()*mom2.P() + mass2*mass2));
321 mom3.SetE(TMath::Sqrt(mom3.P()*mom3.P() + mass3*mass3));
323 // boost to Lab system
324 mom1.Boost(velocityBW);
325 mom2.Boost(velocityBW);
326 mom3.Boost(velocityBW);
328 p1.SetLastMotherPdg(parentBW.Encoding());
329 p1.SetLastMotherDecayCoor(parentBW.Pos());
330 p1.SetLastMotherDecayMom(parentBW.Mom());
331 p2.SetLastMotherPdg(parentBW.Encoding());
332 p2.SetLastMotherDecayCoor(parentBW.Pos());
333 p2.SetLastMotherDecayMom(parentBW.Mom());
334 p3.SetLastMotherPdg(parentBW.Encoding());
335 p3.SetLastMotherDecayCoor(parentBW.Pos());
336 p3.SetLastMotherDecayMom(parentBW.Mom());
338 //set to daughters the same type as has mother
343 // energy conservation check in the lab system
344 Double_t 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()) +
345 (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y()) +
346 (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z()) +
347 (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E()));
348 // if deltaS is too big then repeat the kinematic procedure
350 // cout << "3-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << " + " << p3.Encoding() << endl;
351 // cout << "Mother (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
352 // cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
353 // cout << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << "\t" << p2.Mom().X() << "\t" << p2.Mom().Y() << "\t" << p2.Mom().Z() << endl;
354 // cout << "Daughter3 (e,px,py,pz): " << p3.Mom().E() << "\t" << p3.Mom().X() << "\t" << p3.Mom().Y() << "\t" << p3.Mom().Z() << endl;
355 // cout << "3-body decay delta(sqrtS) = " << deltaS << endl;
356 // cout << "Repeating the decay algorithm..." << endl;
362 Int_t parentIndex = parent.GetIndex();
366 p1.SetMother(parentIndex);
367 p2.SetMother(parentIndex);
368 p3.SetMother(parentIndex);
369 parent.SetFirstDaughterIndex(p1.GetIndex());
370 parent.SetLastDaughterIndex(p3.GetIndex());
371 allocator.AddParticle(p1, output);
372 allocator.AddParticle(p2, output);
373 allocator.AddParticle(p3, output);