R3BROOT
R3B analysis software
Loading...
Searching...
No Matches
IsElastic.cxx
Go to the documentation of this file.
1/******************************************************************************
2 * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH *
3 * Copyright (C) 2019-2025 Members of R3B Collaboration *
4 * *
5 * This software is distributed under the terms of the *
6 * GNU General Public Licence (GPL) version 3, *
7 * copied verbatim in the file "LICENSE". *
8 * *
9 * In applying this license GSI does not waive the privileges and immunities *
10 * granted to it by virtue of its status as an Intergovernmental Organization *
11 * or submit itself to any jurisdiction. *
12 ******************************************************************************/
13
14#include "IsElastic.h"
15#include "TVector3.h"
16
17// Code extracted from R3BNeutronTracker2D. It is absolutely atrocious. Don't look at it. Nothing too see here. Go away.
18
19bool Neuland::IsElastic(const R3BNeulandCluster* cl1, const R3BNeulandCluster* cl2)
20{
21 const R3BNeulandCluster* c1;
22 const R3BNeulandCluster* c2;
23
24 // if (cl1->GetT() < cl2->GetT())
25 // {
26 c1 = cl1;
27 c2 = cl2;
28 /*}
29 else
30 {
31 c1 = cl2;
32 c2 = cl1;
33 }*/
34
35 Double_t c = 29.9792458;
36 Double_t dio = 10.6; // 3 times half the diogonal of a paddle
37 Double_t amu = 931.494028;
38 // Double_t mNeutron = 1.0086649156 * amu;
39
40 Double_t x0 = 0.;
41 Double_t y0 = 0.;
42 Double_t z0 = 0.;
43 Double_t t0 = 0.;
44
45 // Double_t v1xmin,v1ymin,v1zmin,v1xmax,v1ymax,v1zmax;
46 Double_t v4xmin, v4ymin, v4zmin; //,v4xmax,v4ymax,v4zmax;
47 Double_t v3xmin, v3ymin, v3zmin; //,v3xmax,v3ymax,v3zmax;
48 // Double_t v6xmin,v6ymin,v6zmin,v6xmax,v6ymax,v6zmax;
49
50 // search for scattering
51 // elastic scattering of particle 1 (neutron) on a particle 2 (proton) at
52 // rest. Outgoing particles are 3 (scattered neutron) and 4 (recoil
53 // proton). angle of particles after scattering is theta3 and theta4 K is
54 // kinetic energy, p is momentum
55
56 // TVector3 pos1;
57 // TVector3 pos2;
58
59 const TVector3 pos1 = c1->GetFirstHit().GetPosition();
60
61 // incoming particle
62 // vector from previous interaction to present interaction
63 Double_t v1x = pos1.X() - x0;
64 Double_t v1y = pos1.Y() - y0;
65 Double_t v1z = pos1.Z() - z0;
66 // time difference and distance between previous interaction to present
67 // interaction
68 Double_t dt = c1->GetFirstHit().GetT() - t0;
69 Double_t dr = sqrt(v1x * v1x + v1y * v1y + v1z * v1z);
70
71 Double_t beta1 = dr / dt / c;
72 Double_t beta1max = (dr + dio) / dt / c;
73 Double_t beta1min = (dr - dio) / dt / c;
74 if (beta1max > 0.99)
75 beta1max = 0.99;
76 if (beta1min < 0.)
77 beta1min = 0.;
78
79 Double_t gamma1 = 1. / sqrt(1. - beta1 * beta1);
80 Double_t gamma1min = 1. / sqrt(1. - beta1min * beta1min);
81 Double_t gamma1max = 1. / sqrt(1. - beta1max * beta1max);
82 Double_t p1 = beta1 * gamma1 * 1. * amu;
83 Double_t p1min = beta1min * gamma1min * 1. * amu;
84 Double_t p1max = beta1max * gamma1max * 1. * amu;
85 Double_t En1 = sqrt(p1 * p1 + amu * amu);
86 Double_t En1min = sqrt(p1min * p1min + amu * amu);
87 Double_t En1max = sqrt(p1max * p1max + amu * amu);
88 Double_t K1 = En1 - amu;
89 Double_t K1min = En1min - amu;
90 Double_t K1max = En1max - amu;
91
92 // particle 4 is proton
93 const TVector3 pos2 = c1->GetLastHit().GetPosition();
94 Double_t v4x = (pos2 - pos1).X();
95 Double_t v4y = (pos2 - pos1).Y();
96 Double_t v4z = (pos2 - pos1).Z();
97 // vector perpendicular to scattering plane
98 Double_t v5x = v1y * v4z - v1z * v4y;
99 Double_t v5y = v1z * v4x - v1x * v4z;
100 Double_t v5z = v1x * v4y - v1y * v4x;
101
102 Double_t tempAngle;
103 dt = c1->GetLastHit().GetT() - c1->GetFirstHit().GetT();
104 dr = sqrt(v4x * v4x + v4y * v4y + v4z * v4z);
105
106 Double_t theta4Measured = acos((v1x * v4x + v1y * v4y + v1z * v4z) / sqrt(v1x * v1x + v1y * v1y + v1z * v1z) /
107 sqrt(v4x * v4x + v4y * v4y + v4z * v4z));
108 Double_t theta4Measuredmin = theta4Measured;
109 Double_t theta4Measuredmax = theta4Measured;
110 if ((v4x * v4x + v4y * v4y + v4z * v4z) == 0.)
111 {
112 theta4Measured = 1.55;
113 theta4Measuredmin = 0.;
114 theta4Measuredmax = 1.55;
115 }
116 for (Int_t k = 1; k < 3; k++)
117 {
118 if (k == 1)
119 {
120 v4xmin = v4x - dio;
121 }
122 else
123 {
124 v4xmin = v4x + dio;
125 }
126 for (Int_t l = 1; l < 3; l++)
127 {
128 if (l == 1)
129 {
130 v4ymin = v4y - dio;
131 }
132 else
133 {
134 v4ymin = v4y + dio;
135 }
136 for (Int_t m = 1; m < 3; m++)
137 {
138 if (m == 1)
139 {
140 v4zmin = v4z - dio;
141 }
142 else
143 {
144 v4zmin = v4z + dio;
145 }
146 tempAngle =
147 acos((v1x * v4xmin + v1y * v4ymin + v1z * v4zmin) / sqrt(v1x * v1x + v1y * v1y + v1z * v1z) /
148 sqrt(v4xmin * v4xmin + v4ymin * v4ymin + v4zmin * v4zmin));
149 if (tempAngle > theta4Measuredmax)
150 theta4Measuredmax = tempAngle;
151 if (tempAngle < theta4Measuredmin)
152 theta4Measuredmin = tempAngle;
153 }
154 }
155 }
156 if (theta4Measuredmax > 1.55)
157 theta4Measuredmax = 1.55;
158
159 // calculate velocity of neutron after scattering
160 const TVector3 pos3 = c2->GetFirstHit().GetPosition();
161 Double_t v3x = (pos3 - pos1).X();
162 Double_t v3y = (pos3 - pos1).Y();
163 Double_t v3z = (pos3 - pos1).Z();
164
165 Double_t v6x = v1y * v3z - v1z * v3y;
166 Double_t v6y = v1z * v3x - v1x * v3z;
167 Double_t v6z = v1x * v3y - v1y * v3x;
168
169 Double_t theta56 = acos((v5x * v6x + v5y * v6y + v5z * v6z) / sqrt(v5x * v5x + v5y * v5y + v5z * v5z) /
170 sqrt(v6x * v6x + v6y * v6y + v6z * v6z));
171 if ((v5x * v5x + v5y * v5y + v5z * v5z) == 0. || (v6x * v6x + v6y * v6y + v6z * v6z) == 0.)
172 {
173 theta56 = 3.14;
174 }
175
176 Double_t theta3 = acos((v1x * v3x + v1y * v3y + v1z * v3z) / sqrt(v1x * v1x + v1y * v1y + v1z * v1z) /
177 sqrt(v3x * v3x + v3y * v3y + v3z * v3z));
178 Double_t theta3min = theta3;
179 Double_t theta3max = theta3;
180
181 for (Int_t k = 1; k < 3; k++)
182 {
183 if (k == 1)
184 {
185 v3xmin = v3x - dio;
186 }
187 else
188 {
189 v3xmin = v3x + dio;
190 }
191 for (Int_t l = 1; l < 3; l++)
192 {
193 if (l == 1)
194 {
195 v3ymin = v3y - dio;
196 }
197 else
198 {
199 v3ymin = v3y + dio;
200 }
201 for (Int_t m = 1; m < 3; m++)
202 {
203 if (m == 1)
204 {
205 v3zmin = v3z - dio;
206 }
207 else
208 {
209 v3zmin = v3z + dio;
210 }
211 tempAngle =
212 acos((v1x * v3xmin + v1y * v3ymin + v1z * v3zmin) / sqrt(v1x * v1x + v1y * v1y + v1z * v1z) /
213 sqrt(v3xmin * v3xmin + v3ymin * v3ymin + v3zmin * v3zmin));
214 if (tempAngle > theta3max)
215 theta3max = tempAngle;
216 if (tempAngle < theta3min)
217 theta3min = tempAngle;
218 }
219 }
220 }
221
222 if (theta3max > 1.55)
223 theta3max = 1.55;
224
225 dt = c2->GetLastHit().GetT() - c2->GetFirstHit().GetT();
226 dr = sqrt(v3x * v3x + v3y * v3y + v3z * v3z);
227
228 Double_t beta3 = dr / dt / c;
229 Double_t beta3max = (dr + dio) / dt / c;
230 Double_t beta3min = (dr - dio) / dt / c;
231
232 if (beta3max > 0.99)
233 beta3max = 0.99;
234 if (beta3min < 0.)
235 beta3min = 0.;
236
237 Double_t gamma3 = 1. / sqrt(1. - beta3 * beta3);
238 Double_t gamma3min = 1. / sqrt(1. - beta3min * beta3min);
239 Double_t gamma3max = 1. / sqrt(1. - beta3max * beta3max);
240 Double_t p3 = beta3 * gamma3 * 1. * amu;
241 Double_t p3min = beta3min * gamma3min * 1. * amu;
242 Double_t p3max = beta3max * gamma3max * 1. * amu;
243 if (p3 > p1)
244 p3 = p1;
245 // Double_t En3=sqrt(p3*p3+amu*amu);
246 // Double_t En3min=sqrt(p3min*p3min+amu*amu);
247 // Double_t En3max=sqrt(p3max*p3max+amu*amu);
248 // Double_t K3=En3-amu;
249 // Double_t K3min=En3min-amu;
250 // Double_t K3max=En3max-amu;
251
252 Double_t Ma, Mb, Mc, Md, Ka, Thc, Ei, Pa, AA, BB, a, b, cc, Pc1 /*, Pc2*/;
253 Double_t Pd1, Thd1, /*Ec1,*/ Ed1, /*Kc1,*/ Kd1 /*, Qsqr1*/;
254 // Double_t Pd2,Thd2,Ec2,Ed2,Kc2,Kd2,Qsqr2;
255 Ma = 1.0087 * amu;
256 Mb = 1.0073 * amu;
257 Mc = Ma;
258 Md = Mb;
259 // calculate momentum of scattered proton and neutron
260 Ka = K1;
261 Thc = theta3;
262 Ei = Ma + Mb + Ka;
263 Pa = sqrt(Ka * Ka + 2. * Ka * Ma);
264 AA = Ei * Ei - Md * Md + Mc * Mc - Pa * Pa;
265 BB = AA * AA - 4. * Ei * Ei * Mc * Mc;
266 a = 4. * Pa * Pa * cos(Thc) * cos(Thc) - 4. * Ei * Ei;
267 b = 4. * AA * Pa * cos(Thc);
268 cc = BB;
269 Pc1 = -b / (2. * a) + sqrt((b * b) / (4. * a * a) - (cc / a));
270 if (Pc1 < 0.)
271 Pc1 = 0.;
272 // Pc2 = -b / (2. * a) - sqrt((b * b) / (4. * a * a) - (cc / a));
273 Pd1 = sqrt(Pc1 * Pc1 + Pa * Pa - 2. * Pc1 * Pa * cos(Thc));
274 Thd1 = acos((Pc1 * Pc1 - Pd1 * Pd1 - Pa * Pa) / (-2. * Pd1 * Pa));
275 // Ec1 = sqrt(Pc1 * Pc1 + Ma * Ma);
276 Ed1 = sqrt(Pd1 * Pd1 + Mb * Mb);
277 // Kc1 = Ec1 - Ma;
278 Kd1 = Ed1 - Mb;
279 // Qsqr1 = (-(Ka - Kc1) * (Ka - Kc1) + (Pa * Pa + Pc1 * Pc1 - 2. * Pa * Pc1 * cos(Thc))) / 197.327 / 197.327;
280 Double_t p3b = Pc1;
281 // Double_t K3b=Kc1;
282 // Double_t E3b=Ec1;
283 Double_t theta4 = Thd1;
284 Double_t p4b = Pd1;
285 // Double_t K4b=Kd1;
286 // Double_t E4b=Ed1;
287
288 Double_t p3bmin = p3b;
289 Double_t p3bmax = p3b;
290 Double_t p4bmin = p4b;
291 Double_t p4bmax = p4b;
292 Double_t theta4min = theta4;
293 Double_t theta4max = theta4;
294
295 // calculate minimum and maximum(within errors) momentum of scattered proton
296 // and neutron
297 for (Int_t m = 1; m < 5; m++)
298 {
299 if (m == 1)
300 {
301 Ka = K1min;
302 Thc = theta3min;
303 }
304 if (m == 2)
305 {
306 Ka = K1min;
307 Thc = theta3max;
308 }
309 if (m == 3)
310 {
311 Ka = K1max;
312 Thc = theta3min;
313 }
314 if (m == 4)
315 {
316 Ka = K1max;
317 Thc = theta3max;
318 }
319
320 Ei = Ma + Mb + Ka;
321 Pa = sqrt(Ka * Ka + 2. * Ka * Ma);
322 AA = Ei * Ei - Md * Md + Mc * Mc - Pa * Pa;
323 BB = AA * AA - 4. * Ei * Ei * Mc * Mc;
324 a = 4. * Pa * Pa * cos(Thc) * cos(Thc) - 4. * Ei * Ei;
325 b = 4. * AA * Pa * cos(Thc);
326 cc = BB;
327 Pc1 = -b / (2. * a) + sqrt((b * b) / (4. * a * a) - (cc / a));
328 if (((b * b) / (4. * a * a) - (cc / a)) < 0.)
329 Pc1 = 0.;
330 if (Pc1 < 0.)
331 Pc1 = 0.;
332 // Pc2 = -b / (2. * a) - sqrt((b * b) / (4. * a * a) - (cc / a));
333 Pd1 = sqrt(Pc1 * Pc1 + Pa * Pa - 2. * Pc1 * Pa * cos(Thc));
334 Thd1 = acos((Pc1 * Pc1 - Pd1 * Pd1 - Pa * Pa) / (-2. * Pd1 * Pa));
335
336 if (Pc1 < p3bmin)
337 p3bmin = Pc1;
338 if (Pc1 > p3bmax)
339 p3bmax = Pc1;
340 if (Pd1 < p4bmin)
341 p4bmin = Pd1;
342 if (Pd1 > p4bmax)
343 p4bmax = Pd1;
344 if (Thd1 < theta4min)
345 theta4min = Thd1;
346 if (Thd1 > theta4max)
347 theta4max = Thd1;
348 }
349
350 // Ec1 = sqrt(p3bmin * p3bmin + Ma * Ma);
351 // Kc1 = Ec1 - Ma;
352 Ed1 = sqrt(p4bmin * p4bmin + Mb * Mb);
353 Kd1 = Ed1 - Mb;
354 // Qsqr1 = (-(Ka - Kc1) * (Ka - Kc1) + (Pa * Pa + Pc1 * Pc1 - 2. * Pa * Pc1 * cos(Thc))) / 197.327 / 197.327;
355 // Double_t K3bmin=Kc1;
356 // Double_t E3bmin=Ec1;
357 Double_t K4bmin = Kd1;
358 // Double_t E4bmin=Ed1;
359 // Ec1 = sqrt(p3bmax * p3bmax + Ma * Ma);
360 // Kc1 = Ec1 - Ma;
361 Ed1 = sqrt(p4bmax * p4bmax + Mb * Mb);
362 Kd1 = Ed1 - Mb;
363 // Qsqr1 = (-(Ka - Kc1) * (Ka - Kc1) + (Pa * Pa + Pc1 * Pc1 - 2. * Pa * Pc1 * cos(Thc))) / 197.327 / 197.327;
364
365 // Double_t K3bmax=Kc1;
366 // Double_t E3bmax=Ec1;
367 Double_t K4bmax = Kd1;
368 // Double_t E4bmax=Ed1;
369
370 // decide if Cluster comes from scattered neutron or another neutron is
371 // needed!
372 Double_t protonEnergy = 8.76839 + 4.13858 * c1->GetE() - 0.00337368 * c1->GetE() * c1->GetE();
373
374 return p3bmax > p3min && p3bmin < p3max && theta4max > theta4Measuredmin && theta4min < theta4Measuredmax &&
375 K4bmax > protonEnergy && K4bmin < protonEnergy && theta56 * 180. / 3.14 > 120.;
376}
c
static const Double_t c
Definition R3BNeulandNeutron.cxx:18
c2
static const Double_t c2
Definition R3BNeulandNeutronReconstructionMon.cxx:39
R3BNeulandCluster::GetE
Double_t GetE() const
Definition R3BNeulandCluster.cxx:68
R3BNeulandCluster::GetFirstHit
R3BNeulandHit GetFirstHit() const
Definition R3BNeulandCluster.cxx:32
R3BNeulandCluster::GetLastHit
R3BNeulandHit GetLastHit() const
Definition R3BNeulandCluster.cxx:44
R3B::Neuland::IsElastic
bool IsElastic(const R3BNeulandCluster *, const R3BNeulandCluster *)
Definition IsElastic.cxx:19
R3BNeulandHit::GetT
auto GetT() const -> double
Definition R3BNeulandHit.h:57
R3BNeulandHit::GetPosition
auto GetPosition() const -> TVector3
Definition R3BNeulandHit.h:61