GamR  0.0.0
GammaROOT
Loading...
Searching...
No Matches
Coulex.cc
Go to the documentation of this file.
1
8
9#include <iostream>
10#include <fstream>
11#include <cmath>
12#include <vector>
13#include <algorithm>
14#include <string>
15#include <complex>
16
17#include <TMath.h>
18
19#include <angdist/AngDist.hh>
20#include <toolkit/Misc.hh>
21#include <utils/Utilities.hh>
22#include "Eloss.hh"
23#include "WDBSubs.hh"
24#include "Coulex.hh"
25
26namespace GamR {
27 namespace Coulex {
50 int Kinematics(double A1I, double A2I, double A3I, double A4I, double E1I, double QEFFI, double PSII, double &E3I, double &E3AI, double &E4I, double &E4AI, double &ZETAI, double &ZETAAI, double &thetaI, double &thetaaI, double &PHII, double &PHIAI)
51 {
52 /*
53 C
54 C SUBROUTINE TO CALCULATE REACTION KINEMATICS
55 C SEE MARION AND YOUNG
56 C THIS CODE WAS ORIGINALLY WRITTEN BY ANDREW STUCHBERY CIRCA 1980
57 * converted to 'C' by Brendan McCormick 01/2019
58 C
59 C A1=INCIDENT MASS
60 C A2=TARGET MASS
61 C A3="LIGHT" PRODUCT MASS
62 C A4="HEAVY" PRODUCT MASS
63 C E1=INCIDENT BEAM ENERGY
64 C QEFF=EFFECTIVE Q VALUE FOR REACTION (INCLUDE ENERGY LEVEL)
65 C PSI=LAB RECOIL(BACKSCATTER) ANGLE OF "LIGHT" PRODUCT (RADIANS)
66 C E3=LAB ENERGY OF "LIGHT" PRODUCT
67 C E3A=OTHER POSSIBLE E3 FOR B.GT.D (SEE BELOW)
68 C E4=LAB ENERGY OF "HEAVY" PRODUCT
69 C E4A=OTHER POSSIBLE E4
70 C ZETA=LAB RECOIL ANGLE OF HEAVY PRODUCT (RADIANS)
71 C ZETAA=OTHER POSSIBLE ZETA
72 C THETA=CM ANGLE OF LIGHT (I.E. DETECTED) PRODUCT (RADIANS)
73 C THETAA = OTHER POSSIBLE THETA
74 C PHI=CM ANGLE OF HEAVY (I.E. UNDETECTED) PRODUCT (RADIANS)
75 C PHIA = OTHER POSSIBLE PHI
76 C IFLAG (now return value) = +1 NORMAL CASE - UNIQUE SOLUTION
77 C 0 CASE WHERE B>D (E.G. IF A1>A2) - TWO POSSIBLE SOLUTIONS
78 C -1 AS FOR 0 BUT AN IMPOSSIBLE SCATTERING ANGLE (PSI) HAS
79 C BEEN REQUESTED
80 C
81 */
82 int iFlag = 1;
83
84 if(A1I==0 || A2I==0 || A3I==0 || A4I==0) {
85 std::cerr << "!!ERROR!! One or more of the particles in the kinematics routine has no mass." << std::endl;
86 return -1;
87 }
88
89 double A1,A2,A3,A4,E1,QEFF,PSI,E3,E3A,E4,E4A;
90 double ZETA,ZETAA,THETA,THETAA,PHI,PHIA,ET,AD,ED,QD;
91 double A,B,C,D,CPSI,SPSI,DUM1,DUM2,DUM3,E3T,E3TA;
92 double SZETA,SZA,CTHETA,CTHETAA,PSIMX,DUM0;
93
94 // INPUT PARAMETERS
95 A1=A1I;
96 A2=A2I;
97 A3=A3I;
98 A4=A4I;
99 E1=E1I;
100 QEFF=QEFFI;
101 PSI=PSII;
102
103 ET=E1+QEFF;
104 AD=(A1+A2)*(A3+A4);
105 ED=A1*E1/(ET*AD);
106 QD=(A2+A1*QEFF/ET)/AD;
107
108 A=A4*ED;
109 B=A3*ED;
110 C=A3*QD;
111 D=A4*QD;
112
113 SPSI=std::sin(PSI);
114 CPSI=std::cos(PSI);
115 iFlag=1;
116 DUM0=(D/B-SPSI*SPSI);
117 if(DUM0<0)DUM0=0.;
118 DUM1=std::sqrt(DUM0);
119 if(B>D) { // Two solutions, set flag to 0 and check recoil angle is possible.
120 PSIMX=std::asin(std::sqrt(D/B));
121 if(PSI>PSIMX) return -1; // PSI greater than max possible
122 else iFlag = 0;
123 }
124 DUM2=CPSI+DUM1;
125 DUM3=CPSI-DUM1;
126 E3T=B*DUM2*DUM2;
127 E3TA=B*DUM3*DUM3;
128 E3=ET*B*DUM2*DUM2;
129 E3A=ET*B*DUM3*DUM3;
130 E4=ET-E3;
131 E4A=ET-E3A;
132 // CALCULATE ZETA
133 SZETA=std::sqrt(A3*E3/(A4*E4))*SPSI;
134 if(E4A>0) SZA=std::sqrt(A3*E3A/(A4*E4A))*SPSI;
135 if(E4A<0) SZA=0;
136 ZETA=std::asin(SZETA);
137 ZETAA=std::asin(SZA);
138 // CALCULATE THETA
139 CTHETA=(E3T-B-D)/(2.*std::sqrt(A*C));
140 CTHETAA=(E3TA-B-D)/(2.*std::sqrt(A*C));
141 if(CTHETA>1.) CTHETA=1.;
142 if(CTHETAA>1.) CTHETAA=1.;
143 if(CTHETA<-1.) CTHETA=-1.;
144 if(CTHETAA<-1.) CTHETAA=-1.;
145 THETA=std::acos(CTHETA);
146 THETAA=std::acos(CTHETAA);
147 // CALCULATE PHI
148 PHI=TMath::Pi()-THETA;
149 PHIA=TMath::Pi()-THETAA;
150
151 // OUTPUT PARAMETERS
152 E3I=E3;
153 E3AI=E3A;
154 E4I=E4;
155 E4AI=E4A;
156 ZETAI=ZETA;
157 ZETAAI=ZETAA;
158 thetaI=THETA;
159 thetaaI=THETAA;
160 PHII=PHI;
161 PHIAI=PHIA;
162
163 return iFlag;
164 }
165
166 // :: DEPRECIATED See Eloss.cc for new routine :: Simplified interface to the fortran-based eloss routine by A.E. Stuchbery based on Ziegler's stopping powers. Accepts a monoelemental target only. Dfoil can be 0 and will pull a value from SCOEF.dat.
167 /*void Eloss(float Ein, float Abeam, int Zbeam, float Afoil, int Zfoil, float &Dfoil, float FoilThickness, float &Eout, float &transittime, float &range)
168 {
169 int strips = round(FoilThickness * 100 / Dfoil);
170 float racc = 0.001;
171 float eacc = 0.001;
172 int Zf[10];
173 Zf[0] = Zfoil;
174 float Af[10];
175 Af[0] = Afoil;
176 int nsp = 1; // Number of species = 1
177 float form[10];
178 form[0] = 1.; // Molecular forula = 1
179 eloss_(Ein,Abeam,Zbeam,Af,Zf,nsp,form,Dfoil,FoilThickness,strips,racc,eacc,Eout,transittime,range);
180 }*/
181
182 double xcmlr(double anglab, double angCM, double Ap, double At, double Ep, double Q, bool Tgtex)
183 {
184 /* computes solid angle correction for cm to lab transformation
185 c XCLR = (cm cross section)/(lab cross section) (now return value)
186 c Tlab = lab scattering angle
187 c Tcm = centre of mass scattering angle
188 c
189 c Added extra April 2006 to get case where 180 degrees CM correct
190 c
191 c
192 c Ap = projectile mass
193 c At = target mass
194 c Ep = beam energy
195 c Q = Q value in same units as beam energy
196 c pt = p or t character to specify if it is the target or projectile
197 c solid angle that is required
198 c
199 c Reference: Marion and Young Andrew Stuchbery May 2004
200 converted to 'C' by Brendan McCormick 01/2019
201 c
202 c ******************** angles in radians *******************/
203 double tau, taubar, Etp;
204
205 if(angCM==0 && !Tgtex) return 1.;
206 if(anglab==0) {
207 // CM angle for projectile must be 180 degrees. For target it can be zero:
208 // See AW 1975 p 267 for formulas
209 if(Tgtex) {
210 Etp = Ep + Q * (1. + Ap/At);
211 taubar = std::sqrt(Ep/Etp);
212 return 1. / std::pow(1.+taubar,2);
213 }
214 else {
215 Etp = Ep + Q * (1. + Ap/At);
216 tau = (Ap/At) * std::sqrt(Ep/Etp);
217 return 1. / std::pow(1.-tau,2);
218 }
219 }
220 return std::pow(std::sin(anglab)/std::sin(angCM),2) * std::fabs(std::cos(angCM-anglab));
221 }
222
223 /* Creates the data file which FREADER will read to setup the common block variables for the coulex calculation.
224 * Total states is the total number of nuclear states to be considered in all level schemes, Z_P is projectile Z, A_P is projectile mass, Z_T is target Z, A_T is target mass, nucgroups is the number of different nuclear level schemes,
225 * grouptype[nucgroups] specifies whether the group is a projectile(1) or target(2), nuclvls[nucgroups] specifies how many levels are in the level scheme provided for a scheme, nuclvl[nucgroups][nuclvls] is a struct containing level information,
226 * nucEMtrans[nucgroups] is the number of reduced matrix elements between transitions provided, nucEMtran[nucgroups][nucE/Mtrans] is a structure containing sqrt(B(EM/MM)) values, and filename is an optional output name for the data file.
227 * Function returns 0 on success, 1 on a file access error.
228 */
229 int MakeDatafile(int Zbeam, float Abeam, int Ztgt, float Atgt, int grouptype, std::vector<WDB_nuclvl> &nuclvls, std::vector<WDB_nuctrans> &nucEtrans, std::vector<WDB_nuctrans> &nucMtrans, std::string filename /*= ""*/)
230 {
231 int nlvls = nuclvls.size();
232 int nEtrans = nucEtrans.size();
233 int nMtrans = nucMtrans.size();
234
235 std::string strout;
236 std::ofstream datafile;
237
238 if(filename.empty()) filename = "coulex.dat";
239 datafile.open(filename.c_str());
240 if(datafile.fail()) return 1;
241
242 strout = "1 " + std::to_string(nlvls) +"\n";
243 datafile.write(strout.c_str(),strout.size());
244 strout = "7 0\n";
245 datafile.write(strout.c_str(),strout.size());
246 strout = "8 0\n";
247 datafile.write(strout.c_str(),strout.size());
248 strout = "12 1\n";
249 datafile.write(strout.c_str(),strout.size());
250 strout = "17 " + std::to_string(Zbeam) + ' ' + std::to_string(Abeam) + '\n';
251 datafile.write(strout.c_str(),strout.size());
252 strout = "18 " + std::to_string(Ztgt) + ' ' + std::to_string(Atgt) + '\n';
253 datafile.write(strout.c_str(),strout.size());
254 strout = "19 1.0 \n";
255 datafile.write(strout.c_str(),strout.size());
256 strout = "21 1.0 \n";
257 datafile.write(strout.c_str(),strout.size());
258 for(int l=0;l<nlvls;l++) {
259 if(nuclvls[l].par) strout = "22 " + std::to_string(l+1) + ' ' + std::to_string(nuclvls[l].spin) + ' ' + std::to_string(nuclvls[l].E) + " +1 " + std::to_string(nuclvls[l].K_Band) + '\n';
260 else strout = "22 " + std::to_string(l+1) + ' ' + std::to_string(nuclvls[l].spin) + ' ' + std::to_string(nuclvls[l].E) + " -1 " + std::to_string(nuclvls[l].K_Band) + '\n';
261 datafile.write(strout.c_str(),strout.size());
262 }
263 for(int tE=0;tE<nEtrans;tE++) {
264 strout = "23 " + std::to_string(nucEtrans[tE].lvl1) + ' ' + std::to_string(nucEtrans[tE].lvl2) + ' ' + std::to_string(nucEtrans[tE].B) + ' ' + std::to_string(nucEtrans[tE].M) + '\n';
265 datafile.write(strout.c_str(),strout.size());
266 }
267 for(int tM=0;tM<nMtrans;tM++) {
268 strout = "24 " + std::to_string(nucMtrans[tM].lvl1) + ' ' + std::to_string(nucMtrans[tM].lvl2) + ' ' + std::to_string(nucMtrans[tM].B) + ' ' + std::to_string(nucMtrans[tM].M) + '\n';
269 datafile.write(strout.c_str(),strout.size());
270 }
271 strout = "26 1\n";
272 datafile.write(strout.c_str(),strout.size());
273 //strout = "27 1 1 " + std::to_string(nuclvls) + ' ' + std::to_string(nuclvl[nlvls-1].spin) + '\n';
274 //datafile.write(strout.c_str(),strout.size());
275 strout = "29 1 " + std::to_string(grouptype) + '\n';
276 datafile.write(strout.c_str(),strout.size());
277 strout = "0\n";
278 datafile.write(strout.c_str(),strout.size());
279 strout = "1000\n";
280 datafile.write(strout.c_str(),strout.size());
281
282 datafile.close();
283 return 0;
284 }
285
286 // Calls the coulex setup and data file reading routines that are part of the Winther-DeBoer codes. This only needs to be done once. Use Makedatafile to setup the data file.
287 int ReadDatafile(std::string datafile/*=""*/)
288 {
289 // Error code returned by the freader routine. If non-zero, there was a problem with the data file.
290 int ierr;
291
292 WDB::setup_coulex_();
293
294 // Check if user supplied a data file name and call freader to setup values for coulex routine
295 if(datafile.empty()) datafile = "coulex.dat";
296 char *df = GamR::Utils::c_to_f_str(datafile);
297 WDB::freader_(df,ierr);
298 delete[] df;
299
300 return ierr;
301 }
302
303 // This function is a wrapper to handle datafile construction, reading and removal. Returns 0 on success, 1 on file write error, 2 on read error, and 3 on removal error.
304 int setupcoulex(int Zbeam, float Abeam, int Ztgt, float Atgt, int grouptype, std::vector<WDB_nuclvl> &nuclvls, std::vector<WDB_nuctrans> &nucEtrans, std::vector<WDB_nuctrans> &nucMtrans)
305 {
306 int retval = MakeDatafile(Zbeam,Abeam,Ztgt,Atgt,grouptype,nuclvls,nucEtrans,nucMtrans);
307 if(retval) return retval;
308 retval = ReadDatafile();
309 if(retval) return 2;
310 retval = remove("coulex.dat");
311 if(retval) return 3;
312 return 0;
313 }
314
315 // Gives the Winther-DeBoer codes a projectile energy and target scattering angle (CM) which calculates nuclear stat tensors after coulex, then copies the tensors and reaction cross-section for your chosen level into the variables provided (whether the target or projectile tensors are calculated is determined in the MakeDatafile routine).
316 // Please ensure rho is dimensioned [3][5] (k=0,2,4 q<=0 to 2*k)
317 void GP_AC_Tensors(float Ebeam, float CMAng, int lvl, std::complex<double> **rho, float &xsect)
318 {
319 float rhor[3][5];
320
321 WDB::coulparcm_(Ebeam, CMAng);
322
323 WDB::coulex_();
324
325 int lvl1 = lvl+1; //+1 because FORTRAN doesn't 0 index
326 WDB::getcoulex_(lvl1, rhor, xsect);
327
328 for(int k=0; k<3; k++) {
329 for(int q=0; q<=2*k; q++) {
330 rho[k][q] = std::complex<double>(double(rhor[k][q]/rhor[0][0]), 0.);
331 }
332 }
333 }
334
335 // Draws particle detector configuration which tensors will be calculated for
336 void DrawSetup()
337 {
338 std::cout << " _________________________" << std::endl;
339 std::cout << "| --- Top-down view --- |" << std::endl;
340 std::cout << "| Phi angle rotates |" << std::endl;
341 std::cout << "| particle detector |" << std::endl;
342 std::cout << "| counter-clockwise |" << std::endl;
343 std::cout << "| around beam direction |" << std::endl;
344 std::cout << "| in the x-y plane. |" << std::endl;
345 std::cout << "| Theta angle rotates |" << std::endl;
346 std::cout << "| gamma-ray detector |" << std::endl;
347 std::cout << "| clockwise around tgt |" << std::endl;
348 std::cout << "| in x-z plane. |" << std::endl;
349 std::cout << "|_________________________|" << std::endl;
350 std::cout << "| |" << std::endl;
351 std::cout << "|Beam Tgt stop|" << std::endl;
352 std::cout << "|=======|---------------[]|" << std::endl;
353 std::cout << "| < z >^ |" << std::endl;
354 std::cout << "| r |" << std::endl;
355 std::cout << "| v |" << std::endl;
356 std::cout << "| x |\\y |" << std::endl;
357 std::cout << "| | Ptcle | \\ |" << std::endl;
358 std::cout << "| | dtctr | | |" << std::endl;
359 std::cout << "| ----z \\ |x |" << std::endl;
360 std::cout << "| \\ \\| |" << std::endl;
361 std::cout << "| \\ |" << std::endl;
362 std::cout << "| y |" << std::endl;
363 std::cout << "|_________________________|" << std::endl;
364 }
365
366
367
368 /*
369 //NOTE pass Qks vector by reference, nullptr is acceptable and results in no return value.
370 void ExperimentalSetup(int &Zbeam, double &Abeam, int &Ztgt, double &Atgt, bool &ProjectileDetected, bool &TargetExcited, double &Ebeam, double &Qex, double &TargetThickness, double &TargetDensity, double &PDetx, double &PDety, double &PDetz, double &PDetr, int &ngamma, int &nstates, std::vector<WDB_nuclvl> &NucLevels, int &Tenslevel, int &nEtrans, std::vector<WDB_nuctrans> &Etrans, int &nMtrans, std::vector<WDB_nuctrans> &Mtrans, std::vector<GamR::AngDist::SolidAttenuation> *Qks )
371 {
372
373 // Beam and target (charges, masses in AMU, beam energy in MeV and target thickness in mg/cm2, excited-state energy in MeV, excited-state spin, ground-state spin, B(E2) in eb, specify if target excites [or projectile], specify if projectile is detected [or target], excited-state partiy and ground-state parity)
374 int ttrans;
375 float *SE;
376 float *SS;
377 bool *SP;
378 float *B;
379 int *Bl1;
380 int *Bl2;
381 int *M;
382 //Gamma-ray detector info for Qk calculation
383 float *gL, *gR, *gD, *gr;
384 int trans1;
385 // Input handling
386 std::string inp;
387 std::vector<std::string> allinp;
388
389 // Answer save file
390 std::ifstream ansfile;
391 std::string ans;
392 std::string ansstr = "Wtheta.ans";
393
394 // Get all the inputs necessary to calculate the stat tensors
395 std::cout << "This routine will assist you to setup the variables required for the CoulexNucTensors routine for a given target and detector geometry." << std::endl << std::endl;
396 retryansfile:
397 std::cout << "Answer file name [" << ansstr << "]: ";
398 std::getline(std::cin, inp);
399 if(ansstr.empty() && inp.empty()) ansstr = "Wtheta.ans";
400 else {
401 if(!inp.empty()) ansstr = inp;
402 ansfile.open(ansstr.c_str());
403 if(ansfile.fail()) {
404 std::cout << "Could not open answer file! Please try again or leave blank to skip." << std::endl;
405 ansstr = "";
406 goto retryansfile;
407 }
408 }
409 std::cout << "Beam and target details -" << std::endl;
410 ans = readans(ansfile);
411 retryBeam:
412 std::cout << "\tBeam nucleus [" << ans << "]: ";
413 std::getline(std::cin,inp);
414 if(inp.empty()) inp = ans;
415 if(GamR::TK::GetZA(inp,Zbeam,Abeam)) {
416 std::cout << "\tInvalid nucleus. Ensure nucleus is in the format AEl, ElA or El-A e.g. 3He, Li-5." << std::endl;
417 goto retryBeam;
418 }
419 allinp.push_back(inp);
420 ans = readans(ansfile);
421 retryE:
422 std::cout << "\tBeam energy in MeV [" << ans << "]: ";
423 std::getline(std::cin,inp);
424 if(inp.empty()) inp = ans;
425 if(catcherr(inp, Ebeam)) goto retryE;
426 allinp.push_back(inp);
427 ans = readans(ansfile);
428 retryTgt:
429 std::cout << "\tTarget nucleus [" << ans << "]: ";
430 std::getline(std::cin,inp);
431 if(inp.empty()) inp = ans;
432 if(GamR::TK::GetZA(inp,Ztgt,Atgt)) {
433 std::cout << "\tInvalid nucleus. Ensure nucleus is in the format AEl, ElA or El-A e.g. 3He, Li5, C-13." << std::endl;
434 goto retryTgt;
435 }
436 allinp.push_back(inp);
437 ans = readans(ansfile);
438 retryTargetThickness:
439 std::cout << "\tTarget thickness in mg/cm2 [" << ans << "]: ";
440 std::getline(std::cin,inp);
441 if(inp.empty()) inp = ans;
442 if(catcherr(inp, TargetThickness)) goto retryTargetThickness;
443 allinp.push_back(inp);
444 ans = readans(ansfile);
445 retryTargetDensity:
446 std::cout << "\tTarget density in g/cm3 (leave blank to get from SCOEF.DAT) [" << ans << "]: ";
447 std::getline(std::cin,inp);
448 if(inp.empty()) inp = ans;
449 if(inp.empty()) {
450 TargetDensity = 0.;
451 inp = "0";
452 }
453 else {
454 if(catcherr(inp, TargetDensity, true)) goto retryTargetDensity;
455 }
456 allinp.push_back(inp);
457 ans = readans(ansfile);
458 retryEx:
459 std::cout << "\tProjectile(1) or target(2) excitation [" << ans << "]: ";
460 std::getline(std::cin,inp);
461 if(inp.empty()) inp = ans;
462 switch(inp[0]) {
463 case '1':
464 TargetExcited = false;
465 break;
466 case '2':
467 TargetExcited = true;
468 break;
469 default:
470 goto retryEx;
471 }
472 allinp.push_back(inp);
473 ans = readans(ansfile);
474 retryDet:
475 std::cout << "\tProjectile(1) or target(2) detection [" << ans << "]: ";
476 std::getline(std::cin,inp);
477 if(inp.empty()) inp = ans;
478 switch(inp[0]) {
479 case '1':
480 ProjectileDetected = true;
481 break;
482 case '2':
483 ProjectileDetected = false;
484 break;
485 default:
486 goto retryDet;
487 }
488 allinp.push_back(inp);
489 ans = readans(ansfile);
490 int pread = 0;
491 if(!ans.empty()) pread = std::stoi(ans);
492 retryStates:
493 std::cout << "\tTotal states including ground state [" << ans << "]: ";
494 std::getline(std::cin,inp);
495 if(inp.empty()) inp = ans;
496 if(catcherr(inp,nstates)) goto retryStates;
497 if(nstates<2) {
498 std::cout << "\tInvalid entry. Must be >1." << std::endl;
499 goto retryStates;
500 }
501 allinp.push_back(inp);
502 SE = new float[nstates];
503 SS = new float[nstates];
504 SP = new bool[nstates];
505 for(int i=0; i<nstates; i++) {
506 if(pread) ans = readans(ansfile);
507 retrySE:
508 std::cout << "\tState " << i+1 << " energy in MeV [" << ans << "]: ";
509 std::getline(std::cin,inp);
510 if(inp.empty()) inp = ans;
511 try { SE[i] = stof(inp); }
512 catch(const std::invalid_argument& ia) {
513 std::cout << "Invalid entry. Must be numerical." << std::endl;
514 goto retrySE;
515 }
516 if(SE[i]<0) {
517 std::cout << "Invalid entry. Must be >=0." << std::endl;
518 goto retrySE;
519 }
520 allinp.push_back(inp);
521 if(pread) {
522 ans = readans(ansfile);
523 pread--;
524 }
525 retrySS:
526 std::cout << "\tState " << i+1 << " spin and parity [" << ans << "]: ";
527 std::getline(std::cin,inp);
528 if(inp.empty()) inp = ans;
529 if(inp[inp.size()-1]=='+') SP[i] = true;
530 else if(inp[inp.size()-1]=='-') SP[i] = false;
531 else {
532 std::cout << "\tParity not specified. Please specify spin and parity e.g. 0.5+" << std::endl;
533 goto retrySS;
534 }
535 try { SS[i] = std::stof(inp.substr(0,inp.size()-1)); }
536 catch(const std::invalid_argument& ia) {
537 std::cout << "\tInvalid entry. Must be numerical." << std::endl;
538 allinp.pop_back();
539 goto retrySS;
540 }
541 if(SS[i]<0) {
542 std::cout << "\tInvalid entry. Must be >=0." << std::endl;
543 allinp.pop_back();
544 goto retrySS;
545 }
546 allinp.push_back(inp);
547 }
548 if(pread) {
549 for(int i=0; i<pread; i++) {
550 ans = readans(ansfile);
551 ans = readans(ansfile);
552 }
553 }
554 ans = readans(ansfile);
555 retryTenslevel:
556 std::cout << "\tCalculate tensors for state [" << ans << "]: ";
557 std::getline(std::cin,inp);
558 if(inp.empty()) inp = ans;
559 if(catcherr(inp, Tenslevel)) goto retryTenslevel;
560 allinp.push_back(inp);
561 ans = readans(ansfile);
562 if(!ans.empty()) pread = std::stoi(ans);
563 else pread = 0;
564 retryTrans:
565 std::cout << "\tTotal transitions between states [" << ans << "]: ";
566 std::getline(std::cin,inp);
567 if(inp.empty()) inp = ans;
568 if(catcherr(inp,ttrans)) goto retryTrans;
569 allinp.push_back(inp);
570 B = new float[ttrans];
571 Bl1 = new int[ttrans];
572 Bl2 = new int[ttrans];
573 M = new int[ttrans];
574 for(int i=0; i<ttrans; i++) {
575 if(pread) ans = readans(ansfile);
576 retryl1:
577 std::cout << "\tTransition " << i+1 << " lower state number [" << ans << "]: ";
578 std::getline(std::cin,inp);
579 if(inp.empty()) inp = ans;
580 if(catcherr(inp, Bl1[i])) goto retryl1;
581 allinp.push_back(inp);
582 if(pread) ans = readans(ansfile);
583 retryl2:
584 std::cout << "\tTransition " << i+1 << " upper state number [" << ans << "]: ";
585 std::getline(std::cin,inp);
586 if(inp.empty()) inp = ans;
587 if(catcherr(inp, Bl2[i])) goto retryl2;
588 allinp.push_back(inp);
589 if(pread) ans = readans(ansfile);
590 retryB:
591 std::cout << "\tTransition " << i+1 << " B up in eb [" << ans << "]: ";
592 std::getline(std::cin,inp);
593 if(inp.empty()) inp = ans;
594 if(catcherr(inp, B[i])) goto retryB;
595 allinp.push_back(inp);
596 if(pread) {
597 ans = readans(ansfile);
598 pread--;
599 }
600 retryM:
601 std::cout << "\tTransition " << i+1 << " multipolarity [" << ans << "]: ";
602 std::getline(std::cin,inp);
603 if(inp.empty()) inp = ans;
604 if(catcherr(inp, M[i])) goto retryM;
605 allinp.push_back(inp);
606 }
607 if(pread) {
608 for(int i=0; i<pread; i++) {
609 ans = readans(ansfile);
610 ans = readans(ansfile);
611 ans = readans(ansfile);
612 ans = readans(ansfile);
613 }
614 }
615 ans = readans(ansfile);
616 retrytrans1:
617 std::cout << "\tObserved transition [" << ans << "]: ";
618 std::getline(std::cin,inp);
619 if(inp.empty()) inp = ans;
620 if(catcherr(inp, trans1)) goto retrytrans1;
621 allinp.push_back(inp);
622
623 // ==================================================================================================================
624 // ------------------------------------------------------------------------------------------------------------------
625 // ==================================================================================================================
626
627 std::cout << "Enter details for detector layout shown below. Rotate tensors later for phi != 0" << std::endl;
628 std::cout << "NOTE! Enter a negative z for back-scattered particle detection." << std::endl;
629 DrawSetup();
630 std::cout << std::endl << "Particle detector details -" << std::endl;
631 ans = readans(ansfile);
632 retryPDetx:
633 std::cout << "\tDetector width x in mm [" << ans << "]: ";
634 std::getline(std::cin,inp);
635 if(inp.empty()) inp = ans;
636 if(catcherr(inp, PDetx)) goto retryPDetx;
637 allinp.push_back(inp);
638 ans = readans(ansfile);
639 retryPDety:
640 std::cout << "\tDetector height y in mm [" << ans << "]: ";
641 std::getline(std::cin,inp);
642 if(inp.empty()) inp = ans;
643 if(catcherr(inp, PDety)) goto retryPDety;
644 allinp.push_back(inp);
645 ans = readans(ansfile);
646 retryPDetz:
647 std::cout << "\tHorizontal distance from target to detector z in mm [" << ans << "]: ";
648 std::getline(std::cin,inp);
649 if(inp.empty()) inp = ans;
650 try { PDetz = stof(inp); }
651 catch(const std::invalid_argument& ia) {
652 std::cout << "\tInvalid distance. Must be numerical." << std::endl;
653 goto retryPDetz;
654 }
655 if(!PDetz) {
656 std::cout << "\tInvalid distance. Must be non-zero." << std::endl;
657 goto retryPDetz;
658 }
659 allinp.push_back(inp);
660 ans = readans(ansfile);
661 retryPDetr:
662 std::cout << "\tRadial distance from centre to detector edge r in mm [" << ans << "]: ";
663 std::getline(std::cin,inp);
664 if(inp.empty()) inp = ans;
665 try { PDetr = stof(inp); }
666 catch(const std::invalid_argument& ia) {
667 std::cout << "\tInvalid distance. Must be numerical." << std::endl;
668 goto retryPDetr;
669 }
670 allinp.push_back(inp);
671 if(Qks!=nullptr) {
672 std::cout << "Finally, enter details for Ge(Li) gamma-ray detectors." << std::endl;
673 ans = readans(ansfile);
674 retrygtotal:
675 std::cout << "\tTotal number of gamma-ray detectors [" << ans << "]: ";
676 std::getline(std::cin,inp);
677 if(inp.empty()) inp = ans;
678 if(catcherr(inp, ngamma)) goto retrygtotal;
679 allinp.push_back(inp);
680 if(!ans.empty()) pread = stoi(ans);
681 else pread = 0;
682 gD = new float[ngamma];
683 gR = new float[ngamma];
684 gL = new float[ngamma];
685 gr = new float[ngamma];
686 for(int i=0; i<ngamma; i++) {
687 if(pread) ans = readans(ansfile);
688 retrygD:
689 std::cout << "\tGe crystal " << i+1 << " distance from source in mm [" << ans << "]: ";
690 std::getline(std::cin,inp);
691 if(inp.empty()) inp = ans;
692 if(catcherr(inp, gD[i])) goto retrygD;
693 gD[i] *= 0.1;
694 allinp.push_back(inp);
695 if(pread) ans = readans(ansfile);
696 retrygL:
697 std::cout << "\tGe crystal " << i+1 << " length in mm [" << ans << "]: ";
698 std::getline(std::cin,inp);
699 if(inp.empty()) inp = ans;
700 if(catcherr(inp, gL[i])) goto retrygL;
701 gL[i] *= 0.1;
702 allinp.push_back(inp);
703 if(pread) ans = readans(ansfile);
704 retrygR:
705 std::cout << "\tGe crystal " << i+1 << " outer radius in mm [" << ans << "]: ";
706 std::getline(std::cin,inp);
707 if(inp.empty()) inp = ans;
708 if(catcherr(inp, gR[i])) goto retrygR;
709 gR[i] *= 0.1;
710 allinp.push_back(inp);
711 if(pread) {
712 ans = readans(ansfile);
713 pread--;
714 }
715 retrygr:
716 std::cout << "\tGe crystal " << i+1 << " inner radius in mm [" << ans << "]: ";
717 std::getline(std::cin,inp);
718 if(inp.empty()) inp = ans;
719 if(catcherr(inp, gr[i])) goto retrygr;
720 gr[i] *= 0.1;
721 allinp.push_back(inp);
722 }
723 if(pread) {
724 for(int i=0; i<pread; i++) {
725 ans = readans(ansfile);
726 ans = readans(ansfile);
727 ans = readans(ansfile);
728 ans = readans(ansfile);
729 }
730 }
731 }
732
733 // Write user answers to file
734 ansfile.close();
735 std::ofstream ansout;
736 ansout.open(ansstr.c_str());
737 if(!ansout.fail()) {
738 for(auto outp : allinp) ansout << outp << std::endl;
739 }
740 ansout.close();
741
742 for(int i=0; i<nstates; i++) {
743 WDB_nuclvl nlvl;
744 nlvl.spin = SS[i];
745 nlvl.E = SE[i];
746 nlvl.par = SP[i];
747 nlvl.K_Band = 0;
748 NucLevels.push_back(nlvl);
749 }
750
751 nMtrans = 0;
752 nEtrans = 0;
753 for(int i=0; i<ttrans; i++) {
754 int m = 0;
755 int e = 0;
756 bool oddm = M[i] % 2;
757 int l1 = Bl1[i];
758 int l2 = Bl2[i];
759 WDB_nuctrans ntran;
760 ntran.lvl1 = Bl1[i];
761 ntran.lvl2 = Bl2[i];
762 ntran.B = std::sqrt(B[i]);
763 ntran.M = M[i];
764 if(!oddm && SP[l2-1]==SP[l1-1] || oddm && SP[l1-1]!=SP[l2-1]) {
765 Etrans.push_back(ntran);
766 nEtrans++;
767 }
768 else {
769 Mtrans.push_back(ntran);
770 nMtrans++;
771 }
772 }
773
774 // Calculate Qk for each detector
775 if(Qks!=nullptr) {
776 for(int g=0; g<ngamma; g++) {
777 int l2 = Bl2[trans1-1] - 1;
778 int l1 = Bl1[trans1-1] - 1;
779 double Et = NucLevels[l2].E - NucLevels[l1].E;
780 std::vector<double> Qk;
781 for(int k=0; k<=4; k++) {
782 if(k%2) Qk.push_back(0.);
783 else Qk.push_back(getQk(Et, gL[g], gR[g], gr[g], gD[g], double(k)));
784 }
785 GamR::AngDist::SolidAttenuation GamRQk(Qk);
786 Qks->push_back(GamRQk);
787 }
788 }
789
790 delete[] Bl1;
791 delete[] Bl2;
792 delete[] B;
793 delete[] M;
794 if(Qks!=nullptr) {
795 delete[] gL;
796 delete[] gR;
797 delete[] gD;
798 delete[] gr;
799 }
800 }
801*/
802
803 /*
804 // Routine to facilitate stat tensors calculation. Calls ExperimentalSetup and then passes inputs to CoulexNucTensors to return the tensors. Returns a nullptr if there was an error.
805 // Pass getQk by reference, or leave null to not get Qks.
806 GamR::AngDist::StatTensor *DoCoulex(std::vector<GamR::AngDist::SolidAttenuation> *getQk = nullptr)
807 {
808 int Zp, Zt, ngam, nstates, tenslvl, nEt, nMt;
809 double Ap, At, E, Q, Tt, Td, Px, Py, Pz, Pr;
810 std::vector<WDB_nuclvl> nlvl;
811 std::vector<WDB_nuctrans> Etrans, Mtrans;
812 bool ProjDet, TgtEx;
813 if(getQk==nullptr) ExperimentalSetup(Zp, Ap, Zt, At, ProjDet, TgtEx, E, Q, Tt, Td, Px, Py, Pz, Pr, ngam, nstates, nlvl, tenslvl, nEt, Etrans, nMt, Mtrans);
814 else ExperimentalSetup(Zp, Ap, Zt, At, ProjDet, TgtEx, E, Q, Tt, Td, Px, Py, Pz, Pr, ngam, nstates, nlvl, tenslvl, nEt, Etrans, nMt, Mtrans, getQk);
815 return CoulexNucTensors(Zp, Ap, Zt, At, ProjDet, TgtEx, E, Q, Td, Tt, tenslvl-1, nstates, nlvl, nEt, Etrans, nMt, Mtrans, Px, Py, Pz, Pr);
816 }
817 */
818 }
819}
GamR::Coulex::WDB::coulparcm_
void coulparcm_(float &eb, float &tcmdeg)
GamR::Coulex::WDB::freader_
void freader_(char readfile[], int &inperr)
GamR::Coulex::WDB::getcoulex_
void getcoulex_(int &lvl, float rho_cre[3][5], float &xsect)
GamR::Coulex::ReadDatafile
int ReadDatafile(std::string datafile)
Definition Coulex.cc:287
GamR::Coulex::GP_AC_Tensors
void GP_AC_Tensors(float Ebeam, float CMAng, int lvl, std::complex< double > **rho, float &xsect)
Definition Coulex.cc:317
GamR::Coulex::DrawSetup
void DrawSetup()
Definition Coulex.cc:336
GamR::Coulex::Kinematics
int Kinematics(double A1I, double A2I, double A3I, double A4I, double E1I, double QEFFI, double PSII, double &E3I, double &E3AI, double &E4I, double &E4AI, double &ZETAI, double &ZETAAI, double &thetaI, double &thetaaI, double &PHII, double &PHIAI)
Definition Coulex.cc:50
GamR::Coulex::setupcoulex
int setupcoulex(int Zbeam, float Abeam, int Ztgt, float Atgt, int grouptype, std::vector< WDB_nuclvl > &nuclvls, std::vector< WDB_nuctrans > &nucEtrans, std::vector< WDB_nuctrans > &nucMtrans)
Definition Coulex.cc:304
GamR::Coulex::xcmlr
double xcmlr(double anglab, double angCM, double Ap, double At, double Ep, double Q, bool Tgtex)
Definition Coulex.cc:182
GamR::Coulex::MakeDatafile
int MakeDatafile(int Zbeam, float Abeam, int Ztgt, float Atgt, int grouptype, std::vector< WDB_nuclvl > &nuclvls, std::vector< WDB_nuctrans > &nucEtrans, std::vector< WDB_nuctrans > &nucMtrans, std::string filename)
Definition Coulex.cc:229
GamR::Utils::c_to_f_str
char * c_to_f_str(std::string strin)
Definition Utilities.cc:284
GamR
Definition Gain.cc:19