Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc (revision 3247) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc (revision 3247) @@ -0,0 +1,614 @@ +/* ======================================================================== *\ +! +! * +! * 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 09/2003 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MCalibrationChargeCalc +// +// Task to calculate the calibration conversion factors from the FADC +// time slices. The integrated time slices have to be delivered by an +// MExtractedSignalCam. The pedestals by an MPedestalCam. +// +// The output container MCalibrationCam holds one entry of type MCalibrationPix +// for every pixel. It is filled in the following way: +// +// ProProcess: Search for MPedestalCam, MExtractedSignalCam and MExtractedSignalBlindPixel +// Initialize MCalibrationCam +// Initialize pulser light wavelength +// +// ReInit: MCalibrationCam::InitSize(NumPixels) is called which allocates +// memory in a TClonesArray of type MCalibrationPix +// Initialize number of used FADC slices +// Optionally exclude pixels from calibration +// +// Process: Every MCalibrationPix holds a histogram class, +// MHCalibrationPixel which itself hold histograms of type: +// HCharge(npix) (distribution of summed FADC time slice +// entries) +// HTime(npix) (distribution of position of maximum) +// HChargevsN(npix) (distribution of charges vs. event number. +// +// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian +// All histograms HTime(npix) are fitted to a Gaussian +// The histogram HBlindPixelCharge (blind pixel) is fitted to +// a single PhE fit +// +// The histograms of the PIN Diode are fitted to Gaussians +// +// Fits can be excluded via the commands: +// MalibrationCam::SkipBlindPixelFits() (skip all blind +// pixel fits) +// +// Hi-Gain vs. Lo-Gain Calibration (very memory-intensive) +// can be skipped with the command: +// MalibrationCam::SkipHiLoGainCalibration() +// +// Input Containers: +// MRawEvtData +// MPedestalCam +// MExtractedSignalCam +// MExtractedSignalBlindPixel +// +// Output Containers: +// MCalibrationCam +// +////////////////////////////////////////////////////////////////////////////// +#include "MCalibrationChargeCalc.h" + +// FXIME: Usage of fstream is a preliminary workaround! +#include + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MRawRunHeader.h" +#include "MRawEvtPixelIter.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" + +#include "MCalibrationChargeCam.h" +#include "MCalibrationChargePINDiode.h" +#include "MCalibrationPix.h" + +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" +#include "MExtractedSignalBlindPixel.h" + +#include "MCalibrationBlindPix.h" + +ClassImp(MCalibrationChargeCalc); + +using namespace std; + +const UInt_t MCalibrationChargeCalc::fgBlindPixelIdx = 559; +const UInt_t MCalibrationChargeCalc::fgPINDiodeIdx = 9999; +const UInt_t MCalibrationChargeCalc::fgBlindPixelSinglePheCut = 400; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title) + : fPedestals(NULL), fCam(NULL), + fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL), + fSignals(NULL), fBlindPixel(NULL), fPINDiode(NULL) +{ + + fName = name ? name : "MCalibrationChargeCalc"; + fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam "; + + AddToBranchList("MRawEvtData.fHiGainPixId"); + AddToBranchList("MRawEvtData.fLoGainPixId"); + AddToBranchList("MRawEvtData.fHiGainFadcSamples"); + AddToBranchList("MRawEvtData.fLoGainFadcSamples"); + + Clear(); + SetBlindPixelIdx(); + SetPINDiodeIdx(); +} + +void MCalibrationChargeCalc::Clear(const Option_t *o) +{ + + SETBIT(fFlags, kUseBlindPixelFit); + SETBIT(fFlags, kUseQualityChecks); + SETBIT(fFlags, kHiLoGainCalibration); + + CLRBIT(fFlags, kHiGainOverFlow); + CLRBIT(fFlags, kLoGainOverFlow); + + fNumBlindPixelSinglePhe = 0; + fNumBlindPixelPedestal = 0; + + fNumHiGainSamples = 0; + fNumLoGainSamples = 0; + fConversionHiLo = 0; + fNumExcludedPixels = 0; + +} + + +// -------------------------------------------------------------------------- +// +// The PreProcess searches for the following input containers: +// - MRawEvtData +// - MPedestalCam +// +// The following output containers are also searched and created if +// they were not found: +// +// - MHCalibrationBlindPixel +// - MCalibrationCam +// +// The following output containers are only searched, but not created +// +// - MTime +// +Int_t MCalibrationChargeCalc::PreProcess(MParList *pList) +{ + + fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); + if (!fRawEvt) + { + *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl; + return kFALSE; + } + + const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!runheader) + *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; + else + if (runheader->GetRunType() == kRTMonteCarlo) + { + return kTRUE; + } + + fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam"); + if (!fCam) + { + *fLog << err << dbginf << "MCalibrationChargeCam could not be created ... aborting." << endl; + return kFALSE; + } + + fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode"); + if (!fPINDiode) + { + *fLog << err << dbginf << "MCalibrationChargePINDiode could not be created ... aborting." << endl; + return kFALSE; + } + + fCam->SetPINDiode(fPINDiode); + + fEvtTime = (MTime*)pList->FindObject("MTime"); + + fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPedestals) + { + *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl; + return kFALSE; + } + + + fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); + if (!fSignals) + { + *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl; + return kFALSE; + } + + fBlindPixel = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel"); + if (!fBlindPixel) + { + *fLog << err << dbginf << "Cannot find MExtractedSignalBlindPixel ... aborting" << endl; + return kFALSE; + } + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// The ReInit searches for the following input containers: +// - MRawRunHeader +// +Bool_t MCalibrationChargeCalc::ReInit(MParList *pList ) +{ + + fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + + MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!cam) + { + *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl; + return kFALSE; + } + + fCam->SetGeomCam(cam); + + fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices(); + fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices(); + fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples); + + UInt_t npixels = cam->GetNumPixels(); + + for (UInt_t i=0; iGetFirstUsedSliceHiGain(), + fSignals->GetLastUsedSliceHiGain()); + pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(), + fSignals->GetLastUsedSliceLoGain()); + + if (!TESTBIT(fFlags,kUseQualityChecks)) + pix.SetExcludeQualityCheck(); + + // Exclude the blind pixel and the PIN Diode from normal pixel calibration: + if (i == fBlindPixelIdx) + pix.SetExcluded(); + + if (i == fPINDiodeIdx) + pix.SetExcluded(); + + } + + // + // Look for file to exclude pixels from analysis + // + if (!fExcludedPixelsFile.IsNull()) + { + + fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data()); + + // + // Initialize reading the file + // + ifstream in(fExcludedPixelsFile.Data(),ios::in); + + if (in) + { + *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl; + // + // Read the file and count the number of entries + // + UInt_t pixel = 0; + + while (++fNumExcludedPixels) + { + + in >> pixel; + + if (!in.good()) + break; + // + // Check for out of range + // + if (pixel > npixels) + { + *fLog << warn << "WARNING: To be excluded pixel: " << pixel + << " is out of range " << endl; + continue; + } + // + // Exclude pixel + // + MCalibrationPix &pix = (*fCam)[pixel]; + pix.SetExcluded(); + + *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl; + + } + + if (--fNumExcludedPixels == 0) + *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data() + << "'" << " is empty " << endl; + else + fCam->SetNumPixelsExcluded(fNumExcludedPixels); + + } + else + *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl; + } + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Calculate the integral of the FADC time slices and store them as a new +// pixel in the MCerPhotEvt container. +// +Int_t MCalibrationChargeCalc::Process() +{ + + + MRawEvtPixelIter pixel(fRawEvt); + + // + // Create a loop to fill the calibration histograms + // Search for: a signal in MExtractedSignalCam + // Search for: a signal in MExtractedSignalBlindPixel + // Fill histograms with: + // charge + // charge vs. event nr. + // + // + // Initialize pointers to blind pixel and individual pixels + // + // FIXME: filling the bind pixel histograms in this class is preliminary + // and will be replaced soon to fill them with MFillH + // + MCalibrationBlindPix &blindpixel = *(fCam->GetBlindPixel()); + MExtractedSignalBlindPixel &blindsig = (*fBlindPixel); + + const UInt_t signal = blindsig.GetExtractedSignal(); + + if (!blindpixel.FillCharge(signal)) + *fLog << warn << + "Overflow or Underflow occurred filling Blind Pixel sum = " << signal << endl; + + blindpixel.FillGraphs(signal,0); + + TH1I *hist; + + if (signal > fBlindPixelSinglePheCut) + { + hist = (blindpixel.GetHist())->GetHSinglePheFADCSlices(); + fNumBlindPixelSinglePhe++; + } + else + { + hist = (blindpixel.GetHist())->GetHPedestalFADCSlices(); + fNumBlindPixelPedestal++; + } + + pixel.Jump(fBlindPixelIdx); + + const Byte_t *ptr = pixel.GetHiGainSamples(); + + for (Int_t i=1;i<16;i++) + hist->Fill(i,*ptr++); + + ptr = pixel.GetLoGainSamples(); + for (Int_t i=16;i<31;i++) + hist->Fill(i,*ptr++); + + pixel.Reset(); + + while (pixel.Next()) + { + + const UInt_t pixid = pixel.GetPixelId(); + + MCalibrationPix &pix = (*fCam)[pixid]; + MExtractedSignalPix &sig = (*fSignals) [pixid]; + + const Float_t sumhi = sig.GetExtractedSignalHiGain(); + const Float_t sumlo = sig.GetExtractedSignalLoGain(); + + Float_t abstime = 0.; + + if (sig.IsLoGainUsed()) + abstime = (Float_t)pixel.GetIdxMaxLoGainSample(); + else + abstime = (Float_t)pixel.GetIdxMaxHiGainSample(); + + if (pix.IsExcluded()) + continue; + + pix.FillGraphs(sumhi,sumlo); + + if (sig.IsLoGainUsed()) + { + + if (!pix.FillChargeLoGain(sumlo)) + *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid + << " signal = " << sumlo << endl; + + if (!pix.FillAbsTimeLoGain(abstime)) + *fLog << warn << "Could not fill Lo Gain Abs. Time of pixel: " + << pixid << " time = " << abstime << endl; + } + else + { + if (!pix.FillChargeHiGain(sumhi)) + *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid + << " signal = " << sumhi << endl; + + if (!pix.FillAbsTimeHiGain(abstime)) + *fLog << warn << "Could not fill Hi Gain Abs. Time of pixel: " + << pixid << " time = " << abstime << endl; + } + + } /* while (pixel.Next()) */ + + return kTRUE; +} + +Int_t MCalibrationChargeCalc::PostProcess() +{ + *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl; + + // + // Cut edges to make fits and viewing of the hists easier + // + fCam->CutEdges(); + + // + // Fit the blind pixel + // + if (TESTBIT(fFlags,kUseBlindPixelFit)) + { + // + // Get pointer to blind pixel + // + MCalibrationBlindPix &blindpixel = *(fCam->GetBlindPixel()); + + *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl; + + // + // retrieve the histogram containers + // + MHCalibrationBlindPixel *hist = blindpixel.GetHist(); + + // + // retrieve mean and sigma of the blind pixel pedestal, + // so that we can use it for the fit + // + // FIXME: This part has to be migrated into the future MHCalibrationChargeBlindPix class + // + const UInt_t nentries = fPedestals->GetTotalEntries(); + const UInt_t nslices = 12; + const Float_t sqrslice = TMath::Sqrt((Float_t)nslices); + + MPedestalPix &pedpix = (*fPedestals)[fBlindPixelIdx]; + if (&pedpix) + { + const Float_t pedestal = pedpix.GetPedestal()*nslices; + const Float_t pederr = pedpix.GetPedestalRms()*nslices/nentries; + const Float_t pedsigma = pedpix.GetPedestalRms()*sqrslice; + const Float_t pedsigmaerr = pederr/2.; + + hist->SetMeanPedestal(pedestal); + hist->SetMeanPedestalErr(pederr); + hist->SetSigmaPedestal(pedsigma); + hist->SetSigmaPedestalErr(pedsigmaerr); + } + + if (!blindpixel.FitCharge()) + { + *fLog << warn << "Could not fit the blind pixel! " << endl; + *fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl; + fCam->SetBlindPixelMethodValid(kFALSE); + } + else + fCam->SetBlindPixelMethodValid(kTRUE); + + if (blindpixel.CheckOscillations()) + fCam->SetBlindPixelMethodValid(kFALSE); + + TH1I *sphehist = hist->GetHSinglePheFADCSlices(); + TH1I *pedhist = hist->GetHPedestalFADCSlices(); + + if (fNumBlindPixelSinglePhe > 1) + sphehist->Scale(1./fNumBlindPixelSinglePhe); + if (fNumBlindPixelPedestal > 1) + pedhist->Scale(1./fNumBlindPixelPedestal); + + blindpixel.DrawClone(); + } + else + *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl; + + *fLog << inf << "total: " << GetNumExecutions() << " sphe: " << fNumBlindPixelSinglePhe << " ped: " << fNumBlindPixelPedestal << endl; + + + *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl; + + // + // loop over the pedestal events and check if we have calibration + // + for (Int_t pixid=0; pixidGetSize(); pixid++) + { + + MCalibrationPix &pix = (*fCam)[pixid]; + + // + // Check if the pixel has been excluded from the fits + // + if (pix.IsExcluded()) + continue; + + // + // get the pedestals + // + const Float_t ped = (*fPedestals)[pixid].GetPedestal(); + const Float_t prms = (*fPedestals)[pixid].GetPedestalRms(); + + // + // set them in the calibration camera + // + pix.SetPedestal(ped,prms,(Float_t)fNumHiGainSamples,(Float_t)fNumLoGainSamples); + + // + // perform the Gauss fits to the charges + // + pix.FitCharge(); + + // + // check also for oscillations + // + pix.CheckOscillations(); + + // + // calculate the F-Factor method + // + pix.CalcFFactorMethod(); + } + + if (TESTBIT(fFlags,kUseBlindPixelFit) && fCam->IsBlindPixelMethodValid()) + { + if (!fCam->CalcFluxInsidePlexiglass()) + { + *fLog << warn << "Could not calculate the number of photons from the blind pixel " << endl; + *fLog << "You can try to calibrate using the MCalibrationChargeCalc::SkipBlindPixelFit()" << endl; + fCam->SetBlindPixelMethodValid(kFALSE); + } + } + else + *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl; + + if (!fPINDiode->CalcFluxOutsidePlexiglass()) + { + *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl; + fCam->SetPINDiodeMethodValid(kFALSE); + } + else + { + fCam->SetPINDiodeMethodValid(kTRUE); + } + + fCam->SetReadyToSave(); + + return kTRUE; +} + + Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.h (revision 3247) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.h (revision 3247) @@ -0,0 +1,106 @@ +#ifndef MARS_MCalibrationChargeCalc +#define MARS_MCalibrationChargeCalc + +///////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationChargeCalc // +// // +// Integrates the time slices of the all pixels of a calibration event // +// and substract the pedestal value // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#include "TString.h" + +class MRawEvtData; +class MRawRunHeader; +class MPedestalCam; +class MCalibrationChargePINDiode; +class MCalibrationChargeCam; +class MGeomCam; +class MExtractedSignalCam; +class MExtractedSignalBlindPixel; +class MTime; +class MCalibrationChargeCalc : public MTask +{ + +private: + + static const UInt_t fgBlindPixelIdx; // ID of the blind pixel + static const UInt_t fgPINDiodeIdx; // ID of the PIN Diode + static const UInt_t fgBlindPixelSinglePheCut; // FADC sum from which on an event is considered as a S.ph. one. + + MPedestalCam *fPedestals; //! Pedestals of all pixels in the camera + MCalibrationChargeCam *fCam; // Calibration events of all pixels in the camera + MRawEvtData *fRawEvt; //! raw event data (time slices) + MRawRunHeader *fRunHeader; //! RunHeader information + + MTime *fEvtTime; //! Time of the event + + MExtractedSignalCam *fSignals; // Extracted signal of all pixels in the camera + MExtractedSignalBlindPixel *fBlindPixel; // Extracted signal of the blind pixel + MCalibrationChargePINDiode *fPINDiode; // Calibration events of all pixels in the camera + + + UInt_t fBlindPixelIdx; + UInt_t fPINDiodeIdx; + + Byte_t fNumHiGainSamples; + Byte_t fNumLoGainSamples; + Float_t fSqrtHiGainSamples; + + UInt_t fBlindPixelSinglePheCut; + + Int_t fNumBlindPixelSinglePhe; + Int_t fNumBlindPixelPedestal; + + Float_t fConversionHiLo; + Int_t fFlags; // Flag for the fits used + + TString fExcludedPixelsFile; + UInt_t fNumExcludedPixels; + + enum { kUseBlindPixelFit, + kUseQualityChecks, + kHiLoGainCalibration, + kHiGainOverFlow, + kLoGainOverFlow }; + + Bool_t ReInit(MParList *pList); + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + + MCalibrationChargeCalc(const char *name=NULL, const char *title=NULL); + + void Clear(const Option_t *o=""); + + void SkipBlindPixelFit(Bool_t b=kTRUE) + {b ? CLRBIT(fFlags, kUseBlindPixelFit) : SETBIT(fFlags, kUseBlindPixelFit);} + void SkipQualityChecks(Bool_t b=kTRUE) + {b ? CLRBIT(fFlags, kUseQualityChecks) : SETBIT(fFlags, kUseQualityChecks);} + void SkipHiLoGainCalibration(Bool_t b=kTRUE) + {b ? CLRBIT(fFlags, kHiLoGainCalibration) : SETBIT(fFlags, kHiLoGainCalibration);} + + + // Setters + void SetConversionHiLo ( const Float_t conv ) { fConversionHiLo = conv; } + void SetBlindPixelSinglePheCut ( const Int_t cut=fgBlindPixelSinglePheCut) + { fBlindPixelSinglePheCut = cut; } + + void SetPINDiodeIdx( const UInt_t idx=fgPINDiodeIdx ) { fPINDiodeIdx = idx; } + void SetBlindPixelIdx( const UInt_t idx=fgBlindPixelIdx ) { fBlindPixelIdx = idx; } + + // Exclude pixels from configuration file + void ExcludePixelsFromAsciiFile(const char *file) { fExcludedPixelsFile = file; } + + ClassDef(MCalibrationChargeCalc, 1) // Task to fill the Calibration Containers from raw data +}; + +#endif Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCam.cc (revision 3247) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCam.cc (revision 3247) @@ -0,0 +1,1100 @@ +/* ======================================================================== *\ +! +! * +! * 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 11/2003 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MCalibrationChargeCam +// +// Hold the whole Calibration results of the camera: +// +// 1) MCalibrationChargeCam initializes a TClonesArray whose elements are +// pointers to MCalibrationPix Containers +// 2) It initializes a pointer to an MCalibrationBlindPix container +// 3) It initializes a pointer to an MCalibrationPINDiode container +// +// +// The calculated values (types of GetPixelContent) are: +// +// Fitted values: +// ============== +// +// 0: Fitted Charge +// 1: Error of fitted Charge +// 2: Sigma of fitted Charge +// 3: Error of Sigma of fitted Charge +// +// Useful variables derived from the fit results: +// ============================================= +// +// 4: Returned probability of Gauss fit to Charge distribution +// 5: Reduced Sigma of fitted Charge --> sqrt(sigma_Q^2 - PedRMS^2) +// 6: Error Reduced Sigma of fitted Charge +// 7: Reduced Sigma per Charge +// 8: Error of Reduced Sigma per Charge +// +// Results of the different calibration methods: +// ============================================= +// +// 9: Number of Photo-electrons obtained with the F-Factor method +// 10: Error on Number of Photo-electrons obtained with the F-Factor method +// 11: Mean conversion factor obtained with the F-Factor method +// 12: Error on the mean conversion factor obtained with the F-Factor method +// 13: Overall F-Factor of the readout obtained with the F-Factor method +// 14: Error on Overall F-Factor of the readout obtained with the F-Factor method +// 15: Number of Photons inside Plexiglass obtained with the Blind Pixel method +// 16: Error on Number of Photons inside Plexiglass obtained with the Blind Pixel method +// 17: Mean conversion factor obtained with the Blind Pixel method +// 18: Error on the mean conversion factor obtained with the Blind Pixel method +// 19: Overall F-Factor of the readout obtained with the Blind Pixel method +// 20: Error on Overall F-Factor of the readout obtained with the Blind Pixel method +// 21: Number of Photons outside Plexiglass obtained with the PIN Diode method +// 22: Error on Number of Photons outside Plexiglass obtained with the PIN Diode method +// 23: Mean conversion factor obtained with the PIN Diode method +// 24: Error on the mean conversion factor obtained with the PIN Diode method +// 25: Overall F-Factor of the readout obtained with the PIN Diode method +// 26: Error on Overall F-Factor of the readout obtained with the PIN Diode method +// +// Localized defects: +// ================== +// +// 27: Excluded Pixels +// 28: Pixels where the fit did not succeed --> results obtained only from the histograms +// 29: Pixels with succeeded fit, but apparently wrong results +// 30: Pixels with un-expected behavior in the fourier spectrum (e.g. oscillations) +// +// Other classifications of pixels: +// ================================ +// +// 31: Pixels with saturated Hi-Gain +// +// Classification of validity of the calibrations: +// =============================================== +// +// 32: Pixels with valid calibration by the F-Factor-Method +// 33: Pixels with valid calibration by the Blind Pixel-Method +// 34: Pixels with valid calibration by the PIN Diode-Method +// +// Used Pedestals: +// =============== +// +// 35: Mean Pedestal over the entire range of signal extraction +// 36: Error on the Mean Pedestal over the entire range of signal extraction +// 37: Pedestal RMS over the entire range of signal extraction +// 38: Error on the Pedestal RMS over the entire range of signal extraction +// +// Calculated absolute arrival times (very low precision!): +// ======================================================== +// +// 39: Absolute Arrival time of the signal +// 40: Error on the Absolute Arrival time of the signal +// 41: RMS of the Absolute Arrival time of the signal +// 42: Error on the RMS of the Absolute Arrival time of the signal +// +///////////////////////////////////////////////////////////////////////////// +#include "MCalibrationChargeCam.h" + +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MCalibrationPix.h" +#include "MCalibrationBlindPix.h" +#include "MCalibrationChargePINDiode.h" + +#include "MHCalibrationPixel.h" + +ClassImp(MCalibrationChargeCam); + +using namespace std; + +const Int_t MCalibrationChargeCam::gkBlindPixelId = 559; +const Int_t MCalibrationChargeCam::gkPINDiodeId = 9999; +const Float_t MCalibrationChargeCam::gkBlindPixelArea = 100; +// Average QE of Blind Pixel (three colours) +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelQEGreen = 0.154; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelQEBlue = 0.226; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelQEUV = 0.247; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelQECT1 = 0.247; +// Average QE Error of Blind Pixel (three colours) +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelQEGreenErr = 0.015; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelQEBlueErr = 0.02; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelQEUVErr = 0.02; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelQECT1Err = 0.02; +// Attenuation factor Blind Pixel (three colours) +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelAttGreen = 1.97; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelAttBlue = 1.96; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelAttUV = 1.95; +const Float_t MCalibrationChargeCam::gkCalibrationBlindPixelAttCT1 = 1.95; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry +// Later, a call to MCalibrationChargeCam::InitSize(Int_t size) has to be performed +// +// Creates an MCalibrationBlindPix container +// +MCalibrationChargeCam::MCalibrationChargeCam(const char *name, const char *title) + : fBlindPixel(NULL), + fPINDiode(NULL), + fOffsets(NULL), + fSlopes(NULL), + fOffvsSlope(NULL) +{ + fName = name ? name : "MCalibrationChargeCam"; + fTitle = title ? title : "Storage container for the Calibration Information in the camera"; + + fPixels = new TClonesArray("MCalibrationPix",1); + fBlindPixel = new MCalibrationBlindPix(); + + Clear(); +} + +// -------------------------------------------------------------------------- +// +// Delete the TClonesArray of MCalibrationPix containers +// Delete the MCalibrationPINDiode and the MCalibrationBlindPix +// +// Delete the histograms if they exist +// +MCalibrationChargeCam::~MCalibrationChargeCam() +{ + + // + // delete fPixels should delete all Objects stored inside + // + delete fPixels; + delete fBlindPixel; + + if (fOffsets) + delete fOffsets; + if (fSlopes) + delete fSlopes; + if (fOffvsSlope) + delete fOffvsSlope; + +} + +// ------------------------------------------------------------------- +// +// This function simply allocates memory via the ROOT command: +// (TObject**) TStorage::ReAlloc(fCont, newSize * sizeof(TObject*), +// fSize * sizeof(TObject*)); +// newSize corresponds to size in our case +// fSize is the old size (in most cases: 1) +// +void MCalibrationChargeCam::InitSize(const UInt_t i) +{ + + // + // check if we have already initialized to size + // + if (CheckBounds(i)) + return; + + fPixels->ExpandCreate(i); + +} + +// -------------------------------------------------------------------------- +// +// This function returns the current size of the TClonesArray +// independently if the MCalibrationPix is filled with values or not. +// +// It is the size of the array fPixels. +// +Int_t MCalibrationChargeCam::GetSize() const +{ + return fPixels->GetEntriesFast(); +} + +// -------------------------------------------------------------------------- +// +// Check if position i is inside the current bounds of the TClonesArray +// +Bool_t MCalibrationChargeCam::CheckBounds(Int_t i) const +{ + return i < GetSize(); +} + + +// -------------------------------------------------------------------------- +// +// Get i-th pixel (pixel number) +// +MCalibrationPix &MCalibrationChargeCam::operator[](UInt_t i) +{ + return *static_cast(fPixels->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th pixel (pixel number) +// +const MCalibrationPix &MCalibrationChargeCam::operator[](UInt_t i) const +{ + return *static_cast(fPixels->UncheckedAt(i)); +} + + +// -------------------------------------- +// +void MCalibrationChargeCam::Clear(Option_t *o) +{ + + fPixels->ForEach(TObject, Clear)(); + fBlindPixel->Clear(); + + fMeanFluxInsidePlexiglass = -1.; + fMeanFluxErrInsidePlexiglass = -1.; + + fNumExcludedPixels = 0; + + CLRBIT(fFlags,kBlindPixelMethodValid); + CLRBIT(fFlags,kPINDiodeMethodValid); + CLRBIT(fFlags,kFluxInsidePlexiglassAvailable); + + return; +} + +void MCalibrationChargeCam::SetBlindPixelMethodValid(const Bool_t b) +{ + b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid); +} + +void MCalibrationChargeCam::SetPINDiodeMethodValid(const Bool_t b) +{ + b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid); +} + +Bool_t MCalibrationChargeCam::IsBlindPixelMethodValid() const +{ + return TESTBIT(fFlags,kBlindPixelMethodValid); +} + +Bool_t MCalibrationChargeCam::IsPINDiodeMethodValid() const +{ + return TESTBIT(fFlags,kPINDiodeMethodValid); +} + + +Bool_t MCalibrationChargeCam::IsFluxInsidePlexiglassAvailable() const +{ + return TESTBIT(fFlags,kFluxInsidePlexiglassAvailable); +} + + +// -------------------------------------------------------------------------- +// +// Print first the well fitted pixels +// and then the ones which are not FitValid +// +void MCalibrationChargeCam::Print(Option_t *o) const +{ + + *fLog << all << GetDescriptor() << ":" << endl; + int id = 0; + + *fLog << all << "Succesfully calibrated pixels:" << endl; + *fLog << all << endl; + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + + if (pix->IsChargeValid() && !pix->IsExcluded() && !pix->IsOscillating() && pix->IsFitted()) + { + + *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " + << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " + << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() + << " Nr Phe's: " << pix->GetPheFFactorMethod() << endl; + id++; + } + } + + *fLog << all << id << " succesful pixels :-))" << endl; + id = 0; + + *fLog << all << endl; + *fLog << all << "Pixels with errors:" << endl; + *fLog << all << endl; + + TIter Next2(fPixels); + while ((pix=(MCalibrationPix*)Next2())) + { + + if (!pix->IsChargeValid() && !pix->IsExcluded() && !pix->IsFitted()) + { + + *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " + << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " + << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() << endl; + id++; + } + } + *fLog << all << id << " pixels with errors :-((" << endl; + + *fLog << all << endl; + *fLog << all << "Pixels with oscillations:" << endl; + *fLog << all << endl; + + id = 0; + + TIter Next3(fPixels); + while ((pix=(MCalibrationPix*)Next3())) + { + + if (pix->IsOscillating() && !pix->IsExcluded()) + { + + *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " + << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " + << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() << endl; + id++; + } + } + *fLog << all << id << " Oscillating pixels :-((" << endl; + + + *fLog << all << endl; + *fLog << all << "Excluded pixels:" << endl; + *fLog << all << endl; + + TIter Next4(fPixels); + while ((pix=(MCalibrationPix*)Next4())) + if (pix->IsExcluded()) + *fLog << all << pix->GetPixId() << endl; + + *fLog << all << fNumExcludedPixels << " excluded pixels " << endl; +} + +// -------------------------------------------------------------------------- +// +// Return true if pixel is inside bounds of the TClonesArray fPixels +// +Bool_t MCalibrationChargeCam::IsPixelUsed(Int_t idx) const +{ + if (!CheckBounds(idx)) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Return true if pixel has already been fitted once (independent of the result) +// +Bool_t MCalibrationChargeCam::IsPixelFitted(Int_t idx) const +{ + + if (!CheckBounds(idx)) + return kFALSE; + + return (*this)[idx].IsFitted(); +} + +// -------------------------------------------------------------------------- +// +// Sets the user ranges of all histograms such that +// empty bins at the edges are not used. Additionally, it rebins the +// histograms such that in total, 50 bins are used. +// +void MCalibrationChargeCam::CutEdges() +{ + + fBlindPixel->GetHist()->CutAllEdges(); + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + pix->GetHist()->CutAllEdges(); + } + + return; +} + +// -------------------------------------------------------------------------- +// +// The types are as follows: +// +// Fitted values: +// ============== +// +// 0: Fitted Charge +// 1: Error of fitted Charge +// 2: Sigma of fitted Charge +// 3: Error of Sigma of fitted Charge +// +// Useful variables derived from the fit results: +// ============================================= +// +// 4: Returned probability of Gauss fit to Charge distribution +// 5: Reduced Sigma of fitted Charge --> sqrt(sigma_Q^2 - PedRMS^2) +// 6: Error Reduced Sigma of fitted Charge +// 7: Reduced Sigma per Charge +// 8: Error of Reduced Sigma per Charge +// +// Results of the different calibration methods: +// ============================================= +// +// 9: Number of Photo-electrons obtained with the F-Factor method +// 10: Error on Number of Photo-electrons obtained with the F-Factor method +// 11: Mean conversion factor obtained with the F-Factor method +// 12: Error on the mean conversion factor obtained with the F-Factor method +// 13: Overall F-Factor of the readout obtained with the F-Factor method +// 14: Error on Overall F-Factor of the readout obtained with the F-Factor method +// 15: Number of Photons inside Plexiglass obtained with the Blind Pixel method +// 16: Error on Number of Photons inside Plexiglass obtained with the Blind Pixel method +// 17: Mean conversion factor obtained with the Blind Pixel method +// 18: Error on the mean conversion factor obtained with the Blind Pixel method +// 19: Overall F-Factor of the readout obtained with the Blind Pixel method +// 20: Error on Overall F-Factor of the readout obtained with the Blind Pixel method +// 21: Number of Photons outside Plexiglass obtained with the PIN Diode method +// 22: Error on Number of Photons outside Plexiglass obtained with the PIN Diode method +// 23: Mean conversion factor obtained with the PIN Diode method +// 24: Error on the mean conversion factor obtained with the PIN Diode method +// 25: Overall F-Factor of the readout obtained with the PIN Diode method +// 26: Error on Overall F-Factor of the readout obtained with the PIN Diode method +// +// Localized defects: +// ================== +// +// 27: Excluded Pixels +// 28: Pixels where the fit did not succeed --> results obtained only from the histograms +// 29: Pixels with succeeded fit, but apparently wrong results +// 30: Pixels with un-expected behavior in the fourier spectrum (e.g. oscillations) +// +// Other classifications of pixels: +// ================================ +// +// 31: Pixels with saturated Hi-Gain +// +// Classification of validity of the calibrations: +// =============================================== +// +// 32: Pixels with valid calibration by the F-Factor-Method +// 33: Pixels with valid calibration by the Blind Pixel-Method +// 34: Pixels with valid calibration by the PIN Diode-Method +// +// Used Pedestals: +// =============== +// +// 35: Mean Pedestal over the entire range of signal extraction +// 36: Error on the Mean Pedestal over the entire range of signal extraction +// 37: Pedestal RMS over the entire range of signal extraction +// 38: Error on the Pedestal RMS over the entire range of signal extraction +// +// Calculated absolute arrival times (very low precision!): +// ======================================================== +// +// 39: Absolute Arrival time of the signal +// 40: Error on the Absolute Arrival time of the signal +// 41: RMS of the Absolute Arrival time of the signal +// 42: Error on the RMS of the Absolute Arrival time of the signal +// +Bool_t MCalibrationChargeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + + if (idx > GetSize()) + return kFALSE; + + Float_t area = cam[idx].GetA(); + + if (area == 0) + return kFALSE; + + switch (type) + { + case 0: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetCharge(); + break; + case 1: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetChargeErr(); + break; + case 2: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetSigmaCharge(); + break; + case 3: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetSigmaChargeErr(); + break; + case 4: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetChargeProb(); + break; + case 5: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetRSigmaCharge(); + break; + case 6: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetRSigmaChargeErr(); + break; + case 7: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge(); + break; + case 8: + if ((*this)[idx].IsExcluded()) + return kFALSE; + // relative error RsigmaCharge square + val = (*this)[idx].GetRSigmaChargeErr()* (*this)[idx].GetRSigmaChargeErr() + / ((*this)[idx].GetRSigmaCharge() * (*this)[idx].GetRSigmaCharge() ); + // relative error Charge square + val += (*this)[idx].GetChargeErr() * (*this)[idx].GetChargeErr() + / ((*this)[idx].GetCharge() * (*this)[idx].GetCharge() ); + // calculate relative error out of squares + val = TMath::Sqrt(val) ; + // multiply with value to get absolute error + val *= (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge(); + break; + case 9: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetPheFFactorMethod(); + break; + case 10: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetPheFFactorMethodErr(); + break; + case 11: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetMeanConversionFFactorMethod(); + break; + case 12: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetConversionFFactorMethodErr(); + break; + case 13: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorFFactorMethod(); + break; + case 14: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorErrFFactorMethod(); + break; + case 15: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = GetMeanFluxInsidePlexiglass()*area; + break; + case 16: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = GetMeanFluxInsidePlexiglass()*area; + break; + case 17: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetMeanConversionBlindPixelMethod(); + break; + case 18: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetConversionBlindPixelMethodErr(); + break; + case 19: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorBlindPixelMethod(); + break; + case 20: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorErrBlindPixelMethod(); + break; + case 21: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = fPINDiode->GetMeanFluxOutsidePlexiglass()*area; + break; + case 22: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = fPINDiode->GetMeanFluxOutsidePlexiglass()*area; + break; + case 23: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetMeanConversionPINDiodeMethod(); + break; + case 24: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetConversionPINDiodeMethodErr(); + break; + case 25: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorBlindPixelMethod(); + break; + case 26: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorErrBlindPixelMethod(); + break; + case 27: + if ((*this)[idx].IsExcluded()) + val = 1.; + else + return kFALSE; + break; + case 28: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (!(*this)[idx].IsFitted()) + val = 1; + else + return kFALSE; + break; + case 29: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (!(*this)[idx].IsFitted()) + return kFALSE; + if (!(*this)[idx].IsChargeValid()) + val = 1; + else + return kFALSE; + break; + case 30: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsOscillating()) + val = 1; + else + return kFALSE; + break; + case 31: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsHiGainSaturation()) + val = 1; + else + return kFALSE; + break; + case 32: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsFFactorMethodValid()) + val = 1; + else + return kFALSE; + break; + case 33: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsBlindPixelMethodValid()) + val = 1; + else + return kFALSE; + break; + case 34: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsPINDiodeMethodValid()) + val = 1; + else + return kFALSE; + break; + case 35: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetPed(); + break; + case 36: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = 1.; + // val = (*this)[idx].GetPedError(); + break; + case 37: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetPedRms(); + break; + case 38: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = 1.; + // val = (*this)[idx].GetPedRmsError(); + break; + case 39: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetAbsTimeMean(); + break; + case 40: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetAbsTimeMeanErr(); + break; + case 41: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetAbsTimeRms(); + break; + case 42: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetAbsTimeMeanErr()/TMath::Sqrt(2.); + break; + default: + return kFALSE; + } + return val!=-1.; +} + +// -------------------------------------------------------------------------- +// +// What MHCamera needs in order to draw an individual pixel in the camera +// +void MCalibrationChargeCam::DrawPixelContent(Int_t idx) const +{ + (*this)[idx].DrawClone(); +} + + +// -------------------------------------------------------------------------- +// +// +// +Bool_t MCalibrationChargeCam::CalcFluxInsidePlexiglass() +{ + + if (!fBlindPixel->IsFitOK()) + return kFALSE; + + const Float_t mean = fBlindPixel->GetLambda(); + const Float_t merr = fBlindPixel->GetErrLambda(); + + // + // Start calculation of number of photons + // + // The blind pixel has exactly 100 mm^2 area (with negligible error), + // + fMeanFluxInsidePlexiglass = mean*gkBlindPixelArea; + + // Start calculation of number of photons relative Variance (!!) + fMeanFluxErrInsidePlexiglass = merr*merr/mean/mean; + + switch (fColor) + { + case kEGreen: + fMeanFluxInsidePlexiglass /= gkCalibrationBlindPixelQEGreen; + fMeanFluxErrInsidePlexiglass += gkCalibrationBlindPixelQEGreenErr*gkCalibrationBlindPixelQEGreenErr + / gkCalibrationBlindPixelQEGreen / gkCalibrationBlindPixelQEGreen; + + fMeanFluxInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttGreen); // correct for absorption + // attenuation has negligible error + break; + case kEBlue: + fMeanFluxInsidePlexiglass /= gkCalibrationBlindPixelQEBlue; + fMeanFluxErrInsidePlexiglass += gkCalibrationBlindPixelQEBlueErr*gkCalibrationBlindPixelQEBlueErr + / gkCalibrationBlindPixelQEBlue / gkCalibrationBlindPixelQEBlue; + + fMeanFluxInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttBlue); // correct for absorption + // attenuation has negligible error + break; + case kEUV: + fMeanFluxInsidePlexiglass /= gkCalibrationBlindPixelQEUV; + fMeanFluxErrInsidePlexiglass += gkCalibrationBlindPixelQEUVErr*gkCalibrationBlindPixelQEUVErr + / gkCalibrationBlindPixelQEUV / gkCalibrationBlindPixelQEUV; + + fMeanFluxInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttUV); // correct for absorption + // attenuation has negligible error + break; + case kECT1: + default: + fMeanFluxInsidePlexiglass /= gkCalibrationBlindPixelQECT1; + fMeanFluxErrInsidePlexiglass += gkCalibrationBlindPixelQECT1Err*gkCalibrationBlindPixelQECT1Err + / gkCalibrationBlindPixelQECT1 / gkCalibrationBlindPixelQECT1; + + fMeanFluxInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttCT1); // correct for absorption + // attenuation has negligible error + break; + } + + *fLog << inf << endl; + *fLog << inf << " Photon flux [ph/mm^2] inside Plexiglass: " + << fMeanFluxInsidePlexiglass << endl; + + if (fMeanFluxInsidePlexiglass > 0.) + SETBIT(fFlags,kFluxInsidePlexiglassAvailable); + else + { + CLRBIT(fFlags,kFluxInsidePlexiglassAvailable); + return kFALSE; + } + + if (fMeanFluxErrInsidePlexiglass < 0.) + { + *fLog << warn << " Relative Variance on Photon flux inside Plexiglass: " + << fMeanFluxErrInsidePlexiglass << endl; + CLRBIT(fFlags,kFluxInsidePlexiglassAvailable); + return kFALSE; + } + + // Finish calculation of errors -> convert from relative variance to absolute error + fMeanFluxErrInsidePlexiglass = TMath::Sqrt(fMeanFluxErrInsidePlexiglass); + fMeanFluxErrInsidePlexiglass *= fMeanFluxInsidePlexiglass; + + *fLog << inf << " Error on photon flux [ph/mm^2] inside Plexiglass: " + << fMeanFluxErrInsidePlexiglass << endl; + *fLog << inf << endl; + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + + if(pix->IsChargeValid()) + { + + const Int_t idx = pix->GetPixId(); + + const Float_t charge = pix->GetCharge(); + const Float_t area = (*fGeomCam)[idx].GetA(); + const Float_t chargeerr = pix->GetChargeErr(); + + const Float_t nphot = fMeanFluxInsidePlexiglass*area; + const Float_t nphoterr = fMeanFluxErrInsidePlexiglass*area; + const Float_t conversion = nphot/charge; + Float_t conversionerr; + + conversionerr = nphoterr/charge + * nphoterr/charge ; + conversionerr += chargeerr/charge + * chargeerr/charge + * conversion*conversion; + conversionerr = TMath::Sqrt(conversionerr); + + const Float_t conversionsigma = 0.; + + pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma); + + if (conversionerr/conversion < 0.1) + pix->SetBlindPixelMethodValid(); + } + } + return kTRUE; +} + + +void MCalibrationChargeCam::ApplyPINDiodeCalibration() +{ + + Float_t flux = fPINDiode->GetMeanFluxOutsidePlexiglass(); + Float_t fluxerr = fPINDiode->GetMeanFluxErrOutsidePlexiglass(); + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + + if (pix->IsChargeValid()) + { + + const Int_t idx = pix->GetPixId(); + + const Float_t charge = pix->GetCharge(); + const Float_t area = (*fGeomCam)[idx].GetA(); + const Float_t chargeerr = pix->GetChargeErr(); + + const Float_t nphot = flux * area; + const Float_t nphoterr = fluxerr * area; + const Float_t conversion = nphot/charge; + + Float_t conversionerr; + + conversionerr = nphoterr/charge + * nphoterr/charge ; + conversionerr += chargeerr/charge + * chargeerr/charge + * conversion*conversion; + if (conversionerr > 0.) + conversionerr = TMath::Sqrt(conversionerr); + + const Float_t conversionsigma = 0.; + + pix->SetConversionPINDiodeMethod(conversion, conversionerr, conversionsigma); + + if (conversionerr/conversion < 0.1) + pix->SetPINDiodeMethodValid(); + + } + } +} + + + +Bool_t MCalibrationChargeCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) +{ + + if (ipx < 0 || !IsPixelFitted(ipx)) + return kFALSE; + + if (!IsFluxInsidePlexiglassAvailable()) + if (!CalcFluxInsidePlexiglass()) + return kFALSE; + + mean = (*this)[ipx].GetMeanConversionBlindPixelMethod(); + err = (*this)[ipx].GetConversionBlindPixelMethodErr(); + sigma = (*this)[ipx].GetSigmaConversionBlindPixelMethod(); + + return kTRUE; +} + + +Bool_t MCalibrationChargeCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) +{ + + if (ipx < 0 || !IsPixelFitted(ipx)) + return kFALSE; + + Float_t conv = (*this)[ipx].GetMeanConversionFFactorMethod(); + + if (conv < 0.) + return kFALSE; + + mean = conv; + err = (*this)[ipx].GetConversionFFactorMethodErr(); + sigma = (*this)[ipx].GetSigmaConversionFFactorMethod(); + + return kTRUE; +} + + +//----------------------------------------------------------------------------------- +// +// Calculates the conversion factor between the integral of FADCs slices +// (as defined in the signal extractor MExtractSignal.cc) +// and the number of photons reaching the plexiglass for one Inner Pixel +// +// FIXME: The PINDiode is still not working and so is the code +// +Bool_t MCalibrationChargeCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) +{ + + if (ipx < 0 || !IsPixelFitted(ipx)) + return kFALSE; + + mean = (*this)[ipx].GetMeanConversionPINDiodeMethod(); + err = (*this)[ipx].GetConversionPINDiodeMethodErr(); + sigma = (*this)[ipx].GetSigmaConversionPINDiodeMethod(); + + return kFALSE; + +} + +//----------------------------------------------------------------------------------- +// +// Calculates the best combination of the three used methods possible +// between the integral of FADCs slices +// (as defined in the signal extractor MExtractSignal.cc) +// and the number of photons reaching one Inner Pixel. +// The procedure is not yet defined. +// +// FIXME: The PINDiode is still not working and so is the code +// +Bool_t MCalibrationChargeCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) +{ + + if (ipx < 0 || !IsPixelFitted(ipx)) + return kFALSE; + + return kFALSE; + +} + +/* +void MCalibrationChargeCam::DrawHiLoFits() +{ + + if (!fOffsets) + fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-600.,400.); + if (!fSlopes) + fSlopes = new TH1D("mm","Slopes of the HiGain LoGain Fit",100,-2.,2.); + if (!fOffvsSlope) + fOffvsSlope = new TH2D("aa","Slopes vs Offsets of the HiGain LoGain Fit",100,-600.,400.,100,-2.,2.); + + TIter Next(fPixels); + MCalibrationPix *pix; + MHCalibrationPixel *hist; + while ((pix=(MCalibrationPix*)Next())) + { + hist = pix->GetHist(); + hist->FitHiGainvsLoGain(); + fOffsets->Fill(hist->GetOffset(),1.); + fSlopes->Fill(hist->GetSlope(),1.); + fOffvsSlope->Fill(hist->GetOffset(),hist->GetSlope(),1.); + } + + TCanvas *c1 = new TCanvas(); + + c1->Divide(1,3); + c1->cd(1); + fOffsets->Draw(); + gPad->Modified(); + gPad->Update(); + + c1->cd(2); + fSlopes->Draw(); + gPad->Modified(); + gPad->Update(); + + c1->cd(3); + fOffvsSlope->Draw("col1"); + gPad->Modified(); + gPad->Update(); +} + +*/ Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCam.h (revision 3247) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCam.h (revision 3247) @@ -0,0 +1,139 @@ +#ifndef MARS_MCalibrationChargeCam +#define MARS_MCalibrationChargeCam + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +#ifndef MARS_MCamEvent +#include "MCamEvent.h" +#endif + +class TH1D; +class TH2D; +class TClonesArray; + +class MCalibrationPix; +class MCalibrationBlindPix; +class MCalibrationChargePINDiode; +class MCalibrationChargeCam : public MParContainer, public MCamEvent +{ +private: + + static const Int_t gkBlindPixelId; + static const Int_t gkPINDiodeId; + static const Float_t gkBlindPixelArea; // The Blind Pixel area in mm^2 + + static const Float_t gkCalibrationBlindPixelQEGreen; + static const Float_t gkCalibrationBlindPixelQEBlue ; + static const Float_t gkCalibrationBlindPixelQEUV ; + static const Float_t gkCalibrationBlindPixelQECT1 ; + + static const Float_t gkCalibrationBlindPixelQEGreenErr; + static const Float_t gkCalibrationBlindPixelQEBlueErr ; + static const Float_t gkCalibrationBlindPixelQEUVErr ; + static const Float_t gkCalibrationBlindPixelQECT1Err ; + + static const Float_t gkCalibrationBlindPixelAttGreen; + static const Float_t gkCalibrationBlindPixelAttBlue ; + static const Float_t gkCalibrationBlindPixelAttUV ; + static const Float_t gkCalibrationBlindPixelAttCT1 ; + + Int_t fNumPixels; + TClonesArray *fPixels; //-> Array of MCalibrationPix with fit results + + MCalibrationBlindPix *fBlindPixel; //-> Pointer to the Blind Pixel with fit results + const MCalibrationChargePINDiode *fPINDiode; //! Pointer to the PIN Diode with fit results + + const MGeomCam *fGeomCam; //! Need geom cam to know which pixel in inner or outer + + TH1D* fOffsets; //! + TH1D* fSlopes; //! + + TH2D* fOffvsSlope; //! + + Float_t fMeanFluxInsidePlexiglass; // The mean number of photons in an INNER PIXEL inside the plexiglass + Float_t fMeanFluxErrInsidePlexiglass; // The uncertainty about the number of photons in an INNER PIXEL + + UInt_t fNumExcludedPixels; + + Byte_t fFlags; + + enum { kBlindPixelMethodValid, kPINDiodeMethodValid, kCombinedMethodValid, + kFluxInsidePlexiglassAvailable }; + +public: + + enum PulserColor_t { kEBlue, kEGreen, kEUV, kECT1 } ; + +private: + + PulserColor_t fColor; + +public: + + MCalibrationChargeCam(const char *name=NULL, const char *title=NULL); + ~MCalibrationChargeCam(); + + void Clear( Option_t *o="" ); + void InitSize( const UInt_t i ); + + // Setters + void SetColor( const PulserColor_t color ) { fColor = color; } + void SetNumPixelsExcluded( const UInt_t n ) { fNumExcludedPixels = n; } + void SetGeomCam( const MGeomCam *geom) { fGeomCam = geom; } + void SetPINDiode( const MCalibrationChargePINDiode *d ) { fPINDiode = d; } + + // Setters only for MC!! + void SetBlindPixelMethodValid( const Bool_t b = kTRUE ); + void SetPINDiodeMethodValid( const Bool_t b = kTRUE ); + + // Getters + Int_t GetSize() const; + UInt_t GetNumPixels() const { return fNumPixels; } + + MCalibrationBlindPix *GetBlindPixel() { return fBlindPixel; } + const MCalibrationBlindPix *GetBlindPixel() const { return fBlindPixel; } + + Float_t GetMeanFluxInsidePlexiglass() const { return fMeanFluxInsidePlexiglass; } + Float_t GetMeanFluxErrInsidePlexiglass() const { return fMeanFluxErrInsidePlexiglass; } + + Bool_t GetConversionFactorFFactor( Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma ); + Bool_t GetConversionFactorBlindPixel( Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma ); + Bool_t GetConversionFactorPINDiode( Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma ); + Bool_t GetConversionFactorCombined( Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma ); + + Bool_t IsPixelUsed(Int_t idx) const; + Bool_t IsPixelFitted(Int_t idx) const; + + Bool_t IsBlindPixelMethodValid() const; + Bool_t IsPINDiodeMethodValid() const; + + Bool_t IsFluxInsidePlexiglassAvailable() const; + + // Others + MCalibrationPix &operator[](UInt_t i); + const MCalibrationPix &operator[](UInt_t i) const; + + void CutEdges(); + Bool_t CheckBounds(Int_t i) const; + + // Prints + void Print(Option_t *o="") const; + + // Draws + void DrawPixelContent(Int_t num) const; +// void DrawHiLoFits(); + + // Others + Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; + + Bool_t CalcFluxInsidePlexiglass(); + void ApplyPINDiodeCalibration(); + + ClassDef(MCalibrationChargeCam, 1) // Container for calibration information of the camera +}; + +#endif + + + Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.cc (revision 3247) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.cc (revision 3247) @@ -0,0 +1,377 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationChargePINDiode // +// // +// This is the storage container to hold informations about the pedestal // +// (offset) value of one Pixel (PMT). // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MCalibrationChargePINDiode.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MCalibrationChargePINDiode); + +using namespace std; + +const Float_t MCalibrationChargePINDiode::fgChargeLimit = 3.; +const Float_t MCalibrationChargePINDiode::fgChargeErrLimit = 0.; +const Float_t MCalibrationChargePINDiode::fgChargeRelErrLimit = 1.; + +const Float_t MCalibrationChargePINDiode::fgConversionChargePhotons = -1.; +const Float_t MCalibrationChargePINDiode::fgConversionChargePhotonsErr = -1.; +// +// Area of Inner Pixel w.r.t. PIN Diode (which is 1 cm2) +// +// Distance of PIN Diode to pulser D1: 1.5 +- 0.3 m +// Distance of Inner Pixel to pulser D2: 18.0 +- 0.5 m +// +// +// D1*D1 +// conversion C = ------ = 0.0069 +// D2*D2 +// +// Delta C / C = 2 * Sqrt( (Delta D1/D1)2 + (Delta D2/D2)2 ) +// Delta C / C = 0.4 +// +// C = 0.007 +- 0.003 +// +const Float_t MCalibrationChargePINDiode::gkFluxCameravsPINDiode = 0.007; +const Float_t MCalibrationChargePINDiode::gkFluxCameravsPINDiodeErr = 0.003; +// +// Average QE of the PIN Diode +// +const Float_t MCalibrationChargePINDiode::gkPINDiodeQEGreen = -1.0; +const Float_t MCalibrationChargePINDiode::gkPINDiodeQEBlue = -1.0; +const Float_t MCalibrationChargePINDiode::gkPINDiodeQEUV = -1.0; +const Float_t MCalibrationChargePINDiode::gkPINDiodeQECT1 = -1.0; +// +// Average QE of the PIN Diode +// +const Float_t MCalibrationChargePINDiode::gkPINDiodeQEGreenErr = -1.0; +const Float_t MCalibrationChargePINDiode::gkPINDiodeQEBlueErr = -1.0; +const Float_t MCalibrationChargePINDiode::gkPINDiodeQEUVErr = -1.0; +const Float_t MCalibrationChargePINDiode::gkPINDiodeQECT1Err = -1.0; + +const Float_t MCalibrationChargePINDiode::gkPINDiodeArea = 100; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MCalibrationChargePINDiode::MCalibrationChargePINDiode(const char *name, const char *title) + : fFlags(0) +{ + + fName = name ? name : "MCalibrationChargePINDiode"; + fTitle = title ? title : "Container of the fit results of MHCalibrationChargePINDiode"; + + Clear(); + + SetChargeLimit(); + SetChargeErrLimit(); + SetChargeRelErrLimit(); + + SetConversionChargePhotons(); + SetConversionChargePhotonsErr(); + + SetExcluded(kFALSE); + SetExcludeQualityCheck(kFALSE); +} + +// ------------------------------------------------------------------------ +// +// Invalidate values +// +void MCalibrationChargePINDiode::Clear(Option_t *o) +{ + + SetChargeFitValid ( kFALSE ); + SetTimeFitValid ( kFALSE ); + SetMeanTimeInFirstBin ( kFALSE ); + SetMeanTimeInLastBin ( kFALSE ); + + fMeanCharge = -1.; + fMeanChargeErr = -1.; + fSigmaCharge = -1.; + fSigmaChargeErr = -1.; + fChargeProb = -1.; + fPed = -1.; + fPedRms = -1.; + + fRmsChargeMean = -1.; + fRmsChargeMeanErr = -1.; + fRmsChargeSigma = -1.; + fRmsChargeSigmaErr = -1.; + + fAbsTimeMean = -1.; + fAbsTimeRms = -1.; + + fTimeLowerEdge = -1.; + fTimeUpperEdge = -1.; + + fConvertedPhotons = -1.; + fConvertedPhotonsErr = -1.; + + fMeanFluxOutsidePlexiglass = -1.; + fMeanFluxErrOutsidePlexiglass = -1.; + +} + + +// -------------------------------------------------------------------------- +// +// Set the pedestals from outside +// +void MCalibrationChargePINDiode::SetPedestal(Float_t ped, Float_t pedrms) +{ + + fPed = ped; + fPedRms = pedrms; + +} + +// -------------------------------------------------------------------------- +// +// Set the Excluded Bit from outside +// +void MCalibrationChargePINDiode::SetExcluded(Bool_t b ) +{ + b ? SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded); +} + + +// -------------------------------------------------------------------------- +// +// Set the Excluded Bit from outside +// +void MCalibrationChargePINDiode::SetExcludeQualityCheck(Bool_t b ) +{ + b ? SETBIT(fFlags, kExcludeQualityCheck) : CLRBIT(fFlags, kExcludeQualityCheck); +} + +// -------------------------------------------------------------------------- +// +// Set the Excluded Bit from outside +// +void MCalibrationChargePINDiode::SetChargeFitValid(Bool_t b ) +{ + b ? SETBIT(fFlags, kChargeFitValid) : CLRBIT(fFlags, kChargeFitValid); +} + +// -------------------------------------------------------------------------- +// +// Set the Excluded Bit from outside +// +void MCalibrationChargePINDiode::SetTimeFitValid(Bool_t b ) +{ + b ? SETBIT(fFlags, kTimeFitValid) : CLRBIT(fFlags, kTimeFitValid); +} + + +void MCalibrationChargePINDiode::SetMeanTimeInFirstBin(const Bool_t b) +{ + b ? SETBIT(fFlags,kMeanTimeInFirstBin) : CLRBIT(fFlags,kMeanTimeInFirstBin); +} + +void MCalibrationChargePINDiode::SetMeanTimeInLastBin(const Bool_t b) +{ + b ? SETBIT(fFlags,kMeanTimeInLastBin) : CLRBIT(fFlags,kMeanTimeInLastBin); +} + +Bool_t MCalibrationChargePINDiode::IsMeanTimeInFirstBin() const +{ + return TESTBIT(fFlags,kMeanTimeInFirstBin); +} + +Bool_t MCalibrationChargePINDiode::IsMeanTimeInLastBin() const +{ + return TESTBIT(fFlags,kMeanTimeInLastBin); +} + +Bool_t MCalibrationChargePINDiode::IsExcluded() const +{ + return TESTBIT(fFlags,kExcluded); +} + +Bool_t MCalibrationChargePINDiode::IsChargeFitValid() const +{ + return TESTBIT(fFlags, kChargeFitValid); +} + +Bool_t MCalibrationChargePINDiode::IsTimeFitValid() const +{ + return TESTBIT(fFlags, kTimeFitValid); +} + +// +// The check return kTRUE if: +// +// 1) PINDiode has a fitted charge greater than fChargeLimit*PedRMS +// 2) PINDiode has a fit error greater than fChargeErrLimit +// 3) PINDiode has a fitted charge greater its fChargeRelErrLimit times its charge error +// 4) PINDiode has a charge sigma bigger than its Pedestal RMS +// +Bool_t MCalibrationChargePINDiode::CheckChargeFitValidity() +{ + + if (TESTBIT(fFlags,kExcludeQualityCheck)) + return kTRUE; + + if (fMeanCharge < fChargeLimit*GetPedRms()) + { + *fLog << warn << "WARNING: Fitted Charge is smaller than " + << fChargeLimit << " Pedestal RMS in PINDiode " << endl; + return kFALSE; + } + + if (fMeanChargeErr < fChargeErrLimit) + { + *fLog << warn << "WARNING: Error of Fitted Charge is smaller than " + << fChargeErrLimit << " in PINDiode " << endl; + return kFALSE; + } + + if (fMeanCharge < fChargeRelErrLimit*fMeanChargeErr) + { + *fLog << warn << "WARNING: Fitted Charge is smaller than " + << fChargeRelErrLimit << "* its error in PINDiode " << endl; + return kFALSE; + } + + if (fSigmaCharge < GetPedRms()) + { + *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in PINDiode " << endl; + return kFALSE; + } + return kTRUE; +} + +// +// The check returns kTRUE if: +// +// The mean arrival time is at least 1.0 slices from the used edge slices +// +Bool_t MCalibrationChargePINDiode::CheckTimeFitValidity() +{ + + if (TESTBIT(fFlags,kExcludeQualityCheck)) + return kTRUE; + + if ( fAbsTimeMean < fTimeLowerEdge+1.) + { + *fLog << warn << "WARNING: Mean ArrivalTime in first extraction bin of the PINDiode " << endl; + SetMeanTimeInFirstBin(); + return kFALSE; + } + + if ( fAbsTimeMean > fTimeUpperEdge-1.) + { + *fLog << warn << "WARNING: Mean ArrivalTime in last extraction bin of the PINDiode " << endl; + SetMeanTimeInLastBin(); + return kFALSE; + } + + return kTRUE; +} + +Bool_t MCalibrationChargePINDiode::CalcFluxOutsidePlexiglass() +{ + + if (IsChargeFitValid()) + return kFALSE; + + // Start calculation of number of photons per mm^2 on the camera + fMeanFluxOutsidePlexiglass = fConvertedPhotons * gkPINDiodeArea; + // Correct for the distance between camera and PIN Diode and for different areas. + fMeanFluxOutsidePlexiglass *= gkFluxCameravsPINDiode; + + // Start calculation of number of photons relative Variance (!!) + fMeanFluxErrOutsidePlexiglass = fConvertedPhotonsErr * fConvertedPhotonsErr + / fConvertedPhotons / fConvertedPhotons ; + fMeanFluxErrOutsidePlexiglass += gkFluxCameravsPINDiodeErr*gkFluxCameravsPINDiodeErr + / gkFluxCameravsPINDiode/gkFluxCameravsPINDiode; + + switch (fColor) + { + case kEGreen: + fMeanFluxOutsidePlexiglass /= gkPINDiodeQEGreen; + fMeanFluxErrOutsidePlexiglass += gkPINDiodeQEGreenErr*gkPINDiodeQEGreenErr + / gkPINDiodeQEGreen/gkPINDiodeQEGreen; + break; + case kEBlue: + fMeanFluxOutsidePlexiglass /= gkPINDiodeQEBlue; + fMeanFluxErrOutsidePlexiglass += gkPINDiodeQEBlueErr*gkPINDiodeQEBlueErr + / gkPINDiodeQEBlue/gkPINDiodeQEBlue; + break; + case kEUV: + fMeanFluxOutsidePlexiglass /= gkPINDiodeQEUV; + fMeanFluxErrOutsidePlexiglass += gkPINDiodeQEUVErr*gkPINDiodeQEUVErr + / gkPINDiodeQEUV/gkPINDiodeQEUV; + break; + case kECT1: + default: + fMeanFluxOutsidePlexiglass /= gkPINDiodeQECT1; + fMeanFluxErrOutsidePlexiglass += gkPINDiodeQECT1Err*gkPINDiodeQECT1Err + / gkPINDiodeQECT1/gkPINDiodeQECT1; + break; + } + + + *fLog << inf << endl; + *fLog << inf << " Mean Photon flux [ph/mm^2] outside Plexiglass: " + << fMeanFluxOutsidePlexiglass << endl; + + if (fMeanFluxOutsidePlexiglass > 0.) + SETBIT(fFlags,kFluxOutsidePlexiglassAvailable); + else + { + CLRBIT(fFlags,kFluxOutsidePlexiglassAvailable); + return kFALSE; + } + + if (fMeanFluxErrOutsidePlexiglass < 0.) + { + *fLog << warn << "Relative Variance on Photon flux outside Plexiglass: " + << fMeanFluxErrOutsidePlexiglass << endl; + CLRBIT(fFlags,kFluxOutsidePlexiglassAvailable); + return kFALSE; + } + + // Finish calculation of errors -> convert from relative variance to absolute error + fMeanFluxErrOutsidePlexiglass = TMath::Sqrt(fMeanFluxErrOutsidePlexiglass); + fMeanFluxErrOutsidePlexiglass *= fMeanFluxOutsidePlexiglass; + + *fLog << inf << " Error on Photon flux [ph/mm^2] outside Plexiglass: " + << fMeanFluxErrOutsidePlexiglass << endl; + *fLog << inf << endl; + + return kTRUE; +} + + + Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.h (revision 3247) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.h (revision 3247) @@ -0,0 +1,154 @@ +#ifndef MARS_MCalibrationChargePINDiode +#define MARS_MCalibrationChargePINDiode + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MCalibrationChargePINDiode : public MParContainer +{ +private: + + static const Float_t fgChargeLimit; // The default limit (in units of PedRMS) for acceptance of the fitted mean charge + static const Float_t fgChargeErrLimit; // The default limit (in units of PedRMS) for acceptance of the fitted charge sigma + static const Float_t fgChargeRelErrLimit; // The default limit (in units of Error of fitted charge) for acceptance of the fitted mean + + static const Float_t fgConversionChargePhotons; // The default mean conversion from PIN Diode charge to the number of incident photons + static const Float_t fgConversionChargePhotonsErr; // The default error of the mean conversion from PIN Diode charge to the number of incident photons + + static const Float_t gkPINDiodeArea; // The Blind Pixel area in mm^2 + + static const Float_t gkFluxCameravsPINDiode ; + static const Float_t gkFluxCameravsPINDiodeErr ; + + static const Float_t gkPINDiodeQEGreen; + static const Float_t gkPINDiodeQEBlue ; + static const Float_t gkPINDiodeQEUV ; + static const Float_t gkPINDiodeQECT1 ; + + static const Float_t gkPINDiodeQEGreenErr; + static const Float_t gkPINDiodeQEBlueErr ; + static const Float_t gkPINDiodeQEUVErr ; + static const Float_t gkPINDiodeQECT1Err ; + + Float_t fChargeLimit; // The limit (in units of PedRMS) for acceptance of the fitted mean charge + Float_t fChargeErrLimit; // The limit (in units of PedRMS) for acceptance of the fitted charge sigma + Float_t fChargeRelErrLimit; // The limit (in units of Error of fitted charge) for acceptance of the fitted mean + + Float_t fConversionChargePhotons; // The mean conversion from PIN Diode charge to the number of incident photons + Float_t fConversionChargePhotonsErr; // The error of the mean conversion from PIN Diode charge to the number of incident photons + + Float_t fPed; // The mean pedestal (from MPedestalPix) + Float_t fPedRms; // The pedestal RMS (from MPedestalPix) + + Float_t fMeanCharge; // The mean charge after the fit + Float_t fMeanChargeErr; // The error of mean charge after the fit + Float_t fSigmaCharge; // The sigma of the mean charge after the fit + Float_t fSigmaChargeErr; // The error of the sigma of the mean charge after the fit + Float_t fChargeProb; // The probability of the fit function + + Float_t fRmsChargeMean; + Float_t fRmsChargeMeanErr; + Float_t fRmsChargeSigma; + Float_t fRmsChargeSigmaErr; + + Byte_t fFlags; // Flag for the set Bits + + Float_t fAbsTimeMean; + Float_t fAbsTimeRms; + + Float_t fTimeLowerEdge; + Float_t fTimeUpperEdge; + + Float_t fConvertedPhotons; + Float_t fConvertedPhotonsErr; + + Float_t fMeanFluxOutsidePlexiglass; // The mean number of photons in an INNER PIXEL outside the plexiglass + Float_t fMeanFluxErrOutsidePlexiglass; // The uncertainty about the number of photons in an INNER PIXEL + +public: + + enum PulserColor_t { kEBlue, kEGreen, kEUV, kECT1 } ; + +private: + + PulserColor_t fColor; + + enum { kExcluded, kExcludeQualityCheck, + kChargeFitValid, kTimeFitValid, + kMeanTimeInFirstBin, kMeanTimeInLastBin, + kFluxOutsidePlexiglassAvailable }; + + Bool_t CheckChargeFitValidity(); + Bool_t CheckTimeFitValidity(); + +public: + + MCalibrationChargePINDiode(const char *name=NULL, const char *title=NULL); + ~MCalibrationChargePINDiode() {} + + void Clear(Option_t *o=""); + + // Setters + void SetColor( const PulserColor_t color ) { fColor = color; } + + void SetMeanCharge ( const Float_t f ) { fMeanCharge = f; } + void SetMeanChargeErr ( const Float_t f ) { fMeanChargeErr = f; } + void SetSigmaCharge ( const Float_t f ) { fSigmaCharge = f; } + void SetSigmaChargeErr ( const Float_t f ) { fSigmaChargeErr = f; } + + void SetRmsChargeMean ( const Float_t f ) { fRmsChargeMean = f; } + void SetRmsChargeMeanErr ( const Float_t f ) { fRmsChargeMeanErr = f; } + void SetRmsChargeSigma ( const Float_t f ) { fRmsChargeSigma = f; } + void SetRmsChargeSigmaErr( const Float_t f ) { fRmsChargeSigmaErr = f; } + + void SetAbsTimeMean ( const Float_t f ) { fAbsTimeMean = f; } + void SetAbsTimeRms ( const Float_t f ) { fAbsTimeRms = f; } + + void SetChargeLimit ( const Float_t f=fgChargeLimit ) { fChargeLimit = f; } + void SetChargeErrLimit ( const Float_t f=fgChargeErrLimit ) { fChargeErrLimit = f; } + void SetChargeRelErrLimit( const Float_t f=fgChargeRelErrLimit ) { fChargeRelErrLimit = f; } + + void SetConversionChargePhotons ( const Float_t f=fgConversionChargePhotons ) { fConversionChargePhotons = f; } + void SetConversionChargePhotonsErr ( const Float_t f=fgConversionChargePhotonsErr ) { fConversionChargePhotonsErr = f; } + + void SetPedestal(Float_t ped, Float_t pedrms); + + void SetExcluded ( const Bool_t b = kTRUE ); + void SetExcludeQualityCheck( const Bool_t b = kTRUE ); + void SetChargeFitValid ( const Bool_t b = kTRUE ); + void SetTimeFitValid ( const Bool_t b = kTRUE ); + + // Setters ONLY for MC: + void SetMeanTimeInFirstBin( const Bool_t b = kTRUE ); + void SetMeanTimeInLastBin ( const Bool_t b = kTRUE ); + + // Getters + Float_t GetMeanCharge() const { return fMeanCharge; } + Float_t GetMeanChargeErr() const { return fMeanChargeErr; } + Float_t GetSigmaCharge() const { return fSigmaCharge; } + Float_t GetSigmaChargeErr() const { return fSigmaChargeErr; } + Float_t GetChargeProb() const { return fChargeProb; } + + Float_t GetMeanFluxOutsidePlexiglass() const { return fMeanFluxOutsidePlexiglass; } + Float_t GetMeanFluxErrOutsidePlexiglass() const { return fMeanFluxErrOutsidePlexiglass; } + + // Pedestals + Float_t GetPed() const { return fPed; } + Float_t GetPedRms() const { return fPedRms; } + + + Bool_t IsExcluded() const; + Bool_t IsChargeFitValid() const; + Bool_t IsTimeFitValid() const; + + Bool_t IsMeanTimeInFirstBin() const; + Bool_t IsMeanTimeInLastBin() const; + + Bool_t CalcFluxOutsidePlexiglass(); + + ClassDef(MCalibrationChargePINDiode, 1) // Container for Calibration PIN Diode +}; + +#endif /* MARS_MCalibrationChargePINDiode */ +