Merge pull request #5 from hyperk/feature/PileUpApp

Feature/pile up app

.gitignore

@@ -5,6 +5,7 @@
 # binary
 appWCTESingleEvent
 appIWCDSingleEvent
+appGenPileUpSpill
 
 # root file
 *.root

README.md

@@ -29,8 +29,31 @@ Then set up the MDT environment.
 source envMDT.sh
 cd $MDTROOT/cpp; make clean; make all
 cd $MDTROOT/app/utilities/WCRootData; make clean; make all
-cd $MDTROOT/app/application; make clean; make appIWCDSingleEvent
+cd $MDTROOT/app/application; make clean; make all
 cd $MDTROOT
 # edit variables properly in run_test_mdt4wcte.sh
 bash run_test_mdt4iwcd.sh
 ```
+
+## IWCD usage
+Essentially the same usage but with OD support.
+```
+# edit variables properly in run_test_mdt4iwcd.sh
+bash run_test_mdt4iwcd.sh
+```
+
+## IWCD pile-up generation
+Example usage
+```
+$MDTROOT/app/application/appGenPileUpSpill $MDTROOT/example/genPileUpConfig.txt
+```
+The application generates pile-up events by combining ID neutrino interaction events and beam background events in a spill. 
+
+Input variables are:
+- `ListIDNuIntFiles`,`ListBeamBkgFiles`: list of WCSim output files for the ID and background events
+- `OutFileNamePrefix`: output file name will be something like "OutFileNamePrefix".00000.root
+- `MDTParFile`: MDT config file
+- `InitialSeed`: Random seed
+- `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

app/application/CMakeLists.txt

