vpgTLE-tt: new TLE of prompt gamma with time tagging (necessary for t…

…31 of becnhmark) (#662)

* ResetData missing for CT sinograms

* adaptation from gate 8.2 to new gate develop

* discard non null poststep protoninelastic

* refine proton energy before inelastic process by subtracting ionzation losses

* remove secondary protons

* weight tof histo with #PGs

* allocations ImageGamma(s) for Analog similar to vpgTLE

* debug PhaseSpace

* cleaning pTime init

* worldDefaultAir

* worldDefaultAir => set NULL TH1D

* proton=>PG emission time

* 500 t_bins

* Offset error in the 4th dimension of ImageOfHistograms

* tmp debug cout

* with pgTLE TimeBins messengers

* ToF compatible

* mhd user field

* 2ndary protons used

* bin Center correction if SourceOfPG

* correct 1st bin half pgTimeStep

* cleaning stage2

* remove Oreste DataTof class

.github/workflows/main.yml

@@ -278,7 +278,8 @@ jobs:
                         t27_Hits_2SD,
                         t28_Singles_2SD,
                         t29_optical_digi,
-                        t30_dna]
+                        t30_dna,
+                        t31_vpgTLE-tt]
 
     steps:
     - name: Checkout github repo

source/digits_hits/include/GatePhaseSpaceActor.hh

@@ -80,6 +80,8 @@ public:
 
   void SetIsTimeEnabled(bool b) { EnableTime = b; }
 
+  void SetIsIonTimeEnabled(bool b){EnableIonTime = b;}
+
   void SetIsLocalTimeEnabled(bool b) { EnableLocalTime = b; }
 
   void SetIsTimeFromBeginOfEventEnabled(bool b) { EnableTimeFromBeginOfEvent = b; }
@@ -180,6 +182,7 @@ protected:
   bool EnableProdProcess;
   bool EnableWeight;
   bool EnableTime;
+  bool EnableIonTime;
   bool EnableLocalTime;
   bool EnableMass;
   bool EnableSec;
@@ -211,7 +214,7 @@ protected:
   int bSpotID;
   bool bEnableCompact;
   bool bEnablePDGCode;
-  int bPDGCode;
+  long int bPDGCode;
   double trackLength;
 
   bool mMaskIsEnabled;
@@ -250,6 +253,7 @@ protected:
   float tOut;
   float tProd;
   double t; // t is either time or local time.
+  double pt; // Used for proton arrival time recovering from mTime in SourceOfPromptGamma
   double fTimeFromBeginOfEvent;
   double fBeginOfEventTime;
 

source/digits_hits/include/GatePhaseSpaceActorMessenger.hh

@@ -63,6 +63,7 @@ protected:
   G4UIcmdWithABool *pEnableParticleNameCmd;
   G4UIcmdWithABool *pEnableWeightCmd;
   G4UIcmdWithABool *pEnableTimeCmd;
+  G4UIcmdWithABool* pEnableIonTimeCmd;
   G4UIcmdWithABool *pEnableTimeFromBeginOfEventCmd;
   G4UIcmdWithABool *pEnableTrackLengthCmd;
   G4UIcmdWithABool *pEnableMassCmd;

source/digits_hits/include/GatePromptGammaAnalogActor.hh

@@ -31,13 +31,15 @@ public:
   virtual void UserPreTrackActionInVoxel(const int index, const G4Track* t);
   virtual void UserPostTrackActionInVoxel(const int index, const G4Track* t);
   virtual void UserSteppingActionInVoxel(const int index, const G4Step* step);
-
+  virtual void BeginOfEventAction(const G4Event * e);
+
   void SetInputDataFilename(std::string filename);
   virtual void SaveData();
   virtual void ResetData();
 
   void SetOutputCount(bool b) { mSetOutputCount = b; }  //output counts instead of yield
-
+  void SetTimeNbBins(int tbins) { data.SetTimeNbBins(tbins); } // nb of proton time bins
+
 protected:
   GatePromptGammaAnalogActor(G4String name, G4int depth=0);
   GatePromptGammaAnalogActorMessenger * pMessenger;
@@ -49,7 +51,14 @@ protected:
   bool mSetOutputCount;
   bool alreadyHere;
 
+  //helper functions
+  void SetTLEIoH(GateImageOfHistograms*&);
+  void SetTofIoH(GateImageOfHistograms*&);
+
   GateImageOfHistograms * mImageGamma;  //main output (yield)
+  GateImageOfHistograms * mImagetof;
+  double startEvtTime;
+  TH1D * pTime; // TH1D buffer histogram to add a new time value 
 
 };
 //-----------------------------------------------------------------------------

