Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc (revision 2818) +++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc (revision 2821) @@ -18,6 +18,8 @@ ! Author(s): Josep Flix 04/2001 ! Author(s): Thomas Bretz 05/2001 +! Author(s): Sebastian Commichau 12/2003 +! Author(s): Javier Rico 01/2004 ! -! Copyright: MAGIC Software Development, 2000-2001 +! Copyright: MAGIC Software Development, 2000-2004 ! ! @@ -28,4 +30,10 @@ // MPedCalcPedRun // // // +// This task takes a pedestal run file and fills MPedestalCam during // +// the Process() with the pedestal and rms computed in an event basis. // +// In the PostProcess() MPedestalCam is finally filled with the pedestal // +// mean and rms computed in a run basis. // +// More than one run (file) can be merged // +// // // Input Containers: // // MRawEvtData // @@ -33,4 +41,5 @@ // Output Containers: // // MPedestalCam // +// // // // ///////////////////////////////////////////////////////////////////////////// @@ -70,5 +79,5 @@ if (!fRawEvt) { - *fLog << dbginf << "MRawEvtData not found... aborting." << endl; + *fLog << err << "MRawEvtData not found... aborting." << endl; return kFALSE; } @@ -78,69 +87,75 @@ return kFALSE; - MGeomCamMagic magiccam; - - fSumx.Set(magiccam.GetNumPixels()); - fSumx2.Set(magiccam.GetNumPixels()); - - for(UInt_t i=0;iFindObject("MRawRunHeader"); + if (!runheader) + { + *fLog << warn << dbginf; + *fLog << "Warning - cannot check file type, MRawRunHeader not found." << endl; + } + else + if (runheader->GetRunType() == kRTMonteCarlo) + return kTRUE; - fRunheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!fRunheader) - { - *fLog << warn << dbginf << - "Warning - cannot check file type, MRawRunHeader not found." << endl; - } - else - if (fRunheader->GetRunType() == kRTMonteCarlo) - { - return kTRUE; - } + fNumPixels = fPedestals->GetSize(); - fNumHiGainSamples = fRunheader->GetNumSamplesHiGain(); + if(fSumx.GetSize()==0) + { + fSumx.Set(fNumPixels); + fSumx2.Set(fNumPixels); - fPedestals->InitSize(fRunheader->GetNumPixel()); + memset(fSumx.GetArray(), 0, sizeof(Float_t)*fNumPixels); + memset(fSumx2.GetArray(), 0, sizeof(Float_t)*fNumPixels); + } + + // Calculate an even number for the hi gain samples to avoid + // biases due to the fluctuation in pedestal from one slice to + // the other one + fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1; return kTRUE; - } - Int_t MPedCalcPedRun::Process() { - MRawEvtPixelIter pixel(fRawEvt); while (pixel.Next()) { - Byte_t shift=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? 0:1; Byte_t *ptr = pixel.GetHiGainSamples(); - const Byte_t *end = ptr + fRawEvt->GetNumHiGainSamples()-shift; + const Byte_t *end = ptr + fNumHiGainSamples; - const Float_t higainped = CalcHiGainMean(ptr, end); - const Float_t higainrms = CalcHiGainRms(ptr, end, higainped); - - const UInt_t pixid = pixel.GetPixelId(); - MPedestalPix &pix = (*fPedestals)[pixid]; + UInt_t sum = 0; + UInt_t sqr = 0; - // cumulate the sum of pedestals and of pedestal squares - fSumx.AddAt(higainped+fSumx.At(pixid),pixid); - fSumx2.AddAt(GetSumx2(ptr, end)+fSumx2.At(pixid),pixid); - - // set the value of the pedestal and rms computed from the processed event - pix.Set(higainped, higainrms); + do + { + sum += *ptr; + sqr += *ptr * *ptr; + } + while (++ptr != end); + + const Float_t msum = (Float_t)sum; + const Float_t msqr = (Float_t)sqr; + + const Float_t higainped = msum/fNumHiGainSamples; + const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.)); + + const UInt_t idx = pixel.GetPixelId(); + (*fPedestals)[idx].Set(higainped, higainrms); + + fSumx[idx] += msum; + fSumx2[idx] += sqr; } - + fPedestals->SetReadyToSave(); - + fNumSamplesTot+=fNumHiGainSamples; + return kTRUE; } @@ -148,84 +163,24 @@ Int_t MPedCalcPedRun::PostProcess() { -// Compute pedestals and rms from the whole run - - MRawEvtPixelIter pixel(fRawEvt); - - while (pixel.Next()) + // Compute pedestals and rms from the whole run + const Int_t n = fNumSamplesTot; + + MRawEvtPixelIter pixel(fRawEvt); + + while (pixel.Next()) { - const UInt_t pixid = pixel.GetPixelId(); - MPedestalPix &pix = (*fPedestals)[pixid]; - - const Int_t N = GetNumExecutions(); - const Float_t sum = fSumx.At(pixid); - const Float_t sum2 = fSumx2.At(pixid); - const Float_t higainped = sum/N; - const Float_t higainrms = sqrt(1./(N-1.)*(sum2-sum*sum/N)); - pix.Set(higainped,higainrms); - + const UInt_t pixid = pixel.GetPixelId(); + + const Float_t sum = fSumx.At(pixid); + const Float_t sum2 = fSumx2.At(pixid); + + const Float_t higainped = sum/n; + const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.)); + + (*fPedestals)[pixid].Set(higainped, higainrms); } - return kTRUE; - + + return kTRUE; } -Float_t MPedCalcPedRun::CalcHiGainMean(Byte_t *ptr, const Byte_t *end) const -{ - Int_t sum=0; - Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) - ? fNumHiGainSamples : fNumHiGainSamples-1; - - do sum += *ptr; - while (++ptr != end); - - return (Float_t)sum/EvenNumSamples; -} - - - -Float_t MPedCalcPedRun::GetSumx2(Byte_t *ptr, const Byte_t *end) const -{ - Float_t square = 0; - - // Take an even number of time slices to avoid biases due to A/B effect - Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1; - - do - { - const Float_t val = (Float_t)(*ptr); - - square += val*val; - } while (++ptr != end); - - return square/EvenNumSamples; -} - - -Float_t MPedCalcPedRun::CalcHiGainRms(Byte_t *ptr, const Byte_t *end, Float_t higainped) const -{ - Float_t rms = 0; - Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1; - - do - { - const Float_t diff = (Float_t)(*ptr)-higainped; - - rms += diff*diff; - } while (++ptr != end); - - return sqrt(rms/(EvenNumSamples-1)); -} - - - -/* -Float_t MPedCalcPedRun::CalcHiGainMeanErr(Float_t higainrms) const -{ - return higainrms/sqrt((float)fNumHiGainSamples); -} - -Float_t MPedCalcPedRun::CalcHiGainRmsErr(Float_t higainrms) const -{ - return higainrms/sqrt(2.*fNumHiGainSamples); -} -*/ Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h (revision 2818) +++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h (revision 2821) @@ -16,5 +16,4 @@ #include -class MRawRunHeader; class MRawEvtData; class MPedestalCam; @@ -22,7 +21,8 @@ class MPedCalcPedRun : public MTask { - Byte_t fNumHiGainSamples; + Byte_t fNumHiGainSamples; + UShort_t fNumPixels; + ULong_t fNumSamplesTot; - MRawRunHeader *fRunheader; // raw event run header MRawEvtData *fRawEvt; // raw event data (time slices) MPedestalCam *fPedestals; // Pedestals of all pixels in the camera @@ -30,10 +30,4 @@ TArrayF fSumx; // sum of values TArrayF fSumx2; // sum of squared values - - Float_t CalcHiGainMean(Byte_t *ptr, const Byte_t *end) const; - Float_t CalcHiGainRms(Byte_t *ptr, const Byte_t *end, Float_t higainped) const; - Float_t GetSumx2(Byte_t* ptr, const Byte_t* end) const; - //Float_t CalcHiGainMeanErr(Float_t higainrms) const; - //Float_t CalcHiGainRmsErr(Float_t higainrms) const; Bool_t ReInit(MParList *pList); @@ -44,5 +38,4 @@ public: - MPedCalcPedRun(const char *name=NULL, const char *title=NULL);