Changeset 19951 for trunk/Mars

Timestamp:

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

Author:

giangdo

Message:

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

Location:

trunk/Mars/hawc

Files:

• callisto.C (modified) (4 diffs)
• extract_pulse.C (modified) (37 diffs)
• extract_singles.C (modified) (22 diffs)
• fit_singles.C (modified) (7 diffs)
• plot_trace.C (modified) (8 diffs)
• star.C (modified) (2 diffs)
• synchron.C (added)

trunk/Mars/hawc/callisto.C

@@ -27 +27 @@
 
  To run the macro from the command line (assuming you are in a directory
- Mars/build where you have built your Mars environment) ou can do
+ Mars/build where you have built your Mars environment) you can do
 
     root ../hawc/callisto.C\(\"20191002_018.fits.fz\",\"20191002_013.drs.fits\"\)
@@ -54 +54 @@
 
     // mapping file (found in Mars/hawc)
-    const char *mmap = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;
+    const char *mmap = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL;
 
     //const char *pulse_template = "template-pulse.root";
@@ -76 +76 @@
     // Calibration constant (for the moment a single constant to
     // convert the extracted charge to photo-electrons)
-    double scale = 1./22.553;//0.2;
+    //double scale = 1./22.553;//0.2;
+    double scale = 1./17.557;
 
     // ------------------------------------------------------
@@ -317 +318 @@
     write5.AddContainer("MSimSourcePos",       "Events", kFALSE);
     write5.AddContainer("MTime",               "Events", kFALSE);
+    write5.AddContainer("MRawBoardsFACT",      "Events", kFALSE);
     write5.AddContainer("MSrcPosCam",          "Events", kFALSE);
 

trunk/Mars/hawc/extract_pulse.C

@@ -60 +60 @@
 };
 
-//Storage class to kKeep a list of the extracted single
+//Storage class to keep a list of the extracted single
 class MSingles : public MParContainer, public MCamEvent
 {
@@ -68 +68 @@
 
 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";
-        fName = title ? title : "Storeage for vector of singles";
+        fName = title ? title : "Storage for vector of singles";
     }
 
@@ -109 +109 @@
 
 public:
