Index: /trunk/Mars/hawc/extract_pulse.C =================================================================== --- /trunk/Mars/hawc/extract_pulse.C (revision 19876) +++ /trunk/Mars/hawc/extract_pulse.C (revision 19876) @@ -0,0 +1,1427 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "MH.h" +#include "MArrayI.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MDirIter.h" +#include "MFillH.h" +#include "MEvtLoop.h" +#include "MCamEvent.h" +#include "MHCamEvent.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 "MBadPixelsCam.h" +#include "MStatusArray.h" +#include "MStatusDisplay.h" +#include "MTaskInteractive.h" +#include "MPedestalSubtractedEvt.h" +#include "MHCamera.h" +#include "MGeomCamFAMOUS.h" +#include "MRawRunHeader.h" +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MParameters.h" +#include "PixelMap.h" + +using namespace std; +using namespace TMath; + +MHCamera *hgain = 0; +MHCamera *hbase = 0; + +PixelMap pmap; + +// Structure to store the extracted value and the extracted time +struct Single +{ + float fSignal; + float fTime; +}; + +//Storage class to kKeep a list of the extracted single +class MSingles : public MParContainer, public MCamEvent +{ + Int_t fIntegrationWindow; + + vector> fData; + +public: + MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30) + { + fName = name ? name : "MSingles"; + fName = title ? title : "Storeage for vector of singles"; + } + + void InitSize(const UInt_t i) + { + fData.resize(i); + } + + vector &operator[](UInt_t i) { return fData[i]; } + vector &GetVector(UInt_t i) { return fData[i]; } + + UInt_t GetNumPixels() const { return fData.size(); } + + void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; } + Int_t GetIntegrationWindow() const { return fIntegrationWindow; } + + Bool_t GetPixelContent(Double_t &, Int_t , const MGeomCam &, Int_t) const + { + return kTRUE; + } + void DrawPixelContent(Int_t) const { } + + ClassDef(MSingles, 1) +}; + +// Histogram class to extract the baseline +class MHBaseline : public MH +{ + TH2F fBaseline; + + MPedestalCam *fPedestal; + + // The baseline is only extracted where also the signal is extracted + // FIXME: Make sure this is consistent with MExtractSingles + UShort_t fSkipStart; + UShort_t fSkipEnd; + +public: + MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10) + { + fName = "MHBaseline"; + + // Setup the histogram + fBaseline.SetName("Baseline"); + fBaseline.SetTitle("Median spectrum"); + fBaseline.SetXTitle("Pixel [idx]"); + fBaseline.SetYTitle("Median baseline [mV]"); + fBaseline.SetDirectory(NULL); + } + + Bool_t ReInit(MParList *plist) + { + fPedestal = (MPedestalCam*)plist->FindCreateObj("MPedestalCam"); + if (!fPedestal) + return kFALSE; + + const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader"); + if (!header) + { + *fLog << err << "MRawRunHeader not found... abort." << endl; + return kFALSE; + } + + // Bin width should be around 1 dac count which is about 0.5mV + MBinning binsx, binsy; + binsx.SetEdges(64, -0.5, 64-0.5); + binsy.SetEdges(100, -20.5, 29.5); + + // Setup binnning + MH::SetBinning(fBaseline, binsx, binsy); + + return kTRUE; + } + + // Fill the samples into the histogram + Int_t Fill(const MParContainer *par, const Stat_t) + { + const MPedestalSubtractedEvt *evt = dynamic_cast(par); + + const Int_t n = evt->GetNumSamples()-fSkipStart-fSkipEnd; + + // Loop over all pixels + for (int pix=0; pix<64; pix++) + { + // Get samples for each pixel + const Float_t *ptr = evt->GetSamples(pix); + + // Average two consecutive samples + for (int i=0; iInitSize(fBaseline.GetNbinsX()); + fPedestal->SetNumEvents(GetNumExecutions()); + + Int_t cnt = 0; + Double_t avg = 0; + Double_t rms = 0; + + // Loop over all 'pixels' + for (int x=0; x<64; x++) + { + // Get the corresponding slice from the histogram + TH1D *hist = fBaseline.ProjectionY("proj", x+1, x+1); + + // Get the maximum bin + const Int_t bin = hist->GetMaximumBin(); + + // Calculate a parabola through this and the surrounding points + // on logarithmic values (that's a gaussian) + + // + // Quadratic interpolation + // + // calculate the parameters of a parabula such that + // y(i) = a + b*x(i) + c*x(i)^2 + // y'(i) = b + 2*c*x(i) + // + // + + // ------------------------------------------------------------------------- + // a = y2; + // b = (y3-y1)/2; + // c = (y3+y1)/2 - y2; + + const Double_t v1 = hist->GetBinContent(bin-1); + const Double_t v2 = hist->GetBinContent(bin); + const Double_t v3 = hist->GetBinContent(bin+1); + if (v1<=0 || v2<=0 || v3<=0) + continue; + + const Double_t y1 = TMath::Log(v1); + const Double_t y2 = TMath::Log(v2); + const Double_t y3 = TMath::Log(v3); + + //Double_t a = y2; + const Double_t b = (y3-y1)/2; + const Double_t c = (y1+y3)/2 - y2; + if (c>=0) + continue; + + const Double_t w = -1./(2*c); + const Double_t dx = b*w; + + if (dx<-1 || dx>1) + continue; + + // y = exp( - (x-k)^2 / s^2 / 2 ) + // + // -2*s^2 * log(y) = x^2 - 2*k*x + k^2 + // + // a = (k/s0)^2/2 + // b = k/s^2 + // c = -1/(2s^2) --> s = sqrt(-1/2c) + + const Double_t binx = hist->GetBinCenter(bin); + const Double_t binw = hist->GetBinWidth(bin); + + //const Double_t p = hist->GetBinCenter(hist->GetMaximumBin()); + const Double_t p = binx + dx*binw; + + avg += p; + rms += p*p; + + cnt++; + + // Store baseline and sigma + MPedestalPix &pix = (*fPedestal)[x]; + + pix.SetPedestal(p); + pix.SetPedestalRms(TMath::Sqrt(w)*binw); + + delete hist; + } + + avg /= cnt; + rms /= cnt; + + cout << "Baseline(" << cnt << "): " << avg << " +- " << sqrt(rms-avg*avg) << endl; + + return kTRUE; + } + + // Draw histogram + void Draw(Option_t *) + { + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + AppendPad(""); + + pad->SetBorderMode(0); + pad->SetFrameBorderMode(0); + fBaseline.Draw("colz"); + } + + + ClassDef(MHBaseline, 1); +}; + +// Histogram class for the signal and time distribution as +// well as the pulse shape +class MHSingles : public MH +{ + TH2F fSignalAll; + TH2F fSignalOut; + TH2F fSignal1; + TH2F fSignal2; + TH2F fSignal3; + TH2F fSignal4; + TH2F fSignal5; + TH2F fSignal6; + TH2F fSignal7; + TH2F fSignal8; + + TProfile fProfAll; + TProfile fProf1; + TProfile fProf2; + TProfile fProf3; + TProfile fProf4; + TProfile fProf5; + TProfile fProf6; + TProfile fProf7; + TProfile fProf8; + + UInt_t fNumEntries; + + MSingles *fSingles; //! + MPedestalSubtractedEvt *fEvent; //! + MBadPixelsCam *fBadPix; //! + MPedestalCam *fPedestal; //! + +public: + MHSingles() : fNumEntries(0), fSingles(0), fEvent(0), fPedestal(0) + { + fName = "MHSingles"; + + // Setup histograms + fSignalAll.SetName("SignalAll"); + fSignalAll.SetTitle(""); + fSignalAll.SetXTitle("Time [ns]"); + fSignalAll.SetYTitle("Amplitude [a.u.]"); + fSignalAll.GetXaxis()->CenterTitle(); + fSignalAll.GetXaxis()->SetTitleSize(0.045); + fSignalAll.GetYaxis()->SetTitleSize(0.045); + fSignalAll.SetDirectory(NULL); + fSignalAll.SetStats(kFALSE); + + fSignalOut.SetName("SignalOut"); + fSignalOut.SetTitle(""); + fSignalOut.SetXTitle("Time [ns]"); + fSignalOut.SetYTitle("Amplitude [a.u.]"); + fSignalOut.GetXaxis()->CenterTitle(); + fSignalOut.GetXaxis()->SetTitleSize(0.045); + fSignalOut.GetYaxis()->SetTitleSize(0.045); + fSignalOut.SetDirectory(NULL); + fSignalOut.SetStats(kFALSE); + + fSignal1.SetName("Signal1"); + fSignal1.SetTitle(""); + fSignal1.SetXTitle("Time [ns]"); + fSignal1.SetYTitle("Amplitude [a.u.]"); + fSignal1.GetXaxis()->CenterTitle(); + fSignal1.GetXaxis()->SetTitleSize(0.045); + fSignal1.GetYaxis()->SetTitleSize(0.045); + fSignal1.SetDirectory(NULL); + fSignal1.SetStats(kFALSE); + + fSignal2.SetName("Signal2"); + fSignal2.SetTitle(""); + fSignal2.SetXTitle("Time [ns]"); + fSignal2.SetYTitle("Amplitude [a.u.]"); + fSignal2.GetXaxis()->CenterTitle(); + fSignal2.GetXaxis()->SetTitleSize(0.045); + fSignal2.GetYaxis()->SetTitleSize(0.045); + fSignal2.SetDirectory(NULL); + fSignal2.SetStats(kFALSE); + + fSignal3.SetName("Signal3"); + fSignal3.SetTitle(""); + fSignal3.SetXTitle("Time [ns]"); + fSignal3.SetYTitle("Amplitude [a.u.]"); + fSignal3.GetXaxis()->CenterTitle(); + fSignal3.GetXaxis()->SetTitleSize(0.045); + fSignal3.GetYaxis()->SetTitleSize(0.045); + fSignal3.SetDirectory(NULL); + fSignal3.SetStats(kFALSE); + + fSignal4.SetName("Signal4"); + fSignal4.SetTitle(""); + fSignal4.SetXTitle("Time [ns]"); + fSignal4.SetYTitle("Amplitude [a.u.]"); + fSignal4.GetXaxis()->CenterTitle(); + fSignal4.GetXaxis()->SetTitleSize(0.045); + fSignal4.GetYaxis()->SetTitleSize(0.045); + fSignal4.SetDirectory(NULL); + fSignal4.SetStats(kFALSE); + + fSignal5.SetName("Signal5"); + fSignal5.SetTitle(""); + fSignal5.SetXTitle("Time [ns]"); + fSignal5.SetYTitle("Amplitude [a.u.]"); + fSignal5.GetXaxis()->CenterTitle(); + fSignal5.GetXaxis()->SetTitleSize(0.045); + fSignal5.GetYaxis()->SetTitleSize(0.045); + fSignal5.SetDirectory(NULL); + fSignal5.SetStats(kFALSE); + + fSignal6.SetName("Signal6"); + fSignal6.SetTitle(""); + fSignal6.SetXTitle("Time [ns]"); + fSignal6.SetYTitle("Amplitude [a.u.]"); + fSignal6.GetXaxis()->CenterTitle(); + fSignal6.GetXaxis()->SetTitleSize(0.045); + fSignal6.GetYaxis()->SetTitleSize(0.045); + fSignal6.SetDirectory(NULL); + fSignal6.SetStats(kFALSE); + + fSignal7.SetName("Signal7"); + fSignal7.SetTitle(""); + fSignal7.SetXTitle("Time [ns]"); + fSignal7.SetYTitle("Amplitude [a.u.]"); + fSignal7.GetXaxis()->CenterTitle(); + fSignal7.GetXaxis()->SetTitleSize(0.045); + fSignal7.GetYaxis()->SetTitleSize(0.045); + fSignal7.SetDirectory(NULL); + fSignal7.SetStats(kFALSE); + + fSignal8.SetName("Signal8"); + fSignal8.SetTitle(""); + fSignal8.SetXTitle("Time [ns]"); + fSignal8.SetYTitle("Amplitude [a.u.]"); + fSignal8.GetXaxis()->CenterTitle(); + fSignal8.GetXaxis()->SetTitleSize(0.045); + fSignal8.GetYaxis()->SetTitleSize(0.045); + fSignal8.SetDirectory(NULL); + fSignal8.SetStats(kFALSE); + + fProfAll.SetName("ProfAll"); + fProf1.SetName("Prof1"); + fProf2.SetName("Prof2"); + fProf3.SetName("Prof3"); + fProf4.SetName("Prof4"); + fProf5.SetName("Prof5"); + fProf6.SetName("Prof6"); + fProf7.SetName("Prof7"); + fProf8.SetName("Prof8"); + + fProfAll.SetErrorOption("s"); + fProf1.SetErrorOption("s"); + fProf2.SetErrorOption("s"); + fProf3.SetErrorOption("s"); + fProf4.SetErrorOption("s"); + fProf5.SetErrorOption("s"); + fProf6.SetErrorOption("s"); + fProf7.SetErrorOption("s"); + fProf8.SetErrorOption("s"); + + fProfAll.GetYaxis()->SetTitleOffset(0.8); + fProf1.GetYaxis()->SetTitleOffset(0.8); + fProf2.GetYaxis()->SetTitleOffset(0.8); + fProf3.GetYaxis()->SetTitleOffset(0.8); + fProf4.GetYaxis()->SetTitleOffset(0.8); + fProf5.GetYaxis()->SetTitleOffset(0.8); + fProf6.GetYaxis()->SetTitleOffset(0.8); + fProf7.GetYaxis()->SetTitleOffset(0.8); + fProf8.GetYaxis()->SetTitleOffset(0.8); + } + + Bool_t SetupFill(const MParList *plist) + { + fSingles = (MSingles*)plist->FindObject("MSingles"); + if (!fSingles) + { + *fLog << err << "MSingles not found... abort." << endl; + return kFALSE; + } + + fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); + if (!fEvent) + { + *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl; + return kFALSE; + } + + fBadPix = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam"); + if (!fBadPix) + { + *fLog << err << "MBadPixelsCam not found... abort." << endl; + return kFALSE; + } + + fPedestal = (MPedestalCam*)plist->FindObject("MPedestalCam"); + if (!fPedestal) + { + *fLog << err << "MPedestalCam not found... abort." << endl; + return kFALSE; + } + + fNumEntries = 0; + + return kTRUE; + } + + Bool_t ReInit(MParList *plist) + { + const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader"); + if (!header) + { + *fLog << err << "MRawRunHeader not found... abort." << endl; + return kFALSE; + } + + // Setup binning + //const Int_t w = fSingles->GetIntegrationWindow(); + + MBinning binsx, binsy; + binsx.SetEdges(1024, -512, 512); + binsy.SetEdges( 800, -20, 20); + + MH::SetBinning(fSignalAll, binsx, binsy); + MH::SetBinning(fSignalOut, binsx, binsy); + MH::SetBinning(fSignal1, binsx, binsy); + MH::SetBinning(fSignal2, binsx, binsy); + MH::SetBinning(fSignal3, binsx, binsy); + MH::SetBinning(fSignal4, binsx, binsy); + MH::SetBinning(fSignal5, binsx, binsy); + MH::SetBinning(fSignal6, binsx, binsy); + MH::SetBinning(fSignal7, binsx, binsy); + MH::SetBinning(fSignal8, binsx, binsy); + MH::SetBinning(fProfAll, binsx); + MH::SetBinning(fProf1, binsx); + MH::SetBinning(fProf2, binsx); + MH::SetBinning(fProf3, binsx); + MH::SetBinning(fProf4, binsx); + MH::SetBinning(fProf5, binsx); + MH::SetBinning(fProf6, binsx); + MH::SetBinning(fProf7, binsx); + MH::SetBinning(fProf8, binsx); + + fSignal1.GetXaxis()->SetRangeUser(-25, 300); + fSignal1.GetYaxis()->SetRangeUser(-1.5, 2.5); + + fSignal2.GetXaxis()->SetRangeUser(-25, 300); + fSignal2.GetYaxis()->SetRangeUser(-1.5, 3.5); + + fSignal3.GetXaxis()->SetRangeUser(-25, 300); + fSignal3.GetYaxis()->SetRangeUser(-1.5, 4.5); + + fSignal4.GetXaxis()->SetRangeUser(-25, 300); + fSignal4.GetYaxis()->SetRangeUser(-1.5, 5.5); + + fSignal5.GetXaxis()->SetRangeUser(-25, 300); + fSignal5.GetYaxis()->SetRangeUser(-1.5, 6.5); + + fSignal6.GetXaxis()->SetRangeUser(-25, 300); + fSignal6.GetYaxis()->SetRangeUser(-1.5, 7.5); + + fSignal7.GetXaxis()->SetRangeUser(-25, 300); + fSignal7.GetYaxis()->SetRangeUser(-1.5, 8.5); + + fSignal8.GetXaxis()->SetRangeUser(-25, 300); + fSignal8.GetYaxis()->SetRangeUser(-1.5, 9.5); + + return kTRUE; + } + + // Fill singles into histograms + Int_t Fill(const MParContainer *, const Stat_t) + { + // Get pointer to samples to fill samples + const Float_t *ptr = fEvent->GetSamples(0); + const Short_t roi = fEvent->GetNumSamples(); + + if (!ptr) + return kTRUE; + + const Int_t nx = fSignalAll.GetNbinsX(); + const Int_t ny = fSignalAll.GetNbinsY(); + + const double x0 = fSignalAll.GetXaxis()->GetBinCenter(1); + const double y0 = fSignalAll.GetYaxis()->GetBinCenter(1); + + const double wx = fSignalAll.GetXaxis()->GetBinCenter(nx)-x0; + const double wy = fSignalAll.GetYaxis()->GetBinCenter(ny)-y0; + + // Loop over all pixels + for (unsigned int i=0; i<64; i++, ptr+=roi) + { + int idx = i;//pmap.hw(i).index; + + if ((*fBadPix)[i].IsUnsuitable()) + continue; + + if (!hgain->IsUsed(idx)) + continue; + + if (idx!=0) + continue; + + const double gain = hgain->GetBinContent(idx+1)/30; + const double base = hbase->GetBinContent(idx+1); + + const double ped = (*fPedestal)[i].GetPedestal(); + + // loop over all singles + const vector &result = fSingles->GetVector(i); + if (result.size()!=1) + continue; + + // Fill pulse shape + const double tm = result[0].fTime; + const double sig = result[0].fSignal/30 - base; + + bool ok = true;//tm>3 && (ptr[int(tm)-3]-ped-base)/gain > -5; + + double mean = 0; + + for (int s=20; s-100 && x<-10) + mean += y/90; + + if (y<-3) + ok = false; + if (x>5 && x<40 && y<=x/40*-0.5) + ok = false; + if (x>-150 && x<-5 && y>1.5) + ok = false; + if (y>2.5) + ok = false; + } + + // Nothing at or bleow zero from 5 to 40machen + + bool ok1 = sig>gain*0.5 && siggain*1.5 && siggain*2.5 && siggain*3.5 && siggain*4.5 && siggain*5.5 && siggain*6.5 && siggain*7.5 && sig=nx || iy>=ny) + continue; + + int bin = (iy+1)*(nx+2)+(ix+1); + + if (ok) + { + fSignalAll.AddBinContent(bin); + fProfAll.Fill(x, y); + + if (ok1) + { + fSignal1.AddBinContent(bin); + fProf1.Fill(x, y); + } + if (ok2) + { + fSignal2.AddBinContent(bin); + fProf2.Fill(x, y); + } + if (ok3) + { + fSignal3.AddBinContent(bin); + fProf3.Fill(x, y); + } + if (ok4) + { + fSignal4.AddBinContent(bin); + fProf4.Fill(x, y); + } + if (ok5) + { + fSignal5.AddBinContent(bin); + fProf5.Fill(x, y); + } + if (ok6) + { + fSignal6.AddBinContent(bin); + fProf6.Fill(x, y); + } + if (ok7) + { + fSignal7.AddBinContent(bin); + fProf7.Fill(x, y); + } + if (ok8) + { + fSignal8.AddBinContent(bin); + fProf8.Fill(x, y); + } + } + else + fSignalOut.AddBinContent(bin); + } + + if (ok) + { + fSignalAll.SetEntries(fSignalAll.GetEntries()+1); + + if (ok1) + fSignal1.SetEntries(fSignal1.GetEntries()+1); + if (ok2) + fSignal2.SetEntries(fSignal2.GetEntries()+1); + if (ok3) + fSignal3.SetEntries(fSignal3.GetEntries()+1); + if (ok4) + fSignal4.SetEntries(fSignal4.GetEntries()+1); + if (ok5) + fSignal5.SetEntries(fSignal5.GetEntries()+1); + if (ok6) + fSignal6.SetEntries(fSignal6.GetEntries()+1); + if (ok7) + fSignal7.SetEntries(fSignal7.GetEntries()+1); + if (ok8) + fSignal8.SetEntries(fSignal8.GetEntries()+1); + } + else + fSignalOut.SetEntries(fSignalOut.GetEntries()+1); + } + + fNumEntries++; + + return kTRUE; + } + + // Getter for histograms + TH2 *GetSignalAll() { return &fSignalAll; } + TH2 *GetSignalOut() { return &fSignalOut; } + TH2 *GetSignal1() { return &fSignal1; } + TH2 *GetSignal2() { return &fSignal2; } + TH2 *GetSignal3() { return &fSignal3; } + TH2 *GetSignal4() { return &fSignal4; } + TH2 *GetSignal5() { return &fSignal5; } + TH2 *GetSignal6() { return &fSignal6; } + TH2 *GetSignal7() { return &fSignal7; } + TH2 *GetSignal8() { return &fSignal8; } + TH1 *GetProf1() { return &fProf1; } + TH1 *GetProf2() { return &fProf2; } + TH1 *GetProf3() { return &fProf3; } + TH1 *GetProf4() { return &fProf4; } + TH1 *GetProf5() { return &fProf5; } + TH1 *GetProf6() { return &fProf6; } + TH1 *GetProf7() { return &fProf7; } + TH1 *GetProf8() { return &fProf8; } + + UInt_t GetNumEntries() const { return fNumEntries; } + + void Draw(Option_t *) + { + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + AppendPad(""); + + pad->Divide(2,2); + + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignalAll.Draw("colz"); + fProfAll.Draw("same"); + + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignalOut.Draw("colz"); + + pad->cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignal1.Draw("colz"); + fProf1.Draw("same"); + + pad->cd(4); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignal8.Draw("colz"); + //fProf8.Draw("same"); + } + + void DrawCopy(Option_t *) + { + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + AppendPad(""); + + pad->Divide(2,2); + + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignalAll.DrawCopy("colz"); + + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignalOut.DrawCopy("colz"); + + pad->cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignal1.DrawCopy("colz"); + + pad->cd(4); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignal8.DrawCopy("colz"); + } + + ClassDef(MHSingles, 2) +}; + +// Task to extract the singles +class MExtractSingles : public MTask +{ + MSingles *fSingles; + MPedestalCam *fPedestal; + MPedestalSubtractedEvt *fEvent; + + // On time for each pixel in samples + MArrayI fExtractionRange; + + // Number of samples for sliding average + size_t fNumAverage; + + // The range in which the singles have to fit in + // FIXME: Make sure this is consistent with MHBaseline + UInt_t fStartSkip; + UInt_t fEndSkip; + + UInt_t fIntegrationSize; + UInt_t fMaxSearchWindow; + + Int_t fMaxDist; + Double_t fThreshold; + +public: + MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0), + fExtractionRange(64), + fNumAverage(10), fStartSkip(20), fEndSkip(10), + fIntegrationSize(3*10), fMaxSearchWindow(30) + { + } + + void SetMaxDist(Int_t m) { fMaxDist = m; } + void SetThreshold(Int_t th) { fThreshold = th; } + + UInt_t GetIntegrationSize() const { return fIntegrationSize; } + + const MArrayI &GetExtractionRange() const { return fExtractionRange; } + + + Int_t PreProcess(MParList *plist) + { + fSingles = (MSingles*)plist->FindCreateObj("MSingles"); + if (!fSingles) + return kFALSE; + + fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); + if (!fEvent) + { + *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl; + return kFALSE; + } + + fPedestal = (MPedestalCam*)plist->FindObject("MPedestalCam"); + if (!fPedestal) + { + *fLog << err << "MPedestalCam not found... abort." << endl; + return kFALSE; + } + + return kTRUE; + } + + Int_t Process() + { + const UInt_t roi = fEvent->GetNumSamples(); + + const size_t navg = fNumAverage; + const float threshold = fThreshold; + const UInt_t start_skip = fStartSkip; + const UInt_t end_skip = fEndSkip; + const UInt_t integration_size = fIntegrationSize; + const UInt_t max_search_window = fMaxSearchWindow; + const UInt_t max_dist = fMaxDist; + + vector val(roi-navg);//result of first sliding average + for (int pix=0; pix<64; pix++) + { + // Total number of samples as 'on time' + fExtractionRange[pix] += roi-start_skip-navg-end_skip-integration_size; + + // Clear singles for this pixel + vector &result = fSingles->GetVector(pix); + result.clear(); + + // Get pointer to samples + const Float_t *ptr = fEvent->GetSamples(pix); + + // Get Baseline + const Double_t ped = (*fPedestal)[pix].GetPedestal(); + + // Subtract baseline and do a sliding average over + // the samples to remove noise (mainly the 200MHz noise) + for (size_t i=0; ithreshold || + val[i+4]>threshold || + + val[i+5] val[k_max]) + k_max = k; + } + + // Check if the findings make sense + if (k_max == i+5 || k_max == i + max_search_window-1) + continue; + + //search for half maximum before maximum + UInt_t k_half_max = 0; + for (UInt_t k=k_max; k>k_max-25; k--) + { + if (val[k-1] < val[k_max]/2 && + val[k] >= val[k_max]/2 ) + { + k_half_max = k; + break; + } + } + // Check if the finding makes sense + if (k_half_max+integration_size > roi-navg-end_skip) + continue; + if (k_half_max == 0) + continue; + if (k_max - k_half_max > max_dist) + continue; + + // FIXME: Evaluate arrival time more precisely!!! + // FIXME: Get a better integration window + + // Compile "single" + Single single; + single.fSignal = 0; + single.fTime = k_half_max; + + // Crossing of the threshold is at 5 + for (UInt_t j=k_half_max; jimax ? 5+integration_size-(i-imax) : 5+integration_size; + } + } + return kTRUE; + } +}; + +int extract_pulse( + 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_", + 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 *outpath = "." + ) +{ + // ====================================================== + if (!pmap.Read("../hawc/FAMOUSmap171218.txt")) + { + gLog << err << "../hawc/FAMOUSmap171218.txt not found." << endl; + return 1; + } + + // ====================================================== + + + MDirIter iter; + + // built output filename and path + TString outfile = "20191027_006_006-pulse.root"; + TString infile = "20190227_006_006-fit.root"; + + //if (!gSystem->AccessPathName(outfile)) + // return 0; + + // add input files to directory iterator + for (Int_t fileNr = firstRunID; fileNr <= lastRunID; fileNr++) + { + TString filename = runfile; + filename += Form("%03d.fits", fileNr); + iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".fz")); + iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".gz")); + iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename)); + } + + // ====================================================== + + TFile file(infile); + MStatusArray arr; + if (arr.Read()<=0) + { + gLog << "ERROR - Cannot open gain file '" << infile << "'" << endl; + return 1; + } + + hgain = (MHCamera*)arr.FindObjectInCanvas("Gain_copy", "MHCamera", "Cams1"); + hbase = (MHCamera*)arr.FindObjectInCanvas("Baseline_copy", "MHCamera", "Cams1"); + if (!hgain || !hbase) + { + gLog << "ERROR - Cannot find Gain/Baseline" << endl; + return 1; + } + + // ====================================================== + + // true: Display correctly mapped pixels in the camera displays + // but the value-vs-index plot is in software/spiral indices + // false: Display pixels in hardware/linear indices, + // but the order is the camera display is distorted + bool usemap = true; + + // map file to use (get that from La Palma!) + const char *map = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL; + + // ------------------------------------------------------ + + Long_t max0 = 0; + Long_t max1 = 0; + + // ====================================================== + + if (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 1; + } + if (gSystem->AccessPathName(drsfile, kFileExists)) + { + gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl; + return 2; + } + + // -------------------------------------------------------------------------------- + + gLog.Separator("Pulses"); + gLog << all << "Extract pulses with '" << drsfile << "'" << endl; + gLog << endl; + + // ------------------------------------------------------ + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile)) + return 3; + + // ------------------------------------------------------ + + iter.Sort(); + iter.Print(); + + TString title; + title = iter.Next(); + title += "; "; + title += max1; + iter.Reset(); + + // ====================================================== + + MStatusDisplay *d = new MStatusDisplay; + + MBadPixelsCam badpixels; + badpixels.InitSize(64); + badpixels[7].SetUnsuitable(MBadPixelsPix::kUnsuitableRun); + badpixels[17].SetUnsuitable(MBadPixelsPix::kUnsuitableRun); + + MH::SetPalette(); + + MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal"); + + MGeomApply apply; + + MDrsCalibApply drsapply; + drsapply.SetMaxNumPrevEvents(0); // Switched off step-correction -> crashes due to the requirement of N*9 ch + + MRawFitsRead read; + read.LoadMap(map); + read.AddFiles(iter); + + // ====================================================== + + gLog << endl; + gLog.Separator("Calculating baseline"); + + MTaskList tlist0; + + MRawRunHeader header; + MPedestalCam pedcam; + + MParList plist; + plist.AddToList(&tlist0); + plist.AddToList(&drscalib300); + plist.AddToList(&header); + plist.AddToList(&pedcam); + + // ------------------ Setup the tasks --------------- + + MFillH fill0("MHBaseline", "MPedestalSubtractedEvt"); + + drsapply.SetRemoveSpikes(4); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist0.AddToList(&read); + tlist0.AddToList(&apply); + tlist0.AddToList(&drsapply); + tlist0.AddToList(&fill0); + + // ------------------ Setup and run eventloop --------------- + + MEvtLoop loop0(title); + loop0.SetDisplay(d); + loop0.SetParList(&plist); + + if (!loop0.Eventloop(max0)) + return 4; + + if (!loop0.GetDisplay()) + return 5; + + // ---------------------------------------------------------------- + + MGeomCamFAMOUS fact; + MHCamera hped(fact); + hped.SetCamContent(pedcam); + hped.SetCamError(pedcam, 4); + hped.SetAllUsed(); + MHCamEvent display; + display.SetHist(hped); + + d->AddTab("Baseline"); + display.Clone()->Draw(); + + // ====================================================== + + gLog << endl; + gLog.Separator("Extracting singles"); + + MTaskList tlist1; + + MSingles singles; + + plist.Replace(&tlist1); + plist.AddToList(&badpixels); + plist.AddToList(&singles); + + // ------------------ Setup the tasks --------------- + + MExtractSingles extract; + extract.SetMaxDist(max_dist); + extract.SetThreshold(threshold); + + MFillH fill("MHSingles"); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist1.AddToList(&read); + tlist1.AddToList(&apply); + tlist1.AddToList(&drsapply); + tlist1.AddToList(&extract); + tlist1.AddToList(&fill); + + // ------------------ Setup and run eventloop --------------- + + MEvtLoop loop1(title); + loop1.SetDisplay(d); + loop1.SetParList(&plist); + + if (!loop1.Eventloop(max1)) + return 6; + + if (!loop1.GetDisplay()) + return 7; + + if (fill.GetNumExecutions()==0) + return 8; + + // ============================================================= + + MHSingles* h = (MHSingles*) plist.FindObject("MHSingles"); + if (!h) + return 9; + + gStyle->SetOptStat(10); + + TH1 *h1 = h->GetSignal1(); + TH1 *p1 = h->GetProf1(); + + TH1 *h2 = h->GetSignal2(); + TH1 *p2 = h->GetProf2(); + + TH1 *h3 = h->GetSignal3(); + TH1 *p3 = h->GetProf3(); + + TH1 *h4 = h->GetSignal4(); + TH1 *p4 = h->GetProf4(); + + TH1 *h5 = h->GetSignal5(); + TH1 *p5 = h->GetProf5(); + + TH1 *h6 = h->GetSignal6(); + TH1 *p6 = h->GetProf6(); + + TH1 *h7 = h->GetSignal7(); + TH1 *p7 = h->GetProf7(); + + TH1 *h8 = h->GetSignal8(); + 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); + + // ===== + + h1->SetMinimum(h1->GetEntries()/1000); + h2->SetMinimum(h2->GetEntries()/1000); + h3->SetMinimum(h3->GetEntries()/1000); + h4->SetMinimum(h4->GetEntries()/1000); + h5->SetMinimum(h5->GetEntries()/1000); + h6->SetMinimum(h6->GetEntries()/1000); + h7->SetMinimum(h7->GetEntries()/1000); + h8->SetMinimum(h8->GetEntries()/1000); + + // ===== + + d->AddTab("1"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + 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"); + + // ===== + + d->AddTab("2"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + gPad->SetRightMargin(0.15); + h2->DrawCopy("colz"); + p2->DrawCopy("same"); + + pulse.SetParameter(0, 2*1.3); + p2->Fit(&pulse, "N0", "", 1.5, 25); + TF1 *f2 = pulse.DrawCopy("same"); + + // ===== + + d->AddTab("3"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + gPad->SetRightMargin(0.15); + h3->DrawCopy("colz"); + p3->DrawCopy("same"); + + pulse.SetParameter(0, 3*1.3); + p3->Fit(&pulse, "N0", "", 1.5, 25); + TF1 *f3 = pulse.DrawCopy("same"); + + // ===== + + d->AddTab("4"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + gPad->SetRightMargin(0.15); + h4->DrawCopy("colz"); + p4->DrawCopy("same"); + + pulse.SetParameter(0, 4*1.3); + p4->Fit(&pulse, "N0", "", 1.5, 25); + TF1 *f4 = pulse.DrawCopy("same"); + + // ===== + + d->AddTab("5"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + gPad->SetRightMargin(0.15); + h5->DrawCopy("colz"); + p5->DrawCopy("same"); + + pulse.SetParameter(0, 5*1.3); + p5->Fit(&pulse, "N0", "", 1.5, 25); + TF1 *f5 = pulse.DrawCopy("same"); + + // ===== + + d->AddTab("6"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + gPad->SetRightMargin(0.15); + h6->DrawCopy("colz"); + p6->DrawCopy("same"); + + pulse.SetParameter(0, 6*1.3); + p6->Fit(&pulse, "N0", "", 1.5, 25); + TF1 *f6 = pulse.DrawCopy("same"); + + // ===== + + d->AddTab("7"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + gPad->SetRightMargin(0.15); + h7->DrawCopy("colz"); + p7->DrawCopy("same"); + + pulse.SetParameter(0, 7*1.3); + p7->Fit(&pulse, "N0", "", 1.5, 25); + TF1 *f7 = pulse.DrawCopy("same"); + + // ===== + + d->AddTab("8"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + gPad->SetRightMargin(0.15); + h8->DrawCopy("colz"); + p8->DrawCopy("same"); + + pulse.SetParameter(0, 6*1.3); + p8->Fit(&pulse, "N0", "", 1.5, 25); + TF1 *f8 = pulse.DrawCopy("same"); + + // ===== + + d->AddTab("Pulse"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.07); + gPad->SetTopMargin(0.01); + gPad->SetRightMargin(0.15); + p8->SetStats(kFALSE); + p8->GetXaxis()->SetRangeUser(-25, 300); + p8->GetXaxis()->CenterTitle(); + p8->GetXaxis()->SetTitleSize(0.045); + p8->GetYaxis()->SetTitleSize(0.045); + p8->GetYaxis()->SetTitleOffset(0.8); + p8->SetYTitle("Amplitude"); + p8->SetXTitle("Time [ns]"); + p8->DrawCopy(); + p7->DrawCopy("same"); + p6->DrawCopy("same"); + p5->DrawCopy("same"); + p4->DrawCopy("same"); + p3->DrawCopy("same"); + p2->DrawCopy("same"); + p1->DrawCopy("same"); + + f1->Draw("same"); + f2->Draw("same"); + f3->Draw("same"); + f4->Draw("same"); + f5->Draw("same"); + f6->Draw("same"); + f7->Draw("same"); + f8->Draw("same"); + + d->SaveAs("paper_pulse.root"); + + + /* + d->SaveAs(outfile); + + TFile f(outfile, "UPDATE"); + + MParameterI par("NumEvents"); + par.SetVal(fill.GetNumExecutions()); + par.Write(); + + MParameterI win("IntegrationWindow"); + win.SetVal(extract.GetIntegrationSize()); + win.Write(); + + extract.GetExtractionRange().Write("ExtractionRange"); + + if (firstRunID==lastRunID) + header.Write(); + */ + return 0; +}