|
|
|
00001 #include "Validation/Geometry/interface/MaterialBudgetData.h" 00002 00003 #include "FWCore/MessageLogger/interface/MessageLogger.h" 00004 00005 #include "G4Step.hh" 00006 #include "G4Material.hh" 00007 #include "G4EventManager.hh" 00008 #include "G4Event.hh" 00009 00010 //#define TREE_DEBUG 00011 00012 MaterialBudgetData::MaterialBudgetData() 00013 { 00014 //instantiate categorizer to assign an ID to volumes and materials 00015 myMaterialBudgetCategorizer = 0; 00016 allStepsToTree = false; 00017 densityConvertionFactor = 6.24E18; 00018 } 00019 00020 MaterialBudgetData::~MaterialBudgetData() { 00021 } 00022 00023 void MaterialBudgetData::SetAllStepsToTree() 00024 { 00025 allStepsToTree = true; 00026 MAXNUMBERSTEPS = 0; 00027 MAXNUMBERSTEPS = 10000; 00028 theDmb = new float[MAXNUMBERSTEPS]; 00029 theDil = new float[MAXNUMBERSTEPS]; 00030 // rr 00031 theSupportDmb = new float[MAXNUMBERSTEPS]; 00032 theSensitiveDmb = new float[MAXNUMBERSTEPS]; 00033 theCablesDmb = new float[MAXNUMBERSTEPS]; 00034 theCoolingDmb = new float[MAXNUMBERSTEPS]; 00035 theElectronicsDmb = new float[MAXNUMBERSTEPS]; 00036 theOtherDmb = new float[MAXNUMBERSTEPS]; 00037 theAirDmb = new float[MAXNUMBERSTEPS]; 00038 theSupportDil = new float[MAXNUMBERSTEPS]; 00039 theSensitiveDil = new float[MAXNUMBERSTEPS]; 00040 theCablesDil = new float[MAXNUMBERSTEPS]; 00041 theCoolingDil = new float[MAXNUMBERSTEPS]; 00042 theElectronicsDil = new float[MAXNUMBERSTEPS]; 00043 theOtherDil = new float[MAXNUMBERSTEPS]; 00044 theAirDil = new float[MAXNUMBERSTEPS]; 00045 // rr 00046 theInitialX = new double[MAXNUMBERSTEPS]; 00047 theInitialY = new double[MAXNUMBERSTEPS]; 00048 theInitialZ = new double[MAXNUMBERSTEPS]; 00049 theFinalX = new double[MAXNUMBERSTEPS]; 00050 theFinalY = new double[MAXNUMBERSTEPS]; 00051 theFinalZ = new double[MAXNUMBERSTEPS]; 00052 // rr 00053 theVolumeID = new int[MAXNUMBERSTEPS]; 00054 theVolumeName = new std::string[MAXNUMBERSTEPS]; 00055 theVolumeCopy = new int[MAXNUMBERSTEPS]; 00056 theVolumeX = new float[MAXNUMBERSTEPS]; 00057 theVolumeY = new float[MAXNUMBERSTEPS]; 00058 theVolumeZ = new float[MAXNUMBERSTEPS]; 00059 theVolumeXaxis1 = new float[MAXNUMBERSTEPS]; 00060 theVolumeXaxis2 = new float[MAXNUMBERSTEPS]; 00061 theVolumeXaxis3 = new float[MAXNUMBERSTEPS]; 00062 theVolumeYaxis1 = new float[MAXNUMBERSTEPS]; 00063 theVolumeYaxis2 = new float[MAXNUMBERSTEPS]; 00064 theVolumeYaxis3 = new float[MAXNUMBERSTEPS]; 00065 theVolumeZaxis1 = new float[MAXNUMBERSTEPS]; 00066 theVolumeZaxis2 = new float[MAXNUMBERSTEPS]; 00067 theVolumeZaxis3 = new float[MAXNUMBERSTEPS]; 00068 theMaterialID = new int[MAXNUMBERSTEPS]; 00069 theMaterialName = new std::string[MAXNUMBERSTEPS]; 00070 theMaterialX0 = new float[MAXNUMBERSTEPS]; 00071 theMaterialLambda0 = new float[MAXNUMBERSTEPS]; 00072 theMaterialDensity = new float[MAXNUMBERSTEPS]; 00073 theStepID = new int[MAXNUMBERSTEPS]; 00074 theStepInitialPt = new float[MAXNUMBERSTEPS]; 00075 theStepInitialEta = new float[MAXNUMBERSTEPS]; 00076 theStepInitialPhi = new float[MAXNUMBERSTEPS]; 00077 theStepInitialEnergy = new float[MAXNUMBERSTEPS]; 00078 theStepInitialPx = new float[MAXNUMBERSTEPS]; 00079 theStepInitialPy = new float[MAXNUMBERSTEPS]; 00080 theStepInitialPz = new float[MAXNUMBERSTEPS]; 00081 theStepInitialBeta = new float[MAXNUMBERSTEPS]; 00082 theStepInitialGamma = new float[MAXNUMBERSTEPS]; 00083 theStepInitialMass = new float[MAXNUMBERSTEPS]; 00084 theStepFinalPt = new float[MAXNUMBERSTEPS]; 00085 theStepFinalEta = new float[MAXNUMBERSTEPS]; 00086 theStepFinalPhi = new float[MAXNUMBERSTEPS]; 00087 theStepFinalEnergy = new float[MAXNUMBERSTEPS]; 00088 theStepFinalPx = new float[MAXNUMBERSTEPS]; 00089 theStepFinalPy = new float[MAXNUMBERSTEPS]; 00090 theStepFinalPz = new float[MAXNUMBERSTEPS]; 00091 theStepFinalBeta = new float[MAXNUMBERSTEPS]; 00092 theStepFinalGamma = new float[MAXNUMBERSTEPS]; 00093 theStepFinalMass = new float[MAXNUMBERSTEPS]; 00094 theStepPreProcess = new int[MAXNUMBERSTEPS]; 00095 theStepPostProcess = new int[MAXNUMBERSTEPS]; 00096 // rr 00097 } 00098 00099 00100 void MaterialBudgetData::dataStartTrack( const G4Track* aTrack ) 00101 { 00102 // rr 00103 std::cout << "MaterialBudget Analysis of Event #" << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << std::endl; 00104 // rr 00105 00106 const G4ThreeVector& dir = aTrack->GetMomentum() ; 00107 00108 if( myMaterialBudgetCategorizer == 0) myMaterialBudgetCategorizer = new MaterialBudgetCategorizer; 00109 00110 theStepN=0; 00111 theTotalMB=0; 00112 theTotalIL=0; 00113 theEta=0; 00114 thePhi=0; 00115 00116 // rr 00117 theID=0; 00118 thePt=0; 00119 theEnergy=0; 00120 theMass=0; 00121 00122 theSupportMB = 0.; 00123 theSensitiveMB = 0.; 00124 theCablesMB = 0.; 00125 theCoolingMB = 0.; 00126 theElectronicsMB = 0.; 00127 theOtherMB = 0.; 00128 theAirMB = 0.; 00129 theSupportIL = 0.; 00130 theSensitiveIL = 0.; 00131 theCablesIL = 0.; 00132 theCoolingIL = 0.; 00133 theElectronicsIL = 0.; 00134 theOtherIL = 0.; 00135 theAirIL = 0.; 00136 theSupportFractionMB = 0.; 00137 theSensitiveFractionMB = 0.; 00138 theCablesFractionMB = 0.; 00139 theCoolingFractionMB = 0.; 00140 theElectronicsFractionMB = 0.; 00141 theOtherFractionMB = 0.; 00142 theAirFractionMB = 0.; 00143 theSupportFractionIL = 0.; 00144 theSensitiveFractionIL = 0.; 00145 theCablesFractionIL = 0.; 00146 theCoolingFractionIL = 0.; 00147 theElectronicsFractionIL = 0.; 00148 theOtherFractionIL = 0.; 00149 theAirFractionIL = 0.; 00150 // rr 00151 00152 theID = (int)(aTrack->GetDefinition()->GetPDGEncoding()); 00153 thePt = dir.perp(); 00154 if( dir.theta() != 0 ) { 00155 theEta = dir.eta(); 00156 } else { 00157 theEta = -99; 00158 } 00159 // thePhi = dir.phi()/deg; // better not to store in deg 00160 thePhi = dir.phi(); 00161 theEnergy = aTrack->GetTotalEnergy(); 00162 theMass = aTrack->GetDefinition()->GetPDGMass(); 00163 00164 } 00165 00166 00167 void MaterialBudgetData::dataEndTrack( const G4Track* aTrack ) 00168 { 00169 //- std::cout << "[OVAL] MaterialBudget " << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " " << theEta << " " << thePhi << " " << theTotalMB << std::endl; 00170 // rr 00171 std::cout << "Recorded steps " << theStepN << std::endl; 00172 std::cout << " Material Budget: Radiation Length " << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " eta " << theEta << " phi " << thePhi << " total X " << theTotalMB << " SUP " << theSupportMB << " SEN " << theSensitiveMB << " CAB " << theCablesMB << " COL " << theCoolingMB << " ELE " << theElectronicsMB << " other " << theOtherMB << " Air " << theAirMB << std::endl; 00173 std::cout << " Material Budget: Interaction Length " << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " eta " << theEta << " phi " << thePhi << " total L " << theTotalIL << " SUP " << theSupportIL << " SEN " << theSensitiveIL << " CAB " << theCablesIL << " COL " << theCoolingIL << " ELE " << theElectronicsIL << " other " << theOtherIL << " Air " << theAirIL << std::endl; 00174 // rr 00175 } 00176 00177 00178 void MaterialBudgetData::dataPerStep( const G4Step* aStep ) 00179 { 00180 G4Material * theMaterialPre = aStep->GetPreStepPoint()->GetMaterial(); 00181 // G4Material * theMaterialPost = aStep->GetPostStepPoint()->GetMaterial(); 00182 00183 G4StepPoint* prePoint = aStep->GetPreStepPoint(); 00184 G4StepPoint* postPoint = aStep->GetPostStepPoint(); 00185 00186 CLHEP::Hep3Vector prePos = prePoint->GetPosition(); 00187 CLHEP::Hep3Vector postPos = postPoint->GetPosition(); 00188 00189 G4double steplen = aStep->GetStepLength(); 00190 00191 G4double radlen; 00192 G4double intlen; 00193 G4double density; 00194 00195 radlen = theMaterialPre->GetRadlen(); 00196 intlen = theMaterialPre->GetNuclearInterLength(); 00197 density = theMaterialPre->GetDensity() / densityConvertionFactor; // always g/cm3 00198 00199 G4String name = theMaterialPre->GetName(); 00200 // std::cout << " steplen " << steplen << " radlen " << radlen << " mb " << steplen/radlen << " mate " << theMaterialPre->GetName() << std::endl; 00201 00202 G4LogicalVolume* lv = aStep->GetTrack()->GetVolume()->GetLogicalVolume(); 00203 00204 // instantiate the categorizer 00205 int volumeID = myMaterialBudgetCategorizer->volume( lv->GetName() ); 00206 int materialID = myMaterialBudgetCategorizer->material( lv->GetMaterial()->GetName() ); 00207 // rr 00208 std::string volumeName = lv->GetName(); 00209 std::string materialName = lv->GetMaterial()->GetName(); 00210 // rr 00211 00212 // rr 00213 /* 00214 std::cout << " Volume " << lv->GetName() << std::endl; 00215 std::cout << " Material " << lv->GetMaterial()->GetName() << std::endl; 00216 */ 00217 theSupportFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[0]; 00218 theSensitiveFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[1]; 00219 theCablesFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[2]; 00220 theCoolingFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[3]; 00221 theElectronicsFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[4]; 00222 theOtherFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[5]; 00223 theAirFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[6]; 00224 if(theOtherFractionMB!=0) std::cout << " material found with no category " << lv->GetMaterial()->GetName() 00225 << " in volume " << lv->GetName() << std::endl; 00226 theSupportFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[0]; 00227 theSensitiveFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[1]; 00228 theCablesFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[2]; 00229 theCoolingFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[3]; 00230 theElectronicsFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[4]; 00231 theOtherFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[5]; 00232 theAirFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[6]; 00233 if(theOtherFractionIL!=0) std::cout << " material found with no category " << lv->GetMaterial()->GetName() 00234 << " in volume " << lv->GetName() << std::endl; 00235 // if(theOtherFractionMB!=0) LogDebug("MaterialBudgetData") << " material found with no category " << lv->GetMaterial()->GetName() 00236 // << " in volume " << lv->GetName(); 00237 // rr 00238 00239 float dmb = steplen/radlen; 00240 float dil = steplen/intlen; 00241 00242 G4VPhysicalVolume* pv = aStep->GetPreStepPoint()->GetPhysicalVolume(); 00243 const G4VTouchable* t = aStep->GetPreStepPoint()->GetTouchable(); 00244 G4ThreeVector objectTranslation = t->GetTranslation(); 00245 const G4RotationMatrix* objectRotation = t->GetRotation(); 00246 const G4VProcess* interactionPre = prePoint->GetProcessDefinedStep(); 00247 const G4VProcess* interactionPost = postPoint->GetProcessDefinedStep(); 00248 00249 G4Track* track = aStep->GetTrack(); 00250 if(theStepN==0) std::cout << " Simulated Particle " << theID << "\tMass " << theMass << " MeV/c2" 00251 << "\tPt = " << thePt << " MeV/c" << "\tEta = " << theEta << "\tPhi = " << thePhi 00252 << "\tEnergy = " << theEnergy << " MeV" 00253 // << std::endl 00254 // << "\tMagnetic Field at (0,0,0): (" << B000[0] << "," < B000[1] << "," << B000[2] << ")" 00255 << std::endl; 00256 00257 //fill data per step 00258 if( allStepsToTree ){ 00259 if( stepN > MAXNUMBERSTEPS ) stepN = MAXNUMBERSTEPS; 00260 theDmb[theStepN] = dmb; 00261 theDil[theStepN] = dil; 00262 theSupportDmb[theStepN] = (dmb * theSupportFractionMB); 00263 theSensitiveDmb[theStepN] = (dmb * theSensitiveFractionMB); 00264 theCablesDmb[theStepN] = (dmb * theCablesFractionMB); 00265 theCoolingDmb[theStepN] = (dmb * theCoolingFractionMB); 00266 theElectronicsDmb[theStepN] = (dmb * theElectronicsFractionMB); 00267 theOtherDmb[theStepN] = (dmb * theOtherFractionMB); 00268 theAirDmb[theStepN] = (dmb * theAirFractionMB); 00269 theSupportDil[theStepN] = (dil * theSupportFractionIL); 00270 theSensitiveDil[theStepN] = (dil * theSensitiveFractionIL); 00271 theCablesDil[theStepN] = (dil * theCablesFractionIL); 00272 theCoolingDil[theStepN] = (dil * theCoolingFractionIL); 00273 theElectronicsDil[theStepN] = (dil * theElectronicsFractionIL); 00274 theOtherDil[theStepN] = (dil * theOtherFractionIL); 00275 theAirDil[theStepN] = (dil * theAirFractionIL); 00276 theInitialX[theStepN] = prePos.x(); 00277 theInitialY[theStepN] = prePos.y(); 00278 theInitialZ[theStepN] = prePos.z(); 00279 theFinalX[theStepN] = postPos.x(); 00280 theFinalY[theStepN] = postPos.y(); 00281 theFinalZ[theStepN] = postPos.z(); 00282 theVolumeID[theStepN] = volumeID; 00283 theVolumeName[theStepN] = volumeName; 00284 theVolumeCopy[theStepN] = pv->GetCopyNo(); 00285 theVolumeX[theStepN] = objectTranslation.x(); 00286 theVolumeY[theStepN] = objectTranslation.y(); 00287 theVolumeZ[theStepN] = objectTranslation.z(); 00288 theVolumeXaxis1[theStepN] = objectRotation->xx(); 00289 theVolumeXaxis2[theStepN] = objectRotation->xy(); 00290 theVolumeXaxis3[theStepN] = objectRotation->xz(); 00291 theVolumeYaxis1[theStepN] = objectRotation->yx(); 00292 theVolumeYaxis2[theStepN] = objectRotation->yy(); 00293 theVolumeYaxis3[theStepN] = objectRotation->yz(); 00294 theVolumeZaxis1[theStepN] = objectRotation->zx(); 00295 theVolumeZaxis2[theStepN] = objectRotation->zy(); 00296 theVolumeZaxis3[theStepN] = objectRotation->zz(); 00297 theMaterialID[theStepN] = materialID; 00298 theMaterialName[theStepN] = materialName; 00299 theMaterialX0[theStepN] = radlen; 00300 theMaterialLambda0[theStepN] = intlen; 00301 theMaterialDensity[theStepN] = density; 00302 theStepID[theStepN] = track->GetDefinition()->GetPDGEncoding(); 00303 theStepInitialPt[theStepN] = prePoint->GetMomentum().perp(); 00304 theStepInitialEta[theStepN] = prePoint->GetMomentum().eta(); 00305 theStepInitialPhi[theStepN] = prePoint->GetMomentum().phi(); 00306 theStepInitialEnergy[theStepN] = prePoint->GetTotalEnergy(); 00307 theStepInitialPx[theStepN] = prePoint->GetMomentum().x(); 00308 theStepInitialPy[theStepN] = prePoint->GetMomentum().y(); 00309 theStepInitialPz[theStepN] = prePoint->GetMomentum().z(); 00310 theStepInitialBeta[theStepN] = prePoint->GetBeta(); 00311 theStepInitialGamma[theStepN] = prePoint->GetGamma(); 00312 theStepInitialMass[theStepN] = prePoint->GetMass(); 00313 theStepFinalPt[theStepN] = postPoint->GetMomentum().perp(); 00314 theStepFinalEta[theStepN] = postPoint->GetMomentum().eta(); 00315 theStepFinalPhi[theStepN] = postPoint->GetMomentum().phi(); 00316 theStepFinalEnergy[theStepN] = postPoint->GetTotalEnergy(); 00317 theStepFinalPx[theStepN] = postPoint->GetMomentum().x(); 00318 theStepFinalPy[theStepN] = postPoint->GetMomentum().y(); 00319 theStepFinalPz[theStepN] = postPoint->GetMomentum().z(); 00320 theStepFinalBeta[theStepN] = postPoint->GetBeta(); 00321 theStepFinalGamma[theStepN] = postPoint->GetGamma(); 00322 theStepFinalMass[theStepN] = postPoint->GetMass(); 00323 int preProcType = -99; 00324 int postProcType = -99; 00325 if (interactionPre) preProcType = interactionPre->GetProcessType(); 00326 theStepPreProcess[theStepN] = preProcType; 00327 if (interactionPost) postProcType = interactionPost->GetProcessType(); 00328 theStepPostProcess[theStepN] = postProcType; 00329 #ifdef TREE_DEBUG 00330 std::cout << " step " << theStepN 00331 << "\tDelta MB = " << theDmb[theStepN] 00332 << std::endl 00333 << "\t\tSupport " << theSupportDmb[theStepN] 00334 << " Sensitive " << theSensitiveDmb[theStepN] 00335 << " Cables " << theCablesDmb[theStepN] 00336 << " Cooling " << theCoolingDmb[theStepN] 00337 << " Electronics " << theElectronicsDmb[theStepN] 00338 << " Other " << theOtherDmb[theStepN] 00339 << " Air " << theAirDmb[theStepN] 00340 << std::endl 00341 << "\tDelta IL = " << theDil[theStepN] 00342 << std::endl 00343 << "\t\tSupport " << theSupportDil[theStepN] 00344 << " Sensitive " << theSensitiveDil[theStepN] 00345 << " Cables " << theCablesDil[theStepN] 00346 << " Cooling " << theCoolingDil[theStepN] 00347 << " Electronics " << theElectronicsDil[theStepN] 00348 << " Other " << theOtherDil[theStepN] 00349 << " Air " << theAirDil[theStepN] 00350 << std::endl; 00351 if (interactionPre) 00352 std::cout << "\tProcess Pre " << interactionPre->GetProcessName() 00353 << " type " << theStepPreProcess[theStepN] << " " << interactionPre->GetProcessTypeName(G4ProcessType(theStepPreProcess[theStepN])) 00354 << std::endl; 00355 if (interactionPost) 00356 std::cout << "\tProcess Post " << interactionPost->GetProcessName() 00357 << " type " << theStepPostProcess[theStepN] << " " 00358 << interactionPost->GetProcessTypeName(G4ProcessType(theStepPostProcess[theStepN])) 00359 << std::endl; 00360 std::cout << "\tPre x = " << theInitialX[theStepN] 00361 << "\ty = " << theInitialY[theStepN] 00362 << "\tz = " << theInitialZ[theStepN] 00363 << "\tPolar Radius = " << sqrt(prePos.x()*prePos.x()+prePos.y()*prePos.y()) 00364 << "\tPt = " << theStepInitialPt[theStepN] 00365 << "\tEnergy = " << theStepInitialEnergy[theStepN] 00366 // << std::endl 00367 // << "B-field(T) at Pre (cm): " << field->fieldInTesla(GlobalPoint(pos.x()/10.,pos.y()/10.,pos.z()/10.)) 00368 << std::endl; 00369 std::cout << "\tPost x = " << theFinalX[theStepN] 00370 << "\ty = " << theFinalY[theStepN] 00371 << "\tz = " << theFinalZ[theStepN] 00372 << "\tPolar Radius = " << sqrt(postPos.x()*postPos.x()+postPos.y()*postPos.y()) 00373 << "\tPt = " << theStepFinalPt[theStepN] 00374 << "\tEnergy = " << theStepFinalEnergy[theStepN] 00375 << std::endl; 00376 std::cout << "\tvolume " << volumeID << " " << theVolumeName[theStepN] 00377 << " copy number " << theVolumeCopy[theStepN] 00378 << "\tmaterial " << theMaterialID[theStepN] << " " << theMaterialName[theStepN] 00379 << "\tDensity = " << theMaterialDensity[theStepN] << " g/cm3" 00380 << "\tX0 = " << theMaterialX0[theStepN] << " mm" 00381 << "\tLambda0 = " << theMaterialLambda0[theStepN] << " mm" 00382 << std::endl; 00383 std::cout << "\t\tParticle " << theStepID[theStepN] 00384 << " Initial Pt = " << theStepInitialPt[theStepN] << " MeV/c" 00385 << " eta = " << theStepInitialEta[theStepN] 00386 << " phi = " << theStepInitialPhi[theStepN] 00387 << " Energy = " << theStepInitialEnergy[theStepN] << " MeV" 00388 << " Mass = " << theStepInitialMass[theStepN] << " MeV/c2" 00389 << " Beta = " << theStepInitialBeta[theStepN] 00390 << " Gamma = " << theStepInitialGamma[theStepN] 00391 << std::endl 00392 << "\t\tParticle " << theStepID[theStepN] 00393 << " Final Pt = " << theStepFinalPt[theStepN] << " MeV/c" 00394 << " eta = " << theStepFinalEta[theStepN] 00395 << " phi = " << theStepFinalPhi[theStepN] 00396 << " Energy = " << theStepFinalEnergy[theStepN] << " MeV" 00397 << " Mass = " << theStepFinalMass[theStepN] << " MeV/c2" 00398 << " Beta = " << theStepFinalBeta[theStepN] 00399 << " Gamma = " << theStepFinalGamma[theStepN] 00400 << std::endl; 00401 std:: cout << "\tVolume Centre x = " << theVolumeX[theStepN] 00402 << "\ty = " << theVolumeY[theStepN] 00403 << "\tz = " << theVolumeZ[theStepN] 00404 << "\tPolar Radius = " << sqrt( theVolumeX[theStepN]*theVolumeX[theStepN] + 00405 theVolumeY[theStepN]*theVolumeY[theStepN] ) 00406 << std::endl; 00407 std::cout << "\tx axis = (" 00408 << theVolumeXaxis1[theStepN] << "," 00409 << theVolumeXaxis2[theStepN] << "," 00410 << theVolumeXaxis3[theStepN] << ")" 00411 << std::endl; 00412 std::cout << "\ty axis = (" 00413 << theVolumeYaxis1[theStepN] << "," 00414 << theVolumeYaxis2[theStepN] << "," 00415 << theVolumeYaxis3[theStepN] << ")" 00416 << std::endl; 00417 std::cout << "\tz axis = (" 00418 << theVolumeZaxis1[theStepN] << "," 00419 << theVolumeZaxis2[theStepN] << "," 00420 << theVolumeZaxis3[theStepN] << ")" 00421 << std::endl; 00422 std::cout << "\tSecondaries" 00423 << std::endl; 00424 for(G4TrackVector::iterator iSec = aStep->GetSecondary()->begin(); iSec!=aStep->GetSecondary()->end(); iSec++) { 00425 std::cout << "\t\tid " << (*iSec)->GetDefinition()->GetPDGEncoding() 00426 << " created through process " 00427 << " type " << (*iSec)->GetCreatorProcess()->GetProcessType() 00428 << " " << (*iSec)->GetCreatorProcess()->GetProcessTypeName(G4ProcessType((*iSec)->GetCreatorProcess()->GetProcessType())) 00429 << std::endl; 00430 } 00431 #endif 00432 } 00433 00434 theTrkLen = aStep->GetTrack()->GetTrackLength(); 00435 //- std::cout << " theTrkLen " << theTrkLen << " theTrkLen2 " << theTrkLen2 << " postPos " << postPos.mag() << postPos << std::endl; 00436 thePVname = pv->GetName(); 00437 thePVcopyNo = pv->GetCopyNo(); 00438 theRadLen = radlen; 00439 theIntLen = intlen; 00440 theTotalMB += dmb; 00441 theTotalIL += dil; 00442 00443 // rr 00444 theSupportMB += (dmb * theSupportFractionMB); 00445 theSensitiveMB += (dmb * theSensitiveFractionMB); 00446 theCablesMB += (dmb * theCablesFractionMB); 00447 theCoolingMB += (dmb * theCoolingFractionMB); 00448 theElectronicsMB += (dmb * theElectronicsFractionMB); 00449 theOtherMB += (dmb * theOtherFractionMB); 00450 theAirMB += (dmb * theAirFractionMB); 00451 theSupportIL += (dil * theSupportFractionIL); 00452 theSensitiveIL += (dil * theSensitiveFractionIL); 00453 theCablesIL += (dil * theCablesFractionIL); 00454 theCoolingIL += (dil * theCoolingFractionIL); 00455 theElectronicsIL += (dil * theElectronicsFractionIL); 00456 theOtherIL += (dil * theOtherFractionIL); 00457 theAirIL += (dil * theAirFractionIL); 00458 // rr 00459 00460 theStepN++; 00461 00462 } 00463 00464
1.7.3