Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 8164) +++ /trunk/MagicSoft/Mars/Changelog (revision 8165) @@ -18,4 +18,38 @@ -*-*- END OF LINE -*-*- + 2006/10/25 Thomas Bretz + + * mjtrain/MJTrainSeparation.cc: + - print filename into output + + * mbadpixels/MBadPixelsCam.cc: + - turned the warning when merging an empty cam into an info + + * mextralgo/MExtralgoDigitalFilter.cc: + - removed some old comments, added some new ones + + * mextralgo/MExtralgoSpline.[h,cc]: + - added some sanity checks to make sure that the spline + can be initialized and is initialized + - changed the default extraction position for noise extraction to be + in the middle of a 5 slice long spline + + * msignal/MExtractTimeAndCharge.[h,cc]: + - removed a lot of old comment and added some new ones + - changed the behaviour of the determination of the start position + to extract the lo-gains + - added a spline interpolation to estimate the rrival time + of saturating hi-gains + - decreased the random range to the maximum extraction range + + * msignal/MExtractTimeAndChargeDigitalFilter.[h,cc], + msignal/MExtractTimeAndChargeSpline.[h,cc]: + - removed a lot of comments, added some new ones + - removed fraRandomIter (replaced by a real random number) + - replaced some if-conditions by switch + - init the arrays to the maximum possible size + + + 2006/10/25 Daniela Dorner Index: /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.cc (revision 8164) +++ /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.cc (revision 8165) @@ -136,5 +136,5 @@ // -------------------------------------------------------------------------- // -// Merge two MBadPixelsCam together, see also MBadPixelsPix::Merge +// Merge MBadPixelsCam cam into this, see also MBadPixelsPix::Merge // void MBadPixelsCam::Merge(const MBadPixelsCam &cam) @@ -143,5 +143,5 @@ if (n==0) { - *fLog << warn << "MBadPixelsCam::Merge: Container empty." << endl; + *fLog << inf << "MBadPixelsCam::Merge: Container empty." << endl; return; } Index: /trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc =================================================================== --- /trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc (revision 8164) +++ /trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc (revision 8165) @@ -210,28 +210,11 @@ // Take also a possible offset from timefineadjust into account // Sould it be: fTime += fWindowSize/2; ??? + + // HERE we should add the distance from the beginning of the + // extraction window to the leading edge! fTime += 0.5 + 0.5/fWeightsPerBin; - // In the ideal case the would never get <0 - // But timefineadjust can be >0.05 (in 60% of the cases) - // Define in each extractor a lowest and highest extracted value! - // here, the first high-gain slice is already included - // in the fTimeShiftHiGain this shifts the time to the - // start of the rising edge - - // This is from MExtractTimeAndChargeDigitalFilter - // fTime += 0.5 + 1./fWeightsPerBin; - - // Now there is a difference between the rising edge and - // the extracted time. This must be applied after - // the randomization in case of wrog values - const Float_t timefineadjust = sumtime/sumamp; - - // FIXME: Is 3 a good value? - //if (TMath::Abs(timefineadjust)*fWeightsPerBin < 3.) - // fTime -= timefineadjust; - - // if (TMath::FloorNint(timefineadjust+0.5)==0) //if (TMath::Abs(timefineadjust) < 0.2) Index: /trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc (revision 8164) +++ /trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc (revision 8165) @@ -69,4 +69,7 @@ void MExtralgoSpline::InitDerivatives() const { + if (fNum<2) + return; + fDer1[0] = 0.; fDer2[0] = 0.; @@ -211,4 +214,7 @@ Float_t MExtralgoSpline::ExtractNoise(/*Int_t iter*/) { + if (fNum<5) + return 0; + // FIXME: Shell we keep the extraction inside one slice // or randomize it along the extraction window? @@ -216,7 +222,7 @@ if (fExtractionType == kAmplitude) - return Eval(1, nsx); + return Eval(2, nsx); else - return CalcIntegral(2. + nsx); + return CalcIntegral(2 + nsx); } @@ -227,4 +233,7 @@ fSignalDev = -1; fTimeDev = -1; + + if (fNum<2) + return; /* // @@ -373,5 +382,5 @@ Float_t maxpos, maxval; - // FIXME: Check the dfeault if no maximum found!!! + // FIXME: Check the default if no maximum found!!! GetMaxAroundI(maxbin, maxpos, maxval); Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc (revision 8164) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc (revision 8165) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.55 2006-10-24 08:24:52 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.56 2006-10-25 18:36:26 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -69,4 +69,7 @@ #include +#include + +#include "MMath.h" #include "MLog.h" @@ -75,14 +78,6 @@ #include "MParList.h" -//#include "MArrayB.h" -//#include "MArrayF.h" - -//#include "MRawEvtData.h" #include "MRawEvtPixelIter.h" -//#include "MRawRunHeader.h" - -//#include "MPedestalCam.h" -//#include "MPedestalPix.h" -#include "MPedestalSubtractedEvt.h" +#include "MRawRunHeader.h" #include "MArrivalTimeCam.h" @@ -92,9 +87,11 @@ #include "MExtractedSignalPix.h" +#include "MPedestalSubtractedEvt.h" + ClassImp(MExtractTimeAndCharge); using namespace std; -const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -6.0; +const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -2.5; const Byte_t MExtractTimeAndCharge::fgLoGainSwitch = 120; @@ -105,9 +102,9 @@ // Sets: // - fWindowSizeHiGain and fWindowSizeLoGain to 0 -// - fLoGainStartShift to fgLoGainStartShift+fgOffsetLoGain +// - fLoGainStartShift to fgLoGainStartShift // - fLoGainSwitch to fgLoGainSwitch // MExtractTimeAndCharge::MExtractTimeAndCharge(const char *name, const char *title) - : /*fLoGainFirstSave(0),*/ fWindowSizeHiGain(0), fWindowSizeLoGain(0) + : fWindowSizeHiGain(0), fWindowSizeLoGain(0) { fName = name ? name : "MExtractTimeAndCharge"; @@ -133,14 +130,13 @@ Int_t MExtractTimeAndCharge::PreProcess(MParList *pList) { - - if (!MExtractTime::PreProcess(pList)) - return kFALSE; - - fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam)); - if (!fSignals) - return kFALSE; - - *fLog << flush << inf; - return kTRUE; + if (!MExtractTime::PreProcess(pList)) + return kFALSE; + + fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam)); + if (!fSignals) + return kFALSE; + + *fLog << flush << inf; + return kTRUE; } @@ -157,16 +153,14 @@ Bool_t MExtractTimeAndCharge::ReInit(MParList *pList) { - - if (!MExtractTime::ReInit(pList)) - return kFALSE; - - if (!InitArrays()) - return kFALSE; - - if (fSignals) - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples, - fLoGainFirst, fLoGainLast, fNumLoGainSamples); - - return kTRUE; + if (!MExtractTime::ReInit(pList)) + return kFALSE; + + if (!InitArrays(fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain())) + return kFALSE; + + if (fSignals) + fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples, + fLoGainFirst, fLoGainLast, fNumLoGainSamples); + return kTRUE; } @@ -188,101 +182,10 @@ 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 Float_t *sig = fSignal->GetSamples(pixidx); + 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); - } - - // 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 sathi0 = fHiGainFirst; // First slice to extract and first saturating slice + Int_t sathi1 = fHiGainLast; // Last slice to extract and last saturating slice Int_t maxcont; @@ -290,68 +193,46 @@ // 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) + + // Do not even try to extract the hi-gain if we have + // more than one saturating slice + if (numsathi<2) { - - const Int_t rangehi = fHiGainLast - fHiGainFirst + 1/* + fHiLoLast*/; - FindTimeAndChargeHiGain2(sig/*.GetArray()*/+fHiGainFirst, rangehi, - sumhi, deltasumhi, timehi, deltatimehi, - numsathi, maxposhi); - - timehi += fHiGainFirst; + const Int_t rangehi = fHiGainLast - fHiGainFirst + 1; + FindTimeAndChargeHiGain2(sig+fHiGainFirst, rangehi, + sumhi, deltasumhi, timehi, deltatimehi, + numsathi, maxposhi); + + + // Make sure that in cases the time couldn't be correctly determined + // more meaningfull default values are assigned + //if (timehi=fHiGainLast-1) + if (deltatimehi>-0.5 && (timehi<0 || timehi>=fHiGainLast-fHiGainFirst+1)) + timehi = gRandom->Uniform(fHiGainLast-fHiGainFirst+1); + + 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; + + // If we have saturating slices try to get a better estimate + // of the arrival time than timehi or sathi0. This is + // usefull to know where to start lo-gain extraction. + if (numsathi>1) + { + const Int_t p = sathi0>1 ? sathi0-2 : sathi0-1; + if (sathi0>0) + { + // Find the place at which the signal is maxcont/2 + const TVector3 vx(sig[p], sig[p+1], sig[p+2]); + const TVector3 vy(p, p+1, p+2); + timehi=MMath::InterpolParabLin(vx, vy, maxcont/2); + } + else + timehi=0; + } + + Float_t sumlo =0, deltasumlo =-1; // invalidate logain of MExtractedSignalPix + Float_t timelo=0, deltatimelo=-1; // invalidate logain of MArrivalTimePix Int_t numsatlo=0; @@ -363,4 +244,6 @@ // MEANS: Hi gain has saturated, but the signal is to dim // to extract the lo-gain properly + // This could happen because the maxcont was not extracted from + // all slices! // THIS produces pulse positions ~= -1 // The signal might be handled in MCalibrateData, but hwat's about @@ -376,122 +259,75 @@ // FIXME: What to do with the pixel if it saturates too early??? - if (pixel.HasLoGain() && (maxcont > fLoGainSwitch /*|| sathi>0*/) ) + if (pixel.HasLoGain() && maxcont>fLoGainSwitch) { + Int_t first = numh+fLoGainFirst; + Int_t last = numh+fLoGainLast; + // 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; + // + // numh + fLoGainStartShift == 14 / fLoGainStartShift=-1 + // + // The lo-gain is expected to have its raising edge + // at timehi+numh+fOffsetLoGain. We start to search for the + // lo-gain fLoGainStartShift slices earlier. + // + // Instead of fLoGainStartShift the extractor should now how many + // slices before the expected raising edge the start of the + // search must be placed and from there we can step 1.5 slices + // back to be on the safe side. + first = TMath::FloorNint(timehi+numh+fOffsetLoGain+fLoGainStartShift); + if (first<0) + first = 0; + if (first>last) + first=last; + + /* + // Currently we have to stick to this check because at least + // the spline has arrays of this size... + if (first>last) + first = last; + if (firstGetMax(pixidx, satlo0, satlo1, maxlo); + numsatlo = fSignal->GetSaturation(pixidx, fSaturationLimit, satlo0, satlo1); + + FindTimeAndChargeLoGain2(sig+first, last-first+1, + sumlo, deltasumlo, timelo, deltatimelo, + numsatlo, maxposlo); // 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; - } - + //if (timehi=fHiGainLast-1) + if (deltatimelo>-0.5 && (timelo<0 || timelo>=last-first+1)) + timelo = gRandom->Uniform(last-first+1); + + timelo += first-numh; + + // 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; + + //if (TMath::Abs(timelo-fOffsetLoGain - timehi)>1.0) + // deltatimelo = -1; + } + + // Now store the result in the corresponding containers 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()) */ @@ -509,9 +345,5 @@ Int_t MExtractTimeAndCharge::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { - Bool_t rc = MExtractTime::ReadEnv(env, prefix, print); - - if (rc) - SetLoGainStartShift(); if (IsEnvDefined(env, prefix, "LoGainStartShift", print)) @@ -533,17 +365,8 @@ void MExtractTimeAndCharge::Print(Option_t *o) const { - 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 << " LoGainStartShift: " << fLoGainStartShift << endl; - *fLog << " LoGainSwitch: " << (Int_t)fLoGainSwitch << endl; -} + MExtractTime::Print(o); + + *fLog << dec; + *fLog << " LoGainStartShift: " << fLoGainStartShift << endl; + *fLog << " LoGainSwitch: " << (Int_t)fLoGainSwitch << endl; +} Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h (revision 8164) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h (revision 8165) @@ -7,20 +7,15 @@ class MPedestalPix; + class MExtractTimeAndCharge : public MExtractTime { private: - static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to: -2.8) static const Byte_t fgLoGainSwitch; //! Default for fLoGainSwitch (now set to: 100) -private: - -// 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 @@ -31,5 +26,4 @@ public: - MExtractTimeAndCharge(const char *name=NULL, const char *title=NULL); @@ -38,33 +32,24 @@ Byte_t GetLoGainSwitch () const { return fLoGainSwitch; } - void Print(Option_t *o="") const; //*MENU* - void SetLoGainStartShift( const Float_t f=fgLoGainStartShift ) { fLoGainStartShift = f + fOffsetLoGain; } void SetLoGainSwitch ( const Byte_t i=fgLoGainSwitch ) { fLoGainSwitch = i; } - virtual void SetWindowSize(Int_t windowh, Int_t windowl) { fWindowSizeHiGain = windowh; - fWindowSizeLoGain = windowl; } + 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, - Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { } - - virtual 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) { } - */ + virtual Bool_t InitArrays(Int_t n) { return kTRUE; } + 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) { } + Byte_t sat, Int_t maxpos) const { } 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) { } + Byte_t sat, Int_t maxpos) const { } // For MExtractPedestal + Bool_t ReInit(MParList *pList); + + // TObject + void Print(Option_t *o="") const; //*MENU* ClassDef(MExtractTimeAndCharge, 2) // Time And Charge Extractor Base Class Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 8164) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 8165) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.73 2006-10-24 08:24:52 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.74 2006-10-25 18:36:26 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -62,7 +62,5 @@ #include -#include -#include -#include +#include #include "MLog.h" @@ -76,22 +74,17 @@ #include "MExtralgoDigitalFilter.h" -#include "MPedestalPix.h" - ClassImp(MExtractTimeAndChargeDigitalFilter); using namespace std; -const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst = 0; -const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast = 15; -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::fgBinningResolutionHiGain = 10; -const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10; -const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain = 4; -const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain = 4; -const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain = 0.95; -//const Float_t MExtractTimeAndChargeDigitalFilter::fgLoGainStartShift = -1.8; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst = 0; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast = 15; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst = 1; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast = 14; +const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10; +const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10; +const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain = 4; +const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain = 4; +const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain = 0.95; // -------------------------------------------------------------------------- @@ -109,5 +102,5 @@ : fBinningResolutionHiGain(fgBinningResolutionHiGain), fBinningResolutionLoGain(fgBinningResolutionLoGain), - fAutomaticWeights(kTRUE), fRandomIter(0) + fAutomaticWeights(kTRUE) { fName = name ? name : "MExtractTimeAndChargeDigitalFilter"; @@ -116,9 +109,5 @@ SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); SetWindowSize(3, 5); -// SetBinningResolution(); -// SetSignalStartBin(); - - SetOffsetLoGain(fgOffsetLoGain); // Order between both -// SetLoGainStartShift(fgLoGainStartShift); // is important + SetOffsetLoGain(fgOffsetLoGain); } @@ -189,5 +178,4 @@ fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); - } @@ -256,16 +244,8 @@ // Gets called in the ReInit() and initialized the arrays // -Bool_t MExtractTimeAndChargeDigitalFilter::InitArrays() +Bool_t MExtractTimeAndChargeDigitalFilter::InitArrays(Int_t n) { if (!fRunHeader) return kFALSE; - - //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(); @@ -296,5 +276,5 @@ Float_t &sum, Float_t &dsum, Float_t &time, Float_t &dtime, - Byte_t sat, Int_t maxpos) + Byte_t sat, Int_t maxpos) const { // Do some handling if maxpos is last slice! @@ -308,9 +288,5 @@ if (IsNoiseCalculation()) { - if (fRandomIter == fBinningResolutionHiGain) - fRandomIter = 0; - - sum = df.ExtractNoise(fRandomIter); - fRandomIter++; + sum = df.ExtractNoise(gRandom->Integer(fBinningResolutionHiGain)); return; } @@ -324,5 +300,5 @@ Float_t &sum, Float_t &dsum, Float_t &time, Float_t &dtime, - Byte_t sat, Int_t maxpos) + Byte_t sat, Int_t maxpos) const { MExtralgoDigitalFilter df(fBinningResolutionLoGain, fWindowSizeLoGain, @@ -335,8 +311,5 @@ if (IsNoiseCalculation()) { - if (fRandomIter == fBinningResolutionLoGain) - fRandomIter = 0; - - sum = df.ExtractNoise(fRandomIter); + sum = df.ExtractNoise(gRandom->Integer(fBinningResolutionHiGain)); return; } @@ -347,129 +320,8 @@ } -/* + // -------------------------------------------------------------------------- // -// Apply the digital filter algorithm to the high-gain slices. -// -void MExtractTimeAndChargeDigitalFilter::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; - - 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) - { - 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 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++; - } - } - - FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range, - sum, dsum, time, dtime, - sat, 0); -} - -// -------------------------------------------------------------------------- -// -// Apply the digital filter algorithm to the low-gain slices. -// -void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum, - Float_t &time, Float_t &dtime, - 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++; - } - - FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range, - sum, dsum, time, dtime, - sat, 0); - -} -*/ -// -------------------------------------------------------------------------- -// // Read the setup from a TEnv, eg: -// MJPedestal.MExtractor.WindowSizeHiGain: 6 -// MJPedestal.MExtractor.WindowSizeLoGain: 6 -// MJPedestal.MExtractor.BinningResolutionHiGain: 10 -// MJPedestal.MExtractor.BinningResolutionLoGain: 10 // MJPedestal.MExtractor.WeightsFile: filename // MJPedestal.MExtractor.AutomaticWeights: off @@ -485,43 +337,5 @@ rc = kTRUE; } - /* - Byte_t hw = fWindowSizeHiGain; - Byte_t lw = fWindowSizeLoGain; - - if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print)) - { - hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw); - rc = kTRUE; - } - if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print)) - { - lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw); - rc = kTRUE; - } - - if (rc) - SetWindowSize(hw, lw); - - Bool_t rc2 = kFALSE; - Int_t brh = fBinningResolutionHiGain; - Int_t brl = fBinningResolutionLoGain; - - if (IsEnvDefined(env, prefix, "BinningResolutionHiGain", print)) - { - brh = GetEnvValue(env, prefix, brh); - rc2 = kTRUE; - } - if (IsEnvDefined(env, prefix, "BinningResolutionLoGain", print)) - { - brl = GetEnvValue(env, prefix, brl); - rc2 = kTRUE; - } - - if (rc2) - { - SetBinningResolution(brh, brl); - rc = kTRUE; - } - */ + if (IsEnvDefined(env, prefix, "WeightsFile", print)) { @@ -859,276 +673,10 @@ return ReadWeightsFile(name, path); } -/* + //---------------------------------------------------------------------------- // -// Create the weights file -// Beware that the shape-histogram has to contain the pulse starting at bin 1 -// -Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi, - TH1F *shapelo, TH2F *autocorrlo ) -{ - - const Int_t nbinshi = shapehi->GetNbinsX(); - Float_t binwidth = shapehi->GetBinWidth(1); - - TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"), - Form("%s%s",shapehi->GetTitle()," derivative"), - nbinshi, - shapehi->GetBinLowEdge(1), - shapehi->GetBinLowEdge(nbinshi)+binwidth); - - // - // Calculate the derivative of shapehi - // - for (Int_t i = 1; iSetBinContent(i, - ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth)); - derivativehi->SetBinError(i, - (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1) - +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth)); - } - - // - // normalize the shapehi, such that the integral for fWindowSize slices is one! - // - Float_t sum = 0; - Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain); - lasttemp = lasttemp > nbinshi ? nbinshi : lasttemp; - - for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; iGetBinContent(i); - } - sum /= fBinningResolutionHiGain; - - shapehi->Scale(1./sum); - derivativehi->Scale(1./sum); - - // - // read in the noise auto-correlation function: - // - TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain); - - for (Int_t i=0; iGetBinContent(i+1,j+1); //+fSignalStartBinHiGain +fSignalStartBinHiGain - } - } - Bhi.Invert(); - - const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain; - fAmpWeightsHiGain.Set(nsizehi); - fTimeWeightsHiGain.Set(nsizehi); - - // - // Loop over relative time in one BinningResolution interval - // - Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1); - - for (Int_t i = -fBinningResolutionHiGain/2+1; i<=fBinningResolutionHiGain/2; i++) - { - - TMatrix g(fWindowSizeHiGain,1); - TMatrix gT(1,fWindowSizeHiGain); - TMatrix d(fWindowSizeHiGain,1); - TMatrix dT(1,fWindowSizeHiGain); - - for (Int_t count=0; count < fWindowSizeHiGain; count++){ - - g[count][0]=shapehi->GetBinContent(start - +fBinningResolutionHiGain*count+i); - gT[0][count]=shapehi->GetBinContent(start - +fBinningResolutionHiGain*count+i); - d[count][0]=derivativehi->GetBinContent(start - +fBinningResolutionHiGain*count+i); - dT[0][count]=derivativehi->GetBinContent(start - +fBinningResolutionHiGain*count+i); - } - - TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g)); - Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix - - TMatrix m_first = dT*(Bhi*d); // ROOT thinks, m_first is still a matrix - Float_t first = m_first[0][0]/denom; - - TMatrix m_last = gT*(Bhi*d); // ROOT thinks, m_last is still a matrix - Float_t last = m_last[0][0]/denom; - - TMatrix m1 = gT*Bhi; - m1 *= first; - - TMatrix m2 = dT*Bhi; - m2 *=last; - - TMatrix w_amp = m1 - m2; - - TMatrix m_first1 = gT*(Bhi*g); - Float_t first1 = m_first1[0][0]/denom; - - TMatrix m_last1 = gT*(Bhi*d); - Float_t last1 = m_last1 [0][0]/denom; - - TMatrix m11 = dT*Bhi; - m11 *=first1; - - TMatrix m21 = gT*Bhi; - m21 *=last1; - - TMatrix w_time= m11 - m21; - - for (Int_t count=0; count < fWindowSizeHiGain; count++) - { - const Int_t idx = i+fBinningResolutionHiGain/2+fBinningResolutionHiGain*count-1; - fAmpWeightsHiGain [idx] = w_amp [0][count]; - fTimeWeightsHiGain[idx] = w_time[0][count]; - } - - } // end loop over i - - // - // Low Gain histograms - // - TH1F *derivativelo = NULL; - if (shapelo) - { - const Int_t nbinslo = shapelo->GetNbinsX(); - binwidth = shapelo->GetBinWidth(1); - - derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"), - Form("%s%s",shapelo->GetTitle()," derivative"), - nbinslo, - shapelo->GetBinLowEdge(1), - shapelo->GetBinLowEdge(nbinslo)+binwidth); - - // - // Calculate the derivative of shapelo - // - for (Int_t i = 1; iSetBinContent(i, - ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth)); - derivativelo->SetBinError(i, - (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1) - +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth)); - } - - // - // normalize the shapelo, such that the integral for fWindowSize slices is one! - // - sum = 0; - lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain); - lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp; - - for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; iGetBinContent(i); - - sum /= fBinningResolutionLoGain; - - shapelo->Scale(1./sum); - derivativelo->Scale(1./sum); - - // - // read in the noise auto-correlation function: - // - TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain); - - for (Int_t i=0; iGetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain); - } - } - Blo.Invert(); - - const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain; - fAmpWeightsLoGain.Set(nsizelo); - fTimeWeightsLoGain.Set(nsizelo); - - // - // Loop over relative time in one BinningResolution interval - // - Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionLoGain/2; - - for (Int_t i = -fBinningResolutionLoGain/2+1; i<=fBinningResolutionLoGain/2; i++) - { - - TMatrix g(fWindowSizeLoGain,1); - TMatrix gT(1,fWindowSizeLoGain); - TMatrix d(fWindowSizeLoGain,1); - TMatrix dT(1,fWindowSizeLoGain); - - for (Int_t count=0; count < fWindowSizeLoGain; count++){ - - g[count][0] = shapelo->GetBinContent(start - +fBinningResolutionLoGain*count+i); - gT[0][count]= shapelo->GetBinContent(start - +fBinningResolutionLoGain*count+i); - d[count][0] = derivativelo->GetBinContent(start - +fBinningResolutionLoGain*count+i); - dT[0][count]= derivativelo->GetBinContent(start - +fBinningResolutionLoGain*count+i); - } - - TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g)); - Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix - - TMatrix m_first = dT*(Blo*d); // ROOT thinks, m_first is still a matrix - Float_t first = m_first[0][0]/denom; - - TMatrix m_last = gT*(Blo*d); // ROOT thinks, m_last is still a matrix - Float_t last = m_last[0][0]/denom; - - TMatrix m1 = gT*Blo; - m1 *= first; - - TMatrix m2 = dT*Blo; - m2 *=last; - - TMatrix w_amp = m1 - m2; - - TMatrix m_first1 = gT*(Blo*g); - Float_t first1 = m_first1[0][0]/denom; - - TMatrix m_last1 = gT*(Blo*d); - Float_t last1 = m_last1 [0][0]/denom; - - TMatrix m11 = dT*Blo; - m11 *=first1; - - TMatrix m21 = gT*Blo; - m21 *=last1; - - TMatrix w_time= m11 - m21; - - for (Int_t count=0; count < fWindowSizeLoGain; count++) - { - const Int_t idx = i+fBinningResolutionLoGain/2+fBinningResolutionLoGain*count-1; - fAmpWeightsLoGain [idx] = w_amp [0][count]; - fTimeWeightsLoGain[idx] = w_time[0][count]; - } - - } // end loop over i - } - - ofstream fn(filename.Data()); - - fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl; - fn << "# (Amplitude) (Time) " << endl; - - for (Int_t i=0; i maxcont) - { - 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 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++; - } - } - - FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range, - sum, dsum, time, dtime, - sat, maxpos); -} - -// -------------------------------------------------------------------------- -// -// Calculates the arrival time and charge for each pixel -// -void MExtractTimeAndChargeSpline::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) -{ - 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; - } - - if (*p++ >= fSaturationLimit) - sat++; - } - - FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range, - sum, dsum, time, dtime, - sat, maxpos); -} -*/ // -------------------------------------------------------------------------- Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h (revision 8164) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h (revision 8165) @@ -42,8 +42,8 @@ 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 fRandomIter; //! Counter used to randomize weights for noise calculation Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); - Bool_t InitArrays(); + Bool_t InitArrays(Int_t n); public: @@ -94,9 +94,9 @@ 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); + Byte_t sat, Int_t maxpos) const; 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); + Byte_t sat, Int_t maxpos) const; ClassDef(MExtractTimeAndChargeSpline, 4) // Task to Extract Arrival Times and Charges using a Cubic Spline