@@ -28,6 +28,7 @@ pbuilder_install_headers(${HEADERS})
 set(exe_sources
         appIWCDSingleEvent.cc
         appWCTESingleEvent.cc
+        appGenPileUpSpill.cc
 )
 
 pbuilder_executables(

app/application/Makefile

@@ -23,10 +23,18 @@ appIWCDSingleEvent: appIWCDSingleEvent.o PMTResponse3inchR12199_02.o
 	$(LD) $^ $(LDFLAGS) -o $@
 	@echo "$@ done"
 
+appGenPileUpSpill: appGenPileUpSpill.o PMTResponse3inchR12199_02.o
+	$(RM) $@
+	$(LD) $^ $(LDFLAGS) -o $@
+	@echo "$@ done"
+
 %o::%cc
 	$(CXX) $(CXXFLAGS) -c $? -o $@
 
+all: appWCTESingleEvent appIWCDSingleEvent appGenPileUpSpill
+
 clean:
 	$(RM) *.o
 	$(RM) -f appWCTESingleEvent
 	$(RM) -f appIWCDSingleEvent
+	$(RM) -f appGenPileUpSpill

app/application/appGenPileUpSpill.cc

@@ -10,6 +10,9 @@
 #include "WCRootDataPileUpSpill.h"
 #include "BeamTiming.h"
 
+// PMT type used for 3-inch PMTs of IWCD/WCTE
+#include "PMTResponse3inchR12199_02.h"
+
 std::string fInFileTextIDNuInt = "";
 std::string fInFileTextBeamBkg = "";
 std::string fOutFileNamePrefix = "";
@@ -19,11 +22,19 @@ int fNumOfSpillsPerFile = 1000;
 float fIDNuIntRate = 3.5;
 float fBeamBkgRate = 7.5;
 int fTotalNumofSpills = 10000;
+int fUseOD = 0;
 
 void ReadConfiguration(const char*);
 
 int main(int argc, char **argv)
 {
+    if (argc==1)
+    {
+        std::cout<<" ** Please specify input config file ** "<<std::endl;
+        std::cout<<" ** Usage: appGenPileUpSpill config_file ** "<<std::endl;
+        return 0;
+    }
+
     std::string ConfigFileName = argv[1];
     ReadConfiguration(argv[1]);
 
@@ -39,8 +50,20 @@ int main(int argc, char **argv)
     const int seed_beamtiming = rndm->Integer(1000000);
     delete rndm; rndm = 0;
 
+    // IWCD PMT type for ID and OD
+    int NPMTType = !fUseOD ? 1 : 2 ;
+    vector<string> fPMTType(NPMTType);
+    fPMTType[0] = "PMT3inchR12199_02";
+    if (fUseOD>0) fPMTType[1] = "PMT3inchR12199_02_OD";
+
     // Will be manaing marging true hits and digitizing merged true hits
     MDTManager *MDT = new MDTManager(seed_mdt);
+    MDT->RegisterPMTType(fPMTType[0], new PMTResponse3inchR12199_02());
+    if (fUseOD) MDT->RegisterPMTType(fPMTType[1], new PMTResponse3inchR12199_02());
+
+    vector<string> listWCRootEvt(NPMTType);
+    listWCRootEvt[0] = "wcsimrootevent";
+    if (fUseOD) listWCRootEvt[1] = "wcsimrootevent_OD";
 
     // Manage input files for ID interactions 
     WCRootDataIDNuInt *daIDNuInt = new WCRootDataIDNuInt();
@@ -56,12 +79,13 @@ int main(int argc, char **argv)
     daBeamBkg->SetSeed(seed_beambkg);
 
     // Will be used to extract file ID number
-    daIDNuInt->SetNDigitsFileIDNumber(5);
-    daBeamBkg->SetNDigitsFileIDNumber(6);
+    // KMTsui: I don't have it has real usage
+    // daIDNuInt->SetNDigitsFileIDNumber(5);
+    // daBeamBkg->SetNDigitsFileIDNumber(6);
 
     // Load input files from text file
-    daIDNuInt->LoadFiles(fInFileTextIDNuInt.c_str());
-    daBeamBkg->LoadFiles(fInFileTextBeamBkg.c_str());
+    daIDNuInt->LoadFiles(fInFileTextIDNuInt.c_str(),listWCRootEvt);
+    daBeamBkg->LoadFiles(fInFileTextBeamBkg.c_str(),listWCRootEvt);
 
     // Set number of ID interactions per spill
     daIDNuInt->SetInteractionRate(fIDNuIntRate);
@@ -80,7 +104,7 @@ int main(int argc, char **argv)
     daPileUp->SetFileNameForCopyTree(daIDNuInt->GetCurrentInputFileName());
 
     // Create output file. Its name will be something like "OutFileNamePrefix".00000.root
-    daPileUp->CreateTree(fOutFileNamePrefix.c_str());
+    daPileUp->CreateTree(fOutFileNamePrefix.c_str(),listWCRootEvt);
 
     float nuIntTime; // interaction time
     int nuIntBunch;  // bunch id number 
@@ -107,7 +131,8 @@ int main(int argc, char **argv)
             bt->DrawInteractionTime(nuIntTime, nuIntBunch);
 
             // Add true hits of this interaction to MDT
-            daIDNuInt->AddTrueHitsToMDT(MDT, nuIntTime);
+            for(int k=0; k<NPMTType; k++)
+                daIDNuInt->AddTrueHitsToMDT(MDT->GetHitTubeCollection(fPMTType[k]), MDT->GetPMTResponse(fPMTType[k]), nuIntTime, k);
 
             // Add information about this interaction
             daPileUp->AddInteraction(daIDNuInt, nuIntTime, nuIntBunch);
@@ -119,18 +144,24 @@ int main(int argc, char **argv)
             daBeamBkg->Next();
 
             bt->DrawInteractionTime(nuIntTime, nuIntBunch);
-            daBeamBkg->AddTrueHitsToMDT(MDT, nuIntTime);
+            for(int k=0; k<NPMTType; k++)
+                daBeamBkg->AddTrueHitsToMDT(MDT->GetHitTubeCollection(fPMTType[k]), MDT->GetPMTResponse(fPMTType[k]), nuIntTime, k);
             daPileUp->AddInteraction(daBeamBkg, nuIntTime, nuIntBunch);
         }
 
         // Now all the true hits have been merged into one spill
         // Add dark noise hits, and then make digitized hits
-        MDT->DoAddDark();
-        MDT->DoDigitize();
-        MDT->DoTrigger();
+        for(int j=0; j<NPMTType; j++)
+        {
+            MDT->DoAddDark(fPMTType[j]);
+            MDT->DoDigitize(fPMTType[j]);
+            MDT->DoTrigger(fPMTType[j]);
+
+            // Add the resultant digitized hits to the output
+            TriggerInfo *ti = MDT->GetTriggerInfo(fPMTType[j]);
+            daPileUp->AddDigiHits(MDT->GetHitTubeCollection(fPMTType[j]), ti, i, j);
+        }
 
-        // Add the resultant digitized hits to the output
-        daPileUp->AddDigiHits(MDT);
         daPileUp->FillTree();
 
         // Clear all the true and digitized hits of this spill for the next spill
@@ -162,5 +193,6 @@ void ReadConfiguration(const char *filename)
     Conf->GetValue<float>("IDNuIntRate", fIDNuIntRate);
     Conf->GetValue<float>("BeamBkgRate", fBeamBkgRate);
     Conf->GetValue<int>("TotalNumOfSpills", fTotalNumofSpills);
+    Conf->GetValue<int>("UseOD", fUseOD);
     Conf->Finalize();
 }

