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 x/y arrays 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.

<source lang="cpp">

1. include <TTree.h>
2. include <TFile.h>
3. include <TRandom3.h>
4. include <list>
5. 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();

} </source>