Index: trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.cc (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.cc (revision 2525) @@ -0,0 +1,111 @@ +/* ======================================================================== *\ +! +! * +! * 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-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationBlindPix // +// // +// This is the storage container to hold informations about the calibration// +// blind pixel // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MCalibrationBlindPix.h" +#include "MHCalibrationBlindPixel.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MCalibrationBlindPix); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MCalibrationBlindPix::MCalibrationBlindPix(const char *name, const char *title) + : fHist(NULL) +{ + + fName = name ? name : "MCalibrationBlindPix"; + fTitle = title ? title : "Container of the MHCalibrationBlindPixel and the fit results"; + + fHist = new MHCalibrationBlindPixel(); + + if (!fHist) + *fLog << err << dbginf << "Could not create MHCalibrationBlindPixel " << endl; + + fLambda = fMu0 = fMu1 = fSigma0 = fSigma1 = 0; + fErrLambda = fErrMu0 = fErrMu1 = fErrSigma0 = fErrSigma1 = 0; + + fT = fErrT = 0; +} + +MCalibrationBlindPix::~MCalibrationBlindPix() +{ + delete fHist; +} + +// ------------------------------------------------------------------------ +// +// Invalidate values +// +void MCalibrationBlindPix::Clear(Option_t *o) +{ + fHist->Reset(); +} + +Bool_t MCalibrationBlindPix::FitQ() +{ + if (!fHist->FitSinglePhe()) + return kFALSE; + + fLambda = fHist->GetLambda(); + fMu0 = fHist->GetMu0(); + fMu1 = fHist->GetMu1(); + fSigma0 = fHist->GetSigma0(); + fSigma1 = fHist->GetSigma1(); + + fErrLambda = fHist->GetLambdaErr(); + fErrMu0 = fHist->GetMu0Err(); + fErrMu1 = fHist->GetMu1Err(); + fErrSigma0 = fHist->GetSigma0Err(); + fErrSigma1 = fHist->GetSigma1Err(); + + return kTRUE; +} + + + +Bool_t MCalibrationBlindPix::FitT() +{ + + if(!fHist->FitT()) + return kFALSE; + + fT = fHist->GetMeanT(); + fErrT = fHist->GetMeanTErr(); + + return kTRUE; + +} Index: trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.h (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.h (revision 2525) @@ -0,0 +1,69 @@ +#ifndef MARS_MCalibrationBlindPix +#define MARS_MCalibrationBlindPix + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#include "MHCalibrationBlindPixel.h" + +class MCalibrationBlindPix : public MParContainer +{ +private: + + Float_t fLambda; // The mean charge after the fit + Float_t fMu0; // The position of the pedestal-peak + Float_t fMu1; // The position of the first phe-peak + Float_t fSigma0; // The width of the pedestal-peak + Float_t fSigma1; // The width of the first phe-peak + + Float_t fErrLambda; // The error of the mean charge after the fit + Float_t fErrMu0; // The error of the position of the pedestal-peak + Float_t fErrMu1; // The error of the position of the first phe-peak + Float_t fErrSigma0; // The error of the width of the pedestal-peak + Float_t fErrSigma1; // The error of the width of the first phe-peak + + Float_t fT; // The mean arrival time after the fit + Float_t fErrT; // The error of the mean arrival time after the fit + + MHCalibrationBlindPixel *fHist; // Pointer to the histograms performing the fits, etc. + +public: + + MCalibrationBlindPix(const char *name=NULL, const char *title=NULL); + ~MCalibrationBlindPix(); + + void Clear(Option_t *o=""); + + Float_t GetLambda() const { return fLambda; } + Float_t GetMu0() const { return fMu0; } + Float_t GetMu1() const { return fMu1; } + Float_t GetSigma0() const { return fSigma0; } + Float_t GetSigma1() const { return fSigma1; } + + Float_t GetErrLambda() const { return fErrLambda; } + Float_t GetErrMu0() const { return fErrMu0; } + Float_t GetErrMu1() const { return fErrMu1; } + Float_t GetErrSigma0() const { return fErrSigma0; } + Float_t GetErrSigma1() const { return fErrSigma1; } + + Float_t GetT() const { return fT; } + Float_t GetErrT() const { return fErrT; } + + Bool_t FillQ(Int_t q) { return fHist->FillBPQ(q); } + Bool_t FillT(Int_t t) { return fHist->FillBPT(t); } + Bool_t FillRQvsT(Float_t rq, Int_t t) { return fHist->FillBPQvsN(rq,t); } + + Bool_t IsValid() { return fLambda > 0. || fErrLambda > 0.; } + + Bool_t FitQ(); + Bool_t FitT(); + + virtual void Draw(Option_t *opt="") { fHist->Draw(opt); } + MHCalibrationBlindPixel *GetHist() const { return fHist; } + + ClassDef(MCalibrationBlindPix, 1) // Storage Container for Calibration information of one pixel +}; + +#endif + Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.cc (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.cc (revision 2525) @@ -0,0 +1,408 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationCalc // +// // +// This is a task which calculates the number of photons from the FADC // +// time slices. At the moment it integrates simply the FADC values. // +// // +// Input Containers: // +// MRawEvtData // +// // +// Output Containers: // +// MCalibrationCam // +// // +// // +// The class MCalibrationCam hold one entry of type MCalibrationPix for // +// every pixel. It is filled in the following way: // +// PreParocess: MalibrationCam::InitSize(577) is called which allocates // +// memory in an TClonesArray of type MCalibrationPix and // +// all pointers to NULL. // +// // +// Process: The NULL pointer is tested on every pixel via // +// MalibrationCam::IsPixelUsed(npix). // +// // +// In case, IsPixelUsed returns NULL, // +// MalibrationCam::AddPixel(npix) is invoked which creates a // +// new MCalibrationPix(npix) in the npix's entry // +// of the TClonesArray. // +// // +// Every MCalibrationPix holds a histogram class, // +// MHCalibrationPixel which itself hold histograms of type: // +// HQ(npix) (distribution of summed FADC time slice entries) // +// HT(npix) (distribution of position of maximum) // +// HQvsN(npix) (distribution of charges vs. event number. // +// // +// PostProcess: All histograms HQ(npix) are fitted to a Gaussian // +// All histograms HT(npix) are fitted to a Gaussian // +// The histogram HBPQ (blind pixel) is fitted to a single // +// PhE fit // +// The histogram HBPT (blind pixel) is fitted to a Gaussian // +// The histograms of the PIN Diode are fitted to Gaussians // +// // +// Fits can be excluded via the commands: // +// MalibrationCam::SetSkipTFits() (skip all time fits) // +// MalibrationCam::SetSkipBPFits() (skip all blind pixel fits) // +// MalibrationCam::SetSkipPDFits() (skip all PIN Diode fits) // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MCalibrationCalc.h" +#include "MCalibrationConfig.h" +#include "MCalibrationFits.h" + +#include "MCalibrationCam.h" +#include "MCalibrationPix.h" +#include "MCalibrationBlindPix.h" +#include "MCalibrationPINDiode.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MH.h" + +#include "MRawRunHeader.h" +#include "MRawEvtData.h" // MRawEvtData::GetNumPixels +#include "MRawEvtPixelIter.h" + +#include "MTime.h" + +ClassImp(MCalibrationCalc); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default constructor. b is the number of slices before the maximum slice, +// a the number of slices behind the maximum slice which is taken as signal. +// +MCalibrationCalc::MCalibrationCalc(const char *name, const char *title) + : fColor(kEBlue) +{ + + fName = name ? name : "MCalibrationCalc"; + fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam "; + + AddToBranchList("MRawEvtData.fHiGainPixId"); + AddToBranchList("MRawEvtData.fLoGainPixId"); + AddToBranchList("MRawEvtData.fHiGainFadcSamples"); + AddToBranchList("MRawEvtData.fLoGainFadcSamples"); + + SETBIT(fFlags, kUseTFits); + SETBIT(fFlags, kUseBPFit); + SETBIT(fFlags, kUsePDFit); +} + +// -------------------------------------------------------------------------- +// +// The PreProcess searches for the following input containers: +// - MRawRunHeader +// - MRawEvtData +// - MPedestalCam +// +// The following output containers are also searched and created if +// they were not found: +// +// - MHCalibrationBlindPixel +// - MCalibrationCam +// - MTime +// +Int_t MCalibrationCalc::PreProcess(MParList *pList) +{ + + fOverFlow = 0; + fNrEvents = 0; + + fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); + if (!fRawEvt) + { + *fLog << 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; + } + + fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam"); + if (!fCalibrations) + { + *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl; + return kFALSE; + } + + // + // fRawEvt->GetNumPixels() does not work !!!!!!!!!! + // + fCalibrations->InitSize(577); + + switch (fColor) + { + case kEBlue: + fCalibrations->SetColor(MCalibrationCam::kECBlue); + break; + case kEGreen: + fCalibrations->SetColor(MCalibrationCam::kECGreen); + break; + case kEUV: + fCalibrations->SetColor(MCalibrationCam::kECUV); + } + + + + fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPedestals) + { + *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl; + return kFALSE; + } + + + // MTime does not work!!!! + + // fEvtTime = (MTime*)pList->FindObject("MRawEvtTime"); + // if (!fEvtTime) + // { + // *fLog << dbginf << "MTime could not be created ... ¡!¡!¡!¡!" << endl; + // } + + // fEvtTime->SetTime(0.); + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculate the integral of the FADC time slices and store them as a new +// pixel in the MCerPhotEvt container. +// +Int_t MCalibrationCalc::Process() +{ + + fNrEvents++; + + MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); + MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode()); + + MRawEvtPixelIter pixel(fRawEvt); + + while (pixel.Next()) + { + + UShort_t sat = 0; + const Int_t pixid = pixel.GetPixelId(); + + if (!fCalibrations->IsPixelUsed(pixid)) + if (!fCalibrations->AddPixel(pixid)) + *fLog << err << dbginf << "MCalibrationPix(" << pixid << ") could not be created !!" << endl; + + Byte_t *ptr = pixel.GetHiGainSamples(); + Byte_t mid = pixel.GetIdxMaxHiGainSample(); + UInt_t max = pixel.GetMaxHiGainSample(); + + Int_t sum = (max > gkSaturationLimit // overflow ? + ? sat++, pixel.GetSumLoGainSamples() // take Low Gain + : pixel.GetSumHiGainSamples()); // no overflow + + if (sat) + { + + ptr = pixel.GetLoGainSamples(); + max = pixel.GetMaxLoGainSample(); + mid = pixel.GetIdxMaxLoGainSample(); + + sum = (max == gkSaturationLimit // overflow of LoGain ??? -> GimmeABreak!!! + ? fOverFlow++, gkLoGainOverFlow // OUCH (Florian was maybe right) + : sum*gkConversionHiLo ); // OUFF (Florian was wrong) !! + + // *fLog << warn << "Warning: Saturation of HiGain reached in slice " << (int)mid << " !!! " << endl; + // *fLog << warn << "Max = " << max << endl; + + if (fOverFlow) + *fLog << err << dbginf << "Warning: Saturation of LoGain reached! " + << err << dbginf << "sum = " << sum << endl; + + } + + // + // sanity check (old MC files sometimes have pixids>577) + // + if (fPedestals && !fPedestals->CheckBounds(pixid)) + { + *fLog << inf << "Pixel ID larger than camera... skipping event." << endl; + return kCONTINUE; + } + + MPedestalPix &ped = (*fPedestals)[pixid]; + MCalibrationPix &pix = (*fCalibrations)[pixid]; + + Float_t pedes = ped.GetPedestal(); + + // + // FIXME: This is preliminary, we will change to pedestals per slice!!! + // Assume pedestals per time slice ==> multiply with number of slices + // + pedes *= (sat ? fRawEvt->GetNumLoGainSamples() : fRawEvt->GetNumHiGainSamples() ); + + Float_t rsum = (float)sum - pedes; + + switch(pixid) + { + + case gkCalibrationBlindPixelId: + if (!blindpixel.FillQ(sum)) + *fLog << warn << + "Overflow or Underflow occurred filling Blind Pixel means = " << sum << endl; + if (!blindpixel.FillT((int)mid)) + *fLog << warn << + "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl; + if (!blindpixel.FillRQvsT(rsum,fNrEvents)) + *fLog << warn << + "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fNrEvents << endl; + + case gkCalibrationPINDiodeId: + if (!pindiode.FillQ(sum)) + *fLog << warn << + "Overflow or Underflow occurred filling HQ: means = " << sum << endl; + if (!pindiode.FillT((int)mid)) + *fLog << warn << + "Overflow or Underflow occurred filling HT: time = " << (int)mid << endl; + if (!pindiode.FillRQvsT(rsum,fNrEvents)) + *fLog << warn << + "Overflow or Underflow occurred filling HQvsN: eventnr = " << fNrEvents << endl; + + default: + + if (!pix.FillQ(sum)) + *fLog << warn << "Could not fill Q of pixel: " << pixid + << " signal = " << sum << endl; + + // + // Fill the reduced charge into the control histo for better visibility + // + if (!pix.FillRQvsT(rsum,fNrEvents)) + *fLog << warn << "Could not fill red. Q vs. EvtNr of pixel: " << pixid + << " signal = " << rsum << " event Nr: " << fNrEvents << endl; + + + if (!pix.FillT((int)mid)) + *fLog << warn << "Could not fill T of pixel: " << pixid << " time = " << (int)mid << endl; + + + } /* switch(pixid) */ + + } /* while (pixel.Next()) */ + + fCalibrations->SetReadyToSave(); + + return kTRUE; +} + +Int_t MCalibrationCalc::PostProcess() +{ + + *fLog << inf << endl; + *fLog << GetDescriptor() << " Cut Histogram Edges" << endl; + // + // Cut edges to make fits and viewing of the hists easier + // + fCalibrations->CutEdges(); + + // + // Get pointer to blind pixel + // + MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); + + *fLog << GetDescriptor() << " Fitting the Blind Pixel" << endl; + + // + // Fit the blind pixel + // + if (TESTBIT(fFlags,kUseBPFit)) + { + if (!blindpixel.FitQ()) + *fLog << err << dbginf << "Could not fit the blind pixel " << endl; + + if (!blindpixel.FitT()) + *fLog << warn << "Could not the Times of the blind pixel " << endl; + + blindpixel.Draw(); + + } + + *fLog << GetDescriptor() << " Fitting the Normal Pixels" << endl; + // + // loop over the pedestal events and check if we have calibration + // + for (Int_t pixid=0; pixidGetSize(); pixid++) + { + + if (fCalibrations->IsPixelUsed(pixid)) + { + + MCalibrationPix &pix = (*fCalibrations)[pixid]; + + + if (TESTBIT(fFlags,kUseTFits)) + pix.FitT(); + + if (!pix.FitQ()) + continue; + + const Float_t ped = (*fPedestals)[pixid].GetPedestal(); + const Float_t prms = (*fPedestals)[pixid].GetPedestalRms(); + pix.SetPedestal(ped,prms); + } + } + + fCalibrations->SetReadyToSave(); + + if (GetNumExecutions()==0 || fOverFlow==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << setw(7) << fOverFlow << " (" << setw(3) << (int)(fOverFlow*100/GetNumExecutions()) << "%) Evts skipped due to: lo gain saturated." << endl; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.h (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.h (revision 2525) @@ -0,0 +1,77 @@ +#ifndef MARS_MCalibrationCalc +#define MARS_MCalibrationCalc + +///////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationCalc // +// // +// Integrates the time slices of the all pixels of a calibration event // +// and substract the pedestal value // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MRawEvtData; +class MRawRunHeader; + +class MPedestalCam; +class MCalibrationCam; + +class MTime; + +class MCalibrationCalc : public MTask +{ +private: + + MPedestalCam *fPedestals; // Pedestals of all pixels in the camera + MCalibrationCam *fCalibrations; // 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 + + Int_t fNrEvents; // Number of events + Int_t fOverFlow; // Number of events with saturated Low Gain + + Byte_t fFlags; // Flag for the fits used + + enum + { + kUseTFits = 1, + kUseBPFit = 2, + kUsePDFit = 3 + }; + +public: + + enum PulserColor_t { kEGreen, kEBlue, kEUV }; + +private: + + PulserColor_t fColor; + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + + MCalibrationCalc(const char *name=NULL, const char *title=NULL); + + void SetSkipTFits(Bool_t b=kTRUE) + {b ? CLRBIT(fFlags, kUseTFits) : SETBIT(fFlags, kUseTFits);} + void SetSkipBPFit(Bool_t b=kTRUE) + {b ? CLRBIT(fFlags, kUseBPFit) : SETBIT(fFlags, kUseBPFit);} + void SetSkipPDFit(Bool_t b=kTRUE) + {b ? CLRBIT(fFlags, kUsePDFit) : SETBIT(fFlags, kUsePDFit);} + + void SetPulserColor(PulserColor_t color) { fColor = color; } + + ClassDef(MCalibrationCalc, 1) // Task to fill the Calibration Containers from raw data +}; + +#endif Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc (revision 2525) @@ -0,0 +1,503 @@ +/* ======================================================================== *\ +! +! * +! * 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-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationCam // +// // +// Hold the Calibration information for all pixels in the camera // +// // +// The Classes MCalibrationPix's are stored in a TClonesArray // +// The Class MCalibrationBlindPix and the MCalibrationPINDiode // +// are stored in separate pointers // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MCalibrationCam.h" +#include "MCalibrationPix.h" +#include "MCalibrationBlindPix.h" +#include "MCalibrationConfig.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MGeomCam.h" + +ClassImp(MCalibrationCam); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default constructor. +// Creates a TClonesArray of MHCalibrationPix objects, initialized to 0 entries +// Creates an MCalibrationBlindPix and an MCalibrationPINDiode +// +MCalibrationCam::MCalibrationCam(const char *name, const char *title) + : fMeanNrPhotAvailable(kFALSE) +{ + fName = name ? name : "MCalibrationCam"; + fTitle = title ? title : "Storage container for the Calibration Information in the camera"; + + fPixels = new TClonesArray("MCalibrationPix",0); + fBlindPixel = new MCalibrationBlindPix(); + fPINDiode = new MCalibrationPINDiode(); +} + +// -------------------------------------------------------------------------- +// +// Delete the TClonesArray or MCalibrationPix's +// Delete the MCalibrationPINDiode and the MCalibrationBlindPix +// +MCalibrationCam::~MCalibrationCam() +{ + + // + // delete fPixels should delete all Objects stored inside + // + delete fPixels; + delete fBlindPixel; + delete fPINDiode; +} + +// -------------------------------------------------------------------------- +// +// This function return the size of the FILLED MCalibrationCam +// It is NOT the size of the array fPixels !!! +// Note that with every call to AddPixels, GetSize() might change +// +Int_t MCalibrationCam::GetSize() const +{ + + // + // Here it is important to use GetEntriesFast, + // We get the array index of the last "filled" fPixel entry + // (Not the number of "filled" fPixels!!) + // + return fPixels->GetEntriesFast(); +} + +// -------------------------------------------------------------------------- +// +// Get i-th pixel (pixel number) +// +MCalibrationPix &MCalibrationCam::operator[](Int_t i) +{ + return *static_cast(fPixels->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th pixel (pixel number) +// +MCalibrationPix &MCalibrationCam::operator[](Int_t i) const +{ + return *static_cast(fPixels->UncheckedAt(i)); +} + +Bool_t MCalibrationCam::AddPixel(Int_t idx) +{ + + // + // Check bounds first + // + if (!CheckBounds(idx)) + InitSize(idx); + + // + // in case, new cannot allocate memory, + // return FALSE + // + return (new ((*fPixels)[idx]) MCalibrationPix(idx)); + +} + + +// -------------------------------------------------------------------------- +// +// Return a pointer to the pixel with the requested idx. +// NULL if it doesn't exist. +// +MCalibrationPix *MCalibrationCam::GetCalibrationPix(Int_t idx) const +{ + if (idx<0) + return NULL; + + if (!CheckBounds(idx)) + return NULL; + +// return (*this)[idx]; + return (MCalibrationPix*)fPixels->At(idx); +} + + + +// -------------------------------------------------------------------------- +// +// Check if position i is inside bounds +// +Bool_t MCalibrationCam::CheckBounds(Int_t i) const +{ + return i < fPixels->GetEntriesFast(); +} + +Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const +{ + if (!CheckBounds(idx)) + return kFALSE; + + return (&(*this)[idx]); +} + +Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const +{ + if (!IsPixelUsed(idx)) + return kFALSE; + + return ((*this)[idx].GetRQ() > 0. && (*this)[idx].GetErrRQ() > 0.); +} + + + +void MCalibrationCam::Clear(Option_t *o) +{ + fPixels->ForEach(TObject, Clear)(); +} + +// +// Perform the fits of the charges +// If i=-1, then all pixels will be fitted +// Otherwise only the one with index i +// +// The number of succesful fits is returned +// +UShort_t MCalibrationCam::FitQ(Int_t i) +{ + + UShort_t nsuccess = 0; + + // invalidate i if it exceeds the number of entries in fPixels + if (i > fPixels->GetEntriesFast()) + { + *fLog << warn << "Tried to fit pixel out of allowed range " << endl; + return 0; + } + + if (i == -1) // loop over all events + { + + Int_t id = 0; + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + if (!IsPixelUsed(id++)) + continue; + if (pix->FitQ()) + nsuccess++; + } + } + else // fit only the pixel with index i + { + if (!IsPixelUsed(i)) + return 0; + if((*this)[i].FitQ()) + nsuccess++; + } + + return nsuccess; + +} + +// +// Perform the fits of the charges of all pixels +// plus the blind pixel and the PIN Diode +// +// The number of succesful fits is returned +// +UShort_t MCalibrationCam::FitAllQ() +{ + + // FIXME: Once the blind pixel is fully working, + // there must be some penalty in case the + // fit does not succeed + // + UShort_t nsuccess = 0; + + Int_t id = 0; + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + if (!IsPixelUsed(id++)) + continue; + if (pix->FitQ()) + nsuccess++; + } + + if (fBlindPixel->FitQ()) + nsuccess++; + + if (fPINDiode->FitQ()) + nsuccess++; + + return nsuccess; + +} + + + +// +// Perform the fits of the arrival times +// If i=-1, then all pixels will be fitted +// Otherwise only the one with index i +// +// The number of succesful fits is returned +// +UShort_t MCalibrationCam::FitT(Int_t i) +{ + + UShort_t nsuccess = 0; + + // invalidate i if it exceeds the number of entries in fPixels + if (i > fPixels->GetEntriesFast()) + { + *fLog << warn << "Tried to fit pixel out of allowed range " << endl; + return 0; + } + + if (i == -1) // loop over all events + { + + Int_t id = 0; + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + if (!IsPixelUsed(id++)) + continue; + if (pix->FitT()) + nsuccess++; + } + } + else // fit only the pixel with index i + { + if((*this)[i].FitT()) + nsuccess++; + } + + return nsuccess; + +} + + +// +// Perform the fits of the times of all pixels +// plus the blind pixel and the PIN Diode +// +// The number of succesful fits is returned +// +UShort_t MCalibrationCam::FitAllT() +{ + + // FIXME: Once the blind pixel is fully working, + // there must be some penalty in case the + // fit does not succeed + // + UShort_t nsuccess = 0; + + Int_t id = 0; + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + if (!IsPixelUsed(id++)) + continue; + if (pix->FitT()) + nsuccess++; + } + + if (fBlindPixel->FitT()) + nsuccess++; + + if (fPINDiode->FitT()) + nsuccess++; + + return nsuccess; + +} + + + +void MCalibrationCam::CutEdges() +{ + + fBlindPixel->GetHist()->CutAllEdges(); + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + pix->GetHist()->CutAllEdges(); + } + + return; +} + +// --------------------------------------------------------------------- +// +// 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: 0) +// +void MCalibrationCam::InitSize(Int_t size) +{ + + // + // check if we have already initialized to size + // + if (CheckBounds(size)) + return; + + // + // It is important to use Expand and NOT ExpandCreate, because + // we want to keep all pixel not used with a NULL pointer. + // + fPixels->Expand(size); + // + // Set all entries to the null pointer. + // Later we fill the array per pixId with the contruction: new (fPixels[i]) MCalibrationPix(pixid) + // To check, if pixels is already filled, we test the NULL pointer (see IsPixelUsed) + // + for (Int_t i=0; i< size; i++) + { + MCalibrationPix *pix = &(*this)[i]; + pix = NULL; + } +} + +void MCalibrationCam::Print(Option_t *o) const +{ + *fLog << all << GetDescriptor() << ":" << endl; + int id = 0; + + TIter Next(fPixels); + MCalibrationPix *pix; + while ((pix=(MCalibrationPix*)Next())) + { + if (!IsPixelUsed(id)) + continue; + + *fLog << id << ": "; + *fLog << pix->GetQ() << " " << pix->GetRQ() << endl; + + id++; + } +} + + +Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + + if (!IsPixelFitted(idx)) + return kFALSE; + + switch (type) + { + case 0: + val = (*this)[idx].GetRQ(); + break; + case 1: + val = (*this)[idx].GetErrRQ(); + break; + } + return val>=0; +} + +void MCalibrationCam::DrawPixelContent(Int_t num) const +{ + *fLog << warn << "MCalibrationCam::DrawPixelContent - not available." << endl; +} + +Bool_t MCalibrationCam::CalcNrPhotInnerPixel() +{ + if (!fBlindPixel->IsValid()) + return kFALSE; + + const Float_t mean = fBlindPixel->GetLambda(); + // const Float_t merr = fBlindPixel->GetErrLambda(); + + switch (fColor) + { + case kECGreen: + fMeanNrPhotInnerPix = mean / + (gkCalibrationBlindPixelQEGreen + *TMath::Power(10,gkCalibrationBlindPixelAttGreen) + *gkCalibrationBlindPixelArea); + break; + case kECBlue: + fMeanNrPhotInnerPix = mean / + (gkCalibrationBlindPixelQEBlue + *TMath::Power(10,gkCalibrationBlindPixelAttBlue) + *gkCalibrationBlindPixelArea); + break; + case kECUV: + fMeanNrPhotInnerPix = mean / + (gkCalibrationBlindPixelQEUV + *TMath::Power(10,gkCalibrationBlindPixelAttUV) + *gkCalibrationBlindPixelArea); + break; + } + + fMeanNrPhotAvailable = kTRUE; + return kTRUE; +} + + + +Bool_t MCalibrationCam::GetConversionFactor(Int_t ipx, Float_t &mean, Float_t &err) +{ + + if (ipx < 0 || !IsPixelFitted(ipx)) + return kFALSE; + + if (!fMeanNrPhotAvailable) + if (!CalcNrPhotInnerPixel()) + return kFALSE; + + mean = fMeanNrPhotInnerPix / (*this)[ipx].GetRQ(); + + // + // Not yet ready , sorry + // + err = 100000000000.; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h (revision 2525) @@ -0,0 +1,100 @@ +#ifndef MARS_MCalibrationCam +#define MARS_MCalibrationCam + +#ifndef MARS_MCamEvent +#include "MCamEvent.h" +#endif + +#ifndef MARS_MCalibrationPix +#include "MCalibrationPix.h" +#endif + +#ifndef MARS_MCalibrationBlindPix +#include "MCalibrationBlindPix.h" +#endif + +#ifndef MARS_MCalibrationPINDiode +#include "MCalibrationPINDiode.h" +#endif + +#ifndef ROOT_TClonesArray +#include +#endif + + +class MGeomCam; +class MCalibrationPix; +class MCalibrationBlindPix; +class MCalibrationPINDiode; + +class MCalibrationCam : public MCamEvent +{ +private: + + Int_t fNumPixels; + TClonesArray *fPixels; // FIXME: Change TClonesArray away from a pointer? + + MCalibrationBlindPix *fBlindPixel; // containing blind pixel data with fitresults + MCalibrationPINDiode *fPINDiode; // containing PIN Diode data with fit results + + Float_t fMeanNrPhotInnerPix; // The mean number of photons in an inner pixel + Float_t fMeanNrPhotInnerPixErr; // The uncertainty about the number of photons in an inner pixel + Bool_t fMeanNrPhotAvailable; + + Bool_t CalcNrPhotInnerPixel(); + +public: + + enum CalibrationColor_t { kECGreen, kECBlue, kECUV }; + +private: + + CalibrationColor_t fColor; + +public: + + MCalibrationCam(const char *name=NULL, const char *title=NULL); + ~MCalibrationCam(); + + void Clear(Option_t *o=""); + + void InitSize(const Int_t i); + Int_t GetSize() const; + + UInt_t GetNumPixels() const { return fNumPixels; } + + Bool_t IsPixelUsed(Int_t idx) const; + Bool_t IsPixelFitted(Int_t idx) const; + + MCalibrationPix &operator[](Int_t i); + MCalibrationPix &operator[](Int_t i) const; + + UShort_t FitQ(Int_t i=-1); + UShort_t FitAllQ(); + UShort_t FitT(Int_t i=-1); + UShort_t FitAllT(); + + Bool_t CheckBounds(Int_t i) const; + + Bool_t AddPixel(Int_t idx); + + void Print(Option_t *o="") const; + + void CutEdges(); + + Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; + void DrawPixelContent(Int_t num) const; + + MCalibrationPix *GetCalibrationPix(Int_t idx) const; + MCalibrationBlindPix *GetBlindPixel() const { return fBlindPixel; } + MCalibrationPINDiode *GetPINDiode() const { return fPINDiode; } + + void SetColor(CalibrationColor_t color) { fColor = color; } + + Bool_t GetConversionFactor(Int_t ipx, Float_t &mean, Float_t &err); + + ClassDef(MCalibrationCam, 1) // Storage Container for all calibration information of the camera +}; + +#endif + Index: trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h (revision 2525) @@ -0,0 +1,45 @@ +#ifndef MARS_MCalibrationConfig +#define MARS_MCalibrationConfig + +///////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationConfig // +// // +// Contains all configuration data of the Calibration // +// // +///////////////////////////////////////////////////////////////////////////// + + +// The highest value of the High Gain in the linear regime +const Byte_t gkSaturationLimit = 250; + +// The conversion factor between High Gain and Low Gain +const UShort_t gkConversionHiLo = 10; + +// The penalty constant to produce overflow in the histogram +const ULong_t gkLoGainOverFlow = 99999999; + +// ----- BLIND PIXEL ----------------------// + +// Pixel ID of the Calibration Blind Pixel +const UShort_t gkCalibrationBlindPixelId = 559; + +// Average QE of Blind Pixel (three colours) +const Float_t gkCalibrationBlindPixelQEGreen = 15.4; +const Float_t gkCalibrationBlindPixelQEBlue = 22.6; +const Float_t gkCalibrationBlindPixelQEUV = 24.7; + +// Attenuation factor Blind Pixel (three colours) +const Float_t gkCalibrationBlindPixelAttGreen = 1.97; +const Float_t gkCalibrationBlindPixelAttBlue = 1.96; +const Float_t gkCalibrationBlindPixelAttUV = 1.95; + +// Area of Blind Pixel w.r.t. Inner Pixel +const Float_t gkCalibrationBlindPixelArea = 0.25; + +// ----- PIN DIODE ------------------------// + +// Pixel ID of the Calibration PIN Diode +const UShort_t gkCalibrationPINDiodeId = 9999; + +#endif /* MARS_MCalibrationBlindPixelConfig */ Index: trunk/MagicSoft/Mars/manalysis/MCalibrationFits.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationFits.h (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationFits.h (revision 2525) @@ -0,0 +1,310 @@ +#ifndef MARS_MCalibrationFits +#define MARS_MCalibrationFits + +#ifndef ROOT_TMath +#include +#endif + +#define GIMMEABREAK 10000000000.0 + +inline Double_t gfKto5(Double_t *x, Double_t *par) +{ + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return GIMMEABREAK; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma1 < sigma0) + return GIMMEABREAK; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + return par[5]*sum; + +}; + +inline Double_t gfKto6(Double_t *x, Double_t *par) +{ + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return GIMMEABREAK; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma1 < sigma0) + return GIMMEABREAK; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + Double_t mu6 = (6.*mu1)-(5.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + Double_t lambda6 = lambda5*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + // k=6: + arg = (x[0] - mu6)/sigma6; + sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6; + + return par[5]*sum; + +}; + +inline Double_t gfKto7(Double_t *x, Double_t *par) +{ + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return GIMMEABREAK; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma1 < sigma0) + return GIMMEABREAK; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + Double_t mu6 = (6.*mu1)-(5.*mu0); + Double_t mu7 = (7.*mu1)-(6.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0)); + Double_t sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + Double_t lambda6 = lambda5*lambda; + Double_t lambda7 = lambda6*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + // k=6: + arg = (x[0] - mu6)/sigma6; + sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6; + + // k=7: + arg = (x[0] - mu7)/sigma7; + sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7; + + return par[5]*sum; + +}; + + +inline Double_t gfKto8(Double_t *x, Double_t *par) +{ + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return GIMMEABREAK; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma1 < sigma0) + return GIMMEABREAK; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + Double_t mu6 = (6.*mu1)-(5.*mu0); + Double_t mu7 = (7.*mu1)-(6.*mu0); + Double_t mu8 = (8.*mu1)-(7.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0)); + Double_t sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0)); + Double_t sigma8 = TMath::Sqrt((8.*sigma1*sigma1) - (7.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + Double_t lambda6 = lambda5*lambda; + Double_t lambda7 = lambda6*lambda; + Double_t lambda8 = lambda7*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + // k=6: + arg = (x[0] - mu6)/sigma6; + sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6; + + // k=7: + arg = (x[0] - mu7)/sigma7; + sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7; + + // k=8: + arg = (x[0] - mu8)/sigma8; + sum += 0.000024801587315*lambda8*TMath::Exp(-0.5*arg*arg)/sigma7; + + return par[5]*sum; + +}; + +#endif /* MARS_MCalibrationFits */ + Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.cc (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.cc (revision 2525) @@ -0,0 +1,100 @@ +/* ======================================================================== *\ +! +! * +! * 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-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationPINDiode // +// // +// This is the storage container to hold informations about the pedestal // +// (offset) value of one Pixel (PMT). // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MCalibrationPINDiode.h" + +#include "MLog.h" + +ClassImp(MCalibrationPINDiode); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MCalibrationPINDiode::MCalibrationPINDiode(const char *name, const char *title) + : fHist(NULL) +{ + + fName = name ? name : "MCalibrationPINDiode"; + fTitle = title ? title : "Container of the MHCalibrationPINDiode and the fit results"; + + fHist = new MHCalibrationPINDiode(); + + fQ = fErrQ = 0.; + fPed = fPedRms = 0.; + fT = fErrT = 0.; + fRQ = fErrRQ = 0.; + fSigmaQ = fErrSigmaQ = 0.; + +} + +MCalibrationPINDiode::~MCalibrationPINDiode() +{ + delete fHist; +} + +// ------------------------------------------------------------------------ +// +// Invalidate values +// +void MCalibrationPINDiode::Clear(Option_t *o) +{ + fHist->Reset(); +} + +Bool_t MCalibrationPINDiode::FitQ() +{ + if(!fHist->FitQ()) + return kFALSE; + + fQ = fHist->GetQMean(); + fErrQ = fHist->GetQMeanErr(); + fSigmaQ = fHist->GetQSigma(); + fErrSigmaQ = fHist->GetQSigmaErr(); + + return kTRUE; + +} + +Bool_t MCalibrationPINDiode::FitT() +{ + + if(!fHist->FitT()) + return kFALSE; + + fT = fHist->GetT(); + fErrT = fHist->GetErrT(); + + return kTRUE; + +} Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.h (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.h (revision 2525) @@ -0,0 +1,53 @@ +#ifndef MARS_MCalibrationPINDiode +#define MARS_MCalibrationPINDiode + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#include "MHCalibrationPINDiode.h" + +class MCalibrationPINDiode : public MParContainer +{ +private: + + Float_t fQ; // The mean charge after the fit + Float_t fErrQ; // The error of mean charge after the fit + Float_t fPed; // The mean pedestal (from MPedestalPix) + Float_t fPedRms; // The pedestal RMS (from MPedestalPix) + Float_t fSigmaQ; // The sigma of the mean charge after the fit + Float_t fErrSigmaQ; // The error of the sigma of the mean charge after the fit + Float_t fT; // The mean arrival time after the fit + Float_t fErrT; // The error of the mean arrival time after the fit + + Float_t fRQ; // The reduced mean charge after the fit + Float_t fErrRQ; // The error of the reduced mean charge after the fit + + MHCalibrationPINDiode *fHist; // Pointer to the histograms performing the fits, etc. + +public: + + MCalibrationPINDiode(const char *name=NULL, const char *title=NULL); + ~MCalibrationPINDiode(); + + void Clear(Option_t *o=""); + + void SetPed(Float_t ped) { fPed = ped; } + void SetPedRms(Float_t pedrms) { fPedRms = pedrms; } + + Bool_t IsValid() const { return fRQ >=0 || fErrRQ >= 0; } + + Bool_t FillQ(Int_t q) { return fHist->FillQ(q); } + Bool_t FillT(Int_t t) { return fHist->FillT(t); } + Bool_t FillRQvsT(Float_t rq, Int_t t) { return fHist->FillQvsN(rq,t); } + + Bool_t FitQ(); + Bool_t FitT(); + + MHCalibrationPINDiode *GetHist() const { return fHist; } + + ClassDef(MCalibrationPINDiode, 1) // Storage Container for Calibration information of one pixel +}; + +#endif /* MARS_MCalibrationPINDiode */ + Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc (revision 2525) @@ -0,0 +1,139 @@ +/* ======================================================================== *\ +! +! * +! * 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-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MCalibrationPix // +// // +// This is the storage container to hold informations about the pedestal // +// (offset) value of one Pixel (PMT). // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MCalibrationPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MCalibrationPix); + +using namespace std; + +static const TString gsDefHTitle = "Calibration Histograms Pixel "; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MCalibrationPix::MCalibrationPix(Int_t pix, const char *name, const char *title) + : fPixId(pix) +{ + + fName = name ? name : "MCalibrationPixel"; + fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results"; + + fHist = new MHCalibrationPixel(fPixId,"MHCalibrationPixel",gsDefHTitle.Data()+fPixId); + + fQ = fErrQ = 0.; + fPed = fPedRms = 0.; + fT = fErrT = 0.; + fRQ = fErrRQ = 0.; + fSigmaQ = fErrSigmaQ = 0.; + +} + +MCalibrationPix::~MCalibrationPix() +{ + delete fHist; +} + +// ------------------------------------------------------------------------ +// +// Invalidate values +// +void MCalibrationPix::Clear(Option_t *o) +{ + fHist->Reset(); +} + +Bool_t MCalibrationPix::FitQ() +{ + + if (fHist->IsFitted()) + return kTRUE; + + if(!fHist->FitQ()) + { + *fLog << warn << "Could not fit charges of pixel " << fPixId << endl; + fHist->PrintQFitResult(); + return kFALSE; + } + + fQ = fHist->GetQMean(); + fErrQ = fHist->GetQMeanErr(); + fSigmaQ = fHist->GetQSigma(); + fErrSigmaQ = fHist->GetQSigmaErr(); + + if (fPed) + fRQ = fQ - fPed; + if (fPedRms) + fErrRQ = TMath::Sqrt(fErrQ*fErrQ + fPedRms*fPedRms); + + + return kTRUE; + +} + +void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms) +{ + + fPed = ped; + fPedRms = pedrms; + + if ((!fRQ) && fQ) + fRQ = fQ - fPed; + if ((!fErrRQ) && fErrQ) + fErrRQ = TMath::Sqrt(fErrQ*fErrQ + fPedRms*fPedRms); + +} + +Bool_t MCalibrationPix::FitT() +{ + + if(!fHist->FitT()) + { + *fLog << warn << "Could not fit times of pixel " << fPixId << endl; + fHist->PrintTFitResult(); + return kFALSE; + } + + fT = fHist->GetTMean(); + fErrT = fHist->GetTSigma(); + + return kTRUE; + +} + +void MCalibrationPix::Test() +{ + *fLog << "TEST: pixid: " << fPixId << endl; +} Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h (revision 2525) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h (revision 2525) @@ -0,0 +1,71 @@ +#ifndef MARS_MCalibrationPix +#define MARS_MCalibrationPix + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#include "MHCalibrationPixel.h" + +class MCalibrationPix : public MParContainer +{ +private: + + // + // FIXME: Derive class from MCerphotPix ?? + // + Int_t fPixId; // the pixel Id + + Float_t fQ; // The mean charge after the fit + Float_t fErrQ; // The error of mean charge after the fit + Float_t fPed; // The mean pedestal (from MPedestalPix) + Float_t fPedRms; // The pedestal RMS (from MPedestalPix) + Float_t fSigmaQ; // The sigma of the mean charge after the fit + Float_t fErrSigmaQ; // The error of the sigma of the mean charge after the fit + Float_t fT; // The mean arrival time after the fit + Float_t fErrT; // The error of the mean arrival time after the fit + + Float_t fRQ; // The reduced mean charge after the fit + Float_t fErrRQ; // The error of the reduced mean charge after the fit + + MHCalibrationPixel *fHist; // Pointer to the histograms performing the fits, etc. + +public: + + MCalibrationPix(Int_t pix=-1, const char *name=NULL, const char *title=NULL); + ~MCalibrationPix(); + + void Clear(Option_t *o=""); + + Float_t GetQ() const { return fQ; } + Float_t GetRQ() const { return fRQ; } + Float_t GetErrQ() const { return fErrQ; } + Float_t GetErrRQ() const { return fErrRQ; } + + Float_t GetSigmaQ() const { return fSigmaQ; } + Float_t GetErrSigmaQ() const { return fErrSigmaQ; } + Float_t GetT() const { return fT; } + Float_t GetErrT() const { return fErrT; } + + void SetPedestal(Float_t ped, Float_t pedrms); + + Bool_t FillQ(Int_t q) { return fHist->FillQ(q); } + Bool_t FillT(Int_t t) { return fHist->FillT(t); } + Bool_t FillRQvsT(Float_t rq, Int_t t) { return fHist->FillQvsN(rq,t); } + + Bool_t IsValid() const { return fRQ >=0 || fErrRQ >= 0; } + Int_t GetPixId() const { return fPixId; } + + Bool_t FitQ(); + Bool_t FitT(); + + MHCalibrationPixel *GetHist() const { return fHist; } + virtual void Draw(Option_t *opt="") { fHist->Draw(opt); } + + void Test(); + + ClassDef(MCalibrationPix, 1) // Storage Container for Calibration information of one pixel +}; + +#endif +