Cluster convolution for GATE

In her MSc project, Silje Grimstad made a filtered library of the clusters from the summer beamtest in Heidelberg. The result is ~15k Helium clusters stored as a x/y array in a ROOT file:

x_mean	y_mean	size	hit_array	height	width	x_start	y_start	i_event	combined
float	float	int	list<int>	int	int	int	int	int	bool

A random sample of 9 clusters is shown below:

We want to convolute a GATE simulation with clusters from this library. For each hit in the MC data, we pick a random cluster from the library and draw it around the hit position. To do this, use the code at the bottom, with the following result:

The needed cluster database for this is located at GitHib.

#include <TTree.h>
#include <TFile.h>
#include <TRandom3.h>
#include <list>
#include <vector>

using namespace std;

void convert_root_file() {
   TFile         *fCluster = new TFile("database_combined_clusters_filter.root", "READ");
   TFile         *fSimulationIn = new TFile("simulationToBeClustered.root", "READ");
   TFile         *fSimulationOut = new TFile("simulationToBeClustered_clustered.root", "recreate");
   TTree         *treeCluster = (TTree*) fCluster->Get("database");
   TTree         *treeSimulationIn  = (TTree*) fSimulationIn->Get("Hits");
   TTree         *treeSimulationOut = new TTree("Hits", "Clustered GATE tree");

   Double_t       x_mean, y_mean;
   Int_t          randomClusterIdx, binPos, idx_x;
   list<Int_t>   *hit_array = 0;
   TBranch       *b_hit_array;
   TRandom3      *gRand = new TRandom3(0);
   Float_t        pixel_size = 0.02929;
   
   Float_t        x,y,z,outX, outY, outZ;
   Int_t          parentID, eventID, outEventID, outParentID, nClusters, timeInClockSpeed, PDGEncoding;

   treeSimulationIn->SetBranchAddress("posX",&x);
   treeSimulationIn->SetBranchAddress("posY",&y);
   treeSimulationIn->SetBranchAddress("posZ",&z);
   treeSimulationIn->SetBranchAddress("eventID",&eventID);
   treeSimulationIn->SetBranchAddress("parentID",&parentID);
   treeSimulationIn->SetBranchAddress("PDGEncoding", &PDGEncoding);

   treeSimulationOut->Branch("posX", &outX, "posX/F");
   treeSimulationOut->Branch("posY", &outY, "posY/F");
   treeSimulationOut->Branch("posZ", &outZ, "posZ/F");
   treeSimulationOut->Branch("eventID", &outEventID, "eventID/I");
   treeSimulationOut->Branch("parentID", &outParentID, "parentID/I");
   treeSimulationOut->Branch("clockTime", &timeInClockSpeed, "clockTime/I");

   treeCluster->SetBranchAddress("x_mean", &x_mean);
   treeCluster->SetBranchAddress("y_mean", &y_mean);
   treeCluster->SetBranchAddress("hit_array", &hit_array, &b_hit_array);

   nClusters = treeCluster->GetEntriesFast();
   randomClusterIdx = gRand->Integer(nClusters);
   treeCluster->GetEntry(randomClusterIdx);

   for (Int_t i=0, N = treeSimulationIn->GetEntries(); i<N; ++i) {
      treeSimulationIn->GetEntry(i);
      
      timeInClockSpeed = 400 * eventID; // 10 us readout speed
      outParentID = parentID; // =0 for a primary particle, >0 for secondaries
      outEventID = eventID;
      outZ = z;

      randomClusterIdx = gRand->Integer(nClusters);
      treeCluster->GetEntry(randomClusterIdx);

      idx_x = 0;
      if (abs(PDGEncoding) == 2212) { // Only apply this for protons
         for (Int_t n : *hit_array) { // loop x
            for (int binPosPow = 0; binPosPow < 10; binPosPow++) { // loop y
               binPos = pow(2, binPosPow);
               if (binPos & n) {
                  outX = x + (idx_x - x_mean) * pixel_size;
                  outY = y + (binPosPow - y_mean) * pixel_size;
                  treeSimulationOut->Fill();
               }
            }
            idx_x++;
         }
      }
   }

   treeSimulationOut->Write();
}