-    MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10)
+    MHBaseline() : fPedestal(0), fSkipStart(5), fSkipEnd(10)
     {
         fName = "MHBaseline";
@@ -185 +185 @@
         for (int x=0; x<64; x++)
         {
+            //exclude the pixels without cones
+            if (x == 61)
+               continue;
+            if (x == 62)
+               continue;
+            if (x == 63)
+               continue;
+
             // Get the corresponding slice from the histogram
             TH1D *hist = fBaseline.ProjectionY("proj", x+1, x+1);
@@ -286 +294 @@
 class MHSingles : public MH
 {
-    TH2F fSignalAll;
+    TH2F fSignals; // all signals, without normalization
+    TH2F fOsicllation;
+
+    TH2F fSignalAll; // Histogram of all signals
     TH2F fSignalOut;
-    TH2F fSignal1;
-    TH2F fSignal2;
+    TH2F fSignal1; // Histogram of signals with 1pe
+    TH2F fSignal2; // Histogram of signals with 2pe
     TH2F fSignal3;
     TH2F fSignal4;
@@ -320 +331 @@
 
         // Setup histograms
+
+        fSignals.SetName("Signals");
+        fSignals.SetTitle("");
+        fSignals.SetXTitle("Time [ns]");
+        fSignals.SetYTitle("Amplitude [mV]");
+        fSignals.GetXaxis()->CenterTitle();
+        fSignals.GetXaxis()->SetTitleSize(0.045);
+        fSignals.GetYaxis()->SetTitleSize(0.045);
+        fSignals.SetDirectory(NULL);
+        fSignals.SetStats(kFALSE);
+
+        fOsicllation.SetName("Oscillations");
+        fOsicllation.SetTitle("");
+        fOsicllation.SetXTitle("Oscillation's Amplitude [mV]");
+        fOsicllation.SetYTitle("Pulse's Amplitude [mV]");
+        fOsicllation.GetXaxis()->CenterTitle();
+        fOsicllation.GetXaxis()->SetTitleSize(0.045);
+        fOsicllation.GetYaxis()->SetTitleSize(0.045);
+        fOsicllation.SetDirectory(NULL);
+        fOsicllation.SetStats(kFALSE);
+
         fSignalAll.SetName("SignalAll");
         fSignalAll.SetTitle("");
@@ -498 +530 @@
         //const Int_t w = fSingles->GetIntegrationWindow();
 
-        MBinning binsx, binsy;
+        MBinning binsx, binsy, bin_time, bin_pulse, binosc, binpulse;
         binsx.SetEdges(1024,  -512, 512);
         binsy.SetEdges( 800,  -20,   20);
+
+        bin_time.SetEdges(1024, -512, 512);
+        bin_pulse.SetEdges(240, -60, 60); //bin width is about 1 dac count = 0.5mV
+
+        binosc.SetEdges(140, -65, 5);
+        binpulse.SetEdges(140, -5, 65);
+
+        MH::SetBinning(fSignals, bin_time, bin_pulse);
+        MH::SetBinning(fOsicllation, binosc, binpulse);
 
         MH::SetBinning(fSignalAll, binsx, binsy);
@@ -579 +620 @@
                 continue;
 
-            if (idx!=0)
+            if (idx!=43)
                 continue;
 
-            const double gain  = hgain->GetBinContent(idx+1)/30;
+            const double gain  = hgain->GetBinContent(idx+1)/40;
             const double base  = hbase->GetBinContent(idx+1);
 
@@ -594 +635 @@
             // Fill pulse shape
             const double tm  = result[0].fTime;
-            const double sig = result[0].fSignal/30 - base;
+            const double sig = result[0].fSignal/40 - base;
 
             bool ok = true;//tm>3 && (ptr[int(tm)-3]-ped-base)/gain > -5;
 
-            double mean = 0;
-
-            for (int s=20; s<roi-20; s++)
+            double corr = -0.25;
+
+            // Make histograms to check correlation between signals and oscillations
+            for(int s=0; s<roi-10; s++)
+            {   double time = (s-tm)/2;
+                double pulse = ptr[s]-ped-base; // corrected pulse without normalization
+
+                // Fill all extracted signals into one histogram
+                fSignals.Fill(time, pulse);
+            }
+
+            int t_arrival = (int)tm;
+            double max_pulse = ptr[t_arrival]-ped-base;
+            int max_pulse_pos = 0; // position of the maximal pulse in sample
+            // Find pulse's amplitude
+            for(int s = tm; s < tm + 15; s++)
             {
-                double x = (s - tm)/2; // [ns]
+              if (ptr[s]-ped-base > max_pulse)
+                 max_pulse = ptr[s]-ped-base;
+                 max_pulse_pos = s;
+            }
+            //cout << "________Check if the filled histograms are correct________" << endl;
+            //cout << "Arrival time: " << tm << ",Amplitude: " << ptr[t_arrival]-ped-base << endl;
+            //cout << "Pulse's position: "<< max_pulse_pos << ",Max. Pulse's Amplitude: " << ptr[max_pulse_pos]-ped-base << endl;
+
+            // Find maximal oscillation's amplitude
+            double max_osc = 0;
+            for (int osc = max_pulse_pos; osc < max_pulse_pos+40; osc++)
+            {
+                if (ptr[osc]-ped-base < max_osc)
+                {
+                  max_osc = ptr[osc]-ped-base;
+                }
+                else
+                  continue;
+            }
+            fOsicllation.Fill(max_osc, max_pulse);
+
+            // Set conditions for 1pe signals
+            for (int s=5; s<roi-10; s++)
+            {
+                double x = (s - tm)/2; // return the time in [ns]
                 double y = (ptr[s]-ped-base)/gain;
 
-                if (x>-100 && x<-10)
-                    mean += y/90;
-
-                if (y<-3)
+                //if (x>-100 && x<-10)
+                //    mean += y/90;
+
+                if (x>5 && y>2.5+x*-0.01)
+                    ok = false;
+                if (y<-1.5)
                     ok = false;
                 if (x>5 && x<40 && y<=x/40*-0.5)
@@ -614 +694 @@
                 if (x>-150 && x<-5 && y>1.5)
                     ok = false;
-                if (y>2.5)
+                if (x<0 && y > 1.2)
                     ok = false;
+
             }
 
-            // Nothing at or bleow zero from 5 to 40machen
-
+            // Fill nPE signals into histograms
+
+            // Set conditions for signals with n PE
+            // Nothing at or below zero from 5 to 40machen
             bool ok1 = sig>gain*0.5 && sig<gain*1.5;
             bool ok2 = sig>gain*1.5 && sig<gain*2.5;
@@ -630 +713 @@
 
 
-            for (int s=20; s<roi-20; s++)
+            for (int s=5; s<roi-10; s++)
             {
                 double x = (s - tm)/2; // [ns]
-                double y = (ptr[s]-ped-base)/gain -  mean;
+                double y = (ptr[s]-ped-base)/gain - corr;
 
                 int ix = TMath::Nint((x-x0)/wx*nx);
@@ -642 +725 @@
 
                 int bin = (iy+1)*(nx+2)+(ix+1);
+
 
                 if (ok)
@@ -724 +808 @@
 
     // Getter for histograms
+    TH2 *GetSignals() { return &fSignals; }
+    TH2 *GetOscillation() { return &fOsicllation; }
     TH2 *GetSignalAll() { return &fSignalAll; }
     TH2 *GetSignalOut() { return &fSignalOut; }
@@ -774 +860 @@
         gPad->SetFrameBorderMode(0);
         fSignal8.Draw("colz");
-        //fProf8.Draw("same");
+        fProf8.Draw("same");
     }
 
@@ -836 +922 @@
     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)
     {
     }
@@ -877 +963 @@
 
         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;
@@ -934 +1020 @@
                 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
@@ -946 +1041 @@
                     }
                 }
