Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 2127) +++ /trunk/MagicSoft/Mars/Changelog (revision 2128) @@ -1,3 +1,21 @@ -*-*- END OF LINE -*-*- + + 2003/05/21: Wolfgang Wittek + + * mhist/MHBlindPixels.[h,cc] + - change 1D histogram into 2D histogram (pixel Id vs. Theta) + - add 2D histogram : no.of blind pixels vs. Theta + + * mhist/MHSigmaTheta.cc + - correct "BinningPix" + + * manalysis/MPadSchweizer.[h,cc] + - add simulation of blind pixels + + * mhist/MHMatrix.cc + - in DefRefMatrix : allow variable bin size for 'hth' and 'hthd' + + + 2003/05/20: Oscar Blanch Bigas @@ -50,5 +68,4 @@ * manalysis/MFiltercutsCalc.[h,cc]: - added - Index: /trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc (revision 2127) +++ /trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc (revision 2128) @@ -86,4 +86,5 @@ #include "MPedestalCam.h" #include "MPedestalPix.h" +#include "MBlindPixels.h" ClassImp(MPadSchweizer); @@ -101,4 +102,6 @@ fHSigmaPixTheta = NULL; fHDiffPixTheta = NULL; + fHBlindPixIdTheta = NULL; + fHBlindPixNTheta = NULL; fHSigmaPedestal = NULL; @@ -125,15 +128,22 @@ // fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma) // fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) -// -// -void MPadSchweizer::SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff) +// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id) +// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels ) +// +// +void MPadSchweizer::SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff, + TH2D *hist2Pix, TH2D *hist2PixN) { fHSigmaTheta = hist2; fHSigmaPixTheta = hist3; fHDiffPixTheta = hist3Diff; + fHBlindPixIdTheta = hist2Pix; + fHBlindPixNTheta = hist2PixN; fHSigmaTheta->SetDirectory(NULL); fHSigmaPixTheta->SetDirectory(NULL); fHDiffPixTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetDirectory(NULL); Print(); @@ -171,5 +181,6 @@ Bool_t MPadSchweizer::PreProcess(MParList *pList) { - if ( !fHSigmaTheta || !fHSigmaPixTheta || !fHDiffPixTheta) + if ( !fHSigmaTheta || !fHSigmaPixTheta || !fHDiffPixTheta || + !fHBlindPixIdTheta || !fHBlindPixNTheta) { *fLog << err << "At least one of the histograms needed for the padding is not defined ... aborting." << endl; @@ -190,5 +201,5 @@ return kFALSE; } - + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); if (!fCam) @@ -211,4 +222,11 @@ return kFALSE; } + + fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); + if (!fBlinds) + { + *fLog << err << "MBlindPixels not found... aborting." << endl; + return kFALSE; + } @@ -296,4 +314,81 @@ //------------------------------------------- // for the current theta, + // generate blind pixels according to the histograms + // fHBlindPixNTheta and fHBlindPixIDTheta + // + + + Int_t binPix = fHBlindPixNTheta->GetXaxis()->FindBin(theta); + + if ( binPix < 1 || binPix > fHBlindPixNTheta->GetNbinsX() ) + { + //*fLog << "MPadSchweizer::Process(); binNumber out of range : theta, binPix = " + // << theta << ", " << binPix << "; Skip event " << endl; + // event cannot be padded; skip event + + rc = 2; + fErrors[rc]++; + return kCONTINUE; + } + + // numBlind is the number of blind pixels in this event + TH1D *nblind; + UInt_t numBlind; + + nblind = fHBlindPixNTheta->ProjectionY("", binPix, binPix, ""); + if ( nblind->GetEntries() == 0.0 ) + { + *fLog << "MPadSchweizer::Process(); projection for Theta bin " + << binPix << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete nblind; + + rc = 7; + fErrors[rc]++; + return kCONTINUE; + } + else + { + numBlind = (Int_t) (nblind->GetRandom()+0.5); + + //*fLog << "numBlind = " << numBlind << endl; + } + delete nblind; + + + // throw the Id of numBlind different pixels in this event + TH1D *hblind; + UInt_t idBlind; + UInt_t listId[npix]; + UInt_t nlist = 0; + Bool_t equal; + + hblind = fHBlindPixIdTheta->ProjectionY("", binPix, binPix, ""); + if ( hblind->GetEntries() > 0.0 ) + for (UInt_t i=0; iGetRandom()+0.5); + equal = kFALSE; + for (UInt_t j=0; jSetPixelBlind(idBlind); + //*fLog << "idBlind = " << idBlind << endl; + } + delete hblind; + + + //------------------------------------------- + // for the current theta, // generate a sigmabar according to the histogram fHSigmaTheta // @@ -301,13 +396,22 @@ Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(theta); + if (binPix != binNumber) + { + cout << "The binnings of the 2 histograms are not identical; aborting" + << endl; + return kERROR; + } + TH1D *hsigma; - if ( binNumber < 1 || binNumber > fHSigmaTheta->GetNbinsX() ) + hsigma = fHSigmaTheta->ProjectionY("", binNumber, binNumber, ""); + if ( hsigma->GetEntries() == 0.0 ) { - //*fLog << "MPadSchweizer::Process(); binNumber out of range : theta, binNumber = " - // << theta << ", " << binNumber << "; Skip event " << endl; + *fLog << "MPadSchweizer::Process(); projection for Theta bin " + << binNumber << " has no entries; Skip event " << endl; // event cannot be padded; skip event - - rc = 2; + delete hsigma; + + rc = 3; fErrors[rc]++; return kCONTINUE; @@ -315,25 +419,9 @@ else { - hsigma = fHSigmaTheta->ProjectionY("", binNumber, binNumber, ""); - if ( hsigma->GetEntries() == 0.0 ) - { - *fLog << "MPadSchweizer::Process(); projection for Theta bin " - << binNumber << " has no entries; Skip event " << endl; - // event cannot be padded; skip event - delete hsigma; - - rc = 3; - fErrors[rc]++; - return kCONTINUE; - } - else - { - sigmabar = hsigma->GetRandom(); - - //*fLog << "Theta, bin number = " << theta << ", " << binNumber - // << ", sigmabar = " << sigmabar << endl; - } - delete hsigma; - } + sigmabar = hsigma->GetRandom(); + //*fLog << "Theta, bin number = " << theta << ", " << binNumber // << ", sigmabar = " << sigmabar << endl + } + delete hsigma; + const Double_t sigmabar2 = sigmabar*sigmabar; @@ -644,4 +732,8 @@ << "%) Evts skipped due to: No data for generating Sigma" << endl; + *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) + << (int)(fErrors[7]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Blind pixels" << endl; + *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) Index: /trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h (revision 2127) +++ /trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h (revision 2128) @@ -20,4 +20,5 @@ class MSigmabar; class MParList; +class MBlindPixels; class MPadSchweizer : public MTask @@ -29,4 +30,5 @@ MMcEvt *fMcEvt; MPedestalCam *fPed; + MBlindPixels *fBlinds; Int_t fPadFlag; @@ -34,7 +36,9 @@ Int_t fGroup; - Int_t fErrors[7]; + Int_t fErrors[8]; // plots used for the padding + TH2D *fHBlindPixIdTheta; // 2D-histogram (blind pixel Id vs. Theta) + TH2D *fHBlindPixNTheta; // 2D-histogram (no.of blind pixels vs. Theta) TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) @@ -53,5 +57,6 @@ ~MPadSchweizer(); - void SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff); + void SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff, + TH2D *hist2Pix, TH2D *hist2PixN); Bool_t PreProcess(MParList *pList); Index: /trunk/MagicSoft/Mars/mhist/MHBlindPixels.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHBlindPixels.cc (revision 2127) +++ /trunk/MagicSoft/Mars/mhist/MHBlindPixels.cc (revision 2128) @@ -1,2 +1,3 @@ + /* ======================================================================== *\ ! @@ -32,5 +33,12 @@ #include +#include "MMcEvt.hxx" #include "MBlindPixels.h" +#include "MPedestalCam.h" +#include "MParList.h" +#include "MBinning.h" + +#include "MLog.h" +#include "MLogManip.h" ClassImp(MHBlindPixels); @@ -41,20 +49,65 @@ // MHBlindPixels::MHBlindPixels(const char *name, const char *title) - : fHist(fName, fTitle, 577, -.5, 577-.5) { fName = name ? name : "MHBlindPixels"; - fTitle = title ? title : "Histogram for Blind Pixels"; + fTitle = title ? title : "Histogram for Blind Pixels vs. Theta"; // - we initialize the histogram - // - we have to give diferent names for the diferent histograms because + // - we have to give different names for the different histograms because // root don't allow us to have diferent histograms with the same name - fHist.SetName("1D-BlindPixels"); - fHist.SetTitle(fTitle); + fBlindId.SetName("2D-IdBlindPixels"); + fBlindId.SetTitle("2D-IdBlindPixels"); + fBlindId.SetDirectory(NULL); + fBlindId.SetXTitle("\\Theta [\\circ]"); + fBlindId.SetYTitle("pixel Id"); - fHist.SetDirectory(NULL); + fBlindN.SetName("2D-NBlindPixels"); + fBlindN.SetTitle("2D-NBlindPixels"); + fBlindN.SetDirectory(NULL); + fBlindN.SetXTitle("\\Theta [\\circ]"); + fBlindN.SetYTitle("number of blind pixels"); +} - fHist.SetXTitle("Id"); - fHist.SetYTitle("Counts"); +// -------------------------------------------------------------------------- +// +// Set the binnings and prepare the filling of the histogram +// +Bool_t MHBlindPixels::SetupFill(const MParList *plist) +{ + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)plist->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << err << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + + // Get Theta Binning + MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning"); + if (!binstheta) + { + *fLog << err << "Object 'BinningTheta' [MBinning] not found... aborting" << endl; + return kFALSE; + } + + // Get binning for pixel number + const UInt_t npix1 = fPed->GetSize()+1; + + MBinning binspix("BinningPixel"); + binspix.SetEdges(npix1, -0.5, npix1-0.5); + + // Set binnings in histograms + SetBinning(&fBlindId, binstheta, &binspix); + SetBinning(&fBlindN, binstheta, &binspix); + + return kTRUE; } @@ -68,5 +121,33 @@ pad->SetBorderMode(0); - fHist.Draw(option); + pad->Divide(2,2); + + TH1D *h; + + pad->cd(1); + fBlindId.Draw(option); + + pad->cd(2); + fBlindN.Draw(option); + + pad->cd(3); + gPad->SetBorderMode(0); + h = ((TH2*)&fBlindId)->ProjectionY("ProjY-pixId", -1, 9999, ""); + h->SetDirectory(NULL); + h->SetTitle("Distribution of blind pixel Id"); + h->SetXTitle("Id of blind pixel"); + h->SetYTitle("No. of events"); + h->Draw(option); + h->SetBit(kCanDelete); + + pad->cd(4); + gPad->SetBorderMode(0); + h = ((TH2*)&fBlindN)->ProjectionY("ProjY-pixN", -1, 9999, ""); + h->SetDirectory(NULL); + h->SetTitle("Distribution of no.of blind pixels"); + h->SetXTitle("No. of blind pixels"); + h->SetYTitle("No. of events"); + h->Draw(option); + h->SetBit(kCanDelete); pad->Modified(); @@ -79,11 +160,33 @@ return kFALSE; + Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg; const MBlindPixels &bp = *(MBlindPixels*)par; // FIXME: Slow. - for (int i=0; i<577; i++) + const UInt_t npix = fPed->GetSize(); + + UInt_t nb = 0; + for (UInt_t i=0; i +#ifndef ROOT_TH2 +#include #endif + +class MPedestalCam; +class MMcEvt; +class MParList; + class MHBlindPixels : public MH { private: - TH1D fHist; // + MPedestalCam *fPed; //! + MMcEvt *fMcEvt; //! + + TH2D fBlindId; // 2D-histogram : pixel Id vs. Theta + TH2D fBlindN; // 2D-histogram : no.of blind pixels vs. Theta public: MHBlindPixels(const char *name=NULL, const char *title=NULL); - const TH1D *GetHist() { return &fHist; } - const TH1D *GetHist() const { return &fHist; } + const TH2D *GetBlindId() { return &fBlindId; } + const TH2D *GetBlindId() const { return &fBlindId; } - TH1 *GetHistByName(const TString name) { return &fHist; } + const TH2D *GetBlindN() { return &fBlindN; } + const TH2D *GetBlindN() const { return &fBlindN; } + + TH2 *GetBlinIdByName(const TString name) { return &fBlindId; } + TH2 *GetBlinNByName(const TString name) { return &fBlindN; } void Draw(Option_t* option = ""); + Bool_t SetupFill(const MParList *plist); Bool_t Fill(const MParContainer *par, const Stat_t w=1); - ClassDef(MHBlindPixels, 1) // Histogram of blind pixels + ClassDef(MHBlindPixels, 1) // Histogram of blind pixel Id vs. Theta }; #endif + + Index: /trunk/MagicSoft/Mars/mhist/MHMatrix.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 2127) +++ /trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 2128) @@ -788,8 +788,9 @@ // Calculate normalization factors // - const int nbins = thsh.GetNbinsX(); - const double frombin = thsh.GetBinLowEdge(1); - const double tobin = thsh.GetBinLowEdge(nbins+1); - const double dbin = thsh.GetBinWidth(1); + //const int nbins = thsh.GetNbinsX(); + //const double frombin = thsh.GetBinLowEdge(1); + //const double tobin = thsh.GetBinLowEdge(nbins+1); + //const double dbin = thsh.GetBinWidth(1); + const Int_t nrows = fM.GetNrows(); const Int_t ncols = fM.GetNcols(); @@ -798,10 +799,16 @@ // set up the real histogram (distribution before) // - TH1F hth("th", "Distribution before reduction", nbins, frombin, tobin); + //TH1F hth("th", "Distribution before reduction", nbins, frombin, tobin); + TH1F hth; + hth.SetNameTitle("th", "Distribution before reduction"); + SetBinning(&hth, &thsh); hth.SetDirectory(NULL); for (Int_t j=0; jGetSize(); + const UInt_t npix1 = fPed->GetSize()+1; MBinning binspix("BinningPixel"); - binspix.SetEdges(npix, -0.5, 0.5+npix); + binspix.SetEdges(npix1, -0.5, npix1-0.5); // Set binnings in histograms