Implement waveform simulation and hit finding algorithm for mPMT

.gitignore

@@ -7,8 +7,10 @@ appWCTESingleEvent
 appIWCDSingleEvent
 appGenPileUpSpill
 
-# root file
+# output files
 *.root
+*.pdf
+
 
 .idea
 build-*

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
@@ -57,3 +53,25 @@ Input variables are:
 - `IDNuIntRate`,`BeamBkgRate`: Mean number of ID and background events per spill. In each spill, the actual number of interactions are drawn from a Poisson distribution, and interaction timing according to the bunch structure (see `BeamTiming` class under `app/utilities/WCRootData/`)
 - `UseOD`: Process OD hits
 - `NumOfSpillsSavedPerFile`, `TotalNumOfSpills`: output spill setup
+
+## How to simulate and access digitized pulses
+Turn on waveform simulation in the parameter file by `< DigitizerType_PMTType = 1 >`.
+
+For each true hit, a digitized waveform is simulated by sampling the single PE pulse (defined by the `< WaveformFile >` parameter) every 8 ns with 0.1 resolution. If there is another PE arriving within the hit integration window, the waveforms are added. 
+
+To do pulse fitting, the hit finding algorithm [here](https://github.com/hyperk/MDT/issues/8) is implemented to calculate the digitized time and charge.
+
+Optionally, the waveform and digi TQ pulls of the first pulse of each PMT in each event can be saved by setting `< SaveWaveform = 1 >`. To read the pulses,
+```
+// open the file and get the digitzed waveform tree
+TTree* wcsimDigiWFTree = (TTree*)f->Get("wcsimDigiWFTree");
+TClonesArray *arr = new TClonesArray("TH1F");
+wcsimDigiWFTree->GetBranch("wcsimrootevent_waveform")->SetAutoDelete(kFALSE);
+wcsimDigiWFTree->SetBranchAddress("wcsimrootevent_waveform",&arr);
+// In each event, each array index corresponds to PMTId-1 (from 0 to nPMTs-1)
+wcsimDigiWFTree->GetEntry(0); // EvtId
+TH1F* h = (TH1F*)arr->At(0); // PMTId-1
+// Time and charge pulls of digitized hits
+TTree* WCSimDigiPulls = (TTree*)f->Get("WCSimDigiPulls");
+WCSimDigiPulls->Scan("PullT:TrueT:PullQ:TrueQ:EvtId:PMTId","PullT>-99"); // Pull==-99 means no digi hit for this PMT in this event
+```

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/PMTResponse3inchR12199_02.cc

@@ -30,6 +30,8 @@ void PMTResponse3inchR12199_02::Initialize(int seed, const string &pmtname)
     s["TimingResMinimum"] = "TimingResMinimum"; 
     s["ScalFactorTTS"] = "ScalFactorTTS";
     s["SPECDFFile"] = "SPECDFFile";
+    s["PMTDE"] = "PMTDE";
+    s["PMTTime"] = "PMTTime";
     if( fPMTType!="" )
     {
         map<string, string>::iterator i;
@@ -43,7 +45,11 @@ void PMTResponse3inchR12199_02::Initialize(int seed, const string &pmtname)
     Conf->GetValue<float>(s["TimingResMinimum"], fTResMinimum);
     Conf->GetValue<float>(s["ScalFactorTTS"], fSclFacTTS);
     Conf->GetValue<string>(s["SPECDFFile"], fTxtFileSPECDF);
+    Conf->GetValue<string>(s["PMTDE"], fPMTDEFile);
+    Conf->GetValue<string>(s["PMTTime"], fPMTTFile);
     this->LoadCDFOfSPE(fTxtFileSPECDF);
+    this->LoadPMTDE(fPMTDEFile);
+    this->LoadPMTTime(fPMTTFile);
 }
 
 float PMTResponse3inchR12199_02::HitTimeSmearing(float Q)
@@ -53,3 +59,6 @@ float PMTResponse3inchR12199_02::HitTimeSmearing(float Q)
     if( timingResolution<fTResMinimum ){ timingResolution = fTResMinimum; }
     return fRand->Gaus(0.0,timingResolution);
 }
