Index: trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (revision 4628) +++ trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (revision 4629) @@ -394,19 +394,29 @@ const Int_t npixels = fPedestals->GetSize(); + const Int_t areas = fPedestals->GetAverageAreas(); + const Int_t sectors = fPedestals->GetAverageSectors(); if (fSumx.GetSize()==0) - { + { fSumx. Set(npixels); fSumx2.Set(npixels); - + + fAreaSumx. Set(areas); + fAreaSumx2.Set(areas); + fAreaValid.Set(areas); + + fSectorSumx. Set(sectors); + fSectorSumx2.Set(sectors); + fSectorValid.Set(sectors); + fSumx.Reset(); fSumx2.Reset(); - } + } if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0) { - *fLog << err << GetDescriptor() + *fLog << err << GetDescriptor() << ": Number of extracted Slices is 0, abort ... " << endl; - return kFALSE; + return kFALSE; } @@ -429,66 +439,4 @@ Int_t MPedCalcPedRun::Process() { - MRawEvtPixelIter pixel(fRawEvt); - - while (pixel.Next()) - { - const UInt_t idx = pixel.GetPixelId(); - - Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst; - Byte_t *end = ptr + fWindowSizeHiGain; - - UInt_t sum = 0; - UInt_t sqr = 0; - - if (fWindowSizeHiGain != 0) - { - do - { - sum += *ptr; - sqr += *ptr * *ptr; - } - while (++ptr != end); - } - - if (fWindowSizeLoGain != 0) - { - ptr = pixel.GetLoGainSamples() + fLoGainFirst; - end = ptr + fWindowSizeLoGain; - - do - { - sum += *ptr; - sqr += *ptr * *ptr; - } - while (++ptr != end); - } - - const Float_t msum = (Float_t)sum; - - // - // These three lines have been uncommented by Markus Gaug - // If anybody needs them, please contact me!! - // - // const Float_t higainped = msum/fNumHiGainSlices; - // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.)); - // (*fPedestals)[idx].Set(higainped, higainrms); - - fSumx[idx] += msum; - // - // The old version: - // - // const Float_t msqr = (Float_t)sqr; - // fSumx2[idx] += msqr; - // - // The new version: - // - const Float_t sqrsum = msum*msum; - fSumx2[idx] += sqrsum; - } - - fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain; - - return kTRUE; -/* MRawEvtPixelIter pixel(fRawEvt); @@ -497,5 +445,5 @@ const UInt_t idx = pixel.GetPixelId(); const UInt_t aidx = (*fGeom)[idx].GetAidx(); - const UInt_t sector = (*fGeom)[idx].GetSector(); + const UInt_t sector = (*fGeom)[idx].GetSector(); Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst; @@ -542,5 +490,5 @@ fSumx[idx] += msum; fAreaSumx[aidx] += msum; - fSectorSumx[sector] += msum; + fSectorSumx[sector] += msum; // // The old version: @@ -554,10 +502,10 @@ fSumx2[idx] += sqrsum; fAreaSumx2[aidx] += sqrsum; - fSectorSumx2[sector] += sqrsum; + fSectorSumx2[sector] += sqrsum; } fPedestals->SetReadyToSave(); fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain; -*/ + return kTRUE; } @@ -569,82 +517,42 @@ Int_t MPedCalcPedRun::PostProcess() { - const ULong_t nevts = GetNumExecutions(); - if (nevts<1) - { - *fLog << err << "ERROR - Not enough events recorded...abort." << endl; - return kFALSE; - } - - *fLog << flush << inf << "Calculating pedestals..." << flush; - - // Compute pedestals and rms from the whole run - const ULong_t n = fNumSamplesTot; - const ULong_t npix = fGeom->GetNumPixels(); - - for (UInt_t pixidx=0; pixidxSetTotalEntries(fNumSamplesTot); - fPedestals->ReCalc(*fGeom); - fPedestals->SetReadyToSave(); - - return kTRUE; + // Compute pedestals and rms from the whole run + const ULong_t n = fNumSamplesTot; + const ULong_t nevts = GetNumExecutions(); + + MRawEvtPixelIter pixel(fRawEvt); + + while (pixel.Next()) + { + + const Int_t pixid = pixel.GetPixelId(); + const UInt_t aidx = (*fGeom)[pixid].GetAidx(); + const UInt_t sector = (*fGeom)[pixid].GetSector(); + + fAreaValid [aidx]++; + fSectorValid[sector]++; + + const Float_t sum = fSumx.At(pixid); + const Float_t sum2 = fSumx2.At(pixid); + const Float_t higainped = sum/n; + // + // The old version: + // + // const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.)); + // + // The new version: + // + // 1. Calculate the Variance of the sums: + Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.); + // 2. Scale the variance to the number of slices: + higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain); + // 3. Calculate the RMS from the Variance: + (*fPedestals)[pixid].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar)); + + } // // Loop over the (two) area indices to get the averaged pedestal per aidx // - /* - - // Compute pedestals and rms from the whole run - const ULong_t n = fNumSamplesTot; - const ULong_t nevts = GetNumExecutions(); - - MRawEvtPixelIter pixel(fRawEvt); - - while (pixel.Next()) - { - const Int_t pixid = pixel.GetPixelId(); - const UInt_t aidx = (*fGeom)[pixid].GetAidx(); - const UInt_t sector = (*fGeom)[pixid].GetSector(); - - fAreaValid [aidx]++; - fSectorValid[sector]++; - - const Float_t sum = fSumx.At(pixid); - const Float_t sum2 = fSumx2.At(pixid); - const Float_t higainped = sum/n; - // - // The old version: - // - // const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.)); - // - // The new version: - // - // 1. Calculate the Variance of the sums: - Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.); - // 2. Scale the variance to the number of slices: - higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain); - // 3. Calculate the RMS from the Variance: - (*fPedestals)[pixid].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar), 0, nevts); - } for (Int_t aidx=0; aidxGetAverageArea(aidx).Set(higainped, higainrms); } - */ + // // Loop over the (six) sector indices to get the averaged pedestal per sector // - /* for (Int_t sector=0; sectorGetAverageSector(sector).Set(higainped, higainrms); } + fPedestals->SetTotalEntries(fNumSamplesTot); fPedestals->SetReadyToSave(); return kTRUE; - */ } Index: trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.h =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.h (revision 4628) +++ trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.h (revision 4629) @@ -33,5 +33,4 @@ TArrayD fSumx; // sum of values TArrayD fSumx2; // sum of squared values - /* TArrayD fAreaSumx; // averaged sum of values per area idx TArrayD fAreaSumx2; // averaged sum of squared values per area idx @@ -40,5 +39,4 @@ TArrayD fSectorSumx2; // averaged sum of squared values per sector TArrayI fSectorValid; // number of valid pixel with sector idx - */ Int_t PreProcess (MParList *pList);