19 pthread_mutexattr_init(&mutex_attr);
20 pthread_mutexattr_settype(&mutex_attr, PTHREAD_MUTEX_RECURSIVE);
21 pthread_mutex_init(&mutex, &mutex_attr);
23 table_initialized =
false;
26 jcalib = japp->GetJCalibration(run_number);
28 _DBG_<<
"Unable to get JCalibration object!"<<endl;
32 JParameterManager *jparms = japp->GetJParameterManager();
34 _DBG_<<
"Unable to get JParameterManager object!"<<endl;
65 cout<<
"Reading Material map from "<<namepath<<
" ..."<<endl;
66 vector< vector<float> > Mmap;
67 jcalib->Get(namepath, Mmap);
68 cout<<Mmap.size()<<
" entries found (";
81 map<float, int> rvals;
82 map<float, int> zvals;
86 for(
unsigned int i=0; i<Mmap.size(); i++){
87 vector<float> &a = Mmap[i];
102 dr = (rmax-
rmin)/(
double)(
Nr-1);
103 dz = (zmax-
zmin)/(
double)(
Nz-1);
106 cout<<
" ) at 0x"<<hex<<(
unsigned long)
this<<dec<<endl;
111 for(
int ir=0; ir<
Nr; ir++){
112 MatTable[ir]=&buff[ir*
Nz];
116 for(
unsigned int i=0; i<Mmap.size(); i++){
117 vector<float> &a = Mmap[i];
120 int ir = (int)floor((r - r0)/dr);
121 int iz = (int)floor((z - z0)/dz);
122 if(ir<0 || ir>=Nr){
_DBG_<<
"ir out of range: ir="<<ir<<
" Nr="<<Nr<<endl;
continue;}
123 if(iz<0 || iz>=Nz){
_DBG_<<
"iz out of range: iz="<<iz<<
" Nz="<<Nz<<endl;
continue;}
124 VolMat &mat = MatTable[ir][iz];
142 pthread_mutex_lock(&mutex);
144 if(table_initialized){
146 pthread_mutex_unlock(&mutex);
151 cout<<
"Filling material table"; cout.flush();
156 double zmin = -100.0;
160 dr = (rmax-
rmin)/(
double)
Nr;
161 dz = (zmax-
zmin)/(
double)
Nz;
164 for(
int ir=0; ir<
Nr; ir++){
165 double r = r0 + (double)ir*dr;
166 MatTable[ir]=&buff[ir*
Nz];
167 for(
int iz=0; iz<
Nz; iz++){
168 double z = z0 + (double)iz*dz;
174 double d_r = dr/(double)n_r;
175 double d_z = dz/(double)n_z;
176 double d_phi = 2.0*M_PI/(double)n_phi;
177 VolMat avg_mat={0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
179 for(
int i_r=0; i_r<
n_r; i_r++){
180 double my_r = r - dr/2.0 + (double)i_r*d_r;
181 for(
int i_z=0; i_z<
n_z; i_z++){
182 double my_z = z - dz/2.0 + (double)i_z*d_z;
183 for(
int i_phi=0; i_phi<
n_phi; i_phi++){
184 double my_phi = (double)i_phi*d_phi;
186 DVector3 pos(my_r*cos(my_phi), my_r*
sin(my_phi), my_z);
187 double A, Z, density, radlen,LnI;
189 FindMatLL(pos, density, A, Z, radlen,LnI);
191 double rhoZ_overA = density*Z/A;
192 double rhoZ_overA_logI = rhoZ_overA*LnI;
205 avg_mat.
A /= (double)(n_r*n_z*n_phi);
206 avg_mat.
Z /= (double)(n_r*n_z*n_phi);
207 avg_mat.
Density /= (double)(n_r*n_z*n_phi);
208 avg_mat.
RadLen /= (double)(n_r*n_z*n_phi);
209 avg_mat.
rhoZ_overA /= (double)(n_r*n_z*n_phi);
212 MatTable[ir][iz] = avg_mat;
214 cout<<
"\r Filling Material table ... "<<100.0*(double)ir/(
double)Nr<<
"% ";cout.flush();
218 table_initialized=
true;
219 pthread_mutex_unlock(&mutex);
234 #define ROOT_MINOR 28
237 #if ROOT_MAJOR>=5 && ROOT_MINOR>=28
239 cout<<
"Created TGeoManager :"<<gGeoManager<<endl;
241 cout<<
"Skipping explicit TGeoManager creation due to ROOT ver. < 5.28"<<endl;
252 pthread_mutex_lock(&mutex);
254 if(!DRGeom)InitDRGeom();
256 TGeoNode *cnode = DRGeom->FindNode(x[0],x[1],x[2]);
258 pthread_mutex_unlock(&mutex);
270 pthread_mutex_lock(&mutex);
271 if(!DRGeom)InitDRGeom();
272 TGeoNode *cnode = DRGeom->FindNode(x[0],x[1],x[2]);
273 TGeoVolume *cvol = cnode->GetVolume();
274 pthread_mutex_unlock(&mutex);
285 return FindMatLL(pos, density, A, Z, RadLen);
293 return FindMatTable(pos, rhoZ_overA, rhoZ_overA_logI, RadLen);
301 if(!table_initialized)((
DRootGeom*)
this)->InitTable();
305 double r = pos.Perp();
307 int ir = (int)floor((r-r0)/dr);
308 int iz = (int)floor((z-z0)/dz);
309 if(ir<0 || ir>=
Nr || iz<0 || iz>=
Nz){
310 A = Z = density = RadLen = 0.0;
311 return RESOURCE_UNAVAILABLE;
314 const VolMat &mat = MatTable[ir][iz];
329 if(!table_initialized)((
DRootGeom*)
this)->InitTable();
333 double r = pos.Perp();
335 int ir = (int)floor((r-r0)/dr);
336 int iz = (int)floor((z-z0)/dz);
337 if(ir<0 || ir>=
Nr || iz<0 || iz>=
Nz){
338 rhoZ_overA = rhoZ_overA_logI = RadLen = 0.0;
339 return RESOURCE_UNAVAILABLE;
342 const VolMat &mat = MatTable[ir][iz];
352 double &rhoZ_overA_logI,
double &RadLen)
const{
353 pthread_mutex_lock(const_cast<pthread_mutex_t*>(&mutex));
355 if(!DRGeom)((
DRootGeom*)
this)->InitDRGeom();
357 TGeoMaterial *mat=DRGeom->GetMaterial(matname);
359 _DBG_<<
"Missing material " << matname <<endl;
360 pthread_mutex_unlock(const_cast<pthread_mutex_t*>(&mutex));
361 return RESOURCE_UNAVAILABLE;
363 double A=mat->GetA();
364 double Z=mat->GetZ();
365 double density=mat->GetDensity();
366 rhoZ_overA=density*Z/A;
367 RadLen=mat->GetRadLen();
372 if (mat->IsMixture()) {
373 const TGeoMixture * mixt = dynamic_cast <
const TGeoMixture*> (mat);
374 for (
int i=0;i<mixt->GetNelements();i++){
375 double w_i=mixt->GetWmixt()[i];
376 double Z_i=mixt->GetZmixt()[i];
378 double I_i=7.0+12.0*Z_i;
379 if (Z_i>=13) I_i=9.76*Z_i+58.8*pow(Z_i,-0.19);
381 LnI+=w_i*Z_i*log(I_i)/Z;
387 if (Z>=13) I=9.76*Z+58.8*pow(Z,-0.19);
391 rhoZ_overA_logI=rhoZ_overA*LnI;
393 pthread_mutex_unlock(const_cast<pthread_mutex_t*>(&mutex));
407 double density, A, Z,LnI;
408 FindMatLL(pos, density, A, Z, RadLen,LnI);
409 rhoZ_overA = density*Z/A;
410 rhoZ_overA_logI = rhoZ_overA*LnI;
419 double &RadLen)
const{
420 density=RadLen=A=Z=0.;
422 pthread_mutex_lock(const_cast<pthread_mutex_t*>(&mutex));
424 if(!DRGeom)((
DRootGeom*)
this)->InitDRGeom();
426 TGeoNode *cnode = DRGeom->FindNode(pos.X(),pos.Y(),pos.Z());
428 _DBG_<<
"Missing cnode at position (" <<pos.X()<<
","<<pos.Y()<<
","
429 <<pos.Z()<<
")"<<endl;
430 pthread_mutex_unlock(const_cast<pthread_mutex_t*>(&mutex));
431 return RESOURCE_UNAVAILABLE;
433 TGeoVolume *cvol = cnode->GetVolume();
435 _DBG_<<
"Missing cvol" <<endl;
436 pthread_mutex_unlock(const_cast<pthread_mutex_t*>(&mutex));
437 return RESOURCE_UNAVAILABLE;
439 TGeoMaterial *cmat = cvol->GetMedium()->GetMaterial();
441 density=cmat->GetDensity();
442 RadLen=cmat->GetRadLen();
448 pthread_mutex_unlock(const_cast<pthread_mutex_t*>(&mutex));
457 double &RadLen,
double &LnI
459 density=RadLen=A=Z=LnI=0.;
461 pthread_mutex_lock(const_cast<pthread_mutex_t*>(&mutex));
463 if(!DRGeom)((
DRootGeom*)
this)->InitDRGeom();
465 TGeoNode *cnode = DRGeom->FindNode(pos.X(),pos.Y(),pos.Z());
467 _DBG_<<
"Missing cnode at position (" <<pos.X()<<
","<<pos.Y()<<
","
468 <<pos.Z()<<
")"<<endl;
469 pthread_mutex_unlock(const_cast<pthread_mutex_t*>(&mutex));
470 return RESOURCE_UNAVAILABLE;
472 TGeoVolume *cvol = cnode->GetVolume();
474 _DBG_<<
"Missing cvol" <<endl;
475 pthread_mutex_unlock(const_cast<pthread_mutex_t*>(&mutex));
476 return RESOURCE_UNAVAILABLE;
478 TGeoMaterial *cmat = cvol->GetMedium()->GetMaterial();
479 density=cmat->GetDensity();
480 RadLen=cmat->GetRadLen();
489 if (cmat->IsMixture()) {
490 const TGeoMixture * mixt = dynamic_cast <
const TGeoMixture*> (cmat);
492 for (
int i=0;i<mixt->GetNelements();i++){
493 double w_i=mixt->GetWmixt()[i];
494 double Z_i=mixt->GetZmixt()[i];
496 double I_i=7.0+12.0*Z_i;
497 if (Z_i>=13) I_i=9.76*Z_i+58.8*pow(Z_i,-0.19);
499 LnI+=w_i*Z_i*log(I_i)/Z;
505 if (Z>=13) I=9.76*Z+58.8*pow(Z,-0.19);
521 pthread_mutex_lock(&mutex);
522 if(!DRGeom)InitDRGeom();
523 TGeoNode *cnode = DRGeom->FindNode(
x[0],
x[1],
x[2]);
524 TGeoVolume *cvol = cnode->GetVolume();
525 TGeoMaterial *cmat = cvol->GetMedium()->GetMaterial();
527 Current_Node = cnode;
528 Current_Volume = cvol;
529 Current_Material = cmat;
530 Mat_Index = cmat->GetIndex();
532 Mat.A = cmat->GetA();
533 Mat.Z = cmat->GetZ();
534 Mat.Density = cmat->GetDensity();
535 Mat.RadLen = cmat->GetRadLen();
537 pthread_mutex_unlock(&mutex);
struct VolMat FindMat(double *x)
jerror_t FindMatLL(DVector3 pos, double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const
TGeoNode * FindNode(double *x)
DRootGeom(JApplication *japp, unsigned int run_number=1)
int ReadMap(string namepath, int32_t runnumber)
jerror_t FindMatTable(DVector3 pos, double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const
TGeoVolume * FindVolume(double *x)