jpos2gate_converter

src/comptonscattering.cpp

@@ -322,7 +322,7 @@ void ComptonScattering::Scatter(Event* event, int index) const
             //if new_E is within limit -- smear, otherwise fill histogram with new_E
             if((new_E >= fSmearLowLimit_) && (new_E <= fSmearHighLimit_))
             {
-                double Esmear = gRandom->Gaus(new_E, sigmaE(E));
+                double Esmear = gRandom->Gaus(new_E, sigmaE(new_E));
                 fH_electron_E_blur_->Fill(Esmear);
                 event->SetEdepSmearOf(ii, Esmear);
             }

src/comptonscattering.h

@@ -52,7 +52,7 @@ class ComptonScattering
         float fSmearHighLimit_; //higher limit for phenomenologicly derived smearing effect
         static long double KleinNishina_(double* angle, double* energy); //Klein-Nishina function
         static long double KleinNishinaTheta_(double* angle, double* energy); //Klein-Nishina based theta PDF
-        double sigmaE(double E, double coeff=0.0444) const; //calculate std dev for the smearing effevt
+        double sigmaE(double E, double coeff=0.044) const; //calculate std dev for the smearing effevt
 
         static unsigned objectID_;
 

src/event.cpp

@@ -59,6 +59,50 @@ Event::Event(std::vector<TLorentzVector*>* emissionCoordinates, std::vector<TLor
     }
 }
 
+///
+/// \brief Event::Event Developed constructor allowing to create a complete Event.
+/// \param sourcePos Vector of TLorentzVector with entries corresponsing to X,Y,Z,T coordinates of emission point for a given particle.
+/// \param pos Vector of TLorentzVector with entries corresponsing to X,Y,Z,T coordinates of hit point for a given particle.
+/// \param momentum Vector of TLorentzVector with entries corresponsing to pX,pY,pZ,E coordinates of four-momentum for a given particle.
+/// \param phi Vector of azimuthal angles corresponding to given hit points positions.
+/// \param theta Vector of angular angles corresponding to given hit points positions.
+/// \param cutPassing Vector of bool values corresponding, true if particle left a signal in the detector.
+/// \param primary Vector of bool values corresponding, false if particle was scattered in a phantom.
+/// \param edep Vector of deposited energy by particles.
+/// \param edepSmear Vector of deposited energy by particles, smeared according to experimental formula.
+/// \param Id Particle ids.
+/// \param decayType Type of decay.
+///
+Event::Event(std::vector<TLorentzVector> &sourcePos, std::vector<TLorentzVector> &pos, std::vector<TLorentzVector> &momentum,\
+    std::vector<double> &phi, std::vector<double> &theta, std::vector<bool> &cutPassing, std::vector<bool> &primary,\
+    std::vector<double> &edep, std::vector<double> &edepSmear, long Id, int decayType)
+{
+    fCounter_++;
+    fId = Id;
+    fWeight_ = 1.0;
+    fDecayType_ = (DecayType)decayType;
+    fPassFlag_ = true;
+    fEmissionPoint_.resize(sourcePos.size());
+    std::copy(sourcePos.begin(), sourcePos.end(), fEmissionPoint_.begin());
+    fHitPoint_.resize(pos.size());
+    std::copy(pos.begin(), pos.end(), fHitPoint_.begin());
+    fFourMomentum_.resize(momentum.size());
+    std::copy(momentum.begin(), momentum.end(), fFourMomentum_.begin());
+    fCutPassing_.resize(cutPassing.size());
+    std::copy(cutPassing.begin(), cutPassing.end(), fCutPassing_.begin());
+    fHitPhi_.resize(phi.size());
+    std::copy(phi.begin(), phi.end(), fHitPhi_.begin());
+    fHitTheta_.resize(theta.size());
+    std::copy(theta.begin(), theta.end(), fHitTheta_.begin());
+    fPrimaryPhoton_.resize(primary.size());
+    std::copy(primary.begin(), primary.end(), fPrimaryPhoton_.begin());
+    fEdep_.resize(edep.size());
+    std::copy(edep.begin(), edep.end(), fEdep_.begin());
+    fEdepSmear_.resize(edepSmear.size());
+    std::copy(edepSmear.begin(), edepSmear.end(), fEdepSmear_.begin());
+}
+
+
 ///
 /// \brief Event::Event Copy constructor.
 /// \param est An instance of Event class to be copied.

src/event.h

@@ -29,6 +29,9 @@ class Event : public TObject
     public:
         Event();
         Event(std::vector<TLorentzVector*>* emissionCoordinates, std::vector<TLorentzVector*>* fourMomentum, double weight, DecayType type);
