]>
Commit | Line | Data |
---|---|---|
03896fc4 | 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 | ///////////////////////////////////////////////////////////////////////////////////// | |
b1c2e580 | 10 | |
11 | #include <functional> | |
12 | #include <algorithm> | |
13 | #include <vector> | |
14 | #include <iostream> | |
15 | #include <TRandom.h> | |
16 | #include <TError.h> | |
17 | #include <TMath.h> | |
18 | ||
b1c2e580 | 19 | #include "DatabasePDG.h" |
b1c2e580 | 20 | #include "ParticlePDG.h" |
b1c2e580 | 21 | #include "DecayChannel.h" |
b1c2e580 | 22 | #include "UKUtility.h" |
b1c2e580 | 23 | #include "Particle.h" |
b1c2e580 | 24 | #include "HYJET_COMMONS.h" |
03896fc4 | 25 | #include "HadronDecayer.h" |
b1c2e580 | 26 | |
27 | //calculates decay time in fm/c | |
28 | //calculates 1,2 and 3 body decays | |
29 | ||
30 | using std::cout; | |
31 | using std::endl; | |
32 | ||
33 | Double_t GetDecayTime(const Particle &parent, Double_t weakDecayLimit) { | |
03896fc4 | 34 | // |
35 | // return a random decay time according to the particle's width | |
36 | // | |
37 | ||
b1c2e580 | 38 | ParticlePDG *pDef = parent.Def(); |
39 | Double_t width = pDef->GetWidth(); //GeV | |
7b7936e9 | 40 | |
41 | // if particle is set to be stable then return 0 | |
42 | if(pDef->GetStableStatus()) | |
43 | return 0.; | |
44 | ||
45 | if(width > weakDecayLimit && weakDecayLimit>=0.0) { | |
b1c2e580 | 46 | const Double_t slope = parent.E() * 0.1973 / (pDef->GetMass() * width); |
b1c2e580 | 47 | return -slope * TMath::Log(gRandom->Rndm());//in fm/c |
48 | } | |
49 | ||
50 | return 0.; | |
51 | } | |
52 | ||
53 | ||
54 | extern "C" void mydelta_(); | |
55 | extern SERVICEEVCommon SERVICEEV; | |
56 | ||
03896fc4 | 57 | void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, const DatabasePDG* database) { |
58 | // | |
59 | // perform the decay | |
60 | // | |
61 | ||
b1c2e580 | 62 | // Get the PDG properties of the particle |
63 | ParticlePDG *pDef = parent.Def(); | |
64 | ||
65 | // Get the number of posible decay channels | |
66 | Int_t nDecayChannel = pDef->GetNDecayChannels(); | |
67 | ||
68 | // check the full branching of this specie | |
69 | Double_t fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays | |
70 | ||
71 | // return if particle has no branching | |
72 | if(fullBranching < 0.00001) | |
73 | return; | |
74 | ||
75 | // get the PDG mass of the specie | |
03896fc4 | 76 | Double_t pdgMass = pDef->GetMass(); |
77 | Int_t comprCodePyth=0; | |
78 | Float_t delta =0; | |
b1c2e580 | 79 | |
80 | Bool_t success = kFALSE; | |
81 | Int_t iterations = 0; | |
82 | // Try to decay the particle | |
83 | while(!success) { | |
7b7936e9 | 84 | // std::cout << "HadronDecayer::Decay() iteration #" << iterations << std::endl; |
85 | if(iterations>1000) { | |
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; | |
89 | return; | |
90 | } | |
3fa37a65 | 91 | |
b1c2e580 | 92 | // get a random mass using the Breit Wigner distribution |
03896fc4 | 93 | Double_t bwMass = gRandom->BreitWigner(pdgMass, pDef->GetWidth()); |
94 | // bwMass = pdgMass; | |
b1c2e580 | 95 | // Try to cut the Breit Wigner tail of the particle using the cuts from pythia |
03896fc4 | 96 | // The delta variable is obtained from pythia based on the specie |
b1c2e580 | 97 | int encoding =pDef->GetPDG(); |
98 | SERVICEEV.ipdg = encoding; | |
99 | mydelta_(); | |
03896fc4 | 100 | comprCodePyth=SERVICEEV.KC; |
101 | delta = SERVICEEV.delta;// PYDAT2.PMAS[KC][3]; | |
b1c2e580 | 102 | |
03896fc4 | 103 | //if there are no such particle in PYTHIA particle table, we take delta=0.4 |
104 | if(comprCodePyth==0){ | |
105 | bwMass=pdgMass; | |
106 | delta=0.0; | |
b1c2e580 | 107 | } |
108 | ||
109 | //bad delta - an exception | |
03896fc4 | 110 | if(comprCodePyth==254){ |
111 | bwMass=pdgMass; | |
112 | delta=0.0; | |
b1c2e580 | 113 | } |
114 | ||
7b7936e9 | 115 | // K0 decay into K0s or K0l |
116 | if(TMath::Abs(encoding)==311) { | |
03896fc4 | 117 | bwMass=pdgMass; |
118 | delta=0.0; | |
7b7936e9 | 119 | } |
120 | ||
b1c2e580 | 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 |
03896fc4 | 122 | if(comprCodePyth!=0 && delta>0 && (bwMass<pdgMass-delta || bwMass>pdgMass+delta)){ |
7b7936e9 | 123 | iterations++; |
b1c2e580 | 124 | continue; |
125 | } | |
b1c2e580 | 126 | |
127 | // check how many decay channels are allowed with the generated mass | |
03896fc4 | 128 | Int_t nAllowedChannels = database->GetNAllowedChannels(pDef, bwMass); |
b1c2e580 | 129 | // if no decay channels are posible with this mass, then generate another BW mass |
130 | if(nAllowedChannels==0) { | |
131 | iterations++; | |
3fa37a65 | 132 | std::cout << "HadronDecayer::Decay() no decays allowed at this BW mass" << std::endl; |
b1c2e580 | 133 | continue; |
134 | } | |
135 | ||
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; | |
141 | ||
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; | |
147 | while(!found) { | |
3fa37a65 | 148 | for(Int_t nChannel = 0; nChannel < nDecayChannel; nChannel++) { |
b1c2e580 | 149 | randValue -= pDef->GetDecayChannel(nChannel)->GetBranching(); |
03896fc4 | 150 | if(randValue <= 0. && database->IsChannelAllowed(pDef->GetDecayChannel(nChannel), bwMass)) { |
b1c2e580 | 151 | chosenChannel = nChannel; |
152 | found = kTRUE; | |
153 | break; | |
154 | } | |
155 | } | |
156 | channelIterations++; | |
157 | } | |
158 | ||
159 | // get the PDG information for the chosen decay channel | |
160 | DecayChannel *dc = pDef->GetDecayChannel(chosenChannel); | |
161 | Int_t nSec = dc->GetNDaughters(); | |
162 | ||
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()); | |
03896fc4 | 166 | Double_t bwEnergy = TMath::Sqrt(parent.Mom().X()*parent.Mom().X() + |
b1c2e580 | 167 | parent.Mom().Y()*parent.Mom().Y() + |
168 | parent.Mom().Z()*parent.Mom().Z() + | |
03896fc4 | 169 | bwMass*bwMass); |
b1c2e580 | 170 | |
03896fc4 | 171 | Int_t nb = (Int_t)parent.GetType(); //particle from jets |
b1c2e580 | 172 | |
03896fc4 | 173 | TLorentzVector MomparentBW(parent.Mom().X(), parent.Mom().Y(), parent.Mom().Z(), bwEnergy); |
b1c2e580 | 174 | parentBW.Mom(MomparentBW); |
175 | // take into account BW when calculating boost velocity (for wide resonances it matters) | |
176 | TVector3 velocityBW(parentBW.Mom().BoostVector()); | |
177 | ||
178 | // now we have an allowed decay | |
179 | // first case: one daughter particle | |
180 | if(nSec == 1) { | |
181 | // initialize the daughter particle | |
182 | Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0))); | |
7b7936e9 | 183 | p1.Pos(parentBW.Pos()); |
184 | p1.Mom(parent.Mom()); | |
b1c2e580 | 185 | p1.SetLastMotherPdg(parentBW.Encoding()); |
186 | p1.SetLastMotherDecayCoor(parentBW.Pos()); | |
187 | p1.SetLastMotherDecayMom(parentBW.Mom()); | |
03896fc4 | 188 | p1.SetType(nb); |
b1c2e580 | 189 | |
3fa37a65 | 190 | // add the daughter to the list of secondaries |
191 | Int_t parentIndex = parent.GetIndex(); | |
b1c2e580 | 192 | Int_t p1Index = p1.SetIndex(); // set the daughter index |
193 | p1.SetMother(parentIndex); // set the mother index for this daughter | |
3fa37a65 | 194 | parent.SetFirstDaughterIndex(p1Index); |
195 | parent.SetLastDaughterIndex(p1Index); | |
b1c2e580 | 196 | allocator.AddParticle(p1, output); |
197 | success = kTRUE; | |
198 | } | |
199 | // second case: two daughter particles | |
200 | else if(nSec == 2) { | |
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()); | |
206 | ||
207 | // calculate the momenta in rest frame of mother for the two particles (theta and phi are isotropic) | |
03896fc4 | 208 | MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), bwMass); |
b1c2e580 | 209 | |
210 | // boost to the laboratory system (to the mother velocity) | |
211 | p1.Mom().Boost(velocityBW); | |
212 | p2.Mom().Boost(velocityBW); | |
213 | ||
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()); | |
b1c2e580 | 221 | //set to daughters the same type as has mother |
03896fc4 | 222 | p1.SetType(nb); |
223 | p2.SetType(nb); | |
b1c2e580 | 224 | |
225 | ||
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 | |
232 | ||
b1c2e580 | 233 | if(deltaS>0.001) { |
b1c2e580 | 234 | iterations++; |
235 | continue; | |
236 | } | |
237 | // push particles to the list of secondaries | |
3fa37a65 | 238 | Int_t parentIndex = parent.GetIndex(); |
b1c2e580 | 239 | p1.SetIndex(); |
240 | p2.SetIndex(); | |
241 | p1.SetMother(parentIndex); | |
242 | p2.SetMother(parentIndex); | |
3fa37a65 | 243 | parent.SetFirstDaughterIndex(p1.GetIndex()); |
244 | parent.SetLastDaughterIndex(p2.GetIndex()); | |
b1c2e580 | 245 | allocator.AddParticle(p1, output); |
246 | allocator.AddParticle(p2, output); | |
247 | success = kTRUE; | |
3fa37a65 | 248 | |
b1c2e580 | 249 | } |
250 | ||
251 | // third case: three daughter particle | |
252 | else if(nSec == 3) { | |
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(); | |
03896fc4 | 264 | Double_t deltaMass = bwMass - mass1 - mass2 - mass3; |
b1c2e580 | 265 | |
266 | do { | |
267 | Double_t rd1 = gRandom->Rndm(); | |
268 | Double_t rd2 = gRandom->Rndm(); | |
269 | if (rd2 > rd1) | |
270 | std::swap(rd1, rd2); | |
271 | // 1 | |
272 | Double_t e = rd2*deltaMass; | |
273 | pAbs1 = TMath::Sqrt(e*e + 2*e*mass1); | |
274 | sumPabs = pAbs1; | |
275 | maxPabs = sumPabs; | |
276 | // 2 | |
277 | e = (1-rd1)*deltaMass; | |
278 | pAbs2 = TMath::Sqrt(e*e + 2*e*mass2); | |
279 | ||
280 | if(pAbs2 > maxPabs) | |
281 | maxPabs = pAbs2; | |
282 | ||
283 | sumPabs += pAbs2; | |
284 | // 3 | |
285 | e = (rd1-rd2)*deltaMass; | |
286 | pAbs3 = TMath::Sqrt(e*e + 2*e*mass3); | |
287 | ||
288 | if (pAbs3 > maxPabs) | |
289 | maxPabs = pAbs3; | |
290 | sumPabs += pAbs3; | |
291 | } while(maxPabs > sumPabs - maxPabs); | |
292 | ||
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); | |
299 | ||
300 | mom1.SetPxPyPzE(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta, 0); | |
301 | mom1 *= pAbs1; | |
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); | |
308 | ||
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, | |
312 | 0.); | |
313 | ||
314 | mom3 *= pAbs3*mom3.P(); | |
315 | mom2 = mom1; | |
316 | mom2 += mom3; | |
317 | mom2 *= -1.; | |
318 | // calculate energy | |
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)); | |
322 | ||
323 | // boost to Lab system | |
324 | mom1.Boost(velocityBW); | |
325 | mom2.Boost(velocityBW); | |
326 | mom3.Boost(velocityBW); | |
327 | ||
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()); | |
337 | ||
338 | //set to daughters the same type as has mother | |
03896fc4 | 339 | p1.SetType(nb); |
340 | p2.SetType(nb); | |
341 | p3.SetType(nb); | |
b1c2e580 | 342 | |
b1c2e580 | 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 | |
349 | if(deltaS>0.001) { | |
7b7936e9 | 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; | |
b1c2e580 | 357 | |
358 | iterations++; | |
359 | continue; | |
360 | } | |
361 | ||
3fa37a65 | 362 | Int_t parentIndex = parent.GetIndex(); |
b1c2e580 | 363 | p1.SetIndex(); |
364 | p2.SetIndex(); | |
365 | p3.SetIndex(); | |
366 | p1.SetMother(parentIndex); | |
367 | p2.SetMother(parentIndex); | |
368 | p3.SetMother(parentIndex); | |
3fa37a65 | 369 | parent.SetFirstDaughterIndex(p1.GetIndex()); |
370 | parent.SetLastDaughterIndex(p3.GetIndex()); | |
b1c2e580 | 371 | allocator.AddParticle(p1, output); |
372 | allocator.AddParticle(p2, output); | |
373 | allocator.AddParticle(p3, output); | |
374 | success = kTRUE; | |
375 | } | |
376 | } | |
377 | return; | |
378 | } |