/* ======================================================================== *\ ! ! * ! * 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): Thomas Bretz 1/2002 ! Author(s): Wolfgang Wittek 1/2002 ! ! Copyright: MAGIC Software Development, 2000-2002 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // MHEffOnTimeTheta // // // // calculates the effective on time for each bin in Theta // // // ////////////////////////////////////////////////////////////////////////////// #include "MHEffOnTimeTheta.h" #include #include #include #include #include "MTime.h" #include "MBinning.h" #include "MParList.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHEffOnTimeTheta); // -------------------------------------------------------------------------- // // Default Constructor. It sets name and title of the histograms. // MHEffOnTimeTheta::MHEffOnTimeTheta(const char *name, const char *title) : fHEffOn() { // // set the name and title of this object // fName = name ? name : "MHEffOnTimeTheta"; fTitle = title ? title : "1-D histogram of Eff On Time"; // effective on time versus theta fHEffOn.SetName("EffOn"); fHEffOn.SetTitle("Effective On Time vs. Theta"); fHEffOn.SetDirectory(NULL); fHEffOn.SetXTitle("\\Theta [\\circ]"); fHEffOn.SetYTitle("t_{eff} [s]"); // chi2/NDF versus theta fHChi2.SetName("Chi2/NDF"); fHChi2.SetTitle("Chi2/NDF of OnTimeFit vs. Theta"); fHChi2.SetDirectory(NULL); fHChi2.SetXTitle("\\Theta [\\circ]"); fHChi2.SetYTitle("chi^{2}/NDF"); // lambda versus theta fHLambda.SetName("lambda"); fHLambda.SetTitle("lambda of OnTimeFit vs. Theta"); fHLambda.SetDirectory(NULL); fHLambda.SetXTitle("\\Theta [\\circ]"); fHLambda.SetYTitle("\\lambda [Hz]"); // N0del versus theta fHN0del.SetName("N0del"); fHN0del.SetTitle("N0del of OnTimeFit vs. Theta"); fHN0del.SetDirectory(NULL); fHN0del.SetXTitle("\\Theta [\\circ]"); fHN0del.SetYTitle("N0del"); } // ----------------------------------------------------------------------- // // Calculate the effective on time by fitting the distribution of // time differences // void MHEffOnTimeTheta::Calc(TH2D *hist) { // nbins = number of Theta bins const Int_t nbins = hist->GetNbinsY(); for (int i=1; i<=nbins; i++) { // TH1D &h = *hist->ProjectionX("Calc-theta", i, i); TH1D &h = *hist->ProjectionX("Calc-theta", i, i, "E"); // Nmdel = Nm * binwidth, with Nm = number of observed events const Double_t Nmdel = h.Integral("width"); const Double_t Nm = h.Integral(); // Double_t mean = h->GetMean(); //................................................... // determine range (yq[0], yq[1]) of time differences // where fit should be performed; // require a fraction >=xq[0] of all entries to lie below yq[0] // and a fraction <=xq[1] of all entries to lie below yq[1]; // within the range (yq[0], yq[1]) there must be no empty bin; // choose pedestrian approach as long as GetQuantiles is not available Double_t xq[2] = { 0.15, 0.95 }; Double_t yq[2]; // GetQuantiles doesn't seem to be available in root 3.01/06 // h->GetQuantiles(2,yq,xq); const Double_t sumtot = h.Integral(); const Int_t jbins = h.GetNbinsX(); if (sumtot > 0.0) { Double_t sum1 = 0.0; yq[0] = h.GetBinLowEdge(jbins+1); for (int j=1; j<=jbins; j++) { if (sum1 >= xq[0]*sumtot) { yq[0] = h.GetBinLowEdge(j); break; } sum1 += h.GetBinContent(j); } Double_t sum2 = 0.0; yq[1] = h.GetBinLowEdge(jbins+1); for (int j=1; j<=jbins; j++) { Double_t content = h.GetBinContent(j); sum2 += content; if (sum2 >= xq[1]*sumtot || content == 0.0) { yq[1] = h.GetBinLowEdge(j); break; } } //................................................... // parameter 0 = lambda // parameter 1 = N0*del with N0 = ideal number of events // and del = bin width of time difference TF1 func("Poisson", "[1] * [0] * exp(-[0] *x)", yq[0], yq[1]); func.SetParameter(0, 100); // [Hz] func.SetParameter(1, Nmdel); func.SetParLimits(0, 0, 1000); // [Hz] func.SetParLimits(1, 0, 10*Nmdel); func.SetParName(0, "lambda"); func.SetParName(1, "Nmdel"); // options : 0 do not plot the function // I use integral of function in bin rather than value at bin center // R use the range specified in the function range // Q quiet mode h.Fit("Poisson", "0IRQ"); // gPad->SetLogy(); // gStyle->SetOptFit(1011); // h->GetXaxis()->SetTitle("time difference [s]"); // h->GetYaxis()->SetTitle("Counts"); // h->DrawCopy(); // func.SetRange(yq[0], yq[1]); // Range of Drawing // func.DrawCopy("same"); const Double_t lambda = func.GetParameter(0); const Double_t N0del = func.GetParameter(1); const Double_t chi2 = func.GetChisquare(); const Int_t NDF = func.GetNDF(); // was fit successful ? if (NDF>0 && chi2<2.5*NDF) { // the effective on time is Nm/lambda fHEffOn.SetBinContent(i, Nm/lambda); // plot chi2/NDF of fit fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0); // lambda of fit fHLambda.SetBinContent(i, lambda); // N0del of fit fHN0del.SetBinContent(i, N0del); delete &h; continue; } } fHEffOn.SetBinContent (i, 1.e-20); fHChi2.SetBinContent (i, 1.e-20); fHLambda.SetBinContent(i, 1.e-20); fHN0del.SetBinContent (i, 1.e-20); delete &h; } } // ------------------------------------------------------------------------- // // Set the binnings and prepare the filling of the histograms // Bool_t MHEffOnTimeTheta::SetupFill(const MParList *plist) { const MBinning* bins = (MBinning*)plist->FindObject("BinningTheta"); if (!bins) { *fLog << err << dbginf << "BinningTheta [MBinning] not found... aborting." << endl; return kFALSE; } SetBinning(&fHEffOn, bins); SetBinning(&fHChi2, bins); SetBinning(&fHLambda, bins); SetBinning(&fHN0del, bins); fHEffOn.Sumw2(); fHChi2.Sumw2(); fHLambda.Sumw2(); fHN0del.Sumw2(); return kTRUE; } // ------------------------------------------------------------------------- // // Draw the histogram // void MHEffOnTimeTheta::Draw(Option_t *opt) { TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); pad->SetBorderMode(0); pad->Divide(2,2); pad->cd(1); gPad->SetBorderMode(0); fHEffOn.Draw(opt); pad->cd(2); gPad->SetBorderMode(0); fHChi2.Draw(opt); pad->cd(3); gPad->SetBorderMode(0); fHLambda.Draw(opt); pad->cd(4); gPad->SetBorderMode(0); fHN0del.Draw(opt); pad->Modified(); pad->Update(); }