+
+///////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////

app/application/appWCTESingleEvent.cc

@@ -42,12 +42,13 @@ 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"};
+    const vector<string> listWCRootCopyTree{"wcsimGeoT","Settings","wcsimRootOptionsT"};
 
 	// WCRootData is an interface class between MDT and WCSim root file
     WCRootData *inData = new WCRootData();
@@ -85,6 +86,7 @@ int main(int argc, char **argv)
         MDT->DoInitialize();
     }
     outData->WriteTree();
+    for (auto s : listWCRootCopyTree) outData->CopyTree(fInFileName.c_str(),s.c_str());
     inData->CloseFile();
 }
 

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/include/WCRootData.h

@@ -69,6 +69,21 @@ class WCRootData
         float fHitTimeOffset;
         bool fMultDigiHits;
 
+        TTree *fWCSimDigiWFT;
+        std::vector<TClonesArray*> fDigiWF;
+        // TClonesArray* fDigiWF;
+        // TClonesArray* fDigiWF2;
+        // TClonesArray* fDigiWFOD;
+        bool fSaveWF;
+
+        TTree *fWCSimDigiPulls;
+        float fPullQ;
+        float fPullT;
+        float fTrueQ;
+        float fTrueT;
+        int fEvtId;
+        int fPMTId;
+
     private:
 		void SetTubes(HitTubeCollection*, const int);
         TString fOutFileName;

app/utilities/WCRootData/src/WCRootData.cc

@@ -15,19 +15,33 @@ WCRootData::WCRootData()
 
     int mult_flag = 1;
     fHitTimeOffset = 0.;
+    int save_wf = 0;
 
     Configuration *Conf = Configuration::GetInstance();
     Conf->GetValue<float>("TimeOffset", fHitTimeOffset);
     Conf->GetValue<int>("FlagMultDigits", mult_flag);
+    Conf->GetValue<int>("SaveWaveform", save_wf);
 
     fMultDigiHits = bool(mult_flag);
+    fSaveWF = bool(save_wf);
+
+    fWCSimDigiWFT = 0;
+    fWCSimDigiPulls = 0;
+    fPullQ = -99.;
+    fPullT = -99.;
+    fTrueQ = -99.;
+    fTrueT = -99.;
+    fEvtId = -99;
+    fPMTId = -99;
 }
 
 WCRootData::~WCRootData()
 {
     if( fWCGeom ){ delete fWCGeom; fWCGeom = 0; }
     if( fWCSimT ){ delete fWCSimT; fWCSimT = 0; }
     if( fWCSimC ){ delete fWCSimC; fWCSimC = 0; }
+    if( fWCSimDigiWFT ){ delete fWCSimDigiWFT; fWCSimDigiWFT = 0; }
+    if( fWCSimDigiPulls ){ delete fWCSimDigiPulls; fWCSimDigiPulls = 0; }
     fSpEvt.clear();
     fSpEvt.shrink_to_fit();
     isOD.clear();
@@ -91,7 +105,8 @@ 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)); }
-            if( !pr->ApplyDE(th) ){ continue; }
+            th->SetCreatorProcess((int)(aHitTime->GetPhotonCreatorProcess()));
+            if( !pr->ApplyDE(th,&(*hc)[tubeID]) ){ continue; }
 
             (&(*hc)[tubeID])->AddRawPE(th);
         }
@@ -185,7 +200,28 @@ void WCRootData::CreateTree(const char *filename, const vector<string> &list)
             fWCSimT->Branch(list[i].c_str(), bAddress, &fSpEvt[i], bufferSize, 2);
             if ( list[i].find("OD")!=std::string::npos ) isOD[i] = true;
         }
