Merge pull request #67 from rest-for-physics/jgalan_analysis_upgrade

New observable boundingSize inside TRestGeant4AnalysisProcess

CMakeLists.txt

@@ -1,4 +1,4 @@
-set(LibraryVersion "1.8")
+set(LibraryVersion "1.9")
 add_definitions(-DLIBRARY_VERSION="${LibraryVersion}")
 
 if (${REST_DECAY0} MATCHES "ON")

inc/TRestGeant4Event.h

@@ -132,6 +132,8 @@ class TRestGeant4Event : public TRestEvent {
     void SetBoundaries(Double_t xMin, Double_t xMax, Double_t yMin, Double_t yMax, Double_t zMin,
                        Double_t zMax);
 
+    Double_t GetBoundingBoxSize();
+
     inline size_t GetNumberOfPrimaries() const { return fPrimaryParticleNames.size(); }
 
     inline TString GetPrimaryEventParticleName(size_t n = 0) const { return fPrimaryParticleNames[n]; }
@@ -238,6 +240,7 @@ class TRestGeant4Event : public TRestEvent {
 
     // Constructor
     TRestGeant4Event();
+
     // Destructor
     virtual ~TRestGeant4Event();
 

src/TRestGeant4AnalysisProcess.cxx

@@ -145,9 +145,12 @@
 ///
 /// * **thetaPrimary**: polar angle of the primary generated particle.
 /// * **phiPrimary**: azimuth angle of the primary generated particle.
-///
 /// * **energyPrimary**: energy of the primary event generated.
 ///
+/// * **boundingSize**: It stores a value with the event size calculated
+/// as the diagonal distance of a bounding box defined to contain all the
+/// hits that  produced an energy deposit.
+///
 /// The following code ilustrates the addition of a primary event
 /// observable.
 ///
@@ -417,8 +420,6 @@ TRestEvent* TRestGeant4AnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
     fInputG4Event = (TRestGeant4Event*)inputEvent;
     *fOutputG4Event = *((TRestGeant4Event*)inputEvent);
 
-    TString obsName;
-
     Double_t energy = fOutputG4Event->GetSensitiveVolumeEnergy();
 
     if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
@@ -456,19 +457,25 @@ TRestEvent* TRestGeant4AnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
     SetObservableValue((string) "energyPrimary", energyPrimary);
 
     Double_t energyTotal = fOutputG4Event->GetTotalDepositedEnergy();
-    obsName = this->GetName() + (TString) ".totalEdep";
     SetObservableValue((string) "totalEdep", energyTotal);
 
+    Double_t size = fOutputG4Event->GetBoundingBoxSize();
+    SetObservableValue((string) "boundingSize", size);
+
     // process names as named by Geant4
     // processes present here will be added to the list of observables which can be used to see if the event
     // contains the process of interest.
     vector<string> processNames = {"phot", "compt"};
-    for (const auto& processName : processNames) {
+    for (auto& processName : processNames) {
         Int_t containsProcess = 0;
         if (fOutputG4Event->ContainsProcess(fG4Metadata->GetGeant4PhysicsInfo().GetProcessID(processName))) {
             containsProcess = 1;
         }
-        SetObservableValue("ContainsProcess" + processName, containsProcess);
+
+        if (processName.size() > 0) {
+            processName[0] = toupper(processName[0]);
+            SetObservableValue("containsProcess" + processName, containsProcess);
+        }
     }
 
     /*
@@ -547,7 +554,7 @@ TRestEvent* TRestGeant4AnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
         cout << "----------------------------" << endl;
     }
 
-    // These cuts should be in another process or eliminated?!!
+    /// We should use here ApplyCut
     if (energy < fLowEnergyCut) return nullptr;
     if (fHighEnergyCut > 0 && energy > fHighEnergyCut) return nullptr;
 

src/TRestGeant4Event.cxx

@@ -281,6 +281,20 @@ void TRestGeant4Event::SetBoundaries(Double_t xMin, Double_t xMax, Double_t yMin
     fMaxZ = zMax;
 }
 
+///////////////////////////////////////////////
+/// \brief This method returns the event size as the size of the bounding box
+/// enclosing all hits.
+///
+Double_t TRestGeant4Event::GetBoundingBoxSize() {
+    SetBoundaries();
+
+    Double_t dX = fMaxX - fMinX;
+    Double_t dY = fMaxY - fMinY;
+    Double_t dZ = fMaxZ - fMinZ;
+
+    return TMath::Sqrt(dX * dX + dY * dY + dZ * dZ);
+}
+
 void TRestGeant4Event::SetBoundaries() {
     Double_t maxX = -1e10, minX = 1e10, maxZ = -1e10, minZ = 1e10, maxY = -1e10, minY = 1e10;
     Double_t minEnergy = 1e10, maxEnergy = -1e10;
@@ -297,6 +311,8 @@ void TRestGeant4Event::SetBoundaries() {
 
             Double_t en = hits.GetEnergy(nhit);
 
+            if (en <= 0) continue;
+
             if (x > maxX) maxX = x;
             if (x < minX) minX = x;
             if (y > maxY) maxY = y;