Index: trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc (revision 2488) +++ trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc (revision 2502) @@ -165,42 +165,74 @@ // -------------------------------------------------------------------------- // -// Replaces each pixel by the average of its surrounding pixels. +// Replaces each pixel (signal, signal error, pedestal, pedestal rms) +// by the average of its surrounding pixels. // If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also // included. // +// NT: Informations about the interpolated pedestals are added. +// When the option Interpolated is called, the blind pixel with the new +// values of signal and fluttuation is included in the calculation of +// the Image Parameters. +// void MBlindPixelCalc::Interpolate() const { - const UShort_t entries = fEvt->GetNumPixels(); - - Double_t *nphot = new Double_t[entries]; - Double_t *perr = new Double_t[entries]; - Double_t *ped = new Double_t[entries]; - - // - // remove the pixels in fPixelsIdx if they are set to be used, - // (set them to 'unused' state) + const UShort_t entries = fGeomCam->GetNumPixels(); + + // + // Create arrays + // + Double_t *nphot = new Double_t[entries]; + Double_t *perr = new Double_t[entries]; + Double_t *ped = new Double_t[entries]; + Double_t *pedrms = new Double_t[entries]; + + // + // Loop over all pixels // for (UShort_t i=0; iIsBlind(idx)) + // + // Check whether pixel with idx i is blind + // + if (!fPixels->IsBlind(i)) continue; - const MGeomPix &gpix = (*fGeomCam)[idx]; - - Int_t num = TESTBIT(fFlags, kUseCentralPixel) ? 1 : 0; - - nphot[i] = TESTBIT(fFlags, kUseCentralPixel) ? pix.GetNumPhotons() : 0; - perr[i] = TESTBIT(fFlags, kUseCentralPixel) ? (*fPed)[idx].GetPedestalRms() : 0; - ped[i] = TESTBIT(fFlags, kUseCentralPixel) ? (*fPed)[idx].GetPedestal() : 0; - - nphot[i] *= fGeomCam->GetPixRatio(idx); - // FIXME? perr - // FIXME? ped - + // + // Get a pointer to this pixel. If it is not yet existing + // create a new entry for this pixel in MCerPhotEvt + // + const MCerPhotPix *pix = fEvt->GetPixById(i); + if (!pix) + pix = fEvt->AddPixel(i, 0, 0); + + // + // Get the corresponding geometry and pedestal + // + const MGeomPix &gpix = (*fGeomCam)[i]; + const MPedestalPix &ppix = (*fPed)[i]; + + // Do Not-Use-Central-Pixel + const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel); + + Int_t num = nucp ? 0 : 1; + + nphot[i] = nucp ? 0 : pix->GetNumPhotons(); + ped[i] = nucp ? 0 : ppix.GetPedestal(); + perr[i] = nucp ? 0 : Sqr(pix->GetErrorPhot()); + pedrms[i] = nucp ? 0 : Sqr(ppix.GetPedestalRms()); + + // + // The values are rescaled to the small pixels area for the right comparison + // + const Double_t ratio = fGeomCam->GetPixRatio(i); + + nphot[i] *= ratio; + ped[i] *= ratio; + perr[i] *= Sqr(ratio); + pedrms[i] *= Sqr(pedrms[i]); + + // + // Loop over all its neighbors + // const Int_t n = gpix.GetNumNeighbors(); for (int j=0; jIsBlind(nidx)) continue; + // + // Check whether the neighbor has a signal stored + // const MCerPhotPix *evtpix = fEvt->GetPixById(nidx); if (!evtpix) continue; - nphot[i] += evtpix->GetNumPhotons()*fGeomCam->GetPixRatio(nidx); - perr[i] += (*fPed)[nidx].GetPedestalRms(); - ped[i] += (*fPed)[nidx].GetPedestal(); - // FIXME? perr - // FIXME? ped + // + // Get the geometry for the neighbor + // + const Double_t nratio = fGeomCam->GetPixRatio(nidx); + MPedestalPix &nppix = (*fPed)[nidx]; + + // + //The error is calculated as the quadratic sum of the errors + // + ped[i] += nratio*nppix.GetPedestal(); + nphot[i] += nratio*evtpix->GetNumPhotons(); + perr[i] += Sqr(nratio*evtpix->GetErrorPhot()); + pedrms[i] += Sqr(nratio*nppix.GetPedestalRms()); num++; } - nphot[i] /= num*fGeomCam->GetPixRatio(idx); - perr[i] /= num/*FIXME:???*/; - ped[i] /= num/*FIXME:???*/; - } - - if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam) - for (UShort_t i=0; iIsBlind(pix.GetPixId())) - { - pix.Set(nphot[i], -1); - (*fPed)[pix.GetPixId()].Set(ped[i], perr[i]); - } - } - + // + // Now the mean is calculated and the values rescaled back to the pixel area + // + nphot[i] /= num*ratio; + ped[i] /= num*ratio; + perr[i] = TMath::Sqrt(perr[i]/num)*ratio; + pedrms[i] = TMath::Sqrt(pedrms[i]/num)*ratio; + + } + + // + // Now the new pixel values are calculated and can be replaced in + // the corresponding containers + // + for (UShort_t i=0; iIsBlind(i)) + continue; + + // + // It must exist, we have created it in the loop before. + // + fEvt->GetPixById(i)->Set(nphot[i], perr[i]); + (*fPed)[i].Set(ped[i], pedrms[i]); + } + + // + // Delete the intermediat arrays + // delete nphot; delete perr; delete ped; + delete pedrms; } Index: trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.h (revision 2488) +++ trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.h (revision 2502) @@ -35,4 +35,6 @@ }; + static Double_t Sqr(Double_t x) { return x*x; } + void Interpolate() const; void Unmap() const; Index: trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h (revision 2488) +++ trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h (revision 2502) @@ -32,5 +32,5 @@ //void InitSize(UInt_t num) { fPixels->Expand(num); } - void AddPixel(Int_t idx, Float_t nph, Float_t er) + MCerPhotPix *AddPixel(Int_t idx, Float_t nph=0, Float_t er=0) { // @@ -47,5 +47,5 @@ fLut[idx] = fNumPixels; - new ((*fPixels)[fNumPixels++]) MCerPhotPix(idx, nph, er); + return new ((*fPixels)[fNumPixels++]) MCerPhotPix(idx, nph, er); }