+
                 // Check if the finding makes sense
                 if (k_half_max+integration_size > roi-navg-end_skip)
@@ -954 +1050 @@
                     continue;
 
-                // FIXME: Evaluate arrival time more precisely!!!
-                // FIXME: Get a better integration window
+                //cout << "________Check if the finding makes sense________" << endl;
+                //cout << "Pulse's position: " << k_max << ", Max. Pulse's Amplitude: " << val[k_max] << endl;
+                //cout << "Half maximum: " << k_half_max << ", Pulse at half maximum: " << val[k_half_max] << endl;
+
 
                 // Compile "single"
@@ -970 +1068 @@
                 // Add single to list
                 result.push_back(single);
+                break; //accept only one pulse per event
 
                 // skip falling edge
@@ -985 +1084 @@
                     Int_t        max_dist   = 14,
                     Double_t     threshold  =  5,
-                    const char  *runfile    = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_",
+                    const char  *runfile    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_",
                     int          firstRunID =  6,
                     int          lastRunID  =  6,
-                    const char  *drsfile    = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_003.drs.fits",
+                    const char  *drsfile    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_004.drs.fits",
                     const char  *outpath    = "."
                    )
 {
     // ======================================================
-    if (!pmap.Read("../hawc/FAMOUSmap171218.txt"))
-    {
-        gLog << err << "../hawc/FAMOUSmap171218.txt not found." << endl;
+    if (!pmap.Read("/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt"))
+    {
+        gLog << err << "HAWCsEyemap181214.txt not found." << endl;
         return 1;
     }
@@ -1005 +1104 @@
 
     // built output filename and path
-    TString outfile = "20191027_006_006-pulse.root";
-    TString infile  = "20190227_006_006-fit.root";
+    TString outfile = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_pulse.root";
+    TString infile  = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_fit.root";
 
     //if (!gSystem->AccessPathName(outfile))
@@ -1048 +1147 @@
 
     // 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;
 
     // ------------------------------------------------------
 
-    Long_t max0 = 0;
-    Long_t max1 = 0;
+    Long_t max0 = 10000;
+    Long_t max1 = 10000;
 
     // ======================================================
@@ -1097 +1196 @@
     MBadPixelsCam badpixels;
     badpixels.InitSize(64);
-    badpixels[7].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
-    badpixels[17].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
+    badpixels[43].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
 
     MH::SetPalette();
@@ -1219 +1317 @@
     gStyle->SetOptStat(10);
 
+    TH2 *hsignals = h->GetSignals();
+    TH2 *hoscillation = h->GetOscillation();
+
+    TH2 *hall = h->GetSignalAll();
+
     TH1 *h1 = h->GetSignal1();
     TH1 *p1 = h->GetProf1();
 
+    cout << "Entries:\t" << (p1->GetEntries()) << endl;
+
     TH1 *h2 = h->GetSignal2();
     TH1 *p2 = h->GetProf2();
@@ -1243 +1348 @@
     TH1 *p8 = h->GetProf8();
 
-    TF1 pulse("pulse", "[0]*(1-exp(-x*x/[1]))*exp(-x/[2])", 0, 300);
-    pulse.SetLineColor(kBlack);
-    pulse.SetParameters(0, 1./0.094, 1./0.05);
+    // This is the old fit
+    //TF1 pulse("pulse", "[0]*(1-exp(-x*x/[1]))*exp(-x/[2])", 0, 300);
+    //pulse.SetLineColor(kBlack);
+    //pulse.SetParameters(0, 1./0.094, 1./0.05);
+
+    TF1 pulse("pulse", "[0]*(1-1/(1+exp((x-[1])/[2])))*exp(-(x-[1])/[3])", 0, 250);
+    pulse.SetLineColor(kRed);
+    pulse.SetParLimits(1, 1.0, 3);
+    pulse.SetParLimits(2, 0.5, 2.0);
+    pulse.SetParLimits(3, 50, 110);
 
     // =====
@@ -1259 +1371 @@
 
     // =====
-
+    d->AddTab("Extracted Signals");
+    gPad->SetGrid();
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
+    gPad->SetRightMargin(0.15);
+    hsignals->DrawCopy("colz");
+
+    d->AddTab("OscillationAmpl");
+    gPad->SetGrid();
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
+    gPad->SetRightMargin(0.15);
+    hoscillation->DrawCopy("colz");
+
+
+    d->AddTab("All Signals");
+    gPad->SetGrid();
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
+    gPad->SetRightMargin(0.15);
+    hall->DrawCopy("colz");
+
+
+    // =====
     d->AddTab("1");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     h1->DrawCopy("colz");
     p1->DrawCopy("same");
 
-    pulse.SetParameter(0, 1*1.3);
-    p1->Fit(&pulse, "N0", "", 1.5, 25);
-    TF1 *f1 = pulse.DrawCopy("same");
+    //pulse.FixParameter(0, 1*1.5);
+    //p1->Fit(&pulse, "N0", "", 1, 70);
+    //TF1 *f1 = pulse.DrawCopy("same");
 
     // =====
@@ -1276 +1411 @@
     d->AddTab("2");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     h2->DrawCopy("colz");
@@ -1290 +1425 @@
     d->AddTab("3");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     h3->DrawCopy("colz");
@@ -1304 +1439 @@
     d->AddTab("4");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     h4->DrawCopy("colz");
@@ -1318 +1453 @@
     d->AddTab("5");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     h5->DrawCopy("colz");
@@ -1332 +1467 @@
     d->AddTab("6");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     h6->DrawCopy("colz");
@@ -1346 +1481 @@
     d->AddTab("7");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     h7->DrawCopy("colz");
@@ -1360 +1495 @@
     d->AddTab("8");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     h8->DrawCopy("colz");
@@ -1374 +1509 @@
     d->AddTab("Pulse");
     gPad->SetGrid();
-    gPad->SetLeftMargin(0.07);
-    gPad->SetTopMargin(0.01);
+    gPad->SetLeftMargin(0.1);
+    gPad->SetTopMargin(0.1);
     gPad->SetRightMargin(0.15);
     p8->SetStats(kFALSE);
@@ -1394 +1529 @@
     p1->DrawCopy("same");
 
-    f1->Draw("same");
+    //f1->Draw("same");
     f2->Draw("same");
     f3->Draw("same");
@@ -1406 +1541 @@
 
 
+
+    d->SaveAs(outfile);
+
+    TFile f(outfile, "UPDATE");
     /*
-    d->SaveAs(outfile);
-
-    TFile f(outfile, "UPDATE");
-
     MParameterI par("NumEvents");
     par.SetVal(fill.GetNumExecutions());

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);
 

trunk/Mars/hawc/fit_singles.C

@@ -93 +93 @@
 
 // Print function parameters
-void PrintFunc(TF1 &f, float integration_window=30)
+void PrintFunc(TF1 &f, float integration_window=40)
 {
     //cout << "--------------------" << endl;
@@ -113 +113 @@
 // The parameters for the function are the filenam from the gain extraction
 // and the output filename
-int fit_singles(const char *filename = "20191027_006_006.root",
-                const char *outfile  = "20190227_006_006-fit.root",
+int fit_singles(const char *filename = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006.root",
+                const char *outfile  = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_fit.root",
                 bool fixednoise=false)
 {
     // Read the mapping between bias channels and hardware channels
     PixelMap pmap;
-    if (!pmap.Read("../hawc/FAMOUSmap171218.txt"))
-    {
-        cout << "FAMOUSmap171218.txt not found." << endl;
+    if (!pmap.Read("/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt"))
+    {
+        cout << "HAWCsEyemap181214.txt not found." << endl;
         return 1;
     }
 
     // It is more correct to fit the integral, but this is super
-    // slow, especially for 1440 pixel. To get that, one would have
+    // slow. To get that, one would have
     // to analytically integrate and fit this function.
-    // (Currently the integral is switched off for the 1440 individual
-    // spectra in all cases)
     bool fast = false; // Switch off using integral
 
     // Values which should be read from the file but not available atm
-    Int_t integration_window = 30;
+    Int_t integration_window = 40;
 
     // Map for which pixel shall be plotted and which not
@@ -266 +264 @@
     TH1F hCoeffR2    ("CoeffR2",    "Coefficient R",          90, -1,     2);
 
-    // Histigram for the sum of all spectrums
-    TH1F hSum("Sum1", "Signal spectrum of all pixels",
+    // Histogram for the sum of all spectrums
+    TH1F hSum("Sum1", "Signal spectrum of all pixels (excl. blind pixels)",
               hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(),  hsignal->GetYaxis()->GetXmax());
     hSum.SetXTitle("pulse integral [mV sample]");
@@ -275 +273 @@
 
     // Histogram for the sum of all pixels excluding the ones with faulty fits
-    TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (incl TM)",
+    TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (excl. faulty and blind pixels)",
                     hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(),  hsignal->GetYaxis()->GetXmax());
     hSumClear1.SetXTitle("pulse integral [mV sample]");
@@ -303 +301 @@
     hPulse.SetStats(false);
 
-    // Spektrum for the normalized individual spectra
-    TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (incl TM)",  120, 0.05, 12.05);
+    // Spectrum for the normalized individual spectra
+    TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (excl. blind pixels)",  120, 0.05, 12.05);
     hSumScale1.SetXTitle("pulse integral [pe]");
     hSumScale1.SetYTitle("Rate");
@@ -332 +330 @@
         if (usePixel[pixel]==0)
             continue;
+
+        //exclude the pixels without cones
+        if (pixel == 61)
+           continue;
+        if (pixel == 62)
+           continue;
+        if (pixel == 63)
+           continue;
+
+
+
 
         TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
@@ -434 +443 @@
         if (usePixel[pixel]==0)
             continue;
+
+        // Exclude pixels without cones
+        if (pixel == 61)
+           continue;
+        if (pixel == 62)
+           continue;
+        if (pixel == 63)
+           continue;
 
         //Projectipon of a certain Pixel out of Ampl.Spectrum

trunk/Mars/hawc/plot_trace.C

@@ +1 @@
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TSystem.h>
+#include <TF1.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+#include <TMath.h>
+#include <TGraph.h>
+#include <TLine.h>
+#include <TFitResultPtr.h>
+#include <TFitResult.h>
+#include <TFile.h>
+#include <TLine.h>
+
+#include <cstdio>
+#include <stdio.h>
+#include <stdint.h>
+
+#include "Getline.h"
+#include "MH.h"
+#include "MArrayI.h"
+#include "MDirIter.h"
+#include "MFillH.h"
+#include "MEvtLoop.h"
+#include "MCamEvent.h"
+#include "MGeomApply.h"
+#include "MTaskList.h"
+#include "MParList.h"
+#include "MContinue.h"
+#include "MBinning.h"
+#include "MDrsCalibApply.h"
+#include "MDrsCalibration.h"
+#include "MRawFitsRead.h"
+#include "MStatusDisplay.h"
+#include "MTaskInteractive.h"
+#include "MPedestalSubtractedEvt.h"
+#include "MGeomCamFAMOUS.h"
+#include "MRawRunHeader.h"
+#include "MPedestalCam.h"
+#include "MPedestalPix.h"
+#include "MParameters.h"
+
 // ==========================================================================
 // ============ see plot_callisto function at the end of the file ===========
@@ -21 +63 @@
 
 // Create a histogram for the trace
-TH1D hist("Trace", "Waveform", 1024, -0.25, 511.75);
+//TH1D hist("Trace", "Waveform", 1024, -0.25, 511.75); //range = [-0.25,511.75]ns
+TH1D hist("Trace", "Waveform (Pixel 5, Event 1)", 1024, -0.5, 1024-0.5);
 
 // All 'data members' that are required globally
@@ -34 +77 @@
 {
     fEvt = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
-    if (!fEvt)
-    {
-        gfLog << err << "MPedestalSubtractedEvt not found... aborting." << endl;
-        return kALSE;
-    }
+
+    //if (!fEvt)
+    //{
+    //    gfLog << err << "MPedestalSubtractedEvt not found... aborting." << endl;
+    //    return kFALSE;
+    //}
 
     // Create a canvas and plot the histogram into it
-    c = new TCanvas;
+    c = new TCanvas();
     hist.SetStats(kFALSE);
     hist.Draw();
-
-    TF1 constant("zero", "[0]", -0.25, 511.75);
+    hist.SetXTitle("Samples");
+    hist.SetYTitle("Amplitude [a.u.]");
+
+    //TF1 constant("zero", "[0]", -0.25, 511.75);
+    TF1 constant("zero", "[0]", -0.5, 1024);
     constant.SetParameter(0, 0);
     constant.SetLineColor(kBlack);
@@ -59 +106 @@
     constant.DrawCopy("same");
 
+
     return kTRUE;
 }
@@ -64 +112 @@
 Int_t Process()
 {
-    int pixel_index = 0;
+    int pixel_index = 6; //plot trace only of this pixel
+    int skipstart = 5;
+    int skipend = 10;
 
     // Number of samples per pixel (usually 1024)
@@ -73 +123 @@
     hist.Reset();
     for (int i=0; i<numsamples; i++)
-        hist.Fill(i*0.5, cal[i]); // Convert sample to nano-second
+        hist.Fill(i, cal[i]); // 0.5ns = 1 sample
 
     // Reset min/max range
-    hist.SetMinimum();
-    hist.SetMaximum();
-
-    // Set a reasonable range
-    if (hist.GetMinimum()>-50)
-        hist.SetMinimum(-50);
-    if (hist.GetMaximum()<250)
-        hist.SetMaximum(250);
+    hist.SetMinimum(-50);
+    hist.SetMaximum(50);
+
+    //// Set a reasonable range
+    //if (hist.GetMinimum()>-30)
+    //    hist.SetMinimum(-20);
+    //if (hist.GetMaximum()<60)
+    //    hist.SetMaximum(50);
+
+
+
+    TLine *line = new TLine(skipstart, 50, skipstart, -50);
+    line->SetLineStyle(kSolid);
+    line->SetLineColor(kRed);
+    line->Draw();
+
+    TLine *line2 = new TLine(1024-skipend, 50, 1024-skipend, -50);
+    line2->SetLineStyle(kSolid);
+    line2->SetLineColor(kRed);
+    line2->Draw();
 
     // Signal root to update the canvas/pad
     c->Modified();
     c->Update();
+    //c->SaveAs("/home/giangdo/Documents/Master/MA/pulse/data_single_pe/traces/pixel3/event_0.png");
+
 
     bool rc = HandleInput();
@@ -131 +195 @@
 
     // mapping file (found in Mars/hawc)
-    const char *mmap = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;
+    const char *mmap = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL;
 
     // ======================================================
@@ -193 +257 @@
     MEvtLoop loop(gSystem->BaseName(datafile));
     loop.SetParList(&plist);
- 
-    if (!loop.Eventloop())
+
+    if (!loop.Eventloop(100))
         return 4;
 

trunk/Mars/hawc/star.C

@@ -18 +18 @@
 
  To run the macro from the command line (assuming you are in a directory
- Mars/build where you have built your Mars environment) ou can do
+ Mars/build where you have built your Mars environment) you can do
 
     root ../hawc/star.C\(\"00000015.003_Y_MonteCarlo003_Events.root\"\)
@@ -206 +206 @@
     write.AddContainer("MPointingPos",        "Events", kFALSE);
     write.AddContainer("MTime",               "Events", kFALSE);
+    write.AddContainer("MRawBoardsFACT",      "Events", kFALSE);
     write.AddContainer("MSrcPosCam",          "Events", kFALSE);