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

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

nvx fichiers dans CVS
File size: 21.6 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// $Id: G4ComptonTest.cc,v 1.28 2008/04/24 14:14:25 pia Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29// -------------------------------------------------------------------
30// GEANT 4 class file --- Copyright CERN 1998
31// CERN Geneva Switzerland
32//
33//
34// File name: G4ComptonTest
35//
36// Author: Maria Grazia Pia
37// Andreas Pfeiffer
38//
39// Creation date: 2 May 2001
40//
41// Modifications:
42// 08 Mar 2008 MGP Updated to recent Geant4 interface changes
43// 28 Nov 2002 AP update to AIDA 3
44// 14 Sep 2001 AP Moved histograms to Lizard
45// 16 Sep 2001 AP Moved ntuples to Lizard
46// 24 Apr 2008 MGP Upgrate to treat couples correctly from Luciano Pandola's PenelopeComptonTest
47//
48// -------------------------------------------------------------------
49//
50// (MGP) The following is obsolete and should be replaced by iAIDA instructions
51//
52// from: geant4/source/processes/electromagnetic/lowenergy/test/
53//
54// execute the following lines _before_ gmake,
55// source /afs/cern.ch/sw/lhcxx/share/LHCXX/latest/scripts/setupAnaphe
56//
57// or, for [t]csh fans:
58//
59// source /afs/cern.ch/sw/lhcxx/share/LHCXX/latest/scripts/setupAnaphe.csh
60//
61// both assume that you have the correct PATH to the compiler
62//
63// [gmake and run your simulation]
64//
65// to start Lizard:
66// /afs/cern.ch/sw/lhcxx/share/LHCXX/latest/scripts/lizard
67//
68// see also: http://cern.ch/Anaphe
69//
70// ********************************************************************
71
72#include "globals.hh"
73#include "G4ios.hh"
74#include <fstream>
75#include <iomanip>
76
77#include "G4Material.hh"
78#include "G4VContinuousDiscreteProcess.hh"
79#include "G4ProcessManager.hh"
80#include "G4LowEnergyBremsstrahlung.hh"
81#include "G4eBremsstrahlung.hh"
82#include "G4LowEnergyCompton.hh"
83#include "G4PenelopeCompton.hh"
84#include "G4LowEnergyPolarizedCompton.hh"
85#include "G4ComptonScattering.hh"
86#include "G4LowEnergyIonisation.hh"
87#include "G4eIonisation.hh"
88#include "G4VeLowEnergyLoss.hh"
89#include "G4EnergyLossTables.hh"
90#include "G4VParticleChange.hh"
91#include "G4ParticleChange.hh"
92#include "G4DynamicParticle.hh"
93#include "G4Electron.hh"
94#include "G4Positron.hh"
95#include "G4Gamma.hh"
96#include "G4Box.hh"
97#include "G4PVPlacement.hh"
98#include "G4Step.hh"
99#include "G4GRSVolume.hh"
100#include "G4UnitsTable.hh"
101#include "G4ProductionCutsTable.hh"
102#include "G4MaterialCutsCouple.hh"
103#include "G4RunManager.hh"
104#include "G4RegionStore.hh"
105#include "G4StateManager.hh"
106#include "G4ApplicationState.hh"
107
108// New Histogramming (from AIDA and Anaphe):
109#include <memory> // for the auto_ptr(T>
110
111#include "AIDA/AIDA.h"
112
113int main()
114{
115 // G4String fileName;
116 // G4cout << "Enter histogram file name" << G4endl;
117 // G4cin >> fileName;
118
119
120 // Setup
121
122 // G4cout.setf( ios::scientific, ios::floatfield );
123
124 // -------------------------------------------------------------------
125
126 // ---- HBOOK initialization
127
128 // Creating the analysis factory
129 std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() );
130
131 // Creating the tree factory
132 std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() );
133
134 // Creating a tree mapped to a new hbook file.
135 bool readOnly = false;
136 bool createFile = true;
137 std::auto_ptr< AIDA::ITree > tree( tf->create( "comptonhisto.hbook", "hbook", readOnly, createFile ) );
138 std::cout << "Tree store : " << tree->storeName() << std::endl;
139
140 // Next create the nTuples using the factory and open it for writing
141 // Creating a tuple factory, whose tuples will be handled by the tree
142 std::auto_ptr< AIDA::ITupleFactory > tpf( af->createTupleFactory( *tree ) );
143
144
145 // ---- Primary ntuple ------
146 // If using Anaphe HBOOK implementation, there is a limitation on the length of the
147 // variable names in a ntuple
148 AIDA::ITuple* ntuple1 = tpf->create( "1", "Primary tuple",
149 "float e0, e1, dedx, dedxNow, px1, py1, pz1, p1, theta1, neminus, neplus, nphoton" );
150
151 // ---- Secondary ntuple ------
152 AIDA::ITuple* ntuple2 = tpf->create( "2", "Secondary tuple",
153 "float px, py, pz, p, e, ekin, theta, phi, type" );
154
155 // ---- Secondaries histos ----
156 // Creating a histogram factory, whose histograms will be handled by the tree
157 std::auto_ptr< AIDA::IHistogramFactory > hf( af->createHistogramFactory( *tree ) );
158
159 // Creating an 1-dimensional histogram in the root directory of the tree
160
161 AIDA::IHistogram1D* hEKin;
162 hEKin = hf->createHistogram1D("10","Kinetic Energy", 100,0.,10.);
163
164 AIDA::IHistogram1D* hP;
165 hP = hf->createHistogram1D("20","Momentum", 100,0.,10.);
166
167 AIDA::IHistogram1D* hNSec;
168 hNSec = hf->createHistogram1D("30","Number of secondaries", 10,0.,10.);
169
170 AIDA::IHistogram1D* hDeposit;
171 hDeposit = hf->createHistogram1D("40","Local energy deposit", 100,0.,10.);
172
173 AIDA::IHistogram1D* hTheta;
174 hTheta = hf->createHistogram1D("50","Theta", 100,0.,pi);
175
176 AIDA::IHistogram1D* hPhi;
177 hPhi = hf->createHistogram1D("60","Phi", 100,-pi,pi);
178
179 AIDA::IHistogram1D* hE1;
180 hE1 = hf->createHistogram1D("70","Scattered particle energy", 100,0.,10.);
181
182 AIDA::IHistogram1D* hEdiff;
183 hEdiff = hf->createHistogram1D("80","Energy difference initial-scattered particle", 100,0.,10.);
184
185
186 // end NEW
187 // ================================================================================
188
189 // ==================== end of Histogram and NTuple handling
190
191 //--------- Materials definition ---------
192
193 G4Material* Be = new G4Material("Beryllium", 4., 9.01*g/mole, 1.848*g/cm3);
194 G4Material* Graphite = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
195 G4Material* Al = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3);
196 G4Material* Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3);
197 G4Material* LAr = new G4Material("LArgon", 18., 39.95*g/mole, 1.393*g/cm3);
198 G4Material* Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3);
199 G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
200 G4Material* W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
201 G4Material* Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3);
202 G4Material* U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3);
203
204 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
205 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
206 G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
207 G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
208 G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
209 G4Element* N = new G4Element("Nitrogen", "N" , 7., 14.01*g/mole);
210
211 G4Material* maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3);
212
213 G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
214 water->AddElement(H,2);
215 water->AddElement(O,1);
216
217 G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
218 ethane->AddElement(H,6);
219 ethane->AddElement(C,2);
220
221 G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
222 csi->AddElement(Cs,1);
223 csi->AddElement(I,1);
224
225 G4Material* air = new G4Material("Air" , 1.290*mg/cm3, 2);
226 air->AddElement(N,0.7);
227 air->AddElement(O,0.3);
228
229
230 // Interactive set-up
231
232 G4cout << "How many interactions? " << G4endl;
233 G4int nIterations;
234 G4cin >> nIterations;
235 if (nIterations <= 0) G4Exception("Wrong input");
236
237 G4cout << "Enter the initial particle energy E (MeV)" << G4endl;
238 G4double initEnergy;
239 G4cin >> initEnergy ;
240 initEnergy = initEnergy * MeV;
241 G4double initialEnergy = initEnergy;
242 if (initEnergy <= 0.) G4Exception("Wrong input");
243
244
245
246 // Dump the material table
247 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
248 G4int nMaterials = G4Material::GetNumberOfMaterials();
249 G4cout << "Available materials are: " << G4endl;
250 for (G4int mat = 0; mat < nMaterials; mat++)
251 {
252 G4cout << mat << ". "
253 << (*theMaterialTable)[mat]->GetName()
254 << G4endl;
255 }
256
257 G4cout << "Which material? " << G4endl;
258 G4int materialId;
259 G4cin >> materialId;
260
261 G4Material* material = (*theMaterialTable)[materialId] ;
262
263 G4cout << "The selected material is: "
264 << material->GetName()
265 << G4endl;
266 // Geometry
267
268 G4double dimX = 1 * mm;
269 G4double dimY = 1 * mm;
270 G4double dimZ = 1 * mm;
271
272 G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
273 G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,material,
274 "LFrame", 0, 0, 0);
275 logicalFrame->SetMaterial(material);
276 G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
277 "PFrame",logicalFrame,0,false,0);
278
279 // RunManager
280 G4RunManager* rm = new G4RunManager();
281
282 rm->GeometryHasBeenModified();
283 rm->DefineWorldVolume(physicalFrame);
284
285 G4cout << "[OK] World is defined " << G4endl;
286 G4StateManager::GetStateManager()->SetNewState(G4State_Idle);
287 rm->DumpRegion("DefaultRegionForTheWorld"); //this forces the region update!
288
289 // Particle definitions
290
291 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
292 G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
293 G4ParticleDefinition* positron = G4Positron::PositronDefinition();
294
295 G4ProductionCutsTable* cutsTable = G4ProductionCutsTable::GetProductionCutsTable();
296 G4ProductionCuts* cuts = cutsTable->GetDefaultProductionCuts();
297 if (cuts == 0) G4cout << " G4ProductionCuts* cuts = 0" << G4endl;
298
299 G4double cutG = 1*micrometer;
300 G4double cutE = 1*micrometer;
301 cuts->SetProductionCut(cutG, gamma); // photons
302 cuts->SetProductionCut(cutE, electron); // electrons
303 cuts->SetProductionCut(cutE, positron); // positrons
304 cuts->SetProductionCut(cutG, 0); //gammas
305 cuts->SetProductionCut(cutE, 1); //electrons
306 cuts->SetProductionCut(cutE, 2); //positrons
307 //G4double cutAll = 1.*micrometer;
308 //cuts->SetProductionCut(cutAll, -1); //all
309
310 G4cout << "Cuts are defined " << G4endl;
311
312 // MGP 8/3/2008 - There is something wrong with the Couple, to be investigated
313 G4MaterialCutsCouple* couple = new G4MaterialCutsCouple(material,cuts);
314
315 logicalFrame->SetMaterialCutsCouple(couple);
316
317 G4cout << "Recalculation needed: " << couple->IsRecalcNeeded() << G4endl;
318
319 cutsTable->UpdateCoupleTable(physicalFrame);
320 cutsTable->PhysicsTableUpdated();
321 cutsTable->DumpCouples();
322
323 //couple->SetUseFlag(true);
324 //cutsTable->UpdateCoupleTable(world);
325 //cutsTable->DumpCouples();
326
327 //RunManager
328 //G4RunManager* rm = new G4RunManager();
329 //rm->GeometryHasBeenModified();
330 //G4VPhysicalVolume* world(physicalFrame);
331 //rm->DefineWorldVolume(world);
332 //G4cout << "[OK] World is defined " << G4endl;
333
334
335 // Processes
336
337 G4int processType;
338 G4cout
339 << "LowEnergy [1] or Penelope [2] or LowEnergyPolarized [3] or Standard [4]?"
340 << G4endl;
341 G4cin >> processType;
342 if (processType < 1 || processType > 4 ) G4Exception("Wrong input");
343
344 G4VContinuousDiscreteProcess* bremProcess;
345 G4VContinuousDiscreteProcess* ioniProcess;
346
347 if (processType < 4)
348 {
349 bremProcess = new G4LowEnergyBremsstrahlung;
350 ioniProcess = new G4LowEnergyIonisation;
351 }
352 else
353 {
354 bremProcess = new G4eBremsstrahlung;
355 ioniProcess = new G4eIonisation;
356 }
357
358
359 G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
360 electron->SetProcessManager(eProcessManager);
361 eProcessManager->AddProcess(bremProcess);
362 eProcessManager->AddProcess(ioniProcess);
363
364 G4ProcessManager* positronProcessManager = new G4ProcessManager(positron);
365 positron->SetProcessManager(positronProcessManager);
366 positronProcessManager->AddProcess(new G4eBremsstrahlung);
367 positronProcessManager->AddProcess( new G4eIonisation());
368
369 // Initialize the physics tables
370 bremProcess->BuildPhysicsTable(*electron);
371 ioniProcess->BuildPhysicsTable(*electron);
372
373 // Photon process
374 G4VDiscreteProcess* photonProcess = 0;
375 if (processType == 1)
376 {
377 photonProcess = new G4LowEnergyCompton;
378 G4cout << "G4LowEnergyCompton CREATED" << G4endl;
379 }
380 if (processType == 2)
381 {
382 photonProcess = new G4PenelopeCompton;
383 }
384 if (processType == 3)
385 {
386 photonProcess = new G4LowEnergyPolarizedCompton;
387 }
388 if (processType == 4)
389 {
390 photonProcess = new G4ComptonScattering;
391 }
392
393 G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
394 gamma->SetProcessManager(gProcessManager);
395 gProcessManager->AddProcess(photonProcess);
396 photonProcess->BuildPhysicsTable(*gamma);
397
398 // Primary direction
399 G4double initX = 0. * mm;
400 G4double initY = 0. * mm;
401 G4double initZ = 1. * mm;
402 G4ParticleMomentum gDirection(initX,initY,initZ);
403
404 // Check applicability
405
406 if (! (photonProcess->IsApplicable(*gamma))) G4Exception("Not Applicable");
407
408 // Initialize the physics tables (in which material?)
409 photonProcess->BuildPhysicsTable(*gamma);
410
411
412 // --------- Test the DoIt -------------------------------------------------------------------
413
414 G4cout << "DoIt in material " << material->GetName() << G4endl;
415
416
417 G4Track* gTrack = 0;
418 G4GRSVolume* touche = 0;
419 G4Step* step = 0;
420 G4StepPoint* point = 0;
421 G4StepPoint* newPoint = 0;
422 G4double gEnergy = 0.;
423
424 for (G4int iter=0; iter<nIterations; iter++)
425 {
426 // Primary energy
427 // gEnergy = initEnergy*MeV*G4UniformRand();
428 gEnergy = initEnergy*MeV;
429 G4cout << "---- Initial energy = " << gEnergy/MeV << " MeV" << G4endl;
430
431 // Dynamic particle (incident primary)
432 G4DynamicParticle dynamicPhoton(G4Gamma::Gamma(),gDirection,gEnergy);
433
434 // Track (incident)
435 G4ThreeVector position(0.,0.,0.);
436 G4double time = 0. ;
437 gTrack = new G4Track(&dynamicPhoton,time,position);
438
439 // Do I really need this?
440 touche = new G4GRSVolume(physicalFrame, 0, position);
441 //gTrack->SetTouchable(touche);
442
443 // Step
444 step = new G4Step();
445 step->SetTrack(gTrack);
446 //const G4MaterialCutsCouple* theCouple(cutsTable->GetMaterialCutsCouple(material, cutsTable->GetDefaultProductionCuts()));
447 // if (couple == 0) G4cout << "Couple = 0 in setting Step" <<G4endl;
448 //const G4MaterialCutsCouple* couple = new G4MaterialCutsCouple(material,cuts);
449
450 // PreStep point
451 point = new G4StepPoint();
452 point->SetPosition(position);
453 point->SetMaterial(material);
454 point->SetMaterialCutsCouple(couple);
455 G4double safety = 10000.*cm;
456 point->SetSafety(safety);
457 step->SetPreStepPoint(point);
458
459 // PostStep point
460 newPoint = new G4StepPoint();
461 G4ThreeVector newPosition(0.,0.,1.*mm);
462 newPoint->SetPosition(newPosition);
463 newPoint->SetMaterial(material);
464 newPoint->SetMaterialCutsCouple(couple);
465 newPoint->SetSafety(safety);
466 step->SetPostStepPoint(newPoint);
467
468 // Step length
469 step->SetStepLength(1.*micrometer);
470
471 gTrack->SetStep(step);
472
473 // const G4MaterialCutsCouple* couple = gTrack->GetMaterialCutsCouple();
474
475 G4cout << "Iteration = "
476 << iter
477 << " - Step Length = "
478 << step->GetStepLength()/mm
479 << " mm "
480 << G4endl;
481
482 G4ParticleChange* particleChange = (G4ParticleChange*) photonProcess->PostStepDoIt(*gTrack,*step);
483
484 // Primary physical quantities
485
486 G4double energyChange = particleChange->GetEnergy();
487 G4double dedx = gEnergy - energyChange ;
488 G4double dedxNow = dedx / (step->GetStepLength());
489
490 G4ThreeVector eChange = particleChange->CalcMomentum(energyChange,
491 (*particleChange->GetMomentumDirection()),
492 particleChange->GetMass());
493 G4double pxChange = eChange.x();
494 G4double pyChange = eChange.y();
495 G4double pzChange = eChange.z();
496 G4double pChange = std::sqrt(pxChange*pxChange + pyChange*pyChange + pzChange*pzChange);
497
498 G4double xChange = particleChange->GetPosition()->x();
499 G4double yChange = particleChange->GetPosition()->y();
500 G4double zChange = particleChange->GetPosition()->z();
501 G4double thetaChange = particleChange->GetMomentumDirection()->theta();
502
503 G4cout << "---- Energy: "
504 << energyChange/MeV << " MeV, "
505 << "(px,py,pz): ("
506 << pxChange/MeV << ","
507 << pyChange/MeV << ","
508 << pzChange/MeV << ") MeV"
509 << G4endl;
510
511 G4cout << "---- Energy loss (dE) = " << dedx/keV << " keV" << G4endl;
512
513 // Primary
514
515 hNSec->fill(particleChange->GetNumberOfSecondaries());
516 hDeposit->fill(particleChange->GetLocalEnergyDeposit());
517 hE1->fill(energyChange);
518 hEdiff->fill(dedx);
519
520 G4int nElectrons = 0;
521 G4int nPositrons = 0;
522 G4int nPhotons = 0;
523
524
525
526 // Secondaries
527
528 for (G4int i = 0; i < (particleChange->GetNumberOfSecondaries()); i++)
529 {
530 G4Track* finalParticle = particleChange->GetSecondary(i) ;
531
532 G4double e = finalParticle->GetTotalEnergy();
533 G4double eKin = finalParticle->GetKineticEnergy();
534 G4double px = (finalParticle->GetMomentum()).x();
535 G4double py = (finalParticle->GetMomentum()).y();
536 G4double pz = (finalParticle->GetMomentum()).z();
537 G4double theta = (finalParticle->GetMomentum()).theta();
538 G4double phi = (finalParticle->GetMomentum()).phi();
539 G4double p = std::sqrt(px*px+py*py+pz*pz);
540
541 if (eKin > gEnergy)
542 {
543 G4cout << "WARNING: eFinal > eInit " << G4endl;
544 }
545
546
547 G4String particleName = finalParticle->GetDefinition()->GetParticleName();
548 G4ParticleDefinition* def = finalParticle->GetDefinition();
549 G4cout << "==== Final "
550 << particleName << " "
551 << "energy: " << e/MeV << " MeV, "
552 << "eKin: " << eKin/MeV << " MeV, "
553 << "(px,py,pz): ("
554 << px/MeV << ","
555 << py/MeV << ","
556 << pz/MeV << ") MeV "
557 << G4endl;
558
559 hEKin->fill(eKin);
560 hP->fill(p);
561 hTheta->fill(theta);
562 hPhi->fill(phi);
563
564 G4double particleType = -1.;
565 if (def == electron)
566 {
567 particleType = 1.;
568 nElectrons++;
569 }
570 if (def == positron)
571 {
572 particleType = 2.;
573 nPositrons++;
574 }
575
576 if (def == gamma)
577 {
578 particleType = 3.;
579 nPhotons++;
580 }
581
582 // Fill the secondaries ntuple
583 ntuple2->fill( ntuple2->findColumn( "px" ), px );
584 ntuple2->fill( ntuple2->findColumn( "py" ), py );
585 ntuple2->fill( ntuple2->findColumn( "pz" ), pz );
586 ntuple2->fill( ntuple2->findColumn( "p" ), p );
587 ntuple2->fill( ntuple2->findColumn( "e" ), e );
588 ntuple2->fill( ntuple2->findColumn( "ekin" ), eKin );
589 ntuple2->fill( ntuple2->findColumn( "theta" ), theta );
590 ntuple2->fill( ntuple2->findColumn( "phi" ), phi );
591 ntuple2->fill( ntuple2->findColumn( "type" ), particleType );
592
593 // NEW: Values of attributes are prepared; store them to the nTuple:
594 ntuple2->addRow(); // check for returning true ...
595
596 delete particleChange->GetSecondary(i);
597 } // end loop over secondaries
598
599 // Fill the primaries ntuple
600
601 ntuple1->fill( ntuple1->findColumn( "e0" ), gEnergy );
602 ntuple1->fill( ntuple1->findColumn( "e1" ), energyChange );
603 ntuple1->fill( ntuple1->findColumn( "dedx" ), dedx );
604 ntuple1->fill( ntuple1->findColumn( "dedxNow" ), dedxNow );
605 ntuple1->fill( ntuple1->findColumn( "px1" ), pxChange );
606 ntuple1->fill( ntuple1->findColumn( "py1" ), pyChange );
607 ntuple1->fill( ntuple1->findColumn( "pz1" ), pzChange );
608 ntuple1->fill( ntuple1->findColumn( "p1" ), pChange );
609 ntuple1->fill( ntuple1->findColumn( "theta1" ), thetaChange );
610 ntuple1->fill( ntuple1->findColumn( "neminus" ), (G4double) nElectrons );
611 ntuple1->fill( ntuple1->findColumn( "neplus" ), (G4double) nPositrons );
612 ntuple1->fill( ntuple1->findColumn( "nphoton" ), (G4double) nPhotons );
613
614 //NEW: Values of attributes are prepared; store them to the nTuple:
615 ntuple1->addRow();
616
617 particleChange->Clear();
618
619 delete touche;
620 delete step;
621 delete gTrack;
622
623
624 } // end loop over events
625
626 G4cout << "-----------------------------------------------------" << G4endl;
627
628 // Committing the transaction with the tree
629 G4cout << "Committing..." << G4endl;
630 tree->commit();
631 G4cout << "Closing the tree..." << G4endl;
632 tree->close();
633
634 G4cout << "END OF TEST" << G4endl;
635
636
637}
Note: See TracBrowser for help on using the repository browser.