app/utilities/WCRootData/include/WCRootDataBeamBkg.h

@@ -11,7 +11,7 @@ class WCRootDataBeamBkg : public WCRootDataNuInt
         WCRootDataBeamBkg();
         virtual ~WCRootDataBeamBkg();
 
-        void LoadFiles(const char *);
+        void LoadFiles(const char *, const std::vector<string> &v=vector<string>());
         float GetEnergyDepositionInOD() const { return fEdepoOD; }
 
     private :

app/utilities/WCRootData/include/WCRootDataIDNuInt.h

@@ -12,7 +12,7 @@ class WCRootDataIDNuInt : public WCRootDataNuInt
         WCRootDataIDNuInt();
         virtual ~WCRootDataIDNuInt();
 
-        void LoadFiles(const char*);
+        void LoadFiles(const char*, const std::vector<string> &v=vector<string>());
         void SetInteractionInformation(PileUpSpill_t*) const;
 
     private :

app/utilities/WCRootData/include/WCRootDataNuInt.h

@@ -22,7 +22,7 @@ class WCRootDataNuInt : public WCRootData
         virtual ~WCRootDataNuInt();
 
         void SetSeed(const int i){ fRnd->SetSeed(i); }
-        void LoadFiles(const char*); 
+        void LoadFiles(const char*, const std::vector<string> &v=vector<string>()); 
         virtual bool Next();
 
         void SetNDigitsFileIDNumber(int i){ fNDigitsFileIdNum = i; }
@@ -36,6 +36,8 @@ class WCRootDataNuInt : public WCRootData
         virtual float GetEnergyDepositionInOD() const { return 0.; }
         virtual void SetInteractionInformation(PileUpSpill_t*) const {};
 
+        void SetHasFriend(bool val) {fHasFriend=val;}
+
     protected:
         int GetFileIdNumber(const char*, const char *prefix="wcsim.");
 
@@ -58,5 +60,6 @@ class WCRootDataNuInt : public WCRootData
         int fFileIdNum; 
         int fNDigitsFileIdNum;
         float fIntRate;
+        std::vector<string> fWCRootEvtList;
 };
 #endif

app/utilities/WCRootData/include/WCRootDataPileUpSpill.h

@@ -27,7 +27,7 @@ class WCRootDataPileUpSpill : public WCRootData
 
         // For writing 
        // void SetOutFilePrefix(const char *f){ fOutFilePrefix = TString(f); }
-        void CreateTree(const char *f=0);
+        void CreateTree(const char *f=0, const vector<string> &v=vector<string>());
         void FillTree();
         void WriteTree();
         void AddInteraction(const WCRootDataNuInt*, float offset_time=0., int bunch_id=99999);

app/utilities/WCRootData/src/WCRootData.cc

