Cluster convolution for GATE - Revision history

Hpe090: Added code to take edep arguments

2019-04-01T11:56:02Z

Added code to take edep arguments

Hpe090 at 14:09, 25 February 2019

2019-02-25T14:09:03Z

← Older revision		Revision as of 14:09, 25 February 2019
Line 1:		Line 1:
	In her MSc project, Silje Grimstad made a filtered library of the clusters from the summer beamtest in Heidelberg.		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:		The result is ~15k Helium clusters stored as x/y arrays in a ROOT file:

	{\| class="wikitable"		{\| class="wikitable"

Hpe090 at 14:08, 25 February 2019

2019-02-25T14:08:54Z

← Older revision		Revision as of 14:08, 25 February 2019
Line 1:		Line 1:
	In her MSc project, Silje Grimstad made a filtered library of the clusters from the summer beamtest in Heidelberg.		In her MSc project, Silje Grimstad made a filtered library of the clusters from the summer beamtest in Heidelberg.
	The result is ~~~20k~~ Helium clusters stored as a x/y array in a ROOT file:		The result is ~15k Helium clusters stored as a x/y array in a ROOT file:

	{\| class="wikitable"		{\| class="wikitable"

Hpe090 at 14:05, 25 February 2019

2019-02-25T14:05:28Z

← Older revision		Revision as of 14:05, 25 February 2019
Line 41:		Line 41:
	#include <TTree.h>		#include <TTree.h>
	#include <TFile.h>		#include <TFile.h>
	#include <~~TRandom~~.h>		#include <TRandom3.h>
	#include <list>		#include <list>
	#include <vector>		#include <vector>
Line 62:		Line 62:
	Float_t pixel_size = 0.02929;		Float_t pixel_size = 0.02929;

	~~Float_t allZ = 0;~~		Float_t x,y,z,outX, outY, outZ;
	~~Int_t fillN = 0;~~
	Float_t x,y,z~~,edep~~,outX, outY, outZ;
	Int_t parentID, eventID, outEventID, outParentID, nClusters, timeInClockSpeed, PDGEncoding;		Int_t parentID, eventID, outEventID, outParentID, nClusters, timeInClockSpeed, PDGEncoding;
	~~Char_t processName[17];~~

	~~TCanvas *c = new TCanvas();~~
	~~TH2C *hOrig = new TH2C("hOrig", "Original hitmap",1024,-15,15,512,-7.5,7.5);~~
	~~TH2C *hNew = new TH2C("hNew", "New hitmap", 1024,-15,15,512,-7.5,7.5);~~
	~~Int_t palette[2] = {kWhite, kBlack};~~

	~~c->Divide(2,1);~~
	~~gStyle->SetPalette(2, palette);~~

	treeSimulationIn->SetBranchAddress("posX",&x);		treeSimulationIn->SetBranchAddress("posX",&x);
	treeSimulationIn->SetBranchAddress("posY",&y);		treeSimulationIn->SetBranchAddress("posY",&y);
	treeSimulationIn->SetBranchAddress("posZ",&z);		treeSimulationIn->SetBranchAddress("posZ",&z);
	~~treeSimulationIn->SetBranchAddress("edep",&edep);~~
	treeSimulationIn->SetBranchAddress("eventID",&eventID);		treeSimulationIn->SetBranchAddress("eventID",&eventID);
	treeSimulationIn->SetBranchAddress("parentID",&parentID);		treeSimulationIn->SetBranchAddress("parentID",&parentID);
Line 102:		Line 90:
	treeSimulationIn->GetEntry(i);		treeSimulationIn->GetEntry(i);

	timeInClockSpeed = 400 * eventID;		timeInClockSpeed = 400 * eventID; // 10 us readout speed
	outParentID = parentID;		outParentID = parentID; // =0 for a primary particle, >0 for secondaries
	outEventID = eventID;		outEventID = eventID;
	outZ = z;		outZ = z;

	~~if (z<2) hOrig->Fill(x,y);~~

	randomClusterIdx = gRand->Integer(nClusters);		randomClusterIdx = gRand->Integer(nClusters);
Line 113:		Line 99:

	idx_x = 0;		idx_x = 0;
	if (abs(PDGEncoding) == 2212) {		if (abs(PDGEncoding) == 2212) { // Only apply this for protons
	for (Int_t n : *hit_array) { // loop x		for (Int_t n : *hit_array) { // loop x
	for (int binPosPow = 0; binPosPow < 10; binPosPow++) { // loop y		for (int binPosPow = 0; binPosPow < 10; binPosPow++) { // loop y
Line 121:		Line 107:
	outY = y + (binPosPow - y_mean) * pixel_size;		outY = y + (binPosPow - y_mean) * pixel_size;
	treeSimulationOut->Fill();		treeSimulationOut->Fill();
	~~if (z<2) hNew->Fill(outX, outY);~~
	}		}
	}		}
Line 130:		Line 115:

	treeSimulationOut->Write();		treeSimulationOut->Write();

	~~c->cd(1);~~
	~~hOrig->Draw("COL");~~
	~~c->cd(2);~~
	~~hNew->Draw("COL");~~

	/*
	~~delete fSimulationIn;~~
	~~delete fSimulationOut;~~
	~~delete fCluster;~~
	*/
	}		}
	</source>		</source>

Hpe090: Made page

2019-02-25T14:03:02Z

Made page

New page

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

{| class="wikitable"
|-
! 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:

[[File:9_clusters.png|600px]]

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:

[[File:original_vs_convoluted_hitmap.PNG|1000px]]

The needed cluster database for this is located at [https://github.com/HelgeEgil/focal/blob/master/DTCToolkit/Data/ClusterSizes/ALPIDE/database_combined_clusters_filter.root GitHib].

<source lang="cpp">
#include <TTree.h>
#include <TFile.h>
#include <TRandom.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 allZ = 0;
Int_t fillN = 0;
Float_t x,y,z,edep,outX, outY, outZ;
Int_t parentID, eventID, outEventID, outParentID, nClusters, timeInClockSpeed, PDGEncoding;
Char_t processName[17];

TCanvas *c = new TCanvas();
TH2C *hOrig = new TH2C("hOrig", "Original hitmap",1024,-15,15,512,-7.5,7.5);
TH2C *hNew = new TH2C("hNew", "New hitmap", 1024,-15,15,512,-7.5,7.5);
Int_t palette[2] = {kWhite, kBlack};

c->Divide(2,1);
gStyle->SetPalette(2, palette);

treeSimulationIn->SetBranchAddress("posX",&x);
treeSimulationIn->SetBranchAddress("posY",&y);
treeSimulationIn->SetBranchAddress("posZ",&z);
treeSimulationIn->SetBranchAddress("edep",&edep);
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;
outParentID = parentID;
outEventID = eventID;
outZ = z;

if (z<2) hOrig->Fill(x,y);

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

idx_x = 0;
if (abs(PDGEncoding) == 2212) {
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();
if (z<2) hNew->Fill(outX, outY);
}
}
idx_x++;
}
}
}

treeSimulationOut->Write();

c->cd(1);
hOrig->Draw("COL");
c->cd(2);
hNew->Draw("COL");

/*
delete fSimulationIn;
delete fSimulationOut;
delete fCluster;
*/
}
</source>