00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <iostream>
00015 #include <fstream>
00016
00017 #include <CLHEP/Units/PhysicalConstants.h>
00018
00019 #include "GaudiKernel/Bootstrap.h"
00020 #include "GaudiKernel/IService.h"
00021 #include "GaudiKernel/Service.h"
00022 #include "GaudiKernel/ISvcLocator.h"
00023 #include "MdcGeomSvc/MdcGeomSvc.h"
00024
00025 #include "BesMdcGeoParameter.hh"
00026 #include "globals.hh"
00027 #include <cstdlib>
00028 #include "ReadBoostRoot.hh"
00029 #include "G4Svc/IG4Svc.h"
00030 #include "G4Svc/G4Svc.h"
00031
00032 BesMdcGeoParameter * BesMdcGeoParameter::fPointer=0;
00033
00034 BesMdcGeoParameter * BesMdcGeoParameter::GetGeo(void){
00035 if (! fPointer) fPointer = new BesMdcGeoParameter();
00036 return fPointer;
00037 }
00038
00039 BesMdcWire::BesMdcWire(double length, double phi, double r,double rotateAngle){
00040 fLength=length;
00041 if(phi<0){
00042 fPhi = phi + 2*pi;
00043 }else if(phi>=2*pi){
00044 fPhi = phi - 2*pi;
00045 }else{
00046 fPhi=phi;
00047 }
00048 fRadius=r;
00049 fRotateAngle=rotateAngle;
00050
00051 fX=r*cos(phi);
00052 fY=r*sin(phi);
00053 }
00054
00055 double BesMdcWire:: Phi(double z) const{
00056
00057
00058 double OB = R()*sin(RotateAngle());
00059 double OC = OB*z*2./fLength;
00060 double phi=fPhi+RotateAngle()-atan2(OC,R()*cos(RotateAngle()));
00061
00062
00063 if(phi<0){
00064 phi += 2*pi;
00065 }else if(phi>=2*pi){
00066 phi -= 2*pi;
00067 }
00068 return phi;
00069 }
00070
00071 double BesMdcWire::X(double){
00072 return fX;
00073 }
00074 double BesMdcWire::Y(double){
00075 return fY;
00076 }
00077
00078 BesMdcGeoParameter::BesMdcGeoParameter(void){
00079 ISvcLocator* svcLocator = Gaudi::svcLocator();
00080 IG4Svc* tmpSvc;
00081 G4Svc* m_G4Svc;
00082 StatusCode sc=svcLocator->service("G4Svc", tmpSvc);
00083 m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
00084 if (!sc.isSuccess())
00085 std::cout<<"BesMdcGeoParameter::Could not open G4 Service"<<std::endl;
00086 if(m_G4Svc->GetMdcDataInput()== 0){
00087 cout<<"------- get BesMdcGeoParameter from file --------"<<endl;
00088 InitFromFile();
00089 }
00090 if(m_G4Svc->GetMdcDataInput()== 1) {
00091 cout<<"=======get BesMdcGeoParameter from MdcGeomSvc======="<<endl;
00092 InitFromSvc();
00093 }
00094
00095
00096 if(fPointer)
00097 { G4Exception("BesMdcGeoParameter constructed twice."); }
00098 fPointer=this;
00099 }
00100
00101
00102 BesMdcWire BesMdcGeoParameter::Wire(int wireNo){
00103
00104 int i;
00105 for(i=0; i<fLayerNo; i++){
00106 if(fLayer[i].BeginWireNo()<=wireNo && wireNo<fLayer[i].SumWireNo()){
00107
00108 break;
00109 }
00110 }
00111
00112 BesMdcWire temp(fLayer[i].Length(), fWirePhi[wireNo], fLayer[i].R(), fLayer[i].RotateAngle());
00113 return temp;
00114 }
00115
00116
00117
00118 BesMdcWire BesMdcGeoParameter::SignalWire(int signalLayerNo, int wireNo)
00119 {
00120
00121 int i=fSignalLayer[signalLayerNo];
00122 int wireNoInLayer=2*wireNo+1-fLayer[i].FirstWire();
00123 double phi=fLayer[i].Phi();
00124 double shiftPhi=fLayer[i].ShiftPhi();
00125 double wirePhi;
00126 wirePhi= wireNoInLayer*shiftPhi+phi;
00127
00128 BesMdcWire temp(fLayer[i].Length(), fWirePhi[fLayer[i].BeginWireNo()+wireNoInLayer], fLayer[i].R(),fLayer[i].RotateAngle());
00129 return temp;
00130 }
00131
00132
00133 const BesMdcLayer& BesMdcGeoParameter::Layer(int layerNumber) const {
00134 if(layerNumber<0 || layerNumber>89){
00135 cout<<"Error: Wrong layerNo: "<<layerNumber<<endl;
00136 }
00137 return fLayer[layerNumber];
00138 }
00139
00140 const BesMdcLayer& BesMdcGeoParameter::SignalLayer(int layerNumber) const {
00141 if(layerNumber<0 || layerNumber>42){
00142 cout<<"Error: Wrong SignallayerNo: "<<layerNumber<<endl;
00143 }
00144 return fLayer[fSignalLayer[layerNumber]];
00145 }
00146
00147
00148 void BesMdcGeoParameter::InitFromFile(void){
00149 int wireNo, firstWire;
00150 double length, phi, r, rotateCell,rotateAngle;
00151 double innerR, outR, z;
00152 string name, line;
00153
00154 G4String geoPath = ReadBoostRoot::GetBoostRoot();
00155 if(!geoPath){
00156 G4Exception("BOOST environment not set!");
00157 }
00158 geoPath += "/dat/Mdc.txt";
00159
00160 ifstream inFile(geoPath);
00161 if(!inFile.good()){
00162 cout<<"Error, mdc parameters file not exist"<<endl;
00163 return;
00164 }
00165
00166 getline(inFile, line);
00167 inFile>>fLayerNo>>fWireNo>>fSignalLayerNo>>fSignalWireR>>fFieldWireR;
00168
00169 inFile.seekg(1,ios::cur);
00170 getline(inFile, line);
00171 int i,signalLayer;
00172 for(i=0; i<fSignalLayerNo; i++){
00173 inFile>>signalLayer;
00174 fSignalLayer[i]=signalLayer-1;
00175 }
00176
00177 inFile.seekg(1,ios::cur);
00178 getline(inFile, line);
00179 getline(inFile, line);
00180 for( i=0; i<fLayerNo; i++){
00181 inFile>>name>>wireNo>>length>>r>>phi>>firstWire>>rotateCell;
00182 getline(inFile, line);
00183
00184 rotateAngle=2*pi*rotateCell/wireNo;
00185
00186 fLayer[i].SetName(name);fLayer[i].SetRadius(r);
00187 fLayer[i].SetLength(length); fLayer[i].SetRotateCell(rotateCell);
00188 fLayer[i].SetRotateAngle(rotateAngle); fLayer[i].SetWireNo(wireNo);
00189 fLayer[i].SetShiftPhi(twopi/wireNo); fLayer[i].SetFirstWire(firstWire);
00190
00191 phi*=(pi/180);
00192 if(phi<0)phi += fLayer[i].ShiftPhi();
00193 fLayer[i].SetPhi(phi);
00194
00195 if(i==0){
00196 fLayer[i].SetSumWireNo(wireNo); fLayer[i].SetBeginWireNo(0);
00197 }else{
00198 fLayer[i].SetBeginWireNo(fLayer[i-1].SumWireNo());
00199 fLayer[i].SetSumWireNo(fLayer[i-1].SumWireNo()+wireNo);
00200 }
00201
00202 for(int j=0; j<wireNo; j++){
00203 fWirePhi[fLayer[i].BeginWireNo()+j]=j*fLayer[i].ShiftPhi()+phi;
00204 }
00205 }
00206
00207 if(fLayer[fLayerNo-1].SumWireNo()!= fWireNo){
00208 cout<<"Total wire number is not consistant!"<<endl;
00209 }
00210
00211 getline(inFile, line);
00212 inFile>>fSegmentNo;
00213 inFile.seekg(1,ios::cur);
00214 getline(inFile, line);
00215 getline(inFile, line);
00216
00217 for(i=0; i<fSegmentNo; i++){
00218 inFile>>length>>innerR>>outR>>z>>name;
00219 getline(inFile,line);
00220
00221 fMdcSegment[i].SetLength(length); fMdcSegment[i].SetInnerR(innerR);
00222 fMdcSegment[i].SetOutR(outR); fMdcSegment[i].SetZ(z);
00223 fMdcSegment[i].SetName(name);
00224 }
00225
00226 }
00227
00228
00229 void BesMdcGeoParameter::InitFromSvc() {
00230 ISvcLocator* svcLocator = Gaudi::svcLocator();
00231 IMdcGeomSvc* ISvc;
00232 MdcGeomSvc* mdcGeomSvc;
00233 StatusCode sc=svcLocator->service("MdcGeomSvc", ISvc);
00234 mdcGeomSvc=dynamic_cast<MdcGeomSvc *>(ISvc);
00235 if (!sc.isSuccess())
00236 std::cout<<"BesMdcGeoParameter::Could not open Geometry Service"<<std::endl;
00237
00238 fLayerNo= mdcGeomSvc->Misc()->LayerNo();
00239 fWireNo=mdcGeomSvc->Misc()->WireNo();
00240 fSignalLayerNo=mdcGeomSvc->Misc()->SLayerNo();
00241 fSignalWireR=mdcGeomSvc->Misc()->SWireR();
00242 fFieldWireR=mdcGeomSvc->Misc()->FWireR();
00243
00244 int i,signalLayer;
00245 for(i=0; i<fSignalLayerNo; i++){
00246 signalLayer=mdcGeomSvc->Layer(i)->SLayer();
00247 fSignalLayer[i]=signalLayer-1;
00248 }
00249
00250 string name;
00251 int wireNo,firstWire;
00252 double length, r, phi,rotateCell,rotateAngle;
00253 for(i=0;i<fLayerNo;i++){
00254 name=mdcGeomSvc->GeneralLayer(i)->LayerName();
00255 wireNo=mdcGeomSvc->GeneralLayer(i)->NCell()*2;
00256 length= mdcGeomSvc->GeneralLayer(i)->Length();
00257 r= mdcGeomSvc->GeneralLayer(i)->Radius();
00258 phi=mdcGeomSvc->GeneralLayer(i)->nomPhi();
00259 firstWire=mdcGeomSvc->GeneralLayer(i)->First();
00260 rotateCell= mdcGeomSvc->GeneralLayer(i)->nomShift();
00261
00262 rotateAngle=2*pi*rotateCell/wireNo;
00263
00264 fLayer[i].SetName(name);fLayer[i].SetRadius(r);
00265 fLayer[i].SetLength(length); fLayer[i].SetRotateCell(rotateCell);
00266 fLayer[i].SetRotateAngle(rotateAngle); fLayer[i].SetWireNo(wireNo);
00267 fLayer[i].SetShiftPhi(twopi/wireNo); fLayer[i].SetFirstWire(firstWire);
00268 fLayer[i].SetPhi(phi);
00269
00270 if(i==0){
00271 fLayer[i].SetSumWireNo(wireNo); fLayer[i].SetBeginWireNo(0);
00272 }else{
00273 fLayer[i].SetBeginWireNo(fLayer[i-1].SumWireNo());
00274 fLayer[i].SetSumWireNo(fLayer[i-1].SumWireNo()+wireNo);
00275 }
00276
00277 for(int j=0; j<wireNo; j++){
00278 fWirePhi[fLayer[i].BeginWireNo()+j]=j*fLayer[i].ShiftPhi()+phi;
00279 }
00280 }
00281
00282 if(fLayer[fLayerNo-1].SumWireNo()!= fWireNo){
00283 cout<<"Total wire number is not consistant!"<<endl;
00284 }
00285
00286
00287 double innerR,outR,z;
00288 fSegmentNo=mdcGeomSvc->getSegmentNo();
00289 for(i=0;i<fSegmentNo;i++){
00290 length=mdcGeomSvc->End(i)->Length();
00291 innerR=mdcGeomSvc->End(i)->InnerR();
00292 outR=mdcGeomSvc->End(i)->OutR();
00293 z=mdcGeomSvc->End(i)->Z();
00294 name=mdcGeomSvc->End(i)->Name();
00295
00296 fMdcSegment[i].SetLength(length); fMdcSegment[i].SetInnerR(innerR);
00297 fMdcSegment[i].SetOutR(outR); fMdcSegment[i].SetZ(z);
00298 fMdcSegment[i].SetName(name);
00299 }
00300
00301 }
00302
00303 void BesMdcGeoParameter::Dump() {
00304
00305 cout<<" fLayerNo: "<<fLayerNo<<endl;
00306 cout<<" fWireNo: "<<fWireNo<<endl;
00307 cout<<" fSignalLayerNo: "<<fSignalLayerNo<<endl;
00308 cout<<" fSignalWireR: "<<fSignalWireR<<endl;
00309 cout<<" fFieldWireR: "<<fFieldWireR<<endl;
00310
00311 cout<<"fSingalLayer:"<<endl;
00312 for(int i=0; i<fSignalLayerNo; i++){
00313 cout<<fSignalLayer[i]+1<<' '; }
00314 cout<<endl;
00315
00316 for(int i=0;i<fLayerNo;i++){
00317 cout<<"Layer["<<i<<"]: "
00318 <<" name:"<<fLayer[i].Name() <<" wireNo:"<<fLayer[i].WireNo()
00319 <<" length: "<<fLayer[i].Length() <<" r: "<<fLayer[i].R();
00320 if (i<75) cout<<" phi:"<<fLayer[i].Phi()*180/pi;
00321 else cout<<" phi:"<<(fLayer[i].Phi()-fLayer[i].ShiftPhi())*180/pi;
00322 cout<<" firstWire: "<<fLayer[i].FirstWire()
00323 <<" rotateCell: "<<fLayer[i].RotateCell()<<endl;
00324 }
00325
00326 cout<<"fSegmentNo:"<<fSegmentNo<<endl;
00327 for(int j=0;j<fSegmentNo;j++){
00328 cout<<"length:"<<fMdcSegment[j].Length()
00329 <<" innerR:"<<fMdcSegment[j].InnerR()
00330 <<" outR:"<<fMdcSegment[j].OutR()
00331 <<" z:"<<fMdcSegment[j].Z()
00332 <<" name:"<<fMdcSegment[j].Name()<<endl;
00333 }
00334
00335 }
00336