Changeset 19951 for trunk/Mars/hawc/extract_singles.C

Timestamp:

04/29/20 21:48:05 (5 years ago)

Author:

giangdo

Message:

script to synchronise the data from HAWC's Eye and HAWC

File:

• trunk/Mars/hawc/extract_singles.C (modified) (22 diffs)

trunk/Mars/hawc/extract_singles.C

@@ -57 +57 @@
 {
     Int_t fIntegrationWindow;
-
     vector<vector<Single>> fData;
 
 public:
-    MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30)
+    MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(40)
     {
         fName = name ? name : "MSingles";
@@ -88 +87 @@
     ClassDef(MSingles, 1)
 };
-
 // Histogram class to extract the baseline
 class MHBaseline : public MH
@@ -102 +100 @@
 
 public:
-    MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10)
+    MHBaseline() : fPedestal(0), fSkipStart(10), fSkipEnd(10)
     {
         fName = "MHBaseline";
@@ -108 +106 @@
         // Setup the histogram
         fBaseline.SetName("Baseline");
-        fBaseline.SetTitle("Median spectrum");
+        fBaseline.SetTitle("Baseline Spectrum");
         fBaseline.SetXTitle("Pixel [idx]");
-        fBaseline.SetYTitle("Median baseline [mV]");
+        fBaseline.SetYTitle("Amplitude [mV]");
         fBaseline.SetDirectory(NULL);
     }
@@ -282 +280 @@
     TH2F       fTime;
     TProfile2D fPulse;
+    TH2F       fPulse2D; //plot the extracted singles in a 2D plot vs samples(time)
 
     UInt_t fNumEntries;
@@ -296 +295 @@
         // Setup histograms
         fSignal.SetName("Signal");
-        fSignal.SetTitle("pulse integral spectrum");
+        fSignal.SetTitle("Pulse Integral Spectrum");
         fSignal.SetXTitle("pixel [SoftID]");
         fSignal.SetYTitle("time [a.u.]");
@@ -302 +301 @@
 
         fTime.SetName("Time");
-        fTime.SetTitle("arival time spectrum");
+        fTime.SetTitle("Arrival Time Spectrum");
         fTime.SetXTitle("pixel [SoftID]");
         fTime.SetYTitle("time [a.u.]");
@@ -308 +307 @@
 
         fPulse.SetName("Pulse");
-        fPulse.SetTitle("average pulse shape spectrum");
+        fPulse.SetTitle("Average Pulse Shape Spectrum");
         fPulse.SetXTitle("pixel [SoftID]");
         fPulse.SetYTitle("time [a.u.]");
         fPulse.SetDirectory(NULL);
+
+        fPulse2D.SetName("Pulse2D");
+        fPulse2D.SetTitle("Extracted Pulses");
+        fPulse2D.SetXTitle("time [a.u.]");
+        fPulse2D.SetYTitle("amplitude [a.u.]");
+        fPulse2D.SetDirectory(NULL);
     }
 
@@ -363 +368 @@
 
         // Correct for DRS timing!!
+
+        //Setup binning for the average time spectrum
         MBinning binst(roi, -0.5, roi-0.5);
         MH::SetBinning(fTime, binsx, binst);
 
+        //Setup binning for the average pulse spectrum
         MBinning binspy(2*roi, -roi-0.5, roi-0.5);
+        //MBinning binspy(1024, -522-0.5, 522-0.5);
         MH::SetBinning(fPulse, binsx, binspy);
+
+        //Setup binning for the 2D plot amplitude vs time
+        MBinning bins_time(2*roi, -roi-0.5, roi-0.5);
+        //MBinning bins_time(1024, -522-0.5, 522-0.5);
+        MBinning bins_ampl(160, -30, 50); //N_bin is chosen such that 1 DAC = 0.5mV
+        MH::SetBinning(fPulse2D, bins_time, bins_ampl);
 
         return kTRUE;
@@ -402 +417 @@
                 for (int s=0; s<roi; s++)
                     fPulse.Fill(i, s-tm, ptr[s]);
