Index: trunk/MagicSoft/Mars/mhist/MH.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MH.cc (revision 2522) +++ 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 2522) +++ 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