Merge pull request #7 from hyperk/update_to_WCSim_v1.12.16

Update to compile with WCSim v1.12.16 and new WCTE PMT response

Author: kmtsui

README.md

@@ -18,11 +18,7 @@ The library provides C++ classes that manage three tasks:
 Simple runtime procedures. First set up your ROOT and WCSIM:
 ```
 source your_thisroot.sh
-export WCSIMDIR=your_WCSIM_installation
-```
-If `libWCSimRoot.so` is not directly under `WCSIMDIR`, export to `WCSIMROOTDIR`.
-```
-export WCSIMROOTDIR=your_libWCSimRoot.so_installation
+source your_this_wcsim.sh # or just export WCSIM_BUILD_DIR
 ```
 Then set up the MDT environment.
 ```
@@ -32,7 +28,7 @@ cd $MDTROOT/app/utilities/WCRootData; make clean; make all
 cd $MDTROOT/app/application; make clean; make all
 cd $MDTROOT
 # edit variables properly in run_test_mdt4wcte.sh
-bash run_test_mdt4iwcd.sh
+bash run_test_mdt4wcte.sh
 ```
 
 ## IWCD usage

app/application/Makefile

@@ -2,12 +2,12 @@ CXX=g++
 LD=g++
 
 CXXFLAGS += -Wall -std=c++11 -g $(shell root-config --cflags)\
-			-I$(WCSIMDIR)/include\
+			-I$(WCSIM_BUILD_DIR)/include/WCSim\
 			-I${MDTROOT}/cpp/include\
 			-I${WCRDROOT}/include
 
 LDFLAGS += $(shell root-config --ldflags) $(shell root-config --libs) -lTreePlayer\
-		   -L$(WCSIMROOTDIR) -lWCSimRoot\
+		   -L$(WCSIM_BUILD_DIR)/lib -lWCSimRoot\
 		   -L${MDTROOT}/cpp -lMDT\
 		   -L${WCRDROOT} -lWCRData
 

app/application/appWCTESingleEvent.cc

@@ -42,10 +42,10 @@ int main(int argc, char **argv)
 	// WCTE will use single PMT type, so define the corresponding type of 3-inch PMT
     const int NPMTType = 1;
     string fPMTType[NPMTType];
-    fPMTType[0] = "PMT3inchR12199_02";
+    fPMTType[0] = "PMT3inchR14374_WCTE";
 
     MDTManager *MDT = new MDTManager(fSeed);
-    MDT->RegisterPMTType(fPMTType[0], new PMTResponse3inchR12199_02());
+    MDT->RegisterPMTType(fPMTType[0], new Response3inchR14374_WCTE());
 
     const vector<string> listWCRootEvt{"wcsimrootevent"};
 

app/utilities/WCRootData/Makefile

@@ -1,9 +1,9 @@
 CXX=g++
 LD=g++
 
-CXXFLAGS += -Wall -O2 -std=c++11 -g -fPIC -I./include $(shell root-config --cflags) -I$(WCSIMDIR)/include -I${MDTROOT}/cpp/include
+CXXFLAGS += -Wall -O2 -std=c++11 -g -fPIC -I./include $(shell root-config --cflags) -I$(WCSIM_BUILD_DIR)/include/WCSim -I${MDTROOT}/cpp/include
 #LDFLAGS += -shared $(shell root-config --ldflags) $(shell root-config --libs) -lTreePlayer -lMinuit2 -L$(WCSIMDIR) -lWCSimRoot -L${MDTROOT}/cpp -lMDT
-LDFLAGS += -shared $(shell root-config --ldflags) $(shell root-config --libs) -L$(WCSIMROOTDIR) -lWCSimRoot -L${MDTROOT}/cpp -lMDT
+LDFLAGS += -shared $(shell root-config --ldflags) $(shell root-config --libs) -L$(WCSIM_BUILD_DIR)/lib -lWCSimRoot -L${MDTROOT}/cpp -lMDT
 
 TARGET=libWCRData.so
 

app/utilities/WCRootData/src/WCRootData.cc

