UserDecay.C 5.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181
  1. void SetMode1_for_PDG(Int_t pdg, Float_t b0,Int_t m00,Int_t m01)
  2. {
  3. Int_t mode[6][3];
  4. Float_t bratio[6];
  5. for (Int_t kz = 0; kz < 6; kz++) {
  6. bratio[kz] = 0.;
  7. mode[kz][0] = 0;
  8. mode[kz][1] = 0;
  9. mode[kz][2] = 0;
  10. }
  11. bratio[0] = b0; mode[0][0] = m00; mode[0][1] = m01;
  12. gMC->SetDecayMode(pdg,bratio,mode);
  13. }
  14. void SetMode2_for_PDG(Int_t pdg, Float_t b0,Int_t m00,Int_t m01,Int_t m02,Float_t b1,Int_t m10,Int_t m11,Int_t m12)
  15. {
  16. Int_t mode[6][3];
  17. Float_t bratio[6];
  18. for (Int_t kz = 0; kz < 6; kz++) {
  19. bratio[kz] = 0.;
  20. mode[kz][0] = 0;
  21. mode[kz][1] = 0;
  22. mode[kz][2] = 0;
  23. }
  24. bratio[0] = b0; mode[0][0] = m00; mode[0][1] = m01; mode[0][2] = m02;
  25. bratio[0] = b1; mode[1][0] = m10; mode[1][1] = m11; mode[1][2] = m12;
  26. gMC->SetDecayMode(pdg,bratio,mode);
  27. }
  28. void SetMode_for_PDG(Int_t pdg, Float_t* b, Int_t* m0, Int_t* m1, Int_t *m2, Int_t N)
  29. {
  30. Int_t mode[6][3];
  31. Float_t bratio[6];
  32. for (Int_t kz = 0; kz < 6; kz++)
  33. {
  34. bratio[kz] = 0.;
  35. mode[kz][0] = 0;
  36. mode[kz][1] = 0;
  37. mode[kz][2] = 0;
  38. }
  39. if (N < 7)
  40. {
  41. for (Int_t nn = 0; nn < N; nn++)
  42. {
  43. bratio[nn] = b[nn];
  44. mode[nn][0] = m0[nn];
  45. mode[nn][1] = m1[nn];
  46. mode[nn][2] = m2[nn];
  47. }
  48. for (Int_t nn = 0; nn < 6; nn++)
  49. gMC->SetDecayMode(pdg,bratio,mode);
  50. }
  51. }
  52. void UserDecayConfig()
  53. {
  54. const Double_t kAu2Gev=0.9314943228;
  55. const Double_t khSlash = 1.0545726663e-27;
  56. const Double_t kErg2Gev = 1/1.6021773349e-3;
  57. const Double_t khShGev = khSlash*kErg2Gev;
  58. const Double_t kYear2Sec = 3600*24*365.25;
  59. cout << "Loading User Decay Config from macro"<< endl;
  60. TDatabasePDG *db= TDatabasePDG::Instance();
  61. TParticlePDG *p=0;
  62. Int_t m1[6], m2[6], m3[6], ip;
  63. Float_t b[6], totB;
  64. Int_t AlphaPDG, He5PDG, He3, H3L, H4L, prot, deut, neut, H3, He4L;
  65. p = db->GetParticle("proton");
  66. if (p) prot = p->PdgCode();
  67. p = db->GetParticle("Deuteron");
  68. if (p) deut = p->PdgCode();
  69. else { deut = 1000010020; db->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,0,3,"Ion",deut); }
  70. p = db->GetParticle("neutron");
  71. if (p) neut = p->PdgCode();
  72. p = db->GetParticle("Triton");
  73. if (p) H3 = p->PdgCode();
  74. else { H3 = 1000010030; db->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,khShGev/(12.33*kYear2Sec),3,"Ion",H3); }
  75. p = db->GetParticle("He4L");
  76. if (p) He4L = p->PdgCode();
  77. else { He4L = 1010020040; db->AddParticle("He4L","He4L",3.92501,kFALSE,khShGev/(12.33*kYear2Sec),6,"Ion",He4L);
  78. gMC->DefineParticle(He4L, "He4L", kPTHadron, 3.92501 , 2.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);}
  79. p = db->GetParticle("Alpha");
  80. if (!p) {
  81. p = db->GetParticle("alpha");
  82. if (!p) {
  83. p = db->GetParticle("He4");
  84. }
  85. }
  86. if(p) AlphaPDG = p->PdgCode();
  87. else { AlphaPDG = 1000020040; db->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,khShGev/(12.33*kYear2Sec),6,"Ion",AlphaPDG); }
  88. cout << "!!!! 0-0 !!!!!!! " << (db->ParticleList()->GetEntries()) << " " << AlphaPDG << endl;
  89. p = db->GetParticle("He3");
  90. if (!p) p = db->GetParticle("HE3");
  91. if (p) He3 = p->PdgCode();
  92. else { He3=1000020030; db->AddParticle("HE3","HE3", 2.80923,kFALSE,0,6,"Ion",He3); }
  93. cout << "!!!! 0-1 !!!!!!! " << (db->ParticleList()->GetEntries()) << " " << He3 << endl;
  94. if (He3)
  95. {
  96. p = db->GetParticle("H3L");
  97. if (p) H3L = p->PdgCode();
  98. else { H3L = 1010010030; gMC->DefineParticle(H3L, "H3L", kPTHadron, 2.99131 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE); }
  99. // ----- H3L decays -----
  100. //SetMode1_for_PDG(H3L,100,-211,He3);
  101. for (Int_t kz=0; kz<6; kz++) { b[kz] = 0.; m1[kz] = 0; m2[kz] = 0; m3[kz] = 0; }
  102. ip = -1; totB = 0.;
  103. b[++ip] = 24.7; m1[ip] = -211; m2[ip] = He3; m3[ip] = 0; // pi-, He3
  104. b[++ip] = 12.35; m1[ip] = 111; m2[ip] = H3; m3[ip] = 0; // pi0, H3
  105. b[++ip] = 36.7; m1[ip] = -211; m2[ip] = prot; m3[ip] = deut; // pi-, p, d
  106. b[++ip] = 18.35; m1[ip] = 111; m2[ip] = neut; m3[ip] = deut; // pi0, n, d
  107. b[++ip] = 0.2; m1[ip] = neut; m2[ip] = deut; m3[ip] = 0; // n, d
  108. b[++ip] = 1.5; m1[ip] = neut; m2[ip] = neut; m3[ip] = prot; // n, n, p
  109. for (Int_t kz=0; kz<ip+1; kz++) totB += b[kz] / 100.;
  110. for (Int_t kz=0; kz<ip+1; kz++) b[kz] /= totB;
  111. SetMode_for_PDG(H3L,b,m1,m2,m3,++ip);
  112. p = db->GetParticle("H4L");
  113. if (p) H4L = p->PdgCode();
  114. else { H4L = 1010010040; gMC->DefineParticle(H4L, "H4L", kPTHadron, 3.92503 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); }
  115. // ----- H4L decays -----
  116. //SetMode1_for_PDG(H4L,100,-211,AlphaPDG);
  117. for (Int_t kz=0; kz<6; kz++) { b[kz] = 0.; m1[kz] = 0; m2[kz] = 0; m3[kz] = 0; }
  118. ip = -1; totB = 0.;
  119. b[++ip] = 75.; m1[ip] = -211; m2[ip] = AlphaPDG; m3[ip] = 0; // pi-, He4
  120. b[++ip] = 25.; m1[ip] = H3; m2[ip] = prot; m3[ip] = -211; // H3, prot, pi-
  121. for (Int_t kz=0; kz<ip+1; kz++) totB += b[kz] / 100.;
  122. for (Int_t kz=0; kz<ip+1; kz++) b[kz] /= totB;
  123. SetMode_for_PDG(H4L,b,m1,m2,m3,++ip);
  124. }
  125. // ----- He4L decays -----
  126. for (Int_t kz=0; kz<6; kz++) { b[kz] = 0.; m1[kz] = 0; m2[kz] = 0; m3[kz] = 0; }
  127. ip = -1; totB = 0.;
  128. b[++ip] = 32.; m1[ip] = -211; m2[ip] = He3; m3[ip] = prot; // pi-, He3, p
  129. b[++ip] = 35.; m1[ip] = 111; m2[ip] = AlphaPDG; m3[ip] = 0; // pi0, He4
  130. b[++ip] = 14.; m1[ip] = 111; m2[ip] = He3; m3[ip] = neut; // pi0, He3, neut
  131. b[++ip] = 19.; m1[ip] = prot; m2[ip] = H3; m3[ip] = 0; // p, H3
  132. for (Int_t kz=0; kz<ip+1; kz++) totB += b[kz] / 100.;
  133. for (Int_t kz=0; kz<ip+1; kz++) b[kz] /= totB;
  134. SetMode_for_PDG(He4L,b,m1,m2,m3,++ip);
  135. p = db->GetParticle("LN");
  136. if (!p) {
  137. p = db->GetParticle("LambdaNeutron");
  138. if (!p) {
  139. gMC->DefineParticle(1010000020, "LN", kPTNeutron, 2.054 , 0.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
  140. p = db->GetParticle("LN");
  141. }
  142. }
  143. if (p) SetMode1_for_PDG(p->PdgCode(),100,1000010020,-211);
  144. cout << "!!!! 1-0 !!!!!!! " << (db->ParticleList()->GetEntries()) << " " << H3L<< " " << H4L << endl;
  145. }