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
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124 #include "StEmcGeom.h"
00125 #include <assert.h>
00126 #include <strings.h>
00127 #include <stdlib.h>
00128 #include <TROOT.h>
00129 #include "StMaker.h"
00130 #include "StEmcUtil/others/emcInternalDef.h"
00131
00132 ClassImp(StEmcGeom)
00133
00134 StEmcGeom *StEmcGeom::mGeom[8] = {0,0,0,0,0,0,0,0};
00135
00136 const Float_t perr=0.01;
00137 Float_t rmin, rsmdEta, rsmdPhi;
00138
00139 StEmcGeom
00140 *StEmcGeom::instance(const Int_t det)
00141 {
00142 return getEmcGeom(det);
00143 }
00144 StEmcGeom
00145 *StEmcGeom::getEmcGeom(const Int_t det)
00146 {
00147 if(det>=1 && det<=4) {
00148 Int_t indDet = det - 1;
00149 if(mGeom[indDet] == 0) mGeom[indDet] = new StEmcGeom(det);
00150 return mGeom[indDet];
00151 } else return 0;
00152 }
00153
00154
00155 StEmcGeom
00156 *StEmcGeom::instance(const Char_t *cdet)
00157 {
00158 return getEmcGeom(cdet);
00159 }
00160 StEmcGeom
00161 *StEmcGeom::getEmcGeom(const Char_t *cdet)
00162 {
00163 Int_t det=getDetNumFromName(cdet);
00164 return getEmcGeom(det);
00165 }
00166
00167
00168 StEmcGeom
00169 *StEmcGeom::instance(const Int_t det, const Char_t* mode)
00170 {
00171 return getEmcGeom(det,mode);
00172 }
00173 StEmcGeom
00174 *StEmcGeom::getEmcGeom(const Int_t det, const Char_t* mode)
00175 {
00176 if(det>=1 && det<=4) {
00177 Int_t indDet = det - 1;
00178 if(mGeom[indDet] == 0) mGeom[indDet] = new StEmcGeom(det, mode);
00179 return mGeom[indDet];
00180 } else return 0;
00181 }
00182
00183
00184 StEmcGeom::StEmcGeom(const Int_t det)
00185 {
00186 initGeom(det);
00187 }
00188
00189 StEmcGeom::StEmcGeom(const Char_t *cdet)
00190 {
00191 Int_t det=getDetNumFromName(cdet);
00192 if(det) initGeom(det);
00193 }
00194
00195 Int_t
00196 StEmcGeom::getDetNumFromName(const Char_t *cdet)
00197 {
00198 Int_t det=0;
00199 if (!strcmp(cdet,"bemc")) {det=1;}
00200 else if(!strcmp(cdet,"bprs")) {det=2;}
00201 else if(!strcmp(cdet,"bsmde")){det=3;}
00202 else if(!strcmp(cdet,"bsmdp")){det=4;}
00203 else if(!strcmp(cdet,"eemc")) {det=5;}
00204 else if(!strcmp(cdet,"eprs")) {det=6;}
00205 else if(!strcmp(cdet,"esmde")){det=7;}
00206 else if(!strcmp(cdet,"esmdp")){det=8;}
00207 else {LOG_ERROR << Form(" StEmcGeom: Bad value of cdet %s ", cdet) << endm;}
00208 return det;
00209 }
00210
00211 StEmcGeom::~StEmcGeom()
00212 {
00213 mGeom[mDetector] = 0;
00214 }
00215
00216 void StEmcGeom::initGeom(const Int_t det)
00217 {
00218 mMode="default";
00219 mDetector = det;
00220 mGeantGeom = 0;
00221 mCalg = 0;
00222 mCalg_st = 0;
00223 mCalr = 0;
00224 mCalr_st = 0;
00225
00226 mPhiOffset.Set(2); mPhiStep.Set(2); mPhiBound.Set(2);
00227
00228 defineDefaultCommonConstants();
00229 defineModuleGridOnPhi();
00230
00231 mMaxAdc = 1024; if(det==1) mMaxAdc = 4096;
00232
00233
00234 switch (det){
00235 case 1:
00236 case 2:
00237 mRadius = 225.405;
00238 mYWidth = 11.174;
00239 initBEMCorBPRS();
00240 break;
00241 case 3:
00242 mRadius = 230.705;
00243 mYWidth = 11.2014*2;
00244 initBSMDE();
00245 break;
00246 case 4:
00247 mRadius = 232.742;
00248 mYWidth = 22.835;
00249 initBSMDP();
00250 break;
00251 case 5:
00252 case 6:
00253 initEEMCorEPRS();
00254 break;
00255 case 7:
00256 initESMDE();
00257 break;
00258 case 8:
00259 initESMDP();
00260 break;
00261 default:
00262 LOG_FATAL << Form(" StEmcGeom: Bad value of mDetector %i ", mDetector) << endm;
00263 assert(0);
00264 }
00265 mGeom[det-1] = this;
00266 }
00267
00268 StEmcGeom::StEmcGeom(const Int_t det, const Char_t *mode)
00269 {
00270 mMode=mode; mMode.ToLower();
00271 if(mMode == "geant"){
00272 getGeantGeometryTable();
00273 }
00274 else mMode.Append(" : wrong option !!! ");
00275
00276 if(!(mMode=="geant")){
00277 LOG_WARN << Form("<W> Something wrong(%s)=> using default",mMode.Data()) << endm;
00278 initGeom(det);
00279 return;
00280 }
00281
00282 LOG_INFO << Form("<I> Used Geant Geometry for BEMC, version %5.2f ",mCalg_st->version) << endm;
00283 mDetector = det;
00284 mPhiOffset.Set(2); mPhiStep.Set(2); mPhiBound.Set(2);
00285
00286 defineCommonConstants();
00287 defineModuleGridOnPhi();
00288
00289 Float_t currentDepth;
00290 switch (det){
00291 case 1:
00292 case 2:
00293 mRadius = rmin;
00294
00295 currentDepth = mRadius+mCalg_st->scintthk[0]+2.*mCalg_st->abpapthk;
00296 mYWidth = currentDepth*tan(mPhiStepHalf)-mCalg_st->crackwd;
00297 initBEMCorBPRS();
00298 break;
00299 case 3:
00300 mRadius = rsmdEta;
00301 mYWidth = mCalg_st->smalfwdh*2.;
00302 initBSMDE();
00303 break;
00304 case 4:
00305 mRadius = rsmdPhi;
00306 mYWidth = mCalg_st->smalfwdh*2.;
00307 initBSMDP();
00308 break;
00309 case 5:
00310 case 6:
00311 initEEMCorEPRS();
00312 break;
00313 case 7:
00314 initESMDE();
00315 break;
00316 case 8:
00317 initESMDP();
00318 break;
00319 default:
00320 LOG_ERROR << Form(" StEmcGeom: Bad value of mDetector %i ", mDetector) << endm;
00321 }
00322 }
00323
00324 void
00325 StEmcGeom::defineDefaultCommonConstants()
00326 {
00327
00328 mNModule = 120;
00329 mEtaMax = 0.984;
00330 mEtaMin = 0.0035;
00331
00332 mPhiOffset[0] = (75.-3.)/ 180. * C_PI;
00333 mPhiOffset[1] = (105.+3.)/180. * C_PI;
00334 mPhiStep[0] = -C_2PI /(mNModule/2);
00335 mPhiStep[1] = -mPhiStep[0];
00336
00337 mPhiBound[0] = 75. /180. * C_PI;
00338 mPhiBound[1] = 105./180. * C_PI;
00339 mPhiStepHalf = 3. * C_PI/180.;
00340 }
00341
00342 void
00343 StEmcGeom::defineCommonConstants()
00344 {
00345 Float_t lW[2], smdW;
00346 mNModule = 120;
00347 mEtaMax = 0.984;
00348 mEtaMin = 0.0035;
00349
00350 mPhiStepHalf = 360. / (Float_t)mNModule;
00351
00352
00353
00354
00355 mPhiOffset[0] = (75.-3.)/ 180. * C_PI;
00356 mPhiOffset[1] = (105.+3.)/180. * C_PI;
00357
00358 mPhiStep[0] = -toRad(2.*mPhiStepHalf);
00359 mPhiStep[1] = -mPhiStep[0];
00360 mPhiStepHalf = toRad(mPhiStepHalf);
00361
00362 rmin = mCalr_st->rmin;
00363 Int_t nsuper=(Int_t)mCalg_st->nsuper;
00364 Int_t nsmd=(Int_t)mCalg_st->nsmd;
00365 for(Int_t i=0; i<2; i++){
00366 lW[i] = mCalg_st->scintthk[i]+mCalg_st->absorthk+2.*mCalg_st->abpapthk;
00367 }
00368
00369 rsmdEta=rmin+2.*(lW[0]*nsuper+lW[1]*(nsmd-nsuper));
00370
00371 smdW=2.*(mCalg_st->g10sbthk+mCalg_st->smalfthk+mCalg_st->abpapthk);
00372
00373 rsmdPhi=rsmdEta + 2.*smdW;
00374 }
00375
00376 void
00377 StEmcGeom::defineModuleGridOnPhi()
00378 {
00379
00380
00381
00382 Int_t mw,ew,sw;
00383 Float_t etaw=-0.1;
00384 mPhiModule.Set(mNModule);
00385 Int_t im = 0;
00386 for(Int_t i=0; i<mNModule/2; i++){
00387
00388 Double_t phiW = mPhiOffset[im] + mPhiStep[im]*i;
00389 while(phiW >= C_PI) phiW -= C_2PI;
00390 while(phiW < -C_PI) phiW += C_2PI;
00391 if(phiW > (C_PI-0.0001)) phiW = -C_PI;
00392 mPhiModule[i] = phiW;
00393 Int_t cond = getBin(mPhiModule[i], etaw, mw,ew,sw);
00394 if (!cond) mPhiModule[mw-1] = phiW;
00395 else {LOG_WARN << "<W> Something wrong in StEmcGeom::defineModuleGridOnPhi()" << endm;}
00396 }
00397 }
00398
00399 void
00400 StEmcGeom::initBEMCorBPRS()
00401 {
00402 if(mMode.CompareTo("geant") == 0) {
00403 mNEta = (Int_t)mCalg_st->netat;
00404 mNSub = (Int_t)mCalg_st->nsub;
00405 }
00406 else {mNEta = 20; mNSub= 2;}
00407
00408 mNes = mNEta * mNSub;
00409 mNRaw = mNes * mNModule;
00410 mZlocal.Set(mNEta); mEta.Set(mNEta);
00411 mYlocal.Set(mNSub); mPhi.Set(mNSub);
00412
00413
00414 mEtaB.Set(mNEta+1); Int_t i;
00415
00416 for(i=0; i<mNEta; i++) {mEtaB[i] = 0.05*i;} mEtaB[mNEta]=mEtaMax; mEtaB[0]=mEtaMin;
00417
00418 for(i=0; i< mNEta; i++){
00419 mEta[i] = (mEtaB[i+1] + mEtaB[i])/2.;
00420 mZlocal[i] = mRadius * sinh(mEta[i]);
00421 }
00422
00423
00424
00425 mYlocal.Set(mNSub);
00426 mYlocal[0] = -mYWidth/2.; mYlocal[1] = mYWidth/2.;
00427
00428 mPhi.Set(mNSub);
00429 for(Int_t i=0;i<mNSub; i++) mPhi[i] = atan2(mYlocal[i],mRadius);
00430
00431 mYB.Set(mNSub+1); mPhiB.Set(mNSub+1);
00432 mYB[0] = -mYWidth;
00433 mYB[1] = 0.0;
00434 mYB[2] = +mYWidth;
00435 for(Int_t i=0; i<mNSub+1; i++) mPhiB[i] = atan2(mYB[i],mRadius);
00436
00437 }
00438
00439 void StEmcGeom::initBSMDE(){
00440 Float_t smetawdh, seta1wdh, seta2wdh, seta12wdh;
00441 if(mMode.CompareTo("geant") == 0) {
00442 mNEta = (Int_t)mCalg_st->netfirst + (Int_t)mCalg_st->netsecon;
00443 mNSub = 1;
00444 smetawdh = mCalg_st->smetawdh;
00445 seta1wdh = mCalg_st->seta1wdh;
00446 seta2wdh = mCalg_st->seta2wdh;
00447 seta12wdh= mCalg_st->set12wdh;
00448 }
00449 else {
00450 mNEta = 150;
00451 mNSub = 1;
00452 smetawdh=0.9806;
00453 seta1wdh=0.7277;
00454 seta2wdh=0.9398;
00455 seta12wdh=0.0406;
00456 }
00457
00458 mNes = mNEta * mNSub;
00459 mNRaw = mNes * mNModule;
00460 mZlocal.Set(mNEta); mEta.Set(mNEta);
00461 mYlocal.Set(mNSub); mPhi.Set(mNSub);
00462 mEtaB.Set(mNEta+1);
00463 TArrayF zb(mNEta+1);
00464
00465
00466 Int_t i;
00467 Float_t shift1, shift2;
00468 shift1 = 2.*smetawdh + seta1wdh;
00469
00470 shift2 = shift1 + (seta1wdh+seta12wdh)*2*74 + (seta1wdh+seta12wdh)
00471 + (seta2wdh+seta12wdh);
00472 zb[0] = 2.*smetawdh;
00473 for(i=0; i<mNEta; i++) {
00474 if(i<mNEta/2) {
00475 mZlocal[i] = shift1 + (seta1wdh+seta12wdh)*2*i;
00476 zb[i+1] = zb[i] + (seta1wdh+seta12wdh)*2;
00477 }
00478 else {
00479 mZlocal[i] = shift2 + (seta2wdh+seta12wdh)*2*(i-75);
00480 zb[i+1] = zb[i] + (seta2wdh+seta12wdh)*2;
00481 }
00482 mEta[i] = -::log(tan(atan2(mRadius,mZlocal[i])/2.0));
00483 mEtaB[i] = -::log(tan(atan2(mRadius,zb[i])/2.0));
00484 }
00485 mEtaB[mNEta] = -::log(tan(atan2(mRadius,zb[mNEta])/2.0));
00486
00487 mEtaMin = mEtaB[0];
00488 mEtaMax = mEtaB[mNEta];
00489
00490
00491 mYlocal.Set(mNSub); mYlocal[0] = 0.0;
00492 mPhi.Set(mNSub); mPhi[0] = 0.0;
00493
00494 mYB.Set(mNSub+1); mPhiB.Set(mNSub+1);
00495 mYB[0] = -mYWidth/2.;
00496 mYB[1] = mYWidth/2.;
00497 for(Int_t i=0; i<mNSub+1; i++) mPhiB[i] = atan2(mYB[i],mRadius);
00498 }
00499
00500 void StEmcGeom::initBSMDP()
00501 {
00502 Float_t sphiwdh, sphidwdh, shift, smdgaswdh;
00503
00504 if(mMode.CompareTo("geant") == 0) {
00505 mNEta = (Int_t)mCalg_st->netasmdp;
00506 mNSub = (Int_t)mCalg_st->nphistr;
00507 sphiwdh = mCalg_st->sphiwdh;
00508 sphidwdh = mCalg_st->sphidwdh;
00509 smdgaswdh= mCalg_st->smgaswdh;
00510
00511 }
00512 else{
00513 mNEta = 10;
00514 mNSub = 15;
00515 sphiwdh = 0.668;
00516 sphidwdh = 0.07874;
00517 smdgaswdh = 0.295;
00518 }
00519 shift = - mYWidth/2. + smdgaswdh + sphiwdh;
00520
00521 mNes = mNEta * mNSub;
00522 mNRaw = mNes * mNModule;
00523 mZlocal.Set(mNEta); mEta.Set(mNEta);
00524 mYlocal.Set(mNSub); mPhi.Set(mNSub);
00525
00526
00527 mEtaB.Set(mNEta+1); Int_t i;
00528
00529 for(i=0; i<mNEta; i++) {mEtaB[i] = 0.1*i;}
00530 mEtaB[mNEta]=mEtaMax;
00531
00532 for(i=0; i< mNEta; i++){
00533 mEta[i] = (mEtaB[i+1] + mEtaB[i])/2.;
00534 mZlocal[i] = mRadius * sinh(mEta[i]);
00535 }
00536
00537
00538
00539 mYlocal.Set(mNSub); mPhi.Set(mNSub);
00540 mYB.Set(mNSub+1); mPhiB.Set(mNSub+1);
00541 mYB[0] = -mYWidth/2. + smdgaswdh;
00542
00543 if(mMode.CompareTo("geant") == 0){
00544 Int_t n2=mNSub/2;
00545 mYlocal[n2] = 0.0;
00546 for(i=n2+1; i<mNSub; i++){
00547 mYlocal[i] = (sphiwdh+sphidwdh)*2*(i-n2);
00548 mYlocal[mNSub-i-1] = -mYlocal[i];
00549 }
00550 for(i=0; i<mNSub; i++){
00551 mPhi[i] = atan2(mYlocal[i], mRadius);
00552 }
00553 }
00554 else{
00555 for(i=0; i<mNSub; i++){
00556 mYlocal[i] = shift + (sphiwdh+sphidwdh)*2*i;
00557 mPhi[i] = atan2(mYlocal[i], mRadius);
00558 if(i==0 || i==(mNSub-1)) mYB[i+1] = mYB[i] + sphiwdh*2. + sphidwdh;
00559 else mYB[i+1] = mYB[i] + (sphiwdh + sphidwdh)*2.;
00560 mPhiB[i] = atan2(mYB[i], mRadius);
00561 }
00562 mPhiB[mNSub] = atan2(mYB[mNSub], mRadius);
00563 }
00564 }
00565
00566 Int_t StEmcGeom::getVolIdBemc(const Int_t ivid, Int_t &module,Int_t &eta,Int_t &sub, Int_t &detector)
00567 {
00568
00569
00570
00571 static Int_t emcIvid[5]={10000000,100000,100,10,1};
00572 Int_t emcChid[5], i, ividw, rl, phi, dep;
00573
00574 if(mCalg_st == 0) getGeantGeometryTable();
00575
00576 ividw = ivid;
00577 for(i=0; i<5; i++){
00578 emcChid[i] = ividw/emcIvid[i];
00579 ividw = ividw%emcIvid[i];
00580 }
00581 if(ividw == 0){
00582 rl = emcChid[0];
00583 eta = emcChid[1];
00584 phi = emcChid[2];
00585 sub = emcChid[3];
00586 dep = emcChid[4];
00587 switch (dep) {
00588 case 1:
00589 detector = BPRS; break;
00590 case 2:
00591 detector = BEMC; break;
00592 default:
00593 LOG_WARN << Form("<W> StEmcGeom::getVolIdBemc => wrong value of dep %i ",dep) << endm;
00594 }
00595 if (rl==1) {
00596
00597
00598 phi += Int_t((75.-mCalg_st->shift[0])/6.);
00599 while (phi<=0) phi+=60;
00600 while (phi>=61) phi-=60;
00601 module=phi;
00602
00603 }
00604 else if(rl==2) {
00605
00606 phi += Int_t((105.-mCalg_st->shift[1])/6.);
00607 while (phi<=0) phi+=60;
00608 while (phi>=61) phi-=60;
00609 module=phi+60;
00610 sub =(sub+1)%2+1;
00611 }
00612 else{
00613 LOG_ERROR << Form("<E> getVolIdBemc -- error decoding BEMC Geant volume Id %i; rl=%i", ivid, rl) << endm;
00614 return 1;
00615 }
00616 }
00617 else {
00618 LOG_ERROR << Form("<E> getVolIdBemc -- error decoding BEMC Geant volume Id %i=>%i", ivid, ividw) << endm;
00619 return 1;
00620 }
00621
00622
00623 return 0;
00624 }
00625
00626 Int_t StEmcGeom::getVolIdBsmd(const Int_t ivid, Int_t &module,Int_t &eta,Int_t &sub, Int_t &detector)
00627 {
00628
00629
00630 static Int_t smdIvid[5]={100000000,1000000,1000,100,1};
00631 Int_t smdChid[5], i, ividw, rl, phi, t, strip;
00632
00633 if(mCalg_st == 0) getGeantGeometryTable();
00634
00635 ividw = ivid;
00636 for(i=0; i<5; i++){
00637 smdChid[i] = ividw/smdIvid[i];
00638 ividw = ividw%smdIvid[i];
00639 }
00640 if(ividw == 0){
00641 rl = smdChid[0];
00642 eta = smdChid[1];
00643 phi = smdChid[2];
00644 t = smdChid[3];
00645 strip = smdChid[4];
00646 if (rl==1) {
00647
00648 phi += Int_t((75. - mCalg_st->shift[0])/6.);
00649 while (phi<=0) phi+=60;
00650 while (phi>=61) phi-=60;
00651 module=phi;
00652 }
00653 else if(rl==2) {
00654
00655 phi += Int_t((105. - mCalg_st->shift[1])/6.);
00656 while (phi<=0) phi+=60;
00657 while (phi>=61) phi-=60;
00658 module=phi+60;
00659 }
00660 else{
00661 LOG_ERROR << Form("<E> getVolIdBsmd -- error decoding BSMD Geant volume Id %i; rl=%i", ivid, rl) << endm;
00662 return 1;
00663 }
00664 if (t==1){
00665 detector = BSMDE;
00666 eta = strip;
00667 sub = 1;
00668 }
00669 else if(t==2){
00670 detector = BSMDE;
00671 eta = strip + 75;
00672 sub = 1;
00673 }
00674 else if(t==3){
00675 detector = BSMDP;
00676 eta = abs(eta);
00677
00678
00679
00680 switch (rl) {
00681 case 1:
00682 sub = 16 - strip;
00683 break;
00684 case 2:
00685 sub = strip;
00686 break;
00687 }
00688 }
00689 else {
00690 LOG_ERROR << Form("<E> getVolIdBsmd: Type mismatch %i ",t) << endm;
00691 return 1;
00692 }
00693 }
00694 else {
00695 LOG_ERROR << Form("<E> getVolIdBsmd -- error decoding BSMD Geant volume Id %i=>%i", ivid, ividw) << endm;
00696 return 1;
00697 }
00698
00699
00700 return 0;
00701 }
00702
00703 Int_t StEmcGeom::getVolId(const Int_t ivid, Int_t &module,Int_t &eta,Int_t &sub, Int_t &det)
00704 {
00705 if (mDetector==1||mDetector==2) return getVolIdBemc(ivid,module,eta,sub,det);
00706 else if(mDetector==3||mDetector==4) return getVolIdBsmd(ivid,module,eta,sub,det);
00707 else {
00708 LOG_ERROR << Form("<E> getVolId -- wrong detectot number %i ",mDetector) << endm;
00709 return 0;
00710 }
00711 }
00712
00713 void StEmcGeom::initEEMCorEPRS()
00714 {
00715 mNModule = 24;
00716 mNEta = 12;
00717 mNSub = 5;
00718 mNes = mNEta * mNSub;
00719 mNRaw = mNes * mNModule;
00720 mRadius = 225.405;
00721 mYWidth = 11.1716;
00722 mZlocal.Set(mNEta); mEta.Set(mNEta);
00723 mYlocal.Set(mNSub); mPhi.Set(mNSub);
00724
00725
00726 mEtaB.Set(mNEta+1); Int_t i;
00727
00728 for(i=0; i<mNEta; i++) {mEtaB[i] = 0.05*i;} mEtaB[mNEta]=0.99;
00729
00730 for(i=0; i< mNEta; i++){
00731 mEta[i] = (mEtaB[i+1] + mEtaB[i])/2.;
00732 mZlocal[i] = mRadius * sinh(mEta[i]);
00733 }
00734
00735
00736 mYlocal.Set(mNSub);
00737 mYlocal[0] = mYWidth/2; mYlocal[1] = - mYlocal[0];
00738
00739 mPhi.Set(mNSub);
00740 mPhi[0] = atan2(mYWidth/2,mRadius); mPhi[1] = -mPhi[0];
00741
00742
00743 }
00744
00745 void StEmcGeom::initESMDE(){
00746 mNModule = 24;
00747 mNEta = 200;
00748 mNSub = 200;
00749 mNes = mNEta * mNSub;
00750 mNRaw = mNes * mNModule;
00751 mRadius = 230.467;
00752 mYWidth = 11.2014*2;
00753 mZlocal.Set(mNEta); mEta.Set(mNEta);
00754 mYlocal.Set(mNSub); mPhi.Set(mNSub);
00755
00756
00757 Int_t i;
00758 Float_t smetawdh=0.9806;
00759 Float_t seta1wdh=0.7277, seta2wdh=0.9398, seta12wdh=0.0406;
00760 Float_t shift1, shift2;
00761 shift1 = 2.*smetawdh + seta1wdh;
00762
00763 shift2 = shift1 + (seta1wdh+seta12wdh)*2*74 + (seta1wdh+seta12wdh)
00764 + (seta2wdh+seta12wdh);
00765
00766 for(i=0; i<mNEta; i++) {
00767 if(i<mNEta/2) mZlocal[i] = shift1 + (seta1wdh+seta12wdh)*2*i;
00768 else mZlocal[i] = shift2 + (seta2wdh+seta12wdh)*2*(i-75);
00769 mEta[i] = -::log(tan(atan2(mRadius,mZlocal[i])/2.0));
00770 }
00771
00772
00773 mYlocal.Set(mNSub); mYlocal[0] = 0.0;
00774
00775 mPhi.Set(mNSub); mPhi[0] = 0.0;
00776
00777
00778
00779 }
00780
00781 void StEmcGeom::initESMDP()
00782 {
00783 mNModule = 24;
00784 mNEta = 200;
00785 mNSub = 200;
00786 mNes = mNEta * mNSub;
00787 mNRaw = mNes * mNModule;
00788 mRadius = 232.467;
00789 mYWidth = 22.835;
00790 mZlocal.Set(mNEta); mEta.Set(mNEta);
00791 mYlocal.Set(mNSub); mPhi.Set(mNSub);
00792
00793
00794 mEtaB.Set(mNEta+1); Int_t i;
00795
00796 for(i=0; i<mNEta; i++) {mEtaB[i] = 0.1*i;}
00797 mEtaB[mNEta]=0.99;
00798
00799 for(i=0; i< mNEta; i++){
00800 mEta[i] = (mEtaB[i+1] + mEtaB[i])/2.;
00801 mZlocal[i] = mRadius * sinh(mEta[i]);
00802 }
00803
00804
00805
00806 mYlocal.Set(mNSub); mPhi.Set(mNSub);
00807
00808 Float_t sphiwdh = 0.668;
00809 Float_t sphidwdh = 0.07874;
00810 Float_t shift = mYWidth/2. - 0.295 - sphiwdh;
00811
00812 for(i=0; i<mNSub; i++){
00813 mYlocal[i] = shift - (sphiwdh+sphidwdh)*2*i;
00814 mPhi[i] = atan2(mYlocal[i], mRadius);
00815 }
00816 }
00817
00818 void StEmcGeom::printGeom() const
00819 {
00820 cout<<" mMode "<<mMode.Data()<<endl;
00821 cout<<" mDetector "<<mDetector<<endl;
00822 cout<<" mNModule "<<mNModule<<endl;
00823 cout<<" mNEta "<<mNEta<<endl;
00824 cout<<" mNSub "<<mNSub<<endl;
00825 cout<<" mNes "<<mNes<<endl;
00826 cout<<" mNRaw "<<mNRaw<<endl;
00827 cout<<" mRadius "<<mRadius<<endl;
00828 cout<<" mYWidth "<<mYWidth<<endl;
00829 cout<<" mEtaMin "<<mEtaMin<<endl;
00830 cout<<" mEtaMax "<<mEtaMax<<endl;
00831 cout<<" mPhiOffset "<<mPhiOffset[0]<<"("<<toDeg(mPhiOffset[0])<<") "
00832 <<mPhiOffset[1]<<"("<<toDeg(mPhiOffset[0])<<")"<<endl;
00833 cout<<" mPhiStep "<<mPhiStep[0]<<"("<<toDeg(mPhiStep[0])<<") "
00834 <<mPhiStep[1]<<"("<<toDeg(mPhiStep[1])<<")"<<endl;
00835 cout<<" Max ADC "<<mMaxAdc<<endl;
00836
00837 Int_t i;
00838 cout<<"\n Z grid and Eta grid "<<endl;
00839 cout<<" mEtaB[0] "<<mEtaB[0]<<endl;
00840 for(i=0; i<mNEta; i++){
00841 printf(" i %3i Zl %7.3f Eta %8.5f mEtaB %7.5f\n"
00842 , i, mZlocal[i], mEta[i], mEtaB[i+1]);
00843 }
00844
00845 cout<<"\n Y grid and Phi grid "<<endl;
00846
00847 printf(" mYB %7.3f mPhiB %9.5f \n", mYB[0], mPhiB[0]);
00848 for(i=0; i<mNSub; i++){
00849 printf(" i %2i Yl %7.3f Phi %9.5f ", i, mYlocal[i], mPhi[i]);
00850 printf(" mYB %7.3f mPhiB %9.5f \n", mYB[i+1], mPhiB[i+1]);
00851 }
00852 cout<<"\n Phi grid of center of modules in STAR system\n";
00853 cout<<" =============================================\n";
00854 Int_t mw,ew,sw;
00855 Float_t etaw=-0.1;
00856 for(i=0; i<mNModule/2; i++){
00857 printf(" %3i phi %10.7f (%6.1f)",
00858 i+1,mPhiModule[i],mPhiModule[i]*C_DEG_PER_RAD);
00859 Int_t cond = getBin(mPhiModule[i], etaw, mw,ew,sw);
00860 if(!cond) {
00861 printf(" => %3i phi %10.7f (%6.1f)",
00862 mw,mPhiModule[mw-1],mPhiModule[mw-1]*C_DEG_PER_RAD);
00863 }
00864 printf("\n");
00865 }
00866 cout<<"\n == "<<endl;
00867 for(Int_t i=0;i<4; i++){
00868 cout<<" Pointer for det = "<<i+1<<" -> "<<mGeom[i]<<endl;
00869 }
00870 }
00871
00872 void StEmcGeom::compare(const StEmcGeom &g, Bool_t key=kFALSE) const
00873 {
00874 Float_t err;
00875 Int_t i;
00876 if(mDetector==g.Detector()) {
00877 printf(" mMode %10s | %10s \n", mMode.Data(), g.Mode()->Data());
00878 printf("---------------------------------------------------\n");
00879 if(mNModule != g.NModule() || key)
00880 printf(" mNModule %10i | %10i\n",mNModule,g.NModule());
00881 if(mNEta != g.NEta() || key)
00882 printf(" mNEta %10i | %10i\n",mNEta,g.NEta());
00883 if(mNSub != g.NSub() || key)
00884 printf(" mNSub %10i | %10i\n",mNSub,g.NSub());
00885 if(mNes != g.Nes() || key)
00886 printf(" mNes %10i | %10i\n",mNes,g.Nes());
00887 if(mNRaw != g.NRaw() || key)
00888 printf(" mNRaw %10i | %10i\n",mNRaw,g.NRaw());
00889
00890 err=relativeError(mRadius,g.Radius());
00891 if(err>perr || key){
00892 printf(" mRadius %10.3f | %10.3f",mRadius,g.Radius());
00893 printError(err);
00894 }
00895
00896 err=relativeError(mYWidth,g.YWidth());
00897 if(err>perr || key) {
00898 printf(" mYWidth %10.3f | %10.3f",mYWidth,g.YWidth());
00899 printError(err);
00900 }
00901
00902 err=relativeError(mEtaMax,g.EtaMax());
00903 if(err>perr || key) {
00904 printf(" mEtaMax %10.3f | %10.3f",mEtaMax,g.EtaMax());
00905 printError(err);
00906 }
00907
00908 for(i=0; i<2; i++){
00909 err=relativeError(mPhiOffset[i],g.PhiOffset()[i]);
00910 if(err>perr || key){
00911 printf(" mPhiOffset[%1i] %7.3f | %10.3f",i,mPhiOffset[i],g.PhiOffset()[i]);
00912 printError(err);
00913 }
00914 }
00915 for(i=0; i<2; i++){
00916 err=relativeError(mPhiStep[i],g.PhiStep()[i]);
00917 if(err>perr || key){
00918 printf(" mPhiStep[%1i] %7.3f | %10.3f",i,mPhiStep[i],g.PhiStep()[i]);
00919 printError(err);
00920 }
00921 }
00922 printf("\n Phi grid for center of module in STAR system mNModule=%i\n",mNModule);
00923 for(i=0; i<mNModule; i++){
00924 err=relativeError(mPhiModule[i],g.PhiModule()[i]);
00925 if(err>perr || key){
00926 printf(" %3i phi %7.2f | %10.2f",i,toDeg(mPhiModule[i])
00927 ,toDeg(g.PhiModule()[i]));
00928 printError(err);
00929 }
00930 }
00931
00932 printf("\n Z grid and Eta grid => mNEta=%i\n",mNEta);
00933 for(i=0; i<mNEta; i++){
00934 err=relativeError(mZlocal[i], g.Zlocal()[i]);
00935 if(err>perr || key){
00936 printf(" %3i Zl %7.2f | %10.2f || ",i, mZlocal[i], g.Zlocal()[i]);
00937 printf("Eta %7.3f | %7.3f",mEta[i], g.Eta()[i]);
00938 printError(err);
00939 }
00940 }
00941
00942 printf("\n Y grid and Phi grid => mNSub=%i\n",mNSub);
00943 for(i=0; i<mNSub; i++){
00944 err=relativeError(mYlocal[i], g.Ylocal()[i]);
00945 if(err>perr || key){
00946 printf("%3i Yl %7.3f | %10.3f || ",i, mYlocal[i], g.Ylocal()[i]);
00947 printf("Phi %9.6f | %9.6f",mPhi[i], g.Phi()[i]);
00948 printError(err);
00949 }
00950 }
00951 }
00952 else printf("<W> You compare geometry for different detector %i != %i \n",
00953 mDetector, g.mDetector);
00954 }
00955
00956 Float_t StEmcGeom::relativeError(Float_t a, Float_t b) const
00957 {
00958 Float_t sum = fabs(a) + fabs(b), perr;
00959 if(sum == 0.0) return 0.0;
00960 else {
00961 perr = 200.*fabs(a-b)/sum;
00962 return perr;
00963 }
00964 }
00965
00966 void
00967 StEmcGeom::printError(Float_t err) const
00968 {
00969 if(err>perr) {LOG_INFO << Form(" | perr=%6.3f%% ",err) << endm;}
00970 else {LOG_INFO << " " << endm; }
00971 }
00972
00973 void
00974 StEmcGeom::getGeantGeometryTable()
00975 {
00976
00977
00978 mGeantGeom = NULL;
00979 mCalg = 0;
00980 mCalr = 0;
00981 mCalg_st = 0;
00982 mCalr_st = 0;
00983 StMaker maker;
00984 mChain = maker.GetChain();
00985
00986 if(mChain) mGeantGeom = mChain->GetDataSet(".const/geom");
00987
00988 TString line;
00989
00990 if(mGeantGeom) {
00991 mCalg = (St_calb_calg *) mGeantGeom->Find("calb_calg");
00992 if(mCalg && mCalg->GetNRows()) {
00993 mCalg_st = mCalg->GetTable();
00994 printf("calb_calr get from Geant::");
00995 for(Int_t i=0;i<2;i++) printf(" Barrel %i Angle shift %6.0f ", i+1, mCalg_st->shift[i]);
00996 printf("\n");
00997 mCalr = (St_calb_calr *) mGeantGeom->Find("calb_calr");
00998 if(mCalr) mCalr_st = mCalr->GetTable();
00999 }
01000 }
01001 if(!mCalg_st || !mCalr_st) {
01002 mMode.Append(" : No table");
01003 LOG_INFO << Form("StEmcGeom::getGeantGeometryTable() could not find geom") << endm;
01004 LOG_INFO << Form("StEmcGeom::getGeantGeometryTable() create own calb_calg/r") << endm;
01005 mCalg = new St_calb_calg("calg", 1);
01006 mCalr = new St_calb_calr("calr", 1);
01007 mCalg_st = mCalg->GetTable();
01008 mCalr_st = mCalr->GetTable();
01009
01010
01011
01012
01013 mCalg_st[0].shift[0]=75.0;
01014 mCalg_st[0].shift[1]=105.0;
01015 }
01016 }