@@ -91,6 +91,7 @@ void WCRootData::AddTrueHitsToMDT(HitTubeCollection *hc, PMTResponse *pr, float
 
             th->SetStartTime(aHitTime->GetPhotonStartTime()+intTime);
             for(int k=0; k<3; k++){ th->SetStartPosition(k, aHitTime->GetPhotonStartPos(k)); }
+            th->SetCreatorProcess((int)(aHitTime->GetPhotonCreatorProcess()));
             if( !pr->ApplyDE(th) ){ continue; }
 
             (&(*hc)[tubeID])->AddRawPE(th);
@@ -207,6 +208,7 @@ void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID
     std::vector<TVector3> photonEndPos;
     std::vector<TVector3> photonStartDir;
     std::vector<TVector3> photonEndDir;
+    std::vector<ProcessType_t> photonCreatorProcess;
     for(hc->Begin(); !hc->IsEnd(); hc->Next())
     {
         // Get tube ID
@@ -227,6 +229,7 @@ void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID
             photonEndPos.push_back(TVector3(PEs[iPE]->GetPosition(0),PEs[iPE]->GetPosition(1),PEs[iPE]->GetPosition(2)));
             photonStartDir.push_back(TVector3(PEs[iPE]->GetStartDirection(0),PEs[iPE]->GetStartDirection(1),PEs[iPE]->GetStartDirection(2)));
             photonEndDir.push_back(TVector3(PEs[iPE]->GetDirection(0),PEs[iPE]->GetDirection(1),PEs[iPE]->GetDirection(2)));
+            photonCreatorProcess.push_back((ProcessType_t)(PEs[iPE]->GetCreatorProcess()));
         }
 
         anEvent->AddCherenkovHit(tubeID,
@@ -238,7 +241,8 @@ void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID
                                 photonStartPos,
                                 photonEndPos,
                                 photonStartDir,
-                                photonEndDir);
+                                photonEndDir,
+                                photonCreatorProcess);
 
         truetime.clear();
         primaryParentID.clear();
@@ -247,6 +251,7 @@ void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID
         photonEndPos.clear();
         photonStartDir.clear();
         photonEndDir.clear();
+        photonCreatorProcess.clear();
     }
     const int nTriggers = ti->GetNumOfTrigger(); 
     for(int iTrig=0; iTrig<nTriggers; iTrig++) 
@@ -358,7 +363,7 @@ void WCRootData::AddTracks(const WCSimRootTrigger *aEvtIn, float offset_time, in
         Int_t     parenttype = aTrack->GetParenttype();
         Int_t     id = aTrack->GetId();
         Int_t     idPrnt = aTrack->GetParentId();
-
+        ProcessType_t creatorProcess = aTrack->GetCreatorProcess();
 
         Double_t   dir[3];
         Double_t   pdir[3];
@@ -393,6 +398,7 @@ void WCRootData::AddTracks(const WCSimRootTrigger *aEvtIn, float offset_time, in
 			  stop,
 			  start,
 			  parenttype,
+              creatorProcess,
 			  time,
 			  id,
 			  idPrnt,

cpp/Makefile

@@ -1,8 +1,8 @@
 CXX=g++
 LD=g++
 
-CXXFLAGS += -Wall -O2 -std=c++11 -g -fPIC -I./include
-LDFLAGS += -shared
+CXXFLAGS += -Wall -O2 -std=c++11 -g -fPIC -I./include $(shell root-config --cflags)
+LDFLAGS += -shared $(shell root-config --ldflags) $(shell root-config --libs)
 
 TARGET=libMDT.so
 SRCDIR = ./src

cpp/include/PMTResponse.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <string>
+#include <TGraph.h>
 #include "MTRandom.h"
 #include "HitTube.h"
 
@@ -78,3 +79,16 @@ class Response3inchR14374 : public GenericPMTResponse
     private:
         double fTimeResAt1PE;
 };
+
+class Response3inchR14374_WCTE : public GenericPMTResponse
+{
+    public: 
+        Response3inchR14374_WCTE(int, const string &s="");
+        Response3inchR14374_WCTE();
+        virtual ~Response3inchR14374_WCTE();
+        float HitTimeSmearing(float);
+        void Initialize(int, const string &s="");
+        
+    private:
+        TGraph  *gTResol = nullptr; // timing resolution
+};

cpp/include/TrueHit.h

@@ -15,13 +15,15 @@ class TrueHit
         float GetStartPosition(int i) const { return fStartPosition[i]; }
         float GetStartDirection(int i) const { return fStartDirection[i]; }
         int GetPosBin(int i) const { return fBin[i]; }
+        int GetCreatorProcess() const { return fCreatorProcess; }
 
         void SetPosition(int i, float f) { fPosition[i] = f; }
         void SetDirection(int i, float f) { fDirection[i] = f; }
         void SetStartTime(float f) { fStartTime = f; }
         void SetStartPosition(int i, float f) { fStartPosition[i] = f; }
         void SetStartDirection(int i, float f) { fStartDirection[i] = f; }
         void SetPosBin(int i, int b){ fBin[i] = b; }
+        void SetCreatorProcess(int i){ fCreatorProcess = i; }
 
     private:
         double fTime;
@@ -32,4 +34,5 @@ class TrueHit
         float fStartDirection[3]; // Photon track initial direction
         int fParentId;
         int fBin[3];
+        int fCreatorProcess;
 };

cpp/src/PMTResponse.cc

@@ -229,3 +229,78 @@ float Response3inchR14374::HitTimeSmearing(float Q)
     float timingResolution = 0.6*fSclFacTTS;
     return fRand->Gaus(0., timingResolution);
 }