+        Event(std::vector<TLorentzVector> &sourcePos, std::vector<TLorentzVector> &pos, std::vector<TLorentzVector> &momentum,\
+        std::vector<double> &phi, std::vector<double> &theta, std::vector<bool> &cutPassing, std::vector<bool> &primary,\
+        std::vector<double> &edep, std::vector<double> &edepSmear, long Id, int decayType);
         Event(const Event& est);
         Event& operator=(const Event &est);
         virtual ~Event();
@@ -67,7 +70,7 @@ class Event : public TObject
         //number of event
         long fId;
         //ROOT stuff
-        ClassDef(Event, 17)
+        ClassDef(Event, 18)
 
     private:
         static long fCounter_; //static variable incremented with every call of a constructor (but not copy constructor)

tools/jpos2gate_converter/Makefile

@@ -0,0 +1,59 @@
+CXX=g++
+CXXFLAGS= -std=c++11 -Wall `root-config --cflags`
+LDFLAGS= `root-config --ldflags --glibs`
+
+OBJDIR=./obj
+SRCDIR=src
+H_FILES := $(wildcard $(SRCDIR)/*.h) 
+CPP_FILES := $(wildcard $(SRCDIR)/*.cpp) $(SRCDIR)/MyGateOutputDict.cpp $(SRCDIR)/event.cpp 
+OBJ_FILES := $(addprefix $(OBJDIR)/,$(notdir $(CPP_FILES:.cpp=.o)))
+
+EVPATH = "$(shell pwd)/$(SRCDIR)/"
+#checks if a dictionary exists
+DICT_EXISTS=$(shell [ -e "$(shell pwd)/$(OBJDIR)/MyGateOutputDict.o" ] && echo 1 || echo 0 )
+
+EVENT_FILES_PRESENT = 0
+COPY_FILES = $(SRCDIR)/EventDict.cpp $(SRCDIR)/event.h $(SRCDIR)/event.cpp $(SRCDIR)/event_linkdef.h $(SRCDIR)/constants.h 
+
+all:  converter 
+	@echo "COMPILATION COMPLETE!!!"
+
+converter: $(OBJDIR)/MyGateOutputDict.o $(OBJ_FILES) $(OBJDIR)/EventDict.o 
+	@echo "Creating executable: $@"
+	@echo ${CPP_FILES}
+	@(cp $(SRCDIR)/*.pcm . &&  $(CXX) -o converter $^ $(LDFLAGS))
+
+$(SRCDIR)/MyGateOutputDict.cpp: $(SRCDIR)/MyGateOutput.h $(COPY_FILES)
+	@echo "Compiling $@"
+	@(cd src && rootcint -f MyGateOutputDict.cpp -c $(CXXFLAGS) -p  MyGateOutput.h mygateoutput_linkdef.h)
+
+$(OBJDIR)/MyGateOutputDict.o: $(SRCDIR)/MyGateOutputDict.cpp
+	@echo "Compiling $@"
+	@$(CXX) $(SRCDIR)/MyGateOutputDict.cpp -o $(OBJDIR)/MyGateOutputDict.o -c $(CXXFLAGS)
+
+$(OBJDIR)/EventDict.o: ../../$(COPY_FILES)
+	@echo "Compiling EventDict.o"
+	@(cd src && rootcint -f EventDict.cpp -c $(CXXFLAGS) -p  event.h event_linkdef.h)
+	@$(CC) $(SRCDIR)/EventDict.cpp -o $(OBJDIR)/EventDict.o -c $(CXXFLAGS)
+
+$(OBJDIR)/%.o: $(SRCDIR)/%.cpp 
+	@echo "Compiling $@"
+	@$(CXX) $(CXXFLAGS) -c -o $@ $<
+
+$(OBJDIR)/%.o: $(SRCDIR)/%.C
+	@echo "Compiling $@"
+	@$(CXX) $(CXXFLAGS) -c -o $@ $<
+
+$(COPY_FILES):
+	@echo "Copying file: $@"
+	@cp ../../$@ ./src/
+
+clean:
+	@echo "Cleaning..."
+	@rm -f $(SRCDIR)/*.gch $(SRCDIR)/*.d $(SRCDIR)/MyGateOutputDict.cpp $(SRCDIR)/*.so $(SRCDIR)/Auto* $(OBJDIR)/*.o $(SRCDIR)/MyGateOutputDict* MyGateOutputDict* converter
+	@rm $(COPY_FILES)
+	@rm *.pcm $(SRCDIR)/*.pcm
+
+
+
+

tools/jpos2gate_converter/README.md

@@ -0,0 +1,80 @@
+# jpos2gate converter
+
+## Author: Rafał Masełek
+## Email: [email protected]
+
+### About:
+This a program to convert root outputs between JPOS and GATE. 
+!!! It is Beta version: no optimalization, little tests were made.
+
+### Requirements:
++ JPOS
++ pyROOT
++ Python 2.7 or newer
++ bash concole
+
+### Installation:
+To build the application type *make* in the console. Some file from the JPOS src directory will be copied and application will be built.
+
+### Running:
+Instead of directly running builded C++ application, use the provided convert.py Python script:
+>python convert.py [path to input file] [name of output file] [1 for JPOS->GATE conversion, 0 for GATE->JPOS]
+
+Example GATE->JPOS:
+>python convert ~/misie_pysie/gate_output.root "my_jpos_output.root" 0
+
+
+Example JPOS->GATE:
+>python convert ~/pysie_mysie/jpos_output.root "my_gate_output" 1
+
+Note that providing root extension is non-mandatory.
+
+### JPOS->GATE:
+Since it is Beta version, there are some constraints. The following fields in Hits tree (GATE output) will be filled with data:
++ time
++ sourcePosX
++ sourcePosY
++ sourcePosZ
++ posX
++ posY
++ posZ
++ edep
++ eventID
+
+The rest will be filled with some default values, that can be checked in main::jpos2gate function.
+
+If the folder name inside JPOS root is not "0_0_0_0_0_0", see section *"More on input files"*
+
+### GATE->JPOS:
+Since it is Beta version, there are some constraints. All JPOS files created by this converter will contain
+folder named "0_0_0_0_0_0" no matter what was the source in input (GATE) file.
+
+
+### More on input files
+jpos2gate converter uses automatically generated ROOT template classes to read data. It allows to use program for different data files without the need to recompile, as long as they have the same structure as the original ones. Unfortunately, JPOS output files contain subdirectories with names indicating the position and momentum of the source, so if one wants to convert JPOS file with subdirectory other than "0_0_0_0_0_0", one must regenerate template classes and recompile the converter. The instruction (algorithm) of what to do is provided below:
+
+#### You want to convert GATE->JPOS or JPOS->GATE with subdirectory named "0_0_0_0_0_0"? 
+
+Simply use the convert.py script as described in section *"Running"*.
+
+#### You want to convert JPOS->GATE with subdirectory named OTHER than "0_0_0_0_0_0"?
+
+You have to:
+1. Use prepare.py Python script to generate template classes. Type in bash console:
+>python prepare.py [path to JPOS root file] [path to GATE root file] [name of the subdirectory]
+
+**NOTE1:** You can use "data/gate.root" as a second parameter if you want.
+
+**NOTE2:** If your JPOS output file contains only one subdirectory, then providing only one parameter -- the path to JPOS
+root file should be sufficient.
+
+2. Compile the converter:
+
+>make
+
+3. Use convert.py Python script as described in section *"Running"*
+
+**IMPORTANT NOTE:** When you generate template classes for some subdirectory name, conversion 
+JPOS->GATE will be available only for JPOS root files containing subdirectories with the same name.
+To use it for different name, you have to redo the procedure (Do it if you want to be able to convert files
+with subdirectory named "0_0_0_0_0_0").

tools/jpos2gate_converter/convert.py

@@ -0,0 +1,30 @@
+#!/usr/bin/python 
+"""
+@author: Rafal Maselek
+@email: [email protected]
+This script copies input data files into data folder, change their name to meet the requirements
+of the converter application and runs the program.
+"""
+import sys
+import os
+
+if len(sys.argv) < 4:
+	print("Insufficient number of arguments!")
+	print("python convert.py [input file path] [output file name] [1 if goja2jpos, 0 otherwise]")
+	sys.exit(1)
+else:
+	# input file must have specific name
+	if sys.argv[3]:
+		in_name = "gate.root"
+	else:
+		in_name="jpos.root"
+	# check if user provided a name of the output file with or without extension
+	if sys.argv[2][-5:] == ".root":
+		out_name = sys.argv[2]
+	else:
+		out_name = sys.argv[2]+".root"
+	# prepare and execute bash commands
+	cp_command = "cp "+sys.argv[1]+" ./data/"+in_name
+	command = "./converter"+" "+out_name+" "+sys.argv[3]
+	os.system(cp_command)
+	os.system(command)

tools/jpos2gate_converter/obj/.gitignore

@@ -0,0 +1,2 @@
+*
+!.gitignore

tools/jpos2gate_converter/prepare.py

@@ -0,0 +1,87 @@
+#!/usr/bin/python
+#author: Rafal Maselek
+#email: [email protected]
+
+# this script generates EventHits files for a given .root file name, and folder name inside that file.
+
+import sys
+from ROOT import gROOT, TFile, TTree
+import os
+import shutil
+
+# check for command line arguments
+if(len(sys.argv) > 1):
+    nameJPOS = sys.argv[1]
+else:
+    nameJPOS = 'data/jpos.root'
+if(len(sys.argv) > 2):
+    nameGATE = sys.argv[2]
+else:
+    nameGATE = 'data/gate.root'
+if(len(sys.argv) == 4):
+    folder = sys.argv[3]
+else:
+    folder = None
+
+# generate files
+gROOT.Reset()
+gROOT.ProcessLine(".L ../../src/event.cpp")
+
+#########################################################################################################
+#########################################################################################################
+#########################################################################################################
+
+# FIRST DEAL WITH JPOS DATA
+f = TFile(nameJPOS)
+if folder:
+	tree = f.Get(folder+"/tree")
+else:
+	tree = f.FindObjectAny("tree")
+# tree.__class__ = TTree
+
+tree.MakeClass("JPOSOutput")
+
+thisFile = "JPOSOutput.C"
+base = os.path.splitext(thisFile)[0]
+newThisFile = base + ".cpp"
+os.rename(thisFile, newThisFile)
+# read the file
+file = open("./JPOSOutput.h", "r")
+contents = file.readlines()
+file.close()
+
+# add and write necessary content
+contents.insert(14, "#include <vector>\n")
+contents.insert(15, "#include <TLorentzVector.h>\n")
+for ii in range(0, len(contents)):
+    if "vector<" in contents[ii]:
+        contents[ii] = contents[ii].replace("vector<", "std::vector<")
+file = open("./JPOSOutput.h", "w")
+contents = "".join(contents)
+file.write(contents)
+file.close()
+
+path1 = os.path.abspath(newThisFile)
+path2 = file.name
+shutil.move(path1, os.path.join(os.path.dirname(path1), 'src/', newThisFile))
+shutil.move(path2, os.path.join(os.path.dirname(path2), 'src/',  os.path.basename(path2)))
+f.Close()
+#########################################################################################################
+#########################################################################################################
+#########################################################################################################
+
+# SECONDLY DEAL WITH GATE DATA
+f = TFile(nameGATE)
+tree = f.Get("Hits")
+tree.MakeClass("GateOutput")
+
+thisFile = "GateOutput.C"
+base = os.path.splitext(thisFile)[0]
+newThisFile = base + ".cpp"
+os.rename(thisFile, newThisFile)
+
+path1 = os.path.abspath(newThisFile)
+path2 = os.path.abspath("GateOutput.h")
+shutil.move(path1, os.path.join(os.path.dirname(path1), 'src/', newThisFile))
+shutil.move(path2, os.path.join(os.path.dirname(path2), 'src/',  "GateOutput.h"))
+f.Close()

tools/jpos2gate_converter/src/GateOutput.cpp

@@ -0,0 +1,43 @@
+#define GateOutput_cxx
+#include "GateOutput.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+
+void GateOutput::Loop()
+{
+//   In a ROOT session, you can do:
+//      root> .L GateOutput.C
+//      root> GateOutput t
+//      root> t.GetEntry(12); // Fill t data members with entry number 12
+//      root> t.Show();       // Show values of entry 12
+//      root> t.Show(16);     // Read and show values of entry 16
+//      root> t.Loop();       // Loop on all entries
+//
+
+//     This is the loop skeleton where:
+//    jentry is the global entry number in the chain
+//    ientry is the entry number in the current Tree
+//  Note that the argument to GetEntry must be:
+//    jentry for TChain::GetEntry
+//    ientry for TTree::GetEntry and TBranch::GetEntry
+//
+//       To read only selected branches, Insert statements like:
+// METHOD1:
+//    fChain->SetBranchStatus("*",0);  // disable all branches
+//    fChain->SetBranchStatus("branchname",1);  // activate branchname
+// METHOD2: replace line
+//    fChain->GetEntry(jentry);       //read all branches
+//by  b_branchname->GetEntry(ientry); //read only this branch
+   if (fChain == 0) return;
+
+   Long64_t nentries = fChain->GetEntriesFast();
+
+   Long64_t nbytes = 0, nb = 0;
+   for (Long64_t jentry=0; jentry<nentries;jentry++) {
+      Long64_t ientry = LoadTree(jentry);
+      if (ientry < 0) break;
+      nb = fChain->GetEntry(jentry);   nbytes += nb;
+      // if (Cut(ientry) < 0) continue;
+   }
+}