Index: trunk/MagicSoft/Mars/Makefile =================================================================== --- trunk/MagicSoft/Mars/Makefile (revision 8153) +++ trunk/MagicSoft/Mars/Makefile (revision 8154) @@ -45,4 +45,5 @@ mhist \ manalysis \ + mextralgo \ msignal \ mbadpixels \ Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc =================================================================== --- trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc (revision 8153) +++ trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc (revision 8154) @@ -38,6 +38,102 @@ } +// ----------------------------------------------------------------------------- +// +// Calculates the chi2 of the fit, once the weights have been iterated. +// Argument: time, obtained after a call to EvalDigitalFilterHiGain +// +Float_t MExtralgoDigitalFilter::GetChisq(const Int_t maxp, const Int_t frac, const Float_t sum) const +{ + /* + TMatrix g (windowh,1); + TMatrix gt(windowh,1); + TMatrix y (windowh,1); + + const Float_t etau = fFineAdjustHi*sumhi; + const Float_t t_fine = TMath::Abs(fFineAdjustHi)< 1./fBinningResolutionHiGain ? -fFineAdjustHi : 0.; + + // if (t_fine==0.) + // return -1.; + + if (fDebug) + gLog << all << " fMaxPHi: " << fMaxPHi << " fIterPHi " << fIterPHi << " t_fine: " << t_fine << endl; + + // + // Slide with a window of size windowh over the sample + // and calculate the arrays by interpolating the pulse shape using the + // fine-tuned time information. + // + for (Int_t sample=0; sample < windowh; sample++) + { + const Int_t idx = fArrBinningResHalfHiGain[sample] + fIterPHi; + const Int_t ids = fMaxPHi + sample; + + y [sample][0] = fHiGainSignalDF[ids]; + g [sample][0] = t_fine >= 0 + ? (fPulseShapeHiGain[idx] + t_fine*(fPulseShapeHiGain[idx+1] -fPulseShapeHiGain[idx]) )*sumhi + : (fPulseShapeHiGain[idx] + t_fine*(fPulseShapeHiGain[idx] -fPulseShapeHiGain[idx-1]))*sumhi; + gt[sample][0] = t_fine >= 0 + ? (fPulseShapeDotHiGain[idx] + t_fine*(fPulseShapeDotHiGain[idx+1]-fPulseShapeDotHiGain[idx]) )*etau + : (fPulseShapeDotHiGain[idx] + t_fine*(fPulseShapeDotHiGain[idx] -fPulseShapeDotHiGain[idx-1]) )*etau; + } + + TMatrix second = y - g - gt; + TMatrix first(TMatrix::kTransposed,second); + TMatrix chisq = first * ((*fBHiInv)*second); + + return chisq[0][0]/(windowh-2); + */ +/* + + TMatrix S(fWindowSize, 1); // Signal (start==start of window) + for (int i=0; i -void MExtralgoDigitalFilter::Extract() +void MExtralgoDigitalFilter::Extract(Int_t maxpos) { fSignal = 0; // default is: no pulse found @@ -54,6 +150,5 @@ // Calculate the sum of the first fWindowSize slices // - - // For the case of an even numberof weights/bin there is + // For the case of an even number of weights/bin there is // no central bin.So we create an artificial central bin. for (Int_t i=0; i @@ -145,5 +248,5 @@ #include -Bool_t MExtralgoDigitalFilter::CalculateWeights(TH1 &shape, const TH2 &autocorr, TArrayF &weightsamp, TArrayF &weightstime, Int_t wpb) +Int_t MExtralgoDigitalFilter::CalculateWeights(TH1 &shape, const TH2 &autocorr, TArrayF &weightsamp, TArrayF &weightstime, Int_t wpb) { const Int_t weightsperbin = wpb<=0?shape.GetNbinsX()/autocorr.GetNbinsX():wpb; @@ -154,5 +257,5 @@ cout << " Shape: " << shape.GetNbinsX() << endl; cout << " ACorr: " << autocorr.GetNbinsX() << endl; - return kFALSE; + return -1; } @@ -180,4 +283,5 @@ // FIXME: DELETE!!! TH1 &derivative = *static_cast(shape.Clone()); + derivative.SetDirectory(0); derivative.Reset(); @@ -233,5 +337,5 @@ { cout << "ERROR - Division by zero: denom[0][0]==0 for i=" << i << "." << endl; - return kFALSE; + return -1; } @@ -248,8 +352,8 @@ } - return kTRUE; + return first*weightsperbin; } -void MExtralgoDigitalFilter::CalculateWeights2(TH1F &shape, const TH2F &autocorr, TArrayF &weightsamp, TArrayF &weightstime, Int_t wpb) +Int_t MExtralgoDigitalFilter::CalculateWeights2(TH1 &shape, const TH2 &autocorr, TArrayF &weightsamp, TArrayF &weightstime, Int_t wpb) { const Int_t weightsperbin = wpb<=0?shape.GetNbinsX()/autocorr.GetNbinsX():wpb; @@ -260,5 +364,5 @@ cout << " Shape: " << shape.GetNbinsX() << endl; cout << " ACorr: " << autocorr.GetNbinsX() << endl; - return; + return -1; } @@ -282,5 +386,7 @@ TSpline5 val("Signal", &gr); - TH1F derivative(shape); + // FIXME: DELETE!!! + TH1 &derivative = *static_cast(shape.Clone()); + derivative.SetDirectory(0); derivative.Reset(); @@ -331,4 +437,10 @@ const TMatrixD denom = (gT*(B*g))*(dT*(B*d)) - (dT*(B*g))*(dT*(B*g)); + if (denom[0][0]==0) + { + cout << "ERROR - Division by zero: denom[0][0]==0 for i=" << i << "." << endl; + return -1; + } + const TMatrixD w_amp = (dT*(B*d))*(gT*B) - (gT*(B*d))*(dT*B); const TMatrixD w_time = (gT*(B*g))*(dT*B) - (gT*(B*d))*(gT*B); @@ -342,46 +454,4 @@ } } + return first*weightsperbin; } -/* -Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1); -for (Int_t i = -fBinningResolutionHiGain/2+1; i<=fBinningResolutionHiGain/2; i++) -{ - // i = -4...5 - for (Int_t count=0; count < fWindowSizeHiGain; count++) - { - // count=0: pos=10*(start+0) + 10 + i - // count=1: pos=10*(start+1) + 10 + i - // count=2: pos=10*(start+2) + 10 + i - // count=3: pos=10*(start+3) + 10 + i - } - - for (Int_t count=0; count < fWindowSizeHiGain; count++) - { - // count=0: pos = 10*0 + 5 +i-1 - // count=1: pos = 10*1 + 5 +i-1 - // count=2: pos = 10*2 + 5 +i-1 - // count=3: pos = 10*3 + 5 +i-1 - } -} - -Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionLoGain/2; -for (Int_t i = -fBinningResolutionLoGain/2+1; i<=fBinningResolutionLoGain/2; i++) -{ - // i=-4..5 - for (Int_t count=0; count < fWindowSizeLoGain; count++) - { - // count=0: pos = 10*(start+0) + 5 + i - // count=1: pos = 10*(start+1) + 5 + i - // count=2: pos = 10*(start+2) + 5 + i - // count=3: pos = 10*(start+3) + 5 + i - } - - for (Int_t count=0; count < fWindowSizeLoGain; count++) - { - // count=0: pos = 10*0 + 5 +i-1 - // count=1: pos = 10*1 + 5 +i-1 - // count=2: pos = 10*2 + 5 +i-1 - // count=3: pos = 10*3 + 5 +i-1 - } -} -*/ Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h =================================================================== --- trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h (revision 8153) +++ trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h (revision 8154) @@ -2,6 +2,6 @@ #define MARS_MExtralgoDigitalFilter -#ifndef ROOT_TROOT -#include +#ifndef ROOT_TMatrix +#include #endif @@ -12,13 +12,18 @@ class TArrayF; +//#include + class MExtralgoDigitalFilter { private: // Input - Float_t *fVal; - Int_t fNum; + const Float_t *fVal; + Int_t fNum; Float_t const *fWeightsAmp; Float_t const *fWeightsTime; + Float_t const *fPulseShape; + + const TMatrix *fAinv; const Int_t fWeightsPerBin; // Number of weights per data bin @@ -30,4 +35,27 @@ Float_t fSignal; Float_t fSignalDev; + + Float_t GetChisq(const Int_t maxp, const Int_t frac, const Float_t sum) const; + + inline Double_t ChiSq(const Double_t sum, const Int_t startv, const Int_t startw=0) const + { + // + // Slide with a window of size windowsize over the sample + // and multiply the entries with the corresponding weights + // + Double_t chisq = 0; + + // Shift the start of the weight to the center of sample 0 + Float_t const *w = fPulseShape + startw; + + const Float_t *beg = fVal+startv; + for (Float_t const *pex=beg; pex Float_t MExtralgoSpline::ExtractNoise(Int_t iter) { - const Float_t nsx = iter * fResolution; + const Float_t nsx = gRandom->Uniform(); //iter * fResolution; if (fExtractionType == kAmplitude) @@ -223,5 +224,5 @@ fSignalDev = -1; fTimeDev = -1; - +/* // // Allow no saturated slice and @@ -229,5 +230,4 @@ // - /* const Bool_t limlo = maxbin < TMath::Ceil(fRiseTime); const Bool_t limup = maxbin > fNum-TMath::Ceil(fFallTime)-1; @@ -248,7 +248,4 @@ return; } - */ - - fTimeDev = fResolution; // @@ -256,6 +253,4 @@ // - - /* Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one @@ -346,6 +341,12 @@ } } - */ - + + fTime = maxpos; + fTimeDev = fResolution; + fSignal = CalcIntegral(maxpos); + fSignalDev = 0; // means: is valid + + return; +*/ // --- Start NEW --- @@ -369,4 +370,5 @@ Float_t maxpos, maxval; + // FIXME: Check the dfeault if no maximum found!!! GetMaxAroundI(maxbin, maxpos, maxval); @@ -376,4 +378,5 @@ { fTime = maxpos; + fTimeDev = fResolution; fSignal = maxval; fSignalDev = 0; // means: is valid @@ -381,30 +384,4 @@ } - // - // Loop from the beginning of the slice upwards to reach the maxhalf: - // With means of bisection: - // - /* - static const Float_t sqrt2 = TMath::Sqrt(2.); - - step = sqrt2*3*0.061;//fRiseTime; - Float_t x = maxpos-0.86-3*0.061;//fRiseTime*1.25; - -// step = sqrt2*0.5;//fRiseTime; -// Float_t x = maxpos-1.25;//fRiseTime*1.25; - - Int_t cnt =0; - while (cnt++<30) - { - const Float_t y=EvalAt(x); - - if (TMath::Abs(y-maxval/2)maxval/2 ? -step : +step; - } - */ - // Search downwards for maxval/2 // By doing also a search upwards we could extract the pulse width @@ -412,9 +389,33 @@ fTime = x1; + fTimeDev = fResolution; fSignal = CalcIntegral(maxpos); fSignalDev = 0; // means: is valid - //if (fSignal>100) - // cout << "V=" << maxpos-x1 << endl; + // + // Loop from the beginning of the slice upwards to reach the maxhalf: + // With means of bisection: + // + /* + static const Float_t sqrt2 = TMath::Sqrt(2.); + + step = sqrt2*3*0.061;//fRiseTime; + Float_t x = maxpos-0.86-3*0.061;//fRiseTime*1.25; + +// step = sqrt2*0.5;//fRiseTime; +// Float_t x = maxpos-1.25;//fRiseTime*1.25; + + Int_t cnt =0; + while (cnt++<30) + { + const Float_t y=EvalAt(x); + + if (TMath::Abs(y-maxval/2)maxval/2 ? -step : +step; + } + */ // @@ -422,4 +423,5 @@ // // fTime = cnt==31 ? -1 : x; + // fTimeDev = fResolution; // fSignal = cnt==31 ? CalcIntegral(x) : CalcIntegral(maxpos); // fSignalDev = 0; // means: is valid Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h =================================================================== --- trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h (revision 8153) +++ trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h (revision 8154) @@ -21,6 +21,6 @@ // Input - Float_t *fVal; - Int_t fNum; + Float_t const *fVal; + const Int_t fNum; Float_t *fDer1; @@ -238,4 +238,5 @@ // Take a default in case no maximum is found + // FIXME: Check THIS!!! xmax=i; ymax=fVal[i]; @@ -438,5 +439,5 @@ public: - MExtralgoSpline(Float_t *val, Int_t n, Float_t *der1, Float_t *der2) + MExtralgoSpline(const Float_t *val, Int_t n, Float_t *der1, Float_t *der2) : fExtractionType(kIntegral), fVal(val), fNum(n), fDer1(der1), fDer2(der2), fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1) { Index: trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc (revision 8153) +++ trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc (revision 8154) @@ -94,4 +94,5 @@ #include "MFTriggerPattern.h" #include "MGeomApply.h" +#include "MPedestalSubtract.h" //#include "MMcPedestalCopy.h" #include "MPointingPosCalc.h" @@ -450,4 +451,8 @@ MContinue conttp(&ftp, "ContTrigPattern"); + // Create the pedestal subtracted raw-data + MPedestalSubtract pedsub; + pedsub.SetPedestalCam(&pedcamab); + // Do signal and pedestal calculation MPedCalcFromLoGain pedlo1("MPedCalcFundamental"); @@ -471,15 +476,30 @@ extractor1->SetPedestals(&pedcamab); + // Setup to use the hi-gain extraction window in the lo-gain + // range (the start of the lo-gain range is added automatically + // by MPedCalcFromLoGain) + // + // The window size of the extractor is not yet initialized, + // so we have to stick to the extraction range + const Int_t f = extractor1->GetHiGainFirst(); + const Int_t l = extractor1->GetHiGainLast(); + const Int_t w = (l-f+1)&~1; + + pedlo1.SetExtractWindow(f, w); + if (extractor1->InheritsFrom("MExtractTimeAndCharge")) { pedlo2.SetExtractor((MExtractTimeAndCharge*)extractor1); pedlo3.SetExtractor((MExtractTimeAndCharge*)extractor1); + /* const Int_t win = ((MExtractTimeAndCharge*)extractor1)->GetWindowSizeHiGain(); pedlo1.SetExtractWindow(15, win); - pedlo2.SetExtractWindow(15, win/*obsolete*/); - pedlo3.SetExtractWindow(15, win/*obsolete*/); + pedlo2.SetExtractWindow(15, win//obsolete//); + pedlo3.SetExtractWindow(15, win//obsolete//); + */ } else { + /* // FIXME: How to get the fixed value 15 automatically? const Int_t f = (Int_t)(15.5+extractor1->GetHiGainFirst()); @@ -488,4 +508,8 @@ pedlo2.SetExtractWindow(f, n); pedlo3.SetExtractWindow(f, n); + */ + pedlo2.SetExtractWindow(f, w); + pedlo3.SetExtractWindow(f, w); + } } @@ -657,4 +681,5 @@ tlist2.AddToList(&apply); tlist2.AddToList(&merge); + tlist2.AddToList(&pedsub); tlist2.AddToList(&pedlo1); tlist2.AddToList(&pedlo2); Index: trunk/MagicSoft/Mars/mjobs/MJCalibration.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCalibration.cc (revision 8153) +++ trunk/MagicSoft/Mars/mjobs/MJCalibration.cc (revision 8154) @@ -154,4 +154,5 @@ #include "MRawFileRead.h" #include "MGeomApply.h" +#include "MPedestalSubtract.h" #include "MTaskEnv.h" #include "MBadPixelsMerge.h" @@ -1907,6 +1908,21 @@ tlist.AddToList(&apply); - MPedCalcPedRun pedcalc; - pedcalc.SetExtractWindow(fExtractor->GetHiGainFirst(),TMath::Nint(fExtractor->GetNumHiGainSamples())); + // Produce pedestal subtracted raw-data + MPedestalSubtract pedsub; + pedsub.SetPedestalCam(fExtractor->GetPedestals()); + tlist.AddToList(&pedsub); + + // Setup to use the hi-gain extraction window in the lo-gain + // range (the start of the lo-gain range is added automatically + // by MPedCalcFromLoGain) + // + // The window size of the extractor is not yet initialized, + // so we have to stick to the extraction range + const Int_t f = fExtractor->GetHiGainFirst(); + const Int_t l = fExtractor->GetHiGainLast(); + const Int_t w = (l-f+1)&~1; + + MPedCalcPedRun pedcalc; + pedcalc.SetExtractWindow(f, w); if (IsIntensity()) Index: trunk/MagicSoft/Mars/mjobs/MJPedestal.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJPedestal.cc (revision 8153) +++ trunk/MagicSoft/Mars/mjobs/MJPedestal.cc (revision 8154) @@ -82,4 +82,5 @@ #include "MRawEvtData.h" #include "MGeomApply.h" +#include "MPedestalSubtract.h" #include "MTriggerPatternDecode.h" #include "MBadPixelsMerge.h" @@ -1064,9 +1065,9 @@ // if (!fBadPixelsFile.IsNull()) - { + { *fLog << inf << "Excluding: " << fBadPixelsFile << endl; - ifstream fin(fBadPixelsFile.Data()); - fBadPixels.AsciiRead((istream&)fin); - } + ifstream fin(fBadPixelsFile); + fBadPixels.AsciiRead(fin); + } MGeomApply geomapl; @@ -1074,5 +1075,5 @@ MPedCalcPedRun pedcalc; - pedcalc.SetPedestalUpdate(kFALSE); + //pedcalc.SetPedestalUpdate(kFALSE); MPedCalcFromLoGain pedlogain; @@ -1132,4 +1133,12 @@ tlist.AddToList(&fillpul); } + + // produce pedestal subtracted raw-data + MPedestalSubtract pedsub; + if (fExtractor && fExtractionType!=kFundamental) + pedsub.SetPedestalCam(&fPedestalCamIn); + else + pedsub.SetNamePedestalCam(""); // Only copy hi- and lo-gain together! + tlist.AddToList(&pedsub); // ---------------------------------------------------------------------- @@ -1186,38 +1195,35 @@ if (fExtractor) { - fExtractor->SetPedestals(&fPedestalCamIn); - - if (fExtractionType!=kFundamental) + fExtractor->SetPedestals(&fPedestalCamIn); + + if (fExtractionType!=kFundamental) { - pedcalc.SetRandomCalculation(fExtractionType==kWithExtractorRndm); - pedlogain.SetRandomCalculation(fExtractionType==kWithExtractorRndm); - - pedcalc.SetExtractor((MExtractTimeAndCharge*)fExtractor); - pedlogain.SetExtractor((MExtractTimeAndCharge*)fExtractor); + pedcalc.SetRandomCalculation(fExtractionType==kWithExtractorRndm); + pedlogain.SetRandomCalculation(fExtractionType==kWithExtractorRndm); + + pedcalc.SetExtractor((MExtractTimeAndCharge*)fExtractor); + pedlogain.SetExtractor((MExtractTimeAndCharge*)fExtractor); } - - if (fExtractor->InheritsFrom("MExtractTimeAndCharge")) + else { - - const Float_t f = 0.1+fExtractor->GetHiGainFirst(); - const Int_t win = ((MExtractTimeAndCharge*)fExtractor)->GetWindowSizeHiGain(); - pedcalc.SetExtractWindow((Int_t)f, win); - pedlogain.SetExtractWindow((Int_t)(15+f), win); - + // The window size of the extractor is not yet initialized, + // so we have to stick to the extraction range + const Int_t f = fExtractor->GetHiGainFirst(); + const Int_t l = fExtractor->GetHiGainLast(); + const Int_t w = (l-f+1)&~1; + + // Setup to use the hi-gain extraction window in the lo-gain + // range (the start of the lo-gain range is added automatically + // by MPedCalcFromLoGain) + pedcalc.SetExtractWindow(f, w); + pedlogain.SetExtractWindow(f, w); } - else + + if (!fExtractor->InheritsFrom("MExtractTimeAndCharge") && fExtractionType!=kFundamental) { - const Float_t f = 0.1+fExtractor->GetHiGainFirst(); - const Float_t n = 0.1+fExtractor->GetNumHiGainSamples(); - pedcalc.SetExtractWindow((Int_t)f, (Int_t)n); - pedlogain.SetExtractWindow((Int_t)(15+f), (Int_t)n); - - if (fExtractionType!=kFundamental) - { - *fLog << inf; - *fLog << "Signal extractor doesn't inherit from MExtractTimeAndCharge..." << endl; - *fLog << " --> falling back to fundamental pedestal extraction." << endl; - fExtractionType=kFundamental; - } + *fLog << inf; + *fLog << "Signal extractor doesn't inherit from MExtractTimeAndCharge..." << endl; + *fLog << " --> falling back to fundamental pedestal extraction." << endl; + fExtractionType=kFundamental; } } Index: trunk/MagicSoft/Mars/mpedestal/MPedestalSubtract.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedestalSubtract.cc (revision 8153) +++ trunk/MagicSoft/Mars/mpedestal/MPedestalSubtract.cc (revision 8154) @@ -1,3 +1,5 @@ /* ======================================================================== *\ +! $Name: not supported by cvs2svn $:$Id: MPedestalSubtract.cc,v 1.2 2006-10-24 08:18:07 tbretz Exp $ +! -------------------------------------------------------------------------- ! ! * Index: trunk/MagicSoft/Mars/mpedestal/MPedestalSubtract.h =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedestalSubtract.h (revision 8154) +++ trunk/MagicSoft/Mars/mpedestal/MPedestalSubtract.h (revision 8154) @@ -0,0 +1,39 @@ +#ifndef MARS_MPedestalSubtract +#define MARS_MPedestalSubtract + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MRawEvtData; +class MPedestalCam; +class MPedestalSubtractedEvt; + +class MPedestalSubtract : public MTask +{ +private: + static const TString fgNamePedestalCam; //! "MPedestalCam" + static const TString fgNamePedestalSubtractedEvt; //! "MPedestalSubtractedEvt" + + MRawEvtData *fRawEvt; //! Input Raw data + MPedestalCam *fPedestals; //! Pedestals of all pixels in the camera + MPedestalSubtractedEvt *fSignal; //! Output container + + TString fNamePedestalCam; // Name of the 'MPedestalCam' container + TString fNamePedestalSubtractedEvt; // Name of the 'MPedestalSubtractedEvt' container + + Int_t PreProcess(MParList *pList); + Int_t Process(); + +public: + MPedestalSubtract(const char *name=NULL, const char *title=NULL); + + void SetPedestalCam(MPedestalCam *pedcam) { fPedestals = pedcam; } + void SetNamePedestalCam(const char *name=fgNamePedestalCam.Data()) { fNamePedestalCam = name; } + + void SetNamePedestalSubtractedEvt(const char *name=fgNamePedestalSubtractedEvt.Data()) { fNamePedestalSubtractedEvt = name; } + + ClassDef(MPedestalSubtract, 0) // Class to subtract pedestal including AB-noise +}; + +#endif Index: trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.cc (revision 8154) +++ trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.cc (revision 8154) @@ -0,0 +1,170 @@ +/* ======================================================================== *\ +! $Name: not supported by cvs2svn $:$Id: MPedestalSubtractedEvt.cc,v 1.1 2006-10-24 08:22:41 tbretz Exp $ +! -------------------------------------------------------------------------- +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 10/2006 +! +! Copyright: MAGIC Software Development, 2000-2006 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MPedestalSubtractedEvt +// +// Storage container to store the raw FADC values. +// +///////////////////////////////////////////////////////////////////////////// + +#include "MPedestalSubtractedEvt.h" + +ClassImp(MPedestalSubtractedEvt); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Initialize number of samples (per pixel) and number of pixels. +// +// Initialize the correct length of fSamples and fSamplesRaw +// +// And reset its contents to 0. +// +void MPedestalSubtractedEvt::InitSamples(UInt_t pixels, UInt_t samples) +{ + fNumSamples = samples; + fNumPixels = pixels; + + fSamples.Set(pixels*samples); + fSamplesRaw.Set(pixels*samples); + + fSamples.Reset(); + fSamplesRaw.Reset(); +} + +// -------------------------------------------------------------------------- +// +// Return a pointer to the first slice with subtracted pedestal of +// the samples of the pixel with given pixel-index. If the number +// exceeds the number of pixels NULL is returned. +// +// The user is responsible not to exceed the slices for one pixel! +// +Float_t *MPedestalSubtractedEvt::GetSamples(UInt_t pixel) const +{ + return pixel>=fNumPixels ? NULL : fSamples.GetArray()+pixel*fNumSamples; +} + +// -------------------------------------------------------------------------- +// +// Return a pointer to the first slice of the raw-data samples of the pixel +// with given pixel-index. If the number exceeds the number of +// pixels NULL is returned. +// +// The user is responsible not to exceed the slices for one pixel! +// +Byte_t *MPedestalSubtractedEvt::GetSamplesRaw(UInt_t pixel) const +{ + return pixel>=fNumPixels ? NULL : fSamplesRaw.GetArray()+pixel*fNumSamples; +} + +// -------------------------------------------------------------------------- +// +// Return some information about saturation in the raw-data of pixel idx. +// +// The search range is defined by [first,last]. Saturation is considered if +// contents is >= limit. +// +// The number of saturating slices are returned and first/last are filled +// with the first and last saturating slice index w.r.t. the beginning of +// the raw-data of this pixel not first. +// +// Warning: No range checks and no sanity checks are done! +// +Int_t MPedestalSubtractedEvt::GetSaturation(const Int_t idx, Int_t limit, Int_t &first, Int_t &last) const +{ + // Determin saturation of hi-gains + Byte_t *p0 = GetSamplesRaw(idx); + + Byte_t *sathi0 = 0; // first saturating hi-gain slice + Byte_t *sathi1 = 0; // last saturating hi-gain slice + + Int_t num = 0; + + const Byte_t *end = p0+last; + for (Byte_t *ptr=p0+first; ptr<=end; ptr++) + { + if (*ptr>=limit) + { + sathi1 = ptr; + if (!sathi0) + sathi0 = ptr; + num++; + } + } + + last = sathi1 ? sathi1-p0 : -1; + first = sathi0 ? sathi0-p0 : -1; + + return num; +} + +// -------------------------------------------------------------------------- +// +// Get the maximum of the slices [first,last] of pixel index. +// +// The position returned is the index of the position of the pedestal +// subtracted maximum w.r.t. to first. +// The value returned is the maximum of the raw-data. +// +// Warning: No range checks and no sanity checks are done! +// +Int_t MPedestalSubtractedEvt::GetMax(const Int_t idx, const Int_t first, const Int_t last, Int_t &val) const +{ + // Get pointer to first slice to be considered + Byte_t const *samb = GetSamplesRaw(idx); + Float_t const *samf = GetSamples(idx); + + Byte_t const *ptrb = samb+first; + Float_t const *ptrf = samf+first; + + // The best information so far: the first slice is the maximum + const Byte_t *maxb = ptrb; + const Float_t *maxf = ptrf; + + // Store the pointer to the first slice +// const Byte_t *begb = ptrb; + const Float_t *begf = ptrf; + + // Calculate the last slice to be considered + const Byte_t *endb = samb+last; + + while (ptrb*maxb) + maxb = ptrb; + if (*++ptrf>*maxf) + maxf = ptrf; + } + + val = *maxb; + return maxf-begf; +} + Index: trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.h =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.h (revision 8154) +++ trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.h (revision 8154) @@ -0,0 +1,64 @@ +#ifndef MARS_MPedestalSubtractedEvt +#define MARS_MPedestalSubtractedEvt + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +#ifndef MARS_MCamEvent +#include "MCamEvent.h" +#endif +#ifndef MARS_MArrayF +#include "MArrayF.h" +#endif +#ifndef MARS_MArrayB +#include "MArrayB.h" +#endif + +class MPedestalSubtractedEvt : public MParContainer, public MCamEvent +{ +private: + MArrayF fSamples; // list of all samples with pedestal subtracted + MArrayB fSamplesRaw; // list of all samples (raw) + + UInt_t fNumSamples; // number of samples per pixel + UInt_t fNumPixels; // number of pixels + +public: + MPedestalSubtractedEvt(const char *name=NULL, const char *title=NULL) + : fNumSamples(0), fNumPixels(0) + { + } + + void InitSamples(UInt_t pixels, UInt_t samples); + + Float_t *GetSamples(UInt_t pixel) const; + Byte_t *GetSamplesRaw(UInt_t pixel) const; + + UInt_t GetNumSamples() const { return fNumSamples; } + UShort_t GetNumPixels() const { return fNumPixels; } + + Int_t GetSaturation(const Int_t idx, Int_t limit, Int_t &first, Int_t &last) const; + Int_t GetMax(const Int_t pixidx, const Int_t first, const Int_t last, Int_t &val) const; + Int_t GetMax(const Int_t pixidx, const Int_t first, const Int_t last) const + { + Int_t val; + return GetMax(pixidx, first, last, val); + } + Int_t GetMax(const Int_t pixidx, Int_t &val) const + { + return GetMax(pixidx, 0, fNumSamples, val); + } + + Int_t GetMax(const Int_t pixidx) const + { + Int_t val; + return GetMax(pixidx, 0, fNumSamples, val); + } + + 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 { } + + ClassDef(MPedestalSubtractedEvt, 6) //Container to store the raw Event Data +}; + +#endif Index: trunk/MagicSoft/Mars/mpedestal/Makefile =================================================================== --- trunk/MagicSoft/Mars/mpedestal/Makefile (revision 8153) +++ trunk/MagicSoft/Mars/mpedestal/Makefile (revision 8154) @@ -42,6 +42,5 @@ MHPedestalCor.cc \ MPedestalSubtract.cc \ - MPedestalSubtractedEvt.cc \ - MPedCalcFromData.cc + MPedestalSubtractedEvt.cc ############################################################ Index: trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h (revision 8153) +++ trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h (revision 8154) @@ -19,5 +19,4 @@ #pragma link C++ class MPedestalSubtract+; #pragma link C++ class MPedestalSubtractedEvt+; -#pragma link C++ class MPedCalcFromData+; #pragma link C++ class MHPedestalCor+; Index: trunk/MagicSoft/Mars/msignal/MExtractTime.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTime.cc (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractTime.cc (revision 8154) @@ -1,3 +1,5 @@ /* ======================================================================== *\ +! $Name: not supported by cvs2svn $:$Id: MExtractTime.cc,v 1.22 2006-10-24 08:24:52 tbretz Exp $ +! -------------------------------------------------------------------------- ! ! * @@ -17,6 +19,7 @@ ! ! Author(s): Markus Gaug, 04/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 +! Author(s): Thomas Bretz, 04/2004 +! +! Copyright: MAGIC Software Development, 2000-2006 ! ! @@ -137,5 +140,5 @@ // - MArrivalTimeCam::SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples, // fLoGainFirst, fLoGainLast, fNumLoGainSamples); -// +/* Bool_t MExtractTime::ReInit(MParList *pList) { @@ -148,10 +151,10 @@ return kTRUE; } - +*/ // -------------------------------------------------------------------------- // // Calculate the integral of the FADC time slices and store them as a new // pixel in the MArrivalTimeCam container. -// +/* Int_t MExtractTime::Process() { @@ -183,16 +186,16 @@ pix.SetGainSaturation(sathi, satlo); - } /* while (pixel.Next()) */ + } fArrTime->SetReadyToSave(); return kTRUE; -} +}*/ void MExtractTime::Print(Option_t *o) const { - if (IsA()==MExtractTime::Class()) - *fLog << GetDescriptor() << ":" << endl; +// if (IsA()==MExtractTime::Class()) +// *fLog << GetDescriptor() << ":" << endl; + MExtractor::Print(o); *fLog << " Offset Lo-Gain: " << fOffsetLoGain << endl; - MExtractor::Print(o); -} +} Index: trunk/MagicSoft/Mars/msignal/MExtractTime.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTime.h (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractTime.h (revision 8154) @@ -17,13 +17,13 @@ MArrivalTimeCam *fArrTime; //! Container with the photons arrival times - +/* virtual void FindTimeHiGain(Byte_t *firstused, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const {} virtual void FindTimeLoGain(Byte_t *firstused, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const {} - + */ Int_t PreProcess( MParList *pList ); - Bool_t ReInit ( MParList *pList ); - Int_t Process (); +// Bool_t ReInit ( MParList *pList ); +// Int_t Process (); public: Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc (revision 8154) @@ -1,3 +1,5 @@ /* ======================================================================== *\ +! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.55 2006-10-24 08:24:52 tbretz Exp $ +! -------------------------------------------------------------------------- ! ! * @@ -17,6 +19,7 @@ ! ! Author(s): Markus Gaug, 05/2004 +! Author(s): Thomas Bretz, 05/2004 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2006 ! ! @@ -60,4 +63,6 @@ // MExtractedSignalCam // +// For weired events check: Run 94127 Event 672, 1028 +// ////////////////////////////////////////////////////////////////////////////// #include "MExtractTimeAndCharge.h" @@ -70,10 +75,14 @@ #include "MParList.h" -#include "MRawEvtData.h" +//#include "MArrayB.h" +//#include "MArrayF.h" + +//#include "MRawEvtData.h" #include "MRawEvtPixelIter.h" -#include "MRawRunHeader.h" - -#include "MPedestalCam.h" -#include "MPedestalPix.h" +//#include "MRawRunHeader.h" + +//#include "MPedestalCam.h" +//#include "MPedestalPix.h" +#include "MPedestalSubtractedEvt.h" #include "MArrivalTimeCam.h" @@ -87,6 +96,6 @@ using namespace std; -const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -3.5; -const Byte_t MExtractTimeAndCharge::fgLoGainSwitch = 120; +const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -6.0; +const Byte_t MExtractTimeAndCharge::fgLoGainSwitch = 120; // -------------------------------------------------------------------------- @@ -95,5 +104,4 @@ // // Sets: -// - fLoGainFirstSave to 0 // - fWindowSizeHiGain and fWindowSizeLoGain to 0 // - fLoGainStartShift to fgLoGainStartShift+fgOffsetLoGain @@ -101,5 +109,5 @@ // MExtractTimeAndCharge::MExtractTimeAndCharge(const char *name, const char *title) - : fLoGainFirstSave(0), fWindowSizeHiGain(0), fWindowSizeLoGain(0) + : /*fLoGainFirstSave(0),*/ fWindowSizeHiGain(0), fWindowSizeLoGain(0) { fName = name ? name : "MExtractTimeAndCharge"; @@ -157,7 +165,6 @@ if (fSignals) - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples, + fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples, fLoGainFirst, fLoGainLast, fNumLoGainSamples); - return kTRUE; @@ -175,148 +182,322 @@ Int_t MExtractTimeAndCharge::Process() { - - MRawEvtPixelIter pixel(fRawEvt); - - while (pixel.Next()) - { - // COPY HERE PRODUCING ARRAY WITH SAMPLES - - /* - const MPedestalPix &ped = (*fPedestals)[pixidx]; - - const Float_t pedes = ped.GetPedestal(); - const Float_t ABoffs = ped.GetPedestalABoffset(); - - const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; - - const Int_t num = pixel.GetNumHiGainSamples()+pixel.GetNumLoGainSamples(); - - MArrayF sample(num); - - const Int_t abflag = pixel.HasABFlag() ? 1 : 0; - const Int_t ids0 = fHiGainFirst + abflag; - - Int_t ids = ids0; - - Float_t null = 0; // Starting value for maxpos - Float_t *maxpos = &null; // Position of maximum - Float_t *sat1 = 0; // First saturating slice - Float_t *sat2 = 0; // Last saturating slice - - const Float_t *beg = sample.GetArray(); - const Float_t *end = beg + num; - - Float_t *ptr = beg; - while (ptr *maxpos) - // maxpos = ptr; - // - // if (*ptr >= fSaturationLimit) - // { - // sat2 = ptr; - // if (!sat1) - // sat1 = ptr; - // } - // - // ptr++; + MRawEvtPixelIter pixel(fRawEvt); + + while (pixel.Next()) + { + const Int_t pixidx = pixel.GetPixelId(); + + // + // Preparations for the pedestal subtraction (with AB-noise correction) + // + Float_t *sig = fSignal->GetSamples(pixidx); + + // ==================================================================== + // MOVE THIS TO MExtractTimeAndCharge + // ==================================================================== + // number of hi- and lo-gains + const Int_t numh = pixel.GetNumHiGainSamples(); + const Int_t numl = pixel.GetNumLoGainSamples(); +/* + // COPY HERE PRODUCING ARRAY WITH SAMPLES + static MArrayB sample; + sample.Set(numh+numl); + + if (numl>0) + { + memcpy(sample.GetArray(), pixel.GetHiGainSamples(), numh); + memcpy(sample.GetArray()+numh, pixel.GetLoGainSamples(), numl); } - */ - - // - // Find signal in hi- and lo-gain - // - Float_t sumhi =0., deltasumhi =0; // Set hi-gain of MExtractedSignalPix valid - Float_t timehi=0., deltatimehi=0; // Set hi-gain of MArrivalTimePix valid - Byte_t sathi=0; - - // Initialize fMaxBinContent for the case, it gets not set by the derived class - fMaxBinContent = fLoGainSwitch + 1; - - const Int_t pixidx = pixel.GetPixelId(); - const MPedestalPix &ped = (*fPedestals)[pixidx]; - const Bool_t higainabflag = pixel.HasABFlag(); - - FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), - sumhi, deltasumhi, timehi, deltatimehi, - sathi, ped, higainabflag); - - // Make sure that in cases the time couldn't be correctly determined - // more meaningfull default values are assigned - if (timehi<=0 || timehi>pixel.GetNumHiGainSamples()) - timehi = gRandom->Uniform(pixel.GetNumHiGainSamples()); - - Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix - Float_t timelo=0., deltatimelo=-1; // invalidate logain of MArrivalTimePix - Byte_t satlo=0; - - // - // Adapt the low-gain extraction range from the obtained high-gain time - // - - // IN THIS CASE THE PIXEL SHOULD BE MARKED BAD!!!! - // MEANS: Hi gain has saturated, but it is too early to extract - // the lo-gain properly - // THIS produces pulse positions ~= -1 - // The signal might be handled in MCalibrateData, but hwat's about - // the arrival times in MCalibrateRelTime - if (sathi && fMaxBinContent<=fLoGainSwitch) - deltasumlo=deltasumhi=deltatimelo=deltatimehi=-1; - - // FIXME: What to do with the pixel if it saturates too early??? - if (pixel.HasLoGain() && (fMaxBinContent > fLoGainSwitch /*|| sathi>0*/) ) - { - fLoGainFirstSave = fLoGainFirst; - - // sathi is the number (not index!) of the first saturating slice - // 0 indicates that no saturation was found - // FIMXME: Is 3 an accurate value? - - const Float_t pos = sathi==0 ? timehi : (int)(sathi)-3; - - if (pos+fLoGainStartShift>0) - fLoGainFirst = (Byte_t)(pos + fLoGainStartShift); - - if (fLoGainFirst= fWindowSizeLoGain) - { - deltasumlo = 0; // make logain of MExtractedSignalPix valid - deltatimelo = 0; // make logain of MArrivalTimePix valid - - const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1; - FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst, - sumlo, deltasumlo, timelo, deltatimelo, - satlo, ped, logainabflag); - - if (satlo>6) - deltasumlo=deltatimelo=-1; - } - fLoGainFirst = fLoGainFirstSave; - - // Make sure that in cases the time couldn't be correctly determined - // more meaningfull default values are assigned - if (timelo<=0 || timelo>pixel.GetNumLoGainSamples()) - timelo = gRandom->Uniform(pixel.GetNumLoGainSamples()); - } - - MExtractedSignalPix &pix = (*fSignals)[pixidx]; - MArrivalTimePix &tix = (*fArrTime)[pixidx]; - pix.SetExtractedSignal(sumhi, deltasumhi, sumlo, deltasumlo); - pix.SetGainSaturation(sathi, satlo); - - tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo); - tix.SetGainSaturation(sathi, satlo); - + + // get pedestal information + const MPedestalPix &pedpix = (*fPedestals)[pixidx]; + + // pedestal information + const Int_t ab = pixel.HasABFlag() ? 1 : 0; + const Float_t ped = pedpix.GetPedestal(); + + // destination array + static MArrayF sig; + sig.Set(numh+numl); + + // start of destination array, end of hi-gain destination array + // and start of hi-gain samples + Float_t *beg = sig.GetArray(); + Float_t *end = beg + sig.GetSize(); + + // determine with which pedestal (+/- AB offset) to start + const Bool_t noswap = (ab&1)==((beg-beg)&1); + const Float_t offh = noswap ? pedpix.GetPedestalABoffset() : -pedpix.GetPedestalABoffset(); + const Float_t mean[2] = { ped + offh, ped - offh }; + + //const Int_t rangeh = fHiGainLast - fHiGainFirst + 1 + fHiLoLast; + //const Int_t rangel = fLoGainLast - fLoGainFirst + 1; + + const Byte_t *src = numl>0 ? sample.GetArray() : pixel.GetHiGainSamples(); + //const Byte_t *src = sample.GetArray(); + + // Copy hi-gains into array and substract pedestal + // FIXME: Shell we really subtract the pedestal from saturating slices??? + for (Float_t *ptr=beg; ptrGetSamplesRaw(pixidx); //numl>0 ? sample.GetArray() : pixel.GetHiGainSamples(); +// Byte_t *p0 = sample.GetArray(); + + Byte_t maxcont = 0; + Int_t maxposhi = -1; + + Byte_t *sathi0 = 0; // first saturating hi-gain slice + Byte_t *sathi1 = 0; // last saturating hi-gain slice + for (Byte_t *ptr=p0+fHiGainFirst; ptrmaxcont) + { + // ORGINAL: + //Int_t range = (fHiGainLast - fHiGainFirst + 1 + fHiLoLast) - fWindowSizeHiGain + 1; + // maxpos>1 && maxpos<=range + + // This is a strange workaround to fit the previous + // Spline setup: TO BE CHECKED! + const Int_t pos = ptr-p0; + if (pos1 && maxposhi=fSaturationLimit) + { + sathi1 = ptr; + if (!sathi0) + sathi0 = ptr; + } + } +*/ + Int_t sathi0 = fHiGainFirst; // First slice to extract and first saturating slice + Int_t sathi1 = fHiGainLast/*+fHiLoLast*/; // Last slice to extract and last saturating slice + + Int_t maxcont; + Int_t maxposhi = fSignal->GetMax(pixidx, sathi0, sathi1, maxcont); + // Would it be better to take lastsat-firstsat? + Int_t numsathi = fSignal->GetSaturation(pixidx, fSaturationLimit, sathi0, sathi1); +/* + // sathi2 is the number (not the index) of first saturating slice +// const Byte_t sathi2 = sathi0<0 ? 0 : sathi0+1; + + // Number (not index) of first saturating slice + // const Byte_t sathi2 = sathi0==0 ? 0 : sathi0-p0+1; +// const Byte_t numsathi = sathi0==0 ? 0 : sathi1-sathi0+1; + +// Int_t maxb = maxcont; + + // ==================================================================== + // FIXME: Move range out of the extractors... + + // Initialize maxcont for the case, it gets not set by the derived class + Float_t sumhi2 =0., deltasumhi2 =-1; // Set hi-gain of MExtractedSignalPix valid + Float_t timehi2=0., deltatimehi2=-1; // Set hi-gain of MArrivalTimePix valid +*/ + Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid + Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid + if (numsathi<1) + { + + const Int_t rangehi = fHiGainLast - fHiGainFirst + 1/* + fHiLoLast*/; + FindTimeAndChargeHiGain2(sig/*.GetArray()*/+fHiGainFirst, rangehi, + sumhi, deltasumhi, timehi, deltatimehi, + numsathi, maxposhi); + + timehi += fHiGainFirst; + } +/* + Float_t sumhi = sumhi2; + Float_t deltasumhi = deltasumhi2; + Float_t timehi = timehi2; + Float_t deltatimehi=deltatimehi2; +// Byte_t sathi = sathi2; + + + // + // Find signal in hi- and lo-gain + // + Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid + Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid + Byte_t sathi=0; + + // Initialize maxcont for the case, it gets not set by the derived class + maxcont = fLoGainSwitch + 1; + + //const Int_t pixidx = pixel.GetPixelId(); + //const MPedestalPix &ped = (*fPedestals)[pixidx]; + const Bool_t higainabflag = pixel.HasABFlag(); + FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), + sumhi, deltasumhi, timehi, deltatimehi, + sathi, pedpix, higainabflag); + timehi += fHiGainFirst; + + */ + // Make sure that in cases the time couldn't be correctly determined + // more meaningfull default values are assigned + if (timehi<=fHiGainFirst || timehi>=fHiGainLast) + timehi = gRandom->Uniform(fHiGainLast-fHiGainFirst)+fHiGainFirst; + + Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix + Float_t timelo=0., deltatimelo=-1; // invalidate logain of MArrivalTimePix + Byte_t satlo=0; + Int_t numsatlo=0; + + // + // Adapt the low-gain extraction range from the obtained high-gain time + // + + // IN THIS CASE THE PIXEL SHOULD BE MARKED BAD!!!! + // MEANS: Hi gain has saturated, but the signal is to dim + // to extract the lo-gain properly + // THIS produces pulse positions ~= -1 + // The signal might be handled in MCalibrateData, but hwat's about + // the arrival times in MCalibrateRelTime + if (numsathi>0 && maxcont<=fLoGainSwitch) + deltasumlo=deltasumhi=deltatimelo=deltatimehi=-1; + + // If more than 8 hi-gain slices have saturated this is + // no physical event. We just assume that something with + // the extraction is wrong + if (numsathi>8) // FIXME: Should be something like 2? + deltasumhi=deltatimehi=-1; + + // FIXME: What to do with the pixel if it saturates too early??? + if (pixel.HasLoGain() && (maxcont > fLoGainSwitch /*|| sathi>0*/) ) + { + // To determin the window in which the lo-gain is extracted + // clearly more information about the relation between the + // extraction window and the reslting time is necessary. + //fLoGainStartShift = -6.0; + + const Byte_t fLoGainFirstSave = fLoGainFirst; + + // sathi is the number (not index!) of the first saturating slice + // 0 indicates that no saturation was found + // FIXME: Is 0.5 should be expressed by the rise time of + // the hi-gain signal! + const Float_t pos = numsathi==0 ? timehi : sathi0-0.5; + + const Int_t lostart = TMath::FloorNint(pos+6); + + const Int_t newfirst = TMath::FloorNint(pos+fLoGainStartShift); + fLoGainFirst = newfirst>fLoGainFirstSave ? newfirst : fLoGainFirstSave; + +// if (fLoGainLast-fLoGainFirst >= fWindowSizeLoGain) + { +/* + Float_t sumlo2 =0., deltasumlo2 =-1.; // invalidate logain of MExtractedSignalPix + Float_t timelo2=0., deltatimelo2=-1.; // invalidate logain of MArrivalTimePix + + // Determin saturation in lo-gains + Byte_t *satlo0 = 0; // first saturating hi-gain slice + Byte_t *satlo1 = 0; // last saturating hi-gain slice + Int_t maxposlo = -1; + Byte_t maxlo = 0; + + Byte_t *p0 = fSignal->GetSamplesRaw(pixidx); + for (Byte_t *ptr=p0+numh+fLoGainFirst; ptrmaxlo) + { + // This is a strange workaround to fit the previous + // Spline setup: TO BE CHECKED! + const Int_t ipos = ptr-p0-numh; + if (ipos>=fLoGainFirst && ipos<=fLoGainLast) + { + maxposlo = ipos-fLoGainFirst; + maxlo = *ptr; + } + } + if (*ptr>=fSaturationLimit) + { + satlo1 = ptr; + if (!satlo0) + satlo0 = ptr; + } + // FIXME: Do we also have to really COUNT the number + // of saturating slices, eg. for crosschecks? + } + + // Number of saturating lo-gain slices (this is not necessary + // identical to a counter counting the number of saturating + // slices!) + const Byte_t numsatlo = satlo0==0 ? 0 : satlo1-satlo0+1; +*/ + + Int_t satlo0 = numh+fLoGainFirst; // First slice to extract and first saturating slice + Int_t satlo1 = numh+fLoGainLast; // Last slice to extract and last saturating slice + + // Would it be better to take lastsat-firstsat? + Int_t maxlo; + Int_t maxposlo = fSignal->GetMax(pixidx, satlo0, satlo1, maxlo); + numsatlo = fSignal->GetSaturation(pixidx, fSaturationLimit, satlo0, satlo1); + + // const Byte_t satlo2 = numsatlo; + + const Int_t rangelo = fLoGainLast - fLoGainFirst + 1; + FindTimeAndChargeLoGain2(sig/*.GetArray()*/+numh+fLoGainFirst, rangelo, + sumlo, deltasumlo, timelo, deltatimelo, + numsatlo, maxposlo); + timelo += fLoGainFirst; + +// sumlo =sumlo2, deltasumlo=deltasumlo2; // invalidate logain of MExtractedSignalPix +// timelo=timelo2, deltatimelo=deltatimelo2; // invalidate logain of MArrivalTimePix +// satlo = satlo2; + + /* + // THIS IS NOW THE JOB OF THE EXTRACTOR! + //deltasumlo = 0; // make logain of MExtractedSignalPix valid + //deltatimelo = 0; // make logain of MArrivalTimePix valid + + const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1; + FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst, + sumlo, deltasumlo, timelo, deltatimelo, + satlo, pedpix, logainabflag); + timelo += fLoGainFirst; + */ + + + // If more than 6 lo-gain slices have saturated this is + // no physical event. We just assume that something with + // the extraction is wrong + if (numsatlo>6) + deltasumlo=deltatimelo=-1; + } + fLoGainFirst = fLoGainFirstSave; + + // Make sure that in cases the time couldn't be correctly determined + // more meaningfull default values are assigned + if (timelo<=fLoGainFirst || timelo>=fLoGainLast) + timelo = gRandom->Uniform(fLoGainLast-fLoGainFirst)+fLoGainFirst; + } + + MExtractedSignalPix &pix = (*fSignals)[pixidx]; + MArrivalTimePix &tix = (*fArrTime)[pixidx]; + pix.SetExtractedSignal(sumhi, deltasumhi, sumlo, deltasumlo); + pix.SetGainSaturation(numsathi, numsatlo); +// pix.SetGainSaturation(sathi, satlo); + + tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo); + tix.SetGainSaturation(numsathi, numsatlo); +// tix.SetGainSaturation(sathi, satlo); } /* while (pixel.Next()) */ - fArrTime->SetReadyToSave(); - fSignals->SetReadyToSave(); - - return kTRUE; + fArrTime->SetReadyToSave(); + fSignals->SetReadyToSave(); + + return kTRUE; } @@ -352,17 +533,17 @@ void MExtractTimeAndCharge::Print(Option_t *o) const { - if (IsA()==Class()) - *fLog << GetDescriptor() << ":" << endl; + MExtractTime::Print(o); +// if (IsA()==Class()) +// *fLog << GetDescriptor() << ":" << endl; *fLog << dec; - *fLog << " Taking " << fWindowSizeHiGain - << " HiGain samples from slice " << (Int_t)fHiGainFirst - << " to " << (Int_t)(fHiGainLast+fHiLoLast) << " incl" << endl; - *fLog << " Taking " << fWindowSizeLoGain - << " LoGain samples from slice " << (Int_t)fLoGainFirst - << " to " << (Int_t)fLoGainLast << " incl" << endl; - + //*fLog << " Taking " << fWindowSizeHiGain + // << " HiGain samples from slice " << (Int_t)fHiGainFirst + // << " to " << (Int_t)(fHiGainLast/*+fHiLoLast*/) << " incl" << endl; + //*fLog << " Taking " << fWindowSizeLoGain + // << " LoGain samples from slice " << (Int_t)fLoGainFirst + // << " to " << (Int_t)fLoGainLast << " incl" << endl; + // *fLog << " LoGainStartShift: " << fLoGainStartShift << endl; *fLog << " LoGainSwitch: " << (Int_t)fLoGainSwitch << endl; - MExtractTime::Print(o); } Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h (revision 8154) @@ -14,20 +14,18 @@ static const Byte_t fgLoGainSwitch; //! Default for fLoGainSwitch (now set to: 100) -protected: +private: - Byte_t fLoGainFirstSave; //! Temporary variable to store the original position of low-gain start slice +// Byte_t fLoGainFirstSave; //! Temporary variable to store the original position of low-gain start slice Float_t fLoGainStartShift; // Shift to start searching the low-gain signal obtained from the high-gain times. Byte_t fLoGainSwitch; // Limit for max. bin content before the low-gain gets extracted +protected: + Int_t fWindowSizeHiGain; // Window Size High-Gain Int_t fWindowSizeLoGain; // Window Size Low-Gain - Byte_t fMaxBinContent; // Maximum bin content - Int_t PreProcess(MParList *pList); Int_t Process(); - Bool_t ReInit(MParList *pList); - Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); @@ -47,7 +45,9 @@ virtual void SetWindowSize(Int_t windowh, Int_t windowl) { fWindowSizeHiGain = windowh; fWindowSizeLoGain = windowl; } + + Bool_t ReInit(MParList *pList); virtual Bool_t InitArrays() { return kTRUE; } - +/* virtual void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum, Float_t &time, Float_t &dtime, @@ -57,4 +57,14 @@ Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { } + */ + virtual void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t sat, Int_t maxpos) { } + + virtual void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t sat, Int_t maxpos) { } + + // For MExtractPedestal ClassDef(MExtractTimeAndCharge, 2) // Time And Charge Extractor Base Class Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 8154) @@ -1,3 +1,5 @@ /* ======================================================================== *\ +! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.73 2006-10-24 08:24:52 tbretz Exp $ +! -------------------------------------------------------------------------- ! ! * @@ -72,4 +74,5 @@ #include "MCalibrationPattern.h" #include "MExtractedSignalCam.h" +#include "MExtralgoDigitalFilter.h" #include "MPedestalPix.h" @@ -83,6 +86,6 @@ const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst = 1; const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast = 14; -const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain = 4; -const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain = 6; +//const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain = 4; +//const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain = 6; const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10; const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10; @@ -90,5 +93,5 @@ const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain = 4; const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain = 0.95; -const Float_t MExtractTimeAndChargeDigitalFilter::fgLoGainStartShift = -1.8; +//const Float_t MExtractTimeAndChargeDigitalFilter::fgLoGainStartShift = -1.8; // -------------------------------------------------------------------------- @@ -104,5 +107,7 @@ // MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title) - : fTimeShiftHiGain(0.), fTimeShiftLoGain(0.), fAutomaticWeights(kTRUE), fRandomIter(0) + : fBinningResolutionHiGain(fgBinningResolutionHiGain), + fBinningResolutionLoGain(fgBinningResolutionLoGain), + fAutomaticWeights(kTRUE), fRandomIter(0) { fName = name ? name : "MExtractTimeAndChargeDigitalFilter"; @@ -110,10 +115,10 @@ SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); - SetWindowSize(); - SetBinningResolution(); - SetSignalStartBin(); + SetWindowSize(3, 5); +// SetBinningResolution(); +// SetSignalStartBin(); SetOffsetLoGain(fgOffsetLoGain); // Order between both - SetLoGainStartShift(fgLoGainStartShift); // is important +// SetLoGainStartShift(fgLoGainStartShift); // is important } @@ -121,20 +126,16 @@ // // Checks: -// - if a window is bigger than the one defined by the ranges, set it to the available range +// - if a window is bigger than the one defined by the ranges, set it +// to the available range // // Sets: // - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain // - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain -// +// +// This function might be used to turn the digital filter into a +// sliding window extractor by setting the filename to NULL +// void MExtractTimeAndChargeDigitalFilter::SetWindowSize(Int_t windowh, Int_t windowl) { - - if (windowh != fgWindowSizeHiGain) - *fLog << warn << GetDescriptor() - << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default." << endl; - if (windowl != fgWindowSizeLoGain) - *fLog << warn << GetDescriptor() - << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default" << endl; - fWindowSizeHiGain = windowh; fWindowSizeLoGain = windowl; @@ -184,6 +185,6 @@ // size of 6. The exact numbers have to be found still. // - fNumHiGainSamples = fAmpWeightsHiGain.GetSum()/fBinningResolutionHiGain; - fNumLoGainSamples = fAmpWeightsLoGain.GetSum()/fBinningResolutionLoGain; + fNumHiGainSamples = fWindowSizeHiGain; + fNumLoGainSamples = fWindowSizeLoGain; fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); @@ -191,4 +192,5 @@ } + // -------------------------------------------------------------------------- // @@ -213,6 +215,4 @@ // // If new weights are set -// fTimeShiftHiGain -// fTimeShiftLoGain // fNumHiGainSamples // fNumLoGainSamples @@ -234,18 +234,10 @@ // - // shift the times back to the right reference (start counting from 0) + // We need here the effective number of samples. In pricipal the number + // is different depending on the weights used and must be set + // event by event. // - // The high-gain has the extraction offset (fHiGainFirst) already included - // here for speed reason. The low-gain has a changing extraction offset, - // so it will be added at every event (in FindTimeAndChargeLoGain) - fTimeShiftHiGain = 0.5 + 1./fBinningResolutionHiGain; - fTimeShiftLoGain = 0.5 + 1./fBinningResolutionLoGain; - - // - // We need here the effective number of samples which is about 2.5 in the case of a window - // size of 6. The exact numbers have to be found still. - // - fNumHiGainSamples = (Float_t)fWindowSizeHiGain/2.4; - fNumLoGainSamples = (Float_t)fWindowSizeLoGain/2.4; + fNumHiGainSamples = fAmpWeightsHiGain.GetSum()/fBinningResolutionHiGain; + fNumLoGainSamples = fAmpWeightsLoGain.GetSum()/fBinningResolutionLoGain; fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); @@ -253,5 +245,5 @@ // From MExtractTimeAndCharge::ReInit if (fSignals) - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples, + fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples, fLoGainFirst, fLoGainLast, fNumLoGainSamples); return kTRUE; @@ -269,9 +261,11 @@ return kFALSE; - const Int_t rangehi = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast); - const Int_t rangelo = (Int_t)(fLoGainLast - fLoGainFirst + 1); - - fHiGainSignal.Set(rangehi); - fLoGainSignal.Set(rangelo); + //const Int_t rangehi = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast); + //const Int_t rangelo = (Int_t)(fLoGainLast - fLoGainFirst + 1); + + //cout << "ARRAYS INITIALIZED TO : " << rangehi << " " << rangelo << endl; + + //fHiGainSignal.Set(rangehi); + //fLoGainSignal.Set(rangelo); return GetWeights(); @@ -295,220 +289,137 @@ } -void MExtractTimeAndChargeDigitalFilter::CalcBinningResArrays() -{ - - fArrBinningResHiGain.Set(fWindowSizeHiGain); - fArrBinningResHalfHiGain.Set(fWindowSizeHiGain); - - for (int i=0; i fMaxBinContent) - { - Byte_t maxpos = p-ptr; - - // range-fWindowSizeHiGain+1 == fHiLoLast isn't it? - if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1)) - fMaxBinContent = *p; - } - - if (*p++ >= fSaturationLimit) - if (!sat) - sat = ids-fHiGainFirst; - } - - // - // Slide with a window of size fWindowSizeHiGain over the sample - // and multiply the entries with the corresponding weights - // - if (IsNoiseCalculation()) - { - if (fRandomIter == fBinningResolutionHiGain) - fRandomIter = 0; - for (Int_t ids=0; ids < fWindowSizeHiGain; ids++) + Int_t range = fHiGainLast - fHiGainFirst + 1; + const Byte_t *end = first + range; + Byte_t *p = first; + + const Float_t pedes = ped.GetPedestal(); + const Float_t ABoffs = ped.GetPedestalABoffset(); + + const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; + + Byte_t maxcont = 0; + Int_t maxpos = 0; // obsolete for Digital Filter (used in spline!) + sat = 0; + + // + // Check for saturation in all other slices + // + Int_t ids = fHiGainFirst; + Float_t *sample = fHiGainSignal.GetArray(); + + while (p maxcont) { - const Int_t idx = fArrBinningResHiGain[ids] + fRandomIter; - sum += fAmpWeightsHiGain [idx]*fHiGainSignal[ids]; + maxpos = ids-fHiGainFirst-1; + // range-fWindowSizeHiGain+1 == fHiLoLast isn't it? + //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1)) + if (maxpos > 1 && maxpos < (range-1*//* - fWindowSizeHiGain + 1*//*)) + maxcont = *p; } - fRandomIter++; - return; - } - - if (fHiLoLast != 0) - { - - end = logain + fHiLoLast; - - while (logain= fSaturationLimit) + if (!sat) + sat = ids-fHiGainFirst; + } + + if (fHiLoLast != 0) + { + end = logain + fHiLoLast; + + while (logain= fSaturationLimit) - if (!sat) - sat = ids-fHiGainFirst; + *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1]; + + if (*logain > maxcont) + { + maxpos = ids-fHiGainFirst-1; + + // FIXME!!! MUST BE THERE! + // range-fWindowSizeHiGain+1 == fHiLoLast isn't it? + //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1)) + // maxcont = *logain; + } + + if (*logain++ >= fSaturationLimit) + if (!sat) + sat = ids-fHiGainFirst; + + range++; } } - // - // allow no saturated slice - // - // if (sat > 0) - // return; - - Float_t time_sum = 0.; - Float_t fmax = -FLT_MAX; - Float_t ftime_max = 0.; - Int_t max_p = -1; - - // - // Calculate the sum of the first fWindowSize slices - // - for (Int_t i=0;ifmax) - { - fmax = sum; - ftime_max = time_sum; - max_p = i; - } - } /* for (Int_t i=0;i fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain ) - { - if (t_iter > fBinningResolutionHalfHiGain-1) - { - t_iter -= fBinningResolutionHiGain; - max_p--; - sample_iter--; - } - if (t_iter < -fBinningResolutionHalfHiGain) - { - t_iter += fBinningResolutionHiGain; - max_p++; - sample_iter++; - } - } - - sum = 0.; - time_sum = 0.; - if (max_p < 0) - max_p = 0; - if (max_p+fWindowSizeHiGain > range) - max_p = range-fWindowSizeHiGain; - // - // Slide with a window of size fWindowSizeHiGain over the sample - // and multiply the entries with the corresponding weights - // - for (Int_t sample=0; sample < fWindowSizeHiGain; sample++) - { - const Int_t idx = fArrBinningResHalfHiGain[sample] + t_iter; - const Int_t ids = max_p + sample; - const Float_t pex = fHiGainSignal[ids]; - sum += fAmpWeightsHiGain [idx]*pex; - time_sum += fTimeWeightsHiGain[idx]*pex; - } - - if (sum == 0) - { - time = -1; - dsum = 0; - dtime = 0; - return; - } - - // here, the first high-gain slice is already included in the fTimeShiftHiGain -// time = fTimeShiftHiGain + max_p - Float_t(t_iter)/fBinningResolutionHiGain; - time = max_p + fTimeShiftHiGain + (Float_t)fHiGainFirst /* this shifts the time to the start of the rising edge */ - - ((Float_t)t_iter)/fBinningResolutionHiGain; - - const Float_t timefineadjust = time_sum/sum; - - if (TMath::Abs(timefineadjust) < 4./fBinningResolutionHiGain) - time -= timefineadjust; - + FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range, + sum, dsum, time, dtime, + sat, 0); } @@ -521,151 +432,37 @@ Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { - - const Int_t range = fLoGainLast - fLoGainFirst + 1; - - const Byte_t *end = ptr + range; - Byte_t *p = ptr; - // - // Prepare the low-gain pedestal - // - const Float_t pedes = ped.GetPedestal(); - const Float_t ABoffs = ped.GetPedestalABoffset(); - - const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; - - // - // Check for saturation in all other slices - // - Float_t *sample = fLoGainSignal.GetArray(); - Int_t ids = fLoGainFirst; - while (p= fSaturationLimit) - sat++; - } - - // - // Slide with a window of size fWindowSizeLoGain over the sample - // and multiply the entries with the corresponding weights - // - if (IsNoiseCalculation()) - { - if (fRandomIter == fBinningResolutionLoGain) - fRandomIter = 0; - for (Int_t ids=0; ids < fWindowSizeLoGain; ids++) - { - const Int_t idx = fArrBinningResLoGain[ids] + fRandomIter; - sum += fAmpWeightsLoGain [idx]*fLoGainSignal[ids]; - } - return; - } - - Float_t time_sum = 0.; - Float_t fmax = -FLT_MAX; - Float_t ftime_max = 0.; - Int_t max_p = -1; - - // - // Calculate the sum of the first fWindowSize slices - // - for (Int_t i=0;ifmax) - { - fmax = sum; - ftime_max = time_sum; - max_p = i; - } - } /* for (Int_t i=0;i fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain ) - { - if (t_iter > fBinningResolutionHalfLoGain-1) - { - t_iter -= fBinningResolutionLoGain; - max_p--; - sample_iter--; - } - if (t_iter < -fBinningResolutionHalfLoGain) - { - t_iter += fBinningResolutionLoGain; - max_p++; - sample_iter++; - } - } - - sum = 0.; - time_sum = 0.; - - // - // Slide with a window of size fWindowSizeLoGain over the sample - // and multiply the entries with the corresponding weights - // - for (Int_t sample=0; sample < fWindowSizeLoGain; sample++) - { - const Int_t idx = fArrBinningResHalfLoGain[sample] + t_iter; - const Int_t ids = max_p + sample; - const Float_t pex = ids < 0 ? 0. : ( ids >= range ? 0. : fLoGainSignal[ids]); - sum += fAmpWeightsLoGain [idx]*pex; - time_sum += fTimeWeightsLoGain[idx]*pex; - } - - if (sum == 0) - { - time = -1; - dsum = 0; - dtime = 0; - return; - } - - time = max_p + fTimeShiftLoGain + (Float_t)fLoGainFirst /* this shifts the time to the start of the rising edge */ - - ((Float_t)t_iter)/fBinningResolutionLoGain; - - const Float_t timefineadjust = time_sum/sum; - - if (TMath::Abs(timefineadjust) < 4./fBinningResolutionLoGain) - time -= timefineadjust; - -} - + const Int_t range = fLoGainLast - fLoGainFirst + 1; + + const Byte_t *end = ptr + range; + Byte_t *p = ptr; + + // + // Prepare the low-gain pedestal + // + const Float_t pedes = ped.GetPedestal(); + const Float_t ABoffs = ped.GetPedestalABoffset(); + + const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; + + // + // Check for saturation in all other slices + // + Float_t *sample = fLoGainSignal.GetArray(); + Int_t ids = fLoGainFirst; + + while (p= fSaturationLimit) + sat++; + } + + FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range, + sum, dsum, time, dtime, + sat, 0); + +} +*/ // -------------------------------------------------------------------------- // @@ -681,6 +478,4 @@ { - Byte_t hw = fWindowSizeHiGain; - Byte_t lw = fWindowSizeLoGain; Bool_t rc = kFALSE; @@ -690,4 +485,7 @@ rc = kTRUE; } + /* + Byte_t hw = fWindowSizeHiGain; + Byte_t lw = fWindowSizeLoGain; if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print)) @@ -704,5 +502,5 @@ if (rc) SetWindowSize(hw, lw); - + Bool_t rc2 = kFALSE; Int_t brh = fBinningResolutionHiGain; @@ -725,5 +523,5 @@ rc = kTRUE; } - + */ if (IsEnvDefined(env, prefix, "WeightsFile", print)) { @@ -823,5 +621,5 @@ } - if (2!=sscanf(str.Data(), "# High Gain Weights:%2i %2i", &fWindowSizeHiGain, &fBinningResolutionHiGain)) + if (2!=sscanf(str.Data(), "# High Gain Weights: %2i %2i", &fWindowSizeHiGain, &fBinningResolutionHiGain)) { *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl; @@ -833,4 +631,5 @@ fAmpWeightsHiGain .Set(len); fTimeWeightsHiGain.Set(len); + fPulseHiGain.Set(len); hi = kTRUE; continue; @@ -845,5 +644,5 @@ } - if (2!=sscanf(str.Data(),"# Low Gain Weights:%2i %2i", &fWindowSizeLoGain, &fBinningResolutionLoGain)) + if (2!=sscanf(str.Data(),"# Low Gain Weights: %2i %2i", &fWindowSizeLoGain, &fBinningResolutionLoGain)) { *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl; @@ -855,4 +654,5 @@ fAmpWeightsLoGain .Set(len); fTimeWeightsLoGain.Set(len); + fPulseLoGain.Set(len); lo = kTRUE; continue; @@ -867,7 +667,8 @@ continue; - if (2!=sscanf(str.Data(), "%f %f", + if (3!=sscanf(str.Data(), "%f %f %f", lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt], - lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt])) + lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt], + lo ? &fPulseLoGain[cnt] : &fPulseHiGain[cnt])) { *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl; @@ -903,5 +704,5 @@ *fLog << "with a resolution of " << fBinningResolutionLoGain << endl; - CalcBinningResArrays(); + //CalcBinningResArrays(); switch (fWindowSizeHiGain) @@ -1058,5 +859,5 @@ return ReadWeightsFile(name, path); } - +/* //---------------------------------------------------------------------------- // @@ -1125,5 +926,5 @@ Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1); - for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++) + for (Int_t i = -fBinningResolutionHiGain/2+1; i<=fBinningResolutionHiGain/2; i++) { @@ -1178,5 +979,5 @@ for (Int_t count=0; count < fWindowSizeHiGain; count++) { - const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1; + const Int_t idx = i+fBinningResolutionHiGain/2+fBinningResolutionHiGain*count-1; fAmpWeightsHiGain [idx] = w_amp [0][count]; fTimeWeightsHiGain[idx] = w_time[0][count]; @@ -1246,7 +1047,7 @@ // Loop over relative time in one BinningResolution interval // - Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain; - - for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++) + Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionLoGain/2; + + for (Int_t i = -fBinningResolutionLoGain/2+1; i<=fBinningResolutionLoGain/2; i++) { @@ -1301,5 +1102,5 @@ for (Int_t count=0; count < fWindowSizeLoGain; count++) { - const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1; + const Int_t idx = i+fBinningResolutionLoGain/2+fBinningResolutionLoGain*count-1; fAmpWeightsLoGain [idx] = w_amp [0][count]; fTimeWeightsLoGain[idx] = w_time[0][count]; @@ -1329,5 +1130,5 @@ return kTRUE; } - +*/ void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const { @@ -1336,5 +1137,4 @@ MExtractTimeAndCharge::Print(o); - *fLog << " Time Shift HiGain: " << fTimeShiftHiGain << " LoGain: " << fTimeShiftLoGain << endl; *fLog << " Window Size HiGain: " << fWindowSizeHiGain << " LoGain: " << fWindowSizeLoGain << endl; *fLog << " Binning Res HiGain: " << fBinningResolutionHiGain << " LoGain: " << fBinningResolutionHiGain << endl; Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h (revision 8154) @@ -22,106 +22,113 @@ { private: - static const Byte_t fgHiGainFirst; //! Default for fHiGainFirst (now set to: 0) - static const Byte_t fgHiGainLast; //! Default for fHiGainLast (now set to:14) - static const Byte_t fgLoGainFirst; //! Default for fLoGainFirst (now set to: 3) - static const Byte_t fgLoGainLast; //! Default for fLoGainLast (now set to:14) - static const Int_t fgWindowSizeHiGain; //! Default for fWindowSizeHiGain (now set to: 6) - static const Int_t fgWindowSizeLoGain; //! Default for fWindowSizeLoGain (now set to: 6) - static const Int_t fgBinningResolutionHiGain; //! Default for fBinningResolutionHiGain (now set to: 10) - static const Int_t fgBinningResolutionLoGain; //! Default for fBinningResolutionLoGain (now set to: 10) - static const Int_t fgSignalStartBinHiGain; //! Default for fSignalStartBinHiGain (now set to: 4) - static const Int_t fgSignalStartBinLoGain; //! Default for fSignalStartBinLoGain (now set to: 4) - static const TString fgNameWeightsFile; //! "cosmics_weights.dat" - static const Float_t fgOffsetLoGain; //! Default for fOffsetLoGain (now set to 1.7) - static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to -1.8) + static const Byte_t fgHiGainFirst; //! Default for fHiGainFirst (now set to: 0) + static const Byte_t fgHiGainLast; //! Default for fHiGainLast (now set to:14) + static const Byte_t fgLoGainFirst; //! Default for fLoGainFirst (now set to: 3) + static const Byte_t fgLoGainLast; //! Default for fLoGainLast (now set to:14) +// static const Int_t fgWindowSizeHiGain; //! Default for fWindowSizeHiGain (now set to: 6) +// static const Int_t fgWindowSizeLoGain; //! Default for fWindowSizeLoGain (now set to: 6) + static const Int_t fgBinningResolutionHiGain; //! Default for fBinningResolutionHiGain (now set to: 10) + static const Int_t fgBinningResolutionLoGain; //! Default for fBinningResolutionLoGain (now set to: 10) + static const Int_t fgSignalStartBinHiGain; //! Default for fSignalStartBinHiGain (now set to: 4) + static const Int_t fgSignalStartBinLoGain; //! Default for fSignalStartBinLoGain (now set to: 4) + static const TString fgNameWeightsFile; //! "cosmics_weights.dat" + static const Float_t fgOffsetLoGain; //! Default for fOffsetLoGain (now set to 1.7) +// static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to -1.8) - MCalibrationPattern *fCalibPattern; //! Calibration DM pattern - - MArrayF fHiGainSignal; //! Need fast access to the signals in a float way - MArrayF fLoGainSignal; //! Store them in separate arrays + MCalibrationPattern *fCalibPattern; //! Calibration DM pattern - MArrayI fArrBinningResHiGain; //! helping arrays to hold binningres * i - MArrayI fArrBinningResLoGain; //! helping arrays to hold binningres * i - MArrayI fArrBinningResHalfHiGain; //! helping arrays to hold binningres * i + binningreshalf - MArrayI fArrBinningResHalfLoGain; //! helping arrays to hold binningres * i + binningreshalf - - Float_t fTimeShiftHiGain; // Time shift from when on to apply the filter - Float_t fTimeShiftLoGain; // Time shift from when on to apply the filter - - Int_t fSignalStartBinHiGain; //! Start bin from when on to apply weights - Int_t fSignalStartBinLoGain; //! Start bin from when on to apply weights +// MArrayF fHiGainSignal; //! Need fast access to the signals in a float way +// MArrayF fLoGainSignal; //! Store them in separate arrays - Int_t fBinningResolutionHiGain; // Number of weights per bin High-Gain - Int_t fBinningResolutionHalfHiGain; // Half Number of weights per bin High-Gain - Int_t fBinningResolutionLoGain; // Number of weights per bin Low-Gain - Int_t fBinningResolutionHalfLoGain; // Half Number of weights per bin Low-Gain - - MArrayF fAmpWeightsHiGain; //! Amplitude weights High-Gain (from weights file) - MArrayF fTimeWeightsHiGain; //! Time weights High-Gain (from weights file) - MArrayF fAmpWeightsLoGain; //! Amplitude weights Low-Gain (from weights file) - MArrayF fTimeWeightsLoGain; //! Time weights Low-Gain (from weights file) +// Int_t fSignalStartBinHiGain; //! Start bin from when on to apply weights +// Int_t fSignalStartBinLoGain; //! Start bin from when on to apply weights - TString fNameWeightsFile; // Name of the weights file - Bool_t fAutomaticWeights; // Flag whether weight should be determined automatically - TString fNameWeightsFileSet; //! Flag if weights have alreayd been set - Int_t fRandomIter; //! Counter used to randomize weights for noise calculation + Int_t fBinningResolutionHiGain; // Number of weights per bin High-Gain + Int_t fBinningResolutionLoGain; // Number of weights per bin Low-Gain - // MExtractTimeAndChargeDigitalFilter - void CalcBinningResArrays(); - Int_t GetAutomaticWeights(); - Bool_t GetWeights(); - Int_t ReadWeightsFile(TString filename, TString path=""); - TString CompileWeightFileName(TString path, const TString &name) const; + MArrayF fAmpWeightsHiGain; //! Amplitude weights High-Gain (from weights file) + MArrayF fTimeWeightsHiGain; //! Time weights High-Gain (from weights file) + MArrayF fAmpWeightsLoGain; //! Amplitude weights Low-Gain (from weights file) + MArrayF fTimeWeightsLoGain; //! Time weights Low-Gain (from weights file) - // MExtractTimeAndCharge - Bool_t InitArrays(); + MArrayF fPulseHiGain; //! + MArrayF fPulseLoGain; //! - // MTask - Int_t PreProcess(MParList *pList); - Int_t Process(); + TString fNameWeightsFile; // Name of the weights file + Bool_t fAutomaticWeights; // Flag whether weight should be determined automatically + TString fNameWeightsFileSet; //! Flag if weights have alreayd been set + Int_t fRandomIter; //! Counter used to randomize weights for noise calculation + + // MExtractTimeAndChargeDigitalFilter + void CalcBinningResArrays(); + Int_t GetAutomaticWeights(); + Bool_t GetWeights(); + Int_t ReadWeightsFile(TString filename, TString path=""); + TString CompileWeightFileName(TString path, const TString &name) const; + + + // MExtractTimeAndCharge + Bool_t InitArrays(); + + // MTask + Int_t PreProcess(MParList *pList); + Int_t Process(); protected: - // MParContainer - Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); public: - MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL); - ~MExtractTimeAndChargeDigitalFilter() { } + MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL); + ~MExtractTimeAndChargeDigitalFilter() { } +/* + Bool_t WriteWeightsFile(TString filename, + TH1F *shapehi, TH2F *autocorrhi, + TH1F *shapelo=NULL, TH2F *autocorrlo=NULL ); - Bool_t WriteWeightsFile(TString filename, - TH1F *shapehi, TH2F *autocorrhi, - TH1F *shapelo=NULL, TH2F *autocorrlo=NULL ); + */ + void SetNameWeightsFile(TString s="") + { + s.ReplaceAll("\015", ""); // This is a fix for TEnv files edited with windows editors + fNameWeightsFile = s; + fNameWeightsFileSet=""; + } + void EnableAutomaticWeights(Bool_t b=kTRUE) { fAutomaticWeights = b; } - void SetNameWeightsFile(TString s="") { s.ReplaceAll("\015", ""); // This is a fix for TEnv files edited with windows editors - fNameWeightsFile = s; fNameWeightsFileSet=""; } - void EnableAutomaticWeights(Bool_t b=kTRUE) { fAutomaticWeights = b; } + void SetBinningResolution(const Int_t rh=fgBinningResolutionHiGain, const Int_t rl=fgBinningResolutionLoGain) + { + fBinningResolutionHiGain = rh & ~1; + fBinningResolutionLoGain = rl & ~1; + } +/* + void SetSignalStartBin( const Int_t sh=fgSignalStartBinHiGain, const Int_t sl=fgSignalStartBinLoGain) + { + fSignalStartBinHiGain = sh; + fSignalStartBinLoGain = sl; + } - void SetBinningResolution(const Int_t rh=fgBinningResolutionHiGain, const Int_t rl=fgBinningResolutionLoGain) { - fBinningResolutionHiGain = rh & ~1; - fBinningResolutionHalfHiGain = fBinningResolutionHiGain/2; - fBinningResolutionLoGain = rl & ~1; - fBinningResolutionHalfLoGain = fBinningResolutionLoGain/2; - } - - void SetSignalStartBin( const Int_t sh=fgSignalStartBinHiGain, const Int_t sl=fgSignalStartBinLoGain) { - fSignalStartBinHiGain = sh; - fSignalStartBinLoGain = sl; - } + */ + void SetWindowSize( Int_t windowh, Int_t windowl); + const char* GetNameWeightsFile() const { return fNameWeightsFile.Data(); } - void SetWindowSize( Int_t windowh=fgWindowSizeHiGain, Int_t windowl=fgWindowSizeLoGain); + void Print(Option_t *o="") const; //*MENU* +/* + void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); + void FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); + */ + void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t sat, Int_t maxpos); - const char* GetNameWeightsFile() const { return fNameWeightsFile.Data(); } + void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t sat, Int_t maxpos); - void Print(Option_t *o="") const; //*MENU* - - void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum, - Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); - void FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum, Float_t &dsum, - Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); - - ClassDef(MExtractTimeAndChargeDigitalFilter, 2) // Hendrik's digital filter + ClassDef(MExtractTimeAndChargeDigitalFilter, 2) // Hendrik's digital filter }; Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 8154) @@ -1,3 +1,5 @@ /* ======================================================================== *\ +! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeSpline.cc,v 1.61 2006-10-24 08:24:52 tbretz Exp $ +! -------------------------------------------------------------------------- ! ! * @@ -15,7 +17,8 @@ ! * ! -! Author(s): Markus Gaug 09/2004 +! Author(s): Thomas Bretz +! Author(s): Markus Gaug 09/2004 ! -! Copyright: MAGIC Software Development, 2002-2004 +! Copyright: MAGIC Software Development, 2002-2006 ! ! @@ -142,4 +145,6 @@ #include "MExtractTimeAndChargeSpline.h" +#include "MExtralgoSpline.h" + #include "MPedestalPix.h" @@ -156,9 +161,9 @@ const Byte_t MExtractTimeAndChargeSpline::fgLoGainLast = 14; const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.05; -const Float_t MExtractTimeAndChargeSpline::fgRiseTimeHiGain = 0.5; -const Float_t MExtractTimeAndChargeSpline::fgFallTimeHiGain = 0.5; +const Float_t MExtractTimeAndChargeSpline::fgRiseTimeHiGain = 0.7; +const Float_t MExtractTimeAndChargeSpline::fgFallTimeHiGain = 1.0; const Float_t MExtractTimeAndChargeSpline::fgLoGainStretch = 1.5; const Float_t MExtractTimeAndChargeSpline::fgOffsetLoGain = 1.3; -const Float_t MExtractTimeAndChargeSpline::fgLoGainStartShift = -2.4; +//const Float_t MExtractTimeAndChargeSpline::fgLoGainStartShift = -2.4; // -------------------------------------------------------------------------- @@ -177,6 +182,5 @@ // MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title) - : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.), - fRandomIter(0), fExtractionType(kIntegral) + : fRandomIter(0), fExtractionType(kIntegral) { @@ -187,5 +191,5 @@ SetLoGainStretch(); SetOffsetLoGain(fgOffsetLoGain); - SetLoGainStartShift(fgLoGainStartShift); +// SetLoGainStartShift(fgLoGainStartShift); SetRiseTimeHiGain(); @@ -310,7 +314,7 @@ { - Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast; - - fHiGainSignal .Set(range); + Int_t range = fHiGainLast - fHiGainFirst + 1;// + fHiLoLast; + +// fHiGainSignal .Set(range); fHiGainFirstDeriv .Set(range); fHiGainSecondDeriv.Set(range); @@ -318,13 +322,13 @@ range = fLoGainLast - fLoGainFirst + 1; - fLoGainSignal .Set(range); +// fLoGainSignal .Set(range); fLoGainFirstDeriv .Set(range); fLoGainSecondDeriv.Set(range); - fHiGainSignal .Reset(); +// fHiGainSignal .Reset(); fHiGainFirstDeriv .Reset(); fHiGainSecondDeriv.Reset(); - fLoGainSignal .Reset(); +// fLoGainSignal .Reset(); fLoGainFirstDeriv .Reset(); fLoGainSecondDeriv.Reset(); @@ -361,362 +365,132 @@ } +void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain2(const Float_t *ptr, Int_t num, + Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t sat, Int_t maxpos) +{ + // const Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast; + + // Do some handling if maxpos is last slice! + + MExtralgoSpline s(ptr, num, fHiGainFirstDeriv.GetArray(), fHiGainSecondDeriv.GetArray()); + + s.SetRiseFallTime(fRiseTimeHiGain, fFallTimeHiGain); + s.SetResolution(fResolution); + + if (IsNoiseCalculation()) + { + if (fRandomIter == int(1./fResolution)) + fRandomIter = 0; + + sum = s.ExtractNoise(fRandomIter); + fRandomIter++; + return; + } + + s.Extract(sat, maxpos); + s.GetTime(time, dtime); + s.GetSignal(sum, dsum); +} + +void MExtractTimeAndChargeSpline::FindTimeAndChargeLoGain2(const Float_t *ptr, Int_t num, + Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t sat, Int_t maxpos) +{ + MExtralgoSpline s(ptr, num, fLoGainFirstDeriv.GetArray(), fLoGainSecondDeriv.GetArray()); + + s.SetRiseFallTime(fRiseTimeLoGain, fFallTimeLoGain); + s.SetResolution(fResolution); + + if (IsNoiseCalculation()) + { + if (fRandomIter == int(1./fResolution)) + fRandomIter = 0; + + sum = s.ExtractNoise(fRandomIter); + return; + } + + s.Extract(sat, maxpos); + s.GetTime(time, dtime); + s.GetSignal(sum, dsum); +} + // -------------------------------------------------------------------------- // -// Calculates the arrival time and charge for each pixel -// -void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum, +// Calculates the arrival time and charge for each pixel +// +/* +void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { - Int_t range = fHiGainLast - fHiGainFirst + 1; - const Byte_t *end = first + range; - Byte_t *p = first; - - sat = 0; - dsum = 0; // In all cases the extracted signal is valid - - const Float_t pedes = ped.GetPedestal(); - const Float_t ABoffs = ped.GetPedestalABoffset(); - - const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; - - fAbMax = 0.; - fAbMaxPos = 0.; - fHalfMax = 0.; - fMaxBinContent = 0; - Int_t maxpos = 0; - - // - // Check for saturation in all other slices - // - Int_t ids = fHiGainFirst; - Float_t *sample = fHiGainSignal.GetArray(); - while (p fMaxBinContent) + Int_t range = fHiGainLast - fHiGainFirst + 1; + const Byte_t *end = first + range; + Byte_t *p = first; + + const Float_t pedes = ped.GetPedestal(); + const Float_t ABoffs = ped.GetPedestalABoffset(); + + const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; + + Byte_t maxcont = 0; + Int_t maxpos = 0; + sat = 0; + + // + // Check for saturation in all other slices + // + Int_t ids = fHiGainFirst; + Float_t *sample = fHiGainSignal.GetArray(); + while (p maxcont) { - maxpos = ids-fHiGainFirst-1; - // range-fWindowSizeHiGain+1 == fHiLoLast isn't it? - if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1)) - fMaxBinContent = *p; + maxpos = ids-fHiGainFirst-1; + // range-fWindowSizeHiGain+1 == fHiLoLast isn't it? + //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1)) + if (maxpos > 1 && maxpos < (range-1*//* - fWindowSizeHiGain + 1*//*)) + maxcont = *p; } - if (*p++ >= fSaturationLimit) - if (!sat) - sat = ids-fHiGainFirst; - - } - - if (fHiLoLast != 0) - { - - end = logain + fHiLoLast; - - while (logain= fSaturationLimit) + if (!sat) + sat = ids-fHiGainFirst; + } + + if (fHiLoLast != 0) + { + end = logain + fHiLoLast; + + while (logain fMaxBinContent) - { + *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1]; + + if (*logain > maxcont) + { maxpos = ids-fHiGainFirst-1; + + // FIXME!!! MUST BE THERE! // range-fWindowSizeHiGain+1 == fHiLoLast isn't it? //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1)) - // fMaxBinContent = *logain; + // maxcont = *logain; } - - if (*logain++ >= fSaturationLimit) - if (!sat) - sat = ids-fHiGainFirst; - - range++; + + if (*logain++ >= fSaturationLimit) + if (!sat) + sat = ids-fHiGainFirst; + + range++; } } - fAbMax = fHiGainSignal[maxpos]; - - fHiGainSecondDeriv[0] = 0.; - fHiGainFirstDeriv[0] = 0.; - - for (Int_t i=1;i=0;k--) - fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k]; - for (Int_t k=range-2;k>=0;k--) - fHiGainSecondDeriv[k] /= 6.; - - if (IsNoiseCalculation()) - { - - if (fRandomIter == int(1./fResolution)) - fRandomIter = 0; - - const Float_t nsx = fRandomIter * fResolution; - - if (fExtractionType == kAmplitude) - { - const Float_t b = nsx; - const Float_t a = 1. - nsx; - - sum = a*fHiGainSignal[1] - + b*fHiGainSignal[2] - + (a*a*a-a)*fHiGainSecondDeriv[1] - + (b*b*b-b)*fHiGainSecondDeriv[2]; - } - else - sum = CalcIntegralHiGain(2. + nsx, range); - - fRandomIter++; - return; - } - - // - // Allow no saturated slice and - // Don't start if the maxpos is too close to the limits. - // - const Bool_t limlo = maxpos < TMath::Ceil(fRiseTimeHiGain); - const Bool_t limup = maxpos > range-TMath::Ceil(fFallTimeHiGain)-1; - if (sat || limlo || limup) - { - dtime = 1.0; - if (fExtractionType == kAmplitude) - { - sum = fAbMax; - time = (Float_t)(fHiGainFirst + maxpos); - return; - } - - sum = CalcIntegralHiGain(limlo ? 0 : range, range); - time = (Float_t)(fHiGainFirst + maxpos - 1); - return; - } - - dtime = fResolution; - - // - // Now find the maximum - // - Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one - Float_t lower = -1. + maxpos; - Float_t upper = (Float_t)maxpos; - fAbMaxPos = upper; - Float_t x = lower; - Float_t y = 0.; - Float_t a = 1.; - Float_t b = 0.; - Int_t klo = maxpos-1; - Int_t khi = maxpos; - - // - // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2. - // If no maximum is found, go to interval maxpos+1. - // - while ( x < upper - 0.3 ) - { - - x += step; - a -= step; - b += step; - - y = a*fHiGainSignal[klo] - + b*fHiGainSignal[khi] - + (a*a*a-a)*fHiGainSecondDeriv[klo] - + (b*b*b-b)*fHiGainSecondDeriv[khi]; - - if (y > fAbMax) - { - fAbMax = y; - fAbMaxPos = x; - } - - } - - // - // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2 - // - if (fAbMaxPos > upper-0.1) - { - upper = 1. + maxpos; - lower = (Float_t)maxpos; - x = lower; - a = 1.; - b = 0.; - khi = maxpos+1; - klo = maxpos; - - while (x fAbMax) - { - fAbMax = y; - fAbMaxPos = x; - } - } - } - // - // Now, the time, abmax and khicont and klocont are set correctly within the previous precision. - // Try a better precision. - // - const Float_t up = fAbMaxPos+step - 3.0*fResolution; - const Float_t lo = fAbMaxPos-step + 3.0*fResolution; - const Float_t maxpossave = fAbMaxPos; - - x = fAbMaxPos; - a = upper - x; - b = x - lower; - - step = 2.*fResolution; // step size of 0.1 FADC slices - - while (x fAbMax) - { - fAbMax = y; - fAbMaxPos = x; - } - } - - // - // Second, try from time down to time-0.2 in steps of fResolution. - // - x = maxpossave; - - // - // Test the possibility that the absolute maximum has not been found between - // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos - // which requires new setting of klocont and khicont - // - if (x < lower + fResolution) - { - klo--; - khi--; - upper -= 1.; - lower -= 1.; - } - - a = upper - x; - b = x - lower; - - while (x>lo) - { - - x -= step; - a += step; - b -= step; - - y = a*fHiGainSignal[klo] - + b*fHiGainSignal[khi] - + (a*a*a-a)*fHiGainSecondDeriv[klo] - + (b*b*b-b)*fHiGainSecondDeriv[khi]; - - if (y > fAbMax) - { - fAbMax = y; - fAbMaxPos = x; - } - } - - if (fExtractionType == kAmplitude) - { - time = fAbMaxPos + (Int_t)fHiGainFirst; - sum = fAbMax; - return; - } - - fHalfMax = fAbMax/2.; - - // - // Now, loop from the maximum bin leftward down in order to find the position of the half maximum. - // First, find the right FADC slice: - // - klo = maxpos; - while (klo > 0) - { - if (fHiGainSignal[--klo] < fHalfMax) - break; - } - - khi = klo+1; - // - // Loop from the beginning of the slice upwards to reach the fHalfMax: - // With means of bisection: - // - x = (Float_t)klo; - a = 1.; - b = 0.; - - step = 0.5; - Bool_t back = kFALSE; - - Int_t maxcnt = 20; - Int_t cnt = 0; - - while (TMath::Abs(y-fHalfMax) > fResolution) - { - - if (back) - { - x -= step; - a += step; - b -= step; - } - else - { - x += step; - a -= step; - b += step; - } - - y = a*fHiGainSignal[klo] - + b*fHiGainSignal[khi] - + (a*a*a-a)*fHiGainSecondDeriv[klo] - + (b*b*b-b)*fHiGainSecondDeriv[khi]; - - back = y > fHalfMax; - - if (++cnt > maxcnt) - break; - - step /= 2.; - } - - // - // Now integrate the whole thing! - // - time = (Float_t)fHiGainFirst + x; - sum = CalcIntegralHiGain(fAbMaxPos - fRiseTimeHiGain, range); -} - + FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range, + sum, dsum, time, dtime, + sat, maxpos); +} // -------------------------------------------------------------------------- @@ -728,425 +502,40 @@ Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { - Int_t range = fLoGainLast - fLoGainFirst + 1; - const Byte_t *end = first + range; - Byte_t *p = first; - - const Float_t pedes = ped.GetPedestal(); - const Float_t ABoffs = ped.GetPedestalABoffset(); - - const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; - - fAbMax = 0.; - fAbMaxPos = 0.; - Int_t maxpos = 0; - Int_t max = -9999; - - dsum = 0; // In all cases the extracted signal is valid - - // - // Check for saturation in all other slices - // - Int_t ids = fLoGainFirst; - Float_t *sample = fLoGainSignal.GetArray(); - while (p max) + Int_t range = fLoGainLast - fLoGainFirst + 1; + const Byte_t *end = first + range; + Byte_t *p = first; + + const Float_t pedes = ped.GetPedestal(); + const Float_t ABoffs = ped.GetPedestalABoffset(); + + const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; + + Int_t maxpos = 0; + Int_t max = -9999; + + // + // Check for saturation in all other slices + // + Int_t ids = fLoGainFirst; + Float_t *sample = fLoGainSignal.GetArray(); + while (p max) { - maxpos = ids-fLoGainFirst-1; - max = *p; + maxpos = ids-fLoGainFirst-1; + max = *p; } - if (*p++ >= fSaturationLimit) - sat++; - } - - fAbMax = fLoGainSignal[maxpos]; - - fLoGainSecondDeriv[0] = 0.; - fLoGainFirstDeriv[0] = 0.; - - for (Int_t i=1;i=0;k--) - fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k]; - for (Int_t k=range-2;k>=0;k--) - fLoGainSecondDeriv[k] /= 6.; - - if (IsNoiseCalculation()) - { - if (fRandomIter == int(1./fResolution)) - fRandomIter = 0; - - const Float_t nsx = fRandomIter * fResolution; - - if (fExtractionType == kAmplitude) - { - const Float_t b = nsx; - const Float_t a = 1. - nsx; - - sum = a*fLoGainSignal[1] - + b*fLoGainSignal[2] - + (a*a*a-a)*fLoGainSecondDeriv[1] - + (b*b*b-b)*fLoGainSecondDeriv[2]; - } - else - sum = CalcIntegralLoGain(2. + nsx, range); - - fRandomIter++; - return; - } - // - // Allow no saturated slice and - // Don't start if the maxpos is too close to the limits. - // - const Bool_t limlo = maxpos < TMath::Ceil(fRiseTimeLoGain); - const Bool_t limup = maxpos > range-TMath::Ceil(fFallTimeLoGain)-1; - if (sat || limlo || limup) - { - dtime = 1.0; - if (fExtractionType == kAmplitude) - { - time = (Float_t)(fLoGainFirst + maxpos); - sum = fAbMax; - return; - } - - sum = CalcIntegralLoGain(limlo ? 0 : range, range); - time = (Float_t)(fLoGainFirst + maxpos - 1); - return; - } - - dtime = fResolution; - - // - // Now find the maximum - // - Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one - Float_t lower = -1. + maxpos; - Float_t upper = (Float_t)maxpos; - fAbMaxPos = upper; - Float_t x = lower; - Float_t y = 0.; - Float_t a = 1.; - Float_t b = 0.; - Int_t klo = maxpos-1; - Int_t khi = maxpos; - - // - // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2. - // If no maximum is found, go to interval maxpos+1. - // - while ( x < upper - 0.3 ) - { - - x += step; - a -= step; - b += step; - - y = a*fLoGainSignal[klo] - + b*fLoGainSignal[khi] - + (a*a*a-a)*fLoGainSecondDeriv[klo] - + (b*b*b-b)*fLoGainSecondDeriv[khi]; - - if (y > fAbMax) - { - fAbMax = y; - fAbMaxPos = x; - } - - } - - // - // Test the possibility that the absolute maximum has not been found before the - // maxpos and search from maxpos to maxpos+1 in steps of 0.2 - // - if (fAbMaxPos > upper-0.1) - { - - upper = 1. + maxpos; - lower = (Float_t)maxpos; - x = lower; - a = 1.; - b = 0.; - khi = maxpos+1; - klo = maxpos; - - while (x fAbMax) - { - fAbMax = y; - fAbMaxPos = x; - } - } - } - - - // - // Now, the time, abmax and khicont and klocont are set correctly within the previous precision. - // Try a better precision. - // - const Float_t up = fAbMaxPos+step - 3.0*fResolution; - const Float_t lo = fAbMaxPos-step + 3.0*fResolution; - const Float_t maxpossave = fAbMaxPos; - - x = fAbMaxPos; - a = upper - x; - b = x - lower; - - step = 2.*fResolution; // step size of 0.1 FADC slice - - while (x fAbMax) - { - fAbMax = y; - fAbMaxPos = x; - } - } - - // - // Second, try from time down to time-0.2 in steps of 0.025. - // - x = maxpossave; - - // - // Test the possibility that the absolute maximum has not been found between - // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos - // which requires new setting of klocont and khicont - // - if (x < lower + fResolution) - { - klo--; - khi--; - upper -= 1.; - lower -= 1.; - } - - a = upper - x; - b = x - lower; - - while (x>lo) - { - - x -= step; - a += step; - b -= step; - - y = a*fLoGainSignal[klo] - + b*fLoGainSignal[khi] - + (a*a*a-a)*fLoGainSecondDeriv[klo] - + (b*b*b-b)*fLoGainSecondDeriv[khi]; - - if (y > fAbMax) - { - fAbMax = y; - fAbMaxPos = x; - } - } - - if (fExtractionType == kAmplitude) - { - time = fAbMaxPos + (Int_t)fLoGainFirst; - sum = fAbMax; - return; - } - - fHalfMax = fAbMax/2.; - - // - // Now, loop from the maximum bin leftward down in order to find the position of the half maximum. - // First, find the right FADC slice: - // - klo = maxpos; - while (klo > 0) - { - klo--; - if (fLoGainSignal[klo] < fHalfMax) - break; - } - - khi = klo+1; - // - // Loop from the beginning of the slice upwards to reach the fHalfMax: - // With means of bisection: - // - x = (Float_t)klo; - a = 1.; - b = 0.; - - step = 0.5; - Bool_t back = kFALSE; - - Int_t maxcnt = 20; - Int_t cnt = 0; - - while (TMath::Abs(y-fHalfMax) > fResolution) - { - - if (back) - { - x -= step; - a += step; - b -= step; - } - else - { - x += step; - a -= step; - b += step; - } - - y = a*fLoGainSignal[klo] - + b*fLoGainSignal[khi] - + (a*a*a-a)*fLoGainSecondDeriv[klo] - + (b*b*b-b)*fLoGainSecondDeriv[khi]; - - back = y > fHalfMax; - - if (++cnt > maxcnt) - break; - - step /= 2.; - } - - // - // Now integrate the whole thing! - // - time = x + (Int_t)fLoGainFirst; - sum = CalcIntegralLoGain(fAbMaxPos - fRiseTimeLoGain, range); -} - -Float_t MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t start, Float_t range) const -{ - // The number of steps is calculated directly from the integration - // window. This is the only way to ensure we are not dealing with - // numerical rounding uncertanties, because we always get the same - // value under the same conditions -- it might still be different on - // other machines! - const Float_t step = 0.2; - const Float_t width = fRiseTimeHiGain+fFallTimeHiGain; - const Float_t max = range-1 - (width+step); - const Int_t num = TMath::Nint(width/step); - - // The order is important. In some cases (limlo-/limup-check) it can - // happen than max<0. In this case we start at 0 - if (start > max) - start = max; - if (start < 0) - start = 0; - - start += step/2; - - Double_t sum = 0.; - for (Int_t i=0; i max) - start = max; - if (start < 0) - start = 0; - - start += step/2; - - Double_t sum = 0.; - for (Int_t i=0; i= fSaturationLimit) + sat++; + } + + FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range, + sum, dsum, time, dtime, + sat, maxpos); +} +*/ // -------------------------------------------------------------------------- @@ -1198,6 +587,3 @@ return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc; - -} - - +} Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h (revision 8154) @@ -11,99 +11,95 @@ class MPedestalPix; + class MExtractTimeAndChargeSpline : public MExtractTimeAndCharge { +private: + static const Byte_t fgHiGainFirst; //! Default for fHiGainFirst (now set to: 2) + static const Byte_t fgHiGainLast; //! Default for fHiGainLast (now set to: 14) + static const Byte_t fgLoGainFirst; //! Default for fLoGainFirst (now set to: 2) + static const Byte_t fgLoGainLast; //! Default for fLoGainLast (now set to: 14) + static const Float_t fgResolution; //! Default for fResolution (now set to: 0.003) + static const Float_t fgRiseTimeHiGain; //! Default for fRiseTime (now set to: 1.5) + static const Float_t fgFallTimeHiGain; //! Default for fFallTime (now set to: 4.5) + static const Float_t fgLoGainStretch; //! Default for fLoGainStretch (now set to: 1.5) + static const Float_t fgOffsetLoGain; //! Default for fOffsetLoGain (now set to 1.7) +// static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to -1.6) + +// MArrayF fHiGainSignal; //! Need fast access to the signals in a float way +// MArrayF fLoGainSignal; //! Store them in separate arrays + MArrayF fHiGainFirstDeriv; //! High-gain discretized first derivatives + MArrayF fLoGainFirstDeriv; //! Low-gain discretized first derivatives + MArrayF fHiGainSecondDeriv; //! High-gain discretized second derivatives + MArrayF fLoGainSecondDeriv; //! Low-gain discretized second derivatives + + Float_t fResolution; // The time resolution in FADC units + + Float_t fRiseTimeHiGain; // The usual rise time of the pulse in the high-gain + Float_t fFallTimeHiGain; // The usual fall time of the pulse in the high-gain + Float_t fRiseTimeLoGain; // The usual rise time of the pulse in the low-gain + Float_t fFallTimeLoGain; // The usual fall time of the pulse in the low-gain + + Float_t fLoGainStretch; // The stretch of the low-gain w.r.t. the high-gain pulse + + Int_t fRandomIter; //! Counter used to randomize weights for noise calculation + + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + Bool_t InitArrays(); + +public: + enum ExtractionType_t { kAmplitude, kIntegral }; //! Possible time and charge extraction types private: - - static const Byte_t fgHiGainFirst; //! Default for fHiGainFirst (now set to: 2) - static const Byte_t fgHiGainLast; //! Default for fHiGainLast (now set to: 14) - static const Byte_t fgLoGainFirst; //! Default for fLoGainFirst (now set to: 2) - static const Byte_t fgLoGainLast; //! Default for fLoGainLast (now set to: 14) - static const Float_t fgResolution; //! Default for fResolution (now set to: 0.003) - static const Float_t fgRiseTimeHiGain; //! Default for fRiseTime (now set to: 1.5) - static const Float_t fgFallTimeHiGain; //! Default for fFallTime (now set to: 4.5) - static const Float_t fgLoGainStretch; //! Default for fLoGainStretch (now set to: 1.5) - static const Float_t fgOffsetLoGain; //! Default for fOffsetLoGain (now set to 1.7) - static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to -1.6) - - MArrayF fHiGainSignal; //! Need fast access to the signals in a float way - MArrayF fLoGainSignal; //! Store them in separate arrays - MArrayF fHiGainFirstDeriv; //! High-gain discretized first derivatives - MArrayF fLoGainFirstDeriv; //! Low-gain discretized first derivatives - MArrayF fHiGainSecondDeriv; //! High-gain discretized second derivatives - MArrayF fLoGainSecondDeriv; //! Low-gain discretized second derivatives - - Float_t fAbMax; //! Current maximum of the spline - Float_t fAbMaxPos; //! Current position of the maximum of the spline - Float_t fHalfMax; //! Current half maximum of the spline - - Float_t fResolution; // The time resolution in FADC units - - Float_t fRiseTimeHiGain; // The usual rise time of the pulse in the high-gain - Float_t fFallTimeHiGain; // The usual fall time of the pulse in the high-gain - Float_t fRiseTimeLoGain; // The usual rise time of the pulse in the low-gain - Float_t fFallTimeLoGain; // The usual fall time of the pulse in the low-gain - - Float_t fLoGainStretch; // The stretch of the low-gain w.r.t. the high-gain pulse - - Int_t fRandomIter; //! Counter used to randomize weights for noise calculation - - Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); - - Bool_t InitArrays(); - - Float_t CalcIntegralHiGain(Float_t start, Float_t range) const; - Float_t CalcIntegralLoGain(Float_t start, Float_t range) const; - -public: - enum ExtractionType_t { kAmplitude, kIntegral }; //! Possible time and charge extraction types - -private: - - ExtractionType_t fExtractionType; + ExtractionType_t fExtractionType; public: + MExtractTimeAndChargeSpline(const char *name=NULL, const char *title=NULL); + ~MExtractTimeAndChargeSpline() {} - MExtractTimeAndChargeSpline(const char *name=NULL, const char *title=NULL); - ~MExtractTimeAndChargeSpline() {} + Float_t GetRiseTimeHiGain() const { return fRiseTimeHiGain; } + Float_t GetFallTimeHiGain() const { return fFallTimeHiGain; } - Float_t GetRiseTimeHiGain() const { return fRiseTimeHiGain; } - Float_t GetFallTimeHiGain() const { return fFallTimeHiGain; } + void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0 ); - void SetRange ( Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0 ); - void SetResolution ( const Float_t f=fgResolution ) { fResolution = f; } - void SetRiseTimeHiGain( const Float_t f=fgRiseTimeHiGain ) + void SetResolution(const Float_t f=fgResolution) { fResolution = f; } + + void SetRiseTimeHiGain(const Float_t f=fgRiseTimeHiGain) { - fRiseTimeHiGain = f; - fRiseTimeLoGain = f*fLoGainStretch; - fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); - fWindowSizeHiGain = TMath::Nint(fRiseTimeHiGain + fFallTimeHiGain); + fRiseTimeHiGain = f; + fRiseTimeLoGain = f*fLoGainStretch; + fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); + fWindowSizeHiGain = TMath::Nint(fRiseTimeHiGain + fFallTimeHiGain); } - void SetFallTimeHiGain( const Float_t f=fgFallTimeHiGain ) + void SetFallTimeHiGain(const Float_t f=fgFallTimeHiGain) { - fFallTimeHiGain = f; - fFallTimeLoGain = f*fLoGainStretch; - fNumHiGainSamples = fRiseTimeHiGain + fFallTimeHiGain; - fNumLoGainSamples = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.; - fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); - fWindowSizeLoGain = TMath::Nint(fRiseTimeLoGain + fFallTimeLoGain); + fFallTimeHiGain = f; + fFallTimeLoGain = f*fLoGainStretch; + fNumHiGainSamples = fRiseTimeHiGain + fFallTimeHiGain; + fNumLoGainSamples = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.; + fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); + fWindowSizeLoGain = TMath::Nint(fRiseTimeLoGain + fFallTimeLoGain); } - void SetLoGainStretch ( const Float_t f=fgLoGainStretch ) { fLoGainStretch = f; } - - void SetChargeType ( const ExtractionType_t typ=kIntegral); - - void FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum, - Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); - void FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum, Float_t &dsum, - Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); + void SetLoGainStretch(const Float_t f=fgLoGainStretch) { fLoGainStretch = f; } - ClassDef(MExtractTimeAndChargeSpline, 4) // Task to Extract Arrival Times and Charges using a Cubic Spline + void SetChargeType(const ExtractionType_t typ=kIntegral); +/* + void FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); + void FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); +*/ + void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t sat, Int_t maxpos); + + void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t sat, Int_t maxpos); + + ClassDef(MExtractTimeAndChargeSpline, 4) // Task to Extract Arrival Times and Charges using a Cubic Spline }; #endif - - - Index: trunk/MagicSoft/Mars/msignal/MExtractedSignalPix.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractedSignalPix.h (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractedSignalPix.h (revision 8154) @@ -15,6 +15,9 @@ Float_t fExtractedSignalLoGainError; // error of the mean value of the extracted signal - Byte_t fNumHiGainSaturated; // Number of first hi-gain slice which has saturated (could be negative if already the first slice saturates) - Byte_t fNumLoGainSaturated; // Number of first lo-gain slices which have saturated + Byte_t fNumHiGainSaturated; // Number of how many hi-gain slices saturated + Byte_t fNumLoGainSaturated; // Number of how many lo-gain slices saturated + + Byte_t GetNumHiGainSaturated() const { return fNumHiGainSaturated; } + Byte_t GetNumLoGainSaturated() const { return fNumLoGainSaturated; } public: @@ -36,7 +39,4 @@ Float_t GetExtractedSignalLoGainError() const { return fExtractedSignalLoGainError; } - Byte_t GetNumHiGainSaturated() const { return fNumHiGainSaturated; } - Byte_t GetNumLoGainSaturated() const { return fNumLoGainSaturated; } - Bool_t IsHiGainSaturated() const { return fNumHiGainSaturated>0; } Bool_t IsLoGainSaturated() const { return fNumLoGainSaturated>0; } Index: trunk/MagicSoft/Mars/msignal/MExtractor.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractor.cc (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractor.cc (revision 8154) @@ -93,4 +93,5 @@ #include "MPedestalCam.h" #include "MPedestalPix.h" +//#include "MPedestalSubtractEvt.h" #include "MExtractedSignalCam.h" @@ -101,5 +102,5 @@ using namespace std; -const Byte_t MExtractor::fgSaturationLimit = 250; +const Byte_t MExtractor::fgSaturationLimit = 245; const TString MExtractor::fgNamePedestalCam = "MPedestalCam"; const TString MExtractor::fgNameSignalCam = "MExtractedSignalCam"; @@ -123,6 +124,7 @@ MExtractor::MExtractor(const char *name, const char *title) : fResolutionPerPheHiGain(0), fResolutionPerPheLoGain(0), - fPedestals(NULL), fSignals(NULL), fRawEvt(NULL), fRunHeader(NULL), - fHiLoLast(0), fNumHiGainSamples(0), fNumLoGainSamples(0) + fPedestals(NULL), fSignals(NULL), fRawEvt(NULL), fRunHeader(NULL), + fSignal(NULL), + /*fHiLoLast(0),*/ fNumHiGainSamples(0), fNumLoGainSamples(0) { fName = name ? name : "MExtractor"; @@ -172,4 +174,11 @@ { *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl; + return kFALSE; + } + + fSignal = (MPedestalSubtractedEvt*)pList->FindObject(AddSerialNumber("MPedestalSubtractedEvt")); + if (!fSignal) + { + *fLog << err << AddSerialNumber("MPedestalSubtractedEvt") << " not found... aborting." << endl; return kFALSE; } @@ -231,5 +240,20 @@ Bool_t MExtractor::ReInit(MParList *pList) { - + const Int_t numl = fRunHeader->GetNumSamplesLoGain(); + const Int_t numh = fRunHeader->GetNumSamplesHiGain(); + const Int_t num = numh+numl; + + if (fHiGainLast>=num) + { + *fLog << err << "ERROR - Last hi-gain slice must not exceed " << num-1 << endl; + return kFALSE; + } + if (fLoGainLast>=num) + { + *fLog << err << "ERROR - Last lo-gain slice must not exceed " << num-1 << endl; + return kFALSE; + } + +/* const Int_t logainsamples = fRunHeader->GetNumSamplesLoGain(); @@ -289,5 +313,6 @@ } - return kTRUE; + */ + return kTRUE; } @@ -431,5 +456,5 @@ *fLog << " Hi Gain Range: " << Form("%2d %2d", fHiGainFirst, fHiGainLast) << endl; *fLog << " Lo Gain Range: " << Form("%2d %2d", fLoGainFirst, fLoGainLast) << endl; - *fLog << " Gain Overlap to Lo: " << Form("%2d", fHiLoLast) << endl; +// *fLog << " Gain Overlap to Lo: " << Form("%2d", fHiLoLast) << endl; *fLog << " Saturation Lim: " << Form("%3d", fSaturationLimit) << endl; *fLog << " Num Samples Hi/Lo: " << Form("%2.1f %2.1f", fNumHiGainSamples, fNumLoGainSamples) << endl; Index: trunk/MagicSoft/Mars/msignal/MExtractor.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractor.h (revision 8153) +++ trunk/MagicSoft/Mars/msignal/MExtractor.h (revision 8154) @@ -18,4 +18,5 @@ class MPedestalCam; +class MPedestalSubtractedEvt; class MExtractedSignalCam; @@ -42,4 +43,6 @@ MRawEvtData *fRawEvt; //! Raw event data (time slices) MRawRunHeader *fRunHeader; //! RunHeader information + + MPedestalSubtractedEvt *fSignal; //! Byte_t fHiGainFirst; // First FADC slice nr. to extract the High Gain signal @@ -48,5 +51,5 @@ Byte_t fLoGainLast; // Last FADC slice nr. to extract the Low Gain signal - Byte_t fHiLoLast; // Number of slices in fLoGainSamples counted for the High-Gain signal +// Byte_t fHiLoLast; // Number of slices in fLoGainSamples counted for the High-Gain signal Float_t fNumHiGainSamples; // Number High Gain FADC slices used to extract the signal @@ -91,4 +94,5 @@ Float_t GetResolutionPerPheHiGain() const { return fResolutionPerPheHiGain; } Float_t GetResolutionPerPheLoGain() const { return fResolutionPerPheLoGain; } + Byte_t GetSaturationLimit() const { return fSaturationLimit; } Bool_t IsNoiseCalculation () const { return fNoiseCalculation; } @@ -107,9 +111,10 @@ void SetPedestals (MPedestalCam *pedcam) { fPedestals = pedcam; } + MPedestalCam *GetPedestals() { return fPedestals; } // TObject void Clear(Option_t *o="") { - fHiGainFirst = fHiGainLast = fLoGainFirst = fLoGainLast = fHiLoLast = 0; + fHiGainFirst = fHiGainLast = fLoGainFirst = fLoGainLast = /*fHiLoLast =*/ 0; } Index: trunk/MagicSoft/Mars/msignal/Makefile =================================================================== --- trunk/MagicSoft/Mars/msignal/Makefile (revision 8153) +++ trunk/MagicSoft/Mars/msignal/Makefile (revision 8154) @@ -21,5 +21,5 @@ INCLUDES = -I. -I../mbase -I../mgui -I../mraw -I../manalysis \ -I../mgeom -I../mtools -I../mpedestal -I../mbadpixels \ - -I../mcalib + -I../mcalib -I../mextralgo # mgui (MCamEvent): MExtractSignalCam Index: trunk/MagicSoft/Mars/msignal/calibration_weights_UV46.dat =================================================================== --- trunk/MagicSoft/Mars/msignal/calibration_weights_UV46.dat (revision 8153) +++ trunk/MagicSoft/Mars/msignal/calibration_weights_UV46.dat (revision 8154) @@ -1,104 +1,104 @@ # High Gain Weights: 4 10 # (Amplitude) (Time) --0.735581 -0.673266 --0.800386 -0.47606 --0.902801 -0.192719 --0.945945 0.237423 --0.851618 0.770424 --0.598107 1.28444 --0.290663 1.69546 -0.0967702 1.89144 -0.362204 1.98161 -0.610111 2.06506 -0.919902 2.20752 -1.12282 1.94486 -1.46305 1.70309 -1.85728 1.58904 -2.19155 1.21011 -2.37154 0.603168 -2.37187 0.000360131 -2.25639 -0.535921 -2.12495 -0.948388 -1.97328 -1.286 -1.72155 -1.70358 -1.55218 -1.75991 -1.21489 -1.86382 -0.759117 -2.09403 -0.266577 -2.03167 --0.102906 -1.60487 --0.284545 -1.16641 --0.367327 -0.679947 --0.376403 -0.3943 --0.379568 -0.229929 --0.218954 0.0443019 --0.22659 0.163101 --0.195144 0.249649 --0.136182 0.328464 --0.0541723 0.334035 --0.00249719 0.207255 -0.00794832 0.102902 -0.00785676 0.0396167 -0.00459524 0.0111356 -0.00887564 0.0405 +-0.735581 -0.673266 0 +-0.800386 -0.47606 0 +-0.902801 -0.192719 0 +-0.945945 0.237423 0 +-0.851618 0.770424 0 +-0.598107 1.28444 0 +-0.290663 1.69546 0 +0.0967702 1.89144 0 +0.362204 1.98161 0 +0.610111 2.06506 0 +0.919902 2.20752 0 +1.12282 1.94486 0 +1.46305 1.70309 0 +1.85728 1.58904 0 +2.19155 1.21011 0 +2.37154 0.603168 0 +2.37187 0.000360131 0 +2.25639 -0.535921 0 +2.12495 -0.948388 0 +1.97328 -1.286 0 +1.72155 -1.70358 0 +1.55218 -1.75991 0 +1.21489 -1.86382 0 +0.759117 -2.09403 0 +0.266577 -2.03167 0 +-0.102906 -1.60487 0 +-0.284545 -1.16641 0 +-0.367327 -0.679947 0 +-0.376403 -0.3943 0 +-0.379568 -0.229929 0 +-0.218954 0.0443019 0 +-0.22659 0.163101 0 +-0.195144 0.249649 0 +-0.136182 0.328464 0 +-0.0541723 0.334035 0 +-0.00249719 0.207255 0 +0.00794832 0.102902 0 +0.00785676 0.0396167 0 +0.00459524 0.0111356 0 +0.00887564 0.0405 0 # Low Gain Weights: 6 10 # (Amplitude) (Time) -0.0446612 -0.385273 -0.038191 -0.00418687 -0.0386966 0.0212324 -0.0402881 0.0744799 -0.0415794 0.229615 -0.0598731 0.44332 -0.0758477 0.661518 -0.101509 1.10641 -0.159323 1.64997 -0.497256 2.83685 -0.245087 3.27499 -0.140546 2.46177 -0.58086 2.2849 -0.632721 2.45587 -0.72819 2.52835 -0.889583 2.48099 -0.980812 2.50031 -1.09885 2.55892 -1.21374 2.78769 -1.61928 3.08069 -1.38544 1.95583 -1.31998 1.1792 -1.50633 0.591226 -1.50916 0.0793899 -1.5008 -0.33188 -1.47339 -0.575386 -1.45362 -0.915309 -1.40214 -1.31593 -1.34175 -1.77904 -1.0661 -2.05471 -1.31087 -1.49798 -1.33793 -1.34758 -1.10172 -1.21719 -1.08133 -1.09356 -1.04007 -0.981455 -0.976745 -1.08299 -0.930979 -1.14774 -0.874203 -1.18348 -0.816708 -1.20126 -0.587354 -1.92869 -0.783078 -1.89621 -0.792771 -1.03439 -0.622278 -0.781807 -0.61184 -0.745831 -0.578792 -0.683741 -0.537336 -0.596328 -0.51443 -0.592858 -0.482294 -0.560586 -0.462351 -0.827587 -0.317989 -1.05649 -0.459672 -0.775035 -0.468287 -0.619961 -0.374182 -0.31635 -0.376946 -0.225242 -0.367075 -0.347444 -0.340737 -0.393231 -0.321054 -0.187384 -0.320654 -0.225558 -0.302148 -0.399499 -0.232954 -0.607578 +0.0446612 -0.385273 0 +0.038191 -0.00418687 0 +0.0386966 0.0212324 0 +0.0402881 0.0744799 0 +0.0415794 0.229615 0 +0.0598731 0.44332 0 +0.0758477 0.661518 0 +0.101509 1.10641 0 +0.159323 1.64997 0 +0.497256 2.83685 0 +0.245087 3.27499 0 +0.140546 2.46177 0 +0.58086 2.2849 0 +0.632721 2.45587 0 +0.72819 2.52835 0 +0.889583 2.48099 0 +0.980812 2.50031 0 +1.09885 2.55892 0 +1.21374 2.78769 0 +1.61928 3.08069 0 +1.38544 1.95583 0 +1.31998 1.1792 0 +1.50633 0.591226 0 +1.50916 0.0793899 0 +1.5008 -0.33188 0 +1.47339 -0.575386 0 +1.45362 -0.915309 0 +1.40214 -1.31593 0 +1.34175 -1.77904 0 +1.0661 -2.05471 0 +1.31087 -1.49798 0 +1.33793 -1.34758 0 +1.10172 -1.21719 0 +1.08133 -1.09356 0 +1.04007 -0.981455 0 +0.976745 -1.08299 0 +0.930979 -1.14774 0 +0.874203 -1.18348 0 +0.816708 -1.20126 0 +0.587354 -1.92869 0 +0.783078 -1.89621 0 +0.792771 -1.03439 0 +0.622278 -0.781807 0 +0.61184 -0.745831 0 +0.578792 -0.683741 0 +0.537336 -0.596328 0 +0.51443 -0.592858 0 +0.482294 -0.560586 0 +0.462351 -0.827587 0 +0.317989 -1.05649 0 +0.459672 -0.775035 0 +0.468287 -0.619961 0 +0.374182 -0.31635 0 +0.376946 -0.225242 0 +0.367075 -0.347444 0 +0.340737 -0.393231 0 +0.321054 -0.187384 0 +0.320654 -0.225558 0 +0.302148 -0.399499 0 +0.232954 -0.607578 0 Index: trunk/MagicSoft/Mars/msignal/cosmics_weights46.dat =================================================================== --- trunk/MagicSoft/Mars/msignal/cosmics_weights46.dat (revision 8153) +++ trunk/MagicSoft/Mars/msignal/cosmics_weights46.dat (revision 8154) @@ -1,104 +1,104 @@ # High Gain Weights: 4 10 # (Amplitude) (Time) --0.637443 -0.648763 --0.707555 -0.610204 --0.772465 -0.469694 --0.896081 -0.305143 --1.02783 0.00696361 --1.03684 0.731979 --0.73698 1.49428 --0.23419 1.57509 -0.12767 1.83675 -0.441087 1.70526 -0.639658 1.50546 -0.806192 1.58667 -0.966865 1.41151 -1.24058 1.34506 -1.58293 1.39453 -2.02092 1.46038 -2.31508 0.8539 -2.31705 0.051831 -2.18497 -0.493294 -2.01129 -0.749924 -1.8675 -0.959329 -1.75322 -1.25026 -1.60962 -1.37231 -1.3216 -1.58832 -0.880364 -1.9924 -0.219661 -2.50675 --0.35077 -1.96955 --0.553454 -0.855454 --0.491715 -0.387673 --0.38701 -0.0297507 --0.293549 0.092153 --0.278646 0.194742 --0.263716 0.25195 --0.221238 0.323789 --0.121941 0.381038 -0.00379863 0.393009 -0.0805122 0.201843 -0.0443888 -0.0531952 --0.0199606 -0.175625 --0.0809393 -0.21419 +-0.637443 -0.648763 0 +-0.707555 -0.610204 0 +-0.772465 -0.469694 0 +-0.896081 -0.305143 0 +-1.02783 0.00696361 0 +-1.03684 0.731979 0 +-0.73698 1.49428 0 +-0.23419 1.57509 0 +0.12767 1.83675 0 +0.441087 1.70526 0 +0.639658 1.50546 0 +0.806192 1.58667 0 +0.966865 1.41151 0 +1.24058 1.34506 0 +1.58293 1.39453 0 +2.02092 1.46038 0 +2.31508 0.8539 0 +2.31705 0.051831 0 +2.18497 -0.493294 0 +2.01129 -0.749924 0 +1.8675 -0.959329 0 +1.75322 -1.25026 0 +1.60962 -1.37231 0 +1.3216 -1.58832 0 +0.880364 -1.9924 0 +0.219661 -2.50675 0 +-0.35077 -1.96955 0 +-0.553454 -0.855454 0 +-0.491715 -0.387673 0 +-0.38701 -0.0297507 0 +-0.293549 0.092153 0 +-0.278646 0.194742 0 +-0.263716 0.25195 0 +-0.221238 0.323789 0 +-0.121941 0.381038 0 +0.00379863 0.393009 0 +0.0805122 0.201843 0 +0.0443888 -0.0531952 0 +-0.0199606 -0.175625 0 +-0.0809393 -0.21419 0 # Low Gain Weights: 6 10 # (Amplitude) (Time) -0.0446612 -0.385273 -0.038191 -0.00418687 -0.0386966 0.0212324 -0.0402881 0.0744799 -0.0415794 0.229615 -0.0598731 0.44332 -0.0758477 0.661518 -0.101509 1.10641 -0.159323 1.64997 -0.497256 2.83685 -0.245087 3.27499 -0.140546 2.46177 -0.58086 2.2849 -0.632721 2.45587 -0.72819 2.52835 -0.889583 2.48099 -0.980812 2.50031 -1.09885 2.55892 -1.21374 2.78769 -1.61928 3.08069 -1.38544 1.95583 -1.31998 1.1792 -1.50633 0.591226 -1.50916 0.0793899 -1.5008 -0.33188 -1.47339 -0.575386 -1.45362 -0.915309 -1.40214 -1.31593 -1.34175 -1.77904 -1.0661 -2.05471 -1.31087 -1.49798 -1.33793 -1.34758 -1.10172 -1.21719 -1.08133 -1.09356 -1.04007 -0.981455 -0.976745 -1.08299 -0.930979 -1.14774 -0.874203 -1.18348 -0.816708 -1.20126 -0.587354 -1.92869 -0.783078 -1.89621 -0.792771 -1.03439 -0.622278 -0.781807 -0.61184 -0.745831 -0.578792 -0.683741 -0.537336 -0.596328 -0.51443 -0.592858 -0.482294 -0.560586 -0.462351 -0.827587 -0.317989 -1.05649 -0.459672 -0.775035 -0.468287 -0.619961 -0.374182 -0.31635 -0.376946 -0.225242 -0.367075 -0.347444 -0.340737 -0.393231 -0.321054 -0.187384 -0.320654 -0.225558 -0.302148 -0.399499 -0.232954 -0.607578 +0.0446612 -0.385273 0 +0.038191 -0.00418687 0 +0.0386966 0.0212324 0 +0.0402881 0.0744799 0 +0.0415794 0.229615 0 +0.0598731 0.44332 0 +0.0758477 0.661518 0 +0.101509 1.10641 0 +0.159323 1.64997 0 +0.497256 2.83685 0 +0.245087 3.27499 0 +0.140546 2.46177 0 +0.58086 2.2849 0 +0.632721 2.45587 0 +0.72819 2.52835 0 +0.889583 2.48099 0 +0.980812 2.50031 0 +1.09885 2.55892 0 +1.21374 2.78769 0 +1.61928 3.08069 0 +1.38544 1.95583 0 +1.31998 1.1792 0 +1.50633 0.591226 0 +1.50916 0.0793899 0 +1.5008 -0.33188 0 +1.47339 -0.575386 0 +1.45362 -0.915309 0 +1.40214 -1.31593 0 +1.34175 -1.77904 0 +1.0661 -2.05471 0 +1.31087 -1.49798 0 +1.33793 -1.34758 0 +1.10172 -1.21719 0 +1.08133 -1.09356 0 +1.04007 -0.981455 0 +0.976745 -1.08299 0 +0.930979 -1.14774 0 +0.874203 -1.18348 0 +0.816708 -1.20126 0 +0.587354 -1.92869 0 +0.783078 -1.89621 0 +0.792771 -1.03439 0 +0.622278 -0.781807 0 +0.61184 -0.745831 0 +0.578792 -0.683741 0 +0.537336 -0.596328 0 +0.51443 -0.592858 0 +0.482294 -0.560586 0 +0.462351 -0.827587 0 +0.317989 -1.05649 0 +0.459672 -0.775035 0 +0.468287 -0.619961 0 +0.374182 -0.31635 0 +0.376946 -0.225242 0 +0.367075 -0.347444 0 +0.340737 -0.393231 0 +0.321054 -0.187384 0 +0.320654 -0.225558 0 +0.302148 -0.399499 0 +0.232954 -0.607578 0