Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationCam.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationCam.cc (revision 4929) @@ -0,0 +1,1111 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationCam +// +// Base class for camera calibration classes. Incorporates the TObjArray's: +// - fHiGainArray (for calibrated High Gains per pixel) +// - fLoGainArray (for calibrated Low Gains per pixel) +// - fAverageHiGainAreas (for averaged High Gains events per camera area index) +// - fAverageLoGainAreas (for averaged High Gains events per camera area index) +// - fAverageHiGainSectors (for averaged High Gains events per camera sector ) +// - fAverageLoGainSectors (for averaged High Gains events per camera sector ) +// These TObjArray's are called by their default constructors, thus no objects +// are created, until the derived class does so. +// +// The corresponding operators: [],() and the operators GetAverageHiGainArea(), +// GetAverageLoGainArea(), GetAverageHiGainSector() and GetAverageLoGainSector() +// have to be cast to the corresponding class. It is assumed that all classes +// dealing with calibration pixels derive from MHCalibrationPix. +// +// The following flag can be set: +// - SetDebug() for debug output. +// - SetLoGain() for the case that low-gain slices are available, but +// MRawRunHeader::GetNumLoGainSlices() gives still 0. +// +// The flag fLoGain steers if low-gain signal is treated at all or not. +// +///////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationCam.h" + +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MCalibrationPix.h" +#include "MCalibrationCam.h" + +#include "MHCalibrationPix.h" + +#include "MBadPixelsPix.h" +#include "MBadPixelsCam.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MParList.h" + +#include "MRawRunHeader.h" + +ClassImp(MHCalibrationCam); + +using namespace std; + +const Int_t MHCalibrationCam::fgAverageNbins = 2000; +const Int_t MHCalibrationCam::fgPulserFrequency = 500; + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - all pointers to NULL +// +// Initializes and sets owner of: +// - fHiGainArray, fLoGainArray +// - fAverageHiGainAreas, fAverageLoGainAreas +// - fAverageHiGainSectors, fAverageLoGainSectors +// +// Initializes: +// - fPulserFrequency to fgPulserFrequency +// +MHCalibrationCam::MHCalibrationCam(const char *name, const char *title) + : fColor(MCalibrationCam::kNONE), + fBadPixels(NULL), fIntensCam(NULL), fCam(NULL), fGeom(NULL), fRunHeader(NULL) +{ + + fHiGainArray = new TObjArray; + fHiGainArray->SetOwner(); + + fLoGainArray = new TObjArray; + fLoGainArray->SetOwner(); + + fAverageHiGainAreas = new TObjArray; + fAverageHiGainAreas->SetOwner(); + + fAverageLoGainAreas = new TObjArray; + fAverageLoGainAreas->SetOwner(); + + fAverageHiGainSectors = new TObjArray; + fAverageHiGainSectors->SetOwner(); + + fAverageLoGainSectors = new TObjArray; + fAverageLoGainSectors->SetOwner(); + + SetAverageNbins(); + SetPulserFrequency(); + + SetDebug (kFALSE); + SetLoGain (kTRUE); + SetOscillations(kTRUE); +} + +const Int_t MHCalibrationCam::GetSize() const +{ + return fHiGainArray->GetSize(); +} + +// -------------------------------------------------------------------------- +// +// Deletes the TClonesArray of: +// - fHiGainArray, fLoGainArray +// - fAverageHiGainAreas, fAverageLoGainAreas +// - fAverageHiGainSectors, fAverageLoGainSectors +// +MHCalibrationCam::~MHCalibrationCam() +{ + + delete fHiGainArray; + delete fLoGainArray; + + delete fAverageHiGainAreas; + delete fAverageLoGainAreas; + + delete fAverageHiGainSectors; + delete fAverageLoGainSectors; +} + +void MHCalibrationCam::ResetHists() +{ + + if (fHiGainArray) + { fHiGainArray->ForEach(MHCalibrationPix,Reset)(); } + if (fLoGainArray) + { fLoGainArray->ForEach(MHCalibrationPix,Reset)(); } + + if (fAverageHiGainAreas) + { fAverageHiGainAreas->ForEach(MHCalibrationPix,Reset)(); } + if (fAverageLoGainAreas) + { fAverageLoGainAreas->ForEach(MHCalibrationPix,Reset)(); } + if (fAverageHiGainSectors) + { fAverageHiGainSectors->ForEach(MHCalibrationPix,Reset)(); } + if (fAverageLoGainSectors) + { fAverageLoGainSectors->ForEach(MHCalibrationPix,Reset)(); } + +} + + + +// -------------------------------------------------------------------------- +// +// Get i-th High Gain pixel (pixel number) +// +MHCalibrationPix &MHCalibrationCam::operator[](UInt_t i) +{ + return *static_cast(fHiGainArray->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th High Gain pixel (pixel number) +// +const MHCalibrationPix &MHCalibrationCam::operator[](UInt_t i) const +{ + return *static_cast(fHiGainArray->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th Low Gain pixel (pixel number) +// +MHCalibrationPix &MHCalibrationCam::operator()(UInt_t i) +{ + return *static_cast(fLoGainArray->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th Low Gain pixel (pixel number) +// +const MHCalibrationPix &MHCalibrationCam::operator()(UInt_t i) const +{ + return *static_cast(fLoGainArray->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Returns the current size of the TObjArray fAverageHiGainAreas +// independently if the MHCalibrationPix is filled with values or not. +// +const Int_t MHCalibrationCam::GetAverageAreas() const +{ + return fAverageHiGainAreas->GetEntries(); +} + +// -------------------------------------------------------------------------- +// +// Get i-th High Gain pixel Area (area number) +// +MHCalibrationPix &MHCalibrationCam::GetAverageHiGainArea(UInt_t i) +{ + return *static_cast(fAverageHiGainAreas->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th High Gain pixel Area (area number) +// +const MHCalibrationPix &MHCalibrationCam::GetAverageHiGainArea(UInt_t i) const +{ + return *static_cast(fAverageHiGainAreas->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th Low Gain pixel Area (area number) +// +MHCalibrationPix &MHCalibrationCam::GetAverageLoGainArea(UInt_t i) +{ + return *static_cast(fAverageLoGainAreas->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th Low Gain pixel Area (area number) +// +const MHCalibrationPix &MHCalibrationCam::GetAverageLoGainArea(UInt_t i) const +{ + return *static_cast(fAverageLoGainAreas->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Returns the current size of the TObjArray fAverageHiGainSectors +// independently if the MHCalibrationPix is filled with values or not. +// +const Int_t MHCalibrationCam::GetAverageSectors() const +{ + return fAverageHiGainSectors->GetEntries(); +} + +// -------------------------------------------------------------------------- +// +// Get i-th High Gain Sector (sector number) +// +MHCalibrationPix &MHCalibrationCam::GetAverageHiGainSector(UInt_t i) +{ + return *static_cast(fAverageHiGainSectors->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th High Gain Sector (sector number) +// +const MHCalibrationPix &MHCalibrationCam::GetAverageHiGainSector(UInt_t i) const +{ + return *static_cast(fAverageHiGainSectors->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th Low Gain Sector (sector number) +// +MHCalibrationPix &MHCalibrationCam::GetAverageLoGainSector(UInt_t i) +{ + return *static_cast(fAverageLoGainSectors->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th Low Gain Sector (sector number) +// +const MHCalibrationPix &MHCalibrationCam::GetAverageLoGainSector(UInt_t i) const +{ + return *static_cast(fAverageLoGainSectors->UncheckedAt(i)); +} + + +const TArrayI &MHCalibrationCam::GetRunNumbers() const +{ + return fRunNumbers; +} + +// -------------------------------------------------------------------------- +// +// Our own clone function is necessary since root 3.01/06 or Mars 0.4 +// I don't know the reason. +// +// Creates new MHCalibrationCam +// Deletes the TObjArray's and Clones them individually +// Copies the TArray's +// Copies the fPulserFrequency +// +TObject *MHCalibrationCam::Clone(const char *) const +{ + + // const Int_t nhi = fHiGainArray->GetEntries(); + // const Int_t nlo = fLoGainArray->GetEntries(); + const Int_t navhi = fAverageHiGainAreas->GetEntries(); + const Int_t navlo = fAverageLoGainAreas->GetEntries(); + const Int_t nsehi = fAverageHiGainSectors->GetEntries(); + const Int_t nselo = fAverageLoGainSectors->GetEntries(); + + // + // FIXME, this might be done faster and more elegant, by direct copy. + // + MHCalibrationCam *cam = new MHCalibrationCam(); + + // cam->fHiGainArray->Expand(nhi); + // cam->fLoGainArray->Expand(nlo); + cam->fAverageHiGainAreas->Expand(navhi); + cam->fAverageLoGainAreas->Expand(navlo); + cam->fAverageHiGainSectors->Expand(nsehi); + cam->fAverageLoGainSectors->Expand(nselo); + + /* + for (int i=0; ifHiGainArray)[i]; + (*cam->fHiGainArray)[i] = (*fHiGainArray)[i]->Clone(); + } + for (int i=0; ifLoGainArray)[i]; + (*cam->fLoGainArray)[i] = (*fLoGainArray)[i]->Clone(); + } + */ + + for (int i=0; ifAverageHiGainAreas)[i]; + (*cam->fAverageHiGainAreas)[i] = (*fAverageHiGainAreas)[i]->Clone(); + } + for (int i=0; ifAverageLoGainAreas)[i]; + (*cam->fAverageLoGainAreas)[i] = (*fAverageLoGainAreas)[i]->Clone(); + } + for (int i=0; ifAverageHiGainSectors)[i]; + (*cam->fAverageHiGainSectors)[i] = (*fAverageHiGainSectors)[i]->Clone(); + } + for (int i=0; ifAverageLoGainSectors)[i]; + (*cam->fAverageLoGainSectors)[i] = (*fAverageLoGainSectors)[i]->Clone(); + } + + cam->fAverageAreaNum = fAverageAreaNum; + cam->fAverageAreaSat = fAverageAreaSat; + cam->fAverageAreaSigma = fAverageAreaSigma; + cam->fAverageAreaSigmaVar = fAverageAreaSigmaVar; + cam->fAverageAreaRelSigma = fAverageAreaRelSigma; + cam->fAverageAreaRelSigmaVar = fAverageAreaRelSigmaVar; + cam->fAverageSectorNum = fAverageSectorNum; + cam->fRunNumbers = fRunNumbers; + + cam->fPulserFrequency = fPulserFrequency; + cam->fAverageNbins = fAverageNbins; + + return cam; +} + +// -------------------------------------------------------------------------- +// +// Gets the pointers to: +// - MGeomCam +// +// Calls SetupHists(const MParList *pList) +// +// Calls Delete-Function of: +// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +Bool_t MHCalibrationCam::SetupFill(const MParList *pList) +{ + + fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fGeom) + { + *fLog << err << GetDescriptor() + << ": MGeomCam not found... aborting." << endl; + return kFALSE; + } + + fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << warn << GetDescriptor() + << ": MRawRunHeader not found... will not store run numbers." << endl; + } + + fHiGainArray->Delete(); + fLoGainArray->Delete(); + + fAverageHiGainAreas->Delete(); + fAverageLoGainAreas->Delete(); + + fAverageHiGainSectors->Delete(); + fAverageLoGainSectors->Delete(); + + return SetupHists(pList); +} + + +Bool_t MHCalibrationCam::SetupHists(const MParList *pList) +{ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// - MBadPixelsCam +// +// Searches pointer to: +// - MArrivalTimeCam +// +// Initializes, if empty to MArrivalTimeCam::GetSize() for: +// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray +// +// Initializes, if empty to MGeomCam::GetNumAreas() for: +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// +// Initializes, if empty to MGeomCam::GetNumSectors() for: +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +// Initializes TArray's to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively +// Fills with number of valid pixels (if !MBadPixelsPix::IsBad()): +// - MHCalibrationCam::fAverageAreaNum[area index] +// - MHCalibrationCam::fAverageSectorNum[area index] +// +// Calls InitializeHists() for every entry in: +// - MHCalibrationCam::fHiGainArray +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +// Sets Titles and Names for the Histograms +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +// Retrieves the run numbers from MRawRunHeader and stores them in fRunNumbers +// +Bool_t MHCalibrationCam::ReInit(MParList *pList) +{ + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nsectors = fGeom->GetNumSectors(); + const Int_t nareas = fGeom->GetNumAreas(); + + fBadPixels = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam"); + if (!fBadPixels) + { + + fBadPixels = (MBadPixelsCam*)pList->FindCreateObj(AddSerialNumber("MBadPixelsCam")); + if (!fBadPixels) + { + *fLog << err << "Cannot find nor create MBadPixelsCam ... abort." << endl; + return kFALSE; + } + else + fBadPixels->InitSize(npixels); + } + + // + // The function TArrayF::Set() already sets all entries to 0. + // + fAverageAreaNum. Set(nareas); + fAverageAreaSat. Set(nareas); + fAverageAreaSigma. Set(nareas); + fAverageAreaSigmaVar. Set(nareas); + fAverageAreaRelSigma. Set(nareas); + fAverageAreaRelSigmaVar.Set(nareas); + fAverageSectorNum. Set(nsectors); + fRunNumbers. Set(fRunNumbers.GetSize()+1); + + for (Int_t aidx=0; aidxGetRunNumber(); + if (IsLoGain()) + SetLoGain(fRunHeader->GetNumSamplesLoGain()); + } + + if (!ReInitHists(pList)) + return kFALSE; + + if (!fRunHeader) + return kTRUE; + + for (Int_t i=0; iGetEntries(); i++) + { + TH1F *h = (*this)[i].GetHGausHist(); + h->SetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," ")); + } + + if (IsLoGain()) + for (Int_t i=0; iGetEntries(); i++) + { + TH1F *h = (*this)(i).GetHGausHist(); + h->SetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," ")); + } + + for (Int_t j=0; jSetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," ")); + } + + if (IsLoGain()) + for (Int_t j=0; jSetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," ")); + } + + for (Int_t j=0; jSetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," ")); + } + + if (IsLoGain()) + for (Int_t j=0; jSetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," ")); + } + + return kTRUE; +} + + + +Bool_t MHCalibrationCam::ReInitHists(MParList *pList) +{ + return kTRUE; +} + +//-------------------------------------------------------------------------------- +// +// Retrieves from MGeomCam: +// - number of pixels +// - number of pixel areas +// - number of sectors +// +// For all TObjArray's (including the averaged ones), the following steps are performed: +// +// 1) Test size and return kFALSE if not matching +// 2) +// +Bool_t MHCalibrationCam::Fill(const MParContainer *par, const Stat_t w) +{ + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nareas = fGeom->GetNumAreas(); + const Int_t nsectors = fGeom->GetNumSectors(); + + if (fHiGainArray->GetEntries() != npixels) + { + *fLog << err << "ERROR - Size mismatch... abort." << endl; + return kFALSE; + } + + if (IsLoGain()) + if (fLoGainArray->GetEntries() != npixels) + { + *fLog << err << "ERROR - Size mismatch... abort." << endl; + return kFALSE; + } + + if (fAverageHiGainAreas->GetEntries() != nareas) + { + *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; + return kFALSE; + } + + if (IsLoGain()) + if (fAverageLoGainAreas->GetEntries() != nareas) + { + *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; + return kFALSE; + } + + if (fAverageHiGainSectors->GetEntries() != nsectors) + { + *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; + return kFALSE; + } + + if (IsLoGain()) + if (fAverageLoGainSectors->GetEntries() != nsectors) + { + *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; + return kFALSE; + } + + return FillHists(par, w); +} + +Bool_t MHCalibrationCam::FillHists(const MParContainer *par, const Stat_t w) +{ + *fLog << warn << GetDescriptor() << "FillHists not overloaded! Can't be used!" << endl; + return kFALSE; +} + +// -------------------------------------------------------------------------- +// +// 0) Ask if fHiGainArray and fLoGainArray have been initialized, +// otherwise return kFALSE. +// 1) FinalizeHists() +// 2) FinalizeBadPixels() +// 3) CalcAverageSigma() +// +Bool_t MHCalibrationCam::Finalize() +{ + + if (fHiGainArray->GetEntries() == 0 && fLoGainArray->GetEntries() == 0) + { + *fLog << err << GetDescriptor() + << ": ERROR: Both (HiGain and LoGain) histogram arrays have not been initialized... abort." << endl; + return kFALSE; + } + + if (!FinalizeHists()) + return kFALSE; + + FinalizeBadPixels(); + CalcAverageSigma(); + + return kTRUE; +} + +Bool_t MHCalibrationCam::FinalizeHists() +{ + return kTRUE; +} + +void MHCalibrationCam::FinalizeBadPixels() +{ +} + + +// ------------------------------------------------------------- +// +// If MBadPixelsPix::IsBad(): +// - calls MHCalibrationPix::SetExcluded() +// +// Calls: +// - MHGausEvents::InitBins() +// - MHCalibrationPix::ChangeHistId(i) +// - MHCalibrationPix::SetEventFrequency(fPulserFrequency) +// +void MHCalibrationCam::InitHists(MHCalibrationPix &hist, MBadPixelsPix &bad, const Int_t i) +{ + + if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun)) + hist.SetExcluded(); + + hist.InitBins(); + hist.ChangeHistId(i); + hist.SetEventFrequency(fPulserFrequency); + + TH1F *h = hist.GetHGausHist(); + h->SetTitle( Form("%s%s", h->GetTitle()," Runs: ")); +} + +void MHCalibrationCam::FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp) +{ + + for (Int_t i=0; iGetSize(); i++) + { + + MHCalibrationPix &hist = (*this)[i]; + + if (hist.IsExcluded()) + continue; + + MCalibrationPix &pix = calcam[i]; + MBadPixelsPix &bad = badcam[i]; + + FitHiGainHists(hist,pix,bad,fittyp,osctyp); + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationPix &hist = GetAverageHiGainArea(j); + MCalibrationPix &pix = calcam.GetAverageArea(j); + MBadPixelsPix &bad = calcam.GetAverageBadArea(j); + + FitHiGainHists(hist,pix,bad,fittyp,osctyp); + } + + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationPix &hist = GetAverageHiGainSector(j); + MCalibrationPix &pix = calcam.GetAverageSector(j); + MBadPixelsPix &bad = calcam.GetAverageBadSector(j); + + FitHiGainHists(hist,pix,bad,fittyp,osctyp); + } + +} + +void MHCalibrationCam::FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp) +{ + + for (Int_t i=0; iGetSize(); i++) + { + + MHCalibrationPix &hist = (*this)(i); + + if (hist.IsExcluded()) + continue; + + MCalibrationPix &pix = calcam[i]; + MBadPixelsPix &bad = badcam[i]; + + FitLoGainHists(hist,pix,bad,fittyp,osctyp); + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationPix &hist = GetAverageLoGainArea(j); + MCalibrationPix &pix = calcam.GetAverageArea(j); + MBadPixelsPix &bad = calcam.GetAverageBadArea(j); + + FitLoGainHists(hist,pix,bad,fittyp,osctyp); + } + + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationPix &hist = GetAverageLoGainSector(j); + MCalibrationPix &pix = calcam.GetAverageSector(j); + MBadPixelsPix &bad = calcam.GetAverageBadSector(j); + + FitLoGainHists(hist,pix,bad,fittyp,osctyp); + } +} + +//------------------------------------------------------------ +// +// For all averaged areas, the fitted sigma is multiplied with the square root of +// the number involved pixels +// +void MHCalibrationCam::CalcAverageSigma() +{ + + if (!fCam) + return; + + for (UInt_t j=0; jGetNumAreas(); j++) + { + + MCalibrationPix &pix = fCam->GetAverageArea(j); + + const Float_t numsqr = TMath::Sqrt((Float_t)fAverageAreaNum[j]); + fAverageAreaSigma[j] = pix.GetSigma () * numsqr; + fAverageAreaSigmaVar[j] = pix.GetSigmaErr () * pix.GetSigmaErr() * numsqr; + + pix.SetSigma (fAverageAreaSigma[j]); + pix.SetSigmaVar(fAverageAreaSigmaVar[j]); + + fAverageAreaRelSigma [j] = fAverageAreaSigma[j] / pix.GetMean(); + fAverageAreaRelSigmaVar[j] = fAverageAreaSigmaVar[j] / (fAverageAreaSigma[j]*fAverageAreaSigma[j]); + fAverageAreaRelSigmaVar[j] += pix.GetMeanRelVar(); + fAverageAreaRelSigmaVar[j] *= fAverageAreaRelSigma[j]; + } +} + +// --------------------------------------------------------------------------- +// +// Returns if the histogram is empty and sets the following flag: +// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) +// +// Fits the histograms with a Gaussian, in case of failure +// calls MHCalibrationPix::RepeatFit(), in case of repeated failure +// calls MHCalibrationPix::BypassFit() and sets the following flags: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp ) +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). +// In case no, sets the following flags: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp ) +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Retrieves the results and store them in MCalibrationPix +// +void MHCalibrationCam::FitHiGainHists(MHCalibrationPix &hist, + MCalibrationPix &pix, + MBadPixelsPix &bad, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp) +{ + + + if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow()) + return; + + // + // 2) Fit the Hi Gain histograms with a Gaussian + // + if (!hist.FitGaus()) + // + // 3) In case of failure set the bit Fitted to false and take histogram means and RMS + // + if (!hist.RepeatFit()) + { + hist.BypassFit(); + bad.SetUncalibrated( fittyp ); + } + + hist.Renorm(); + // + // 4) Check for oscillations + // + if (IsOscillations()) + { + hist.CreateFourierSpectrum(); + + if (!hist.IsFourierSpectrumOK()) + bad.SetUncalibrated( osctyp ); + } + + // + // 5) Retrieve the results and store them in this class + // + pix.SetHiGainMean ( hist.GetMean() ); + pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() ); + pix.SetHiGainSigma ( hist.GetSigma() ); + pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() ); + pix.SetHiGainProb ( hist.GetProb() ); + pix.SetHiGainNumBlackout( hist.GetBlackout() ); + pix.SetHiGainNumPickup ( hist.GetPickup() ); + + if (IsDebug()) + { + *fLog << dbginf << GetDescriptor() << ": ID " << hist.GetPixId() + << " HiGainSaturation: " << pix.IsHiGainSaturation() + << " HiGainMean: " << hist.GetMean () + << " HiGainMeanErr: " << hist.GetMeanErr () + << " HiGainMeanSigma: " << hist.GetSigma () + << " HiGainMeanSigmaErr: " << hist.GetSigmaErr() + << " HiGainMeanProb: " << hist.GetProb () + << " HiGainNumBlackout: " << hist.GetBlackout() + << " HiGainNumPickup : " << hist.GetPickup () + << endl; + } + +} + + +// --------------------------------------------------------------------------- +// +// Returns if the histogram is empty and sets the following flag: +// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) +// +// Fits the histograms with a Gaussian, in case of failure +// calls MHCalibrationPix::RepeatFit(), in case of repeated failure +// calls MHCalibrationPix::BypassFit() and sets the following flags: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp ) +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). +// In case no, sets the following flags: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp ) +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Retrieves the results and store them in MCalibrationPix +// +void MHCalibrationCam::FitLoGainHists(MHCalibrationPix &hist, + MCalibrationPix &pix, + MBadPixelsPix &bad, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp) +{ + + if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow()) + return; + + + // + // 2) Fit the Hi Gain histograms with a Gaussian + // + if (!hist.FitGaus()) + // + // 3) In case of failure set the bit Fitted to false and take histogram means and RMS + // + if (!hist.RepeatFit()) + { + hist.BypassFit(); + bad.SetUncalibrated( fittyp ); + } + + // + // 4) Check for oscillations + // + if (IsOscillations()) + { + hist.CreateFourierSpectrum(); + + if (!hist.IsFourierSpectrumOK()) + bad.SetUncalibrated( osctyp ); + } + + // + // 5) Retrieve the results and store them in this class + // + pix.SetLoGainMean ( hist.GetMean() ); + pix.SetLoGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() ); + pix.SetLoGainSigma ( hist.GetSigma() ); + pix.SetLoGainSigmaVar ( hist.GetSigmaErr() * hist.GetSigmaErr() ); + pix.SetLoGainProb ( hist.GetProb() ); + pix.SetLoGainNumBlackout( hist.GetBlackout() ); + pix.SetLoGainNumPickup ( hist.GetPickup() ); + + if (IsDebug()) + { + *fLog << dbginf << GetDescriptor() << "ID: " << hist.GetPixId() + << " HiGainSaturation: " << pix.IsHiGainSaturation() + << " LoGainMean: " << hist.GetMean () + << " LoGainMeanErr: " << hist.GetMeanErr () + << " LoGainMeanSigma: " << hist.GetSigma () + << " LoGainMeanSigmaErr: " << hist.GetSigmaErr() + << " LoGainMeanProb: " << hist.GetProb () + << " LoGainNumBlackout: " << hist.GetBlackout() + << " LoGainNumPickup : " << hist.GetPickup () + << endl; + } + +} + + + +// -------------------------------------------------------------------------- +// +// Dummy, needed by MCamEvent +// +Bool_t MHCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// What MHCamera needs in order to draw an individual pixel in the camera +// +void MHCalibrationCam::DrawPixelContent(Int_t idx) const +{ +} + +// ----------------------------------------------------------------------------- +// +// Default draw: +// +// Displays the averaged areas, both High Gain and Low Gain +// +// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options +// +void MHCalibrationCam::Draw(const Option_t *opt) +{ + + const Int_t nareas = fAverageHiGainAreas->GetEntries(); + if (nareas == 0) + return; + + TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this); + pad->SetBorderMode(0); + + pad->Divide(IsLoGain() ? 2 : 1,nareas); + + for (Int_t i=0; icd(IsLoGain() ? 2*(i+1)-1 : i+1); + GetAverageHiGainArea(i).Draw(opt); + + if (!fAverageAreaSat[i]) + DrawAverageSigma(fAverageAreaSat[i], i, + fAverageAreaSigma[i], fAverageAreaSigmaVar[i], + fAverageAreaRelSigma[i], fAverageAreaRelSigmaVar[i]); + + if (IsLoGain()) + { + pad->cd(2*(i+1)); + GetAverageLoGainArea(i).Draw(opt); + } + + if (fAverageAreaSat[i]) + DrawAverageSigma(fAverageAreaSat[i], i, + fAverageAreaSigma[i], fAverageAreaSigmaVar[i], + fAverageAreaRelSigma[i], fAverageAreaRelSigmaVar[i]); + } + +} + +// ----------------------------------------------------------------------------- +// +// Default draw: +// +// Displays a TPaveText with the re-normalized sigmas of the average area +// +void MHCalibrationCam::DrawAverageSigma(Bool_t sat, Bool_t inner, + Float_t sigma, Float_t sigmavar, + Float_t relsigma, Float_t relsigmavar) const +{ + + if (sigma != 0 && sigmavar >= 0 && relsigmavar >= 0.) + { + + TPad *newpad = new TPad("newpad","transparent",0,0,1,1); + newpad->SetFillStyle(4000); + newpad->Draw(); + newpad->cd(); + + TPaveText *text = new TPaveText(sat? 0.1 : 0.35,0.7,sat ? 0.4 : 0.7,1.0); + text->SetTextSize(0.07); + const TString line1 = Form("%s%s%s",inner ? "Outer" : "Inner", + " Pixels ", sat ? "Low Gain" : "High Gain"); + TText *txt1 = text->AddText(line1.Data()); + const TString line2 = Form("#sigma per pix: %2.2f #pm %2.2f",sigma,TMath::Sqrt(sigmavar)); + TText *txt2 = text->AddText(line2.Data()); + const TString line3 = Form("Rel. #sigma per pix: %2.2f #pm %2.2f",relsigma,TMath::Sqrt(relsigmavar)); + TText *txt3 = text->AddText(line3.Data()); + text->Draw(""); + + text->SetBit(kCanDelete); + txt1->SetBit(kCanDelete); + txt2->SetBit(kCanDelete); + txt3->SetBit(kCanDelete); + newpad->SetBit(kCanDelete); + } +} + +Int_t MHCalibrationCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "Debug", print)) + { + SetDebug(GetEnvValue(env, prefix, "Debug", IsDebug())); + rc = kTRUE; + } + if (IsEnvDefined(env, prefix, "LoGain", print)) + { + SetDebug(GetEnvValue(env, prefix, "LoGain", IsLoGain())); + rc = kTRUE; + } + if (IsEnvDefined(env, prefix, "Oscillations", print)) + { + SetDebug(GetEnvValue(env, prefix, "Oscillations", IsOscillations())); + rc = kTRUE; + } + + + return rc; +} Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationCam.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationCam.h (revision 4929) @@ -0,0 +1,192 @@ +#ifndef MARS_MHCalibrationCam +#define MARS_MHCalibrationCam + +#ifndef ROOT_TArrayI +#include +#endif + +#ifndef ROOT_TArrayF +#include +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef MARS_MCamEvent +#include "MCamEvent.h" +#endif + +#ifndef MARS_MBadPixelsPix +#include "MBadPixelsPix.h" +#endif + +#ifndef MARS_MCalibrationCam +#include "MCalibrationCam.h" +#endif + +class TText; +class TObjArray; + +class MHCalibrationPix; +class MGeomCam; +class MRawRunHeader; +class MCalibrationIntensityCam; +class MCalibrationCam; +class MCalibrationPix; +class MBadPixelsCam; +class MBadPixelsPix; + +class MHCalibrationCam : public MH, public MCamEvent +{ + +private: + + static const Int_t fgAverageNbins; //! The default for fAverageNbins (now set to: 2000) + static const Int_t fgPulserFrequency; //! The default for fPulserFrequency (now set to: 500) + +protected: + + Float_t fNumHiGainSaturationLimit; // Rel. amount sat. higain FADC slices until pixel is called saturated + Float_t fNumLoGainSaturationLimit; // Rel. amount sat. logain FADC slices until pixel is called saturated + + TArrayI fAverageAreaNum; // Number of pixels in average pixels per area + TArrayF fAverageAreaRelSigma; // Re-normalized relative sigmas in average pixels per area + TArrayF fAverageAreaRelSigmaVar; // Variance Re-normalized relative sigmas in average pixels per area + TArrayI fAverageAreaSat; // Number of saturated slices in average pixels per area + TArrayF fAverageAreaSigma; // Re-normalized sigmas in average pixels per area + TArrayF fAverageAreaSigmaVar; // Variance Re-normalized sigmas in average pixels per area + Int_t fAverageNbins; // Number of bins for the average histograms + TObjArray *fAverageHiGainAreas; //-> Array of calibration pixels, one per pixel area + TObjArray *fAverageHiGainSectors; //-> Array of calibration pixels, one per camera sector + TObjArray *fAverageLoGainAreas; //-> Array of calibration pixels, one per pixel area + TObjArray *fAverageLoGainSectors; //-> Array of calibration pixels, one per camera sector + TArrayI fAverageSectorNum; // Number of pixels in average pixels per sector + TArrayI fRunNumbers; // Numbers of runs used + + MCalibrationCam::PulserColor_t fColor; // Colour of the pulsed LEDs + + MBadPixelsCam *fBadPixels; //! Bad Pixels storage container + MCalibrationIntensityCam *fIntensCam; //! Intensity Calibration Cam with the results + MCalibrationCam *fCam; //! Calibration Cam with the results + MGeomCam *fGeom; //! Camera geometry + MRawRunHeader *fRunHeader; //! Run Header + + TObjArray *fHiGainArray; //-> Array of calibration pixels, one per pixel + TObjArray *fLoGainArray; //-> Array of calibration pixels, one per pixel + + Int_t fPulserFrequency; // Light pulser frequency + + enum { kDebug, kLoGain, kOscillations }; // Possible flags + + Byte_t fFlags; // Bit-field to hold the flags + + virtual Bool_t SetupHists(const MParList *pList); + virtual Bool_t ReInitHists(MParList *pList); + virtual Bool_t FillHists(const MParContainer *par, const Stat_t w=1); + virtual Bool_t FinalizeHists(); + virtual void FinalizeBadPixels(); + + virtual void CalcAverageSigma(); + + void DrawAverageSigma(Bool_t sat, Bool_t inner, + Float_t sigma, Float_t sigmaerr, + Float_t relsigma, Float_t relsigmaerr) const; + + void FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp); + + void FitHiGainHists(MHCalibrationPix &hist, + MCalibrationPix &pix, + MBadPixelsPix &bad, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp); + + void FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp); + + void FitLoGainHists(MHCalibrationPix &hist, + MCalibrationPix &pix, + MBadPixelsPix &bad, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp); + + void InitHists(MHCalibrationPix &hist, MBadPixelsPix &bad, const Int_t i); + + Bool_t IsDebug () const { return TESTBIT(fFlags,kDebug); } + Bool_t IsOscillations() const { return TESTBIT(fFlags,kOscillations); } + Bool_t IsLoGain () const { return TESTBIT(fFlags,kLoGain); } + + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + +public: + + MHCalibrationCam(const char *name=NULL, const char *title=NULL); + virtual ~MHCalibrationCam(); + + Bool_t SetupFill(const MParList *pList); + Bool_t ReInit ( MParList *pList); + Bool_t Fill (const MParContainer *par, const Stat_t w=1); + Bool_t Finalize ( ); + + virtual void ResetHists(); + + // Clone + TObject *Clone(const char *name="") const; + + // Draw + void Draw(const Option_t *opt); + + Bool_t GetPixelContent ( Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; + void DrawPixelContent( Int_t num ) const; + + const Int_t GetAverageAreas () const; + MHCalibrationPix &GetAverageHiGainArea (UInt_t i); + const MHCalibrationPix &GetAverageHiGainArea (UInt_t i) const; + MHCalibrationPix &GetAverageLoGainArea (UInt_t i); + const MHCalibrationPix &GetAverageLoGainArea (UInt_t i) const; + MHCalibrationPix &GetAverageHiGainSector(UInt_t i); + const MHCalibrationPix &GetAverageHiGainSector(UInt_t i) const; + MHCalibrationPix &GetAverageLoGainSector(UInt_t i); + const MHCalibrationPix &GetAverageLoGainSector(UInt_t i) const; + const Int_t GetAverageSectors () const; + const MCalibrationCam::PulserColor_t GetColor () const { return fColor; } + const Float_t GetNumHiGainSaturationLimit() const { return fNumHiGainSaturationLimit; } + const Float_t GetNumLoGainSaturationLimit() const { return fNumLoGainSaturationLimit; } + const TArrayI &GetRunNumbers () const; + const Int_t GetSize () const; + + MHCalibrationPix &operator[] (UInt_t i); + const MHCalibrationPix &operator[] (UInt_t i) const; + MHCalibrationPix &operator() (UInt_t i); + const MHCalibrationPix &operator() (UInt_t i) const; + + void SetColor ( const MCalibrationCam::PulserColor_t color) { fColor = color; } + void SetDebug ( const Bool_t b=kTRUE ) { b + ? SETBIT(fFlags,kDebug) + : CLRBIT(fFlags,kDebug); } + void SetLoGain ( const Bool_t b=kTRUE ) { b + ? SETBIT(fFlags,kLoGain) + : CLRBIT(fFlags,kLoGain); } + void SetOscillations ( const Bool_t b=kTRUE ) { b + ? SETBIT(fFlags,kOscillations) + : CLRBIT(fFlags,kOscillations); } + void SetAverageNbins ( const Int_t bins=fgAverageNbins) { fAverageNbins = bins; } + void SetNumLoGainSaturationLimit( const Float_t lim ) { fNumLoGainSaturationLimit = lim; } + void SetNumHiGainSaturationLimit( const Float_t lim ) { fNumHiGainSaturationLimit = lim; } + void SetPulserFrequency ( const Int_t f=fgPulserFrequency) { fPulserFrequency = f; } + + ClassDef(MHCalibrationCam, 1) // Base Histogram class for Calibration Camera +}; + +#endif + + + + + + + + + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindCam.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindCam.cc (revision 4929) @@ -0,0 +1,333 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationChargeBlindCam +// +// Histogram class for blind pixels in the camera. Incorporates the TObjArray's: +// - fBlindPixelsArray (for calibrated High Gains per pixel) +// +///////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationChargeBlindCam.h" +#include "MHCalibrationChargeBlindPix.h" + +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MCalibrationChargeBlindPix.h" +#include "MCalibrationChargeBlindCam.h" + +#include "MExtractedSignalBlindPixel.h" + +#include "MParList.h" + +#include "MRawRunHeader.h" + +ClassImp(MHCalibrationChargeBlindCam); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - all pointers to NULL +// +// Initializes and sets owner of: +// - fBlindPixelsArray +// +// Sets fFitFunc to kEPoisson4 +// +MHCalibrationChargeBlindCam::MHCalibrationChargeBlindCam(const char *name, const char *title) + : fCam(NULL), fRunHeader(NULL) +{ + + fName = name ? name : "MHCalibrationChargeBlindCam"; + fTitle = title ? title : "Class to fille the blind pixel histograms"; + + fBlindPixelsArray = new TObjArray; + fBlindPixelsArray->SetOwner(); + + fFitFunc = MHCalibrationChargeBlindPix::kEPoisson4; + +} + +// -------------------------------------------------------------------------- +// +// Deletes the TClonesArray of: +// - fBlindPixelsArray +// +MHCalibrationChargeBlindCam::~MHCalibrationChargeBlindCam() +{ + delete fBlindPixelsArray; +} + +// -------------------------------------------------------------------------- +// +// Get i-th High Gain pixel (pixel number) +// +MHCalibrationChargeBlindPix &MHCalibrationChargeBlindCam::operator[](UInt_t i) +{ + return *static_cast(fBlindPixelsArray->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th High Gain pixel (pixel number) +// +const MHCalibrationChargeBlindPix &MHCalibrationChargeBlindCam::operator[](UInt_t i) const +{ + return *static_cast(fBlindPixelsArray->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Our own clone function is necessary since root 3.01/06 or Mars 0.4 +// I don't know the reason. +// +// Creates new MHCalibrationChargeBlindCam +// Deletes the TObjArray's and Clones them individually +// +TObject *MHCalibrationChargeBlindCam::Clone(const char *name) const +{ + + const Int_t nhi = fBlindPixelsArray->GetEntries(); + + // + // FIXME, this might be done faster and more elegant, by direct copy. + // + MHCalibrationChargeBlindCam *cam = new MHCalibrationChargeBlindCam(); + + cam->fBlindPixelsArray->Expand(nhi); + + for (int i=0; ifBlindPixelsArray)[i] = (*fBlindPixelsArray)[i]->Clone(); + + return cam; +} + +// -------------------------------------------------------------------------- +// +// Gets the pointers to: +// - MRunHeader +// - MExtractedSignalBlindPix +// +// Calls Delete-Function of: +// - MHCalibrationChargeBlindCam::fBlindPixelsArray +// +Bool_t MHCalibrationChargeBlindCam::SetupFill(const MParList *pList) +{ + + fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << warn << GetDescriptor() + << ": MRawRunHeader not found... will not store run numbers." << endl; + return kFALSE; + } + + fSignal = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel"); + if (!fSignal) + { + *fLog << err << "MExtractedSignalBlindPixel not found... aborting " << endl; + return kFALSE; + } + + fBlindPixelsArray->Delete(); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Initializes, if empty to MExtractedSignalCam::GetSize() for: +// - MHCalibrationChargeBlindCam::fBlindPixelsArray +// +// Calls InitializeHists() for every entry in: +// - MHCalibrationChargeBlindCam::fBlindPixelsArray +// +// Retrieves the run numbers from MRawRunHeader and stores them in fRunNumbers +// +Bool_t MHCalibrationChargeBlindCam::ReInit(MParList *pList) +{ + + const UInt_t nblindpixels = fSignal->GetNumBlindPixels(); + + Int_t runnr = 0; + + if (fRunHeader) + runnr = fRunHeader->GetRunNumber(); + + fCam = (MCalibrationChargeBlindCam*)pList->FindCreateObj("MCalibrationChargeBlindCam"); + if (!fCam) + { + *fLog << err << "Cannot find nor create MCalibrationChargeBlindCam ... abort." << endl; + return kFALSE; + } + + if (fCam->GetNumBlindPixels() != nblindpixels) + { + *fLog << err; + *fLog << "Size mismatch in MCalibrationChargeBlindCam ... abort." << endl; + *fLog << " Size of MCalibrationChargeBlindCam: " << fCam->GetNumBlindPixels() << endl; + *fLog << " Size of MExtractedSignalBlindPixel: " << nblindpixels << endl; + return kFALSE; + } + + + const Int_t samples = fSignal->GetNumFADCSamples(); + const Int_t integ = fSignal->IsExtractionType( MExtractBlindPixel::kIntegral ); + + if (fBlindPixelsArray->GetEntries()==0) + { + + fBlindPixelsArray->Expand(nblindpixels); + + for (UInt_t i=0; iSetTitle( Form("%s%s", h->GetTitle()," Runs: ")); + (*this)[i].ChangeFitFunc(fFitFunc); + (*this)[i].SetupFill(pList); + (*this)[i].SetCalibrationChargeBlindPix(&(*fCam)[i]); + } + } + + for (UInt_t i=0; iSetTitle( Form("%s%i%s", h->GetTitle(),runnr," ")); + } + + return kTRUE; +} + + +//-------------------------------------------------------------------------------- +// +// Retrieves from MExtractedSignalBlindPixel: +// - number of blind pixels +// +// For all TObjArray's, the following steps are performed: +// +// 1) Test size and return kFALSE if not matching +// 2) +// +Bool_t MHCalibrationChargeBlindCam::Fill(const MParContainer *par, const Stat_t w) +{ + + const Int_t nblindpixels = fSignal->GetNumBlindPixels(); + + if (GetSize() != nblindpixels) + { + gLog << err << "ERROR - Size mismatch... abort." << endl; + return kFALSE; + } + + for (Int_t i=0; iGetEntries(); + + if (size == 0) + return; + + TString option(opt); + option.ToLower(); + + TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this); + pad->SetBorderMode(0); + + switch (size) + { + case 1: + break; + case 2: + pad->Divide(2,1); + break; + case 3: + case 4: + pad->Divide(2,2); + break; + default: + pad->Divide(size/2+1,size/2+1); + break; + } + + for (Int_t i=0; icd(i+1); + (*this)[i].Draw(option); + } + + pad->Modified(); + pad->Update(); + +} Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindCam.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindCam.h (revision 4929) @@ -0,0 +1,65 @@ +#ifndef MARS_MHCalibrationChargeBlindCam +#define MARS_MHCalibrationChargeBlindCam + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef MARS_MHCalibrationChargeBlindPix +#include "MHCalibrationChargeBlindPix.h" +#endif + +class TObjArray; + +class MRawRunHeader; +class MExtractedSignalBlindPixel; +class MCalibrationChargeBlindCam; +class MHCalibrationChargeBlindPix; + +class MHCalibrationChargeBlindCam : public MH +{ +private: + + MExtractedSignalBlindPixel *fSignal; //! Extracted Signal class + MCalibrationChargeBlindCam *fCam; //! Calibration Cam with the results + MRawRunHeader *fRunHeader; //! Run Header + + TObjArray *fBlindPixelsArray; //-> Array of calibration pixels, one per pixel + + MHCalibrationChargeBlindPix::FitFunc_t fFitFunc; + +public: + + MHCalibrationChargeBlindCam(const char *name=NULL, const char *title=NULL); + ~MHCalibrationChargeBlindCam(); + + Bool_t SetupFill( const MParList *pList); + Bool_t ReInit ( MParList *pList); + Bool_t Fill ( const MParContainer *par, const Stat_t w=1); + Bool_t Finalize ( ); + + // Clone + TObject *Clone(const char *name="") const; + + // Draw + void Draw(Option_t *opt=""); + + const Int_t GetSize() const { return fBlindPixelsArray->GetSize(); } + + MHCalibrationChargeBlindPix &operator[] (UInt_t i); + const MHCalibrationChargeBlindPix &operator[] (UInt_t i) const; + + void SetFitFunc(const MHCalibrationChargeBlindPix::FitFunc_t func) { fFitFunc = func; } + + ClassDef(MHCalibrationChargeBlindCam, 1) // Histogram class for Blind Pixel Calibration +}; + +#endif + + + + + + + + + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc (revision 4929) @@ -0,0 +1,1178 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationChargeBlindPix +// +// Histogram class for the charge calibration of the Blind Pixel. +// Stores and fits the charges and stores the averaged assumed pedestal and +// single-phe FADC slice entries. Charges are taken from MExtractedSignalBlindPix. +// Performs the Single Photo-electron fit to extract the Poisson mean and its errors. +// +// Different fits can be chosen with the function ChangeFitFunc(). +// +// The fit result is accepted under the condition that: +// 1) the Probability is greater than fProbLimit (default 0.001 == 99.7%) +// 2) at least fNumSinglePheLimit events are found in the single Photo-electron peak +// +// The single FADC slice entries are averaged and stored in fASinglePheFADCSlices, if +// their sum exceeds fSinglePheCut, otherwise in fAPedestalFADCSlices. +// +// Used numbers are the following: +// +// Electronic conversion factor: +// Assume, we have N_e electrons at the anode, +// thus a charge of N_e*e (e = electron charge) Coulomb. +// +// This charge is AC coupled and runs into a R_pre = 50 Ohm resistency. +// The corresponding current is amplified by a gain factor G_pre = 400 +// (the precision of this value still has to be checked !!!) and again AC coupled to +// the output. +// The corresponding signal goes through the whole transmission and +// amplification chain and is digitized in the FADCs. +// The conversion Signal Area to FADC counts (Conv_trans) has been measured +// by David and Oscar to be approx. 3.9 pVs^-1 +// +// Thus: Conversion FADC counts to Number of Electrons at Anode: +// FADC counts = (1/Conv_tran) * G_pre * R_pre * e * N_e = 8 * 10^-4 N_e. +// +// Also: FADC counts = 8*10^-4 * GAIN * N_phe +// +// In the blind pixel, there is an additional pre-amplifier with an amplification of +// about 10. Therefore, we have for the blind pixel: +// +// +// FADC counts (Blind Pixel) = 8*10^-3 * GAIN * N_phe +// +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationChargeBlindPix.h" +#include "MExtractBlindPixel.h" + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MRawEvtData.h" +#include "MRawEvtPixelIter.h" + +#include "MExtractedSignalBlindPixel.h" +#include "MCalibrationChargeBlindPix.h" + +ClassImp(MHCalibrationChargeBlindPix); + +using namespace std; + +const Double_t MHCalibrationChargeBlindPix::gkElectronicAmp = 0.008; +const Double_t MHCalibrationChargeBlindPix::gkElectronicAmpErr = 0.002; +const Float_t MHCalibrationChargeBlindPix::gkSignalInitializer = -9999.; + +const Int_t MHCalibrationChargeBlindPix::fgChargeNbins = 128; +const Axis_t MHCalibrationChargeBlindPix::fgChargeFirst = -0.5; +const Axis_t MHCalibrationChargeBlindPix::fgChargeLast = 511.5; +const Float_t MHCalibrationChargeBlindPix::fgSinglePheCut = 20.; +const Float_t MHCalibrationChargeBlindPix::fgNumSinglePheLimit = 50.; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - the default number for fNbins (fgChargeNbins) +// - the default number for fFirst (fgChargeFirst) +// - the default number for fLast (fgChargeLast) +// - the default number for fSinglePheCut (fgSingePheCut) +// - the default number for fNumSinglePheLimit (fgNumSinglePheLimit) +// - the default number of bins after stripping (30) +// +// - the default name of the fHGausHist ("HCalibrationChargeBlindPix") +// - the default title of the fHGausHist ("Distribution of Summed FADC slices Blind Pixel ") +// - the default x-axis title for fHGausHist ("Sum FADC Slices") +// - the default y-axis title for fHGausHist ("Nr. of events") +// +// Initializes: +// - all pointers to NULL +// - fASinglePheFADCSlices(0); +// - fAPedestalFADCSlices(0); +// - fPixId to 0 +// +// Calls: +// - Clear() +// +MHCalibrationChargeBlindPix::MHCalibrationChargeBlindPix(const char *name, const char *title) + : fBlindPix(NULL), fSignal(NULL), fRawEvt(NULL), + fSinglePheFit(NULL), + fFitLegend(NULL), + fHSinglePheFADCSlices(NULL), fHPedestalFADCSlices(NULL) +{ + + fName = name ? name : "MHCalibrationChargeBlindPix"; + fTitle = title ? title : "Statistics of the FADC sums of Blind Pixel calibration events"; + + SetNbins( fgChargeNbins ); + SetFirst( fgChargeFirst ); + SetLast ( fgChargeLast ); + + fASinglePheFADCSlices.ResizeTo(1); + fAPedestalFADCSlices.ResizeTo(1); + + SetSinglePheCut(); + SetNumSinglePheLimit(); + SetProbLimit(0.001); + SetBinsAfterStripping(0); + + fHGausHist.SetName("HCalibrationChargeBlindPix"); + fHGausHist.SetTitle("Distribution of Summed FADC slices Blind Pixel"); + fHGausHist.SetXTitle("Signal Amplitude"); + fHGausHist.SetYTitle("Nr. of events"); + + fPixId = 0; + + Clear(); +} + +// -------------------------------------------------------------------------- +// +// Default Destructor. +// +// Deletes (if Pointer is not NULL): +// +// - fSinglePheFit +// - fFitLegend +// - fHSinglePheFADCSlices +// - fHPedestalFADCSlices +// +MHCalibrationChargeBlindPix::~MHCalibrationChargeBlindPix() +{ + + if (fSinglePheFit) + delete fSinglePheFit; + + if (fFitLegend) + delete fFitLegend; + + if (fHSinglePheFADCSlices) + delete fHSinglePheFADCSlices; + + if (fHPedestalFADCSlices) + delete fHPedestalFADCSlices; + +} + +// -------------------------------------------------------------------------- +// +// Sets: +// - all variables to 0., except the fit result variables to gkSignalInitializer +// - all flags to kFALSE +// - all pointers to NULL +// - the default fit function (kEPoisson5) +// +// Deletes: +// - all pointers unequal NULL +// +// Calls: +// - MHCalibrationChargePix::Clear() +// +void MHCalibrationChargeBlindPix::Clear(Option_t *o) +{ + + fLambda = gkSignalInitializer; + fMu0 = gkSignalInitializer; + fMu1 = gkSignalInitializer; + fSigma0 = gkSignalInitializer; + fSigma1 = gkSignalInitializer; + fLambdaErr = gkSignalInitializer; + fMu0Err = gkSignalInitializer; + fMu1Err = gkSignalInitializer; + fSigma0Err = gkSignalInitializer; + fSigma1Err = gkSignalInitializer; + + fLambdaCheck = gkSignalInitializer; + fLambdaCheckErr = gkSignalInitializer; + + // fFitFunc = kEMichele; + fFitFunc = kEPoisson5; + + fNumSinglePhes = 0; + fNumPedestals = 0; + + fChisquare = 0.; + fNDF = 0 ; + fProb = 0.; + + SetSinglePheFitOK ( kFALSE ); + SetPedestalFitOK ( kFALSE ); + + if (fFitLegend) + { + delete fFitLegend; + fFitLegend = NULL; + } + + if (fSinglePheFit) + { + delete fSinglePheFit; + fSinglePheFit = NULL; + } + + if (fHSinglePheFADCSlices) + { + delete fHSinglePheFADCSlices; + fHSinglePheFADCSlices = NULL; + } + + if (fHPedestalFADCSlices) + { + delete fHPedestalFADCSlices; + fHPedestalFADCSlices = NULL; + } + + MHCalibrationPix::Clear(); + return; +} + +// -------------------------------------------------------------------------- +// +// Empty function to overload MHCalibrationChargePix::Reset() +// +void MHCalibrationChargeBlindPix::Reset() +{ +} + +/* +// -------------------------------------------------------------------------- +// +// Our own clone function is necessary since root 3.01/06 or Mars 0.4 +// I don't know the reason. +// +// Creates new MHCalibrationCam +// +TObject *MHCalibrationChargeBlindPix::Clone(const char *) const +{ + + MHCalibrationChargeBlindPix *pix = new MHCalibrationChargeBlindPix(); + this->Copy(*pix); + + this->fHGausHist.Copy(pix->fHGausHist); + this->fSinglePheFit->Copy(*(pix->fSinglePheFit)); + this->fHSinglePheFADCSlices->Copy(*(pix->fHSinglePheFADCSlices)); + this->fHPedestalFADCSlices->Copy(*(pix->fHPedestalFADCSlices)); + + + return pix; +} +*/ + +// -------------------------------------------------------------------------- +// +// Set bit kSinglePheFitOK from outside +// +void MHCalibrationChargeBlindPix::SetSinglePheFitOK (const Bool_t b ) +{ + b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK); +} + +// -------------------------------------------------------------------------- +// +// Set bit kPedestalFitOK from outside +// +void MHCalibrationChargeBlindPix::SetPedestalFitOK(const Bool_t b) +{ + b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK); +} + +// -------------------------------------------------------------------------- +// +// Ask for status of bit kSinglePheFitOK +// +const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const +{ + return TESTBIT(fFlags,kSinglePheFitOK); +} + +// -------------------------------------------------------------------------- +// +// Ask for status of bit kPedestalFitOK +// +const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const +{ + return TESTBIT(fFlags,kPedestalFitOK); +} + +// -------------------------------------------------------------------------- +// +// Gets the pointers to: +// - MRawEvtData +// - MExtractedSignalBlindPixel +// +Bool_t MHCalibrationChargeBlindPix::SetupFill(const MParList *pList) +{ + + fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); + if (!fRawEvt) + { + *fLog << err << "MRawEvtData not found... aborting." << endl; + return kFALSE; + } + + fSignal = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel"); + if (!fSignal) + { + *fLog << err << "MExtractedSignalBlindPixel not found... aborting " << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// - MCalibrationChargeBlindPix +// +// Calls: +// - MHGausHist::InitBins() +// +Bool_t MHCalibrationChargeBlindPix::ReInit(MParList *pList) +{ + + fBlindPix = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix"); + if (!fBlindPix) + return kFALSE; + + const Int_t samples = fSignal->GetNumFADCSamples(); + const Int_t integ = fSignal->IsExtractionType( MExtractBlindPixel::kIntegral ); + + // + // Modify the histogram size in case, integrals have been used + // + if ( fLast < samples*integ*fgChargeLast ) + { + SetLast ( samples * (fgChargeLast+0.5) - 0.5 ); + SetSinglePheCut( samples * fgSinglePheCut ); + } + + MHCalibrationPix::InitBins(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Retrieves from MExtractedSignalBlindPixel: +// - number of FADC samples +// - extracted signal +// - blind Pixel ID +// +// Resizes (if necessary): +// - fASinglePheFADCSlices to sum of HiGain and LoGain samples +// - fAPedestalFADCSlices to sum of HiGain and LoGain samples +// +// Fills the following histograms: +// - MHGausEvents::FillHistAndArray(signal) +// +// Creates MRawEvtPixelIter, jumps to blind pixel ID, +// fills the vectors fASinglePheFADCSlices and fAPedestalFADCSlices +// with the full FADC slices, depending on the size of the signal w.r.t. fSinglePheCut +// +Bool_t MHCalibrationChargeBlindPix::Fill(const MParContainer *par, const Stat_t w) +{ + + const Int_t samples = (Int_t)fRawEvt->GetNumHiGainSamples() + +(Int_t)fRawEvt->GetNumLoGainSamples(); + + if (!fASinglePheFADCSlices.IsValid()) + { + fASinglePheFADCSlices.ResizeTo(samples); + fAPedestalFADCSlices.ResizeTo(samples); + } + + if (fASinglePheFADCSlices.GetNrows() != samples) + { + fASinglePheFADCSlices.ResizeTo(samples); + fAPedestalFADCSlices.ResizeTo(samples); + } + + Float_t slices = (Float_t)fSignal->GetNumFADCSamples(); + + if (slices == 0.) + { + *fLog << err + << "Number of used signal slices in MExtractedSignalBlindPix " + << "is zero ... abort." + << endl; + return kFALSE; + } + + // + // Signal extraction and histogram filling + // + const Float_t signal = fSignal->GetExtractedSignal(fPixId); + + if (signal > -0.5) + FillHist(signal); + else + return kTRUE; + + // + // In order to study the single-phe posistion, we extract the slices + // + const Int_t blindpixIdx = fSignal->GetBlindPixelIdx(fPixId); + + MRawEvtPixelIter pixel(fRawEvt); + pixel.Jump(blindpixIdx); + + if (signal > fSinglePheCut) + FillSinglePheFADCSlices(pixel); + else + FillPedestalFADCSlices(pixel); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Returns kFALSE, if empty +// +// - Creates the fourier spectrum and sets bit MHGausEvents::IsFourierSpectrumOK() +// - Retrieves the pedestals from MExtractedSignalBlindPixel +// - Normalizes fASinglePheFADCSlices and fAPedestalFADCSlices +// - Executes FitPedestal() +// - Executes FitSinglePhe() +// - Retrieves fit results and stores them in MCalibrationChargeBlindPix +// +Bool_t MHCalibrationChargeBlindPix::Finalize() +{ + + if (IsEmpty()) + { + *fLog << err << GetDescriptor() << " ID: " << fPixId + << " My histogram has not been filled !! " << endl; + return kFALSE; + } + + fMeanPedestal = fSignal->GetPed(); + fMeanPedestalErr = fSignal->GetPedErr(); + fSigmaPedestal = fSignal->GetPedRms(); + fSigmaPedestalErr = fSignal->GetPedRmsErr(); + + if (fNumSinglePhes > 1) + for (Int_t i=0;i 1) + for (Int_t i=0;iSetValid(kTRUE); + + if (FitSinglePhe()) + fBlindPix->SetSinglePheFitOK(); + else + fBlindPix->SetValid(IsPedestalFitOK()); + + fBlindPix->SetLambda ( fLambda ); + fBlindPix->SetLambdaVar ( fLambdaErr*fLambdaErr ); + fBlindPix->SetMu0 ( fMu0 ); + fBlindPix->SetMu0Err ( fMu0Err ); + fBlindPix->SetMu1 ( fMu1 ); + fBlindPix->SetMu1Err ( fMu1Err ); + fBlindPix->SetSigma0 ( fSigma0 ); + fBlindPix->SetSigma0Err ( fSigma0Err ); + fBlindPix->SetSigma1 ( fSigma1 ); + fBlindPix->SetSigma1Err ( fSigma1Err ); + fBlindPix->SetProb ( fProb ); + + fBlindPix->SetLambdaCheck ( fLambdaCheck ); + fBlindPix->SetLambdaCheckErr ( fLambdaCheckErr ); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Checks again for the size and fills fASinglePheFADCSlices with the FADC slice entries +// +void MHCalibrationChargeBlindPix::FillSinglePheFADCSlices(const MRawEvtPixelIter &iter) +{ + + const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples(); + + if (fASinglePheFADCSlices.GetNrows() < n) + fASinglePheFADCSlices.ResizeTo(n); + + Int_t i=0; + + Byte_t *start = iter.GetHiGainSamples(); + Byte_t *end = start + iter.GetNumHiGainSamples(); + + for (Byte_t *ptr = start; ptr < end; ptr++, i++) + fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr; + + start = iter.GetLoGainSamples(); + end = start + iter.GetNumLoGainSamples(); + + for (Byte_t *ptr = start; ptr < end; ptr++, i++) + fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr; + + fNumSinglePhes++; +} + +// -------------------------------------------------------------------------- +// +// Checks again for the size and fills fAPedestalFADCSlices with the FADC slice entries +// +void MHCalibrationChargeBlindPix::FillPedestalFADCSlices(const MRawEvtPixelIter &iter) +{ + + const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples(); + + if (fAPedestalFADCSlices.GetNrows() < n) + fAPedestalFADCSlices.ResizeTo(n); + + Int_t i = 0; + Byte_t *start = iter.GetHiGainSamples(); + Byte_t *end = start + iter.GetNumHiGainSamples(); + + for (Byte_t *ptr = start; ptr < end; ptr++, i++) + fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr; + + start = iter.GetLoGainSamples(); + end = start + iter.GetNumLoGainSamples(); + + for (Byte_t *ptr = start; ptr < end; ptr++, i++) + fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr; + + fNumPedestals++; +} + + +// -------------------------------------------------------------------------- +// +// Task to simulate single phe spectrum with the given parameters +// +Bool_t MHCalibrationChargeBlindPix::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1) +{ + + gRandom->SetSeed(); + + if (fHGausHist.GetIntegral() != 0) + { + *fLog << err << "Histogram " << fHGausHist.GetTitle() << " is already filled. " << endl; + *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl; + return kFALSE; + } + + if (!InitFit()) + return kFALSE; + + for (Int_t i=0;i<10000; i++) + fHGausHist.Fill(fSinglePheFit->GetRandom()); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// - Get the ranges from the stripped histogram +// - choose reasonable start values for the fit +// - initialize the fit function depending on fFitFunc +// - initialize parameter names and limits depending on fFitFunc +// +Bool_t MHCalibrationChargeBlindPix::InitFit() +{ + + // + // Get the fitting ranges + // + Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()); + Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()); + + if (rmin < 0.) + rmin = 0.; + + // + // First guesses for the fit (should be as close to reality as possible, + // otherwise the fit goes gaga because of high number of dimensions ... + // + const Stat_t entries = fHGausHist.Integral("width"); + const Double_t lambda_guess = 0.5; + //const Double_t maximum_bin = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + const Double_t norm = entries/TMath::Sqrt(TMath::TwoPi()); + + // + // Initialize the fit function + // + switch (fFitFunc) + { + case kEPoisson4: + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6); + rmin += 6.5; + break; + case kEPoisson5: + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6); + rmin = 0.; + break; + case kEPoisson6: + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6); + break; + case kEPolya: + fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8); + break; + case kEMichele: + fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,9); + break; + default: + *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; + return kFALSE; + break; + } + + if (!fSinglePheFit) + { + *fLog << warn << dbginf << "WARNING: Could not create fit function for Single Phe fit" << endl; + return kFALSE; + } + + const Double_t mu_0_guess = 13.5; + const Double_t si_0_guess = 2.5; + const Double_t mu_1_guess = 30.; + const Double_t si_1_guess = si_0_guess + si_0_guess; + // Michele + const Double_t lambda_1cat_guess = 1.00; + const Double_t lambda_1dyn_guess = lambda_1cat_guess/10.; + const Double_t mu_1cat_guess = 50.; + const Double_t mu_1dyn_guess = 17.; + const Double_t si_1cat_guess = si_0_guess + si_0_guess; + const Double_t si_1dyn_guess = si_0_guess + si_0_guess/2.; + // Polya + const Double_t excessPoisson_guess = 0.5; + const Double_t delta1_guess = 8.; + const Double_t delta2_guess = 5.; + const Double_t electronicAmp_guess = gkElectronicAmp; + const Double_t electronicAmp_limit = gkElectronicAmpErr; + + // + // Initialize boundaries and start parameters + // + switch (fFitFunc) + { + + case kEPoisson4: + fSinglePheFit->SetParNames( "#lambda", "#mu_{0}", "#mu_{1}", "#sigma_{0}", "#sigma_{1}","Area"); + fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); + fSinglePheFit->SetParLimits(0,0.,2.); + fSinglePheFit->SetParLimits(1,10.,17.); + fSinglePheFit->SetParLimits(2,17.,50.); + fSinglePheFit->SetParLimits(3,1.,5.); + fSinglePheFit->SetParLimits(4,5.,30.); + fSinglePheFit->SetParLimits(5,norm-(0.5*norm),norm+(0.7*norm)); + break; + case kEPoisson5: + case kEPoisson6: + fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area"); + fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,800.,si_0_guess,500.,norm); + fSinglePheFit->SetParLimits(0,0.,2.); + fSinglePheFit->SetParLimits(1,0.,100.); + fSinglePheFit->SetParLimits(2,300.,1500.); + fSinglePheFit->SetParLimits(3,30.,250.); + fSinglePheFit->SetParLimits(4,100.,1000.); + fSinglePheFit->SetParLimits(5,norm/1.5,norm*1.5); + break; + + case kEPolya: + fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess, + delta1_guess,delta2_guess, + electronicAmp_guess, + fSigmaPedestal, + norm, + fMeanPedestal); + fSinglePheFit->SetParNames("#lambda","b_{tot}", + "#delta_{1}","#delta_{2}", + "amp_{e}","#sigma_{0}", + "Area", "#mu_{0}"); + fSinglePheFit->SetParLimits(0,0.,1.); + fSinglePheFit->SetParLimits(1,0.,1.); + fSinglePheFit->SetParLimits(2,6.,12.); + fSinglePheFit->SetParLimits(3,3.,8.); + fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit, + electronicAmp_guess+electronicAmp_limit); + fSinglePheFit->SetParLimits(5, + fSigmaPedestal-3.*fSigmaPedestalErr, + fSigmaPedestal+3.*fSigmaPedestalErr); + fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1); + fSinglePheFit->SetParLimits(7, + fMeanPedestal-3.*fMeanPedestalErr, + fMeanPedestal+3.*fMeanPedestalErr); + break; + case kEMichele: + fSinglePheFit->SetParNames("#lambda_{cat}","#lambda_{dyn}", + "#mu_{0}","#mu_{1cat}","#mu_{1dyn}", + "#sigma_{0}","#sigma_{1cat}","#sigma_{1dyn}", + "Area"); + fSinglePheFit->SetParameters(lambda_1cat_guess, lambda_1dyn_guess, + mu_0_guess, mu_1cat_guess,mu_1dyn_guess, + si_0_guess, si_1cat_guess,si_1dyn_guess, + norm); + fSinglePheFit->SetParLimits(0,0.01,2.0); + fSinglePheFit->SetParLimits(1,0.,0.5); + fSinglePheFit->SetParLimits(2,10.,16.); + fSinglePheFit->SetParLimits(3,25.,50.); + fSinglePheFit->SetParLimits(4,16.,18.5); + fSinglePheFit->SetParLimits(5,1.,5.); + fSinglePheFit->SetParLimits(6,10.,50.); + fSinglePheFit->SetParLimits(7,5.,10.); + fSinglePheFit->SetParLimits(8,norm/2.,norm*2.5); + break; + + default: + *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; + return kFALSE; + break; + } + + fSinglePheFit->SetRange(rmin,rmax); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// - Retrieve the parameters depending on fFitFunc +// - Retrieve probability, Chisquare and NDF +// +void MHCalibrationChargeBlindPix::ExitFit() +{ + + + // + // Finalize + // + switch (fFitFunc) + { + + case kEPoisson4: + case kEPoisson5: + case kEPoisson6: + case kEPoisson7: + fLambda = fSinglePheFit->GetParameter(0); + fMu0 = fSinglePheFit->GetParameter(1); + fMu1 = fSinglePheFit->GetParameter(2); + fSigma0 = fSinglePheFit->GetParameter(3); + fSigma1 = fSinglePheFit->GetParameter(4); + + fLambdaErr = fSinglePheFit->GetParError(0); + fMu0Err = fSinglePheFit->GetParError(1); + fMu1Err = fSinglePheFit->GetParError(2); + fSigma0Err = fSinglePheFit->GetParError(3); + fSigma1Err = fSinglePheFit->GetParError(4); + break; + case kEPolya: + fLambda = fSinglePheFit->GetParameter(0); + fMu0 = fSinglePheFit->GetParameter(7); + fMu1 = 0.; + fSigma0 = fSinglePheFit->GetParameter(5); + fSigma1 = 0.; + + fLambdaErr = fSinglePheFit->GetParError(0); + fMu0Err = fSinglePheFit->GetParError(7); + fMu1Err = 0.; + fSigma0Err = fSinglePheFit->GetParError(5); + fSigma1Err = 0.; + case kEMichele: + fLambda = fSinglePheFit->GetParameter(0); + fMu0 = fSinglePheFit->GetParameter(2); + fMu1 = fSinglePheFit->GetParameter(3); + fSigma0 = fSinglePheFit->GetParameter(5); + fSigma1 = fSinglePheFit->GetParameter(6); + + fLambdaErr = fSinglePheFit->GetParError(0); + fMu0Err = fSinglePheFit->GetParError(2); + fMu1Err = fSinglePheFit->GetParError(3); + fSigma0Err = fSinglePheFit->GetParError(5); + fSigma1Err = fSinglePheFit->GetParError(6); + break; + default: + break; + } + + fProb = fSinglePheFit->GetProb(); + fChisquare = fSinglePheFit->GetChisquare(); + fNDF = fSinglePheFit->GetNDF(); + + *fLog << all << "Results of the Blind Pixel Fit: " << endl; + *fLog << all << "Chisquare: " << fChisquare << endl; + *fLog << all << "DoF: " << fNDF << endl; + *fLog << all << "Probability: " << fProb << endl; + +} + +// -------------------------------------------------------------------------- +// +// - Executes InitFit() +// - Fits the fHGausHist with fSinglePheFit +// - Executes ExitFit() +// +// The fit result is accepted under condition: +// 1) The results are not nan's +// 2) The NDF is not smaller than fNDFLimit (5) +// 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%) +// 4) at least fNumSinglePheLimit events are in the single Photo-electron peak +// +Bool_t MHCalibrationChargeBlindPix::FitSinglePhe(Option_t *opt) +{ + + if (!InitFit()) + return kFALSE; + + fHGausHist.Fit(fSinglePheFit,opt); + + ExitFit(); + + // + // The fit result is accepted under condition: + // 1) The results are not nan's + // 2) The NDF is not smaller than fNDFLimit (5) + // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%) + // 4) at least fNumSinglePheLimit events are in the single Photo-electron peak + // + if ( TMath::IsNaN(fLambda) + || TMath::IsNaN(fLambdaErr) + || TMath::IsNaN(fProb) + || TMath::IsNaN(fMu0) + || TMath::IsNaN(fMu0Err) + || TMath::IsNaN(fMu1) + || TMath::IsNaN(fMu1Err) + || TMath::IsNaN(fSigma0) + || TMath::IsNaN(fSigma0Err) + || TMath::IsNaN(fSigma1) + || TMath::IsNaN(fSigma1Err) + || fNDF < fNDFLimit + || fProb < fProbLimit ) + return kFALSE; + + const Stat_t entries = fHGausHist.Integral("width"); + const Float_t numSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries; + + if (numSinglePhe < fNumSinglePheLimit) + { + *fLog << warn << "WARNING - Statistics is too low: Only " << numSinglePhe + << " in the Single Photo-Electron peak " << endl; + return kFALSE; + } + else + *fLog << all << numSinglePhe << " in Single Photo-Electron peak " << endl; + + SetSinglePheFitOK(); + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// - Retrieves limits for the fit +// - Fits the fHGausHist with Gauss +// - Retrieves the results to fLambdaCheck and fLambdaCheckErr +// - Sets a flag IsPedestalFitOK() +// +void MHCalibrationChargeBlindPix::FitPedestal (Option_t *opt) +{ + + // Perform the cross-check fitting only the pedestal: + const Axis_t rmin = 0.; + // const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + const Axis_t rmax = fSinglePheCut; + + FitGaus(opt, rmin, rmax); + + const Stat_t entries = fHGausHist.Integral("width"); + const Double_t pedarea = fFGausFit->Integral(0.,fSinglePheCut); + + fLambdaCheck = TMath::Log(entries/pedarea); + // estimate the error by the error of the obtained area from the Gauss-function: + fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0); + + SetPedestalFitOK(IsGausFitOK()); + return; +} + + +// ------------------------------------------------------------------------- +// +// Draw a legend with the fit results +// +void MHCalibrationChargeBlindPix::DrawLegend(Option_t *opt) +{ + + TString option(opt); + + if (!fFitLegend) + { + fFitLegend = new TPaveText(0.05,0.05,0.95,0.95); + fFitLegend->SetLabel(Form("%s%s", "Results of the single PhE Fit (", + (fFitFunc == kEPoisson4) ? "Poisson(k=4))" : + (fFitFunc == kEPoisson5) ? "Poisson(k=5))" : + (fFitFunc == kEPoisson6) ? "Poisson(k=6))" : + (fFitFunc == kEPolya ) ? "Polya(k=4))" : + (fFitFunc == kEMichele ) ? "Michele)" : " none )" )); + fFitLegend->SetTextSize(0.05); + } + else + fFitLegend->Clear(); + + const TString line1 = + Form("Mean: #lambda = %2.2f #pm %2.2f",fLambda,fLambdaErr); + TText *t1 = fFitLegend->AddText(line1.Data()); + t1->SetBit(kCanDelete); + + const TString line6 = + Form("Mean #lambda_{check} = %2.2f #pm %2.2f",fLambdaCheck,fLambdaCheckErr); + TText *t2 = fFitLegend->AddText(line6.Data()); + t2->SetBit(kCanDelete); + + if (option.Contains("datacheck")) + { + if (fLambda + 3.*fLambdaErr < fLambdaCheck - 3.*fLambdaCheckErr + || + fLambda - 3.*fLambdaErr > fLambdaCheck + 3.*fLambdaCheckErr ) + { + TText *t = fFitLegend->AddText("#lambda and #lambda_{check} more than 3#sigma apart!"); + t->SetBit(kCanDelete); + } + } + else + { + + const TString line2 = + Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",fMu0,fMu0Err); + TText *t3 = fFitLegend->AddText(line2.Data()); + t3->SetBit(kCanDelete); + + const TString line3 = + Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",fSigma0,fSigma0Err); + TText *t4 = fFitLegend->AddText(line3.Data()); + t4->SetBit(kCanDelete); + + const TString line4 = + Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",fMu1,fMu1Err); + TText *t5 = fFitLegend->AddText(line4.Data()); + t5->SetBit(kCanDelete); + + const TString line5 = + Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",fSigma1,fSigma1Err); + TText *t6 = fFitLegend->AddText(line5.Data()); + t6->SetBit(kCanDelete); + } + + const TString line7 = + Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChisquare,fNDF); + TText *t7 = fFitLegend->AddText(line7.Data()); + t7->SetBit(kCanDelete); + + const TString line8 = + Form("Probability: %6.4f ",fProb); + TText *t8 = fFitLegend->AddText(line8.Data()); + t8->SetBit(kCanDelete); + + if (IsSinglePheFitOK()) + { + TText *t = fFitLegend->AddText("Result of the Fit: OK"); + t->SetBit(kCanDelete); + } + else + { + TText *t = fFitLegend->AddText("Result of the Fit: NOT OK"); + t->SetBit(kCanDelete); + } + + fFitLegend->SetFillColor(IsSinglePheFitOK() ? 80 : 2); + fFitLegend->Draw(); + + return; +} + + +// ------------------------------------------------------------------------- +// +// Draw the histogram +// +// The following options can be chosen: +// +// "": displays the fHGausHist, the fits, the legend and fASinglePheFADCSlices and fAPedestalFADCSlices +// "all": executes additionally MHGausEvents::Draw(), with option "fourierevents" +// "datacheck" display the fHGausHist, the fits and the legend +// +void MHCalibrationChargeBlindPix::Draw(Option_t *opt) +{ + + TString option(opt); + option.ToLower(); + + Int_t win = 1; + + TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,900, 600); + TVirtualPad *pad = NULL; + + if (option.Contains("all")) + { + option.ReplaceAll("all",""); + oldpad->Divide(2,1); + win = 2; + oldpad->cd(1); + TVirtualPad *newpad = gPad; + pad = newpad; + pad->Divide(2,2); + pad->cd(1); + } + else if (option.Contains("datacheck")) + { + pad = oldpad; + pad->Divide(1,2); + pad->cd(1); + fHGausHist.SetStats(0); + } + else + { + pad = oldpad; + pad->Divide(2,2); + pad->cd(1); + } + + if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow()) + gPad->SetLogy(); + + gPad->SetTicks(); + + fHGausHist.Draw(); + if (fFGausFit ) + { + fFGausFit->SetLineColor(kBlue); + fFGausFit->Draw("same"); + if (!option.Contains("datacheck")) + { + TLine *line = new TLine(fSinglePheCut, 0., fSinglePheCut, 10.); + line->SetBit(kCanDelete); + line->SetLineColor(kBlue); + line->SetLineWidth(3); + line->DrawLine(fSinglePheCut, 0., fSinglePheCut, 2.); + } + } + + if (fSinglePheFit) + { + fSinglePheFit->SetFillStyle(0); + fSinglePheFit->SetLineWidth(3); + fSinglePheFit->SetLineColor(IsSinglePheFitOK() ? kGreen : kRed); + fSinglePheFit->Draw("same"); + } + + pad->cd(2); + DrawLegend(option.Data()); + + if (option.Contains("datacheck")) + return; + + pad->cd(3); + + if (fASinglePheFADCSlices.GetNrows()!=1) + { + if (fHSinglePheFADCSlices) + delete fHSinglePheFADCSlices; + + fHSinglePheFADCSlices = new TH1F(fASinglePheFADCSlices); + fHSinglePheFADCSlices->SetName("SinglePheFADCSlices"); + fHSinglePheFADCSlices->SetTitle(Form("%s%4.1f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut)); + fHSinglePheFADCSlices->SetXTitle("FADC slice number"); + fHSinglePheFADCSlices->SetYTitle("FADC counts"); + const Int_t nbins = fHSinglePheFADCSlices->GetNbinsX(); + TH2D *nulls = new TH2D("Nulls",fHSinglePheFADCSlices->GetTitle(),nbins,0., + fHSinglePheFADCSlices->GetXaxis()->GetBinCenter(nbins), + 100,0.,50.); + nulls->SetDirectory(NULL); + nulls->SetBit(kCanDelete); + nulls->GetXaxis()->SetTitle(fHSinglePheFADCSlices->GetXaxis()->GetTitle()); + nulls->GetYaxis()->SetTitle(fHSinglePheFADCSlices->GetYaxis()->GetTitle()); + nulls->GetXaxis()->CenterTitle(); + nulls->GetYaxis()->CenterTitle(); + nulls->SetStats(0); + nulls->Draw(); + fHSinglePheFADCSlices->Draw("same"); + } + + pad->cd(4); + if (fAPedestalFADCSlices.GetNrows()!=1) + { + + if (fHPedestalFADCSlices) + delete fHPedestalFADCSlices; + + fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices); + fHPedestalFADCSlices->SetName("PedestalFADCSlices"); + fHPedestalFADCSlices->SetTitle(Form("%s%4.1f","Pedestal FADC Slices, Sum < ",fSinglePheCut)); + fHPedestalFADCSlices->SetXTitle("FADC slice number"); + fHPedestalFADCSlices->SetYTitle("FADC counts"); + const Int_t nbins = fHPedestalFADCSlices->GetNbinsX(); + TH2D *nullp = new TH2D("Nullp",fHPedestalFADCSlices->GetTitle(),nbins,0., + fHPedestalFADCSlices->GetXaxis()->GetBinCenter(nbins), + 100,0.,50.); + nullp->SetDirectory(NULL); + nullp->SetBit(kCanDelete); + nullp->GetXaxis()->SetTitle(fHPedestalFADCSlices->GetXaxis()->GetTitle()); + nullp->GetYaxis()->SetTitle(fHPedestalFADCSlices->GetYaxis()->GetTitle()); + nullp->GetXaxis()->CenterTitle(); + nullp->GetYaxis()->CenterTitle(); + nullp->SetStats(0); + nullp->Draw(); + fHPedestalFADCSlices->Draw("same"); + } + + if (win < 2) + return; + + oldpad->cd(2); + MHCalibrationPix::Draw("fourierevents"); +} + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h (revision 4929) @@ -0,0 +1,522 @@ +#ifndef MARS_MHCalibrationChargeBlindPix +#define MARS_MHCalibrationChargeBlindPix + + +#ifndef MARS_MHCalibrationPix +#include "MHCalibrationPix.h" +#endif + +#ifndef ROOT_TF1 +#include +#endif + +#ifndef ROOT_TVector +#include +#endif + +class TH1F; +class TPaveText; +class TText; +class MRawEvtData; +class MRawEvtPixelIter; +class MCalibrationChargeBlindPix; +class MExtractBlindPixel; +class MExtractedSignalBlindPixel; + +class MHCalibrationChargeBlindPix : public MHCalibrationPix +{ +private: + + static const Int_t fgChargeNbins; //! Default for fNBins (now set to: 5300 ) + static const Axis_t fgChargeFirst; //! Default for fFirst (now set to: -100.5 ) + static const Axis_t fgChargeLast; //! Default for fLast (now set to: 5199.5 ) + static const Float_t fgSinglePheCut; //! Default for fSinglePheCut (now set to: 200 ) + static const Float_t fgNumSinglePheLimit; //! Default for fNumSinglePheLimit (now set to: 50) + static const Float_t gkSignalInitializer; //! Signal initializer (-9999.) + + static const Double_t gkElectronicAmp; // Electronic Amplification after the PMT (in FADC counts/N_e) + static const Double_t gkElectronicAmpErr; // Error of the electronic amplification + + Float_t fSinglePheCut; // Value of summed FADC slices upon which event considered as single-phe + Float_t fNumSinglePheLimit; // Minimum number of single-phe events + + MCalibrationChargeBlindPix *fBlindPix; //! Storage container results + MExtractedSignalBlindPixel *fSignal; //! Storage container extracted signal + MRawEvtData *fRawEvt; //! Storage container raw data + + TVector fASinglePheFADCSlices; // Averaged FADC slice entries supposed single-phe events + TVector fAPedestalFADCSlices; // Averaged FADC slice entries supposed pedestal events + + TF1 *fSinglePheFit; // Single Phe Fit (Gaussians convoluted with Poisson) + + UInt_t fNumSinglePhes; // Number of entries in fASinglePheFADCSlices + UInt_t fNumPedestals; // Number of entries in fAPedestalFADCSlices + + Double_t fLambda; // Poisson mean from Single-phe fit + Double_t fLambdaCheck; // Poisson mean from Pedestal fit alone + Double_t fMu0; // Mean of the pedestal + Double_t fMu1; // Mean of single-phe peak + Double_t fSigma0; // Sigma of the pedestal + Double_t fSigma1; // Sigma of single-phe peak + Double_t fLambdaErr; // Error of Poisson mean from Single-phe fit + Double_t fLambdaCheckErr; // Error of Poisson mean from Pedestal fit alone + Double_t fMu0Err; // Error of Mean of the pedestal + Double_t fMu1Err; // Error of Mean of single-phe peak + Double_t fSigma0Err; // Error of Sigma of the pedestal + Double_t fSigma1Err; // Error of Sigma of single-phe peak + Double_t fChisquare; // Chisquare of single-phe fit + Int_t fNDF; // Ndof of single-phe fit + Double_t fProb; // Probability of singleo-phe fit + Double_t fMeanPedestal; // Mean pedestal from pedestal run + Double_t fSigmaPedestal; // Sigma pedestal from pedestal run + Double_t fMeanPedestalErr; // Error of Mean pedestal from pedestal run + Double_t fSigmaPedestalErr; // Error of Sigma pedestal from pedestal run + + Byte_t fFlags; // Bit-field for the flags + enum { kSinglePheFitOK, kPedestalFitOK }; // Possible bits to be set + + TPaveText *fFitLegend; //! Some legend to display the fit results + TH1F *fHSinglePheFADCSlices; //! A histogram created and deleted only in Draw() + TH1F *fHPedestalFADCSlices; //! A histogram created and deleted only in Draw() + + // Fill histos + void FillSinglePheFADCSlices(const MRawEvtPixelIter &iter); + void FillPedestalFADCSlices( const MRawEvtPixelIter &iter); + + // Fit + Bool_t InitFit(); + void ExitFit(); + +public: + + MHCalibrationChargeBlindPix(const char *name=NULL, const char *title=NULL); + ~MHCalibrationChargeBlindPix(); + + void Clear(Option_t *o=""); + void Reset(); + +// TObject *Clone(const char *) const; + + Bool_t SetupFill(const MParList *pList); + Bool_t ReInit ( MParList *pList); + Bool_t Fill (const MParContainer *par, const Stat_t w=1); + Bool_t Finalize(); + + // Getters + const Double_t GetLambda () const { return fLambda; } + const Double_t GetLambdaCheck () const { return fLambdaCheck; } + const Double_t GetMu0 () const { return fMu0; } + const Double_t GetMu1 () const { return fMu1; } + const Double_t GetSigma0 () const { return fSigma0; } + const Double_t GetSigma1 () const { return fSigma1; } + const Double_t GetLambdaErr () const { return fLambdaErr; } + const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; } + const Double_t GetMu0Err () const { return fMu0Err; } + const Double_t GetMu1Err () const { return fMu1Err; } + const Double_t GetSigma0Err () const { return fSigma0Err; } + const Double_t GetSigma1Err () const { return fSigma1Err; } + const Float_t GetSinglePheCut () const { return fSinglePheCut; } + + TVector &GetASinglePheFADCSlices() { return fASinglePheFADCSlices; } + const TVector &GetASinglePheFADCSlices() const { return fASinglePheFADCSlices; } + + TVector &GetAPedestalFADCSlices() { return fAPedestalFADCSlices; } + const TVector &GetAPedestalFADCSlices() const { return fAPedestalFADCSlices; } + + const Bool_t IsSinglePheFitOK() const; + const Bool_t IsPedestalFitOK() const; + + // Setters + void SetCalibrationChargeBlindPix ( MCalibrationChargeBlindPix *pix) { fBlindPix = pix; } + void SetSinglePheCut ( const Float_t cut =fgSinglePheCut ) { fSinglePheCut = cut; } + void SetNumSinglePheLimit ( const Float_t lim =fgNumSinglePheLimit ) { fNumSinglePheLimit = lim; } + + void SetMeanPedestal ( const Float_t f ) { fMeanPedestal = f; } + void SetMeanPedestalErr ( const Float_t f ) { fMeanPedestalErr = f; } + void SetSigmaPedestal ( const Float_t f ) { fSigmaPedestal = f; } + void SetSigmaPedestalErr ( const Float_t f ) { fSigmaPedestalErr = f; } + + void SetSinglePheFitOK ( const Bool_t b=kTRUE); + void SetPedestalFitOK ( const Bool_t b=kTRUE); + + // Draws + void Draw(Option_t *opt=""); + +private: + void DrawLegend(Option_t *opt=""); + + // Fits +public: + enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele }; // The possible fit functions + +private: + FitFunc_t fFitFunc; + +public: + Bool_t FitSinglePhe (Option_t *opt="RL0+Q"); + void FitPedestal (Option_t *opt="RL0+Q"); + + void ChangeFitFunc(const FitFunc_t func) { fFitFunc = func; } + + // Simulation + Bool_t SimulateSinglePhe(const Double_t lambda, + const Double_t mu0, const Double_t mu1, + const Double_t sigma0, const Double_t sigma1); + +private: + + inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par) + { + + Double_t lambda1cat = par[0]; + Double_t lambda1dyn = par[1]; + Double_t mu0 = par[2]; + Double_t mu1cat = par[3]; + Double_t mu1dyn = par[4]; + Double_t sigma0 = par[5]; + Double_t sigma1cat = par[6]; + Double_t sigma1dyn = par[7]; + + Double_t sumcat = 0.; + Double_t sumdyn = 0.; + Double_t arg = 0.; + + if (lambda1cat < lambda1dyn) + return FLT_MAX; + + if (mu1cat < mu0) + return FLT_MAX; + + if (mu1dyn < mu0) + return FLT_MAX; + + if (mu1cat < mu1dyn) + return FLT_MAX; + + if (sigma0 < 0.0001) + return FLT_MAX; + + if (sigma1cat < sigma0) + return FLT_MAX; + + if (sigma1dyn < sigma0) + return FLT_MAX; + + Double_t mu2cat = (2.*mu1cat)-mu0; + Double_t mu2dyn = (2.*mu1dyn)-mu0; + Double_t mu3cat = (3.*mu1cat)-(2.*mu0); + Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0); + + Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0)); + Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0)); + Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0)); + Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0)); + + Double_t lambda2cat = lambda1cat*lambda1cat; + Double_t lambda2dyn = lambda1dyn*lambda1dyn; + Double_t lambda3cat = lambda2cat*lambda1cat; + Double_t lambda3dyn = lambda2dyn*lambda1dyn; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sumcat = TMath::Exp(-0.5*arg*arg)/sigma0; + sumdyn = sumcat; + + // k=1cat: + arg = (x[0] - mu1cat)/sigma1cat; + sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat; + // k=1dyn: + arg = (x[0] - mu1dyn)/sigma1dyn; + sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn; + + // k=2cat: + arg = (x[0] - mu2cat)/sigma2cat; + sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat; + // k=2dyn: + arg = (x[0] - mu2dyn)/sigma2dyn; + sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn; + + + // k=3cat: + arg = (x[0] - mu3cat)/sigma3cat; + sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat; + // k=3dyn: + arg = (x[0] - mu3dyn)/sigma3dyn; + sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn; + + sumcat = TMath::Exp(-1.*lambda1cat)*sumcat; + sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn; + + return par[8]*(sumcat+sumdyn)/2.; + + } + + inline static Double_t fPoissonKto4(Double_t *x, Double_t *par) + { + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return FLT_MAX; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma0 < 0.0001) + return FLT_MAX; + + if (sigma1 < sigma0) + return FLT_MAX; + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + return TMath::Exp(-1.*lambda)*par[5]*sum; + + } + + + inline static Double_t fPoissonKto5(Double_t *x, Double_t *par) + { + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return FLT_MAX; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma0 < 0.0001) + return FLT_MAX; + + if (sigma1 < sigma0) + return FLT_MAX; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + return TMath::Exp(-1.*lambda)*par[5]*sum; + + } + + + inline static Double_t fPoissonKto6(Double_t *x, Double_t *par) + { + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return FLT_MAX; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma0 < 0.0001) + return FLT_MAX; + + if (sigma1 < sigma0) + return FLT_MAX; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + Double_t mu6 = (6.*mu1)-(5.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + Double_t lambda6 = lambda5*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + // k=6: + arg = (x[0] - mu6)/sigma6; + sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6; + + return TMath::Exp(-1.*lambda)*par[5]*sum; + + } + + inline static Double_t fPolya(Double_t *x, Double_t *par) + { + + const Double_t QEcat = 0.247; // mean quantum efficiency + const Double_t sqrt2 = 1.4142135623731; + const Double_t sqrt3 = 1.7320508075689; + const Double_t sqrt4 = 2.; + + const Double_t lambda = par[0]; // mean number of photons + + const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution + const Double_t delta1 = par[2]; // amplification first dynode + const Double_t delta2 = par[3]; // amplification subsequent dynodes + + const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges + + const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain + const Double_t A = 1. + excessPoisson - QEcat + + 1./delta1 + + 1./delta1/delta2 + + 1./delta1/delta2/delta2; // variance contributions from PMT and QE + + const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion + + const Double_t mu0 = par[7]; // pedestal + const Double_t mu1 = totAmpl; // single phe position + const Double_t mu2 = 2*totAmpl; // double phe position + const Double_t mu3 = 3*totAmpl; // triple phe position + const Double_t mu4 = 4*totAmpl; // quadruple phe position + + const Double_t sigma0 = par[5]; + const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A); + const Double_t sigma2 = sqrt2*sigma1; + const Double_t sigma3 = sqrt3*sigma1; + const Double_t sigma4 = sqrt4*sigma1; + + const Double_t lambda2 = lambda*lambda; + const Double_t lambda3 = lambda2*lambda; + const Double_t lambda4 = lambda3*lambda; + + //-- calculate the area---- + Double_t arg = (x[0] - mu0)/sigma0; + Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + return TMath::Exp(-1.*lambda)*par[6]*sum; + } + + ClassDef(MHCalibrationChargeBlindPix, 1) // Histogram class for Charge Blind Pixel Calibration +}; + +#endif /* MARS_MHCalibrationChargeBlindPix */ Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeCam.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeCam.cc (revision 4929) @@ -0,0 +1,1339 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationChargeCam +// +// Fills the extracted signals of MExtractedSignalCam into the MHCalibrationPix-classes +// MHCalibrationChargeHiGainPix and MHCalibrationChargeLoGainPix for every: +// +// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray and +// MHCalibrationCam::fLoGainArray +// +// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), +// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and +// MHCalibrationCam::fAverageLoGainAreas +// +// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), +// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors and +// MHCalibrationCam::fAverageLoGainSectors +// +// Every signal is taken from MExtractedSignalCam and 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. +// +// Additionally, the (FADC slice) position of the maximum is stored in an Absolute +// Arrival Time histogram. This histogram serves for a rough cross-check if the +// signal does not lie at or outside the edges of the extraction window. +// +// The Charge 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 +// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit()) +// In case this does not make the fit valid, the histogram means and RMS's are +// taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) or +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas +// from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup) +// +// Unless more than fNumHiGainSaturationLimit (default: 1%) of the overall FADC +// slices show saturation, the following flag is set: +// - MCalibrationChargePix::SetHiGainSaturation(); +// In that case, the calibration constants are derived from the low-gain results. +// +// If more than fNumLoGainSaturationLimit (default: 1%) of the overall +// low-gain FADC slices saturate, the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnsuitableRun ) +// +// The class also fills arrays with the signal vs. event number, creates a fourier +// spectrum 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 ) or +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// This same procedure is performed for the average pixels. +// +// The following results are written into MCalibrationChargeCam: +// +// - MCalibrationPix::SetHiGainSaturation() +// - MCalibrationPix::SetHiGainMean() +// - MCalibrationPix::SetHiGainMeanErr() +// - MCalibrationPix::SetHiGainSigma() +// - MCalibrationPix::SetHiGainSigmaErr() +// - MCalibrationPix::SetHiGainProb() +// - MCalibrationPix::SetHiGainNumPickup() +// +// - MCalibrationPix::SetLoGainMean() +// - MCalibrationPix::SetLoGainMeanErr() +// - MCalibrationPix::SetLoGainSigma() +// - MCalibrationPix::SetLoGainSigmaErr() +// - MCalibrationPix::SetLoGainProb() +// - MCalibrationPix::SetLoGainNumPickup() +// +// - MCalibrationChargePix::SetAbsTimeMean() +// - MCalibrationChargePix::SetAbsTimeRms() +// +// 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 "MHCalibrationChargeCam.h" +#include "MHCalibrationCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MHCalibrationChargePix.h" +#include "MHCalibrationPix.h" + +#include "MCalibrationIntensityCam.h" +#include "MCalibrationChargeCam.h" +#include "MCalibrationChargePix.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" + +#include "MRawEvtData.h" +#include "MRawRunHeader.h" +#include "MRawEvtPixelIter.h" + +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" + +#include "MArrayI.h" +#include "MArrayD.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +ClassImp(MHCalibrationChargeCam); + +using namespace std; + +const Axis_t MHCalibrationChargeCam::fgChargeHiGainFirst = -100.5; +const Axis_t MHCalibrationChargeCam::fgChargeHiGainLast = 999.5; +const Axis_t MHCalibrationChargeCam::fgChargeLoGainFirst = -150.5; +const Axis_t MHCalibrationChargeCam::fgChargeLoGainLast = 499.5; + +const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.01; +const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005; +const Float_t MHCalibrationChargeCam::fgTimeLowerLimit = 1.; +const Float_t MHCalibrationChargeCam::fgTimeUpperLimit = 2.; +// 1Led Green, 1 LED blue, 5 LEDs blue, 10 LEDs blue, 10 LEDs UV, CT1, 5Leds Green +const Float_t MHCalibrationChargeCam::gkHiGainInnerRefLines[7] = { 245., 323. , 1065., 1467., 180., 211. , 533.5}; +const Float_t MHCalibrationChargeCam::gkHiGainOuterRefLines[7] = { 217., 307.5, 932. , 1405., 167., 183.5, 405.5}; +const Float_t MHCalibrationChargeCam::gkLoGainInnerRefLines[7] = { 20.8, 28.0 , 121. , 200.2, 16.5, 13.5 , 41.7 }; +const Float_t MHCalibrationChargeCam::gkLoGainOuterRefLines[7] = { 18.9, 26.0 , 108.3, 198. , 14.0, 11. , 42. }; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - all pointers to NULL +// +// Initializes: +// - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit +// - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit +// - fTimeLowerLimit to fgTimeLowerLimit +// - fTimeUpperLimit to fgTimeUpperLimit +// +MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title) + : fRawEvt(NULL) +{ + + fName = name ? name : "MHCalibrationChargeCam"; + fTitle = title ? title : "Class to fill the calibration histograms "; + + SetNumHiGainSaturationLimit(fgNumHiGainSaturationLimit); + SetNumLoGainSaturationLimit(fgNumLoGainSaturationLimit); + SetTimeLowerLimit(); + SetTimeUpperLimit(); + +} + +// -------------------------------------------------------------------------- +// +// Gets the pointers to: +// - MRawEvtData +// +Bool_t MHCalibrationChargeCam::SetupHists(const MParList *pList) +{ + + fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); + if (!fRawEvt) + { + *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// - MExtractedSignalCam +// - MCalibrationChargeCam +// - MBadPixelsCam +// +// Initializes the number of used FADC slices from MExtractedSignalCam +// into MCalibrationChargeCam and test for changes in that variable +// +// Initializes, if empty to MGeomCam::GetNumPixels(): +// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray +// - MHCalibrationCam::fHiGainOverflow, MHCalibrationCam::fLoGainOverflow +// +// 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::fLoGainArray +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +// Sets Titles and Names for the Charge Histograms: +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +Bool_t MHCalibrationChargeCam::ReInitHists(MParList *pList) +{ + + MExtractedSignalCam *signal = + (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam")); + if (!signal) + { + *fLog << err << "MExtractedSignalCam not found... abort." << endl; + return kFALSE; + } + + fIntensCam = (MCalibrationIntensityCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityChargeCam")); + if (fIntensCam) + *fLog << inf << "Found MCalibrationIntensityChargeCam ... " << endl; + else + { + fCam = (MCalibrationCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam")); + if (!fCam) + { + fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationChargeCam")); + if (!fCam) + { + *fLog << err << "Cannot find nor create MCalibrationChargeCam ... abort." << endl; + return kFALSE; + } + fCam->Init(*fGeom); + } + } + + fFirstHiGain = signal->GetFirstUsedSliceHiGain(); + fLastHiGain = signal->GetLastUsedSliceHiGain(); + fFirstLoGain = signal->GetFirstUsedSliceLoGain(); + fLastLoGain = signal->GetLastUsedSliceLoGain(); + + /* + const Float_t numhigain = signal->GetNumUsedHiGainFADCSlices(); + const Float_t numlogain = signal->GetNumUsedLoGainFADCSlices(); + + + if (fCam) + { + if (fCam->GetNumHiGainFADCSlices() == 0.) + fCam->SetNumHiGainFADCSlices ( numhigain ); + else if (fCam->GetNumHiGainFADCSlices() != numhigain) + { + *fLog << err << GetDescriptor() + << ": Number of High Gain FADC extraction slices has changed, abort..." << endl; + return kFALSE; + } + + if (fCam->GetNumLoGainFADCSlices() == 0.) + fCam->SetNumLoGainFADCSlices ( numlogain ); + else if (fCam->GetNumLoGainFADCSlices() != numlogain) + { + *fLog << err << GetDescriptor() + << ": Number of Low Gain FADC extraction slices has changes, abort..." << endl; + return kFALSE; + } + } + */ + + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nsectors = fGeom->GetNumSectors(); + const Int_t nareas = fGeom->GetNumAreas(); + + const Int_t higainsamples = fRunHeader->GetNumSamplesHiGain(); + const Int_t logainsamples = fRunHeader->GetNumSamplesLoGain(); + + if (fHiGainArray->GetEntries()==0) + { + fHiGainArray->Expand(npixels); + for (Int_t i=0; iSetName ("HCalibrationChargeHiGainPix"); + pix.GetHGausHist()->SetTitle("Summed Hi Gain FADC slices Pixel "); + + pix.GetHAbsTime()->SetName ("HAbsTimeHiGainPix"); + pix.GetHAbsTime()->SetTitle("Absolute Arrival Times Hi Gain Pixel "); + + InitHists((*this)[i],(*fBadPixels)[i],i); + } + } + + if (fLoGainArray->GetEntries()==0 && IsLoGain()) + { + fLoGainArray->Expand(npixels); + + for (Int_t i=0; iSetName ("HCalibrationChargeLoGainPix"); + pix.GetHGausHist()->SetTitle("Summed Lo Gain FADC slices Pixel "); + + pix.GetHAbsTime()->SetName ("HAbsTimeLoGainPix"); + pix.GetHAbsTime()->SetTitle("Absolute Arrival Times Lo Gain Pixel "); + + // + // Adapt the range for the case, the intense blue is used: + // FIXME: this is a nasty workaround, but for the moment necessary + // in order to avoid default memory space. + // + if (fGeom->InheritsFrom("MGeomCamMagic")) + { + if ( fColor == MCalibrationCam::kBLUE) + { + pix.SetLast(999.5); + pix.SetNbins((Int_t)(999.5-fgChargeLoGainFirst)); + } + } + InitHists(pix,(*fBadPixels)[i],i); + } + } + + if (fAverageHiGainAreas->GetEntries()==0) + { + fAverageHiGainAreas->Expand(nareas); + + for (Int_t j=0; jSetName ("HCalibrationChargeHiGainAreaIdx"); + pix.GetHGausHist()->SetTitle("Signal averaged on event-by-event basis High Gain Area Idx "); + + pix.GetHAbsTime()->SetName ("HAbsTimeHiGainAreaIdx"); + pix.GetHAbsTime()->SetTitle("Absolute Arrival Times average Hi Gain Area Idx "); + + if (fGeom->InheritsFrom("MGeomCamMagic")) + { + pix.GetHGausHist()->SetTitle(Form("%s%s%s","Signal averaged on event-by-event basis ", + j==0 ? "Inner Pixels " : "Outer Pixels ","High Gain Runs: ")); + pix.InitBins(); + pix.SetEventFrequency(fPulserFrequency); + } + else + { + InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j); + } + } + } + + + if (fAverageLoGainAreas->GetEntries()==0 && IsLoGain()) + { + fAverageLoGainAreas->Expand(nareas); + + for (Int_t j=0; jSetName ("HChargeLoGainAreaIdx"); + pix.GetHGausHist()->SetTitle("Signal averaged on event-by-event basis Low Gain Area Idx "); + + pix.GetHAbsTime()->SetName ("HAbsTimeLoGainAreaIdx"); + pix.GetHAbsTime()->SetTitle("Absolute Arrival Times average Lo Gain Area Idx "); + + // + // Adapt the range for the case, the intense blue is used: + // FIXME: this is a nasty workaround, but for the moment necessary + // in order to avoid default memory space. + // + if (fGeom->InheritsFrom("MGeomCamMagic")) + { + if ( fColor == MCalibrationCam::kBLUE) + { + pix.SetFirst(-10.5); + pix.SetLast(999.5); + pix.SetNbins(3030); + } + } + + if (fGeom->InheritsFrom("MGeomCamMagic")) + { + pix.GetHGausHist()->SetTitle(Form("%s%s%s","Signal averaged on event-by-event basis ", + j==0 ? "Inner Pixels " : "Outer Pixels ","Low Gain Runs: ")); + pix.InitBins(); + pix.SetEventFrequency(fPulserFrequency); + } + else + { + InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j); + } + } + } + + if (fAverageHiGainSectors->GetEntries()==0) + { + fAverageHiGainSectors->Expand(nsectors); + + for (Int_t j=0; jSetName ("HCalibrationChargeHiGainPix"); + pix.GetHGausHist()->SetTitle("Signals averaged on event-by-event basis HiGain Sector "); + + pix.GetHAbsTime()->SetName ("HAbsTimeHiGainPix"); + pix.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Sector "); + + InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j); + } + } + + if (fAverageLoGainSectors->GetEntries()==0 && IsLoGain()) + { + fAverageLoGainSectors->Expand(nsectors); + + for (Int_t j=0; jSetName ("HCalibrationChargeLoGainPix"); + pix.GetHGausHist()->SetTitle("Signals averaged on event-by-event basis LoGain Sector "); + + pix.GetHAbsTime()->SetName ("HAbsTimeLoGainPix"); + pix.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Sector "); + + // + // Adapt the range for the case, the intense blue is used: + // FIXME: this is a nasty workaround, but for the moment necessary + // in order to avoid default memory space. + // + if (fGeom->InheritsFrom("MGeomCamMagic")) + { + if ( fColor == MCalibrationCam::kBLUE) + { + pix.SetFirst(-10.5); + pix.SetLast(999.5); + pix.SetNbins(3030); + } + } + + InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j); + } + } + + fSumhiarea .Set(nareas); + fSumloarea .Set(nareas); + fTimehiarea .Set(nareas); + fTimeloarea .Set(nareas); + fSumhisector.Set(nsectors); + fSumlosector.Set(nsectors); + fTimehisector.Set(nsectors); + fTimelosector.Set(nsectors); + + fSathiarea .Set(nareas); + fSatloarea .Set(nareas); + fSathisector.Set(nsectors); + fSatlosector.Set(nsectors); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Retrieves from MExtractedSignalCam: +// - first used LoGain FADC slice +// +// Retrieves from MGeomCam: +// - number of pixels +// - number of pixel areas +// - number of sectors +// +// For all TObjArray's (including the averaged ones), the following steps are performed: +// +// 1) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with: +// - MExtractedSignalPix::GetExtractedSignalHiGain(); +// - MExtractedSignalPix::GetExtractedSignalLoGain(); +// +// 2) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with: +// - MExtractedSignalPix::GetNumHiGainSaturated(); +// - MExtractedSignalPix::GetNumLoGainSaturated(); +// +// 3) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with: +// - MRawEvtPixelIter::GetIdxMaxHiGainSample(); +// - MRawEvtPixelIter::GetIdxMaxLoGainSample(first slice); +// +Bool_t MHCalibrationChargeCam::FillHists(const MParContainer *par, const Stat_t w) +{ + + MExtractedSignalCam *signal = (MExtractedSignalCam*)par; + if (!signal) + { + *fLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl; + return kFALSE; + } + + const UInt_t npixels = fGeom->GetNumPixels(); + const UInt_t nareas = fGeom->GetNumAreas(); + const UInt_t nsectors = fGeom->GetNumSectors(); + const UInt_t lofirst = signal->GetFirstUsedSliceLoGain(); + + fSumhiarea .Reset(); + fSumloarea .Reset(); + fTimehiarea .Reset(); + fTimeloarea .Reset(); + fSumhisector.Reset(); + fSumlosector.Reset(); + fTimehisector.Reset(); + fTimelosector.Reset(); + + fSathiarea .Reset(); + fSatloarea .Reset(); + fSathisector.Reset(); + fSatlosector.Reset(); + + for (UInt_t i=0; iGetSize(); i++) + { + + MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i]; + + if (histhi.IsExcluded()) + continue; + + MCalibrationChargePix &pix = fIntensCam + ? (MCalibrationChargePix&)(*fIntensCam)[i] + : (MCalibrationChargePix&)(*fCam)[i]; + + if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + pix.SetHiGainSaturation(); + if (IsOscillations()) + histhi.CreateFourierSpectrum(); + continue; + } + + MBadPixelsPix &bad = (*fBadPixels)[i]; + + Stat_t overflow = histhi.GetHGausHist()->GetBinContent(histhi.GetHGausHist()->GetNbinsX()+1); + if (overflow > 0.1) + { + *fLog << warn << GetDescriptor() + << ": HiGain Histogram Overflow occurred " << overflow + << " times in pixel: " << i << " (without saturation!) " << endl; + bad.SetUncalibrated( MBadPixelsPix::kHiGainOverFlow ); + } + + overflow = histhi.GetHGausHist()->GetBinContent(0); + if (overflow > 0.1) + { + *fLog << warn << GetDescriptor() + << ": HiGain Histogram Underflow occurred " << overflow + << " times in pixel: " << i << " (without saturation!) " << endl; + bad.SetUncalibrated( MBadPixelsPix::kHiGainOverFlow ); + } + + FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain); + } + + if (IsLoGain()) + for (Int_t i=0; iGetSize(); i++) + { + + MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i); + MBadPixelsPix &bad = (*fBadPixels)[i]; + + if (histlo.IsExcluded()) + continue; + + if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) + { + *fLog << warn << "Saturated Lo Gain histogram in pixel: " << i << endl; + bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); + if (IsOscillations()) + histlo.CreateFourierSpectrum(); + continue; + } + + Stat_t overflow = histlo.GetHGausHist()->GetBinContent(histlo.GetHGausHist()->GetNbinsX()+1); + if (overflow > 0.1) + { + *fLog << warn << GetDescriptor() + << ": Lo-Gain Histogram Overflow occurred " << overflow + << " times in pixel: " << i << " (without saturation!) " << endl; + bad.SetUncalibrated( MBadPixelsPix::kLoGainOverFlow ); + } + + overflow = histlo.GetHGausHist()->GetBinContent(0); + if (overflow > 0.1) + { + *fLog << warn << GetDescriptor() + << ": Lo-Gain Histogram Underflow occurred " << overflow + << " times in pixel: " << i << " (without saturation!) " << endl; + bad.SetUncalibrated( MBadPixelsPix::kLoGainOverFlow ); + } + + MCalibrationChargePix &pix = fIntensCam + ? (MCalibrationChargePix&)(*fIntensCam)[i] + : (MCalibrationChargePix&)(*fCam)[i]; + + if (pix.IsHiGainSaturation()) + FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain); + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainArea(j); + MCalibrationChargePix &pix = fIntensCam + ? (MCalibrationChargePix&)fIntensCam->GetAverageArea(j) + : (MCalibrationChargePix&)fCam->GetAverageArea(j); + + if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + pix.SetHiGainSaturation(); + if (IsOscillations()) + histhi.CreateFourierSpectrum(); + continue; + } + + MBadPixelsPix &bad = fIntensCam + ? fIntensCam->GetAverageBadArea(j) + : fCam->GetAverageBadArea(j); + + FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain); + } + + if (IsLoGain()) + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainArea(j); + + if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) + { + *fLog << warn << "Saturated Lo Gain histogram in area idx: " << j << endl; + histlo.CreateFourierSpectrum(); + continue; + } + + MCalibrationChargePix &pix = fIntensCam + ? (MCalibrationChargePix&)fIntensCam->GetAverageArea(j) + : (MCalibrationChargePix&)fCam->GetAverageArea(j) ; + + if (pix.IsHiGainSaturation()) + { + MBadPixelsPix &bad = fIntensCam + ? fIntensCam->GetAverageBadArea(j) + : fCam->GetAverageBadArea(j); + FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain); + } + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainSector(j); + MCalibrationChargePix &pix = fIntensCam + ? (MCalibrationChargePix&)fIntensCam->GetAverageSector(j) + : (MCalibrationChargePix&)fCam->GetAverageSector(j); + + if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + pix.SetHiGainSaturation(); + if (IsOscillations()) + histhi.CreateFourierSpectrum(); + continue; + } + + MBadPixelsPix &bad = fIntensCam + ? fIntensCam->GetAverageBadSector(j) + : fCam->GetAverageBadSector(j); + + FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain); + } + + if (IsLoGain()) + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j); + MBadPixelsPix &bad = fIntensCam + ? fIntensCam->GetAverageBadSector(j) + : fCam->GetAverageBadSector(j); + + if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) + { + *fLog << warn << "Saturated Lo Gain histogram in sector: " << j << endl; + bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); + if (IsOscillations()) + histlo.CreateFourierSpectrum(); + continue; + } + + MCalibrationChargePix &pix = fIntensCam + ? (MCalibrationChargePix&)fIntensCam->GetAverageSector(j) + : (MCalibrationChargePix&)fCam->GetAverageSector(j); + + if (pix.IsHiGainSaturation()) + FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain); + } + + // + // Perform the fitting for the High Gain (done in MHCalibrationCam) + // + FitHiGainArrays(fIntensCam ? (MCalibrationCam&)(*fIntensCam->GetCam()) : (MCalibrationCam&)(*fCam), + *fBadPixels, + MBadPixelsPix::kHiGainNotFitted, + MBadPixelsPix::kHiGainOscillating); + + // + // Perform the fitting for the Low Gain (done in MHCalibrationCam) + // + if (IsLoGain()) + FitLoGainArrays(fIntensCam ? (MCalibrationCam&)(*fIntensCam->GetCam()) : (MCalibrationCam&)(*fCam), + *fBadPixels, + MBadPixelsPix::kLoGainNotFitted, + MBadPixelsPix::kLoGainOscillating); + + return kTRUE; +} + +// -------------------------------------------------------------------------------- +// +// Fill the absolute time results into MCalibrationChargePix +// +// Check absolute time validity: +// - Mean arrival time is at least fTimeLowerLimit slices from the lower edge +// - Mean arrival time is at least fUpperLimit slices from the upper edge +// +void MHCalibrationChargeCam::FinalizeAbsTimes(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad, + Byte_t first, Byte_t last) +{ + + const Float_t mean = hist.GetAbsTimeMean(); + const Float_t rms = hist.GetAbsTimeRms(); + + pix.SetAbsTimeMean ( mean ); + pix.SetAbsTimeRms ( rms ); + + const Float_t lowerlimit = (Float_t)first + fTimeLowerLimit; + const Float_t upperlimit = (Float_t)last + fTimeUpperLimit; + + if ( mean < lowerlimit) + { + *fLog << warn << GetDescriptor() + << Form("%s%3.1f%s%2.1f%s%4i",": Mean ArrivalTime: ",mean," smaller than ",fTimeLowerLimit, + " FADC slices from lower edge in pixel ",hist.GetPixId()) << endl; + bad.SetUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin ); + } + + if ( mean > upperlimit ) + { + *fLog << warn << GetDescriptor() + << Form("%s%3.1f%s%2.1f%s%4i",": Mean ArrivalTime: ",mean," greater than ",fTimeUpperLimit, + " FADC slices from upper edge in pixel ",hist.GetPixId()) << endl; + bad.SetUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins ); + } +} + +// -------------------------------------------------------------------------- +// +// Sets all pixels to MBadPixelsPix::kUnsuitableRun, if following flags are set: +// - MBadPixelsPix::kLoGainSaturation +// +// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set: +// - if MBadPixelsPix::kHiGainNotFitted and !MCalibrationPix::IsHiGainSaturation() +// - if MBadPixelsPix::kHiGainOscillating and !MCalibrationPix::IsHiGainSaturation() +// - if MBadPixelsPix::kLoGainNotFitted and MCalibrationPix::IsLoGainSaturation() +// - if MBadPixelsPix::kLoGainOscillating and MCalibrationPix::IsLoGainSaturation() +// +void MHCalibrationChargeCam::FinalizeBadPixels() +{ + + for (Int_t i=0; iGetSize(); i++) + { + + MBadPixelsPix &bad = (*fBadPixels)[i]; + MCalibrationPix &pix = fIntensCam ? (*fIntensCam)[i] : (*fCam)[i]; + + if (bad.IsUncalibrated( MBadPixelsPix::kHiGainNotFitted )) + if (!pix.IsHiGainSaturation()) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted )) + if (pix.IsHiGainSaturation()) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation )) + bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); + + if (IsOscillations()) + { + if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOscillating )) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating )) + if (pix.IsHiGainSaturation()) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + } + } +} + +// -------------------------------------------------------------------------- +// +// Dummy, needed by MCamEvent +// +Bool_t MHCalibrationChargeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calls MHCalibrationPix::DrawClone() for pixel idx +// +void MHCalibrationChargeCam::DrawPixelContent(Int_t idx) const +{ + (*this)[idx].DrawClone(); +} + + +// ----------------------------------------------------------------------------- +// +// Default draw: +// +// Displays the averaged areas, both High Gain and Low Gain +// +// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options +// +void MHCalibrationChargeCam::Draw(const Option_t *opt) +{ + + const Int_t nareas = fAverageHiGainAreas->GetEntries(); + if (nareas == 0) + return; + + TString option(opt); + option.ToLower(); + + if (!option.Contains("datacheck")) + { + MHCalibrationCam::Draw(opt); + return; + } + + // + // From here on , the datacheck - Draw + // + TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this); + pad->SetBorderMode(0); + pad->Divide(1,nareas); + + // + // Loop over inner and outer pixels + // + for (Int_t i=0; icd(i+1); + + MHCalibrationChargePix &hipix = (MHCalibrationChargePix&)GetAverageHiGainArea(i); + // + // Ask for Hi-Gain saturation + // + if (hipix.GetSaturated() > fNumHiGainSaturationLimit*hipix.GetHGausHist()->GetEntries() && IsLoGain()) + { + MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainArea(i); + DrawDataCheckPixel(lopix,i ? gkLoGainOuterRefLines : gkLoGainInnerRefLines); + } + else + DrawDataCheckPixel(hipix,i ? gkHiGainOuterRefLines : gkHiGainInnerRefLines); + } +} + + +// -------------------------------------------------------------------------- +// +// Our own clone function is necessary since root 3.01/06 or Mars 0.4 +// I don't know the reason. +// +// Creates new MHCalibrationCam +// +TObject *MHCalibrationChargeCam::Clone(const char *name) const +{ + + const Int_t navhi = fAverageHiGainAreas->GetEntries(); + const Int_t navlo = fAverageLoGainAreas->GetEntries(); + const Int_t nsehi = fAverageHiGainSectors->GetEntries(); + const Int_t nselo = fAverageLoGainSectors->GetEntries(); + + // + // FIXME, this might be done faster and more elegant, by direct copy. + // + MHCalibrationChargeCam *cam = new MHCalibrationChargeCam(); + + cam->fAverageHiGainAreas->Expand(navhi); + cam->fAverageLoGainAreas->Expand(navlo); + cam->fAverageHiGainSectors->Expand(nsehi); + cam->fAverageLoGainSectors->Expand(nselo); + + cam->fAverageHiGainAreas->Expand(navhi); + cam->fAverageLoGainAreas->Expand(navlo); + cam->fAverageHiGainSectors->Expand(nsehi); + cam->fAverageLoGainSectors->Expand(nselo); + + for (int i=0; ifAverageHiGainAreas) [i] = (*fAverageHiGainAreas) [i]->Clone(); + for (int i=0; ifAverageLoGainAreas) [i] = (*fAverageLoGainAreas) [i]->Clone(); + for (int i=0; ifAverageHiGainSectors)[i] = (*fAverageHiGainSectors)[i]->Clone(); + for (int i=0; ifAverageLoGainSectors)[i] = (*fAverageLoGainSectors)[i]->Clone(); + + cam->fAverageAreaNum = fAverageAreaNum; + cam->fAverageAreaSat = fAverageAreaSat; + cam->fAverageAreaSigma = fAverageAreaSigma; + cam->fAverageAreaSigmaVar = fAverageAreaSigmaVar; + cam->fAverageAreaRelSigma = fAverageAreaRelSigma; + cam->fAverageAreaRelSigmaVar = fAverageAreaRelSigmaVar; + cam->fAverageSectorNum = fAverageSectorNum; + cam->fRunNumbers = fRunNumbers; + + cam->fColor = fColor; + cam->fPulserFrequency = fPulserFrequency; + + return cam; + +} + +void MHCalibrationChargeCam::DrawDataCheckPixel(MHCalibrationChargePix &pix, const Float_t refline[]) +{ + + TVirtualPad *newpad = gPad; + newpad->Divide(1,2); + newpad->cd(1); + + gPad->SetTicks(); + if (!pix.IsEmpty() && !pix.IsOnlyOverflow() && !pix.IsOnlyUnderflow()) + gPad->SetLogy(); + + TH1F *hist = pix.GetHGausHist(); + + TH2D *null = new TH2D("Null",hist->GetTitle(),100,pix.GetFirst() > -1. ? 0. : 100.,pix.GetLast()/2., + 100,0.,hist->GetEntries()/10.); + + null->SetDirectory(NULL); + null->SetBit(kCanDelete); + null->SetStats(kFALSE); + // + // set the labels bigger + // + TAxis *xaxe = null->GetXaxis(); + TAxis *yaxe = null->GetYaxis(); + xaxe->CenterTitle(); + yaxe->CenterTitle(); + xaxe->SetTitleSize(0.07); + yaxe->SetTitleSize(0.07); + xaxe->SetTitleOffset(0.7); + yaxe->SetTitleOffset(0.55); + xaxe->SetLabelSize(0.06); + yaxe->SetLabelSize(0.06); + + xaxe->SetTitle(hist->GetXaxis()->GetTitle()); + yaxe->SetTitle(hist->GetYaxis()->GetTitle()); + null->Draw(); + hist->Draw("same"); + + gStyle->SetOptFit(); + + if (pix.GetFGausFit()) + { + switch ( fColor ) + { + case MCalibrationCam::kGREEN: + pix.GetFGausFit()->SetLineColor(kGreen); + break; + case MCalibrationCam::kBLUE: + pix.GetFGausFit()->SetLineColor(kBlue); + break; + case MCalibrationCam::kUV: + pix.GetFGausFit()->SetLineColor(106); + break; + case MCalibrationCam::kCT1: + pix.GetFGausFit()->SetLineColor(006); + break; + default: + pix.GetFGausFit()->SetLineColor(kRed); + } + pix.GetFGausFit()->Draw("same"); + } + + DisplayRefLines(null,refline); + + newpad->cd(2); + gPad->SetTicks(); + + pix.DrawEvents(); + return; + +} + + +void MHCalibrationChargeCam::DisplayRefLines(const TH2D *hist, const Float_t refline[]) const +{ + + TLine *green1 = new TLine(refline[0],0.,refline[0],hist->GetYaxis()->GetXmax()); + green1->SetBit(kCanDelete); + green1->SetLineColor(kGreen); + green1->SetLineStyle(2); + green1->SetLineWidth(3); + green1->Draw(); + + TLine *green5 = new TLine(refline[6],0.,refline[6],hist->GetYaxis()->GetXmax()); + green5->SetBit(kCanDelete); + green5->SetLineColor(8); + green5->SetLineStyle(2); + green5->SetLineWidth(3); + green5->Draw(); + + TLine *blue1 = new TLine(refline[1],0.,refline[1],hist->GetYaxis()->GetXmax()); + blue1->SetBit(kCanDelete); + blue1->SetLineColor(227); + blue1->SetLineStyle(2); + blue1->SetLineWidth(3); + blue1->Draw(); + + TLine *blue5 = new TLine(refline[2],0.,refline[2],hist->GetYaxis()->GetXmax()); + blue5->SetBit(kCanDelete); + blue5->SetLineColor(68); + blue5->SetLineStyle(2); + blue5->SetLineWidth(3); + blue5->Draw(); + + TLine *blue10 = new TLine(refline[3],0.,refline[3],hist->GetYaxis()->GetXmax()); + blue10->SetBit(kCanDelete); + blue10->SetLineColor(4); + blue10->SetLineStyle(2); + blue10->SetLineWidth(3); + blue10->Draw(); + + TLine *uv10 = new TLine(refline[4],0.,refline[4],hist->GetYaxis()->GetXmax()); + uv10->SetBit(kCanDelete); + uv10->SetLineColor(106); + uv10->SetLineStyle(2); + uv10->SetLineWidth(3); + uv10->Draw(); + + TLine *ct1 = new TLine(refline[5],0.,refline[5],hist->GetYaxis()->GetXmax()); + ct1->SetBit(kCanDelete); + ct1->SetLineColor(6); + ct1->SetLineStyle(2); + ct1->SetLineWidth(3); + ct1->Draw(); + + TLegend *leg = new TLegend(0.7,0.35,0.9,0.99); + leg->SetBit(kCanDelete); + leg->AddEntry(green1,"1 Led GREEN","l"); + leg->AddEntry(green5,"5 Leds GREEN","l"); + leg->AddEntry(blue1,"1 Led BLUE","l"); + leg->AddEntry(blue5,"5 Leds BLUE","l"); + leg->AddEntry(blue10,"10 Leds BLUE","l"); + leg->AddEntry(uv10,"10 Leds UV","l"); + leg->AddEntry(ct1,"CT1-Pulser","l"); + + leg->Draw(); +} Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeCam.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeCam.h (revision 4929) @@ -0,0 +1,104 @@ +#ifndef MARS_MHCalibrationChargeCam +#define MARS_MHCalibrationChargeCam + +#ifndef MARS_MHCalibrationCam +#include "MHCalibrationCam.h" +#endif + +#ifndef MARS_MArrayI +#include "MArrayI.h" +#endif + +#ifndef MARS_MArrayD +#include "MArrayD.h" +#endif + +class TH2D; +class MRawEvtData; +class MExtractedSignalCam; +class MCalibrationChargePix; +class MHCalibrationChargePix; + +class MHCalibrationChargeCam : public MHCalibrationCam +{ +private: + + static const Axis_t fgChargeHiGainFirst; //! First Bin of HiGain Histograms (nw set to: -100.5) + static const Axis_t fgChargeHiGainLast; //! Last Bin of HiGain Histograms (nw set to: 999.5) + static const Axis_t fgChargeLoGainFirst; //! First Bin of LoGain Histograms (nw set to: -150.5) + static const Axis_t fgChargeLoGainLast; //! Last Bin of LoGain Histograms (nw set to: 499.5) + + static const Float_t fgNumHiGainSaturationLimit; //! The default for fNumHiGainSaturationLimit (now at: 0.01) + static const Float_t fgNumLoGainSaturationLimit; //! The default for fNumLoGainSaturationLimit (now at: 0.005) + static const Float_t fgTimeLowerLimit; //! Default for fTimeLowerLimit (now set to: 1.) + static const Float_t fgTimeUpperLimit; //! Default for fTimeUpperLimit (now set to: 2.) + + static const Float_t gkHiGainInnerRefLines[7]; //! + static const Float_t gkHiGainOuterRefLines[7]; //! + static const Float_t gkLoGainInnerRefLines[7]; //! + static const Float_t gkLoGainOuterRefLines[7]; //! + + MArrayD fSumhiarea ; //! + MArrayD fSumloarea ; //! + MArrayD fTimehiarea ; //! + MArrayD fTimeloarea ; //! + MArrayD fSumhisector; //! + MArrayD fSumlosector; //! + MArrayD fTimehisector; //! + MArrayD fTimelosector; //! + + MArrayI fSathiarea ; //! + MArrayI fSatloarea ; //! + MArrayI fSathisector; //! + MArrayI fSatlosector; //! + + Float_t fTimeLowerLimit; // Limit dist. to first signal slice (in units of FADC slices) + Float_t fTimeUpperLimit; // Limit dist. to last signal slice (in units of FADC slices) + Byte_t fFirstHiGain; // First used slice High Gain + Byte_t fLastHiGain; // Last used slice High Gain + Byte_t fFirstLoGain; // First used slice Low Gain + Byte_t fLastLoGain; // Last used slice Low Gain + + MRawEvtData *fRawEvt; //! Raw event data + + Bool_t SetupHists(const MParList *pList); + Bool_t ReInitHists(MParList *pList); + Bool_t FillHists(const MParContainer *par, const Stat_t w=1); + + void FinalizeAbsTimes (MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad, + Byte_t first, Byte_t last); + Bool_t FinalizeHists(); + void FinalizeBadPixels(); + + void DrawDataCheckPixel(MHCalibrationChargePix &pix, const Float_t refline[]); + void DisplayRefLines(const TH2D *hist, const Float_t refline[]) const; + +public: + + MHCalibrationChargeCam(const char *name=NULL, const char *title=NULL); + ~MHCalibrationChargeCam() {} + + TObject *Clone(const char *name="") const; + + void SetTimeLowerLimit ( const Float_t f=fgTimeLowerLimit ) { fTimeLowerLimit = f; } + void SetTimeUpperLimit ( const Float_t f=fgTimeUpperLimit ) { fTimeUpperLimit = f; } + + Bool_t GetPixelContent ( Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; + void DrawPixelContent( Int_t num ) const; + + // Draw + void Draw(const Option_t *opt); + + ClassDef(MHCalibrationChargeCam, 1) // Histogram class for Charge Camera Calibration +}; + +#endif + + + + + + + + + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePINDiode.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePINDiode.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePINDiode.cc (revision 4929) @@ -0,0 +1,398 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationChargePINDiode +// +// Histogram class for the charge calibration of the PIN Diode. +// Stores and fits the charges, the RMS of the charges and stores the +// location of the maximum FADC slice. Charges are taken from MExtractedSignalPINDiode. +// +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationChargePINDiode.h" + +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MExtractedSignalPINDiode.h" +#include "MCalibrationChargePINDiode.h" + +ClassImp(MHCalibrationChargePINDiode); + +using namespace std; + +const Axis_t MHCalibrationChargePINDiode::fgAbsTimeFirst = -0.5; +const Axis_t MHCalibrationChargePINDiode::fgAbsTimeLast = 29.5; +const Int_t MHCalibrationChargePINDiode::fgAbsTimeNbins = 30; +const Axis_t MHCalibrationChargePINDiode::fgChargeFirst = -0.5; +const Axis_t MHCalibrationChargePINDiode::fgChargeLast = 1999.5; +const Int_t MHCalibrationChargePINDiode::fgChargeNbins = 2000; +const Int_t MHCalibrationChargePINDiode::fgRmsChargeNbins = 200; +const Axis_t MHCalibrationChargePINDiode::fgRmsChargeFirst = 0.; +const Axis_t MHCalibrationChargePINDiode::fgRmsChargeLast = 200.; +const Float_t MHCalibrationChargePINDiode::fgTimeLowerLimit = 3.; +const Float_t MHCalibrationChargePINDiode::fgTimeUpperLimit = 4.; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - the default number for fAbsTimeFirst (fgAbsTimeFirst) +// - the default number for fAbsTimeLast (fgAbsTimeLast) +// - the default number for fAbsTimeNbins (fgAbsTimeNbins) +// - the default number for MHGausEvents::fNbins (fgChargeNbins) +// - the default number for MHGausEvents::fFirst (fgChargeFirst) +// - the default number for MHGausEvents::fLast (fgChargeLast) +// - the default number for fRmsChargeNbins (fgRmsChargeNbins) +// - the default number for fRmsChargeFirst (fgRmsChargeFirst) +// - the default number for fRmsChargeLast (fgRmsChargeLast) +// - the default number for fTimeLowerLimit (fgTimeLowerLimit) +// - the default number for fTimeUpperLimit (fgTimeUpperLimit) +// +// - the default name of the fHGausHist ("HCalibrationChargePINDiode") +// - the default title of the fHGausHist ("Distribution of Summed FADC slices PIN Diode") +// - the default x-axis title for fHGausHist ("Sum FADC Slices") +// - the default y-axis title for fHGausHist ("Nr. of events") +// - the default name of the fHAbsTime ("HAbsTimePINDiode") +// - the default title of the fHAbsTime ("Distribution of Absolute Arrival Times PIN Diode") +// - the default x-axis title for fHAbsTime ("Absolute Arrival Time [FADC slice nr]") +// - the default y-axis title for fHAbsTime ("Nr. of events") +// - the default name of the fHRmsCharge ("HRmsChargePINDiode") +// - the default title of the fHRmsCharge ("Distribution of Variances of summed FADC slices PIN Diode") +// - the default x-axis title for fHRmsCharge ("RMS (sum) [FADC slices]") +// - the default y-axis title for fHRmsCharge ("Nr. of events") +// - the default directory of the fHRmsCharge (NULL) +// - the current style for fHRmsCharge (NULL) +// +// Initializes: +// - fHRmsCharge() +// - all pointers to NULL +// +// Calls: +// - Clear() +// +MHCalibrationChargePINDiode::MHCalibrationChargePINDiode(const char *name, const char *title) + : fPINDiode(NULL), fSigPIN(NULL), fHRmsCharge() +{ + + fName = name ? name : "MHCalibrationChargePINDiode"; + fTitle = title ? title : "Fill the FADC sums of the PINDiode events and perform the fits"; + + SetAbsTimeFirst(); + SetAbsTimeLast(); + SetAbsTimeNbins(); + + SetNbins( fgChargeNbins ); + SetFirst( fgChargeFirst ); + SetLast ( fgChargeLast ); + + SetRmsChargeNbins(); + SetRmsChargeFirst(); + SetRmsChargeLast(); + + SetTimeLowerLimit(); + SetTimeUpperLimit(); + + fHGausHist.SetName("HCalibrationChargePINDiode"); + fHGausHist.SetTitle("Distribution of Summed FADC slices PIN Diode"); + fHGausHist.SetXTitle("Sum FADC Slices"); + fHGausHist.SetYTitle("Nr. of events"); + + fHAbsTime.SetName("HAbsTimePINDiode"); + fHAbsTime.SetTitle("Distribution of Absolute Arrival Times PIN Diode"); + fHAbsTime.SetXTitle("Absolute Arrival Time [FADC slice nr]"); + fHAbsTime.SetYTitle("Nr. of events"); + + fHRmsCharge.SetName("HRmsChargePINDiode"); + fHRmsCharge.SetTitle("Distribution of Variances of summed FADC slices PIN Diode"); + fHRmsCharge.SetXTitle("RMS (sum) [FADC slices]"); + fHRmsCharge.SetYTitle("Nr. of events"); + fHRmsCharge.UseCurrentStyle(); + fHRmsCharge.SetDirectory(NULL); + + Clear(); +} + +// -------------------------------------------------------------------------- +// +// Initializes Binning of the following histograms: +// - fHGausHist.SetBins(fNbins,fFirst,fLast); +// - fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast); +// - fHRmsCharge.SetBins(fRmsChargeNbins,fRmsChargeFirst,fRmsChargeLast); +// +Bool_t MHCalibrationChargePINDiode::SetupFill(const MParList *pList) +{ + + MHCalibrationPix::InitBins(); + + fHAbsTime. SetBins(fAbsTimeNbins, fAbsTimeFirst, fAbsTimeLast); + fHRmsCharge.SetBins(fRmsChargeNbins,fRmsChargeFirst,fRmsChargeLast); + + return kTRUE; + +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// - MExtractedSignalPINDiode +// - MCalibrationChargePINDiode +// +Bool_t MHCalibrationChargePINDiode::ReInit(MParList *pList) +{ + + fSigPIN = (MExtractedSignalPINDiode*)pList->FindCreateObj("MExtractedSignalPINDiode"); + if (!fSigPIN) + { + *fLog << err << "MExtractedSignalPINDiode not found... aborting " << endl; + return kFALSE; + } + + fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode"); + if (!fPINDiode) + { + *fLog << err << "MCalibrationChargePINDiode not found... aborting " << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Retrieves from MExtractedSignalPINDiode: +// - Number of used FADC samples via MExtractedSignalPINDiode::GetNumFADCSamples() +// - Extracted signal via MExtractedSignalPINDiode::GetExtractedSignal() +// - Signal Rms MExtractedSignalPINDiode::GetExtractedRms() +// - Arrival Time MExtractedSignalPINDiode::GetExtractedTime() +// +// Fills the following histograms: +// - MHGausEvents::FillHistAndArray(signal) +// - MHCalibrationChargePix::FillAbsTime(time); +// - FillRmsCharge(rms); +// +Bool_t MHCalibrationChargePINDiode::Fill(const MParContainer *par, const Stat_t w) +{ + + MExtractedSignalPINDiode *extractor = (MExtractedSignalPINDiode*)par; + + if (!extractor) + { + *fLog << err << "No argument in MExtractedSignalPINDiode::Fill... abort." << endl; + return kFALSE; + } + + Float_t slices = (Float_t)extractor->GetNumFADCSamples(); + + if (slices == 0.) + { + *fLog << err << "Number of used signal slices in MExtractedSignalPINDiode is zero ... abort." + << endl; + return kFALSE; + } + + const Float_t signal = (float)extractor->GetExtractedSignal(); + const Float_t time = extractor->GetExtractedTime(); + const Float_t rms = extractor->GetExtractedRms(); + + FillHistAndArray(signal); + FillAbsTime(time); + FillRmsCharge(rms); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Returns kTRUE, if empty +// +// Performs the following fits: +// - MHGausEvents::FitGaus() +// - FitRmsCharge() +// +// Creates the fourier spectrum (MHGausEvents::CreateFourierSpectrum() +// and sets bit MCalibrationChargePINDiode::SetOscillating( MHGausEvents::IsFourierSpectrumOK() ) +// Retrieves the results of the following fits and stores them in MCalibrationChargePINDiode: +// - Mean Charge and Error +// - Sigma Charge and Error +// - Fit Probability +// - Abs Time Mean +// - Abs Time Rms +// - Rms Charge Mean and Error +// - Rms Charge Sigma and Error +// +// Performs one consistency check on the arrival time: +// The check returns kFALSE if: +// +// -The mean arrival time is in fTimeLowerLimit slices from the lower edge +// and fUpperLimit slices from the upper edge +// +Bool_t MHCalibrationChargePINDiode::Finalize() +{ + + if (IsGausFitOK() || IsEmpty()) + return kTRUE; + + FitGaus(); + FitRmsCharge(); + + CreateFourierSpectrum(); + fPINDiode->SetOscillating ( !IsFourierSpectrumOK() ); + + fPINDiode->SetMean ( fMean ); + fPINDiode->SetMeanVar ( fMeanErr * fMeanErr ); + fPINDiode->SetSigma ( fSigma ); + fPINDiode->SetSigmaVar ( fSigmaErr * fMeanErr ); + fPINDiode->SetProb ( fProb ); + + fPINDiode->SetAbsTimeMean( GetAbsTimeMean() ); + fPINDiode->SetAbsTimeRms( GetAbsTimeRms() ); + + fPINDiode->SetRmsChargeMean( GetRmsChargeMean() ); + fPINDiode->SetRmsChargeMeanErr( GetRmsChargeMeanErr() ); + fPINDiode->SetRmsChargeSigma( GetRmsChargeSigma() ); + fPINDiode->SetRmsChargeSigmaErr( GetRmsChargeSigmaErr() ); + + fPINDiode->SetValid(kTRUE); + + const Byte_t loweredge = fSigPIN->GetFirstUsedSlice(); + const Byte_t upperedge = fSigPIN->GetLastUsedSlice(); + const Float_t lowerlimit = (Float_t)loweredge + fTimeLowerLimit; + const Float_t upperlimit = (Float_t)upperedge + fTimeUpperLimit; + + if (GetAbsTimeMean() < lowerlimit) + { + *fLog << warn << GetDescriptor() + << Form("%s%3.1f%s%2.1f%s",": Mean ArrivalTime: ",GetAbsTimeMean()," smaller than ", + lowerlimit," FADC slices from lower edge in PIN Diode") << endl; + *fLog << warn << GetDescriptor() << ": No PIN Diode calibration!! " << endl; + fPINDiode->SetValid(kFALSE); + } + + if ( GetAbsTimeMean() > upperlimit ) + { + *fLog << warn << GetDescriptor() + << Form("%s%3.1f%s%2.1f%s",": Mean ArrivalTime: ",GetAbsTimeMean()," bigger than ", + upperlimit," FADC slices from upper edge in PIN Diode") << endl; + *fLog << warn << GetDescriptor() << ": No PIN Diode calibration!! " << endl; + fPINDiode->SetValid(kFALSE); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fills fHRmsCharge with q +// Returns kFALSE, if overflow or underflow occurred, else kTRUE +// +Bool_t MHCalibrationChargePINDiode::FillRmsCharge(const Float_t q) +{ + return fHRmsCharge.Fill(q) > -1; +} + +// ----------------------------------------------------------- +// +// Fits -- not yet implemented +// +Bool_t MHCalibrationChargePINDiode::FitRmsCharge(Option_t *option) +{ + return 1; +} + + +// ------------------------------------------------------------------------- +// +// Draw the histogram +// +// The following options can be chosen: +// +// "": displays the fHGausHist with fits and fHRmsCharge +// "all": executes additionally MHCalibrationPix::Draw(), with option "fourierevents" +// +void MHCalibrationChargePINDiode::Draw(const Option_t *opt) +{ + + TString option(opt); + option.ToLower(); + + Int_t win = 1; + + TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,900, 600); + TVirtualPad *pad = NULL; + + oldpad->SetBorderMode(0); + + if (option.Contains("all")) + { + option.ReplaceAll("all",""); + oldpad->Divide(2,1); + win = 2; + oldpad->cd(1); + TVirtualPad *newpad = gPad; + pad = newpad; + pad->Divide(1,2); + pad->cd(1); + } + else + { + pad = oldpad; + pad->Divide(1,2); + pad->cd(1); + } + + if (IsEmpty()) + return; + + if (!IsOnlyOverflow() && !IsOnlyUnderflow()) + gPad->SetLogy(); + + gPad->SetTicks(); + + fHGausHist.Draw(opt); + if (fFGausFit) + { + fFGausFit->SetLineColor(IsGausFitOK() ? kGreen : kRed); + fFGausFit->Draw("same"); + } + + pad->cd(2); + fHRmsCharge.Draw(opt); + + oldpad->cd(2); + MHCalibrationPix::Draw("fourierevents"); +} + + + + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePINDiode.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePINDiode.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePINDiode.h (revision 4929) @@ -0,0 +1,84 @@ +#ifndef MARS_MHCalibrationChargePINDiode +#define MARS_MHCalibrationChargePINDiode + + +#ifndef MARS_MHCalibrationChargePix +#include "MHCalibrationChargePix.h" +#endif + +class TH1F; +class MExtractedSignalPINDiode; +class MCalibrationChargePINDiode; +class MHCalibrationChargePINDiode : public MHCalibrationChargePix +{ +private: + + static const Axis_t fgAbsTimeFirst; //! Default for fAbsTimeFirst (now set to: -0.5 ) + static const Axis_t fgAbsTimeLast; //! Default for fAbsTimeLast (now set to: 29.5 ) + static const Int_t fgAbsTimeNbins; //! Default for fAbsTimeNBins (now set to: 30 ) + static const Int_t fgChargeNbins; //! Default for fNBins (now set to: 200 ) + static const Axis_t fgChargeFirst; //! Default for fFirst (now set to: -0.5 ) + static const Axis_t fgChargeLast; //! Default for fLast (now set to: 199.5 ) + static const Int_t fgRmsChargeNbins; //! Default for fRmsChargeNBins (now set to: 100 ) + static const Axis_t fgRmsChargeFirst; //! Default for fRmsChargeFirst (now set to: 0. ) + static const Axis_t fgRmsChargeLast; //! Default for fRmsChargeLast (now set to: 100. ) + static const Float_t fgTimeLowerLimit; //! Default for fTimeLowerLimit (now set to: 1.) + static const Float_t fgTimeUpperLimit; //! Default for fTimeUpperLimit (now set to: 2.) + + MCalibrationChargePINDiode *fPINDiode; //! Storage container of the results + MExtractedSignalPINDiode *fSigPIN; //! Storage container of extracted signal + + TH1F fHRmsCharge; // Histogram containing Variance of summed FADC slices + + Axis_t fRmsChargeFirst; // Lower bound bin used for the fHRmsCharge + Axis_t fRmsChargeLast; // Upper bound bin used for the fHRmsCharge + Int_t fRmsChargeNbins; // Number of bins used for the fHRmsCharge + Float_t fRmsChargeMean; // Mean of the Gauss fit + Float_t fRmsChargeMeanErr; // Error of the mean of the Gauss fit + Float_t fRmsChargeSigma; // Sigma of the Gauss fit + Float_t fRmsChargeSigmaErr; // Error of the sigma of the Gauss fit + Float_t fTimeLowerLimit; // Limit dist. to first signal slice (units: FADC slices) + Float_t fTimeUpperLimit; // Limit dist. to last signal slice (units: FADC slices) + +public: + + MHCalibrationChargePINDiode(const char *name=NULL, const char *title=NULL); + ~MHCalibrationChargePINDiode(){} + + Bool_t SetupFill(const MParList *pList); + Bool_t ReInit ( MParList *pList); + Bool_t Fill (const MParContainer *par, const Stat_t w=1); + Bool_t Finalize (); + + // Draw + void Draw(Option_t *opt=""); + + // Getters + TH1F *GetHRmsCharge() { return &fHRmsCharge; } + const TH1F *GetHRmsCharge() const { return &fHRmsCharge; } + Float_t GetRmsChargeMean() const { return fRmsChargeMean; } + Float_t GetRmsChargeMeanErr() const { return fRmsChargeMeanErr; } + Float_t GetRmsChargeSigma() const { return fRmsChargeSigma; } + Float_t GetRmsChargeSigmaErr() const { return fRmsChargeSigmaErr; } + + // Fill histos + Bool_t FillRmsCharge(const Float_t q); + + // Fits + Bool_t FitRmsCharge(Option_t *option="RQ0"); + + // Setters + void SetAbsTimeNbins ( const Int_t bins =fgAbsTimeNbins ) { fAbsTimeNbins = bins; } + void SetAbsTimeFirst ( const Axis_t first=fgAbsTimeFirst ) { fAbsTimeFirst = first; } + void SetAbsTimeLast ( const Axis_t last =fgAbsTimeLast ) { fAbsTimeLast = last; } + void SetRmsChargeNbins ( const Int_t bins =fgRmsChargeNbins ) { fRmsChargeNbins = bins; } + void SetRmsChargeFirst ( const Axis_t first=fgRmsChargeFirst ) { fRmsChargeFirst = first; } + void SetRmsChargeLast ( const Axis_t last =fgRmsChargeLast ) { fRmsChargeLast = last; } + void SetTimeLowerLimit ( const Float_t f=fgTimeLowerLimit ) { fTimeLowerLimit = f; } + void SetTimeUpperLimit ( const Float_t f=fgTimeUpperLimit ) { fTimeUpperLimit = f; } + + ClassDef(MHCalibrationChargePINDiode, 1) // Histogram class for Charge PIN Diode Calibration +}; + +#endif + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePix.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePix.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePix.cc (revision 4929) @@ -0,0 +1,309 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +////////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationChargePix +// +// Histogram class for the charge calibration. +// Stores and fits the charges and stores the location of the maximum FADC +// slice. Charges are taken from MExtractedSignalPix. +// +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationChargePix.h" + +#include +#include + +#include +#include +#include +#include + +#include "MH.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHCalibrationChargePix); + +using namespace std; + +const Int_t MHCalibrationChargePix::fgChargeNbins = 2000; +const Axis_t MHCalibrationChargePix::fgChargeFirst = -0.5; +const Axis_t MHCalibrationChargePix::fgChargeLast = 1999.5; +const Int_t MHCalibrationChargePix::fgAbsTimeNbins = 15; +const Axis_t MHCalibrationChargePix::fgAbsTimeFirst = -0.5; +const Axis_t MHCalibrationChargePix::fgAbsTimeLast = 14.5; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - the default number for fNbins (fgChargeNbins) +// - the default number for fFirst (fgChargeFirst) +// - the default number for fLast (fgChargeLast) +// - the default number for fAbsTimeNbins (fgAbsTimeNbins) +// - the default number for fAbsTimeFirst (fgAbsTimeFirst) +// - the default number for fAbsTimeLast (fgAbsTimeLast) +// +// - the default name of the fHGausHist ("HCalibrationCharge") +// - the default title of the fHGausHist ("Distribution of Summed FADC slices Pixel ") +// - the default x-axis title for fHGausHist ("Sum FADC Slices") +// - the default y-axis title for fHGausHist ("Nr. of events") +// +// - the default name of the fHAbsTime ("HAbsTimePixel") +// - the default title of the fHAbsTime ("Distribution of Absolute Arrival Times Pixel ") +// - the default x-axis title for fHAbsTime ("Absolute Arrival Time [FADC slice nr]") +// - the default y-axis title for fHAbsTime ("Nr. of events"); +// - the default directory of the fHAbsTime (NULL) +// - the current style for fHAbsTime +// +// Initializes: +// - fHAbsTime() +// +// Calls: +// - Clear(); +// +MHCalibrationChargePix::MHCalibrationChargePix(const char *name, const char *title) + : fHAbsTime() +{ + + fName = name ? name : "MHCalibrationChargePix"; + fTitle = title ? title : "Statistics of the FADC sums of calibration events"; + + SetNbins ( fgChargeNbins ); + SetFirst ( fgChargeFirst ); + SetLast ( fgChargeLast ); + + SetAbsTimeNbins(); + SetAbsTimeFirst(); + SetAbsTimeLast(); + + fHGausHist.SetName("HCalibrationCharge"); + fHGausHist.SetTitle("Distribution of Summed FADC slices Pixel"); + fHGausHist.SetXTitle("Sum FADC Slices"); + fHGausHist.SetYTitle("Nr. of events"); + + fHAbsTime.SetName("HAbsTimePixel"); + fHAbsTime.SetTitle("Distribution of Absolute Arrival Times Pixel "); + fHAbsTime.SetXTitle("Absolute Arrival Time [FADC slice nr]"); + fHAbsTime.SetYTitle("Nr. of events"); + + fHAbsTime.UseCurrentStyle(); + fHAbsTime.SetDirectory(NULL); + + Clear(); + +} + +// -------------------------------------------------------------------------- +// +// Sets: +// - fHGausHist.SetBins(fNbins,fFirst,fLast); +// - fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast); +// +void MHCalibrationChargePix::InitBins() +{ + MHGausEvents::InitBins(); + fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast); +} + +// -------------------------------------------------------------------------- +// +// Empty function to overload MHGausEvents::Reset() +// +void MHCalibrationChargePix::Reset() +{ + MHGausEvents::Reset(); + fHAbsTime.Reset(); +} + +// -------------------------------------------------------------------------- +// +// Calls MHCalibrationPix::ChangeHistId() +// +// Add id to names and titles of: +// - fHAbsTime +// +void MHCalibrationChargePix::ChangeHistId(Int_t id) +{ + + MHCalibrationPix::ChangeHistId(id); + + fHAbsTime.SetName (Form("%s%d", fHAbsTime.GetName(), id)); + fHAbsTime.SetTitle(Form("%s%d", fHAbsTime.GetTitle(), id)); + +} + +// -------------------------------------------------------------------------- +// +// returns fHGausHist.Integral("width") +// +const Float_t MHCalibrationChargePix::GetIntegral() const +{ + return fHGausHist.Integral("width"); +} + +// -------------------------------------------------------------------------- +// +// returns fHAbsTime.GetMean() +// +const Float_t MHCalibrationChargePix::GetAbsTimeMean() const +{ + return fHAbsTime.GetMean(); +} + +// -------------------------------------------------------------------------- +// +// returns fHAbsTime.GetRMS() +// +const Float_t MHCalibrationChargePix::GetAbsTimeRms() const +{ + return fHAbsTime.GetRMS(); +} + +// -------------------------------------------------------------------------- +// +// Fills fHAbsTime with t +// Returns kFALSE, if overflow or underflow occurred, else kTRUE +// +Bool_t MHCalibrationChargePix::FillAbsTime(Float_t t) +{ + return fHAbsTime.Fill(t) > -1; +} + +// ----------------------------------------------------------------------------- +// +// Default draw: +// +// The following options can be chosen: +// +// "": displays the fHGausHist and the fHAbsTime +// "all": executes additionally MHCalibrationPix::Draw(), with options +// +// The following picture shows a typical outcome of call to Draw("all"): +// One the left side: +// - The distribution of the values with the Gauss fit are shown (points connected +// with each other). The line is green, thus the Gauss fit has a probability higher +// than 0.5%. +// - The distribution of the positions of the maximum (abs. arrival times) +// is displayed. Most of the events have an arrival time of slice 7 (==hardware:8) +// +// On the right side: +// - The first plot shows the distribution of the values with the Gauss fit +// with error bars +// - The second plot shows the TGraph with the events vs. time +// - The third plot shows the fourier transform and a peak at 100 Hz. +// - The fourth plot shows the projection of the fourier components and an exponential +// fit, with the result that the observed deviation is not statistical, but signficant with a +// probability smaller than 0.5%. +// +//Begin_Html +/* + +*/ +//End_Html +// +void MHCalibrationChargePix::Draw(const Option_t *opt) +{ + + TString option(opt); + option.ToLower(); + + Int_t win = 1; + + TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,600, 600); + TVirtualPad *pad = NULL; + + if (option.Contains("all")) + { + option.ReplaceAll("all",""); + oldpad->Divide(2,1); + win = 2; + oldpad->cd(1); + TVirtualPad *newpad = gPad; + pad = newpad; + pad->Divide(1,2); + pad->cd(1); + } + else if (option.Contains("datacheck")) + { + MHCalibrationPix::Draw("events"); + return; + } + else + { + pad = oldpad; + pad->Divide(1,2); + pad->cd(1); + } + /* + else + { + option.ReplaceAll("time",""); + pad = oldpad; + pad->Divide(1,2); + pad->cd(1); + } + */ + if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow()) + gPad->SetLogy(); + + gPad->SetTicks(); + + fHGausHist.GetXaxis()->SetLabelSize(0.06); + fHGausHist.GetYaxis()->SetLabelSize(0.07); + fHGausHist.GetXaxis()->SetLabelOffset(0.01); + fHGausHist.GetYaxis()->SetLabelOffset(0.01); + fHGausHist.GetXaxis()->SetTitleSize(0.065); + fHGausHist.GetYaxis()->SetTitleSize(0.07); + fHGausHist.GetXaxis()->SetTitleOffset(0.6); + fHGausHist.GetYaxis()->SetTitleOffset(0.6); + fHGausHist.Draw(); + if (fFGausFit) + { + fFGausFit->SetLineColor(IsGausFitOK() ? kGreen : kRed); + fFGausFit->Draw("same"); + } + + pad->cd(2); + gPad->SetTicks(); + + fHAbsTime.GetXaxis()->SetLabelSize(0.06); + fHAbsTime.GetYaxis()->SetLabelSize(0.07); + fHAbsTime.GetXaxis()->SetLabelOffset(0.01); + fHAbsTime.GetYaxis()->SetLabelOffset(0.01); + fHAbsTime.GetXaxis()->SetTitleSize(0.065); + fHAbsTime.GetYaxis()->SetTitleSize(0.07); + fHAbsTime.GetXaxis()->SetTitleOffset(0.6); + fHAbsTime.GetYaxis()->SetTitleOffset(0.6); + fHAbsTime.Draw(); + + if (win < 2) + return; + + oldpad->cd(2); + MHCalibrationPix::Draw("fourierevents"); + +} Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePix.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePix.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePix.h (revision 4929) @@ -0,0 +1,61 @@ +#ifndef MARS_MHCalibrationChargePix +#define MARS_MHCalibrationChargePix + +#ifndef MARS_MHCalibrationPix +#include "MHCalibrationPix.h" +#endif + +class MHCalibrationChargePix : public MHCalibrationPix +{ + +private: + + static const Int_t fgChargeNbins; // Default for fNBins (now set to: 2000 ) + static const Axis_t fgChargeFirst; // Default for fFirst (now set to: -0.5 ) + static const Axis_t fgChargeLast; // Default for fLast (now set to: 1999.5) + static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 ) + static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) + static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 ) + +protected: + + TH1F fHAbsTime; // Histogram containing the absolute arrival times + + Int_t fAbsTimeNbins; // Number of bins used for the fHAbsTime + Axis_t fAbsTimeFirst; // Lower bound bin used for the fHAbsTime + Axis_t fAbsTimeLast; // Upper bound bin used for the fHAbsTime + +public: + + MHCalibrationChargePix(const char *name=NULL, const char *title=NULL); + ~MHCalibrationChargePix() {} + + void Reset(); + void InitBins(); + + // Setters + virtual void SetAbsTimeNbins(const Int_t bins =fgAbsTimeNbins) { fAbsTimeNbins = bins; } + virtual void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst) { fAbsTimeFirst = first; } + virtual void SetAbsTimeLast( const Axis_t last =fgAbsTimeLast) { fAbsTimeLast = last; } + + // Getters + TH1F *GetHAbsTime() { return &fHAbsTime; } + const TH1F *GetHAbsTime() const { return &fHAbsTime; } + + const Float_t GetAbsTimeMean( ) const; + const Float_t GetAbsTimeRms() const; + const Float_t GetIntegral() const; + + // Fill histos + Bool_t FillAbsTime(const Float_t t); + + // Draws + virtual void Draw(Option_t *opt=""); + + // Miscelleaneous + void ChangeHistId(Int_t id); + + ClassDef(MHCalibrationChargePix, 1) // Base Histogram class for Charge Pixel Calibration +}; + +#endif /* MARS_MHCalibrationChargePix */ Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationPix.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationPix.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationPix.cc (revision 4929) @@ -0,0 +1,327 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationPix +// +// A base class for events which are believed to follow a Gaussian distribution +// with time, e.g. calibration events, observables containing white noise, ... +// +// MHCalibrationPix derives from MHGausEvents, thus all features of +// MHGausEvents can be used by a class deriving from MHCalibrationPix +// +// See also: MHGausEvents +// +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationPix.h" + +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHCalibrationPix); + +using namespace std; + +const Float_t MHCalibrationPix::fgBlackoutLimit = 5.; +const Float_t MHCalibrationPix::fgPickupLimit = 5.; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// Sets: +// - the default number for fPickupLimit (fgPickupLimit) +// - the default number for fBlackoutLimit (fgBlackoutLimit) +// +// Initializes: +// - all variables to 0. +// +MHCalibrationPix::MHCalibrationPix(const char *name, const char *title) + : fEventFrequency(0), fPixId(-1) +{ + + fName = name ? name : "MHCalibrationPix"; + fTitle = title ? title : "Calibration histogram events"; + + Clear(); + + SetBlackoutLimit(); + SetPickupLimit(); +} + + + +// -------------------------------------------------------------------------- +// +// Default Clear(), can be overloaded. +// +// Sets: +// - all other pointers to NULL +// - all variables to 0., except fPixId to -1 and keep fEventFrequency +// - all flags to kFALSE +// +// Deletes (if not NULL): +// - all pointers +// +void MHCalibrationPix::Clear(Option_t *o) +{ + + MHGausEvents::Clear(); + fSaturated = 0; +} + + +// ----------------------------------------------------------------------------- +// +// Bypasses the Gauss fit by taking mean and RMS from the histogram +// +// Errors are determined in the following way: +// MeanErr = RMS / Sqrt(entries) +// SigmaErr = RMS / (2.*Sqrt(entries) ) +// +void MHCalibrationPix::BypassFit() +{ + + const Stat_t entries = fHGausHist.GetEntries(); + + if (entries <= 0.) + { + *fLog << warn << GetDescriptor() + << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl; + return; + } + + fMean = fHGausHist.GetMean(); + fMeanErr = fHGausHist.GetRMS() / TMath::Sqrt(entries); + fSigma = fHGausHist.GetRMS() ; + fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.; +} + +// -------------------------------------------------------------------------- +// +// - Set fPixId to id +// +// Add id to names and titles of: +// - fHGausHist +// +void MHCalibrationPix::ChangeHistId(const Int_t id) +{ + + fPixId = id; + + fHGausHist.SetName( Form("%s%d", fHGausHist.GetName(), id)); + fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id)); + + fName = Form("%s%d", fName.Data(), id); + fTitle = Form("%s%d", fTitle.Data(), id); + +} + + +// ----------------------------------------------------------------------------- +// +// Create the x-axis for the event graph +// +Float_t *MHCalibrationPix::CreateEventXaxis(Int_t n) +{ + + Float_t *xaxis = new Float_t[n]; + + if (fEventFrequency) + for (Int_t i=0;iGetXaxis()->SetTitle((fEventFrequency) ? "Time [s]" : "Event Nr."); +} + +// ---------------------------------------------------------------------------------- +// +// Create a graph to display the array fPowerSpectrum +// If the variable fEventFrequency is set, the x-axis is transformed into real frequency. +// +void MHCalibrationPix::CreateGraphPowerSpectrum() +{ + + MHGausEvents::CreateGraphPowerSpectrum(); + + fGraphPowerSpectrum->GetXaxis()->SetTitle((fEventFrequency) ? "Frequency [Hz]" : "Frequency"); +} + +// ------------------------------------------------------------------------------- +// +// Return the number of "blackout" events, which are events with values higher +// than fBlackoutLimit sigmas from the mean +// +// +const Double_t MHCalibrationPix::GetBlackout() const +{ + + if ((fMean == 0.) && (fSigma == 0.)) + return -1.; + + const Int_t first = fHGausHist.GetXaxis()->GetFirst(); + const Int_t last = fHGausHist.GetXaxis()->FindBin(fMean-fBlackoutLimit*fSigma); + + if (first >= last) + return 0.; + + return fHGausHist.Integral(first, last, "width"); +} + + +// ------------------------------------------------------------------------------- +// +// Return the number of "pickup" events, which are events with values higher +// than fPickupLimit sigmas from the mean +// +// +const Double_t MHCalibrationPix::GetPickup() const +{ + + if ((fMean == 0.) && (fSigma == 0.)) + return -1.; + + const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma); + const Int_t last = fHGausHist.GetXaxis()->GetLast(); + + if (first >= last) + return 0.; + + return fHGausHist.Integral(first, last, "width"); +} + + +// -------------------------------------------------------------------------- +// +// Re-normalize the results, has to be overloaded +// +void MHCalibrationPix::Renorm() +{ +} + +// ----------------------------------------------------------------------------- +// +// If flag IsGausFitOK() is set (histogram already successfully fitted), +// returns kTRUE +// +// If both fMean and fSigma are still zero, call FitGaus() +// +// Repeats the Gauss fit in a smaller range, defined by: +// +// min = GetMean() - fBlackoutLimit * GetSigma(); +// max = GetMean() + fPickupLimit * GetSigma(); +// +// The fit results are retrieved and stored in class-own variables. +// +// A flag IsGausFitOK() is set according to whether the fit probability +// is smaller or bigger than fProbLimit, whether the NDF is bigger than +// fNDFLimit and whether results are NaNs. +// +Bool_t MHCalibrationPix::RepeatFit(const Option_t *option) +{ + + if (IsGausFitOK()) + return kTRUE; + + if ((fMean == 0.) && (fSigma == 0.)) + return FitGaus(); + + // + // Get new fitting ranges + // + Axis_t rmin = fMean - fBlackoutLimit * fSigma; + Axis_t rmax = fMean + fPickupLimit * fSigma; + + Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()); + Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ; + + fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax); + + fHGausHist.Fit(fFGausFit,option); + + fMean = fFGausFit->GetParameter(1); + fSigma = fFGausFit->GetParameter(2); + fMeanErr = fFGausFit->GetParError(1) ; + fSigmaErr = fFGausFit->GetParError(2) ; + fProb = fFGausFit->GetProb() ; + + // + // The fit result is accepted under condition: + // 1) The results are not nan's + // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit) + // 3) The Probability is greater than fProbLimit (default: fgProbLimit) + // + if ( TMath::IsNaN ( fMean ) + || TMath::IsNaN ( fMeanErr ) + || TMath::IsNaN ( fProb ) + || TMath::IsNaN ( fSigma ) + || TMath::IsNaN ( fSigmaErr ) + || fFGausFit->GetNDF() < fNDFLimit + || fProb < fProbLimit ) + return kFALSE; + + SetGausFitOK(kTRUE); + return kTRUE; + +} + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationPix.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationPix.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationPix.h (revision 4929) @@ -0,0 +1,60 @@ +#ifndef MARS_MHCalibrationPix +#define MARS_MHCalibrationPix + +#ifndef MARS_MHGausEvents +#include "MHGausEvents.h" +#endif + +class MHCalibrationPix : public MHGausEvents +{ +private: + + const static Float_t fgBlackoutLimit; //! Default for fBlackoutLimit (now set to: 5. ) + const static Float_t fgPickupLimit; //! Default for fPickupLimit (now set to: 5. ) + + Float_t fEventFrequency; // Event frequency in Hertz (to be set) + + Float_t *CreateEventXaxis(Int_t n); // Create an x-axis for the Event TGraphs + Float_t *CreatePSDXaxis(Int_t n); // Create an x-axis for the PSD TGraphs + void CreateGraphEvents(); // Create the TGraph fGraphEvents of fEvents + void CreateGraphPowerSpectrum(); // Create the TGraph fGraphPowerSpectrum out of fPowerSpectrum + +protected: + + Float_t fBlackoutLimit; // Lower nr sigmas from mean until event is considered blackout + Int_t fSaturated; // Number of events classified as saturated + Float_t fPickupLimit; // Upper nr sigmas from mean until event is considered pickup + Int_t fPixId; // Pixel ID + +public: + + MHCalibrationPix(const char* name=NULL, const char* title=NULL); + ~MHCalibrationPix() {} + + void Clear(Option_t *o=""); + + // Getters + const Double_t GetBlackout() const; + const Double_t GetPickup() const; + const Int_t GetPixId() const { return fPixId; } + const Float_t GetSaturated() const { return fSaturated; } + + // Fits + Bool_t RepeatFit(const Option_t *option="RQ0"); // Repeat fit within limits defined by fPickupLimit + void BypassFit(); // Take mean and RMS from the histogram + + // Setters + void SetBlackoutLimit ( const Float_t lim=fgBlackoutLimit ) { fBlackoutLimit = lim; } + void SetEventFrequency ( const Float_t f ) { fEventFrequency = f; } + void SetPickupLimit ( const Float_t lim=fgPickupLimit ) { fPickupLimit = lim; } + void SetPixId ( const Int_t i ) { fPixId = i; } + void SetSaturated ( const Int_t i ) { fSaturated += i; } + + // Miscelleaneous + virtual void ChangeHistId(const Int_t id); // Changes names and titles of the histogram + virtual void Renorm(); // Re-normalize the results + + ClassDef(MHCalibrationPix, 1) // Base class for calibration events +}; + +#endif Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimeCam.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimeCam.cc (revision 4929) @@ -0,0 +1,625 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationRelTimeCam +// +// Fills the extracted relative arrival times of MArrivalTimeCam into +// the MHCalibrationPix-classes MHCalibrationRelTimePix for every: +// +// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray +// or MHCalibrationCam::fHiGainArray, respectively, depending if +// MArrivalTimePix::IsLoGainUsed() is set. +// +// - 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 +// +// Every relative time is calculated as the difference between the individual +// pixel arrival time and the one of pixel 1 (hardware number: 2). +// The relative times 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 +// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit()) +// In case this does not make the fit valid, the histogram means and RMS's are +// taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kRelTimeNotFitted ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas +// from the mean are counted as Pickup events (stored in MHCalibrationPix::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::kRelTimeOscillating ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// This same procedure is performed for the average pixels. +// +// The following results are written into MCalibrationRelTimeCam: +// +// - 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 "MHCalibrationRelTimeCam.h" +#include "MHCalibrationRelTimePix.h" + +#include "MHCalibrationPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MCalibrationIntensityRelTimeCam.h" + +#include "MCalibrationRelTimeCam.h" +#include "MCalibrationRelTimePix.h" +#include "MCalibrationPix.h" + +#include "MArrivalTimeCam.h" +#include "MArrivalTimePix.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" + +ClassImp(MHCalibrationRelTimeCam); + +using namespace std; + +const Float_t MHCalibrationRelTimeCam::fgNumHiGainSaturationLimit = 0.25; +const UInt_t MHCalibrationRelTimeCam::fgReferencePixel = 1; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - fReferencePixel to fgReferencePixel +// +MHCalibrationRelTimeCam::MHCalibrationRelTimeCam(const char *name, const char *title) +{ + + fName = name ? name : "MHCalibrationRelTimeCam"; + fTitle = title ? title : "Histogram class for the relative time calibration of the camera"; + + SetNumHiGainSaturationLimit(fgNumHiGainSaturationLimit); + SetReferencePixel(); +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// - MCalibrationRelTimeCam +// +// Searches pointer to: +// - MArrivalTimeCam +// +// Initializes, if empty to MGeomCam::GetNumPixels(): +// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray +// +// 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::fLoGainArray +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +// Sets Titles and Names for the Histograms +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +// Sets number of bins to MHCalibrationCam::fAverageNbins for: +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +Bool_t MHCalibrationRelTimeCam::ReInitHists(MParList *pList) +{ + + fIntensCam = (MCalibrationIntensityCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityRelTimeCam")); + if (fIntensCam) + *fLog << inf << "Found MCalibrationIntensityRelTimeCam ... " << endl; + else + { + fCam = (MCalibrationCam*)pList->FindObject(AddSerialNumber("MCalibrationRelTimeCam")); + if (!fCam) + { + fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationRelTimeCam")); + if (!fCam) + { + *fLog << err << "Cannot find nor create MCalibrationRelTimeCam ... abort." << endl; + return kFALSE; + } + fCam->Init(*fGeom); + } + } + + MArrivalTimeCam *signal = (MArrivalTimeCam*)pList->FindObject("MArrivalTimeCam"); + if (!signal) + { + *fLog << err << "MArrivalTimeCam 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 && IsLoGain() ) + { + fLoGainArray->Expand(npixels); + for (Int_t i=0; iGetEntries()==0) + { + fAverageHiGainAreas->Expand(nareas); + + for (Int_t j=0; jSetTitle("Rel. Arr. Times average HiGain Area Idx "); + + if (fGeom->InheritsFrom("MGeomCamMagic")) + { + hist.GetHGausHist()->SetTitle(Form("%s%s%s","Signal averaged on event-by-event basis ", + j==0 ? "Inner Pixels " : "Outer Pixels ","High Gain Runs: ")); + hist.InitBins(); + hist.SetEventFrequency(fPulserFrequency); + } + else + { + hist.GetHGausHist()->SetTitle("Signal averaged on event-by-event basis High Gain Area Idx "); + InitHists(hist,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j); + } + } + } + + if (fAverageLoGainAreas->GetEntries()==0 && IsLoGain() ) + { + fAverageLoGainAreas->Expand(nareas); + + for (Int_t j=0; jInheritsFrom("MGeomCamMagic")) + { + hist.GetHGausHist()->SetTitle(Form("%s%s%s","Rel. Arr. Times averaged on event-by-event basis ", + j==0 ? "Inner Pixels " : "Outer Pixels ","Low Gain Runs: ")); + hist.InitBins(); + hist.SetEventFrequency(fPulserFrequency); + } + else + { + hist.GetHGausHist()->SetTitle("Rel. Arr. Times averaged on event-by-event basis Low Gain Area Idx "); + InitHists(hist,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j); + } + } + } + + if (fAverageHiGainSectors->GetEntries()==0) + { + fAverageHiGainSectors->Expand(nsectors); + + for (Int_t j=0; jSetTitle("Rel. Arr. Times averaged on event-by-event basis HiGain Sector "); + hist.SetNbins(fAverageNbins); + + InitHists(hist,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j); + } + } + + if (fAverageLoGainSectors->GetEntries()==0 && IsLoGain() ) + { + fAverageLoGainSectors->Expand(nsectors); + + for (Int_t j=0; jSetTitle("Rel. Arr. Times averaged on event-by-event basis LoGain Sector "); + hist.SetNbins(fAverageNbins); + + InitHists(hist,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j); + } + } + + fSumareahi .Set(nareas); + fSumarealo .Set(nareas); + fSumsectorhi.Set(nsectors); + fSumsectorlo.Set(nsectors); + fNumareahi .Set(nareas); + fNumarealo .Set(nareas); + fNumsectorhi.Set(nsectors); + fNumsectorlo.Set(nsectors); + + return kTRUE; +} + + +// ------------------------------------------------------------------------------- +// +// Retrieves pointer to MArrivalTimeCam: +// +// Retrieves from MGeomCam: +// - number of pixels +// - number of pixel areas +// - number of sectors +// +// Fills HiGain or LoGain histograms (MHGausEvents::FillHistAndArray()), respectively +// depending on MArrivalTimePix::IsLoGainUsed(), with: +// - MArrivalTimePix::GetArrivalTime(pixid) - MArrivalTimePix::GetArrivalTime(1); +// (i.e. the time difference between pixel i and pixel 1 (hardware number: 2) ) +// +Bool_t MHCalibrationRelTimeCam::FillHists(const MParContainer *par, const Stat_t w) +{ + + MArrivalTimeCam *arrtime = (MArrivalTimeCam*)par; + if (!arrtime) + { + gLog << err << "No argument in MArrivalTime::Fill... abort." << endl; + return kFALSE; + } + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nareas = fGeom->GetNumAreas(); + const Int_t nsectors = fGeom->GetNumSectors(); + + fSumareahi .Reset(); + fSumarealo .Reset(); + fSumsectorhi.Reset(); + fSumsectorlo.Reset(); + fNumareahi .Reset(); + fNumarealo .Reset(); + fNumsectorhi.Reset(); + fNumsectorlo.Reset(); + + const MArrivalTimePix &refpix = (*arrtime)[fReferencePixel]; + const Float_t reftime = refpix.IsLoGainUsed() + ? refpix.GetArrivalTimeLoGain() : refpix.GetArrivalTimeHiGain(); + + for (Int_t i=0; iGetSize(); i++) + { + + MHCalibrationRelTimePix &histhi = (MHCalibrationRelTimePix&)(*this)[i]; + + if (histhi.IsExcluded()) + continue; + + MCalibrationRelTimePix &pix = fIntensCam + ? (MCalibrationRelTimePix&)(*fIntensCam)[i] + : (MCalibrationRelTimePix&)(*fCam)[i]; + + if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + pix.SetHiGainSaturation(); + histhi.SetExcluded(); + } + else + if (IsLoGain()) + (*this)(i).SetExcluded(); + + Stat_t overflow = histhi.GetHGausHist()->GetBinContent(histhi.GetHGausHist()->GetNbinsX()+1); + if (overflow > 0.1) + { + *fLog << warn << GetDescriptor() + << ": HiGain Histogram Overflow occurred " << overflow + << " times in pixel: " << i << " (without saturation!) " << endl; + // bad.SetUncalibrated( MBadPixelsPix::kHiGainOverFlow ); + } + + overflow = histhi.GetHGausHist()->GetBinContent(0); + if (overflow > 0.1) + { + *fLog << warn << GetDescriptor() + << ": HiGain Histogram Underflow occurred " << overflow + << " times in pixel: " << i << " (without saturation!) " << endl; + // bad.SetUncalibrated( MBadPixelsPix::kHiGainOverFlow ); + } + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationRelTimePix &histhi = (MHCalibrationRelTimePix&)GetAverageHiGainArea(j); + + if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + MCalibrationRelTimePix &pix = fIntensCam + ? (MCalibrationRelTimePix&)fIntensCam->GetAverageArea(j) + : (MCalibrationRelTimePix&)fCam->GetAverageArea(j); + pix.SetHiGainSaturation(); + histhi.SetExcluded(); + } + else + if (IsLoGain()) + GetAverageLoGainArea(j).SetExcluded(); + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationRelTimePix &histhi = (MHCalibrationRelTimePix&)GetAverageHiGainSector(j); + MCalibrationRelTimePix &pix = fIntensCam + ? (MCalibrationRelTimePix&)fIntensCam->GetAverageSector(j) + : (MCalibrationRelTimePix&)fCam->GetAverageSector(j); + + if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + pix.SetHiGainSaturation(); + histhi.SetExcluded(); + } + else + if (IsLoGain()) + GetAverageLoGainSector(j).SetExcluded(); + } + + FitHiGainArrays(fIntensCam ? (MCalibrationCam&)(*fIntensCam->GetCam()) : (MCalibrationCam&)(*fCam), + *fBadPixels, + MBadPixelsPix::kRelTimeNotFitted, + MBadPixelsPix::kRelTimeOscillating); + + if (IsLoGain()) + FitLoGainArrays(fIntensCam ? (MCalibrationCam&)(*fIntensCam->GetCam()) : (MCalibrationCam&)(*fCam), + *fBadPixels, + MBadPixelsPix::kRelTimeNotFitted, + MBadPixelsPix::kRelTimeOscillating); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set: +// - MBadPixelsPix::kRelTimeNotFitted +// - MBadPixelsPix::kRelTimeOscillating +// +void MHCalibrationRelTimeCam::FinalizeBadPixels() +{ + + for (Int_t i=0; iGetSize(); i++) + { + + MBadPixelsPix &bad = (*fBadPixels)[i]; + + if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeNotFitted )) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeOscillating)) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + } +} + +// -------------------------------------------------------------------------- +// +// The types are as follows: +// +// Fitted values: +// ============== +// +// 0: Fitted Mean Relative Arrival Time in FADC slices (MHGausEvents::GetMean()*MHCalibrationRelTimePix::GetFADCSliceWidth()) +// 1: Error Mean Relative Arrival Time in FADC slices (MHGausEvents::GetMeanErr()*MHCalibrationRelTimePix::GetFADCSliceWidth()) +// 2: Sigma fitted Relative Arrival Time in FADC slices (MHGausEvents::GetSigma()*MHCalibrationRelTimePix::GetFADCSliceWidth()) +// 3: Error Sigma Relative Arrival Time in FADC slices (MHGausEvents::GetSigmaErr()*MHCalibrationRelTimePix::GetFADCSliceWidth()) +// +// 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 MHCalibrationRelTimeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + + if (fHiGainArray->GetSize() <= idx) + return kFALSE; + + const MHCalibrationRelTimePix &pix = (MHCalibrationRelTimePix&)(*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 MHCalibrationPix::DrawClone() for pixel idx +// +void MHCalibrationRelTimeCam::DrawPixelContent(Int_t idx) const +{ + (*this)[idx].DrawClone(); +} Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimeCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimeCam.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimeCam.h (revision 4929) @@ -0,0 +1,57 @@ +#ifndef MARS_MHCalibrationRelTimeCam +#define MARS_MHCalibrationRelTimeCam + +#ifndef MARS_MHCalibrationCam +#include "MHCalibrationCam.h" +#endif + +#ifndef MARS_MArrayI +#include "MArrayI.h" +#endif + +#ifndef MARS_MArrayD +#include "MArrayD.h" +#endif + +class MGeomCam; +class MHCalibrationRelTimeCam : public MHCalibrationCam +{ + +private: + + static const Float_t fgNumHiGainSaturationLimit; //! The default for fNumHiGainSaturationLimit (now at: 0.25) + static const UInt_t fgReferencePixel; //! Default for fReferencePixel (now set to: 1) + + MArrayD fSumareahi ; //! + MArrayD fSumarealo ; //! + MArrayD fSumsectorhi; //! + MArrayD fSumsectorlo; //! + MArrayI fNumareahi ; //! + MArrayI fNumarealo ; //! + MArrayI fNumsectorhi; //! + MArrayI fNumsectorlo; //! + + UInt_t fReferencePixel; // The reference pixel for rel. times + + Bool_t ReInitHists(MParList *pList); + Bool_t FillHists(const MParContainer *par, const Stat_t w=1); + Bool_t FinalizeHists(); + void FinalizeBadPixels(); + +public: + + MHCalibrationRelTimeCam(const char *name=NULL, const char *title=NULL); + ~MHCalibrationRelTimeCam() {} + + UInt_t GetReferencePixel() const { return fReferencePixel; } + void SetReferencePixel(const UInt_t i=fgReferencePixel) { fReferencePixel = i; } + + + Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; + void DrawPixelContent(Int_t idx) const; + + ClassDef(MHCalibrationRelTimeCam, 1) // Histogram class for Relative Time Camera Calibration +}; + +#endif + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimePix.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimePix.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimePix.cc (revision 4929) @@ -0,0 +1,80 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationRelTimePix +// +// Histogram class for the relative arrival time calibration. +// Stores and fits the relative arrival times between pixel fPixId and +// pixel number 1 (hardware index: 2). Times are taken from MArrivalTimePix +// +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationRelTimePix.h" + +#include + +ClassImp(MHCalibrationRelTimePix); + +using namespace std; +// +const Int_t MHCalibrationRelTimePix::fgRelTimeNbins = 900; +const Axis_t MHCalibrationRelTimePix::fgRelTimeFirst = -5.; +const Axis_t MHCalibrationRelTimePix::fgRelTimeLast = 5.; +const Float_t MHCalibrationRelTimePix::fgFADCSliceWidth = 3.3333; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - the default number for fNbins (fgRelTimeNbins) +// - the default number for fFirst (fgRelTimeFirst) +// - the default number for fLast (fgRelTimeLast) +// +// - the default name of the fHGausHist ("HCalibrationRelTime") +// - the default title of the fHGausHist ("Distribution of Relative Arrival Times Pixel ") +// - the default x-axis title for fHGausHist ("FADC Slice") +// - the default y-axis title for fHGausHist ("Nr. of events") +// +// - the default number for fFADCSliceWidth (fgFADCSliceWidth) +// +MHCalibrationRelTimePix::MHCalibrationRelTimePix(const char *name, const char *title) +{ + + fName = name ? name : "MHCalibrationRelTimePix"; + fTitle = title ? title : "Histogrammed Calibration Relative Arrival Time events"; + + SetNbins ( fgRelTimeNbins ); + SetFirst ( fgRelTimeFirst ); + SetLast ( fgRelTimeLast ); + SetFADCSliceWidth(); + + // Create a large number of bins, later we will rebin + fHGausHist.SetName("HCalibrationRelTime"); + fHGausHist.SetTitle("Distribution of Relative Arrival Times Pixel "); + fHGausHist.SetXTitle("FADC Slice"); + fHGausHist.SetYTitle("Nr. of events"); + +} + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimePix.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimePix.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationRelTimePix.h (revision 4929) @@ -0,0 +1,34 @@ +#ifndef MARS_MHCalibrationRelTimePix +#define MARS_MHCalibrationRelTimePix + +#ifndef MARS_MHCalibrationPix +#include "MHCalibrationPix.h" +#endif + +class MHCalibrationRelTimePix : public MHCalibrationPix +{ + +private: + + static const Int_t fgRelTimeNbins; //! Default for fNbins (now set to: 900 ) + static const Axis_t fgRelTimeFirst; //! Default for fFirst (now set to: -13.5 ) + static const Axis_t fgRelTimeLast; //! Default for fLast (now set to: 13.5 ) + static const Float_t fgFADCSliceWidth; //! Default for fFADCSliceWidth (now set to: 3.333) + + Float_t fFADCSliceWidth; // Time FADC Slice Width in ns. + +public: + + MHCalibrationRelTimePix(const char *name=NULL, const char *title=NULL); + ~MHCalibrationRelTimePix() {} + + // Getters + Float_t GetFADCSliceWidth() const { return fFADCSliceWidth; } + + // Setters + void SetFADCSliceWidth( const Float_t f=fgFADCSliceWidth ) { fFADCSliceWidth = f; } + + ClassDef(MHCalibrationRelTimePix, 1) // Histogram class for Relative Time Pixel Calibration +}; + +#endif Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestCam.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestCam.cc (revision 4929) @@ -0,0 +1,510 @@ +/* ======================================================================== *\ +! +! * +! * 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 +// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit()) +// In case this does not make the fit valid, the histogram means and RMS's are +// taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas +// from the mean are counted as Pickup events (stored in MHCalibrationPix::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 "MHCalibrationPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MCalibrationTestCam.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"; + + SetAverageNbins(5000); + +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// - MCalibrationTestCam +// +// Searches pointer to: +// - MCerPhotEvt +// +// Initializes, if empty to MGeomCam::GetNumAreas() for: +// - MHCalibrationCam::fAverageHiGainAreas +// +// Initializes, if empty to MGeomCam::GetNumSectors() for: +// - MHCalibrationCam::fAverageHiGainSectors +// +// 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) + { + 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); + + TH1F *h = GetAverageHiGainArea(j).GetHGausHist(); + h->SetTitle( Form("%s%s", h->GetTitle()," Runs: ")); + } + } + + + 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); + + TH1F *h = GetAverageHiGainSector(j).GetHGausHist(); + h->SetTitle( Form("%s%s", h->GetTitle()," Runs: ")); + + } + } + + + fMeanMeanPhotPerArea.Set(nareas); + fRmsMeanPhotPerArea .Set(nareas); + fMeanSigmaPhotPerArea.Set(nareas); + fRmsSigmaPhotPerArea.Set(nareas); + + SetLoGain(kFALSE); + + return kTRUE; +} + + +// ------------------------------------------------------------------------------- +// +// Retrieves pointer to MCerPhotEvt: +// +// Retrieves from MGeomCam: +// - number of pixels +// - number of pixel areas +// - number of sectors +// +// Fills HiGain histograms (MHGausEvents::FillHistAndArray()) +// with: +// - MCerPhotPix::GetNumPhotons(pixid); +// +Bool_t MHCalibrationTestCam::FillHists(const MParContainer *par, const Stat_t w) +{ + + MCerPhotEvt *calibration = (MCerPhotEvt*)par; + if (!calibration) + { + gLog << err << "No argument in MHCalibrationTestCam::Fill... abort." << endl; + return kFALSE; + } + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nareas = fGeom->GetNumAreas(); + const Int_t nsectors = fGeom->GetNumSectors(); + + TArrayF sumareahi (nareas); + TArrayF sumsectorhi(nsectors); + TArrayI numareahi (nareas); + TArrayI numsectorhi(nsectors); + + for (Int_t i=0; iGetPixById(i); + if (!pix) + continue; + + + const Float_t signal = pix->GetNumPhotons(); + + if (signal < 0.0001) + continue; + + 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; jGetNumAreas()); + + for (Int_t i=0; iGetSize(); i++) + { + + MHCalibrationPix &hist = (*this)[i]; + + if (hist.IsEmpty()) + continue; + + if (!hist.FitGaus()) + if (!hist.RepeatFit()) + hist.BypassFit(); + + hist.CreateFourierSpectrum(); + + const Float_t area = (*fGeom)[i].GetA(); + const Int_t aidx = (*fGeom)[i].GetAidx(); + + fMeanMeanPhotPerArea[aidx] += hist.GetMean() / area; + fRmsMeanPhotPerArea [aidx] += hist.GetMean() / area * hist.GetMean() / area; + fMeanSigmaPhotPerArea[aidx] += hist.GetSigma()/ area; + fRmsSigmaPhotPerArea [aidx] += hist.GetSigma()/ area * hist.GetSigma() / area; + numaidx[aidx]++; + } + + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationPix &hist = GetAverageHiGainArea(j); + if (hist.IsEmpty()) + continue; + + if (!hist.FitGaus()) + if (!hist.RepeatFit()) + hist.BypassFit(); + + hist.CreateFourierSpectrum(); + + fRmsMeanPhotPerArea [j] -= fMeanMeanPhotPerArea [j]*fMeanMeanPhotPerArea [j]/numaidx[j]; + fRmsSigmaPhotPerArea[j] -= fMeanSigmaPhotPerArea[j]*fMeanSigmaPhotPerArea[j]/numaidx[j]; + + fMeanMeanPhotPerArea [j] /= numaidx[j]; + fMeanSigmaPhotPerArea[j] /= numaidx[j]; + fRmsMeanPhotPerArea [j] /= numaidx[j]-1.; + fRmsSigmaPhotPerArea [j] /= numaidx[j]-1.; + + if (fRmsMeanPhotPerArea [j] > 0.) + fRmsMeanPhotPerArea [j] = TMath::Sqrt(fRmsMeanPhotPerArea [j]); + if (fRmsSigmaPhotPerArea [j] > 0.) + fRmsSigmaPhotPerArea [j] = TMath::Sqrt(fRmsSigmaPhotPerArea [j]); + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationPix &hist = GetAverageHiGainSector(j); + if (hist.IsEmpty()) + continue; + + if (!hist.FitGaus()) + if (!hist.RepeatFit()) + hist.BypassFit(); + + hist.CreateFourierSpectrum(); + } + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// 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()) +// +// Converted values: +// ================= +// +// 7: Fitted Mean Test Calibration (MHGausEvents::GetMean()) by MGeomPix::GetA() +// 8: Fitted Mean Error Calibration (MHGausEvents::GetMeanErr()) by MGeomPix::GetA() +// 9: Fitted Sigma Test Calibration (MHGausEvents::GetSigma()) by MGeomPix::GetA() +// 10: Fitted Sigma Error Calibration (MHGausEvents::GetSigmaErr()) by MGeomPix::GetA() +// +Bool_t MHCalibrationTestCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + + if (fHiGainArray->GetSize() <= idx) + return kFALSE; + + const MHCalibrationPix &pix = (*this)[idx]; + + if (pix.IsEmpty()) + return kFALSE; + + 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; + case 7: + val = pix.GetMean()/cam[idx].GetA(); + break; + case 8: + val = pix.GetMeanErr()/cam[idx].GetA(); + break; + case 9: + val = pix.GetSigma()/cam[idx].GetA(); + break; + case 10: + val = pix.GetSigmaErr()/cam[idx].GetA(); + break; + default: + return kFALSE; + } + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calls MHCalibrationPix::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++) + { + + MHCalibrationPix &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/mhcalib/MHCalibrationTestCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestCam.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestCam.h (revision 4929) @@ -0,0 +1,48 @@ +#ifndef MARS_MHCalibrationTestCam +#define MARS_MHCalibrationTestCam + +#ifndef MARS_MHCalibrationCam +#include "MHCalibrationCam.h" +#endif + +#ifndef ROOT_TArrayF +#include +#endif +#ifndef ROOT_TArrayI +#include +#endif + +class MHCalibrationTestCam : public MHCalibrationCam +{ + +private: + + TArrayF fMeanMeanPhotPerArea; + TArrayF fRmsMeanPhotPerArea ; + TArrayF fMeanSigmaPhotPerArea; + TArrayF fRmsSigmaPhotPerArea ; + + Bool_t ReInitHists(MParList *pList); + Bool_t FillHists(const MParContainer *par, const Stat_t w=1); + Bool_t FinalizeHists(); + +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; + + const Float_t GetMeanMeanPhotPerArea ( const Int_t aidx ) const { return fMeanMeanPhotPerArea [aidx]; } + const Float_t GetMeanSigmaPhotPerArea ( const Int_t aidx ) const { return fMeanSigmaPhotPerArea [aidx]; } + const Float_t GetRmsMeanPhotPerArea ( const Int_t aidx ) const { return fRmsMeanPhotPerArea [aidx]; } + const Float_t GetRmsSigmaPhotPerArea ( const Int_t aidx ) const { return fRmsSigmaPhotPerArea [aidx]; } + + void CalcAverageSigma(); + + ClassDef(MHCalibrationTestCam, 1) // Histogram class for Relative Time Camera Calibration +}; + +#endif + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestPix.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestPix.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestPix.cc (revision 4929) @@ -0,0 +1,98 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +////////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationTestPix +// +// Histogram class for the charge calibration. +// Stores and fits the charges and stores the location of the maximum FADC +// slice. Tests are taken from MExtractedSignalPix. +// +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationTestPix.h" + +#include +#include + +#include +#include +#include +#include + +#include "MH.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHCalibrationTestPix); + +using namespace std; + +const Int_t MHCalibrationTestPix::fgChargeNbins = 4000; +const Axis_t MHCalibrationTestPix::fgChargeFirst = -0.5; +const Axis_t MHCalibrationTestPix::fgChargeLast = 39999.5; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - the default number for fNbins (fgChargeNbins) +// - the default number for fFirst (fgChargeFirst) +// - the default number for fLast (fgChargeLast) +// +// - the default name of the fHGausHist ("HCalibrationTest") +// - the default title of the fHGausHist ("Distribution of calibrated FADC slices Pixel ") +// - the default x-axis title for fHGausHist ("Sum FADC Slices") +// - the default y-axis title for fHGausHist ("Nr. of events") +// +// Calls: +// - Clear(); +// +MHCalibrationTestPix::MHCalibrationTestPix(const char *name, const char *title) +{ + + fName = name ? name : "MHCalibrationTestPix"; + fTitle = title ? title : "Statistics of the calibrated FADC sums of calibration events"; + + SetNbins ( fgChargeNbins ); + SetFirst ( fgChargeFirst ); + SetLast ( fgChargeLast ); + + fHGausHist.SetName("HCalibrationTest"); + fHGausHist.SetTitle("Distribution of calibrated Photons Pixel "); + fHGausHist.SetXTitle("Nr. Photons"); + fHGausHist.SetYTitle("Nr. of events"); + +} + + +// -------------------------------------------------------------------------- +// +// returns fHGausHist.Integral("width") +// +const Float_t MHCalibrationTestPix::GetIntegral() const +{ + return fHGausHist.Integral("width"); +} + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestPix.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestPix.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestPix.h (revision 4929) @@ -0,0 +1,27 @@ +#ifndef MARS_MHCalibrationTestPix +#define MARS_MHCalibrationTestPix + +#ifndef MARS_MHCalibrationPix +#include "MHCalibrationPix.h" +#endif + +class MHCalibrationTestPix : public MHCalibrationPix +{ + +private: + + static const Int_t fgChargeNbins; // Default for fNBins (now set to: 2000 ) + static const Axis_t fgChargeFirst; // Default for fFirst (now set to: -0.5 ) + static const Axis_t fgChargeLast; // Default for fLast (now set to: 1999.5) + +public: + + MHCalibrationTestPix(const char *name=NULL, const char *title=NULL); + ~MHCalibrationTestPix() {} + + const Float_t GetIntegral() const; + + ClassDef(MHCalibrationTestPix, 1) // Base Histogram class for Test Pixel Calibration +}; + +#endif Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimeCam.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimeCam.cc (revision 4929) @@ -0,0 +1,462 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationTestTimeCam +// +// Fills the calibrated signal from an MArrivalTime into +// MHCalibrationTestTimePix 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 +// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit()) +// In case this does not make the fit valid, the histogram means and RMS's are +// taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas +// from the mean are counted as Pickup events (stored in MHCalibrationPix::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 "MHCalibrationTestTimeCam.h" +#include "MHCalibrationTestTimePix.h" + +#include "MHCalibrationPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MCalibrationCam.h" +#include "MCalibrationPix.h" + +#include "MArrivalTime.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" + +ClassImp(MHCalibrationTestTimeCam); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHCalibrationTestTimeCam::MHCalibrationTestTimeCam(const char *name, const char *title) +{ + + fName = name ? name : "MHCalibrationTestTimeCam"; + fTitle = title ? title : "Histogram class for testing the calibration of arrival times"; + + SetAverageNbins(5000); +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// +// Searches pointer to: +// - MArrivalTime +// +// Initializes, if empty to MGeomCam::GetNumAreas() for: +// - MHCalibrationCam::fAverageHiGainAreas +// +// Initializes, if empty to MGeomCam::GetNumSectors() for: +// - MHCalibrationCam::fAverageHiGainSectors +// +// 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 MHCalibrationTestTimeCam::ReInitHists(MParList *pList) +{ + + MArrivalTime *signal = (MArrivalTime*)pList->FindObject("MArrivalTime"); + if (!signal) + { + *fLog << err << "MArrivalTime 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) + { + fAverageHiGainAreas->Expand(nareas); + + for (Int_t j=0; jSetTitle("TestTime Calibrations Area Idx "); + GetAverageHiGainArea(j).SetNbins(fAverageNbins); + GetAverageHiGainArea(j).InitBins(); + GetAverageHiGainArea(j).ChangeHistId(j); + GetAverageHiGainArea(j).SetEventFrequency(fPulserFrequency); + + TH1F *h = GetAverageHiGainArea(j).GetHGausHist(); + h->SetTitle( Form("%s%s", h->GetTitle()," Runs: ")); + + + } + } + + + if (fAverageHiGainSectors->GetEntries()==0) + { + fAverageHiGainSectors->Expand(nsectors); + + for (Int_t j=0; jSetTitle("TestTime Calibrations Sector "); + GetAverageHiGainSector(j).SetNbins(fAverageNbins); + GetAverageHiGainSector(j).InitBins(); + GetAverageHiGainSector(j).ChangeHistId(j); + GetAverageHiGainSector(j).SetEventFrequency(fPulserFrequency); + TH1F *h = GetAverageHiGainSector(j).GetHGausHist(); + h->SetTitle( Form("%s%s", h->GetTitle()," Runs: ")); + + + } + } + + SetLoGain(kFALSE); + + return kTRUE; +} + + +// ------------------------------------------------------------------------------- +// +// Retrieves pointer to MArrivalTime: +// +// Retrieves from MGeomCam: +// - number of pixels +// - number of pixel areas +// - number of sectors +// +// Fills HiGain histograms (MHGausEvents::FillHistAndArray()) +// with: +// - MArrivalTime::GetArrivalTime(pixid) - MArrivalTime::GetArrivalTime(1); +// (i.e. the time difference between pixel i and pixel 1 (hardware number: 2) ) +// +Bool_t MHCalibrationTestTimeCam::FillHists(const MParContainer *par, const Stat_t w) +{ + + MArrivalTime *calibration = (MArrivalTime*)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(); + + TArrayF sumareahi (nareas); + TArrayF sumsectorhi(nsectors); + TArrayI numareahi (nareas); + TArrayI numsectorhi(nsectors); + + for (Int_t i=0; iGetSize(); i++) + { + + MHCalibrationPix &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++) + { + + MHCalibrationPix &hist = GetAverageHiGainArea(j); + if (hist.IsEmpty()) + continue; + + if (!hist.FitGaus()) + if (!hist.RepeatFit()) + { + hist.BypassFit(); + } + + hist.CreateFourierSpectrum(); + + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationPix &hist = GetAverageHiGainSector(j); + if (hist.IsEmpty()) + continue; + + if (!hist.FitGaus()) + if (!hist.RepeatFit()) + { + hist.BypassFit(); + } + hist.CreateFourierSpectrum(); + + + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +void MHCalibrationTestTimeCam::FinalizeBadPixels() +{ + +} + +// -------------------------------------------------------------------------- +// +// The types are as follows: +// +// Fitted values: +// ============== +// +// 0: Fitted Mean Time Calibration (MHGausEvents::GetMean()) +// 1: Error Mean Time Calibration (MHGausEvents::GetMeanErr()) +// 2: Sigma fitted Time Calibration (MHGausEvents::GetSigma()) +// 3: Error Sigma Time 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 MHCalibrationTestTimeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + + if (fHiGainArray->GetSize() <= idx) + return kFALSE; + + const MHCalibrationPix &pix = (*this)[idx]; + + if (pix.IsExcluded()) + return kFALSE; + + 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 MHCalibrationPix::DrawClone() for pixel idx +// +void MHCalibrationTestTimeCam::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 MHCalibrationTestTimeCam::CalcAverageSigma() +{ + + for (UInt_t j=0; jGetNumAreas(); j++) + { + + MHCalibrationPix &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/mhcalib/MHCalibrationTestTimeCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimeCam.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimeCam.h (revision 4929) @@ -0,0 +1,32 @@ +#ifndef MARS_MHCalibrationTestTimeCam +#define MARS_MHCalibrationTestTimeCam + +#ifndef MARS_MHCalibrationCam +#include "MHCalibrationCam.h" +#endif + +class MHCalibrationTestTimeCam : 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: + + MHCalibrationTestTimeCam(const char *name=NULL, const char *title=NULL); + ~MHCalibrationTestTimeCam() {} + + 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(MHCalibrationTestTimeCam, 0) // Histogram class for Relative Time Camera Calibration +}; + +#endif + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimePix.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimePix.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimePix.cc (revision 4929) @@ -0,0 +1,78 @@ +/* ======================================================================== *\ +! +! * +! * 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 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ +////////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationTestTimePix +// +// Histogram class for the charge calibration. +// Stores and fits the charges and stores the location of the maximum FADC +// slice. TestTimes are taken from MExtractedSignalPix. +// +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationTestTimePix.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHCalibrationTestTimePix); + +using namespace std; + +const Int_t MHCalibrationTestTimePix::fgChargeNbins = 600; +const Axis_t MHCalibrationTestTimePix::fgChargeFirst = -0.5; +const Axis_t MHCalibrationTestTimePix::fgChargeLast = 29.5; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - the default number for fNbins (fgChargeNbins) +// - the default number for fFirst (fgChargeFirst) +// - the default number for fLast (fgChargeLast) +// +// - the default name of the fHGausHist ("HCalibrationTestTime") +// - the default title of the fHGausHist ("Distribution of calibrated FADC slices Pixel ") +// - the default x-axis title for fHGausHist ("Sum FADC Slices") +// - the default y-axis title for fHGausHist ("Nr. of events") +// +// Calls: +// - Clear(); +// +MHCalibrationTestTimePix::MHCalibrationTestTimePix(const char *name, const char *title) +{ + + fName = name ? name : "MHCalibrationTestTimePix"; + fTitle = title ? title : "Statistics of the calibrated FADC sums of calibration events"; + + SetNbins ( fgChargeNbins ); + SetFirst ( fgChargeFirst ); + SetLast ( fgChargeLast ); + + fHGausHist.SetName("HCalibrationTestTime"); + fHGausHist.SetTitle("Distribution of calibrated Arrival Times Pixel"); + fHGausHist.SetXTitle("Arrival Time [FADC units]"); + fHGausHist.SetYTitle("Nr. of events"); + +} + Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimePix.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimePix.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimePix.h (revision 4929) @@ -0,0 +1,26 @@ +#ifndef MARS_MHCalibrationTestTimePix +#define MARS_MHCalibrationTestTimePix + +#ifndef MARS_MHCalibrationPix +#include "MHCalibrationPix.h" +#endif + +class TH1F; +class MHCalibrationTestTimePix : public MHCalibrationPix +{ + +private: + + static const Int_t fgChargeNbins; // Default for fNBins (now set to: 2000 ) + static const Axis_t fgChargeFirst; // Default for fFirst (now set to: -0.5 ) + static const Axis_t fgChargeLast; // Default for fLast (now set to: 1999.5) + +public: + + MHCalibrationTestTimePix(const char *name=NULL, const char *title=NULL); + ~MHCalibrationTestTimePix() {} + + ClassDef(MHCalibrationTestTimePix, 1) // Base Histogram class for TestTime Pixel Calibration +}; + +#endif Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTimeTestCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTimeTestCam.cc (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTimeTestCam.cc (revision 4929) @@ -0,0 +1,512 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationTimeTestCam +// +// Fills the calibrated signal from an MArrivalTime into +// MHCalibrationTimeTestPix 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 "MHCalibrationTimeTestCam.h" +#include "MHCalibrationTimeTestPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MCalibrationCam.h" +#include "MCalibrationPix.h" + +#include "MArrivalTime.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" + +ClassImp(MHCalibrationTimeTestCam); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHCalibrationTimeTestCam::MHCalibrationTimeTestCam(const char *name, const char *title) +{ + + fName = name ? name : "MHCalibrationTimeTestCam"; + fTitle = title ? title : "Histogram class for testing the calibration of arrival times"; + +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// +// Searches pointer to: +// - MArrivalTime +// +// 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 MHCalibrationTimeTestCam::ReInitHists(MParList *pList) +{ + + MArrivalTime *signal = (MArrivalTime*)pList->FindObject("MArrivalTime"); + if (!signal) + { + *fLog << err << "MArrivalTime 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("TimeTest 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("TimeTest 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("TimeTest 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("TimeTest Calibrations Sector "); + GetAverageLoGainSector(j).SetNbins(fAverageNbins); + GetAverageLoGainSector(j).InitBins(); + GetAverageLoGainSector(j).ChangeHistId(j); + GetAverageLoGainSector(j).SetEventFrequency(fPulserFrequency); + } + } + + + return kTRUE; +} + + +// ------------------------------------------------------------------------------- +// +// Retrieves pointer to MArrivalTime: +// +// 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 MHCalibrationTimeTestCam::FillHists(const MParContainer *par, const Stat_t w) +{ + + MArrivalTime *calibration = (MArrivalTime*)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; iGetSize(); 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::kTimeTestNotFitted +// - MBadPixelsPix::kTimeTestOscillating +// +void MHCalibrationTimeTestCam::FinalizeBadPixels() +{ + +} + +// -------------------------------------------------------------------------- +// +// The types are as follows: +// +// Fitted values: +// ============== +// +// 0: Fitted Mean Time Calibration (MHGausEvents::GetMean()) +// 1: Error Mean Time Calibration (MHGausEvents::GetMeanErr()) +// 2: Sigma fitted Time Calibration (MHGausEvents::GetSigma()) +// 3: Error Sigma Time 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 MHCalibrationTimeTestCam::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]; + + if (pix.IsExcluded()) + return kFALSE; + + 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 MHCalibrationTimeTestCam::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 MHCalibrationTimeTestCam::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/mhcalib/MHCalibrationTimeTestCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTimeTestCam.h (revision 4929) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTimeTestCam.h (revision 4929) @@ -0,0 +1,33 @@ +#ifndef MARS_MHCalibrationTimeTestCam +#define MARS_MHCalibrationTimeTestCam + +#ifndef MARS_MHCalibrationCam +#include "MHCalibrationCam.h" +#endif + +class MGeomCam; +class MHCalibrationTimeTestCam : 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: + + MHCalibrationTimeTestCam(const char *name=NULL, const char *title=NULL); + ~MHCalibrationTimeTestCam() {} + + 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(MHCalibrationTimeTestCam, 0) // Histogram class for Relative Time Camera Calibration +}; + +#endif +