Index: /trunk/MagicSoft/Mars/macros/dohtml.C =================================================================== --- /trunk/MagicSoft/Mars/macros/dohtml.C (revision 4929) +++ /trunk/MagicSoft/Mars/macros/dohtml.C (revision 4930) @@ -71,4 +71,5 @@ sourcedir += "mhist:"; sourcedir += "mhistmc:"; + sourcedir += "mhcalib:"; sourcedir += "mhvstime:"; sourcedir += "mimage:"; Index: /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h (revision 4929) +++ /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h (revision 4930) @@ -37,18 +37,4 @@ #pragma link C++ class MCalibrationChargePINDiode+; -#pragma link C++ class MHCalibrationCam+; -#pragma link C++ class MHCalibrationPix+; -#pragma link C++ class MHCalibrationChargeCam+; -#pragma link C++ class MHCalibrationChargePix+; -#pragma link C++ class MHCalibrationChargeBlindCam+; -#pragma link C++ class MHCalibrationChargeBlindPix+; -#pragma link C++ class MHCalibrationChargePINDiode+; -#pragma link C++ class MHCalibrationRelTimePix+; -#pragma link C++ class MHCalibrationRelTimeCam+; -#pragma link C++ class MHCalibrationTestCam+; -#pragma link C++ class MHCalibrationTestPix+; -#pragma link C++ class MHCalibrationTestTimeCam+; -#pragma link C++ class MHCalibrationTestTimePix+; - #pragma link C++ class MMcCalibrationCalc++; Index: unk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc (revision 4929) +++ (revision ) @@ -1,1094 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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), - fDebug(kFALSE), fLoGain(kFALSE) -{ - - 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(); - -} - -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(); - fLoGain = 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 (fLoGain) - 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 (fLoGain) - 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 (fLoGain) - 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 (fLoGain) - 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 (fLoGain) - 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 (fLoGain) - 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 - // - 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 - // - 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(fLoGain ? 2 : 1,nareas); - - for (Int_t i=0; icd(fLoGain ? 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 (fLoGain) - { - 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; - } - - return rc; -} Index: unk/MagicSoft/Mars/mcalib/MHCalibrationCam.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h (revision 4929) +++ (revision ) @@ -1,181 +1,0 @@ -#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 - Bool_t fDebug; // Debug option is used? - Bool_t fLoGain; // Is Low-Gain available - - 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); - - Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); - -public: - - MHCalibrationCam(const char *name=NULL, const char *title=NULL); - virtual ~MHCalibrationCam(); - - virtual Bool_t SetupFill(const MParList *pList); - virtual Bool_t ReInit ( MParList *pList); - virtual Bool_t Fill (const MParContainer *par, const Stat_t w=1); - virtual Bool_t Finalize ( ); - - virtual void ResetHists(); - - // Clone - virtual TObject *Clone(const char *name="") const; - - // Draw - virtual void Draw(const Option_t *opt); - - virtual Bool_t GetPixelContent ( Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; - virtual 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; - - const Bool_t IsDebug () const { return fDebug; } - - 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 ) { fDebug = b; } - void SetLoGain ( const Bool_t b=kTRUE ) { fLoGain = b; } - 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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindCam.cc (revision 4929) +++ (revision ) @@ -1,333 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindCam.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindCam.h (revision 4929) +++ (revision ) @@ -1,65 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h (revision 4929) +++ (revision ) @@ -1,522 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc (revision 4929) +++ (revision ) @@ -1,1307 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 && fLoGain) - { - 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 && fLoGain) - { - 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 && fLoGain) - { - 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(); - 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 (fLoGain) - 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 ); - 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(); - histhi.CreateFourierSpectrum(); - continue; - } - - MBadPixelsPix &bad = fIntensCam - ? fIntensCam->GetAverageBadArea(j) - : fCam->GetAverageBadArea(j); - - FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain); - } - - if (fLoGain) - 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(); - histhi.CreateFourierSpectrum(); - continue; - } - - MBadPixelsPix &bad = fIntensCam - ? fIntensCam->GetAverageBadSector(j) - : fCam->GetAverageBadSector(j); - - FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain); - } - - if (fLoGain) - 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 ); - 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 (fLoGain) - 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::kHiGainOscillating )) - bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); - - if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted )) - if (pix.IsHiGainSaturation()) - bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); - - if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating )) - if (pix.IsHiGainSaturation()) - bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); - - if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation )) - bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); - } -} - -// -------------------------------------------------------------------------- -// -// 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() && fLoGain) - { - 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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.h (revision 4929) +++ (revision ) @@ -1,104 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargePINDiode.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePINDiode.cc (revision 4929) +++ (revision ) @@ -1,398 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargePINDiode.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePINDiode.h (revision 4929) +++ (revision ) @@ -1,84 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.cc (revision 4929) +++ (revision ) @@ -1,309 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.h (revision 4929) +++ (revision ) @@ -1,61 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationPix.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationPix.cc (revision 4929) +++ (revision ) @@ -1,327 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mcalib/MHCalibrationPix.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationPix.h (revision 4929) +++ (revision ) @@ -1,60 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc (revision 4929) +++ (revision ) @@ -1,625 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 && fLoGain ) - { - 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 && fLoGain ) - { - 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 && fLoGain ) - { - 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 (fLoGain) - (*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 (fLoGain) - 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 (fLoGain) - GetAverageLoGainSector(j).SetExcluded(); - } - - FitHiGainArrays(fIntensCam ? (MCalibrationCam&)(*fIntensCam->GetCam()) : (MCalibrationCam&)(*fCam), - *fBadPixels, - MBadPixelsPix::kRelTimeNotFitted, - MBadPixelsPix::kRelTimeOscillating); - - if (fLoGain) - 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: unk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.h (revision 4929) +++ (revision ) @@ -1,57 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.cc (revision 4929) +++ (revision ) @@ -1,80 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.h (revision 4929) +++ (revision ) @@ -1,34 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.cc (revision 4929) +++ (revision ) @@ -1,510 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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); - - fLoGain = 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: unk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.h (revision 4929) +++ (revision ) @@ -1,48 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationTestTimeCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationTestTimeCam.cc (revision 4929) +++ (revision ) @@ -1,462 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: ")); - - - } - } - - fLoGain = 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: unk/MagicSoft/Mars/mcalib/MHCalibrationTestTimeCam.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationTestTimeCam.h (revision 4929) +++ (revision ) @@ -1,32 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHCalibrationTestTimePix.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationTestTimePix.cc (revision 4929) +++ (revision ) @@ -1,78 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mcalib/MHCalibrationTestTimePix.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationTestTimePix.h (revision 4929) +++ (revision ) @@ -1,26 +1,0 @@ -#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: unk/MagicSoft/Mars/mcalib/MHGausEvents.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc (revision 4929) +++ (revision ) @@ -1,863 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MHGausEvents -// -// A base class for events which are believed to follow a Gaussian distribution -// with time, e.g. calibration events, observables containing white noise, ... -// -// MHGausEvents derives from MH, thus all features of MH can be used by a class -// deriving from MHGausEvents, especially the filling functions. -// -// The central objects are: -// -// 1) The TH1F fHGausHist: -// ==================== -// -// It is created with a default name and title and resides in directory NULL. -// - Any class deriving from MHGausEvents needs to apply a binning to fHGausHist -// (e.g. by setting the variables fNbins, fFirst, fLast and calling the function -// InitBins() or by directly calling fHGausHist.SetBins(..) ) -// - The histogram is filled with the functions FillHist() or FillHistAndArray(). -// - The histogram can be fitted with the function FitGaus(). This involves stripping -// of all zeros at the lower and upper end of the histogram and re-binning to -// a new number of bins, specified in the variable fBinsAfterStripping. -// - The fit result's probability is compared to a reference probability fProbLimit -// The NDF is compared to fNDFLimit and a check is made whether results are NaNs. -// Accordingly, a flag IsGausFitOK() is set. -// - One can repeat the fit within a given amount of sigmas from the previous mean -// (specified by the variables fPickupLimit and fBlackoutLimit) with the function RepeatFit() -// - One can completely skip the fitting to set mean, sigma and its errors directly -// from the histograms with the function BypassFit() -// -// 2) The TArrayF fEvents: -// ========================== -// -// It is created with 0 entries and not expanded unless FillArray() or FillHistAndArray() -// are called. -// - A first call to FillArray() or FillHistAndArray() initializes fEvents by default -// to 512 entries. -// - Any later call to FillArray() or FillHistAndArray() fills up the array. -// Reaching the limit, the array is expanded by a factor 2. -// - The array can be fourier-transformed into the array fPowerSpectrum. -// Note that any FFT accepts only number of events which are a power of 2. -// Thus, fEvents is cut to the next power of 2 smaller than its actual number of entries. -// Be aware that you might lose information at the end of your analysis. -// - Calling the function CreateFourierSpectrum() creates the array fPowerSpectrum -// and its projection fHPowerProbability which in turn is fit to an exponential. -// - The fit result's probability is compared to a referenc probability fProbLimit -// and accordingly the flag IsExpFitOK() is set. -// - The flag IsFourierSpectrumOK() is set accordingly to IsExpFitOK(). -// Later, a closer check will be installed. -// - You can display all arrays by calls to: CreateGraphEvents() and -// CreateGraphPowerSpectrum() and successive calls to the corresponding Getters. -// -// To see an example, have a look at: Draw() -// -////////////////////////////////////////////////////////////////////////////// -#include "MHGausEvents.h" - -#include -#include -#include -#include -#include -#include -#include - -#include "MFFT.h" -#include "MArray.h" - -#include "MH.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHGausEvents); - -using namespace std; - -const Int_t MHGausEvents::fgNDFLimit = 2; -const Float_t MHGausEvents::fgProbLimit = 0.001; -const Int_t MHGausEvents::fgPowerProbabilityBins = 30; -// -------------------------------------------------------------------------- -// -// Default Constructor. -// Sets: -// - the default Probability Bins for fPowerProbabilityBins (fgPowerProbabilityBins) -// - the default Probability Limit for fProbLimit (fgProbLimit) -// - the default NDF Limit for fNDFLimit (fgNDFLimit) -// - the default number of bins after stripping for fBinsAfterStipping (fgBinsAfterStipping) -// - the default name of the fHGausHist ("HGausHist") -// - the default title of the fHGausHist ("Histogram of Events with Gaussian Distribution") -// - the default directory of the fHGausHist (NULL) -// - the default for fNbins (100) -// - the default for fFirst (0.) -// - the default for fLast (100.) -// -// Initializes: -// - fEvents to 0 entries -// - fHGausHist() -// - all other pointers to NULL -// - all variables to 0. -// - all flags to kFALSE -// -MHGausEvents::MHGausEvents(const char *name, const char *title) - : fHPowerProbability(NULL), - fPowerSpectrum(NULL), - fFGausFit(NULL), fFExpFit(NULL), - fFirst(0.), - fGraphEvents(NULL), fGraphPowerSpectrum(NULL), - fLast(100.), fNbins(100) -{ - - fName = name ? name : "MHGausEvents"; - fTitle = title ? title : "Events with expected Gaussian distributions"; - - Clear(); - - SetBinsAfterStripping(); - SetNDFLimit(); - SetPowerProbabilityBins(); - SetProbLimit(); - - fHGausHist.SetName("HGausHist"); - fHGausHist.SetTitle("Histogram of Events with Gaussian Distribution"); - // important, other ROOT will not draw the axes: - fHGausHist.UseCurrentStyle(); - fHGausHist.SetDirectory(NULL); - fHGausHist.GetYaxis()->CenterTitle(); -} - - - -// -------------------------------------------------------------------------- -// -// Default Destructor. -// -// Deletes (if Pointer is not NULL): -// -// - fHPowerProbability -// - fFGausFit -// - fFExpFit -// - fPowerSpectrum -// - fGraphEvents -// - fGraphPowerSpectrum -// -MHGausEvents::~MHGausEvents() -{ - - // delete histograms -// if (fHPowerProbability) -// if (gROOT->FindObject(fHPowerProbability->GetName())) -// delete fHPowerProbability; - - // delete fits - if (fFGausFit) - delete fFGausFit; - if (fFExpFit) - delete fFExpFit; - - // delete arrays - if (fPowerSpectrum) - delete fPowerSpectrum; - - // delete graphs - if (fGraphEvents) - delete fGraphEvents; - if (fGraphPowerSpectrum) - delete fGraphPowerSpectrum; -} - -// -------------------------------------------------------------------------- -// -// Default Clear(), can be overloaded. -// -// Sets: -// - all other pointers to NULL -// - all variables to 0. -// - all flags to kFALSE -// -// Deletes (if not NULL): -// - all pointers -// -void MHGausEvents::Clear(Option_t *o) -{ - - SetGausFitOK ( kFALSE ); - SetExpFitOK ( kFALSE ); - SetFourierSpectrumOK( kFALSE ); - SetExcluded ( kFALSE ); - - fMean = 0.; - fSigma = 0.; - fMeanErr = 0.; - fSigmaErr = 0.; - fProb = 0.; - - fCurrentSize = 0; - - if (fHPowerProbability) - { - if (gROOT->FindObject(fHPowerProbability->GetName())) - delete fHPowerProbability; - fHPowerProbability = NULL; - } - - // delete fits - if (fFGausFit) - { - delete fFGausFit; - fFGausFit = NULL; - } - - if (fFExpFit) - { - delete fFExpFit; - fFExpFit = NULL; - } - - // delete arrays - if (fPowerSpectrum) - { - delete fPowerSpectrum; - fPowerSpectrum = NULL; - } - - // delete graphs - if (fGraphEvents) - { - delete fGraphEvents; - fGraphEvents = NULL; - } - - if (fGraphPowerSpectrum) - { - delete fGraphPowerSpectrum; - fGraphPowerSpectrum = NULL; - } -} - - - -// ----------------------------------------------------------------------------- -// -// Create the x-axis for the event graph -// -Float_t *MHGausEvents::CreateEventXaxis(Int_t n) -{ - - Float_t *xaxis = new Float_t[n]; - - for (Int_t i=0;iSetXTitle("P(f)"); - fHPowerProbability->SetYTitle("Counts"); - fHPowerProbability->GetYaxis()->CenterTitle(); - fHPowerProbability->SetDirectory(NULL); - fHPowerProbability->SetBit(kCanDelete); - // - // First guesses for the fit (should be as close to reality as possible, - // - const Double_t xmax = fHPowerProbability->GetXaxis()->GetXmax(); - - fFExpFit = new TF1("FExpFit","exp([0]-[1]*x)",0.,xmax); - - const Double_t slope_guess = (TMath::Log(fHPowerProbability->GetEntries())+1.)/xmax; - const Double_t offset_guess = slope_guess*xmax; - - fFExpFit->SetParameters(offset_guess, slope_guess); - fFExpFit->SetParNames("Offset","Slope"); - fFExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess); - fFExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess); - fFExpFit->SetRange(0.,xmax); - - fHPowerProbability->Fit(fFExpFit,"RQL0"); - - if (GetExpProb() > fProbLimit) - SetExpFitOK(kTRUE); - - // - // For the moment, this is the only check, later we can add more... - // - SetFourierSpectrumOK(IsExpFitOK()); - - return; -} - -// ---------------------------------------------------------------------------------- -// -// Create a graph to display the array fEvents -// If the variable fEventFrequency is set, the x-axis is transformed into real time. -// -void MHGausEvents::CreateGraphEvents() -{ - - MArray::StripZeros(fEvents); - - const Int_t n = fEvents.GetSize(); - - if (n==0) - return; - - fGraphEvents = new TGraph(n,CreateEventXaxis(n),fEvents.GetArray()); - fGraphEvents->SetTitle("Evolution of Events with time"); - fGraphEvents->GetYaxis()->SetTitle(fHGausHist.GetXaxis()->GetTitle()); - fGraphEvents->GetYaxis()->CenterTitle(); -} - -// ---------------------------------------------------------------------------------- -// -// Create a graph to display the array fPowerSpectrum -// If the variable fEventFrequency is set, the x-axis is transformed into real frequency. -// -void MHGausEvents::CreateGraphPowerSpectrum() -{ - - MArray::StripZeros(*fPowerSpectrum); - - const Int_t n = fPowerSpectrum->GetSize(); - - fGraphPowerSpectrum = new TGraph(n,CreatePSDXaxis(n),fPowerSpectrum->GetArray()); - fGraphPowerSpectrum->SetTitle("Power Spectrum Density"); - fGraphPowerSpectrum->GetYaxis()->SetTitle("P(f)"); - fGraphPowerSpectrum->GetYaxis()->CenterTitle(); - -} - - -// ----------------------------------------------------------------------------- -// -// Create the x-axis for the event graph -// -Float_t *MHGausEvents::CreatePSDXaxis(Int_t n) -{ - - Float_t *xaxis = new Float_t[n]; - - for (Int_t i=0;i -*/ -//End_Html -// -void MHGausEvents::Draw(const Option_t *opt) -{ - - TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this,600, 900); - - TString option(opt); - option.ToLower(); - - Int_t win = 1; - - if (option.Contains("events")) - { - option.ReplaceAll("events",""); - win += 1; - } - if (option.Contains("fourier")) - { - option.ReplaceAll("fourier",""); - win += 2; - } - - pad->SetBorderMode(0); - pad->Divide(1,win); - pad->cd(1); - - if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow()) - gPad->SetLogy(); - - gPad->SetTicks(); - - fHGausHist.Draw(option); - - if (fFGausFit) - { - fFGausFit->SetLineColor(IsGausFitOK() ? kGreen : kRed); - fFGausFit->Draw("same"); - } - switch (win) - { - case 2: - pad->cd(2); - DrawEvents(); - break; - case 3: - pad->cd(2); - DrawPowerSpectrum(*pad,3); - break; - case 4: - pad->cd(2); - DrawEvents(); - pad->cd(3); - DrawPowerSpectrum(*pad,4); - break; - } -} - -void MHGausEvents::DrawEvents() -{ - - if (!fGraphEvents) - CreateGraphEvents(); - - if (!fGraphEvents) - return; - - fGraphEvents->SetBit(kCanDelete); - fGraphEvents->SetTitle("Events with time"); - fGraphEvents->Draw("AL"); - -} - - -void MHGausEvents::DrawPowerSpectrum(TVirtualPad &pad, Int_t i) -{ - - if (fPowerSpectrum) - { - if (!fGraphPowerSpectrum) - CreateGraphPowerSpectrum(); - - fGraphPowerSpectrum->Draw("AL"); - fGraphPowerSpectrum->SetBit(kCanDelete); - } - - pad.cd(i); - - if (fHPowerProbability && fHPowerProbability->GetEntries() > 0) - { - gPad->SetLogy(); - fHPowerProbability->Draw(); - if (fFExpFit) - { - fFExpFit->SetLineColor(IsExpFitOK() ? kGreen : kRed); - fFExpFit->Draw("same"); - } - } -} - - -// -------------------------------------------------------------------------- -// -// Fill fEvents with f -// If size of fEvents is 0, initializes it to 512 -// If size of fEvents is smaller than fCurrentSize, double the size -// Increase fCurrentSize by 1 -// -void MHGausEvents::FillArray(const Float_t f) -{ - if (fEvents.GetSize() == 0) - fEvents.Set(512); - - if (fCurrentSize >= fEvents.GetSize()) - fEvents.Set(fEvents.GetSize()*2); - - fEvents.AddAt(f,fCurrentSize++); -} - - -// -------------------------------------------------------------------------- -// -// Fills fHGausHist with f -// Returns kFALSE, if overflow or underflow occurred, else kTRUE -// -Bool_t MHGausEvents::FillHist(const Float_t f) -{ - return fHGausHist.Fill(f) > -1; -} - -// -------------------------------------------------------------------------- -// -// Executes: -// - FillArray() -// - FillHist() -// -Bool_t MHGausEvents::FillHistAndArray(const Float_t f) -{ - - FillArray(f); - return FillHist(f); -} - -// ------------------------------------------------------------------- -// -// Fit fGausHist with a Gaussian after stripping zeros from both ends -// and rebinned to the number of bins specified in fBinsAfterStripping -// -// 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 MHGausEvents::FitGaus(Option_t *option, const Double_t xmin, const Double_t xmax) -{ - - if (IsGausFitOK()) - return kTRUE; - - // - // First, strip the zeros from the edges which contain only zeros and rebin - // to fBinsAfterStripping bins. If fBinsAfterStripping is 0, reduce bins only by - // a factor 10. - // - // (ATTENTION: The Chisquare method is more sensitive, - // the _less_ bins, you have!) - // - StripZeros(&fHGausHist, - fBinsAfterStripping ? fBinsAfterStripping - : (fNbins > 1000 ? fNbins/10 : 0)); - - TAxis *axe = fHGausHist.GetXaxis(); - // - // Get the fitting ranges - // - Axis_t rmin = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetFirst()) : xmin; - Axis_t rmax = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetLast()) : xmax; - - // - // First guesses for the fit (should be as close to reality as possible, - // - const Stat_t entries = fHGausHist.Integral(axe->FindBin(rmin),axe->FindBin(rmax),"width"); - const Double_t mu_guess = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); - const Double_t sigma_guess = fHGausHist.GetRMS(); - const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess; - - fFGausFit = new TF1("GausFit","gaus",rmin,rmax); - - if (!fFGausFit) - { - *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit " - << "in: " << fName << endl; - return kFALSE; - } - - fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess); - fFGausFit->SetParNames("Area","#mu","#sigma"); - fFGausFit->SetParLimits(0,0.,area_guess*25.); - fFGausFit->SetParLimits(1,rmin,rmax); - fFGausFit->SetParLimits(2,0.,rmax-rmin); - fFGausFit->SetRange(rmin,rmax); - - 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; -} - - -const Double_t MHGausEvents::GetChiSquare() const -{ - return ( fFGausFit ? fFGausFit->GetChisquare() : 0.); -} - - -const Double_t MHGausEvents::GetExpChiSquare() const -{ - return ( fFExpFit ? fFExpFit->GetChisquare() : 0.); -} - - -const Int_t MHGausEvents::GetExpNdf() const -{ - return ( fFExpFit ? fFExpFit->GetNDF() : 0); -} - - -const Double_t MHGausEvents::GetExpProb() const -{ - return ( fFExpFit ? fFExpFit->GetProb() : 0.); -} - - -const Int_t MHGausEvents::GetNdf() const -{ - return ( fFGausFit ? fFGausFit->GetNDF() : 0); -} - -const Double_t MHGausEvents::GetOffset() const -{ - return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.); -} - - - -// -------------------------------------------------------------------------- -// -// If fFExpFit exists, returns fit parameter 1 (Slope of Exponential fit), -// otherwise 0. -// -const Double_t MHGausEvents::GetSlope() const -{ - return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.); -} - -// -------------------------------------------------------------------------- -// -// Default InitBins, can be overloaded. -// -// Executes: -// - fHGausHist.SetBins(fNbins,fFirst,fLast) -// -void MHGausEvents::InitBins() -{ - fHGausHist.SetBins(fNbins,fFirst,fLast); -} - -// -------------------------------------------------------------------------- -// -// Return kFALSE if number of entries is 0 -// -const Bool_t MHGausEvents::IsEmpty() const -{ - return !(fHGausHist.GetEntries()); -} - -// -------------------------------------------------------------------------- -// -// Return kTRUE if number of entries is bin content of fNbins+1 -// -const Bool_t MHGausEvents::IsOnlyOverflow() const -{ - return fHGausHist.GetEntries() == fHGausHist.GetBinContent(fNbins+1); -} - -// -------------------------------------------------------------------------- -// -// Return kTRUE if number of entries is bin content of 0 -// -const Bool_t MHGausEvents::IsOnlyUnderflow() const -{ - return fHGausHist.GetEntries() == fHGausHist.GetBinContent(0); -} - - -const Bool_t MHGausEvents::IsExcluded() const -{ - return TESTBIT(fFlags,kExcluded); -} - - -const Bool_t MHGausEvents::IsExpFitOK() const -{ - return TESTBIT(fFlags,kExpFitOK); -} - -const Bool_t MHGausEvents::IsFourierSpectrumOK() const -{ - return TESTBIT(fFlags,kFourierSpectrumOK); -} - - -const Bool_t MHGausEvents::IsGausFitOK() const -{ - return TESTBIT(fFlags,kGausFitOK); -} - - -// ----------------------------------------------------------------------------------- -// -// A default print -// -void MHGausEvents::Print(const Option_t *o) const -{ - - *fLog << all << endl; - *fLog << all << "Results of the Gauss Fit in: " << fName << endl; - *fLog << all << "Mean: " << GetMean() << endl; - *fLog << all << "Sigma: " << GetSigma() << endl; - *fLog << all << "Chisquare: " << GetChiSquare() << endl; - *fLog << all << "DoF: " << GetNdf() << endl; - *fLog << all << "Probability: " << GetProb() << endl; - *fLog << all << endl; - -} - - -// -------------------------------------------------------------------------- -// -// Default Reset(), can be overloaded. -// -// Executes: -// - Clear() -// - fHGausHist.Reset() -// - fEvents.Set(0) -// -void MHGausEvents::Reset() -{ - - Clear(); - fHGausHist.Reset(); - fEvents.Set(0); - -} - -// -------------------------------------------------------------------------- -// -// Set Excluded bit from outside -// -void MHGausEvents::SetExcluded(const Bool_t b) -{ - b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded); -} - - -// ------------------------------------------------------------------- -// -// The flag setters are to be used ONLY for Monte-Carlo!! -// -void MHGausEvents::SetExpFitOK(const Bool_t b) -{ - - b ? SETBIT(fFlags,kExpFitOK) : CLRBIT(fFlags,kExpFitOK); -} - -// ------------------------------------------------------------------- -// -// The flag setters are to be used ONLY for Monte-Carlo!! -// -void MHGausEvents::SetFourierSpectrumOK(const Bool_t b) -{ - - b ? SETBIT(fFlags,kFourierSpectrumOK) : CLRBIT(fFlags,kFourierSpectrumOK); -} - - -// ------------------------------------------------------------------- -// -// The flag setters are to be used ONLY for Monte-Carlo!! -// -void MHGausEvents::SetGausFitOK(const Bool_t b) -{ - b ? SETBIT(fFlags,kGausFitOK) : CLRBIT(fFlags,kGausFitOK); - -} - Index: unk/MagicSoft/Mars/mcalib/MHGausEvents.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHGausEvents.h (revision 4929) +++ (revision ) @@ -1,157 +1,0 @@ -#ifndef MARS_MHGausEvents -#define MARS_MHGausEvents - -#ifndef ROOT_TH1 -#include -#endif - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TVirtualPad; -class TGraph; -class TArrayF; -class TH1F; -class TH1I; -class TF1; - -class MHGausEvents : public MH -{ -private: - - const static Int_t fgNDFLimit; //! Default for fNDFLimit (now set to: 2) - const static Float_t fgProbLimit; //! Default for fProbLimit (now set to: 0.001) - const static Int_t fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20) - - Int_t fBinsAfterStripping; // Bins for the Gauss Histogram after stripping off the zeros at both ends - Int_t fCurrentSize; // Current size of the array fEvents - Byte_t fFlags; // Bit field for the fit result bits - Int_t fPowerProbabilityBins; // Bins for the projected power spectrum - - TH1I *fHPowerProbability; // Fourier transform of fEvents projected on y-axis - TArrayF *fPowerSpectrum; // Fourier transform of fEvents - - enum { kGausFitOK, - kExpFitOK, - kFourierSpectrumOK, - kExcluded }; // Bits for information about fit results - -protected: - - TArrayF fEvents; // Array which holds the entries of GausHist - TF1 *fFGausFit; // Gauss fit for fHGausHist - TF1 *fFExpFit; // Exponential fit for FHPowerProbability - Axis_t fFirst; // Lower histogram edge for fHGausHist (used by InitBins()) - TGraph *fGraphEvents; //! TGraph to display the event array (will not be cloned!!) - TGraph *fGraphPowerSpectrum; //! TGraph to display the power spectrum array (will not be cloned!!) - TH1F fHGausHist; // Histogram to hold the Gaussian distribution - Axis_t fLast; // Upper histogram edge for fHGausHist (used by InitBins()) - Double_t fMean; // Mean of the Gauss fit - Double_t fMeanErr; // Error of the mean of the Gauss fit - Int_t fNbins; // Number histogram bins for fHGausHist (used by InitBins()) - Int_t fNDFLimit; // NDF limit for judgement if fit is OK - Double_t fSigma; // Sigma of the Gauss fit - Double_t fSigmaErr; // Error of the sigma of the Gauss fit - Double_t fProb; // Probability of the Gauss fit - Float_t fProbLimit; // Probability limit for judgement if fit is OK - - virtual Float_t *CreateEventXaxis(Int_t n); // Create an x-axis for the Event TGraphs - virtual Float_t *CreatePSDXaxis(Int_t n); // Create an x-axis for the PSD TGraphs - virtual void CreateGraphEvents(); // Create the TGraph fGraphEvents of fEvents - virtual void CreateGraphPowerSpectrum(); // Create the TGraph fGraphPowerSpectrum out of fPowerSpectrum - - void DrawPowerSpectrum(TVirtualPad &pad, Int_t i); // Draw graph of fPowerSpectrum and fHPowerProbability - - // Setters - void SetBinsAfterStripping ( const Int_t nbins=0 ) { fBinsAfterStripping =nbins; } - void SetPowerProbabilityBins ( const Int_t nbins=fgPowerProbabilityBins ) { fPowerProbabilityBins=nbins; } - -public: - - MHGausEvents(const char* name=NULL, const char* title=NULL); - ~MHGausEvents(); - - virtual void Clear(Option_t *o=""); - virtual void Reset(); - virtual void InitBins(); - - // Draws - virtual void Draw(Option_t *option=""); // Default Draw - void DrawEvents(); // Draw graph of fEvents - - // Getters - const Double_t GetChiSquare() const; - const Double_t GetExpChiSquare() const; - const Int_t GetExpNdf() const; - const Double_t GetExpProb() const; - TArrayF *GetEvents() { return &fEvents; } - const TArrayF *GetEvents() const { return &fEvents; } - TF1 *GetFExpFit() { return fFExpFit; } - const TF1 *GetFExpFit() const { return fFExpFit; } - TF1 *GetFGausFit() { return fFGausFit; } - const TF1 *GetFGausFit() const { return fFGausFit; } - TGraph *GetGraphEvents() { return fGraphEvents; } - const Double_t GetFirst() const { return fFirst; } - const Double_t GetLast () const { return fLast ; } - const TGraph *GetGraphEvents() const { return fGraphEvents; } - TGraph *GetGraphPowerSpectrum() { return fGraphPowerSpectrum; } - const TGraph *GetGraphPowerSpectrum() const { return fGraphPowerSpectrum; } - TH1F *GetHGausHist() { return &fHGausHist; } - const TH1F *GetHGausHist() const { return &fHGausHist; } - TH1I *GetHPowerProbability() { return fHPowerProbability; } - const TH1I *GetHPowerProbability() const { return fHPowerProbability; } - const Double_t GetMean() const { return fMean; } - const Double_t GetMeanErr() const { return fMeanErr; } - const Int_t GetNdf() const; - const Double_t GetOffset() const; - TArrayF *GetPowerSpectrum() { return fPowerSpectrum; } - const TArrayF *GetPowerSpectrum() const { return fPowerSpectrum; } - const Double_t GetProb() const { return fProb; } - const Double_t GetSigma() const { return fSigma; } - const Double_t GetSigmaErr() const { return fSigmaErr; } - const Double_t GetSlope() const; - - const Bool_t IsExcluded() const; - const Bool_t IsExpFitOK() const; - const Bool_t IsEmpty() const; - const Bool_t IsFourierSpectrumOK() const; - const Bool_t IsGausFitOK() const; - const Bool_t IsOnlyOverflow() const; - const Bool_t IsOnlyUnderflow() const; - - // Fill - void FillArray ( const Float_t f ); // Fill only the array fEvents - Bool_t FillHist ( const Float_t f ); // Fill only the histogram HGausHist - Bool_t FillHistAndArray( const Float_t f ); // Fill bothe the array fEvents and the histogram HGausHist - - // Fits - Bool_t FitGaus( Option_t *option="RQ0", - const Double_t xmin=0., - const Double_t xmax=0.); // Fit the histogram HGausHist with a Gaussian - - // Prints - virtual void Print(const Option_t *o="") const; // Default Print - - // Setters - void SetExcluded ( const Bool_t b=kTRUE ); - void SetExpFitOK ( const Bool_t b=kTRUE ); - void SetFourierSpectrumOK( const Bool_t b=kTRUE ); - void SetGausFitOK ( const Bool_t b=kTRUE ); - void SetLast ( const Double_t d ) { fLast = d; } - void SetFirst ( const Double_t d ) { fFirst = d; } - void SetMean ( const Double_t d ) { fMean = d; } - void SetMeanErr ( const Double_t d ) { fMeanErr = d; } - void SetNbins ( const Int_t i ) { fNbins = i; } - void SetNDFLimit ( const Int_t lim=fgNDFLimit ) { fNDFLimit = lim; } - void SetProb ( const Double_t d ) { fProb = d; } - void SetProbLimit ( const Float_t lim=fgProbLimit ) { fProbLimit = lim; } - void SetSigma ( const Double_t d ) { fSigma = d; } - void SetSigmaErr ( const Double_t d ) { fSigmaErr = d; } - - void CreateFourierSpectrum(); // Create the fourier spectrum out of fEvents - - ClassDef(MHGausEvents, 1) // Base class for events with Gaussian distributed values -}; - -#endif Index: /trunk/MagicSoft/Mars/mcalib/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mcalib/Makefile (revision 4929) +++ /trunk/MagicSoft/Mars/mcalib/Makefile (revision 4930) @@ -19,5 +19,5 @@ # connect the include files defined in the config.mk file # -INCLUDES = -I. -I../mbase -I../mhbase -I../mgui -I../mgeom \ +INCLUDES = -I. -I../mbase -I../mgui -I../mgeom \ -I../manalysis -I../mraw -I../mtools -I../mmc -I../mhist \ -I../mimage -I../msignal -I../mbadpixels -I../mpedestal @@ -60,19 +60,5 @@ MCalibrationChargeBlindPix.cc \ MCalibrationChargePINDiode.cc \ - MHCalibrationChargeBlindPix.cc \ - MHCalibrationChargeBlindCam.cc \ - MHCalibrationChargePix.cc \ - MHCalibrationCam.cc \ - MHCalibrationPix.cc \ - MHCalibrationChargeCam.cc \ - MHCalibrationChargePINDiode.cc \ - MHCalibrationRelTimeCam.cc \ - MHCalibrationRelTimePix.cc \ - MHCalibrationTestCam.cc \ - MHCalibrationTestPix.cc \ - MHCalibrationTestTimeCam.cc \ - MHCalibrationTestTimePix.cc \ - MMcCalibrationCalc.cc \ - MHGausEvents.cc + MMcCalibrationCalc.cc ############################################################ Index: /trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h (revision 4930) +++ /trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h (revision 4930) @@ -0,0 +1,23 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MHCalibrationCam+; +#pragma link C++ class MHCalibrationPix+; +#pragma link C++ class MHCalibrationChargeCam+; +#pragma link C++ class MHCalibrationChargePix+; +#pragma link C++ class MHCalibrationChargeBlindCam+; +#pragma link C++ class MHCalibrationChargeBlindPix+; +#pragma link C++ class MHCalibrationChargePINDiode+; +#pragma link C++ class MHCalibrationRelTimePix+; +#pragma link C++ class MHCalibrationRelTimeCam+; +#pragma link C++ class MHCalibrationTestCam+; +#pragma link C++ class MHCalibrationTestPix+; +#pragma link C++ class MHCalibrationTestTimeCam+; +#pragma link C++ class MHCalibrationTestTimePix+; + +#pragma link C++ class MHGausEvents++; + +#endif Index: /trunk/MagicSoft/Mars/mhcalib/MHGausEvents.cc =================================================================== --- /trunk/MagicSoft/Mars/mhcalib/MHGausEvents.cc (revision 4930) +++ /trunk/MagicSoft/Mars/mhcalib/MHGausEvents.cc (revision 4930) @@ -0,0 +1,912 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHGausEvents +// +// A base class for events which are believed to follow a Gaussian distribution +// with time, e.g. calibration events, observables containing white noise, ... +// +// MHGausEvents derives from MH, thus all features of MH can be used by a class +// deriving from MHGausEvents, especially the filling functions. +// +// The central objects are: +// +// 1) The TH1F fHGausHist: +// ==================== +// +// It is created with a default name and title and resides in directory NULL. +// - Any class deriving from MHGausEvents needs to apply a binning to fHGausHist +// (e.g. by setting the variables fNbins, fFirst, fLast and calling the function +// InitBins() or by directly calling fHGausHist.SetBins(..) ) +// - The histogram is filled with the functions FillHist() or FillHistAndArray(). +// - The histogram can be fitted with the function FitGaus(). This involves stripping +// of all zeros at the lower and upper end of the histogram and re-binning to +// a new number of bins, specified in the variable fBinsAfterStripping. +// - The fit result's probability is compared to a reference probability fProbLimit +// The NDF is compared to fNDFLimit and a check is made whether results are NaNs. +// Accordingly, a flag IsGausFitOK() is set. +// - One can repeat the fit within a given amount of sigmas from the previous mean +// (specified by the variables fPickupLimit and fBlackoutLimit) with the function RepeatFit() +// - One can completely skip the fitting to set mean, sigma and its errors directly +// from the histograms with the function BypassFit() +// +// 2) The TArrayF fEvents: +// ========================== +// +// It is created with 0 entries and not expanded unless FillArray() or FillHistAndArray() +// are called. +// - A first call to FillArray() or FillHistAndArray() initializes fEvents by default +// to 512 entries. +// - Any later call to FillArray() or FillHistAndArray() fills up the array. +// Reaching the limit, the array is expanded by a factor 2. +// - The array can be fourier-transformed into the array fPowerSpectrum. +// Note that any FFT accepts only number of events which are a power of 2. +// Thus, fEvents is cut to the next power of 2 smaller than its actual number of entries. +// Be aware that you might lose information at the end of your analysis. +// - Calling the function CreateFourierSpectrum() creates the array fPowerSpectrum +// and its projection fHPowerProbability which in turn is fit to an exponential. +// - The fit result's probability is compared to a referenc probability fProbLimit +// and accordingly the flag IsExpFitOK() is set. +// - The flag IsFourierSpectrumOK() is set accordingly to IsExpFitOK(). +// Later, a closer check will be installed. +// - You can display all arrays by calls to: CreateGraphEvents() and +// CreateGraphPowerSpectrum() and successive calls to the corresponding Getters. +// +// To see an example, have a look at: Draw() +// +////////////////////////////////////////////////////////////////////////////// +#include "MHGausEvents.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "MFFT.h" +#include "MArray.h" + +#include "MH.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHGausEvents); + +using namespace std; + +const Int_t MHGausEvents::fgNDFLimit = 2; +const Float_t MHGausEvents::fgProbLimit = 0.001; +const Int_t MHGausEvents::fgPowerProbabilityBins = 30; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// Sets: +// - the default Probability Bins for fPowerProbabilityBins (fgPowerProbabilityBins) +// - the default Probability Limit for fProbLimit (fgProbLimit) +// - the default NDF Limit for fNDFLimit (fgNDFLimit) +// - the default number of bins after stripping for fBinsAfterStipping (fgBinsAfterStipping) +// - the default name of the fHGausHist ("HGausHist") +// - the default title of the fHGausHist ("Histogram of Events with Gaussian Distribution") +// - the default directory of the fHGausHist (NULL) +// - the default for fNbins (100) +// - the default for fFirst (0.) +// - the default for fLast (100.) +// +// Initializes: +// - fEvents to 0 entries +// - fHGausHist() +// - all other pointers to NULL +// - all variables to 0. +// - all flags to kFALSE +// +MHGausEvents::MHGausEvents(const char *name, const char *title) + : fHPowerProbability(NULL), + fPowerSpectrum(NULL), + fFGausFit(NULL), fFExpFit(NULL), + fFirst(0.), + fGraphEvents(NULL), fGraphPowerSpectrum(NULL), + fLast(100.), fNbins(100) +{ + + fName = name ? name : "MHGausEvents"; + fTitle = title ? title : "Events with expected Gaussian distributions"; + + Clear(); + + SetBinsAfterStripping(); + SetNDFLimit(); + SetPowerProbabilityBins(); + SetProbLimit(); + + fHGausHist.SetName("HGausHist"); + fHGausHist.SetTitle("Histogram of Events with Gaussian Distribution"); + // important, other ROOT will not draw the axes: + fHGausHist.UseCurrentStyle(); + fHGausHist.SetDirectory(NULL); + fHGausHist.GetYaxis()->CenterTitle(); +} + + + +// -------------------------------------------------------------------------- +// +// Default Destructor. +// +// Deletes (if Pointer is not NULL): +// +// - fHPowerProbability +// - fFGausFit +// - fFExpFit +// - fPowerSpectrum +// - fGraphEvents +// - fGraphPowerSpectrum +// +MHGausEvents::~MHGausEvents() +{ + + // delete histograms +// if (fHPowerProbability) +// if (gROOT->FindObject(fHPowerProbability->GetName())) +// delete fHPowerProbability; + + // delete fits + if (fFGausFit) + delete fFGausFit; + if (fFExpFit) + delete fFExpFit; + + // delete arrays + if (fPowerSpectrum) + delete fPowerSpectrum; + + // delete graphs + if (fGraphEvents) + delete fGraphEvents; + if (fGraphPowerSpectrum) + delete fGraphPowerSpectrum; +} + +// -------------------------------------------------------------------------- +// +// Default Clear(), can be overloaded. +// +// Sets: +// - all other pointers to NULL +// - all variables to 0. +// - all flags to kFALSE +// +// Deletes (if not NULL): +// - all pointers +// +void MHGausEvents::Clear(Option_t *o) +{ + + SetGausFitOK ( kFALSE ); + SetExpFitOK ( kFALSE ); + SetFourierSpectrumOK( kFALSE ); + SetExcluded ( kFALSE ); + + fMean = 0.; + fSigma = 0.; + fMeanErr = 0.; + fSigmaErr = 0.; + fProb = 0.; + + fCurrentSize = 0; + + if (fHPowerProbability) + { + if (gROOT->FindObject(fHPowerProbability->GetName())) + delete fHPowerProbability; + fHPowerProbability = NULL; + } + + // delete fits + if (fFGausFit) + { + delete fFGausFit; + fFGausFit = NULL; + } + + if (fFExpFit) + { + delete fFExpFit; + fFExpFit = NULL; + } + + // delete arrays + if (fPowerSpectrum) + { + delete fPowerSpectrum; + fPowerSpectrum = NULL; + } + + // delete graphs + if (fGraphEvents) + { + delete fGraphEvents; + fGraphEvents = NULL; + } + + if (fGraphPowerSpectrum) + { + delete fGraphPowerSpectrum; + fGraphPowerSpectrum = NULL; + } +} + +TObject *MHGausEvents::Clone(const char *name) const +{ + + MHGausEvents &gaus = *new MHGausEvents; + + gaus.fBinsAfterStripping = fBinsAfterStripping; + gaus.fCurrentSize = fCurrentSize; + gaus.fFlags = fFlags; + gaus.fPowerProbabilityBins = fPowerProbabilityBins; + /* + if (fHPowerProbability) + gaus.fHPowerProbability=(TH1I*)fHPowerProbability->Clone(); + + if (fPowerSpectrum) + gaus.fPowerSpectrum = new TArrayF(*fPowerSpectrum); + */ + + gaus.fEvents.Set(fEvents.GetSize()); + gaus.fEvents = fEvents; + /* + if (fFGausFit) + gaus.fFGausFit=(TF1*)fFGausFit->Clone(); + if (fFExpFit) + gaus.fFExpFit=(TF1*)fFExpFit->Clone(); + */ + + gaus.fFirst = fFirst; + + /* + if (fGraphEvents) + gaus.fGraphEvents=(TGraph*)fGraphEvents->Clone(); + if (fGraphPowerSpectrum) + gaus.fGraphPowerSpectrum=(TGraph*)fGraphPowerSpectrum->Clone(); + */ + + Copy(gaus.fHGausHist); + + gaus.fLast = fLast; + gaus.fMean = fMean; + gaus.fMeanErr = fMeanErr; + gaus.fNbins = fNbins; + gaus.fNDFLimit = fNDFLimit; + gaus.fSigma = fSigma; + gaus.fSigma = fSigmaErr; + gaus.fProb = fProb; + gaus.fProbLimit = fProbLimit; + + return &gaus; +} + + +// ----------------------------------------------------------------------------- +// +// Create the x-axis for the event graph +// +Float_t *MHGausEvents::CreateEventXaxis(Int_t n) +{ + + Float_t *xaxis = new Float_t[n]; + + for (Int_t i=0;iSetXTitle("P(f)"); + fHPowerProbability->SetYTitle("Counts"); + fHPowerProbability->GetYaxis()->CenterTitle(); + fHPowerProbability->SetDirectory(NULL); + fHPowerProbability->SetBit(kCanDelete); + // + // First guesses for the fit (should be as close to reality as possible, + // + const Double_t xmax = fHPowerProbability->GetXaxis()->GetXmax(); + + fFExpFit = new TF1("FExpFit","exp([0]-[1]*x)",0.,xmax); + + const Double_t slope_guess = (TMath::Log(fHPowerProbability->GetEntries())+1.)/xmax; + const Double_t offset_guess = slope_guess*xmax; + + fFExpFit->SetParameters(offset_guess, slope_guess); + fFExpFit->SetParNames("Offset","Slope"); + fFExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess); + fFExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess); + fFExpFit->SetRange(0.,xmax); + + fHPowerProbability->Fit(fFExpFit,"RQL0"); + + if (GetExpProb() > fProbLimit) + SetExpFitOK(kTRUE); + + // + // For the moment, this is the only check, later we can add more... + // + SetFourierSpectrumOK(IsExpFitOK()); + + return; +} + +// ---------------------------------------------------------------------------------- +// +// Create a graph to display the array fEvents +// If the variable fEventFrequency is set, the x-axis is transformed into real time. +// +void MHGausEvents::CreateGraphEvents() +{ + + MArray::StripZeros(fEvents); + + const Int_t n = fEvents.GetSize(); + + if (n==0) + return; + + fGraphEvents = new TGraph(n,CreateEventXaxis(n),fEvents.GetArray()); + fGraphEvents->SetTitle("Evolution of Events with time"); + fGraphEvents->GetYaxis()->SetTitle(fHGausHist.GetXaxis()->GetTitle()); + fGraphEvents->GetYaxis()->CenterTitle(); +} + +// ---------------------------------------------------------------------------------- +// +// Create a graph to display the array fPowerSpectrum +// If the variable fEventFrequency is set, the x-axis is transformed into real frequency. +// +void MHGausEvents::CreateGraphPowerSpectrum() +{ + + MArray::StripZeros(*fPowerSpectrum); + + const Int_t n = fPowerSpectrum->GetSize(); + + fGraphPowerSpectrum = new TGraph(n,CreatePSDXaxis(n),fPowerSpectrum->GetArray()); + fGraphPowerSpectrum->SetTitle("Power Spectrum Density"); + fGraphPowerSpectrum->GetYaxis()->SetTitle("P(f)"); + fGraphPowerSpectrum->GetYaxis()->CenterTitle(); + +} + + +// ----------------------------------------------------------------------------- +// +// Create the x-axis for the event graph +// +Float_t *MHGausEvents::CreatePSDXaxis(Int_t n) +{ + + Float_t *xaxis = new Float_t[n]; + + for (Int_t i=0;i +*/ +//End_Html +// +void MHGausEvents::Draw(const Option_t *opt) +{ + + TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this,600, 900); + + TString option(opt); + option.ToLower(); + + Int_t win = 1; + + if (option.Contains("events")) + { + option.ReplaceAll("events",""); + win += 1; + } + if (option.Contains("fourier")) + { + option.ReplaceAll("fourier",""); + win += 2; + } + + pad->SetBorderMode(0); + pad->Divide(1,win); + pad->cd(1); + + if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow()) + gPad->SetLogy(); + + gPad->SetTicks(); + + fHGausHist.Draw(option); + + if (fFGausFit) + { + fFGausFit->SetLineColor(IsGausFitOK() ? kGreen : kRed); + fFGausFit->Draw("same"); + } + switch (win) + { + case 2: + pad->cd(2); + DrawEvents(); + break; + case 3: + pad->cd(2); + DrawPowerSpectrum(*pad,3); + break; + case 4: + pad->cd(2); + DrawEvents(); + pad->cd(3); + DrawPowerSpectrum(*pad,4); + break; + } +} + +void MHGausEvents::DrawEvents() +{ + + if (!fGraphEvents) + CreateGraphEvents(); + + if (!fGraphEvents) + return; + + fGraphEvents->SetBit(kCanDelete); + fGraphEvents->SetTitle("Events with time"); + fGraphEvents->Draw("AL"); + +} + + +void MHGausEvents::DrawPowerSpectrum(TVirtualPad &pad, Int_t i) +{ + + if (fPowerSpectrum) + { + if (!fGraphPowerSpectrum) + CreateGraphPowerSpectrum(); + + fGraphPowerSpectrum->Draw("AL"); + fGraphPowerSpectrum->SetBit(kCanDelete); + } + + pad.cd(i); + + if (fHPowerProbability && fHPowerProbability->GetEntries() > 0) + { + gPad->SetLogy(); + fHPowerProbability->Draw(); + if (fFExpFit) + { + fFExpFit->SetLineColor(IsExpFitOK() ? kGreen : kRed); + fFExpFit->Draw("same"); + } + } +} + + +// -------------------------------------------------------------------------- +// +// Fill fEvents with f +// If size of fEvents is 0, initializes it to 512 +// If size of fEvents is smaller than fCurrentSize, double the size +// Increase fCurrentSize by 1 +// +void MHGausEvents::FillArray(const Float_t f) +{ + if (fEvents.GetSize() == 0) + fEvents.Set(512); + + if (fCurrentSize >= fEvents.GetSize()) + fEvents.Set(fEvents.GetSize()*2); + + fEvents.AddAt(f,fCurrentSize++); +} + + +// -------------------------------------------------------------------------- +// +// Fills fHGausHist with f +// Returns kFALSE, if overflow or underflow occurred, else kTRUE +// +Bool_t MHGausEvents::FillHist(const Float_t f) +{ + return fHGausHist.Fill(f) > -1; +} + +// -------------------------------------------------------------------------- +// +// Executes: +// - FillArray() +// - FillHist() +// +Bool_t MHGausEvents::FillHistAndArray(const Float_t f) +{ + + FillArray(f); + return FillHist(f); +} + +// ------------------------------------------------------------------- +// +// Fit fGausHist with a Gaussian after stripping zeros from both ends +// and rebinned to the number of bins specified in fBinsAfterStripping +// +// 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 MHGausEvents::FitGaus(Option_t *option, const Double_t xmin, const Double_t xmax) +{ + + if (IsGausFitOK()) + return kTRUE; + + // + // First, strip the zeros from the edges which contain only zeros and rebin + // to fBinsAfterStripping bins. If fBinsAfterStripping is 0, reduce bins only by + // a factor 10. + // + // (ATTENTION: The Chisquare method is more sensitive, + // the _less_ bins, you have!) + // + StripZeros(&fHGausHist, + fBinsAfterStripping ? fBinsAfterStripping + : (fNbins > 1000 ? fNbins/10 : 0)); + + TAxis *axe = fHGausHist.GetXaxis(); + // + // Get the fitting ranges + // + Axis_t rmin = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetFirst()) : xmin; + Axis_t rmax = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetLast()) : xmax; + + // + // First guesses for the fit (should be as close to reality as possible, + // + const Stat_t entries = fHGausHist.Integral(axe->FindBin(rmin),axe->FindBin(rmax),"width"); + const Double_t mu_guess = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + const Double_t sigma_guess = fHGausHist.GetRMS(); + const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess; + + fFGausFit = new TF1("GausFit","gaus",rmin,rmax); + + if (!fFGausFit) + { + *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit " + << "in: " << fName << endl; + return kFALSE; + } + + fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess); + fFGausFit->SetParNames("Area","#mu","#sigma"); + fFGausFit->SetParLimits(0,0.,area_guess*25.); + fFGausFit->SetParLimits(1,rmin,rmax); + fFGausFit->SetParLimits(2,0.,rmax-rmin); + fFGausFit->SetRange(rmin,rmax); + + 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; +} + + +const Double_t MHGausEvents::GetChiSquare() const +{ + return ( fFGausFit ? fFGausFit->GetChisquare() : 0.); +} + + +const Double_t MHGausEvents::GetExpChiSquare() const +{ + return ( fFExpFit ? fFExpFit->GetChisquare() : 0.); +} + + +const Int_t MHGausEvents::GetExpNdf() const +{ + return ( fFExpFit ? fFExpFit->GetNDF() : 0); +} + + +const Double_t MHGausEvents::GetExpProb() const +{ + return ( fFExpFit ? fFExpFit->GetProb() : 0.); +} + + +const Int_t MHGausEvents::GetNdf() const +{ + return ( fFGausFit ? fFGausFit->GetNDF() : 0); +} + +const Double_t MHGausEvents::GetOffset() const +{ + return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.); +} + + + +// -------------------------------------------------------------------------- +// +// If fFExpFit exists, returns fit parameter 1 (Slope of Exponential fit), +// otherwise 0. +// +const Double_t MHGausEvents::GetSlope() const +{ + return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.); +} + +// -------------------------------------------------------------------------- +// +// Default InitBins, can be overloaded. +// +// Executes: +// - fHGausHist.SetBins(fNbins,fFirst,fLast) +// +void MHGausEvents::InitBins() +{ + fHGausHist.SetBins(fNbins,fFirst,fLast); +} + +// -------------------------------------------------------------------------- +// +// Return kFALSE if number of entries is 0 +// +const Bool_t MHGausEvents::IsEmpty() const +{ + return !(fHGausHist.GetEntries()); +} + +// -------------------------------------------------------------------------- +// +// Return kTRUE if number of entries is bin content of fNbins+1 +// +const Bool_t MHGausEvents::IsOnlyOverflow() const +{ + return fHGausHist.GetEntries() == fHGausHist.GetBinContent(fNbins+1); +} + +// -------------------------------------------------------------------------- +// +// Return kTRUE if number of entries is bin content of 0 +// +const Bool_t MHGausEvents::IsOnlyUnderflow() const +{ + return fHGausHist.GetEntries() == fHGausHist.GetBinContent(0); +} + + +const Bool_t MHGausEvents::IsExcluded() const +{ + return TESTBIT(fFlags,kExcluded); +} + + +const Bool_t MHGausEvents::IsExpFitOK() const +{ + return TESTBIT(fFlags,kExpFitOK); +} + +const Bool_t MHGausEvents::IsFourierSpectrumOK() const +{ + return TESTBIT(fFlags,kFourierSpectrumOK); +} + + +const Bool_t MHGausEvents::IsGausFitOK() const +{ + return TESTBIT(fFlags,kGausFitOK); +} + + +// ----------------------------------------------------------------------------------- +// +// A default print +// +void MHGausEvents::Print(const Option_t *o) const +{ + + *fLog << all << endl; + *fLog << all << "Results of the Gauss Fit in: " << fName << endl; + *fLog << all << "Mean: " << GetMean() << endl; + *fLog << all << "Sigma: " << GetSigma() << endl; + *fLog << all << "Chisquare: " << GetChiSquare() << endl; + *fLog << all << "DoF: " << GetNdf() << endl; + *fLog << all << "Probability: " << GetProb() << endl; + *fLog << all << endl; + +} + + +// -------------------------------------------------------------------------- +// +// Default Reset(), can be overloaded. +// +// Executes: +// - Clear() +// - fHGausHist.Reset() +// - fEvents.Set(0) +// +void MHGausEvents::Reset() +{ + + Clear(); + fHGausHist.Reset(); + fEvents.Set(0); + +} + +// -------------------------------------------------------------------------- +// +// Set Excluded bit from outside +// +void MHGausEvents::SetExcluded(const Bool_t b) +{ + b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded); +} + + +// ------------------------------------------------------------------- +// +// The flag setters are to be used ONLY for Monte-Carlo!! +// +void MHGausEvents::SetExpFitOK(const Bool_t b) +{ + + b ? SETBIT(fFlags,kExpFitOK) : CLRBIT(fFlags,kExpFitOK); +} + +// ------------------------------------------------------------------- +// +// The flag setters are to be used ONLY for Monte-Carlo!! +// +void MHGausEvents::SetFourierSpectrumOK(const Bool_t b) +{ + + b ? SETBIT(fFlags,kFourierSpectrumOK) : CLRBIT(fFlags,kFourierSpectrumOK); +} + + +// ------------------------------------------------------------------- +// +// The flag setters are to be used ONLY for Monte-Carlo!! +// +void MHGausEvents::SetGausFitOK(const Bool_t b) +{ + b ? SETBIT(fFlags,kGausFitOK) : CLRBIT(fFlags,kGausFitOK); + +} + Index: /trunk/MagicSoft/Mars/mhcalib/MHGausEvents.h =================================================================== --- /trunk/MagicSoft/Mars/mhcalib/MHGausEvents.h (revision 4930) +++ /trunk/MagicSoft/Mars/mhcalib/MHGausEvents.h (revision 4930) @@ -0,0 +1,161 @@ +#ifndef MARS_MHGausEvents +#define MARS_MHGausEvents + +#ifndef ROOT_TH1 +#include +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TVirtualPad; +class TGraph; +class TArrayF; +class TH1F; +class TH1I; +class TF1; + +class MHGausEvents : public MH +{ +private: + + const static Int_t fgNDFLimit; //! Default for fNDFLimit (now set to: 2) + const static Float_t fgProbLimit; //! Default for fProbLimit (now set to: 0.001) + const static Int_t fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20) + + Int_t fBinsAfterStripping; // Bins for the Gauss Histogram after stripping off the zeros at both ends + Int_t fCurrentSize; // Current size of the array fEvents + Byte_t fFlags; // Bit field for the fit result bits + Int_t fPowerProbabilityBins; // Bins for the projected power spectrum + + TH1I *fHPowerProbability; // Fourier transform of fEvents projected on y-axis + TArrayF *fPowerSpectrum; // Fourier transform of fEvents + + enum { kGausFitOK, + kExpFitOK, + kFourierSpectrumOK, + kExcluded }; // Bits for information about fit results + +protected: + + TArrayF fEvents; // Array which holds the entries of GausHist + TF1 *fFGausFit; // Gauss fit for fHGausHist + TF1 *fFExpFit; // Exponential fit for FHPowerProbability + Axis_t fFirst; // Lower histogram edge for fHGausHist (used by InitBins()) + TGraph *fGraphEvents; //! TGraph to display the event array (will not be cloned!!) + TGraph *fGraphPowerSpectrum; //! TGraph to display the power spectrum array (will not be cloned!!) + TH1F fHGausHist; // Histogram to hold the Gaussian distribution + Axis_t fLast; // Upper histogram edge for fHGausHist (used by InitBins()) + Double_t fMean; // Mean of the Gauss fit + Double_t fMeanErr; // Error of the mean of the Gauss fit + Int_t fNbins; // Number histogram bins for fHGausHist (used by InitBins()) + Int_t fNDFLimit; // NDF limit for judgement if fit is OK + Double_t fSigma; // Sigma of the Gauss fit + Double_t fSigmaErr; // Error of the sigma of the Gauss fit + Double_t fProb; // Probability of the Gauss fit + Float_t fProbLimit; // Probability limit for judgement if fit is OK + + virtual Float_t *CreateEventXaxis(Int_t n); // Create an x-axis for the Event TGraphs + virtual Float_t *CreatePSDXaxis(Int_t n); // Create an x-axis for the PSD TGraphs + virtual void CreateGraphEvents(); // Create the TGraph fGraphEvents of fEvents + virtual void CreateGraphPowerSpectrum(); // Create the TGraph fGraphPowerSpectrum out of fPowerSpectrum + + void DrawPowerSpectrum(TVirtualPad &pad, Int_t i); // Draw graph of fPowerSpectrum and fHPowerProbability + + // Setters + void SetBinsAfterStripping ( const Int_t nbins=0 ) { fBinsAfterStripping =nbins; } + void SetPowerProbabilityBins ( const Int_t nbins=fgPowerProbabilityBins ) { fPowerProbabilityBins=nbins; } + +public: + + MHGausEvents(const char* name=NULL, const char* title=NULL); + ~MHGausEvents(); + + TObject *Clone(const char* name="") const; + + void Clear(Option_t *o=""); + void Reset(); + + // Draws + void Draw(Option_t *option=""); // Default Draw + void DrawEvents(); // Draw graph of fEvents + + // Inits + virtual void InitBins(); + + // Getters + const Double_t GetChiSquare() const; + const Double_t GetExpChiSquare() const; + const Int_t GetExpNdf() const; + const Double_t GetExpProb() const; + TArrayF *GetEvents() { return &fEvents; } + const TArrayF *GetEvents() const { return &fEvents; } + TF1 *GetFExpFit() { return fFExpFit; } + const TF1 *GetFExpFit() const { return fFExpFit; } + TF1 *GetFGausFit() { return fFGausFit; } + const TF1 *GetFGausFit() const { return fFGausFit; } + TGraph *GetGraphEvents() { return fGraphEvents; } + const Double_t GetFirst() const { return fFirst; } + const Double_t GetLast () const { return fLast ; } + const TGraph *GetGraphEvents() const { return fGraphEvents; } + TGraph *GetGraphPowerSpectrum() { return fGraphPowerSpectrum; } + const TGraph *GetGraphPowerSpectrum() const { return fGraphPowerSpectrum; } + TH1F *GetHGausHist() { return &fHGausHist; } + const TH1F *GetHGausHist() const { return &fHGausHist; } + TH1I *GetHPowerProbability() { return fHPowerProbability; } + const TH1I *GetHPowerProbability() const { return fHPowerProbability; } + const Double_t GetMean() const { return fMean; } + const Double_t GetMeanErr() const { return fMeanErr; } + const Int_t GetNdf() const; + const Double_t GetOffset() const; + TArrayF *GetPowerSpectrum() { return fPowerSpectrum; } + const TArrayF *GetPowerSpectrum() const { return fPowerSpectrum; } + const Double_t GetProb() const { return fProb; } + const Double_t GetSigma() const { return fSigma; } + const Double_t GetSigmaErr() const { return fSigmaErr; } + const Double_t GetSlope() const; + + const Bool_t IsExcluded() const; + const Bool_t IsExpFitOK() const; + const Bool_t IsEmpty() const; + const Bool_t IsFourierSpectrumOK() const; + const Bool_t IsGausFitOK() const; + const Bool_t IsOnlyOverflow() const; + const Bool_t IsOnlyUnderflow() const; + + // Fill + void FillArray ( const Float_t f ); // Fill only the array fEvents + Bool_t FillHist ( const Float_t f ); // Fill only the histogram HGausHist + Bool_t FillHistAndArray( const Float_t f ); // Fill bothe the array fEvents and the histogram HGausHist + + // Fits + Bool_t FitGaus( Option_t *option="RQ0", + const Double_t xmin=0., + const Double_t xmax=0.); // Fit the histogram HGausHist with a Gaussian + + // Prints + virtual void Print(const Option_t *o="") const; // Default Print + + // Setters + void SetExcluded ( const Bool_t b=kTRUE ); + void SetExpFitOK ( const Bool_t b=kTRUE ); + void SetFourierSpectrumOK( const Bool_t b=kTRUE ); + void SetGausFitOK ( const Bool_t b=kTRUE ); + void SetLast ( const Double_t d ) { fLast = d; } + void SetFirst ( const Double_t d ) { fFirst = d; } + void SetMean ( const Double_t d ) { fMean = d; } + void SetMeanErr ( const Double_t d ) { fMeanErr = d; } + void SetNbins ( const Int_t i ) { fNbins = i; } + void SetNDFLimit ( const Int_t lim=fgNDFLimit ) { fNDFLimit = lim; } + void SetProb ( const Double_t d ) { fProb = d; } + void SetProbLimit ( const Float_t lim=fgProbLimit ) { fProbLimit = lim; } + void SetSigma ( const Double_t d ) { fSigma = d; } + void SetSigmaErr ( const Double_t d ) { fSigmaErr = d; } + + void CreateFourierSpectrum(); // Create the fourier spectrum out of fEvents + + ClassDef(MHGausEvents, 1) // Base class for events with Gaussian distributed values +}; + +#endif Index: /trunk/MagicSoft/Mars/mhcalib/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mhcalib/Makefile (revision 4930) +++ /trunk/MagicSoft/Mars/mhcalib/Makefile (revision 4930) @@ -0,0 +1,52 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +#------------------------------------------------------------------------------ + +# +# Handling name of the Root Dictionary Files +# +CINT = HCalib + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mhbase \ + -I../manalysis -I../mraw -I../mtools -I../mmc -I../mhist \ + -I../mimage -I../msignal -I../mbadpixels -I../mpedestal + +# mhbase: MBinning MH +# manalysis: MExtractedSignal, MCerPhotEvt (move to mcalib?) +# mraw: MRawRunHeader, MRawEvtHeader, MRawEvtPixelIter (3xMCalibrationCalc) +# mmc: MMcFadcHeader, MMcEvt +# mimage MHillas + +SRCFILES = MHCalibrationChargeBlindPix.cc \ + MHCalibrationChargeBlindCam.cc \ + MHCalibrationChargePix.cc \ + MHCalibrationCam.cc \ + MHCalibrationPix.cc \ + MHCalibrationChargeCam.cc \ + MHCalibrationChargePINDiode.cc \ + MHCalibrationRelTimeCam.cc \ + MHCalibrationRelTimePix.cc \ + MHCalibrationTestCam.cc \ + MHCalibrationTestPix.cc \ + MHCalibrationTestTimeCam.cc \ + MHCalibrationTestTimePix.cc \ + MHGausEvents.cc + +############################################################ + +all: $(OBJS) + +include ../Makefile.rules + +mrproper: clean rmbak