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