Count the time and bounce in Integrating Sphere

Hello to everyone,

I am new to geant4 and am trying to create a simulation of an integrating sphere with liquid (currently water) and a detector inside it.

The structure of the geometry is almost fully completed. My idea is to fire a certain number of particles (fParticleGun already implemented) from the edge of one of the two conduits and count the number of bounces of the optical photon on the surface of the sphere, the number of photons reaching the yellow detector and the time it takes for the photon to reach that detector.

As far as the two surfaces (sphere and detector) are concerned, I think I have already made them sensitive:

  G4MultiFunctionalDetector* SferaIntegratrice_SD = new G4MultiFunctionalDetector("SferaIntegratrice_SD");
  G4SDManager::GetSDMpointer()->AddNewDetector(SferaIntegratrice_SD);
  //G4VPrimitiveScorer* primitiv1 = new G4PSDoseDeposit("dose");
  //SferaIntegratrice_SD->RegisterPrimitive(primitiv1);
  SetSensitiveDetector("SferaIntegratriceLV",SferaIntegratrice_SD);
  
  
  G4MultiFunctionalDetector* SiPM_SD = new G4MultiFunctionalDetector("SiPM_SD");
  G4SDManager::GetSDMpointer()->AddNewDetector(SiPM_SD);
  //G4VPrimitiveScorer* primitiv2 = new G4PSDoseDeposit("dose");
  //SiPM_SD->RegisterPrimitive(primitiv2);
  SetSensitiveDetector("SiPMLV",SiPM_SD);

Correct me if I am wrong. To count the number of bounces, for simplicity’s sake let’s consider firing a single optical photon, I think we need to use the class G4PSSphereSurfaceCurrent.hh (reading from the developer’s manual). I couldn’t find any examples on it unfortunately…

So I ask, how can I count the number of bounces, the number of photons that reach the yellow detector and then disappear, and the time it takes for the photon to reach that detector once it has been fired?

I also upload the “DetectorConstruction.cc” DetectorConstruction.cc (12.3 KB)

Thank you for your help.

1 Like

I was able to implement the right scorer:

  G4MultiFunctionalDetector* SferaIntegratrice_SD = new G4MultiFunctionalDetector("SferaIntegratrice_SD");
  G4SDManager::GetSDMpointer()->AddNewDetector(SferaIntegratrice_SD);
  SetSensitiveDetector("SferaIntegratriceLV",SferaIntegratrice_SD);
  G4VPrimitiveScorer* primitiv1 = new G4PSSphereSurfaceCurrent("number_of_bounce", fCurrent_Out);
  SferaIntegratrice_SD->RegisterPrimitive(primitiv1);
  
  G4cout<<"Created G4MultiFunctionalDetector named "

        <<SferaIntegratrice_SD->GetName()<<", and a G4PSSphereSurfaceCurrent scorer " <<"named "<<primitiv1->GetName()<<G4endl;
  
  
 
  G4MultiFunctionalDetector* SiPM_SD = new G4MultiFunctionalDetector("SiPM_SD");
  G4SDManager::GetSDMpointer()->AddNewDetector(SiPM_SD);
  SetSensitiveDetector("SiPMLV",SiPM_SD);
  G4VPrimitiveScorer* primitiv2 = new G4PSFlatSurfaceFlux("number_photon_SiPM", fFlux_In);
  SiPM_SD->RegisterPrimitive(primitiv2);
  
  
  G4cout<<"Created G4MultiFunctionalDetector named "

        <<SiPM_SD->GetName()<<", and a G4PSFlatSurfaceFlux scorer " <<"named "<<primitiv2->GetName()<<G4endl; 
  

Now the question is: how can i obtain the information about primitiv1 and primitiv2 after each run?

alternatively, you can attach the scorers via commands/macro. there you can also dump results after each run. think you will find suitable examples for that

Edit: sorry, missed that it is a sphere… not so sure about that then :frowning:

How can I do that? What changes if they are inside a sphere?

I tried with this command:

/score/create/boxMesh OpticalBoxMesh
/score/mesh/boxSize 10 10 10 cm
/score/mesh/nBin 1 1 1
/score/quantity/flatSurfaceFlux number_photon_SiPM
/run/beamOn 1000
/score/quantity/flatSurfaceFlux OpticalBoxMesh number_photon_SiPM Risult_Sim.txt
/score/close

but i get this error: G4WT0 > number_photon_SiPM: integer expected.

Of course the number that i expected is an integer but i don’t understand why this not happen…

/score/close before /run/beamOn 1000
after the run, you can dump:
/score/dumpAllQuantitiesToFile scorer filename.csv

or you make use of the histogram filler: Command-based scoring — Book For Application Developers 11.0 documentation

but as mentioned, the SphereSurfaceCurrent is not in the list of available primitive scorers for command based scoring, as you will see when you study that page :frowning:

an idea could be to define an additional hollow sphere, and attach the volumeFlux scorer to it?

1 Like

Mmm… Maybe i can… maybe, I’m not sure.Thanks for the reply.

By the way, I made this mac file:

