Hi Tobias,
You can find in attachment at this mail an example of an including ROOT
Class in a Geant 4 development for a PET simulation.
See You
Seb.
---------------------------------------------
Sebastien JAN
Institut des Sciences Nucleaires - CNRS/IN2P3
Groupe ATLAS
53, Avenue des Martyrs
38026 Grenoble CEDEX
Tel : 04 76 28 41 87
Piece : 329
---------------------------------------------
#include "GepetoHisto.hh"
#include "GepetoHistoMessenger.hh"
// Include the CLHEP Histogram classes
//#include "CLHEP/Hist/Histogram.h"
//#include "CLHEP/Hist/Tuple.h"
//#include "CLHEP/Hist/HBookFile.h"
// Include files for ROOT.
#include "Rtypes.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TNtuple.h"
#include "TTree.h"
#include "TBranch.h"
#include "TRandom.h"
#include "globals.hh"
#include "G4Event.hh"
#include "G4Run.hh"
#include "G4Step.hh"
#include "RecorderBase.hh"
#include "G4ios.hh"
#include <iomanip.h>
#include "G4SDManager.hh"
#include "Randomize.hh"
#include "G4VVisManager.hh"
#include "G4UImanager.hh"
#include "G4VProcess.hh"
#include "G4HCofThisEvent.hh"
#include "GepetoDetectorMessenger.hh"
#include "GepetoHit.hh"
#include "CLHEP/Random/RanecuEngine.h"
GepetoHistograms::GepetoHistograms()
: histName("histfile"),DOI("YES"),Echant("YES")
{
runMessenger = new GepetoHistoMessenger(this);
saveRndm = 1;
hfile = NULL ;
sig511 = 0.0;
Rpet = 0.0;
DeltaX = 0.0;
DeltaY = 0.0;
DeltaZ = 0.0;
NumberRing1 = 0.;
EspaceEntreRing1 = 0.;
pet1Fov = 0.;
NumberCrystalTransaxial = 0;
NbrCoincTotale = 0.;
static TROOT rootBase("simple","Test of histogramming and I/O");
xpos1 = new Float_t[dimOfHitVector] ;
xpos2 = new Float_t[dimOfHitVector] ;
ypos1 = new Float_t[dimOfHitVector] ;
ypos2 = new Float_t[dimOfHitVector] ;
zpos1 = new Float_t[dimOfHitVector] ;
zpos2 = new Float_t[dimOfHitVector] ;
Edep1 = new Float_t[dimOfHitVector] ;
Edep2 = new Float_t[dimOfHitVector] ;
}
GepetoHistograms::~GepetoHistograms(){
delete runMessenger;
if (hfile != NULL) {
delete hfile;
}
delete tree;
delete[] xpos1 ;
delete[] xpos2 ;
delete[] ypos1 ;
delete[] ypos2 ;
delete[] zpos1 ;
delete[] zpos2 ;
delete[] Edep1 ;
delete[] Edep2 ;
}
void GepetoHistograms::RecordBeginOfRun(const G4Run* r)
{
hfile = new TFile(histName,"RECREATE");
// Arbre
tree = new TTree("T", "Tree de Simu");
tree->Branch("numberHits",&numberHits,"numberHits/I");
tree->Branch("Edep1",Edep1,"Edep1[numberHits]/F");
tree->Branch("xpos1",xpos1,"xpos1[numberHits]/F");
tree->Branch("ypos1",ypos1,"ypos1[numberHits]/F");
tree->Branch("zpos1",zpos1,"zpos1[numberHits]/F");
tree->Branch("numberHits1",&numberHits1,"numberHits1/I");
tree->Branch("Edep2",Edep2,"Edep2[numberHits1]/F");
tree->Branch("xpos2",xpos2,"xpos2[numberHits1]/F");
tree->Branch("ypos2",ypos2,"ypos2[numberHits1]/F");
tree->Branch("zpos2",zpos2,"zpos2[numberHits1]/F");
tree->Branch("Ndw1",&Ndw1,"Ndw1/I");
tree->Branch("Ndw2",&Ndw2,"Ndw2/I");
tree->Branch("Ee",&Ee,"Ee/I");
tree->Branch("xe_start",&xe_start,"xe_start/F");
tree->Branch("ye_start",&ye_start,"ye_start/F");
tree->Branch("ze_start",&ze_start,"ze_start/F");
tree->Branch("xe_stop",&xe_stop,"xe_stop/F");
tree->Branch("ye_stop",&ye_stop,"ye_stop/F");
tree->Branch("ze_stop",&ze_stop,"ze_stop/F");
tree->Branch("Erec1",&Erec1,"Erec1/F");
tree->Branch("Erec2",&Erec2,"Erec2/F");
tree->Branch("Xrec1",&Xrec1,"Xrec1/F");
tree->Branch("Xrec2",&Xrec2,"Xrec2/F");
tree->Branch("Yrec1",&Yrec1,"Yrec1/F");
tree->Branch("Yrec2",&Yrec2,"Yrec2/F");
tree->Branch("Zrec1",&Zrec1,"Zrec1/F");
tree->Branch("Zrec2",&Zrec2,"Zrec2/F");
tree->Branch("timeGlobal1",&timeLocal1,"timeLocal1/F");
tree->Branch("timeLocal1",&timeGlobal1,"timeGlobal1/F");
tree->Branch("timeGlobal2",&timeLocal2,"timeLocal2/F");
tree->Branch("timeLocal2",&timeGlobal2,"timeGlobal2/F");
tree->Branch("DeltaTemps",&DeltaTemps,"DeltaTemps/F");
tree->Branch("T_evts",&T_evts,"T_evts/F");
tree->Branch("N_t",&N_t,"N_t/F");
tree->Branch("A_t",&A_t,"A_t/F");
// save Rndm status
if (saveRndm > 0)
{ HepRandom::showEngineStatus();
HepRandom::saveEngineStatus("beginOfRun.rndm");
HepRandom::restoreEngineStatus("endOfRun.rndm");
}
G4UImanager* UI = G4UImanager::GetUIpointer();
G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
if(pVVisManager)
{
UI->ApplyCommand("/vis/scene/notifyHandlers");
}
// Initialisation Single
Nsingle1 = 0.;
Nsingle2 = 0.;
}
void GepetoHistograms::RecordEndOfRun(const G4Run* r)
{
hfile->Write();
hfile->Close();
// save Rndm status
if (saveRndm == 1)
{
HepRandom::showEngineStatus();
HepRandom::saveEngineStatus("endOfRun.rndm");
}
cout<<" ------------------"<<endl;
cout<<"| Activite Simulee | "<< Lambda*StartActivity <<" Bq"<< endl;
cout<<" ------------------"<<endl;
G4double DeltaT = 50.e-9;
cout<<" ------------------------"<<endl;
cout<<"| Coincidences Fortuites | "<< ceil(2*DeltaT*Lambda*StartActivity/(1+2*DeltaT*Lambda*StartActivity)*NbrCoincTotale) <<endl;
cout<<" ------------------------"<<endl;
cout<<" ------------------------"<<endl;
cout<<"| Coincidences Totales | "<< NbrCoincTotale <<endl;
cout<<" ------------------------"<<endl;
cout <<"$$$$$$ Single Evt gamma1 == " <<Nsingle1<<endl;
cout <<"$$$$$$ Single Evt gamma2 == " <<Nsingle2<<endl;
}
void GepetoHistograms::RecordBeginOfEvent(const G4Event* event)
{
for (G4int k = 0 ; k < dimOfHitVector ; k++)
{
Edep1[k] = 0.;
xpos1[k] = 0.;
ypos1[k] = 0.;
zpos1[k] = 0.;
Edep2[k] = 0.;
xpos2[k] = 0.;
ypos2[k] = 0.;
zpos2[k] = 0.;
}
// initialisation des variables Step/Step
Ndw1 = 0; // Nombre de compton dans le fantome pour gamma1 et gamma2
Ndw2 = 0;
Ee = 0;
E_positron = 0.; // Energie cinetique du positron au vertex
xe_start = 0.;
ye_start = 0.;
ze_start = 0.;
xe_stop = 0.;
ye_stop = 0.;
ze_stop = 0.;
// variables utilisees pour le calcul des parcourts de e+
libreX = 0.;
libreY = 0.;
libreZ = 0.;
longueur = 0.;
timeLocal1 = 0.;
timeGlobal1 = 0.;
timeLocal2 = 0.;
timeGlobal2 = 0.;
DeltaTemps = 0.;
// Variables reconstruites :
Erec1 = 0.;
Erec2 = 0.;
erec1 = 0.;
erec2 = 0.;
xrec1 = 0.;
xrec2 = 0.;
yrec1 = 0.;
yrec2 = 0.;
zrec1 = 0.;
zrec2 = 0.;
Xrec1 = 0.;
Xrec2 = 0.;
Yrec1 = 0.;
Yrec2 = 0.;
Zrec1 = 0.;
Zrec2 = 0.;
// variables de temps
Nevts = 0.;
Nt = 0.;
T_evts = 0.;
N_t = 0.;
A_t = 0.;
}
void GepetoHistograms::RecordEndOfEvent(const G4Event* event )
{
// Parametres temps de la simu
Nevts = event->GetEventID();
Nt = StartActivity - Nevts ; // Nt : nombre de noyaux restant
Lambda = log(2.0)/DemiVie ;
T_evts = 1/Lambda*log(StartActivity/Nt);
N_t = StartActivity*exp(-Lambda*T_evts);
A_t = Lambda*N_t;
//
G4SDManager* fSDM = G4SDManager::GetSDMpointer();
G4int collectionID = fSDM->GetCollectionID("CalCollection");
G4int collectionID1 = fSDM->GetCollectionID("CalCollection1");
G4HCofThisEvent* HCofEvent = event->GetHCofThisEvent();
G4HCofThisEvent* HCofEvent1 = event->GetHCofThisEvent();
GepetoHitsCollection* trackerHC =
(GepetoHitsCollection*) (HCofEvent->GetHC(collectionID));
GepetoHitsCollection* trackerHC1 =
(GepetoHitsCollection*) (HCofEvent1->GetHC(collectionID1));
numberHits = trackerHC->entries();
numberHits1 = trackerHC1->entries();
/////////// Parametres de Resolution /////////////
G4double Resolution= (sig511/100.0)*511.0;
G4double a = Resolution/(2.35*sqrt(511.0));
//////////////////////////////////////////////////
if( numberHits > dimOfHitVector ) { numberHits = dimOfHitVector ; }
for (G4int i = 0; i < numberHits ; i++)
{
GepetoHit* aHit = (*trackerHC)[i];
Edep1[i] = aHit->GetEdep() / keV;
G4double sigma1 = a*sqrt(Edep1[i]);
G4double DeltaEdep1 = RandGauss::shoot(Edep1[i],sigma1);
erec1 += DeltaEdep1 ;
xpos1[i] = aHit->GetPos().x() / cm;
ypos1[i] = aHit->GetPos().y() / cm;
zpos1[i] = aHit->GetPos().z() / cm;
xrec1 += xpos1[i]*Edep1[i];
yrec1 += ypos1[i]*Edep1[i];
zrec1 += zpos1[i]*Edep1[i];
}
if( numberHits1 > dimOfHitVector ) { numberHits1 = dimOfHitVector ; }
for (G4int j = 0; j < numberHits1; j++)
{
GepetoHit* aHit1 = (*trackerHC1)[j];
Edep2[j] = aHit1->GetEdep() / keV;
G4double sigma2 = a*sqrt(Edep2[j]);
G4double DeltaEdep2 = RandGauss::shoot(Edep2[j],sigma2);
erec2 += DeltaEdep2 ;
xpos2[j] = aHit1->GetPos().x() / cm;
ypos2[j] = aHit1->GetPos().y() / cm;
zpos2[j] = aHit1->GetPos().z() / cm;
xrec2 += xpos2[j]*Edep2[j];
yrec2 += ypos2[j]*Edep2[j];
zrec2 += zpos2[j]*Edep2[j];
}
// ****** On va mettre une condition pour ne remplir l'arbre que si il y a Hit pour gamma1 && gamma2 ************
// calcul du nombre de single
if (numberHits == 0. && numberHits1 != 0.)
Nsingle1++;
if (numberHits != 0. && numberHits1 == 0.)
Nsingle2++;
//////////////////////////////
if (numberHits != 0 && numberHits1 != 0)
{
G4double PI = 3.141592;
Erec1 = erec1 ;
Erec2 = erec2 ;
Rpet = (pet1Rmin*NumberModule1)/(2*PI);
if (DOI == "YES")
{
Xrec1 = xrec1 / Erec1;
Yrec1 = yrec1 / Erec1;
Zrec1 = zrec1 / Erec1;
Xrec2 = xrec2 / Erec2;
Yrec2 = yrec2 / Erec2;
Zrec2 = zrec2 / Erec2;
}
if (DOI == "NO" )
{
if (yrec1 >= 0.)
{
if (xrec1 >= 0.)
{
Xrec1 = Rpet*cos(atan(yrec1/xrec1));
Yrec1 = Rpet*sin(atan(yrec1/xrec1));
}
else
{
Xrec1 = -Rpet*cos(atan(yrec1/xrec1));
Yrec1 = Rpet*sin(atan(-yrec1/xrec1));
}
}
else
{
if (xrec1 >= 0.)
{
Xrec1 = Rpet*cos(atan(yrec1/xrec1));
Yrec1 = Rpet*sin(atan(yrec1/xrec1));
}
else
{
Xrec1 = -Rpet*cos(atan(yrec1/xrec1));
Yrec1 = Rpet*sin(atan(-yrec1/xrec1));
}
}
if (yrec2 >= 0.)
{
if (xrec2 >= 0.)
{
Xrec2 = Rpet*cos(atan(yrec2/xrec2));
Yrec2 = Rpet*sin(atan(yrec2/xrec2));
}
else
{
Xrec2 = -Rpet*cos(atan(yrec2/xrec2));
Yrec2 = Rpet*sin(atan(-yrec2/xrec2));
}
}
else
{
if (xrec2 >= 0.)
{
Xrec2 = Rpet*cos(atan(yrec2/xrec2));
Yrec2 = Rpet*sin(atan(yrec2/xrec2));
}
else
{
Xrec2 = -Rpet*cos(atan(yrec2/xrec2));
Yrec2 = Rpet*sin(atan(-yrec2/xrec2));
}
}
Zrec2 = zrec2 / Erec2;
Zrec1 = zrec1 / Erec1;
}
// ooOO0OOoo Calcul de Xrec et Yrec avec Echantillonnage ooOO0OOoo
if (Echant == "YES")
{
teta = 0.;
//// Calcul TetaEvt1
if (Xrec1 >= 0.)
{
if (Yrec1 >= 0.) TetaEvt1 = atan(Yrec1/Xrec1);
else
TetaEvt1 = 2*PI + atan(Yrec1/Xrec1);
}
else
{
TetaEvt1 = PI + atan(Yrec1/Xrec1);
}
//// Calcul TetaEvt2
if (Xrec2 >= 0.)
{
if (Yrec2 >= 0.) TetaEvt2 = atan(Yrec2/Xrec2);
else
TetaEvt2 = 2*PI + atan(Yrec2/Xrec2);
}
else
{
TetaEvt2 = PI + atan(Yrec2/Xrec2);
}
///////
for(G4int i = 0 ; i < NumberCrystalTransaxial ; i++)
{
if(TetaEvt1 >= teta && TetaEvt1 <= (teta + 2*PI/NumberCrystalTransaxial))
{
TetaEvt1 = teta + PI/NumberCrystalTransaxial;
}
if(TetaEvt2 >= teta && TetaEvt2 <= (teta + 2*PI/NumberCrystalTransaxial))
{
TetaEvt2 = teta + PI/NumberCrystalTransaxial;
}
teta = teta + 2*PI/NumberCrystalTransaxial;
}
Revt1 = sqrt(Xrec1*Xrec1 + Yrec1*Yrec1);
Revt2 = sqrt(Xrec2*Xrec2 + Yrec2*Yrec2);
Xrec1 = Revt1*cos(TetaEvt1);
Yrec1 = Revt1*sin(TetaEvt1);
Xrec2 = Revt2*cos(TetaEvt2);
Yrec2 = Revt2*sin(TetaEvt2);
// ooOO0OOoo Fin du Calcul de Xrec et Yrec avec Echantillonnage ooOO0OOoo
// ooOO0OOoo Calcul de Zrec1 avec Echantillonnage ooOO0OOoo
for(G4int j = 1 ; j < NumberRing1 + 1; j++)
{
if (Zrec1 >= -(NumberRing1*pet1Fov+(NumberRing1-1)*EspaceEntreRing1)/2 +(j-1)*pet1Fov +
(j-1)*EspaceEntreRing1 && Zrec1 <= -(NumberRing1*pet1Fov+(NumberRing1-1)*EspaceEntreRing1)/2 +j*pet1Fov +
(j-1)*EspaceEntreRing1)
{
Zrec1 = -(NumberRing1*pet1Fov+(NumberRing1-1)*EspaceEntreRing1)/2+(j-1)*pet1Fov+(j-1)*EspaceEntreRing1+pet1Fov/2 ;
goto Exit_Zrec1;
}
}
Exit_Zrec1:
// ooOO0OOoo Fin du Calcul de Zrec1 avec Echantillonnage ooOO0OOoo
// ooOO0OOoo Calcul de Zrec2 avec Echantillonnage ooOO0OOoo
for(G4int k = 1 ; k < NumberRing1 + 1; k++)
{
if (Zrec2 >= -(NumberRing1*pet1Fov+(NumberRing1-1)*EspaceEntreRing1)/2 +(k-1)*pet1Fov +
(k-1)*EspaceEntreRing1 && Zrec2 <= -(NumberRing1*pet1Fov+(NumberRing1-1)*EspaceEntreRing1)/2 +k*pet1Fov + (k-1)*EspaceEntreRing1)
{
Zrec2 = -(NumberRing1*pet1Fov+(NumberRing1-1)*EspaceEntreRing1)/2 +(k-1)*pet1Fov+(k-1)*EspaceEntreRing1+pet1Fov/2;
goto Exit_Zrec2;
}
}
// ooOO0OOoo Fin du Calcul de Zrec2 avec Echantillonnage ooOO0OOoo
}
Exit_Zrec2:
// ooOO0OOoo Remplissage du Tree ooOO0OOoo
// Condition d'exclusion d'evenements sur Zrec1 et Zrec2
if (abs(Zrec1) <= (NumberRing1*pet1Fov+(NumberRing1-1)*EspaceEntreRing1)/2 && abs(Zrec2) <= (NumberRing1*pet1Fov+(NumberRing1-1)*EspaceEntreRing1)/2)
{
// Calcul de la fenetre temporelle entre gamma1 et gamm2
DeltaTemps = abs(timeGlobal1 - timeGlobal2) ;
// Calcul du nombre de coincidences Totales -> evenement "double gamma" dans le detecteur
NbrCoincTotale ++ ;
// ooOO0OOoo Remplissage du Tree ooOO0OOoo
tree->Fill();
}
}
}
void GepetoHistograms::RecordStep(const G4Step* s)
{
G4Track* track = s->GetTrack();
G4double EdepStep = s->GetDeltaEnergy();
G4Material* material = track->GetMaterial();
G4String name = material->GetName();
G4int nbr_photon = track -> GetTrackID();
G4ParticleDefinition* particule = track ->GetDefinition();
G4String Nom = particule->GetParticleName();
G4StepPoint* sp = s -> GetPostStepPoint();
const G4VProcess* pro = sp -> GetProcessDefinedStep();
G4String nompro = pro -> GetProcessName();
G4ThreeVector vertg1,vertg2;
if (Nom == "e+") {
E_positron = track -> GetVertexKineticEnergy() / keV;
Ee = E_positron;
G4ThreeVector vertex_positron = track -> GetVertexPosition();
xe_start = vertex_positron.x() / cm ;
ye_start = vertex_positron.y() / cm ;
ze_start = vertex_positron.z() / cm ;
//cout << "xe = "<<xe<<endl;
//cout << "ye = "<<ye<<endl;
//cout << "ze = "<<ze<<endl;
//////////// parcourt des e+ dans le fantome/////////
delta = s->GetDeltaPosition();
delta_X = delta.x();
delta_Y = delta.y();
delta_Z = delta.z();
libreX += delta_X;
libreY += delta_Y;
libreZ += delta_Z;
leng = s -> GetStepLength();
longueur += leng;
vertg_e = track -> GetPosition();
xe_stop = vertg_e.x() / cm ;
ye_stop = vertg_e.y() / cm ; // coordonnees au photopic
ze_stop = vertg_e.z() / cm ;
}
//////////////////////////////////////////
// tout le reste pour les gammas !!!!!!!!!
//////////////////////////////////////////
if (Nom == "gamma" )
{
if (nbr_photon == 2) // 1er gamma
{
if (name == "Water")
{if (nompro == "LowEnCompton" || nompro == "LowEnRayleigh") Ndw1 += 1;} //compton dans l'eau
timeGlobal1 = track -> GetGlobalTime() / ns;
timeLocal1 = track -> GetLocalTime();
}
if (nbr_photon == 3) // 2eme gamma
{
if (name == "Water")
{if (nompro == "LowEnCompton" || nompro == "LowEnRayleigh") Ndw2 += 1;}//compton dans l'eau
timeGlobal2 = track -> GetGlobalTime() / ns;
timeLocal2 = track -> GetLocalTime();
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SethistName(G4String name)
{
histName = name;
G4cout << " hist file = " << histName << G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SetresoEner(G4double SigE)
{
sig511 = SigE;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SetresoX(G4double resox)
{
DeltaX = resox;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SetresoY(G4double resoy)
{
DeltaY = resoy;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SetresoZ(G4double resoz)
{
DeltaZ = resoz;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SetDOI(G4String doi)
{
DOI = doi;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SetEchant(G4String echant)
{
Echant = echant;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::Set1NumberRing(G4double ring)
{
NumberRing1 = ring;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::Set1AxialFov(G4double axialmod)
{
pet1Fov = axialmod;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::Set1InterRing(G4double espacering)
{
EspaceEntreRing1 = espacering;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SetNumberCrystalTransaxial(G4int val)
{
NumberCrystalTransaxial = val;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::Set1NumberModule(G4double val)
{
NumberModule1 = val;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::Set1Rmin(G4double val)
{
pet1Rmin = val;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
// Fonctions timing
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GepetoHistograms::SetDemiVie(G4double val)
{
DemiVie = val;
}
void GepetoHistograms::SetStartActivity(G4double val)
{
StartActivity = val;
}
#ifndef GepetoHistograms_h
#define GepetoHistograms_h 1
#include "RecorderBase.hh"
#include "G4UserRunAction.hh"
//#include "CLHEP/Hist/Histogram.h"
//#include "CLHEP/Hist/Tuple.h"
//#include "CLHEP/Hist/HBookFile.h"
// Include files for ROOT.
#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TNtuple.h"
#include "TTree.h"
#include "TBranch.h"
#include "TRandom.h"
#include "g4std/iostream"
#include "globals.hh"
#include "G4Step.hh"
#include "G4Event.hh"
#include "G4VProcess.hh"
class GepetoHit;
class GepetoHistoMessenger;
class G4Run;
static const int dimOfHitVector = 500 ;
class GepetoHistograms: public RecorderBase
{
public:
GepetoHistograms();
~GepetoHistograms();
void RecordBeginOfRun(const G4Run *);
void RecordEndOfRun(const G4Run *);
void RecordBeginOfEvent(const G4Event *);
void RecordEndOfEvent(const G4Event * );
void RecordStep(const G4Step *);
void SethistName(G4String name);
void SetDOI(G4String doi);
G4String GetDOI() {return DOI;};
void SetEchant(G4String echant);
G4String GetEchant() {return Echant;};
void SetresoEner(G4double SigE);
G4double GetresoEner() {return sig511;};
void SetresoX(G4double resox);
G4double GetresoX() {return DeltaX;};
void SetresoY(G4double resoy);
G4double GetresoY() {return DeltaY;};
void SetresoZ(G4double resoz);
G4double GetresoZ() {return DeltaZ;};
void SetRndmFreq(G4int val) {saveRndm = val;}
G4int GetRndmFreq() {return saveRndm;}
void Set1NumberRing(G4double val);
G4double Get1NumberRing() {return NumberRing1;};
void Set1Rmin(G4double val);
G4double Get1Rmin() {return pet1Rmin;};
void Set1NumberModule(G4double val);
G4double Get1NumberModule() {return NumberModule1;};
void Set1InterRing(G4double);
G4double Get1InterRing() {return EspaceEntreRing1;};
void SetDemiVie(G4double val);
G4double GetDemiVie() {return DemiVie;};
void SetStartActivity(G4double);
G4double GetStartActivity() {return StartActivity;};
void Set1AxialFov(G4double);
G4double Get1AxialFov() {return pet1Fov;};
void SetNumberCrystalTransaxial(G4int);
G4int GetNumberCrystalTransaxial() {return NumberCrystalTransaxial;};
private:
TFile * hfile;
G4String histName;
G4String DOI;
G4String Echant;
Float_t sig511;
Float_t Rpet;
Float_t DeltaX;
Float_t DeltaY;
Float_t DeltaZ;
Int_t Ndw1,Ndw2;
Int_t Ee;
G4double E_positron;
Float_t xe_start;
Float_t ye_start;
Float_t ze_start;
Float_t xe_stop;
Float_t ye_stop;
Float_t ze_stop;
Float_t timeLocal1;
Float_t timeGlobal1;
Float_t timeLocal2;
Float_t timeGlobal2;
Float_t DeltaTemps;
Float_t DemiVie;
Float_t StartActivity;
Float_t Nevts;
Float_t Nt;
Float_t T_evts;
Float_t Lambda;
Float_t A_t;
Float_t N_t;
Float_t NbrCoincTotale;
Float_t Lp,longueur,leng,delta_X,delta_Y,delta_Z,libreX,libreY,libreZ;
G4ThreeVector delta;
Float_t Erec1,Erec2,erec1,erec2;
Float_t xrec1,yrec1,zrec1;
Float_t xrec2,yrec2,zrec2;
Float_t Xrec1,Yrec1,Zrec1;
Float_t Xrec2,Yrec2,Zrec2;
Float_t * xpos1 ;
Float_t * ypos1 ;
Float_t * zpos1 ;
Float_t * xpos2 ;
Float_t * ypos2 ;
Float_t * zpos2 ;
Float_t * Edep1 ;
Float_t * Edep2 ;
G4ThreeVector vertg_e;
TTree *tree;
Int_t numberHits;
Int_t numberHits1;
G4int Ne;
GepetoHistoMessenger* runMessenger;
G4int saveRndm;
//////// Variables de l'echantillonnage
G4double teta;
G4double TetaEvt1,TetaEvt2;
G4double Revt1,Revt2;
G4int NumberCrystalTransaxial;
G4double NumberRing1;
G4double EspaceEntreRing1;
G4double pet1Fov;
G4double pet1Rmin;
G4double NumberModule1;
G4double Nsingle1,Nsingle2;
};
#endif
This archive was generated by hypermail 2b29 : Sat Jan 04 2003 - 23:50:53 MET