FairShip
Loading...
Searching...
No Matches
MeanMaterialBudget.cxx
Go to the documentation of this file.
1// SPDX-License-Identifier: BSD-3-Clause
2// SPDX-FileCopyrightText: ALICE Experiment at CERN, All rights reserved
3
4#include "MeanMaterialBudget.h"
5
6#include "FairLogger.h"
7#include "TGeoManager.h"
8#include "TGeoMaterial.h"
9#include "TGeoNode.h"
10#include "TGeoShape.h"
11#include "TGeoVolume.h"
12#include "TMath.h"
13
14namespace shipgen {
15
16double MeanMaterialBudget(std::span<const double, 3> start,
17 std::span<const double, 3> end,
18 std::span<double, 10> mparam) {
19 //
20 // Calculate mean material budget and material properties between
21 // the points "start" and "end".
22 //
23 // "mparam" - parameters used for the energy and multiple scattering
24 // corrections:
25 //
26 // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
27 // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [dimensionless]
28 // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [dimensionless]
29 // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [dimensionless]
30 // mparam[4] - length: sum(x_i) [cm]
31 // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [dimensionless]
32 // mparam[6] - number of boundary crosses
33 // mparam[7] - maximum density encountered (g/cm^3)
34 // mparam[8] - equivalent interaction length fraction: sum(x_i/I0_i)
35 // [dimensionless] mparam[9] - maximum cross section encountered (mbarn)
36 //
37 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
38 //
39 // Original ALICE improvements by
40 // Andrea Dainese, Andrea.Dainese@lnl.infn.it
41 // Andrei Gheata, Andrei.Gheata@cern.ch
42 //
43 // SHiP enhancements:
44 // Thomas Ruf, Thomas.Ruf@cern.ch (interaction length support, Dec
45 // 2016) Anupama Reghunath, anupama.reghunath@cern.ch (error logging,
46 // Nov 2024)
47 //
48
49 mparam[0] = 0;
50 mparam[1] = 1;
51 mparam[2] = 0;
52 mparam[3] = 0;
53 mparam[4] = 0;
54 mparam[5] = 0;
55 mparam[6] = 0;
56 mparam[7] = 0;
57 mparam[8] = 0;
58 mparam[9] = 0;
59 //
60 double bparam[7]; // total parameters
61 double lparam[7]; // local parameters
62 double mbarn = 1E-3 * 1E-24 * TMath::Na(); // cm^2 * Avogadro
63
64 for (Int_t i = 0; i < 7; i++) bparam[i] = 0;
65
66 if (!gGeoManager) {
67 return 0.;
68 }
69 //
70 double length;
71 double dir[3];
72 length = TMath::Sqrt((end[0] - start[0]) * (end[0] - start[0]) +
73 (end[1] - start[1]) * (end[1] - start[1]) +
74 (end[2] - start[2]) * (end[2] - start[2]));
75 mparam[4] = length;
76 if (length < TGeoShape::Tolerance()) return 0.0;
77 double invlen = 1. / length;
78 dir[0] = (end[0] - start[0]) * invlen;
79 dir[1] = (end[1] - start[1]) * invlen;
80 dir[2] = (end[2] - start[2]) * invlen;
81
82 // Initialise start point and direction
83 TGeoNode* currentnode = 0;
84 TGeoNode* startnode = gGeoManager->InitTrack(start.data(), dir);
85 if (!startnode) {
86 LOG(error) << "MeanMaterialBudget: start point out of geometry: x "
87 << start[0] << ", y " << start[1] << ", z " << start[2];
88 return 0.0;
89 }
90 TGeoMaterial* material = startnode->GetVolume()->GetMedium()->GetMaterial();
91 lparam[0] = material->GetDensity();
92 if (lparam[0] > mparam[7]) mparam[7] = lparam[0];
93 lparam[1] = material->GetRadLen();
94 lparam[2] = material->GetA();
95 lparam[3] = material->GetZ();
96 lparam[4] = length;
97 lparam[5] = lparam[3] / lparam[2];
98 lparam[6] = material->GetIntLen();
99 double n = lparam[0] / lparam[2];
100 double sigma = 1. / (n * lparam[6]) / mbarn;
101 if (sigma > mparam[9]) mparam[9] = sigma;
102 if (material->IsMixture()) {
103 TGeoMixture* mixture = dynamic_cast<TGeoMixture*>(material);
104 lparam[5] = 0;
105 double sum = 0;
106 for (Int_t iel = 0; iel < mixture->GetNelements(); iel++) {
107 sum += mixture->GetWmixt()[iel];
108 lparam[5] += mixture->GetZmixt()[iel] * mixture->GetWmixt()[iel] /
109 mixture->GetAmixt()[iel];
110 }
111 lparam[5] /= sum;
112 }
113
114 // Locate next boundary within length without computing safety.
115 // Propagate either with length (if no boundary found) or just cross boundary
116 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
117 double step = 0.0; // Step made
118 double snext = gGeoManager->GetStep();
119 // If no boundary within proposed length, return current density
120 if (!gGeoManager->IsOnBoundary()) {
121 mparam[0] = lparam[0];
122 mparam[1] = lparam[4] / lparam[1];
123 mparam[2] = lparam[2];
124 mparam[3] = lparam[3];
125 mparam[4] = lparam[4];
126 mparam[8] = lparam[4] / lparam[6];
127 return lparam[0];
128 }
129 // Try to cross the boundary and see what is next
130 Int_t nzero = 0;
131 while (length > TGeoShape::Tolerance()) {
132 currentnode = gGeoManager->GetCurrentNode();
133 if (snext < 2. * TGeoShape::Tolerance())
134 nzero++;
135 else
136 nzero = 0;
137 if (nzero > 3) {
138 // This means navigation has problems on one boundary
139 // Try to cross by making a small step
140 mparam[0] = bparam[0] / step;
141 mparam[1] = bparam[1];
142 mparam[2] = bparam[2] / step;
143 mparam[3] = bparam[3] / step;
144 mparam[5] = bparam[5] / step;
145 mparam[8] = bparam[6];
146 mparam[4] = step;
147 mparam[0] = 0.; // if crash of navigation take mean density 0
148 mparam[1] = 1000000; // and infinite rad length
149 return bparam[0] / step;
150 }
151 mparam[6] += 1.;
152 step += snext;
153 bparam[1] += snext / lparam[1];
154 bparam[2] += snext * lparam[2];
155 bparam[3] += snext * lparam[3];
156 bparam[5] += snext * lparam[5];
157 bparam[6] += snext / lparam[6];
158 bparam[0] += snext * lparam[0];
159
160 if (snext >= length) break;
161 if (!currentnode) break;
162 length -= snext;
163 material = currentnode->GetVolume()->GetMedium()->GetMaterial();
164 lparam[0] = material->GetDensity();
165 if (lparam[0] > mparam[7]) mparam[7] = lparam[0];
166 lparam[1] = material->GetRadLen();
167 lparam[2] = material->GetA();
168 lparam[3] = material->GetZ();
169 lparam[5] = lparam[3] / lparam[2];
170 lparam[6] = material->GetIntLen();
171 n = lparam[0] / lparam[2];
172 sigma = 1. / (n * lparam[6]) / mbarn;
173 if (sigma > mparam[9]) mparam[9] = sigma;
174 if (material->IsMixture()) {
175 TGeoMixture* mixture = dynamic_cast<TGeoMixture*>(material);
176 lparam[5] = 0;
177 double sum = 0;
178 for (Int_t iel = 0; iel < mixture->GetNelements(); iel++) {
179 sum += mixture->GetWmixt()[iel];
180 lparam[5] += mixture->GetZmixt()[iel] * mixture->GetWmixt()[iel] /
181 mixture->GetAmixt()[iel];
182 }
183 lparam[5] /= sum;
184 }
185 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
186 snext = gGeoManager->GetStep();
187 }
188 mparam[0] = bparam[0] / step;
189 mparam[1] = bparam[1];
190 mparam[2] = bparam[2] / step;
191 mparam[3] = bparam[3] / step;
192 mparam[5] = bparam[5] / step;
193 mparam[8] = bparam[6];
194 return bparam[0] / step;
195}
196
197} // namespace shipgen
mbarn
const Double_t mbarn
Definition: FixedTargetGenerator.cxx:30
E
Definition: diagrams_e.h:4
shipgen::MeanMaterialBudget
double MeanMaterialBudget(std::span< const double, 3 > start, std::span< const double, 3 > end, std::span< double, 10 > mparam)
Definition: MeanMaterialBudget.cxx:16