source/digits_hits/include/GatePromptGammaAnalogActorMessenger.hh

@@ -31,6 +31,7 @@ protected:
 
   G4UIcmdWithAString * pSetInputDataFileCmd;
   G4UIcmdWithABool * pSetOutputCountCmd;
+  G4UIcmdWithAnInteger * pTimeNbBinsCmd;
 };
 //-----------------------------------------------------------------------------
 

source/digits_hits/include/GatePromptGammaData.hh

@@ -33,13 +33,19 @@ public:
   void SetProtonNbBins(int x);
   void SetGammaNbBins(int x);
 
+  void SetTimeNbBins(int x);
+  void SetTimeTMax(double x);
+
   double GetProtonEMin() { return min_proton_energy; }
   double GetProtonEMax() { return max_proton_energy; }
   double GetGammaEMin()  { return min_gamma_energy; }
   double GetGammaEMax()  { return max_gamma_energy; }
   int GetProtonNbBins()  { return proton_bin; }
   int GetGammaNbBins()   { return gamma_bin; }
 
+  double GetTimeTMax() { return max_gamma_time;}
+  int GetTimeNbBins() { return time_bin;}
+
   void Initialize(std::string & filename, const G4Material * m);
   void InitializeMaterial(bool DebugOutputEnabled);
   void Read(std::string & filename);
@@ -75,6 +81,9 @@ protected:
   unsigned int proton_bin;
   unsigned int gamma_bin;
 
+  double max_gamma_time;
+  unsigned int time_bin;
+
   void SetCurrentPointerForThisElement(const G4Element * elem);
 
   //final outputs used to create tle source.

source/digits_hits/include/GatePromptGammaTLEActor.hh

@@ -42,6 +42,8 @@ public:
   virtual void SaveData();
   virtual void ResetData();
 
+  void SetTimeNbBins(int tbins) { data.SetTimeNbBins(tbins); } // nb of proton time bins
+
   void EnableDebugOutput(bool b) { mIsDebugOutputEnabled = b; }
   void EnableOutputMatch(bool b) { mIsOutputMatchEnabled = b; }
   //void EnableSysVarianceImage(bool b) { mIsSysVarianceImageEnabled = b; }
@@ -61,6 +63,7 @@ protected:
   //helper functions
   void SetTrackIoH(GateImageOfHistograms*&);
   void SetTLEIoH(GateImageOfHistograms*&);
+  void SetTofIoH(GateImageOfHistograms*&);
   GateVImageVolume* GetPhantom();
   void BuildVarianceOutput(); //converts trackl,tracklsq into mImageGamma and tlevar per voxel. Not used.
   //void BuildSysVarianceOutput(); //converts trackl into mImageGamma and tlesysvarv. Not used.
@@ -73,13 +76,18 @@ protected:
   GateImageOfHistograms * tracklsq;     //L_i^2. also intermediate output: track squared length per voxel per E_proton
 
   //output, calculated at end of simu
+  GateImageOfHistograms * mImagetof;
   GateImageOfHistograms * mImageGamma;  //oldstyle main output,
   GateImageOfHistograms * tle;          //main output (yield)
   GateImageOfHistograms * tlesysvar;    //systematic variance per voxel, per E_gamma. Not used.
   GateImageOfHistograms * tlevariance;  //uncertainty per voxel, per E_gamma. Not used.
 
   GateImageInt mLastHitEventImage;      //store eventID when last updated.
   int mCurrentEvent;                    //monitor event. TODO: not sure if necesary
+  int mCurrentIndex;                    //monitor the input and output of the particule in a voxel /** Modif Oreste **/
+  double startEvtTime;
+
+  TH1D * pTime; // TH1D buffer histogram to add a new time value 
 };
 //-----------------------------------------------------------------------------
 

source/digits_hits/include/GatePromptGammaTLEActorMessenger.hh

@@ -27,11 +27,12 @@ public:
 
 protected:
   void BuildCommands(G4String base);
-  GatePromptGammaTLEActor* pTLEActor;
+  GatePromptGammaTLEActor* pPGTLEActor;
 
   G4UIcmdWithAString * pSetInputDataFileCmd;
   G4UIcmdWithABool * pEnableDebugOutputCmd;
   G4UIcmdWithABool * pEnableOutputMatchCmd;
+  G4UIcmdWithAnInteger * pTimeNbBinsCmd;
 };
 //-----------------------------------------------------------------------------
 

source/digits_hits/src/GatePhaseSpaceActor.cc

@@ -51,6 +51,7 @@ GatePhaseSpaceActor::GatePhaseSpaceActor(G4String name, G4int depth) : GateVActo
     EnableProdProcess = true;
     EnableWeight = true;
     EnableTime = false;
