source: trunk/source/processes/electromagnetic/lowenergy/test/G4LowEnergyPolarizedComptonTest.cc@ 1314

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

nvx fichiers dans CVS
File size: 22.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: G4LowEnergyPolarizedComptonTest.cc,v 1.8 2006/06/29 19:44:05 gunter 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: G4ComptonScatteringTest.cc
36//
37// Author: Francesco Longo & Gerardo Depaola
38//
39// Creation date: 23 january 2001
40//
41// Modifications:
42//
43// -------------------------------------------------------------------
44
45#include "globals.hh"
46#include "G4ios.hh"
47#include <fstream>
48#include <iomanip>
49
50#include "G4ParticleDefinition.hh"
51#include "G4ParticleTypes.hh"
52#include "G4ParticleTable.hh"
53#include "G4Material.hh"
54#include "G4MaterialTable.hh"
55#include "G4VDiscreteProcess.hh"
56#include "G4VProcess.hh"
57#include "G4ProcessManager.hh"
58
59#include "G4ComptonScattering.hh"
60#include "G4PolarizedComptonScattering.hh"
61#include "G4LowEnergyCompton.hh"
62#include "G4LowEnergyPolarizedCompton.hh"
63
64#include "G4EnergyLossTables.hh"
65#include "G4VParticleChange.hh"
66#include "G4ParticleChange.hh"
67#include "G4DynamicParticle.hh"
68
69#include "G4LowEnergyBremsstrahlung.hh"
70#include "G4LowEnergyIonisation.hh"
71#include "G4eIonisation.hh"
72#include "G4MultipleScattering.hh"
73#include "G4eIonisation.hh"
74#include "G4eBremsstrahlung.hh"
75#include "G4eplusAnnihilation.hh"
76
77#include "G4Electron.hh"
78#include "G4Positron.hh"
79#include "G4Gamma.hh"
80
81#include "G4GRSVolume.hh"
82#include "G4Box.hh"
83#include "G4PVPlacement.hh"
84#include "G4Step.hh"
85
86#include "G4UnitsTable.hh"
87
88#include "CLHEP/Hist/TupleManager.h"
89#include "CLHEP/Hist/HBookFile.h"
90#include "CLHEP/Hist/Histogram.h"
91#include "CLHEP/Hist/Tuple.h"
92
93
94HepTupleManager* hbookManager;
95
96G4int main()
97{
98
99 // Setup
100
101 G4int nIterations = 100000;
102 G4int materialId = 3;
103
104 // G4cout.setf( ios::scientific, ios::floatfield );
105
106 // -------------------------------------------------------------------
107
108 // ---- HBOOK initialization
109
110
111 hbookManager = new HBookFile("comptontest.hbook", 58);
112 assert (hbookManager != 0);
113
114 // ---- Book a histogram and ntuples
115
116 G4cout<<"Hbook file name: "<<((HBookFile*) hbookManager)->filename()<<G4endl;
117
118
119 G4double initEnergy = 1*MeV;
120 G4double initX = 0.;
121 G4double initY = 0.;
122 G4double initZ = 1.;
123
124 G4cout << "Enter the initial particle energy E (keV)" << G4endl;
125 G4cin >> initEnergy ;
126 initEnergy = initEnergy * keV;
127 G4double limit = initEnergy/keV;
128
129 G4cout << limit << G4endl;
130
131 if (initEnergy <= 0.) G4Exception("Wrong input");
132
133 // ---- primary ntuple ------
134 HepTuple* ntuple1 = hbookManager->ntuple("Primary Ntuple");
135 assert (ntuple1 != 0);
136
137 // ---- secondary ntuple ------
138 HepTuple* ntuple2 = hbookManager->ntuple("Secondary Ntuple");
139 assert (ntuple2 != 0);
140
141 /*
142 // ---- table ntuple ------
143 HepTuple* ntuple3 = hbookManager->ntuple("Mean Free Path Ntuple");
144 assert (ntuple3 != 0);
145 */
146
147 // ---- secondaries histos ----
148
149 HepHistogram* heETot;
150 heETot = hbookManager->histogram("Electron Total Energy", 100,0.,limit);
151 assert (heETot != 0);
152
153 HepHistogram* heP;
154 heP = hbookManager->histogram("Electron Momentum", 100,0.,limit);
155 assert (heP != 0);
156
157 HepHistogram* hgETot;
158 hgETot = hbookManager->histogram("Gamma Total Energy", 100,0.,limit);
159 assert (hgETot != 0);
160
161 HepHistogram* hgP;
162 hgP = hbookManager->histogram("Gamma Momentum", 100,0.,limit);
163 assert (hgP != 0);
164
165 HepHistogram* hgTheta;
166 hgTheta = hbookManager->histogram("Theta Scattered Gamma ", 100,0.,4.);
167 assert (hgTheta != 0);
168
169 HepHistogram* hgPhi;
170 hgPhi = hbookManager->histogram("Phi Scattered Gamma ", 100,-4.,4.);
171 assert (hgPhi != 0);
172
173 HepHistogram* hSumE;
174 hSumE = hbookManager->histogram("Total Energy", 100,0.,2*limit);
175 assert (hSumE != 0);
176
177 HepHistogram* hgRapp;
178 hgRapp = hbookManager->histogram("Energy Theta Relation", 100,0.,2.);
179 assert (hgRapp != 0);
180
181 HepHistogram* hNSec;
182 hNSec = hbookManager->histogram("Number of secondaries", 100,0.,10.);
183 assert (hNSec != 0);
184
185 HepHistogram* hDebug;
186 hDebug = hbookManager->histogram("Debug", 100,0.,limit);
187 assert (hDebug != 0);
188
189
190 //--------- Materials definition ---------
191
192 G4Material* Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3);
193 G4Material* Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3);
194 G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
195 G4Material* W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
196 G4Material* Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3);
197 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
198 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
199 G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
200 G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
201 G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
202 G4Material* maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3);
203
204 G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
205 csi->AddElement(Cs,1);
206 csi->AddElement(I,1);
207
208
209 // Interactive set-up
210
211 G4cout << "How many interactions? " << G4endl;
212 G4cin >> nIterations;
213
214 if (nIterations <= 0) G4Exception("Wrong input");
215
216
217
218 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
219
220 G4int nMaterials = G4Material::GetNumberOfMaterials();
221
222 G4cout << "Available materials are: " << G4endl;
223 for (G4int mat = 0; mat < nMaterials; mat++)
224 {
225 G4cout << mat << ") "
226 << (*theMaterialTable)[mat]->GetName()
227 << G4endl;
228 }
229
230 G4cout << "Which material? " << G4endl;
231 G4cin >> materialId;
232
233 G4Material* material = (*theMaterialTable)[materialId] ;
234
235 G4cout << "The selected material is: "
236 << material->GetName()
237 << G4endl;
238
239 G4double dimX = 1*mm;
240 G4double dimY = 1*mm;
241 G4double dimZ = 1*mm;
242
243 // Geometry
244
245 G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
246
247 G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
248 (*theMaterialTable)[materialId],
249 "LFrame", 0, 0, 0);
250 logicalFrame->SetMaterial(material);
251
252 G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
253 "PFrame",logicalFrame,0,false,0);
254
255 // Particle definitions
256
257 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
258 G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
259 G4ParticleDefinition* positron = G4Positron::PositronDefinition();
260
261 gamma->SetCuts(0.1*micrometer);
262 electron->SetCuts(0.1*micrometer);
263 positron->SetCuts(0.1*micrometer);
264
265 G4Gamma::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
266 G4Electron::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
267 G4Positron::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
268
269 G4cout<<"the cut in energy for gamma in: "<<
270 (*theMaterialTable)[materialId]->GetName()
271 <<" is: "<<(G4Gamma::GetCutsInEnergy()[materialId])/keV
272 <<" keV" <<G4endl;
273 G4cout<<"the cut in energy for e- in: "<<
274 (*theMaterialTable)[materialId]->GetName()
275 <<" is: "<<(G4Electron::GetCutsInEnergy()[materialId])/keV
276 <<" keV" <<G4endl;
277
278 // Processes
279
280
281 G4int processType;
282 G4cout << "LowEnergy [1] or Standard [2] Compton or Standard PolarizedCompton[3] or LowEnergyPolarizedCompton [4]" << G4endl;
283 G4cin >> processType;
284 if ( !(processType == 1 || processType == 2
285 || processType == 3 || processType == 4))
286 {
287 G4Exception("Wrong input");
288 }
289
290 G4VDiscreteProcess* gammaProcess;
291
292
293 if (processType == 1)
294 {
295 gammaProcess = new G4LowEnergyCompton();
296 }
297 else if (processType == 2)
298 {
299 gammaProcess = new G4ComptonScattering();
300 }
301 else if (processType == 3)
302 {
303 gammaProcess = new G4PolarizedComptonScattering();
304 }
305 else if (processType == 4)
306 {
307 gammaProcess = new G4LowEnergyPolarizedCompton();
308 }
309
310 G4VProcess* theeminusMultipleScattering = new G4MultipleScattering();
311 G4VProcess* theeminusIonisation = new G4eIonisation();
312 G4VProcess* theeminusBremsstrahlung = new G4eBremsstrahlung();
313 G4VProcess* theeplusMultipleScattering = new G4MultipleScattering();
314 G4VProcess* theeplusIonisation = new G4eIonisation();
315 G4VProcess* theeplusBremsstrahlung = new G4eBremsstrahlung();
316 G4VProcess* theeplusAnnihilation = new G4eplusAnnihilation();
317
318 //----------------
319 // process manager
320 //----------------
321
322 // gamma
323
324 G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
325 gamma->SetProcessManager(gProcessManager);
326 gProcessManager->AddDiscreteProcess(gammaProcess);
327 G4ForceCondition* condition;
328
329 //electron
330
331 G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
332 electron->SetProcessManager(eProcessManager);
333 eProcessManager->AddProcess(theeminusMultipleScattering);
334 eProcessManager->AddProcess(theeminusIonisation);
335 eProcessManager->AddProcess(theeminusBremsstrahlung);
336
337 //positron
338
339 G4ProcessManager* pProcessManager = new G4ProcessManager(positron);
340 positron->SetProcessManager(pProcessManager);
341 pProcessManager->AddProcess(theeplusMultipleScattering);
342 pProcessManager->AddProcess(theeplusIonisation);
343 pProcessManager->AddProcess(theeplusBremsstrahlung);
344 pProcessManager->AddProcess(theeplusAnnihilation);
345
346 //--------------
347 // set ordering
348 //--------------
349
350
351 eProcessManager->
352 SetProcessOrdering(theeminusMultipleScattering, idxAlongStep,1);
353 eProcessManager->
354 SetProcessOrdering(theeminusIonisation, idxAlongStep,2);
355
356 eProcessManager->
357 SetProcessOrdering(theeminusMultipleScattering, idxPostStep,1);
358 eProcessManager->
359 SetProcessOrdering(theeminusIonisation, idxPostStep,2);
360 eProcessManager->
361 SetProcessOrdering(theeminusBremsstrahlung, idxPostStep,3);
362
363 pProcessManager->SetProcessOrderingToFirst(theeplusAnnihilation, idxAtRest);
364 pProcessManager->
365 SetProcessOrdering(theeplusMultipleScattering, idxAlongStep,1);
366 pProcessManager->
367 SetProcessOrdering(theeplusIonisation, idxAlongStep,2);
368
369 pProcessManager->
370 SetProcessOrdering(theeplusMultipleScattering, idxPostStep,1);
371 pProcessManager->
372 SetProcessOrdering(theeplusIonisation, idxPostStep,2);
373 pProcessManager->
374 SetProcessOrdering(theeplusBremsstrahlung, idxPostStep,3);
375 pProcessManager->
376 SetProcessOrdering(theeplusAnnihilation, idxPostStep,4);
377
378 // Create a DynamicParticle
379
380 // G4double eEnergy = initEnergy*keV;
381 G4double eEnergy = initEnergy;
382
383 // G4cout << eEnergy/keV << " INIT ENERGY (keV)" << G4endl;
384
385 G4ParticleMomentum eDirection(initX,initY,initZ);
386 G4DynamicParticle dynamicGamma(G4Gamma::Gamma(),eDirection,eEnergy);
387
388 G4cout << eDirection << " Direction" << G4endl;
389
390
391 // if (processType == 3 || processType == 4)
392 // {
393 G4double PolX, PolY, PolZ;
394 G4cout << "Polarization Vector" << G4endl;
395 G4cin >> PolX >> PolY >> PolZ;
396 dynamicGamma.SetPolarization(PolX, PolY, PolZ);
397 //G4cout << "polarization" << dynamicGamma.GetPolarization() << G4endl;
398 // }
399
400
401 dynamicGamma.DumpInfo();
402
403 // Track
404
405 G4ThreeVector aPosition(0.,0.,0.);
406 G4double aTime = 0. ;
407
408 G4Track* gTrack = new G4Track(&dynamicGamma,aTime,aPosition);
409
410 G4GRSVolume* touche = new G4GRSVolume(physicalFrame, NULL, aPosition);
411 gTrack->SetTouchable(touche);
412
413 // Step
414
415 G4Step* step = new G4Step();
416 step->SetTrack(gTrack);
417
418 G4StepPoint* aPoint = new G4StepPoint();
419 aPoint->SetPosition(aPosition);
420 aPoint->SetMaterial(material);
421 G4double safety = 10000.*cm;
422 aPoint->SetSafety(safety);
423 step->SetPreStepPoint(aPoint);
424
425 // Check applicability
426
427 if (! (gammaProcess->IsApplicable(*gamma)))
428 {
429 G4Exception("Not Applicable");
430 }
431 else
432 {
433 G4cout<< "applicability OK" << G4endl;
434 }
435
436 // Initialize the physics tables (in which material?)
437
438 gammaProcess->BuildPhysicsTable(*gamma);
439
440 theeminusMultipleScattering->BuildPhysicsTable(*electron);
441 theeminusIonisation->BuildPhysicsTable(*electron);
442 theeminusBremsstrahlung->BuildPhysicsTable(*electron);
443 theeplusMultipleScattering->BuildPhysicsTable(*positron);
444 theeplusIonisation->BuildPhysicsTable(*positron);
445 theeplusBremsstrahlung->BuildPhysicsTable(*positron);
446 theeplusAnnihilation->BuildPhysicsTable(*positron) ;
447
448 // G4cout<< "table OK" << endl;
449
450 /*
451
452 // Test GetMeanFreePath()
453
454 G4Material* apttoMaterial ;
455 G4String MaterialName ;
456
457 G4double minArg = 100*eV,maxArg = 100*GeV, argStp;
458 const G4int pntNum = 300;
459 G4double Tkin[pntNum+1];
460 G4double meanFreePath=0. ;
461
462 argStp = (std::log10(maxArg)-std::log10(minArg))/pntNum;
463
464 for(G4int d = 0; d < pntNum+1; d++)
465 {
466 Tkin[d] = std::pow(10,(std::log10(minArg) + d*argStp));
467 }
468
469 G4double sti = 1.*mm;
470 step->SetStepLength(sti);
471
472 // for ( G4int J = 0 ; J < nMaterials ; J++ )
473 // {
474 apttoMaterial = (*theMaterialTable)[materialId] ;
475 MaterialName = apttoMaterial->GetName() ;
476 logicalFrame->SetMaterial(apttoMaterial);
477
478 gTrack->SetStep(step);
479
480 G4LowEnergyCompton* gammaLowEProcess =
481 (G4LowEnergyCompton*) gammaProcess;
482 G4ComptonScattering* gammaStdProcess =
483 (G4ComptonScattering*) gammaProcess;
484
485
486 for (G4int i=0 ; i<pntNum; i++)
487 {
488 dynamicGamma.SetKineticEnergy(Tkin[i]);
489 if (processType == 1)
490 {
491 meanFreePath=gammaLowEProcess
492 ->GetMeanFreePath(*gTrack, sti, condition);
493 }
494 else
495 {
496 meanFreePath=gammaStdProcess
497 ->GetMeanFreePath(*gTrack, sti, condition);
498 }
499
500 ntuple3->column("kinen",Tkin[i]);
501 ntuple3->column("mfp",meanFreePath/cm);
502 ntuple3->dumpData();
503
504 // G4cout << meanFreePath/cm << G4endl;
505
506 }
507 G4cout << "Mean Free Path OK" << G4endl;
508 */
509
510 // --------- Test the DoIt
511
512 G4cout << "DoIt in " << material->GetName() << G4endl;
513
514 dynamicGamma.SetKineticEnergy(eEnergy);
515 dynamicGamma.SetMomentumDirection(initX,initY,initZ);
516
517 for (G4int iter=0; iter<nIterations; iter++)
518 {
519
520 step->SetStepLength(1*micrometer);
521
522
523
524 G4cout << "Iteration = " << iter
525 << " - Step Length = "
526 << step->GetStepLength()/mm << " mm "
527 << G4endl;
528
529 gTrack->SetStep(step);
530
531 G4StepPoint* preStep = step->GetPreStepPoint();
532 G4StepPoint* postStep = step->GetPostStepPoint();
533 G4ThreeVector prePosition = preStep->GetPosition();
534 G4ThreeVector postPosition = postStep->GetPosition();
535
536 //G4cout << prePosition << "pre step point "<< G4endl;
537 //G4cout << postPosition << "post step point "<< G4endl;
538
539 G4ThreeVector polInitial=dynamicGamma.GetPolarization();
540
541 G4cout << polInitial << " Initial Polarization" << G4endl;
542
543 G4VParticleChange* dummy;
544 dummy = gammaProcess->PostStepDoIt(*gTrack, *step);
545 G4ParticleChange* particleChange = (G4ParticleChange*) dummy;
546
547
548 // Primary physical quantities
549
550
551 // particleChange->DumpInfo();
552
553 G4double energyChange = particleChange->GetEnergyChange();
554
555 G4double dedx = initEnergy - energyChange ;
556
557 G4double dedxNow = dedx / (step->GetStepLength());
558
559 G4ThreeVector eChange =
560 particleChange->CalcMomentum(energyChange,
561 (*particleChange->GetMomentumChange()),
562 particleChange->GetMassChange());
563
564 G4double pxChange = eChange.x();
565 G4double pyChange = eChange.y();
566 G4double pzChange = eChange.z();
567 G4double pChange =
568 std::sqrt(pxChange*pxChange + pyChange*pyChange + pzChange*pzChange);
569 G4double thetaChange = eChange.theta();
570
571
572 const G4ThreeVector* momChange =particleChange->GetMomentumDirectionChange();
573
574 G4cout << (momChange->x()) << " " << (momChange->y()) << " " << (momChange->z()) << " " << G4endl;
575
576 G4cout << eChange << "newdir" << G4endl;
577 G4double phiChange = eChange.phi();
578
579 G4double xChange = particleChange->GetPositionChange()->x();
580 G4double yChange = particleChange->GetPositionChange()->y();
581 G4double zChange = particleChange->GetPositionChange()->z();
582
583 //G4cout << "Theta " << thetaChange << G4endl;
584 //G4cout << "Phi " << phiChange << G4endl;
585
586 G4cout << "---- Primary after the step ---- " << G4endl;
587
588 G4cout << "Position (x,y,z) = "
589 << xChange << " "
590 << yChange << " "
591 << zChange << " "
592 << G4endl;
593
594 G4cout << " Initial Energy " << initEnergy/keV << " keV" << G4endl;
595 G4cout << "---- Energy: " << energyChange/MeV << " MeV, "
596 << "(px,py,pz): ("
597 << pxChange/keV << ","
598 << pyChange/keV << ","
599 << pzChange/keV << ") keV"
600 << G4endl;
601 /* G4cout << "---- Energy loss (dE) = " << dedx/keV << " keV" << G4endl;
602 G4cout << "Stopping power (dE/dx)=" << dedxNow << G4endl;
603 */
604
605
606 G4double electronMass = 511.22*keV; // da inserire la definizione
607
608 G4double Ratio = energyChange/
609 (initEnergy/(1 + (initEnergy*(1-std::cos(thetaChange))/electronMass)));
610 // testenergy
611
612 //G4cout << Ratio << "RATIO" << G4endl;
613 //G4cout << energyChange/keV << "ENERGY (keV)" << G4endl;
614
615
616 const G4ThreeVector* polChange=particleChange->GetPolarizationChange();
617
618
619 //G4cout << pxChange/pChange << "X" << G4endl;
620 //G4cout << pyChange/pChange << "Y" << G4endl;
621 //G4cout << pzChange/pChange << "Z" << G4endl;
622
623 //G4cout << polChange->x() << "pol X" << G4endl;
624 //G4cout << polChange->y() << "pol Y" << G4endl;
625 //G4cout << polChange->z() << "pol Z" << G4endl;
626 //G4cout << polChange->mag() << "pol mag" << G4endl;
627
628
629 G4double ScalarProduct = (polChange->x())*(pxChange/pChange)+
630 (pyChange/pChange)*(polChange->y())+
631 (pzChange/pChange)*(polChange->z());
632
633 //G4cout << ScalarProduct << "scalar product" << G4endl;
634
635 hgETot->accumulate(energyChange/keV);
636 hgP->accumulate(pChange/keV);
637 hgTheta->accumulate(thetaChange);
638 hgPhi->accumulate(phiChange);
639 hgRapp->accumulate(Ratio);
640
641 // Secondaries
642
643 ntuple1->column("eprimary", initEnergy/keV);
644 ntuple1->column("energyf", energyChange/keV);
645 ntuple1->column("de", dedx/keV);
646 ntuple1->column("dedx", dedxNow/keV);
647 ntuple1->column("pxch", pxChange);
648 ntuple1->column("pych", pyChange);
649 ntuple1->column("pzch", pzChange);
650 ntuple1->column("pch", pChange);
651 ntuple1->column("polx",(polInitial.x()));
652 ntuple1->column("poly",(polInitial.y()));
653 ntuple1->column("polz",(polInitial.z()));
654 ntuple1->column("polchx",(polChange->x()));
655 ntuple1->column("polchy",(polChange->y()));
656 ntuple1->column("polchz",(polChange->z()));
657 ntuple1->column("thetach", thetaChange);
658 ntuple1->column("phich", phiChange);
659 ntuple1->dumpData();
660
661 // Secondaries physical quantities
662
663 hNSec->accumulate(particleChange->GetNumberOfSecondaries());
664 hDebug->accumulate(particleChange->GetLocalEnergyDeposit());
665
666 G4cout << " secondaries " <<
667 particleChange->GetNumberOfSecondaries() << G4endl;
668
669 G4double Etotal = 0.;
670 Etotal += energyChange;
671
672 //G4cout << " Total energy" << Etotal << G4endl;
673
674 for (G4int i = 0; i < (particleChange->GetNumberOfSecondaries()); i++)
675 {
676 // The following two items should be filled per event, not
677 // per secondary; filled here just for convenience, to avoid
678 // complicated logic to dump ntuple when there are no secondaries
679
680 G4Track* finalParticle = particleChange->GetSecondary(i) ;
681
682 G4double e = finalParticle->GetTotalEnergy();
683 G4double eKin = finalParticle->GetKineticEnergy();
684 G4double px = (finalParticle->GetMomentum()).x();
685 G4double py = (finalParticle->GetMomentum()).y();
686 G4double pz = (finalParticle->GetMomentum()).z();
687 G4double theta = (finalParticle->GetMomentum()).theta();
688 G4double p = std::sqrt(px*px+py*py+pz*pz);
689
690 if (eKin > initEnergy)
691 {
692 G4cout << "WARNING: eKinFinal > eKinInit " << G4endl;
693 // << e
694 // << " > " initEnergy
695
696 }
697
698 G4String particleName =
699 finalParticle->GetDefinition()->GetParticleName();
700 G4cout << "==== Final "
701 << particleName << " "
702 << "energy: " << e/keV << " keV, "
703 << "eKin: " << eKin/keV << " keV, "
704 << "(px,py,pz): ("
705 << px/keV << ","
706 << py/keV << ","
707 << pz/keV << ") keV "
708 << G4endl;
709
710 // G4cout << " energia secondaria" << e << G4endl;
711
712 heETot->accumulate(eKin/keV);
713 heP->accumulate(p/keV);
714
715 Etotal += eKin;
716 //G4cout << " energia totale" << Etotal << G4endl;
717
718 G4int partType;
719 if (particleName == "e-") partType = 1;
720 else if (particleName == "e+") partType = 2;
721 else if (particleName == "gamma") partType = 3;
722
723 // Fill the secondaries ntuple
724
725 ntuple2->column("event",iter);
726 ntuple2->column("eprimary",initEnergy/keV);
727 ntuple2->column("px", px);
728 ntuple2->column("py", py);
729 ntuple2->column("pz", pz);
730 ntuple2->column("p", p);
731 ntuple2->column("e", e/keV);
732 ntuple2->column("theta", theta);
733 ntuple2->column("ekin", eKin/keV);
734 ntuple2->column("type", partType);
735
736 ntuple2->dumpData();
737
738 delete particleChange->GetSecondary(i);
739 }
740
741 // G4cout << Etotal/keV << " E total (keV) " << G4endl;
742 hSumE->accumulate(Etotal/keV);
743 particleChange->Clear();
744
745 }
746
747
748 // G4cout << "Iteration number: " << iter << G4endl;
749 hbookManager->write();
750 delete hbookManager;
751
752 delete step;
753
754 G4cout << "END OF THE MAIN PROGRAM" << G4endl;
755}
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
Note: See TracBrowser for help on using the repository browser.