00001
00002 #include "StiVMCToolKit.h"
00003 #ifdef __ROOT__
00004 #include "StMaker.h"
00005 #endif
00006 #if !defined(__CINT__) || defined(__MAKECINT__)
00007 #include <string.h>
00008 #include <assert.h>
00009 #include "Riostream.h"
00010 #include <stdio.h>
00011 #include "TSystem.h"
00012 #include "TMath.h"
00013 #include "TFile.h"
00014 #else
00015 #define BIT(n) (1 << (n))
00016 #endif
00017 static Int_t __addelement__(Int_t NElem,const TGeoMaterial *mat, Elem_t *ElementList) {
00018 return StiVMCToolKit::Add2ElementList(NElem, mat, ElementList);
00019 }
00020
00021 struct Elem_t {
00022 Int_t index;
00023 Double_t W;
00024 Double_t A;
00025 Double_t Z;
00026 Elem_t() : index(0),W(0),A(0),Z(0) {}
00027 };
00028 struct ElemV_t {
00029 const Char_t *name;
00030 Double_t A;
00031 Double_t Z;
00032 Double_t V;
00033 Double_t W;
00034 };
00035 static Int_t m_Debug = 0;
00036 static const Int_t NoElemMax = 10000;
00037 static VolumeMap_t VolumesToBeAveraged[] = {
00038 {"PIPE","the STAR beam pipe mother volume","HALL_1/CAVE_1/PIPE_1-2/*","",""},
00039 {"PIPC","the Central Beam PIPe Volum","HALL_1/CAVE_1/PIPE_1-2/PIPC_1/*","",""},
00040 {"PVAC","the Vacuum Volume of Be section of pipe","HALL_1/CAVE_1/PIPE_1-2/PIPC_1/PVAC_1","",""},
00041 {"PIPO","Steel pipe from Be to 1st flanges","HALL_1/CAVE_1/PIPE_1-2/PIPO_1/PVAO_1","",""},
00042 {"PVAO","its cavity","HALL_1/CAVE_1/PIPE_1-2/PIPO_1/PVAO_1","",""},
00043
00044 {"PXBX","Extrenal Berillium tube","HALL_1/CAVE_1/PXBX_1","",""},
00045
00046 {"SOUM", "Outer shileding structure","HALL_1/CAVE_1/SVTT_1/SOUM_1/*","",""},
00047 {"SXRL", "Circular water feeds","HALL_1/CAVE_1/SVTT_1/SXRL_1-2/*","",""},
00048 {"SXR1", "Circular water feeds","HALL_1/CAVE_1/SVTT_1/SXR1_3-4/*","",""},
00049 {"SXR2", "Circular water feeds","HALL_1/CAVE_1/SVTT_1/SXR2_5-6/*","",""},
00050 {"SCBM", "Mother of All Cables","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/*","",""},
00051 {"SCBL", "The bundles of cables connecting PCBs with the transition boards","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SCBL_1","",""},
00052 {"SCB1", "The bundles of cables connecting PCBs with the transition boards","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SCB1_2","",""},
00053 {"SCB2", "The bundles of cables connecting PCBs with the transition boards","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SCB2_3","",""},
00054 {"SCB3", "The bundles of cables connecting PCBs with the transition boards","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SCB3_4","",""},
00055 {"SFED", "bundles of water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SFED_1","",""},
00056 {"SFE1", "bundles of water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SFE1_2","",""},
00057 {"SFE2", "bundles of water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SFE2_3","",""},
00058 {"SPLS", "plastic of the water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SPLS_1","",""},
00059 {"SPL1", "plastic of the water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SPL1_2","",""},
00060 {"SPL2", "plastic of the water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SPL2_3","",""},
00061 {"SALM", "aluminum shield mesh","HALL_1/CAVE_1/SVTT_1/SALM_1-2","",""},
00062 {"SOSH", "SVT outer shield","HALL_1/CAVE_1/SVTT_1/SOSH_1","",""},
00063 {"SISH", "SVT inner shield","HALL_1/CAVE_1/SVTT_1/SISH_1","",""},
00064 {"SLYD", "layer mother","HALL_1/CAVE_1/SVTT_1/SLYD_1/*","",""},
00065 {"SLY1", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY1_2/*","",""},
00066 {"SLY2", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY2_3/*","",""},
00067 {"SLY3", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY3_4/*","",""},
00068 {"SLY4", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY4_5/*","",""},
00069 {"SLY5", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY5_6/*","",""},
00070 {"SVTD", "an active wafer volume","HALL_1/CAVE_1/SVTT_1/SLY*/SLS*/SLD*/STL*/STS*/SVTD_1","svt","SVTD"},
00071 {"SROD", "Support rod","HALL_1/CAVE_1/SVTT_1/SROD_1-2","",""},
00072 {"SBSP", "Beampipe support mother","HALL_1/CAVE_1/SVTT_1/SBSP_1-2","",""},
00073
00074 {"SBWC", "water manifold to support cone bracket mother","HALL_1/CAVE_1/SVTT_1/SBWC_1-2/*","",""},
00075 {"SWMM", "water manifold mother","HALL_1/CAVE_1/SVTT_1/SWMM_1-2/*","",""},
00076 {"SIES", "Volume to hold inner endring screws","HALL_1/CAVE_1/SVTT_1/SIES_1-2/*","",""},
00077 {"SOES", "Volume to hold outer endring screws","HALL_1/CAVE_1/SVTT_1/SOES_1-2/*","",""},
00078 {"SBRG", "Bracket joining the end rungs","HALL_1/CAVE_1/SVTT_1/SBRG_1-2/*","",""},
00079 {"SOER", "outer end ring","HALL_1/CAVE_1/SVTT_1/SOER_1-2/*","",""},
00080 {"SIRT", "inner end ring tube piece ","HALL_1/CAVE_1/SVTT_1/SIRT_1-2","",""},
00081 {"SIRP", "inner end ring polygon piece ","HALL_1/CAVE_1/SVTT_1/SIRP_1-2","",""},
00082 {"STAC", "twinax cable approximation, copper","HALL_1/CAVE_1/SVTT_1/SCON_1/STAC_1-2","",""},
00083
00084
00085 {"SCMP","SSD mounting plate inserted in the cone","HALL_1/CAVE_1/SVTT_1/SFMO_1/SCMP_1-8","",""},
00086 {"SCVM","SSD V-shape mouting piece","HALL_1/CAVE_1/SVTT_1/SFMO_1/SCVM_1-8/*","",""},
00087 {"SSLT","the linking (sector to the cone) tube","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSLT_1-8","",""},
00088 {"SSLB","the linking (sector to the cone)","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSLB_1-8","",""},
00089 {"SSRS","the side of the small rib","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSRS_1-4","",""},
00090 {"SSRT","the top of the side rib","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSRT_1-4","",""},
00091 {"SSSS","Side parts of the small sectors","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSSS_1-4","",""},
00092 {"SSST","Top parts of the small sectors","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSST_1-4","",""},
00093
00094 {"SFSM","the structure mother volume","HALL_1/CAVE_1/SVTT_1/SFMO_1/SFLM_1-20/SFSM_1/*","",""},
00095 {"SFDM","the detectors and adcs mother volume","HALL_1/CAVE_1/SVTT_1/SFMO_1/SFLM_1-20/SFDM_1/*","",""},
00096 {"SFSD","the strip detector", "HALL_1/CAVE_1/SVTT_1/SFMO_1/SFLM_1-20/SFDM_1/SFSW_1-16/SFSD_1","ssd",""},
00097
00098
00099 {"TPCW","the TPC supporting endcap Wheel","HALL_1/CAVE_1/TPCE_1/TPCW_1-2/*","",""},
00100 {"TPEA","one endcap placed in TPC","HALL_1/CAVE_1/TPCE_1/TPEA_1-2/*","",""},
00101 {"TPCM","the Central Membrane placed in TPC","HALL_1/CAVE_1/TPCE_1/TPCM_1","",""},
00102 {"TOFC","defines outer field cage - fill it with insulating gas already","HALL_1/CAVE_1/TPCE_1/TOFC_1/*","",""},
00103 {"TIFC","defines the Inner Field Cage placed in TPC","HALL_1/CAVE_1/TPCE_1/TIFC_1/*","",""},
00104 {"TPGV","the Gas Volume placed in TPC","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/*","",""},
00105 {"TPSS","a division of gas volume corresponding to a supersectors","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/TPSS_1-12/*","",""},
00106 {"TPAD","(inner) real padrow with dimensions defined at positioning time","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/TPSS_1-12/TPAD_1-39","tpc",""},
00107 {"TPA1","(outer) real padrow with dimensions defined at positioning time","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/TPSS_1-12/TPA1_40-73","tpc",""}
00108 };
00109 static Int_t NofVolToBEAveraged = sizeof(VolumesToBeAveraged)/sizeof(VolumeMap_t);
00110 static VolumeMap_t TopVolumes[] = {
00111 {"BBCM","Beam Beam Counter Modules Geometry","","",""},
00112 {"BTOF","the whole CTF system envelope","","",""},
00113 {"CALB","the geometry of the Barrel EM Calorimeter","","",""},
00114 {"ECAL","the EM EndCap Calorimeter GEOmetry","","",""},
00115 {"FBOX","one Pb-Glass fpd detector","","",""},
00116 {"FBO1","an other Pb-Glass fpd detector","","",""},
00117 {"FTMO","the mother of the single FTPC RO barrel","","",""},
00118 {"IBEM","the IBeam structure beneath the Bell reducer cone","","",""},
00119 {"MAGP","the geometry of the STAR magnet","","",""},
00120 {"PHMD","the Photon Multiplicity Detector","","",""},
00121 {"SUPO","the geometry of the Forward TPC supports in STAR","","",""},
00122
00123 {"SVTT","the SVT geometry for STAR","","",""},
00124
00125 {"SFMO","the mother of all Silicon Strip Detector volumes (inside of SVTT)","","",""},
00126
00127 {"FTPC","the geometry of the Forward TPC in STAR (inside of SVTT)","","",""},
00128
00129 {"UPST","the geometry of the UPSTREAM AreA.","","",""},
00130 {"VPDD","the Pseudo Vertex Position Detector of STAR","","",""},
00131 {"ZCAL","the geometry of the Zero deg. Quartz Calorimeter","","",""}
00132 };
00133 static Int_t nTopVol = sizeof(TopVolumes)/sizeof(VolumeMap_t);
00134
00135
00136
00137
00138
00139
00140
00141 static ElemV_t Air[] = {
00142 {"N_2", 14.00674, 7, 2*78.084, 0},
00143 {"O_2", 15.9994 , 8, 2*(20.946 + 0.033), 0},
00144 {"C", 12.011, 6, 0.033, 0},
00145 {"Ar", 39.948, 18, 0.934, 0}
00146 };
00147 static Int_t NAir = 0;
00148 static ElemV_t P10[] = {
00149 {"Ar", 40,18,0.95744680,0},
00150 {"C" , 12, 6,0.03191489,0},
00151 {"H" , 1, 1,0.01063830,0}
00152 };
00153 static Int_t NP10 = 3;
00154
00155
00156
00157
00158
00159
00160
00161
00162 void StiVMCToolKit::SetDebug(Int_t m) {m_Debug = m;}
00163
00164 Int_t StiVMCToolKit::Debug() {return m_Debug;}
00165
00166 void StiVMCToolKit::PrintShape(TGeoShape *shape) {
00167 TGeoBBox *box = 0, *Box = 0;
00168 TGeoTrd1 *trd1 = 0;
00169 TGeoTrd2 *trd2 = 0;
00170 TGeoTube *tube = 0;
00171 TGeoTubeSeg *tubs = 0;
00172 TGeoPcon *pcon = 0;
00173 TGeoPgon *pgon = 0;
00174 TGeoCone *cone = 0;
00175 TGeoConeSeg *cons = 0;
00176 TGeoArb8 *arb8 = 0;
00177 TGeoEltu *eltu = 0;
00178 Double_t *XY;
00179 Double_t dZ;
00180
00181 Double_t paramsBC[3];
00182 Int_t i, j;
00183 Int_t Nz = 0;
00184 shape->GetBoundingCylinder(paramsBC);
00185 shape->ComputeBBox();
00186 Box = (TGeoBBox *) shape;
00187 for (Int_t bit = 24; bit >= 9; bit--) {
00188 if (shape->TestShapeBit(BIT(bit))) {
00189 switch (BIT(bit)) {
00190 case TGeoShape::kGeoBox:
00191 box = (TGeoBBox *) shape;
00192 cout << "Box \tdX\t" << box->GetDX() << "\tdY\t" << box->GetDY() << "\tdZ\t" << box->GetDZ()
00193 << endl;
00194 break;
00195 case TGeoShape::kGeoTrd1:
00196 trd1 = (TGeoTrd1 *) shape;
00197 cout << "Trd1\tdX1\t" << trd1->GetDx1() << "\tdX2\t" << trd1->GetDx2()
00198 << "\tdY\t" << trd1->GetDy() << "\tdZ\t" << trd1->GetDz()
00199 << endl;
00200 break;
00201 case TGeoShape::kGeoTrd2:
00202 trd2 = (TGeoTrd2 *) shape;
00203 cout << "Trd2\tdX1\t" << trd2->GetDx1() << "\tdX2\t" << trd2->GetDx2()
00204 << "\tdY1\t" << trd2->GetDy1() << "\tdY2\t" << trd2->GetDy2()
00205 << "\tdZ\t" << trd2->GetDz() << endl;
00206 break;
00207 case TGeoShape::kGeoTubeSeg:
00208 tubs = (TGeoTubeSeg *) shape;
00209 cout << "Tubs\tRmin\t" << tubs->GetRmin() << "\tRmax\t" << tubs->GetRmax() << "\tdZ\t" << tubs->GetDz()
00210 << "\tPhi1\t" << tubs->GetPhi1() << "\tPhi2\t" << tubs->GetPhi2()
00211 << endl;
00212 break;
00213 case TGeoShape::kGeoTube:
00214 tube = (TGeoTube *) shape;
00215 cout << "Tube\tRmin\t" << tube->GetRmin() << "\tRmax\t" << tube->GetRmax() << "\tdZ\t" << tube->GetDz()
00216 << endl;
00217 break;
00218 case TGeoShape::kGeoPcon:
00219 pcon = (TGeoPcon *) shape;
00220 Nz = pcon->GetNz();
00221 cout << "Pcon"
00222 << "\tPhi1\t" << pcon->GetPhi1() << "\tDphi\t" << pcon->GetDphi() << "\tNz\t" << Nz << endl;
00223 for (i = 0; i < Nz; i++) {
00224 cout << i << "\tZ\t" << pcon->GetZ(i) << "\tRmin\t" << pcon->GetRmin(i) << "\tRmax\t" << pcon->GetRmax(i) << endl;
00225 }
00226 cout << endl;
00227 break;
00228 case TGeoShape::kGeoPgon:
00229 pgon = (TGeoPgon *) shape;
00230 Nz = pgon->GetNz();
00231
00232 cout << "Pgon\tNedges\t" << pgon->GetNedges()
00233 << "\tPhi1\t" << pgon->GetPhi1() << "\tDphi\t" << pgon->GetDphi() << "\tNz\t" <<Nz << endl;
00234 for (i = 0; i <Nz; i++) {
00235 cout << i << "\tZ\t" << pgon->GetZ(i) << "\tRmin\t" << pgon->GetRmin(i) << "\tRmax\t" << pgon->GetRmax(i) << endl;
00236 }
00237 cout << endl;
00238 break;
00239 case TGeoShape::kGeoCone:
00240 cone = (TGeoCone *) shape;
00241 cout << "Cone\tdZ\t" << cone->GetDz()
00242 << "\tRmin1\t" << cone->GetRmin1() << "\tRmax1\t" << cone->GetRmax1()
00243 << "\tRmin2\t" << cone->GetRmin2() << "\tRmax2\t" << cone->GetRmax2()
00244 << endl;
00245 break;
00246 case TGeoShape::kGeoConeSeg:
00247 cons = (TGeoConeSeg *) shape;
00248 cout << "Cons\tdZ\t" << cons->GetDz()
00249 << "\tPhi1\t" << cons->GetPhi1() << "\tPhi2\t" << cons->GetPhi2()
00250 << "\tRmin1\t" << cons->GetRmin1() << "\tRmax1\t" << cons->GetRmax1()
00251 << "\tRmin2\t" << cons->GetRmin2() << "\tRmax2\t" << cons->GetRmax2()
00252 << endl;
00253 break;
00254 case TGeoShape::kGeoArb8:
00255 case TGeoShape::kGeoTrap:
00256 arb8 = (TGeoArb8 *) shape;
00257 XY = arb8->GetVertices();
00258 dZ = arb8->GetDz();
00259 for (j = 0; j < 2; j++) {
00260 cout << "Trap/Arb8\tdZ\t";
00261 if (j == 0) cout << -dZ;
00262 else cout << dZ;
00263 for (i = 4*j; i < 4*j+4; i++)
00264 cout << "\t(" << XY[2*i] << "," << XY[2*i+1] << ")";
00265 cout << endl;
00266 }
00267 break;
00268 case TGeoShape::kGeoEltu:
00269 eltu = (TGeoEltu *) shape;
00270 cout << "Eltu\tdZ\t" << eltu->GetDz()
00271 << "\tA\t" << eltu->GetA() << "\tB\t" << eltu->GetB()
00272 << endl;
00273 break;
00274 case TGeoShape::kGeoTorus:
00275 case TGeoShape::kGeoPara:
00276 case TGeoShape::kGeoSph:
00277 case TGeoShape::kGeoCtub:
00278 default:
00279 cout << bit << "\t has not yet implemented for " << shape->GetName() << endl;
00280 break;
00281 }
00282 break;
00283 }
00284 }
00285 }
00286
00287 Double_t StiVMCToolKit::GetShapeVolume(TGeoShape *shape) {
00288 TGeoBBox *box = 0, *Box = 0;
00289 TGeoTrd1 *trd1 = 0;
00290 TGeoTrd2 *trd2 = 0;
00291 TGeoTube *tube = 0;
00292 TGeoTubeSeg *tubs = 0;
00293 TGeoPcon *pcon = 0;
00294 TGeoPgon *pgon = 0;
00295 TGeoCone *cone = 0;
00296 TGeoConeSeg *cons = 0;
00297 TGeoEltu *eltu = 0;
00298 Double_t volume = 0;
00299 Double_t paramsBC[3];
00300 Double_t volBB = 0;
00301 Double_t volBC = 0;
00302 Int_t i;
00303 Int_t Nz = 0;
00304 shape->GetBoundingCylinder(paramsBC);
00305 shape->ComputeBBox();
00306 Box = (TGeoBBox *) shape;
00307 volBB = 8*Box->GetDX()*Box->GetDY()*Box->GetDZ();
00308 volBC = 2*TMath::Pi()*(paramsBC[1] - paramsBC[0])*Box->GetDZ();;
00309 for (Int_t bit = 24; bit >= 9; bit--) {
00310 if (shape->TestShapeBit(BIT(bit))) {
00311 Int_t kbit = BIT(bit);
00312 switch (kbit) {
00313 case TGeoShape::kGeoBox:
00314 box = (TGeoBBox *) shape;
00315 volume = 8*box->GetDX()*box->GetDY()*box->GetDZ();
00316 break;
00317 case TGeoShape::kGeoTrd1:
00318 trd1 = (TGeoTrd1 *) shape;
00319 volume = 4*(trd1->GetDx1()+trd1->GetDx2())*trd1->GetDy()*trd1->GetDz();
00320 break;
00321 case TGeoShape::kGeoTrd2:
00322 trd2 = (TGeoTrd2 *) shape;
00323 volume = 2*(trd2->GetDx1() + trd2->GetDx2())*(trd2->GetDy1() + trd2->GetDy2())*trd2->GetDz();
00324 break;
00325 case TGeoShape::kGeoTubeSeg:
00326 case TGeoShape::kGeoTube:
00327 tube = (TGeoTube *) shape;
00328 volume = 2*TMath::Pi()*(tube->GetRmax()* tube->GetRmax() - tube->GetRmin()*tube->GetRmin())*
00329 tube->GetDz();
00330 if (kbit == TGeoShape::kGeoTubeSeg) {
00331 tubs = (TGeoTubeSeg *) shape;
00332 volume *= TMath::Abs(tubs->GetPhi2()-tubs->GetPhi1())/360.;
00333 }
00334 break;
00335 case TGeoShape::kGeoPcon:
00336 case TGeoShape::kGeoPgon:
00337 pcon = (TGeoPcon *) shape;
00338 Nz = pcon->GetNz();
00339 volume = 0;
00340 for (i = 1; i < Nz; i++) {
00341 volume +=
00342 ((pcon->GetRmax(i)*pcon->GetRmax(i) + pcon->GetRmax(i-1)*pcon->GetRmax(i-1) + pcon->GetRmax(i)*pcon->GetRmax(i-1)) -
00343 (pcon->GetRmin(i)*pcon->GetRmin(i) + pcon->GetRmin(i-1)*pcon->GetRmin(i-1) + pcon->GetRmin(i)*pcon->GetRmin(i-1)))*
00344 TMath::Abs((pcon->GetZ(i) - pcon->GetZ(i-1)));
00345 }
00346 if (kbit == TGeoShape::kGeoPgon) {
00347 pgon = (TGeoPgon *) shape;
00348 volume *= TMath::Tan(TMath::Pi()/180.*TMath::Abs(pgon->GetDphi()/pgon->GetNedges()/2.))/3.*pgon->GetNedges();
00349 }
00350 else {
00351 volume *= TMath::Pi()/3.*TMath::Abs(pcon->GetDphi())/360.;
00352 }
00353 break;
00354 case TGeoShape::kGeoCone:
00355 case TGeoShape::kGeoConeSeg:
00356 cone = (TGeoCone *) shape;
00357 volume = TMath::Pi()*2./3.*
00358 ((cone->GetRmax2()*cone->GetRmax2() + cone->GetRmax1()*cone->GetRmax1() + cone->GetRmax2()*cone->GetRmax1()) -
00359 (cone->GetRmin2()*cone->GetRmin2() + cone->GetRmin1()*cone->GetRmin1() + cone->GetRmin2()*cone->GetRmin1()))*
00360 cone->GetDz();
00361 if (kbit == TGeoShape::kGeoConeSeg) {
00362 cons = (TGeoConeSeg *) shape;
00363 volume *= TMath::Abs(cons->GetPhi2() - cons->GetPhi1())/360.;
00364 }
00365 break;
00366 case TGeoShape::kGeoEltu:
00367 eltu = (TGeoEltu *) shape;
00368 volume = 2*TMath::Pi()*eltu->GetA()*eltu->GetB()*eltu->GetDz();
00369 break;
00370 case TGeoShape::kGeoArb8:
00371 case TGeoShape::kGeoTrap:
00372 case TGeoShape::kGeoTorus:
00373 case TGeoShape::kGeoPara:
00374 case TGeoShape::kGeoSph:
00375 case TGeoShape::kGeoCtub:
00376 default:
00377 if (Debug())
00378 cout << bit << "\t has not yet implemented for " << shape->GetName() << endl;
00379 volume = TMath::Min(volBB,volBC);
00380 break;
00381 }
00382 break;
00383 }
00384 }
00385 if (Debug())
00386 cout << " =============> Volume = " << volume
00387 << "\t" << shape->GetName() << "\tBB\t" << volBB << "\tBC\t" << volBC << endl;
00388 assert(volume >= 0);
00389 if (Debug() && ! (volBB - volume > -1e-7 && volBC - volume > -1.e-7) ) {
00390 PrintShape(shape);
00391 cout << "volume\t" << volume << "\tvolBB\t" << volBB << "\tvolBC\t" << volBC << endl;
00392 cout << "\tBoundingBox dX\t" << Box->GetDX() << "\tdY\t" << Box->GetDY() << "\tdZ\t" << Box->GetDZ() << endl;
00393 cout << "\tBoundingCylinder\trMin\t" << TMath::Sqrt(paramsBC[0]) << "\trMax\t" << TMath::Sqrt(paramsBC[1])
00394 << "\tdZ\t" << Box->GetDZ() << endl;
00395
00396 assert(volBB - volume > -1e-7 && volBC - volume > -1.e-7);
00397 }
00398 return volume;
00399 }
00400
00401 Int_t StiVMCToolKit::Add2ElementList(Int_t NElem,const TGeoMaterial *mat, Elem_t *ElementList) {
00402 assert(NElem>=0 && NElem<NoElemMax);
00403 if (! NAir) {
00404 NAir = sizeof(Air)/sizeof(ElemV_t);
00405 Double_t W = 0;
00406 for (Int_t i = 0; i < NAir; i++) W += Air[i].A*Air[i].V;
00407 for (Int_t i = 0; i < NAir; i++) Air[i].W = Air[i].A*Air[i].V/W;
00408 }
00409 if (mat->InheritsFrom("TGeoMixture")) {
00410 TGeoMixture *mix = (TGeoMixture *) mat;
00411 Int_t N = mix->GetNelements();
00412 Double_t *A = mix->GetAmixt();
00413 Double_t *Z = mix->GetZmixt();
00414 Double_t *W = mix->GetWmixt();
00415 for (Int_t i = 0; i < N; i++) {
00416 if (TMath::Abs(Z[i] - 7.30) < 1.e-3 &&
00417 TMath::Abs(A[i] - 14.61) < 1.e-3) {
00418 for (Int_t j = 0; j < NAir; j++) {
00419 ElementList[NElem].index = NElem+1;
00420 ElementList[NElem].W = W[i]*Air[j].W;
00421 ElementList[NElem].A = Air[j].A;
00422 ElementList[NElem].Z = Air[j].Z;
00423 NElem++;
00424 }
00425 } else {
00426 if (TMath::Abs(Z[i] - 17.436 ) < 1.e-3 &&
00427 TMath::Abs(A[i] - 38.691) < 1.e-3) {
00428 for (Int_t j = 0; j < NP10; j++) {
00429 ElementList[NElem].index = NElem+1;
00430 ElementList[NElem].W = W[i]*P10[j].W;
00431 ElementList[NElem].A = P10[j].A;
00432 ElementList[NElem].Z = P10[j].Z;
00433 NElem++;
00434 }
00435 } else {
00436 ElementList[NElem].index = NElem+1;
00437 ElementList[NElem].W = W[i];
00438 ElementList[NElem].A = A[i];
00439 ElementList[NElem].Z = Z[i];
00440 NElem++;
00441 }
00442 }
00443 }
00444 } else {
00445 Double_t A = mat->GetA();
00446 Double_t Z = mat->GetZ();
00447 if (TMath::Abs(Z - 7.30) < 1.e-3 &&
00448 TMath::Abs(A - 14.61) < 1.e-3) {
00449 for (Int_t j = 0; j < NAir; j++) {
00450 ElementList[NElem].index = NElem+1;
00451 ElementList[NElem].W = Air[j].W;
00452 ElementList[NElem].A = Air[j].A;
00453 ElementList[NElem].Z = Air[j].Z;
00454 NElem++;
00455 }
00456 } else {
00457 if (TMath::Abs(Z - 17.436 ) < 1.e-3 &&
00458 TMath::Abs(A - 38.691) < 1.e-3) {
00459 for (Int_t j = 0; j < NP10; j++) {
00460 ElementList[NElem].index = NElem+1;
00461 ElementList[NElem].W = P10[j].W;
00462 ElementList[NElem].A = P10[j].A;
00463 ElementList[NElem].Z = P10[j].Z;
00464 NElem++;
00465 }
00466 } else {
00467 ElementList[NElem].index = NElem+1;
00468 ElementList[NElem].W = 1.;
00469 ElementList[NElem].A = A;
00470 ElementList[NElem].Z = Z;
00471 NElem++;
00472 }
00473 }
00474 }
00475 assert(NElem < NoElemMax);
00476 return NElem;
00477 }
00478
00479 Int_t StiVMCToolKit::Merge2ElementList(Int_t NElem, Elem_t *ElementList,
00480 Int_t NElemD, Elem_t *ElementListD, Double_t weight) {
00481 assert(NElem>=0 && NElem<NoElemMax);
00482 for (Int_t i = 0; i < NElemD; i++) {
00483 ElementList[NElem].index = NElem+1;
00484 ElementList[NElem].W = weight*ElementListD[i].W;
00485 ElementList[NElem].A = ElementListD[i].A;
00486 ElementList[NElem].Z = ElementListD[i].Z;
00487 NElem++;
00488 assert(NElem < NoElemMax);
00489 }
00490 return NElem;
00491 }
00492
00493 Int_t StiVMCToolKit::NormolizeElementList(Int_t NElem, Elem_t *ElementList){
00494 assert(NElem>=0 && NElem<NoElemMax);
00495 Double_t W = 0;
00496 for (Int_t i = 0; i < NElem; i++) W += ElementList[i].W ;
00497 for (Int_t i = 0; i < NElem; i++) ElementList[i].W /= W;
00498 for (Int_t k = 0; k < NElem; k++) {
00499 if (ElementList[k].W > 0) {
00500 for (Int_t l = k+1; l < NElem; l++) {
00501 if (ElementList[l].W > 0 &&
00502 ElementList[k].A == ElementList[l].A &&
00503 ElementList[k].Z == ElementList[l].Z) {
00504 ElementList[k].W += ElementList[l].W;
00505 ElementList[l].W = -1;
00506 }
00507 }
00508 }
00509 }
00510 Int_t N = 0;
00511
00512 for (Int_t k = 0; k < NElem; k++) {
00513 if (ElementList[k].W <= 0) continue;
00514 if (k != N)
00515 ElementList[N] = ElementList[k];
00516 N++;
00517 }
00518 return N;
00519 }
00520
00521 Double_t StiVMCToolKit::GetWeight(TGeoNode *nodeT, const TString &pathT,
00522 Int_t *NElem, Elem_t *ElementList) {
00523 Double_t Weight = 0;
00524 Double_t WeightT = 0;
00525 if (! nodeT) {
00526 if (! gGeoManager) return Weight;
00527 gGeoManager->RestoreMasterVolume();
00528 gGeoManager->CdTop();
00529 gGeoManager->cd(pathT.Data());
00530 nodeT = gGeoManager->GetCurrentNode();
00531 if (! nodeT) return Weight;
00532 }
00533 TGeoVolume *volT = nodeT->GetVolume();
00534 const TGeoMaterial *mat = volT->GetMaterial();
00535
00536 Double_t dens = mat->GetDensity();
00537 TGeoShape *shapeT = (TGeoShape *) volT->GetShape();
00538 if (! shapeT) return Weight;
00539 Double_t volume = GetShapeVolume(shapeT);
00540 WeightT = Weight = dens*volume;
00541 #if 0
00542 if (! ElementList ) {
00543 ElementList = new Elem_t[NoElemMax];
00544 NElem = new Int_t[1]; NElem[0] = 0;
00545 }
00546 #else
00547 assert(ElementList && NElem);
00548 #endif
00549 Int_t im1 = NElem[0];
00550 NElem[0] = __addelement__(NElem[0], mat, ElementList);
00551 Int_t im2 = NElem[0];
00552
00553
00554 TObjArray *nodes = volT->GetNodes();
00555 Int_t nd = nodeT->GetNdaughters();
00556
00557 if (nd > 0) {
00558 vector<Double_t> weights(nd);
00559 for (Int_t id = 0; id < nd; id++) {
00560 TGeoNode *node = (TGeoNode*)nodes->UncheckedAt(id);
00561 if (! node) continue;
00562
00563 TGeoVolume *vol = node->GetVolume();
00564 if (! vol) continue;
00565 TString name(TString(vol->GetName()));
00566
00567 TGeoShape *shapeD = (TGeoShape *) vol->GetShape();
00568 if (! shapeD ) continue;
00569 Double_t weight = dens*GetShapeVolume(shapeD);
00570 WeightT -= weight;
00571 Weight -= weight;
00572 TString path = pathT;
00573 if (path != "") path += "/";
00574 path += node->GetName();
00575 Int_t NElemD[1] = {0};
00576 Elem_t ElementListD[NoElemMax];
00577 weights[id] = GetWeight(node,path,NElemD,ElementListD);
00578 Weight += weights[id];
00579 NElem[0] = Merge2ElementList(NElem[0], ElementList, NElemD[0], ElementListD, weights[id]);
00580 }
00581 }
00582 for (Int_t i = im1; i < im2; i++) {
00583 ElementList[i].W *= WeightT;
00584 }
00585 NElem[0] = NormolizeElementList(NElem[0],ElementList);
00586 return Weight;
00587 }
00588
00589
00590 Double_t StiVMCToolKit::GetVolumeWeight(TGeoVolume *volT, Int_t *NElem, Elem_t *ElementList)
00591 {
00592 if (! volT || !volT->GetShape()) return 0;
00593
00594 Double_t Weight = 0;
00595 Double_t WeightT = 0;
00596 const TGeoMaterial *mat = volT->GetMaterial();
00597 TGeoShape *shapeT = (TGeoShape *) volT->GetShape();
00598 if (! shapeT) return 0;
00599
00600 Double_t dens = mat->GetDensity();
00601 if (Debug())
00602 cout << "volT \t" << volT->GetName() << "\t" << volT->GetTitle() << endl;
00603 Double_t volume = GetShapeVolume(shapeT);
00604 Weight = dens*volume;
00605 WeightT = Weight;
00606 #if 0
00607 if (! ElementList ) {
00608
00609 ElementList = new Elem_t[NoElemMax];
00610 NElem = new Int_t[1]; NElem[0] = 0;
00611 }
00612 #else
00613 assert(ElementList && "StiVMCToolKit::GetWeight: Create the the array first!");
00614 #endif
00615 Int_t im1 = *NElem;
00616 Int_t im2 = __addelement__(im1, mat, ElementList);
00617 *NElem = im2;
00618
00619
00620 TObjArray *nodes = volT->GetNodes();
00621 Int_t nd = volT->GetNdaughters();
00622
00623 if (nd > 0) {
00624 vector<Double_t> weights(nd);
00625 for (Int_t id = 0; id < nd; id++) {
00626 TGeoNode *node = (TGeoNode*) nodes->UncheckedAt(id);
00627 if (! node) continue;
00628
00629 TGeoVolume *vol = node->GetVolume();
00630 if (! vol) continue;
00631 TString name(TString(vol->GetName()));
00632
00633 TGeoShape *shapeD = (TGeoShape *) vol->GetShape();
00634 if (! shapeD ) continue;
00635 Double_t weight = dens*GetShapeVolume(shapeD);
00636 WeightT -= weight;
00637 Weight -= weight;
00638 Int_t NElemD[1] = {0};
00639 Elem_t ElementListD[NoElemMax];
00640 weights[id] = GetVolumeWeight(vol,NElemD,ElementListD);
00641 Weight += weights[id];
00642 Int_t im =*NElem;
00643 Int_t imm = Merge2ElementList(im, ElementList, NElemD[0], ElementListD, weights[id]);
00644 *NElem = imm;
00645 }
00646 }
00647 for (Int_t i = im1; i < im2; i++) {
00648 ElementList[i].W *= WeightT;
00649 }
00650 Int_t in = *NElem;
00651 Int_t inn= NormolizeElementList(in,ElementList);
00652 *NElem = inn;
00653 return Weight;
00654 }
00655
00656 void StiVMCToolKit::MakeAverageVolume(TGeoVolume *volT, TGeoShape *&newshape, TGeoMedium *&newmed, Double_t *xyzM) {
00657 if (! volT) return;
00658 vector<Elem_t> ElementList(NoElemMax);
00659 Int_t NElem = 0;
00660
00661 Int_t *al = &NElem;
00662 Elem_t *el = &ElementList[0];
00663 TGeoVolume *v = volT;
00664
00665 Double_t Weight =
00666 GetVolumeWeight(v
00667 , al
00668 , el
00669 );
00670 const TGeoMedium *med = volT->GetMedium();
00671 const TGeoMaterial *mat = volT->GetMaterial();
00672 if (Debug()) {
00673 if (med->GetParam(0)) cout << "===================" << endl;
00674 else cout << "+++++++++++++++++++" << endl;
00675 cout << volT->GetName() << " === "
00676 << " Weight " << Weight << "[g]\t";
00677 cout << "material\t" << mat->GetName() << "\tmedium\t" << med->GetName() << endl;
00678 }
00679 TGeoShape *shapeT = (TGeoShape *) volT->GetShape();
00680 Double_t volume = GetShapeVolume(shapeT);
00681 if (Debug()) PrintShape(shapeT);
00682 shapeT->ComputeBBox();
00683 TGeoBBox * Box = (TGeoBBox *) shapeT;
00684 Double_t dx=-2001, dy=-2002, dz=-2003, rmin=-2009, rmax=-2010;
00685 Double_t paramsBC[4];
00686 shapeT->GetBoundingCylinder(paramsBC);
00687 Double_t volBB = 8*Box->GetDX()*Box->GetDY()*Box->GetDZ();
00688 Double_t volBC = 2*TMath::Pi()*(paramsBC[1] - paramsBC[0])*Box->GetDZ();
00689 Double_t volB = 0;
00690 if (Debug()) cout << shapeT->GetName();
00691 enum ShapeId {kBox = 1, kTube, kKeep};
00692 Int_t iShape = 3;
00693 TString name(volT->GetName());
00694 if (name == "SROD") {
00695 iShape = 1;
00696 Double_t S = TMath::Pi()*(paramsBC[1] - paramsBC[0]);
00697 dy = TMath::Sqrt(paramsBC[1]) - TMath::Sqrt(paramsBC[0]);
00698 dx = S/(4*dy);
00699 dz = Box->GetDZ();
00700 Box->SetBoxDimensions(dx,dy,dz);
00701 volB = volBB = volBC;
00702 } else {
00703 if ((shapeT->TestShapeBit(TGeoShape::kGeoTube) ||
00704 shapeT->TestShapeBit(TGeoShape::kGeoTubeSeg)) &&
00705 xyzM && xyzM[0]*xyzM[0] + xyzM[1]*xyzM[1] < 1.e-3) {
00706 iShape = 3;
00707 volB = volume;
00708 if (Debug()) cout << "\tKeep shape ===========\t";
00709 }
00710 else {
00711 if ( (volBB <= volBC) || (xyzM && xyzM[0]*xyzM[0] + xyzM[1]*xyzM[1] > 1.e-3) ) {
00712 iShape = 1;
00713 volB = volBB;
00714 dx = Box->GetDX();
00715 dy = Box->GetDY();
00716 dz = Box->GetDZ();
00717 if (Debug()) cout << "\tBoundingBox dX\t" << dx << "\tdY\t" << dy << "\tdZ\t" << dz;
00718 } else {
00719 iShape = 2;
00720 volB = volBC;
00721 rmin = TMath::Sqrt(paramsBC[0]);
00722 rmax = TMath::Sqrt(paramsBC[1]);
00723 dz = Box->GetDZ();
00724 if (Debug()) cout << "\tBoundingCylinder\trMin\t" << rmin << "\trMax\t" << rmax << "\tdZ\t" << dz;
00725 }
00726 }
00727 }
00728 Double_t Dens = Weight/volB;
00729 if (Debug()) cout << "\tvolume\t" << volume << "\tvolB\t" << volB << "\tDens\t" << Dens << endl;
00730 TString newmatname(mat->GetName());
00731 newmatname += "_";
00732 newmatname += name;
00733 TGeoMixture *newmat = new TGeoMixture(newmatname.Data(),NElem, Dens);
00734 for (Int_t i = 0; i < NElem; i++) {
00735 if (Debug())
00736 cout << Form("id%3i\tW%10.3g\tA%5.1f\t%5.1f\n",
00737 ElementList[i].index,ElementList[i].W,ElementList[i].A,ElementList[i].Z);
00738
00739
00740 newmat->DefineElement(i,ElementList[i].A,ElementList[i].Z,ElementList[i].W);
00741 }
00742 newmat->SetUniqueID(gGeoManager->AddMaterial(newmat));
00743 Int_t matid = newmat->GetUniqueID();
00744
00745 Double_t params[10];
00746 for (Int_t i = 0; i < 10; i++) params[i] = med->GetParam(i);
00747 newmed = new TGeoMedium(newmatname.Data(),matid,newmat,params);
00748 if (Debug() && xyzM) {
00749 cout << "\txyzM x\t" << xyzM[0] << "\ty\t" << xyzM[1] << "\tz\t" << xyzM[2] << endl;
00750 }
00751 newshape = shapeT;
00752 if (iShape == 1) newshape = new TGeoBBox(dx, dy, dz);
00753 else if (iShape == 2) {
00754 newshape = new TGeoTube(rmin, rmax, dz);
00755 }
00756 }
00757
00758 TGeoManager * StiVMCToolKit::GetVMC() {
00759 TGeoManager *gGeo = 0;
00760 #ifndef __NOVMC__
00761
00763 gGeo = gGeoManager;
00764 if (gGeo) return gGeo;
00765 cout << "StiVMCToolKit::GetVMC() -I- Get VMC geometry" <<endl;
00766 if (StMaker::GetChain()) {
00767 StMaker::GetChain()->GetDataBase("VmcGeometry");
00768 }
00769 if (! gGeoManager)
00770 cout << "StiVMCToolKit::GetVMC() -E- Can't get VMC geometry" <<endl;
00771 gGeo = gGeoManager;
00772 #endif
00773 return gGeo;
00774 }
00775
00776 void StiVMCToolKit::TestVMC4Reconstruction(){
00779 GetVMC();
00780
00781 Int_t NV[2] = {nTopVol,NofVolToBEAveraged};
00782 VolumeMap_t *list = 0;
00783 TGeoVolume *volT = 0;
00784 cout << "<table>" << endl;
00785 cout << "<tr><td>name</td><td>Comment </td><td>Volume[cm**3]</td><td>Weight[g]</td><td>Av.Dens.[g/cm**3]</td></tr>" << endl;
00786 for (Int_t i = 0; i < 2; i++) {
00787 if (i == 0) list = TopVolumes;
00788 else list = VolumesToBeAveraged;
00789 for (Int_t j = 0; j < NV[i]; j++) {
00790 volT = gGeoManager->GetVolume(list[j].name);
00791
00792 if (! volT) {cout << "Can't find " << list[j].name << "\t" << list[j].comment << endl; continue;}
00793 TGeoShape *shapeT = volT->GetShape();
00794 Double_t volume = GetShapeVolume(shapeT);
00795 Double_t weight = GetVolumeWeight(volT, 0, 0);
00796 cout << "<tr><td>"<< list[j].name << "</td><td>" << list[j].comment << "</td>"
00797 << "<td>" << Form("%10.3g",volume) << "</td>"
00798 << "<td>" << Form("%10.3g",weight) << "</td>"
00799 << "<td>" << Form("%10.3g",weight/volume) << "</td></tr>"
00800 << endl;
00801 }
00802 }
00803 cout << "</table>" << endl;
00804 }
00805
00806 TGeoPhysicalNode *StiVMCToolKit::Alignment(const TGeoNode *nodeT, const Char_t *pathT,
00807 TGeoVolume *volT, TGeoShape *newshape, TGeoMedium* newmed) {
00808
00809 TObjArray *listP = gGeoManager->GetListOfPhysicalNodes();
00810 TGeoPhysicalNode *nodeP = 0;
00811 TGeoCombiTrans *trP = 0;
00812 if (listP) {
00813 Int_t N = listP->GetEntries();
00814 for (Int_t i = 0; i < N; i++) {
00815 TGeoPhysicalNode *nod = dynamic_cast<TGeoPhysicalNode *> (listP->At(i));
00816 if (nod && TString(nod->GetName()) == TString(pathT)) {
00817 nodeP = nod; break;
00818 }
00819 }
00820 }
00821 Double_t local[3] = {0,0,0};
00822 TGeoShape *shapeT = volT->GetShape();
00823 if (Debug()) {
00824 cout << "StiVMCToolKit::Alignment -I node\t" << nodeT->GetName()
00825 << "\tvolume\t" << volT->GetName() << "\t:" << volT->GetTitle() << endl;
00826 cout << "\tshape\t" << shapeT->GetName() << "\tmaterial\t" << volT->GetMaterial()->GetName();
00827 if (newshape) cout << "\tnewshape\t" << newshape->GetName();
00828 cout << endl;
00829 }
00830 if (shapeT->TestShapeBit(TGeoShape::kGeoPcon) || shapeT->TestShapeBit(TGeoShape::kGeoPgon)) {
00831 TGeoPcon *pcon = (TGeoPcon *) shapeT;
00832 Int_t Nz = pcon->GetNz();
00833 local[2] = 0.5*(pcon->GetZ(0) + pcon->GetZ(Nz-1));
00834 }
00835 Double_t master[3];
00836 nodeT->LocalToMaster(local,master);
00837 if (Debug())
00838 cout << "\tmaster x\t" << master[0] << "\ty\t" << master[1] << "\tz\t" << master[2] << endl;
00839 if (!nodeP) nodeP = gGeoManager->MakePhysicalNode(pathT);
00840 if (!nodeP) return nodeP;
00841 #if 0
00842 if (nodeP->IsAligned()) {
00843 trP = nodeP->GetNode()->GetMatrix();
00844 trP->SetTranslation(master[0], master[1], master[2]);
00845 } else {
00846 trP = new TGeoCombiTrans(master[0], master[1], master[2], nodeP->GetNode()->GetMatrix());
00847 }
00848 #endif
00849 nodeP->Align(trP,newshape);
00850 TGeoVolume *newvol = nodeP->GetNode(-1)->GetVolume();
00851 if (Debug())
00852 cout << "newvol\t" << newvol->GetName() << "\tmed\t" << newvol->GetMedium()->GetName() << endl;
00853 newvol->SetMedium(newmed);
00854 newvol->SetTitle("AVERAGED");
00855 newvol->SetLineColor(1);
00856 newvol->SetVisibility(1);
00857 TObjArray *nodes = newvol->GetNodes();
00858 if (nodes && ! newvol->TObject::TestBit(TGeoVolume::kVolumeImportNodes)) {
00859 nodes->Delete();
00860 delete nodes;
00861 newvol->ClearNodes();
00862 }
00863 return nodeP;
00864 }
00865
00866 TGeoPhysicalNode *StiVMCToolKit::LoopOverNodes(const TGeoNode *nodeT, const Char_t *pathT,
00867 const Char_t *name,void ( *callback)(TGeoPhysicalNode *)){
00868 TGeoVolume *volT = nodeT->GetVolume();
00869 const Char_t *nameT = volT->GetName();
00870 TGeoPhysicalNode *nodeP = 0;
00871 TGeoShape *newshape = 0;
00872 TGeoMedium* newmed = 0;
00873 if (nameT && name && ! strncmp(nameT,name,4)) {
00874
00875 Double_t local[3] = {0, 0, 0};
00876 TGeoShape *shapeT = volT->GetShape();
00877 if (shapeT->TestShapeBit(TGeoShape::kGeoPcon) || shapeT->TestShapeBit(TGeoShape::kGeoPgon)) {
00878 TGeoPcon *pcon = (TGeoPcon *) shapeT;
00879 Int_t Nz = pcon->GetNz();
00880 local[2] = 0.5*(pcon->GetZ(0) + pcon->GetZ(Nz-1));
00881 }
00882 Double_t master[3];
00883 nodeT->LocalToMaster(local,master);
00884 MakeAverageVolume(volT, newshape, newmed, master);
00885 nodeP = Alignment(nodeT,pathT, volT, newshape, newmed);
00886 if (nodeP) {
00887 if (! callback) MakeVolume(nodeP);
00888 else callback(nodeP);
00889 }
00890 return nodeP;
00891 }
00892 TObjArray *nodes = volT->GetNodes();
00893 Int_t nd = volT->GetNdaughters();
00894 for (Int_t id = 0; id < nd; id++) {
00895 TGeoNode *node = (TGeoNode*) nodes->UncheckedAt(id);
00896 if (! node) continue;
00897 TString path = pathT;
00898 if (path != "") path += "/";
00899 path += node->GetName();
00900 if (! name && Debug()) {
00901 cout << path;
00902 TGeoVolume *vol = node->GetVolume();
00903 const TGeoMedium *med = vol->GetMedium();
00904 if (med->GetParam(0)) {cout << "\t===================" << endl; continue;}
00905 else cout << endl;
00906 }
00907 LoopOverNodes(node, path, name, callback);
00908 }
00909 return nodeP;
00910 }
00911
00912 void StiVMCToolKit::MakeVolume(TGeoPhysicalNode *nodeP) {
00913 if (Debug()) {
00914 cout << "StiVMCToolKit::MakeVolume -I TGeoPhysicalNode\t" << nodeP->GetName() << endl;
00915 TGeoVolume *volP = nodeP->GetVolume();
00916 TGeoMaterial *matP = volP->GetMaterial(); matP->Print("");
00917 TGeoShape *shapeP = nodeP->GetShape(); cout << "New Shape\t"; PrintShape(shapeP);
00918 TGeoHMatrix *hmat = nodeP->GetMatrix(); hmat->Print("");
00919 }
00920 }
00921
00922 Double_t StiVMCToolKit::GetPotI(const TGeoMaterial *mat) {
00923 Double_t PotI = 0;
00924 if (mat) {
00925 Double_t s1 = 0, s2 = 0;
00926 if (mat->InheritsFrom("TGeoMixture")) {
00927 TGeoMixture *mix = (TGeoMixture *) mat;
00928 Int_t N = mix->GetNelements();
00929 assert(N);
00930 Double_t *A = mix->GetAmixt();
00931 Double_t *Z = mix->GetZmixt();
00932 Double_t *W = mix->GetWmixt();
00933 for (Int_t i = 0; i < N; i++) {
00934 s1 += W[i]*Z[i]/A[i];
00935 s2 += W[i]*Z[i]*TMath::Log(Z[i])/A[i];
00936 }
00937 PotI=16.e-9*TMath::Exp(0.9*s2/s1);
00938 } else {
00939 Double_t Z = mat->GetZ();
00940 PotI=16.e-9*TMath::Power(Z,0.9);
00941 }
00942 }
00943 return PotI;
00944 }
00945
00946 void StiVMCToolKit::GetVMC4Reconstruction(const Char_t *pathT, const Char_t *nameT){
00952 GetVMC();
00953 gGeoManager->RestoreMasterVolume();
00954 gGeoManager->CdTop();
00955 TString path("");
00956 if (pathT) {gGeoManager->cd(pathT); path = pathT;}
00957 else {path = gGeoManager->GetCurrentNode()->GetName();}
00958 TGeoNode *nodeT = gGeoManager->GetCurrentNode();
00959 if (! nodeT) return;
00960 #if 1
00961 LoopOverNodes(nodeT, path, nameT);
00962 #else
00963 for (Int_t i = 0; i < NofVolToBEAveraged; i++) {
00964 LoopOverNodes(nodeT, path, VolumesToBeAveraged[i].name);
00965 }
00966 #endif
00967 }
00968 #ifndef __ROOT__
00969
00970 void TestVMCTK() {
00973
00974 TString path("HALL_1/CAVE_1/SVTT_1");
00975 Double_t weight = 1.e-3*StiVMCToolKit::GetWeight(0,path, 0, 0);
00976 cout << "StiVMCToolKit::TestVMCTK() -I- total weight for "
00977 << path.Data() << "\t" << weight << "[kg]" << endl;
00978 StiVMCToolKit::TestVMC4Reconstruction();
00979 }
00980 #endif
00981
00982 Double_t StiVMCToolKit::Nice(Double_t phi) {
00983 while (phi < 2*TMath::Pi()) phi += 2*TMath::Pi();
00984 while (phi >= 2*TMath::Pi()) phi -= 2*TMath::Pi();
00985 return phi;
00986 }
