Merge pull request #9 from hyperk/merge_waveform_simulation

Implement waveform simulation and hit finding algorithm for mPMT

Author: kmtsui

.gitignore

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

README.md

@@ -53,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/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

@@ -48,6 +48,7 @@ int main(int argc, char **argv)
     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/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();
@@ -92,7 +106,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; }
+            if( !pr->ApplyDE(th,&(*hc)[tubeID]) ){ continue; }
 
             (&(*hc)[tubeID])->AddRawPE(th);
         }
@@ -186,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)
@@ -333,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()
@@ -345,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();
+        }
+    }
 }
 
 
@@ -415,6 +493,11 @@ void WCRootData::WriteTree()
     TFile *f = fWCSimT->GetCurrentFile();
     f->cd();
     fWCSimT->Write();
+    if (fSaveWF) 
+    {
+        fWCSimDigiWFT->Write();
+        fWCSimDigiPulls->Write();
+    }
     f->Close();
 }
 

cpp/include/HitDigitizer.h

@@ -6,6 +6,8 @@
 #include "PMTResponse.h"
 #include "MTRandom.h"
 
+#include "TH1F.h"
+
 using std::cout;
 using std::endl;
 using std::vector;
@@ -15,13 +17,13 @@ class HitDigitizer
     public:
         HitDigitizer(int s=67823);
         virtual ~HitDigitizer();
-        void Digitize(HitTubeCollection*, PMTResponse*);
-        void DigitizeTube(HitTube*, PMTResponse*);
+        virtual void Digitize(HitTubeCollection*, PMTResponse*);
+        virtual void DigitizeTube(HitTube*, PMTResponse*);
 
-        void ApplyThreshold(double&, bool&);
-        double DoTruncate(const double, const double);
+        virtual void ApplyThreshold(double&, bool&);
+        virtual double DoTruncate(const double, const double);
 
-    private:
+    protected:
         float fPrecisionCharge;
         float fPrecisionTiming;
         float fEfficiency;
@@ -30,3 +32,31 @@ class HitDigitizer
 
         MTRandom *fRand;
 };
+
+// mPMT specific digitizer
+class HitDigitizer_mPMT : public HitDigitizer
+{
+    public:
+        HitDigitizer_mPMT(int s=67823);
+        virtual ~HitDigitizer_mPMT();
+        void LoadWaveform(const string &filename);
+        void DigitizeTube(HitTube*, PMTResponse*);
+        void ApplyThreshold(double&, bool&);
+        TH1F BuildWavetrain(const vector<TrueHit*> PEs, double waveform_window);
+        void FitWavetrain(TH1F hist, vector<double>& vDigiT, vector<double>& vDigiQ);
+
+    private:
+        TH1F* hWF;
+        float fAmplitudeSigma;
+        float fDt;
+        float fDv;
+        float fWaveformOffset;
+        float fADCToPE;
+        int fIntegralPreceding;
+        int fIntegralFollowing;
+        int fChargeWindowBefore;
+        int fChargeWindowAfter;
+        float fDigiTimeOffset;
+        int fHitInsensitivityPeriod;
+        float fAmplitudeThreshold;
+};

cpp/include/HitTube.h

@@ -22,6 +22,8 @@
 #include "TrueHit.h"
 using std::vector;
 
+#include "TH1F.h"
+
 class HitTube
 {
     public: 
@@ -74,6 +76,22 @@ class HitTube
         float GetChargeDigi(const int i) const { return fChargeDigi[i]; }
         const vector<int>& GetParentCompositionDigi(const int i) const { return fParentCompDigi[i]; }
 
+        void SetDigiWF(const TH1F& hist)
+        {
+            if(hDigiWF != nullptr) delete hDigiWF;
+            hDigiWF = new TH1F(hist);
+            hDigiWF->SetDirectory(0);
+        }
+        TH1F* GetDigiWF()
+        {
+            return hDigiWF;
+        }
+        void SetDigiPulls(float pullT, float pullQ) { fPullT = pullT; fPullQ = pullQ;}
+        void SetTrueTQ(float trueT, float trueQ) { fTrueT = trueT; fTrueQ = trueQ;}
+        float GetPullT() {return fPullT;}
+        float GetPullQ() {return fPullQ;}
+        float GetTrueT() {return fTrueT;}
+        float GetTrueQ() {return fTrueQ;}
 
     private:
         int fNRawPE;
@@ -88,4 +106,10 @@ class HitTube
         vector<float> fTimeDigi;
         vector<float> fChargeDigi;
         vector<vector<int>> fParentCompDigi;
+
+        TH1F* hDigiWF;
+        float fPullT; 
+        float fPullQ;
+        float fTrueT; 
+        float fTrueQ;
 };

cpp/include/MDTManager.h

@@ -37,8 +37,7 @@ class MDTManager
         PMTResponse* GetPMTResponse(const string &s="Def") { return fPMTResp[s]; }
 
     private:
-        // TriggerAlgo *fTrigAlgo;
-        HitDigitizer *fDgtzr;
+        map<string, HitDigitizer*> fDgtzr;
         MTRandom *fRndm;
 
         map<string, TriggerAlgo*> fTrigAlgo;

cpp/include/PMTResponse.h

@@ -18,6 +18,7 @@ class PMTResponse
         virtual ~PMTResponse() {};
         virtual double GetRawSPE(const TrueHit* th=NULL, const HitTube* ht=NULL) = 0; 
         virtual float HitTimeSmearing(float) = 0;
+        virtual float HitTimeSmearing(float,int) = 0;
         virtual void Initialize(int, const string &s="") = 0;
         virtual bool ApplyDE(const TrueHit* th=NULL, const HitTube* ht=NULL) = 0;
 
@@ -38,12 +39,21 @@ class GenericPMTResponse : public PMTResponse
         void Initialize(int, const string &s="");
         virtual double GetRawSPE(const TrueHit* th=NULL, const HitTube* ht=NULL);
         virtual float HitTimeSmearing(float);
+        virtual float HitTimeSmearing(float,int);
         virtual bool ApplyDE(const TrueHit* th=NULL, const HitTube* ht=NULL);
 
     protected:
         void LoadCDFOfSPE(const string &s);
         float fqpe0[501];
         string fTxtFileSPECDF;
+        void LoadPMTDE(const string &s);
+        int fLoadDE;
+        std::vector<double> fDE;
+        string fPMTDEFile;
+        void LoadPMTTime(const string &s);
+        int fLoadT;
+        std::vector<double> fT;
+        string fPMTTFile;
 
     //private:
     protected: