Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TPythia6Decayer.cxx
Go to the documentation of this file.
1 // @(#)root/pythia6:$Id$
2 // Author: Christian Holm Christensen 22/04/06
3 // Much of this code has been lifted from AliROOT.
4 
5 /*************************************************************************
6  * Copyright (C) 1995-2006, Rene Brun and Fons Rademakers. *
7  * All rights reserved. *
8  * *
9  * For the licensing terms see $ROOTSYS/LICENSE. *
10  * For the list of contributors see $ROOTSYS/README/CREDITS. *
11  *************************************************************************/
12 
13 /** \class TPythia6Decayer
14  \ingroup pythia6
15 
16 This class implements the TVirtualMCDecayer interface.
17 
18 The TPythia6 singleton object is used to decay particles. Note, that since
19 this class modifies common blocks (global variables) it is defined as a
20 singleton.
21 */
22 
23 #include "TPythia6.h"
24 #include "TPythia6Decayer.h"
25 #include "TPDGCode.h"
26 #include "TLorentzVector.h"
27 #include "TClonesArray.h"
28 
29 ClassImp(TPythia6Decayer);
30 
31 TPythia6Decayer* TPythia6Decayer::fgInstance = 0;
32 
33 ////////////////////////////////////////////////////////////////////////////////
34 /// Get the singleton object.
35 
36 TPythia6Decayer* TPythia6Decayer::Instance()
37 {
38  if (!fgInstance) fgInstance = new TPythia6Decayer;
39  return fgInstance;
40 }
41 
42 ////////////////////////////////////////////////////////////////////////////////
43 /// Constructor
44 
45 TPythia6Decayer::TPythia6Decayer()
46  : fDecay(kMaxDecay),
47  fBraPart(501)
48 {
49  fBraPart.Reset(1);
50 }
51 
52 ////////////////////////////////////////////////////////////////////////////////
53 /// Initialize the decayer
54 
55 void TPythia6Decayer::Init()
56 {
57  static Bool_t init = kFALSE;
58  if (init) return;
59  init = kTRUE;
60  ForceDecay();
61 }
62 
63 ////////////////////////////////////////////////////////////////////////////////
64 /// Decay a particle of type IDPART (PDG code) and momentum P.
65 
66 void TPythia6Decayer::Decay(Int_t idpart, TLorentzVector* p)
67 {
68  if (!p) return;
69  TPythia6::Instance()->Py1ent(0, idpart, p->Energy(), p->Theta(), p->Phi());
70  TPythia6::Instance()->GetPrimaries();
71 }
72 
73 ////////////////////////////////////////////////////////////////////////////////
74 /// Get the decay products into the passed PARTICLES TClonesArray of
75 /// TParticles
76 
77 Int_t TPythia6Decayer::ImportParticles(TClonesArray *particles)
78 {
79  return TPythia6::Instance()->ImportParticles(particles,"All");
80 }
81 
82 ////////////////////////////////////////////////////////////////////////////////
83 /// Force a particular decay type
84 
85 void TPythia6Decayer::SetForceDecay(Int_t type)
86 {
87  if (type > kMaxDecay) {
88  Warning("SetForceDecay", "Invalid decay mode: %d", type);
89  return;
90  }
91  fDecay = EDecayType(type);
92 }
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 /// Force a particle decay mode
96 
97 void TPythia6Decayer::ForceDecay()
98 {
99  EDecayType decay=fDecay;
100  TPythia6::Instance()->SetMSTJ(21,2);
101  if (decay == kNoDecayHeavy) return;
102 
103  //
104  // select mode
105  Int_t products[3];
106  Int_t mult[3];
107 
108  switch (decay) {
109  case kHardMuons:
110  products[0] = 13;
111  products[1] = 443;
112  products[2] = 100443;
113  mult[0] = 1;
114  mult[1] = 1;
115  mult[2] = 1;
116  ForceParticleDecay( 511, products, mult, 3);
117  ForceParticleDecay( 521, products, mult, 3);
118  ForceParticleDecay( 531, products, mult, 3);
119  ForceParticleDecay( 5122, products, mult, 3);
120  ForceParticleDecay( 5132, products, mult, 3);
121  ForceParticleDecay( 5232, products, mult, 3);
122  ForceParticleDecay( 5332, products, mult, 3);
123  ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
124  ForceParticleDecay( 443, 13, 2); // J/Psi -> mu+ mu-
125 
126  ForceParticleDecay( 411,13,1); // D+/-
127  ForceParticleDecay( 421,13,1); // D0
128  ForceParticleDecay( 431,13,1); // D_s
129  ForceParticleDecay( 4122,13,1); // Lambda_c
130  ForceParticleDecay( 4132,13,1); // Xsi_c
131  ForceParticleDecay( 4232,13,1); // Sigma_c
132  ForceParticleDecay( 4332,13,1); // Omega_c
133  break;
134  case kSemiMuonic:
135  ForceParticleDecay( 411,13,1); // D+/-
136  ForceParticleDecay( 421,13,1); // D0
137  ForceParticleDecay( 431,13,1); // D_s
138  ForceParticleDecay( 4122,13,1); // Lambda_c
139  ForceParticleDecay( 4132,13,1); // Xsi_c
140  ForceParticleDecay( 4232,13,1); // Sigma_c
141  ForceParticleDecay( 4332,13,1); // Omega_c
142  ForceParticleDecay( 511,13,1); // B0
143  ForceParticleDecay( 521,13,1); // B+/-
144  ForceParticleDecay( 531,13,1); // B_s
145  ForceParticleDecay( 5122,13,1); // Lambda_b
146  ForceParticleDecay( 5132,13,1); // Xsi_b
147  ForceParticleDecay( 5232,13,1); // Sigma_b
148  ForceParticleDecay( 5332,13,1); // Omega_b
149  break;
150  case kDiMuon:
151  ForceParticleDecay( 113,13,2); // rho
152  ForceParticleDecay( 221,13,2); // eta
153  ForceParticleDecay( 223,13,2); // omega
154  ForceParticleDecay( 333,13,2); // phi
155  ForceParticleDecay( 443,13,2); // J/Psi
156  ForceParticleDecay(100443,13,2);// Psi'
157  ForceParticleDecay( 553,13,2); // Upsilon
158  ForceParticleDecay(100553,13,2);// Upsilon'
159  ForceParticleDecay(200553,13,2);// Upsilon''
160  break;
161  case kSemiElectronic:
162  ForceParticleDecay( 411,11,1); // D+/-
163  ForceParticleDecay( 421,11,1); // D0
164  ForceParticleDecay( 431,11,1); // D_s
165  ForceParticleDecay( 4122,11,1); // Lambda_c
166  ForceParticleDecay( 4132,11,1); // Xsi_c
167  ForceParticleDecay( 4232,11,1); // Sigma_c
168  ForceParticleDecay( 4332,11,1); // Omega_c
169  ForceParticleDecay( 511,11,1); // B0
170  ForceParticleDecay( 521,11,1); // B+/-
171  ForceParticleDecay( 531,11,1); // B_s
172  ForceParticleDecay( 5122,11,1); // Lambda_b
173  ForceParticleDecay( 5132,11,1); // Xsi_b
174  ForceParticleDecay( 5232,11,1); // Sigma_b
175  ForceParticleDecay( 5332,11,1); // Omega_b
176  break;
177  case kDiElectron:
178  ForceParticleDecay( 113,11,2); // rho
179  ForceParticleDecay( 333,11,2); // phi
180  ForceParticleDecay( 221,11,2); // eta
181  ForceParticleDecay( 223,11,2); // omega
182  ForceParticleDecay( 443,11,2); // J/Psi
183  ForceParticleDecay(100443,11,2);// Psi'
184  ForceParticleDecay( 553,11,2); // Upsilon
185  ForceParticleDecay(100553,11,2);// Upsilon'
186  ForceParticleDecay(200553,11,2);// Upsilon''
187  break;
188  case kBJpsiDiMuon:
189 
190  products[0] = 443;
191  products[1] = 100443;
192  mult[0] = 1;
193  mult[1] = 1;
194 
195  ForceParticleDecay( 511, products, mult, 2); // B0 -> J/Psi (Psi') X
196  ForceParticleDecay( 521, products, mult, 2); // B+/- -> J/Psi (Psi') X
197  ForceParticleDecay( 531, products, mult, 2); // B_s -> J/Psi (Psi') X
198  ForceParticleDecay( 5122, products, mult, 2); // Lambda_b -> J/Psi (Psi') X
199  ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
200  ForceParticleDecay( 443,13,2); // J/Psi -> mu+ mu-
201  break;
202  case kBPsiPrimeDiMuon:
203  ForceParticleDecay( 511,100443,1); // B0
204  ForceParticleDecay( 521,100443,1); // B+/-
205  ForceParticleDecay( 531,100443,1); // B_s
206  ForceParticleDecay( 5122,100443,1); // Lambda_b
207  ForceParticleDecay(100443,13,2); // Psi'
208  break;
209  case kBJpsiDiElectron:
210  ForceParticleDecay( 511,443,1); // B0
211  ForceParticleDecay( 521,443,1); // B+/-
212  ForceParticleDecay( 531,443,1); // B_s
213  ForceParticleDecay( 5122,443,1); // Lambda_b
214  ForceParticleDecay( 443,11,2); // J/Psi
215  break;
216  case kBJpsi:
217  ForceParticleDecay( 511,443,1); // B0
218  ForceParticleDecay( 521,443,1); // B+/-
219  ForceParticleDecay( 531,443,1); // B_s
220  ForceParticleDecay( 5122,443,1); // Lambda_b
221  break;
222  case kBPsiPrimeDiElectron:
223  ForceParticleDecay( 511,100443,1); // B0
224  ForceParticleDecay( 521,100443,1); // B+/-
225  ForceParticleDecay( 531,100443,1); // B_s
226  ForceParticleDecay( 5122,100443,1); // Lambda_b
227  ForceParticleDecay(100443,11,2); // Psi'
228  break;
229  case kPiToMu:
230  ForceParticleDecay(211,13,1); // pi->mu
231  break;
232  case kKaToMu:
233  ForceParticleDecay(321,13,1); // K->mu
234  break;
235  case kWToMuon:
236  ForceParticleDecay( 24, 13,1); // W -> mu
237  break;
238  case kWToCharm:
239  ForceParticleDecay( 24, 4,1); // W -> c
240  break;
241  case kWToCharmToMuon:
242  ForceParticleDecay( 24, 4,1); // W -> c
243  ForceParticleDecay( 411,13,1); // D+/- -> mu
244  ForceParticleDecay( 421,13,1); // D0 -> mu
245  ForceParticleDecay( 431,13,1); // D_s -> mu
246  ForceParticleDecay( 4122,13,1); // Lambda_c
247  ForceParticleDecay( 4132,13,1); // Xsi_c
248  ForceParticleDecay( 4232,13,1); // Sigma_c
249  ForceParticleDecay( 4332,13,1); // Omega_c
250  break;
251  case kZDiMuon:
252  ForceParticleDecay( 23, 13,2); // Z -> mu+ mu-
253  break;
254  case kHadronicD:
255  ForceHadronicD();
256  break;
257  case kPhiKK:
258  ForceParticleDecay(333,321,2); // Phi->K+K-
259  break;
260  case kOmega:
261  ForceOmega();
262  case kAll:
263  break;
264  case kNoDecay:
265  TPythia6::Instance()->SetMSTJ(21,0);
266  break;
267  case kNoDecayHeavy: break; // cannot get here: early return above
268  case kMaxDecay: break;
269  }
270 }
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 /// Get the partial branching ratio for a particle of type IPART (a
274 /// PDG code).
275 
276 Float_t TPythia6Decayer::GetPartialBranchingRatio(Int_t ipart)
277 {
278  Int_t kc = TPythia6::Instance()->Pycomp(TMath::Abs(ipart));
279  // return TPythia6::Instance()->GetBRAT(kc);
280  return fBraPart[kc];
281 }
282 
283 ////////////////////////////////////////////////////////////////////////////////
284 /// Get the life-time of a particle of type KF (a PDG code).
285 
286 Float_t TPythia6Decayer::GetLifetime(Int_t kf)
287 {
288  Int_t kc=TPythia6::Instance()->Pycomp(TMath::Abs(kf));
289  return TPythia6::Instance()->GetPMAS(kc,4) * 3.3333e-12;
290 }
291 
292 ////////////////////////////////////////////////////////////////////////////////
293 /// Read in particle data from an ASCII file. The file name must
294 /// previously have been set using the member function
295 /// SetDecayTableFile.
296 
297 void TPythia6Decayer::ReadDecayTable()
298 {
299  if (fDecayTableFile.IsNull()) {
300  Warning("ReadDecayTable", "No file set");
301  return;
302  }
303  Int_t lun = 15;
304  TPythia6::Instance()->OpenFortranFile(lun,
305  const_cast<char*>(fDecayTableFile.Data()));
306  TPythia6::Instance()->Pyupda(3,lun);
307  TPythia6::Instance()->CloseFortranFile(lun);
308 }
309 
310 // ===================================================================
311 // BEGIN COMMENT
312 //
313 // It would be better if the particle and decay information could be
314 // read from the current TDatabasePDG instance.
315 //
316 // However, it seems to me that some information is missing. In
317 // particular
318 //
319 // - The broadning cut-off,
320 // - Resonance width
321 // - Color charge
322 // - MWID (?)
323 //
324 // Further more, it's not clear to me at least, what all the
325 // parameters Pythia needs are.
326 //
327 // Code like the below could be used to make a temporary file that
328 // Pythia could then read in. Ofcourse, one could also manipulate
329 // the data structures directly, but that's propably more dangerous.
330 //
331 #if 0
332 void PrintPDG(TParticlePDG* pdg)
333 {
334  TParticlePDG* anti = pdg->AntiParticle();
335  const char* antiName = (anti ? anti->GetName() : "");
336  Int_t color = 0;
337  switch (TMath::Abs(pdg->PdgCode())) {
338  case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: // Quarks
339  color = 1; break;
340  case 21: // Gluon
341  color = 2; break;
342  case 1103:
343  case 2101: case 2103: case 2203:
344  case 3101: case 3103: case 3201: case 3203: case 3303:
345  case 4101: case 4103: case 4201: case 4203: case 4301: case 4303: case 4403:
346  case 5101: case 5103: case 5201: case 5203: case 5301: case 5303: case 5401:
347  case 5403: case 5503:
348  // Quark combinations
349  color = -1; break;
350  case 1000001: case 1000002: case 1000003: case 1000004: case 1000005:
351  case 1000006: // super symmetric partners to quars
352  color = 1; break;
353  case 1000021: // ~g
354  color = 2; break;
355  case 2000001: case 2000002: case 2000003: case 2000004: case 2000005:
356  case 2000006: // R hadrons
357  color = 1; break;
358  case 3000331: case 3100021: case 3200111: case 3100113: case 3200113:
359  case 3300113: case 3400113:
360  // Technicolor
361  color = 2; break;
362  case 4000001: case 4000002:
363  color = 1; break;
364  case 9900443: case 9900441: case 9910441: case 9900553: case 9900551:
365  case 9910551:
366  color = 2; break;
367  }
368  std::cout << std::right
369  << " " << std::setw(9) << pdg->PdgCode()
370  << " " << std::left << std::setw(16) << pdg->GetName()
371  << " " << std::setw(16) << antiName
372  << std::right
373  << std::setw(3) << Int_t(pdg->Charge())
374  << std::setw(3) << color
375  << std::setw(3) << (anti ? 1 : 0)
376  << std::fixed << std::setprecision(5)
377  << std::setw(12) << pdg->Mass()
378  << std::setw(12) << pdg->Width()
379  << std::setw(12) << 0 // Broad
380  << std::scientific
381  << " " << std::setw(13) << pdg->Lifetime()
382  << std::setw(3) << 0 // MWID
383  << std::setw(3) << pdg->Stable()
384  << std::endl;
385 }
386 
387 void MakeDecayList()
388 {
389  TDatabasePDG* pdgDB = TDatabasePDG::Instance();
390  pdgDB->ReadPDGTable();
391  const THashList* pdgs = pdgDB->ParticleList();
392  TParticlePDG* pdg = 0;
393  TIter nextPDG(pdgs);
394  while ((pdg = static_cast<TParticlePDG*>(nextPDG()))) {
395  // std::cout << "Processing " << pdg->GetName() << std::endl;
396  PrintPDG(pdg);
397 
398  TObjArray* decays = pdg->DecayList();
399  TDecayChannel* decay = 0;
400  TIter nextDecay(decays);
401  while ((decay = static_cast<TDecayChannel*>(nextDecay()))) {
402  // std::cout << "Processing decay number " << decay->Number() << std::endl;
403  }
404  }
405 }
406 #endif
407 // END COMMENT
408 // ===================================================================
409 
410 ////////////////////////////////////////////////////////////////////////////////
411 /// write particle data to an ASCII file. The file name must
412 /// previously have been set using the member function
413 /// SetDecayTableFile.
414 ///
415 /// Users can use this function to make an initial decay list file,
416 /// which then can be edited by hand, and re-loaded into the decayer
417 /// using ReadDecayTable.
418 ///
419 /// The file syntax is
420 ///
421 /// particle_list : partcle_data
422 /// | particle_list particle_data
423 /// ;
424 /// particle_data : particle_info
425 /// | particle_info '\n' decay_list
426 /// ;
427 /// particle_info : See below
428 /// ;
429 /// decay_list : decay_entry
430 /// | decay_list decay_entry
431 /// ;
432 /// decay_entry : See below
433 ///
434 /// The particle_info consists of 13 fields:
435 ///
436 /// PDG code int
437 /// Name string
438 /// Anti-particle name string if there's no anti-particle,
439 /// then this field must be the
440 /// empty string
441 /// Electic charge int in units of |e|/3
442 /// Color charge int in units of quark color charges
443 /// Have anti-particle int 1 of there's an anti-particle
444 /// to this particle, or 0
445 /// otherwise
446 /// Mass float in units of GeV
447 /// Resonance width float
448 /// Max broadning float
449 /// Lifetime float
450 /// MWID int ??? (some sort of flag)
451 /// Decay int 1 if it decays. 0 otherwise
452 ///
453 /// The format to write these entries in are
454 ///
455 /// " %9 %-16s %-16s%3d%3d%3d%12.5f%12.5f%12.5f%13.gf%3d%d\n"
456 ///
457 /// The decay_entry consists of 8 fields:
458 ///
459 /// On/Off int 1 for on, -1 for off
460 /// Matrix element type int
461 /// Branching ratio float
462 /// Product 1 int PDG code of decay product 1
463 /// Product 2 int PDG code of decay product 2
464 /// Product 3 int PDG code of decay product 3
465 /// Product 4 int PDG code of decay product 4
466 /// Product 5 int PDG code of decay product 5
467 ///
468 /// The format for these lines are
469 ///
470 /// " %5d%5d%12.5f%10d%10d%10d%10d%10d\n"
471 ///
472 
473 void TPythia6Decayer::WriteDecayTable()
474 {
475  if (fDecayTableFile.IsNull()) {
476  Warning("ReadDecayTable", "No file set");
477  return;
478  }
479  Int_t lun = 15;
480  TPythia6::Instance()->OpenFortranFile(lun,
481  const_cast<char*>(fDecayTableFile.Data()));
482  TPythia6::Instance()->Pyupda(1,lun);
483  TPythia6::Instance()->CloseFortranFile(lun);
484 }
485 
486 ////////////////////////////////////////////////////////////////////////////////
487 /// Count number of decay products
488 
489 Int_t TPythia6Decayer::CountProducts(Int_t channel, Int_t particle)
490 {
491  Int_t np = 0;
492  for (Int_t i = 1; i <= 5; i++)
493  if (TMath::Abs(TPythia6::Instance()->GetKFDP(channel,i)) == particle) np++;
494  return np;
495 }
496 
497 ////////////////////////////////////////////////////////////////////////////////
498 /// Force golden D decay modes
499 
500 void TPythia6Decayer::ForceHadronicD()
501 {
502  const Int_t kNHadrons = 4;
503  Int_t channel;
504  Int_t hadron[kNHadrons] = {411, 421, 431, 4112};
505 
506  // for D+ -> K0* (-> K- pi+) pi+
507  Int_t iKstar0 = 313;
508  Int_t iKstarbar0 = -313;
509  Int_t products[2] = {kKPlus, kPiMinus}, mult[2] = {1, 1};
510  ForceParticleDecay(iKstar0, products, mult, 2);
511 
512  // for Ds -> Phi pi+
513  Int_t iPhi = 333;
514  ForceParticleDecay(iPhi,kKPlus,2); // Phi->K+K-
515 
516  Int_t decayP1[kNHadrons][3] = {
517  {kKMinus, kPiPlus, kPiPlus},
518  {kKMinus, kPiPlus, 0 },
519  {kKPlus , iKstarbar0, 0 },
520  {-1 , -1 , -1 }
521  };
522  Int_t decayP2[kNHadrons][3] = {
523  {iKstarbar0, kPiPlus, 0 },
524  {-1 , -1 , -1 },
525  {iPhi , kPiPlus, 0 },
526  {-1 , -1 , -1 }
527  };
528 
529  TPythia6* pyth = TPythia6::Instance();
530  for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++) {
531  Int_t kc = pyth->Pycomp(hadron[ihadron]);
532  pyth->SetMDCY(kc,1,1);
533  Int_t ifirst = pyth->GetMDCY(kc,2);
534  Int_t ilast = ifirst + pyth->GetMDCY(kc,3)-1;
535 
536  for (channel = ifirst; channel <= ilast; channel++) {
537  if ((pyth->GetKFDP(channel,1) == decayP1[ihadron][0] &&
538  pyth->GetKFDP(channel,2) == decayP1[ihadron][1] &&
539  pyth->GetKFDP(channel,3) == decayP1[ihadron][2] &&
540  pyth->GetKFDP(channel,4) == 0) ||
541  (pyth->GetKFDP(channel,1) == decayP2[ihadron][0] &&
542  pyth->GetKFDP(channel,2) == decayP2[ihadron][1] &&
543  pyth->GetKFDP(channel,3) == decayP2[ihadron][2] &&
544  pyth->GetKFDP(channel,4) == 0)) {
545  pyth->SetMDME(channel,1,1);
546  } else {
547  pyth->SetMDME(channel,1,0);
548  fBraPart[kc] -= pyth->GetBRAT(channel);
549  } // selected channel ?
550  } // decay channels
551  } // hadrons
552 }
553 
554 ////////////////////////////////////////////////////////////////////////////////
555 ///
556 /// Force decay of particle into products with multiplicity mult
557 
558 void TPythia6Decayer::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
559 {
560  TPythia6* pyth = TPythia6::Instance();
561 
562  Int_t kc = pyth->Pycomp(particle);
563  pyth->SetMDCY(kc,1,1);
564 
565  Int_t ifirst = pyth->GetMDCY(kc,2);
566  Int_t ilast = ifirst + pyth->GetMDCY(kc,3)-1;
567  fBraPart[kc] = 1;
568 
569  //
570  // Loop over decay channels
571  for (Int_t channel= ifirst; channel <= ilast; channel++) {
572  if (CountProducts(channel,product) >= mult) {
573  pyth->SetMDME(channel,1,1);
574  } else {
575  pyth->SetMDME(channel,1,0);
576  fBraPart[kc]-=pyth->GetBRAT(channel);
577  }
578  }
579 }
580 
581 ////////////////////////////////////////////////////////////////////////////////
582 ///
583 /// Force decay of particle into products with multiplicity mult
584 
585 void TPythia6Decayer::ForceParticleDecay(Int_t particle, Int_t* products,
586  Int_t* mult, Int_t npart)
587 {
588  TPythia6* pyth = TPythia6::Instance();
589 
590  Int_t kc = pyth->Pycomp(particle);
591  pyth->SetMDCY(kc,1,1);
592  Int_t ifirst = pyth->GetMDCY(kc,2);
593  Int_t ilast = ifirst+pyth->GetMDCY(kc,3)-1;
594  fBraPart[kc] = 1;
595  //
596  // Loop over decay channels
597  for (Int_t channel = ifirst; channel <= ilast; channel++) {
598  Int_t nprod = 0;
599  for (Int_t i = 0; i < npart; i++)
600  nprod += (CountProducts(channel, products[i]) >= mult[i]);
601  if (nprod)
602  pyth->SetMDME(channel,1,1);
603  else {
604  pyth->SetMDME(channel,1,0);
605  fBraPart[kc] -= pyth->GetBRAT(channel);
606  }
607  }
608 }
609 
610 ////////////////////////////////////////////////////////////////////////////////
611 /// Force Omega -> Lambda K- Decay
612 
613 void TPythia6Decayer::ForceOmega()
614 {
615  TPythia6* pyth = TPythia6::Instance();
616 
617  Int_t kc = pyth->Pycomp(3334);
618  pyth->SetMDCY(kc,1,1);
619  Int_t ifirst = pyth->GetMDCY(kc,2);
620  Int_t ilast = ifirst + pyth->GetMDCY(kc,3)-1;
621  for (Int_t channel = ifirst; channel <= ilast; channel++) {
622  if (pyth->GetKFDP(channel,1) == kLambda0 &&
623  pyth->GetKFDP(channel,2) == kKMinus &&
624  pyth->GetKFDP(channel,3) == 0)
625  pyth->SetMDME(channel,1,1);
626  else
627  pyth->SetMDME(channel,1,0);
628  // selected channel ?
629  } // decay channels
630 }