Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationPulseTimeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationPulseTimeCam.cc (revision 6593) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationPulseTimeCam.cc (revision 6593) @@ -0,0 +1,654 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationPulseTimeCam +// +// Fills the extracted signals of MExtractedSignalCam into the MHCalibrationPix-classes +// MHCalibrationPulseTimeHiGainPix and MHCalibrationPulseTimeLoGainPix for every: +// +// - Pixel, stored in the TOrdCollection's MHCalibrationCam::fHiGainArray and +// MHCalibrationCam::fLoGainArray +// +// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), +// stored in the TOrdCollection's MHCalibrationCam::fAverageHiGainAreas and +// MHCalibrationCam::fAverageLoGainAreas +// +// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), +// stored in the TOrdCollection'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 PulseTime 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: +// - MCalibrationPulseTimePix::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 MCalibrationPulseTimeCam: +// +// - MCalibrationPix::SetHiGainSaturation() +// - MCalibrationPix::SetHiGainMean() +// - MCalibrationPix::SetHiGainMeanErr() +// - MCalibrationPix::SetHiGainSigma() +// - MCalibrationPix::SetHiGainSigmaErr() +// - MCalibrationPix::SetHiGainProb() +// - MCalibrationPix::SetHiGainNumPickup() +// +// 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 "MHCalibrationPulseTimeCam.h" +#include "MHCalibrationCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MHCalibrationPix.h" + +#include "MCalibrationIntensityCam.h" +#include "MCalibrationChargeCam.h" +#include "MCalibrationChargePix.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MBadPixelsIntensityCam.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 +#include + +ClassImp(MHCalibrationPulseTimeCam); + +using namespace std; + +const Int_t MHCalibrationPulseTimeCam::fgHiGainNbins = 40; +const Axis_t MHCalibrationPulseTimeCam::fgHiGainFirst = -0.5; +const Axis_t MHCalibrationPulseTimeCam::fgHiGainLast = 19.5; +const Float_t MHCalibrationPulseTimeCam::fgProbLimit = 0.0001; +const Float_t MHCalibrationPulseTimeCam::fgLowerSignalLimit = 100.; +const TString MHCalibrationPulseTimeCam::gsHistName = "PulseTime"; +const TString MHCalibrationPulseTimeCam::gsHistTitle = "Extracted Times"; +const TString MHCalibrationPulseTimeCam::gsHistXTitle = "Time [FADC slices]"; +const TString MHCalibrationPulseTimeCam::gsHistYTitle = "Nr. events"; +const TString MHCalibrationPulseTimeCam::fgReferenceFile = "mjobs/calibrationref.rc"; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// Sets: +// - all pointers to NULL +// +// - fNbins to fgHiGainNbins +// - fFirst to fgHiGainFirst +// - fLast to fgHiGainLast +// +// - fHistName to gsHistName +// - fHistTitle to gsHistTitle +// - fHistXTitle to gsHistXTitle +// - fHistYTitle to gsHistYTitle +// +MHCalibrationPulseTimeCam::MHCalibrationPulseTimeCam(const char *name, const char *title) + : fRawEvt(NULL) +{ + + fName = name ? name : "MHCalibrationPulseTimeCam"; + fTitle = title ? title : "Class to fill the extracted pulse times for cosmics "; + + SetNbins(fgHiGainNbins); + SetFirst(fgHiGainFirst); + SetLast (fgHiGainLast ); + + SetProbLimit(fgProbLimit); + + SetHistName (gsHistName .Data()); + SetHistTitle (gsHistTitle .Data()); + SetHistXTitle(gsHistXTitle.Data()); + SetHistYTitle(gsHistYTitle.Data()); + + SetReferenceFile(); + SetLoGain(kFALSE); + + SetLowerSignalLimit(); + + fInnerRefTime = 5.; + fOuterRefTime = 5.; +} +// -------------------------------------------------------------------------- +// +// Creates new MHCalibrationPulseTimeCam only with the averaged areas: +// the rest has to be retrieved directly, e.g. via: +// MHCalibrationPulseTimeCam *cam = MParList::FindObject("MHCalibrationPulseTimeCam"); +// - cam->GetAverageSector(5).DrawClone(); +// - (*cam)[100].DrawClone() +// +TObject *MHCalibrationPulseTimeCam::Clone(const char *) const +{ + + MHCalibrationPulseTimeCam *cam = new MHCalibrationPulseTimeCam(); + + // + // Copy the data members + // + cam->fColor = fColor; + cam->fRunNumbers = fRunNumbers; + cam->fPulserFrequency = fPulserFrequency; + cam->fFlags = fFlags; + cam->fNbins = fNbins; + cam->fFirst = fFirst; + cam->fLast = fLast; + cam->fReferenceFile = fReferenceFile; + + if (!IsAverageing()) + return cam; + + const Int_t navhi = fAverageHiGainAreas->GetSize(); + + for (int i=0; ifAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i); + + return cam; +} + +// -------------------------------------------------------------------------- +// +// Gets the pointers to: +// - MRawEvtData +// +Bool_t MHCalibrationPulseTimeCam::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 +// - MCalibrationPulseTimeCam or MCalibrationIntensityPulseTimeCam +// - MBadPixelsCam +// +// Initializes the number of used FADC slices from MExtractedSignalCam +// into MCalibrationPulseTimeCam and test for changes in that variable +// +// Calls: +// - InitHiGainArrays() +// +// Sets: +// - fSumhiarea to nareas +// - fSumloarea to nareas +// - fSumhisector to nsectors +// - fSumlosector to nsectors +// +Bool_t MHCalibrationPulseTimeCam::ReInitHists(MParList *pList) +{ + + MExtractedSignalCam *signal = + (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam")); + if (!signal) + { + *fLog << err << "MExtractedSignalCam not found... abort." << endl; + return kFALSE; + } + + if (!InitCams(pList,"PulseTime")) + return kFALSE; + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nsectors = fGeom->GetNumSectors(); + const Int_t nareas = fGeom->GetNumAreas(); + + InitHiGainArrays(npixels,nareas,nsectors); + + fSumhiarea .Set(nareas); + fSumloarea .Set(nareas); + fSumhisector.Set(nsectors); + fSumlosector.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 TOrdCollection's (including the averaged ones), the following steps are performed: +// +// 1) Fill PulseTimes histograms (MHGausEvents::FillHistAndArray()) with: +// - MExtractedSignalPix::GetExtractedSignalHiGain(); +// - MExtractedSignalPix::GetExtractedSignalLoGain(); +// +Bool_t MHCalibrationPulseTimeCam::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(); + + fSumhiarea .Reset(); + fSumhisector.Reset(); + + MRawEvtPixelIter pixel(fRawEvt); + while (pixel.Next()) + { + + const Int_t i = pixel.GetPixelId(); + + MHCalibrationPix &histhi = (*this)[i]; + + if (histhi.IsExcluded()) + continue; + + const MExtractedSignalPix &pix = (*signal)[i]; + + const Int_t sathi = (Int_t)pix.GetNumHiGainSaturated(); + + if (sathi) + continue; + + const Float_t sighi = pix.GetExtractedSignalHiGain(); + + if (sighi < fLowerSignalLimit) + continue; + + const Float_t sighi = (Float_t)pixel.GetIdxMaxHiGainSample(); + + histhi.FillHist(sighi); + + const Int_t aidx = (*fGeom)[i].GetAidx(); + const Int_t sector = (*fGeom)[i].GetSector(); + + fSumhiarea[aidx] += sighi; + fSumhisector[sector] += sighi; + + } + + for (UInt_t j=0; jGetCam() : fCam; + MBadPixelsCam *badcam = fIntensBad ? fIntensBad->GetCam() : fBadPixels; + + // + // Perform the fitting for the High Gain (done in MHCalibrationCam) + // + FitHiGainArrays(*pulsecam, *badcam, + MBadPixelsPix::kHiGainNotFitted, + MBadPixelsPix::kHiGainOscillating); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calls MHCalibrationPix::DrawClone() for pixel idx +// +void MHCalibrationPulseTimeCam::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 MHCalibrationPulseTimeCam::Draw(const Option_t *opt) +{ + + const Int_t nareas = fAverageHiGainAreas->GetSize(); + 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); + + MHCalibrationPix &hipix = GetAverageHiGainArea(i); + DrawDataCheckPixel(hipix,i ? fOuterRefPulseTime : fInnerRefPulseTime); + } +} + +// ----------------------------------------------------------------------------- +// +// Draw the average pixel for the datacheck: +// +// Displays the averaged areas, both High Gain and Low Gain +// +// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options +// +void MHCalibrationPulseTimeCam::DrawDataCheckPixel(MHCalibrationPix &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(); + + TH1F *null = new TH1F("Null",hist->GetTitle(),100, + pix.GetFirst() > 0. ? pix.GetFirst() : 0., + pix.GetLast() > pix.GetFirst() + ? ( pix.GetLast() > 450. + ? ( fColor == MCalibrationCam::kBLUE ? 800. : 450. ) + : pix.GetLast() ) + : pix.GetFirst()*2.); + + null->SetMaximum(1.1*hist->GetMaximum()); + 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(); + + TF1 *fit = pix.GetFGausFit(); + + if (fit) + { + switch ( fColor ) + { + case MCalibrationCam::kGREEN: + fit->SetLineColor(kGreen); + break; + case MCalibrationCam::kBLUE: + fit->SetLineColor(kBlue); + break; + case MCalibrationCam::kUV: + fit->SetLineColor(106); + break; + case MCalibrationCam::kCT1: + fit->SetLineColor(006); + break; + default: + fit->SetLineColor(kRed); + } + fit->Draw("same"); + } + + DisplayRefLines(null,refline); + + newpad->cd(2); + gPad->SetTicks(); + + TH1F *null2 = new TH1F("Null2",hist->GetTitle(),100,0.,pix.GetEvents()->GetSize()/pix.GetEventFrequency()); + + null2->SetMinimum(pix.GetMean()-10.*pix.GetSigma()); + null2->SetMaximum(pix.GetMean()+10.*pix.GetSigma()); + null2->SetDirectory(NULL); + null2->SetBit(kCanDelete); + null2->SetStats(kFALSE); + // + // set the labels bigger + // + TAxis *xaxe2 = null2->GetXaxis(); + TAxis *yaxe2 = null2->GetYaxis(); + xaxe2->CenterTitle(); + yaxe2->CenterTitle(); + xaxe2->SetTitleSize(0.07); + yaxe2->SetTitleSize(0.07); + xaxe2->SetTitleOffset(0.7); + yaxe2->SetTitleOffset(0.55); + xaxe2->SetLabelSize(0.06); + yaxe2->SetLabelSize(0.06); + + pix.CreateGraphEvents(); + TGraph *gr = pix.GetGraphEvents(); + + xaxe2->SetTitle(gr->GetXaxis()->GetTitle()); + yaxe2->SetTitle(gr->GetYaxis()->GetTitle()); + + null2->Draw(); + + pix.DrawEvents("same"); + + // newpad->cd(3); + // pix.DrawPowerSpectrum(*newpad,4); + + return; + +} + +void MHCalibrationPulseTimeCam::DisplayRefLines(const TH1F *hist, const Float_t refline) const +{ + + TGraph *uv10 = new TGraph(2); + uv10->SetPoint(0,refline,0.); + uv10->SetPoint(1,refline,hist->GetMaximum()); + uv10->SetBit(kCanDelete); + uv10->SetLineColor(106); + uv10->SetLineStyle(2); + uv10->SetLineWidth(3); + uv10->Draw("L"); + + TLegend *leg = new TLegend(0.8,0.55,0.99,0.99); + leg->SetBit(kCanDelete); + leg->AddEntry(uv10,"10 Leds UV","l"); + + leg->Draw(); +} + +Int_t MHCalibrationPulseTimeCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + + Bool_t rc = kFALSE; + + if (MHCalibrationCam::ReadEnv(env,prefix,print)) + rc = kTRUE; + + if (IsEnvDefined(env, prefix, "LowerSignalLimit", print)) + { + SetLowerSignalLimit(GetEnvValue(env,prefix,"LowerSignalLimit",fLowerSignalLimit)); + rc = kTRUE; + } + + if (IsEnvDefined(env, prefix, "ReferenceFile", print)) + { + SetReferenceFile(GetEnvValue(env,prefix,"ReferenceFile",fReferenceFile.Data())); + rc = kTRUE; + } + + TEnv refenv(fReferenceFile); + + fInnerRefPulseTime = refenv.GetValue("InnerRefPulseTime",fInnerRefPulseTime); + fOuterRefPulseTime = refenv.GetValue("OuterRefPulseTime",fOuterRefPulseTime); + + return rc; +} Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationPulseTimeCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationPulseTimeCam.h (revision 6592) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationPulseTimeCam.h (revision 6593) @@ -15,9 +15,10 @@ class MRawEvtData; +class MHCalibrationChargePix; class MHCalibrationPulseTimeCam : public MHCalibrationCam { private: - static const Float_t fgLowerSignalLimit //! Default for fLowerSignalLimit + static const Float_t fgLowerSignalLimit; //! Default for fLowerSignalLimit static const Int_t fgHiGainNbins; //! Nr. bins of HiGain Histograms (now set to: 550 )