Changeset 2128 for trunk/MagicSoft/Mars

Timestamp:

05/21/03 12:21:58 (22 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• manalysis/MPadSchweizer.cc (modified) (10 diffs)
• manalysis/MPadSchweizer.h (modified) (4 diffs)
• mhist/MHBlindPixels.cc (modified) (5 diffs)
• mhist/MHBlindPixels.h (modified) (1 diff)
• mhist/MHMatrix.cc (modified) (4 diffs)
• mhist/MHSigmaTheta.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- 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 +68 @@
    * manalysis/MFiltercutsCalc.[h,cc]:
      - added
-
 
 

trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc

@@ -86 +86 @@
 #include "MPedestalCam.h"
 #include "MPedestalPix.h"
+#include "MBlindPixels.h"
 
 ClassImp(MPadSchweizer);
@@ -101 +102 @@
   fHSigmaPixTheta = NULL;
   fHDiffPixTheta  = NULL;
+  fHBlindPixIdTheta = NULL;
+  fHBlindPixNTheta  = NULL;
 
   fHSigmaPedestal = NULL;
@@ -125 +128 @@
 // 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 +181 @@
 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 +201 @@
        return kFALSE;
      }
-  
+
    fCam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!fCam)
@@ -211 +222 @@
        return kFALSE;
      }
+
+   fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
+   if (!fBlinds)
+     {
+       *fLog << err << "MBlindPixels not found... aborting." << endl;
+       return kFALSE;
+     }
    
 
@@ -296 +314 @@
   //-------------------------------------------
   // 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; i<numBlind; i++)
+    {
+      while(1)
+      {
+        idBlind = (Int_t) (hblind->GetRandom()+0.5);
+        equal = kFALSE;
+        for (UInt_t j=0; j<nlist; j++)
+          if (idBlind == listId[j])
+          { 
+            equal = kTRUE;
+            break;
+          }
+        if (!equal) break;
+      }
+      listId[nlist] = idBlind;
+      nlist++;
+
+      fBlinds->SetPixelBlind(idBlind);
+      //*fLog << "idBlind = " << idBlind << endl;
+    }
+  delete hblind;
+
+
+  //-------------------------------------------
+  // for the current theta, 
   // generate a sigmabar according to the histogram fHSigmaTheta
   //
@@ -301 +396 @@
   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 +419 @@
   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 +732 @@
           << "%) 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()) 

trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h

@@ -20 +20 @@
 class MSigmabar;
 class MParList;
+class MBlindPixels;
 
 class MPadSchweizer : public MTask
@@ -29 +30 @@
     MMcEvt         *fMcEvt;
     MPedestalCam   *fPed;
+    MBlindPixels   *fBlinds;
 
     Int_t          fPadFlag;
@@ -34 +36 @@
     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 +57 @@
     ~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);

trunk/MagicSoft/Mars/mhist/MHBlindPixels.cc

@@ +1 @@
+
 /* ======================================================================== *\
 !
@@ -32 +33 @@
 #include <TCanvas.h>
 
+#include "MMcEvt.hxx"
 #include "MBlindPixels.h"
+#include "MPedestalCam.h"
+#include "MParList.h"
+#include "MBinning.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
 
 ClassImp(MHBlindPixels);
@@ -41 +49 @@
 //
 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 +121 @@
     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 +160 @@
         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<npix; i++)
+    {
         if (bp.IsBlind(i))
-            fHist.Fill(i, w);
+        {
+          fBlindId.Fill(theta, i, w);
+          nb++;
+        }   
+    }
+    fBlindN.Fill(theta, nb, w);
 
     return kTRUE;
 }
+
+
+
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/mhist/MHBlindPixels.h

@@ -5 +5 @@
 #include "MH.h"
 #endif
-#ifndef ROOT_TH1
-#include <TH1.h>
+#ifndef ROOT_TH2
+#include <TH2.h>
 #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
+
+

trunk/MagicSoft/Mars/mhist/MHMatrix.cc

@@ -788 +788 @@
     // 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 +799 @@
     // 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; j<nrows; j++)
         hth.Fill(fM(j, refcolumn));
 
-    TH1F hthd("thd", "Correction factors", nbins, frombin, tobin);
+    //TH1F hthd("thd", "Correction factors", nbins, frombin, tobin);
+    TH1F hthd;
+    hthd.SetNameTitle("thd", "Correction factors");
+    SetBinning(&hthd, &thsh);
     hthd.SetDirectory(NULL);
     hthd.Divide((TH1F*)&thsh, &hth, 1, 1);
@@ -840 +847 @@
     TVector v(fM.GetNrows());
     v = TMatrixColumn(fM, refcolumn);
-    v += -frombin;
-    v *= 1/dbin;
+    //v += -frombin;
+    //v *= 1/dbin;
 
     //
@@ -850 +857 @@
     for (ir=0; ir<nrows; ir++)
     {
-        const Int_t indref = (Int_t)v(ind[ir]);
-
+        // const Int_t indref = (Int_t)v(ind[ir]);
+        const Int_t indref = hthd.FindBin(v(ind[ir])) - 1;
         cumulweight[indref] += hthd.GetBinContent(indref+1);
         if (cumulweight[indref]<=0.5)

trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc

@@ -100 +100 @@
 
     MBinning binspix("BinningPixel");
-    binspix.SetEdges(577, -0.5, 577.5);
+    binspix.SetEdges(578, -0.5, 577.5);
 
     SetBinning(&fSigmaTheta,    &binst, &binsb);
@@ -169 +169 @@
 
     // Get binning for pixel number
-    const UInt_t npix = fPed->GetSize();
+    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