source: trunk/source/processes/electromagnetic/lowenergy/test/G4PenelopeComptonTest.cc@ 1201

Last change on this file since 1201 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 23.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4PenelopeComptonTest.cc,v 1.10 2008/03/27 13:44:04 pandola Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// -------------------------------------------------------------------
31// GEANT 4 class file --- Copyright CERN 1998
32// CERN Geneva Switzerland
33//
34//
35// File name: G4PenelopeComptonTest.cc
36//
37// Author: Francesco Longo
38//
39// Creation date: 04 january 2001
40//
41// Modifications: Luciano Pandola (27 november 2002)
42// Adapted in order to test G4PenelopeCompton
43// Minor modification in n-tuple filling
44// Updated analysis to AIDA 3.0
45//
46// -------------------------------------------------------------------
47
48#include "globals.hh"
49#include "G4ios.hh"
50#include <fstream>
51#include <iomanip>
52
53#include "G4ParticleDefinition.hh"
54#include "G4ParticleTypes.hh"
55#include "G4ParticleTable.hh"
56#include "G4Material.hh"
57#include "G4MaterialTable.hh"
58#include "G4VDiscreteProcess.hh"
59#include "G4VProcess.hh"
60#include "G4ProcessManager.hh"
61#include "G4LowEnergyCompton.hh"
62#include "G4PenelopeCompton.hh"
63#include "G4EnergyLossTables.hh"
64#include "G4VParticleChange.hh"
65#include "G4ParticleChange.hh"
66#include "G4DynamicParticle.hh"
67#include "G4ForceCondition.hh"
68#include "G4RunManager.hh"
69#include "G4RegionStore.hh"
70#include "G4LowEnergyBremsstrahlung.hh"
71#include "G4LowEnergyIonisation.hh"
72#include "G4eIonisation.hh"
73#include "G4MultipleScattering.hh"
74#include "G4eIonisation.hh"
75#include "G4eBremsstrahlung.hh"
76#include "G4eplusAnnihilation.hh"
77#include "G4StateManager.hh"
78#include "G4ApplicationState.hh"
79#include "G4ComptonScattering52.hh"
80#include "G4PhotoElectricEffect.hh"
81
82#include "G4Electron.hh"
83#include "G4Positron.hh"
84#include "G4Gamma.hh"
85
86#include "G4GRSVolume.hh"
87#include "G4Box.hh"
88#include "G4PVPlacement.hh"
89#include "G4Step.hh"
90#include "G4ProductionCutsTable.hh"
91#include "G4MaterialCutsCouple.hh"
92
93#include "G4UnitsTable.hh"
94#include "AIDA/IManagedObject.h"
95
96#include <memory>
97#include "AIDA/IAnalysisFactory.h"
98#include "AIDA/ITreeFactory.h"
99#include "AIDA/ITree.h"
100#include "AIDA/IHistogramFactory.h"
101#include "AIDA/IHistogram1D.h"
102#include "AIDA/IHistogram2D.h"
103#include "AIDA/IHistogram3D.h"
104#include "AIDA/ITupleFactory.h"
105#include "AIDA/ITuple.h"
106
107
108G4int main()
109{
110
111 // Setup
112
113 G4int nIterations = 100000;
114 G4int materialId = 3;
115
116 //G4cout.setf(std::ios::scientific,std::ios::floatfield );
117
118 // -------------------------------------------------------------------
119
120 // ---- HBOOK initialization
121
122 G4String fileName;
123 G4cout << "Filename?" << G4endl;
124 G4cin >> fileName;
125
126 //build up the factories
127 AIDA::IAnalysisFactory *af = AIDA_createAnalysisFactory();
128 G4cout << "Analysis factory created" << G4endl;
129 //parameters for the TreeFactory
130 G4bool fileExists = true;
131 G4bool readOnly = false;
132 AIDA::ITreeFactory *tf = af->createTreeFactory();
133 AIDA::ITree *tree = tf->create(fileName,"hbook",readOnly,fileExists);
134 G4cout << "Tree store : " << tree->storeName() << G4endl;
135 G4cout << " Booked Hbook File " << G4endl;
136
137 //AIDA::IHistogramFactory *hf = af->createHistogramFactory(*tree);
138 AIDA::ITupleFactory *tpf = af->createTupleFactory(*tree );
139
140 // ---- primary ntuple ------
141 AIDA::ITuple* ntuple1 = tpf->create("1",
142 "Primary Ntuple","double eprimary,energyf,de,eloc,dedx,pxch,pych,pzch,pch,thetach");
143
144 // ---- secondary ntuple ------
145 AIDA::ITuple* ntuple2 = tpf->create("2",
146 "Secondary Ntuple","double eprimary,p_el,e_el,theta_el,costheta_el,ekin_el,ekin_rel");
147
148 // ---- table ntuple ------
149 AIDA::ITuple* ntuple3 = tpf->create("3","Mean Free Path Ntuple","double kinen,mfp");
150
151 // ---- fluorescence ntuple -------
152 AIDA::ITuple* ntuple4 = tpf->create("4","Fluorescence Ntuple","double eprimary,px,py,pz,e,partType");
153
154 //--------- Materials definition ---------
155 G4Element* SiEl = new G4Element ("SiEl","Si",14.,28.055*g/mole);
156 G4Element* FeEl = new G4Element ("FeEl","Fe",26.,55.58*g/mole);
157 G4Element* CuEl = new G4Element ("CuEl","Cu",29.,63.55*g/mole);
158 G4Element* WEl = new G4Element ("WEl","W",74.,183.85*g/mole);
159 G4Element* PbEl = new G4Element ("PbEl","Pb",82.,207.19*g/mole);
160 G4Element* UEl = new G4Element ("UEl","U",92.,238.03*g/mole);
161 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
162 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
163 G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
164 G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
165 G4Element* I = new G4Element ("Iodine" , "I", 53. , 126.9044*g/mole);
166 G4Element* N = new G4Element ("Nitrogen","N",7.,14.0*g/mole);
167 G4Element* AuEl = new G4Element ("Gold","Au",79.,196.97*g/mole);
168 G4Element* AlEl = new G4Element ("Aluminum","Al",13,26.981*g/mole);
169
170 G4Material* Si = new G4Material("Silicon", 2.33*g/cm3,1);
171 Si->AddElement(SiEl,1);
172 G4Material* Fe = new G4Material("Iron",7.87*g/cm3,1);
173 Fe->AddElement(FeEl,1);
174 G4Material* Cu = new G4Material("Copper", 8.96*g/cm3,1);
175 Cu->AddElement(CuEl,1);
176 G4Material* W = new G4Material("Tungsten",19.30*g/cm3,1);
177 W->AddElement(WEl,1);
178 G4Material* Pb = new G4Material("Lead",11.35*g/cm3,1);
179 Pb->AddElement(PbEl,1);
180 G4Material* U = new G4Material("Uranium",18.95*g/cm3,1);
181 U->AddElement(UEl,1);
182 G4Material* maO = new G4Material("Oxygen", 1.1*g/cm3,1);
183 maO->AddElement(O,2);
184 G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
185 water->AddElement(H,2);
186 water->AddElement(O,1);
187 G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
188 ethane->AddElement(H,6);
189 ethane->AddElement(C,2);
190 G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
191 csi->AddElement(Cs,1);
192 csi->AddElement(I,1);
193 G4Material* Air = new G4Material
194 ("Air", 1.2929*kg/m3, 2, kStateGas, 300.00*kelvin, 1.0*atmosphere);
195 Air->AddElement(N,0.8);
196 Air->AddElement(O,0.2);
197 G4Material* Gold = new G4Material("Gold",19.3*g/cm3,1);
198 Gold->AddElement(AuEl,1);
199 G4Material* Al = new G4Material("Aluminum",2.70*g/cm3,1);
200 Al->AddElement(AlEl,1);
201
202
203 // Interactive set-up
204
205 G4cout << "How many interactions? " << G4endl;
206 G4cin >> nIterations;
207
208 if (nIterations <= 0) G4Exception("Wrong input");
209
210 G4double initEnergy = 10*keV;
211 G4double initX = 0.;
212 G4double initY = 0.;
213 G4double initZ = 1.;
214
215 G4cout << "Enter the initial particle energy E (MeV)" << G4endl;
216 G4cin >> initEnergy ;
217
218 initEnergy = initEnergy*MeV;
219
220 if (initEnergy <= 0.) G4Exception("Wrong input");
221
222 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
223
224 G4int nMaterials = G4Material::GetNumberOfMaterials();
225
226 G4cout << "Available materials are: " << G4endl;
227 for (G4int mat = 0; mat < nMaterials; mat++)
228 {
229 G4cout << mat << ") "
230 << (*theMaterialTable)[mat]->GetName()
231 << G4endl;
232 }
233
234 G4cout << "Which material? " << G4endl;
235 G4cin >> materialId;
236
237 G4Material* material = (*theMaterialTable)[materialId] ;
238
239 G4cout << "The selected material is: "
240 << material->GetName()
241 << G4endl;
242
243 G4double dimX = 1*mm;
244 G4double dimY = 1*mm;
245 G4double dimZ = 1*mm;
246
247 // Geometry
248
249 G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
250
251 G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
252 (*theMaterialTable)[materialId],
253 "LFrame", 0, 0, 0);
254 logicalFrame->SetMaterial(material);
255
256 G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
257 "PFrame",logicalFrame,0,false,0);
258
259 G4RunManager* rm = new G4RunManager();
260
261 rm->GeometryHasBeenModified();
262 rm->DefineWorldVolume(physicalFrame);
263
264 G4cout << "[OK] World is defined " << G4endl;
265 G4StateManager::GetStateManager()->SetNewState(G4State_Idle);
266 rm->DumpRegion("DefaultRegionForTheWorld"); //this forces the region update!
267
268 // Particle definitions
269
270 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
271 G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
272 G4ParticleDefinition* positron = G4Positron::PositronDefinition();
273
274 G4ProductionCutsTable* cutsTable = G4ProductionCutsTable::GetProductionCutsTable();
275 G4ProductionCuts* cuts = cutsTable->GetDefaultProductionCuts();
276 G4double cutG=1*nanometer;
277 G4double cutE=1*nanometer;
278 cuts->SetProductionCut(cutG, 0); //gammas
279 cuts->SetProductionCut(cutE, 1); //electrons
280 cuts->SetProductionCut(cutE, 2); //positrons
281 G4cout << "Cuts are defined " << G4endl;
282
283 //G4Gamma::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
284 //G4Electron::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
285 //G4Positron::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
286
287 G4MaterialCutsCouple* theCouple = new G4MaterialCutsCouple(material,cuts);
288
289 logicalFrame->SetMaterialCutsCouple(theCouple);
290
291 //G4cout << "La coppia: " << theCouple << G4endl;
292 G4cout << "Recalculation needed: " << theCouple->IsRecalcNeeded() << G4endl;
293
294 cutsTable->UpdateCoupleTable(physicalFrame);
295 cutsTable->PhysicsTableUpdated();
296 cutsTable->DumpCouples();
297
298
299 // Processes
300
301 G4int processType;
302 G4cout << "LowEnergy[1], Penelope [2] or Standard-5.2 [3] Compton?" << G4endl;
303 G4cin >> processType;
304 if ( !(processType == 1 || processType == 2 || processType == 3))
305 {
306 G4Exception("Wrong input");
307 }
308
309 G4VDiscreteProcess* gammaProcess = NULL;
310
311 if (processType == 1)
312 {
313 gammaProcess = new G4LowEnergyCompton();
314 G4cout << "The selected model is Low Energy" << G4endl;
315 }
316 else if (processType == 2)
317 {
318 G4PenelopeCompton* thePenelopeCompton = new G4PenelopeCompton();
319 //thePenelopeCompton->SetUseAtomicDeexcitation(false);
320 gammaProcess = thePenelopeCompton;
321 G4cout << "The selected model is Penelope" << G4endl;
322 //G4cout << "Fluorescence is disabled" << G4endl;
323 }
324 else if (processType == 3)
325 {
326 gammaProcess = new G4ComptonScattering52();
327 G4cout << "The selected model is Standard" << G4endl;
328 }
329
330 G4VProcess* theeminusMultipleScattering = new G4MultipleScattering();
331 G4VProcess* theeminusIonisation = new G4eIonisation();
332 G4VProcess* theeminusBremsstrahlung = new G4eBremsstrahlung();
333 G4VProcess* theeplusMultipleScattering = new G4MultipleScattering();
334 G4VProcess* theeplusIonisation = new G4eIonisation();
335 G4VProcess* theeplusBremsstrahlung = new G4eBremsstrahlung();
336 G4VProcess* theeplusAnnihilation = new G4eplusAnnihilation();
337
338 //----------------
339 // process manager
340 //----------------
341
342 // gamma
343
344 G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
345 gamma->SetProcessManager(gProcessManager);
346 gProcessManager->AddProcess(gammaProcess);
347 G4ForceCondition* condition; //l'ho fissata a zero! E' onesto??
348
349 //electron
350
351 G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
352 electron->SetProcessManager(eProcessManager);
353 eProcessManager->AddProcess(theeminusMultipleScattering);
354 eProcessManager->AddProcess(theeminusIonisation);
355 eProcessManager->AddProcess(theeminusBremsstrahlung);
356
357 //positron
358
359 G4ProcessManager* pProcessManager = new G4ProcessManager(positron);
360 positron->SetProcessManager(pProcessManager);
361 pProcessManager->AddProcess(theeplusMultipleScattering);
362 pProcessManager->AddProcess(theeplusIonisation);
363 pProcessManager->AddProcess(theeplusBremsstrahlung);
364 pProcessManager->AddProcess(theeplusAnnihilation);
365
366 //--------------
367 // set ordering
368 //--------------
369
370
371 eProcessManager->
372 SetProcessOrdering(theeminusMultipleScattering, idxAlongStep,1);
373 eProcessManager->
374 SetProcessOrdering(theeminusIonisation, idxAlongStep,2);
375
376 eProcessManager->
377 SetProcessOrdering(theeminusMultipleScattering, idxPostStep,1);
378 eProcessManager->
379 SetProcessOrdering(theeminusIonisation, idxPostStep,2);
380 eProcessManager->
381 SetProcessOrdering(theeminusBremsstrahlung, idxPostStep,3);
382
383
384
385 pProcessManager->SetProcessOrderingToFirst(theeplusAnnihilation, idxAtRest);
386 pProcessManager->
387 SetProcessOrdering(theeplusMultipleScattering, idxAlongStep,1);
388 pProcessManager->
389 SetProcessOrdering(theeplusIonisation, idxAlongStep,2);
390
391 pProcessManager->
392 SetProcessOrdering(theeplusMultipleScattering, idxPostStep,1);
393 pProcessManager->
394 SetProcessOrdering(theeplusIonisation, idxPostStep,2);
395 pProcessManager->
396 SetProcessOrdering(theeplusBremsstrahlung, idxPostStep,3);
397 pProcessManager->
398 SetProcessOrdering(theeplusAnnihilation, idxPostStep,4);
399
400 // G4LowEnergyIonisation IonisationProcess;
401 // eProcessManager->AddProcess(&IonisationProcess);
402 // eProcessManager->SetProcessOrdering(&IonisationProcess,idxAlongStep,1);
403 // eProcessManager->SetProcessOrdering(&IonisationProcess,idxPostStep, 1);
404
405 // G4LowEnergyBremsstrahlung BremstrahlungProcess;
406 // eProcessManager->AddProcess(&BremstrahlungProcess);
407 // eProcessManager->SetProcessOrdering(&BremstrahlungProcess,idxAlongStep,1);
408 // eProcessManager->SetProcessOrdering(&BremstrahlungProcess,idxPostStep, 1);
409
410 // G4eIonisation IonisationPlusProcess;
411 // pPositronProcessManager->AddProcess(&IonisationPlusProcess);
412 // pProcessManager->
413 // SetProcessOrdering(&IonisationPlusProcess,idxAlongStep,1);
414 // pProcessManager->SetProcessOrdering(&IonisationPlusProcess,idxPostStep,1);
415
416
417
418 // Create a DynamicParticle
419
420 G4double eEnergy = initEnergy*MeV;
421 G4ParticleMomentum eDirection(initX,initY,initZ);
422 G4DynamicParticle dynamicGamma(G4Gamma::Gamma(),eDirection,eEnergy);
423
424 //dynamicGamma.DumpInfo(0);
425
426 // Track
427
428 G4ThreeVector aPosition(0.,0.,0.);
429 G4double aTime = 0. ;
430
431 G4Track* gTrack = new G4Track(&dynamicGamma,aTime,aPosition);
432
433 G4GRSVolume* touche = new G4GRSVolume(physicalFrame, NULL, aPosition);
434 //gTrack->SetTouchableHandle(touche); //verificare!!!!!!!!!!!!
435
436
437 // Step
438
439 G4Step* step = new G4Step();
440 step->SetTrack(gTrack);
441
442 G4StepPoint* aPoint = new G4StepPoint();
443 aPoint->SetPosition(aPosition);
444 aPoint->SetMaterial(material);
445 aPoint->SetMaterialCutsCouple(theCouple);
446 G4double safety = 10000.*cm;
447 aPoint->SetSafety(safety);
448 step->SetPreStepPoint(aPoint);
449
450 // PostStep point
451 G4StepPoint* newPoint = new G4StepPoint();
452 G4ThreeVector newPosition(0.,0.,1.*mm);
453 newPoint->SetPosition(newPosition);
454 newPoint->SetMaterial(material);
455 newPoint->SetMaterialCutsCouple(theCouple);
456 newPoint->SetSafety(safety);
457 step->SetPostStepPoint(newPoint);
458
459 // Check applicability
460
461 if (! (gammaProcess->IsApplicable(*gamma)))
462 {
463 G4Exception("Not Applicable");
464 }
465 else
466 {
467 G4cout<< "applicability OK" << G4endl;
468 }
469
470 // Initialize the physics tables (in which material?)
471 gammaProcess->BuildPhysicsTable(*gamma);
472
473 /*
474 theeminusMultipleScattering->BuildPhysicsTable(*electron);
475 theeminusIonisation->BuildPhysicsTable(*electron);
476 theeminusBremsstrahlung->BuildPhysicsTable(*electron);
477 theeplusMultipleScattering->BuildPhysicsTable(*positron);
478 theeplusIonisation->BuildPhysicsTable(*positron);
479 theeplusBremsstrahlung->BuildPhysicsTable(*positron);
480 theeplusAnnihilation->BuildPhysicsTable(*positron) ;
481 */
482
483 G4cout<< "table OK" << G4endl;
484
485 G4Material* apttoMaterial ;
486 G4String MaterialName ;
487
488 G4double minArg = 100*eV,maxArg = 100*GeV, argStp;
489 const G4int pntNum = 300;
490 G4double Tkin[pntNum+1];
491 G4double meanFreePath=0. ;
492
493 argStp = (std::log10(maxArg)-std::log10(minArg))/pntNum;
494
495 for(G4int d = 0; d < pntNum+1; d++)
496 {
497 Tkin[d] = std::pow(10,(std::log10(minArg) + d*argStp));
498 }
499
500 G4double sti = 1.*mm;
501 step->SetStepLength(sti);
502
503 apttoMaterial = (*theMaterialTable)[materialId] ;
504 MaterialName = apttoMaterial->GetName() ;
505 logicalFrame->SetMaterial(apttoMaterial);
506
507 gTrack->SetStep(step);
508
509 G4LowEnergyCompton* gammaLowEProcess2 = (G4LowEnergyCompton*) gammaProcess;
510 G4PenelopeCompton* gammaLowEProcess = (G4PenelopeCompton*) gammaProcess;
511
512 for (G4int i=0 ; i<pntNum; i++)
513 {
514 dynamicGamma.SetKineticEnergy(Tkin[i]);
515 if (processType == 1)
516 {
517 meanFreePath=gammaLowEProcess2
518 ->DumpMeanFreePath(*gTrack, sti, condition);
519 }
520 else if (processType == 2)
521 {
522 meanFreePath=gammaLowEProcess
523 ->DumpMeanFreePath(*gTrack, sti, condition);
524 }
525
526 ntuple3->fill(ntuple3->findColumn("kinen"),std::log10(Tkin[i]));
527 ntuple3->fill(ntuple3->findColumn("mfp"),std::log10(meanFreePath/cm));
528 ntuple3->addRow();
529 }
530 G4cout << "Mean Free Path OK" << G4endl;
531
532 // --------- Test the DoIt
533
534 G4cout << "DoIt in " << material->GetName() << G4endl;
535
536
537 dynamicGamma.SetKineticEnergy(eEnergy);
538 G4int iter;
539 for (iter=0; iter<nIterations; iter++)
540 {
541
542 step->SetStepLength(1*nanometer);
543
544 G4cout << "Iteration = " << iter
545 << " - Step Length = "
546 << step->GetStepLength()/mm << " mm "
547 << G4endl;
548
549
550 gTrack->SetStep(step);
551
552 G4ParticleChange* particleChange = (G4ParticleChange*) gammaProcess->PostStepDoIt(*gTrack,*step);
553
554 // Primary physical quantities
555
556 G4double energyChange = particleChange->GetEnergy();
557
558 G4double dedx = initEnergy - energyChange ;
559 G4double dedxNow = dedx / (step->GetStepLength());
560
561 G4ThreeVector eChange =
562 particleChange->CalcMomentum(energyChange,
563 (*particleChange->GetMomentumDirection()),
564 particleChange->GetMass());
565
566 G4double pxChange = eChange.x();
567 G4double pyChange = eChange.y();
568 G4double pzChange = eChange.z();
569 G4double pChange =
570 std::sqrt(pxChange*pxChange + pyChange*pyChange + pzChange*pzChange);
571 G4double thetaChange = eChange.theta();
572 thetaChange = thetaChange/deg; //conversion in degrees
573 G4double xChange = particleChange->GetPosition()->x();
574 G4double yChange = particleChange->GetPosition()->y();
575 G4double zChange = particleChange->GetPosition()->z();
576 G4double localEnergyDeposit = particleChange->GetLocalEnergyDeposit();
577
578 //G4double thetaChange = particleChange->GetPositionChange()->theta();
579
580 G4cout << "---- Primary after the step ---- " << G4endl;
581
582 G4cout << "Position (x,y,z) = "
583 << xChange << " "
584 << yChange << " "
585 << zChange << " "
586 << G4endl;
587
588 G4cout << "---- Energy: " << energyChange/MeV << " MeV, "
589 << "(px,py,pz): ("
590 << pxChange/MeV << ","
591 << pyChange/MeV << ","
592 << pzChange/MeV << ") MeV"
593 << G4endl;
594
595 G4cout << "---- Energy loss (dE) = " << dedx/keV << " keV" << G4endl;
596 // G4cout << "Stopping power (dE/dx)=" << dedxNow << G4endl;
597
598 ntuple1->fill(ntuple1->findColumn("eprimary"),initEnergy/MeV);
599 ntuple1->fill(ntuple1->findColumn("energyf"),energyChange/initEnergy);
600 ntuple1->fill(ntuple1->findColumn("de"),dedx/MeV);
601 ntuple1->fill(ntuple1->findColumn("eloc"),localEnergyDeposit/MeV);
602 ntuple1->fill(ntuple1->findColumn("dedx"),dedxNow/(MeV/cm));
603 ntuple1->fill(ntuple1->findColumn("pxch"),pxChange/initEnergy);
604 ntuple1->fill(ntuple1->findColumn("pych"),pyChange/initEnergy);
605 ntuple1->fill(ntuple1->findColumn("pzch"),pzChange/initEnergy);
606 ntuple1->fill(ntuple1->findColumn("pch"),pChange/initEnergy);
607 ntuple1->fill(ntuple1->findColumn("thetach"),thetaChange);
608 ntuple1->addRow();
609
610 // Secondaries physical quantities
611
612 // Secondaries
613 G4cout << " Secondaries " <<
614 particleChange->GetNumberOfSecondaries() << G4endl;
615
616 G4double px_el=0,py_el=0,pz_el=0,p_el=0,e_el=0,theta_el=0,eKin_el=0;
617 G4double costheta_el = 0;
618
619 for (G4int i = 0; i < (particleChange->GetNumberOfSecondaries()); i++)
620 {
621 // The following two items should be filled per event, not
622 // per secondary; filled here just for convenience, to avoid
623 // complicated logic to dump ntuple when there are no secondaries
624
625 G4Track* finalParticle = particleChange->GetSecondary(i) ;
626
627 G4double e = finalParticle->GetTotalEnergy();
628 G4double eKin = finalParticle->GetKineticEnergy();
629 G4double px = (finalParticle->GetMomentum()).x();
630 G4double py = (finalParticle->GetMomentum()).y();
631 G4double pz = (finalParticle->GetMomentum()).z();
632 G4double theta = (finalParticle->GetMomentum()).theta();
633 G4double cosTheta = (finalParticle->GetMomentum()).cosTheta();
634 G4double p = std::sqrt(px*px+py*py+pz*pz);
635 theta = theta/deg; //conversion in degrees
636 if (eKin > initEnergy)
637 {
638 G4cout << "WARNING: eFinal > eInit " << G4endl;
639 }
640
641 G4String particleName =
642 finalParticle->GetDefinition()->GetParticleName();
643 G4cout << "==== Final "
644 << particleName << " "
645 << "energy: " << e/MeV << " MeV, "
646 << "eKin: " << eKin/MeV << " MeV, "
647 << "(px,py,pz): ("
648 << px/MeV << ","
649 << py/MeV << ","
650 << pz/MeV << ") MeV "
651 << G4endl;
652
653 G4int partType=-1;
654 if (particleName == "e-") {
655 partType = 1;
656 px_el=px;
657 py_el=py;
658 pz_el=pz;
659 p_el=p;
660 e_el=e;
661 theta_el=theta;
662 costheta_el = cosTheta;
663 eKin_el=eKin;
664 }
665 else if (particleName == "gamma") partType = 2;
666
667
668 delete particleChange->GetSecondary(i);
669
670
671 // Fill the secondaries ntuple
672
673 // Normalize all to the energy of primary
674 // for gammas initEnergy=initP
675 if (i==0) {
676 ntuple2->fill(ntuple2->findColumn("eprimary"),initEnergy);
677 ntuple2->fill(ntuple2->findColumn("p_el"),p_el/initEnergy);
678 ntuple2->fill(ntuple2->findColumn("e_el"),e_el/(initEnergy+electron_mass_c2));
679 ntuple2->fill(ntuple2->findColumn("theta_el"),theta_el);
680 ntuple2->fill(ntuple2->findColumn("costheta_el"),costheta_el);
681 ntuple2->fill(ntuple2->findColumn("ekin_el"),eKin_el);
682 ntuple2->fill(ntuple2->findColumn("ekin_rel"),eKin_el/initEnergy);
683 ntuple2->addRow();
684 }
685 else //fluorescence
686 {
687 ntuple4->fill(ntuple4->findColumn("eprimary"),initEnergy);
688 ntuple4->fill(ntuple4->findColumn("px"),px/keV);
689 ntuple4->fill(ntuple4->findColumn("py"),py/keV);
690 ntuple4->fill(ntuple4->findColumn("pz"),pz/keV);
691 ntuple4->fill(ntuple4->findColumn("e"),eKin/keV);
692 ntuple4->fill(ntuple4->findColumn("partType"),(G4double) partType);
693 ntuple4->addRow();
694 }
695
696 }
697 particleChange->Clear();
698
699 }
700
701
702
703 G4cout << "Iteration number: " << iter << G4endl;
704
705 G4cout << "Committing.............." << G4endl;
706 tree->commit();
707 G4cout << "Closing the tree........" << G4endl;
708 tree->close();
709
710 delete step;
711 delete tpf;
712 delete tree;
713 //delete af;
714
715
716 G4cout << "END OF THE MAIN PROGRAM" << G4endl;
717 return 0;
718}
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
Note: See TracBrowser for help on using the repository browser.