+    EnableIonTime = false;
     EnableLocalTime = false;
     EnableTimeFromBeginOfEvent = false;
     EnableMass = true;
@@ -170,7 +171,7 @@ void GatePhaseSpaceActor::Construct()
             pIAEARecordType->iw = 1;
         if (EnableWeight)
             pIAEARecordType->iweight = 1;
-        if (EnableTime || EnableLocalTime)
+        if (EnableTime || EnableLocalTime || EnableIonTime)
         {
             GateWarning("'Time' is not available in IAEA phase space.");
         }
@@ -239,6 +240,8 @@ void GatePhaseSpaceActor::InitTree()
         mFile->write_variable("Weight", &w);
     if (EnableTime || EnableLocalTime)
         mFile->write_variable("Time", &t);
+    if (EnableIonTime)
+        mFile->write_variable("IonTime", &pt);
     if (EnableMass)
         mFile->write_variable("Mass", &m); // in MeV/c2
     if (EnableXPosition)
@@ -377,6 +380,7 @@ void GatePhaseSpaceActor::BeginOfEventAction(const G4Event *e)
 {
     // Set Primary Energy
     bPrimaryEnergy = e->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy(); // GetInitialEnergy oid.
+    if (EnableIonTime) pt = e->GetPrimaryVertex()->GetT0();
 
     // Store the application time of the event
     auto app = GateApplicationMgr::GetInstance();
@@ -815,10 +819,15 @@ void GatePhaseSpaceActor::UserSteppingAction(const GateVVolume *, const G4Step *
         pIAEARecordType->write_particle();
 
         pIAEAheader->update_counters(pIAEARecordType);
-    }
+    } 
     else
     {
         mFile->fill();
+//#define CORRELATION_TIME_ENERGY_SPACE
+#ifdef CORRELATION_TIME_ENERGY_SPACE
+	// std::cout << e << " " << tProd << std::endl;
+	std::cout << e << " " << tProd << " " << x << " " << y << " " << z << std::endl;
+#endif
     }
     mIsFirstStep = false;
 
@@ -853,7 +862,7 @@ void GatePhaseSpaceActor::SaveData()
         fclose(pIAEAheader->fheader);
         fclose(pIAEARecordType->p_file);
     }
-    else
+    else 
     {
         if (!this->mOverWriteFilesFlag)
         {

source/digits_hits/src/GatePhaseSpaceActorMessenger.cc

@@ -171,6 +171,12 @@ void GatePhaseSpaceActorMessenger::BuildCommands(G4String base)
     pEnableTimeCmd->SetGuidance(guidance);
     pEnableTimeCmd->SetParameterName("State", false);
 
+    bb = base+"/enableIonTime";
+    pEnableIonTimeCmd = new G4UIcmdWithABool(bb,this);
+    guidance = "Save the time of primary particles in the phase space file.";
+    pEnableIonTimeCmd->SetGuidance(guidance);
+    pEnableIonTimeCmd->SetParameterName("State",false);
+
     bb = base + "/enableTimeFromBeginOfEvent";
     pEnableTimeFromBeginOfEventCmd = new G4UIcmdWithABool(bb, this);
     guidance = "Save the time of particles since the begin of event time.";
@@ -346,6 +352,8 @@ void GatePhaseSpaceActorMessenger::SetNewValue(G4UIcommand *command, G4String pa
         pActor->SetIsWeightEnabled(pEnableWeightCmd->GetNewBoolValue(param));
     if (command == pEnableTimeCmd)
         pActor->SetIsTimeEnabled(pEnableTimeCmd->GetNewBoolValue(param));
+    if (command == pEnableIonTimeCmd)
+        pActor->SetIsIonTimeEnabled(pEnableIonTimeCmd->GetNewBoolValue(param));
     if (command == pEnableTimeFromBeginOfEventCmd)
         pActor->SetIsTimeFromBeginOfEventEnabled(pEnableTimeFromBeginOfEventCmd->GetNewBoolValue(param));
     if (command == pEnableTrackLengthCmd)
@@ -401,6 +409,8 @@ void GatePhaseSpaceActorMessenger::SetNewValue(G4UIcommand *command, G4String pa
 
     if (command == pEnableTOutCmd)
         pActor->SetIsTOutEnabled(pEnableTOutCmd->GetNewBoolValue(param));
+    if(command == pEnableTProdCmd)
+        pActor->SetIsTProdEnabled(pEnableTProdCmd->GetNewBoolValue(param));
     if (command == pEnableTProdCmd)
         pActor->SetIsTProdEnabled(pEnableTProdCmd->GetNewBoolValue(param));
     if (command == pUseMaskCmd)