How to calculate standard error for RE02 example

[hungbt1908]

Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!

Geant4 Version: 10.07.p4
Operating System: Window 11
Compiler/Version:
CMake Version:

---

Dear all,
I’m trying to explore the “examples\extended\runAndEvent\RE02” example to learn the way to use G4MultiFunctionalDetector. I ran the example and the results are shown in the attached image.

My question is how to calculate the standard error for each voxel with a certain quality, the G4PSEnergyDeposit3D for example.

I saw in the B1 example show the way to calculate the rms but it has to know the eDep^2. I don’t know how to apply it for the RE02 example.

  // Compute dose = total energy deposit in a run and its variance
  //
  G4double edep  = fEdep.GetValue();
  G4double edep2 = fEdep2.GetValue();
  
  G4double rms = edep2 - edep*edep/nofEvents;
  if (rms > 0.) rms = std::sqrt(rms); else rms = 0.;  

  const B1DetectorConstruction* detectorConstruction
   = static_cast<const B1DetectorConstruction*>
     (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
  G4double mass = detectorConstruction->GetScoringVolume()->GetMass();
  G4double dose = edep/mass;
  G4double rmsDose = rms/mass;

[mcortesg]

Dear Bui Tien Hung,

Sorry for delay in response. You need to work on class RE02Run as follows:

1. Create this member variable: std::vector<G4THitsMap<G4double>*> fRunMap2;
2. Then, the idea is to copy whatever is done with fRunMap, so that fRunMap2 does the same. For instance, in RE02Run constructor you should add a line stating fRunMap2.push_back(new G4THitsMap<G4double>(detName,collectionName)) below the push_back done for fRunMap. And so on.
3. Of course, the key difference is in RecordEvent() function. Here, you instead add up the square of the quantities stored in *evtMap, i.e.: *fRunMap2[i] += (*evtMap)*(*evtMap);

Then, once you verify that the accumulation of the squared weights is working, you need to adapt the output accordingly, which is done in RunAction::EndOfRunAction() method.

I know that these instructions are brief and may not enough to have this working completely, but I hope they help you to understand the approach. If you encounter technical problems while doing this, please let us know.

Regards,
Miguel

[hungbt1908]

Dear Miguel,
I want to thank you for your response. I followed your instructions step by step but when I complied the source code, It showed me an error like “error C2676: binary ‘*’: ‘G4THitsMap’ does not define this operator or a conversion to a type acceptable to the predefined operator”.
I hope you can explain the error to me. Thank you again.
Best regards,
Hung

[mcortesg]

Oops, my fault, I forgot that detail, my sincere apologies. This happens because the operation ‘*’ is not defined in G4THitsMap, thus my instruction #3 needs to be prepared beforehand.

I am out of office for some days, but let me check if I can access my repo. I did this in the past, so I’ll get back to you with more accurate info as soon as I can.

Regards
Miguel

[mcortesg]

Hi again,

A possible workaround to get step #3 working is:

a. Create this method in Run class to produce a map storing the squares of quantities accumulated at end of event (please don’t forget to declare it in header file):

// Get the squares of all the quantities stored
// in the EvtMap in order to keep information
// about the uncertainties during all the run.
G4THitsMap<G4double>* 
Run::ComputeEvtHitsMap2( G4THitsMap<G4double>* evtHitsMap ) const
{
  G4THitsMap<G4double>* evtHitsMap2 = evtHitsMap;
  std::map<G4int, G4double*>* evtMap = evtHitsMap->GetMap();

  std::map<G4int, G4double*>::iterator it;
  for(it = evtMap->begin(); it != evtMap->end(); it++)  {
    G4int key = it->first;
    G4double evtX2 = (*(it->second))*(*(it->second));
    evtHitsMap2->set(key, evtX2);
  }

  return evtHitsMap2;
}

b. Then, go back to previous step #3, i.e. in method RecordEvent(), and write this instead of my previous hint:

//--- Create the evtMap2 to score squares.
G4THitsMap<G4double>* evtMap2 = ComputeEvtHitsMap2( evtMap );
// And sum it up to HitsMap of squares persistent in RUN
*fRunMap2[i] += *evtMap2;

I hope this works. Please check carefully variable names, as I adapted on the fly my naming criterion to this one.

Regards,
Miguel