40 ClassImp(TGeoMCGeometry);
 
   42 TGeoMCGeometry* TGeoMCGeometry::fgInstance=0;
 
   49 TGeoMCGeometry::TGeoMCGeometry(
const char *name, 
const char *title,
 
   50                                Bool_t g3CompatibleVolumeNames)
 
   51   : TVirtualMCGeometry(name, title),
 
   52     fG3CompatibleVolumeNames(g3CompatibleVolumeNames)
 
   61 TGeoMCGeometry::TGeoMCGeometry()
 
   62   : TVirtualMCGeometry(),
 
   63     fG3CompatibleVolumeNames(kFALSE)
 
   72 TGeoMCGeometry::~TGeoMCGeometry()
 
   86 TGeoManager* TGeoMCGeometry::GetTGeoManager()
 const 
   88   if ( ! gGeoManager ) 
new TGeoManager(
"TGeo", 
"Root geometry manager");
 
   97 Double_t* TGeoMCGeometry::CreateDoubleArray(Float_t* array, Int_t size)
 const 
   99    Double_t* doubleArray;
 
  101       doubleArray = 
new Double_t[size];
 
  102       for (Int_t i=0; i<size; i++) doubleArray[i] = array[i];
 
  105       doubleArray = 
new Double_t[1];
 
  115 void TGeoMCGeometry::Vname(
const char *name, 
char *vname)
 const 
  117    if (fG3CompatibleVolumeNames) {
 
  118       Int_t l = strlen(name);
 
  121       for (i=0;i<l;i++) vname[i] = toupper(name[i]);
 
  122       for (i=l;i<4;i++) vname[i] = 
' ';
 
  125       Int_t l = strlen(name);
 
  127       for (Int_t i=0;i<l;i++) vname[i] = name[i];
 
  154 void TGeoMCGeometry::Material(Int_t& kmat, 
const char* name, Double_t a, Double_t z,
 
  155                        Double_t dens, Double_t radl, Double_t absl, Float_t* buf,
 
  158    Double_t* dbuf = CreateDoubleArray(buf, nwbuf);
 
  159    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
 
  181 void TGeoMCGeometry::Material(Int_t& kmat, 
const char* name, Double_t a, Double_t z,
 
  182                      Double_t dens, Double_t radl, Double_t absl, Double_t* ,
 
  185    GetTGeoManager()->Material(name, a, z, dens, kmat, radl, absl);
 
  202 void TGeoMCGeometry::Mixture(Int_t& kmat, 
const char* name, Float_t* a, Float_t* z,
 
  203                     Double_t dens, Int_t nlmat, Float_t* wmat)
 
  205    Double_t* da = CreateDoubleArray(a, TMath::Abs(nlmat));
 
  206    Double_t* dz = CreateDoubleArray(z, TMath::Abs(nlmat));
 
  207    Double_t* dwmat = CreateDoubleArray(wmat, TMath::Abs(nlmat));
 
  209    Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
 
  210    for (Int_t i=0; i< TMath::Abs(nlmat); i++) {
 
  211       a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
 
  233 void TGeoMCGeometry::Mixture(Int_t& kmat, 
const char* name, Double_t* a, Double_t* z,
 
  234                     Double_t dens, Int_t nlmat, Double_t* wmat)
 
  240       for (i=0;i<nlmat;i++) {
 
  241          amol += a[i]*wmat[i];
 
  243       for (i=0;i<nlmat;i++) {
 
  244          wmat[i] *= a[i]/amol;
 
  247    GetTGeoManager()->Mixture(name, a, z, dens, nlmat, wmat, kmat);
 
  273 void TGeoMCGeometry::Medium(Int_t& kmed, 
const char* name, Int_t nmat, Int_t isvol,
 
  274                    Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 
  275                    Double_t stemax, Double_t deemax, Double_t epsil,
 
  276                    Double_t stmin, Float_t* ubuf, Int_t nbuf)
 
  279    Double_t* dubuf = CreateDoubleArray(ubuf, nbuf);
 
  280    Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil,
 
  308 void TGeoMCGeometry::Medium(Int_t& kmed, 
const char* name, Int_t nmat, Int_t isvol,
 
  309                    Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 
  310                    Double_t stemax, Double_t deemax, Double_t epsil,
 
  311                    Double_t stmin, Double_t* , Int_t )
 
  313    GetTGeoManager()->Medium(name,kmed,nmat, isvol, ifield, fieldm, tmaxfd, stemax,deemax, epsil, stmin);
 
  328 void TGeoMCGeometry::Matrix(Int_t& krot, Double_t thex, Double_t phix, Double_t they,
 
  329                    Double_t phiy, Double_t thez, Double_t phiz)
 
  331    krot = GetTGeoManager()->GetListOfMatrices()->GetEntriesFast();
 
  332    GetTGeoManager()->Matrix(krot, thex, phix, they, phiy, thez, phiz);
 
  345 Int_t TGeoMCGeometry::Gsvolu(
const char *name, 
const char *shape, Int_t nmed,
 
  346                     Float_t *upar, Int_t npar)
 
  348    Double_t* dupar = CreateDoubleArray(upar, npar);
 
  349    Int_t 
id = Gsvolu(name, shape, nmed, dupar, npar);
 
  364 Int_t TGeoMCGeometry::Gsvolu(
const char *name, 
const char *shape, Int_t nmed,
 
  365                     Double_t *upar, Int_t npar)
 
  372    TGeoVolume* vol = GetTGeoManager()->Volume(vname, vshape, nmed, upar, npar);
 
  374       Fatal(
"Gsvolu", 
"Could not create volume %s", name);
 
  377    return vol->GetNumber();
 
  391 void  TGeoMCGeometry::Gsdvn(
const char *name, 
const char *mother, Int_t ndiv,
 
  397    Vname(mother,vmother);
 
  399    GetTGeoManager()->Division(vname, vmother, iaxis, ndiv, 0, 0, 0, 
"n");
 
  410 void  TGeoMCGeometry::Gsdvn2(
const char *name, 
const char *mother, Int_t ndiv,
 
  411                     Int_t iaxis, Double_t c0i, Int_t numed)
 
  416    Vname(mother,vmother);
 
  418    GetTGeoManager()->Division(vname, vmother, iaxis, ndiv, c0i, 0, numed, 
"nx");
 
  432 void  TGeoMCGeometry::Gsdvt(
const char *name, 
const char *mother, Double_t step,
 
  433                    Int_t iaxis, Int_t numed, Int_t )
 
  438    Vname(mother,vmother);
 
  440    GetTGeoManager()->Division(vname, vmother, iaxis, 0, 0, step, numed, 
"s");
 
  455 void  TGeoMCGeometry::Gsdvt2(
const char *name, 
const char *mother, Double_t step,
 
  456                     Int_t iaxis, Double_t c0, Int_t numed, Int_t )
 
  461    Vname(mother,vmother);
 
  463    GetTGeoManager()->Division(vname, vmother, iaxis, 0, c0, step, numed, 
"sx");
 
  482 void  TGeoMCGeometry::Gsord(
const char * , Int_t )
 
  502 void  TGeoMCGeometry::Gspos(
const char *name, Int_t nr, 
const char *mother, Double_t x,
 
  503                    Double_t y, Double_t z, Int_t irot, 
const char *konly)
 
  505    TString only = konly;
 
  507    Bool_t isOnly = kFALSE;
 
  508    if (only.Contains(
"only")) isOnly = kTRUE;
 
  512    Vname(mother,vmother);
 
  515    GetTGeoManager()->Node(vname, nr, vmother, x, y, z, irot, isOnly, upar);
 
  523 void  TGeoMCGeometry::Gsposp(
const char *name, Int_t nr, 
const char *mother,
 
  524                     Double_t x, Double_t y, Double_t z, Int_t irot,
 
  525                     const char *konly, Float_t *upar, Int_t np )
 
  527    Double_t* dupar = CreateDoubleArray(upar, np);
 
  528    Gsposp(name, nr, mother, x, y, z, irot, konly, dupar, np);
 
  537 void  TGeoMCGeometry::Gsposp(
const char *name, Int_t nr, 
const char *mother,
 
  538                     Double_t x, Double_t y, Double_t z, Int_t irot,
 
  539                     const char *konly, Double_t *upar, Int_t np )
 
  541    TString only = konly;
 
  543    Bool_t isOnly = kFALSE;
 
  544    if (only.Contains(
"only")) isOnly = kTRUE;
 
  548    Vname(mother,vmother);
 
  550    GetTGeoManager()->Node(vname,nr,vmother, x,y,z,irot,isOnly,upar,np);
 
  557 Int_t TGeoMCGeometry::VolId(
const char *name)
 const 
  559    Int_t uid = GetTGeoManager()->GetUID(name);
 
  561       printf(
"VolId: Volume %s not found\n",name);
 
  571 Int_t TGeoMCGeometry::MediumId(
const char *name)
 const 
  573    TGeoMedium* medium = GetTGeoManager()->GetMedium(name);
 
  574    if (medium) 
return medium->GetId();
 
  576    printf(
"MediumId: Medium %s not found\n",name);
 
  584 const char* TGeoMCGeometry::VolName(Int_t 
id)
 const 
  586    TGeoVolume *volume = GetTGeoManager()->GetVolume(
id);
 
  588       Error(
"VolName",
"volume with id=%d does not exist",
id);
 
  591    return volume->GetName();
 
  598 Int_t TGeoMCGeometry::NofVolumes()
 const 
  600    return GetTGeoManager()->GetListOfUVolumes()->GetEntriesFast()-1;
 
  607 Int_t TGeoMCGeometry::NofVolDaughters(
const char* volName)
 const 
  609    TGeoVolume* volume = GetTGeoManager()->GetVolume(volName);
 
  612       Error(
"NofVolDaughters", 
"Volume %s not found.", volName);
 
  616    return volume->GetNdaughters();
 
  623 const char*  TGeoMCGeometry::VolDaughterName(
const char* volName, Int_t i)
 const 
  626    TGeoVolume* volume = GetTGeoManager()->GetVolume(volName);
 
  628       Error(
"VolDaughterName", 
"Volume %s not found.", volName);
 
  633    if (i<0 || i>=volume->GetNdaughters()) {
 
  634       Error(
"VolDaughterName", 
"Volume %s Index out of limits", volName);
 
  639    return volume->GetNode(i)->GetVolume()->GetName();
 
  646 Int_t TGeoMCGeometry::VolDaughterCopyNo(
const char* volName, Int_t i)
 const 
  649    TGeoVolume* volume = GetTGeoManager()->GetVolume(volName);
 
  651       Error(
"VolDaughterName", 
"Volume %s not found.", volName);
 
  656    if (i<0 || i>=volume->GetNdaughters()) {
 
  657       Error(
"VolDaughterName", 
"Volume %s Index out of limits", volName);
 
  662    return volume->GetNode(i)->GetNumber();
 
  669 Int_t TGeoMCGeometry::VolId2Mate(Int_t 
id)
 const 
  671    TGeoVolume *volume = GetTGeoManager()->GetVolume(
id);
 
  673       Error(
"VolId2Mate",
"volume with id=%d does not exist",
id);
 
  676    TGeoMedium *med = volume->GetMedium();
 
  702 Bool_t TGeoMCGeometry::GetTransformation(
const TString &volumePath,TGeoHMatrix &mat)
 
  705    GetTGeoManager()->PushPath();
 
  706    if (!GetTGeoManager()->cd(volumePath.Data())) {
 
  707       GetTGeoManager()->PopPath();
 
  710    mat = *GetTGeoManager()->GetCurrentMatrix();
 
  711    GetTGeoManager()->PopPath();
 
  727 Bool_t TGeoMCGeometry::GetShape(
const TString &volumePath,TString &shapeType,
 
  731    GetTGeoManager()->PushPath();
 
  732    if (!GetTGeoManager()->cd(volumePath.Data())) {
 
  733       GetTGeoManager()->PopPath();
 
  736    TGeoVolume * vol = GetTGeoManager()->GetCurrentVolume();
 
  737    GetTGeoManager()->PopPath();
 
  738    if (!vol) 
return kFALSE;
 
  739    TGeoShape *shape = vol->GetShape();
 
  740    TClass *class_type = shape->IsA();
 
  741    if (class_type==TGeoBBox::Class()) {
 
  745       TGeoBBox *box = (TGeoBBox*)shape;
 
  746       par.AddAt(box->GetDX(),0);
 
  747       par.AddAt(box->GetDY(),1);
 
  748       par.AddAt(box->GetDZ(),2);
 
  751    if (class_type==TGeoTrd1::Class()) {
 
  755       TGeoTrd1 *trd1 = (TGeoTrd1*)shape;
 
  756       par.AddAt(trd1->GetDx1(),0);
 
  757       par.AddAt(trd1->GetDx2(),1);
 
  758       par.AddAt(trd1->GetDy(), 2);
 
  759       par.AddAt(trd1->GetDz(), 3);
 
  762    if (class_type==TGeoTrd2::Class()) {
 
  766       TGeoTrd2 *trd2 = (TGeoTrd2*)shape;
 
  767       par.AddAt(trd2->GetDx1(),0);
 
  768       par.AddAt(trd2->GetDx2(),1);
 
  769       par.AddAt(trd2->GetDy1(),2);
 
  770       par.AddAt(trd2->GetDy2(),3);
 
  771       par.AddAt(trd2->GetDz(), 4);
 
  774    if (class_type==TGeoTrap::Class()) {
 
  778       TGeoTrap *trap = (TGeoTrap*)shape;
 
  779       Double_t tth = TMath::Tan(trap->GetTheta()*TMath::DegToRad());
 
  780       par.AddAt(trap->GetDz(),0);
 
  781       par.AddAt(tth*TMath::Cos(trap->GetPhi()*TMath::DegToRad()),1);
 
  782       par.AddAt(tth*TMath::Sin(trap->GetPhi()*TMath::DegToRad()),2);
 
  783       par.AddAt(trap->GetH1(),3);
 
  784       par.AddAt(trap->GetBl1(),4);
 
  785       par.AddAt(trap->GetTl1(),5);
 
  786       par.AddAt(TMath::Tan(trap->GetAlpha1()*TMath::DegToRad()),6);
 
  787       par.AddAt(trap->GetH2(),7);
 
  788       par.AddAt(trap->GetBl2(),8);
 
  789       par.AddAt(trap->GetTl2(),9);
 
  790       par.AddAt(TMath::Tan(trap->GetAlpha2()*TMath::DegToRad()),10);
 
  793    if (class_type==TGeoTube::Class()) {
 
  797       TGeoTube *tube = (TGeoTube*)shape;
 
  798       par.AddAt(tube->GetRmin(),0);
 
  799       par.AddAt(tube->GetRmax(),1);
 
  800       par.AddAt(tube->GetDz(),2);
 
  803    if (class_type==TGeoTubeSeg::Class()) {
 
  807       TGeoTubeSeg *tubs = (TGeoTubeSeg*)shape;
 
  808       par.AddAt(tubs->GetRmin(),0);
 
  809       par.AddAt(tubs->GetRmax(),1);
 
  810       par.AddAt(tubs->GetDz(),2);
 
  811       par.AddAt(tubs->GetPhi1(),3);
 
  812       par.AddAt(tubs->GetPhi2(),4);
 
  815    if (class_type==TGeoCone::Class()) {
 
  819       TGeoCone *cone = (TGeoCone*)shape;
 
  820       par.AddAt(cone->GetDz(),0);
 
  821       par.AddAt(cone->GetRmin1(),1);
 
  822       par.AddAt(cone->GetRmax1(),2);
 
  823       par.AddAt(cone->GetRmin2(),3);
 
  824       par.AddAt(cone->GetRmax2(),4);
 
  827    if (class_type==TGeoConeSeg::Class()) {
 
  831       TGeoConeSeg *cons = (TGeoConeSeg*)shape;
 
  832       par.AddAt(cons->GetDz(),0);
 
  833       par.AddAt(cons->GetRmin1(),1);
 
  834       par.AddAt(cons->GetRmax1(),2);
 
  835       par.AddAt(cons->GetRmin2(),3);
 
  836       par.AddAt(cons->GetRmax2(),4);
 
  837       par.AddAt(cons->GetPhi1(),5);
 
  838       par.AddAt(cons->GetPhi2(),6);
 
  841    if (class_type==TGeoSphere::Class()) {
 
  845       TGeoSphere *sphe = (TGeoSphere*)shape;
 
  846       par.AddAt(sphe->GetRmin(),0);
 
  847       par.AddAt(sphe->GetRmax(),1);
 
  848       par.AddAt(sphe->GetTheta1(),2);
 
  849       par.AddAt(sphe->GetTheta2(),3);
 
  850       par.AddAt(sphe->GetPhi1(),4);
 
  851       par.AddAt(sphe->GetPhi2(),5);
 
  854    if (class_type==TGeoPara::Class()) {
 
  858       TGeoPara *para = (TGeoPara*)shape;
 
  859       par.AddAt(para->GetX(),0);
 
  860       par.AddAt(para->GetY(),1);
 
  861       par.AddAt(para->GetZ(),2);
 
  862       par.AddAt(para->GetTxy(),3);
 
  863       par.AddAt(para->GetTxz(),4);
 
  864       par.AddAt(para->GetTyz(),5);
 
  867    if (class_type==TGeoPgon::Class()) {
 
  869       TGeoPgon *pgon = (TGeoPgon*)shape;
 
  870       Int_t nz = pgon->GetNz();
 
  871       const Double_t *rmin = pgon->GetRmin();
 
  872       const Double_t *rmax = pgon->GetRmax();
 
  873       const Double_t *z = pgon->GetZ();
 
  876       par.AddAt(pgon->GetPhi1(),0);
 
  877       par.AddAt(pgon->GetDphi(),1);
 
  878       par.AddAt(pgon->GetNedges(),2);
 
  879       par.AddAt(pgon->GetNz(),3);
 
  880       for (Int_t i=0; i<nz; i++) {
 
  881          par.AddAt(z[i], 4+3*i);
 
  882          par.AddAt(rmin[i], 4+3*i+1);
 
  883          par.AddAt(rmax[i], 4+3*i+2);
 
  887    if (class_type==TGeoPcon::Class()) {
 
  889       TGeoPcon *pcon = (TGeoPcon*)shape;
 
  890       Int_t nz = pcon->GetNz();
 
  891       const Double_t *rmin = pcon->GetRmin();
 
  892       const Double_t *rmax = pcon->GetRmax();
 
  893       const Double_t *z = pcon->GetZ();
 
  896       par.AddAt(pcon->GetPhi1(),0);
 
  897       par.AddAt(pcon->GetDphi(),1);
 
  898       par.AddAt(pcon->GetNz(),2);
 
  899       for (Int_t i=0; i<nz; i++) {
 
  900          par.AddAt(z[i], 3+3*i);
 
  901          par.AddAt(rmin[i], 3+3*i+1);
 
  902          par.AddAt(rmax[i], 3+3*i+2);
 
  906    if (class_type==TGeoEltu::Class()) {
 
  910       TGeoEltu *eltu = (TGeoEltu*)shape;
 
  911       par.AddAt(eltu->GetA(),0);
 
  912       par.AddAt(eltu->GetB(),1);
 
  913       par.AddAt(eltu->GetDz(),2);
 
  916    if (class_type==TGeoHype::Class()) {
 
  920       TGeoHype *hype = (TGeoHype*)shape;
 
  921       par.AddAt(TMath::Sqrt(hype->RadiusHypeSq(0.,kTRUE)),0);
 
  922       par.AddAt(TMath::Sqrt(hype->RadiusHypeSq(0.,kFALSE)),1);
 
  923       par.AddAt(hype->GetDZ(),2);
 
  924       par.AddAt(hype->GetStIn(),3);
 
  925       par.AddAt(hype->GetStOut(),4);
 
  928    if (class_type==TGeoGtra::Class()) {
 
  932       TGeoGtra *trap = (TGeoGtra*)shape;
 
  933       Double_t tth = TMath::Tan(trap->GetTheta()*TMath::DegToRad());
 
  934       par.AddAt(trap->GetDz(),0);
 
  935       par.AddAt(tth*TMath::Cos(trap->GetPhi()*TMath::DegToRad()),1);
 
  936       par.AddAt(tth*TMath::Sin(trap->GetPhi()*TMath::DegToRad()),2);
 
  937       par.AddAt(trap->GetH1(),3);
 
  938       par.AddAt(trap->GetBl1(),4);
 
  939       par.AddAt(trap->GetTl1(),5);
 
  940       par.AddAt(TMath::Tan(trap->GetAlpha1()*TMath::DegToRad()),6);
 
  941       par.AddAt(trap->GetH2(),7);
 
  942       par.AddAt(trap->GetBl2(),8);
 
  943       par.AddAt(trap->GetTl2(),9);
 
  944       par.AddAt(TMath::Tan(trap->GetAlpha2()*TMath::DegToRad()),10);
 
  945       par.AddAt(trap->GetTwistAngle(),11);
 
  948    if (class_type==TGeoCtub::Class()) {
 
  952       TGeoCtub *ctub = (TGeoCtub*)shape;
 
  953       const Double_t *lx = ctub->GetNlow();
 
  954       const Double_t *tx = ctub->GetNhigh();
 
  955       par.AddAt(ctub->GetRmin(),0);
 
  956       par.AddAt(ctub->GetRmax(),1);
 
  957       par.AddAt(ctub->GetDz(),2);
 
  958       par.AddAt(ctub->GetPhi1(),3);
 
  959       par.AddAt(ctub->GetPhi2(),4);
 
  969    Error(
"GetShape",
"Getting shape parameters for shape %s not implemented", shape->ClassName());
 
  994 Bool_t TGeoMCGeometry::GetMaterial(
const TString &volumeName,
 
  995                             TString &name,Int_t &imat,
 
  996                             Double_t &a,Double_t &z,Double_t &dens,
 
  997                             Double_t &radl,Double_t &inter,TArrayD &par)
 
  999    TGeoVolume *vol = GetTGeoManager()->GetVolume(volumeName.Data());
 
 1000    if (!vol) 
return kFALSE;
 
 1001    TGeoMedium *med = vol->GetMedium();
 
 1002    if (!med) 
return kFALSE;
 
 1003    TGeoMaterial *mat = med->GetMaterial();
 
 1004    imat = mat->GetUniqueID();
 
 1005    name = mat->GetName();
 
 1006    name = name.Strip(TString::kTrailing, 
'$');
 
 1009    dens   = mat->GetDensity();
 
 1010    radl   = mat->GetRadLen();
 
 1011    inter  = mat->GetIntLen(); 
 
 1040 Bool_t TGeoMCGeometry::GetMedium(
const TString &volumeName,TString &name,
 
 1041                           Int_t &imed,Int_t &nmat,Int_t &isvol,Int_t &ifield,
 
 1042                           Double_t &fieldm,Double_t &tmaxfd,Double_t &stemax,
 
 1043                           Double_t &deemax,Double_t &epsil, Double_t &stmin,
 
 1046    TGeoVolume *vol = GetTGeoManager()->GetVolume(volumeName.Data());
 
 1047    if (!vol) 
return kFALSE;
 
 1048    TGeoMedium *med = vol->GetMedium();
 
 1049    if (!med) 
return kFALSE;
 
 1050    TGeoMaterial *mat = med->GetMaterial();
 
 1051    nmat = mat->GetUniqueID();
 
 1052    imed = med->GetId();
 
 1053    name = med->GetName();
 
 1054    name = name.Strip(TString::kTrailing, 
'$');
 
 1056    isvol  = (Int_t)med->GetParam(0);
 
 1057    ifield = (Int_t)med->GetParam(1);
 
 1058    fieldm = med->GetParam(2);
 
 1059    tmaxfd = med->GetParam(3);
 
 1060    stemax = med->GetParam(4);
 
 1061    deemax = med->GetParam(5);
 
 1062    epsil  = med->GetParam(6);
 
 1063    stmin  = med->GetParam(7);