+
+                for (int s=0; s<roi; s++)
+                    fPulse2D.Fill(s-tm, ptr[s]);
             }
 
@@ -416 +434 @@
     const TH2 *GetTime() const   { return &fTime; }
     const TH2 *GetPulse() const  { return &fPulse; }
+    const TH2F *GetPulse2D() const {return &fPulse2D;}
 
     UInt_t GetNumEntries() const { return fNumEntries; }
@@ -441 +460 @@
         gPad->SetFrameBorderMode(0);
         fPulse.Draw("colz");
+
+        pad->cd(4);
+        gPad->SetBorderMode(0);
+        gPad->SetFrameBorderMode(0);
+        fPulse2D.Draw("colz");
     }
 
@@ -465 +489 @@
         gPad->SetFrameBorderMode(0);
         fPulse.DrawCopy("colz");
+
+        pad->cd(4);
+        gPad->SetBorderMode(0);
+        gPad->SetFrameBorderMode(0);
+        fPulse2D.Draw("colz");
     }
 
@@ -497 +526 @@
     MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0),
         fExtractionRange(64),
-        fNumAverage(10), fStartSkip(20), fEndSkip(10),
-        fIntegrationSize(3*10), fMaxSearchWindow(30)
+        fNumAverage(10), fStartSkip(10), fEndSkip(10),
+        fIntegrationSize(4*10), fMaxSearchWindow(20)
     {
     }
@@ -538 +567 @@
 
         const size_t navg              = fNumAverage;
-        const float  threshold         = fThreshold; 
+        const float  threshold         = fThreshold;
         const UInt_t start_skip        = fStartSkip;
         const UInt_t end_skip          = fEndSkip;
@@ -595 +624 @@
                 if (k_max == i+5 || k_max == i + max_search_window-1)
                     continue;
+
+                //Make sure the pulse is isolated (area after the maximum is empty)
+                UInt_t k_falling = k_max+200;
+                for (UInt_t k_after = k_falling; k_after < k_max +250; k_after++)
+                {
+                    if (val[k_after] > val[k_falling] &&
+                        val[k_after + fNumAverage/2] > val[k_falling])
+                    continue;
+                }
 
                 //search for half maximum before maximum
@@ -623 +661 @@
                 single.fTime   = k_half_max;
 
-                // Crossing of the threshold is at 5
+                // Crossing of the threshold is at 7
                 for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
                     single.fSignal += ptr[j];
@@ -631 +669 @@
                 // Add single to list
                 result.push_back(single);
+                break; //accept only one pulse per event
 
                 // skip falling edge
@@ -644 +683 @@
 
 int extract_singles(
-                    Int_t        max_dist   = 14,
-                    Double_t     threshold  =  5,
-                    const char  *runfile    = "raw/2019/10/27/20191027_",
-                    int          firstRunID =  6,
-                    int          lastRunID  =  6,
-                    const char  *drsfile    = "raw/2019/10/27/20191027_003.drs.fits",
-                    const char  *outpath    = "."
+                      Int_t        max_dist   = 14,
+                      Double_t     threshold  = 5,
+                      const char  *runfile    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_",
+                      int          firstRunID =  6,
+                      int          lastRunID  =  6,
+                      const char  *drsfile    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_004.drs.fits",
+                      const char  *outpath    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted"
                    )
 {
@@ -688 +727 @@
 
     // map file to use (get that from La Palma!)
-    const char *map = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;
+    const char *map = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL;
 
     // ------------------------------------------------------
@@ -856 +895 @@
     const TH2 *htime   = h->GetTime();
     const TH2 *hpulse  = h->GetPulse();
+    const TH2F *ppulse2d = h->GetPulse2D();
 
     d->AddTab("Time");
@@ -869 +909 @@
     hp->Draw();
 
+    d->AddTab("2DPulse");
+    ppulse2d->DrawCopy("colz");
+
     d->SaveAs(outfile);