@@ -83,13 +83,13 @@ void WCRootData::AddTrueHitsToMDT(HitTubeCollection *hc, PMTResponse *pr, float
             if( truetime<0. ){ continue; }
             if( aHitTime->GetParentID()<0 ){ continue; }
 
-            TrueHit *th = new TrueHit(truetime, aHitTime->GetParentID());
+            TrueHit *th = new TrueHit(truetime+intTime, aHitTime->GetParentID());
 
             for(int k=0; k<3; k++){ th->SetPosition(k, aHitTime->GetPhotonEndPos(k)); }
             for(int k=0; k<3; k++){ th->SetDirection(k, aHitTime->GetPhotonEndDir(k)); }
             for(int k=0; k<3; k++){ th->SetStartDirection(k, aHitTime->GetPhotonStartDir(k)); }
 
-            th->SetStartTime(aHitTime->GetPhotonStartTime());
+            th->SetStartTime(aHitTime->GetPhotonStartTime()+intTime);
             for(int k=0; k<3; k++){ th->SetStartPosition(k, aHitTime->GetPhotonStartPos(k)); }
             if( !pr->ApplyDE(th) ){ continue; }
 
@@ -198,6 +198,56 @@ void WCRootData::AddDigiHits(MDTManager *mdt, int eventID, int iPMT)
 void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID, int iPMT)
 {
     WCSimRootTrigger* anEvent = fSpEvt[iPMT]->GetTrigger(0);
+    // Save raw hits
+    // container for photon info
+    std::vector<double> truetime;
+    std::vector<int>   primaryParentID;
+    std::vector<float> photonStartTime;
+    std::vector<TVector3> photonStartPos;
+    std::vector<TVector3> photonEndPos;
+    std::vector<TVector3> photonStartDir;
+    std::vector<TVector3> photonEndDir;
+    for(hc->Begin(); !hc->IsEnd(); hc->Next())
+    {
+        // Get tube ID
+        HitTube *aPH = &(*hc)();
+        int tubeID = aPH->GetTubeID();
+        int mPMTID = aPH->GetmPMTID();
+        int mPMT_PMTID = aPH->GetmPMT_PMTID();
+
+
+        const int NPE = aPH->GetNRawPE();
+        const vector<TrueHit*> PEs = aPH->GetPhotoElectrons();
+        for(int iPE=0; iPE<NPE; iPE++)
+        {
+            truetime.push_back(PEs[iPE]->GetTime());
+            primaryParentID.push_back(PEs[iPE]->GetParentId());
+            photonStartTime.push_back(PEs[iPE]->GetStartTime());
+            photonStartPos.push_back(TVector3(PEs[iPE]->GetStartPosition(0),PEs[iPE]->GetStartPosition(1),PEs[iPE]->GetStartPosition(2)));
+            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)));
+        }
+
+        anEvent->AddCherenkovHit(tubeID,
+                                mPMTID,
+                                mPMT_PMTID,
+                                truetime,
+                                primaryParentID,
+                                photonStartTime,
+                                photonStartPos,
+                                photonEndPos,
+                                photonStartDir,
+                                photonEndDir);
+
+        truetime.clear();
+        primaryParentID.clear();
+        photonStartTime.clear();
+        photonStartPos.clear();
+        photonEndPos.clear();
+        photonStartDir.clear();
+        photonEndDir.clear();
+    }
     const int nTriggers = ti->GetNumOfTrigger(); 
     for(int iTrig=0; iTrig<nTriggers; iTrig++) 
     {
@@ -237,14 +287,13 @@ void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID
         {
             HitTube *aPH = &(*hc)();
             int tubeID = aPH->GetTubeID();
+            int mPMTID = aPH->GetmPMTID();
+            int mPMT_PMTID = aPH->GetmPMT_PMTID();
             int nCount = 0;
             for(int i=0; i<aPH->GetNDigiHits(); i++)
             {
                 float t = aPH->GetTimeDigi(i);
 
-                // Need to be updated 
-                int mPMT_module_id = 0;
-                int mPMT_pmt_id = 0;
                 if( t>=tWinLowEdge && t<=tWinUpEdge )
                 {
                     bool doFill = false;
@@ -262,8 +311,8 @@ void WCRootData::AddDigiHits(HitTubeCollection *hc, TriggerInfo *ti, int eventID
                         anEvent->AddCherenkovDigiHit(q,
                                                      t,
                                                      tubeID,
-                                                     mPMT_module_id,
-                                                     mPMT_pmt_id,
+                                                     mPMTID,
+                                                     mPMT_PMTID,
                                                      true_pe_comp);
                         nHits += 1;
                         sumQ += q;
@@ -477,11 +526,15 @@ void WCRootData::SetTubes(HitTubeCollection *hc, const int iPMT)
 		const WCSimRootPMT *tube = !isOD.at(iPMT) ? fWCGeom->GetPMTPtr(i, bool(iPMT)) : fWCGeom->GetODPMTPtr(i) ;
 
 		const int tubeID = tube->GetTubeNo();
+        const int mPMTID = tube->GetmPMTNo();
+        const int mPMT_PMTID = tube->GetmPMT_PMTNo();
 		hc->AddHitTube(tubeID);
 		for(int j=0; j<3; j++)
 		{
 			(&(*hc)[tubeID])->SetPosition(j, tube->GetPosition(j));
 			(&(*hc)[tubeID])->SetOrientation(j, tube->GetOrientation(j));
 		}
+        (&(*hc)[tubeID])->SetmPMTID(mPMTID);
+        (&(*hc)[tubeID])->SetmPMT_PMTID(mPMT_PMTID);
 	}
 }

app/utilities/WCRootData/src/WCRootDataBeamBkg.cc

@@ -2,18 +2,20 @@
 
 WCRootDataBeamBkg::WCRootDataBeamBkg()
 {
-    fHasFriend = true;
+    fHasFriend = false;
     fNuIntType = eBeamBkg;
 };
 
 WCRootDataBeamBkg::~WCRootDataBeamBkg()
 {
 }
 
-void WCRootDataBeamBkg::LoadFiles(const char *filename)
+void WCRootDataBeamBkg::LoadFiles(const char *filename, const vector<string> &list)
 {
     // Read input file list of WCSim first
-    this->WCRootDataNuInt::LoadFiles(filename);
+    this->WCRootDataNuInt::LoadFiles(filename,list);
+
+    if (!fHasFriend) return;
 
     ifstream fin(filename);
     string aLine;