source: trunk/examples/advanced/xray_fluorescence/src/XrayFluoAnalysisManager.cc@ 1276

Last change on this file since 1276 was 1230, checked in by garnier, 16 years ago

update to geant4.9.3
File size: 25.4 KB
RevLine 
[807]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// -----------
[1230]34// 20 Jul 2007 A.Mantero, code cleaning and update. Replaced histos with tuple
[807]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#ifdef G4ANALYSIS_USE
43
44#include "G4VProcess.hh"
45#include <fstream>
46#include "G4ios.hh"
47#include "XrayFluoAnalysisManager.hh"
48#include "G4Step.hh"
[1230]49#include "XrayFluoDetectorConstruction.hh"
50#include "G4VPhysicalVolume.hh"
[807]51
52XrayFluoAnalysisManager* XrayFluoAnalysisManager::instance = 0;
53
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55
56XrayFluoAnalysisManager::XrayFluoAnalysisManager()
57 :outputFileName("xrayfluo"), visPlotter(false), phaseSpaceFlag(false), physicFlag (false), persistencyType("xml"),
[1230]58 deletePersistencyFile(true), gunParticleEnergies(0), gunParticleTypes(0), analysisFactory(0), tree(0), treeDet(0), histogramFactory(0), plotter(0)
[807]59{
60 //creating the messenger
61 analisysMessenger = new XrayFluoAnalysisMessenger(this);
62
63 // Hooking an AIDA compliant analysis system.
64 // creating Analysis factroy, necessary to create/manage analysis
65 analysisFactory = AIDA_createAnalysisFactory();
66
67 CreatePersistency(outputFileName,persistencyType);
68
69 if(analysisFactory) {
70
71
72 // Creating the plotter factory
73
74 if (visPlotter){
75 plotterFactory = analysisFactory->createPlotterFactory();
76 }
77 }
78
79 G4cout << "XrayFluoAnalysisManager created" << G4endl;
80}
81
82//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83
84XrayFluoAnalysisManager::~XrayFluoAnalysisManager()
85{
86 delete histogramFactory;
87 histogramFactory=0;
88
89 delete analysisFactory;
90 analysisFactory = 0;
91
92 delete plotterFactory;
93 plotterFactory=0;
94
95 delete tree;
96
97 if ( gunParticleEnergies ) delete gunParticleEnergies;
98 gunParticleEnergies = 0;
99 if ( gunParticleTypes ) delete gunParticleTypes;
100 gunParticleTypes = 0;
101
102 delete instance;
103
104 G4cout << "XrayFluoAnalysisManager delete" << G4endl;
105
106}
107
108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109
110XrayFluoAnalysisManager* XrayFluoAnalysisManager::getInstance()
111
112{
113 if (instance == 0) {instance = new XrayFluoAnalysisManager;}
114 return instance;
115}
116
117//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118
119void XrayFluoAnalysisManager::CreatePersistency(G4String fileName,G4String persistencyType,
120 G4bool readOnly, G4bool createNew)
121{
122
123 if (tree) {
124 delete tree;
125 tree = 0;
126 }
127 if (treeDet) {
128 delete treeDet;
129 treeDet =0;
130 }
131 if (histogramFactory) {
132 delete histogramFactory;
133 histogramFactory = 0;
134 }
135 if (tupleFactory) {
136 delete tupleFactory;
137 tupleFactory = 0;
138 }
139
140 if (tupleDetFactory) {
141 delete tupleDetFactory;
142 tupleDetFactory = 0;
143 }
144 G4String fileNameDet = fileName+"Detector";
145
146 AIDA::ITreeFactory* treeFactory = analysisFactory->createTreeFactory();
147 if(treeFactory) {
148 if (persistencyType == "hbook") {
149 fileName = fileName + ".hbk";
150 fileNameDet = fileNameDet + ".hbk";
151 }
152 else if (persistencyType == "xml"){
153 fileName = fileName + ".xml";
154 fileNameDet = fileNameDet + ".xml";
155 }
156
157 if (phaseSpaceFlag) {
158
159 tree = treeFactory->create(fileName,persistencyType,readOnly,createNew); // output file
160 treeDet = treeFactory->create(fileNameDet,persistencyType,readOnly,createNew); // output file
161 if(analysisFactory) {
162
163 tupleDetFactory = analysisFactory->createTupleFactory(*treeDet);
164 tupleFactory = analysisFactory->createTupleFactory(*tree);
165
166 }
167 }
168
169 // trees to be stored in case of phase-space production
170 else {
171
172 tree = treeFactory->create(fileName,persistencyType,readOnly,createNew); // output file
173
174 if(analysisFactory) {
175
176 histogramFactory = analysisFactory->createHistogramFactory(*tree);
177
178 }
179 }
180
181 delete treeFactory; // Will not delete the ITree.
182 }
183}
184
185
186void XrayFluoAnalysisManager::book()
187{
188
189 if (phaseSpaceFlag) {
190 /*
191 // Book clouds
192
193 cloud_1 = histogramFactory->createCloud1D("Gamma Exting Sample","ciao!",-1);
194 cloud_2 = histogramFactory->createCloud1D("Gamma Incident on the detector","ciao!",-1);
195 cloud_3 = histogramFactory->createCloud1D("Electrons Exiting the Sample","ciao!",-1);
196 beamCloud = histogramFactory->createCloud1D("Incident Radiation Spectrum","ciao!",-1);
197 */
198 // Book output Tuple
199
200 // Book tuple
201 std::vector<std::string> columnNames;
[1230]202 columnNames.push_back("particle");
203 columnNames.push_back("energies");
[807]204 columnNames.push_back("momentumTheta");
205 columnNames.push_back("momentumPhi");
[1230]206 columnNames.push_back("processes");
207 columnNames.push_back("material");
208 columnNames.push_back("origin");
209 columnNames.push_back("depth");
[807]210
211 std::vector<std::string> columnTypes;
212 columnTypes.push_back("int");
213 columnTypes.push_back("double");
214 columnTypes.push_back("double");
215 columnTypes.push_back("double");
216 columnTypes.push_back("int"); // useful for hbk
[1230]217 columnTypes.push_back("int"); // useful for hbk
218 columnTypes.push_back("int");
219 columnTypes.push_back("double");
[807]220
221 tupleFluo = tupleFactory->create("10", "Total Tuple", columnNames, columnTypes, "");
222 assert(tupleFluo);
223
224 }
225
226 else {
227 // Book histograms
228
229 histo_1 = histogramFactory->createHistogram1D("1","Energy Deposit", 500,0.,10.); //20eV def.
230 histo_2 = histogramFactory->createHistogram1D("2","Gamma born in the sample", 100,0.,10.);
231 histo_3 = histogramFactory->createHistogram1D("3","Electrons born in the sample", 100,0.,10.);
232 histo_4 = histogramFactory->createHistogram1D("4","Gammas leaving the sample", 300,0.,10.);
233 histo_5 = histogramFactory->createHistogram1D("5","Electrons leaving the sample ",200000 ,0.,10.0); // .05 eV def.
234 histo_6 = histogramFactory->createHistogram1D("6","Gammas reaching the detector", 100,0.,10.);
235 histo_7 = histogramFactory->createHistogram1D("7","Spectrum of the incident particles", 100,0.,10.);
236 histo_8 = histogramFactory->createHistogram1D("8","Protons reaching the detector", 100,0.,10.);
237 histo_9 = histogramFactory->createHistogram1D("9","Protons leaving the sample", 100,0.,10.);
238
239 // Debugging-purpose Histos
240
241 //histo_10 = histogramFactory->createHistogram1D("10","Photon Origin", 4,0.,3.);
242 //histo_11 = histogramFactory->createHistogram1D("11","Spectrum from LowEnPhotoELectric", 300,0.,10.);
243 //histo_12 = histogramFactory->createHistogram1D("12","Spectrum From the other processes (unknown)", 300,0.,10.);
244
245 // IHistogram2D* histo_20 = histogramFactory->
246 //create2D("20","Phi, Theta",80 ,-3.14,3.14,80,0.,3.14);
247 }
248}
249
250//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251
252void XrayFluoAnalysisManager::LoadGunData(G4String fileName, G4bool raileighFlag) {
253
254 // gunParticleEnergies = new std::vector<G4double>;
255 // gunParticleTypes = new std::vector<G4String>;
256
257 G4String ext = fileName.substr(fileName.size()-3,fileName.size()-1);
258 G4String persistencyType;
259
260 if (ext == "xml") {
261
262 persistencyType = "xml";
263 }
264 else if (ext == "hbk") {
265
266 persistencyType = "hbook";
267 }
268
269 gunParticleEnergies = new std::vector<G4double>;
270 gunParticleTypes = new std::vector<G4String>;
271
[1230]272 // dal punto di vista della programmazione e' un po' una porcata.....
273
[807]274 AIDA::ITreeFactory* treeDataFactory = analysisFactory->createTreeFactory();
[1230]275 AIDA::ITree* treeData = treeDataFactory->create(fileName,persistencyType,false, false); // input file
[807]276 AIDA::IManagedObject* mo = treeData->find("10");
[1230]277 assert(mo);
278 G4cout <<"mo Name:" << mo->name()<< G4endl;
[807]279 AIDA::ITuple* tupleData = dynamic_cast<AIDA::ITuple*>(mo);
[1230]280
281 //AIDA::ITuple* tupleData = (AIDA::ITuple*) treeData->find("10");
282
283 assert(tupleData);
284 // da riscrivere un pochino melgio. Del resto serve per pochissmo, poi non serve piu'
285
286
[807]287 tupleData->start();
[1230]288 G4int processesIndex = tupleData->findColumn("processes");
289 G4int energyIndex = tupleData->findColumn("energies");
290 G4int particleIndex = tupleData->findColumn("particle");
[807]291
292 while (tupleData->next()) {
[1230]293 if (raileighFlag ^ (!raileighFlag && (tupleData->getInt(processesIndex) == 1 ||
294 tupleData->getInt(processesIndex) == 2)) ) {
295 gunParticleEnergies->push_back(tupleData->getDouble(energyIndex)*MeV);
296 if (tupleData->getInt(particleIndex) == 1 ) gunParticleTypes->push_back("gamma");
297 if (tupleData->getInt(particleIndex) == 0 ) gunParticleTypes->push_back("e-");
[807]298 }
299
300 }
301 G4cout << "Maximum mumber of events: "<< gunParticleEnergies->size() <<G4endl;
302
303
304}
305
306std::vector<G4double>* XrayFluoAnalysisManager::GetEmittedParticleEnergies() {
307
308 return gunParticleEnergies;
309
310}
311
312std::vector<G4String>* XrayFluoAnalysisManager::GetEmittedParticleTypes() {
313
314 return gunParticleTypes;
315
316}
317
318
319void XrayFluoAnalysisManager::finish()
320{
321
322 if (tupleFluo) {ExtractData();};
323
324 if(plotter) {
325 // Wait for the keyboard return to avoid destroying the plotter window too quickly.
326 G4cout << "Press <ENTER> to exit" << G4endl;
327 G4cin.get();
328 plotter->hide(); //hide plotter windows, but doesn't delete plotter
329 }
330
331 if(tree) {
332 tree->commit(); // Write histos in file.
333 tree->close();
334 }
335 if (treeDet) {
336 treeDet->commit();
337 treeDet->close();
338 }
339
340 deletePersistencyFile = false;
341
342}
343
344//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345
346void XrayFluoAnalysisManager::InitializePlotter()
347{
348
349 // If speciefied (visPlotter=true),
350 // a window for the visulizaton of partial results is created
351 if(plotterFactory && visPlotter && deletePersistencyFile)
352 {
353 plotter = plotterFactory->create();
354 // Set the page title :
355 plotter->setParameter("pageTitle","XrayFluo");
356 }
357
358 if(plotter && visPlotter) {
359 plotter->show(); // shows plotter window
360 }
361
362}
363//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
364
365void XrayFluoAnalysisManager::PlotCurrentResults()
366{
367 if(plotter) {
368 if (phaseSpaceFlag){
[1230]369
370
371 // altra porcata di programmazione e cmq mi serve un istogramma. Da rivedere il tutto.
[807]372 // Plotting the spectrum of Gamma exiting the sample - cloud1
373 AIDA::ICloud1D& cloud = *cloud_1;
374 AIDA::IFilter* filterGamma = tupleFactory->createFilter(
375 " Particle == std::string(\"gamma\") ");
376 AIDA::IEvaluator* evaluatorEnergy = tupleFactory->createEvaluator("Energies");
[1230]377
378 // di nuovo, del resto serve solo per un attimo.. da correggere cmq.
379
[807]380 filterGamma->initialize(*tupleFluo);
381 evaluatorEnergy->initialize(*tupleFluo);
382 tupleFluo->project(cloud,*evaluatorEnergy,*filterGamma);
383
384 plotter->currentRegion().plot( cloud, "Exiting Gammas " );
385 plotter->refresh();
386 }
387
388 else{
389 // Plotting the Detector Energy Deposit - histo_1
390 AIDA::IHistogram1D& histo1p = *histo_1;
391 plotter->currentRegion().plot( histo1p, "Detector Energy Deposition" );
392 plotter->refresh();
393 }
394
395 }
396}
397
398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
399G4bool XrayFluoAnalysisManager::GetDeletePersistencyFileFlag()
400{
401 return deletePersistencyFile;
402}
403
404//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
405
406void XrayFluoAnalysisManager::SetPhysicFlag(G4bool val)
407{
408 physicFlag = val;
409}
410
411//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
412
413
414
415
416
417void XrayFluoAnalysisManager::analyseStepping(const G4Step* aStep)
418{
419
420 if (phaseSpaceFlag){
421
422
423 G4String particleType="";
424 G4String parentProcess="";
425 G4ThreeVector momentum=0;
426 G4double particleEnergy=0;
[1230]427 G4String sampleMaterial="";
428 G4double particleDepth=0;
429 G4int isBornInTheSample=0;
430 XrayFluoDetectorConstruction* detector = XrayFluoDetectorConstruction::GetInstance();
431
432 // Select volume from which the step starts
[807]433 if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
[1230]434 G4ThreeVector creationPos = aStep->GetTrack()->GetVertexPosition();
435 // qua serve ancora una selezione: deve essere interno al sample
436 //codice "a mano" per controllare il volume di orgine della traccia
437 /*
438 G4Double sampleXplus = detector->GetSampleXplusFaceCoord();
439 G4Double sampleXminus = detector->GetSampleXminusFaceCoord();
440 G4Double sampleYplus = detector->GetSampleYplusFaceCoord();
441 G4Double sampleYminus = detector->GetSampleYminusFaceCoord();
442 G4Double sampleZplus = detector->GetSampleYplusFaceCoord();
443 G4Double sampleZminus = detector->GetSampleYminusFaceCoord();
444
445 G4bool ZsampleCheck = (creationPos.z()<=sampleZminus) && (creationPos.z() >= sampleZplus);
446 G4bool XsampleCheck = (creationPos.x()<=sampleZminus) && (creationPos.x() >= sampleXplus);
447 G4bool YsampleCheck = (creationPos.y()<=sampleZminus) && (creationPos.y() >= sampleYplus);
448 */
449 G4VPhysicalVolume* creationPosVolume = detector->GetGeometryNavigator()->LocateGlobalPointAndSetup(creationPos);
450
451 // if physicflag is on, I record all the data of all the particle born in the sample.
452 // If it is off (but phase space production is on) I record data
453 // only for particles creted inside the sample and exiting it.
454 if (physicFlag ^ (!physicFlag &&
455 (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
456 //ZsampleCheck && XsampleCheck && YsampleCheck)) {
457 ) ) {
[807]458 //
459 //
460 // G4cout << "physicFlag: "<< physicFlag << G4endl
461 // << "NextVolume: "<< aStep->GetTrack()->GetNextVolume()->GetName() << G4endl;
462
463 // extracting information needed
464 particleType = aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetParticleName();
465 momentum = aStep->GetTrack()->GetDynamicParticle()->GetMomentum();
466 particleEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();
[1230]467 if (creationPosVolume->GetName() == "Sample" ) {
468 isBornInTheSample = 1;
469 particleDepth = creationPosVolume->GetLogicalVolume()->GetSolid()->DistanceToOut(creationPos, G4ThreeVector(0,0,-1));
470 }
471 else {
472 particleDepth = -1;
473 }
474 // metodo per ottenere la profondita' "a mano"
475 // particleDepth = std::abs(creationPos.z() - detector->GetSampleIlluminatedFaceCoord());
476
[807]477 G4int parent;
478 if(aStep->GetTrack()->GetCreatorProcess()){
479 parentProcess = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
[1230]480 parent = 5;
481 // hack for HBK
482 if (parentProcess == "") parent = 0;
483 if (parentProcess == "IONI") parent = 1;
484 if (parentProcess == "LowEnPhotoElec") parent = 2;
485 if (parentProcess == "Transportation") parent = 3;// should never happen
486 if (parentProcess == "initStep") parent = 4;// should never happen
[807]487 }
488 else {
[1230]489 parentProcess = "Not Known -- (primary generator + Rayleigh)";
490 parent = 5;
[807]491 }
[1230]492 G4int sampleMat=0;
493 if(aStep->GetTrack()){
494 sampleMaterial = aStep->GetTrack()->GetMaterial()->GetName();
495 if (sampleMaterial == ("Dolorite" || "Anorthosite" || "Mars1" || "IceBasalt" || "HPGe")) sampleMat=1;
496 }
[807]497 // filling tuple
498
499 // G4cout<< particleType << G4endl;
[1230]500 G4int part = -1 ;
[807]501 if (particleType == "gamma") part =1;
502 if (particleType == "e-") part = 0;
[1230]503 if (particleType == "proton") part = 2;
[807]504
[1230]505
[807]506 tupleFluo->fill(0,part);
507 tupleFluo->fill(1,particleEnergy);
508 tupleFluo->fill(2,momentum.theta());
509 tupleFluo->fill(3,momentum.phi());
[1230]510 tupleFluo->fill(4,parent); //hacked to be used with hbk
511 tupleFluo->fill(5,sampleMat); //hacked to be used with hbk
512 tupleFluo->fill(6,isBornInTheSample);
513 tupleFluo->fill(7,particleDepth);
[807]514
515 tupleFluo->addRow();
516 }
517 }
518 }
519
520
521 // Normal behaviour, without creation of phase space
522 else {
523
524 G4double gammaAtTheDetPre=0;
525 G4double protonsAtTheDetPre=0;
526 G4double gammaLeavingSample=0;
527 //G4double gammaLeavingSamplePhi=0;
528 //G4double gammaLeavingSampleTheta=0;
529
530 G4double eleLeavingSample=0;
531 G4double protonsLeavSam=0;
532 G4double gammaBornInSample=0;
533 G4double eleBornInSample=0;
534
535 // Filling the histograms with data, passing thru stepping.
536
537
538 // Select volume from wich the step starts
539 if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
540 // Select volume from wich the step ends
541 if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" ) {
542 // Select the particle type
543 if ((aStep->GetTrack()->GetDynamicParticle()
544 ->GetDefinition()-> GetParticleName()) == "gamma" )
545
546 {
547
548
549 // Control histos, used for debugging purpose.
550
551 // if(aStep->GetTrack()->GetCreatorProcess()){
552 // G4String process = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
553 // G4cout << "Il processo di origine e' " << process << G4endl;
554 // if(histo_10) {
555 // histo_10->fill(1.);
556 // gammaLeavingSample =
557 // (aStep->GetPreStepPoint()->GetKineticEnergy());
558 // if(histo_11) {
559 // histo_11->fill(gammaLeavingSample/keV);
560
561 // }
562 // }
563 // }
564 // else {
565 // //G4cout << "Sembra che non ci sia un processo d'origine" << G4endl;
566 // if(histo_10) {
567 // histo_10->fill(2.);
568
569 // gammaLeavingSample =
570 // (aStep->GetPreStepPoint()->GetKineticEnergy());
571 // if(histo_12) {
572 // histo_12->fill(gammaLeavingSample/keV);
573
574 // }
575
576 // }
577 // }
578
579
580 // //
581 // Filling Histos //
582 // //
583
584
585 gammaLeavingSample =
586 (aStep->GetPreStepPoint()->GetKineticEnergy());
587 if(histo_4) {
588 histo_4->fill(gammaLeavingSample/keV);
589
590 }
591
592 // Other debugging purpose histos
593 /* gammaLeavingSamplePhi =
594 (aStep->GetPreStepPoint()->GetMomentumDirection().phi());
595 G4cout << "questo e' Phi: " << gammaLeavingSamplePhi << G4endl;
596 gammaLeavingSampleTheta =
597 (aStep->GetPreStepPoint()->GetMomentumDirection().theta());
598 G4cout << "questo e' Theta: " << gammaLeavingSampleTheta << G4endl;
599 if(histo_20) {
600 // G4cout << "histo_20 esiste" << G4endl;
601 histo_20->fill(gammaLeavingSamplePhi,gammaLeavingSampleTheta,1.);
602 } */
603
604
605 }
606 }
607 }
608
609
610
611 if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
612
613 if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" )
614 {
615 if ((aStep->GetTrack()->GetDynamicParticle()
616 ->GetDefinition()-> GetParticleName()) == "e-" )
617 {
618 eleLeavingSample = (aStep->GetPreStepPoint()->GetKineticEnergy());
619 if(histo_5) {
620 histo_5->fill(eleLeavingSample/keV);
621 }
622 }
623 else if ((aStep->GetTrack()->GetDynamicParticle()
624 ->GetDefinition()-> GetParticleName()) == "proton" )
625 {
626 protonsLeavSam = (aStep->GetPreStepPoint()->GetKineticEnergy());
627 if(histo_9) {
628 histo_9->fill(protonsLeavSam/keV);
629 }
630 }
631
632 }
633 }
634
635 // electrons from the detector -- Debugging
636
637// if((aStep->GetTrack()->GetDynamicParticle()
638// ->GetDefinition()-> GetParticleName()) == "e-" ){
639
640// if(1== (aStep->GetTrack()->GetCurrentStepNumber())){
641
642// if(0 != aStep->GetTrack()->GetParentID()){
643// if(aStep->GetTrack()->GetVolume()->GetName() == "HPGeDetector")
644
645// if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" ){
646
647// if(aStep->GetTrack()->GetCreatorProcess()){
648// G4String process = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
649// G4cout << "Origin Porcess Name: " << process << G4endl;
650// if(histo_10) {
651// histo_10->fill(1.);
652// gammaLeavingSample =
653// (aStep->GetPreStepPoint()->GetKineticEnergy());
654// if(histo_11) {
655// histo_11->fill(gammaLeavingSample/keV);
656
657// }
658// }
659// }
660// else {
661// G4cout << "No Origin Process Name" << G4endl;
662// if(histo_10) {
663// histo_10->fill(2.);
664
665// gammaLeavingSample =
666// (aStep->GetPreStepPoint()->GetKineticEnergy());
667// if(histo_12) {
668// histo_12->fill(gammaLeavingSample/keV);
669
670// }
671
672// }
673// }
674// }
675// }
676
677// }
678// }
679
680
681
682
683 if((aStep->GetTrack()->GetDynamicParticle()
684 ->GetDefinition()-> GetParticleName()) == "gamma" )
685
686 {if(1== (aStep->GetTrack()->GetCurrentStepNumber()))
687
688 {if(0 != aStep->GetTrack()->GetParentID())
689
690 {if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
691 {
692 gammaBornInSample = (aStep->GetPreStepPoint()->GetKineticEnergy());
693 if(histo_2) {
694 histo_2->fill(gammaBornInSample/keV);
695 }
696 }
697 }
698 }
699 }
700 if((aStep->GetTrack()->GetDynamicParticle()
701 ->GetDefinition()-> GetParticleName()) == "e-" )
702
703 {if(1== (aStep->GetTrack()->GetCurrentStepNumber()))
704
705 {if(0 != aStep->GetTrack()->GetParentID())
706
707 {if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
708{
709 eleBornInSample = (aStep->GetPreStepPoint()->GetKineticEnergy());
710 if(histo_3) {
711 histo_3->fill(eleBornInSample/keV);
712 }
713 }
714 }
715 }
716 }
717
718 if(aStep->GetTrack()->GetNextVolume()){
719
720 //if(aStep->GetTrack()->GetVolume()->GetName() == "World"){
721
722 if(aStep->GetTrack()->GetNextVolume()->GetName() == "HPGeDetector")
723
724 {
725 if ((aStep->GetTrack()->GetDynamicParticle()
726 ->GetDefinition()-> GetParticleName()) == "gamma" )
727 {
728
729 gammaAtTheDetPre =
730 (aStep->GetPreStepPoint()->GetKineticEnergy());
731 if(histo_6) {
732 histo_6->fill( gammaAtTheDetPre/keV);
733 }
734 }
735 else if ((aStep->GetTrack()->GetDynamicParticle()
736 ->GetDefinition()-> GetParticleName()) == "proton" )
737 {
738 protonsAtTheDetPre =
739 (aStep->GetPreStepPoint()->GetKineticEnergy());
740 if(histo_8) {
741 histo_8->fill( protonsAtTheDetPre/keV);
742 }
743 }
744 }
745 //} close of if(aStep->GetTrack()->GetVolume()->GetName() == "World"){
746 }
747 }
748}
749
750//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751
752void XrayFluoAnalysisManager::ExtractData(){
753
754 if (tupleFluo->rows()) {
755
756
757 // AIDA::IFilter* filterGamma = tupleFactory->createFilter(" Particle == std::string(\"gamma\")");
758 // AIDA::IFilter* filterEminus = tupleFactory->createFilter(" Particle == std::string(\"e-\")");
759
760
[1230]761 AIDA::IFilter* filterAngle = tupleFactory->createFilter(
762
763"(momentumPhi >= (220. * (3.1415926/180.) )) && (momentumPhi <= (230. * (3.1415926/180.) )) && (momentumTheta >= (130. * (3.1415926/180.) )) && (momentumTheta <= (140. * (3.1415926/180.) ))" );
[807]764
765
766 // filterGamma ->initialize(*tupleFluo);
767 // filterEminus ->initialize(*tupleFluo);
768 filterAngle->initialize(*tupleFluo);
769
770 // Create IEvaluator and initialize it to this ITuple
771 // AIDA::IEvaluator* evaluatorEnergy = tupleFactory->createEvaluator("Energies");
772 // evaluatorEnergy->initialize(*tupleFluo);
773
774 // tupleFluo->project(*cloud_1,*evaluatorEnergy,*filterGamma);
775 // tupleFluo->project(*cloud_2,*evaluatorEnergy,*filterAngle);
776 // tupleFluo->project(*cloud_3,*evaluatorEnergy,*filterEminus);
777
778 tupleDetFluo = tupleDetFactory->createFiltered("1", *tupleFluo, *filterAngle);
779 assert(tupleDetFluo);
780 }
781}
782//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
783
784
785void XrayFluoAnalysisManager::analyseEnergyDep(G4double energyDep)
786{
787
788 // Filling of Energy Deposition, called by XrayFluoEventAction
789
790 histo_1->fill(energyDep/keV);
791
792}
793
794//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
795
796void XrayFluoAnalysisManager::analysePrimaryGenerator(G4double energy)
797{
798
799 // Filling of energy spectrum histogram of the primary generator
800
801 if (phaseSpaceFlag){
802 // beamCloud->fill(energy/keV);
803 }
804 else {
805 histo_7->fill(energy/keV);
806 }
807}
808
809
810// not used -- Created for future development
811
812void XrayFluoAnalysisManager::SetOutputFileName(G4String newName)
813{
814
815 outputFileName = newName;
816}
817
818void XrayFluoAnalysisManager::SetOutputFileType(G4String newType)
819{
820
821 persistencyType = newType;
822}
823
824#endif
Note: See TracBrowser for help on using the repository browser.