20 void Usage(JApplication &app);
52 string mkstr(
const T &val,
unsigned int width)
56 string s = ssval.str();
58 int Nspaces = (int)width - (
int)s.size();
59 if(Nspaces>0)ss<<
string(Nspaces,
' ');
69 int main(
int narg,
char *argv[])
79 cout<<
"Filling material table ..."; cout.flush();
87 for(
int ir=0; ir<
Nr; ir++){
88 double r = r0 + dr/2.0 + (double)ir*dr;
89 MatTable[ir]=&buff[ir*
Nz];
90 for(
int iz=0; iz<
Nz; iz++){
91 double z = z0 + dz/2.0 + (double)iz*dz;
94 double d_r = dr/(double)
n_r;
95 double d_z = dz/(double)
n_z;
96 double d_phi = 2.0*M_PI/(double)
n_phi;
97 Material avg_mat={0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
98 for(
int i_r=0; i_r<
n_r; i_r++){
99 double my_r = r - dr/2.0 + (double)i_r*d_r;
100 for(
int i_z=0; i_z<
n_z; i_z++){
101 double my_z = z - dz/2.0 + (double)i_z*d_z;
102 for(
int i_phi=0; i_phi<
n_phi; i_phi++){
103 double my_phi = (double)i_phi*d_phi;
105 DVector3 pos(my_r*cos(my_phi), my_r*
sin(my_phi), my_z);
106 double A, Z, density, radlen, rhoZ_overA, rhoZ_overA_logI;
107 jerror_t err = g->
FindMatLL(pos, density, A, Z, radlen);
109 rhoZ_overA = density*Z/A;
110 double I = (Z*12.0 + 7.0)*1.0E-9;
111 rhoZ_overA_logI = rhoZ_overA*log(I);
113 A = Z = density = radlen = rhoZ_overA = rhoZ_overA_logI = 0.0;
119 avg_mat.
RadLen += 1.0/radlen;
143 MatTable[ir][iz] = avg_mat;
145 cout<<
"\r Filling Material table ... "<<100.0*(double)ir/(
double)Nr<<
"% ";cout.flush();
155 cols.push_back(
"density");
156 cols.push_back(
"radlen");
157 cols.push_back(
"rhoZ_overA");
158 cols.push_back(
"rhoZ_overA_logI");
159 cols.push_back(
"chi2c_factor");
160 cols.push_back(
"chi2a_factor");
161 cols.push_back(
"chi2a_corr");
164 vector<unsigned int> max_width;
165 for(
unsigned int i=0; i<cols.size(); i++)max_width.push_back(cols[i].size()+2);
166 for(
int ir=0; ir<
Nr; ir++){
167 double r = r0 + dr/2.0 + (double)ir*dr;
168 for(
int iz=0; iz<
Nz; iz++){
169 double z = z0 + dz/2.0 + (double)iz*dz;
174 ss.str(
""); ss<<r; strs.push_back(ss.str());
175 ss.str(
""); ss<<z; strs.push_back(ss.str());
176 ss.str(
""); ss<<mat.
A; strs.push_back(ss.str());
177 ss.str(
""); ss<<mat.
Z; strs.push_back(ss.str());
178 ss.str(
""); ss<<mat.
Density; strs.push_back(ss.str());
179 ss.str(
""); ss<<mat.
RadLen; strs.push_back(ss.str());
180 ss.str(
""); ss<<mat.
rhoZ_overA; strs.push_back(ss.str());
182 ss.str(
""); ss<<mat.
chi2c_factor; strs.push_back(ss.str());
183 ss.str(
""); ss<<mat.
chi2a_factor; strs.push_back(ss.str());
184 ss.str(
""); ss<<mat.
chi2a_corr; strs.push_back(ss.str());
186 for(
unsigned int i=0; i<strs.size(); i++){
187 if(strs[i].
size()+2 > max_width[i])max_width[i] = strs[i].
size()+2;
194 stringstream colheader;
196 for(
unsigned int i=0; i<cols.size(); i++)colheader<<
mkstr(cols[i], max_width[i] + (i==0 ? -1:0));
200 cout<<
"Writing material table to file ..."<<endl;
201 ofstream of(
"material_map");
202 time_t t = time(NULL);
204 of<<
"# Material map generated with src/programs/Utilities/mkMaterialMap"<<endl;
205 of<<
"# generated: "<<ctime(&t);
207 of<<
"# Generated with the following parameters:"<<endl;
208 of<<
"# Nr = "<<Nr<<endl;
209 of<<
"# Nz = "<<
Nz<<endl;
210 of<<
"# rmin = "<<
rmin<<endl;
211 of<<
"# rmax = "<<
rmax<<endl;
212 of<<
"# zmin = "<<
zmin<<endl;
213 of<<
"# zmax = "<<
zmax<<endl;
215 of<<
"# sampling points per cell:"<<endl;
216 of<<
"# n_r = "<<
n_r<<endl;
217 of<<
"# n_z = "<<
n_z<<endl;
218 of<<
"# n_phi = "<<
n_phi<<endl;
223 unsigned int entries_written = 0;
224 for(
int ir=0; ir<
Nr; ir++){
225 double r = r0 + dr/2.0 + (double)ir*dr;
226 for(
int iz=0; iz<
Nz; iz++){
227 double z = z0 + dz/2.0 + (double)iz*dz;
231 if(entries_written%50 == 0)of<<colheader.str()<<endl;
232 if(entries_written==0) of <<
"#% 00 01 02 03 04 05 06 07 08 09 10" <<endl;
234 of<<
mkstr(r, max_width[0]);
235 of<<
mkstr(z, max_width[1]);
236 of<<
mkstr(mat.
A, max_width[2]);
237 of<<
mkstr(mat.
Z, max_width[3]);
264 for(
int i=1; i<narg; i++){
265 string arg = argv[i];
266 string arg2 = ((i+1)<narg) ? argv[i+1]:
"";
271 else if(arg==
"-Nr" ){
Nr = atoi(arg2.c_str()); i++;}
272 else if(arg==
"-Nz" ){
Nz = atoi(arg2.c_str()); i++;}
273 else if(arg==
"-rmin" ){
rmin = atof(arg2.c_str()); i++;}
274 else if(arg==
"-rmax" ){
rmax = atof(arg2.c_str()); i++;}
275 else if(arg==
"-zmin" ){
zmin = atof(arg2.c_str()); i++;}
276 else if(arg==
"-zmax" ){
zmax = atof(arg2.c_str()); i++;}
277 else if(arg==
"-n_r" ){
n_r = atoi(arg2.c_str()); i++;}
278 else if(arg==
"-n_z" ){
n_z = atoi(arg2.c_str()); i++;}
279 else if(arg==
"-n_phi"){
n_phi = atoi(arg2.c_str()); i++;}
281 cout<<
"Unknown argument \""<<arg<<
"\""<<endl;
285 cout<<
"option \""<<arg<<
"\" requires a value!"<<endl;
294 void Usage(JApplication &app)
297 cout<<
"Usage:"<<endl;
298 cout<<
" mkMaterialMap [options] source1 source2 source3 ..."<<endl;
300 cout<<
" -Nr # Set number of grid points in r (def:"<<
Nr<<
")"<<endl;
301 cout<<
" -Nz # Set number of grid points in z (def:"<<
Nz<<
")"<<endl;
302 cout<<
" -rmin # Set lower boundary of map in r (def:"<<
rmin<<
")"<<endl;
303 cout<<
" -rmax # Set upper boundary of map in r (def:"<<
rmax<<
")"<<endl;
304 cout<<
" -zmin # Set lower boundary of map in z (def:"<<
zmin<<
")"<<endl;
305 cout<<
" -zmax # Set upper boundary of map in z (def:"<<
zmax<<
")"<<endl;
306 cout<<
" -n_r # Set number of sampling points in r (def:"<<
n_r<<
")"<<endl;
307 cout<<
" -n_z # Set number of sampling points in z (def:"<<
n_r<<
")"<<endl;
308 cout<<
" -n_phi # Set number of sampling points in phi (def:"<<
n_phi<<
")"<<endl;
string mkstr(const T &val, unsigned int width)
jerror_t FindMatLL(DVector3 pos, double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const
void ParseCommandLineArguments(int &narg, char *argv[])
void Usage(JApplication &app)
int main(int argc, char *argv[])