BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
UserDecay.C
Go to the documentation of this file.
1void 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
6 for (Int_t kz = 0; kz < 6; kz++) {
7 bratio[kz] = 0.;
8 mode[kz][0] = 0;
9 mode[kz][1] = 0;
10 mode[kz][2] = 0;
11 }
12
13 bratio[0] = b0;
14 mode[0][0] = m00;
15 mode[0][1] = m01;
16
17 gMC->SetDecayMode(pdg, bratio, mode);
18}
19
20void SetMode2_for_PDG(Int_t pdg,
21 Float_t b0,
22 Int_t m00,
23 Int_t m01,
24 Int_t m02,
25 Float_t b1,
26 Int_t m10,
27 Int_t m11,
28 Int_t m12)
29{
30 Int_t mode[6][3];
31 Float_t bratio[6];
32
33 for (Int_t kz = 0; kz < 6; kz++) {
34 bratio[kz] = 0.;
35 mode[kz][0] = 0;
36 mode[kz][1] = 0;
37 mode[kz][2] = 0;
38 }
39
40 bratio[0] = b0;
41 mode[0][0] = m00;
42 mode[0][1] = m01;
43 mode[0][2] = m02;
44 bratio[0] = b1;
45 mode[1][0] = m10;
46 mode[1][1] = m11;
47 mode[1][2] = m12;
48
49 gMC->SetDecayMode(pdg, bratio, mode);
50}
51
52void SetMode_for_PDG(Int_t pdg, Float_t* b, Int_t* m0, Int_t* m1, Int_t* m2, Int_t N)
53{
54 Int_t mode[6][3];
55 Float_t bratio[6];
56
57 for (Int_t kz = 0; kz < 6; kz++) {
58 bratio[kz] = 0.;
59 mode[kz][0] = 0;
60 mode[kz][1] = 0;
61 mode[kz][2] = 0;
62 }
63 if (N < 7) {
64 for (Int_t nn = 0; nn < N; nn++) {
65 bratio[nn] = b[nn];
66 mode[nn][0] = m0[nn];
67 mode[nn][1] = m1[nn];
68 mode[nn][2] = m2[nn];
69 }
70 for (Int_t nn = 0; nn < 6; nn++)
71 gMC->SetDecayMode(pdg, bratio, mode);
72 }
73}
74
75void UserDecayConfig()
76{
77 const Double_t kAu2Gev = 0.9314943228;
78 const Double_t khSlash = 1.0545726663e-27;
79 const Double_t kErg2Gev = 1 / 1.6021773349e-3;
80 const Double_t khShGev = khSlash * kErg2Gev;
81 const Double_t kYear2Sec = 3600 * 24 * 365.25;
82 cout << "Loading User Decay Config from macro" << endl;
83 TDatabasePDG* db = TDatabasePDG::Instance();
84 TParticlePDG* p = 0;
85
86 Int_t m1[6], m2[6], m3[6], ip;
87 Float_t b[6], totB;
88 Int_t AlphaPDG, He5PDG, He3, H3L, H4L, prot, deut, neut, H3, He4L;
89
90 p = db->GetParticle("proton");
91 if (p)
92 prot = p->PdgCode();
93
94 p = db->GetParticle("Deuteron");
95 if (p)
96 deut = p->PdgCode();
97 else {
98 deut = 1000010020;
99 db->AddParticle("Deuteron", "Deuteron", 2 * kAu2Gev + 8.071e-3, kTRUE, 0, 3, "Ion", deut);
100 }
101
102 p = db->GetParticle("neutron");
103 if (p)
104 neut = p->PdgCode();
105
106 p = db->GetParticle("Triton");
107 if (p)
108 H3 = p->PdgCode();
109 else {
110 H3 = 1000010030;
111 db->AddParticle("Triton", "Triton", 3 * kAu2Gev + 14.931e-3, kFALSE, khShGev / (12.33 * kYear2Sec), 3, "Ion",
112 H3);
113 }
114
115 p = db->GetParticle("He4L");
116 if (p)
117 He4L = p->PdgCode();
118 else {
119 He4L = 1010020040;
120 db->AddParticle("He4L", "He4L", 3.92501, kFALSE, khShGev / (12.33 * kYear2Sec), 6, "Ion", He4L);
121 gMC->DefineParticle(He4L, "He4L", kPTHadron, 3.92501, 2.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4,
122 kFALSE);
123 }
124
125 p = db->GetParticle("Alpha");
126 if (!p) {
127 p = db->GetParticle("alpha");
128 if (!p) {
129 p = db->GetParticle("He4");
130 }
131 }
132
133 if (p)
134 AlphaPDG = p->PdgCode();
135 else {
136 AlphaPDG = 1000020040;
137 db->AddParticle("Alpha", "Alpha", 4 * kAu2Gev + 2.424e-3, kTRUE, khShGev / (12.33 * kYear2Sec), 6, "Ion",
138 AlphaPDG);
139 }
140
141 cout << "!!!! 0-0 !!!!!!! " << (db->ParticleList()->GetEntries()) << " " << AlphaPDG << endl;
142
143 p = db->GetParticle("He3");
144 if (!p)
145 p = db->GetParticle("HE3");
146 if (p)
147 He3 = p->PdgCode();
148 else {
149 He3 = 1000020030;
150 db->AddParticle("HE3", "HE3", 2.80923, kFALSE, 0, 6, "Ion", He3);
151 }
152
153 cout << "!!!! 0-1 !!!!!!! " << (db->ParticleList()->GetEntries()) << " " << He3 << endl;
154
155 if (He3) {
156 p = db->GetParticle("H3L");
157 if (p)
158 H3L = p->PdgCode();
159 else {
160 H3L = 1010010030;
161 gMC->DefineParticle(H3L, "H3L", kPTHadron, 2.99131, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3,
162 kFALSE);
163 }
164
165 // ----- H3L decays -----
166 // SetMode1_for_PDG(H3L,100,-211,He3);
167 for (Int_t kz = 0; kz < 6; kz++) {
168 b[kz] = 0.;
169 m1[kz] = 0;
170 m2[kz] = 0;
171 m3[kz] = 0;
172 }
173 ip = -1;
174 totB = 0.;
175 b[++ip] = 24.7;
176 m1[ip] = -211;
177 m2[ip] = He3;
178 m3[ip] = 0; // pi-, He3
179 b[++ip] = 12.35;
180 m1[ip] = 111;
181 m2[ip] = H3;
182 m3[ip] = 0; // pi0, H3
183 b[++ip] = 36.7;
184 m1[ip] = -211;
185 m2[ip] = prot;
186 m3[ip] = deut; // pi-, p, d
187 b[++ip] = 18.35;
188 m1[ip] = 111;
189 m2[ip] = neut;
190 m3[ip] = deut; // pi0, n, d
191 b[++ip] = 0.2;
192 m1[ip] = neut;
193 m2[ip] = deut;
194 m3[ip] = 0; // n, d
195 b[++ip] = 1.5;
196 m1[ip] = neut;
197 m2[ip] = neut;
198 m3[ip] = prot; // n, n, p
199
200 for (Int_t kz = 0; kz < ip + 1; kz++)
201 totB += b[kz] / 100.;
202 for (Int_t kz = 0; kz < ip + 1; kz++)
203 b[kz] /= totB;
204
205 SetMode_for_PDG(H3L, b, m1, m2, m3, ++ip);
206
207 p = db->GetParticle("H4L");
208 if (p)
209 H4L = p->PdgCode();
210 else {
211 H4L = 1010010040;
212 gMC->DefineParticle(H4L, "H4L", kPTHadron, 3.92503, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4,
213 kFALSE);
214 }
215
216 // ----- H4L decays -----
217 // SetMode1_for_PDG(H4L,100,-211,AlphaPDG);
218 for (Int_t kz = 0; kz < 6; kz++) {
219 b[kz] = 0.;
220 m1[kz] = 0;
221 m2[kz] = 0;
222 m3[kz] = 0;
223 }
224 ip = -1;
225 totB = 0.;
226 b[++ip] = 75.;
227 m1[ip] = -211;
228 m2[ip] = AlphaPDG;
229 m3[ip] = 0; // pi-, He4
230 b[++ip] = 25.;
231 m1[ip] = H3;
232 m2[ip] = prot;
233 m3[ip] = -211; // H3, prot, pi-
234
235 for (Int_t kz = 0; kz < ip + 1; kz++)
236 totB += b[kz] / 100.;
237 for (Int_t kz = 0; kz < ip + 1; kz++)
238 b[kz] /= totB;
239
240 SetMode_for_PDG(H4L, b, m1, m2, m3, ++ip);
241 }
242
243 // ----- He4L decays -----
244 for (Int_t kz = 0; kz < 6; kz++) {
245 b[kz] = 0.;
246 m1[kz] = 0;
247 m2[kz] = 0;
248 m3[kz] = 0;
249 }
250 ip = -1;
251 totB = 0.;
252 b[++ip] = 32.;
253 m1[ip] = -211;
254 m2[ip] = He3;
255 m3[ip] = prot; // pi-, He3, p
256 b[++ip] = 35.;
257 m1[ip] = 111;
258 m2[ip] = AlphaPDG;
259 m3[ip] = 0; // pi0, He4
260 b[++ip] = 14.;
261 m1[ip] = 111;
262 m2[ip] = He3;
263 m3[ip] = neut; // pi0, He3, neut
264 b[++ip] = 19.;
265 m1[ip] = prot;
266 m2[ip] = H3;
267 m3[ip] = 0; // p, H3
268
269 for (Int_t kz = 0; kz < ip + 1; kz++)
270 totB += b[kz] / 100.;
271 for (Int_t kz = 0; kz < ip + 1; kz++)
272 b[kz] /= totB;
273
274 SetMode_for_PDG(He4L, b, m1, m2, m3, ++ip);
275
276 p = db->GetParticle("LN");
277 if (!p) {
278 p = db->GetParticle("LambdaNeutron");
279 if (!p) {
280 gMC->DefineParticle(1010000020, "LN", kPTNeutron, 2.054, 0.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 2,
281 kFALSE);
282 p = db->GetParticle("LN");
283 }
284 }
285 if (p)
286 SetMode1_for_PDG(p->PdgCode(), 100, 1000010020, -211);
287
288 cout << "!!!! 1-0 !!!!!!! " << (db->ParticleList()->GetEntries()) << " " << H3L << " " << H4L << endl;
289}
Int_t
int Int_t
Float_t
float Float_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
TDatabasePDG::Instance
static TDatabasePDG * Instance()
TDatabasePDG::ParticleList
const THashList * ParticleList() const
TDatabasePDG::GetParticle
TParticlePDG * GetParticle(const char *name) const
TDatabasePDG::AddParticle
virtual TParticlePDG * AddParticle(const char *Name, const char *Title, Double_t Mass, Bool_t Stable, Double_t DecayWidth, Double_t Charge, const char *ParticleClass, Int_t PdgCode, Int_t Anti=-1, Int_t TrackingCode=0)
TParticlePDG::PdgCode
Int_t PdgCode() const
SetMode_for_PDG
void SetMode_for_PDG(Int_t pdg, Float_t *b, Int_t *m0, Int_t *m1, Int_t *m2, Int_t N)
Definition UserDecay.C:52
UserDecayConfig
void UserDecayConfig()
Definition UserDecay.C:75
SetMode2_for_PDG
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)
Definition UserDecay.C:20
SetMode1_for_PDG
void SetMode1_for_PDG(Int_t pdg, Float_t b0, Int_t m00, Int_t m01)
Definition UserDecay.C:1