/run/numberOfThreads 4
# Initialize kernel
/run/initialize
#
/control/verbose 2
/tracking/verbose 2
#
/score/create/boxMesh OpticalBoxMesh
/score/mesh/boxSize 10 10 10 cm
/score/mesh/nBin 1 1 1
/score/quantity/flatSurfaceFlux number_photon_SiPM
/run/beamOn 1000
/score/close
#/score/quantity/flatSurfaceFlux OpticalBoxMesh number_photon_SiPM Risult_Sim.txt
/score/dumpAllQuantitiesToFile OpticalBoxMesh Result_Sim.csv

and i get this file Result_Sim.csv:

# mesh name: OpticalBoxMesh
# primitive scorer name: number_photon_SiPM
# iX, iY, iZ, total(value) [percm2], total(val^2), entry
0,0,0,0,0,0

Seems that no info was stored…

1 Like

it seems that you do not close the scorer before you switch on the beam :wink:

I did it but nothing change…

# Macro file of "exampleB3.cc"
# Change the default number of workers (in multi-threading mode) 
/run/numberOfThreads 4
# Initialize kernel
/run/initialize
#
/control/verbose 2
/tracking/verbose 2
#
/score/create/boxMesh OpticalBoxMesh
/score/mesh/boxSize 10 10 10 cm
/score/mesh/nBin 1 1 1
/score/quantity/flatSurfaceFlux number_photon_SiPM
/score/close
/run/beamOn 5000
#/score/quantity/flatSurfaceFlux OpticalBoxMesh number_photon_SiPM Risult_Sim.txt
/score/dumpAllQuantitiesToFile OpticalBoxMesh Result_Sim.csv

Why?

1 Like

flatsurfaceflux is only measured in the -z plane of the cube. maybe there are no photons passing through? try to rotate the mesh accordingly

https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Detector/commandScore.html?highlight=flatsurfaceflux#id9

1 Like

Thanks @weller .

I was able to resolve part of my problem in this way:

