80 #include "StMessMgr.h"
81 #include "St_DataSetIter.h"
82 #include "PhysicalConstants.h"
83 #include "tables/St_tofSlatGeom_Table.h"
84 #include "StThreeVectorD.hh"
85 #include "StPhysicalHelixD.hh"
87 #include "StTofUtil/StTofGeometry.h"
88 #include "StTofUtil/tofPathLength.hh"
129 LOG_INFO <<
"StTofGeometry: unable to initialize without maker" << endm;
132 LOG_INFO <<
"StTofGeometry: eventtime = ";
133 maker->GetDateTime().Print();
136 mDbDataSet= maker->GetDataBase(
"Geometry/tof");
139 St_tofSlatGeom* mSlatGeom =
static_cast<St_tofSlatGeom*
>(mDbDataSet->
Find(
"tofSlatGeom"));
141 tofSlatGeom_st *mslats=
static_cast<tofSlatGeom_st*
>(mSlatGeom->GetArray());
142 int numRows=mSlatGeom->GetNRows();
143 LOG_INFO <<
"StTofGeometry: numRows = " << numRows << endm;
144 for (
int i=0;i<numRows;i++) mTofSlatVec.push_back(mslats[i]);
147 mTofParam.i_eta_max = 10;
148 mTofParam.i_eta_min = 1;
149 mTofParam.i_phi_max = 5;
150 mTofParam.i_phi_min = 1;
151 mTofParam.n_counter_eta = 10;
152 mTofParam.n_counter_phi = 5;
153 mTofParam.n_tray_eta = 1;
154 mTofParam.n_tray_phi = 1;
155 mTofParam.counter_thickness = 2;
156 mTofParam.counter_width = 4;
157 mTofParam.r = 213.95;
158 mTofParam.tray_height = 4.7;
159 mTofParam.tray_width = 10.795;
160 mTofParam.tray_length = 120.81;
161 mTofParam.tray_phi_zero = 0;
168 int StTofGeometry::calcSlatId(
const int iphi,
const int ieta)
const {
169 return ((ieta - 1) * 4 + iphi);
176 tofSlatGeom_st thisSlat;
178 thisSlat = mTofSlatVec[slatId-1];
180 LOG_INFO <<
"StTofGeometry: Warning: slatId ("<< slatId <<
") seriously out of range" << endm;
188 LOG_INFO <<
"StTofGeometry: Initializing default DAQ and SlatId mappings" << endm;
191 mTofDaqMap[ 0]=37; mTofDaqMap[ 1]=38; mTofDaqMap[ 2]=39; mTofDaqMap[ 3]=40; mTofDaqMap[ 4]=41;
192 mTofDaqMap[ 5]=33; mTofDaqMap[ 6]=34; mTofDaqMap[ 7]=35; mTofDaqMap[ 8]=36;
193 mTofDaqMap[ 9]=29; mTofDaqMap[10]=30; mTofDaqMap[11]=31; mTofDaqMap[12]=32;
194 mTofDaqMap[13]=25; mTofDaqMap[14]=26; mTofDaqMap[15]=27; mTofDaqMap[16]=28;
195 mTofDaqMap[17]=21; mTofDaqMap[18]=22; mTofDaqMap[19]=23; mTofDaqMap[20]=24;
196 mTofDaqMap[21]=17; mTofDaqMap[22]=18; mTofDaqMap[23]=19; mTofDaqMap[24]=20;
197 mTofDaqMap[25]=13; mTofDaqMap[26]=14; mTofDaqMap[27]=15; mTofDaqMap[28]=16;
198 mTofDaqMap[29]= 9; mTofDaqMap[30]=10; mTofDaqMap[31]=11; mTofDaqMap[32]=12;
199 mTofDaqMap[33]= 5; mTofDaqMap[34]= 6; mTofDaqMap[35]= 7; mTofDaqMap[36]= 8;
200 mTofDaqMap[37]= 1; mTofDaqMap[38]= 2; mTofDaqMap[39]= 3; mTofDaqMap[40]= 4;
212 for (
int i=0;i<41;i++) mTofSlatMap[mTofDaqMap[i]]=i+1;
220 tofSlatGeom_st thisSlat =
tofSlat(slatId);
221 double cosAng = thisSlat.cosang;
222 double sinAng = ::sqrt(1.0 - cosAng*cosAng);
223 double tanAng = fabs(sinAng/cosAng);
224 double r = (fabs(thisSlat.z) + thisSlat.r/tanAng) * sinAng;
225 double x = r * fabs(cosAng) * cos(thisSlat.phi);
226 double y = r * fabs(cosAng) * sin(thisSlat.phi);
227 double z = r * sinAng * cosAng/fabs(cosAng)*(-1.0);
229 return slatNormPoint;
241 tofSlatGeom_st thisSlat =
tofSlat(slatId);
242 int iEta = thisSlat.ieta;
243 int iPhi, centerSlatId;
246 centerSlatId=calcSlatId(iPhi,iEta);
250 for (iPhi=2;iPhi<4;iPhi++){
251 centerSlatId=calcSlatId(iPhi,iEta);
256 return planeNormPoint;
263 os <<
"------StTofGeometry::printGeo()------" << endm;
264 os <<
"eta id max & min = " << mTofParam.i_eta_max <<
" "
265 << mTofParam.i_eta_min << endm;
266 os <<
"phi id max & min = " << mTofParam.i_phi_max <<
" "
267 << mTofParam.i_phi_min << endm;
268 os <<
"counters/trays/eta(phi) = " << mTofParam.n_counter_eta <<
" "
269 << mTofParam.n_counter_phi << endm;
270 os <<
"trays in eta & phi = " << mTofParam.n_tray_eta <<
" "
271 << mTofParam.n_tray_phi << endm;
272 os <<
"slat thickness & width = " << mTofParam.counter_thickness <<
" "
273 << mTofParam.counter_width << endm;
274 os <<
"mean counter radius = " << mTofParam.r << endm;
275 os <<
"tray height & width = " << mTofParam.tray_height <<
" "
276 << mTofParam.tray_width << endm;
277 os <<
"tray length & phi0 = " << mTofParam.tray_length <<
" "
278 << mTofParam.tray_phi_zero << endm;
279 LOG_INFO <<
"---------------------------------------" << endm;
286 os <<
"------StTofGeometry::printSlat()------" << endm;
287 os <<
"Slat: id, tray, eta, phi = " <<
" " << slatId <<
" "
289 <<
" " <<
tofSlat(slatId).iphi << endm;
290 os <<
"ieta, cosang = "<<
tofSlat(slatId).ieta<<
" "
291 <<
tofSlat(slatId).cosang<< endm;
292 os <<
"eta, eta_max, etamin = " <<
tofSlat(slatId).eta<<
" "
293 <<
tofSlat(slatId).eta_max<<
" "
294 <<
tofSlat(slatId).eta_min<< endm;
295 os <<
"r, z, z_max, z_min = " <<
tofSlat(slatId).r<<
" "
298 <<
tofSlat(slatId).z_min<< endm;
299 os <<
"iphi, phi, phi_max, phi_min = "<<
tofSlat(slatId).iphi<<
" "
301 <<
tofSlat(slatId).phi_max <<
" "
302 <<
tofSlat(slatId).phi_min<< endm;
303 LOG_INFO <<
"------------------------------------" << endm;
311 float phi = point.phi();
314 float etaMin, etaMax;
315 if(tofSlat.eta < 0) {
316 etaMin = tofSlat.eta_min;
317 etaMax = tofSlat.eta_max;
320 etaMin = tofSlat.eta_max;
321 etaMax = tofSlat.eta_min;
325 if(point.pseudoRapidity() >= etaMin &&
326 point.pseudoRapidity() <= etaMax) {
328 if (phi >= tofSlat.phi_min && phi <= tofSlat.phi_max)
341 int trayEta = volumeId/100000;
342 int trayPhi =
static_cast<int>(fmod(volumeId,100000.)/1000) ;
343 int counterPhi =
static_cast<int>(fmod(volumeId,1000.)/100) ;
344 int counterEta =
static_cast<int>(fmod(volumeId,100.)) ;
347 phiId = 14 - trayPhi ;
348 if (phiId<1) phiId = phiId + 60 ;
350 phiId = phiId * 5 - counterPhi + 1;
352 phiId = phiId * 4 - counterPhi + 1;
353 etaId = counterEta + 10 ;
354 LOG_INFO <<
"StTofGeometry: WARNING TOFp tray not in EAST barrel" << endm;
356 else if (trayEta==2) {
357 phiId = trayPhi - 42 ;
358 if (phiId<1) phiId = phiId + 60 ;
360 phiId = phiId * 5 + counterPhi - 5;
362 phiId = phiId * 4 + counterPhi - 4;
363 etaId = 11 - counterEta ;
366 LOG_INFO<<
" StTofGeometry::tofSlatCrossId unknown trayId "<<trayEta<<endm ;
368 int slatId = calcSlatId(counterPhi,etaId);
379 for (
unsigned int i=0;i<mTofSlatVec.size();i++){
380 if(point.z() >= mTofSlatVec[i].z_min &&
381 point.z() <= mTofSlatVec[i].z_max &&
382 point.phi() >= mTofSlatVec[i].phi_min &&
383 point.phi() <= mTofSlatVec[i].phi_max) {
384 etaId = mTofSlatVec[i].ieta;
385 phiId = mTofSlatVec[i].iphi;
388 if (etaId!=-1 && phiId!=-1)
break;
393 if(etaId > 0 && phiId > 0) slatId = calcSlatId(phiId,etaId);
403 idVector slatIdVec) {
404 idVector idErasedVec = slatIdVec;
405 idVectorIter slatIdIter, idErasedIter;
410 tofSlatHitVector slatHitVec;
414 while (slatIdVec.size() != 0) {
417 slatIdIter = slatIdVec.begin();
419 int trayId = this->
tofSlat(*slatIdIter).trayId;
420 float cosang = this->tofSlat(*slatIdIter).cosang;
421 int iEta = this->tofSlat(*slatIdIter).ieta;
425 unsigned int idMiddleLayer = (mMaxSlatLayers-1)/2;
426 float layerSeperation = mTofParam.counter_thickness/(mMaxSlatLayers-1);
429 StThreeVectorD slatNormalVec = slatNormLayer[idMiddleLayer]/slatNormLayer[idMiddleLayer].mag();
430 for (
unsigned int ll=1;ll<(idMiddleLayer+1);ll++){
431 slatNormLayer[idMiddleLayer+ll] = slatNormLayer[idMiddleLayer]
432 *((slatNormLayer[idMiddleLayer].mag() - ll*layerSeperation )
433 /slatNormLayer[idMiddleLayer].mag());
434 slatNormLayer[idMiddleLayer-ll] = slatNormLayer[idMiddleLayer]
435 *((slatNormLayer[idMiddleLayer].mag() + ll*layerSeperation )
436 /slatNormLayer[idMiddleLayer].mag());
439 for (
unsigned int ll=0;ll<mMaxSlatLayers;ll++){
440 pathLength = helix.
pathLength(slatNormLayer[ll], slatNormalVec);
441 pathLength = (pathLength>0) ? pathLength : 1.0e+33;
442 hitAtLayer[ll] = helix.at(pathLength);
446 idErasedIter = idErasedVec.begin();
447 while (idErasedIter != idErasedVec.end()) {
450 if(this->
tofSlat(*idErasedIter).cosang == cosang &&
451 this->tofSlat(*idErasedIter).ieta == iEta &&
452 this->tofSlat(*idErasedIter).trayId == trayId) {
454 bool layer[mMaxSlatLayers];
455 int numberOfHitLayers(0);
456 for (
unsigned int ll=0;ll<mMaxSlatLayers;ll++){
458 if (layer[ll]) numberOfHitLayers++;
462 if (numberOfHitLayers>0) {
463 slatHit.hitProfile = 0;
464 if (Debug()) LOG_INFO <<
"L(0-"<<mMaxSlatLayers-1<<
"): ";
465 for (
unsigned int ll=0;ll<mMaxSlatLayers;ll++){
466 slatHit.hitProfile = 2*slatHit.hitProfile + layer[ll]*1;
467 if (Debug()) LOG_INFO <<layer[ll];
469 if (Debug()) LOG_INFO << endm;
471 unsigned int innerLayer(0), outerLayer(mMaxSlatLayers-1);
473 while (!layer[innerLayer] && (innerLayer<=mMaxSlatLayers)) innerLayer++;
474 while (!layer[outerLayer] && (outerLayer>= 0 )) outerLayer--;
475 if (Debug()) LOG_INFO <<
"Li="<<innerLayer<<
" Lo="<<outerLayer<<
" nLx="<<numberOfHitLayers<<endm;
481 if (innerLayer!=outerLayer){
482 s = tofPathLength(&hitAtLayer[innerLayer],&hitAtLayer[outerLayer],helix.curvature());
483 distance = hitAtLayer[outerLayer] - hitAtLayer[innerLayer];
486 float theta_xy = atan2(distance.x(),distance.y());
487 float theta_zr = atan2(distance.z(),distance.perp());
490 theta_zr -= acos(this->
tofSlat(*idErasedIter).cosang);
491 theta_xy += pi/2 + this->
tofSlat(*idErasedIter).phi;
494 slatHit.theta_xy = theta_xy;
495 slatHit.theta_zr = theta_zr;
497 slatHit.slatIndex = *idErasedIter;
501 for (
unsigned int i=0;i<mMaxSlatLayers;i++)
503 averageHitPos += hitAtLayer[i]/numberOfHitLayers;
504 slatHit.layerHitPositions.push_back(hitAtLayer[i]);
506 slatHit.hitPosition = averageHitPos;
508 slatHitVec.push_back(slatHit);
509 slatHit.layerHitPositions.clear();
512 idErasedVec.erase(idErasedIter);
517 slatIdVec = idErasedVec;
527 float zmin = this->
tofSlat(slatId).z_min;
528 float zmax = this->
tofSlat(slatId).z_max;
529 float cosang = this->
tofSlat(slatId).cosang;
532 float length = (zmax-hitPoint->z())/cosang ;
533 float max_distance = (zmax-zmin)/cosang ;
535 if (length>max_distance || length<0){
536 LOG_INFO <<
"HitPositionCorrection: length="<<length<<
" max="<<max_distance
537 <<
" zmin="<<zmin<<
" zmax="<<zmax<<
" cosang="<<cosang<<endm;
545 LOG_INFO <<
"slatHitPosition: hit point out of range;"
546 <<
" phi=" << hitPoint->phi() <<
" z=" << hitPoint->z() << endm;
553 tofSlatGeom_st middleSlat = this->
tofSlat(slatId);
554 int slatPhi = middleSlat.iphi;
555 int slatEta = middleSlat.ieta;
559 int thisSlatId, ieta, iphi;
563 for (ieta= slatEta-1;ieta<=slatEta+1;ieta++)
564 for (iphi= slatPhi-1;iphi<=slatPhi+1;iphi++){
565 if ((ieta==slatEta) && (iphi==slatPhi))
continue;
566 if ((ieta <1) || (iphi<1) || (iphi>4))
continue;
567 thisSlatId = calcSlatId(iphi,ieta);
568 neighbours.push_back(thisSlatId);
572 else if (slatEta==9){
573 for (ieta= slatEta-1;ieta<=slatEta;ieta++)
574 for (iphi= slatPhi-1;iphi<=slatPhi+1;iphi++){
575 if ((ieta==slatEta) && (iphi==slatPhi))
continue;
576 if ((iphi<1) || (iphi>4))
continue;
577 thisSlatId = calcSlatId(iphi,ieta);
578 neighbours.push_back(thisSlatId);
580 ieta = slatEta+1; iphi = slatPhi;
581 thisSlatId = calcSlatId(iphi,ieta);
582 neighbours.push_back(thisSlatId);
584 thisSlatId = calcSlatId(iphi,ieta);
585 neighbours.push_back(thisSlatId);
588 else if (slatEta==10){
589 for (iphi= slatPhi-1;iphi<=slatPhi+1;iphi++){
590 if (iphi==slatPhi)
continue;
591 if ((iphi<1) || (iphi>5))
continue;
593 thisSlatId = calcSlatId(iphi,ieta);
594 neighbours.push_back(thisSlatId);
597 for (iphi=slatPhi-1;iphi<=slatPhi;iphi++){
598 if ((iphi<1) || (iphi>4))
continue;
599 thisSlatId = calcSlatId(iphi,ieta);
600 neighbours.push_back(thisSlatId);
611 tofSlatGeom_st middleSlat = this->
tofSlat(slatId);
612 int slatPhi = middleSlat.iphi;
613 int slatEta = middleSlat.ieta;
617 int thisSlatId, ieta, iphi;
619 for (ieta=slatEta-2;ieta<=slatEta+2;ieta++){
620 if ((ieta<1) || (ieta>10))
continue;
621 int iphiMax = (ieta==10)?5:4;
622 for (iphi = 1; iphi<=iphiMax; iphi++){
623 if ((ieta==slatEta) && (iphi==slatPhi))
continue;
624 thisSlatId = calcSlatId(iphi,ieta);
625 neighbours.push_back(thisSlatId);
637 float phimin = this->
tofSlat(slatId).phi_min;
638 float phimax = this->
tofSlat(slatId).phi_max;
640 float philoc = hitPoint->phi() - phimin;
641 float max_deltaphi = (phimax-phimin);
643 if (philoc>max_deltaphi || philoc<0){
644 LOG_INFO <<
"slatHitPhiPosition: phi="<<philoc<<
" max="<<max_deltaphi
645 <<
" phimin="<<phimin<<
" phimax="<<phimax<<endm;
653 LOG_INFO <<
"slatHitPhiPosition: hit point out of range;"
654 <<
" phi=" << hitPoint->phi() <<
" z=" << hitPoint->z() << endm;
660 Bool_t StTofGeometry::projTrayVector(
const StHelixD &helix, idVector &trayVec)
const {
666 if(s1<0.) s1 = helix.
pathLength(R_tof).second;
668 double phi = point.phi()*180/3.14159;
671 int itray_east = (255+(int)phi)%360/6+61;
672 trayVec.push_back(itray_east);
674 int itray_east1 = (255+(int)(phi+res))%360/6+61;
675 int itray_east2 = (255+(int)(phi-res))%360/6+61;
676 if(itray_east1!=itray_east) {
677 trayVec.push_back(itray_east1);
679 if(itray_east2!=itray_east&&itray_east2!=itray_east1) {
680 trayVec.push_back(itray_east2);
684 int itray_west = (435-(int)phi)%360/6+1;
685 trayVec.push_back(itray_west);
687 int itray_west1 = (435-(int)(phi+res))%360/6+1;
688 int itray_west2 = (435-(int)(phi-res))%360/6+1;
689 if(itray_west1!=itray_west) {
690 trayVec.push_back(itray_west1);
692 if(itray_west2!=itray_west&&itray_west2!=itray_west1) {
693 trayVec.push_back(itray_west2);
702 if(trayVec.size()>0)
return kTRUE;
tofSlatGeom_st tofSlat(const Int_t slatId) const
return slat geometry structure for slatId
idVector slatNeighbours(const int)
returns idVector of 3x3 (max) neighbouring slatIds
void printGeo(ostream &os=cout) const
print global geometry parameters
idVector slatNeighboursWide(const int)
returns idVector of 5x5 (max) neighbouring slatIds
void initDaqMap()
set-up the default Daq-to-SlatId and Slat-to-DaqId mappings
void initGeomFromXdf(const Char_t *="/afs/rhic.bnl.gov/star/users/geurts/public/dbase/ctg_pars.xdf")
initialize TOF Slat parameters from XDF file
tofSlatHitVector tofHelixToArray(const StPhysicalHelixD &helix, idVector slatIdVec)
finds slats in an array of trays which are crossed by a track-helix.
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
float slatPhiPosition(StThreeVectorD *)
returns the local Phi angle of a track inside a slat
float slatHitPosition(StThreeVectorD *)
returns 1-D hit position on the TOFp slat (Zhit)
StThreeVectorD tofSlatNormPoint(const Int_t slatId) const
calculate the normal vector <r> to a slat
~StTofGeometry()
default empty destructor
StThreeVectorD tofPlaneNormPoint(const Int_t slatId) const
calculate the normal vector to a slats-plane
int tofSlatCross(const StThreeVectorD &point, const tofSlatGeom_st tofSlat) const
check if a point is in a slat
void printSlat(const Int_t slatId, ostream &os=cout) const
print slat-specific geometry parameters
int tofSlatCrossId(const StThreeVectorD &point) const
return the index of a slat if the point is in the slat
void initGeomFromDbase(StMaker *)
initialize TOF Slat parameters from STAR dBase
StTofGeometry()
defaulty constructor
virtual TDataSet * Find(const char *path) const
void init()
initialize geometry class from XDF file and set-up DAQ/Slat mappings