//*-- AUTHOR : Ilse Koenig //*-- Created : 14/11/03 by Ilse Koenig /////////////////////////////////////////////////////////////////////////////// // HGeomGeant3Builder // // Class to create the geometry in ROOT // /////////////////////////////////////////////////////////////////////////////// #include "hgeomgeant3builder.h" #include "hgeommedium.h" #include "hgeomnode.h" #include "hgeomtransform.h" #include "hgeomhit.h" #include "TMath.h" #include "geantdef.h" #include "cfortran.h" #include "geant321.h" #include "geafuncdec.h" // MDC Common Block for GEANT position of sensitive volume in the reference // coordinate system of the MDC and for shape of volume struct MDCDISTCONST_DEF { float mdcdx[4][7]; float mdcdy[4][7]; char mdcshape[5]; }; #define MDCDCONST COMMON_BLOCK(MDCDCONST,mdcdconst) COMMON_BLOCK_DEF(MDCDISTCONST_DEF,MDCDCONST); // SHOWER Common Block for GEANT position of sensitive volume in the reference // coordinate system of the Shower and for shape of volume struct SHOWERDISTCONST_DEF { float showerdx[3]; float showerdy[3]; char showershape[5]; }; #define SHOWERDCONST COMMON_BLOCK(SHOWERDCONST,showerdconst) COMMON_BLOCK_DEF(SHOWERDISTCONST_DEF,SHOWERDCONST); // TARGET Common Block #define SEGMAXNO 99 struct GEANTTARGET_DEF { int nSegments; int targetMednum[SEGMAXNO]; float size[SEGMAXNO][3]; float position[SEGMAXNO][3]; char shape[5]; }; #define GEATARGET COMMON_BLOCK(GEATARGET,geatarget) COMMON_BLOCK_DEF(GEANTTARGET_DEF,GEATARGET); // MEDIA Common Block #define NMAXCHAR 120 #define MEDMAXNO 100 struct GEANTMEDIA_DEF { int ndefmedia; char mediumname[MEDMAXNO][NMAXCHAR]; int mediumindex[MEDMAXNO]; }; #define GEAMEDIA COMMON_BLOCK(GEAMEDIA,geamedia) COMMON_BLOCK_DEF(GEANTMEDIA_DEF,GEAMEDIA); HGeomGeant3Builder::HGeomGeant3Builder() { // default constructor GEATARGET.nSegments=0; GEAMEDIA.ndefmedia=0; MDCDCONST.mdcshape[0]='\0'; SHOWERDCONST.showershape[0]='\0'; debugfile=0; } HGeomGeant3Builder::HGeomGeant3Builder(const char* name,const char* title) : HGeomBuilder(name,title) { // constructor with name and title GEATARGET.nSegments=0; GEAMEDIA.ndefmedia=0; MDCDCONST.mdcshape[0]='\0'; SHOWERDCONST.showershape[0]='\0'; debugfile=0; } HGeomGeant3Builder::~HGeomGeant3Builder() { // destructor closeDebugFile(); } void HGeomGeant3Builder::openDebugFile(const char* fname) { // opens an output file for debugging closeDebugFile(); if (strlen(fname)>0) { debugfile=new fstream(); debugfile->open(fname,ios::out); } } void HGeomGeant3Builder::closeDebugFile() { if (debugfile) { if(debugfile->rdbuf()->is_open()==1) debugfile->close(); delete debugfile; debugfile=0; } } Bool_t HGeomGeant3Builder::createNode(HGeomNode* volu) { // creates and positions a volume if (!volu) return kFALSE; TString nodeName(volu->GetName(),4); HGeomNode* mother=volu->getMotherNode(); if (!mother&&!volu->isTopNode()) { Error("createNode","Mother volume of %s not found\n",volu->GetName()); return kFALSE; } TString nodeNameSubStr(volu->GetName(),2); HGeomNode* cv=volu->getCopyNode(); TArrayD* dPar=0; if (!cv||!cv->isCreated()) { HGeomMedium* medium=volu->getMedium(); Int_t nMed=medium->getMediumIndex(); if (nMed<=0) nMed=createMedium(medium); if (nMed<=0) return kFALSE; dPar=volu->getParameters(); Int_t nPar=dPar->GetSize(); Float_t* fPar=new Float_t[nPar]; Int_t ivolu; for (Int_t i=0;iAt(i)); GSVOLU((char*)nodeName.Data(),(char*)volu->getShape().Data(),nMed,fPar,nPar,ivolu); if (debugfile) { *debugfile<<"GSVOLU: "<getShape()<<" "<setCreated(); if (cv) cv->setCreated(); } if (!volu->isTopNode()) { char cTmp[]="ONLY"; HGeomTransform* trans=volu->getPosition(); const HGeomRotation& rot=trans->getRotMatrix(); const HGeomVector& pos=trans->getTransVector(); Int_t nR=0; if (!((HGeomRotation&)rot).isUnitMatrix()) { Float_t rowmat[3], theta[3], phi[3], sinth, costh, sinph, cosph; Int_t rotate; const Double_t deg2rad=TMath::DegToRad(); for (int i=0;i<3;i++) { for (int j=0;j<3;j++) rowmat[j]=rot(i+j*3); GFANG(rowmat,costh,sinth,cosph,sinph,rotate); theta[i]=(float)(atan2((sinth*deg2rad),(costh*deg2rad))/deg2rad); if (theta[i]<0.0F) theta[i]+=180.0F; theta[i]=((int)(theta[i]*10000.F+0.5F))/10000.F; phi[i]=(float)(atan2((sinph*deg2rad),(cosph*deg2rad))/deg2rad); if (phi[i]<0.0F) phi[i]+=360.0F; phi[i]=((int)(phi[i]*10000.F+0.5F))/10000.F; } nR=++nRot; GSROTM(nR,theta[0],phi[0],theta[1],phi[1],theta[2],phi[2]); if (debugfile) { *debugfile<<"GSROTM: "<getCopyNo(),(char*)(volu->getMother().Data()), (Float_t)pos.getX(),(Float_t)pos.getY(),(Float_t)pos.getZ(),nR,cTmp); if (debugfile) { *debugfile<<"GSPOS: "<getCopyNo()<<" "<getMother() <<" "<<(Float_t)pos.getX() <<" "<<(Float_t)pos.getY()<<" "<<(Float_t)pos.getZ()<<" "<getParameters(); TString shape(volu->getShape()); Int_t nSeg=GEATARGET.nSegments; if (volu->getCopyNo()==1) { strcpy(GEATARGET.shape,shape.Data()); } GEATARGET.targetMednum[nSeg]=nMed; if (shape.CompareTo("BOX ")==0||shape.CompareTo("TUBE")==0 ||shape.CompareTo("ELTU")==0) { GEATARGET.size[nSeg][0]=(Float_t)(dPar->At(0)*10.); GEATARGET.size[nSeg][1]=(Float_t)(dPar->At(1)*10.); GEATARGET.size[nSeg][2]=(Float_t)(dPar->At(2)*10.); } else if (shape.CompareTo("CONE")==0) { strcpy(GEATARGET.shape,"TUBE"); // no change of COMMON block and rantarg.F GEATARGET.size[nSeg][0]=0.F; if (dPar->At(2)<=dPar->At(4)) { GEATARGET.size[nSeg][1]=(Float_t)(dPar->At(2)*10.); } else { GEATARGET.size[nSeg][1]=(Float_t)(dPar->At(4)*10.); } GEATARGET.size[nSeg][2]=(Float_t)(dPar->At(0)*10.); } else { Error("createNode","wrong shape %s for target",shape.Data()); return kFALSE; } HGeomTransform* center=volu->getShapePointer()->getCenterPosition(); HGeomTransform* tTrans=volu->getLabTransform(); if (!tTrans) { Error("createNode","No Lab transform for target %s",nodeName.Data()); return kFALSE; } HGeomVector tPos=tTrans->getTransVector(); tPos+=center->getTransVector(); GEATARGET.position[nSeg][0]=(Float_t)tPos.getX(); // in mm GEATARGET.position[nSeg][1]=(Float_t)tPos.getY(); GEATARGET.position[nSeg][2]=(Float_t)tPos.getZ(); GEATARGET.nSegments++; } return kTRUE; } Int_t HGeomGeant3Builder::createMedium(HGeomMedium* med) { // creates a material and medium if (!med) return 0; Int_t nComp=med->getNComponents(); Int_t nCompW=nComp*med->getWeightFac(); Double_t p[3]; if (nComp==1) { med->getComponent(0,p); GSMATE(++nMed,(char*)med->GetName(),(Float_t)p[0],(Float_t)p[1], (Float_t)med->getDensity(),(Float_t)med->getRadiationLength(), -1.,0,0); } else { Float_t* a=new Float_t[nComp]; Float_t* z=new Float_t[nComp]; Float_t* w=new Float_t[nComp]; for(Int_t i=0;igetComponent(i,p); a[i]=(Float_t)p[0]; z[i]=(Float_t)p[1]; w[i]=(Float_t)p[2]; } GSMIXT(++nMed,(char*)med->GetName(),a,z,med->getDensity(),nCompW,w); delete [] a; delete [] z; delete [] w; } Double_t mp[10]; med->getMediumPar(mp); GSTMED(nMed,(char*)med->GetName(),nMed,(int)mp[0],(int)mp[1],(float)mp[2], (float)mp[3],(float)mp[4],(float)mp[5],(float)mp[6], (float)mp[7],0,0); med->setMediumIndex(nMed); if (GEAMEDIA.ndefmediaGetName()); GEAMEDIA.mediumindex[GEAMEDIA.ndefmedia]=nMed; GEAMEDIA.ndefmedia++; } else Warning("createMedium", "Number of media exceeds max number of %i.\nList of media is truncated.", MEDMAXNO); Int_t nPckov=med->getNpckov(); if (nPckov>0) { Double_t pc[4]; Float_t* ppckov=new Float_t[nPckov]; Float_t* absco=new Float_t[nPckov]; Float_t* effic=new Float_t[nPckov]; Float_t* rindex=new Float_t[nPckov]; for(Int_t i=0;igetCerenkovPar(i,pc); ppckov[i]=(Float_t)pc[0]*1.e-9; absco[i]=(Float_t)pc[1]; effic[i]=(Float_t)pc[2]; rindex[i]=(Float_t)pc[3]; } GSCKOV(nMed,nPckov,ppckov,absco,effic,rindex); delete [] ppckov; delete [] absco; delete [] effic; delete [] rindex; } return nMed; } Bool_t HGeomGeant3Builder::createHit(HGeomHit* hit) { if (!hit) return kFALSE; char* cS=(char*)(hit->GetName()); char* cN=(char*)(hit->getComponentName()); Int_t nv=hit->getNv(); Int_t nh=hit->getNh(); Int_t iset, idet; if (nv<=0||nv>20||nh<=0) return kFALSE; if (nh>50) { Error("createHit(HGeomHit*)", "number of hit components too large, only 50 allowed!\n"); return kFALSE; } char* cv=hit->getChnmsv(); char chnmsv[20][5]; for (Int_t i=0;igetChnamh(); char chnamh[50][5]; for (Int_t i=0;igetNbitsv(), hit->getIdType(),1000,1000,iset,idet); GSDETH(cS,cN,nh,chnamh,hit->getNbitsh(), hit->getOrig(),hit->getFact()); TString setName(hit->getDetectorName()); if (setName.Contains("mdc")) return fillMdcCommonBlock(hit); if (setName.Contains("shower")) return fillShowerCommonBlock(hit); return kTRUE; } Bool_t HGeomGeant3Builder::fillMdcCommonBlock(HGeomHit* hit) { TString compName(hit->getComponentName()); Int_t l1=(Int_t)(compName[1]-'0')-1; Int_t l2=(Int_t)(compName[3]-'0')-1; if(compName[2]=='C'&&compName[3]=='4') l2 = 6; TString s; HGeomVector v; Bool_t rc=hit->calcRefPos(v,s); if (rc) { if (strlen(MDCDCONST.mdcshape)==0) { strcpy(MDCDCONST.mdcshape,s.Data()); } else { if (strcmp(MDCDCONST.mdcshape,s.Data())!=0) { Error("fillMdcCommonBlock", "Different shapes for sensitive volumes: %s, %s", MDCDCONST.mdcshape,s.Data()); return kFALSE; } } MDCDCONST.mdcdx[l1][l2]=(Float_t)(v.getX()/10.); MDCDCONST.mdcdy[l1][l2]=(Float_t)(v.getY()/10.); } return rc; } Bool_t HGeomGeant3Builder::fillShowerCommonBlock(HGeomHit* hit) { TString compName(hit->getComponentName()); Int_t l1=(Int_t)(compName[1]-'0')-1; TString s; HGeomVector v; Bool_t rc=hit->calcRefPos(v,s); if (rc) { if (strlen(SHOWERDCONST.showershape)==0) { strcpy(SHOWERDCONST.showershape,s.Data()); } else { if (strcmp(SHOWERDCONST.showershape,s.Data())!=0) { Error("fillShowerCommonBlock", "Different shapes for sensitive volumes: %s, %s", SHOWERDCONST.showershape,s.Data()); return kFALSE; } } SHOWERDCONST.showerdx[l1]=(Float_t)(v.getX()/10.); SHOWERDCONST.showerdy[l1]=(Float_t)(v.getY()/10.); } return rc; }