+
+        if (fSaveWF) 
+        {
+            fWCSimDigiWFT = new TTree("wcsimDigiWFTree","Digitized waveform for each PMT");
+            fDigiWF.clear();
+            fDigiWF = vector<TClonesArray*>(list.size(), 0);
+            for(unsigned int i=0; i<list.size(); i++)
+            {
+                fDigiWF[i] = new TClonesArray("TH1F");
+                fWCSimDigiWFT->Branch(Form("%s_waveform",list[i].c_str()),&fDigiWF[i]);
+            }
+
+            fWCSimDigiPulls = new TTree("WCSimDigiPulls","Time and charge pulls of digitized hits");
+            fWCSimDigiPulls->Branch("PullQ",&fPullQ);
+            fWCSimDigiPulls->Branch("PullT",&fPullT);
+            fWCSimDigiPulls->Branch("TrueQ",&fTrueQ);
+            fWCSimDigiPulls->Branch("TrueT",&fTrueT);
+            fWCSimDigiPulls->Branch("EvtId",&fEvtId);
+            fWCSimDigiPulls->Branch("PMTId",&fPMTId);
+        }
     }
+
 }
 
 void WCRootData::AddDigiHits(MDTManager *mdt, int eventID, int iPMT)
@@ -207,6 +243,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 +264,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 +276,8 @@ void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID
                                 photonStartPos,
                                 photonEndPos,
                                 photonStartDir,
-                                photonEndDir);
+                                photonEndDir,
+                                photonCreatorProcess);
 
         truetime.clear();
         primaryParentID.clear();
@@ -247,6 +286,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++) 
@@ -328,6 +368,41 @@ void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID
         WCSimRootEventHeader *eh = anEvent->GetHeader();
         eh->SetDate( int(triggerTime) );
     }
+
+    if (fSaveWF)
+    {
+        TClonesArray &fDigiWFarray = *(fDigiWF[iPMT]);
+        for(hc->Begin(); !hc->IsEnd(); hc->Next())
+        {
+            HitTube *aPH = &(*hc)();
+
+            // if (aPH->GetDigiWF())
+            //     new(fDigiWFarray[aPH->GetTubeID()-1]) TH1F(*(aPH->GetDigiWF()));
+            // else new(fDigiWFarray[aPH->GetTubeID()-1]) TH1F();
+            TH1F* h = (TH1F*)fDigiWFarray.ConstructedAt(aPH->GetTubeID()-1);
+            if (aPH->GetDigiWF()) 
+            {
+                int nbins = aPH->GetDigiWF()->GetNbinsX();
+                h->SetBins(nbins,aPH->GetDigiWF()->GetBinLowEdge(1),aPH->GetDigiWF()->GetBinLowEdge(nbins+1));
+                for (int i=1;i<=nbins;i++) h->SetBinContent(i,aPH->GetDigiWF()->GetBinContent(i));
+            }
+            else
+            {
+                h->Reset();
+            }
+
+
+            fPullQ = aPH->GetPullQ();
+            fPullT = aPH->GetPullT();
+            fTrueQ = aPH->GetTrueQ();
+            fTrueT = aPH->GetTrueT();
+            fEvtId = eventID;
+            fPMTId = aPH->GetTubeID();
+            fWCSimDigiPulls->Fill();
+
+            aPH = NULL;
+        } // PMT loop
+    }
 }
 
 void WCRootData::FillTree()
@@ -340,6 +415,14 @@ void WCRootData::FillTree()
     {
         fSpEvt[i]->ReInitialize();
     }
+    if (fSaveWF) 
+    {
+        fWCSimDigiWFT->Fill();
+        for(unsigned int i=0; i<fDigiWF.size(); i++)
+        {
+            fDigiWF[i]->Clear();
+        }
+    }
 }
 
 
@@ -358,7 +441,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 +476,7 @@ void WCRootData::AddTracks(const WCSimRootTrigger *aEvtIn, float offset_time, in
 			  stop,
 			  start,
 			  parenttype,
+              creatorProcess,
 			  time,
 			  id,
 			  idPrnt,
@@ -409,6 +493,11 @@ void WCRootData::WriteTree()
     TFile *f = fWCSimT->GetCurrentFile();
     f->cd();
     fWCSimT->Write();
+    if (fSaveWF) 
+    {
+        fWCSimDigiWFT->Write();
+        fWCSimDigiPulls->Write();
+    }
     f->Close();
 }
 

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