00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00057
00076
00077 #include "StMessMgr.h"
00078 #include "St_DataSetIter.h"
00079 #include "PhysicalConstants.h"
00080 #include "ctf/St_ctg_Module.h"
00081 #include "tables/St_tofSlatGeom_Table.h"
00082 #include "StThreeVectorD.hh"
00083 #include "StPhysicalHelixD.hh"
00084 #include "StMaker.h"
00085 #include "StTofUtil/StTofGeometry.h"
00086 #include "StTofUtil/tofPathLength.hh"
00087
00088 #include <utility>
00089 using std::pair;
00090
00091
00093 StTofGeometry::StTofGeometry(){
00094 mDebug = false;
00095 }
00096
00098 StTofGeometry::~StTofGeometry(){ };
00099
00100
00101
00103 void StTofGeometry::init(){
00104 initGeomFromXdf();
00105 initDaqMap();
00106 }
00107
00108
00109
00111 void StTofGeometry::init(StMaker* maker){
00112 initGeomFromDbase(maker);
00113 initDaqMap();
00114 }
00115
00116
00117
00119 void StTofGeometry::initGeomFromXdf(const Char_t* InputXdfFile){
00120 #if 0
00121 LOG_INFO << "StTofGeometry: loading dBase from " << InputXdfFile << endm;
00122 St_XDFFile xdf(InputXdfFile);
00123 St_DataSet *ctfg = xdf.NextEventGet();
00124 St_DataSetIter gime(ctfg);
00125 if (!ctfg){
00126 LOG_INFO << " ERROR: unable read from file" << endm;
00127 assert(0);
00128 return;
00129 }
00130
00131 St_ctg_geo* stafTofParam(0);
00132 St_ctg_slat_phi* stafSlatPhi(0);
00133 St_ctg_slat_eta* stafSlatEta(0);
00134 St_ctg_slat* stafSlatParam(0);
00135
00136 stafTofParam = static_cast<St_ctg_geo*>(gime("tof"));
00137 stafSlatPhi = static_cast<St_ctg_slat_phi*>(gime("tof_slat_phi"));
00138 stafSlatEta = static_cast<St_ctg_slat_eta*>(gime("tof_slat_eta"));
00139 stafSlatParam = static_cast<St_ctg_slat*>(gime("tof_slat"));
00140
00141 ctg_geo_st *geo = stafTofParam->GetTable();
00142 if (geo){
00143
00144 mTofParam.i_eta_max = geo->i_eta_max;
00145 mTofParam.i_eta_min = geo->i_eta_min;
00146 mTofParam.i_phi_max = geo->i_phi_max;
00147 mTofParam.i_phi_min = geo->i_phi_min;
00148 mTofParam.n_counter_eta = geo->n_counter_eta;
00149 mTofParam.n_counter_phi = geo->n_counter_phi;
00150 mTofParam.n_tray_eta = geo->n_tray_eta;
00151 mTofParam.n_tray_phi = geo->n_tray_phi;
00152 mTofParam.counter_thickness = geo->counter_thickness;
00153 mTofParam.counter_width = geo->counter_width;
00154 mTofParam.r = geo->r;
00155 mTofParam.tray_height = geo->tray_height;
00156 mTofParam.tray_width = geo->tray_width;
00157 mTofParam.tray_length = geo->tray_length;
00158 mTofParam.tray_phi_zero = geo->tray_phi_zero;
00159 }
00160 else {
00161 LOG_INFO << " ERROR: unable to read TOF param table" << endm;
00162 return;
00163 }
00164
00165
00166 ctg_slat_eta_st *eta = stafSlatEta->GetTable();
00167 if (eta){
00168 for(int iRow=0; iRow<stafSlatEta->GetNRows(); iRow++,eta++){
00169 StructTofSlatEta tofSlatEta;
00170 tofSlatEta.ieta = eta->ieta;
00171 tofSlatEta.cosang = eta->cosang;
00172 tofSlatEta.eta = eta->eta;
00173 tofSlatEta.eta_max = eta->eta_max;
00174 tofSlatEta.eta_min = eta->eta_min;
00175 tofSlatEta.r = eta->r;
00176 tofSlatEta.z = eta->z;
00177 tofSlatEta.z_max = eta->z_max;
00178 tofSlatEta.z_min = eta->z_min;
00179 mTofSlatEtaVec.push_back(tofSlatEta);
00180 }
00181 }
00182 else {
00183 LOG_INFO << " ERROR: unable to read TOF eta table" << endm;
00184 return;
00185 }
00186
00187 ctg_slat_phi_st *phi = stafSlatPhi->GetTable();
00188 if (phi){
00189 for(int iRow=0; iRow<stafSlatPhi->GetNRows(); iRow++,phi++){
00190 StructTofSlatPhi tofSlatPhi;
00191 tofSlatPhi.iphi = phi->iphi;
00192
00193 tofSlatPhi.phi = ((phi->phi >180)?phi->phi -360:phi>phi )*degree;
00194 tofSlatPhi.phi_max = ((phi->phi_max>180)?phi->phi_max-360:phi->phi_max)*degree;
00195 tofSlatPhi.phi_min = ((phi->phi_min>180)?phi->phi_min-360:phi->phi_min)*degree;
00196 mTofSlatPhiVec.push_back(tofSlatPhi);
00197 }
00198 }
00199 else {
00200 LOG_INFO << " ERROR: unable to read TOF phi table" << endm;
00201 return;
00202 }
00203
00204
00205 for (int i=0; i<stafSlatEta->GetNRows();i++){
00206 int iPhiMin, iPhiMax;
00207 int iEta = mTofSlatEtaVec[i].ieta;
00208 if (iEta==10) {iPhiMin= 1 ; iPhiMax= 5;}
00209 else if (iEta== 9) {iPhiMin= 6 ; iPhiMax= 9;}
00210 else if (iEta < 9) {iPhiMin=10 ; iPhiMax=13;}
00211 else {
00212 LOG_INFO << "StTofGeometry: slat eta out of range " << iEta << endm;
00213 iPhiMin=0; iPhiMax=-1;}
00214 for (int j=iPhiMin-1;j<iPhiMax;j++){
00215 tofSlatGeom_st* tofSlat = new tofSlatGeom_st;
00216 tofSlat->ieta = mTofSlatEtaVec[i].ieta;
00217 tofSlat->z = mTofSlatEtaVec[i].z;
00218 tofSlat->z_min = mTofSlatEtaVec[i].z_min;
00219 tofSlat->z_max = mTofSlatEtaVec[i].z_max;
00220 tofSlat->cosang = mTofSlatEtaVec[i].cosang;
00221 tofSlat->r = mTofSlatEtaVec[i].r;
00222 tofSlat->eta = mTofSlatEtaVec[i].eta;
00223 tofSlat->eta_min = mTofSlatEtaVec[i].eta_min;
00224 tofSlat->eta_max = mTofSlatEtaVec[i].eta_max;
00225 tofSlat->iphi = mTofSlatPhiVec[j].iphi;
00226 tofSlat->phi = mTofSlatPhiVec[j].phi;
00227 tofSlat->phi_min = mTofSlatPhiVec[j].phi_min;
00228 tofSlat->phi_max = mTofSlatPhiVec[j].phi_max;
00229 tofSlat->trayId = 32;
00230 mTofSlatVec.push_back(*tofSlat);
00231 delete tofSlat;
00232 }
00233 }
00234 #endif //0
00235 }
00236
00237
00238
00240 void StTofGeometry::initGeomFromDbase(StMaker *maker){
00241 if (!maker){
00242 LOG_INFO << "StTofGeometry: unable to initialize without maker" << endm;
00243 return;
00244 }
00245 LOG_INFO << "StTofGeometry: eventtime = ";
00246 maker->GetDateTime().Print();
00247
00248 TDataSet* mDbDataSet;
00249 mDbDataSet= maker->GetDataBase("Geometry/tof");
00250
00251 assert(mDbDataSet);
00252 St_tofSlatGeom* mSlatGeom = static_cast<St_tofSlatGeom*>(mDbDataSet->Find("tofSlatGeom"));
00253 assert(mSlatGeom);
00254 tofSlatGeom_st *mslats= static_cast<tofSlatGeom_st*>(mSlatGeom->GetArray());
00255 int numRows=mSlatGeom->GetNRows();
00256 LOG_INFO << "StTofGeometry: numRows = " << numRows << endm;
00257 for (int i=0;i<numRows;i++) mTofSlatVec.push_back(mslats[i]);
00258
00259
00260 mTofParam.i_eta_max = 10;
00261 mTofParam.i_eta_min = 1;
00262 mTofParam.i_phi_max = 5;
00263 mTofParam.i_phi_min = 1;
00264 mTofParam.n_counter_eta = 10;
00265 mTofParam.n_counter_phi = 5;
00266 mTofParam.n_tray_eta = 1;
00267 mTofParam.n_tray_phi = 1;
00268 mTofParam.counter_thickness = 2;
00269 mTofParam.counter_width = 4;
00270 mTofParam.r = 213.95;
00271 mTofParam.tray_height = 4.7;
00272 mTofParam.tray_width = 10.795;
00273 mTofParam.tray_length = 120.81;
00274 mTofParam.tray_phi_zero = 0;
00275
00276 }
00277
00278
00279
00281 int StTofGeometry::calcSlatId(const int iphi, const int ieta) const {
00282 return ((ieta - 1) * 4 + iphi);
00283 }
00284
00285
00286
00288 tofSlatGeom_st StTofGeometry::tofSlat(const Int_t slatId) const {
00289 tofSlatGeom_st thisSlat;
00290 if(slatId > 0)
00291 thisSlat = mTofSlatVec[slatId-1];
00292 else
00293 LOG_INFO << "StTofGeometry: Warning: slatId ("<< slatId <<") seriously out of range" << endm;
00294 return thisSlat;
00295 }
00296
00297
00298
00300 void StTofGeometry::initDaqMap(){
00301 LOG_INFO << "StTofGeometry: Initializing default DAQ and SlatId mappings" << endm;
00302
00303
00304 mTofDaqMap[ 0]=37; mTofDaqMap[ 1]=38; mTofDaqMap[ 2]=39; mTofDaqMap[ 3]=40; mTofDaqMap[ 4]=41;
00305 mTofDaqMap[ 5]=33; mTofDaqMap[ 6]=34; mTofDaqMap[ 7]=35; mTofDaqMap[ 8]=36;
00306 mTofDaqMap[ 9]=29; mTofDaqMap[10]=30; mTofDaqMap[11]=31; mTofDaqMap[12]=32;
00307 mTofDaqMap[13]=25; mTofDaqMap[14]=26; mTofDaqMap[15]=27; mTofDaqMap[16]=28;
00308 mTofDaqMap[17]=21; mTofDaqMap[18]=22; mTofDaqMap[19]=23; mTofDaqMap[20]=24;
00309 mTofDaqMap[21]=17; mTofDaqMap[22]=18; mTofDaqMap[23]=19; mTofDaqMap[24]=20;
00310 mTofDaqMap[25]=13; mTofDaqMap[26]=14; mTofDaqMap[27]=15; mTofDaqMap[28]=16;
00311 mTofDaqMap[29]= 9; mTofDaqMap[30]=10; mTofDaqMap[31]=11; mTofDaqMap[32]=12;
00312 mTofDaqMap[33]= 5; mTofDaqMap[34]= 6; mTofDaqMap[35]= 7; mTofDaqMap[36]= 8;
00313 mTofDaqMap[37]= 1; mTofDaqMap[38]= 2; mTofDaqMap[39]= 3; mTofDaqMap[40]= 4;
00314
00315
00316 mTofDaqMap[41]=99;
00317 mTofDaqMap[42]=51;
00318 mTofDaqMap[43]=52;
00319 mTofDaqMap[44]=53;
00320 mTofDaqMap[45]=54;
00321 mTofDaqMap[46]=55;
00322 mTofDaqMap[47]=56;
00323
00324
00325 for (int i=0;i<41;i++) mTofSlatMap[mTofDaqMap[i]]=i+1;
00326 mTofSlatMap[0]=99;
00327 }
00328
00329
00330
00332 StThreeVectorD StTofGeometry::tofSlatNormPoint(const Int_t slatId) const {
00333 tofSlatGeom_st thisSlat = tofSlat(slatId);
00334 double cosAng = thisSlat.cosang;
00335 double sinAng = ::sqrt(1.0 - cosAng*cosAng);
00336 double tanAng = fabs(sinAng/cosAng);
00337 double r = (fabs(thisSlat.z) + thisSlat.r/tanAng) * sinAng;
00338 double x = r * fabs(cosAng) * cos(thisSlat.phi);
00339 double y = r * fabs(cosAng) * sin(thisSlat.phi);
00340 double z = r * sinAng * cosAng/fabs(cosAng)*(-1.0);
00341 StThreeVectorD slatNormPoint = StThreeVectorD(x,y,z);
00342 return slatNormPoint;
00343 }
00344
00345
00346
00348
00352 StThreeVectorD StTofGeometry::tofPlaneNormPoint(Int_t slatId) const {
00353 StThreeVectorD planeNormPoint(0.0, 0.0, 0.0);
00354 tofSlatGeom_st thisSlat = tofSlat(slatId);
00355 int iEta = thisSlat.ieta;
00356 int iPhi, centerSlatId;
00357 if (iEta==10){
00358 iPhi=3;
00359 centerSlatId=calcSlatId(iPhi,iEta);
00360 planeNormPoint = tofSlatNormPoint(centerSlatId);
00361 }
00362 else {
00363 for (iPhi=2;iPhi<4;iPhi++){
00364 centerSlatId=calcSlatId(iPhi,iEta);
00365 planeNormPoint += tofSlatNormPoint(centerSlatId)/2;
00366 }
00367
00368 }
00369 return planeNormPoint;
00370 }
00371
00372
00373
00375 void StTofGeometry::printGeo(ostream& os) const {
00376 os << "------StTofGeometry::printGeo()------" << endm;
00377 os << "eta id max & min = " << mTofParam.i_eta_max << " "
00378 << mTofParam.i_eta_min << endm;
00379 os << "phi id max & min = " << mTofParam.i_phi_max << " "
00380 << mTofParam.i_phi_min << endm;
00381 os << "counters/trays/eta(phi) = " << mTofParam.n_counter_eta << " "
00382 << mTofParam.n_counter_phi << endm;
00383 os << "trays in eta & phi = " << mTofParam.n_tray_eta << " "
00384 << mTofParam.n_tray_phi << endm;
00385 os << "slat thickness & width = " << mTofParam.counter_thickness << " "
00386 << mTofParam.counter_width << endm;
00387 os << "mean counter radius = " << mTofParam.r << endm;
00388 os << "tray height & width = " << mTofParam.tray_height << " "
00389 << mTofParam.tray_width << endm;
00390 os << "tray length & phi0 = " << mTofParam.tray_length << " "
00391 << mTofParam.tray_phi_zero << endm;
00392 LOG_INFO << "---------------------------------------" << endm;
00393 }
00394
00395
00396
00398 void StTofGeometry::printSlat(const Int_t slatId, ostream& os) const {
00399 os << "------StTofGeometry::printSlat()------" << endm;
00400 os << "Slat: id, tray, eta, phi = " << " " << slatId << " "
00401 << tofSlat(slatId).trayId << " " << tofSlat(slatId).ieta
00402 << " " << tofSlat(slatId).iphi << endm;
00403 os << "ieta, cosang = "<< tofSlat(slatId).ieta<< " "
00404 << tofSlat(slatId).cosang<< endm;
00405 os << "eta, eta_max, etamin = " << tofSlat(slatId).eta<< " "
00406 << tofSlat(slatId).eta_max<< " "
00407 << tofSlat(slatId).eta_min<< endm;
00408 os << "r, z, z_max, z_min = " << tofSlat(slatId).r<< " "
00409 << tofSlat(slatId).z<< " "
00410 << tofSlat(slatId).z_max<< " "
00411 << tofSlat(slatId).z_min<< endm;
00412 os << "iphi, phi, phi_max, phi_min = "<< tofSlat(slatId).iphi<< " "
00413 << tofSlat(slatId).phi<< " "
00414 << tofSlat(slatId).phi_max << " "
00415 << tofSlat(slatId).phi_min<< endm;
00416 LOG_INFO << "------------------------------------" << endm;
00417 }
00418
00419
00420
00422 int StTofGeometry::tofSlatCross(const StThreeVectorD& point, const tofSlatGeom_st tofSlat) const {
00423 int slatCross=0;
00424 float phi = point.phi();
00425
00426
00427 float etaMin, etaMax;
00428 if(tofSlat.eta < 0) {
00429 etaMin = tofSlat.eta_min;
00430 etaMax = tofSlat.eta_max;
00431 }
00432 else {
00433 etaMin = tofSlat.eta_max;
00434 etaMax = tofSlat.eta_min;
00435 }
00436
00437
00438 if(point.pseudoRapidity() >= etaMin &&
00439 point.pseudoRapidity() <= etaMax) {
00440
00441 if (phi >= tofSlat.phi_min && phi <= tofSlat.phi_max)
00442 slatCross = 1;
00443 }
00444 return slatCross;
00445 }
00446
00447
00448
00450 int StTofGeometry::tofSlatCrossId(const int volumeId) const {
00451 int phiId = -1;
00452 int etaId = -1;
00453
00454 int trayEta = volumeId/100000;
00455 int trayPhi = static_cast<int>(fmod(volumeId,100000.)/1000) ;
00456 int counterPhi = static_cast<int>(fmod(volumeId,1000.)/100) ;
00457 int counterEta = static_cast<int>(fmod(volumeId,100.)) ;
00458
00459 if (trayEta==1) {
00460 phiId = 14 - trayPhi ;
00461 if (phiId<1) phiId = phiId + 60 ;
00462 if (counterEta==1)
00463 phiId = phiId * 5 - counterPhi + 1;
00464 else
00465 phiId = phiId * 4 - counterPhi + 1;
00466 etaId = counterEta + 10 ;
00467 LOG_INFO << "StTofGeometry: WARNING TOFp tray not in EAST barrel" << endm;
00468 }
00469 else if (trayEta==2) {
00470 phiId = trayPhi - 42 ;
00471 if (phiId<1) phiId = phiId + 60 ;
00472 if (counterEta==1)
00473 phiId = phiId * 5 + counterPhi - 5;
00474 else
00475 phiId = phiId * 4 + counterPhi - 4;
00476 etaId = 11 - counterEta ;
00477 }
00478 else
00479 LOG_INFO<<" StTofGeometry::tofSlatCrossId unknown trayId "<<trayEta<<endm ;
00480
00481 int slatId = calcSlatId(counterPhi,etaId);
00482 return slatId;
00483 }
00484
00485
00486
00488 int StTofGeometry::tofSlatCrossId(const StThreeVectorD& point) const {
00489 int etaId = -1;
00490 int phiId = -1;
00491
00492 for (unsigned int i=0;i<mTofSlatVec.size();i++){
00493 if(point.z() >= mTofSlatVec[i].z_min &&
00494 point.z() <= mTofSlatVec[i].z_max &&
00495 point.phi() >= mTofSlatVec[i].phi_min &&
00496 point.phi() <= mTofSlatVec[i].phi_max) {
00497 etaId = mTofSlatVec[i].ieta;
00498 phiId = mTofSlatVec[i].iphi;
00499 break;
00500 }
00501 if (etaId!=-1 && phiId!=-1) break;
00502 }
00503
00504
00505 int slatId;
00506 if(etaId > 0 && phiId > 0) slatId = calcSlatId(phiId,etaId);
00507 else slatId = -1;
00508
00509 return slatId;
00510 }
00511
00512
00513
00515 tofSlatHitVector StTofGeometry::tofHelixToArray(const StPhysicalHelixD& helix,
00516 idVector slatIdVec) {
00517 idVector idErasedVec = slatIdVec;
00518 idVectorIter slatIdIter, idErasedIter;
00519
00520 double pathLength;
00521
00522 StructSlatHit slatHit;
00523 tofSlatHitVector slatHitVec;
00524 slatHitVec.clear();
00525
00526
00527 while (slatIdVec.size() != 0) {
00528
00529
00530 slatIdIter = slatIdVec.begin();
00531
00532 int trayId = this->tofSlat(*slatIdIter).trayId;
00533 float cosang = this->tofSlat(*slatIdIter).cosang;
00534 int iEta = this->tofSlat(*slatIdIter).ieta;
00535
00536
00537
00538 unsigned int idMiddleLayer = (mMaxSlatLayers-1)/2;
00539 float layerSeperation = mTofParam.counter_thickness/(mMaxSlatLayers-1);
00540 StThreeVectorD slatNormLayer[mMaxSlatLayers];
00541 slatNormLayer[idMiddleLayer] = this->tofPlaneNormPoint(*slatIdIter);
00542 StThreeVectorD slatNormalVec = slatNormLayer[idMiddleLayer]/slatNormLayer[idMiddleLayer].mag();
00543 for (unsigned int ll=1;ll<(idMiddleLayer+1);ll++){
00544 slatNormLayer[idMiddleLayer+ll] = slatNormLayer[idMiddleLayer]
00545 *((slatNormLayer[idMiddleLayer].mag() - ll*layerSeperation )
00546 /slatNormLayer[idMiddleLayer].mag());
00547 slatNormLayer[idMiddleLayer-ll] = slatNormLayer[idMiddleLayer]
00548 *((slatNormLayer[idMiddleLayer].mag() + ll*layerSeperation )
00549 /slatNormLayer[idMiddleLayer].mag());
00550 }
00551 StThreeVectorD hitAtLayer[mMaxSlatLayers];
00552 for (unsigned int ll=0;ll<mMaxSlatLayers;ll++){
00553 pathLength = helix.pathLength(slatNormLayer[ll], slatNormalVec);
00554 pathLength = (pathLength>0) ? pathLength : 1.0e+33;
00555 hitAtLayer[ll] = helix.at(pathLength);
00556 }
00557
00558
00559 idErasedIter = idErasedVec.begin();
00560 while (idErasedIter != idErasedVec.end()) {
00561
00562
00563 if(this->tofSlat(*idErasedIter).cosang == cosang &&
00564 this->tofSlat(*idErasedIter).ieta == iEta &&
00565 this->tofSlat(*idErasedIter).trayId == trayId) {
00566
00567 bool layer[mMaxSlatLayers];
00568 int numberOfHitLayers(0);
00569 for (unsigned int ll=0;ll<mMaxSlatLayers;ll++){
00570 layer[ll] = tofSlatCross(hitAtLayer[ll], this->tofSlat(*idErasedIter));
00571 if (layer[ll]) numberOfHitLayers++;
00572 }
00573
00574
00575 if (numberOfHitLayers>0) {
00576 slatHit.hitProfile = 0;
00577 if (Debug()) LOG_INFO << "L(0-"<<mMaxSlatLayers-1<<"): ";
00578 for (unsigned int ll=0;ll<mMaxSlatLayers;ll++){
00579 slatHit.hitProfile = 2*slatHit.hitProfile + layer[ll]*1;
00580 if (Debug()) LOG_INFO <<layer[ll];
00581 }
00582 if (Debug()) LOG_INFO << endm;
00583
00584 unsigned int innerLayer(0), outerLayer(mMaxSlatLayers-1);
00585
00586 while (!layer[innerLayer] && (innerLayer<=mMaxSlatLayers)) innerLayer++;
00587 while (!layer[outerLayer] && (outerLayer>= 0 )) outerLayer--;
00588 if (Debug()) LOG_INFO << "Li="<<innerLayer<< " Lo="<<outerLayer<<" nLx="<<numberOfHitLayers<<endm;
00589
00590
00591 float s(0);
00592 StThreeVectorD distance(0,0,0);
00593
00594 if (innerLayer!=outerLayer){
00595 s = tofPathLength(&hitAtLayer[innerLayer],&hitAtLayer[outerLayer],helix.curvature());
00596 distance = hitAtLayer[outerLayer] - hitAtLayer[innerLayer];
00597 }
00598
00599 float theta_xy = atan2(distance.x(),distance.y());
00600 float theta_zr = atan2(distance.z(),distance.perp());
00601
00602
00603 theta_zr -= acos(this->tofSlat(*idErasedIter).cosang);
00604 theta_xy += pi/2 + this->tofSlat(*idErasedIter).phi;
00605
00606 slatHit.s = s;
00607 slatHit.theta_xy = theta_xy;
00608 slatHit.theta_zr = theta_zr;
00609
00610 slatHit.slatIndex = *idErasedIter;
00611
00612
00613 StThreeVectorD averageHitPos(0,0,0);
00614 for (unsigned int i=0;i<mMaxSlatLayers;i++)
00615 if (layer[i]){
00616 averageHitPos += hitAtLayer[i]/numberOfHitLayers;
00617 slatHit.layerHitPositions.push_back(hitAtLayer[i]);
00618 }
00619 slatHit.hitPosition = averageHitPos;
00620
00621 slatHitVec.push_back(slatHit);
00622 slatHit.layerHitPositions.clear();
00623 }
00624
00625 idErasedVec.erase(idErasedIter);
00626 idErasedIter--;
00627 }
00628 idErasedIter++;
00629 }
00630 slatIdVec = idErasedVec;
00631 }
00632 return slatHitVec;
00633 }
00634
00635
00637 float StTofGeometry::slatHitPosition(StThreeVectorD* hitPoint){
00638 int slatId = this->tofSlatCrossId(*hitPoint);
00639 if (slatId>0){
00640 float zmin = this->tofSlat(slatId).z_min;
00641 float zmax = this->tofSlat(slatId).z_max;
00642 float cosang = this->tofSlat(slatId).cosang;
00643
00644
00645 float length = (zmax-hitPoint->z())/cosang ;
00646 float max_distance = (zmax-zmin)/cosang ;
00647
00648 if (length>max_distance || length<0){
00649 LOG_INFO << "HitPositionCorrection: length="<<length<<" max="<<max_distance
00650 << " zmin="<<zmin<<" zmax="<<zmax<<" cosang="<<cosang<<endm;
00651 this->printSlat(slatId);
00652 }
00653
00654 return length;
00655 }
00656
00657
00658 LOG_INFO << "slatHitPosition: hit point out of range;"
00659 << " phi=" << hitPoint->phi() << " z=" << hitPoint->z() << endm;
00660 return -100.;
00661 }
00662
00663
00665 idVector StTofGeometry::slatNeighbours(int slatId){
00666 tofSlatGeom_st middleSlat = this->tofSlat(slatId);
00667 int slatPhi = middleSlat.iphi;
00668 int slatEta = middleSlat.ieta;
00669
00670 idVector neighbours;
00671 neighbours.clear();
00672 int thisSlatId, ieta, iphi;
00673
00674
00675 if (slatEta<9){
00676 for (ieta= slatEta-1;ieta<=slatEta+1;ieta++)
00677 for (iphi= slatPhi-1;iphi<=slatPhi+1;iphi++){
00678 if ((ieta==slatEta) && (iphi==slatPhi)) continue;
00679 if ((ieta <1) || (iphi<1) || (iphi>4)) continue;
00680 thisSlatId = calcSlatId(iphi,ieta);
00681 neighbours.push_back(thisSlatId);
00682 }
00683 }
00684
00685 else if (slatEta==9){
00686 for (ieta= slatEta-1;ieta<=slatEta;ieta++)
00687 for (iphi= slatPhi-1;iphi<=slatPhi+1;iphi++){
00688 if ((ieta==slatEta) && (iphi==slatPhi)) continue;
00689 if ((iphi<1) || (iphi>4)) continue;
00690 thisSlatId = calcSlatId(iphi,ieta);
00691 neighbours.push_back(thisSlatId);
00692 }
00693 ieta = slatEta+1; iphi = slatPhi;
00694 thisSlatId = calcSlatId(iphi,ieta);
00695 neighbours.push_back(thisSlatId);
00696 iphi = slatPhi+1;
00697 thisSlatId = calcSlatId(iphi,ieta);
00698 neighbours.push_back(thisSlatId);
00699 }
00700
00701 else if (slatEta==10){
00702 for (iphi= slatPhi-1;iphi<=slatPhi+1;iphi++){
00703 if (iphi==slatPhi) continue;
00704 if ((iphi<1) || (iphi>5)) continue;
00705 ieta = slatEta;
00706 thisSlatId = calcSlatId(iphi,ieta);
00707 neighbours.push_back(thisSlatId);
00708 }
00709 ieta=slatEta-1;
00710 for (iphi=slatPhi-1;iphi<=slatPhi;iphi++){
00711 if ((iphi<1) || (iphi>4)) continue;
00712 thisSlatId = calcSlatId(iphi,ieta);
00713 neighbours.push_back(thisSlatId);
00714 }
00715 }
00716
00717 return neighbours;
00718 }
00719
00720
00722 idVector StTofGeometry::slatNeighboursWide(int slatId){
00723
00724 tofSlatGeom_st middleSlat = this->tofSlat(slatId);
00725 int slatPhi = middleSlat.iphi;
00726 int slatEta = middleSlat.ieta;
00727
00728 idVector neighbours;
00729 neighbours.clear();
00730 int thisSlatId, ieta, iphi;
00731
00732 for (ieta=slatEta-2;ieta<=slatEta+2;ieta++){
00733 if ((ieta<1) || (ieta>10)) continue;
00734 int iphiMax = (ieta==10)?5:4;
00735 for (iphi = 1; iphi<=iphiMax; iphi++){
00736 if ((ieta==slatEta) && (iphi==slatPhi)) continue;
00737 thisSlatId = calcSlatId(iphi,ieta);
00738 neighbours.push_back(thisSlatId);
00739 }
00740 }
00741
00742 return neighbours;
00743 }
00744
00745
00747 float StTofGeometry::slatPhiPosition(StThreeVectorD* hitPoint){
00748 int slatId = this->tofSlatCrossId(*hitPoint);
00749 if (slatId>0){
00750 float phimin = this->tofSlat(slatId).phi_min;
00751 float phimax = this->tofSlat(slatId).phi_max;
00752
00753 float philoc = hitPoint->phi() - phimin;
00754 float max_deltaphi = (phimax-phimin);
00755
00756 if (philoc>max_deltaphi || philoc<0){
00757 LOG_INFO << "slatHitPhiPosition: phi="<<philoc<<" max="<<max_deltaphi
00758 << " phimin="<<phimin<<" phimax="<<phimax<<endm;
00759 this->printSlat(slatId);
00760 }
00761
00762 return philoc;
00763 }
00764
00765
00766 LOG_INFO << "slatHitPhiPosition: hit point out of range;"
00767 << " phi=" << hitPoint->phi() << " z=" << hitPoint->z() << endm;
00768 return -100.;
00769 }
00770
00771
00772
00773 Bool_t StTofGeometry::projTrayVector(const StHelixD &helix, idVector &trayVec) const {
00774
00775 trayVec.clear();
00776 double R_tof = 220.;
00777 double res = 5.0;
00778 double s1 = helix.pathLength(R_tof).first;
00779 if(s1<0.) s1 = helix.pathLength(R_tof).second;
00780 StThreeVectorD point = helix.at(s1);
00781 double phi = point.phi()*180/3.14159;
00782
00783
00784 int itray_east = (255+(int)phi)%360/6+61;
00785 trayVec.push_back(itray_east);
00786
00787 int itray_east1 = (255+(int)(phi+res))%360/6+61;
00788 int itray_east2 = (255+(int)(phi-res))%360/6+61;
00789 if(itray_east1!=itray_east) {
00790 trayVec.push_back(itray_east1);
00791 }
00792 if(itray_east2!=itray_east&&itray_east2!=itray_east1) {
00793 trayVec.push_back(itray_east2);
00794 }
00795
00796
00797 int itray_west = (435-(int)phi)%360/6+1;
00798 trayVec.push_back(itray_west);
00799
00800 int itray_west1 = (435-(int)(phi+res))%360/6+1;
00801 int itray_west2 = (435-(int)(phi-res))%360/6+1;
00802 if(itray_west1!=itray_west) {
00803 trayVec.push_back(itray_west1);
00804 }
00805 if(itray_west2!=itray_west&&itray_west2!=itray_west1) {
00806 trayVec.push_back(itray_west2);
00807 }
00808
00809
00810
00811
00812
00813
00814
00815 if(trayVec.size()>0) return kTRUE;
00816 else return kFALSE;
00817 }
