Index: /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 5430) +++ /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 5431) @@ -51,4 +51,5 @@ #pragma link C++ class MPad+; #pragma link C++ class MMcCalibrationUpdate+; +#pragma link C++ class MMakePadHistograms+; #endif Index: /trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.cc (revision 5431) +++ /trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.cc (revision 5431) @@ -0,0 +1,350 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek, 09/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MMakePadHistograms +// +// reads calibrated data of a particular type (ON, OFF or MC) +// and produces the histograms for the padding +// +// these histograms are further treated in MPad +// +///////////////////////////////////////////////////////////////////////////// +#include "MMakePadHistograms.h" + +#include +#include + +#include +#include +#include +#include + +#include "MBadPixelsCalc.h" +#include "MBadPixelsTreat.h" +#include "MContinue.h" +#include "MEvtLoop.h" +#include "MFillH.h" +#include "MFSelBasic.h" +#include "MGeomApply.h" +#include "MHBadPixels.h" +#include "MHSigmaTheta.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MProgressBar.h" +#include "MReadMarsFile.h" +#include "MReadReports.h" +#include "MTaskList.h" +//#include "MExtrapolatePointingPos.h" +#include "MPointingPosCalc.h" + +ClassImp(MMakePadHistograms); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MMakePadHistograms::MMakePadHistograms(const char *name, const char *title) +{ + fName = name ? name : "MMakePadHistograms"; + fTitle = title ? title : "Class to make Pad histograms"; + + fHSigmaTheta = NULL; + fHSigmaThetaOuter = NULL; + fHDiffPixTheta = NULL; + fHSigmaPixTheta = NULL; + + fHBadIdTheta = NULL; + fHBadNTheta = NULL; + + // default values + fType = ""; + fMaxEvents = -1; + fNamePedPhotCam = "MPedPhotCamFromData"; + fNameInputFile = ""; + fNameOutputFile = ""; + fPedestalLevel = 3.0; + fUseInterpolation = kTRUE; + fProcessPedestal = kTRUE; + fProcessTime = kFALSE; +} + + +// -------------------------------------------------------------------------- +// +// SetDataType. +// +// the type may be ON, OFF or MC +// +void MMakePadHistograms::SetDataType(const char *type) +{ + fType = type; +} + +// -------------------------------------------------------------------------- +// +// SetMaxEvents. +// +// set the maximum number of events to be read +// +void MMakePadHistograms::SetMaxEvents(Int_t maxev) +{ + fMaxEvents = maxev; +} + +// -------------------------------------------------------------------------- +// +// SetNameInputFile. +// +// the input file contains the calibrated data +// +void MMakePadHistograms::SetNameInputFile(const char *input) +{ + fNameInputFile = input; +} + +// -------------------------------------------------------------------------- +// +// SetNameOutputFile. +// +// the output file contains the histgrams for the padding +// +void MMakePadHistograms::SetNameOutputFile(const char *output) +{ + fNameOutputFile = output; +} + +// -------------------------------------------------------------------------- +// +// SetNamePedPhotCam. +// +// set the name of the MPedPhotCam container; +// it is used in : MBadPixelsCalc, MBadPixelsTreat, MHSigmaTheta +// +void MMakePadHistograms::SetNamePedPhotCam(const char *name) +{ + fNamePedPhotCam = name; +} + +// -------------------------------------------------------------------------- +// +// SetPedestalLevel. +// +// set the pedestal level; +// it is used in : MBadPixelsCalc +// +void MMakePadHistograms::SetPedestalLevel(Float_t pedlevel) +{ + fPedestalLevel = pedlevel; +} + +// -------------------------------------------------------------------------- +// +// SetUseInterpolation. +// +// set the option "UseInterpolation"; +// it is used in : MBadPixelsTreat +// +void MMakePadHistograms::SetUseInterpolation(Bool_t useinter) +{ + fUseInterpolation = useinter; +} + +// -------------------------------------------------------------------------- +// +// SetProcessPedestal. +// +// set the option "ProcessPedestal"; +// it is used in : MBadPixelsTreat +// +void MMakePadHistograms::SetProcessPedestal(Bool_t procped) +{ + fProcessPedestal = procped; +} + +// -------------------------------------------------------------------------- +// +// SetProcessTime. +// +// set the option "ProcessTime"; +// it is used in : MBadPixelsTreat +// +void MMakePadHistograms::SetProcessTime(Bool_t proctim) +{ + fProcessTime = proctim; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MMakePadHistograms::MakeHistograms() +{ + if (fType == "") + { + *fLog << err << "MMakePadHistograms::MakeHistograms; type of input data is not defined" + << endl; + return kFALSE; + } + + if (fNameInputFile == "") + { + *fLog << err << "MMakePadHistograms::MakeHistograms; name of input file is not defined" + << endl; + return kFALSE; + } + + if (fNameOutputFile == "") + { + *fLog << err << "MMakePadHistograms::MakeHistograms; name of output file is not defined" + << endl; + return kFALSE; + } + + *fLog << "------------" << endl; + *fLog << fType << " events :" << endl; + *fLog << "------------" << endl; + *fLog << "MMakePadHistograms : calibrated " << fType + << " data are read from file " << fNameInputFile << endl; + + + MTaskList tlist; + MParList plist; + + MReadReports read; + read.AddTree("Drive"); + read.AddTree("Events","MTime.",kTRUE); + static_cast(read).AddFile(fNameInputFile); + + read.AddToBranchList("MReportDrive.*"); + read.AddToBranchList("MRawEvtHeader.*"); + +// MReadMarsFile read("Events", fNameInputFile); +// read.DisableAutoScheme(); + + + MGeomApply apply; + +// MExtrapolatePointingPos extraPointing(fNameInputFile); + MPointingPosCalc pointposcalc; + + MBadPixelsCalc badcalc; + badcalc.SetNamePedPhotContainer(fNamePedPhotCam); + badcalc.SetPedestalLevel(fPedestalLevel); + + MBadPixelsTreat badtreat; + badtreat.SetNamePedPhotCam(fNamePedPhotCam); + badtreat.SetUseInterpolation(fUseInterpolation); + badtreat.SetProcessPedestal(fProcessPedestal); + badtreat.SetProcessTimes(fProcessTime); + + MFSelBasic selbasic; + MContinue contbasic(&selbasic); + contbasic.SetName("SelBasic"); + + MHBadPixels bad("BadPixels"); + bad.SetNamePedPhotCam(fNamePedPhotCam); + MFillH fillbad("BadPixels[MHBadPixels]", "MBadPixelsCam"); + fillbad.SetName("FillBad"); + + MHSigmaTheta sigth("SigmaTheta"); + sigth.SetNamePedPhotCam(fNamePedPhotCam); + MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", ""); + fillsigtheta.SetName("FillSigTheta"); + + //***************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(&bad); + plist.AddToList(&sigth); + + //***************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&apply); + tlist.AddToList(&pointposcalc, "Events"); + +// tlist.AddToList(&badcalc); done in callisto +// tlist.AddToList(&badtreat); done in callisto + + tlist.AddToList(&contbasic, "Events"); + tlist.AddToList(&fillbad, "Events"); + tlist.AddToList(&fillsigtheta, "Events"); + //***************************** + + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetProgressBar(&bar); + + // if ( !evtloop.Eventloop(fMaxEvents) ) + // return kFALSE; + evtloop.Eventloop(fMaxEvents); + + tlist.PrintStatistics(0, kTRUE); + + bad.DrawClone(); + sigth.DrawClone(); + + // save the histograms for the padding + fHSigmaTheta = (TH2D*)sigth.GetSigmaTheta(); + fHSigmaThetaOuter= (TH2D*)sigth.GetSigmaThetaOuter(); + fHSigmaPixTheta = (TH3D*)sigth.GetSigmaPixTheta(); + fHDiffPixTheta = (TH3D*)sigth.GetDiffPixTheta(); + + fHBadIdTheta = (TH2D*)bad.GetBadId(); + fHBadNTheta = (TH2D*)bad.GetBadN(); + + TFile filewrite(fNameOutputFile, "RECREATE", ""); + fHSigmaTheta->Write(); + fHSigmaThetaOuter->Write(); + //fHSigmaPixTheta->Write(); + fHDiffPixTheta->Write(); + //fHBadIdTheta->Write(); + //fHBadNTheta->Write(); + filewrite.Close(); + + *fLog << "" << endl; + *fLog << "MMakePadHistograms : padding plots for " << fType + << " data written onto file " << fNameOutputFile << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.h (revision 5431) +++ /trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.h (revision 5431) @@ -0,0 +1,73 @@ +#ifndef MARS_MMakePadHistograms +#define MARS_MMakePadHistograms + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH2D; +class TH3D; + +class MMakePadHistograms : public MH +{ +private: + + TString fType; // type of data to be padded (ON, OFF, MC) + Int_t fMaxEvents; // maximum no.of events to be read + TString fNamePedPhotCam; // name of the 'MPedPhotCam' container + TString fNameInputFile; // name of input file (calibrated data) + TString fNameOutputFile; // name of output file (Pad histograms) + Float_t fPedestalLevel; // pedestal level for MBadPixelsCalc + Bool_t fUseInterpolation; // option for MBadPixelsTreat + Bool_t fProcessPedestal; // option for MBadPixelsTreat + Bool_t fProcessTime; // option for MBadPixelsTreat + + + //---------------------------------- + // plots needed for the padding + + TH2D *fHSigmaTheta; // 2D-histogram (Theta, sigmabar(inner)) + TH2D *fHSigmaThetaOuter; // 2D-histogram (Theta, sigmabar(outer)) + TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, (sigma^2-sigmabar^2)/Area) + TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) + + TH2D *fHBadIdTheta; // 2D-histogram (Theta, Id of bad pixel) + TH2D *fHBadNTheta; // 2D-histogram (Theta, no.of bad pixels) + + +public: + MMakePadHistograms(const char *name=NULL, const char *title=NULL); + + void SetDataType(const char *); // type of data (ON, OFF, MC) + void SetMaxEvents(Int_t); // mx. no.of events to be read + void SetNameInputFile(const char *); // Calibrated data + void SetNameOutputFile(const char *); // Pad histograms + void SetNamePedPhotCam(const char *); // "MPedPhotCamFromData" + + // parameters for MBadPixelsCalc + void SetPedestalLevel(Float_t); + + // parameters for MBadPixelsTreat + void SetUseInterpolation(Bool_t); + void SetProcessPedestal(Bool_t); + void SetProcessTime(Bool_t); + + Bool_t MakeHistograms(); + + ClassDef(MMakePadHistograms, 0) // class for making the padding histograms +}; + +#endif + + + + + + + + + Index: unk/MagicSoft/Mars/manalysis/MPad.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 5430) +++ (revision ) @@ -1,2282 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Wolfgang Wittek, 06/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MPad -// -// The aim of the padding is to make the noise level (NSB + electronic noise) -// equal for the MC, ON and OFF data samples -// -// The target noise level is determined by 'merging' the (sigmabar vs. Theta) -// distributions of MC, ON and OFF data. -// -// The prescription on which fraction of events has to padded from -// bin k of sigmabar to bin j is stored in the matrices -// fHgMC, fHgON and fHgOFF. -// -// By the padding the number of photons, its error and the pedestal sigmas -// are altered. On average, the number of photons added is zero. -// -// The formulas used can be found in Thomas Schweizer's Thesis, -// Section 2.2.1 -// -// The padding is done as follows : -// -// Increase the sigmabar according to the matrices fHg.... Then generate -// a pedestal sigma for each pixel using the 3D-histogram Theta, pixel no., -// Sigma^2-Sigmabar^2 (fHDiffPixTheta). -// -// The padding has to be done before the image cleaning because the -// image cleaning depends on the pedestal sigmas. -// -///////////////////////////////////////////////////////////////////////////// -#include "MPad.h" - -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include "MBinning.h" -#include "MPointingPos.h" -#include "MLog.h" -#include "MLogManip.h" -#include "MParList.h" -#include "MGeomCam.h" - -#include "MCerPhotPix.h" -#include "MCerPhotEvt.h" - -#include "MH.h" - -#include "MPedPhotCam.h" -#include "MPedPhotPix.h" - -#include "MBadPixelsCam.h" -#include "MBadPixelsPix.h" - -#include "MRead.h" -#include "MFilterList.h" -#include "MTaskList.h" -#include "MHBadPixels.h" -#include "MFillH.h" -#include "MHSigmaTheta.h" -#include "MEvtLoop.h" - -ClassImp(MPad); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MPad::MPad(const char *name, const char *title) -{ - fName = name ? name : "MPad"; - fTitle = title ? title : "Task for the padding"; - - fType = ""; - - fHSigmaThetaMC = NULL; - fHSigmaThetaON = NULL; - fHSigmaThetaOFF = NULL; - - fHDiffPixThetaMC = NULL; - fHDiffPixThetaON = NULL; - fHDiffPixThetaOFF = NULL; - - fHSigmaPixThetaMC = NULL; - fHSigmaPixThetaON = NULL; - fHSigmaPixThetaOFF = NULL; - - fHgMC = NULL; - fHgON = NULL; - fHgOFF = NULL; - - fHSigmaPedestal = NULL; - fHPhotons = NULL; - fHNSB = NULL; - - fNamePedPhotCam = "MPedPhotCamFromData"; -} - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MPad::~MPad() -{ - delete fHgMC; - delete fHgON; - delete fHgOFF; - - delete fHSigmaTheta; - delete fHSigmaPixTheta; - delete fHDiffPixTheta; - - delete fHSigmaPedestal; - delete fHPhotons; - delete fHNSB; - - delete fInfile; -} -// -------------------------------------------------------------------------- -// -// Merge the sigmabar distributions of 2 samples (samples A and B) -// -// input : the original distributions for samples A and B (hista, histb) -// -// output : the prescription which fraction of events should be padded -// from the sigmabar bin k to bin j (fHgMC, fHgON, fHgOFF) -// - -Bool_t MPad::Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap, - TH2D *fHga, TH2D *fHgb, Int_t nbinssig ) -{ - // hista is the normalized 1D histogram of sigmabar for sample A - // histb is the normalized 1D histogram of sigmabar for sample B - // histc is the difference A-B - - // at the beginning, histap is a copy of hista - // at the end, it will be the 1D histogram for sample A after padding - - // at the beginning, histbp is a copy of histb - // at the end, it will be the 1D histogram for sample B after padding - - // at the beginning, histcp is a copy of histc - // at the end, it should be the difference histap-histbp, - // which should be zero - - // fHga[k][j] tells which fraction of events from sample A should be padded - // from sigma_k to sigma_j - - - // fHgb[k][j] tells which fraction of events from sample B should be padded - // from sigma_k to sigma_j - - //*fLog << all << "MPad::Merge2Distributions(); Entry" << endl; - - Double_t eps = 1.e-10; - - TH1D *histbp = new TH1D( (const TH1D&)*histb ); - - TH1D *histc = new TH1D( (const TH1D&)*hista ); - histc->Add(histb, -1.0); - - TH1D *histcp = new TH1D( (const TH1D&)*histc ); - - - //-------- start j loop ------------------------------------------------ - // loop over bins in histc, - // starting from the end (i.e. at the largest sigmabar) - Double_t v, s, w, t, x, u, a, b, arg; - - //*fLog << "Content of histcp : " << endl; - - for (Int_t j=nbinssig; j >= 1; j--) - { - //*fLog << "j, hista, histb , histap : " << j << ", " - // << hista->GetBinContent(j) << ", " - // << histb->GetBinContent(j) << ", " - // << histap->GetBinContent(j) << endl; - - v = histcp->GetBinContent(j); - //*fLog << "cp : j, v = " << j << ", " << v << endl; - - if ( fabs(v) < eps ) continue; - if (v >= 0.0) - s = 1.0; - else - s = -1.0; - - //................ start k loop ...................................... - // look for a bin k which may compensate the content of bin j - for (Int_t k=j-1; k >= 1; k--) - { - w = histcp->GetBinContent(k); - if ( fabs(w) < eps ) continue; - if (w >= 0.0) - t = 1.0; - else - t = -1.0; - - // if s==t bin k cannot compensate, go to next k bin - if (s == t) continue; - - x = v + w; - if (x >= 0.0) - u = 1.0; - else - u = -1.0; - - // if u==s bin k will partly compensate : pad the whole fraction - // w from bin k to bin j - if (u == s) - arg = -w; - - // if u!=s bin k would overcompensate : pad only the fraction - // v from bin k to bin j - else - arg = v; - - if (arg <=0.0) - { - fHga->SetBinContent(k, j, -arg); - fHgb->SetBinContent(k, j, 0.0); - } - else - { - fHga->SetBinContent(k, j, 0.0); - fHgb->SetBinContent(k, j, arg); - } - - //*fLog << all << "k, j, arg = " << k << ", " << j - // << ", " << arg << endl; - - //...................................... - // this is for checking the procedure - if (arg < 0.0) - { - a = histap->GetBinContent(k); - histap->SetBinContent(k, a+arg); - a = histap->GetBinContent(j); - histap->SetBinContent(j, a-arg); - } - else - { - b = histbp->GetBinContent(k); - histbp->SetBinContent(k, b-arg); - b = histbp->GetBinContent(j); - histbp->SetBinContent(j, b+arg); - } - //...................................... - - if (u == s) - { - histcp->SetBinContent(k, 0.0); - histcp->SetBinContent(j, x); - - v = histcp->GetBinContent(j); - if ( fabs(v) < eps ) break; - - // bin j was only partly compensated : go to next k bin - continue; - } - else - { - histcp->SetBinContent(k, x); - histcp->SetBinContent(j, 0.0); - - // bin j was completely compensated : go to next j bin - break; - } - - } - //................ end k loop ...................................... - } - //-------- end j loop ------------------------------------------------ - - // check results - - //*fLog << "Content of histap, histbp, histcp : " << endl; - - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = histap->GetBinContent(j); - Double_t b = histbp->GetBinContent(j); - Double_t c = histcp->GetBinContent(j); - - //*fLog << "j, a, b, c = " << j << ": " << a << ", " << b << ", " - // << c << endl; - - if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps ) - *fLog << err << "MPad::Merge2Distributions(); inconsistency in results; j, a, b, c = " - << j << ", " << a << ", " << b << ", " << c << endl; - } - - //--------------------------------------------------------------- - TCanvas *pad = MH::MakeDefCanvas("Merging", "", 600, 600); - gROOT->SetSelectedPad(NULL); - - pad->Divide(2, 2); - - pad->cd(1); - gPad->SetBorderMode(0); - hista->SetDirectory(NULL); - hista->Draw(); - hista->SetBit(kCanDelete); - - pad->cd(2); - gPad->SetBorderMode(0); - histb->SetDirectory(NULL); - histb->Draw(); - histb->SetBit(kCanDelete); - - pad->cd(3); - gPad->SetBorderMode(0); - histap->SetDirectory(NULL); - histap->Draw(); - histap->SetBit(kCanDelete); - - pad->DrawClone(); - - //-------------------------------------------------------------------- - - delete histc; - delete histbp; - delete histcp; - - //*fLog << all << "MPad::Merge2Distributions(); Exit" << endl; - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Merge the distributions -// fHSigmaTheta 2D-histogram (Theta, sigmabar) -// -// ===> of ON, OFF and MC data <============== -// -// and define the matrices fHgMC, fHgON and fHgOFF -// -// These matrices tell which fraction of events should be padded -// from bin k of sigmabar to bin j -// - -Bool_t MPad::MergeONOFFMC( - TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, - TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, - TH2D& sigthoff, TH3D& diffpixthoff, TH3D& sigmapixthoff) -{ - //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - // for testing - /* - TAxis *xa = sigthon.GetXaxis(); - Int_t nbitheta = xa->GetNbins(); - - TAxis *ya = sigthon.GetYaxis(); - Int_t nbisig = ya->GetNbins(); - for (Int_t i=1; i<=nbitheta; i++) - { - for (Int_t j=1; j<=nbisig; j++) - { - if (i>0) - { - sigthmc.SetBinContent( i, j, (Float_t) (625000.0+(nbisig*nbisig*nbisig-j*j*j)) ); - sigthon.SetBinContent( i, j, (Float_t)(1.0) ); - sigthoff.SetBinContent( i, j, (Float_t)( (0.5*nbisig-j)*(0.5*nbisig-j) ) ); - } - else - { - sigthmc.SetBinContent( i, j, 0.0); - sigthon.SetBinContent( i, j, 0.0); - sigthoff.SetBinContent(i, j, 0.0); - } - } - } - */ - //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - - *fLog << all << "----------------------------------------------------------------------------------" << endl; - *fLog << all << "MPad::MergeONOFFMC(); Merge the ON, OFF and MC histograms to obtain the target distributions for the padding" - << endl; - - fHSigmaTheta = new TH2D( (const TH2D&) sigthon ); - fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); - - *fLog << "after booking fHSigmaTheta" << endl; - - fHDiffPixTheta = new TH3D( (const TH3D&) diffpixthon ); - fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); - - *fLog << "after booking fHDiffPixTheta" << endl; - - fHSigmaPixTheta = new TH3D( (const TH3D&) sigmapixthon ); - fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); - *fLog << "after booking fHSigmaPixTheta" << endl; - - *fLog << all << "Histograms for the merged padding plots were booked" - << endl; - - //-------------------------- - fHSigmaThetaMC = &sigthmc; - fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); - fHSigmaThetaON = &sigthon; - fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); - fHSigmaThetaOFF = &sigthoff; - fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (orig.)"); - - //*fLog << all << "addresses of SigmaTheta padding plots were copied" - // << endl; - - //-------------------------- - fHDiffPixThetaMC = &diffpixthmc; - fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); - fHDiffPixThetaON = &diffpixthon; - fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); - fHDiffPixThetaOFF = &diffpixthoff; - fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (orig.)"); - - //*fLog << all << "addresses of DiffPixTheta padding plots were copied" - // << endl; - - //-------------------------- - fHSigmaPixThetaMC = &sigmapixthmc; - fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); - fHSigmaPixThetaON = &sigmapixthon; - fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); - fHSigmaPixThetaOFF = &sigmapixthoff; - fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (orig.)"); - - //*fLog << all << "addresses of SigmaPixTheta padding plots were copied" - // << endl; - - //-------------------------- - - - // get binning for fHgON, fHgOFF and fHgMC from sigthon - // binning in Theta - TAxis *ax = sigthon.GetXaxis(); - Int_t nbinstheta = ax->GetNbins(); - TArrayD edgesx; - edgesx.Set(nbinstheta+1); - for (Int_t i=0; i<=nbinstheta; i++) - { - edgesx[i] = ax->GetBinLowEdge(i+1); - //*fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] - // << endl; - } - MBinning binth; - binth.SetEdges(edgesx); - - // binning in sigmabar - TAxis *ay = sigthon.GetYaxis(); - Int_t nbinssig = ay->GetNbins(); - TArrayD edgesy; - edgesy.Set(nbinssig+1); - for (Int_t i=0; i<=nbinssig; i++) - { - edgesy[i] = ay->GetBinLowEdge(i+1); - //*fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] - // << endl; - } - MBinning binsg; - binsg.SetEdges(edgesy); - - - fHgMC = new TH3D; - MH::SetBinning(fHgMC, &binth, &binsg, &binsg); - fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC"); - - fHgON = new TH3D; - MH::SetBinning(fHgON, &binth, &binsg, &binsg); - fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); - - fHgOFF = new TH3D; - MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); - fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); - - //*fLog << all << "fHg histograms were booked" << endl; - - //............ loop over Theta bins ........................... - - - - //*fLog << all << "MPad::MergeONOFFMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; - for (Int_t l=1; l<=nbinstheta; l++) - { - //------------------------------------------- - // merge ON and OFF distributions - // input : hista is the normalized 1D distr. of sigmabar for ON data - // histb is the normalized 1D distr. of sigmabar for OFF data - // output : histap will be the merged distribution (ON-OFF) - // fHga(k,j) will tell which fraction of ON events should be - // padded from bin k to bin j - // fHgb(k,j) will tell which fraction of OFF events should be - // padded from bin k to bin j - - TH2D * fHga = new TH2D; - MH::SetBinning(fHga, &binsg, &binsg); - TH2D * fHgb = new TH2D; - MH::SetBinning(fHgb, &binsg, &binsg); - - TH1D *hista = sigthon.ProjectionY("sigon_y", l, l, ""); - hista->SetNameTitle("ON-orig", "ON (orig)"); - Stat_t suma = hista->Integral(); - if (suma != 0.0) hista->Scale(1./suma); - - TH1D *histap = new TH1D( (const TH1D&)*hista ); - histap->SetNameTitle("ON-OFF-merged)", "ON-OFF (merged)"); - - TH1D *histb = sigthoff.ProjectionY("sigoff_y", l, l, ""); - histb->SetNameTitle("OFF-orig", "OFF (orig)"); - Stat_t sumb = histb->Integral(); - if (sumb != 0.0) histb->Scale(1./sumb); - - if (suma == 0.0 || sumb == 0.0) - { - delete hista; - delete histb; - delete fHga; - delete fHgb; - delete histap; - - *fLog << err - << "cannot call Merge2Distributions(ON, OFF) for theta bin l = " - << l << "; NevON, NevOFF = " << suma << ", " << sumb << endl; - - // go to next theta bin (l) - continue; - } - - //*fLog << "call Merge2Distributions(ON, OFF) for theta bin l = " - // << l << endl; - - Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); - - // fill fHgON and fHgOFF - for (Int_t k=1; k<=nbinssig; k++) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = fHga->GetBinContent(k,j); - fHgON->SetBinContent(l, k, j, a); - - Double_t b = fHgb->GetBinContent(k,j); - fHgOFF->SetBinContent(l, k, j, b); - } - - //------------------------------------------- - // now merge ON-OFF and MC distributions - // input : histe is the result of the merging of ON and OFF distributions - // histf is the normalized 1D distr. of sigmabar for MC data - // output : histep will be the merged distribution (target distribution) - // fHge(k,j) will tell which fraction of ON, OFF events should be - // padded from bin k to bin j - // fHgf(k,j) will tell which fraction of MC events should be - // padded from bin k to bin j - - TH2D * fHge = new TH2D; - MH::SetBinning(fHge, &binsg, &binsg); - TH2D * fHgf = new TH2D; - MH::SetBinning(fHgf, &binsg, &binsg); - - TH1D *histe = new TH1D( (const TH1D&)*histap); - histe->SetNameTitle("ON-OFF-merged", "ON-OFF (merged)"); - - TH1D *histep = new TH1D( (const TH1D&)*histe); - histep->SetNameTitle("ON-OFF-MC-merged)", "ON-OFF-MC (merged)"); - - TH1D *histf = sigthmc.ProjectionY("sigmc_y", l, l, ""); - histf->SetNameTitle("MC-orig", "MC (orig)"); - Double_t sumf = histf->Integral(); - if (sumf != 0.0) histf->Scale(1./sumf); - - if (sumf == 0.0) - { - delete hista; - delete histb; - delete fHga; - delete fHgb; - delete histap; - - delete histe; - delete histf; - delete fHge; - delete fHgf; - delete histep; - - *fLog << err - << "cannot call Merge2Distributions(ON-OFF,MC) for theta bin l = " - << l << "; NevMC = " << sumf << endl; - - // go to next theta bin (l) - continue; - } - - //*fLog << "call Merge2Distributions(ON-OFF, MC) for theta bin l = " - // << l << endl; - - Merge2Distributions(histe, histf, histep, fHge, fHgf, nbinssig); - - // update fHgON and fHgOFF - UpdateHg(fHga, histap, fHge, fHgON, nbinssig, l); - UpdateHg(fHgb, histap, fHge, fHgOFF, nbinssig, l); - - // fill fHgMC - for (Int_t k=1; k<=nbinssig; k++) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t f = fHgf->GetBinContent(k,j); - fHgMC->SetBinContent(l, k, j, f); - } - - // fill target distribution fHSigmaTheta - // (only for plotting, it will not be used in the padding) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t ep = histep->GetBinContent(j); - fHSigmaTheta->SetBinContent(l, j, ep); - } - //------------------------------------------- - - delete hista; - delete histb; - delete fHga; - delete fHgb; - delete histap; - - delete histe; - delete histf; - delete fHge; - delete fHgf; - delete histep; - } - //............ end of loop over Theta bins .................... - - - // Define the target distributions - // fHDiffPixTheta, fHSigmaPixTheta - // (they are calculated as - // averages of the ON and OFF distributions) - // (fHDiffPixTheta will be used in the padding of MC events) - - //............ start of new loop over Theta bins .................... - for (Int_t j=1; j<=nbinstheta; j++) - { - // fraction of ON/OFF/MC events to be padded : fracON, fracOFF, fracMC - - Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); - Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); - Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); - - Double_t normON; - Double_t normOFF; - - // set ranges for 2D-projections of 3D-histograms - TAxis *ax = diffpixthon.GetXaxis(); - ax->SetRange(j, j); - ax = diffpixthoff.GetXaxis(); - ax->SetRange(j, j); - - TH2D *hist; - TH2D *histOFF; - - - // weights for ON and OFF distrubtions when - // calculating the weighted average - Double_t facON = 0.5 * (1. - fracON + fracOFF); - Double_t facOFF = 0.5 * (1. + fracON - fracOFF); - - *fLog << all << "Theta bin j : fracON, fracOFF, fracMC, facON, facOFF = " - << j << ", " << fracON << ", " << fracOFF << ", " - << fracMC << ", " << facON << ", " << facOFF << endl; - - //------------------------------------------------------------------ - // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' - ay = diffpixthon.GetYaxis(); - Int_t nbinspixel = ay->GetNbins(); - - TAxis *az = diffpixthon.GetZaxis(); - Int_t nbinsdiff = az->GetNbins(); - - hist = (TH2D*)diffpixthon.Project3D("zy"); - hist->SetName("dummy"); - histOFF = (TH2D*)diffpixthoff.Project3D("zy"); - - normON = hist->Integral(); - normOFF = histOFF->Integral(); - if (normON != 0.0) hist->Scale(1.0/normON); - if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); - - // weighted average of ON and OFF distributions - hist->Scale(facON); - hist->Add(histOFF, facOFF); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinsdiff; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHDiffPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histOFF; - - //------------------------------------------------------------------ - // define target distribution 'sigma vs. pixel number' - ay = sigmapixthon.GetYaxis(); - nbinspixel = ay->GetNbins(); - - az = sigmapixthon.GetZaxis(); - Int_t nbinssigma = az->GetNbins(); - - hist = (TH2D*)sigmapixthon.Project3D("zy"); - hist->SetName("dummy"); - histOFF = (TH2D*)sigmapixthoff.Project3D("zy"); - - normON = hist->Integral(); - normOFF = histOFF->Integral(); - if (normON != 0.0) hist->Scale(1.0/normON); - if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); - - // weighted average of ON and OFF distributions - hist->Scale(facON); - hist->Add(histOFF, facOFF); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinssigma; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHSigmaPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histOFF; - } - //............ end of new loop over Theta bins .................... - - *fLog << all - << "MPad::MergeONOFFMC(); The target distributions for the padding have been created" - << endl; - *fLog << all << "----------------------------------------------------------" - << endl; - //-------------------------------------------- - - fHSigmaThetaMC->SetDirectory(NULL); - fHSigmaThetaON->SetDirectory(NULL); - fHSigmaThetaOFF->SetDirectory(NULL); - - fHDiffPixThetaMC->SetDirectory(NULL); - fHDiffPixThetaON->SetDirectory(NULL); - fHDiffPixThetaOFF->SetDirectory(NULL); - - fHSigmaPixThetaMC->SetDirectory(NULL); - fHSigmaPixThetaON->SetDirectory(NULL); - fHSigmaPixThetaOFF->SetDirectory(NULL); - - fHgMC->SetDirectory(NULL); - fHgON->SetDirectory(NULL); - fHgOFF->SetDirectory(NULL); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Update the matrix fHgA -// which contains the final padding prescriptions -// - -Bool_t MPad::UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA, - Int_t nbinssig, Int_t l) -{ - // histap target distribution after ON-OFF merging - // fHga padding prescription for ON (or OFF) data to get to histap - // fHge padding prescription to get from histap to the target - // distribution after the ON-OFF-MC merging - // fHgA updated padding prescription for ON (or OFF) data to get - // from the original ON (or OFF) histogram to the target - // distribution after the ON-OFF-MC merging - // l Theta bin - - Double_t eps = 1.e-10; - - for (Int_t k=1; k<=nbinssig; k++) - { - Double_t na = fHga->Integral(1, nbinssig, k, k, " "); - Double_t ne = fHge->Integral(k, k, 1, nbinssig, " "); - Double_t nap = histap->GetBinContent(k); - - if (ne <= eps) - { - // go to next k - continue; - } - - Double_t r = nap - na; - - if (r >= ne-eps) - { - for (Int_t j=k+1; j<=nbinssig; j++) - { - Double_t e = fHge->GetBinContent(k,j); - Double_t A = fHgA->GetBinContent(l,k,j); - fHgA->SetBinContent(l, k, j, A + e); - } - // go to next k - continue; - } - - for (Int_t j=k+1; j<=nbinssig; j++) - { - Double_t e = fHge->GetBinContent(k,j); - Double_t A = fHgA->GetBinContent(l,k,j); - fHgA->SetBinContent(l, k, j, A + r*e/ne); - } - - if (na <= eps) - { - // go to next k - continue; - } - - for (Int_t i=1; iGetBinContent(i,k); - Double_t A = fHgA->GetBinContent(l,i,k); - fHgA->SetBinContent(l, i, k, A - (ne-r)*a/na); - } - - for (Int_t i=1; i<=nbinssig; i++) - { - if (i == k) continue; - for (Int_t j=i+1; j<=nbinssig; j++) - { - if (j == k) continue; - Double_t a = fHga->GetBinContent(i,k); - Double_t e = fHge->GetBinContent(k,j); - Double_t A = fHgA->GetBinContent(l,i,j); - fHgA->SetBinContent(l, i, j, A + (ne-r)*a*e/(na*ne) ); - } - } - } - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Merge the distributions -// fHSigmaTheta 2D-histogram (Theta, sigmabar) -// -// ===> of ON and MC data <============== -// -// and define the matrices fHgMC and fHgON -// -// These matrices tell which fraction of events should be padded -// from bin k of sigmabar to bin j -// - -Bool_t MPad::MergeONMC( - TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, - TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon) -{ - *fLog << all << "----------------------------------------------------------------------------------" << endl; - *fLog << all << "MPad::MergeONMC(); Merge the ON and MC histograms to obtain the target distributions for the padding" - << endl; - - fHSigmaTheta = new TH2D( (TH2D&) sigthon ); - fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); - - fHDiffPixTheta = new TH3D( (TH3D&) diffpixthon ); - fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); - - fHSigmaPixTheta = new TH3D( (TH3D&) sigmapixthon ); - fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); - - //-------------------------- - fHSigmaThetaMC = &sigthmc; - fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); - fHSigmaThetaON = &sigthon; - fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); - - //-------------------------- - fHDiffPixThetaMC = &diffpixthmc; - fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); - fHDiffPixThetaON = &diffpixthon; - fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); - - //-------------------------- - fHSigmaPixThetaMC = &sigmapixthmc; - fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); - fHSigmaPixThetaON = &sigmapixthon; - fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); - - //-------------------------- - - - // get binning for fHgON and fHgMC from sigthon - // binning in Theta - TAxis *ax = sigthon.GetXaxis(); - Int_t nbinstheta = ax->GetNbins(); - TArrayD edgesx; - edgesx.Set(nbinstheta+1); - for (Int_t i=0; i<=nbinstheta; i++) - { - edgesx[i] = ax->GetBinLowEdge(i+1); - *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] - << endl; - } - MBinning binth; - binth.SetEdges(edgesx); - - // binning in sigmabar - TAxis *ay = sigthon.GetYaxis(); - Int_t nbinssig = ay->GetNbins(); - TArrayD edgesy; - edgesy.Set(nbinssig+1); - for (Int_t i=0; i<=nbinssig; i++) - { - edgesy[i] = ay->GetBinLowEdge(i+1); - *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] - << endl; - } - MBinning binsg; - binsg.SetEdges(edgesy); - - - fHgMC = new TH3D; - MH::SetBinning(fHgMC, &binth, &binsg, &binsg); - fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC"); - - fHgON = new TH3D; - MH::SetBinning(fHgON, &binth, &binsg, &binsg); - fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); - - - //............ loop over Theta bins ........................... - - - - *fLog << all << "MPad::MergeONMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; - for (Int_t l=1; l<=nbinstheta; l++) - { - TH2D * fHga = new TH2D; - MH::SetBinning(fHga, &binsg, &binsg); - TH2D * fHgb = new TH2D; - MH::SetBinning(fHgb, &binsg, &binsg); - - //------------------------------------------- - // merge ON and MC distributions - // input : hista is the normalized 1D distr. of sigmabar for ON data - // histb is the normalized 1D distr. of sigmabar for MC data - // output : histap will be the merged distribution (ON-MC) - // fHga(k,j) will tell which fraction of ON events should be - // padded from bin k to bin j - // fHgb(k,j) will tell which fraction of MC events should be - // padded from bin k to bin j - - TH1D *hista = sigthon.ProjectionY("sigon_y", l, l, ""); - Stat_t suma = hista->Integral(); - if (suma != 0.0) hista->Scale(1./suma); - - TH1D *histap = new TH1D( (const TH1D&)*hista ); - - TH1D *histb = sigthmc.ProjectionY("sigmc_y", l, l, ""); - Stat_t sumb = histb->Integral(); - if (sumb != 0.0) histb->Scale(1./sumb); - - Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); - delete hista; - delete histb; - - // fill fHgON and fHgMC - for (Int_t k=1; k<=nbinssig; k++) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = fHga->GetBinContent(k,j); - fHga->SetBinContent (k, j, 0.0); - fHgON->SetBinContent(l, k, j, a); - - Double_t b = fHgb->GetBinContent(k,j); - fHgb->SetBinContent (k, j, 0.0); - fHgMC->SetBinContent(l, k, j, b); - } - - - // fill target distribution fHSigmaTheta - // (only for plotting, it will not be used in the padding) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = histap->GetBinContent(j); - fHSigmaTheta->SetBinContent(l, j, a); - } - - delete fHga; - delete fHgb; - - delete histap; - } - //............ end of loop over Theta bins .................... - - - // Define the target distributions - // fHDiffPixTheta, fHSigmaPixTheta - // (they are calculated as - // averages of the ON and MCdistributions) - // (fHDiffPixTheta will be used in the padding of MC events) - - //............ start of new loop over Theta bins .................... - for (Int_t j=1; j<=nbinstheta; j++) - { - // fraction of ON/MC events to be padded : fracON, fracMC - - Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); - Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); - - Double_t normON; - Double_t normMC; - - // set ranges for 2D-projections of 3D-histograms - TAxis *ax = diffpixthon.GetXaxis(); - ax->SetRange(j, j); - ax = diffpixthmc.GetXaxis(); - ax->SetRange(j, j); - - TH2D *hist; - TH2D *histMC; - - - // weights for ON and MC distrubtions when - // calculating the weighted average - Double_t facON = 0.5 * (1. - fracON + fracMC); - Double_t facMC = 0.5 * (1. + fracON - fracMC); - - *fLog << all << "Theta bin j : fracON, fracMC, facON, facMC = " - << j << ", " << fracON << ", " << fracMC << ", " - << facON << ", " << facMC << endl; - - //------------------------------------------------------------------ - // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' - ay = diffpixthon.GetYaxis(); - Int_t nbinspixel = ay->GetNbins(); - - TAxis *az = diffpixthon.GetZaxis(); - Int_t nbinsdiff = az->GetNbins(); - - hist = (TH2D*)diffpixthon.Project3D("zy"); - hist->SetName("dummy"); - histMC = (TH2D*)diffpixthmc.Project3D("zy"); - - normON = hist->Integral(); - normMC = histMC->Integral(); - if (normON != 0.0) hist->Scale(1.0/normON); - if (normMC != 0.0) histMC->Scale(1.0/normMC); - - // weighted average of ON and MC distributions - hist->Scale(facON); - hist->Add(histMC, facMC); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinsdiff; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHDiffPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histMC; - - //------------------------------------------------------------------ - // define target distribution 'sigma vs. pixel number' - ay = sigmapixthon.GetYaxis(); - nbinspixel = ay->GetNbins(); - - az = sigmapixthon.GetZaxis(); - Int_t nbinssigma = az->GetNbins(); - - hist = (TH2D*)sigmapixthon.Project3D("zy"); - hist->SetName("dummy"); - histMC = (TH2D*)sigmapixthmc.Project3D("zy"); - - normON = hist->Integral(); - normMC = histMC->Integral(); - if (normON != 0.0) hist->Scale(1.0/normON); - if (normMC != 0.0) histMC->Scale(1.0/normMC); - - // weighted average of ON and MC distributions - hist->Scale(facON); - hist->Add(histMC, facMC); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinssigma; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHSigmaPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histMC; - //------------------------------------------------------------------ - } - //............ end of new loop over Theta bins .................... - - *fLog << all - << "MPad::MergeONMC(); The target distributions for the padding have been created" - << endl; - *fLog << all << "----------------------------------------------------------" - << endl; - //-------------------------------------------- - - fHSigmaThetaMC->SetDirectory(NULL); - fHSigmaThetaON->SetDirectory(NULL); - - fHDiffPixThetaMC->SetDirectory(NULL); - fHDiffPixThetaON->SetDirectory(NULL); - - fHSigmaPixThetaMC->SetDirectory(NULL); - fHSigmaPixThetaON->SetDirectory(NULL); - - fHgMC->SetDirectory(NULL); - fHgON->SetDirectory(NULL); - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Read padding distributions from a file -// -// -Bool_t MPad::ReadPaddingDist(const char* namefilein) -{ - *fLog << all << "MPad : Read padding histograms from file " - << namefilein << endl; - - fInfile = new TFile(namefilein); - fInfile->ls(); - - //------------------------------------ - - /* - fHSigmaThetaMC = - (TH2D*) gROOT->FindObject("2D-ThetaSigmabarMC"); - if (!fHSigmaThetaMC) - { - *fLog << all - << "MPad : Object '2D-ThetaSigmabarMC' not found on root file" - << endl; - return kFALSE; - } - */ - fHSigmaThetaMC = new TH2D; - fHSigmaThetaMC->Read("2D-ThetaSigmabarMC"); - *fLog << all - << "MPad : Object '2D-ThetaSigmabarMC' was read in" << endl; - - /* - fHSigmaThetaON = - (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON"); - if (!fHSigmaThetaON) - { - *fLog << all - << "MPad : Object '2D-ThetaSigmabarON' not found on root file" - << endl; - return kFALSE; - } - */ - fHSigmaThetaON = new TH2D; - fHSigmaThetaON->Read("2D-ThetaSigmabarON"); - *fLog << all - << "MPad : Object '2D-ThetaSigmabarON' was read in" << endl; - - /* - fHSigmaThetaOFF = - (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF"); - if (!fHSigmaThetaOFF) - { - *fLog << all - << "MPad : Object '2D-ThetaSigmabarOFF' not found on root file" - << endl; - return kFALSE; - } - */ - fHSigmaThetaOFF = new TH2D; - fHSigmaThetaOFF->Read("2D-ThetaSigmabarOFF"); - *fLog << all - << "MPad:Object '2D-ThetaSigmabarOFF' was read in" << endl; - - //------------------------------------ - /* - fHDiffPixTheta = - (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); - if (!fHDiffPixTheta) - { - *fLog << all - << "MPad : Object '3D-ThetaPixDiff' not found on root file" - << endl; - return kFALSE; - } - */ - fHDiffPixTheta = new TH3D; - fHDiffPixTheta->Read("3D-ThetaPixDiff"); - *fLog << all - << "MPad : Object '3D-ThetaPixDiff' was read in" << endl; - - /* - fHDiffPixThetaMC = - (TH3D*) gROOT->FindObject("3D-ThetaPixDiffMC"); - if (!fHDiffPixThetaMC) - { - *fLog << all - << "MPad : Object '3D-ThetaPixDiffMC' not found on root file" - << endl; - return kFALSE; - } - */ - fHDiffPixThetaMC = new TH3D; - fHDiffPixThetaMC->Read("3D-ThetaPixDiffMC"); - *fLog << all - << "MPad : Object '3D-ThetaPixDiffMC' was read in" << endl; - - /* - fHDiffPixThetaON = - (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON"); - if (!fHDiffPixThetaON) - { - *fLog << all - << "MPad : Object '3D-ThetaPixDiffON' not found on root file" - << endl; - return kFALSE; - } - */ - fHDiffPixThetaON = new TH3D; - fHDiffPixThetaON->Read("3D-ThetaPixDiffON"); - *fLog << all - << "MPad : Object '3D-ThetaPixDiffON' was read in" << endl; - - /* - fHDiffPixThetaOFF = - (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF"); - if (!fHDiffPixThetaOFF) - { - *fLog << all - << "MPad : Object '3D-ThetaPixDiffOFF' not found on root file" - << endl; - return kFALSE; - } - */ - fHDiffPixThetaOFF = new TH3D; - fHDiffPixThetaOFF->Read("3D-ThetaPixDiffOFF"); - *fLog << all - << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl; - - - //------------------------------------ - /* - fHSigmaPixTheta = - (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); - if (!fHSigmaPixTheta) - { - *fLog << all - << "MPad : Object '3D-ThetaPixSigma' not found on root file" - << endl; - return kFALSE; - } - */ - fHSigmaPixTheta = new TH3D; - fHSigmaPixTheta->Read("3D-ThetaPixSigma"); - *fLog << all - << "MPad : Object '3D-ThetaPixSigma' was read in" << endl; - - /* - fHSigmaPixThetaMC = - (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaMC"); - if (!fHSigmaPixThetaMC) - { - *fLog << all - << "MPad : Object '3D-ThetaPixSigmaMC' not found on root file" - << endl; - return kFALSE; - } - */ - fHSigmaPixThetaMC = new TH3D; - fHSigmaPixThetaMC->Read("3D-ThetaPixSigmaMC"); - *fLog << all - << "MPad : Object '3D-ThetaPixSigmaMC' was read in" << endl; - - /* - fHSigmaPixThetaON = - (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON"); - if (!fHSigmaPixThetaON) - { - *fLog << all - << "MPad : Object '3D-ThetaPixSigmaON' not found on root file" - << endl; - return kFALSE; - } - */ - fHSigmaPixThetaON = new TH3D; - fHSigmaPixThetaON->Read("3D-ThetaPixSigmaON"); - *fLog << all - << "MPad : Object '3D-ThetaPixSigmaON' was read in" << endl; - - /* - fHSigmaPixThetaOFF = - (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF"); - if (!fHSigmaPixThetaOFF) - { - *fLog << all - << "MPad : Object '3D-ThetaPixSigmaOFF' not found on root file" - << endl; - return kFALSE; - } - */ - fHSigmaPixThetaOFF = new TH3D; - fHSigmaPixThetaOFF->Read("3D-ThetaPixSigmaOFF"); - *fLog << all - << "MPad : Object '3D-ThetaPixSigmaOFF' was read in" << endl; - - //------------------------------------ - /* - fHgMC = - (TH3D*) gROOT->FindObject("3D-PaddingMatrixMC"); - if (!fHgMC) - { - *fLog << all - << "MPad : Object '3D-PaddingMatrixMC' not found on root file" - << endl; - return kFALSE; - } - */ - fHgMC = new TH3D; - fHgMC->Read("3D-PaddingMatrixMC"); - *fLog << all - << "MPad : Object '3D-PaddingMatrixMC' was read in" << endl; - - /* - fHgON = - (TH3D*) gROOT->FindObject("3D-PaddingMatrixON"); - if (!fHgON) - { - *fLog << all - << "MPad : Object '3D-PaddingMatrixON' not found on root file" - << endl; - return kFALSE; - } - */ - fHgON = new TH3D; - fHgON->Read("3D-PaddingMatrixON"); - *fLog << all - << "MPad : Object '3D-PaddingMatrixON' was read in" << endl; - - /* - fHgOFF = - (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF"); - if (!fHgOFF) - { - *fLog << all - << "MPad : Object '3D-PaddingMatrixOFF' not found on root file" - << endl; - return kFALSE; - } - */ - fHgOFF = new TH3D; - fHgOFF->Read("3D-PaddingMatrixOFF"); - *fLog << all - << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl; - - //------------------------------------ - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Write padding distributions onto a file -// -// -Bool_t MPad::WritePaddingDist(const char* namefileout) -{ - *fLog << all << "MPad : Write padding histograms onto file " - << namefileout << endl; - - TFile outfile(namefileout, "RECREATE"); - - fHSigmaThetaMC->Write(); - fHSigmaThetaON->Write(); - fHSigmaThetaOFF->Write(); - - fHDiffPixTheta->Write(); - fHDiffPixThetaMC->Write(); - fHDiffPixThetaON->Write(); - fHDiffPixThetaOFF->Write(); - - fHSigmaPixTheta->Write(); - fHSigmaPixThetaMC->Write(); - fHSigmaPixThetaON->Write(); - fHSigmaPixThetaOFF->Write(); - - fHgMC->Write(); - fHgON->Write(); - fHgOFF->Write(); - - *fLog << all - << "MPad::WriteTargetDist(); padding histograms written onto file " - << namefileout << endl; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Set type of data to be padded -// -// this is not necessary if the type of the data can be recognized -// directly from the input -// -// -void MPad::SetDataType(const char* type) -{ - fType = type; - *fLog << all << "MPad::SetDataType(); type of data to be padded : " - << fType << endl; - - return; -} - -// -------------------------------------------------------------------------- -// -// Set the pointers and prepare the histograms -// -Int_t MPad::PreProcess(MParList *pList) -{ - if ( !fHSigmaThetaMC || !fHSigmaThetaON || !fHSigmaThetaOFF || - !fHDiffPixThetaMC || !fHDiffPixThetaON || !fHDiffPixThetaOFF || - !fHgMC || !fHgON || !fHgOFF ) - { - *fLog << err - << "MPad : At least one of the histograms needed for the padding is not defined ... aborting." - << endl; - return kFALSE; - } - - fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); - if (!fPointPos) - { - *fLog << err << dbginf << "MPad : MPointingPos not found... aborting." - << endl; - return kFALSE; - } - - fPed = (MPedPhotCam*)pList->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam"); - if (!fPed) - { - *fLog << err << AddSerialNumber(fNamePedPhotCam) - << "[MPedPhotCam] not found... aborting." << endl; - return kFALSE; - } - - - fCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fCam) - { - *fLog << err - << "MPad : MGeomCam not found (no geometry information available)... aborting." - << endl; - return kFALSE; - } - - fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); - if (!fEvt) - { - *fLog << err << "MPad : MCerPhotEvt not found... aborting." - << endl; - return kFALSE; - } - - - fBad = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam"); - if (!fBad) - { - *fLog << inf - << "MPad : MBadPixelsCam container not present... continue." - << endl; - } - - - if (fType !="ON" && fType !="OFF" && fType !="MC") - { - *fLog << err - << "MPad : Type of data to be padded not defined... aborting." - << endl; - return kFALSE; - } - - - //-------------------------------------------------------------------- - // histograms for checking the padding - // - // plot pedestal sigmas - fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", - 100, 0.0, 75.0, 100, 0.0, 75.0); - fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); - fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); - - // plot no.of photons (before vs. after padding) - fHPhotons = new TH2D("Photons","Photons: after vs.before padding", - 100, -100.0, 300.0, 100, -100, 300); - fHPhotons->SetXTitle("no.of photons before padding"); - fHPhotons->SetYTitle("no.of photons after padding"); - - // plot of added NSB - fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -100.0, 200.0); - fHNSB->SetXTitle("no.of added NSB photons"); - fHNSB->SetYTitle("no.of pixels"); - - - //-------------------------------------------------------------------- - - fIter = 10; - - memset(fErrors, 0, sizeof(fErrors)); - memset(fWarnings, 0, sizeof(fWarnings)); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Do the Padding (noise adjustment) -// -// input for the padding : -// - the matrices fHgON, fHgOFF, fHgMC -// - the original distributions fHSigmaTheta, fHDiffPixTheta -// - -Int_t MPad::Process() -{ - //*fLog << all << "Entry MPad::Process();" << endl; - - // this is the conversion factor from the number of photons - // to the number of photo electrons - // later to be taken from MCalibration - // fPEperPhoton = PW * LG * QE * 1D - Double_t fPEperPhoton = 0.92 * 0.94 * 0.23 * 0.9; - - //------------------------------------------------ - // Add pixels to MCerPhotEvt which are not yet in; - // set their number of photons equal to zero. - // Purpose : by the padding the number of photons is changed - // so that cleaning and cuts may be affected. - // However, no big effect is expectec unless the padding is strong - // (strong increase of the noise level) - // At present comment out this code - - /* - UInt_t imaxnumpix = fCam->GetNumPixels(); - - for (UInt_t i=0; iGetNumPixels(); - for (UInt_t j=0; jAddPixel(i, 0.0, (*fPed)[i].GetRms()); - } - */ - - - //----------------------------------------- - Int_t rc=0; - Int_t rw=0; - - const UInt_t npix = fEvt->GetNumPixels(); - - //------------------------------------------- - // Calculate average sigma of the event - // - - //*fLog << all << "MPad::Process(); before padding : " << endl; - - // inner sigmabar/sqrt(area) - Double_t ratioinn = fCam->GetPixRatio(0); - Double_t sigbarInnerold_phot = (fPed->GetArea(0)).GetRms() * sqrt(ratioinn); - Double_t sigbarInnerold = sigbarInnerold_phot * fPEperPhoton; - Double_t sigbarInnerold2 = sigbarInnerold*sigbarInnerold; - - // outer sigmabar/sqrt(area) - Double_t ratioout = fCam->GetPixRatio(500); - Double_t sigbarOuterold_phot = (fPed->GetArea(1)).GetRms() * sqrt(ratioout); - Double_t sigbarOuterold = sigbarOuterold_phot * fPEperPhoton; - Double_t sigbarOuterold2 = sigbarOuterold*sigbarOuterold; - - const Double_t theta = fPointPos->GetZd(); - - //*fLog << all << "theta = " << theta << endl; - - Int_t binTheta = fHSigmaThetaON->GetXaxis()->FindBin(theta); - if ( binTheta < 1 || binTheta > fHSigmaThetaON->GetNbinsX() ) - { - *fLog << warn - << "MPad::Process(); binNumber out of range : theta, binTheta = " - << theta << ", " << binTheta << endl; - - rc = 2; - fErrors[rc]++; - return kCONTINUE; - } - - //*fLog << all << "binTheta = " << binTheta << endl; - - - //------------------------------------------- - // get number of events in this theta bin - // and number of events in this sigbarInnerold_phot bin - - Double_t nTheta; - Double_t nSigma; - - TH2D *st = NULL; - if (fType == "MC") st = fHSigmaThetaMC; - else if (fType == "ON") st = fHSigmaThetaON; - else if (fType == "OFF") st = fHSigmaThetaOFF; - else - { - *fLog << err << "MPad : illegal data type '" << fType - << "'" << endl; - return kERROR; - } - - TH1D *hn; - hn = st->ProjectionY("", binTheta, binTheta, ""); - - //*fLog << "Title of histogram : " << st->GetTitle() << endl; - //for (Int_t i=1; i<=hn->GetNbinsX(); i++) - //{ - // *fLog << "bin " << i << " : no.of entries = " << hn->GetBinContent(i) - // << endl; - //} - - nTheta = hn->Integral(); - Int_t binSigma = hn->FindBin(sigbarInnerold_phot); - nSigma = hn->GetBinContent(binSigma); - - //*fLog << all - // << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " - // << theta << ", " << sigbarInnerold_phot << ", " << binTheta - // << ", " - // << binSigma << ", " << nTheta << ", " << nSigma << endl; - - delete hn; - - //------------------------------------------- - - //****************************************************************** - // has the event to be padded ? - // if yes : to which sigmabar should it be padded ? - // - - Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarInnerold_phot); - //*fLog << all << "binSig, sigbarInnerold_phot = " << binSig << ", " - // << sigbarInnerold_phot << endl; - - Double_t prob; - TH1D *hpad = NULL; - - TH3D *Hg = NULL; - if (fType == "MC") Hg = fHgMC; - else if (fType == "ON") Hg = fHgON; - else if (fType == "OFF") Hg = fHgOFF; - else - { - *fLog << err << "MPad : illegal type of data '" << fType << "'" - << endl; - return kERROR; - } - - hpad = Hg->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, ""); - - //Int_t nb = hpad->GetNbinsX(); - //for (Int_t i=1; i<=nb; i++) - //{ - // if (hpad->GetBinContent(i) != 0.0) - // *fLog << all << "i, fHgON = " << i << ", " - // << hpad->GetBinContent(i) << endl; - //} - - if (nSigma != 0.0) - prob = hpad->Integral() * nTheta/nSigma; - else - prob = 0.0; - - //*fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma - // << ", " << prob << endl; - - if ( prob <= 0.0 || gRandom->Uniform() > prob ) - { - delete hpad; - // event should not be padded - //*fLog << all << "event is not padded" << endl; - - rc = 8; - fErrors[rc]++; - return kTRUE; - } - // event should be padded - //*fLog << all << "event will be padded" << endl; - - - //------------------------------------------- - // for the current theta, generate a sigmabar : - // for MC/ON/OFF data according to the matrix fHgMC/ON/OFF - // - Double_t sigmabarInner_phot = 0; - Double_t sigmabarInner = 0; - - - //Int_t nbinsx = hpad->GetNbinsX(); - //for (Int_t i=1; i<=nbinsx; i++) - //{ - // if (hpad->GetBinContent(i) != 0.0) - // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i) - // << endl; - //} - - sigmabarInner_phot = hpad->GetRandom(); - sigmabarInner = sigmabarInner_phot * fPEperPhoton; - - //*fLog << all << "sigmabarInner_phot = " << sigmabarInner_phot << endl; - - delete hpad; - - // new inner sigmabar2/area - const Double_t sigmabarInner2 = sigmabarInner*sigmabarInner; - - //------------------------------------------- - - //*fLog << all << "MPad::Process(); sigbarInnerold, sigmabarInner = " - // << sigbarInnerold << ", "<< sigmabarInner << endl; - - // Stop if target sigmabar is <= sigbarold - if (sigmabarInner <= sigbarInnerold) - { - *fLog << all << "MPad::Process(); target sigmabarInner is less than sigbarInnerold : " - << sigmabarInner << ", " << sigbarInnerold << ", aborting" - << endl; - - rc = 4; - fErrors[rc]++; - return kCONTINUE; - } - - - //------------------------------------------- - // - // Calculate average number of NSB photo electrons to be added (lambdabar) - // from the value of sigmabar, - // - making assumptions about the average electronic noise (elNoise2) and - // - using a fixed value (F2excess) for the excess noise factor - - Double_t elNoise2; // [photo electrons] - Double_t F2excess = 1.3; - Double_t lambdabar; // [photo electrons] - - - Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta); - if (binTheta != bincheck) - { - *fLog << err - << "MPad::Process(); The binnings of the 2 histograms are not identical; aborting" - << endl; - return kERROR; - } - - // In this Theta bin, get RMS of (Sigma^2-sigmabar^2) for inner pixels - // The average electronic noise (to be added) has to be in the order of - // this RMS, otherwise the electronic noise of an individual pixel - // (elNoise2Pix) may become negative - - // Attention : maximum ID of inner pixels hard coded !!! - Int_t idmaxpixinner = 396; - Int_t binmaxpixinner = - fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)idmaxpixinner ); - - TH1D *hnoise; - if (fType == "MC") - hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); - else if (fType == "ON") - hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); - else if (fType == "OFF") - hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); - else - { - *fLog << err << "MPad::Process(); illegal data type... aborting" - << endl; - return kERROR; - } - - Double_t RMS_phot = hnoise->GetRMS(1); - Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton; - delete hnoise; - - elNoise2 = TMath::Min(2.0*RMS, sigmabarInner2 - sigbarInnerold2); - //*fLog << all << "elNoise2 = " << elNoise2 << endl; - - lambdabar = (sigmabarInner2 - sigbarInnerold2 - elNoise2) / F2excess; - - // This value of lambdabar is the same for all pixels; - // note that lambdabar is normalized to the area of pixel 0 - - //---------- start loop over pixels --------------------------------- - // do the padding for each pixel - // - // pad only pixels - which are used and not bad (before image cleaning) - // - - Double_t sigma2 = 0; - - Double_t diff_phot = 0; - Double_t diff = 0; - - Double_t addSig2_phot= 0; - Double_t addSig2 = 0; - - Double_t elNoise2Pix = 0; - - - // throw the Sigma for the pixels from the distribution fHDiffPixTheta - // MC : from fHDiffPixTheta - // ON : from fHDiffPixThetaOFF - // OFF : from fHDiffPixThetaON - - TH3D *sp = NULL; - - if (fType == "MC") sp = fHDiffPixTheta; - else if (fType == "ON") sp = fHDiffPixThetaOFF; - else if (fType == "OFF") sp = fHDiffPixThetaON; - else - { - *fLog << err << "MPad::Process(); illegal data type... aborting" - << endl; - return kERROR; - } - - Double_t sigbarold2; - Double_t sigmabar2; - Double_t sigmabarOuter2; - - - for (UInt_t i=0; iGetPixRatio(j); - - if (ratioArea < 1.5) - { - sigbarold2 = sigbarInnerold2; - sigmabar2 = sigmabarInner2; - } - else - { - sigbarold2 = sigbarOuterold2; - //sigbarOuter2 = sigmabarInner2 * sigbarOuterold2 / sigbarInnerold2; - sigmabarOuter2 = sigmabarInner2 + (sigbarOuterold2 - sigbarInnerold2); - sigmabar2 = sigmabarOuter2; - } - - MPedPhotPix &ppix = (*fPed)[j]; - - if ( fBad != NULL && (*fBad)[j].IsUnsuitable() ) - { - // this should never occur, because bad pixels should have - // been set unused - *fLog << all << "MPad::Process; bad pixel found which is used, j = " - << j << "... go to next pixel." << endl; - continue; - } - - // count number of pixels treated - fWarnings[0]++; - - - Double_t oldsigma_phot = ppix.GetRms(); - Double_t oldsigma = oldsigma_phot * fPEperPhoton; - Double_t oldsigma2 = oldsigma*oldsigma / ratioArea; - - //--------------------------------- - // throw the Sigma for this pixel - // - Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j ); - - Int_t count; - Bool_t ok; - - TH1D *hdiff; - - - hdiff = sp->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - Double_t integral = hdiff->Integral(); - // if there are no entries in hdiff, diff cannot be thrown - // in this case diff will be set to zero - if ( integral == 0 ) - { - //*fLog << warn << "MPad::Process(); fType = " << fType - // << ", projection of '" - // << sp->GetName() << "' for Theta bin " - // << binTheta << " and pixel " << j - // << " has no entries; set diff equal to the old difference " - // << endl; - - diff = TMath::Max(oldsigma2 - sigbarold2, - lambdabar*F2excess + oldsigma2 - sigmabar2); - - rc = 2; - fWarnings[rc]++; - } - - // start of else ------------------- - else - { - count = 0; - ok = kFALSE; - for (Int_t m=0; mGetRandom(); - - //*fLog << "after GetRandom : j, m, diff_phot = " << j << " : " - // << m << ", " << diff_phot << endl; - - diff = diff_phot * fPEperPhoton * fPEperPhoton; - - // the following condition ensures that elNoise2Pix > 0.0 - if ( (diff + sigmabar2 - oldsigma2 - - lambdabar*F2excess) > 0.0 ) - { - ok = kTRUE; - break; - } - } - - if (!ok) - { - //*fLog << all << "theta, j, count, sigmabar2, diff, oldsigma2, ratioArea, lambdabar = " - // << theta << ", " - // << j << ", " << count << ", " << sigmabar2 << ", " - // << diff << ", " << oldsigma2 << ", " << ratioArea << ", " - // << lambdabar << endl; - diff = lambdabar*F2excess + oldsigma2 - sigmabar2; - - rw = 1; - fWarnings[rw]++; - } - } - // end of else -------------------- - - delete hdiff; - sigma2 = diff + sigmabar2; - - - //--------------------------------- - // get the additional sigma^2 for this pixel (due to the padding) - - addSig2 = (sigma2 - oldsigma2) * ratioArea; - addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); - - //--------------------------------- - // get the additional electronic noise for this pixel - - elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea; - - - //--------------------------------- - // throw actual number of additional NSB photo electrons (NSB) - // and its RMS (sigmaNSB) - - Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea); - Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix; - Double_t sigmaNSB0; - - if (arg >= 0.0) - { - sigmaNSB0 = sqrt( arg ); - } - else - { - sigmaNSB0 = 0.0000001; - if (arg < -1.e-10) - { - *fLog << warn << "MPad::Process(); argument of sqrt < 0 : " - << arg << endl; - } - } - - - //--------------------------------- - // smear NSB0 according to sigmaNSB0 - // and subtract lambdabar because of AC coupling - - Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea; - Double_t NSB_phot = NSB / fPEperPhoton; - - //--------------------------------- - - // add additional NSB to the number of photons - Double_t oldphotons_phot = pix.GetNumPhotons(); - Double_t oldphotons = oldphotons_phot * fPEperPhoton; - Double_t newphotons = oldphotons + NSB; - Double_t newphotons_phot = newphotons / fPEperPhoton; - pix.SetNumPhotons( newphotons_phot ); - - - fHNSB->Fill( NSB_phot/ratioArea ); - fHPhotons->Fill( oldphotons_phot/ratioArea, - newphotons_phot/ratioArea ); - - - // error: add sigma of padded noise quadratically - Double_t olderror_phot = pix.GetErrorPhot(); - Double_t newerror_phot = - sqrt( olderror_phot*olderror_phot + addSig2_phot ); - pix.SetErrorPhot( newerror_phot ); - - - Double_t newsigma = sqrt( sigma2 * ratioArea ); - Double_t newsigma_phot = newsigma / fPEperPhoton; - ppix.SetRms( newsigma_phot ); - - fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot ); - } - //---------- end of loop over pixels --------------------------------- - - // Calculate sigmabar again and crosscheck - - //*fLog << all << "MPad::Process(); after padding : " << endl; - - //*fLog << all << "Exit MPad::Process();" << endl; - - rc = 0; - fErrors[rc]++; - return kTRUE; - //****************************************************************** -} - -// -------------------------------------------------------------------------- -// -// -Int_t MPad::PostProcess() -{ - if (GetNumExecutions() != 0) - { - - *fLog << inf << endl; - *fLog << GetDescriptor() << " execution statistics:" << endl; - *fLog << dec << setfill(' '); - - if (fWarnings[0] == 0 ) fWarnings[0] = 1; - - *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) - << (int)(fWarnings[1]*100/fWarnings[0]) - << "%) Pixel: iteration to find acceptable sigma failed" - << ", fIter = " << fIter << endl; - - *fLog << " " << setw(7) << fWarnings[2] << " (" << setw(3) - << (int)(fWarnings[2]*100/fWarnings[0]) - << "%) Pixel: No data for generating Sigma^2-Sigmabar^2" << endl; - - - *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) - << (int)(fErrors[2]*100/GetNumExecutions()) - << "%) Evts skipped due to: Zenith angle out of range" << endl; - - *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) - << (int)(fErrors[4]*100/GetNumExecutions()) - << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; - - *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3) - << (int)(fErrors[8]*100/GetNumExecutions()) - << "%) Evts didn't have to be padded" << endl; - - *fLog << " " << fErrors[0] << " (" - << (int)(fErrors[0]*100/GetNumExecutions()) - << "%) Evts were successfully padded!" << endl; - *fLog << endl; - - } - - //--------------------------------------------------------------- - TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 600)); - c.Divide(3, 2); - gROOT->SetSelectedPad(NULL); - - c.cd(1); - fHNSB->SetDirectory(NULL); - fHNSB->DrawCopy(); - fHNSB->SetBit(kCanDelete); - - c.cd(2); - fHSigmaPedestal->SetDirectory(NULL); - fHSigmaPedestal->DrawCopy(); - fHSigmaPedestal->SetBit(kCanDelete); - - c.cd(3); - fHPhotons->SetDirectory(NULL); - fHPhotons->DrawCopy(); - fHPhotons->SetBit(kCanDelete); - - //-------------------------------------------------------------------- - - /* - c.cd(4); - fHSigmaTheta->SetDirectory(NULL); - fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta"); - fHSigmaTheta->DrawCopy(); - fHSigmaTheta->SetBit(kCanDelete); - */ - - //-------------------------------------------------------------------- - // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2 - - /* - c.cd(5); - TH2D *l1 = NULL; - l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx"); - l1->SetDirectory(NULL); - l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); - l1->SetXTitle("\\Theta [\\circ]"); - l1->SetYTitle("Sigma^2-Sigmabar^2"); - - l1->DrawCopy("box"); - l1->SetBit(kCanDelete);; - - c.cd(8); - TH2D *l2 = NULL; - l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy"); - l2->SetDirectory(NULL); - l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); - l2->SetXTitle("pixel"); - l2->SetYTitle("Sigma^2-Sigmabar^2"); - - l2->DrawCopy("box"); - l2->SetBit(kCanDelete);; - */ - - //-------------------------------------------------------------------- - // draw the 3D histogram (target): Theta, pixel, Sigma - - /* - c.cd(6); - TH2D *k1 = NULL; - k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx"); - k1->SetDirectory(NULL); - k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)"); - k1->SetXTitle("\\Theta [\\circ]"); - k1->SetYTitle("Sigma"); - - k1->DrawCopy("box"); - k1->SetBit(kCanDelete); - - c.cd(9); - TH2D *k2 = NULL; - k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy"); - k2->SetDirectory(NULL); - k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)"); - k2->SetXTitle("pixel"); - k2->SetYTitle("Sigma"); - - k2->DrawCopy("box"); - k2->SetBit(kCanDelete);; - */ - - //-------------------------------------------------------------------- - - - c.cd(4); - TH2D *m1; - m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); - m1->SetDirectory(NULL); - m1->SetTitle("(fHgON) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); - m1->SetXTitle("Sigmabar-old"); - m1->SetYTitle("Sigmabar-new"); - - m1->DrawCopy("box"); - m1->SetBit(kCanDelete);; - - c.cd(5); - TH2D *m2; - m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); - m2->SetDirectory(NULL); - m2->SetTitle("(fHgOFF) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); - m2->SetXTitle("Sigmabar-old"); - m2->SetYTitle("Sigmabar-new"); - - m2->DrawCopy("box"); - m2->SetBit(kCanDelete);; - - c.cd(6); - TH2D *m3; - m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy"); - m3->SetDirectory(NULL); - m3->SetTitle("(fHgMC) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)"); - m3->SetXTitle("Sigmabar-old"); - m3->SetYTitle("Sigmabar-new"); - - m3->DrawCopy("box"); - m3->SetBit(kCanDelete);; - - - return kTRUE; -} - -// -------------------------------------------------------------------------- - - - - - - - - - - Index: unk/MagicSoft/Mars/manalysis/MPad.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MPad.h (revision 5430) +++ (revision ) @@ -1,121 +1,0 @@ -#ifndef MARS_MPad -#define MARS_MPad - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH1D; -class TH2D; -class TH3D; - -class MGeomCam; -class MCerPhotEvt; -class MPedPhotCam; -class MPointingPos; -class MSigmabar; -class MParList; -class MBadPixelsCam; -class MRead; -class MFilterList; - - -class MPad : public MTask -{ -private: - MGeomCam *fCam; - MCerPhotEvt *fEvt; - MSigmabar *fSigmabar; - MPointingPos *fPointPos; - MPedPhotCam *fPed; - MBadPixelsCam *fBad; - - TString fNamePedPhotCam; // name of the 'MPedPhotCam' container - TString fType; // type of data to be padded - TFile *fInfile; // input file containing padding histograms - - Int_t fIter; - - Int_t fErrors[9]; - Int_t fWarnings[3]; - - //---------------------------------- - // plots used for the padding - // for all plots it is assumed that the pedestal RMS is given in units of "number of photons" - - TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) - TH2D *fHSigmaThetaMC; // 2D-histogram (sigmabar vs. Theta) - TH2D *fHSigmaThetaON; // 2D-histogram (sigmabar vs. Theta) - TH2D *fHSigmaThetaOFF; // 2D-histogram (sigmabar vs. Theta) - - //--------------------- - TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - TH3D *fHDiffPixThetaMC; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - TH3D *fHDiffPixThetaON; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - TH3D *fHDiffPixThetaOFF; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - - //--------------------- - TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHSigmaPixThetaMC; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHSigmaPixThetaON; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHSigmaPixThetaOFF; // 3D-histogram (Theta, pixel, sigma) - - //--------------------- - TH3D *fHgMC; // matrix (Theta, sigbarold, sigbarnew) for MC data - TH3D *fHgON; // matrix (Theta, sigbarold, sigbarnew) for ON data - TH3D *fHgOFF; // matrix (Theta, sigbarold, sigbarnew) for OFF data - - //------------------------------- - // plots for checking the padding - TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after - // versus before padding - TH2D *fHPhotons; // 2D-histogram : no.of photons after - // versus before padding - TH1D *fHNSB; // 1D-histogram : additional NSB - - - Bool_t Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap, - TH2D *fHga, TH2D *fHgb, Int_t nbinssig); - - Bool_t UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA, - Int_t nbinssig, Int_t l); - -public: - MPad(const char *name=NULL, const char *title=NULL); - ~MPad(); - - void SetNamePedPhotCam(const char *name) { fNamePedPhotCam = name; } - - Bool_t MergeONOFFMC( - TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, - TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, - TH2D& sigthoff, TH3D& diffpixthoff,TH3D& sigmapixthoff); - - Bool_t MergeONMC( - TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, - TH2D& sigthon, TH3D&diffpixthon, TH3D& sigmapixthon); - - Bool_t ReadPaddingDist(const char *filein); - Bool_t WritePaddingDist(const char *fileout); - - void SetDataType(const char *type); // type of data to be padded - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - ClassDef(MPad, 0) // task for the padding -}; - -#endif - - - - - - - Index: /trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- /trunk/MagicSoft/Mars/manalysis/Makefile (revision 5430) +++ /trunk/MagicSoft/Mars/manalysis/Makefile (revision 5431) @@ -23,5 +23,5 @@ -I../mhistmc -I../mfileio -I../mmain -I../mhcalib -I../mcalib \ -I../msignal -I../mpointing -I../mtools -I../mfbase \ - -I../mbadpixels -I../mastro -I../mpedestal + -I../mbadpixels -I../mastro -I../mpedestal -I../mtemp/mucm/classes # mcalib: MGeomApply (MCalibrationCam) @@ -55,5 +55,6 @@ MFiltercutsCalc.cc \ MPad.cc \ - MMcCalibrationUpdate.cc + MMcCalibrationUpdate.cc \ + MMakePadHistograms.cc ############################################################ @@ -66,2 +67,3 @@ mrproper: clean rmbak + Index: /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc (revision 5430) +++ /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc (revision 5431) @@ -106,5 +106,5 @@ if (!fPedPhotCam) { - *fLog << err << "MPedPhotCam not found... aborting." << endl; + *fLog << err << fNamePedPhotContainer << "[MPedPhotCam] not found... aborting." << endl; return kFALSE; } @@ -117,4 +117,8 @@ } } + + *fLog << inf << "Name of MPedPhotCam container : " << fNamePedPhotContainer + << endl; + return kTRUE; } @@ -257,4 +261,5 @@ fPedPhotCam->ReCalc(*fGeomCam, fBadPixels); } + return kTRUE; } Index: /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h (revision 5430) +++ /trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h (revision 5431) @@ -37,2 +37,6 @@ #endif + + + + Index: /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.cc =================================================================== --- /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.cc (revision 5430) +++ /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.cc (revision 5431) @@ -35,4 +35,5 @@ #include "MBadPixelsCam.h" #include "MGeomCam.h" +#include "MPedPhotCam.h" #include "MParList.h" #include "MBinning.h" @@ -47,24 +48,55 @@ // ------------------------------------------------------------------------- // -// Default Constructor. +// Constructor. // MHBadPixels::MHBadPixels(const char *name, const char *title) - : fNamePedPhotCam("MPedPhotCam") { fName = name ? name : "MHBadPixels"; fTitle = title ? title : "Histogram for Bad Pixels vs. Theta"; - fBadId.SetName("2D-IdBadPixels"); - fBadId.SetTitle("2D-IdBadPixels"); - fBadId.SetDirectory(NULL); - fBadId.SetXTitle("\\Theta [\\circ]"); - fBadId.SetYTitle("pixel Id"); - - fBadN.SetName("2D-NBadPixels"); - fBadN.SetTitle("2D-NBadPixels"); - fBadN.SetDirectory(NULL); - fBadN.SetXTitle("\\Theta [\\circ]"); - fBadN.SetYTitle("number of bad pixels"); - + fBadId = new TH2D; + fBadId->SetName("2D-IdBadPixels"); + fBadId->SetTitle("2D-IdBadPixels"); + fBadId->SetDirectory(NULL); + fBadId->UseCurrentStyle(); + fBadId->SetXTitle("\\Theta [\\circ]"); + fBadId->SetYTitle("pixel Id"); + fBadId->SetTitleOffset(1.2, "Y"); + + fBadN = new TH2D; + fBadN->SetName("2D-NBadPixels"); + fBadN->SetTitle("2D-NBadPixels"); + fBadN->SetDirectory(NULL); + fBadN->UseCurrentStyle(); + fBadN->SetXTitle("\\Theta [\\circ]"); + fBadN->SetYTitle("number of bad pixels"); + fBadN->SetTitleOffset(1.2, "Y"); + + // define default binnings + fBinsTheta = new MBinning; + fBinsTheta->SetEdgesCos(10, 0.0, 90.0); // theta + + fBinsPix = new MBinning; + fBinsPix->SetEdges(578, -0.5, 577.5); // pixel id + + SetBinning(fBadId, fBinsTheta, fBinsPix); + SetBinning(fBadN, fBinsTheta, fBinsPix); + + //----------------------------------------- + fNamePedPhotCam = "MPedPhotCamFromData"; +} + + +// ------------------------------------------------------------------------- +// +// Destructor. +// +MHBadPixels::~MHBadPixels() +{ + delete fBadId; + delete fBadN; + + delete fBinsTheta; + delete fBinsPix; } @@ -75,10 +107,12 @@ Bool_t MHBadPixels::SetupFill(const MParList *plist) { - MGeomCam *fCam = (MGeomCam*)plist->FindObject(AddSerialNumber("MGeomCam")); + fCam = (MGeomCam*)plist->FindObject(AddSerialNumber("MGeomCam")); if (!fCam) { - *fLog << err << "MGeomCam not found... aborting." << endl; - return kFALSE; - } + *fLog << err << "MHBadPixels::SetupFill; MGeomCam not found... aborting." << endl; + return kFALSE; + } + *fLog << "MHBadPixels::SetupFill; fCam = " << fCam << endl; + fPointPos = (MPointingPos*)plist->FindObject("MPointingPos"); if (!fPointPos) @@ -87,23 +121,44 @@ return kFALSE; } + + fPedPhot = (MPedPhotCam*)plist->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam"); + if (!fPedPhot) + { + *fLog << err << AddSerialNumber(fNamePedPhotCam) + << "[MPedPhotCam] not found... aborting." << endl; + return kFALSE; + } + fPedPhot->InitSize(fCam->GetNumPixels()); + //---------------------------------------------------- + *fLog << inf << "Name of MPedPhotCam container : " << fNamePedPhotCam + << endl; + + + //---------------------------------------------------- + // redefine the binnings + // Get Theta Binning - MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning"); + const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning"); if (!binstheta) { - *fLog << err << "BinningTheta [MBinning] not found... aborting" << endl; - return kFALSE; - } - - // Get binning for pixel number - const UInt_t npix1 = fCam->GetNumPixels()+1; - + *fLog << err << "Object 'BinningTheta' [MBinning] not found... use default binning." << endl; + binstheta = fBinsTheta; + } + + // Define binning for pixel number + const UInt_t npix1 = fPedPhot->GetSize()+1; + //*fLog << "MHBadPixels::SetupFill(); npix1 = " << npix1 << endl; MBinning binspix("BinningPixel"); binspix.SetEdges(npix1, -0.5, npix1-0.5); // Set binnings in histograms - SetBinning(&fBadId, binstheta, &binspix); - SetBinning(&fBadN, binstheta, &binspix); + SetBinning(fBadId, binstheta, &binspix); + SetBinning(fBadN, binstheta, &binspix); + + *fLog << "MHBadPixels::SetupFill(); binnings were set" << endl; + + return kTRUE; @@ -122,18 +177,21 @@ Double_t theta = fPointPos->GetZd(); - const MBadPixelsCam *badpixels = (MBadPixelsCam*)par; - - const UInt_t entries = badpixels->GetSize(); - + const MBadPixelsCam *fBadPixels = (MBadPixelsCam*)par; + + const UInt_t entries = fPedPhot->GetSize(); UInt_t nb = 0; for (UInt_t i=0; iFill(theta, i, w); + nb++; + } + } + fBadN->Fill(theta, nb, w); return kTRUE; @@ -154,16 +212,17 @@ pad->cd(1); - fBadId.Draw(option); + fBadId->Draw(option); pad->cd(2); - fBadN.Draw(option); + fBadN->Draw(option); pad->cd(3); gPad->SetBorderMode(0); - h = ((TH2*)&fBadId)->ProjectionY("ProjY-pixId", -1, 9999, ""); + h = ((TH2*)fBadId)->ProjectionY("ProjY-pixId", -1, 9999, ""); h->SetDirectory(NULL); h->SetTitle("Distribution of bad pixel Id"); h->SetXTitle("Id of bad pixel"); h->SetYTitle("No. of events"); + h->SetTitleOffset(1.2, "Y"); h->Draw(option); h->SetBit(kCanDelete); @@ -171,9 +230,10 @@ pad->cd(4); gPad->SetBorderMode(0); - h = ((TH2*)&fBadN)->ProjectionY("ProjY-pixN", -1, 9999, ""); + h = ((TH2*)fBadN)->ProjectionY("ProjY-pixN", -1, 9999, ""); h->SetDirectory(NULL); h->SetTitle("Distribution of no.of bad pixels"); h->SetXTitle("No. of bad pixels"); h->SetYTitle("No. of events"); + h->SetTitleOffset(1.2, "Y"); h->Draw(option); h->SetBit(kCanDelete); @@ -182,2 +242,14 @@ pad->Update(); } +//========================================================================== + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.h =================================================================== --- /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.h (revision 5430) +++ /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.h (revision 5431) @@ -18,21 +18,28 @@ private: MGeomCam *fCam; //! + MPedPhotCam *fPedPhot; //! MPointingPos *fPointPos; //! TString fNamePedPhotCam; //! name of the 'MPedPhotCam' container + MBinning *fBinsTheta; //! + MBinning *fBinsPix; //! - TH2D fBadId; // 2D-histogram : pixel Id vs. Theta - TH2D fBadN; // 2D-histogram : no.of bad pixels vs. Theta + TH2D *fBadId; // 2D-histogram : pixel Id vs. Theta + TH2D *fBadN; // 2D-histogram : no.of bad pixels vs. Theta public: MHBadPixels(const char *name=NULL, const char *title=NULL); + ~MHBadPixels(); void SetNamePedPhotCam(const char *name) { fNamePedPhotCam = name; } - const TH2D *GetBadId() { return &fBadId; } - const TH2D *GetBadId() const { return &fBadId; } + const TH2D *GetBadId() { return fBadId; } + const TH2D *GetBadId() const { return fBadId; } - const TH2D *GetBadN() { return &fBadN; } - const TH2D *GetBadN() const { return &fBadN; } + const TH2D *GetBadN() { return fBadN; } + const TH2D *GetBadN() const { return fBadN; } + + TH2 *GetBadIdByName(const TString name) { return fBadId; } + TH2 *GetBadNByName(const TString name) { return fBadN; } void Draw(Option_t* option = ""); Index: /trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc (revision 5430) +++ /trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc (revision 5431) @@ -86,44 +86,47 @@ fThetaMax = thetamax; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers +// +// +Int_t MFSelBasic::PreProcess(MParList *pList) +{ + fRawRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRawRun) + { + *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); + if (!fPointPos) + { + *fLog << dbginf << "MPointingPos not found... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + + memset(fCut, 0, sizeof(fCut)); + + //------------------------- *fLog << inf << "MFSelBasic cut values : fMinPhotons, fThetaMin, fThetaMax = "; *fLog << fMinPhotons <<", " << fThetaMin << ", " << fThetaMax << endl; -} - -// -------------------------------------------------------------------------- -// -// Set the pointers -// -// -Int_t MFSelBasic::PreProcess(MParList *pList) -{ - fRawRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!fRawRun) - { - *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; - return kFALSE; - } - - fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); - if (!fPointPos) - { - *fLog << dbginf << "MPointingPos not found... aborting." << endl; - return kFALSE; - } - - fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); - if (!fEvt) - { - *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; - return kFALSE; - } - - fCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fCam) - { - *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; - return kFALSE; - } - - memset(fCut, 0, sizeof(fCut)); + //------------------------- return kTRUE; Index: /trunk/MagicSoft/Mars/mfilter/MFSelFinal.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFSelFinal.cc (revision 5430) +++ /trunk/MagicSoft/Mars/mfilter/MFSelFinal.cc (revision 5431) @@ -86,6 +86,4 @@ fDistMax = distmax; - *fLog << inf << "MFSelFinal cut values : fHadronnessMax, fAlphaMax, fDistMax = "; - *fLog << fHadronnessMax << ", " << fAlphaMax << ", " << fDistMax << endl; } @@ -120,4 +118,9 @@ memset(fCut, 0, sizeof(fCut)); + + //---------------------- + *fLog << inf << "MFSelFinal cut values : fHadronnessMax, fAlphaMax, fDistMax = "; + *fLog << fHadronnessMax << ", " << fAlphaMax << ", " << fDistMax << endl; + //---------------------- return kTRUE; Index: /trunk/MagicSoft/Mars/mfilter/MFSelStandard.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFSelStandard.cc (revision 5430) +++ /trunk/MagicSoft/Mars/mfilter/MFSelStandard.cc (revision 5431) @@ -92,4 +92,45 @@ fWidthMin = widthmin; +} + +// -------------------------------------------------------------------------- +// +// MISSING +// +Int_t MFSelStandard::PreProcess(MParList *pList) +{ + fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); + if (!fHil) + { + *fLog << err << fHilName << " [MHillas] not found... aborting." << endl; + return kFALSE; + } + + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); + if (!fHilSrc) + { + *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; + return kFALSE; + } + + fNewImgPar = (MNewImagePar*)pList->FindObject(fImgParName, "MNewImagePar"); + if (!fNewImgPar) + { + *fLog << err << fImgParName << " [MNewImagePar] not found... aborting." << endl; + return kFALSE; + } + + MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!cam) + { + *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = cam->GetConvMm2Deg(); + + memset(fCut, 0, sizeof(fCut)); + + //-------------------- *fLog << inf << "MFSelStandard cut values : fUsedPixelsMax, fCorePixelsMin = "; *fLog << fUsedPixelsMax << ", " << fCorePixelsMin << endl; @@ -98,43 +139,5 @@ *fLog << inf << " fLengthMin, fWidthMin = " << fLengthMin ; *fLog << ", " << fWidthMin << endl; -} - -// -------------------------------------------------------------------------- -// -// MISSING -// -Int_t MFSelStandard::PreProcess(MParList *pList) -{ - fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); - if (!fHil) - { - *fLog << err << fHilName << " [MHillas] not found... aborting." << endl; - return kFALSE; - } - - fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); - if (!fHilSrc) - { - *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; - return kFALSE; - } - - fNewImgPar = (MNewImagePar*)pList->FindObject(fImgParName, "MNewImagePar"); - if (!fNewImgPar) - { - *fLog << err << fImgParName << " [MNewImagePar] not found... aborting." << endl; - return kFALSE; - } - - MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!cam) - { - *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; - return kFALSE; - } - - fMm2Deg = cam->GetConvMm2Deg(); - - memset(fCut, 0, sizeof(fCut)); + //-------------------- return kTRUE; Index: /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 5430) +++ /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 5431) @@ -32,4 +32,5 @@ // - the 3D-histogram sigma, pixel no., Theta // // - the 3D-histogram (sigma^2-sigmabar^2), pixel no., Theta // +// - the 2D-histogram Theta vs. Phi // // // ////////////////////////////////////////////////////////////////////////////// @@ -64,5 +65,5 @@ // -------------------------------------------------------------------------- // -// Default Constructor. It sets name and title of the histogram. +// Constructor. // MHSigmaTheta::MHSigmaTheta(const char *name, const char *title) @@ -71,49 +72,104 @@ fTitle = title ? title : "2D histogram sigmabar vs. Theta"; - fSigmaTheta.SetDirectory(NULL); - fSigmaTheta.SetName("2D-ThetaSigmabar(Inner)"); - fSigmaTheta.SetTitle("2D: \\bar{\\sigma}, \\Theta"); - fSigmaTheta.SetXTitle("\\Theta [\\circ]"); - fSigmaTheta.SetYTitle("Sigmabar(Inner) / SQRT(Area)"); - - fSigmaThetaOuter.SetDirectory(NULL); - fSigmaThetaOuter.SetName("2D-ThetaSigmabar(Outer)"); - fSigmaThetaOuter.SetTitle("2D: \\bar{\\sigma}, \\Theta"); - fSigmaThetaOuter.SetXTitle("\\Theta [\\circ]"); - fSigmaThetaOuter.SetYTitle("Sigmabar(Outer) / SQRT(Area)"); - - fSigmaPixTheta.SetDirectory(NULL); - fSigmaPixTheta.SetName("3D-ThetaPixSigma"); - fSigmaPixTheta.SetTitle("3D: \\Theta, Pixel Id, \\sigma"); - fSigmaPixTheta.SetXTitle("\\Theta [\\circ]"); - fSigmaPixTheta.SetYTitle("Pixel Id"); - fSigmaPixTheta.SetZTitle("Sigma"); - - fDiffPixTheta.SetDirectory(NULL); - fDiffPixTheta.SetName("3D-ThetaPixDiff"); - fDiffPixTheta.SetTitle("3D: \\Theta, Pixel Id, {\\sigma}^{2}-\\bar{\\sigma}^{2}"); - fDiffPixTheta.SetXTitle("\\Theta [\\circ]"); - fDiffPixTheta.SetYTitle("Pixel Id"); - fDiffPixTheta.SetZTitle("(Sigma2 - Sigmabar2)/Area"); - - // Set default binning - // FIXME: Maybe ist's necessary to adapt the value to the - // Magic default values - MBinning binsb; - MBinning binst; - MBinning binsd; - binsd.SetEdges(100, -10, 20); - binsb.SetEdges(100, 0, 10); - binst.SetEdgesCos(10, 0, 90); - - MBinning binspix("BinningPixel"); - binspix.SetEdges(578, -0.5, 577.5); - - SetBinning(&fSigmaTheta, &binst, &binsb); - SetBinning(&fSigmaThetaOuter, &binst, &binsb); - SetBinning(&fSigmaPixTheta, &binst, &binspix, &binsb); - SetBinning(&fDiffPixTheta, &binst, &binspix, &binsd); - + fThetaPhi = new TH2D; + fThetaPhi->SetDirectory(NULL); + fThetaPhi->UseCurrentStyle(); + fThetaPhi->SetName("2D-ThetaPhi"); + fThetaPhi->SetTitle("2D: \\Theta vs. \\Phi"); + fThetaPhi->SetXTitle("\\phi [\\circ]"); + fThetaPhi->SetYTitle("\\Theta [\\circ]"); + fThetaPhi->SetTitleOffset(1.2,"Y"); + + fSigmaTheta = new TH2D; + fSigmaTheta->SetDirectory(NULL); + fSigmaTheta->UseCurrentStyle(); + fSigmaTheta->SetName("2D-ThetaSigmabar(Inner)"); + fSigmaTheta->SetTitle("2D: \\bar{\\sigma}, \\Theta"); + fSigmaTheta->SetXTitle("\\Theta [\\circ]"); + fSigmaTheta->SetYTitle("Sigmabar(Inner)"); + fSigmaTheta->SetTitleOffset(1.2,"Y"); + + fSigmaThetaOuter = new TH2D; + fSigmaThetaOuter->SetDirectory(NULL); + fSigmaThetaOuter->UseCurrentStyle(); + fSigmaThetaOuter->SetName("2D-ThetaSigmabar(Outer)"); + fSigmaThetaOuter->SetTitle("2D: \\bar{\\sigma}, \\Theta"); + fSigmaThetaOuter->SetXTitle("\\Theta [\\circ]"); + fSigmaThetaOuter->SetYTitle("Sigmabar(Outer)"); + fSigmaThetaOuter->SetTitleOffset(1.2,"Y"); + + fSigmaPixTheta = new TH3D; + fSigmaPixTheta->SetDirectory(NULL); + fSigmaPixTheta->UseCurrentStyle(); + fSigmaPixTheta->SetName("3D-ThetaPixSigma"); + fSigmaPixTheta->SetTitle("3D: \\Theta, Pixel Id, \\sigma"); + fSigmaPixTheta->SetXTitle("\\Theta [\\circ]"); + fSigmaPixTheta->SetYTitle("Pixel Id"); + fSigmaPixTheta->SetZTitle("Sigma"); + + fDiffPixTheta = new TH3D; + fDiffPixTheta->SetDirectory(NULL); + fDiffPixTheta->UseCurrentStyle(); + fDiffPixTheta->SetName("3D-ThetaPixDiff"); + //fDiffPixTheta->SetTitle("3D: \\Theta, Pixel Id, {\\sigma}^{2}-\\bar{\\sigma}^{2}"); + fDiffPixTheta->SetTitle("3D: \\Theta, Pixel Id, (Sigma2-Sigmabar2)/Area"); + fDiffPixTheta->SetXTitle("\\Theta [\\circ]"); + fDiffPixTheta->SetYTitle("Pixel Id"); + fDiffPixTheta->SetZTitle("(Sigma2 - Sigmabar2)/Area"); + + // Define default binning + fBinsPhi = new MBinning; + fBinsPhi->SetEdges( 20, 0.0, 360.0); + + Double_t fThetaLo = 0.0; + Double_t fThetaHi = 90.0; + fBinsTheta = new MBinning; + fBinsTheta->SetEdgesCos( 10, fThetaLo, fThetaHi); // theta + //fBinsTheta->SetEdges( 1, fThetaLo, fThetaHi); // theta + + fBinsSigma = new MBinning; + fBinsSigma->SetEdges( 100, 0.0, 120.0); // sigma + + fBinsSigmabarIn = new MBinning; + fBinsSigmabarOut = new MBinning;; + fBinsSigmabarIn->SetEdges(100, 0.0, 40.0); // sigmabar (inner) + fBinsSigmabarOut->SetEdges(100, 0.0, 80.0); // sigmabar (outer) + + fBinsPix = new MBinning; + fBinsPix->SetEdges(578, -0.5, 577.5); + + fBinsDiff = new MBinning; + fBinsDiff->SetEdges( 100, -500.0, 1500.0); // (sigma2-sigmabar2)/area + + SetBinning(fThetaPhi, fBinsPhi, fBinsTheta); + SetBinning(fSigmaTheta, fBinsTheta, fBinsSigmabarIn); + SetBinning(fSigmaThetaOuter, fBinsTheta, fBinsSigmabarOut); + SetBinning(fSigmaPixTheta, fBinsTheta, fBinsPix, fBinsSigma); + SetBinning(fDiffPixTheta, fBinsTheta, fBinsPix, fBinsDiff); + + //-------------------------------------------- fNamePedPhotCam = "MPedPhotCamFromData"; +} + + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MHSigmaTheta::~MHSigmaTheta() +{ + delete fThetaPhi; + delete fSigmaTheta; + delete fSigmaThetaOuter; + delete fSigmaPixTheta; + delete fDiffPixTheta; + + delete fBinsPhi; + delete fBinsTheta; + delete fBinsSigma; + delete fBinsSigmabarIn; + delete fBinsSigmabarOut; + delete fBinsPix; + delete fBinsDiff; } @@ -133,7 +189,8 @@ fPointPos = (MPointingPos*)plist->FindObject("MPointingPos"); if (!fPointPos) - *fLog << warn << "MPointingPos not found... aborting." << endl; - - + { + *fLog << err << "MPointingPos not found... aborting." << endl; + return kFALSE; + } fPed = (MPedPhotCam*)plist->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam"); @@ -144,5 +201,5 @@ return kFALSE; } - //fPed->InitSize(fCam->GetNumPixels()); + fPed->InitSize(fCam->GetNumPixels()); @@ -162,40 +219,76 @@ + //--------------------------------------------------------------- + *fLog << inf << "Name of MPedPhotCam container : " + << fNamePedPhotCam << endl; + + + //--------------------------------------------------------------- + + // Get Phi Binning + MBinning* binsphi = (MBinning*)plist->FindObject("BinningPhi", "MBinning"); + if (!binsphi) + { + *fLog << warn << "Object 'BinningPhi' [MBinning] not found... use default binning." << endl; + binsphi = fBinsPhi; + } + // Get Theta Binning MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning"); if (!binstheta) { - *fLog << warn << "Object 'BinningTheta' [MBinning] not found... no binning applied." << endl; - return kTRUE; - } - - // Get Sigmabar binning - MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar", "MBinning"); - if (!binssigma) - *fLog << warn << "Object 'BinningSigmabar' [MBinning] not found... no binning applied." << endl; + *fLog << warn << "Object 'BinningTheta' [MBinning] not found... use default binning." << endl; + binstheta = fBinsTheta; + } + + // Get Sigma binning + MBinning* binssig = (MBinning*)plist->FindObject("BinningSigma", "MBinning"); + if (!binssig) + { + *fLog << warn << "Object 'BinningSigma' [MBinning] not found... use default binning." << endl; + binssig = fBinsSigma; + } + + // Get SigmabarIn binning + MBinning* binssigmain = (MBinning*)plist->FindObject("BinningSigmabarIn", "MBinning"); + if (!binssigmain) + { + *fLog << warn << "Object 'BinningSigmabarIn' [MBinning] not found... use default binning." << endl; + binssigmain = fBinsSigmabarIn; + } + + // Get SigmabarOut binning + MBinning* binssigmaout = (MBinning*)plist->FindObject("BinningSigmabarOut", "MBinning"); + if (!binssigmaout) + { + *fLog << warn << "Object 'BinningSigmabarOut' [MBinning] not found... use default binning." << endl; + binssigmaout = fBinsSigmabarOut; + } // Get binning for (sigma^2-sigmabar^2) MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2", "MBinning"); if (!binsdiff) - *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... no binning applied." << endl; - - //FIXME: Missing: Apply new binning to only one axis! + { + *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... use default binning." << endl; + binsdiff = fBinsDiff; + } + + + //--------------------------------------------------------------- // Get binning for pixel number const UInt_t npix1 = fPed->GetSize()+1; - *fLog << "npix1 = " << npix1 << endl; + //*fLog << "MHSigmaTheata::SetupFill(); npix1 = " << npix1 << endl; MBinning binspix("BinningPixel"); binspix.SetEdges(npix1, -0.5, npix1-0.5); // Set binnings in histograms - if (binssigma) - { - SetBinning(&fSigmaTheta, binstheta, binssigma); - SetBinning(&fSigmaThetaOuter, binstheta, binssigma); - SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma); - } - - if (binsdiff) - SetBinning(&fDiffPixTheta, binstheta, &binspix, binsdiff); + SetBinning(fThetaPhi, binsphi, binstheta); + SetBinning(fSigmaTheta, binstheta, binssigmain); + SetBinning(fSigmaThetaOuter, binstheta, binssigmaout); + SetBinning(fSigmaPixTheta, binstheta, &binspix, binssig); + SetBinning(fDiffPixTheta, binstheta, &binspix, binsdiff); + + *fLog << "MHSigmaTheta::SetupFill(); binnings were set" << endl; return kTRUE; @@ -211,22 +304,25 @@ { Double_t theta = fPointPos->GetZd(); + Double_t phi = fPointPos->GetAz(); + + fPed->ReCalc(*fCam, fBad); Double_t mysig = (fPed->GetArea(0)).GetRms(); Double_t mysigouter = (fPed->GetArea(1)).GetRms(); - *fLog << "theta, mysig, mysigouter = " << theta << ", " << mysig - << ", " << mysigouter << endl; - - fSigmaTheta.Fill(theta, mysig); - fSigmaThetaOuter.Fill(theta, mysigouter); - - const UInt_t npix = fEvt->GetNumPixels(); - - for (UInt_t i=0; iFill(phi, theta); + fSigmaTheta->Fill(theta, mysig); + fSigmaThetaOuter->Fill(theta, mysigouter); + + MCerPhotPix *cerpix = 0; + TIter Next(*fEvt); + while ( (cerpix=(MCerPhotPix*)Next()) ) + { + const Int_t id = cerpix->GetPixId(); + + if (!cerpix->IsPixelUsed()) { //*fLog << all << "MHSigmaTheta::Fill; unused pixel found, id = " @@ -237,13 +333,4 @@ const MPedPhotPix &pix = (*fPed)[id]; - if ( fBad != NULL && ((*fBad)[id]).IsUnsuitable() ) - { - // this should never occur, because bad pixels should have - // been set unused - *fLog << all << "MHSigmaTheta::Fill; bad pixel found which is used, id = " - << id << "... go to next pixel." << endl; - continue; - } - // ratio is the area of pixel 0 // divided by the area of the current pixel @@ -251,5 +338,5 @@ const Double_t sigma = pix.GetRms(); - fSigmaPixTheta.Fill(theta, (Double_t)id, sigma*sqrt(ratio)); + fSigmaPixTheta->Fill(theta, (Double_t)id, sigma); Double_t diff; @@ -258,13 +345,13 @@ { // inner pixel - diff = sigma*sigma*ratio - mysig*mysig; + diff = (sigma*sigma - mysig*mysig) * ratio; } else { // outer pixel - diff = sigma*sigma*ratio - mysigouter*mysigouter; + diff = (sigma*sigma - mysigouter*mysigouter) * ratio; } - fDiffPixTheta.Fill(theta, (Double_t)id, diff); + fDiffPixTheta->Fill(theta, (Double_t)id, diff); } @@ -285,5 +372,5 @@ if ((h = (TH1D*)gPad->FindObject("ProjX-Theta"))) { - ProjectionX(*h, fSigmaTheta); + ProjectionX(*h, *fSigmaTheta); if (h->GetEntries()!=0) gPad->SetLogy(); @@ -292,5 +379,5 @@ padsave->cd(4); if ((h = (TH1D*)gPad->FindObject("ProjY-sigma"))) - ProjectionY(*h, fSigmaTheta); + ProjectionY(*h, *fSigmaTheta); gPad = padsave; @@ -305,5 +392,4 @@ TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); pad->SetBorderMode(0); - AppendPad(""); @@ -315,9 +401,11 @@ pad->cd(1); gPad->SetBorderMode(0); - h = fSigmaTheta.ProjectionX("ProjX-Theta", -1, 9999, "E"); + h = fSigmaTheta->ProjectionX("ProjX-Theta", -1, 9999, "E"); h->SetDirectory(NULL); + h->UseCurrentStyle(); h->SetTitle("Distribution of \\Theta"); h->SetXTitle("\\Theta [\\circ]"); h->SetYTitle("No.of events"); + h->SetTitleOffset(1.2,"Y"); h->Draw(opt); h->SetBit(kCanDelete); @@ -325,9 +413,12 @@ pad->cd(2); gPad->SetBorderMode(0); - h = fDiffPixTheta.Project3D("zx"); + h = fDiffPixTheta->Project3D("zx"); h->SetDirectory(NULL); - h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)"); + h->UseCurrentStyle(); + //h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)"); + h->SetTitle("(Sigma2-Sigmabar2)/Area vs. \\Theta (all pixels)"); h->SetXTitle("\\Theta [\\circ]"); h->SetYTitle("(Sigma2 - Sigmabar2) / Area"); + h->SetTitleOffset(1.2,"Y"); h->Draw("box"); h->SetBit(kCanDelete); @@ -335,9 +426,11 @@ pad->cd(3); gPad->SetBorderMode(0); - h = fSigmaPixTheta.Project3D("zx"); + h = fSigmaPixTheta->Project3D("zx"); h->SetDirectory(NULL); + h->UseCurrentStyle(); h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)"); h->SetXTitle("\\Theta [\\circ]"); h->SetYTitle("Sigma"); + h->SetTitleOffset(1.2,"Y"); h->Draw("box"); h->SetBit(kCanDelete); @@ -345,9 +438,11 @@ //pad->cd(7); //gPad->SetBorderMode(0); - //h = fSigmaTheta.ProjectionY("ProjY-sigma", -1, 9999, "E"); + //h = fSigmaTheta->ProjectionY("ProjY-sigma", -1, 9999, "E"); //h->SetDirectory(NULL); + //h->UseCurrentStyle(); //h->SetTitle("Distribution of \\bar{\\sigma}_{ped}"); //h->SetXTitle("\\bar{\\sigma}_{ped}"); //h->SetYTitle("No.of events"); + //h->SetTitleOffset(1.2,"Y"); //h->Draw(opt); //h->SetBit(kCanDelete); @@ -355,9 +450,12 @@ pad->cd(5); gPad->SetBorderMode(0); - h = fDiffPixTheta.Project3D("zy"); + h = fDiffPixTheta->Project3D("zy"); h->SetDirectory(NULL); - h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)"); + h->UseCurrentStyle(); + //h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)"); + h->SetTitle("(Sigma2-Sigmabar2)/Area vs. pixel Id (all \\Theta)"); h->SetXTitle("Pixel Id"); - h->SetYTitle("Sigma2 - sigmabar2) / Area"); + h->SetYTitle("(Sigma2 - sigmabar2) / Area"); + h->SetTitleOffset(1.2,"Y"); h->Draw("box"); h->SetBit(kCanDelete); @@ -365,23 +463,37 @@ pad->cd(6); gPad->SetBorderMode(0); - h = fSigmaPixTheta.Project3D("zy"); + h = fSigmaPixTheta->Project3D("zy"); h->SetDirectory(NULL); + h->UseCurrentStyle(); h->SetTitle("\\sigma_{ped} vs. pixel Id (all \\Theta)"); h->SetXTitle("Pixel Id"); h->SetYTitle("Sigma"); + h->SetTitleOffset(1.2,"Y"); h->Draw("box"); h->SetBit(kCanDelete); pad->cd(4); - fSigmaTheta.Draw(opt); + gPad->SetBorderMode(0); + fSigmaTheta->Draw(); pad->cd(7); - fSigmaThetaOuter.Draw(opt); - - //pad->cd(8); - //fDiffPixTheta.Draw(opt); - - //pad->cd(9); - //fSigmaPixTheta.Draw(opt); + gPad->SetBorderMode(0); + fSigmaThetaOuter->Draw(); + + pad->cd(8); + gPad->SetBorderMode(0); + fThetaPhi->Draw(); + + pad->cd(9); + gPad->SetBorderMode(0); + h = fThetaPhi->ProjectionX("ProjX-Phi", -1, 9999, "E"); + h->SetDirectory(NULL); + h->UseCurrentStyle(); + h->SetTitle("Distribution of \\Phi"); + h->SetXTitle("\\phi [\\circ]"); + h->SetYTitle("No. of events"); + h->SetTitleOffset(1.2,"Y"); + h->Draw(opt); + h->SetBit(kCanDelete); } Index: /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 5430) +++ /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 5431) @@ -35,9 +35,19 @@ // sigmabar is the average pedestal sigma - TH2D fSigmaTheta; // 2D-distribution sigmabar versus Theta (Inner) - TH2D fSigmaThetaOuter; // 2D-distribution sigmabar versus Theta (Outer) + TH2D *fSigmaTheta; // 2D-distribution sigmabar versus Theta (Inner) + TH2D *fSigmaThetaOuter; // 2D-distribution sigmabar versus Theta (Outer) - TH3D fSigmaPixTheta; // 3D-distr.:Theta, pixel, pedestal sigma - TH3D fDiffPixTheta; // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2 + TH3D *fSigmaPixTheta; // 3D-distr.:Theta, pixel, pedestal sigma + TH3D *fDiffPixTheta; // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2 + + TH2D *fThetaPhi; // 2D-distribution Theta versus Phi + + MBinning *fBinsPhi; //! + MBinning *fBinsTheta; //! + MBinning *fBinsSigma; //! + MBinning *fBinsSigmabarIn; //! + MBinning *fBinsSigmabarOut; //! + MBinning *fBinsPix; //! + MBinning *fBinsDiff; //! void Paint(Option_t *opt=""); @@ -45,4 +55,5 @@ public: MHSigmaTheta(const char *name=NULL, const char *title=NULL); + ~MHSigmaTheta(); void SetNamePedPhotCam(const char *name) { fNamePedPhotCam = name; } @@ -51,15 +62,18 @@ Bool_t Fill(const MParContainer *par, const Stat_t w=1); - const TH2D *GetSigmaTheta() { return &fSigmaTheta; } - const TH2D *GetSigmaTheta() const { return &fSigmaTheta; } + const TH2D *GetThetaPhi() { return fThetaPhi; } + const TH2D *GetThetaPhi() const { return fThetaPhi; } - const TH2D *GetSigmaThetaOuter() { return &fSigmaThetaOuter; } - const TH2D *GetSigmaThetaOuter() const { return &fSigmaThetaOuter; } + const TH2D *GetSigmaTheta() { return fSigmaTheta; } + const TH2D *GetSigmaTheta() const { return fSigmaTheta; } - const TH3D *GetSigmaPixTheta() { return &fSigmaPixTheta; } - const TH3D *GetSigmaPixTheta() const { return &fSigmaPixTheta; } + const TH2D *GetSigmaThetaOuter() { return fSigmaThetaOuter; } + const TH2D *GetSigmaThetaOuter() const { return fSigmaThetaOuter; } - const TH3D *GetDiffPixTheta() { return &fDiffPixTheta; } - const TH3D *GetDiffPixTheta() const { return &fDiffPixTheta; } + const TH3D *GetSigmaPixTheta() { return fSigmaPixTheta; } + const TH3D *GetSigmaPixTheta() const { return fSigmaPixTheta; } + + const TH3D *GetDiffPixTheta() { return fDiffPixTheta; } + const TH3D *GetDiffPixTheta() const { return fDiffPixTheta; } void Draw(Option_t *option=""); Index: /trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc =================================================================== --- /trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc (revision 5430) +++ /trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc (revision 5431) @@ -560,4 +560,6 @@ return kFALSE; } + *fLog << all << "Name of MPedPhotCam container = " << fNamePedPhotCam + << endl; fTime = (MArrivalTime*)pList->FindObject(AddSerialNumber("MArrivalTime"));