Index: /trunk/MagicSoft/Mars/macros/calibration.C =================================================================== --- /trunk/MagicSoft/Mars/macros/calibration.C (revision 2525) +++ /trunk/MagicSoft/Mars/macros/calibration.C (revision 2525) @@ -0,0 +1,197 @@ +/* ======================================================================== *\ +! +! * +! * 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-2003 +! +! +\* ======================================================================== */ + +void calibration(TString pedname="./20031102_02399_P_Unavailable_E.root", + TString calname="./20031102_02400_D_Flip500Hz_E.root") +{ + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + // + // Now setup the tasks and tasklist for the pedestals: + // --------------------------------------------------- + // + + MReadMarsFile read("Events", pedname); + read.DisableAutoScheme(); + + MGeomApply geomapl; + MPedCalcPedRun pedcalc; + // MPedestalCam pedcam; + + tlist.AddToList(&read); + tlist.AddToList(&geomapl); + tlist.AddToList(&pedcalc); + + // plist.AddToList(&pedcam); + + MHCamEvent hist; + hist.SetType(1); + plist.AddToList(&hist); + MFillH fill("MHCamEvent", "MPedestalCam"); + + tlist.AddToList(&fill); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + // + // Execute first analysis + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + // plist.FindObject("MPedestalCam")->Print(); + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist2; + + MTaskList tlist2; + plist2.AddToList(&tlist2); + + plist2.AddToList((MPedestalCam*)plist.FindObject("MPedestalCam")); + + MGeomApply geomapl2; + tlist2.AddToList(&geomapl2); + + // + // Now setup the new tasks and tasklist for the calibration + // --------------------------------------------------- + // + + MReadMarsFile read2("Events", calname); + read2.DisableAutoScheme(); + + MCalibrationCalc calcalc; + calcalc.SetSkipTFits(); + // + // As long, as we don't have digital modules, + // we have to set the color by hand + // + calcalc.SetPulserColor(MCalibrationCalc::kEBlue); + + tlist2.AddToList(&read2); + tlist2.AddToList(&calcalc); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop2; + evtloop2.SetParList(&plist2); + + // + // Execute second analysis + // + if (!evtloop2.Eventloop()) + return; + + tlist2.PrintStatistics(); + + plist2.FindObject("MCalibrationCam")->Print(); + // + // just one example how to get the plots of individual pixels + // + MCalibrationCam *cam = plist2.FindObject("MCalibrationCam"); + MCalibrationPix *pix = cam->GetCalibrationPix(560); + pix->Draw(); + + /* + TCanvas *c = MH::MakeDefCanvas(); + c->Divide(3, 2); + + MHCamEvent &h = *(MHCamEvent*)plist->FindObject("MHCamEvent"); + MHCamera *disp0 = h.GetHistByName(); + MHCamera *disp1 = new MHCamera(geom, "MPedestalCam;avg", "Pedestals Avarage"); + MHCamera *disp2 = new MHCamera(geom, "MPedestalCam;rms", "Pedestals Avarage Error"); + MHCamera *disp3 = new MHCamera(geom, "MPedestalCam;rel", "Pedestals Avarage ERR/VAL"); + disp1->SetBit(kCanDelete); + disp2->SetBit(kCanDelete); + disp3->SetBit(kCanDelete); + + disp1->SetCamContent(*disp0, 0); + disp2->SetCamContent(*disp0, 1); + disp3->SetCamContent(*disp0, 2); + + disp1->SetYTitle("P [au]"); + disp2->SetYTitle("\\sigma_{P} [au]"); + disp3->SetYTitle("\\sigma_{P} [%]"); + + TText text(0.1, 0.5, &fname[fname.Last('/')+1]); + text.SetTextSize(0.015); + text.DrawClone(); + + c->cd(1); + gStyle->SetOptStat(1111); + disp1->Draw("hist"); + gPad->Update(); + + c->cd(2); + gStyle->SetOptStat(1101); + disp2->Draw("hist"); + gPad->Update(); + + c->cd(3); + gStyle->SetOptStat(1101); + disp3->Draw("hist"); + gPad->Update(); + + c->cd(4); + gPad->SetBorderMode(0); + gPad->Divide(1,1); + gPad->cd(1); + disp1->Draw(); + + c->cd(5); + gPad->SetBorderMode(0); + gPad->Divide(1,1); + gPad->cd(1); + disp2->Draw(); + + c->cd(6); + gPad->SetBorderMode(0); + gPad->Divide(1,1); + gPad->cd(1); + disp3->Draw(); + + c->SaveAs(fname(0, fname.Last('.')+1) + "ps"); + //c->SaveAs(fname(0, fname.Last('.')+1) + "root"); + + */ + +} + 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 + Index: /trunk/MagicSoft/Mars/mhist/MH.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MH.cc (revision 2524) +++ /trunk/MagicSoft/Mars/mhist/MH.cc (revision 2525) @@ -882,4 +882,83 @@ } +// -------------------------------------------------------------------------- +// +// Cuts the bins containing only zeros at the edges. +// +// A new number of bins can be defined with nbins != 0 +// In the case of nbins == 0, no rebinning will take place +// +// Returns the new (real) number of bins +// +Int_t MH::CutEdges(TH1 *h, Int_t nbins) +{ + + TAxis* axe = h->GetXaxis(); + + Int_t min1 = axe->GetFirst(); + Int_t max1 = axe->GetLast(); + Int_t range1 = max1-min1; + + // + // Check for useless zeros + // + if (range1 == 0) + return 0; + + Int_t min2 = 0; + Int_t max2 = 0; + + for (int i=min1; i<=max1; i++) + { + Double_t x = h->GetBinContent(i); + if (x != 0.) + { + min2 = i; + break; + } + } + + // + // If the histogram consists of zeros only + // + if (min2 == max1) + return 0; + + for (int i=max1; i>=min2; i--) + { + Double_t x = h->GetBinContent(i); + if (x != 0.) + { + max2 = i; + break; + } + } + + // + // Check for rebinning + // + if (nbins < 1) + { + axe->SetRange(min2,max2); + return axe->GetLast()-axe->GetFirst(); + } + + // + // Appying TAxis->SetRange before ReBin does not work ... + // But this workaround helps quite fine + // + Axis_t min = h->GetBinLowEdge(min2); + Axis_t max = h->GetBinLowEdge(max2)+h->GetBinWidth(max2); + + Int_t ngroup = (int)((max2-min2)*h->GetNbinsX()/nbins/(max1-min1)); + + if (ngroup > 1) + h->Rebin(ngroup); + + axe->SetRangeUser(min,max); + + return axe->GetLast()-axe->GetFirst(); +} + void MH::ProjectionX(TH1D &dest, const TH2 &src, Int_t firstybin, Int_t lastybin) { Index: /trunk/MagicSoft/Mars/mhist/MH.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MH.h (revision 2524) +++ /trunk/MagicSoft/Mars/mhist/MH.h (revision 2525) @@ -88,5 +88,6 @@ static void FindGoodLimits(Int_t nbins, Int_t &newbins, Double_t &xmin, Double_t &xmax, Bool_t isInteger); static Double_t GetMinimumGT(const TH1 &h, Double_t gt=0); - + static Int_t CutEdges(TH1 *h, Int_t nbins); + static void ProjectionX(TH1D &dest, const TH2 &src, Int_t firstybin=-1, Int_t lastybin=9999); static void ProjectionY(TH1D &dest, const TH2 &src, Int_t firstxbin=-1, Int_t lastxbin=9999); Index: /trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.cc (revision 2525) +++ /trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.cc (revision 2525) @@ -0,0 +1,465 @@ +/* ======================================================================== *\ +! +! * +! * 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-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHCalibrationBlindPixel // +// // +// Performs all the necessary fits to extract the mean number of photons // +// out of the derived light flux // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationBlindPixel.h" +#include "MHCalibrationConfig.h" +#include "MCalibrationFits.h" + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHCalibrationBlindPixel); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title) + : fSinglePheFit(NULL), fTimeGausFit(NULL) +{ + + fName = name ? name : "MHCalibrationBlindPixel"; + fTitle = title ? title : "Fill the accumulated charges and times all Blind Pixel events and perform fits"; + + // Create a large number of bins, later we will rebin + fBPQfirst = 0; + fBPQlast = gkStartBlindPixelBinNr; + fBPQnbins = gkStartBlindPixelBinNr; + + fHBPQ = new TH1I("HBPQ","Distribution of Summed FADC Slices",fBPQnbins,fBPQfirst,fBPQlast); + fHBPQ->SetXTitle("Sum FADC Slices"); + fHBPQ->SetYTitle("Nr. of events"); + fHBPQ->Sumw2(); + + fErrBPQfirst = 0.; + fErrBPQlast = gkStartBlindPixelBinNr; + fErrBPQnbins = gkStartBlindPixelBinNr; + + fHBPErrQ = new TH1F("HBPErrQ","Distribution of Variances of Summed FADC Slices", + fErrBPQnbins,fErrBPQfirst,fErrBPQlast); + fHBPErrQ->SetXTitle("Variance Summed FADC Slices"); + fHBPErrQ->SetYTitle("Nr. of events"); + fHBPErrQ->Sumw2(); + + Axis_t tfirst = -0.5; + Axis_t tlast = 15.5; + Int_t nbins = 16; + + fHBPT = new TH1I("HBPT","Distribution of Mean Arrival Times",nbins,tfirst,tlast); + fHBPT->SetXTitle("Mean Arrival Times [FADC slice nr]"); + fHBPT->SetYTitle("Nr. of events"); + fHBPT->Sumw2(); + + // We define a reasonable number and later enlarge it if necessary + nbins = 20000; + Axis_t nfirst = -0.5; + Axis_t nlast = (Axis_t)nbins - 0.5; + + fHBPQvsN = new TH1I("HBPQvsN","Sum of Charges vs. Event Number",nbins,nfirst,nlast); + fHBPQvsN->SetXTitle("Event Nr."); + fHBPQvsN->SetYTitle("Sum of FADC slices"); + + fgSinglePheFitFunc = &gfKto8; + fgSinglePheFitNPar = 5; +} + +MHCalibrationBlindPixel::~MHCalibrationBlindPixel() +{ + + delete fHBPQ; + delete fHBPT; + delete fHBPErrQ; + + if (fSinglePheFit) + delete fSinglePheFit; + if (fTimeGausFit) + delete fTimeGausFit; +} + + + +void MHCalibrationBlindPixel::ResetBin(Int_t i) +{ + fHBPQ->SetBinContent (i, 1.e-20); + fHBPErrQ->SetBinContent (i, 1.e-20); + fHBPT->SetBinContent(i, 1.e-20); +} + + +// ------------------------------------------------------------------------- +// +// Draw a legend with the fit results +// +void MHCalibrationBlindPixel::DrawLegend() +{ + + fFitLegend = new TPaveText(0.05,0.05,0.95,0.95); + + if (fFitOK) + fFitLegend->SetFillColor(80); + else + fFitLegend->SetFillColor(2); + + fFitLegend->SetLabel("Results of the single PhE Fit (to k=6):"); + fFitLegend->SetTextSize(0.05); + + char line1[32]; + sprintf(line1,"Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr()); + fFitLegend->AddText(line1); + + char line2[32]; + sprintf(line2,"Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err()); + fFitLegend->AddText(line2); + + char line3[32]; + sprintf(line3,"Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err()); + fFitLegend->AddText(line3); + + char line4[32]; + sprintf(line4,"1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err()); + fFitLegend->AddText(line4); + + char line5[32]; + sprintf(line5,"Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err()); + fFitLegend->AddText(line5); + + char line7[32]; + sprintf(line7,"#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf()); + fFitLegend->AddText(line7); + + char line8[32]; + sprintf(line8,"Probability: %4.2f ",GetProb()); + fFitLegend->AddText(line8); + + if (fFitOK) + fFitLegend->AddText("Result of the Fit: OK"); + else + fFitLegend->AddText("Result of the Fit: NOT OK"); + + fFitLegend->SetBit(kCanDelete); + fFitLegend->Draw(); + +} + + +// ------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHCalibrationBlindPixel::Draw(Option_t *opt) +{ + + gStyle->SetOptFit(0); + gStyle->SetOptStat(1111111); + + TCanvas *c = MakeDefCanvas(this,550,700); + + c->Divide(2,2); + + gROOT->SetSelectedPad(NULL); + + c->cd(1); + gPad->SetLogy(1); + gPad->SetTicks(); + + fHBPQ->DrawCopy(opt); + + if (fSinglePheFit) + { + if (fFitOK) + fSinglePheFit->SetLineColor(kGreen); + else + fSinglePheFit->SetLineColor(kRed); + + fSinglePheFit->DrawCopy("same"); + c->Modified(); + c->Update(); + } + + c->cd(2); + DrawLegend(); + c->Modified(); + c->Update(); + + c->cd(3); + gPad->SetLogy(1); + gPad->SetBorderMode(0); + fHBPT->DrawCopy(opt); + + if (fHBPT->GetFunction("GausTime")) + { + TF1 *tfit = fHBPT->GetFunction("GausTime"); + if (tfit->GetProb() < 0.01) + tfit->SetLineColor(kRed); + else + tfit->SetLineColor(kGreen); + + tfit->DrawCopy("same"); + c->Modified(); + c->Update(); + } + + c->cd(4); + + fHBPQvsN->DrawCopy(opt); + + c->Modified(); + c->Update(); +} + + + +Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1) +{ + gRandom->SetSeed(); + + if (fHBPQ->GetEntries() != 0) + { + *fLog << err << "Histogram " << fHBPQ->GetTitle() << " is already filled. " << endl; + *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl; + return kFALSE; + } + + TF1 *simulateSinglePhe = new TF1("simulateSinglePhe",fgSinglePheFitFunc, + fBPQfirst,fBPQlast,fgSinglePheFitNPar); + + simulateSinglePhe->SetParameters(lambda,mu0,mu1,sigma0,sigma1); + simulateSinglePhe->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1"); + simulateSinglePhe->SetNpx(fBPQnbins); + + for (Int_t i=0;i<10000; i++) + { + fHBPQ->Fill(simulateSinglePhe->GetRandom()); + } + + return kTRUE; +} + + +void MHCalibrationBlindPixel::ChangeFitFunc(BPFitFunc fitfunc, Int_t par) +{ + + fgSinglePheFitFunc = fitfunc; + fgSinglePheFitNPar = par; + +} + + + +Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt) +{ + + if (fSinglePheFit) + return kFALSE; + + + // + // Get the fitting ranges + // + rmin = (rmin != 0.) ? rmin : fBPQfirst; + rmax = (rmax != 0.) ? rmax : fBPQlast; + + // + // First guesses for the fit (should be as close to reality as possible, + // otherwise the fit goes gaga because of high number of dimensions ... + // + const Stat_t entries = fHBPQ->GetSumOfWeights(); + const Double_t lambda_guess = 0.2; + const Double_t mu_0_guess = fHBPQ->GetBinCenter(fHBPQ->GetMaximumBin()); + const Double_t si_0_guess = mu_0_guess/500.; + const Double_t mu_1_guess = mu_0_guess + 2.; + const Double_t si_1_guess = si_0_guess + si_0_guess; + + fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar+1); + fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,entries); + fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","area"); + fSinglePheFit->SetParLimits(0,0.,5.); + fSinglePheFit->SetParLimits(1,rmin,rmax); + fSinglePheFit->SetParLimits(2,rmin,rmax); + fSinglePheFit->SetParLimits(3,1.0,rmax-rmin); + fSinglePheFit->SetParLimits(4,1.7,rmax-rmin); + fSinglePheFit->SetParLimits(5,0.,2.*entries); + // + // Normalize the histogram to facilitate faster fitting of the area + // For speed reasons, FKto8 is normalized to Sqrt(2 pi). + // Therefore also normalize the histogram to that value + // + // ROOT gives us another nice example of user-unfriendly behavior: + // Although the normalization of the function fSinglePhe and the + // Histogram fHBPQ agree (!!), the fit does not normalize correctly INTERNALLY + // in the fitting procedure !!! + // + // This has to do with the fact that the internal function histogramming + // uses 100 bins and does not adapt to the binning of the fitted histogram, unlike PAW does + // (very important if you use Sumw2, see e.g. ROOTTALK: Mon May 26 1997 - 09:56:03 MEST) + // + // So, WE have to adapt to that internal flaw of ROOT: + // + const Int_t npx = fSinglePheFit->GetNpx(); + const Int_t bins = fHBPQ->GetXaxis()->GetLast()-fHBPQ->GetXaxis()->GetFirst(); + // fHBPQ->Scale(gkSq2Pi*(float)bins/npx/entries); + + // + // we need this, otherwise, ROOT does not calculate the area correctly + // don't ask me why it does not behave correctly, it's one of the nasty + // mysteries of ROOT which takes you a whole day to find out :-) + // + // fSinglePheFit->SetNpx(fQnbins); + + fHBPQ->Fit("SinglePheFit",opt); + + fLambda = fSinglePheFit->GetParameter(0); + fMu0 = fSinglePheFit->GetParameter(1); + fMu1 = fSinglePheFit->GetParameter(2); + fSigma0 = fSinglePheFit->GetParameter(3); + fSigma1 = fSinglePheFit->GetParameter(4); + + fLambdaErr = fSinglePheFit->GetParError(0); + fMu0Err = fSinglePheFit->GetParError(1); + fMu1Err = fSinglePheFit->GetParError(2); + fSigma0Err = fSinglePheFit->GetParError(3); + fSigma1Err = fSinglePheFit->GetParError(4); + + fProb = fSinglePheFit->GetProb(); + fChisquare = fSinglePheFit->GetChisquare(); + fNdf = fSinglePheFit->GetNDF(); + + *fLog << "Results of the Blind Pixel Fit: " << endl; + *fLog << "Chisquare: " << fChisquare << endl; + *fLog << "DoF: " << fNdf << endl; + *fLog << "Probability: " << fProb << endl; + *fLog << "Integral: " << fSinglePheFit->Integral(rmin,rmax); + + // + // The fit result is accepted under condition + // The Probability is greater than gkProbLimit (default 0.01 == 99%) + // + if (fProb < gkProbLimit) + { + *fLog << warn << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl; + fFitOK = kFALSE; + return kFALSE; + } + + + fFitOK = kTRUE; + + return kTRUE; +} + + +void MHCalibrationBlindPixel::CutAllEdges() +{ + + // + // The number 100 is necessary because it is the internal binning + // of ROOT functions. A call to SetNpx() does NOT help + // If you find another solution which WORKS!!, please tell me!! + // + Int_t nbins = 100; + + *fLog << "New number of bins in HSinQ: " << CutEdges(fHBPQ,nbins) << endl; + + fBPQfirst = fHBPQ->GetBinLowEdge(fHBPQ->GetXaxis()->GetFirst()); + fBPQlast = fHBPQ->GetBinLowEdge(fHBPQ->GetXaxis()->GetLast())+fHBPQ->GetBinWidth(0); + fBPQnbins = nbins; + + *fLog << "New number of bins in HErrQ: " << CutEdges(fHBPErrQ,30) << endl; + fErrBPQfirst = fHBPErrQ->GetBinLowEdge(fHBPErrQ->GetXaxis()->GetFirst()); + fErrBPQlast = fHBPErrQ->GetBinLowEdge(fHBPErrQ->GetXaxis()->GetLast())+fHBPErrQ->GetBinWidth(0); + fErrBPQnbins = nbins; + + CutEdges(fHBPQvsN,0); + +} + +Bool_t MHCalibrationBlindPixel::FitT(Axis_t rmin, Axis_t rmax, Option_t *opt) +{ + + if (fTimeGausFit) + return kFALSE; + + rmin = (rmin != 0.) ? rmin : 0.; + rmax = (rmax != 0.) ? rmax : 16.; + + const Stat_t entries = fHBPT->GetEntries(); + const Double_t mu_guess = fHBPT->GetBinCenter(fHBPT->GetMaximumBin()); + const Double_t sigma_guess = (rmax - rmin)/2.; + const Double_t area_guess = entries/gkSq2Pi; + + fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax); + fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess); + fTimeGausFit->SetParNames("Area","#mu","#sigma"); + fTimeGausFit->SetParLimits(0,0.,entries); + fTimeGausFit->SetParLimits(1,rmin,rmax); + fTimeGausFit->SetParLimits(2,0.,rmax-rmin); + + fHBPT->Fit("GausTime",opt); + + fMeanT = fTimeGausFit->GetParameter(2); + fSigmaT = fTimeGausFit->GetParameter(3); + fMeanTErr = fTimeGausFit->GetParError(2); + fSigmaTErr = fTimeGausFit->GetParError(3); + + Float_t prob = fTimeGausFit->GetProb(); + + *fLog << "Results of the Times Fit: " << endl; + *fLog << "Chisquare: " << fTimeGausFit->GetChisquare() << endl; + *fLog << "Ndf: " << fTimeGausFit->GetNDF() << endl; + *fLog << "Probability: " << prob << endl; + + if (prob < gkProbLimit) + { + *fLog << warn << "Fit of the Arrival times failed ! " << endl; + return kFALSE; + } + + return kTRUE; + +} Index: /trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.h (revision 2525) +++ /trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.h (revision 2525) @@ -0,0 +1,132 @@ +#ifndef MARS_MHCalibrationBlindPixel +#define MARS_MHCalibrationBlindPixel + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef MARS_MHCalibrationConfig +#include "MHCalibrationConfig.h" +#endif + +#ifndef ROOT_TH1 +#include "TH1.h" +#endif + +#ifndef ROOT_TH1F +#include "TH1F.h" +#endif + +#ifndef ROOT_TF1 +#include "TF1.h" +#endif + +#ifndef ROOT_TPaveText +#include "TPaveText.h" +#endif + +class TMath; +class MParList; +class MHCalibrationBlindPixel : public MH +{ +private: + + TH1I* fHBPQ; //-> Histogram with the single Phe spectrum + TH1F* fHBPErrQ; //-> Variance of summed FADC slices + TH1I* fHBPT; //-> Variance of summed FADC slices + TH1I* fHBPQvsN; //-> Summed Charge vs. Event Nr. + + TF1 *fSinglePheFit; + TF1 *fTimeGausFit; + + Axis_t fBPQfirst; + Axis_t fBPQlast; + Int_t fBPQnbins; + + Axis_t fErrBPQfirst; + Axis_t fErrBPQlast; + Int_t fErrBPQnbins; + + void ResetBin(Int_t i); + void DrawLegend(); + + TPaveText *fFitLegend; + Bool_t fFitOK; + + BPFitFunc fgSinglePheFitFunc; // In the beginning, + Int_t fgSinglePheFitNPar; // we want to be flexible using different functions + + Double_t fLambda; + Double_t fMu0; + Double_t fMu1; + Double_t fSigma0; + Double_t fSigma1; + + Double_t fLambdaErr; + Double_t fMu0Err; + Double_t fMu1Err; + Double_t fSigma0Err; + Double_t fSigma1Err; + + Double_t fChisquare; + Double_t fProb; + Int_t fNdf; + + Double_t fMeanT; + Double_t fMeanTErr; + Double_t fSigmaT; + Double_t fSigmaTErr; + +public: + + MHCalibrationBlindPixel(const char *name=NULL, const char *title=NULL); + ~MHCalibrationBlindPixel(); + + Bool_t FillBPQ(Int_t q) { return fHBPQ->Fill(q) > -1; } + Bool_t FillErrBPQ(Float_t errq) { return fHBPErrQ->Fill(errq) > -1; } + Bool_t FillBPT(Int_t t) { return fHBPT->Fill(t) > -1; } + Bool_t FillBPQvsN(Stat_t rq, Int_t t) { return fHBPQvsN->Fill(t,rq) > -1; } + + const Double_t GetLambda() const { return fLambda; } + const Double_t GetMu0() const { return fMu0; } + const Double_t GetMu1() const { return fMu1; } + const Double_t GetSigma0() const { return fSigma0; } + const Double_t GetSigma1() const { return fSigma1; } + + const Double_t GetLambdaErr() const { return fLambdaErr; } + const Double_t GetMu0Err() const { return fMu0Err; } + const Double_t GetMu1Err() const { return fMu1Err; } + const Double_t GetSigma0Err() const { return fSigma0Err; } + const Double_t GetSigma1Err() const { return fSigma1Err; } + + const Double_t GetChiSquare() const { return fChisquare; } + const Double_t GetProb() const { return fProb; } + const Int_t GetNdf() const { return fNdf; } + + const Double_t GetMeanT() const { return fMeanT; } + const Double_t GetMeanTErr() const { return fMeanTErr; } + const Double_t GetSigmaT() const { return fSigmaT; } + const Double_t GetSigmaTErr() const { return fSigmaTErr; } + + const TH1F *GetHErrQ() { return fHBPErrQ; } + const TH1F *GetHErrQ() const { return fHBPErrQ; } + + Bool_t SimulateSinglePhe(Double_t lambda, + Double_t mu0, + Double_t mu1, + Double_t sigma0, + Double_t sigma1); + + Bool_t FitSinglePhe(Axis_t rmin=0, Axis_t rmax=0, Option_t *opt="R0+"); + Bool_t FitT(Axis_t rmin=0., Axis_t rmax=0.,Option_t *opt="R0+"); + + void ChangeFitFunc(BPFitFunc fitfunc, Int_t par=5); + + + void CutAllEdges(); + void Draw(Option_t *option=""); + + ClassDef(MHCalibrationBlindPixel, 0) +}; + +#endif /* MARS_MHCalibrationBlindPixel */ Index: /trunk/MagicSoft/Mars/mhist/MHCalibrationConfig.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCalibrationConfig.h (revision 2525) +++ /trunk/MagicSoft/Mars/mhist/MHCalibrationConfig.h (revision 2525) @@ -0,0 +1,34 @@ +#ifndef MARS_MHCalibrationConfig +#define MARS_MHCalibrationConfig + +///////////////////////////////////////////////////////////////////////////// +// // +// MHCalibrationConfig // +// // +// Contains all configuration data of the Calibration // +// // +///////////////////////////////////////////////////////////////////////////// + + +// Global rejection criteria for the acceptance of a fit: Prob=0.01 == 99% Probability +const Float_t gkProbLimit = 0.01; + +// Starting number of bins for the histo: +const Int_t gkStartBlindPixelBinNr = 10000; + +// Starting number of bins for the histo: (maximum by hardware: 4080) +const Int_t gkStartPINDiodeBinNr = 4000; + +// Starting number of bins for the histo: (maximum possible by hardware = 40800) +const Int_t gkStartPixelBinNr = 20000; + +// Starting number for the highest value of the Q-histo: +const Axis_t gkStartQlast = 5000.; + +// Square root of 2 pi: +const Float_t gkSq2Pi = 2.506628274631; + +// typedef to the fitting functions for the blind pixel +typedef Double_t (*BPFitFunc)(Double_t *, Double_t *); + +#endif /* MARS_MHCalibrationBlindPixelConfig */ Index: /trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.cc (revision 2525) +++ /trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.cc (revision 2525) @@ -0,0 +1,110 @@ +/* ======================================================================== *\ +! +! * +! * 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-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHCalibrationPINDiode // +// // +// Performs all the necessary fits to extract the mean number of photons // +// out of the derived light flux // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MHCalibrationPINDiode.h" +#include "MHCalibrationConfig.h" + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHCalibrationPINDiode); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHCalibrationPINDiode::MHCalibrationPINDiode(const char *name, const char *title) + : fVarGausFit(NULL) +{ + + fName = name ? name : "MHCalibrationPINDiode"; + fTitle = title ? title : "Fill the accumulated charges and times all PINDiode events and perform fits"; + + // Create a large number of bins, later we will rebin + fQfirst = 0; + fQlast = gkStartPINDiodeBinNr; + fQnbins = gkStartPINDiodeBinNr; + + fHPQ = new TH1I("HPQ","Distribution of Summed FADC Slices",fQnbins,fQfirst,fQlast); + fHPQ->SetXTitle("Sum FADC Slices"); + fHPQ->SetYTitle("Nr. of events"); + fHPQ->Sumw2(); + + fErrQfirst = 0.; + fErrQlast = gkStartPINDiodeBinNr; + fErrQnbins = gkStartPINDiodeBinNr; + + fHErrQ = new TH1F("HErrQ","Distribution of Variances of Summed FADC Slices",fErrQnbins,fErrQfirst,fErrQlast); + fHErrQ->SetXTitle("Variance Summed FADC Slices"); + fHErrQ->SetYTitle("Nr. of events"); + fHErrQ->Sumw2(); + + Int_t tfirst = 0; + Int_t tlast = 31; + Int_t nbins = 32; + + fHPT = new TH1I("HPT","Distribution of Mean Arrival Times",nbins,tfirst,tlast); + fHPT->SetXTitle("Mean Arrival Times [FADC slice nr]"); + fHPT->SetYTitle("Nr. of events"); + fHPT->Sumw2(); + +} + +MHCalibrationPINDiode::~MHCalibrationPINDiode() +{ + + delete fHPQ; + delete fHErrQ; + + if (fVarGausFit) + delete fVarGausFit; + + delete fHPT; + +} + Index: /trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.h (revision 2525) +++ /trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.h (revision 2525) @@ -0,0 +1,51 @@ +#ifndef MARS_MHCalibrationPINDiode +#define MARS_MHCalibrationPINDiode + +#ifndef MARS_MHCalibrationPixel +#include "MHCalibrationPixel.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH1 +#include "TH1.h" +#endif + +#ifndef ROOT_TH1F +#include "TH1F.h" +#endif + +#ifndef ROOT_TF1 +#include "TF1.h" +#endif + +class MHCalibrationPINDiode : public MHCalibrationPixel +{ +private: + + TH1I* fHPQ; //-> Histogram containing the summed 32 PINDiode slices + TH1F* fHErrQ; //-> Variance of summed FADC slices + TH1I* fHPT; //-> Histogram with time evolution of summed charges + + + TF1 *fVarGausFit; + + Float_t fErrQfirst; + Float_t fErrQlast; + UShort_t fErrQnbins; + +public: + + MHCalibrationPINDiode(const char *name=NULL, const char *title=NULL); + ~MHCalibrationPINDiode(); + + const Double_t GetT() const { return fVarGausFit->GetParameter(2); } + const Double_t GetErrT() const { return fVarGausFit->GetParameter(3); } + + ClassDef(MHCalibrationPINDiode, 0) +}; + +#endif + Index: /trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.cc (revision 2525) +++ /trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.cc (revision 2525) @@ -0,0 +1,433 @@ +/* ======================================================================== *\ +! +! * +! * 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-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHCalibrationPixel // +// // +// Performs all the necessary fits to extract the mean number of photons // +// out of the derived light flux // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHCalibrationPixel.h" +#include "MHCalibrationConfig.h" +#include "MCalibrationFits.h" + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHCalibrationPixel); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHCalibrationPixel::MHCalibrationPixel(Int_t pix, const char *name, const char *title) + : fFitOK(kFALSE), fPixId(pix), fTGausFit(NULL), fQGausFit(NULL), fFitLegend(NULL) +{ + + fName = name ? name : "MHCalibrationPixel"; + fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits"; + + TString qname = "HQ"; + TString qtitle = "Distribution of Summed FADC Slices Pixel "; + qname += pix; + qtitle += pix; + + // Create a large number of bins, later we will rebin + fQfirst = -0.5; + fQlast = gkStartQlast - 0.5; + fQnbins = gkStartPixelBinNr; + + fHQ = new TH1I( qname.Data(),qtitle.Data(), + fQnbins,fQfirst,fQlast); + fHQ->SetXTitle("Sum FADC Slices"); + fHQ->SetYTitle("Nr. of events"); + fHQ->Sumw2(); + + TString tname = "HT"; + TString ttitle = "Distribution of Mean Arrival Times Pixel "; + tname += pix; + ttitle += pix; + + Axis_t tfirst = -0.5; + Axis_t tlast = 15.5; + Int_t nbins = 16; + + fHT = new TH1I(tname.Data(),ttitle.Data(), + nbins,tfirst,tlast); + fHT->SetXTitle("Mean Arrival Times [FADC slice nr]"); + fHT->SetYTitle("Nr. of events"); + fHT->Sumw2(); + + TString qvsnname = "HQvsN"; + TString qvsntitle = "Sum of Charges vs. Event Number Pixel "; + qvsnname += pix; + qvsntitle += pix; + + // We define a reasonable number and later enlarge it if necessary + nbins = 20000; + Axis_t nfirst = -0.5; + Axis_t nlast = (Axis_t)nbins - 0.5; + + fHQvsN = new TH1I(qvsnname.Data(),qvsntitle.Data(), + nbins,nfirst,nlast); + fHQvsN->SetXTitle("Event Nr."); + fHQvsN->SetYTitle("Sum of FADC slices"); + + fQChisquare = -1.; + fQProb = -1.; + fQNdf = -1; + + fTChisquare = -1.; + fTProb = -1.; + fTNdf = -1; + +} + +MHCalibrationPixel::~MHCalibrationPixel() +{ + + delete fHQ; + delete fHT; + delete fHQvsN; + + if (fQGausFit) + delete fQGausFit; + if (fTGausFit) + delete fTGausFit; + if (fFitLegend) + delete fFitLegend; + +} + +void MHCalibrationPixel::Reset() +{ + + for (Int_t i = fHQ->FindBin(fQfirst); i <= fHQ->FindBin(fQlast); i++) + fHQ->SetBinContent(i, 1.e-20); + + for (Int_t i = 0; i < 16; i++) + fHT->SetBinContent(i, 1.e-20); + + fQlast = gkStartQlast; + + fHQ->GetXaxis()->SetRangeUser(0.,fQlast); + + return; +} + + +// ------------------------------------------------------------------------- +// +// Set the binnings and prepare the filling of the histograms +// +Bool_t MHCalibrationPixel::SetupFill(const MParList *plist) +{ + + fHQ->Reset(); + fHT->Reset(); + + return kTRUE; +} + + + +// ------------------------------------------------------------------------- +// +// Draw a legend with the fit results +// +void MHCalibrationPixel::DrawLegend() +{ + + fFitLegend = new TPaveText(0.05,0.05,0.95,0.95); + + if (fFitOK) + fFitLegend->SetFillColor(80); + else + fFitLegend->SetFillColor(2); + + fFitLegend->SetLabel("Results of the Gauss Fit:"); + fFitLegend->SetTextSize(0.05); + + char line1[32]; + sprintf(line1,"Mean: Q_{#mu} = %2.2f #pm %2.2f",fQMean,fQMeanErr); + fFitLegend->AddText(line1); + + char line4[32]; + sprintf(line4,"Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fQSigma,fQSigmaErr); + fFitLegend->AddText(line4); + + char line7[32]; + sprintf(line7,"#chi^{2} / N_{dof}: %4.2f / %3i",fQChisquare,fQNdf); + fFitLegend->AddText(line7); + + char line8[32]; + sprintf(line8,"Probability: %4.2f ",fQProb); + fFitLegend->AddText(line8); + + if (fFitOK) + fFitLegend->AddText("Result of the Fit: OK"); + else + fFitLegend->AddText("Result of the Fit: NOT OK"); + + fFitLegend->SetBit(kCanDelete); + fFitLegend->Draw(); + +} + + +// ------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHCalibrationPixel::Draw(Option_t *opt) +{ + + gStyle->SetOptFit(0); + gStyle->SetOptStat(1111111); + + TCanvas *c = MakeDefCanvas(this,600,900); + + gROOT->SetSelectedPad(NULL); + + c->Divide(2,2); + + c->cd(1); + gPad->SetBorderMode(0); + gPad->SetLogy(1); + gPad->SetTicks(); + + fHQ->DrawCopy(opt); + + if (fQGausFit) + { + if (fFitOK) + fQGausFit->SetLineColor(kGreen); + else + fQGausFit->SetLineColor(kRed); + + fQGausFit->DrawCopy("same"); + c->Modified(); + c->Update(); + } + + + c->cd(2); + DrawLegend(); + c->Update(); + + c->cd(3); + gStyle->SetOptStat(1111111); + + gPad->SetLogy(1); + fHT->DrawCopy(opt); + + if (fHT->GetFunction("GausTime")) + { + TF1 *tfit = fHT->GetFunction("GausTime"); + if (tfit->GetProb() < 0.01) + tfit->SetLineColor(kRed); + else + tfit->SetLineColor(kGreen); + + tfit->DrawCopy("same"); + c->Modified(); + c->Update(); + } + + c->Modified(); + c->Update(); + + c->cd(4); + fHQvsN->DrawCopy(opt); +} + + + +Bool_t MHCalibrationPixel::FitT(Axis_t rmin, Axis_t rmax, Option_t *option) +{ + + if (fTGausFit) + return kFALSE; + + rmin = (rmin != 0.) ? rmin : -0.5; + rmax = (rmax != 0.) ? rmax : 15.5; + + const Stat_t entries = fHT->GetEntries(); + const Double_t mu_guess = fHT->GetBinCenter(fHT->GetMaximumBin()); + const Double_t sigma_guess = (rmax - rmin)/2.; + const Double_t area_guess = entries/gkSq2Pi; + + fTGausFit = new TF1("GausTime","gaus",rmin,rmax); + + if (!fTGausFit) + { + *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl; + return kFALSE; + } + + fTGausFit->SetParameters(area_guess,mu_guess,sigma_guess); + fTGausFit->SetParNames("Area","#mu","#sigma"); + fTGausFit->SetParLimits(0,0.,entries); + fTGausFit->SetParLimits(1,rmin,rmax); + fTGausFit->SetParLimits(2,0.,rmax-rmin); + + fHT->Fit("GausTime",option); + + fTChisquare = fTGausFit->GetChisquare(); + fTNdf = fTGausFit->GetNDF(); + fTProb = fTGausFit->GetProb(); + fTMean = fTGausFit->GetParameter(1); + fTSigma = fTGausFit->GetParameter(2); + + if (fTProb < gkProbLimit) + { + *fLog << warn << "Fit of the Arrival times failed ! " << endl; + return kFALSE; + } + + return kTRUE; + +} + +Bool_t MHCalibrationPixel::FitQ(Axis_t rmin, Axis_t rmax, Option_t *option) +{ + + if (fQGausFit) + return kFALSE; + + // + // Get the fitting ranges + // + rmin = (rmin != 0.) ? rmin : fQfirst; + rmax = (rmax != 0.) ? rmax : fQlast; + + // + // First guesses for the fit (should be as close to reality as possible, + // otherwise the fit goes gaga because of high number of dimensions ... + // + const Stat_t entries = fHQ->GetEntries(); + const Double_t ar_guess = entries/gkSq2Pi; + const Double_t mu_guess = fHQ->GetBinCenter(fHQ->GetMaximumBin()); + const Double_t si_guess = mu_guess/500.; + + fQGausFit = new TF1("QGausFit","gaus",rmin,rmax); + + if (!fQGausFit) + { + *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl; + return kFALSE; + } + + fQGausFit->SetParameters(ar_guess,mu_guess,si_guess); + fQGausFit->SetParNames("Area","#mu","#sigma"); + fQGausFit->SetParLimits(0,0.,entries); + fQGausFit->SetParLimits(1,rmin,rmax); + fQGausFit->SetParLimits(2,0.,rmax-rmin); + + fHQ->Fit("QGausFit",option); + + fQChisquare = fQGausFit->GetChisquare(); + fQNdf = fQGausFit->GetNDF(); + fQProb = fQGausFit->GetProb(); + fQMean = fQGausFit->GetParameter(1); + fQMeanErr = fQGausFit->GetParError(1); + fQSigma = fQGausFit->GetParameter(2); + fQSigmaErr = fQGausFit->GetParError(2); + + // + // The fit result is accepted under condition + // The Probability is greater than gkProbLimit (default 0.01 == 99%) + // + if (fQProb < gkProbLimit) + { + *fLog << warn << "Prob: " << fQProb << " is smaller than the allowed value: " << gkProbLimit << endl; + fFitOK = kFALSE; + return kFALSE; + } + + + fFitOK = kTRUE; + + return kTRUE; +} + + +void MHCalibrationPixel::CutAllEdges() +{ + + // + // The number 100 is necessary because it is the internal binning + // of ROOT functions. A call to SetNpx() does NOT help + // If you find another solution which WORKS!!, please tell me!! + // + Int_t nbins = 100; + + // *fLog << "New number of bins in HQ: " << CutEdges(fHQ,nbins) << endl; + + fQfirst = fHQ->GetBinLowEdge(fHQ->GetXaxis()->GetFirst()); + fQlast = fHQ->GetBinLowEdge(fHQ->GetXaxis()->GetLast())+fHQ->GetBinWidth(0); + fQnbins = nbins; + + CutEdges(fHQvsN,0); + +} + +void MHCalibrationPixel::PrintQFitResult() +{ + + *fLog << "Results of the Summed Charges Fit: " << endl; + *fLog << "Chisquare: " << fQChisquare << endl; + *fLog << "DoF: " << fQNdf << endl; + *fLog << "Probability: " << fQProb << endl; + *fLog << endl; + +} + +void MHCalibrationPixel::PrintTFitResult() +{ + + *fLog << "Results of the Arrival Time Slices Fit: " << endl; + *fLog << "Chisquare: " << fTChisquare << endl; + *fLog << "Ndf: " << fTNdf << endl; + *fLog << "Probability: " << fTProb << endl; + *fLog << endl; + +} Index: /trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.h (revision 2525) +++ /trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.h (revision 2525) @@ -0,0 +1,117 @@ +#ifndef MARS_MHCalibrationPixel +#define MARS_MHCalibrationPixel + +#ifndef ROOT_TH1 +#include "TH1.h" +#endif + +#ifndef ROOT_TH1 +#include "TH1F.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TF1 +#include "TF1.h" +#endif + +class TPaveText; +class TMath; +class MParList; + +class MHCalibrationPixel : public MH +{ + +private: + + Int_t fPixId; //-> Pixel Nr + +protected: + + TH1I* fHQ; //-> Summed FADC slices + TH1I* fHT; //-> Mean arrival time in number of FADC sice + + TH1I* fHQvsN; //-> Summed Charge vs. Event Nr. + + TF1* fQGausFit; + TF1* fTGausFit; + + TPaveText *fFitLegend; + + Axis_t fQfirst; + Axis_t fQlast; + Int_t fQnbins; + + Bool_t fFitOK; + + Double_t fQChisquare; + Double_t fQProb; + Int_t fQNdf; + + Double_t fQMean; + Double_t fQMeanErr; + Double_t fQSigma; + Double_t fQSigmaErr; + + Double_t fTChisquare; + Double_t fTProb; + Int_t fTNdf; + + Double_t fTMean; + Double_t fTSigma; + + virtual void DrawLegend(); + +public: + + MHCalibrationPixel(Int_t pix=-1, const char *name=NULL, const char *title=NULL); + ~MHCalibrationPixel(); + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *, const Stat_t w=1) { return kTRUE; } + + Bool_t FillQ(Int_t q) { return fHQ->Fill(q) > -1; } + Bool_t FillT(Int_t t) { return fHT->Fill(t) > -1; } + Bool_t FillQvsN(Float_t q, Int_t n) { return fHQvsN->Fill(n,q) > -1; } + + const TH1I *GetHQ() { return fHQ; } + const TH1I *GetHQ() const { return fHQ; } + + const Double_t GetQMean() const { return fQMean; } + const Double_t GetQMeanErr() const { return fQMeanErr; } + const Double_t GetQSigma() const { return fQSigma; } + const Double_t GetQSigmaErr() const { return fQSigmaErr; } + const Double_t GetArea() const { return fQGausFit->GetParameter(0); } + const Double_t GetAreaErr() const { return fQGausFit->GetParError(0); } + + const Double_t GetQChiSquare() const { return fQChisquare; } + const Double_t GetQProb() const { return fQProb; } + const Int_t GetQNdf() const { return fQNdf; } + + const Double_t GetTMean() const { return fTMean; } + const Double_t GetTSigma() const { return fTSigma; } + + const TH1I *GetHT() { return fHT; } + const TH1I *GetHT() const { return fHT; } + + const TH1I *GetHQvsN() { return fHQvsN; } + const TH1I *GetHQvsN() const { return fHQvsN; } + + Bool_t FitQ(Axis_t rmin=0, Axis_t rmax=0, Option_t *option="RQ0+"); + Bool_t FitT(Axis_t rmin=0, Axis_t rmax=0, Option_t *option="RQ0+"); + + virtual void Draw(Option_t *option=""); + virtual void CutAllEdges(); + virtual void Reset(); + + void PrintQFitResult(); + void PrintTFitResult(); + + Bool_t IsFitted() { return fFitOK; } + + ClassDef(MHCalibrationPixel, 1) +}; + +#endif