+
+///////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////
+Response3inchR14374_WCTE::Response3inchR14374_WCTE(int seed, const string &pmtname)
+{
+    double charge[14] = 
+    {
+        0.2, 0.4, 0.6, 0.8, 1.0,
+        1.2, 1.4, 1.6, 1.8, 2.0,
+        2.5, 3.0, 3.5, 4.0
+    };
+    double resol[14] =
+    {
+        1.1654, 0.61088, 0.4186, 0.32532, 0.26484,
+        0.23084, 0.20969, 0.19297, 0.17716, 0.17046,
+        0.15455, 0.1427, 0.13699, 0.13229
+    };
+    gTResol = new TGraph(14,charge,resol);
+
+    this->Initialize(seed, pmtname);
+}
+
+Response3inchR14374_WCTE::Response3inchR14374_WCTE()
+{
+    double charge[14] = 
+    {
+        0.2, 0.4, 0.6, 0.8, 1.0,
+        1.2, 1.4, 1.6, 1.8, 2.0,
+        2.5, 3.0, 3.5, 4.0
+    };
+    double resol[14] =
+    {
+        1.1654, 0.61088, 0.4186, 0.32532, 0.26484,
+        0.23084, 0.20969, 0.19297, 0.17716, 0.17046,
+        0.15455, 0.1427, 0.13699, 0.13229
+    };
+    gTResol = new TGraph(14,charge,resol);
+}
+
+Response3inchR14374_WCTE::~Response3inchR14374_WCTE()
+{
+    delete gTResol;
+}
+
+void Response3inchR14374_WCTE::Initialize(int seed, const string &pmtname)
+{
+    fPMTType = pmtname;
+    fRand = new MTRandom(seed);
+
+    map<string, string> s;
+    s["ScalFactorTTS"] = "ScalFactorTTS";
+    s["SPECDFFile"] = "SPECDFFile";
+    if( fPMTType!="" )
+    {
+        map<string, string>::iterator i;
+        for(i=s.begin(); i!=s.end(); i++)
+        {
+            i->second += "_" + fPMTType;
+        }
+    }
+    Configuration *Conf = Configuration::GetInstance();
+    Conf->GetValue<float>(s["ScalFactorTTS"], fSclFacTTS);
+    Conf->GetValue<string>(s["SPECDFFile"], fTxtFileSPECDF);
+    this->LoadCDFOfSPE(fTxtFileSPECDF);
+}
+
+float Response3inchR14374_WCTE::HitTimeSmearing(float Q)
+{
+    float pmt_tts = 1.5;
+    if (Q>4.0) Q = 4.0; // limit Q to valid range
+    float val = gTResol->Eval(Q,0,"S");
+    float timingResolution = sqrt(pmt_tts*pmt_tts+val*val)/2.355; // conversion from FWHM to sigma
+    timingResolution *= fSclFacTTS;
+    return fRand->Gaus(0.0,timingResolution);
+}

envMDT.sh

@@ -27,7 +27,7 @@ if [ -n "${MDTROOT}" ] ; then
 		drop_from_path "$LD_LIBRARY_PATH" "${old_mdtroot}/cpp"
 		LD_LIBRARY_PATH=$newpath
 
-		drop_from_path "$LD_LIBRARY_PATH" "${WCSIMROOTDIR}"
+		drop_from_path "$LD_LIBRARY_PATH" "${WCSIM_BUILD_DIR}/lib"
 		LD_LIBRARY_PATH=$newpath
 
 		drop_from_path "$LD_LIBRARY_PATH" "${WCRDROOT}"
@@ -59,18 +59,15 @@ else
 fi
 
 #--- WCSim
-if [[ -z "${WCSIMDIR}" ]]; then
-	echo "WCSIMDIR not set. Exit!"
+if [[ -z "${WCSIM_BUILD_DIR}" ]]; then
+	echo "WCSIM_BUILD_DIR not set. Exit!"
 	return 
 else
-	echo "Using WCSIM installed in $WCSIMDIR"
+	echo "Using WCSIM installed in $WCSIM_BUILD_DIR"
 	#export LD_LIBRARY_PATH=${WCSIMDIR}:$LD_LIBRARY_PATH
 
-	if [[ -z "${WCSIMROOTDIR}" ]]; then
-		WCSIMROOTDIR=$WCSIMDIR
-	fi
-	echo "Using libWCSimRoot.so installed in $WCSIMROOTDIR"
-	export LD_LIBRARY_PATH=${WCSIMROOTDIR}:$LD_LIBRARY_PATH
+	echo "Using libWCSimRoot.so installed in $WCSIM_BUILD_DIR/lib"
+	export LD_LIBRARY_PATH=${WCSIM_BUILD_DIR}/lib:$LD_LIBRARY_PATH
 fi
 
 #--- MDT's utility