Geant4  10.01
XrayFluoAnalysisManager.cc
Go to the documentation of this file.
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: XrayFluoAnalysisManager.cc
28 // GEANT4 tag $Name: xray_fluo-V03-02-00
29 //
30 // Author: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
31 //
32 // History:
33 // -----------
34 // 20 Jul 2007 A.Mantero, code cleaning and update. Replaced histos with tuple
35 // 11 Jul 2003 A.Mantero, code cleaning / Plotter-XML addiction
36 // Sep 2002 A.Mantero, AIDA3.0 Migration
37 // 06 Dec 2001 A.Pfeiffer updated for singleton
38 // 30 Nov 2001 Guy Barrand : migrate to AIDA-2.2.
39 // 28 Nov 2001 Elena Guardincerri Created
40 //
41 // -------------------------------------------------------------------
42 #include "g4root.hh"
43 
44 #include "G4VProcess.hh"
46 #include "G4Step.hh"
48 #include "G4VPhysicalVolume.hh"
49 #include "G4Gamma.hh"
50 #include "G4Electron.hh"
51 #include "G4Proton.hh"
52 #include "G4SystemOfUnits.hh"
53 
54 
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58 
59 namespace {
60  //Mutex to acquire access to singleton instance check/creation
61  G4Mutex instanceMutex = G4MUTEX_INITIALIZER;
62  //Mutex to acquire accss to histograms creation/access
63  //It is also used to control all operations related to histos
64  //File writing and check analysis
65  G4Mutex dataManipulationMutex = G4MUTEX_INITIALIZER;
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 
71  :outputFileName("xrayfluo"), phaseSpaceFlag(false), physicFlag (false),
72  gunParticleEnergies(0), gunParticleTypes(0)
73 {
74  dataLoaded = false;
76 
77  //creating the messenger
79 
80  //Instantiate the analysis manager
81  G4AnalysisManager::Instance();
82 
83  G4cout << "XrayFluoAnalysisManager created" << G4endl;
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 
89 {
92  if ( gunParticleTypes ) delete gunParticleTypes;
93  gunParticleTypes = 0;
94 
95  delete instance;
96  instance = 0;
97 
98  delete G4AnalysisManager::Instance();
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
104 
105 {
106  G4AutoLock l(&instanceMutex);
107  if (instance == 0) {instance = new XrayFluoAnalysisManager;}
108  return instance;
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 
114 {
115  G4AutoLock l(&dataManipulationMutex);
116  // Get analysis manager
117  G4AnalysisManager* man = G4AnalysisManager::Instance();
118  // Open an output file
119  man->OpenFile(outputFileName);
120  man->SetVerboseLevel(1);
121  man->SetFirstHistoId(1);
122  man->SetFirstNtupleId(1);
123 
124  G4cout << "Open output file: " << outputFileName << G4endl;
125 
126  if (phaseSpaceFlag)
127  {
128  // Book output Tuple (ID = 1)
129  man->CreateNtuple("101","OutputNTuple");
130  man->CreateNtupleIColumn("particle");
131  man->CreateNtupleDColumn("energies");
132  man->CreateNtupleDColumn("momentumTheta");
133  man->CreateNtupleDColumn("momentumPhi");
134  man->CreateNtupleIColumn("processes");
135  man->CreateNtupleIColumn("material");
136  man->CreateNtupleIColumn("origin");
137  man->CreateNtupleDColumn("depth");
138  man->FinishNtuple();
139  G4cout << "Created ntuple for phase space" << G4endl;
140  }
141  else {
142  // Book histograms
143  man->CreateH1("h1","Energy Deposit", 500,0.,10.); //20eV def.
144  man->CreateH1("h2","Gamma born in the sample", 100,0.,10.);
145  man->CreateH1("h3","Electrons born in the sample", 100,0.,10.);
146  man->CreateH1("h4","Gammas leaving the sample", 300,0.,10.);
147  man->CreateH1("h5","Electrons leaving the sample ",200000 ,0.,10.0); // .05 eV def.
148  man->CreateH1("h6","Gammas reaching the detector", 100,0.,10.);
149  man->CreateH1("h7","Spectrum of the incident particles", 100,0.,10.);
150  man->CreateH1("h8","Protons reaching the detector", 100,0.,10.);
151  man->CreateH1("h9","Protons leaving the sample", 100,0.,10.);
152  G4cout << "Created histos" << G4endl;
153  }
154 }
155 
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 
159 {
160  G4AutoLock l(&dataManipulationMutex);
161 
162  //Do not allow more than one thread to trigger the file reading
163  if (dataLoaded)
164  return;
165 
166  // Get analysis reader manager
167  G4AnalysisReader* analysisReader = G4AnalysisReader::Instance();
168  analysisReader->SetVerboseLevel(1);
169 
170  //This is for testing purposes
171  G4int ntupleId = analysisReader->GetNtuple("101",fileName);
172  G4cout << "Got ntupleId: " << ntupleId << G4endl;
173 
174  gunParticleEnergies = new std::vector<G4double>;
175  gunParticleTypes = new std::vector<G4String>;
176 
177  G4int particleID, processID;
179  analysisReader->SetNtupleIColumn("processes", processID);
180  analysisReader->SetNtupleDColumn("energies", energy);
181  analysisReader->SetNtupleIColumn("particles", particleID);
182 
183  while (analysisReader->GetNtupleRow() )
184  {
185  if (raileighFlag ^ (!raileighFlag && (processID == 1 ||
186  processID == 2)) )
187  {
188  gunParticleEnergies->push_back(energy*MeV);
189  if (particleID == 1)
190  gunParticleTypes->push_back("gamma");
191  else if (particleID == 0)
192  gunParticleTypes->push_back("e-");
193  }
194  }
195 
196  G4cout << "Maximum mumber of events: "<< gunParticleEnergies->size() <<G4endl;
197 
198  dataLoaded= true;
199 }
200 
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 
203 const std::pair<G4double,G4String> XrayFluoAnalysisManager::GetEmittedParticleEnergyAndType()
204 {
205  G4AutoLock l(&dataManipulationMutex);
206  std::pair<G4double,G4String> result;
207 
209  {
212  result.first = energy;
213  result.second = name;
214  }
215 
217  return result;
218 }
219 
220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
221 
222 
224 {
225  G4AutoLock l(&dataManipulationMutex);
226  G4cout << "Going to save histograms" << G4endl;
227  // Save histograms
228  G4AnalysisManager* man = G4AnalysisManager::Instance();
229  man->Write();
230  man->CloseFile();
231 }
232 
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
234 
235 
237 {
238  physicFlag = val;
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
242 
244 {
245  G4AutoLock l(&dataManipulationMutex);
246  G4AnalysisManager* man = G4AnalysisManager::Instance();
247 
248  if (phaseSpaceFlag){
249  G4ParticleDefinition* particleType= 0;
250  G4String parentProcess="";
251  G4ThreeVector momentum(0.,0.,0.);
252  G4double particleEnergy=0;
253  G4String sampleMaterial="";
254  G4double particleDepth=0;
255  G4int isBornInTheSample=0;
257 
258  // Select volume from which the step starts
259  if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
260  G4ThreeVector creationPos = aStep->GetTrack()->GetVertexPosition();
261  // qua serve ancora una selezione: deve essere interno al sample
262  //codice "a mano" per controllare il volume di orgine della traccia
263 
264  G4VPhysicalVolume* creationPosVolume = detector->GetGeometryNavigator()->LocateGlobalPointAndSetup(creationPos);
265 
266  // if physicflag is on, I record all the data of all the particle born in the sample.
267  // If it is off (but phase space production is on) I record data
268  // only for particles creted inside the sample and exiting it.
269  if (physicFlag ^ (!physicFlag &&
271  ))
272  {
273  // extracting information needed
274  particleType = aStep->GetTrack()->GetDynamicParticle()->GetDefinition();
275  momentum = aStep->GetTrack()->GetDynamicParticle()->GetMomentum();
276  particleEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();
277  if (creationPosVolume->GetName() == "Sample" ) {
278  isBornInTheSample = 1;
279  particleDepth = creationPosVolume->GetLogicalVolume()->GetSolid()
280  ->DistanceToOut(creationPos, G4ThreeVector(0,0,-1));
281  }
282  else {
283  particleDepth = -1;
284  }
285  // metodo per ottenere la profondita' "a mano"
286  // particleDepth = std::abs(creationPos.z() - detector->GetSampleIlluminatedFaceCoord());
287 
288  G4int parent=0;
289  if(aStep->GetTrack()->GetCreatorProcess())
290  {
291  parentProcess = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
292  parent = 5;
293  // hack for HBK
294  if (parentProcess == "") parent = 0;
295  if (parentProcess == "ioni") parent = 1;
296  if (parentProcess == "LowEnPhotoElec") parent = 2;
297  if (parentProcess == "Transportation") parent = 3;// should never happen
298  if (parentProcess == "initStep") parent = 4;// should never happen
299  }
300  else {
301  parentProcess = "Not Known -- (primary generator + Rayleigh)";
302  parent = 5;
303  }
304  G4int sampleMat=0;
305  if(aStep->GetTrack()){
306  sampleMaterial = aStep->GetTrack()->GetMaterial()->GetName();
307  if (sampleMaterial == ("Dolorite" || "Anorthosite" || "Mars1" || "IceBasalt" || "HPGe")) sampleMat=1;
308  }
309 
310  G4int part = -1 ;
311  if (particleType == G4Gamma::Definition()) part =1;
312  if (particleType == G4Electron::Definition()) part = 0;
313  if (particleType == G4Proton::Definition()) part = 2;
314 
315  //Fill ntuple
316  man->FillNtupleIColumn(1,0, part);
317  man->FillNtupleDColumn(1,1,particleEnergy);
318  man->FillNtupleDColumn(1,2,momentum.theta());
319  man->FillNtupleDColumn(1,3,momentum.phi());
320  man->FillNtupleIColumn(1,4,parent);
321  man->FillNtupleIColumn(1,5,sampleMat);
322  man->FillNtupleIColumn(1,6,isBornInTheSample);
323  man->FillNtupleDColumn(1,7,particleDepth);
324  man->AddNtupleRow(1);
325 
326  }
327  }
328  }
329 
330  // Normal behaviour, without creation of phase space
331  else
332  {
333 
334  // Filling the histograms with data, passing thru stepping.
335 
336  // Select volume from wich the step starts
337  if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
338  // Select volume from wich the step ends
339  if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" ) {
340  // Select the particle type
341  if ((aStep->GetTrack()->GetDynamicParticle()
343 
344  {
345  G4double gammaLeavingSample =
346  aStep->GetPreStepPoint()->GetKineticEnergy();
347  man->FillH1(4,gammaLeavingSample/keV);
348 
349  }
350  else if ((aStep->GetTrack()->GetDynamicParticle()
352  {
353  G4double eleLeavingSample =
354  aStep->GetPreStepPoint()->GetKineticEnergy();
355  man->FillH1(5,eleLeavingSample/keV);
356  }
357  else if ((aStep->GetTrack()->GetDynamicParticle()
359  {
360  G4double
361  protonsLeavSam = aStep->GetPreStepPoint()->GetKineticEnergy();
362  man->FillH1(9,protonsLeavSam/keV);
363  }
364  }
365  }
366 
367  if((aStep->GetTrack()->GetDynamicParticle()
369  {
370  if(aStep->GetTrack()->GetCurrentStepNumber() == 1)
371  {
372  if(aStep->GetTrack()->GetParentID() != 0)
373  {
374  if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
375  {
376  G4double gammaBornInSample =
377  aStep->GetPreStepPoint()->GetKineticEnergy();
378  man->FillH1(2,gammaBornInSample/keV);
379  }
380  }
381  }
382  }
383  else if ((aStep->GetTrack()->GetDynamicParticle()
385  {
386  if(aStep->GetTrack()->GetCurrentStepNumber() == 1)
387  {
388  if(aStep->GetTrack()->GetParentID() != 0)
389  {
390  if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
391  {
392  G4double eleBornInSample =
393  aStep->GetPreStepPoint()->GetKineticEnergy();
394  man->FillH1(3,eleBornInSample/keV);
395  }
396  }
397  }
398  }
399 
400 
401  if(aStep->GetTrack()->GetNextVolume())
402  {
403 
404 
405  if(aStep->GetTrack()->GetNextVolume()->GetName() == "HPGeDetector")
406 
407  {
408  if ((aStep->GetTrack()->GetDynamicParticle()
410  {
411 
412  G4double gammaAtTheDetPre =
413  aStep->GetPreStepPoint()->GetKineticEnergy();
414  man->FillH1(6,gammaAtTheDetPre/keV);
415  }
416  else if ((aStep->GetTrack()->GetDynamicParticle()
417  ->GetDefinition() ) == G4Proton::Definition() )
418  {
419  G4double protonsAtTheDetPre =
420  aStep->GetPreStepPoint()->GetKineticEnergy();
421  man->FillH1(8,protonsAtTheDetPre);
422  }
423  }
424  }
425  }
426 }
427 
428 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
429 
431 {
432  G4AutoLock l(&dataManipulationMutex);
433  // Filling of Energy Deposition, called by XrayFluoEventAction
434  G4AnalysisManager* man = G4AnalysisManager::Instance();
435 
436  if (!phaseSpaceFlag)
437  man->FillH1(1,energyDep/keV);
438 
439 }
440 
441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
442 
444 {
445  G4AutoLock l(&dataManipulationMutex);
446  // Filling of energy spectrum histogram of the primary generator
447  G4AnalysisManager* man = G4AnalysisManager::Instance();
448  if (!phaseSpaceFlag)
449  man->FillH1(7,energy/keV);
450 
451 }
452 
453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
454 
456 {
457 
458  outputFileName = newName;
459 }
460 
461 
G4int GetParentID() const
static const double MeV
Definition: G4SIunits.hh:193
G4int CreateNtupleIColumn(const G4String &name)
G4bool SetFirstHistoId(G4int firstId)
G4int CreateH1(const G4String &name, const G4String &title, G4int nbins, G4double xmin, G4double xmax, const G4String &unitName="none", const G4String &fcnName="none", const G4String &binSchemeName="linear")
static XrayFluoAnalysisManager * instance
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
void analyseStepping(const G4Step *aStep)
G4String name
Definition: TRTMaterials.hh:40
G4StepStatus GetStepStatus() const
const G4String & GetName() const
Definition: G4Material.hh:176
void SetVerboseLevel(G4int verboseLevel)
static G4Electron * Definition()
Definition: G4Electron.cc:49
G4int CreateNtuple(const G4String &name, const G4String &title)
G4bool SetNtupleDColumn(const G4String &columnName, G4double &value)
G4ParticleDefinition * GetDefinition() const
const std::pair< G4double, G4String > GetEmittedParticleEnergyAndType()
int G4int
Definition: G4Types.hh:78
static XrayFluoDetectorConstruction * GetInstance()
G4bool OpenFile(const G4String &fileName="")
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:163
void LoadGunData(G4String, G4bool)
G4bool SetNtupleIColumn(const G4String &columnName, G4int &value)
G4VPhysicalVolume * GetNextVolume() const
void analysePrimaryGenerator(G4double energy)
static G4Proton * Definition()
Definition: G4Proton.cc:49
G4StepPoint * GetPreStepPoint() const
const G4VProcess * GetCreatorProcess() const
G4bool FillNtupleIColumn(G4int id, G4int value)
G4GLOB_DLL std::ostream G4cout
G4int GetCurrentStepNumber() const
G4VPhysicalVolume * GetPhysicalVolume() const
const G4String & GetName() const
G4bool FillNtupleDColumn(G4int id, G4double value)
static XrayFluoAnalysisManager * getInstance()
bool G4bool
Definition: G4Types.hh:79
Definition: G4Step.hh:76
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
std::vector< G4double > * gunParticleEnergies
const G4ThreeVector & GetVertexPosition() const
G4Material * GetMaterial() const
G4bool SetFirstNtupleId(G4int firstId)
G4int G4Mutex
Definition: G4Threading.hh:161
tools::rroot::ntuple * GetNtuple() const
G4LogicalVolume * GetLogicalVolume() const
G4double energy(const ThreeVector &p, const G4double m)
G4bool FillH1(G4int id, G4double value, G4double weight=1.0)
G4StepPoint * GetPostStepPoint() const
std::vector< G4String > * gunParticleTypes
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
Definition: G4Navigator.cc:122
G4VPhysicalVolume * GetVolume() const
void SetVerboseLevel(G4int verboseLevel)
#define G4endl
Definition: G4ios.hh:61
XrayFluoAnalysisMessenger * analisysMessenger
static const double keV
Definition: G4SIunits.hh:195
G4int CreateNtupleDColumn(const G4String &name)
void analyseEnergyDep(G4double eDep)
G4double GetKineticEnergy() const
double G4double
Definition: G4Types.hh:76
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4Track * GetTrack() const
G4VSolid * GetSolid() const
static G4Gamma * Definition()
Definition: G4Gamma.cc:49
G4ThreeVector GetMomentum() const