XrayFluoAnalysisManager.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: XrayFluoAnalysisManager.cc
// GEANT4 tag $Name: xray_fluo-V03-02-00
//
// Author: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
//
// History:
// -----------
// 20 Jul 2007 A.Mantero, code cleaning and update. Replaced histos with tuple
// 11 Jul 2003 A.Mantero, code cleaning / Plotter-XML addiction
//    Sep 2002 A.Mantero, AIDA3.0 Migration
// 06 Dec 2001 A.Pfeiffer updated for singleton
// 30 Nov 2001 Guy Barrand : migrate to AIDA-2.2.
// 28 Nov 2001 Elena Guardincerri     Created
//
// -------------------------------------------------------------------
#include "g4root.hh"

#include "G4VProcess.hh"
#include "XrayFluoAnalysisManager.hh"
#include "G4Step.hh"
#include "XrayFluoDetectorConstruction.hh"
#include "G4VPhysicalVolume.hh"
#include "G4Gamma.hh"
#include "G4Electron.hh"
#include "G4Proton.hh"
#include "G4SystemOfUnits.hh"


XrayFluoAnalysisManager* XrayFluoAnalysisManager::instance = 0;

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

namespace { 
  //Mutex to acquire access to singleton instance check/creation
  G4Mutex instanceMutex = G4MUTEX_INITIALIZER;
  //Mutex to acquire accss to histograms creation/access
  //It is also used to control all operations related to histos 
  //File writing and check analysis
  G4Mutex dataManipulationMutex = G4MUTEX_INITIALIZER;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

XrayFluoAnalysisManager::XrayFluoAnalysisManager()
  :outputFileName("xrayfluo"), phaseSpaceFlag(false), physicFlag (false), 
   gunParticleEnergies(0), gunParticleTypes(0)
{
  dataLoaded = false;
  fParticleEnergyAndTypeIndex = 0;

  //creating the messenger
  analisysMessenger = new XrayFluoAnalysisMessenger(this);
  
  //Instantiate the analysis manager
  G4AnalysisManager::Instance();

  G4cout << "XrayFluoAnalysisManager created" << G4endl;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

XrayFluoAnalysisManager::~XrayFluoAnalysisManager() 
{
  if ( gunParticleEnergies ) delete gunParticleEnergies;
  gunParticleEnergies = 0;
  if ( gunParticleTypes ) delete gunParticleTypes;
  gunParticleTypes = 0;

  delete instance;
  instance = 0;
  
  delete G4AnalysisManager::Instance();  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

XrayFluoAnalysisManager* XrayFluoAnalysisManager::getInstance()

{
  G4AutoLock l(&instanceMutex);
  if (instance == 0) {instance = new XrayFluoAnalysisManager;}
  return instance;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void XrayFluoAnalysisManager::book()
{
  G4AutoLock l(&dataManipulationMutex);
  // Get analysis manager
  G4AnalysisManager* man = G4AnalysisManager::Instance();
  // Open an output file
  man->OpenFile(outputFileName);
  man->SetVerboseLevel(1);
  man->SetFirstHistoId(1);
  man->SetFirstNtupleId(1);

  G4cout << "Open output file: " << outputFileName << G4endl;

  if (phaseSpaceFlag) 
    {   
      // Book output Tuple (ID = 1)
      man->CreateNtuple("101","OutputNTuple");
      man->CreateNtupleIColumn("particle");      
      man->CreateNtupleDColumn("energies");
      man->CreateNtupleDColumn("momentumTheta");
      man->CreateNtupleDColumn("momentumPhi");
      man->CreateNtupleIColumn("processes");
      man->CreateNtupleIColumn("material");
      man->CreateNtupleIColumn("origin");
      man->CreateNtupleDColumn("depth");
      man->FinishNtuple();
      G4cout << "Created ntuple for phase space" << G4endl;
    }
  else {
    // Book histograms
    man->CreateH1("h1","Energy Deposit", 500,0.,10.); //20eV def.
    man->CreateH1("h2","Gamma born in the sample", 100,0.,10.);
    man->CreateH1("h3","Electrons  born in the sample", 100,0.,10.);
    man->CreateH1("h4","Gammas leaving the sample", 300,0.,10.);
    man->CreateH1("h5","Electrons leaving the sample ",200000 ,0.,10.0); // .05 eV def.
    man->CreateH1("h6","Gammas reaching the detector", 100,0.,10.);
    man->CreateH1("h7","Spectrum of the incident particles", 100,0.,10.);
    man->CreateH1("h8","Protons reaching the detector", 100,0.,10.);
    man->CreateH1("h9","Protons leaving the sample", 100,0.,10.);
    G4cout << "Created histos" << G4endl;
  }
} 

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void XrayFluoAnalysisManager::LoadGunData(G4String fileName, G4bool raileighFlag) 
{
  G4AutoLock l(&dataManipulationMutex);

  //Do not allow more than one thread to trigger the file reading
  if (dataLoaded)
    return;

  // Get analysis reader manager
  G4AnalysisReader* analysisReader = G4AnalysisReader::Instance();
  analysisReader->SetVerboseLevel(1);

  //This is for testing purposes
  G4int ntupleId = analysisReader->GetNtuple("101",fileName);
  G4cout << "Got ntupleId: " << ntupleId << G4endl;

  gunParticleEnergies = new std::vector<G4double>;
  gunParticleTypes = new std::vector<G4String>;

  G4int particleID, processID;
  G4double energy;
  analysisReader->SetNtupleIColumn("processes", processID);
  analysisReader->SetNtupleDColumn("energies", energy);
  analysisReader->SetNtupleIColumn("particles", particleID);

  while (analysisReader->GetNtupleRow() ) 
    {
      if (raileighFlag ^ (!raileighFlag && (processID == 1 || 
                        processID == 2)) )  
    {
      gunParticleEnergies->push_back(energy*MeV);
      if (particleID == 1)
        gunParticleTypes->push_back("gamma");
      else if (particleID == 0)
        gunParticleTypes->push_back("e-");
    }
    }

  G4cout << "Maximum mumber of events: "<< gunParticleEnergies->size() <<G4endl;  

  dataLoaded= true;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

const std::pair<G4double,G4String> XrayFluoAnalysisManager::GetEmittedParticleEnergyAndType() 
{
  G4AutoLock l(&dataManipulationMutex);
  std::pair<G4double,G4String> result;
  
  if(fParticleEnergyAndTypeIndex < (G4int) gunParticleEnergies->size())
    {
      G4double energy = gunParticleEnergies->at(fParticleEnergyAndTypeIndex);
      G4String name = gunParticleTypes->at(fParticleEnergyAndTypeIndex);
      result.first = energy;
      result.second = name;
    }

  fParticleEnergyAndTypeIndex++;
  return result;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....


void XrayFluoAnalysisManager::finish()
{ 
  G4AutoLock l(&dataManipulationMutex);
  G4cout << "Going to save histograms" << G4endl;
  // Save histograms
  G4AnalysisManager* man = G4AnalysisManager::Instance();
  man->Write();
  man->CloseFile();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....


void XrayFluoAnalysisManager::SetPhysicFlag(G4bool val)
{
  physicFlag = val;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void XrayFluoAnalysisManager::analyseStepping(const G4Step* aStep)
{
  G4AutoLock l(&dataManipulationMutex);
  G4AnalysisManager* man = G4AnalysisManager::Instance();

  if (phaseSpaceFlag){       
    G4ParticleDefinition* particleType= 0;
    G4String parentProcess="";
    G4ThreeVector momentum(0.,0.,0.);
    G4double particleEnergy=0;
    G4String sampleMaterial="";
    G4double particleDepth=0;
    G4int isBornInTheSample=0;
    XrayFluoDetectorConstruction* detector = XrayFluoDetectorConstruction::GetInstance();    

    // Select volume from which the step starts
    if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
      G4ThreeVector creationPos = aStep->GetTrack()->GetVertexPosition();
      // qua serve ancora una selezione: deve essere interno al sample
      //codice "a mano" per controllare il volume di orgine della traccia
     
      G4VPhysicalVolume* creationPosVolume = detector->GetGeometryNavigator()->LocateGlobalPointAndSetup(creationPos);

      // if physicflag is on, I record all the data of all the particle born in the sample. 
      // If it is off (but phase space production is on) I record data 
      // only for particles creted inside the sample and exiting it.  
      if (physicFlag ^ (!physicFlag && 
            (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) 
            )) 
    {   
      // extracting information needed
      particleType = aStep->GetTrack()->GetDynamicParticle()->GetDefinition();
      momentum = aStep->GetTrack()->GetDynamicParticle()->GetMomentum();
      particleEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();
      if (creationPosVolume->GetName() == "Sample" ) {
        isBornInTheSample = 1;
        particleDepth = creationPosVolume->GetLogicalVolume()->GetSolid()
          ->DistanceToOut(creationPos, G4ThreeVector(0,0,-1));
      }
      else {
        particleDepth = -1;
      }
      // metodo per ottenere la profondita' "a mano"
      //    particleDepth = std::abs(creationPos.z() - detector->GetSampleIlluminatedFaceCoord()); 

      G4int parent=0;
      if(aStep->GetTrack()->GetCreatorProcess())
        {
          parentProcess = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
          parent = 5;
          // hack for HBK
          if (parentProcess == "") parent = 0;
          if (parentProcess == "ioni") parent = 1;
          if (parentProcess == "LowEnPhotoElec") parent = 2;
          if (parentProcess == "Transportation") parent = 3;// should never happen
          if (parentProcess == "initStep") parent = 4;// should never happen
        }   
      else {
        parentProcess = "Not Known -- (primary generator + Rayleigh)";
        parent = 5;
      }
      G4int sampleMat=0;
      if(aStep->GetTrack()){
        sampleMaterial = aStep->GetTrack()->GetMaterial()->GetName();
        if (sampleMaterial == ("Dolorite" || "Anorthosite" || "Mars1" || "IceBasalt" || "HPGe")) sampleMat=1;
      }
    
      G4int part = -1 ;
      if (particleType == G4Gamma::Definition()) part =1; 
      if (particleType == G4Electron::Definition()) part = 0;
      if (particleType == G4Proton::Definition()) part = 2;
      
      //Fill ntuple  
      man->FillNtupleIColumn(1,0, part);
      man->FillNtupleDColumn(1,1,particleEnergy);
      man->FillNtupleDColumn(1,2,momentum.theta());
      man->FillNtupleDColumn(1,3,momentum.phi());
      man->FillNtupleIColumn(1,4,parent); 
      man->FillNtupleIColumn(1,5,sampleMat);
      man->FillNtupleIColumn(1,6,isBornInTheSample);
      man->FillNtupleDColumn(1,7,particleDepth);
      man->AddNtupleRow(1);
      
    }
    }
  }
    
  // Normal behaviour, without creation of phase space
  else 
    {
   
      // Filling the histograms with data, passing thru stepping.

      // Select volume from wich the step starts
      if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
    // Select volume from wich the step ends
    if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" ) { 
      // Select the particle type
      if ((aStep->GetTrack()->GetDynamicParticle()
           ->GetDefinition()) == G4Gamma::Definition() )
    
        {         
          G4double gammaLeavingSample = 
        aStep->GetPreStepPoint()->GetKineticEnergy();
          man->FillH1(4,gammaLeavingSample/keV);

        }   
      else if ((aStep->GetTrack()->GetDynamicParticle()
            ->GetDefinition()) == G4Electron::Definition() ) 
        {
          G4double eleLeavingSample = 
        aStep->GetPreStepPoint()->GetKineticEnergy();
          man->FillH1(5,eleLeavingSample/keV);
        }
      else if ((aStep->GetTrack()->GetDynamicParticle()
            ->GetDefinition()) == G4Proton::Definition() )
        {
          G4double 
        protonsLeavSam = aStep->GetPreStepPoint()->GetKineticEnergy();
          man->FillH1(9,protonsLeavSam/keV);
        }
    }   
      }
    
      if((aStep->GetTrack()->GetDynamicParticle()
      ->GetDefinition()) == G4Gamma::Definition() )
        {
      if(aStep->GetTrack()->GetCurrentStepNumber() == 1)
        {         
          if(aStep->GetTrack()->GetParentID() != 0)
        {   
          if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
            {
              G4double gammaBornInSample = 
            aStep->GetPreStepPoint()->GetKineticEnergy();
              man->FillH1(2,gammaBornInSample/keV);
            }
        }
        }
    }    
      else if ((aStep->GetTrack()->GetDynamicParticle()
        ->GetDefinition() ) == G4Electron::Definition() )
    {
      if(aStep->GetTrack()->GetCurrentStepNumber() == 1)
        {         
          if(aStep->GetTrack()->GetParentID() != 0)
        {   
          if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
            {
              G4double eleBornInSample = 
            aStep->GetPreStepPoint()->GetKineticEnergy();
              man->FillH1(3,eleBornInSample/keV);
            }
        }
        }
    }
      
  
      if(aStep->GetTrack()->GetNextVolume())
    {

      
      if(aStep->GetTrack()->GetNextVolume()->GetName() == "HPGeDetector")
        
        { 
          if ((aStep->GetTrack()->GetDynamicParticle()
           ->GetDefinition()) == G4Gamma::Definition() ) 
        {

          G4double gammaAtTheDetPre = 
            aStep->GetPreStepPoint()->GetKineticEnergy();
          man->FillH1(6,gammaAtTheDetPre/keV);
        }
          else if ((aStep->GetTrack()->GetDynamicParticle()
            ->GetDefinition() ) == G4Proton::Definition() ) 
        {
          G4double protonsAtTheDetPre = 
            aStep->GetPreStepPoint()->GetKineticEnergy();
          man->FillH1(8,protonsAtTheDetPre);
        }
        }
    }
    }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void XrayFluoAnalysisManager::analyseEnergyDep(G4double energyDep)
{
  G4AutoLock l(&dataManipulationMutex);
  // Filling of Energy Deposition, called by XrayFluoEventAction  
  G4AnalysisManager* man = G4AnalysisManager::Instance();

  if (!phaseSpaceFlag)
    man->FillH1(1,energyDep/keV);
  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void XrayFluoAnalysisManager::analysePrimaryGenerator(G4double energy)
{
  G4AutoLock l(&dataManipulationMutex);
  // Filling of energy spectrum histogram of the primary generator
  G4AnalysisManager* man = G4AnalysisManager::Instance();
  if (!phaseSpaceFlag)
    man->FillH1(7,energy/keV);
  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void XrayFluoAnalysisManager::SetOutputFileName(G4String newName)
{
  
  outputFileName = newName;
}