source: trunk/source/processes/electromagnetic/lowenergy/test/processTest/src/G4ProcessTestAnalysis.cc@ 1201

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

nvx fichiers dans CVS
File size: 10.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//
27// $Id: G4ProcessTestAnalysis.cc,v 1.11 2006/06/29 19:48:52 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// Author: A. Pfeiffer (Andreas.Pfeiffer@cern.ch)
31// (copy of his UserAnalyser class)
32//
33// History:
34// -----------
35// 5 Nov 2001 MGP Implementation of process test analysis
36//
37// -------------------------------------------------------------------
38
39#include "G4ProcessTestAnalysis.hh"
40#include "globals.hh"
41#include "G4ParticleChange.hh"
42#include "G4Track.hh"
43
44G4ProcessTestAnalysis* G4ProcessTestAnalysis::instance = 0;
45
46G4ProcessTestAnalysis::G4ProcessTestAnalysis()
47{
48 histoManager = createIHistoManager();
49 ntFactory = Lizard::createNTupleFactory();
50}
51
52G4ProcessTestAnalysis::~G4ProcessTestAnalysis()
53{
54 delete ntFactory;
55 ntFactory = 0;
56
57 // std::map< G4String,IHistogram1D*,std::less<G4String> >::iterator pos1;
58 // for (pos1 = histo1D.begin(); pos1 != histo1D.end(); ++pos1)
59 // {
60 // IHistogram* h = pos1->second;
61 // delete h;
62 // }
63
64 delete histoManager;
65 histoManager = 0;
66}
67
68void G4ProcessTestAnalysis::book(const G4String& storeName)
69{
70 const char* nameStore = (storeName + ".hbook").c_str();
71 histoManager->selectStore(nameStore);
72
73 // Book histograms
74 histoManager->create1D("1","Number of secondaries", 10,0.,10.);
75 histoManager->create1D("2","Local energy deposit", 100,0.,10.);
76 histoManager->create1D("3","Kinetic Energy", 100,0.,10.);
77 histoManager->create1D("4","Theta", 100,0.,pi);
78 histoManager->create1D("5","Phi", 100,-pi,pi);
79
80 /*
81 IHistogram1D* hNSec = histoManager->create1D("1","Number of secondaries", 10,0.,10.);
82 histo1D["nSec"] = hNSec;
83 IHistogram1D* hDeposit = histoManager->create1D("2","Local energy deposit", 100,0.,10.);
84 histo1D["eDeposit"] = hDeposit;
85 IHistogram1D* hEKin = histoManager->create1D("3","Kinetic Energy", 100,0.,10.);
86 histo1D["eKin"] = hEKin;
87 IHistogram1D* hTheta = histoManager->create1D("4","Theta", 100,0.,pi);
88 histo1D["theta"] = hTheta;
89 IHistogram1D* hPhi = histoManager->create1D("5","Phi", 100,-pi,pi);
90 histo1D["phi"] = hPhi;
91 */
92
93 // Book ntuples
94
95 G4String ntFileName = storeName + ".hbook::1";
96 const char* name = ntFileName.c_str();
97 ntuple1 = ntFactory->createC(name);
98
99 // Add and bind the attributes to the general final state nTuple
100
101 if( !( ntuple1->addAndBind( "eprimary" , initialEnergy) &&
102 ntuple1->addAndBind( "pxch" , pxChange ) &&
103 ntuple1->addAndBind( "pych" , pyChange ) &&
104 ntuple1->addAndBind( "pzch" , pzChange ) &&
105 ntuple1->addAndBind( "thetach" , thetaChange ) &&
106 ntuple1->addAndBind( "phich" , phiChange ) &&
107 ntuple1->addAndBind( "deposit" , eDeposit ) &&
108 ntuple1->addAndBind( "eminus" , nElectrons ) &&
109 ntuple1->addAndBind( "eplus" , nPositrons ) &&
110 ntuple1->addAndBind( "photons" , nPhotons ) ) )
111 {
112 G4cerr << "Error: unable to add attribute to nTuple1." << G4endl;
113 // Must be cleaned up properly before any exit.
114 delete ntuple1;
115 G4Exception("Could not addAndBind ntuple1");
116 }
117 // ntuples["primary"] = ntuple1;
118
119 ntFileName = storeName + ".hbook::2";
120 name = ntFileName.c_str();
121 ntuple2 = ntFactory->createC(name);
122
123 // Add and bind the attributes to the secondary nTuple
124 if ( !( ntuple2->addAndBind( "e0" , initialEnergy) &&
125 ntuple2->addAndBind( "px" , px ) &&
126 ntuple2->addAndBind( "py" , py ) &&
127 ntuple2->addAndBind( "pz" , pz ) &&
128 ntuple2->addAndBind( "p" , p ) &&
129 ntuple2->addAndBind( "e" , e ) &&
130 ntuple2->addAndBind( "ekin" , eKin ) &&
131 ntuple2->addAndBind( "theta" , theta ) &&
132 ntuple2->addAndBind( "phi" , phi ) &&
133 ntuple2->addAndBind( "type" , partType ) ) )
134 {
135 G4cerr << "Error: unable to add attribute to nTuple2" << G4endl;
136 // Must be cleaned up properly before any exit
137 delete ntuple2;
138 G4Exception("Could not addAndBind ntuple2");
139 }
140 // ntuples["secondaries"] = ntuple2;
141}
142
143void G4ProcessTestAnalysis::finish()
144{
145 // Because of a Lizard feature, ntuples must be deleted at this stage,
146 // not in the destructor (otherwise the ntuples are not stored)
147
148 delete ntuple1;
149 delete ntuple2;
150
151 /*
152 std::map< G4String,Lizard::NTuple*,std::less<G4String> >::iterator pos;
153 Lizard::NTuple* ntuple;
154 pos = ntuples.find("primary");
155 if (pos != ntuples.end())
156 {
157 ntuple = pos->second;
158 delete ntuple;
159 ntuples["primary"] = 0;
160 }
161
162 pos = ntuples.find("secondaries");
163 if (pos != ntuples.end())
164 {
165 ntuple = pos->second;
166 delete ntuple;
167 ntuples["secondary"] = 0;
168 }
169 */
170 histoManager->store("1");
171 histoManager->store("2");
172 histoManager->store("3");
173 histoManager->store("4");
174 histoManager->store("5");
175}
176
177void G4ProcessTestAnalysis::analyseSecondaries(const G4ParticleChange* particleChange)
178{
179 G4int nSecondaries = particleChange->GetNumberOfSecondaries();
180 for (G4int i = 0; i < nSecondaries; i++)
181 {
182 G4Track* finalParticle = particleChange->GetSecondary(i) ;
183
184 // The quantities bound to ntuple2
185 px = (finalParticle->GetMomentum()).x();
186 py = (finalParticle->GetMomentum()).y();
187 pz = (finalParticle->GetMomentum()).z();
188 p = std::sqrt(px*px + py*py + pz*pz);
189 e = finalParticle->GetTotalEnergy();
190 eKin = finalParticle->GetKineticEnergy();
191 theta = (finalParticle->GetMomentum()).theta();
192 phi = (finalParticle->GetMomentum()).phi();
193
194 partType = 0;
195 G4String particleName = finalParticle->GetDefinition()->GetParticleName();
196
197 if (particleName == "e-") partType = 1;
198 else if (particleName == "e+") partType = 2;
199 else if (particleName == "gamma") partType = 3;
200
201 // Fill ntuple
202 ntuple2->addRow();
203
204 // Fill histograms
205 IHistogram1D* h = histoManager->retrieveHisto1D("3");
206 h->fill(eKin);
207 h = histoManager->retrieveHisto1D("4");
208 h->fill(theta);
209 h = histoManager->retrieveHisto1D("5");
210 h->fill(phi);
211
212 /*
213 std::map< G4String,Lizard::NTuple*,std::less<G4String> >::iterator pos;
214 pos = ntuples.find("secondaries");
215 if (pos != ntuples.end())
216 {
217 Lizard::NTuple* ntuple2 = pos->second;
218 ntuple2->addRow();
219 }
220
221 std::map< G4String,IHistogram1D*,std::less<G4String> >::iterator pos1;
222
223 pos1 = histo1D.find("eKin");
224 if (pos1 != histo1D.end())
225 {
226 IHistogram1D* h = pos1->second;
227 h->fill(eKin);
228 }
229
230 pos1 = histo1D.find("theta");
231 if (pos1 != histo1D.end())
232 {
233 IHistogram1D* h = pos1->second;
234 h->fill(theta);
235 }
236
237 pos1 = histo1D.find("phi");
238 if (pos1 != histo1D.end())
239 {
240 IHistogram1D* h = pos1->second;
241 h->fill(phi);
242 }
243 */
244
245 G4cout << "==== Final "
246 << particleName << " "
247 << "energy: " << e/MeV << " MeV, "
248 << "eKin: " << eKin/MeV << " MeV, "
249 << "(px,py,pz): ("
250 << px/MeV << ","
251 << py/MeV << ","
252 << pz/MeV << ") MeV "
253 << G4endl;
254 }
255}
256
257void G4ProcessTestAnalysis::analyseGeneral(const G4Track& track,
258 const G4ParticleChange* particleChange)
259{
260 // Primary physical quantities
261
262 energyChange = particleChange->GetEnergyChange();
263 initialEnergy = track.GetKineticEnergy();
264
265 G4ThreeVector eChange = particleChange->CalcMomentum(energyChange,
266 (*particleChange->GetMomentumChange()),
267 particleChange->GetMassChange());
268
269 pxChange = eChange.x();
270 pyChange = eChange.y();
271 pzChange = eChange.z();
272 thetaChange = particleChange->GetPositionChange()->theta();
273 phiChange = particleChange->GetPositionChange()->phi();
274 eDeposit = particleChange->GetLocalEnergyDeposit();
275
276 nElectrons = 0;
277 nPositrons = 0;
278 nPhotons = 0;
279
280 G4int nSecondaries = particleChange->GetNumberOfSecondaries();
281
282 for (G4int i = 0; i < nSecondaries; i++)
283 {
284 G4Track* finalParticle = particleChange->GetSecondary(i) ;
285 G4String particleName = finalParticle->GetDefinition()->GetParticleName();
286
287 if (particleName == "e-") nElectrons++;
288 else if (particleName == "e+") nPositrons++;
289 else if (particleName == "gamma") nPhotons++;
290 }
291
292 // Fill ntuple
293 ntuple1->addRow();
294
295 // Fill histograms
296 IHistogram1D* h = histoManager->retrieveHisto1D("1");
297 h ->fill(float(nSecondaries));
298 h = histoManager->retrieveHisto1D("2");
299 h->fill(eDeposit);
300
301 /*
302 std::map< G4String,Lizard::NTuple*,std::less<G4String> >::iterator pos;
303 pos = ntuples.find("primary");
304 if (pos != ntuples.end())
305 {
306 Lizard::NTuple* ntuple1 = pos->second;
307 ntuple1->addRow();
308 }
309 std::map< G4String,IHistogram1D*,std::less<G4String> >::iterator pos1;
310 pos1 = histo1D.find("nSec");
311 if (pos1 != histo1D.end())
312 {
313 IHistogram1D* h2 = pos1->second;
314 h2->fill(float(nSecondaries));
315 }
316
317 pos1 = histo1D.find("eDeposit");
318 if (pos1 != histo1D.end())
319 {
320 IHistogram1D* h2 = pos1->second;
321 h2->fill(eDeposit);
322 }
323 */
324}
325
326G4ProcessTestAnalysis* G4ProcessTestAnalysis::getInstance()
327{
328 if (instance == 0) instance = new G4ProcessTestAnalysis;
329 return instance;
330}
Note: See TracBrowser for help on using the repository browser.