Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3844) +++ trunk/MagicSoft/Mars/Changelog (revision 3845) @@ -20,4 +20,9 @@ 2004/04/26: Markus Gaug + + * mcalib/MCalibrationQEPix.cc + - found a bug in the setting of validity flags of colours. This + made the QE used for calibration uncorrect by about a factor 1.5 + for the outer pixels if only one colour is used for calibration. * mjobs/MJCalibration.[h,cc] Index: trunk/MagicSoft/Mars/mcalib/MCalibrate.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrate.cc (revision 3844) +++ trunk/MagicSoft/Mars/mcalib/MCalibrate.cc (revision 3845) @@ -335,5 +335,5 @@ if (pix.IsFFactorMethodValid()) { - calibQE = qe.GetQECascadesFFactor ( zenith ); + calibQE = qe.GetQECascadesFFactor ( zenith ); calibQEVar = qe.GetQECascadesFFactorVar( zenith ); } Index: trunk/MagicSoft/Mars/mcalib/MCalibrationQEPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationQEPix.cc (revision 3844) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationQEPix.cc (revision 3845) @@ -204,6 +204,7 @@ const Float_t newavqe = GetQEBlindPixel(col) / GetDefaultQE (col) ; const Float_t newavqevar = ( GetQEBlindPixelRelVar(col) + GetDefaultQERelVar(col) ) * newavqe * newavqe; - wav += newavqe / newavqevar; - sumw += 1. / newavqevar ; + const Float_t weight = 1./newavqevar; + wav += newavqe * weight; + sumw += weight; } } @@ -224,6 +225,7 @@ const Float_t newavqe = GetQECombined(col) / GetDefaultQE (col) ; const Float_t newavqevar = ( GetQECombinedRelVar(col) + GetDefaultQERelVar(col) ) * newavqe * newavqe; - wav += newavqe / newavqevar; - sumw += 1./ newavqevar; + const Float_t weight = 1./newavqevar; + wav += newavqe * weight; + sumw += weight; } @@ -243,8 +245,9 @@ if (IsFFactorMethodValid (col)) { - const Float_t newavqe = GetQEFFactor(col) / GetDefaultQE (col) ; + const Float_t newavqe = GetQEFFactor(col) / GetDefaultQE (col) ; const Float_t newavqevar = ( GetQEFFactorRelVar(col) + GetDefaultQERelVar(col) ) * newavqe * newavqe; - wav += newavqe / newavqevar; - sumw += 1./ newavqevar; + const Float_t weight = 1./newavqevar; + wav += newavqe *weight; + sumw += weight; } @@ -268,6 +271,7 @@ const Float_t newavqe = GetQEPINDiode(col) / GetDefaultQE (col) ; const Float_t newavqevar = ( GetQEPINDiodeRelVar(col) + GetDefaultQERelVar(col) ) * newavqe * newavqe; - wav += newavqe / newavqevar; - sumw += 1./ newavqevar; + const Float_t weight = 1./newavqevar; + wav += newavqe *weight; + sumw += weight; } } @@ -932,5 +936,5 @@ Bool_t MCalibrationQEPix::IsBlindPixelMethodValid (MCalibrationCam::PulserColor_t col) const { - return TESTBIT(fValidFlags[ MCalibrationCam::kGREEN ],kBlindPixelMethodValid); + return TESTBIT(fValidFlags[ col ],kBlindPixelMethodValid); } @@ -942,5 +946,5 @@ Bool_t MCalibrationQEPix::IsCombinedMethodValid (MCalibrationCam::PulserColor_t col) const { - return TESTBIT(fValidFlags[ MCalibrationCam::kGREEN ],kCombinedMethodValid); + return TESTBIT(fValidFlags[ col ],kCombinedMethodValid); } @@ -951,5 +955,5 @@ Bool_t MCalibrationQEPix::IsFFactorMethodValid (MCalibrationCam::PulserColor_t col) const { - return TESTBIT(fValidFlags[ MCalibrationCam::kGREEN ],kFFactorMethodValid); + return TESTBIT(fValidFlags[ col ],kFFactorMethodValid); } @@ -960,5 +964,5 @@ Bool_t MCalibrationQEPix::IsPINDiodeMethodValid (MCalibrationCam::PulserColor_t col) const { - return TESTBIT(fValidFlags[ MCalibrationCam::kGREEN ],kPINDiodeMethodValid); + return TESTBIT(fValidFlags[ col ],kPINDiodeMethodValid); } @@ -1018,7 +1022,7 @@ { if (b) - SETBIT(fValidFlags[ MCalibrationCam::kGREEN ],kBlindPixelMethodValid); + SETBIT(fValidFlags[ col ],kBlindPixelMethodValid); else - CLRBIT(fValidFlags[ MCalibrationCam::kGREEN ],kBlindPixelMethodValid); + CLRBIT(fValidFlags[ col ],kBlindPixelMethodValid); } @@ -1030,7 +1034,7 @@ { if (b) - SETBIT(fValidFlags[ MCalibrationCam::kGREEN ],kCombinedMethodValid); + SETBIT(fValidFlags[ col ],kCombinedMethodValid); else - CLRBIT(fValidFlags[ MCalibrationCam::kGREEN ],kCombinedMethodValid); + CLRBIT(fValidFlags[ col ],kCombinedMethodValid); } @@ -1042,7 +1046,7 @@ { if (b) - SETBIT(fValidFlags[ MCalibrationCam::kGREEN ],kFFactorMethodValid); + SETBIT(fValidFlags[ col ],kFFactorMethodValid); else - CLRBIT(fValidFlags[ MCalibrationCam::kGREEN ],kFFactorMethodValid); + CLRBIT(fValidFlags[ col ],kFFactorMethodValid); } @@ -1054,7 +1058,7 @@ { if (b) - SETBIT(fValidFlags[ MCalibrationCam::kGREEN ],kPINDiodeMethodValid); + SETBIT(fValidFlags[ col ],kPINDiodeMethodValid); else - CLRBIT(fValidFlags[ MCalibrationCam::kGREEN ],kPINDiodeMethodValid); + CLRBIT(fValidFlags[ col ],kPINDiodeMethodValid); } @@ -1128,5 +1132,5 @@ AddAverageFFactorQEs(MCalibrationCam::kUV , weightedav, sumweights); AddAverageFFactorQEs(MCalibrationCam::kCT1 , weightedav, sumweights); - + if (weightedav == 0. || sumweights == 0.) return kFALSE; Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h (revision 3844) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h (revision 3845) @@ -72,5 +72,5 @@ virtual void FinalizeBadPixels(); - void CalcAverageSigma(); + virtual void CalcAverageSigma(); void DrawAverageSigma(Bool_t sat, Bool_t inner, Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.cc (revision 3845) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.cc (revision 3845) @@ -0,0 +1,526 @@ +/* ======================================================================== *\ +! +! * +! * 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): Markus Gaug 02/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationTestCam +// +// Fills the calibrated signal from an MCerPhotEvt into +// MHCalibrationTestPix for every: +// +// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray +// or MHCalibrationCam::fHiGainArray, respectively. +// +// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), +// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and +// MHCalibrationCam::fAverageHiGainAreas +// +// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), +// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors +// and MHCalibrationCam::fAverageHiGainSectors +// +// The signals are filled into a histogram and an array, in order to perform +// a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an +// event-by-event basis and written into the corresponding average pixels. +// +// The histograms are fitted to a Gaussian, mean and sigma with its errors +// and the fit probability are extracted. If none of these values are NaN's and +// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%), +// the fit is declared valid. +// Otherwise, the fit is repeated within ranges of the previous mean +// +- MHGausEvents::fPickupLimit (default: 5) sigma (see MHGausEvents::RepeatFit()) +// In case this does not make the fit valid, the histogram means and RMS's are +// taken directly (see MHGausEvents::BypassFit()) and the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Outliers of more than MHGausEvents::fPickupLimit (default: 5) sigmas +// from the mean are counted as Pickup events (stored in MHGausEvents::fPickup) +// +// The class also fills arrays with the signal vs. event number, creates a fourier +// spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the +// projected fourier components follow an exponential distribution. +// In case that the probability of the exponential fit is less than +// MHGausEvents::fProbLimit (default: 0.5%), the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// This same procedure is performed for the average pixels. +// +// The following results are written into an MCalibrationCam: +// +// - MCalibrationPix::SetMean() +// - MCalibrationPix::SetMeanErr() +// - MCalibrationPix::SetSigma() +// - MCalibrationPix::SetSigmaErr() +// - MCalibrationPix::SetProb() +// - MCalibrationPix::SetNumPickup() +// +// For all averaged areas, the fitted sigma is multiplied with the square root of +// the number involved pixels in order to be able to compare it to the average of +// sigmas in the camera. +// +///////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationTestCam.h" +#include "MHCalibrationTestPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MCalibrationCam.h" +#include "MCalibrationPix.h" + +#include "MCerPhotEvt.h" +#include "MCerPhotPix.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" + +ClassImp(MHCalibrationTestCam); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHCalibrationTestCam::MHCalibrationTestCam(const char *name, const char *title) +{ + + fName = name ? name : "MHCalibrationTestCam"; + fTitle = title ? title : "Histogram class for testing the calibration"; + +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// - MCalibrationTestCam +// +// Searches pointer to: +// - MCerPhotEvt +// +// Initializes, if empty to MGeomCam::GetNumAreas() for: +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// +// Initializes, if empty to MGeomCam::GetNumSectors() for: +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +// Calls MHCalibrationCam::InitHists() for every entry in: +// - MHCalibrationCam::fHiGainArray +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +// Sets Titles and Names for the Histograms +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +// Sets number of bins to MHCalibrationCam::fAverageNbins for: +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +Bool_t MHCalibrationTestCam::ReInitHists(MParList *pList) +{ + + MCerPhotEvt *signal = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!signal) + { + *fLog << err << "MCerPhotEvt not found... abort." << endl; + return kFALSE; + } + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nsectors = fGeom->GetNumSectors(); + const Int_t nareas = fGeom->GetNumAreas(); + + if (fHiGainArray->GetEntries()==0) + { + fHiGainArray->Expand(npixels); + for (Int_t i=0; iGetEntries()==0) + { + fLoGainArray->Expand(npixels); + for (Int_t i=0; iGetEntries()==0) + { + fAverageHiGainAreas->Expand(nareas); + + for (Int_t j=0; jSetTitle("Test Calibrations Area Idx "); + GetAverageHiGainArea(j).SetNbins(fAverageNbins); + GetAverageHiGainArea(j).InitBins(); + GetAverageHiGainArea(j).ChangeHistId(j); + GetAverageHiGainArea(j).SetEventFrequency(fPulserFrequency); + + } + } + + if (fAverageLoGainAreas->GetEntries()==0) + { + fAverageLoGainAreas->Expand(nareas); + + for (Int_t j=0; jSetTitle("Test Calibrations Area Idx "); + GetAverageLoGainArea(j).SetNbins(fAverageNbins); + GetAverageLoGainArea(j).InitBins(); + GetAverageLoGainArea(j).ChangeHistId(j); + GetAverageLoGainArea(j).SetEventFrequency(fPulserFrequency); + + } + } + + + if (fAverageHiGainSectors->GetEntries()==0) + { + fAverageHiGainSectors->Expand(nsectors); + + for (Int_t j=0; jSetTitle("Test Calibrations Sector "); + GetAverageHiGainSector(j).SetNbins(fAverageNbins); + GetAverageHiGainSector(j).InitBins(); + GetAverageHiGainSector(j).ChangeHistId(j); + GetAverageHiGainSector(j).SetEventFrequency(fPulserFrequency); + } + } + + + if (fAverageLoGainSectors->GetEntries()==0) + { + fAverageLoGainSectors->Expand(nsectors); + + for (Int_t j=0; jSetTitle("Test Calibrations Sector "); + GetAverageLoGainSector(j).SetNbins(fAverageNbins); + GetAverageLoGainSector(j).InitBins(); + GetAverageLoGainSector(j).ChangeHistId(j); + GetAverageLoGainSector(j).SetEventFrequency(fPulserFrequency); + } + } + + + return kTRUE; +} + + +// ------------------------------------------------------------------------------- +// +// Retrieves pointer to MCerPhotEvt: +// +// Retrieves from MGeomCam: +// - number of pixels +// - number of pixel areas +// - number of sectors +// +// Fills HiGain or LoGain histograms (MHGausEvents::FillHistAndArray()), respectively +// depending on MCerPhotPix::IsLoGainUsed(), with: +// - MCerPhotPix::GetArrivalTime(pixid) - MCerPhotPix::GetArrivalTime(1); +// (i.e. the time difference between pixel i and pixel 1 (hardware number: 2) ) +// +Bool_t MHCalibrationTestCam::FillHists(const MParContainer *par, const Stat_t w) +{ + + MCerPhotEvt *calibration = (MCerPhotEvt*)par; + if (!calibration) + { + gLog << err << "No argument in MHCalibrationRelTimeCam::Fill... abort." << endl; + return kFALSE; + } + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nareas = fGeom->GetNumAreas(); + const Int_t nsectors = fGeom->GetNumSectors(); + + Float_t sumareahi [nareas]; + Float_t sumsectorhi[nsectors]; + Int_t numareahi [nareas]; + Int_t numsectorhi[nsectors]; + + memset(sumareahi, 0, nareas * sizeof(Float_t)); + memset(sumsectorhi, 0, nsectors*sizeof(Float_t)); + + for (Int_t i=0; iGetPixById(i); + if (!pix) + continue; + + const Float_t signal = pix->GetNumPhotons(); + + const Int_t aidx = (*fGeom)[i].GetAidx(); + const Int_t sector = (*fGeom)[i].GetSector(); + + histhi.FillHistAndArray(signal) ; + sumareahi [aidx] += signal; + numareahi [aidx] ++; + sumsectorhi[sector] += signal; + numsectorhi[sector] ++; + } + + for (Int_t j=0; jGetSize(); i++) + { + + MHGausEvents &hist = (*this)[i]; + + if (hist.IsExcluded()) + continue; + + if (hist.IsEmpty()) + continue; + + if (!hist.FitGaus()) + if (!hist.RepeatFit()) + { + hist.BypassFit(); + } + + hist.CreateFourierSpectrum(); + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHGausEvents &hist = GetAverageHiGainArea(j); + if (hist.IsEmpty()) + continue; + + if (!hist.FitGaus()) + if (!hist.RepeatFit()) + { + hist.BypassFit(); + } + + hist.CreateFourierSpectrum(); + + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHGausEvents &hist = GetAverageHiGainSector(j); + if (hist.IsEmpty()) + continue; + + if (!hist.FitGaus()) + if (!hist.RepeatFit()) + { + hist.BypassFit(); + } + hist.CreateFourierSpectrum(); + + + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set: +// - MBadPixelsPix::kTestNotFitted +// - MBadPixelsPix::kTestOscillating +// +void MHCalibrationTestCam::FinalizeBadPixels() +{ + +} + +// -------------------------------------------------------------------------- +// +// The types are as follows: +// +// Fitted values: +// ============== +// +// 0: Fitted Mean Test Calibration (MHGausEvents::GetMean()) +// 1: Error Mean Test Calibration (MHGausEvents::GetMeanErr()) +// 2: Sigma fitted Test Calibration (MHGausEvents::GetSigma()) +// 3: Error Sigma Test Calibration (MHGausEvents::GetSigmaErr()) +// +// Useful variables derived from the fit results: +// ============================================= +// +// 4: Returned probability of Gauss fit (calls: MHGausEvents::GetProb()) +// +// Localized defects: +// ================== +// +// 5: Gaus fit not OK (calls: MHGausEvents::IsGausFitOK()) +// 6: Fourier spectrum not OK (calls: MHGausEvents::IsFourierSpectrumOK()) +// +Bool_t MHCalibrationTestCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + + if (fHiGainArray->GetSize() <= idx) + return kFALSE; + + const MHGausEvents &pix = (*this)[idx]; + + switch (type) + { + case 0: + val = pix.GetMean(); + break; + case 1: + val = pix.GetMeanErr(); + break; + case 2: + val = pix.GetSigma(); + break; + case 3: + val = pix.GetSigmaErr(); + break; + case 4: + val = pix.GetProb(); + break; + case 5: + if (!pix.IsGausFitOK()) + val = 1.; + break; + case 6: + if (!pix.IsFourierSpectrumOK()) + val = 1.; + break; + default: + return kFALSE; + } + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calls MHGausEvents::DrawClone() for pixel idx +// +void MHCalibrationTestCam::DrawPixelContent(Int_t idx) const +{ + (*this)[idx].DrawClone(); +} + + +//------------------------------------------------------------ +// +// For all averaged areas, the fitted sigma is multiplied with the square root of +// the number involved pixels +// +void MHCalibrationTestCam::CalcAverageSigma() +{ + + for (UInt_t j=0; jGetNumAreas(); j++) + { + + MHGausEvents &hist = GetAverageHiGainArea(j); + + const Float_t numsqr = TMath::Sqrt((Float_t)fAverageAreaNum[j]); + fAverageAreaSigma[j] = hist.GetSigma () * numsqr; + fAverageAreaSigmaVar[j] = hist.GetSigmaErr () * hist.GetSigmaErr() * numsqr; + + fAverageAreaRelSigma [j] = fAverageAreaSigma[j] / hist.GetMean(); + fAverageAreaRelSigmaVar[j] = fAverageAreaSigmaVar[j] / (fAverageAreaSigma[j]*fAverageAreaSigma[j]); + fAverageAreaRelSigmaVar[j] += hist.GetMeanErr()*hist.GetMeanErr()/hist.GetMean()/hist.GetMean(); + fAverageAreaRelSigmaVar[j] *= fAverageAreaRelSigma[j]; + } +} Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.h (revision 3845) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.h (revision 3845) @@ -0,0 +1,33 @@ +#ifndef MARS_MHCalibrationTestCam +#define MARS_MHCalibrationTestCam + +#ifndef MARS_MHCalibrationCam +#include "MHCalibrationCam.h" +#endif + +class MGeomCam; +class MHCalibrationTestCam : public MHCalibrationCam +{ + +private: + + Bool_t ReInitHists(MParList *pList); + Bool_t FillHists(const MParContainer *par, const Stat_t w=1); + Bool_t FinalizeHists(); + void FinalizeBadPixels(); + +public: + + MHCalibrationTestCam(const char *name=NULL, const char *title=NULL); + ~MHCalibrationTestCam() {} + + Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; + void DrawPixelContent(Int_t idx) const; + + void CalcAverageSigma(); + + ClassDef(MHCalibrationTestCam, 0) // Histogram class for Relative Time Camera Calibration +}; + +#endif +