Index: fact/tools/marsmacros/singlepe.C =================================================================== --- fact/tools/marsmacros/singlepe.C (revision 14165) +++ fact/tools/marsmacros/singlepe.C (revision 14166) @@ -3,4 +3,6 @@ #include #include +#include +#include #include #include @@ -28,4 +30,5 @@ #include "MHCamera.h" #include "MGeomCamFACT.h" +#include "MRawRunHeader.h" using namespace std; @@ -36,14 +39,16 @@ struct Single { - float fSignal; - float fTime; + float fSignal; + float fTime; }; class MSingles : public MParContainer, public MCamEvent { + Int_t fIntegrationWindow; + vector> fData; public: - MSingles(const char *name=NULL, const char *title=NULL) + MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0) { fName = name ? name : "MSingles"; @@ -60,6 +65,10 @@ 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 &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const { + return kTRUE; } void DrawPixelContent(Int_t num) const { } @@ -68,21 +77,91 @@ }; +class MH2F : public TH2F +{ +public: + Int_t Fill(Double_t x,Double_t y) + { + Int_t binx, biny, bin; + fEntries++; + binx = TMath::Nint(x+1); + biny = fYaxis.FindBin(y); + bin = biny*(fXaxis.GetNbins()+2) + binx; + AddBinContent(bin); + if (fSumw2.fN) ++fSumw2.fArray[bin]; + if (biny == 0 || biny > fYaxis.GetNbins()) { + if (!fgStatOverflows) return -1; + } + ++fTsumw; + ++fTsumw2; + fTsumwy += y; + fTsumwy2 += y*y; + return bin; + } + ClassDef(MH2F, 1); +}; + +class MProfile2D : public TProfile2D +{ +public: + Int_t Fill(Double_t x, Double_t y, Double_t z) + { + Int_t bin,binx,biny; + fEntries++; + binx =TMath::Nint(x+1); + biny =fYaxis.FindBin(y); + bin = GetBin(binx, biny); + fArray[bin] += z; + fSumw2.fArray[bin] += z*z; + fBinEntries.fArray[bin] += 1; + if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1; + if (biny == 0 || biny > fYaxis.GetNbins()) { + if (!fgStatOverflows) return -1; + } + ++fTsumw; + ++fTsumw2; + fTsumwy += y; + fTsumwy2 += y*y; + fTsumwz += z; + fTsumwz2 += z*z; + return bin; + } + ClassDef(MProfile2D, 1); +}; + + class MHSingles : public MH { - MSingles *fSingles; - - TH2F fHist; + MH2F fSignal; + MH2F fTime; + MProfile2D fPulse; + + UInt_t fNumEntries; + + MSingles *fSingles; //! + MPedestalSubtractedEvt *fEvent; //! public: - MHSingles() : fSingles(0) + MHSingles() : fNumEntries(0), fSingles(0), fEvent(0) { fName = "MHSingles"; - fHist.SetName("Name"); - fHist.SetTitle("Title"); - fHist.SetXTitle("X title"); - fHist.SetYTitle("Y title"); - fHist.SetDirectory(NULL); + fSignal.SetName("Signal"); + fSignal.SetTitle("Title"); + fSignal.SetXTitle("X title"); + fSignal.SetYTitle("Y title"); + fSignal.SetDirectory(NULL); + + fTime.SetName("Time"); + fTime.SetTitle("Title"); + fTime.SetXTitle("X title"); + fTime.SetYTitle("Y title"); + fTime.SetDirectory(NULL); + + fPulse.SetName("Pulse"); + fPulse.SetTitle("Title"); + fPulse.SetXTitle("X title"); + fPulse.SetYTitle("Y title"); + fPulse.SetDirectory(NULL); } @@ -96,13 +175,40 @@ } + fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); + if (!fEvent) + { + *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl; + return kFALSE; + } + + fNumEntries = 0; + return kTRUE; } - Bool_t ReInit(MParList *) - { + Bool_t ReInit(MParList *plist) + { + const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader"); + if (!header) + { + *fLog << /* err << */ "MRawRunHeader not found... abort." << endl; + return kFALSE; + } + + const Int_t w = fSingles->GetIntegrationWindow(); + MBinning binsx, binsy; binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5); - binsy.SetEdges(300, -25, 600-25); - - MH::SetBinning(fHist, binsx, binsy); + binsy.SetEdges(300, -7.5*w, (60-7.5)*w); + + MH::SetBinning(fSignal, binsx, binsy); + + const UShort_t roi = header->GetNumSamples(); + + // Correct for DRS timing!! + MBinning binst(roi, -0.5, roi-0.5); + MH::SetBinning(fTime, binsx, binst); + + MBinning binspy(2*roi, -roi-0.5, roi-0.5); + MH::SetBinning(fPulse, binsx, binspy); return kTRUE; @@ -111,4 +217,7 @@ Int_t Fill(const MParContainer *par, const Stat_t weight=1) { + const Float_t *ptr = fEvent->GetSamples(0); + const Short_t roi = fEvent->GetNumSamples(); + for (unsigned int i=0; iGetNumPixels(); i++) { @@ -116,11 +225,30 @@ for (unsigned int j=0; jDivide(1,1); + pad->Divide(2,2); pad->cd(1); - gPad->SetBorderMode(0); gPad->SetFrameBorderMode(0); - - fHist.Draw("colz"); + fSignal.Draw("colz"); + + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fTime.Draw("colz"); + + pad->cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fPulse.Draw("colz"); } @@ -143,5 +279,12 @@ }; +MStatusDisplay *d = new MStatusDisplay; + MSingles *fSingles; + +TH1F *fluct1 = new TH1F("", "", 100, -50, 50); +TH1F *fluct2 = new TH1F("", "", 100, -50, 50); + +TGraph weights("weights.txt"); Int_t PreProcess(MParList *plist) @@ -158,4 +301,13 @@ } + d->AddTab("Fluct"); + fluct2->Draw(); + fluct1->Draw("same"); + + fSingles->SetIntegrationWindow(20); + + //if (weights.GetN()==0) + // return kFALSE; + return kTRUE; } @@ -163,64 +315,54 @@ Int_t Process() { - - for (int i=0; iGetNumPixels(); i++) - { - vector &result = fSingles->GetVector(i); + const UInt_t roi = fEvent->GetNumSamples(); + + const size_t navg = 10; + const float threshold = 5; + const UInt_t start_skip = 20; + const UInt_t end_skip = 10; + const UInt_t integration_size = 10*3; + + vector val(roi-navg);//result of first sliding average + for (int pix=0; pixGetNumPixels(); pix++) + { + vector &result = fSingles->GetVector(pix); result.clear(); - const UInt_t roi = fEvent->GetNumSamples(); - - const Float_t *ptr = fEvent->GetSamples(i); - - //--------------------------- - // ...Extract pulses here.... - //--------------------------- - - //sliding average parameter (center weight = 3) - int avg1 = 8; - vector Vslide(roi-2*avg1);//result of first sliding average + const Float_t *ptr = fEvent->GetSamples(pix); //Sliding window - for (UInt_t k=avg1; kthreshold) - &&(Vslide[k+5]>threshold) - ) - { - //average baseline before crossing - //double baseline = (Vslide[k-start_sample]+Vslide[k-start_sample+1]+Vslide[k-start_sample-1])/3.; - single.fSignal = 0; - single.fTime = k+integration_delay; - for(UInt_t l=integration_delay; lthreshold || + val[i+4]>threshold || + + val[i+5]SetPalette(1,0); + // ------------------------------------------------------ @@ -295,5 +433,5 @@ // ====================================================== - MStatusDisplay *d = new MStatusDisplay; + //MStatusDisplay *d = new MStatusDisplay; MBadPixelsCam badpixels; @@ -316,5 +454,5 @@ //MDrsCalibrationTime timecam; - gStyle->SetOptFit(kTRUE); + MH::SetPalette(); // ====================================================== @@ -326,4 +464,5 @@ MSingles singles; + MRawRunHeader header; MParList plist1; @@ -332,4 +471,5 @@ plist1.AddToList(&badpixels); plist1.AddToList(&singles); + plist1.AddToList(&header); MEvtLoop loop1("DetermineCalConst"); @@ -371,9 +511,24 @@ return 11; + gStyle->SetOptFit(kTRUE); + MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles"); if (!h) return 0; - TH2 *hist = h->GetHist(); + TH2 *hsignal = h->GetSignal(); + TH2 *htime = h->GetTime(); + TH2 *hpulse = h->GetPulse(); + + d->AddTab("Time"); + TH1 *ht = htime->ProjectionY(); + ht->Scale(1./1440); + ht->Draw(); + + d->AddTab("Pulse"); + TH1 *hp = hpulse->ProjectionY(); + hp->Scale(1./1440); + hp->Draw(); + // ------------------ Spectrum Fit --------------- @@ -381,220 +536,138 @@ // access the spektrumhistogram by the pointer *hist - hist->SetTitle("Amplitude Spectrum vs. Pixel ID"); - hist->SetXTitle("Pixel SoftID"); - hist->SetYTitle("Amplitude [mV]"); - hist->SetZTitle("counts"); - MGeomCamFACT fact; MHCamera hcam(fact); - struct fit - { - float fAmplitude; - float fGain; - float f2GainRMS; - float fCrosstlk; - float fOffset; - }; - fit fitParameters[singles.GetNumPixels()]; - - TCanvas &canv1 = d->AddTab("Fits"); - TCanvas &canv2 = d->AddTab("Distributions"); - canv2.Divide(3,2); - TCanvas &gain_canv = d->AddTab("Gain Distribution"); - //built an array of Amplitude spectra - TH1D* hAmplSpek[singles.GetNumPixels()]; - TString histotitle; - - TF1 gaus4("g4","gaus"); - gaus4.SetLineColor(2); - - TH1F hAmplitude( - "hAmplitude", - "Amplitude of Single Preak", - 200, - 50, - 250 - ); - hAmplitude.GetXaxis()->SetTitle("Amplitude (from Integral 5 Slices) [mV]"); - hAmplitude.GetYaxis()->SetTitle("counts"); - - TH1F hGain( - "hGain", - "Gain", - 40*5, - 30, - 70 - ); - hGain.GetXaxis()->SetTitle("Gain (from Integral 5 Slices) [mV]"); - hGain.GetYaxis()->SetTitle("counts"); - - TH1F hGainRMS( - "hSigma", - "Sigma of gauss peak in gainfit", - 20*5, - 0, - 40 - ); - hGainRMS.GetXaxis()->SetTitle("Sigma (from Integral 5 Slices) [mV]"); - hGainRMS.GetYaxis()->SetTitle("counts"); - - TH1F hCrosstlk( - "hCrosstlk", - "Parameter proportional to Crosstalk", - 30*5, - 0, - 4 - ); - hCrosstlk.GetXaxis()->SetTitle("~Crosstalk (from Integral 5 Slices) [mV]"); - hCrosstlk.GetYaxis()->SetTitle("counts"); - - TH1F hOffset( - "hOffset", - "X-Offset of Spectrum", - 100*2, - 0, - 4 - ); - hOffset.GetXaxis()->SetTitle(" X-Offset of Spectrum (from Integral 5 Slices) [mV]"); - hOffset.GetYaxis()->SetTitle("counts"); - - TH1F hFitProb( - "hFitProb", - "Fit Probability", - 100*2, - 0, - 100 - ); - hFitProb.GetXaxis()->SetTitle("fit-probability [%]"); - hFitProb.GetYaxis()->SetTitle("counts"); + TH1F hAmplitude("Rate", "Average number of dark counts per event", + 100, 0, 10); + TH1F hGain ("Gain", "Gain distribution", + 50, 100, 250); + TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma", + 100, 0, 30); + TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability", + 100, 0, 25); + TH1F hOffset ("Offset", "Baseline offset distribution", + 50, -7.5, 2.5); + TH1F hFitProb ("FitProb", "FitProb distribution", + 50, 0, 100); + TH1F hSum1 ("Sum1", "Sum", + 100, 0, 2000); + TH1F hSum2 ("Sum2", "Sum", + 100, 0, 10); // define fit function for Amplitudespectrum - TF1 fspek_fit("fspek_fit", spek_fit, 10, 600, 5); - fspek_fit.SetParNames("Ampl 1.Max", "Gain", "Sigma", "Ampl 2.Max", "Offset"); - - // fit parameters - Double_t par[5]; + TF1 func("spektrum", fcn, 0, 1000, 5); + func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset"); + func.SetLineColor(kBlue); // Map for which pixel shall be plotted and which not TArrayC usePixel(1440); + memset(usePixel.GetArray(), 1, 1440); // List of Pixel that should be ignored in camera view - Int_t doNotUse[12]= { - 424, - 923, - 1208, - 583, - 830, - 1399, - 113, - 115, - 354, - 423, - 1195, - 1393 - }; + usePixel[424] = 0; + usePixel[923] = 0; + usePixel[1208] = 0; + usePixel[583] = 0; + usePixel[830] = 0; + usePixel[1399] = 0; + usePixel[113] = 0; + usePixel[115] = 0; + usePixel[354] = 0; + usePixel[423] = 0; + usePixel[1195] = 0; + usePixel[1393] = 0; + + int count_ok = 0; // Begin of Loop over Pixels for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++) { - usePixel[pixel] = 0; - - //define projection histogram title - histotitle = "hPixelPro_"; - histotitle += pixel; + if (usePixel[pixel]==0) + continue; //Projectipon of a certain Pixel out of Ampl.Spectrum - hAmplSpek[pixel] = hist->ProjectionY( histotitle , pixel+1, pixel+1); - - //set histogram and Axis title - histotitle = "Amplitude Spectrum of Pixel #"; - histotitle += pixel; - hAmplSpek[pixel]->SetTitle(histotitle); - - // define appearance of fit - fspek_fit.SetLineColor(kBlue); - fspek_fit.SetLineStyle(2); - - par[0] = hAmplSpek[pixel]->GetBinContent(hAmplSpek[pixel]->GetMaximumBin()); + TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); + hist->SetDirectory(0); + + const Int_t maxbin = hist->GetMaximumBin(); + const Double_t binwidth = hist->GetBinWidth(maxbin); + const Double_t ampl = hist->GetBinContent(maxbin); + + double fwhm = 0; + for (int i=1; iGetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)GetBinCenter(hAmplSpek[pixel]->GetMaximumBin()); - par[2] = 0.2*par[1]; - par[3] = 0.1*par[0]; - par[4] = 0.2*hAmplSpek[pixel]->GetBinCenter(hAmplSpek[pixel]->GetMaximumBin()); - fspek_fit.SetParameters(par); + Double_t par[5] = + { + ampl, hist->GetBinCenter(maxbin), fwhm*1.4, 0.1, -10 + }; + + func.SetParameters(par); + + const double min = par[1]-fwhm*1.4; + const double max = par[1]*5; //Fit Pixels spectrum - canv1.cd(); - TFitResultPtr fitStatus(hAmplSpek[ pixel ]->Fit(&fspek_fit,"+S","",0.8*par[1],5*par[1])); - double fit_prob = fitStatus->Prob(); - - //cout << "pixel: " << pixel << "\t status: " << fitStatus << endl; - - fitParameters[pixel].fAmplitude = fspek_fit.GetParameter(0); - fitParameters[pixel].fGain = fspek_fit.GetParameter(1); - fitParameters[pixel].f2GainRMS = fspek_fit.GetParameter(2); - fitParameters[pixel].fCrosstlk = fspek_fit.GetParameter(3)/fspek_fit.GetParameter(0); - fitParameters[pixel].fOffset = fspek_fit.GetParameter(4); - - gain_canv.cd(); - gPad->SetLogy(1); - hGain.Fill(fitParameters[pixel].fGain); - hGain.Fit(&gaus4); - hGain.DrawCopy(); - - canv2.cd(2); - gPad->SetLogy(1); - hGain.DrawCopy(); - - canv2.cd(3); - gPad->SetLogy(1); - hGainRMS.Fill(fitParameters[pixel].f2GainRMS); - hGainRMS.DrawCopy(); - - canv2.cd(1); - gPad->SetLogy(1); - hAmplitude.Fill(fitParameters[pixel].fAmplitude); - hAmplitude.DrawCopy(); - - canv2.cd(4); - gPad->SetLogy(1); - hCrosstlk.Fill(fitParameters[pixel].fCrosstlk); - hCrosstlk.DrawCopy(); - - canv2.cd(5); - gPad->SetLogy(1); - hOffset.Fill(fitParameters[pixel].fOffset); - hOffset.DrawCopy(); - - canv2.cd(6); - gPad->SetLogy(1); + const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", min, max); + + const bool ok = int(rc)==0; + + const double fit_prob = rc->Prob(); + const float fGain = func.GetParameter(1); + const float fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0); + const float f2GainRMS = func.GetParameter(2); + const float fCrosstlk = func.GetParameter(3); + const float fOffset = func.GetParameter(4)/singles.GetIntegrationWindow(); + + hAmplitude.Fill(fAmplitude); + hGain.Fill(fGain); + //hGainRMS.Fill(f2GainRMS); + hCrosstlk.Fill(fCrosstlk*100); + hOffset.Fill(fOffset); hFitProb.Fill(100*fit_prob); - hFitProb.DrawCopy(); - - - hcam.SetBinContent(pixel+1, fitParameters[pixel].fGain); - if (!fitStatus) - { - usePixel[pixel] = 1; - cout << "Pixel: " << pixel << " success!" << endl << endl; - } - else - { - cout << "Pixel: " << pixel << " failed!" << endl << endl; - } - - } -// TCanvas &canvBroken[12]; -// TString name; - for (int i = 0; i < 12; i++) - { -// name = "Pixel"; -// name += doNotUse[i]; -// canvBroken[i] = d->AddTab(name); -// hAmplSpek[doNotUse[ i ]]->Draw(); - - usePixel[doNotUse[ i ]] = 0; + hGainRMS.Fill(100*f2GainRMS/fGain); + + hcam.SetBinContent(pixel+1, fGain); + + usePixel[pixel] = ok; + + if (!ok) + { + cout << "Pixel " << setw(4) << pixel << ": failed!" << endl; + continue; + } + + for (int b=1; b<=hist->GetNbinsX(); b++) + { + hSum1.Fill(hist->GetBinCenter(b), hist->GetBinContent(b)); + hSum2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b)); + } + + if (count_ok==0) + { + TCanvas &c = d->AddTab(Form("Pix%d", pixel)); + c.cd(); + hist->SetName(Form("Pix%d", pixel)); + hist->DrawCopy()->SetDirectory(0); + //func.DrawCopy("same"); + } + + count_ok++; + + delete hist; } @@ -606,61 +679,70 @@ hcam.DrawCopy(); - - canv1.cd(); + TCanvas &c11 = d->AddTab("SumHist"); + gPad->SetLogy(); + hSum1.DrawCopy(); + + TCanvas &c12 = d->AddTab("GainHist"); + gPad->SetLogy(); + hSum2.DrawCopy(); + + TCanvas &c2 = d->AddTab("Result"); + c2.Divide(3,2); + c2.cd(1); + gPad->SetLogy(); + hAmplitude.DrawCopy(); + + c2.cd(2); + gPad->SetLogy(); + hGain.DrawCopy(); + + c2.cd(3); + gPad->SetLogy(); + hGainRMS.DrawCopy(); + + c2.cd(4); + gPad->SetLogy(); + hCrosstlk.DrawCopy(); + + c2.cd(5); + gPad->SetLogy(); + hOffset.DrawCopy(); + + c2.cd(6); + gPad->SetLogy(); + hFitProb.DrawCopy(); /* - TString title = "-- Calibrated signal #"; - title += seq.GetSequence(); - title += " ("; - title += drsfile; - title += ") --"; - d->SetTitle(title, kFALSE); - - TString path; - path += Form("%s/20%6d_%03d-calibration.root", outpath, - seq.GetSequence()/1000, seq.GetSequence()%1000); - - d->SaveAs(path); - */ - + TCanvas &c3 = d->AddTab("Separation"); + gPad->SetLogy(); + hSep.DrawCopy(); + */ return 0; } -Double_t -spek_fit( - Double_t* x, - Double_t* par - ) -{ - Double_t arg = 0; - Double_t cross = 1; - Double_t peak = 0; -// Double_t arg2 = 0; -// Double_t arg3 = 0; -// Double_t arg4 = 0; - -// if (par[0] != 0) cross = par[3]/par[0]; -// if (par[2] != 0) arg = (x[0] - par[1] + par[4])/par[2]; -// if (par[2] != 0) arg2 = (x[0] - 2*par[1] + par[4])/par[2]; -// if (par[2] != 0) arg3 = (x[0] - 3*par[1] + par[4])/par[2]; -// if (par[2] != 0) arg4 = (x[0] - 4*par[1] + par[4])/par[2]; -// if (par[0] = 0) return 0; -// if (par[2] = 0) return 0; - - if (par[0] != 0) cross = par[3]/par[0]; +Double_t fcn(Double_t *xx, Double_t *par) +{ + Double_t x = xx[0]; + + Double_t ampl = par[0]; + Double_t gain = par[1]; + Double_t sigma = par[2]; + Double_t cross = par[3]; + Double_t shift = par[4]; + + Double_t xTalk = 1; + + Double_t y = 0; for (int order = 1; order < 5; order++) { - arg = (x[0] - order*par[1] - par[4] )/par[2]; - Double_t xTalk = TMath::Power(cross, order-1); -// Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); - Double_t gauss = TMath::Exp(-0.5 * arg * arg); - peak += xTalk*par[0]*gauss; - } - -// peak += par[0]*TMath::Exp(-0.5 * arg * arg); -// peak += cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2); -// peak += cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3); -// peak += cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4); - - return peak; + Double_t arg = (x - order*gain - shift)/sigma; + + Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); + + y += xTalk*gauss; + + xTalk *= cross; + } + + return ampl*y; }