source: ZHANGProjects/ICOSIM/CPP/trunk/source/Collimator.cc @ 2

Last change on this file since 2 was 2, checked in by zhangj, 10 years ago

Initial import
File size: 16.7 KB
Line 
1//Parent class for the collimators
2
3#include <iostream>
4#include <fstream>
5#include <sstream>
6#include <vector>
7#include <string>
8#include <cmath>
9#include "Collimator.h"
10using namespace std;
11
12Collimator::Collimator(const double& ALFX, const double& ALFY, const double& APER_1, const double& APER_2, const double& APER_3, const double& APER_4, const string& APERTYPE, const double& BETX, const double& BETY, const double& DPX, const double& DPY, const double& DX, const double& DY, const string& KEYWORD, const double& L, const double& MUX, const double& MUY, const string& NAME, const double& PTC, const double& PXC, const double& PYC, const double& S, const double& TC, const double& XC, const double& YC, const double& K0L, const double& K0SL, const double& K1L, const double& K1SL, const double& K2L, const double& K2SL, const string& PARENT, const string& meth, const long double& hgap, const long double& hgap2, const double& collang, const long double& pdepth, const long double& pdepth2, const double& tcang, const double& nsig)
13    : Element(ALFX, ALFY, APER_1, APER_2, APER_3, APER_4, APERTYPE, BETX, BETY, DPX, DPY, DX, DY, KEYWORD, L, MUX, MUY, NAME, PTC, PXC, PYC, S, TC, XC, YC, K0L, K0SL, K1L, K1SL, K2L, K2SL, PARENT),
14      method(meth),
15      hgap(hgap),
16      hgap2(hgap2),
17      phi(collang),
18      pdepth(pdepth),
19      pdepth2(pdepth2),
20      tcang(tcang),
21      nsig(nsig)
22{};
23
24Collimator::Collimator(const double& ALFX, const double& ALFY, const double& APER_1, const double& APER_2, const double& APER_3, const double& APER_4, const string& APERTYPE, const double& BETX, const double& BETY, const double& DPX, const double& DPY, const double& DX, const double& DY, const string& KEYWORD, const double& L, const double& MUX, const double& MUY, const string& NAME, const double& PTC, const double& PXC, const double& PYC, const double& S, const double& TC, const double& XC, const double& YC, const double& K0L, const double& K0SL, const double& K1L, const double& K1SL, const double& K2L, const double& K2SL, const string& PARENT, const string& meth, const long double& hgap, const long double& hgap2, const double& collang, const long double& pdepth, const long double& pdepth2, const double& tcang, const double& nsig, const long double& Bmax, const long double& thicknessMagneticField, const double& energyPerIon, const double& mass)
25    : Element(ALFX, ALFY, APER_1, APER_2, APER_3, APER_4, APERTYPE, BETX, BETY, DPX, DPY, DX, DY, KEYWORD, L, MUX, MUY, NAME, PTC, PXC, PYC, S, TC, XC, YC, K0L, K0SL, K1L, K1SL, K2L, K2SL, PARENT),
26      method(meth),
27      hgap(hgap),
28      hgap2(hgap2),
29      phi(collang),
30      pdepth(pdepth),
31      pdepth2(pdepth2),
32      tcang(tcang),
33      nsig(nsig),
34      Bmax(Bmax),
35      thicknessMagneticField(thicknessMagneticField),
36      energyPerIon(energyPerIon),
37      mass(mass)
38{};
39
40Collimator::Collimator(Element elt, const double& tcang, const double& nsig, const string& meth)
41    : Element(elt),
42      method(meth),
43      tcang(tcang),
44      nsig(nsig)
45{};
46
47Collimator::Collimator(Element elt, const double& tcang, const double& nsig, const string& meth, const string& material)
48    : Element(elt),
49      method(meth),
50      tcang(tcang),
51      nsig(nsig),
52      material(material)
53{};
54
55Collimator::Collimator(const Collimator& obj)
56    : Element(obj),
57      method(obj.method),
58      hgap(obj.hgap),
59      hgap2(obj.hgap2),
60      phi(obj.phi),
61      pdepth(obj.pdepth),
62      pdepth2(obj.pdepth2),
63      tcang(obj.tcang),
64      nsig(obj.nsig),
65      material(obj.material),
66      crossSectionPath(obj.crossSectionPath),
67      Bmax(obj.Bmax),
68      thicknessMagneticField(obj.thicknessMagneticField),
69      energyPerIon(obj.energyPerIon),
70      mass(obj.mass)
71{}
72
73
74void Collimator::affiche()
75{
76
77    this->Element::affiche();
78    cout << "Method: " << method << endl;
79
80};
81
82void Collimator::collipassInteraction(Particle& p1, double Apr, double Zpr, double betgam, double lcoll)
83{
84
85    double Navo(6.022137e23);
86    double s(0);
87    double atarget, rho, stre;
88    string abbreviation;
89
90    ifstream entree;
91    string filename;
92    filename = this->crossSectionPath + "materialinfo.csv";;
93    entree.open(filename.c_str());
94
95    if (entree.fail()) {
96        cerr << "Error: impossible to read the file " << filename << endl;
97    } else {
98
99        string word;
100        double nbre;
101
102        getline(entree, word);
103
104        while (!entree.eof()) {
105
106            entree >> ws;
107
108            getline(entree, word, ',');
109
110            if (word == "atarget") {
111                entree >> atarget;//mass number of collimator
112            } else if (word == "rho") {
113                entree >> rho;//density of collimator material
114            } else if (word == "stre") {
115                entree >> stre;//total EM cross-section
116            } else if (word == "abbreviation") {
117                entree >> abbreviation;//chemical symbol for collimator material
118            } else {
119                entree >> word;
120            }
121
122        }//end while(!entree.eof())
123
124    }
125    entree.close();
126
127    double Atarget;
128
129    Atarget = atarget * 0.001;
130
131
132    while (s < lcoll) { //loop until particle is lost or exits the collimator
133
134        string file;
135        ifstream entree1;
136        double sn, wx, wy, dpopdx;
137
138        std::stringstream apzp;
139        apzp << p1.Ap0 << p1.Zp0;
140        file = this->crossSectionPath + abbreviation + apzp.str() + "nuclEM.dat";
141        entree1.open(file.c_str());
142
143        if (entree1.fail()) {
144            cerr << "Warning!! The file " << file << " does not have a valid fomat or does not exist." << endl;
145
146            //we calculate analytical approximations to the EMD and hadronic cross section using a BCV parametrization of an impact-parameter cutoff model
147
148            double stee, Bbcv, steh, ste, intlengthe;
149            stee = stre * p1.Zp0 * (p1.Ap0 - p1.Zp0) / p1.Ap0 / (Zpr * (Apr - Zpr) / Apr); //approximation of electromagnetic dissociation cross section
150            Bbcv = 1.34 * (pow(p1.Ap0, 1. / 3.) + pow(atarget, 1. / 3.) - 0.75 * (pow(p1.Ap0, -1. / 3.) + pow(atarget, -1. / 3.))); //radius of the overlap between two "balls"
151            Bbcv = Bbcv / 1e15;
152            steh = M_PI * Bbcv * Bbcv * 1e31; //area of overlap = hadronic cross section
153            ste = steh + stee;//total cross section
154            intlengthe = Atarget / (Navo * rho * ste * 1e-31); //approximation to the interaction length
155
156            sn = s - intlengthe * log(-p1.random()); //distance between the beginning of the collimator and the next interaction
157
158            if ((sn > lcoll) && (p1.Ap0 > 1) && (p1.Zp0 <= p1.Ap0)) { //particle is going out of the collimator
159                BeBloMuSca(wx, wy, dpopdx, betgam, filename, rho, lcoll - s, atarget, p1); //calculating change in dx/ds, dy/ds and momentum
160                p1.coordonnees[1][0] = p1.coordonnees[0][0];//x coordinate when particle hits
161                p1.coordonnees[1][1] = p1.coordonnees[0][1] + wx;//xs=dx/ds when particle hits
162                p1.coordonnees[1][2] = p1.coordonnees[0][2];//y coordinate when particle hits
163                p1.coordonnees[1][3] = p1.coordonnees[0][3] + wy;//ys=dy/ds when particle hits
164                p1.coordonnees[1][4] = p1.coordonnees[0][4] - dpopdx;//dpop = momentum when the particle hits
165                return;
166            } else {//particle is absorbed
167                p1.Ap0 = 0;
168                p1.Zp0 = -1;
169                p1.coordonnees[1][0] = 0;
170                p1.coordonnees[1][1] = 0;
171                p1.coordonnees[1][2] = 0;
172                p1.coordonnees[1][3] = 0;
173                p1.coordonnees[1][4] = 0;
174                return;
175            }
176
177        } else {
178
179            double intlength, sigt;
180            vector <double> da, dz, isig;
181
182            genipsfastnewxc(sigt, isig, da, dz, file, p1.Ap0, p1.Zp0);
183
184            if (da.size() == 0) {
185                return;
186            }
187
188            intlength = Atarget / (Navo * rho * sigt * 1e-31);
189
190            sn = s - intlength * log(-p1.random()); //Sample from exponential distribution. sn is the position for the next interaction
191
192            if (sn > lcoll) { //particle exits collimator
193                BeBloMuSca(wx, wy, dpopdx, betgam, filename, rho, lcoll - s, atarget, p1);
194                p1.coordonnees[1][0] = p1.coordonnees[0][0];
195                p1.coordonnees[1][1] = p1.coordonnees[0][1] + wx;
196                p1.coordonnees[1][2] = p1.coordonnees[0][2];
197                p1.coordonnees[1][3] = p1.coordonnees[0][3] + wy;
198                p1.coordonnees[1][4] = p1.coordonnees[0][4] - dpopdx;
199                return;
200            } else if (lcoll - sn > 10 * intlength) { //particle is lost if it's moving more than 10 interaction lengths inside the collimator
201                p1.Ap0 = 0;
202                p1.Zp0 = -3;
203                p1.coordonnees[1][0] = 0;
204                p1.coordonnees[1][1] = 0;
205                p1.coordonnees[1][2] = 0;
206                p1.coordonnees[1][3] = 0;
207                p1.coordonnees[1][4] = 0;
208                return;
209            }
210
211            double dao, dzo, r;
212            int chan(0);
213
214            r = -p1.random();
215            int temp(0);
216
217            //deciding which new isotope is created
218            for (int i(0); i < isig.size(); ++i) {
219                if (r > isig[i]) {
220                    temp = temp + 1;
221                } else {
222                    chan = temp;
223                }
224            }
225
226            dao = da[chan - 1]; //decrease in mass number
227            dzo = dz[chan - 1]; //decrease in charge
228
229            if (dao == -999) { //corresponding to the omitted reactions, see genipsfastnewc. Particle assumed lost
230                p1.Ap0 = 0;
231                p1.Zp0 = -2;
232                p1.coordonnees[1][0] = 0;
233                p1.coordonnees[1][1] = 0;
234                p1.coordonnees[1][2] = 0;
235                p1.coordonnees[1][3] = 0;
236                p1.coordonnees[1][4] = 0;
237                return;
238            }
239
240            BeBloMuSca(wx, wy, dpopdx, betgam, filename, rho, sn - s, atarget, p1); //calculating change in dx/ds, dy/ds and momentum
241
242            p1.coordonnees[0][1] = p1.coordonnees[0][1] + wx;
243            p1.coordonnees[0][3] = p1.coordonnees[0][3] + wy;
244            p1.coordonnees[0][4] = p1.coordonnees[0][4] - dpopdx;
245            p1.Ap0 = p1.Ap0 - dao;
246            p1.Zp0 = p1.Zp0 - dzo;
247
248            s = sn;//'moving' particle to the place for the next interaction
249        }
250        entree1.close();
251    }//end while(s < this->lcoll)
252
253};
254
255
256void Collimator::genipsfastnewxc(double& sigt, vector <double>& isig, vector <double>& da, vector <double>& dz, string path, double a, double z)
257{
258
259    ifstream entree;
260    vector <double> temp1, temp2, temp3, sig, temp11, temp22, temp33;
261    double maximum(0);
262
263    entree.open(path.c_str());
264
265    if (entree.fail()) {
266        cerr << "We do not have cross-section information for the isotope with mass number " << a << " and charge " << z << " , i.e. the file " << path << "does not exist or has a wrong format."  << endl;
267    } else {
268
269        double nber1, nber2, nber3;
270        vector <double> sig;
271
272
273        entree >> nber1 >> nber2 >> nber3;
274
275        sigt = nber3;
276
277        while (!entree.eof()) {
278
279            entree >> nber1 >> nber2 >> nber3;
280
281            temp1.push_back(nber1);
282            temp2.push_back(nber2);
283            temp3.push_back(nber3);
284        }
285    }
286    entree.close();
287
288    double sum(0), msig;
289
290    for (int j(0); j < temp1.size(); ++j) {
291        if ((temp1[j] < a / 2) && (temp1[j] > 0)) { //only considering reactions where the mass decreases and the particle doesn't loose more than half its mass
292
293            temp11.push_back(temp1[j]);
294            temp22.push_back(temp2[j]);
295            temp33.push_back(temp3[j]);
296        }
297    }
298
299    for (int k(0); k < temp33.size(); ++k) {
300        if (temp33[k] > maximum) {
301            maximum = temp33[k];
302        }
303    }
304
305    for (int j(0); j < temp33.size(); ++j) {
306        if (temp33[j] > maximum / 100) { //excluding reactions with very small cross-section
307            sig.push_back(temp33[j]);//cross sections
308            da.push_back(temp11[j]);//change in mass number
309            dz.push_back(temp22[j]);//change in charge
310            sum = sum + temp33[j];
311        }
312    }
313    temp1.clear();
314    temp2.clear();
315    temp3.clear();
316    temp11.clear();
317    temp22.clear();
318    temp33.clear();
319
320    if (sig.size() == 0) {
321        cout << "No valid change in the cross-section." << endl;
322        return;
323    }
324    msig = sigt - sum; //sum of the excluded cross sections
325
326    sig.push_back(sig[sig.size() - 1]);
327    da.push_back(da[da.size() - 1]);
328    dz.push_back(dz[dz.size() - 1]);
329
330    for (int k(sig.size() - 2); k > 0; --k) {
331        sig[k] = sig[k - 1];
332        da[k] = da[k - 1];
333        dz[k] = dz[k - 1];
334    }
335
336    sig[0] = msig;
337    da[0] = -999;
338    da.push_back(-998);
339    dz[0] = sigt;
340    dz.push_back(-998);
341
342    vector <double> hsign;
343
344    for (int m(0); m < sig.size(); ++m) {
345        hsign.push_back(sig[m] / sigt); //the fraction of the total cross section for each channel, with the missing cross sections as first element
346    }
347
348    //add up the fractional partial cross sections in a vector isig. First element is 0.
349    //the different channels will then correspond to different intervals between 0 and 1 with lengths corresponding to their probability.
350
351    isig.push_back(0);
352
353    for (int p(1); p <= sig.size(); ++p) {
354        isig.push_back(hsign[p - 1] + isig[p - 1]);
355    }
356
357};
358
359
360void Collimator::BeBloMuSca(double& wx, double& wy, double& dpopdx, double betgam, string path, double rho, double L, double At, Particle p1)
361{
362
363    double me(510998.9);//electron mass [eV/c2]
364    double K(0.307075e5);//K/A in table 23.1, PDG page 163
365    double beta(sqrt(betgam * betgam / (betgam * betgam + 1))); //relativistic beta
366    double gamma(betgam / beta); //relativistic gamma
367    double X = log10(betgam);
368    double Zt, X0, X1, m, a, C, M0, X0MS, Tmax, delta, wrms, Ii, pr, w, dEdx;
369
370    ifstream entree;
371
372    entree.open(path.c_str());
373
374    if (entree.fail()) {
375        cerr << "Warning: problem with the file '" << path << "'." << endl;
376    } else {
377
378        string word;
379        double nbre;
380
381        while (!entree.eof()) {
382
383            getline(entree, word, ',');
384
385            if (word == "ztarget") {
386                getline(entree, word);
387                Zt = atof(word.c_str());//atom number of collimator
388            } else if (word == "X0") {
389                getline(entree, word);
390                X0 = atof(word.c_str());//Sternheimer parameter needed for density correction (not checked)
391            } else if (word == "X1") {
392                getline(entree, word);
393                X1 = atof(word.c_str());//Sternheimer parameter needed for density correction (not checked)
394            } else if (word == "m") {
395                getline(entree, word);
396                m = atof(word.c_str());//Sternheimer parameter needed for density correction (not checked)
397            } else if (word == "a") {
398                getline(entree, word);
399                a = atof(word.c_str());//Sternheimer parameter needed for density correction (not checked)
400            } else if (word == "C") {
401                getline(entree, word);
402                C = atof(word.c_str());//Sternheimer parameter needed for density correction (not checked)
403            } else {
404                getline(entree, word);
405            }
406        }//end while(!entree.eof())
407        entree.close();
408    }
409
410    M0 = 931.494e6 * p1.Ap0; //mass of ion [MeV/c2]
411    X0MS = 716.4 * At / (Zt * (Zt + 1) * log(287 / sqrt(Zt))) / (rho / 1000);
412    if (Zt < 13) { //mean excitation potential [eV]
413        Ii = Zt * (12 + 7 / Zt);
414    } else {
415        Ii = Zt * (9.76 + 58.8 * pow(Zt, -1.19));
416    }
417
418    Tmax = 2 * me * betgam * betgam / (1 + 2 * gamma * me / M0 + (me / M0) * (me / M0)); //Tmax (in eV) in Bethe-Bloch, PDG page 164
419
420    //DENSITY CORRECTION
421
422    if (X < X0) {
423        delta = 0;//valid for non-conductors
424    } else if ((X0 < X) && (X < X1)) {
425        delta = 4.6052 * X + C + a * pow(X1 - X, m);
426    } else {
427        delta = 4.6052 * X + C;
428    }
429
430    L = L * 100;
431    wrms = 19.2333e6 / (beta * betgam * M0) * p1.Zp0 * sqrt(L / X0MS) * (1 + 0.038 * log(L / X0MS)); //rms angular deviation from mult.scatt. Compare PDG page 166 - where does the prefactor come from? --from the fact that it's 3D not plane.
432
433    L = L / 100;
434
435
436    pr = -p1.random() * 2 * M_PI; //This line and the next samples from a 2D gaussian distribution the angular deviation in x and y.
437    w = wrms * sqrt(-log(1 + p1.random()));
438
439    //In order to make it analytically possible to sample, it is converted to polar coordinates, and the radius w and angle phi are sampled
440
441    wx = w * cos(pr); //change in dx/ds
442    wy = w * sin(pr); //change in dy/ds
443
444    dEdx = rho * K * p1.Zp0 * p1.Zp0 * Zt / At / (beta * beta) * (0.5 * log(2 * me * betgam * betgam * Tmax / (Ii * Ii)) - beta * beta - delta / 2); //energy loss acc. to Bete-Bloch on page 163 in PDG [eV/Meter]
445
446    dpopdx = L * dEdx / (M0 * gamma) / (beta * beta); //delta p/p caused by energy change due to Bethe-Bloch
447};
Note: See TracBrowser for help on using the repository browser.