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

Last change on this file since 1321 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
File size: 25.4 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: 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#ifdef G4ANALYSIS_USE
43
44#include "G4VProcess.hh"
45#include <fstream>
46#include "G4ios.hh"
47#include "XrayFluoAnalysisManager.hh"
48#include "G4Step.hh"
49#include "XrayFluoDetectorConstruction.hh"
50#include "G4VPhysicalVolume.hh"
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"), 
58   deletePersistencyFile(true), gunParticleEnergies(0), gunParticleTypes(0), analysisFactory(0), tree(0), treeDet(0), histogramFactory(0), plotter(0)
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;
202    columnNames.push_back("particle");
203    columnNames.push_back("energies");
204    columnNames.push_back("momentumTheta");
205    columnNames.push_back("momentumPhi");
206    columnNames.push_back("processes");
207    columnNames.push_back("material");
208    columnNames.push_back("origin");
209    columnNames.push_back("depth");
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
217    columnTypes.push_back("int"); // useful for hbk
218    columnTypes.push_back("int"); 
219    columnTypes.push_back("double");
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
272  // dal punto di vista della programmazione e' un po' una porcata.....
273
274  AIDA::ITreeFactory* treeDataFactory = analysisFactory->createTreeFactory();
275  AIDA::ITree* treeData = treeDataFactory->create(fileName,persistencyType,false, false); // input file
276  AIDA::IManagedObject* mo = treeData->find("10");
277  assert(mo);
278  G4cout <<"mo Name:" << mo->name()<< G4endl;
279  AIDA::ITuple* tupleData = dynamic_cast<AIDA::ITuple*>(mo);
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
287  tupleData->start();
288  G4int processesIndex = tupleData->findColumn("processes");
289  G4int energyIndex = tupleData->findColumn("energies");
290  G4int particleIndex = tupleData->findColumn("particle");
291
292  while (tupleData->next()) {
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-");
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){
369
370
371      // altra porcata di programmazione e cmq mi serve un istogramma. Da rivedere il tutto.
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");
377
378      // di nuovo, del resto serve solo per un attimo.. da correggere cmq.
379
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;
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
433    if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
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                        ) ) {
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();
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
477        G4int parent;
478        if(aStep->GetTrack()->GetCreatorProcess()){
479          parentProcess = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
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
487        }       
488        else {
489          parentProcess = "Not Known -- (primary generator + Rayleigh)";
490          parent = 5;
491        }
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            }
497        // filling tuple
498       
499        //      G4cout<< particleType << G4endl;
500        G4int part = -1 ;
501        if (particleType == "gamma") part =1; 
502        if (particleType == "e-") part = 0;
503        if (particleType == "proton") part = 2;
504       
505
506        tupleFluo->fill(0,part);
507        tupleFluo->fill(1,particleEnergy);
508        tupleFluo->fill(2,momentum.theta());
509        tupleFluo->fill(3,momentum.phi());
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); 
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   
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.) ))" );
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.