void SteppingAction::UserSteppingAction(const G4Step* step)
{
// Collect energy and track length step by step

  // get volume of the current step
  auto volume_pre = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
  
  auto volume_post = step->GetPostStepPoint()->GetPhysicalVolume()->GetName();

  G4double stepLength = 0.;
  G4double stepTime = 0.;
  
  if ( step->GetTrack()->GetTrackID() == 1. ) {
    stepLength = step->GetStepLength()/CLHEP::cm;
    stepTime = step->GetTrack()->GetGlobalTime()/CLHEP::ns;
    
    if ( stepLength == 0 ) {
     
     ++n_bounce; 
     
    }
    
    if (step->GetTrack()->GetTrackStatus() == fStopAndKill) {
    
    	if (volume_post == "SiPMPhys"){
    
    G4cout << "Particle is died after " << stepTime << " [ns], the # of bouce before to exit was " << n_bounce << "\t" << G4endl; 
    
    	}
    
    }
    
  } 

Now i want to save in a .root file the info about stepTime and n_bounce in order to do two different histo.

How can I do that? Can you help me?

like so?
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Analysis/managers.html?#histograms

1 Like

Works @weller! For example:

But i don’t understand because if I shot, hard coded, a number of photon > 1, the .root file is not generated:

/* //with this part of code, the simulation doesn't create the .root file
  G4int rate = 220;
  G4int time = 60;
  G4double energy = 45*10^(-6);
  G4int n_particle = rate*time*(energy/Energia_fotone);
*/
  
  fParticleGun  = new G4ParticleGun(1);

But works with the command: /run/beamOn 10000000. Why this happens?

    G4double energyA = 45*10^(-6);
    G4double energyB = 45e-6;
    G4cerr << energyA << G4endl << energyB << G4endl;

output:

-456
4.5e-05

¯\_(ツ)_/¯

1 Like

@samuelemillesoli In C++, ^ is the bitwise XOR operator. So in the energyA expression, the OP is computing `(45*10) ^ (-6) = 0x000001C2 xor 0xFFFFFFFA = 0xFFFFFE38 = -456.

C++ doesn’t have an exponential operator. You could (foolishly :slight_smile: ) use std::pow(10.,-6) for this, but 45e-6 is the correct notation.

1 Like

Thanks so much guys for the help and explanation! :slight_smile:

Is there a way to change the default run/beamOn 1 command to the command of my interest /run/beamOn n_particle ?

Moreover, by setting:


  G4int rate = 220;
  G4int time = 60;
  G4double energy = 45e-6;
  G4int n_particle = rate*time*(energy/photon_energy);

  
  fParticleGun = new G4ParticleGun(n_particle);

the programme does not generate the root file where the nTuples are contained, I do not understand the reason. If i put:

fParticleGun = new G4ParticleGun(n_particle);

and then i do the command /run/beamOn 5000 everything is fine!

what is the actual value of n_particle during execution (G4cerr << n_particle << G4endl;)? does it work if you put in fParticleGun = new G4ParticleGun(5000);
?

btw: the particle gun that way will produce 5000 identical primaries, at least same origin and same direction (not sure about energy), so if you expect a distribution in any of these properties there you need to be cautious.
also i think it is favorable to have multiple events with few primaries instead of few events with massively many primaries, as each event is independent from the others and can be computed in parallel. /run/beamOn 5000 on the other hand creates 5k events with the number of primaries you specify in the respective action

regarding the gui buttons: you find instruction to append your own button in the documentation, and also in the examples. not sure if it is as straight forward to adjust the default buttons, but please report back if you find out :slight_smile:

1 Like

I understand, thanks for the explanation!

Running the simulation I realised that the reflectivity parameter seems to have no weight, in fact if I enter 0.01 or 0.99, the number of bounces for the same number of photons sent is the same:

G4OpticalSurface* OpSurface = new G4OpticalSurface("Superfice_Sfera_Integratrice");
  
  
  G4LogicalBorderSurface* Surface = new
    G4LogicalBorderSurface("Superfice_Sfera_Integratrice", physSferaIntegratrice, physCilindro, OpSurface);


	OpSurface->SetType(dielectric_dielectric);
	OpSurface->SetModel(unified);
	OpSurface->SetFinish(groundfrontpainted);
	//OpSurface->SetSigmaAlpha(0.1);
  
 
  #include "Energy.hh"
  G4MaterialPropertiesTable* TEFLON_MPT = new G4MaterialPropertiesTable();
  
  std::vector<G4double> energy = {Energia_fotone*eV};
  std::vector<G4double> rindex_teflon = {1.30};
  std::vector<G4double> transmission = {0.9};
  std::vector<G4double> reflectivity = {0.1}; 
  
  //
  //properties of teflon that depend on energy
  //
  
  TEFLON_MPT->AddProperty("RINDEX", energy, rindex_teflon);
  TEFLON_MPT->AddProperty("TRANSMITTANCE", energy, transmission);
  TEFLON_MPT->AddProperty("REFLECTIVITY", energy, reflectivity);
  
  teflon_mat->SetMaterialPropertiesTable(TEFLON_MPT);
  
  OpSurface->SetMaterialPropertiesTable(TEFLON_MPT);
  
  G4MaterialPropertiesTable* WATER_MPT = new G4MaterialPropertiesTable();
  
  //std::vector<G4double> energy1 = {Energia_fotone*eV};
  std::vector<G4double> rindex_water = {1.34};
  std::vector<G4double> absorption1 = {344.8*cm};
  //std::vector<G4double> reflectivity1 = {0.}; 
  
  
  //
  //properties of water that depend on energy
  //
  
  WATER_MPT->AddProperty("RINDEX", energy, rindex_water);
  WATER_MPT->AddProperty("ABSLENGTH", energy, absorption1);
  //WATER_MPT->AddProperty("REFLECTIVITY", energy1, reflectivity1);
  
  water_mat->SetMaterialPropertiesTable(WATER_MPT);
  
  OpSurface->SetMaterialPropertiesTable(WATER_MPT);

Yet I followed the developer’s guide, following the same pattern in the code described in pag. 257 (listing 5.9).

Ok, I resolved it. I post here the right code:

  if ( lambda_fotone >= 250 && lambda_fotone < 269 ) { /*see Dicken - Accelerated UV weathering device based on integrating sphere technology*/

  std::vector<G4double> reflectivity = {0.965};
  std::vector<G4double> trasmittance = {0.035};
  std::vector<G4double> rindex_teflon = {1.30};

  TEFLON_MPT->AddProperty("RINDEX", energy, rindex_teflon);
  TEFLON_MPT->AddProperty("REFLECTIVITY", energy, reflectivity);
  TEFLON_MPT->AddProperty("TRANSMITTANCE", energy, reflectivity);
  

      } else if ( lambda_fotone >= 269 && lambda_fotone < 279 ) {

        std::vector<G4double> reflectivity = {0.975};
        std::vector<G4double> trasmittance = {0.025};
        std::vector<G4double> rindex_teflon = {1.30};

        TEFLON_MPT->AddProperty("RINDEX", energy, rindex_teflon);
        TEFLON_MPT->AddProperty("REFLECTIVITY", energy, reflectivity);
        TEFLON_MPT->AddProperty("TRANSMITTANCE", energy, reflectivity);

          } else if ( lambda_fotone >= 279 && lambda_fotone < 300s ) {

            std::vector<G4double> reflectivity = {0.985};
            std::vector<G4double> trasmittance = {0.015};
            std::vector<G4double> rindex_teflon = {1.30};

            TEFLON_MPT->AddProperty("RINDEX", energy, rindex_teflon);
            TEFLON_MPT->AddProperty("REFLECTIVITY", energy, reflectivity);
            TEFLON_MPT->AddProperty("TRANSMITTANCE", energy, reflectivity);

            } else { 

              std::vector<G4double> reflectivity = {0.99}; 
              std::vector<G4double> trasmittance = {0.01};
              std::vector<G4double> rindex_teflon = {1.30};

              TEFLON_MPT->AddProperty("RINDEX", energy, rindex_teflon);
              TEFLON_MPT->AddProperty("REFLECTIVITY", energy, reflectivity);
              TEFLON_MPT->AddProperty("TRANSMITTANCE", energy, reflectivity);

  }

Thanks for all the help!!

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.