00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "EmcRecGeoSvc/EmcRecROOTGeo.h"
00010 #include "EmcRecGeoSvc/EmcRecCrystal.h"
00011 #include "ROOTGeo/EmcROOTGeo.h"
00012 #include "CLHEP/Vector/Rotation.h"
00013
00014 #include <TGeoManager.h>
00015 #include "TGeoArb8.h"
00016
00017 using namespace std;
00018 using namespace CLHEP;
00019
00020 EmcRecROOTGeo::EmcRecROOTGeo()
00021 {
00022 m_crystalMap.clear();
00023 }
00024
00025 EmcRecROOTGeo::~EmcRecROOTGeo()
00026 {
00027 m_crystalMap.clear();
00028 }
00029
00030 void EmcRecROOTGeo::InitFromXML()
00031 {
00032 if(!gGeoManager) gGeoManager = new TGeoManager("BesGeo", "Bes geometry");
00033
00034 EmcROOTGeo *fEmcROOTGeo = new EmcROOTGeo();
00035 fEmcROOTGeo->SetPhysicalNode();
00036
00037
00038 for (int part = 0; part < fEmcROOTGeo->GetPartNb(); part++) {
00039 for (int phi = 0; phi < fEmcROOTGeo->GetPhiNb(part); phi++) {
00040 for (int theta = 0; theta < fEmcROOTGeo->GetThetaNb(part); theta++) {
00041 TGeoPhysicalNode *crystalPhysicalNode
00042 = fEmcROOTGeo->GetPhysicalCrystal(part,phi,theta);
00043
00044 TGeoArb8 *arb = (TGeoArb8*)crystalPhysicalNode->GetShape();
00045
00046
00047 double *pRot;
00048 pRot = crystalPhysicalNode->GetMatrix()->GetRotationMatrix();
00049 Hep3Vector rotX(pRot[0], pRot[3], pRot[6]);
00050 Hep3Vector rotY(pRot[1], pRot[4], pRot[7]);
00051 Hep3Vector rotZ(pRot[2], pRot[5], pRot[8]);
00052 HepRotation rotateCrystalToGlobal(rotX, rotY, rotZ);
00053
00054
00055 Hep3Vector rotTX(pRot[0], pRot[1], pRot[2]);
00056 Hep3Vector rotTY(pRot[3], pRot[4], pRot[5]);
00057 Hep3Vector rotTZ(pRot[6], pRot[7], pRot[8]);
00058 HepRotation rotateGlobalToCrystal(rotTX, rotTY, rotTZ);
00059
00060 double *pTrans;
00061 pTrans = crystalPhysicalNode->GetMatrix()->GetTranslation();
00062 Hep3Vector translateCrystalToGlobal(pTrans[0], pTrans[1], pTrans[2]);
00063
00064 Hep3Vector fPnt[10];
00065 for(int i=0;i<8;i++) {
00066 if(i<4) {
00067 fPnt[i] = Hep3Vector(arb->GetVertices()[i*2],
00068 arb->GetVertices()[i*2+1],
00069 -arb->GetDz());
00070 } else {
00071 fPnt[i] = Hep3Vector(arb->GetVertices()[i*2],
00072 arb->GetVertices()[i*2+1],
00073 arb->GetDz());
00074 }
00075
00076 fPnt[i] *= rotateCrystalToGlobal;
00077 fPnt[i] += translateCrystalToGlobal;
00078 }
00079
00080 EmcRecCrystal aCrystal;
00081 aCrystal.Type(EmcRecCrystal::SixPlane());
00082
00083 Hep3Vector tempVec;
00084 if(part==1) {
00085 tempVec = fPnt[3];
00086 fPnt[3] = fPnt[0];
00087 fPnt[0] = fPnt[1];
00088 fPnt[1] = fPnt[2];
00089 fPnt[2] = tempVec;
00090 tempVec = fPnt[7];
00091 fPnt[7] = fPnt[4];
00092 fPnt[4] = fPnt[5];
00093 fPnt[5] = fPnt[6];
00094 fPnt[6] = tempVec;
00095
00096 for(int i=0;i<8;i++) {
00097 aCrystal.Set(i,fPnt[i]);
00098 }
00099 FillCrystalMap(aCrystal,part,theta,phi);
00100
00101 } else {
00102 for(int i=0;i<8;i++) {
00103 if(i%2==0) {
00104 tempVec = fPnt[i];
00105 fPnt[i] = fPnt[i+1];
00106 fPnt[i+1] = tempVec;
00107 }
00108 aCrystal.Set(i,fPnt[i]);
00109 }
00110
00111 int sector, nb;
00112 if(theta<30) {
00113 sector = phi;
00114 nb = theta;
00115 } else {
00116 sector = phi;
00117 aCrystal.Type(EmcRecCrystal::SevenPlane());
00118 if(theta>=30&&theta<32) {
00119 nb = theta-25;
00120 } else {
00121 nb = theta-18;
00122 }
00123 int newTheta, newPhi;
00124 ComputeThetaPhi(part,nb,sector,newTheta,newPhi);
00125 Identifier id = EmcID::crystal_id(part,newTheta,newPhi);
00126 EmcRecCrystal tempCrystal = m_crystalMap[id];
00127
00128 aCrystal.Set(0,tempCrystal.Get(0));
00129 aCrystal.Set(1,fPnt[0]);
00130 aCrystal.Set(2,fPnt[1]);
00131 aCrystal.Set(3,tempCrystal.Get(2));
00132 aCrystal.Set(4,tempCrystal.Get(3));
00133 aCrystal.Set(5,tempCrystal.Get(4));
00134 aCrystal.Set(6,fPnt[4]);
00135 aCrystal.Set(7,fPnt[5]);
00136 aCrystal.Set(8,tempCrystal.Get(6));
00137 aCrystal.Set(9,tempCrystal.Get(7));
00138 }
00139
00140 FillCrystalMap(aCrystal,part,nb,sector);
00141 }
00142 }
00143 }
00144 }
00145
00146 delete fEmcROOTGeo;
00147 }
00148
00149 void EmcRecROOTGeo::FillCrystalMap(EmcRecCrystal& aCrystal,
00150 const int part, const int theta, const int phi)
00151 {
00152 Identifier id;
00153 if(part==1) {
00154 id = EmcID::crystal_id(EmcID::getBARREL(),theta,phi);
00155 m_crystalMap[id]=aCrystal;
00156 } else {
00157 int newTheta, newPhi;
00158 ComputeThetaPhi(part,theta,phi,newTheta,newPhi);
00159 id = EmcID::crystal_id(part,newTheta,newPhi);
00160 m_crystalMap[id]=aCrystal;
00161 }
00162 }
00163
00164 EmcRecCrystal EmcRecROOTGeo::GetCrystal(const Identifier& id) const
00165 {
00166 return m_crystalMap.find(id)->second;
00167 }
00168
00169 HepPoint3D EmcRecROOTGeo::GetCCenter(const Identifier& id) const
00170 {
00171 return GetCrystal(id).Center();
00172 }
00173
00174 HepPoint3D EmcRecROOTGeo::GetCFrontCenter(const Identifier& id) const
00175 {
00176 return GetCrystal(id).FrontCenter();
00177 }
00178
00179 void EmcRecROOTGeo::ComputeThetaPhi(const int part,
00180 const int theta, const int phi, int& newTheta, int& newPhi)
00181 {
00182 int sector = phi;
00183 if((sector>=0)&&(sector<3))
00184 sector+=16;
00185
00186 {
00187 if((theta>=0)&&(theta<4))
00188 {
00189 newTheta = 0;
00190 newPhi = (sector-3)*4+theta;
00191 }
00192 else if((theta>=4)&&(theta<8))
00193 {
00194 newTheta = 1;
00195 newPhi = (sector-3)*4+theta-4;
00196 }
00197 else if((theta>=8)&&(theta<13))
00198 {
00199 newTheta = 2;
00200 newPhi = (sector-3)*5+theta-8;
00201 }
00202 else if((theta>=13)&&(theta<18))
00203 {
00204 newTheta = 3;
00205 newPhi = (sector-3)*5+theta-13;
00206 }
00207 else if((theta>=18)&&(theta<24))
00208 {
00209 newTheta = 4;
00210 newPhi = (sector-3)*6+theta-18;
00211 }
00212 else if((theta>=24)&&(theta<30))
00213 {
00214 newTheta = 5;
00215 newPhi = (sector-3)*6+theta-24;
00216 }
00217 }
00218
00219 if(part==2)
00220 {
00221 switch(newTheta)
00222 {
00223 case 0:
00224 if(newPhi<32)
00225 newPhi = 31-newPhi;
00226 else
00227 newPhi = 95-newPhi;
00228 break;
00229 case 1:
00230 if(newPhi<32)
00231 newPhi = 31-newPhi;
00232 else
00233 newPhi = 95-newPhi;
00234 break;
00235 case 2:
00236 if(newPhi<40)
00237 newPhi = 39-newPhi;
00238 else
00239 newPhi = 119-newPhi;
00240 break;
00241 case 3:
00242 if(newPhi<40)
00243 newPhi = 39-newPhi;
00244 else
00245 newPhi = 119-newPhi;
00246 break;
00247 case 4:
00248 if(newPhi<48)
00249 newPhi = 47-newPhi;
00250 else
00251 newPhi = 143-newPhi;
00252 break;
00253 case 5:
00254 if(newPhi<48)
00255 newPhi = 47-newPhi;
00256 else
00257 newPhi = 143-newPhi;
00258 default:
00259 ;
00260 }
00261 }
00262
00263 }