ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtTPC20MgDecay.cxx
Go to the documentation of this file.
1 #include "AtTPC20MgDecay.h"
2 
3 #include <FairLogger.h> // for Logger, LOG
4 #include <FairPrimaryGenerator.h>
5 
6 #include <TF1.h>
7 #include <TH1.h>
8 #include <TMath.h>
9 #include <TRandom.h>
10 
11 #include <cmath> // for acos
12 #include <iostream> // for operator<<, endl, basic_ostream, cout
13 #include <memory> // for make_unique, unique_ptr
14 
15 Bool_t AtTPC20MgDecay::Init()
16 {
17  // Initialize generatorTH1F*h1 = new TH1F("h1", "h1", 1000,0,11236.8);
18  return true;
19 }
20 
21 // ----- Public method ReadEvent --------------------------------------
22 Bool_t AtTPC20MgDecay::ReadEvent(FairPrimaryGenerator *primGen)
23 {
24 
25  if (fBoxVtxIsSet) {
26  fX = gRandom->Uniform(fX1, fX2);
27  fY = gRandom->Uniform(fY1, fY2);
28  fZ = gRandom->Uniform(fZ1, fZ2);
29  }
30 
31  // Proton of 1210keV and alpha of 506keV
32  Int_t protonPDGID = 2212;
33  Int_t alphaPDGID = 1000020040;
34  Int_t gammaPDGID = 22;
35  Int_t betaPDGID = 11;
36  // Check for particle type
37 
38  // Protons from the decay of 20Na to 19Ne
39  // Double32_t kinEneProton = 0.001210; //GeV
40  // Double32_t kinEneAlpha = 0.000506; //GeV
41  // Double32_t kinEneGamma =0.004033; //GeV and it has zero rest mass
42  Double32_t ptProton = 0, pxProton = 0, pyProton = 0, pzProton = 0;
43  Double32_t pabsProton = 0.0469; // GeV/c , 1.2 MeV
44  Double32_t thetaProton = acos(gRandom->Uniform(-1, 1));
45  Double32_t brp = 0;
46  Double32_t phiProton = gRandom->Uniform(0, 360) * TMath::DegToRad();
47  pzProton = pabsProton * TMath::Cos(thetaProton);
48  ptProton = pabsProton * TMath::Sin(thetaProton);
49  pxProton = ptProton * TMath::Cos(phiProton);
50  pyProton = ptProton * TMath::Sin(phiProton);
51 
52  Double32_t ptAlpha = 0, pxAlpha = 0, pyAlpha = 0, pzAlpha = 0;
53  // Double32_t bra=0;
54  Double32_t pabsAlpha = 0.06162; // GeV/c, 506 keV from decay of 19Ne to 15O
55  Double32_t thetaAlpha = acos(gRandom->Uniform(-1, 1));
56  Double32_t phiAlpha = gRandom->Uniform(0, 360) * TMath::DegToRad();
57  pzAlpha = pabsAlpha * TMath::Cos(thetaAlpha);
58  ptAlpha = pabsAlpha * TMath::Sin(thetaAlpha);
59  pxAlpha = ptAlpha * TMath::Cos(phiAlpha);
60  pyAlpha = ptAlpha * TMath::Sin(phiAlpha);
61 
62  Double32_t ptProton1 = 0, pxProton1 = 0, pyProton1 = 0, pzProton1 = 0;
63  Double32_t pabsProton1 = 0.0389; // GeV/c , 806 keV
64  Double32_t thetaProton1 = acos(gRandom->Uniform(-1, 1));
65  Double32_t phiProton1 = gRandom->Uniform(0, 360) * TMath::DegToRad();
66  pzProton1 = pabsProton1 * TMath::Cos(thetaProton1);
67  ptProton1 = pabsProton1 * TMath::Sin(thetaProton1);
68  pxProton1 = ptProton1 * TMath::Cos(phiProton1);
69  pyProton1 = ptProton1 * TMath::Sin(phiProton1);
70 
71  Double32_t ptProton2 = 0, pxProton2 = 0, pyProton2 = 0, pzProton2 = 0;
72  Double32_t pabsProton2 = 0.04476; // GeV/c, 1056 keV
73  Double32_t thetaProton2 = acos(gRandom->Uniform(-1, 1));
74  Double32_t phiProton2 = gRandom->Uniform(0, 360) * TMath::DegToRad();
75  pzProton2 = pabsProton2 * TMath::Cos(thetaProton2);
76  ptProton2 = pabsProton2 * TMath::Sin(thetaProton2);
77  pxProton2 = ptProton2 * TMath::Cos(phiProton2);
78  pyProton2 = ptProton2 * TMath::Sin(phiProton2);
79 
80  Double32_t ptProton3 = 0, pxProton3 = 0, pyProton3 = 0, pzProton3 = 0;
81  Double32_t pabsProton3 = 0.05169; // GeV/c, 1416 keV
82  Double32_t thetaProton3 = acos(gRandom->Uniform(-1, 1));
83  Double32_t phiProton3 = gRandom->Uniform(0, 360) * TMath::DegToRad();
84  pzProton3 = pabsProton3 * TMath::Cos(thetaProton3);
85  ptProton3 = pabsProton3 * TMath::Sin(thetaProton3);
86  pxProton3 = ptProton3 * TMath::Cos(phiProton3);
87  pyProton3 = ptProton3 * TMath::Sin(phiProton3);
88 
89  Double32_t ptProton4 = 0, pxProton4 = 0, pyProton4 = 0, pzProton4 = 0;
90  Double32_t pabsProton4 = 0.05636; // GeV/c, 1679 keV
91  Double32_t thetaProton4 = acos(gRandom->Uniform(-1, 1));
92  Double32_t phiProton4 = gRandom->Uniform(0, 360) * TMath::DegToRad();
93  pzProton4 = pabsProton4 * TMath::Cos(thetaProton4);
94  ptProton4 = pabsProton4 * TMath::Sin(thetaProton4);
95  pxProton4 = ptProton4 * TMath::Cos(phiProton4);
96  pyProton4 = ptProton4 * TMath::Sin(phiProton4);
97 
98  Double32_t ptProton5 = 0, pxProton5 = 0, pyProton5 = 0, pzProton5 = 0;
99  Double32_t pabsProton5 = 0.06018; // GeV/c, 1928 keV
100  Double32_t thetaProton5 = acos(gRandom->Uniform(-1, 1));
101  Double32_t phiProton5 = gRandom->Uniform(0, 360) * TMath::DegToRad();
102  pzProton5 = pabsProton5 * TMath::Cos(thetaProton5);
103  ptProton5 = pabsProton5 * TMath::Sin(thetaProton5);
104  pxProton5 = ptProton5 * TMath::Cos(phiProton5);
105  pyProton5 = ptProton5 * TMath::Sin(phiProton5);
106 
107  Double32_t ptProton6 = 0, pxProton6 = 0, pyProton6 = 0, pzProton6 = 0;
108  Double32_t pabsProton6 = 0.0651; // GeV/c, 2256 keV
109  Double32_t thetaProton6 = acos(gRandom->Uniform(-1, 1));
110  Double32_t phiProton6 = gRandom->Uniform(0, 360) * TMath::DegToRad();
111  pzProton6 = pabsProton6 * TMath::Cos(thetaProton6);
112  ptProton6 = pabsProton6 * TMath::Sin(thetaProton6);
113  pxProton6 = ptProton6 * TMath::Cos(phiProton6);
114  pyProton6 = ptProton6 * TMath::Sin(phiProton6);
115 
116  Double32_t ptProton7 = 0, pxProton7 = 0, pyProton7 = 0, pzProton7 = 0;
117  Double32_t pabsProton7 = 0.06636; // GeV/c, 2344 keV
118  Double32_t thetaProton7 = acos(gRandom->Uniform(-1, 1));
119  Double32_t phiProton7 = gRandom->Uniform(0, 360) * TMath::DegToRad();
120  pzProton7 = pabsProton7 * TMath::Cos(thetaProton7);
121  ptProton7 = pabsProton7 * TMath::Sin(thetaProton7);
122  pxProton7 = ptProton7 * TMath::Cos(phiProton7);
123  pyProton7 = ptProton7 * TMath::Sin(phiProton7);
124 
125  Double32_t ptProton8 = 0, pxProton8 = 0, pyProton8 = 0, pzProton8 = 0;
126  Double32_t pabsProton8 = 0.06934; // GeV/c, 2559 keV
127  Double32_t thetaProton8 = acos(gRandom->Uniform(-1, 1));
128  Double32_t phiProton8 = gRandom->Uniform(0, 360) * TMath::DegToRad();
129  pzProton8 = pabsProton8 * TMath::Cos(thetaProton8);
130  ptProton8 = pabsProton8 * TMath::Sin(thetaProton8);
131  pxProton8 = ptProton8 * TMath::Cos(phiProton8);
132  pyProton8 = ptProton8 * TMath::Sin(phiProton8);
133 
134  Double32_t ptProton9 = 0, pxProton9 = 0, pyProton9 = 0, pzProton9 = 0;
135  Double32_t pabsProton9 = 0.07513; // GeV/c, 3003 keV
136  Double32_t thetaProton9 = acos(gRandom->Uniform(-1, 1));
137  Double32_t phiProton9 = gRandom->Uniform(0, 360) * TMath::DegToRad();
138  pzProton9 = pabsProton9 * TMath::Cos(thetaProton9);
139  ptProton9 = pabsProton9 * TMath::Sin(thetaProton9);
140  pxProton9 = ptProton9 * TMath::Cos(phiProton9);
141  pyProton9 = ptProton9 * TMath::Sin(phiProton9);
142 
143  Double32_t ptProton10 = 0, pxProton10 = 0, pyProton10 = 0, pzProton10 = 0;
144  Double32_t pabsProton10 = 0.07982; // GeV/c, 3389 keV
145  Double32_t thetaProton10 = acos(gRandom->Uniform(-1, 1));
146  Double32_t phiProton10 = gRandom->Uniform(0, 360) * TMath::DegToRad();
147  pzProton10 = pabsProton10 * TMath::Cos(thetaProton10);
148  ptProton10 = pabsProton10 * TMath::Sin(thetaProton10);
149  pxProton10 = ptProton10 * TMath::Cos(phiProton10);
150  pyProton10 = ptProton10 * TMath::Sin(phiProton10);
151 
152  Double32_t ptProton11 = 0, pxProton11 = 0, pyProton11 = 0, pzProton11 = 0;
153  Double32_t pabsProton11 = 0.08475; // GeV/c, 3820 keV
154  Double32_t thetaProton11 = acos(gRandom->Uniform(-1, 1));
155  Double32_t phiProton11 = gRandom->Uniform(0, 360) * TMath::DegToRad();
156  pzProton11 = pabsProton11 * TMath::Cos(thetaProton11);
157  ptProton11 = pabsProton11 * TMath::Sin(thetaProton11);
158  pxProton11 = ptProton11 * TMath::Cos(phiProton11);
159  pyProton11 = ptProton11 * TMath::Sin(phiProton11);
160 
161  Double32_t ptProton12 = 0, pxProton12 = 0, pyProton12 = 0, pzProton12 = 0;
162  Double32_t pabsProton12 = 0.0875; // GeV/c, 4071 keV
163  Double32_t thetaProton12 = acos(gRandom->Uniform(-1, 1));
164  Double32_t phiProton12 = gRandom->Uniform(0, 360) * TMath::DegToRad();
165  pzProton12 = pabsProton12 * TMath::Cos(thetaProton12);
166  ptProton12 = pabsProton12 * TMath::Sin(thetaProton12);
167  pxProton12 = ptProton12 * TMath::Cos(phiProton12);
168  pyProton12 = ptProton12 * TMath::Sin(phiProton12);
169 
170  Double32_t ptProton13 = 0, pxProton13 = 0, pyProton13 = 0, pzProton13 = 0;
171  Double32_t pabsProton13 = 0.0902; // GeV/c, 4326 keV
172  Double32_t thetaProton13 = acos(gRandom->Uniform(-1, 1));
173  Double32_t phiProton13 = gRandom->Uniform(0, 360) * TMath::DegToRad();
174  pzProton13 = pabsProton13 * TMath::Cos(thetaProton13);
175  ptProton13 = pabsProton13 * TMath::Sin(thetaProton13);
176  pxProton13 = ptProton13 * TMath::Cos(phiProton13);
177  pyProton13 = ptProton13 * TMath::Sin(phiProton13);
178 
179  // beta particles emitted from 20Mg decay to 20Na which will undergo proton emission
180  auto h1 = std::make_unique<TH1F>("h1", "h1", 1000, 0, 8650); // 984 keV
181  auto f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 8650);
182  f1->SetParameters(1.4373e-6, 2.8193e-8, 1.7812e-11, -5.31546e-15, 3.35539e-19, -3.31743e-25);
183  h1->FillRandom("f1");
184  h1->Draw();
185  Double32_t r1 = f1->GetRandom();
186  // std::cout<<"r1="<<r1<<std::endl;
187  Double32_t ptBeta1 = 0, pxBeta1 = 0, pyBeta1 = 0, pzBeta1 = 0;
188  Double32_t pabsBeta1 = (TMath::Sqrt((r1 + 511) * (r1 + 511) - (511 * 511))) * 1e-6;
189 
190  // std::cout<<"pabsBeta1="<<pabsBeta1<<std::endl;
191 
192  Double32_t brbeta = 0;
193  Double32_t brb = 0;
194  Double32_t thetaBeta1 = acos(gRandom->Uniform(-1, 1));
195  Double32_t phiBeta1 = gRandom->Uniform(0, 360) * TMath::DegToRad();
196  pzBeta1 = pabsBeta1 * TMath::Cos(thetaBeta1);
197  ptBeta1 = pabsBeta1 * TMath::Sin(thetaBeta1);
198  pxBeta1 = ptBeta1 * TMath::Cos(phiBeta1);
199  pyBeta1 = ptBeta1 * TMath::Sin(phiBeta1);
200 
201  h1 = std::make_unique<TH1F>("h2", "h2", 1000, 0, 7000); // 2645 keV
202  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 7000);
203  f1->SetParameters(3.21795e-9, 7.31673e-11, 4.05604e-14, -1.61734e-17, 1.25431e-21, 2.00144e-27);
204  h1->FillRandom("f1");
205  h1->Draw();
206  Double32_t r2 = f1->GetRandom();
207  // std::cout<<"r2="<<r2<<std::endl;
208 
209  Double32_t ptBeta2 = 0, pxBeta2 = 0, pyBeta2 = 0, pzBeta2 = 0;
210  Double32_t pabsBeta2 = (TMath::Sqrt((r2 + 511) * (r2 + 511) - (511 * 511))) * 1e-6;
211 
212  // std::cout<<"pabsBeta2="<<pabsBeta2<<std::endl;
213 
214  Double32_t thetaBeta2 = acos(gRandom->Uniform(-1, 1));
215  Double32_t phiBeta2 = gRandom->Uniform(0, 360) * TMath::DegToRad();
216  pzBeta2 = pabsBeta2 * TMath::Cos(thetaBeta2);
217  ptBeta2 = pabsBeta2 * TMath::Sin(thetaBeta2);
218  pxBeta2 = ptBeta2 * TMath::Cos(phiBeta2);
219  pyBeta2 = ptBeta2 * TMath::Sin(phiBeta2);
220 
221  h1 = std::make_unique<TH1F>("h3", "h3", 1000, 0, 6600); // 3001 kev
222  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 6600);
223  f1->SetParameters(4.12126e-7, 9.72611e-9, 5.16976e-12, -2.22576e-15, 1.81306e-19, 4.8888e-25);
224  h1->FillRandom("f1");
225  h1->Draw();
226  Double32_t r3 = f1->GetRandom();
227  // std::cout<<"r3="<<r3<<std::endl;
228 
229  Double32_t ptBeta3 = 0, pxBeta3 = 0, pyBeta3 = 0, pzBeta3 = 0;
230  Double32_t pabsBeta3 = (TMath::Sqrt((r3 + 511) * (r3 + 511) - (511 * 511))) * 1e-6;
231 
232  // std::cout<<"pabsBeta3="<<pabsBeta3<<std::endl;
233 
234  Double32_t thetaBeta3 = acos(gRandom->Uniform(-1, 1));
235  Double32_t phiBeta3 = gRandom->Uniform(0, 360) * TMath::DegToRad();
236  pzBeta3 = pabsBeta3 * TMath::Cos(thetaBeta3);
237  ptBeta3 = pabsBeta3 * TMath::Sin(thetaBeta3);
238  pxBeta3 = ptBeta3 * TMath::Cos(phiBeta3);
239  pyBeta3 = ptBeta3 * TMath::Sin(phiBeta3);
240 
241  h1 = std::make_unique<TH1F>("h4", "h4", 1000, 0, 5750); // 3874 keV
242  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 5750);
243  f1->SetParameters(2.28812e-7, 5.99596e-9, 2.76505e-12, -1.49003e-15, 1.37999e-19, 1.00668e-24);
244  h1->FillRandom("f1");
245  h1->Draw();
246  Double32_t r4 = f1->GetRandom();
247  // std::cout<<"r4="<<r4<<std::endl;
248 
249  Double32_t ptBeta4 = 0, pxBeta4 = 0, pyBeta4 = 0, pzBeta4 = 0;
250  Double32_t pabsBeta4 = (TMath::Sqrt((r4 + 511) * (r4 + 511) - (511 * 511))) * 1e-6;
251 
252  // std::cout<<"pabsBeta4="<<pabsBeta4<<std::endl;
253 
254  Double32_t thetaBeta4 = acos(gRandom->Uniform(-1, 1));
255  Double32_t phiBeta4 = gRandom->Uniform(0, 360) * TMath::DegToRad();
256  pzBeta4 = pabsBeta4 * TMath::Cos(thetaBeta4);
257  ptBeta4 = pabsBeta4 * TMath::Sin(thetaBeta4);
258  pxBeta4 = ptBeta4 * TMath::Cos(phiBeta4);
259  pyBeta4 = ptBeta4 * TMath::Sin(phiBeta4);
260 
261  h1 = std::make_unique<TH1F>("h5", "h5", 1000, 0, 5500); // 4123 keV
262  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 5500);
263  f1->SetParameters(1.40485e-7, 3.80994e-9, 1.66204e-12, -9.67707e-16, 9.31417e-20, 8.78505e-25);
264  h1->FillRandom("f1");
265  h1->Draw();
266  Double32_t r5 = f1->GetRandom();
267  // std::cout<<"r5="<<r5<<std::endl;
268 
269  Double32_t ptBeta5 = 0, pxBeta5 = 0, pyBeta5 = 0, pzBeta5 = 0;
270  Double32_t pabsBeta5 = (TMath::Sqrt((r5 + 511) * (r5 + 511) - (511 * 511))) * 1e-6;
271 
272  // std::cout<<"pabsBeta5="<<pabsBeta5<<std::endl;
273 
274  Double32_t thetaBeta5 = acos(gRandom->Uniform(-1, 1));
275  Double32_t phiBeta5 = gRandom->Uniform(0, 360) * TMath::DegToRad();
276  pzBeta5 = pabsBeta5 * TMath::Cos(thetaBeta5);
277  ptBeta5 = pabsBeta5 * TMath::Sin(thetaBeta5);
278  pxBeta5 = ptBeta5 * TMath::Cos(phiBeta5);
279  pyBeta5 = ptBeta5 * TMath::Sin(phiBeta5);
280 
281  h1 = std::make_unique<TH1F>("h6", "h6", 1000, 0, 4810); // 4800 keV
282  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 4810);
283  f1->SetParameters(1.27496e-7, 3.84163e-9, 1.35989e-12, -1.02747e-15, 1.10166e-19, 2.06696e-24);
284  h1->FillRandom("f1");
285  h1->Draw();
286  Double32_t r6 = f1->GetRandom();
287  // std::cout<<"r6="<<r6<<std::endl;
288 
289  Double32_t ptBeta6 = 0, pxBeta6 = 0, pyBeta6 = 0, pzBeta6 = 0;
290  Double32_t pabsBeta6 = (TMath::Sqrt((r6 + 511) * (r6 + 511) - (511 * 511))) * 1e-6;
291 
292  // std::cout<<"pabsBeta6="<<pabsBeta6<<std::endl;
293  Double32_t thetaBeta6 = acos(gRandom->Uniform(-1, 1));
294  Double32_t phiBeta6 = gRandom->Uniform(0, 360) * TMath::DegToRad();
295  pzBeta6 = pabsBeta6 * TMath::Cos(thetaBeta6);
296  ptBeta6 = pabsBeta6 * TMath::Sin(thetaBeta6);
297  pxBeta6 = ptBeta6 * TMath::Cos(phiBeta6);
298  pyBeta6 = ptBeta6 * TMath::Sin(phiBeta6);
299 
300  h1 = std::make_unique<TH1F>("h7", "h7", 1000, 0, 4010); // 5600
301  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 4010);
302  f1->SetParameters(1.40707e-7, 4.9678e-9, 1.07582e-12, -1.36341e-15, 1.65065e-19, 7.62092e-24);
303  h1->FillRandom("f1");
304  h1->Draw();
305  Double32_t r7 = f1->GetRandom();
306  // std::cout<<"r7="<<r7<<std::endl;
307 
308  Double32_t ptBeta7 = 0, pxBeta7 = 0, pyBeta7 = 0, pzBeta7 = 0;
309  Double32_t pabsBeta7 = (TMath::Sqrt((r7 + 511) * (r7 + 511) - (511 * 511))) * 1e-6;
310 
311  // std::cout<<"pabsBeta7="<<pabsBeta7<<std::endl;
312 
313  Double32_t thetaBeta7 = acos(gRandom->Uniform(-1, 1));
314  Double32_t phiBeta7 = gRandom->Uniform(0, 360) * TMath::DegToRad();
315  pzBeta7 = pabsBeta7 * TMath::Cos(thetaBeta7);
316  ptBeta7 = pabsBeta7 * TMath::Sin(thetaBeta7);
317  pxBeta7 = ptBeta7 * TMath::Cos(phiBeta7);
318  pyBeta7 = ptBeta7 * TMath::Sin(phiBeta7);
319 
320  h1 = std::make_unique<TH1F>("h8", "h8", 1000, 0, 3770.); // 5836 keV
321  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 3770);
322  f1->SetParameters(7.89154e-7, 2.94918e-8, 4.8237e-12, -8.02751e-15, 9.98793e-19, 6.29715e-23);
323  h1->FillRandom("f1");
324  h1->Draw();
325  Double32_t r8 = f1->GetRandom();
326  // std::cout<<"r8="<<r8<<std::endl;
327 
328  Double32_t ptBeta8 = 0, pxBeta8 = 0, pyBeta8 = 0, pzBeta8 = 0;
329  Double32_t pabsBeta8 = (TMath::Sqrt((r8 + 511) * (r8 + 511) - (511 * 511))) * 1e-6;
330 
331  // std::cout<<"pabsBeta8="<<pabsBeta8<<std::endl;
332  Double32_t thetaBeta8 = acos(gRandom->Uniform(-1, 1));
333  Double32_t phiBeta8 = gRandom->Uniform(0, 360) * TMath::DegToRad();
334  pzBeta8 = pabsBeta8 * TMath::Cos(thetaBeta8);
335  ptBeta8 = pabsBeta8 * TMath::Sin(thetaBeta8);
336  pxBeta8 = ptBeta8 * TMath::Cos(phiBeta8);
337  pyBeta8 = ptBeta8 * TMath::Sin(phiBeta8);
338 
339  h1 = std::make_unique<TH1F>("h9", "h9", 1000, 0, 3340); // 6266 keV
340  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 3340);
341  f1->SetParameters(1.5398e-7, 6.4777e-9, 3.04221e-13, -1.67316e-15, 2.11676e-19, 2.6079e-23);
342  h1->FillRandom("f1");
343  h1->Draw();
344  Double32_t r9 = f1->GetRandom();
345  // std::cout<<"r9="<<r9<<std::endl;
346 
347  Double32_t ptBeta9 = 0, pxBeta9 = 0, pyBeta9 = 0, pzBeta9 = 0;
348  Double32_t pabsBeta9 = (TMath::Sqrt((r9 + 511) * (r9 + 511) - (511 * 511))) * 1e-6;
349 
350  // std::cout<<"pabsBeta9="<<pabsBeta9<<std::endl;
351 
352  Double32_t thetaBeta9 = acos(gRandom->Uniform(-1, 1));
353  Double32_t phiBeta9 = gRandom->Uniform(0, 360) * TMath::DegToRad();
354  pzBeta9 = pabsBeta9 * TMath::Cos(thetaBeta9);
355  ptBeta9 = pabsBeta9 * TMath::Sin(thetaBeta9);
356  pxBeta9 = ptBeta9 * TMath::Cos(phiBeta9);
357  pyBeta9 = ptBeta9 * TMath::Sin(phiBeta9);
358 
359  h1 = std::make_unique<TH1F>("h10", "h10", 1000, 0, 3071); // 6534
360  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 3071);
361  f1->SetParameters(4.84941e-7, 2.22562e-8, -9.262068e-13, -5.32068e-15, 6.40919e-19, 1.36991e-22);
362  h1->FillRandom("f1");
363  h1->Draw();
364 
365  Double32_t r10 = f1->GetRandom();
366  // std::cout<<"r10="<<r10<<std::endl;
367 
368  Double32_t ptBeta10 = 0, pxBeta10 = 0, pyBeta10 = 0, pzBeta10 = 0;
369  Double32_t pabsBeta10 = (TMath::Sqrt((r10 + 511) * (r10 + 511) - (511 * 511))) * 1e-6;
370 
371  // std::cout<<"pabsBeta10="<<pabsBeta1<<std::endl;
372  Double32_t thetaBeta10 = acos(gRandom->Uniform(-1, 1));
373  Double32_t phiBeta10 = gRandom->Uniform(0, 360) * TMath::DegToRad();
374  pzBeta10 = pabsBeta10 * TMath::Cos(thetaBeta10);
375  ptBeta10 = pabsBeta10 * TMath::Sin(thetaBeta10);
376  pxBeta10 = ptBeta10 * TMath::Cos(phiBeta10);
377  pyBeta10 = ptBeta10 * TMath::Sin(phiBeta10);
378 
379  h1 = std::make_unique<TH1F>("h11", "h11", 1000, 0, 2835); // 6770 keV
380  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 2835);
381  f1->SetParameters(1.33516e-6, 6.68463e-8, -9.17265e-12, -1.41377e-14, 1.4453e-18, 6.10901e-22);
382  h1->FillRandom("f1");
383  h1->Draw();
384  Double32_t r11 = f1->GetRandom();
385  // std::cout<<"r11="<<r11<<std::endl;
386 
387  Double32_t ptBeta11 = 0, pxBeta11 = 0, pyBeta11 = 0, pzBeta11 = 0;
388  Double32_t pabsBeta11 = (TMath::Sqrt((r11 + 511) * (r11 + 511) - (511 * 511))) * 1e-6;
389 
390  // std::cout<<"pabsBeta11="<<pabsBeta11<<std::endl;
391  Double32_t thetaBeta11 = acos(gRandom->Uniform(-1, 1));
392  Double32_t phiBeta11 = gRandom->Uniform(0, 360) * TMath::DegToRad();
393  pzBeta11 = pabsBeta11 * TMath::Cos(thetaBeta11);
394  ptBeta11 = pabsBeta11 * TMath::Sin(thetaBeta11);
395  pxBeta11 = ptBeta11 * TMath::Cos(phiBeta11);
396  pyBeta11 = ptBeta11 * TMath::Sin(phiBeta11);
397 
398  h1 = std::make_unique<TH1F>("h12", "h12", 1000, 0, 2685); // 6920 keV
399  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 2685);
400  f1->SetParameters(9.03068e-7, 4.81333e-8, -9.93371e-12, -8.9499e-15, 6.66662e-19, 5.74898e-22);
401  h1->FillRandom("f1");
402  h1->Draw();
403  Double32_t r12 = f1->GetRandom();
404  // std::cout<<"r12="<<r12<<std::endl;
405 
406  Double32_t ptBeta12 = 0, pxBeta12 = 0, pyBeta12 = 0, pzBeta12 = 0;
407  Double32_t pabsBeta12 = (TMath::Sqrt((r12 + 511) * (r12 + 511) - (511 * 511))) * 1e-6;
408 
409  // std::cout<<"pabsBeta12="<<pabsBeta12<<std::endl;
410 
411  Double32_t thetaBeta12 = acos(gRandom->Uniform(-1, 1));
412  Double32_t phiBeta12 = gRandom->Uniform(0, 360) * TMath::DegToRad();
413  pzBeta12 = pabsBeta12 * TMath::Cos(thetaBeta12);
414  ptBeta12 = pabsBeta12 * TMath::Sin(thetaBeta12);
415  pxBeta12 = ptBeta12 * TMath::Cos(phiBeta12);
416  pyBeta12 = ptBeta12 * TMath::Sin(phiBeta12);
417 
418  h1 = std::make_unique<TH1F>("h13", "h13", 1000, 0, 2165); // 7440 keV
419  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 2165);
420  f1->SetParameters(2.42189e-9, 1.68134e-10, -8.74415e-14, -3.26775e-18, -8.54988e-21, 5.56934e-24);
421  h1->FillRandom("f1");
422  h1->Draw();
423  Double32_t r13 = f1->GetRandom();
424  // std::cout<<"r13="<<r13<<std::endl;
425 
426  Double32_t ptBeta13 = 0, pxBeta13 = 0, pyBeta13 = 0, pzBeta13 = 0;
427  Double32_t pabsBeta13 = (TMath::Sqrt((r13 + 511) * (r13 + 511) - (511 * 511))) * 1e-6;
428  // std::cout<<"pabsBeta13="<<pabsBeta13<<std::endl;
429 
430  Double32_t thetaBeta13 = acos(gRandom->Uniform(-1, 1));
431  Double32_t phiBeta13 = gRandom->Uniform(0, 360) * TMath::DegToRad();
432  pzBeta13 = pabsBeta13 * TMath::Cos(thetaBeta13);
433  ptBeta13 = pabsBeta13 * TMath::Sin(thetaBeta13);
434  pxBeta13 = ptBeta13 * TMath::Cos(phiBeta13);
435  pyBeta13 = ptBeta13 * TMath::Sin(phiBeta13);
436 
437  // Beta particles from decay of 20Na to 20Ne
438 
439  h1 = std::make_unique<TH1F>("h14", "h14", 1000, 0, 11236.8);
440  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 11236.8);
441  f1->SetParameters(1.59711e-7, 1.68487e-8, 1.0013e-11, -2.1646e-15, 9.92165e-20, 2.01015e-25);
442  h1->FillRandom("f1");
443  h1->Draw();
444  Double32_t r14 = f1->GetRandom();
445  // std::cout<<"r1="<<r1<<std::endl;
446 
447  Double32_t ptBeta14 = 0, pxBeta14 = 0, pyBeta14 = 0, pzBeta14 = 0;
448  Double32_t pabsBeta14 = (TMath::Sqrt((r14 + 511) * (r14 + 511) - (511 * 511))) * 1e-6;
449  // std::cout<<"pabsBeta14="<<pabsBeta14<<std::endl;
450 
451  Double32_t thetaBeta14 = acos(gRandom->Uniform(-1, 1));
452  Double32_t phiBeta14 = gRandom->Uniform(0, 360) * TMath::DegToRad();
453  pzBeta14 = pabsBeta14 * TMath::Cos(thetaBeta14);
454  ptBeta14 = pabsBeta14 * TMath::Sin(thetaBeta14);
455  pxBeta14 = ptBeta14 * TMath::Cos(phiBeta14);
456  pyBeta14 = ptBeta14 * TMath::Sin(phiBeta14);
457 
458  h1 = std::make_unique<TH1F>("h15", "h15", 1000, 0, 5448.6); // 11320 keV
459  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 5448.6);
460  f1->SetParameters(3.1588e-8, 1.25472e-9, 6.2571e-14, -3.16054e-16, 3.84227e-20, 5.00163e-24);
461  h1->FillRandom("f1");
462  h1->Draw();
463  Double32_t r15 = f1->GetRandom();
464 
465  // std::cout<<"r15="<<r15<<std::endl;
466  Double32_t ptBeta15 = 0, pxBeta15 = 0, pyBeta15 = 0, pzBeta15 = 0;
467  Double32_t pabsBeta15 = (TMath::Sqrt((r15 + 511) * (r15 + 511) - (511 * 511))) * 1e-6;
468  // std::cout<<"pabsBeta15="<<pabsBeta15<<std::endl;
469 
470  Double32_t thetaBeta15 = acos(gRandom->Uniform(-1, 1));
471  Double32_t phiBeta15 = gRandom->Uniform(0, 360) * TMath::DegToRad();
472  pzBeta15 = pabsBeta15 * TMath::Cos(thetaBeta15);
473  ptBeta15 = pabsBeta15 * TMath::Sin(thetaBeta15);
474  pxBeta15 = ptBeta15 * TMath::Cos(phiBeta15);
475  pyBeta15 = ptBeta15 * TMath::Sin(phiBeta15);
476 
477  h1 = std::make_unique<TH1F>("h16", "h16", 1000, 0, 5041.5);
478  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 5041.5);
479  f1->SetParameters(4.55e-9, 2.01854e-10, -1.63887e-14, -4.4892e-17, 4.87227e-21, 1.46295e-24);
480  h1->FillRandom("f1");
481  h1->Draw();
482  Double32_t r16 = f1->GetRandom();
483  // std::cout<<"r16="<<r16<<std::endl;
484 
485  Double32_t ptBeta16 = 0, pxBeta16 = 0, pyBeta16 = 0, pzBeta16 = 0;
486  Double32_t pabsBeta16 = (TMath::Sqrt((r16 + 511) * (r16 + 511) - (511 * 511))) * 1e-6;
487  // std::cout<<"pabsBeta6="<<pabsBeta6<<std::endl;
488  Double32_t thetaBeta16 = acos(gRandom->Uniform(-1, 1));
489  Double32_t phiBeta16 = gRandom->Uniform(0, 360) * TMath::DegToRad();
490  pzBeta16 = pabsBeta16 * TMath::Cos(thetaBeta16);
491  ptBeta16 = pabsBeta16 * TMath::Sin(thetaBeta16);
492  pxBeta16 = ptBeta16 * TMath::Cos(phiBeta16);
493  pyBeta16 = ptBeta16 * TMath::Sin(phiBeta16);
494 
495  h1 = std::make_unique<TH1F>("h17", "h17", 1000, 0, 4070.5);
496  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 4070.5);
497  f1->SetParameters(5.85765e-7, 3.02486e-8, -7.97958e-12, -4.63986e-15, 7.78501e-20, 4.44419e-22);
498  h1->FillRandom("f1");
499  h1->Draw();
500  Double32_t r17 = f1->GetRandom();
501  // std::cout<<"r17="<<r17<<std::endl;
502 
503  Double32_t ptBeta17 = 0, pxBeta17 = 0, pyBeta17 = 0, pzBeta17 = 0;
504  Double32_t pabsBeta17 = (TMath::Sqrt((r17 + 511) * (r17 + 511) - (511 * 511))) * 1e-6;
505  // std::cout<<"pabsBeta17="<<pabsBeta17<<std::endl;
506 
507  Double32_t thetaBeta17 = acos(gRandom->Uniform(-1, 1));
508  Double32_t phiBeta17 = gRandom->Uniform(0, 360) * TMath::DegToRad();
509  pzBeta17 = pabsBeta17 * TMath::Cos(thetaBeta17);
510  ptBeta17 = pabsBeta17 * TMath::Sin(thetaBeta17);
511  pxBeta17 = ptBeta17 * TMath::Cos(phiBeta17);
512  pyBeta17 = ptBeta17 * TMath::Sin(phiBeta17);
513 
514  h1 = std::make_unique<TH1F>("h18", "h18", 1000, 0, 6000.);
515  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 3383.5);
516  f1->SetParameters(2.16857e-8, 1.29089e-9, -5.80327e-13, -7.02664e-17, -4.61969e-20, 3.46765e-23);
517  h1->FillRandom("f1");
518  h1->Draw();
519  Double32_t r18 = f1->GetRandom();
520  // std::cout<<"r18="<<r18<<std::endl;
521 
522  Double32_t ptBeta18 = 0, pxBeta18 = 0, pyBeta18 = 0, pzBeta18 = 0;
523  Double32_t pabsBeta18 = (TMath::Sqrt((r18 + 511) * (r18 + 511) - (511 * 511))) * 1e-6;
524  // std::cout<<"pabsBeta81="<<pabsBeta8<<std::endl;
525  Double32_t thetaBeta18 = acos(gRandom->Uniform(-1, 1));
526  Double32_t phiBeta18 = gRandom->Uniform(0, 360) * TMath::DegToRad();
527  pzBeta18 = pabsBeta18 * TMath::Cos(thetaBeta18);
528  ptBeta18 = pabsBeta18 * TMath::Sin(thetaBeta18);
529  pxBeta18 = ptBeta18 * TMath::Cos(phiBeta18);
530  pyBeta18 = ptBeta18 * TMath::Sin(phiBeta18);
531 
532  h1 = std::make_unique<TH1F>("h19", "h19", 1000, 0, 2997.5);
533  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 2997.5);
534  f1->SetParameters(5.06675e-8, 3.48106e-9, -2.25971e-12, 2.97006e-16, -3.48781e-19, 1.63705e-22);
535  h1->FillRandom("f1");
536  h1->Draw();
537  Double32_t r19 = f1->GetRandom();
538  // std::cout<<"r19="<<r19<<std::endl;
539 
540  Double32_t ptBeta19 = 0, pxBeta19 = 0, pyBeta19 = 0, pzBeta19 = 0;
541  Double32_t pabsBeta19 = (TMath::Sqrt((r19 + 511) * (r19 + 511) - (511 * 511))) * 1e-6;
542  // std::cout<<"pabsBeta19="<<pabsBeta19<<std::endl;
543 
544  Double32_t thetaBeta19 = acos(gRandom->Uniform(-1, 1));
545  Double32_t phiBeta19 = gRandom->Uniform(0, 360) * TMath::DegToRad();
546  pzBeta19 = pabsBeta19 * TMath::Cos(thetaBeta19);
547  ptBeta19 = pabsBeta9 * TMath::Sin(thetaBeta19);
548  pxBeta19 = ptBeta19 * TMath::Cos(phiBeta19);
549  pyBeta19 = ptBeta19 * TMath::Sin(phiBeta19);
550 
551  h1 = std::make_unique<TH1F>("h20", "h20", 1000, 0, 2596.5);
552  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 2596.5);
553  f1->SetParameters(3.50494e-8, 2.46542e-9, -1.68544e-12, 2.73518e-16, -2.77506e-19, 1.26138e-22);
554  h1->FillRandom("f1");
555  h1->Draw();
556  Double32_t r20 = f1->GetRandom();
557  // std::cout<<"r20="<<r20<<std::endl;
558 
559  Double32_t ptBeta20 = 0, pxBeta20 = 0, pyBeta20 = 0, pzBeta20 = 0;
560  Double32_t pabsBeta20 = (TMath::Sqrt((r20 + 511) * (r20 + 511) - (511 * 511))) * 1e-6;
561  // std::cout<<"pabsBeta20="<<pabsBeta20<<std::endl;
562  Double32_t thetaBeta20 = acos(gRandom->Uniform(-1, 1));
563  Double32_t phiBeta20 = gRandom->Uniform(0, 360) * TMath::DegToRad();
564  pzBeta20 = pabsBeta20 * TMath::Cos(thetaBeta20);
565  ptBeta20 = pabsBeta20 * TMath::Sin(thetaBeta20);
566  pxBeta20 = ptBeta20 * TMath::Cos(phiBeta20);
567  pyBeta20 = ptBeta20 * TMath::Sin(phiBeta20);
568 
569  h1 = std::make_unique<TH1F>("h21", "h21", 1000, 0, 2286.5);
570  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 2286.5);
571  f1->SetParameters(7.84682e-8, 7.34482e-9, -8.13364e-12, 3.91808e-15, -2.64262e-18, 9.8004e-22);
572  h1->FillRandom("f1");
573  h1->Draw();
574  Double32_t r21 = f1->GetRandom();
575  // std::cout<<"r21="<<r21<<std::endl;
576 
577  Double32_t ptBeta21 = 0, pxBeta21 = 0, pyBeta21 = 0, pzBeta21 = 0;
578  Double32_t pabsBeta21 = (TMath::Sqrt((r21 + 511) * (r21 + 511) - (511 * 511))) * 1e-6;
579  // std::cout<<"pabsBeta21="<<pabsBeta21<<std::endl;
580  Double32_t thetaBeta21 = acos(gRandom->Uniform(-1, 1));
581  Double32_t phiBeta21 = gRandom->Uniform(0, 360) * TMath::DegToRad();
582  pzBeta21 = pabsBeta21 * TMath::Cos(thetaBeta21);
583  ptBeta21 = pabsBeta21 * TMath::Sin(thetaBeta21);
584  pxBeta21 = ptBeta21 * TMath::Cos(phiBeta21);
585  pyBeta21 = ptBeta21 * TMath::Sin(phiBeta21);
586 
587  h1 = std::make_unique<TH1F>("h22", "h22", 1000, 0, 2027.5);
588  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 2027.5);
589  f1->SetParameters(1.1966e-8, 1.17991e-9, -1.40285e-12, 7.42593e-16, -4.98985e-19, 1.84071e-22);
590  h1->FillRandom("f1");
591  h1->Draw();
592  Double32_t r22 = f1->GetRandom();
593  // std::cout<<"r22="<<r22<<std::endl;
594 
595  Double32_t ptBeta22 = 0, pxBeta22 = 0, pyBeta22 = 0, pzBeta22 = 0;
596  Double32_t pabsBeta22 = (TMath::Sqrt((r22 + 511) * (r22 + 511) - (511 * 511))) * 1e-6;
597 
598  // std::cout<<"pabsBeta22="<<pabsBetav<<std::endl;
599 
600  Double32_t thetaBeta22 = acos(gRandom->Uniform(-1, 1));
601  Double32_t phiBeta22 = gRandom->Uniform(0, 360) * TMath::DegToRad();
602  pzBeta22 = pabsBeta22 * TMath::Cos(thetaBeta22);
603  ptBeta22 = pabsBeta22 * TMath::Sin(thetaBeta22);
604  pxBeta22 = ptBeta22 * TMath::Cos(phiBeta22);
605  pyBeta22 = ptBeta22 * TMath::Sin(phiBeta22);
606 
607  h1 = std::make_unique<TH1F>("h23", "h23", 1000, 0, 1986.5);
608  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 1986.5);
609  f1->SetParameters(8.82998e-10, 1.87941e-10, -4.67582e-13, 5.54255e-16, -4.80667e-19, 2.05681e-22);
610  h1->FillRandom("f1");
611  h1->Draw();
612  Double32_t r23 = f1->GetRandom();
613  // std::cout<<"r23="<<r13<<std::endl;
614 
615  Double32_t ptBeta23 = 0, pxBeta23 = 0, pyBeta23 = 0, pzBeta23 = 0;
616  Double32_t pabsBeta23 = (TMath::Sqrt((r23 + 511) * (r23 + 511) - (511 * 511))) * 1e-6;
617 
618  // std::cout<<"pabsBeta23="<<pabsBeta23<<std::endl;
619 
620  Double32_t thetaBeta23 = acos(gRandom->Uniform(-1, 1));
621  Double32_t phiBeta23 = gRandom->Uniform(0, 360) * TMath::DegToRad();
622  pzBeta23 = pabsBeta23 * TMath::Cos(thetaBeta23);
623  ptBeta23 = pabsBeta23 * TMath::Sin(thetaBeta23);
624  pxBeta23 = ptBeta23 * TMath::Cos(phiBeta23);
625  pyBeta23 = ptBeta23 * TMath::Sin(phiBeta23);
626 
627  h1 = std::make_unique<TH1F>("h24", "h24", 1000, 0, 1605);
628  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 1605);
629  f1->SetParameters(7.8520e-8, 7.34997e-9, -8.14208e-12, 3.92149e-15, -2.6442e-18, 9.80965e-22);
630  h1->FillRandom("f1");
631  h1->Draw();
632  Double32_t r24 = f1->GetRandom();
633  // std::cout<<"r24="<<r24<<std::endl;
634 
635  Double32_t ptBeta24 = 0, pxBeta24 = 0, pyBeta24 = 0, pzBeta24 = 0;
636  Double32_t pabsBeta24 = (TMath::Sqrt((r24 + 511) * (r24 + 511) - (511 * 511))) * 1e-6;
637  // std::cout<<"pabsBeta24="<<pabsBeta24<<std::endl;
638  Double32_t thetaBeta24 = acos(gRandom->Uniform(-1, 1));
639  Double32_t phiBeta24 = gRandom->Uniform(0, 360) * TMath::DegToRad();
640  pzBeta24 = pabsBeta24 * TMath::Cos(thetaBeta24);
641  ptBeta24 = pabsBeta24 * TMath::Sin(thetaBeta24);
642  pxBeta24 = ptBeta24 * TMath::Cos(phiBeta24);
643  pyBeta24 = ptBeta24 * TMath::Sin(phiBeta24);
644 
645  h1 = std::make_unique<TH1F>("h25", "h25", 1000, 0, 1550);
646  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 1550);
647  f1->SetParameters(1.43596e-8, 1.41593e-9, -1.68346e-12, 8.91141e-16, -5.98803e-19, 2.20892e-22);
648  h1->FillRandom("f1");
649  h1->Draw();
650  Double32_t r25 = f1->GetRandom();
651  // std::cout<<"r25="<<r25<<std::endl;
652 
653  Double32_t ptBeta25 = 0, pxBeta25 = 0, pyBeta25 = 0, pzBeta25 = 0;
654  Double32_t pabsBeta25 = (TMath::Sqrt((r25 + 511) * (r25 + 511) - (511 * 511))) * 1e-6;
655  // std::cout<<"pabsBeta25="<<pabsBeta25<<std::endl;
656 
657  Double32_t thetaBeta25 = acos(gRandom->Uniform(-1, 1));
658  Double32_t phiBeta25 = gRandom->Uniform(0, 360) * TMath::DegToRad();
659  pzBeta25 = pabsBeta25 * TMath::Cos(thetaBeta25);
660  ptBeta25 = pabsBeta25 * TMath::Sin(thetaBeta25);
661  pxBeta25 = ptBeta25 * TMath::Cos(phiBeta25);
662  pyBeta25 = ptBeta25 * TMath::Sin(phiBeta25);
663 
664  h1 = std::make_unique<TH1F>("h26", "h26", 900, 0, 984);
665  f1 = std::make_unique<TF1>("f1", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x", 0, 984);
666  f1->SetParameters(8.82998e-10, 1.87941e-10, -4.67582e-13, 5.54255e-16, -4.80667e-19, 2.05681e-22);
667  h1->FillRandom("f1");
668  h1->Draw();
669  Double32_t r26 = f1->GetRandom();
670  // std::cout<<"r26="<<r26<<std::endl;
671  //
672  Double32_t ptBeta26 = 0, pxBeta26 = 0, pyBeta26 = 0, pzBeta26 = 0;
673  Double32_t pabsBeta26 = (TMath::Sqrt((r26 + 511) * (r26 + 511) - (511 * 511))) * 1e-6;
674  //::cout<<"pabsBeta26="<<pabsBeta26<<std::endl;
675 
676  Double32_t thetaBeta26 = acos(gRandom->Uniform(-1, 1));
677  Double32_t phiBeta26 = gRandom->Uniform(0, 360) * TMath::DegToRad();
678  pzBeta26 = pabsBeta26 * TMath::Cos(thetaBeta26);
679  ptBeta26 = pabsBeta26 * TMath::Sin(thetaBeta26);
680  pxBeta26 = ptBeta26 * TMath::Cos(phiBeta26);
681  pyBeta26 = ptBeta26 * TMath::Sin(phiBeta26);
682 
683  // Alpha particles from decay of 20Ne to 16O
684 
685  Double32_t ptAlpha1 = 0, pxAlpha1 = 0, pyAlpha1 = 0, pzAlpha1 = 0;
686  Double32_t pabsAlpha1 = 0.2328; // GeV/c, 7.26 MeV
687  Double32_t thetaAlpha1 = acos(gRandom->Uniform(-1, 1));
688  Double32_t phiAlpha1 = gRandom->Uniform(0, 360) * TMath::DegToRad();
689  pzAlpha1 = pabsAlpha1 * TMath::Cos(thetaAlpha1);
690  ptAlpha1 = pabsAlpha1 * TMath::Sin(thetaAlpha1);
691  pxAlpha1 = ptAlpha1 * TMath::Cos(phiAlpha1);
692  pyAlpha1 = ptAlpha1 * TMath::Sin(phiAlpha1);
693 
694  Double32_t ptAlpha2 = 0, pxAlpha2 = 0, pyAlpha2 = 0, pzAlpha2 = 0;
695  Double32_t pabsAlpha2 = 0.2204; // GeV/c, 6.561 MeV
696  Double32_t thetaAlpha2 = acos(gRandom->Uniform(-1, 1));
697  Double32_t phiAlpha2 = gRandom->Uniform(0, 360) * TMath::DegToRad();
698  pzAlpha2 = pabsAlpha2 * TMath::Cos(thetaAlpha2);
699  ptAlpha2 = pabsAlpha2 * TMath::Sin(thetaAlpha2);
700  pxAlpha2 = ptAlpha2 * TMath::Cos(phiAlpha2);
701  pyAlpha2 = ptAlpha2 * TMath::Sin(phiAlpha2);
702 
703  Double32_t ptAlpha3 = 0, pxAlpha3 = 0, pyAlpha3 = 0, pzAlpha3 = 0;
704  Double32_t pabsAlpha3 = 0.2134; // GeV/c, 6.106 Mev
705  Double32_t thetaAlpha3 = acos(gRandom->Uniform(-1, 1));
706  Double32_t phiAlpha3 = gRandom->Uniform(0, 360) * TMath::DegToRad();
707  pzAlpha3 = pabsAlpha3 * TMath::Cos(thetaAlpha3);
708  ptAlpha3 = pabsAlpha3 * TMath::Sin(thetaAlpha3);
709  pxAlpha3 = ptAlpha3 * TMath::Cos(phiAlpha3);
710  pyAlpha3 = ptAlpha3 * TMath::Sin(phiAlpha3);
711 
712  Double32_t ptAlpha4 = 0, pxAlpha4 = 0, pyAlpha4 = 0, pzAlpha4 = 0;
713  Double32_t pabsAlpha4 = 0.2088; // GeV/c, 5.844 MeV
714  Double32_t thetaAlpha4 = acos(gRandom->Uniform(-1, 1));
715  Double32_t phiAlpha4 = gRandom->Uniform(0, 360) * TMath::DegToRad();
716  pzAlpha4 = pabsAlpha4 * TMath::Cos(thetaAlpha4);
717  ptAlpha4 = pabsAlpha4 * TMath::Sin(thetaAlpha4);
718  pxAlpha4 = ptAlpha4 * TMath::Cos(phiAlpha4);
719  pyAlpha4 = ptAlpha4 * TMath::Sin(phiAlpha4);
720 
721  Double32_t ptAlpha5 = 0, pxAlpha5 = 0, pyAlpha5 = 0, pzAlpha5 = 0;
722  Double32_t pabsAlpha5 = 0.2033; // GeV/c, 5.540
723  Double32_t thetaAlpha5 = acos(gRandom->Uniform(-1, 1));
724  Double32_t phiAlpha5 = gRandom->Uniform(0, 360) * TMath::DegToRad();
725  pzAlpha5 = pabsAlpha5 * TMath::Cos(thetaAlpha5);
726  ptAlpha5 = pabsAlpha5 * TMath::Sin(thetaAlpha5);
727  pxAlpha5 = ptAlpha5 * TMath::Cos(phiAlpha5);
728  pyAlpha5 = ptAlpha5 * TMath::Sin(phiAlpha5);
729 
730  Double32_t ptAlpha6 = 0, pxAlpha6 = 0, pyAlpha6 = 0, pzAlpha6 = 0;
731  Double32_t pabsAlpha6 = 0.1882; // GeV/c, 4.749 MeV
732  Double32_t thetaAlpha6 = acos(gRandom->Uniform(-1, 1));
733  Double32_t phiAlpha6 = gRandom->Uniform(0, 360) * TMath::DegToRad();
734  pzAlpha6 = pabsAlpha6 * TMath::Cos(thetaAlpha6);
735  ptAlpha6 = pabsAlpha6 * TMath::Sin(thetaAlpha6);
736  pxAlpha6 = ptAlpha6 * TMath::Cos(phiAlpha6);
737  pyAlpha6 = ptAlpha6 * TMath::Sin(phiAlpha6);
738 
739  Double32_t ptAlpha7 = 0, pxAlpha7 = 0, pyAlpha7 = 0, pzAlpha7 = 0;
740  Double32_t pabsAlpha7 = 0.1757; // GeV/c, 4.140 MeV
741  Double32_t thetaAlpha7 = acos(gRandom->Uniform(-1, 1));
742  Double32_t phiAlpha7 = gRandom->Uniform(0, 360) * TMath::DegToRad();
743  pzAlpha7 = pabsAlpha7 * TMath::Cos(thetaAlpha7);
744  ptAlpha7 = pabsAlpha7 * TMath::Sin(thetaAlpha7);
745  pxAlpha7 = ptAlpha7 * TMath::Cos(phiAlpha7);
746  pyAlpha7 = ptAlpha7 * TMath::Sin(phiAlpha7);
747 
748  Double32_t ptAlpha8 = 0, pxAlpha8 = 0, pyAlpha8 = 0, pzAlpha8 = 0;
749  Double32_t pabsAlpha8 = 0.152; // GeV/c, 3.099 MeV
750  Double32_t thetaAlpha8 = acos(gRandom->Uniform(-1, 1));
751  Double32_t phiAlpha8 = gRandom->Uniform(0, 360) * TMath::DegToRad();
752  pzAlpha8 = pabsAlpha8 * TMath::Cos(thetaAlpha8);
753  ptAlpha8 = pabsAlpha8 * TMath::Sin(thetaAlpha8);
754  pxAlpha8 = ptAlpha8 * TMath::Cos(phiAlpha8);
755  pyAlpha8 = ptAlpha8 * TMath::Sin(phiAlpha8);
756 
757  Double32_t ptAlpha9 = 0, pxAlpha9 = 0, pyAlpha9 = 0, pzAlpha9 = 0;
758  Double32_t pabsAlpha9 = 0.1417; // GeV/c, 2.692 MeV
759  Double32_t thetaAlpha9 = acos(gRandom->Uniform(-1, 1));
760  Double32_t phiAlpha9 = gRandom->Uniform(0, 360) * TMath::DegToRad();
761  pzAlpha9 = pabsAlpha9 * TMath::Cos(thetaAlpha9);
762  ptAlpha9 = pabsAlpha9 * TMath::Sin(thetaAlpha9);
763  pxAlpha9 = ptAlpha9 * TMath::Cos(phiAlpha9);
764  pyAlpha9 = ptAlpha9 * TMath::Sin(phiAlpha9);
765 
766  Double32_t ran = 0;
767  Double32_t bra1 = 0;
768  Double32_t bra = 0;
769  if (fNuclearDecayChainIsSet) {
770 
771  if (!(protonPDGID == 2212))
772  LOG(fatal) << "AtTPC20MgDecayGenerator:PDG code" << protonPDGID << "is not a proton!";
773  brp = gRandom->Uniform(0, 1); // uniform random number for proton
774  bra = gRandom->Uniform(0, 1); // uniform random number for alpha (20Ne decay to 16O)
775  brbeta = gRandom->Uniform(0, 1); // uniform random number for beta (20Mg decay)
776  brb = gRandom->Uniform(0, 1); // uniform random number for beta (20Mgdecay to 20Na followed by 19Ne)
777  ran = gRandom->Uniform(0, 1); // uniform random number for beta(20Na decay to 20Ne)
778  bra1 = gRandom->Uniform(0, 1); // uniform random number for alpha (19Ne decay to 15O)
779 
780  for (Int_t i = 0; i < fParticlesDefinedInNuclearDecay; i++) {
781 
782  if ((0 < brbeta) && (brbeta <= 0.72)) {
783  std::cout << "proton branch along with alpha decay to 15O" << std::endl;
784  if ((0 < brb) && (brb <= 0.7142)) {
785  Double32_t BetaMomentum1 = TMath::Sqrt(pxBeta1 * pxBeta1 + pyBeta1 * pyBeta1 + pzBeta1 * pzBeta1);
786  pxBeta1 = pxBeta1 * r1 * 1e-6 * fParticleEnergies[i] / BetaMomentum1 * fParticleEnergies[i];
787  pyBeta1 = pyBeta1 * fParticleEnergies[i] * r1 * 1e-6 / BetaMomentum1 * fParticleEnergies[i];
788  pzBeta1 = pzBeta1 * r1 * 1e-6 * fParticleEnergies[i] / BetaMomentum1 * fParticleEnergies[i];
789  }
790 
791  else if ((0.7142 < brb) && (brb <= 0.71523)) {
792  Double32_t BetaMomentum2 = TMath::Sqrt(pxBeta2 * pxBeta2 + pyBeta2 * pyBeta2 + pzBeta2 * pzBeta2);
793  pxBeta2 = pxBeta2 * r2 * 1e-6 * fParticleEnergies[i] / BetaMomentum2 * fParticleEnergies[i];
794  pyBeta2 = pyBeta2 * r2 * 1e-6 * fParticleEnergies[i] / BetaMomentum2 * fParticleEnergies[i];
795  pzBeta2 = pzBeta2 * r2 * 1e-6 * fParticleEnergies[i] / BetaMomentum2 * fParticleEnergies[i];
796  }
797 
798  else if ((0.999733 < brb) && (brb <= 1)) {
799  {
800  Double32_t BetaMomentum13 =
801  TMath::Sqrt(pxBeta13 * pxBeta13 + pyBeta13 * pyBeta13 + pzBeta13 * pzBeta13);
802  pxBeta13 = pxBeta13 * r13 * 1e-6 * fParticleEnergies[i] / BetaMomentum13 * fParticleEnergies[i];
803  pyBeta13 = pyBeta13 * r13 * 1e-6 * fParticleEnergies[i] / BetaMomentum13 * fParticleEnergies[i];
804  pzBeta13 = pzBeta13 * r13 * 1e-6 * fParticleEnergies[i] / BetaMomentum13 * fParticleEnergies[i];
805  }
806  if (brp == 1) {
807  {
808  Double32_t ProtonMomentum =
809  TMath::Sqrt(pxProton * pxProton + pyProton * pyProton + pzProton * pzProton);
810  pxProton = pxProton * fParticleEnergies[i] / ProtonMomentum;
811  pyProton = pyProton * fParticleEnergies[i] / ProtonMomentum;
812  pzProton = pzProton * fParticleEnergies[i] / ProtonMomentum;
813  }
814 
815  if (bra1 == 1) {
816  Double32_t AlphaMomentum = TMath::Sqrt(pxAlpha * pxAlpha + pyAlpha * pyAlpha + pzAlpha * pzAlpha);
817 
818  pxAlpha = pxAlpha * fParticleEnergies[i] / AlphaMomentum;
819  pyAlpha = pyAlpha * fParticleEnergies[i] / AlphaMomentum;
820  pzAlpha = pzAlpha * fParticleEnergies[i] / AlphaMomentum;
821  }
822  }
823  }
824 
825  else if ((0.71523 < brb) && (brb <= 0.83426)) {
826  {
827  Double32_t BetaMomentum3 = TMath::Sqrt(pxBeta3 * pxBeta3 + pyBeta3 * pyBeta3 + pzBeta3 * pzBeta3);
828  pxBeta3 = pxBeta3 * r3 * 1e-6 * fParticleEnergies[i] / BetaMomentum3 * fParticleEnergies[i];
829  pyBeta3 = pyBeta3 * r3 * 1e-6 * fParticleEnergies[i] / BetaMomentum3 * fParticleEnergies[i];
830  pzBeta3 = pzBeta3 * r3 * 1e-6 * fParticleEnergies[i] / BetaMomentum3 * fParticleEnergies[i];
831  }
832  if (brp == 1) {
833 
834  Double32_t ProtonMomentum1 =
835  TMath::Sqrt(pxProton1 * pxProton1 + pyProton1 * pyProton1 + pzProton1 * pzProton1);
836  pxProton1 = pxProton1 * fParticleEnergies[i] / ProtonMomentum1;
837  pyProton1 = pyProton1 * fParticleEnergies[i] / ProtonMomentum1;
838  pzProton1 = pzProton1 * fParticleEnergies[i] / ProtonMomentum1;
839  }
840  }
841 
842  else if ((0.9118 < brb) && (brb <= 0.93146)) {
843  {
844  Double32_t BetaMomentum6 = TMath::Sqrt(pxBeta6 * pxBeta6 + pyBeta6 * pyBeta6 + pzBeta6 * pzBeta6);
845  pxBeta6 = pxBeta6 * r6 * 1e-6 * fParticleEnergies[i] / BetaMomentum6 * fParticleEnergies[i];
846  pyBeta6 = pyBeta6 * r6 * 1e-6 * fParticleEnergies[i] / BetaMomentum6 * fParticleEnergies[i];
847  pzBeta6 = pzBeta6 * r6 * 1e-6 * fParticleEnergies[i] / BetaMomentum6 * fParticleEnergies[i];
848  }
849  if ((0 < brp) && (brp <= 0.4375)) {
850  Double32_t ProtonMomentum2 =
851  TMath::Sqrt(pxProton2 * pxProton2 + pyProton2 * pyProton2 + pzProton2 * pzProton2);
852  pxProton2 = pxProton2 * fParticleEnergies[i] / ProtonMomentum2;
853  pyProton2 = pyProton2 * fParticleEnergies[i] / ProtonMomentum2;
854  pzProton2 = pzProton2 * fParticleEnergies[i] / ProtonMomentum2;
855  } else if ((0.4375 < brp) && (brp <= 0.6250)) {
856  Double32_t ProtonMomentum6 =
857  TMath::Sqrt(pxProton6 * pxProton6 + pyProton6 * pyProton6 + pzProton6 * pzProton6);
858  pxProton6 = pxProton6 * fParticleEnergies[i] / ProtonMomentum6;
859  pyProton6 = pyProton6 * fParticleEnergies[i] / ProtonMomentum6;
860  pzProton6 = pzProton6 * fParticleEnergies[i] / ProtonMomentum6;
861  } else if ((0.625 < brp) && (brp <= 0.875)) {
862  Double32_t ProtonMomentum7 =
863  TMath::Sqrt(pxProton7 * pxProton7 + pyProton7 * pyProton7 + pzProton7 * pzProton7);
864  pxProton7 = pxProton7 * fParticleEnergies[i] / ProtonMomentum7;
865  pyProton7 = pyProton7 * fParticleEnergies[i] / ProtonMomentum7;
866  pzProton7 = pzProton7 * fParticleEnergies[i] / ProtonMomentum7;
867  } else if ((0.875 < brp) && (brp <= 1)) {
868  Double32_t ProtonMomentum8 =
869  TMath::Sqrt(pxProton8 * pxProton8 + pyProton8 * pyProton8 + pzProton8 * pzProton8);
870  pxProton8 = pxProton8 * fParticleEnergies[i] / ProtonMomentum8;
871  pyProton8 = pyProton8 * fParticleEnergies[i] / ProtonMomentum8;
872  pzProton8 = pzProton8 * fParticleEnergies[i] / ProtonMomentum8;
873  }
874  }
875 
876  else if ((0.83426 < brb) && (brb <= 0.88386)) {
877  {
878  Double32_t BetaMomentum4 = TMath::Sqrt(pxBeta4 * pxBeta4 + pyBeta4 * pyBeta4 + pzBeta4 * pzBeta4);
879  pxBeta4 = pxBeta4 * r4 * 1e-6 * fParticleEnergies[i] / BetaMomentum4 * fParticleEnergies[i];
880  pyBeta4 = pyBeta4 * r4 * 1e-6 * fParticleEnergies[i] / BetaMomentum4 * fParticleEnergies[i];
881  pzBeta4 = pzBeta4 * r4 * 1e-6 * fParticleEnergies[i] / BetaMomentum4 * fParticleEnergies[i];
882  }
883  if ((0 < brp) && (brp <= 0.0769)) {
884  Double32_t ProtonMomentum3 =
885  TMath::Sqrt(pxProton3 * pxProton3 + pyProton3 * pyProton3 + pzProton3 * pzProton3);
886  pxProton3 = pxProton3 * fParticleEnergies[i] / ProtonMomentum3;
887  pyProton3 = pyProton3 * fParticleEnergies[i] / ProtonMomentum3;
888  pzProton3 = pzProton3 * fParticleEnergies[i] / ProtonMomentum3;
889  }
890 
891  else if ((0.0769 < brp) && (brp <= 1)) {
892  Double32_t ProtonMomentum4 =
893  TMath::Sqrt(pxProton4 * pxProton4 + pyProton4 * pyProton4 + pzProton4 * pzProton4);
894  pxProton4 = pxProton4 * fParticleEnergies[i] / ProtonMomentum4;
895  pyProton4 = pyProton4 * fParticleEnergies[i] / ProtonMomentum4;
896  pzProton4 = pzProton4 * fParticleEnergies[i] / ProtonMomentum4;
897  }
898  }
899 
900  else if ((0.88386 < brb) && (brb <= 0.9118)) {
901  {
902  Double32_t BetaMomentum5 = TMath::Sqrt(pxBeta5 * pxBeta5 + pyBeta5 * pyBeta5 + pzBeta5 * pzBeta5);
903  pxBeta5 = pxBeta5 * r5 * 1e-6 * fParticleEnergies[i] / BetaMomentum5 * fParticleEnergies[i];
904  pyBeta5 = pyBeta5 * r5 * 1e-6 * fParticleEnergies[i] / BetaMomentum5 * fParticleEnergies[i];
905  pzBeta5 = pzBeta5 * r5 * 1e-6 * fParticleEnergies[i] / BetaMomentum5 * fParticleEnergies[i];
906  }
907  if (brp == 1) {
908  Double32_t ProtonMomentum5 =
909  TMath::Sqrt(pxProton5 * pxProton5 + pyProton5 * pyProton5 + pzProton5 * pzProton5);
910  pxProton5 = pxProton5 * fParticleEnergies[i] / ProtonMomentum5;
911  pyProton5 = pyProton5 * fParticleEnergies[i] / ProtonMomentum5;
912  pzProton5 = pzProton5 * fParticleEnergies[i] / ProtonMomentum5;
913  }
914  }
915 
916  else if ((0.93146 < brb) && (brb <= 0.94698)) {
917  {
918  Double32_t BetaMomentum7 = TMath::Sqrt(pxBeta7 * pxBeta7 + pyBeta7 * pyBeta7 + pzBeta7 * pzBeta7);
919  pxBeta7 = pxBeta7 * r7 * 1e-6 * fParticleEnergies[i] / BetaMomentum7 * fParticleEnergies[i];
920  pyBeta7 = pyBeta7 * r7 * 1e-6 * fParticleEnergies[i] / BetaMomentum7 * fParticleEnergies[i];
921  pzBeta7 = pzBeta7 * r7 * 1e-6 * fParticleEnergies[i] / BetaMomentum7 * fParticleEnergies[i];
922  }
923  if (brp == 1) {
924  Double32_t ProtonMomentum9 =
925  TMath::Sqrt(pxProton9 * pxProton9 + pyProton9 * pyProton9 + pzProton9 * pzProton9);
926  pxProton9 = pxProton9 * fParticleEnergies[i] / ProtonMomentum9;
927  pyProton9 = pyProton9 * fParticleEnergies[i] / ProtonMomentum9;
928  pzProton9 = pzProton9 * fParticleEnergies[i] / ProtonMomentum9;
929  }
930  }
931 
932  else if ((0.94698 < brb) && (brb <= 0.95277)) {
933  {
934  Double32_t BetaMomentum8 = TMath::Sqrt(pxBeta8 * pxBeta8 + pyBeta8 * pyBeta8 + pzBeta8 * pzBeta8);
935  pxBeta8 = pxBeta8 * r8 * 1e-6 * fParticleEnergies[i] / BetaMomentum8 * fParticleEnergies[i];
936  pyBeta8 = pyBeta8 * r8 * 1e-6 * fParticleEnergies[i] / BetaMomentum8 * fParticleEnergies[i];
937  pzBeta8 = pzBeta8 * r8 * 1e-6 * fParticleEnergies[i] / BetaMomentum8 * fParticleEnergies[i];
938  }
939  if (brp == 1) {
940  Double32_t ProtonMomentum10 =
941  TMath::Sqrt(pxProton10 * pxProton10 + pyProton10 * pyProton10 + pzProton10 * pzProton10);
942  pxProton10 = pxProton10 * fParticleEnergies[i] / ProtonMomentum10;
943  pyProton10 = pyProton10 * fParticleEnergies[i] / ProtonMomentum10;
944  pzProton10 = pzProton10 * fParticleEnergies[i] / ProtonMomentum10;
945  }
946  }
947 
948  else if ((0.95277 < brb) && (brb <= 0.96517)) {
949  {
950  Double32_t BetaMomentum9 = TMath::Sqrt(pxBeta9 * pxBeta9 + pyBeta9 * pyBeta9 + pzBeta9 * pzBeta9);
951  pxBeta9 = pxBeta9 * r9 * 1e-6 * fParticleEnergies[i] / BetaMomentum9 * fParticleEnergies[i];
952  pyBeta9 = pyBeta9 * r9 * 1e-6 * fParticleEnergies[i] / BetaMomentum9 * fParticleEnergies[i];
953  pzBeta9 = pzBeta9 * r9 * 1e-6 * fParticleEnergies[i] / BetaMomentum9 * fParticleEnergies[i];
954  }
955  if (brp == 1) {
956  Double32_t ProtonMomentum11 =
957  TMath::Sqrt(pxProton11 * pxProton11 + pyProton11 * pyProton11 + pzProton11 * pzProton11);
958  pxProton11 = pxProton11 * fParticleEnergies[i] / ProtonMomentum11;
959  pyProton11 = pyProton11 * fParticleEnergies[i] / ProtonMomentum11;
960  pzProton11 = pzProton11 * fParticleEnergies[i] / ProtonMomentum11;
961  }
962  }
963 
964  else if ((0.96517 < brb) && (brb <= 0.99932)) {
965  {
966  Double32_t BetaMomentum10 =
967  TMath::Sqrt(pxBeta10 * pxBeta10 + pyBeta10 * pyBeta10 + pzBeta10 * pzBeta10);
968  pxBeta10 = pxBeta10 * r10 * 1e-6 * fParticleEnergies[i] / BetaMomentum10 * fParticleEnergies[i];
969  pyBeta10 = pyBeta10 * r10 * 1e-6 * fParticleEnergies[i] / BetaMomentum10 * fParticleEnergies[i];
970  pzBeta10 = pzBeta10 * r10 * 1e-6 * fParticleEnergies[i] / BetaMomentum10 * fParticleEnergies[i];
971  }
972  if ((0 < brp) && (brp <= 0.333)) {
973  Double32_t ProtonMomentum12 =
974  TMath::Sqrt(pxProton12 * pxProton12 + pyProton12 * pyProton12 + pzProton12 * pzProton12);
975  pxProton12 = pxProton12 * fParticleEnergies[i] / ProtonMomentum12;
976  pyProton12 = pyProton12 * fParticleEnergies[i] / ProtonMomentum12;
977  pzProton12 = pzProton12 * fParticleEnergies[i] / ProtonMomentum12;
978  }
979 
980  else if ((0.333 < brp) && (brp <= 1)) {
981  Double32_t ProtonMomentum13 =
982  TMath::Sqrt(pxProton13 * pxProton13 + pyProton13 * pyProton13 + pzProton13 * pzProton13);
983  pxProton13 = pxProton13 * fParticleEnergies[i] / ProtonMomentum13;
984  pyProton13 = pyProton13 * fParticleEnergies[i] / ProtonMomentum13;
985  pzProton13 = pzProton13 * fParticleEnergies[i] / ProtonMomentum13;
986  }
987  }
988 
989  else if ((0.99932 < brb) && (brb <= 0.99963)) {
990  Double32_t BetaMomentum11 = TMath::Sqrt(pxBeta11 * pxBeta11 + pyBeta11 * pyBeta11 + pzBeta11 * pzBeta11);
991  pxBeta11 = pxBeta11 * r11 * 1e-6 * fParticleEnergies[i] / BetaMomentum11 * fParticleEnergies[i];
992  pyBeta11 = pyBeta11 * r11 * 1e-6 * fParticleEnergies[i] / BetaMomentum11 * fParticleEnergies[i];
993  pzBeta11 = pzBeta11 * r11 * 1e-6 * fParticleEnergies[i] / BetaMomentum11 * fParticleEnergies[i];
994  }
995 
996  else if ((0.99963 < brb) && (brb <= 0.999733)) {
997  Double32_t BetaMomentum12 = TMath::Sqrt(pxBeta12 * pxBeta12 + pyBeta12 * pyBeta12 + pzBeta12 * pzBeta12);
998  pxBeta12 = pxBeta12 * r12 * 1e-6 * fParticleEnergies[i] / BetaMomentum12 * fParticleEnergies[i];
999  pyBeta12 = pyBeta12 * r12 * 1e-6 * fParticleEnergies[i] / BetaMomentum12 * fParticleEnergies[i];
1000  pzBeta12 = pzBeta12 * r12 * 1e-6 * fParticleEnergies[i] / BetaMomentum12 * fParticleEnergies[i];
1001  }
1002  } else if ((0.72 < brbeta) && (brbeta <= 1)) {
1003  // std::cout<<"beta branch starts"<<std::endl;
1004  if (ran == 1) {
1005  {
1006  Double32_t BetaMomentum14 =
1007  TMath::Sqrt(pxBeta14 * pxBeta14 + pyBeta14 * pyBeta14 + pzBeta14 * pzBeta14);
1008  pxBeta14 = pxBeta14 * r14 * 1e-6 * fParticleEnergies[i] / BetaMomentum14 * fParticleEnergies[i];
1009  pyBeta14 = pyBeta14 * r14 * 1e-6 * fParticleEnergies[i] / BetaMomentum14 * fParticleEnergies[i];
1010  pzBeta14 = pzBeta14 * r14 * 1e-6 * fParticleEnergies[i] / BetaMomentum14 * fParticleEnergies[i];
1011  }
1012  if (bra == 1) {
1013  Double32_t AlphaMomentum1 =
1014  TMath::Sqrt(pxAlpha1 * pxAlpha1 + pyAlpha1 * pyAlpha1 + pzAlpha1 * pzAlpha1);
1015 
1016  pxAlpha1 = pxAlpha1 * fParticleEnergies[i] / AlphaMomentum1;
1017  pyAlpha1 = pyAlpha1 * fParticleEnergies[i] / AlphaMomentum1;
1018  pzAlpha1 = pzAlpha1 * fParticleEnergies[i] / AlphaMomentum1;
1019  }
1020  }
1021 
1022  else if ((0.000016 < ran) && (ran <= 0.002646)) {
1023  {
1024  Double32_t BetaMomentum15 =
1025  TMath::Sqrt(pxBeta15 * pxBeta15 + pyBeta15 * pyBeta15 + pzBeta15 * pzBeta15);
1026  pxBeta15 = pxBeta15 * r15 * 1e-6 * fParticleEnergies[i] / BetaMomentum15 * fParticleEnergies[i];
1027  pyBeta15 = pyBeta15 * r15 * 1e-6 * fParticleEnergies[i] / BetaMomentum15 * fParticleEnergies[i];
1028  pzBeta15 = pzBeta15 * r15 * 1e-6 * fParticleEnergies[i] / BetaMomentum15 * fParticleEnergies[i];
1029  }
1030  if (bra == 1) {
1031  Double32_t AlphaMomentum2 =
1032  TMath::Sqrt(pxAlpha2 * pxAlpha2 + pyAlpha2 * pyAlpha2 + pzAlpha2 * pzAlpha2);
1033 
1034  pxAlpha2 = pxAlpha2 * fParticleEnergies[i] / AlphaMomentum2;
1035  pyAlpha2 = pyAlpha2 * fParticleEnergies[i] / AlphaMomentum2;
1036  pzAlpha2 = pzAlpha2 * fParticleEnergies[i] / AlphaMomentum2;
1037  }
1038  }
1039 
1040  else if ((0.002646 < ran) && (ran <= 0.004696)) {
1041  Double32_t BetaMomentum16 = TMath::Sqrt(pxBeta16 * pxBeta16 + pyBeta16 * pyBeta16 + pzBeta16 * pzBeta16);
1042  pxBeta16 = pxBeta16 * r16 * 1e-6 * fParticleEnergies[i] / BetaMomentum16 * fParticleEnergies[i];
1043  pyBeta16 = pyBeta16 * r16 * 1e-6 * fParticleEnergies[i] / BetaMomentum16 * fParticleEnergies[i];
1044  pzBeta16 = pzBeta16 * r16 * 1e-6 * fParticleEnergies[i] / BetaMomentum16 * fParticleEnergies[i];
1045  }
1046 
1047  else if ((0.004696 < ran) && (ran <= 0.005866)) {
1048  Double32_t BetaMomentum17 = TMath::Sqrt(pxBeta17 * pxBeta17 + pyBeta17 * pyBeta17 + pzBeta17 * pzBeta17);
1049  pxBeta17 = pxBeta1 * r17 * 1e-6 * fParticleEnergies[i] / BetaMomentum17 * fParticleEnergies[i];
1050  pyBeta17 = pyBeta1 * r17 * 1e-6 * fParticleEnergies[i] / BetaMomentum17 * fParticleEnergies[i];
1051  pzBeta17 = pzBeta1 * r17 * 1e-6 * fParticleEnergies[i] / BetaMomentum17 * fParticleEnergies[i];
1052  }
1053 
1054  else if ((0.005866 < ran) && (ran <= 0.007606)) {
1055  {
1056  Double32_t BetaMomentum18 =
1057  TMath::Sqrt(pxBeta18 * pxBeta18 + pyBeta18 * pyBeta18 + pzBeta18 * pzBeta18);
1058  pxBeta18 = pxBeta18 * r18 * 1e-6 * fParticleEnergies[i] / BetaMomentum18 * fParticleEnergies[i];
1059  pyBeta18 = pyBeta18 * r18 * 1e-6 * fParticleEnergies[i] / BetaMomentum18 * fParticleEnergies[i];
1060  pzBeta18 = pzBeta18 * r18 * 1e-6 * fParticleEnergies[i] / BetaMomentum18 * fParticleEnergies[i];
1061  }
1062  if (bra == 1) {
1063  Double32_t AlphaMomentum3 =
1064  TMath::Sqrt(pxAlpha3 * pxAlpha3 + pyAlpha3 * pyAlpha3 + pzAlpha3 * pzAlpha3);
1065 
1066  pxAlpha3 = pxAlpha3 * fParticleEnergies[i] / AlphaMomentum3;
1067  pyAlpha3 = pyAlpha3 * fParticleEnergies[i] / AlphaMomentum3;
1068  pzAlpha3 = pzAlpha3 * fParticleEnergies[i] / AlphaMomentum3;
1069  }
1070  }
1071 
1072  else if ((0.007606 < ran) && (ran <= 0.008489)) {
1073  {
1074  Double32_t BetaMomentum19 =
1075  TMath::Sqrt(pxBeta19 * pxBeta19 + pyBeta19 * pyBeta19 + pzBeta19 * pzBeta19);
1076  pxBeta19 = pxBeta19 * r19 * 1e-6 * fParticleEnergies[i] / BetaMomentum19 * fParticleEnergies[i];
1077  pyBeta19 = pyBeta19 * r19 * 1e-6 * fParticleEnergies[i] / BetaMomentum19 * fParticleEnergies[i];
1078  pzBeta19 = pzBeta19 * r19 * 1e-6 * fParticleEnergies[i] / BetaMomentum19 * fParticleEnergies[i];
1079  }
1080  if (bra == 1) {
1081  Double32_t AlphaMomentum4 =
1082  TMath::Sqrt(pxAlpha4 * pxAlpha4 + pyAlpha4 * pyAlpha4 + pzAlpha4 * pzAlpha4);
1083 
1084  pxAlpha4 = pxAlpha4 * fParticleEnergies[i] / AlphaMomentum4;
1085  pyAlpha4 = pyAlpha4 * fParticleEnergies[i] / AlphaMomentum4;
1086  pzAlpha4 = pzAlpha4 * fParticleEnergies[i] / AlphaMomentum4;
1087  }
1088  }
1089 
1090  else if ((0.008489 < ran) && (ran <= 0.037259)) {
1091  {
1092  Double32_t BetaMomentum20 =
1093  TMath::Sqrt(pxBeta20 * pxBeta20 + pyBeta20 * pyBeta20 + pzBeta20 * pzBeta20);
1094  pxBeta20 = pxBeta20 * r20 * 1e-6 * fParticleEnergies[i] / BetaMomentum20 * fParticleEnergies[i];
1095  pyBeta20 = pyBeta20 * r20 * 1e-6 * fParticleEnergies[i] / BetaMomentum20 * fParticleEnergies[i];
1096  pzBeta20 = pzBeta20 * r20 * 1e-6 * fParticleEnergies[i] / BetaMomentum20 * fParticleEnergies[i];
1097  }
1098  if (bra == 1) {
1099  Double32_t AlphaMomentum5 =
1100  TMath::Sqrt(pxAlpha5 * pxAlpha5 + pyAlpha5 * pyAlpha5 + pzAlpha5 * pzAlpha5);
1101 
1102  pxAlpha5 = pxAlpha5 * fParticleEnergies[i] / AlphaMomentum5;
1103  pyAlpha5 = pyAlpha5 * fParticleEnergies[i] / AlphaMomentum5;
1104  pzAlpha5 = pzAlpha5 * fParticleEnergies[i] / AlphaMomentum5;
1105  }
1106  }
1107 
1108  else if ((0.037259 < ran) && (ran <= 0.037539)) {
1109  Double32_t BetaMomentum21 = TMath::Sqrt(pxBeta21 * pxBeta21 + pyBeta21 * pyBeta21 + pzBeta21 * pzBeta21);
1110  pxBeta21 = pxBeta21 * r21 * 1e-6 * fParticleEnergies[i] / BetaMomentum21 * fParticleEnergies[i];
1111  pyBeta21 = pyBeta21 * r21 * 1e-6 * fParticleEnergies[i] / BetaMomentum21 * fParticleEnergies[i];
1112  pzBeta21 = pzBeta21 * r21 * 1e-6 * fParticleEnergies[i] / BetaMomentum21 * fParticleEnergies[i];
1113  }
1114 
1115  else if ((0.037539 < ran) && (ran <= 0.039949)) {
1116  {
1117  Double32_t BetaMomentum22 =
1118  TMath::Sqrt(pxBeta22 * pxBeta22 + pyBeta22 * pyBeta22 + pzBeta22 * pzBeta22);
1119  pxBeta22 = pxBeta22 * r22 * 1e-6 * fParticleEnergies[i] / BetaMomentum22 * fParticleEnergies[i];
1120  pyBeta22 = pyBeta22 * r22 * 1e-6 * fParticleEnergies[i] / BetaMomentum22 * fParticleEnergies[i];
1121  pzBeta22 = pzBeta22 * r22 * 1e-6 * fParticleEnergies[i] / BetaMomentum22 * fParticleEnergies[i];
1122  }
1123  if (bra == 1) {
1124  Double32_t AlphaMomentum6 =
1125  TMath::Sqrt(pxAlpha6 * pxAlpha6 + pyAlpha6 * pyAlpha6 + pzAlpha6 * pzAlpha6);
1126 
1127  pxAlpha6 = pxAlpha6 * fParticleEnergies[i] / AlphaMomentum6;
1128  pyAlpha6 = pyAlpha6 * fParticleEnergies[i] / AlphaMomentum6;
1129  pzAlpha6 = pzAlpha6 * fParticleEnergies[i] / AlphaMomentum6;
1130  }
1131  }
1132 
1133  else if ((0.039949 < ran) && (ran <= 0.040574)) {
1134  {
1135  Double32_t BetaMomentum23 =
1136  TMath::Sqrt(pxBeta23 * pxBeta23 + pyBeta23 * pyBeta23 + pzBeta23 * pzBeta23);
1137  pxBeta23 = pxBeta23 * r23 * 1e-6 * fParticleEnergies[i] / BetaMomentum23 * fParticleEnergies[i];
1138  pyBeta23 = pyBeta23 * r23 * 1e-6 * fParticleEnergies[i] / BetaMomentum23 * fParticleEnergies[i];
1139  pzBeta23 = pzBeta23 * r23 * 1e-6 * fParticleEnergies[i] / BetaMomentum23 * fParticleEnergies[i];
1140  }
1141 
1142  if (bra == 1) {
1143  Double32_t AlphaMomentum7 =
1144  TMath::Sqrt(pxAlpha7 * pxAlpha7 + pyAlpha7 * pyAlpha7 + pzAlpha7 * pzAlpha7);
1145 
1146  pxAlpha7 = pxAlpha7 * fParticleEnergies[i] / AlphaMomentum7;
1147  pyAlpha7 = pyAlpha7 * fParticleEnergies[i] / AlphaMomentum7;
1148  pzAlpha7 = pzAlpha7 * fParticleEnergies[i] / AlphaMomentum7;
1149  }
1150  }
1151 
1152  else if ((0.040574 < ran) && (ran <= 0.046404)) {
1153  {
1154  Double32_t BetaMomentum24 =
1155  TMath::Sqrt(pxBeta24 * pxBeta24 + pyBeta24 * pyBeta24 + pzBeta24 * pzBeta24);
1156  pxBeta24 = pxBeta24 * r24 * 1e-6 * fParticleEnergies[i] / BetaMomentum24 * fParticleEnergies[i];
1157  pyBeta24 = pyBeta24 * r24 * 1e-6 * fParticleEnergies[i] / BetaMomentum24 * fParticleEnergies[i];
1158  pzBeta24 = pzBeta24 * r24 * 1e-6 * fParticleEnergies[i] / BetaMomentum24 * fParticleEnergies[i];
1159  }
1160  if (bra == 1) {
1161  Double32_t AlphaMomentum8 =
1162  TMath::Sqrt(pxAlpha8 * pxAlpha8 + pyAlpha8 * pyAlpha8 + pzAlpha8 * pzAlpha8);
1163 
1164  pxAlpha8 = pxAlpha8 * fParticleEnergies[i] / AlphaMomentum8;
1165  pyAlpha8 = pyAlpha8 * fParticleEnergies[i] / AlphaMomentum8;
1166  pzAlpha8 = pzAlpha8 * fParticleEnergies[i] / AlphaMomentum8;
1167  }
1168  }
1169 
1170  else if ((0.046404 < ran) && (ran <= 0.206004)) {
1171  {
1172  Double32_t BetaMomentum25 =
1173  TMath::Sqrt(pxBeta25 * pxBeta25 + pyBeta25 * pyBeta25 + pzBeta25 * pzBeta25);
1174  pxBeta25 = pxBeta25 * r25 * 1e-6 * fParticleEnergies[i] / BetaMomentum25 * fParticleEnergies[i];
1175  pyBeta25 = pyBeta25 * r25 * 1e-6 * fParticleEnergies[i] / BetaMomentum25 * fParticleEnergies[i];
1176  pzBeta25 = pzBeta25 * r25 * 1e-6 * fParticleEnergies[i] / BetaMomentum25 * fParticleEnergies[i];
1177  }
1178  if (bra == 1) {
1179  Double32_t AlphaMomentum9 =
1180  TMath::Sqrt(pxAlpha9 * pxAlpha9 + pyAlpha9 * pyAlpha9 + pzAlpha9 * pzAlpha9);
1181 
1182  pxAlpha9 = pxAlpha9 * fParticleEnergies[i] / AlphaMomentum9;
1183  pyAlpha9 = pyAlpha9 * fParticleEnergies[i] / AlphaMomentum9;
1184  pzAlpha9 = pzAlpha9 * fParticleEnergies[i] / AlphaMomentum9;
1185  }
1186  }
1187 
1188  else if ((0.206004 < ran) && (ran <= 1)) {
1189  Double32_t BetaMomentum26 = TMath::Sqrt(pxBeta26 * pxBeta26 + pyBeta26 * pyBeta26 + pzBeta26 * pzBeta26);
1190  pxBeta26 = pxBeta26 * r26 * 1e-6 * fParticleEnergies[i] / BetaMomentum26 * fParticleEnergies[i];
1191  pyBeta26 = pyBeta26 * r26 * 1e-6 * fParticleEnergies[i] / BetaMomentum26 * fParticleEnergies[i];
1192  pzBeta26 = pzBeta26 * r26 * 1e-6 * fParticleEnergies[i] / BetaMomentum26 * fParticleEnergies[i];
1193  }
1194  }
1195  }
1196 
1197  primGen->AddTrack(22, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); // dummy photon for track ID 0
1198  if ((0 < brb) && (brb <= 0.7142)) {
1199  primGen->AddTrack(betaPDGID, pxBeta1, pyBeta1, pzBeta1, fX, fY, fZ);
1200  } else if ((0.7142 < brb) && (brb <= 0.71523)) {
1201  primGen->AddTrack(betaPDGID, pxBeta2, pyBeta2, pzBeta2, fX, fY, fZ);
1202  } else if ((0.71523 < brb) && (brb <= 0.83426)) {
1203  {
1204  primGen->AddTrack(betaPDGID, pxBeta3, pyBeta3, pzBeta3, fX, fY, fZ);
1205  }
1206  if (brp == 1) {
1207  primGen->AddTrack(protonPDGID, pxProton1, pyProton1, pzProton1, fX, fY, fZ);
1208  }
1209  } else if ((0.83426 < brb) && (brb <= 0.88386)) {
1210  {
1211  primGen->AddTrack(betaPDGID, pxBeta4, pyBeta4, pzBeta4, fX, fY, fZ);
1212  }
1213  if ((0 < brp) && (brp <= 0.0769)) {
1214  primGen->AddTrack(protonPDGID, pxProton3, pyProton3, pzProton3, fX, fY, fZ);
1215  } else if ((0.0769 < brp) && (brp <= 1)) {
1216  primGen->AddTrack(protonPDGID, pxProton4, pyProton4, pzProton4, fX, fY, fZ);
1217  }
1218  } else if ((0.88386 < brb) && (brb <= 0.9118)) {
1219  {
1220  primGen->AddTrack(betaPDGID, pxBeta5, pyBeta5, pzBeta5, fX, fY, fZ);
1221  }
1222  if (brp == 1) {
1223  primGen->AddTrack(protonPDGID, pxProton5, pyProton5, pzProton5, fX, fY, fZ);
1224  }
1225  } else if ((0.9118 < brb) && (brb <= 0.93146)) {
1226  {
1227  primGen->AddTrack(betaPDGID, pxBeta6, pyBeta6, pzBeta6, fX, fY, fZ);
1228  }
1229  if ((0 < brp) && (brp <= 0.4375)) {
1230  primGen->AddTrack(protonPDGID, pxProton2, pyProton2, pzProton2, fX, fY, fZ);
1231  } else if ((0.4375 < brp) && (brp <= 0.6250)) {
1232  primGen->AddTrack(protonPDGID, pxProton6, pyProton6, pzProton6, fX, fY, fZ);
1233  } else if ((0.625 < brp) && (brp <= 0.875)) {
1234  primGen->AddTrack(protonPDGID, pxProton7, pyProton7, pzProton7, fX, fY, fZ);
1235  } else if ((0.875 < brp) && (brp <= 1)) {
1236  primGen->AddTrack(protonPDGID, pxProton8, pyProton8, pzProton8, fX, fY, fZ);
1237  }
1238  } else if ((0.93146 < brb) && (brb <= 0.94698)) {
1239  {
1240  primGen->AddTrack(betaPDGID, pxBeta7, pyBeta7, pzBeta7, fX, fY, fZ);
1241  }
1242  if (brp == 1) {
1243  primGen->AddTrack(protonPDGID, pxProton9, pyProton9, pzProton9, fX, fY, fZ);
1244  }
1245  } else if ((0.94698 < brb) && (brb <= 0.95277)) {
1246  {
1247  primGen->AddTrack(betaPDGID, pxBeta8, pyBeta8, pzBeta8, fX, fY, fZ);
1248  }
1249  if (brp == 1) {
1250  primGen->AddTrack(protonPDGID, pxProton10, pyProton10, pzProton10, fX, fY, fZ);
1251  }
1252  } else if ((0.95277 < brb) && (brb <= 0.96517)) {
1253  {
1254  primGen->AddTrack(betaPDGID, pxBeta9, pyBeta9, pzBeta9, fX, fY, fZ);
1255  }
1256  if (brp == 1) {
1257  primGen->AddTrack(protonPDGID, pxProton11, pyProton11, pzProton11, fX, fY, fZ);
1258  }
1259  } else if ((0.96517 < brb) && (brb <= 0.99932)) {
1260  {
1261  primGen->AddTrack(betaPDGID, pxBeta10, pyBeta10, pzBeta10, fX, fY, fZ);
1262  }
1263  if ((0 < brp) && (brp <= 0.333)) {
1264  primGen->AddTrack(protonPDGID, pxProton12, pyProton12, pzProton12, fX, fY, fZ);
1265  } else if ((0.333 < brp) && (brp <= 1)) {
1266  primGen->AddTrack(protonPDGID, pxProton13, pyProton13, pzProton13, fX, fY, fZ);
1267  }
1268  } else if ((0.99932 < brb) && (brb <= 0.99963)) {
1269  primGen->AddTrack(betaPDGID, pxBeta11, pyBeta11, pzBeta11, fX, fY, fZ);
1270  } else if ((0.99963 < brb) && (brb <= 0.999733)) {
1271  primGen->AddTrack(betaPDGID, pxBeta12, pyBeta12, pzBeta12, fX, fY, fZ);
1272  } else if ((0.999733 < brb) && (brb <= 1)) {
1273  {
1274  primGen->AddTrack(betaPDGID, pxBeta13, pyBeta13, pzBeta13, fX, fY, fZ);
1275  }
1276  if (brp == 1) {
1277  {
1278  primGen->AddTrack(protonPDGID, pxProton, pyProton, pzProton, fX, fY, fZ);
1279  }
1280  if (bra1 <= 1) {
1281  primGen->AddTrack(alphaPDGID, pxAlpha, pyAlpha, pzAlpha, fX, fY, fZ);
1282  }
1283  }
1284  }
1285 
1286  if ((0 < ran) && (ran <= 0.000016)) {
1287  {
1288  primGen->AddTrack(betaPDGID, pxBeta14, pyBeta14, pzBeta14, fX, fY, fZ);
1289  }
1290  if (bra1 == 1) {
1291  primGen->AddTrack(alphaPDGID, pxAlpha1, pyAlpha1, pzAlpha1, fX, fY, fZ);
1292  }
1293  } else if ((0.000016 < ran) && (ran <= 0.002646)) {
1294  {
1295  primGen->AddTrack(betaPDGID, pxBeta15, pyBeta15, pzBeta15, fX, fY, fZ);
1296  }
1297  if (bra == 1) {
1298  primGen->AddTrack(alphaPDGID, pxAlpha2, pyAlpha2, pzAlpha2, fX, fY, fZ);
1299  }
1300  } else if ((0.002646 < ran) && (ran <= 0.004696)) {
1301  primGen->AddTrack(betaPDGID, pxBeta16, pyBeta16, pzBeta16, fX, fY, fZ);
1302  } else if ((0.004696 < ran) && (ran <= 0.005866)) {
1303  primGen->AddTrack(betaPDGID, pxBeta17, pyBeta17, pzBeta17, fX, fY, fZ);
1304  } else if ((0.005866 < ran) && (ran <= 0.007606)) {
1305  {
1306  primGen->AddTrack(betaPDGID, pxBeta18, pyBeta18, pzBeta18, fX, fY, fZ);
1307  }
1308  if (bra == 1) {
1309  primGen->AddTrack(alphaPDGID, pxAlpha3, pyAlpha3, pzAlpha3, fX, fY, fZ);
1310  }
1311  } else if ((0.007606 < ran) && (ran <= 0.008489)) {
1312  {
1313  primGen->AddTrack(betaPDGID, pxBeta19, pyBeta19, pzBeta19, fX, fY, fZ);
1314  }
1315  if (bra == 1) {
1316  primGen->AddTrack(alphaPDGID, pxAlpha4, pyAlpha4, pzAlpha4, fX, fY, fZ);
1317  }
1318  } else if ((0.008489 < ran) && (ran <= 0.037259)) {
1319  {
1320  primGen->AddTrack(betaPDGID, pxBeta20, pyBeta20, pzBeta20, fX, fY, fZ);
1321  }
1322  if (bra == 1) {
1323  primGen->AddTrack(alphaPDGID, pxAlpha5, pyAlpha5, pzAlpha5, fX, fY, fZ);
1324  }
1325  } else if ((0.037259 < ran) && (ran <= 0.037539)) {
1326  primGen->AddTrack(betaPDGID, pxBeta21, pyBeta21, pzBeta21, fX, fY, fZ);
1327  } else if ((0.037539 < ran) && (ran <= 0.039949)) {
1328  {
1329  primGen->AddTrack(betaPDGID, pxBeta22, pyBeta22, pzBeta22, fX, fY, fZ);
1330  }
1331  if (bra == 1) {
1332  primGen->AddTrack(alphaPDGID, pxAlpha6, pyAlpha6, pzAlpha6, fX, fY, fZ);
1333  }
1334  } else if ((0.039949 < ran) && (ran <= 0.040574)) {
1335  {
1336  primGen->AddTrack(betaPDGID, pxBeta23, pyBeta23, pzBeta23, fX, fY, fZ);
1337  }
1338  if (bra == 1) {
1339  primGen->AddTrack(alphaPDGID, pxAlpha7, pyAlpha7, pzAlpha7, fX, fY, fZ);
1340  }
1341  } else if ((0.040574 < ran) && (ran <= 0.046404)) {
1342  {
1343  primGen->AddTrack(betaPDGID, pxBeta24, pyBeta24, pzBeta24, fX, fY, fZ);
1344  }
1345  if (bra == 1) {
1346  primGen->AddTrack(alphaPDGID, pxAlpha8, pyAlpha8, pzAlpha8, fX, fY, fZ);
1347  }
1348  } else if ((0.046404 < ran) && (ran <= 0.206004)) {
1349  {
1350  primGen->AddTrack(betaPDGID, pxBeta25, pyBeta25, pzBeta25, fX, fY, fZ);
1351  }
1352  if (bra == 1) {
1353  primGen->AddTrack(alphaPDGID, pxAlpha9, pyAlpha9, pzAlpha9, fX, fY, fZ);
1354  }
1355  } else if ((0.206004 < ran) && (ran <= 1)) {
1356  primGen->AddTrack(betaPDGID, pxBeta26, pyBeta26, pzBeta26, fX, fY, fZ);
1357  }
1358  }
1359  return kTRUE;
1360 }
1361 void AtTPC20MgDecay::SetDecayChainPoint(Double32_t ParticleEnergy, Double32_t ParticleBranchingRatio)
1362 {
1363 
1364  for (Int_t i = 0; i < fParticlesDefinedInNuclearDecay; i++) {
1365  fParticleEnergies[i] = ParticleEnergy;
1366  fParticleBranchingRatios[i] = ParticleBranchingRatio;
1367  }
1368 }
1369 
1370 ClassImp(AtTPC20MgDecay)
AtTPC20MgDecay::ReadEvent
virtual Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: AtTPC20MgDecay-Rnalpha.cxx:35
ClassImp
ClassImp(AtFindVertex)
AtTPC20MgDecay::Init
virtual Bool_t Init()
Definition: AtTPC20MgDecay-Rnalpha.cxx:29
AtTPC20MgDecay.h
AtTPC20MgDecay::SetDecayChainPoint
void SetDecayChainPoint(Double32_t ParticleEnergy=0, Double32_t ParticleBranchingRatio=0)
Definition: AtTPC20MgDecay-Rnalpha.cxx:152
AtTPC20MgDecay
Definition: AtTPC20MgDecay.h:14
Generated on Sat Nov 11 2023 15:30:37 for ATTPCROOT by   doxygen 1.8.18