Changeset 4584 for trunk/MagicSoft/Mars/manalysis

Timestamp:

08/12/04 06:58:34 (20 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/manalysis

Files:

• MPad.cc (modified) (87 diffs)
• MPad.h (modified) (5 diffs)
• MSigmabar.cc (modified) (7 diffs)
• MSigmabar.h (modified) (1 diff)
• MSigmabarCalc.cc (modified) (5 diffs)
• MSigmabarCalc.h (modified) (2 diffs)

trunk/MagicSoft/Mars/manalysis/MPad.cc

@@ -64 +64 @@
 #include <TCanvas.h>
 #include <TFile.h>
+#include <TVirtualPad.h>
 
 #include "MBinning.h"
 #include "MSigmabar.h"
-#include "MMcEvt.hxx"
+#include "MPointingPos.h"
 #include "MLog.h"
 #include "MLogManip.h"
@@ -75 +76 @@
 #include "MCerPhotPix.h"
 #include "MCerPhotEvt.h"
+
+#include "MH.h"
 
 #include "MPedPhotCam.h"
@@ -190 +193 @@
   //             from sigma_k to sigma_j
 
+  //*fLog << all << "MPad::Merge2Distributions(); Entry" << endl;
+
   Double_t eps = 1.e-10;
 
@@ -205 +210 @@
   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) 
@@ -255 +269 @@
       }
 
-      *fLog << all << "k, j, arg = " << k << ",  " << j 
-            << ",  " << arg << endl;
+      //*fLog << all << "k, j, arg = " << k << ",  " << j 
+      //      << ",  " << arg << endl;
 
       //......................................
@@ -302 +316 @@
 
   // check results 
+
+  //*fLog << "Content of histap, histbp, histcp : " << endl;
+
   for (Int_t j=1; j<=nbinssig; j++)
   {
@@ -308 +325 @@
     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 = "
@@ -314 +334 @@
 
   //---------------------------------------------------------------
-  TCanvas &c = *(MH::MakeDefCanvas("Merging", "", 600, 600)); 
-  c.Divide(2, 2);
+  TCanvas *pad = MH::MakeDefCanvas("Merging", "", 600, 600); 
   gROOT->SetSelectedPad(NULL);
 
-  c.cd(1);
+  pad->Divide(2, 2);
+
+  pad->cd(1);
+  gPad->SetBorderMode(0);
   hista->SetDirectory(NULL);
-  hista->DrawCopy();
+  hista->Draw();
   hista->SetBit(kCanDelete);    
 
-  c.cd(2);
+  pad->cd(2);
+  gPad->SetBorderMode(0);
   histb->SetDirectory(NULL);
-  histb->DrawCopy();
+  histb->Draw();
   histb->SetBit(kCanDelete);    
 
-  c.cd(3);
+  pad->cd(3);
+  gPad->SetBorderMode(0);
   histap->SetDirectory(NULL);
-  histap->DrawCopy();
+  histap->Draw();
   histap->SetBit(kCanDelete);    
+
+  pad->DrawClone();
 
   //--------------------------------------------------------------------
@@ -339 +365 @@
   delete histcp;
 
+  //*fLog << all << "MPad::Merge2Distributions(); Exit" << endl;
+
   return kTRUE;
 }
@@ -357 +385 @@
 
 Bool_t MPad::MergeONOFFMC(
-   TH2D *sigthmc,  TH3D *diffpixthmc,  TH3D *sigmapixthmc,
-   TH2D *blindidthmc,  TH2D *blindnthmc,
-   TH2D *sigthon,  TH3D *diffpixthon,  TH3D *sigmapixthon,
-   TH2D *blindidthon,  TH2D *blindnthon,
-   TH2D *sigthoff, TH3D *diffpixthoff, TH3D *sigmapixthoff,
-   TH2D *blindidthoff, TH2D *blindnthoff)
+   TH2D& sigthmc,  TH3D& diffpixthmc,  TH3D& sigmapixthmc,
+   TH2D& blindidthmc,  TH2D& blindnthmc,
+   TH2D& sigthon,  TH3D& diffpixthon,  TH3D& sigmapixthon,
+   TH2D& blindidthon,  TH2D& blindnthon,
+   TH2D& sigthoff, TH3D& diffpixthoff, TH3D& sigmapixthoff,
+   TH2D& blindidthoff, TH2D& blindnthoff)
 {
+  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+  // 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( (TH2D&)*sigthon );
+  fHSigmaTheta = new TH2D( (const TH2D&) sigthon );
   fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)");
 
-  fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon );
+  *fLog << "after booking fHSigmaTheta" << endl;
+
+  fHDiffPixTheta = new TH3D( (const TH3D&) diffpixthon );
   fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)");
 
-  fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon );
+  *fLog << "after booking fHDiffPixTheta" << endl;
+
+  fHSigmaPixTheta = new TH3D( (const TH3D&) sigmapixthon );
   fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)");
-
-  fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon );
+  *fLog << "after booking fHSigmaPixTheta" << endl;
+
+  fHBlindPixNTheta = new TH2D( (const TH2D&) blindnthon );
   fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)");
 
-  fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon );
+  *fLog << "after booking fHBlindPixNTheta" << endl;
+
+  fHBlindPixIdTheta = new TH2D( (const TH2D&) blindidthon );
   fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)");
 
+  *fLog << "after booking fHBlindPixIdTheta" << endl;
+  *fLog << all << "Histograms for the merged padding plots were booked" 
+        << endl;
+
   //--------------------------
-  fHSigmaThetaMC = sigthmc;
+  fHSigmaThetaMC = &sigthmc;
   fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)");
-  fHSigmaThetaON = sigthon;
+  fHSigmaThetaON = &sigthon;
   fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)");
-  fHSigmaThetaOFF = sigthoff;
+  fHSigmaThetaOFF = &sigthoff;
   fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (orig.)");
 
+  //*fLog << all << "addresses of SigmaTheta padding plots were copied" 
+  //      << endl;
+
   //--------------------------
-  fHBlindPixNThetaMC = blindnthmc;
+  fHBlindPixNThetaMC = &blindnthmc;
   fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)");
-  fHBlindPixNThetaON = blindnthon;
+  fHBlindPixNThetaON = &blindnthon;
   fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)");
-  fHBlindPixNThetaOFF = blindnthoff;
+  fHBlindPixNThetaOFF = &blindnthoff;
   fHBlindPixNThetaOFF->SetNameTitle("2D-ThetaBlindNOFF", "2D-ThetaBlindNOFF (orig.)");
 
+  //*fLog << all << "addresses of BlindPixTheta padding plots were copied" 
+  //      << endl;
+
   //--------------------------
-  fHBlindPixIdThetaMC = blindidthmc;
+  fHBlindPixIdThetaMC = &blindidthmc;
   fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)");
-  fHBlindPixIdThetaON = blindidthon;
+  fHBlindPixIdThetaON = &blindidthon;
   fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)");
-  fHBlindPixIdThetaOFF = blindidthoff;
+  fHBlindPixIdThetaOFF = &blindidthoff;
   fHBlindPixIdThetaOFF->SetNameTitle("2D-ThetaBlindIdOFF", "2D-ThetaBlindIdOFF (orig.)");
 
+  //*fLog << all << "addresses of BlindPixIdTheta padding plots were copied" 
+  //      << endl;
+
   //--------------------------
-  fHDiffPixThetaMC = diffpixthmc;
+  fHDiffPixThetaMC = &diffpixthmc;
   fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)");
-  fHDiffPixThetaON = diffpixthon;
+  fHDiffPixThetaON = &diffpixthon;
   fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)");
-  fHDiffPixThetaOFF = diffpixthoff;
+  fHDiffPixThetaOFF = &diffpixthoff;
   fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (orig.)");
 
+  //*fLog << all << "addresses of DiffPixTheta padding plots were copied" 
+  //      << endl;
+
   //--------------------------
-  fHSigmaPixThetaMC = sigmapixthmc;
+  fHSigmaPixThetaMC = &sigmapixthmc;
   fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)");
-  fHSigmaPixThetaON = sigmapixthon;
+  fHSigmaPixThetaON = &sigmapixthon;
   fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)");
-  fHSigmaPixThetaOFF = sigmapixthoff;
+  fHSigmaPixThetaOFF = &sigmapixthoff;
   fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (orig.)");
+
+  //*fLog << all << "addresses of SigmaPixTheta padding plots were copied" 
+  //      << endl;
+
   //--------------------------
 
@@ -427 +511 @@
   // get binning for fHgON, fHgOFF and fHgMC  from sigthon
   // binning in Theta
-  TAxis *ax = sigthon->GetXaxis();
+  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)
+  //        fHBlindPixNTheta, fHBlindPixIdTheta
+  //               (they are calculated as 
+  //               the sum 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;
+
+    TH1D *hist1;
+    TH1D *hist1OFF;
+
+    // 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;
+
+
+    //------------------------------------------------------------------
+    // define target distribution  'number of blind pixels per event'
+    ay = blindnthon.GetYaxis();
+    Int_t nbinsn = ay->GetNbins();
+
+    hist1    = fHBlindPixNThetaON->ProjectionY("", j, j, "");
+    hist1->SetName("dummy");
+    hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, "");
+
+    normON  = hist1->Integral();
+    normOFF = hist1OFF->Integral();
+    if (normON  != 0.0) hist1->Scale(1.0/normON);
+    if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF);
+
+    // sum of ON and OFF distributions
+    hist1->Add(hist1OFF, 1.0); 
+    Stat_t sum1 = hist1->Integral();
+    if (sum1 != 0.0) hist1->Scale( 1.0/sum1 );
+
+    for (Int_t k=1; k<=nbinsn; k++)
+      {
+        Double_t cont = hist1->GetBinContent(k);
+        fHBlindPixNTheta->SetBinContent(j, k, cont);  
+      }
+
+    delete hist1;
+    delete hist1OFF;
+
+    //------------------------------------------------------------------
+    // define target distribution  'id of blind pixel'
+    ay = blindidthon.GetYaxis();
+    Int_t nbinsid = ay->GetNbins();
+
+    hist1    = fHBlindPixIdThetaON->ProjectionY("", j, j, "");
+    hist1->SetName("dummy");
+    hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, "");
+
+    // divide by the number of events (from fHBlindPixNTheta)
+    if (normON  != 0.0) hist1->Scale(1.0/normON);
+    if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF);
+
+    // sum of ON and OFF distributions
+    hist1->Add(hist1OFF, 1.0); 
+
+    for (Int_t k=1; k<=nbinsid; k++)
+      {
+        Double_t cont = hist1->GetBinContent(k);
+        fHBlindPixIdTheta->SetBinContent(j, k, cont);  
+      }
+
+    delete hist1;
+    delete hist1OFF;
+  }
+  //............   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);
+
+  fHBlindPixIdThetaMC->SetDirectory(NULL);
+  fHBlindPixIdThetaON->SetDirectory(NULL);
+  fHBlindPixIdThetaOFF->SetDirectory(NULL);
+
+  fHBlindPixNThetaMC->SetDirectory(NULL);
+  fHBlindPixNThetaON->SetDirectory(NULL);
+  fHBlindPixNThetaOFF->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; i<k; i++)
+    {
+      Double_t a = fHga->GetBinContent(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& blindidthmc,  TH2D& blindnthmc,
+   TH2D& sigthon,  TH3D& diffpixthon,  TH3D& sigmapixthon,
+   TH2D& blindidthon,  TH2D& blindnthon)
+{
+  *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)");
+
+  fHBlindPixNTheta = new TH2D( (TH2D&) blindnthon );
+  fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)");
+
+  fHBlindPixIdTheta = new TH2D( (TH2D&) blindidthon );
+  fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)");
+
+  //--------------------------
+  fHSigmaThetaMC = &sigthmc;
+  fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)");
+  fHSigmaThetaON = &sigthon;
+  fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)");
+
+  //--------------------------
+  fHBlindPixNThetaMC = &blindnthmc;
+  fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)");
+  fHBlindPixNThetaON = &blindnthon;
+  fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)");
+
+  //--------------------------
+  fHBlindPixIdThetaMC = &blindidthmc;
+  fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)");
+  fHBlindPixIdThetaON = &blindidthon;
+  fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (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;
@@ -441 +1073 @@
 
   // binning in sigmabar
-  TAxis *ay = sigthon->GetYaxis();
+  TAxis *ay = sigthon.GetYaxis();
   Int_t nbinssig = ay->GetNbins();
   TArrayD edgesy;
@@ -463 +1095 @@
   fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON");
 
-  fHgOFF = new TH3D;
-  MH::SetBinning(fHgOFF, &binth, &binsg, &binsg);
-  fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF");
-
-  //............   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++)
-  {
-    TH2D * fHga = new TH2D;
-    MH::SetBinning(fHga, &binsg, &binsg);
-    TH2D * fHgb = new TH2D;
-    MH::SetBinning(fHgb, &binsg, &binsg);
-
-    //-------------------------------------------
-    // 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
-
-    TH1D *hista  = sigthon->ProjectionY("sigon_y", l, l, "");
-    Stat_t suma = hista->Integral();
-    hista->Scale(1./suma);
-
-    TH1D *histap  = new TH1D( (const TH1D&)*hista );
-
-    TH1D *histb  = sigthoff->ProjectionY("sigoff_y", l, l, "");
-    Stat_t sumb = histb->Integral();
-    histb->Scale(1./sumb);
-
-    Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig);
-    delete hista;
-    delete histb;
-
-    // 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);
-        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);
-        fHgOFF->SetBinContent(l, k, j, b);
-      }
-
-    //-------------------------------------------
-    // now merge ON-OFF and MC distributions
-    // input : hista is the result of the merging of ON and OFF distributions  
-    //         histb is the normalized 1D distr. of sigmabar for MC data
-    // output : histap    will be the merged distribution (target distribution)
-    //          fHga(k,j) will tell which fraction of ON, OFF 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
-
-    hista  = new TH1D( (const TH1D&)*histap);
-
-    histb  = sigthmc->ProjectionY("sigmc_y", l, l, "");
-    sumb = histb->Integral();
-    histb->Scale(1./sumb);
-
-    Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig);
-    delete hista;
-    delete histb;
-
-    // update fHgON and fHgOFF and fill fHgMC
-    for (Int_t k=1; k<=nbinssig; k++)
-      for (Int_t j=1; j<=nbinssig; j++)
-      {
-        Double_t a   =  fHga->GetBinContent  (k,j);
-
-        Double_t aON = fHgON->GetBinContent(l,k,j);
-        fHgON->SetBinContent(l, k, j, aON+a);
-
-        Double_t aOFF = fHgOFF->GetBinContent(l,k,j);
-        fHgOFF->SetBinContent(l, k, j, aOFF+a);
-
-        Double_t b = fHgb->GetBinContent(k,j);
-        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 OFF distributions)
-  //        fHBlindPixNTheta, fHBlindPixIdTheta
-  //               (they are calculated as 
-  //               the sum 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;
-
-    TH1D *hist1;
-    TH1D *hist1OFF;
-
-    // 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, 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;
-
-
-    //------------------------------------------------------------------
-    // define target distribution  'number of blind pixels per event'
-    ay = blindnthon->GetYaxis();
-    Int_t nbinsn = ay->GetNbins();
-
-    hist1    = fHBlindPixNThetaON->ProjectionY("", j, j, "");
-    hist1->SetName("dummy");
-    hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, "");
-
-    normON  = hist1->Integral();
-    normOFF = hist1OFF->Integral();
-    if (normON  != 0.0) hist1->Scale(1.0/normON);
-    if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF);
-
-    // sum of ON and OFF distributions
-    hist1->Add(hist1OFF, 1.0); 
-    hist1->Scale( 1.0/hist1->Integral() );
-
-    for (Int_t k=1; k<=nbinsn; k++)
-      {
-        Double_t cont = hist1->GetBinContent(k);
-        fHBlindPixNTheta->SetBinContent(j, k, cont);  
-      }
-
-    delete hist1;
-    delete hist1OFF;
-
-    //------------------------------------------------------------------
-    // define target distribution  'id of blind pixel'
-    ay = blindidthon->GetYaxis();
-    Int_t nbinsid = ay->GetNbins();
-
-    hist1    = fHBlindPixIdThetaON->ProjectionY("", j, j, "");
-    hist1->SetName("dummy");
-    hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, "");
-
-    // divide by the number of events (from fHBlindPixNTheta)
-    if (normON  != 0.0) hist1->Scale(1.0/normON);
-    if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF);
-
-    // sum of ON and OFF distributions
-    hist1->Add(hist1OFF, 1.0); 
-
-    for (Int_t k=1; k<=nbinsid; k++)
-      {
-        Double_t cont = hist1->GetBinContent(k);
-        fHBlindPixIdTheta->SetBinContent(j, k, cont);  
-      }
-
-    delete hist1;
-    delete hist1OFF;
-  }
-  //............   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);
-
-  fHBlindPixIdThetaMC->SetDirectory(NULL);
-  fHBlindPixIdThetaON->SetDirectory(NULL);
-  fHBlindPixIdThetaOFF->SetDirectory(NULL);
-
-  fHBlindPixNThetaMC->SetDirectory(NULL);
-  fHBlindPixNThetaON->SetDirectory(NULL);
-  fHBlindPixNThetaOFF->SetDirectory(NULL);
-
-  fHgMC->SetDirectory(NULL);
-  fHgON->SetDirectory(NULL);
-  fHgOFF->SetDirectory(NULL);
-
-  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 *blindidthmc,  TH2D *blindnthmc,
-   TH2D *sigthon,  TH3D *diffpixthon,  TH3D *sigmapixthon,
-   TH2D *blindidthon,  TH2D *blindnthon)
-{
-  *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)");
-
-  fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon );
-  fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)");
-
-  fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon );
-  fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)");
-
-  //--------------------------
-  fHSigmaThetaMC = sigthmc;
-  fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)");
-  fHSigmaThetaON = sigthon;
-  fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)");
-
-  //--------------------------
-  fHBlindPixNThetaMC = blindnthmc;
-  fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)");
-  fHBlindPixNThetaON = blindnthon;
-  fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)");
-
-  //--------------------------
-  fHBlindPixIdThetaMC = blindidthmc;
-  fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)");
-  fHBlindPixIdThetaON = blindidthon;
-  fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (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   ...........................
@@ -893 +1118 @@
     //                    padded from bin k to bin j
 
-    TH1D *hista  = sigthon->ProjectionY("sigon_y", l, l, "");
+    TH1D *hista  = sigthon.ProjectionY("sigon_y", l, l, "");
     Stat_t suma = hista->Integral();
-    hista->Scale(1./suma);
+    if (suma != 0.0) hista->Scale(1./suma);
 
     TH1D *histap  = new TH1D( (const TH1D&)*hista );
 
-    TH1D *histb  = sigthmc->ProjectionY("sigmc_y", l, l, "");
+    TH1D *histb  = sigthmc.ProjectionY("sigmc_y", l, l, "");
     Stat_t sumb = histb->Integral();
-    histb->Scale(1./sumb);
+    if (sumb != 0.0) histb->Scale(1./sumb);
 
     Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig);
@@ -958 +1183 @@
 
     // set ranges for 2D-projections of 3D-histograms
-    TAxis *ax = diffpixthon->GetXaxis();
+    TAxis *ax = diffpixthon.GetXaxis();
     ax->SetRange(j, j);
-    ax = diffpixthmc->GetXaxis();
+    ax = diffpixthmc.GetXaxis();
     ax->SetRange(j, j);
 
@@ -980 +1205 @@
     //------------------------------------------------------------------
     // define target distribution  'sigma^2-sigmavar^2 vs. pixel number'
-    ay = diffpixthon->GetYaxis();
+    ay = diffpixthon.GetYaxis();
     Int_t nbinspixel = ay->GetNbins();
 
-    TAxis *az = diffpixthon->GetZaxis();
+    TAxis *az = diffpixthon.GetZaxis();
     Int_t nbinsdiff = az->GetNbins();
 
-    hist    = (TH2D*)diffpixthon->Project3D("zy");
+    hist    = (TH2D*)diffpixthon.Project3D("zy");
     hist->SetName("dummy");
-    histMC = (TH2D*)diffpixthmc->Project3D("zy");
+    histMC = (TH2D*)diffpixthmc.Project3D("zy");
 
     normON  = hist->Integral();
@@ -1011 +1236 @@
     //------------------------------------------------------------------
     // define target distribution  'sigma vs. pixel number'
-    ay = sigmapixthon->GetYaxis();
+    ay = sigmapixthon.GetYaxis();
     nbinspixel = ay->GetNbins();
 
-    az = sigmapixthon->GetZaxis();
+    az = sigmapixthon.GetZaxis();
     Int_t nbinssigma = az->GetNbins();
 
-    hist    = (TH2D*)sigmapixthon->Project3D("zy");
+    hist    = (TH2D*)sigmapixthon.Project3D("zy");
     hist->SetName("dummy");
-    histMC = (TH2D*)sigmapixthmc->Project3D("zy");
+    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);
+    if (normMC  != 0.0) histMC->Scale(1.0/normMC);
 
     // weighted average of ON and MC distributions
@@ -1043 +1268 @@
     //------------------------------------------------------------------
     // define target distribution  'number of blind pixels per event'
-    ay = blindnthon->GetYaxis();
+    ay = blindnthon.GetYaxis();
     Int_t nbinsn = ay->GetNbins();
 
@@ -1053 +1278 @@
     normMC  = hist1MC->Integral();
     if (normON  != 0.0) hist1->Scale(1.0/normON);
-    if (normMC != 0.0)  hist1MC->Scale(1.0/normMC);
+    if (normMC  != 0.0) hist1MC->Scale(1.0/normMC);
 
     // sum of ON and MC distributions
     hist1->Add(hist1MC, 1.0); 
-    hist1->Scale( 1.0/hist1->Integral() );
+    Stat_t sum1 = hist1->Integral();
+    if (sum1 != 0.0) hist1->Scale( 1.0/sum1 );
 
     for (Int_t k=1; k<=nbinsn; k++)
@@ -1070 +1296 @@
     //------------------------------------------------------------------
     // define target distribution  'id of blind pixel'
-    ay = blindidthon->GetYaxis();
+    ay = blindidthon.GetYaxis();
     Int_t nbinsid = ay->GetNbins();
 
@@ -1079 +1305 @@
     // divide by the number of events (from fHBlindPixNTheta)
     if (normON  != 0.0) hist1->Scale(1.0/normON);
-    if (normMC != 0.0)  hist1MC->Scale(1.0/normMC);
+    if (normMC  != 0.0) hist1MC->Scale(1.0/normMC);
 
     // sum of ON and MC distributions
@@ -1138 +1364 @@
 
     //------------------------------------
+  
+      /*
       fHSigmaThetaMC = 
       (TH2D*) gROOT->FindObject("2D-ThetaSigmabarMC");
@@ -1147 +1375 @@
           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");
@@ -1159 +1391 @@
           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");
@@ -1171 +1407 @@
           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");
@@ -1184 +1424 @@
           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");
@@ -1196 +1440 @@
           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");
@@ -1208 +1456 @@
           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");
@@ -1220 +1472 @@
           return kFALSE;
         }
+      */
+      fHDiffPixThetaOFF = new TH3D;
+      fHDiffPixThetaOFF->Read("3D-ThetaPixDiffOFF");
       *fLog << all 
             << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl;
@@ -1225 +1480 @@
 
     //------------------------------------
+      /*
        fHSigmaPixTheta = 
       (TH3D*) gROOT->FindObject("3D-ThetaPixSigma");
@@ -1234 +1490 @@
           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");
@@ -1246 +1506 @@
           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");
@@ -1258 +1522 @@
           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");
@@ -1270 +1538 @@
           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");
@@ -1283 +1555 @@
           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");
@@ -1295 +1571 @@
           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");
@@ -1307 +1587 @@
           return kFALSE;
         }
+      */
+      fHgOFF = new TH3D;
+      fHgOFF->Read("3D-PaddingMatrixOFF");
       *fLog << all 
             << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl;
 
     //------------------------------------
+      /*
+      fHBlindPixIdTheta = 
+      (TH2D*) gROOT->FindObject("2D-ThetaBlindId");
+      if (!fHBlindPixIdTheta)
+        {
+          *fLog << all 
+                << "MPad : Object '2D-ThetaBlindId' not found on root file" 
+                << endl;
+          return kFALSE;
+        }
+      */
+      fHBlindPixIdTheta = new TH2D;
+      fHBlindPixIdTheta->Read("2D-ThetaBlindId");
+      *fLog << all 
+            << "MPad : Object '2D-ThetaBlindId' was read in" << endl;
+
+      /*
       fHBlindPixIdThetaMC = 
       (TH2D*) gROOT->FindObject("2D-ThetaBlindIdMC");
@@ -1320 +1620 @@
           return kFALSE;
         }
+      */
+      fHBlindPixIdThetaMC = new TH2D;
+      fHBlindPixIdThetaMC->Read("2D-ThetaBlindIdMC");
       *fLog << all 
             << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl;
 
+      /*
       fHBlindPixIdThetaON = 
       (TH2D*) gROOT->FindObject("2D-ThetaBlindIdON");
@@ -1332 +1636 @@
           return kFALSE;
         }
+      */
+      fHBlindPixIdThetaON = new TH2D;
+      fHBlindPixIdThetaON->Read("2D-ThetaBlindIdON");
       *fLog << all 
             << "MPad : Object '2D-ThetaBlindIdON' was read in" << endl;
 
+      /*
       fHBlindPixIdThetaOFF = 
       (TH2D*) gROOT->FindObject("2D-ThetaBlindIdOFF");
@@ -1344 +1652 @@
           return kFALSE;
         }
+      */
+      fHBlindPixIdThetaOFF = new TH2D;
+      fHBlindPixIdThetaOFF->Read("2D-ThetaBlindIdOFF");
       *fLog << all 
-            << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl;
+            << "MPad : Object '2D-ThetaBlindIdOFF' was read in" << endl;
 
     //------------------------------------
+      /*
+      fHBlindPixNTheta = 
+      (TH2D*) gROOT->FindObject("2D-ThetaBlindN");
+      if (!fHBlindPixNTheta)
+        {
+          *fLog << all 
+                << "MPad : Object '2D-ThetaBlindN' not found on root file" 
+                << endl;
+          return kFALSE;
+        }
+      */
+      fHBlindPixNTheta = new TH2D;
+      fHBlindPixNTheta->Read("2D-ThetaBlindN");
+      *fLog << all 
+            << "MPad : Object '2D-ThetaBlindN' was read in" << endl;
+
+      /*
       fHBlindPixNThetaMC = 
       (TH2D*) gROOT->FindObject("2D-ThetaBlindNMC");
@@ -1357 +1685 @@
           return kFALSE;
         }
+      */
+      fHBlindPixNThetaMC = new TH2D;
+      fHBlindPixNThetaMC->Read("2D-ThetaBlindNMC");
       *fLog << all 
             << "MPad : Object '2D-ThetaBlindNMC' was read in" << endl;
 
+      /*
       fHBlindPixNThetaON = 
       (TH2D*) gROOT->FindObject("2D-ThetaBlindNON");
@@ -1369 +1701 @@
           return kFALSE;
         }
+      */
+      fHBlindPixNThetaON = new TH2D;
+      fHBlindPixNThetaON->Read("2D-ThetaBlindNON");
       *fLog << all 
             << "MPad : Object '2D-ThetaBlindNON' was read in" << endl;
 
+      /*
       fHBlindPixNThetaOFF = 
       (TH2D*) gROOT->FindObject("2D-ThetaBlindNOFF");
@@ -1381 +1717 @@
           return kFALSE;
         }
+      */
+      fHBlindPixNThetaOFF = new TH2D;
+      fHBlindPixNThetaOFF->Read("2D-ThetaBlindNOFF");
       *fLog << all 
             << "MPad : Object '2D-ThetaBlindNOFF' was read in" << endl;
@@ -1402 +1741 @@
   TFile outfile(namefileout, "RECREATE");
 
+  fHBlindPixNTheta->Write();
   fHBlindPixNThetaMC->Write();
   fHBlindPixNThetaON->Write();
   fHBlindPixNThetaOFF->Write();
 
+  fHBlindPixIdTheta->Write();
   fHBlindPixIdThetaMC->Write();
   fHBlindPixIdThetaON->Write();
@@ -1470 +1811 @@
   }
 
-  fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
-  if (!fMcEvt)
+  fPointPos = (MPointingPos*)pList->FindObject("MPointingPos");
+  if (!fPointPos)
     {
-       *fLog << err << dbginf << "MPad : MMcEvt not found... aborting." 
+       *fLog << err << dbginf << "MPad : MPointingPos not found... aborting." 
              << endl;
        return kFALSE;
-     }
+    }
   
    fPed = (MPedPhotCam*)pList->FindObject("MPedPhotCam");
@@ -1510 +1851 @@
      }
 
-   fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
-   if (!fBlinds)
+   fBlindPix = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
+   if (!fBlindPix)
      {
        *fLog << err << "MPad : MBlindPixels not found... aborting." 
@@ -1532 +1873 @@
    // plot pedestal sigmas
    fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", 
-                     100, 0.0, 5.0, 100, 0.0, 5.0);
+                     100, 0.0, 75.0, 100, 0.0, 75.0);
    fHSigmaPedestal->SetXTitle("Pedestal sigma before padding");
    fHSigmaPedestal->SetYTitle("Pedestal sigma after padding");
@@ -1538 +1879 @@
    // plot no.of photons (before vs. after padding) 
    fHPhotons = new TH2D("Photons","Photons: after vs.before padding", 
-                        100, -10.0, 90.0, 100, -10, 90);
+                        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, -10.0, 10.0);
+   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");
@@ -1550 +1891 @@
    //--------------------------------------------------------------------
 
-   fIter = 20;
+   fIter = 10;
 
    memset(fErrors,   0, sizeof(fErrors));
@@ -1570 +1911 @@
 Int_t MPad::Process()
 {
-  *fLog << all << "Entry MPad::Process();" << endl;
+  //*fLog << all << "Entry MPad::Process();" << endl;
 
   // this is the conversion factor from the number of photons
@@ -1576 +1917 @@
   // later to be taken from MCalibration
   // fPEperPhoton = PW * LG * QE * 1D
-  Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9;
+  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
-
+  // 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();
 
@@ -1603 +1950 @@
     fEvt->AddPixel(i, 0.0, (*fPed)[i].GetRms());
   }
-
+  */
 
 
@@ -1615 +1962 @@
   // Calculate average sigma of the event
   //
-  Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt);
-  *fLog << all << "MPad::Process(); before padding : " << endl;
-  fSigmabar->Print("");
-  Double_t sigbarold  = sigbarold_phot * fPEperPhoton;
-  Double_t sigbarold2 = sigbarold*sigbarold;
-
-  const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
-  *fLog << all << "theta = " << theta << endl;
+
+  fSigmabar->Calc(*fCam, *fPed, *fEvt);
+  //*fLog << all << "MPad::Process(); before padding : " << endl;
+  //fSigmabar->Print("");
+
+  // inner sigmabar/sqrt(area)
+  Double_t sigbarInnerold_phot = fSigmabar->GetSigmabarInner();
+  Double_t sigbarInnerold  = sigbarInnerold_phot * fPEperPhoton;
+  Double_t sigbarInnerold2 = sigbarInnerold*sigbarInnerold;
+
+  // outer sigmabar/sqrt(area)
+  Double_t sigbarOuterold_phot = fSigmabar->GetSigmabarOuter();
+  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 = fHBlindPixNTheta->GetXaxis()->FindBin(theta);
@@ -1629 +1986 @@
     *fLog << warn 
           << "MPad::Process(); binNumber out of range : theta, binTheta = "
-          << theta << ",  " << binTheta << ";  Skip event " << endl;
-    // event cannot be padded; skip event
+          << theta << ",  " << binTheta << ";  aborting " << endl;
 
     rc = 2;
@@ -1637 +1993 @@
   }
 
+  //*fLog << all << "binTheta = " << binTheta << endl;
+  
+
   //-------------------------------------------
   // get number of events in this theta bin
-  // and number of events in this sigbarold_phot bin
+  // and number of events in this sigbarInnerold_phot bin
 
   Double_t nTheta;
@@ -1657 +2016 @@
   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(sigbarold_phot);
+  Int_t binSigma = hn->FindBin(sigbarInnerold_phot);
   nSigma = hn->GetBinContent(binSigma);
 
-   *fLog << all           
-         << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
-         << theta << ",  " << sigbarold_phot << ",  " << binTheta << ",  "
-         << binSigma << ",  " << nTheta << ",  " << nSigma   << endl;      
+  //*fLog << all           
+  //       << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
+  //       << theta << ",  " << sigbarInnerold_phot << ",  " << binTheta 
+  //       << ",  "
+  //       << binSigma << ",  " << nTheta << ",  " << nSigma   << endl;      
 
   delete hn;
-  
-
+ 
   //-------------------------------------------
   // for the current theta :
@@ -1688 +2055 @@
     if ( nblind->Integral() == 0.0 )
     {
-      *fLog << warn << "MPad::Process(); projection for Theta bin " 
-            << binTheta << " has no entries; Skip event " << endl;
+      *fLog << warn << "MPad::Process(); projection of '"
+            << fHBlindPixNThetaOFF->GetName() << "' for Theta bin " 
+            << binTheta << " has no entries; aborting " << endl;
       // event cannot be padded; skip event
       delete nblind;
@@ -1732 +2100 @@
           nlist++;
 
-          fBlinds->SetPixelBlind(idBlind);
-          *fLog << all << "idBlind = " << idBlind << endl;
+          fBlindPix->SetPixelBlind(idBlind);
+          //*fLog << all << "idBlind = " << idBlind << endl;
         }
-      fBlinds->SetReadyToSave();
+      fBlindPix->SetReadyToSave();
 
       delete hblind;
@@ -1752 +2120 @@
     if ( nblind->Integral() == 0.0 )
     {
-      *fLog << warn << "MPad::Process(); projection for Theta bin " 
+      *fLog << warn << "MPad::Process(); projection of '" 
+            << fHBlindPixNThetaON->GetName() << "' for Theta bin " 
             << binTheta << " has no entries; Skip event " << endl;
-      // event cannot be padded; skip event
       delete nblind;
 
@@ -1796 +2164 @@
           nlist++;
 
-          fBlinds->SetPixelBlind(idBlind);
-          *fLog << all << "idBlind = " << idBlind << endl;
+          fBlindPix->SetPixelBlind(idBlind);
+          //*fLog << all << "idBlind = " << idBlind << endl;
         }
-      fBlinds->SetReadyToSave();
+      fBlindPix->SetReadyToSave();
 
       delete hblind;
@@ -1810 +2178 @@
   //
 
-  Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot);
-  *fLog << all << "binSig, sigbarold_phot = " << binSig << ",  " 
-          << sigbarold_phot << endl;
+  Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarInnerold_phot);
+  //*fLog << all << "binSig, sigbarInnerold_phot = " << binSig << ",  " 
+  //        << sigbarInnerold_phot << endl;
 
   Double_t prob;
@@ -1843 +2211 @@
     prob = 0.0;
 
-   *fLog << all << "nTheta, nSigma, prob = " << nTheta << ",  " << nSigma 
-         << ",  " << prob << endl;
+  //*fLog << all << "nTheta, nSigma, prob = " << nTheta << ",  " << nSigma 
+  //       << ",  " << prob << endl;
 
   if ( prob <= 0.0  ||  gRandom->Uniform() > prob )
@@ -1850 +2218 @@
     delete hpad;
     // event should not be padded
-    *fLog << all << "event is not padded" << endl;
+    //*fLog << all << "event is not padded" << endl;
 
     rc = 8;
@@ -1857 +2225 @@
   }
   // event should be padded
-  *fLog << all << "event will be padded" << endl;  
+  //*fLog << all << "event will be padded" << endl;  
 
 
@@ -1864 +2232 @@
   //     for MC/ON/OFF data according to the matrix fHgMC/ON/OFF
   //
-  Double_t sigmabar_phot = 0;
-  Double_t sigmabar      = 0;
+  Double_t sigmabarInner_phot = 0;
+  Double_t sigmabarInner      = 0;
 
   
@@ -1876 +2244 @@
   //}
 
-  sigmabar_phot = hpad->GetRandom();
-  sigmabar = sigmabar_phot * fPEperPhoton;
-
-  *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl;
+  sigmabarInner_phot = hpad->GetRandom();
+  sigmabarInner = sigmabarInner_phot * fPEperPhoton;
+
+  //*fLog << all << "sigmabarInner_phot = " << sigmabarInner_phot << endl;
 
   delete hpad;
   
-  const Double_t sigmabar2 = sigmabar*sigmabar;
+  // new inner sigmabar2/area
+  const Double_t sigmabarInner2 = sigmabarInner*sigmabarInner;
 
   //-------------------------------------------
 
-  *fLog << all << "MPad::Process(); sigbarold, sigmabar = " 
-        << sigbarold << ",  "<< sigmabar << endl;
-
-  // Skip event if target sigmabar is <= sigbarold
-  if (sigmabar <= sigbarold)
+  //*fLog << all << "MPad::Process(); sigbarInnerold, sigmabarInner = " 
+  //      << sigbarInnerold << ",  "<< sigmabarInner << endl;
+
+  // Stop if target sigmabar is <= sigbarold
+  if (sigmabarInner <= sigbarInnerold)
   {
-    *fLog << all << "MPad::Process(); target sigmabar is less than sigbarold : "
-          << sigmabar << ",  " << sigbarold << ",   Skip event" << endl;
+    *fLog << all << "MPad::Process(); target sigmabarInner is less than sigbarInnerold : "
+          << sigmabarInner << ",  " << sigbarInnerold << ",   aborting" 
+          << endl;
 
     rc = 4;
@@ -1923 +2293 @@
   }
 
-  // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin.
-  // The average electronic noise (to be added) has to be well above this RMS,
-  // otherwise the electronic noise of an individual pixel (elNoise2Pix)
-  // may become negative
+  // 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, 9999, "");
+    hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, "");
   else if (fType == "ON")
-    hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
+    hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, "");
   else if (fType == "OFF")
-    hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
+    hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, "");
   else
   {
@@ -1946 +2321 @@
   delete hnoise;
 
-  elNoise2 = TMath::Min(RMS,  sigmabar2 - sigbarold2);
-  *fLog << all << "elNoise2 = " << elNoise2 << endl; 
-
-  lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess;  
+  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;
@@ -1957 +2332 @@
   // do the padding for each pixel
   //
-  // pad only pixels   - which are used (before image cleaning)
+  // pad only pixels   - which are used and not blind (before image cleaning)
   //
 
@@ -1971 +2346 @@
 
 
+  // 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; i<npix; i++) 
   {   
+
     MCerPhotPix &pix = (*fEvt)[i];
     if ( !pix.IsPixelUsed() )
@@ -1990 +2388 @@
     Double_t ratioArea = 1.0 / fCam->GetPixRatio(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 ( fBlindPix != NULL  &&  fBlindPix->IsBlind(j) )
+    {
+       // this should never occur, because blind pixels should have
+       // been set unused by MBlindPixelsCalc2::UnMap()
+       //*fLog << all << "MPad::Process; blind 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;
+    Double_t oldsigma2 = oldsigma*oldsigma / ratioArea;
 
     //---------------------------------
@@ -2005 +2430 @@
     TH1D *hdiff;
 
-    // throw the Sigma for this pixel 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;
-    }
-
-    hdiff = sp->ProjectionZ("", binTheta, binTheta,
-                                binPixel, binPixel, "");
-
-    if ( hdiff->GetEntries() == 0 )
-    {
-      *fLog << warn << "MPad::Process(); projection for Theta bin " 
-            << binTheta << "  and pixel bin " << binPixel  
-            << " has no entries;  aborting " << endl;
-      delete hdiff;
-
-      rc = 5;
-      fErrors[rc]++;
-      return kCONTINUE;     
-    }
-
-    count = 0;
-    ok = kFALSE;
-    for (Int_t m=0; m<fIter; m++)
-    {
-      count += 1;
-      diff_phot = hdiff->GetRandom();
-      diff = diff_phot * fPEperPhoton * fPEperPhoton;
+
+     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; m<fIter; m++)
+       {
+         count += 1;
+         diff_phot = hdiff->GetRandom();
+
+         //*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/ratioArea
-                               - lambdabar*F2excess) > 0.0 )
-      {
-        ok = kTRUE;
-        break;
-      }
-    }
-
-    if (!ok)
-    {
-      *fLog << all << "theta, j, count, sigmabar, diff = " << theta 
-            << ",  " 
-            << j << ",  " << count << ",  " << sigmabar << ",  " 
-            << diff << endl;
-      diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2;
-
-      rw = 1;
-      fWarnings[rw]++;
-    }
-    else
-    {
-      rw = 0;
-      fWarnings[rw]++;
-    }
+         // 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;
@@ -2078 +2498 @@
     // get the additional sigma^2 for this pixel (due to the padding)
 
-    addSig2 = sigma2*ratioArea - oldsigma2;
+    addSig2 = (sigma2 - oldsigma2) * ratioArea;
     addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton);
 
@@ -2139 +2559 @@
 
 
-    Double_t newsigma = sqrt( oldsigma2 + addSig2 ); 
+    Double_t newsigma = sqrt( sigma2 * ratioArea ); 
     Double_t newsigma_phot = newsigma / fPEperPhoton; 
     ppix.SetRms( newsigma_phot );
@@ -2150 +2570 @@
 
   fSigmabar->Calc(*fCam, *fPed, *fEvt);
-  *fLog << all << "MPad::Process(); after padding : " << endl;
-  fSigmabar->Print("");
-
-  *fLog << all << "Exit MPad::Process();" << endl;
+  //*fLog << all << "MPad::Process(); after padding : " << endl;
+  //fSigmabar->Print("");
+
+  //*fLog << all << "Exit MPad::Process();" << endl;
 
   rc = 0;
@@ -2173 +2593 @@
     *fLog << dec << setfill(' ');
 
-    int temp;
-    if (fWarnings[0] != 0.0)    
-      temp = (int)(fWarnings[1]*100/fWarnings[0]);
-    else
-      temp = 100;
+    if (fWarnings[0] == 0 ) fWarnings[0] = 1;
 
     *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) 
-          << temp 
-          << "%) Warning : iteration to find acceptable sigma failed" 
+          << (int)(fWarnings[1]*100/fWarnings[0])
+          << "%) Pixel: iteration to find acceptable sigma failed" 
           << ", fIter = " << fIter << endl;
 
-    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
-          << (int)(fErrors[1]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: Sigmabar_old > 0" << 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) 
@@ -2192 +2609 @@
           << "%) Evts skipped due to: Zenith angle out of range" << endl;
 
-    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) 
-          << (int)(fErrors[3]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: No data for generating Sigmabar" << 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[5] << " (" << setw(3) 
-          << (int)(fErrors[5]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
-
-    *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) 
-          << (int)(fErrors[6]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: No data for generating Sigma" << endl;
 
     *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) 
@@ -2224 +2629 @@
 
     //---------------------------------------------------------------
-    TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 1500)); 
-    c.Divide(3, 5);
+    TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 600)); 
+    c.Divide(3, 2);
     gROOT->SetSelectedPad(NULL);
 
@@ -2245 +2650 @@
     //--------------------------------------------------------------------
 
-
+    /*
     c.cd(4);
     fHSigmaTheta->SetDirectory(NULL);
@@ -2251 +2656 @@
     fHSigmaTheta->DrawCopy();     
     fHSigmaTheta->SetBit(kCanDelete);     
-
+    */
 
     //--------------------------------------------------------------------
 
-
+    /*
     c.cd(7);
     fHBlindPixNTheta->SetDirectory(NULL);
@@ -2261 +2666 @@
     fHBlindPixNTheta->DrawCopy();     
     fHBlindPixNTheta->SetBit(kCanDelete);     
+    */
 
     //--------------------------------------------------------------------
 
-
+    /*
     c.cd(10);
     fHBlindPixIdTheta->SetDirectory(NULL);
@@ -2270 +2676 @@
     fHBlindPixIdTheta->DrawCopy();     
     fHBlindPixIdTheta->SetBit(kCanDelete);     
-
+    */
 
     //--------------------------------------------------------------------
     // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2
 
+    /*
     c.cd(5);
     TH2D *l1 = NULL;
@@ -2296 +2703 @@
     l2->DrawCopy("box");
     l2->SetBit(kCanDelete);;
+    */
 
     //--------------------------------------------------------------------
     // draw the 3D histogram (target): Theta, pixel, Sigma
 
+    /*
     c.cd(6);
     TH2D *k1 = NULL;
@@ -2321 +2730 @@
     k2->DrawCopy("box");
     k2->SetBit(kCanDelete);;
-
+    */
 
     //--------------------------------------------------------------------
 
-    c.cd(13);
+    
+    c.cd(4);
     TH2D *m1;
     m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy");
     m1->SetDirectory(NULL);
-    m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all  \\Theta)");
+    m1->SetTitle("(fHgON) Sigmabar-new vs. Sigmabar-old (ON, all  \\Theta)");
     m1->SetXTitle("Sigmabar-old");
     m1->SetYTitle("Sigmabar-new");
@@ -2336 +2746 @@
     m1->SetBit(kCanDelete);;
 
-    c.cd(14);
+    c.cd(5);
     TH2D *m2;
     m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy");
     m2->SetDirectory(NULL);
-    m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all  \\Theta)");
+    m2->SetTitle("(fHgOFF) Sigmabar-new vs. Sigmabar-old (OFF, all  \\Theta)");
     m2->SetXTitle("Sigmabar-old");
     m2->SetYTitle("Sigmabar-new");
@@ -2347 +2757 @@
     m2->SetBit(kCanDelete);;
 
-    c.cd(15);
+    c.cd(6);
     TH2D *m3;
     m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy");
     m3->SetDirectory(NULL);
-    m3->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (MC, all  \\Theta)");
+    m3->SetTitle("(fHgMC) Sigmabar-new vs. Sigmabar-old (MC, all  \\Theta)");
     m3->SetXTitle("Sigmabar-old");
     m3->SetYTitle("Sigmabar-new");
@@ -2369 +2779 @@
 
 
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MPad.h

@@ -17 +17 @@
 class MCerPhotEvt;
 class MPedPhotCam;
-class MMcEvt;
+class MPointingPos;
 class MSigmabar;
 class MParList;
@@ -31 +31 @@
     MCerPhotEvt    *fEvt; 
     MSigmabar      *fSigmabar;
-    MMcEvt         *fMcEvt;
+    MPointingPos   *fPointPos;
     MPedPhotCam    *fPed;
-    MBlindPixels   *fBlinds;
+    MBlindPixels   *fBlindPix;
 
     TString        fType;           // type of data to be padded
@@ -41 +41 @@
 
     Int_t          fErrors[9];
-    Int_t          fWarnings[2];
+    Int_t          fWarnings[3];
 
     //----------------------------------
@@ -89 +89 @@
 
 
-    Bool_t Merge2Distributions(TH1D *hista, TH1D * histb, TH1D * histap,
-                               TH2D *fHga,  TH2D * fHgb, Int_t nbinssig);
+    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:
@@ -98 +100 @@
 
     Bool_t MergeONOFFMC(
-      TH2D *sigthmc,  TH3D *diffpixthmc, TH3D *sigmapixthmc,
-      TH2D *blindidthmc,  TH2D *blindnthmc,
-      TH2D *sigthon,  TH3D *diffpixthon, TH3D *sigmapixthon,
-      TH2D *blindidthon,  TH2D *blindnthon,
-      TH2D *sigthoff=NULL, TH3D *diffpixthoff=NULL,TH3D *sigmapixthoff=NULL,
-      TH2D *blindidthoff=NULL, TH2D *blindnthoff=NULL);
+      TH2D& sigthmc,  TH3D& diffpixthmc, TH3D& sigmapixthmc,
+      TH2D& blindidthmc,  TH2D& blindnthmc,
+      TH2D& sigthon,  TH3D& diffpixthon, TH3D& sigmapixthon,
+      TH2D& blindidthon,  TH2D& blindnthon,
+      TH2D& sigthoff, TH3D& diffpixthoff,TH3D& sigmapixthoff,
+      TH2D& blindidthoff, TH2D& blindnthoff);
 
     Bool_t MergeONMC(
-      TH2D *sigthmc,  TH3D *diffpixthmc, TH3D *sigmapixthmc,
-      TH2D *blindidthmc,  TH2D *blindnthmc,
-      TH2D *sigthon,  TH3D *diffpixthon, TH3D *sigmapixthon,
-      TH2D *blindidthon,  TH2D *blindnthon);
+      TH2D& sigthmc,  TH3D& diffpixthmc, TH3D& sigmapixthmc,
+      TH2D& blindidthmc,  TH2D& blindnthmc,
+      TH2D& sigthon,  TH3D&diffpixthon, TH3D& sigmapixthon,
+      TH2D& blindidthon,  TH2D& blindnthon);
 
     Bool_t ReadPaddingDist(const char *filein);

trunk/MagicSoft/Mars/manalysis/MSigmabar.cc

@@ -30 +30 @@
 //                                                                         //
 // This is the storage container to hold information about                 //
-// "average" pedestal sigmas                                               //
-//                                                                         //
-// In calculating averages all sigmas are normalized to the area of pixel 0//
+// "average" pedestal RMS                                                  //
+//                                                                         //
+// The "average" pedestal RMS is calculated as                             //
+//     <pedRMS> = sqrt( sum_i( (pedRMS_i)^2/area_i ) / no.of pixels )      //
+//                                                                         //
+//     which is the sqrt of the average (pedRMS^2 per area)                //
 //                                                                         //
 /////////////////////////////////////////////////////////////////////////////
@@ -109 +112 @@
 
     //
-    // sum up sigma/sqrt(area) for each sector, 
+    // sum up sigma^2/area for each sector, 
     // separately for inner and outer region;
     //
@@ -147 +150 @@
         const Double_t ratio = geom.GetPixRatio(idx);
 
+        //*fLog << "pixel id, ratio = " << idx << ",  " << ratio << endl;
+
         const MGeomPix &gpix = geom[idx];
 
@@ -171 +176 @@
         {
             innerPixels[sector]++;
-            innerSum[sector]+= sigma * sqrt(ratio);
+            innerSum[sector]+= sigma*sigma * ratio;
         }
         else
         {
             outerPixels[sector]++;
-            outerSum[sector]+= sigma * sqrt(ratio);
+            outerSum[sector]+= sigma*sigma * ratio;
         }
     }
@@ -192 +197 @@
     }
 
-    if (fInnerPixels > 0) fSigmabarInner = fSumInner / fInnerPixels;
-    if (fOuterPixels > 0) fSigmabarOuter = fSumOuter / fOuterPixels;
-
-    //
-    // this is the average sigma/sqrt(area) (over all pixels)
+    if (fInnerPixels > 0) fSigmabarInner = sqrt(fSumInner / fInnerPixels);
+    if (fOuterPixels > 0) fSigmabarOuter = sqrt(fSumOuter / fOuterPixels);
+
+    //
+    // this is the sqrt of the average sigma^2/area 
     //
     fSigmabar = (fInnerPixels+fOuterPixels)<=0 ? 0:
-                (fSumInner+fSumOuter)/(fInnerPixels+fOuterPixels);
+                sqrt( (fSumInner+fSumOuter)/(fInnerPixels+fOuterPixels) );
 
     for (UInt_t i=0; i<6; i++)
@@ -209 +214 @@
 
         const Double_t sum = ip + op;
-        fSigmabarInnerSector[i] = ip <=0 ? 0 :       iss/ip;
-        fSigmabarOuterSector[i] = op <=0 ? 0 :       oss/op;
-        fSigmabarSector[i]      = sum<=0 ? 0 : (iss+oss)/sum;
+        fSigmabarInnerSector[i] = ip <=0 ? 0 :        sqrt(iss/ip);
+        fSigmabarOuterSector[i] = op <=0 ? 0 :        sqrt(oss/op);
+        fSigmabarSector[i]      = sum<=0 ? 0 : sqrt((iss+oss)/sum);
     }
 
@@ -217 +222 @@
     //Print(opt);
 
-  return fSigmabar;
+  return fSigmabarInner;
 }
 

trunk/MagicSoft/Mars/manalysis/MSigmabar.h

@@ -14 +14 @@
 {
 private:
-    Float_t fSigmabar;          // Sigmabar ("average" RMS) of pedestal
+    Float_t fSigmabar;          // Sigmabar ( sqrt(average pedestalRMS^2) ) 
     Float_t fSigmabarInner;     // --only for inner pixels
     Float_t fSigmabarOuter;     // --only for outer pixels  

trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc

@@ -37 +37 @@
 //   MPedPhotCam
 //   MRawRunHeader
-//   MMcEvt  (FIXME: Must be replaced by a 'real-data' container)
+//   MPointingPos
 //   MCerPhotEvt
 //
@@ -58 +58 @@
 #include "MSigmabarParam.h"
 
-#include "MMcEvt.hxx"
+#include "MPointingPos.h"
 
 ClassImp(MSigmabarCalc);
@@ -105 +105 @@
 
     // This is needed for determining min/max Theta
-    fMcEvt = (MMcEvt*)pList->FindObject(AddSerialNumber("MMcEvt"));
-    if (!fMcEvt)
+    fPointPos = (MPointingPos*)pList->FindObject(AddSerialNumber("MPointingPos"));
+    if (!fPointPos)
     {
-        *fLog << warn << "MMcEvt not found... aborting." << endl;
+        *fLog << warn << "MPointingPos not found... aborting." << endl;
         fThetaMin =  0;
         fThetaMax = 90;
@@ -148 +148 @@
     if (rc<fSigmabarMin) fSigmabarMin=rc;
 
-    if (!fMcEvt)
+    if (!fPointPos)
         return kTRUE;
 
-    const Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;
+    const Double_t theta = fPointPos->GetZd();
 
-    if (theta>fSigmabarMax) fSigmabarMax=theta;
-    if (theta<fSigmabarMax) fSigmabarMin=theta;
+    if (theta>fThetaMax) fThetaMax=theta;
+    if (theta<fThetaMin) fThetaMin=theta;
 
     return kTRUE;
@@ -176 +176 @@
 void MSigmabarCalc::Reset()
 {
-    if (!fMcEvt)
+    if (!fPointPos)
         return;
 

trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.h

@@ -10 +10 @@
 #endif
 
-#ifndef MARS_MMcEvt
-#include "MMcEvt.hxx"
+#ifndef MARS_MPointingPos
+#include "MPointingPos.h"
 #endif
 
@@ -33 +33 @@
 {
 private:
-    MMcEvt         *fMcEvt;
+    MPointingPos   *fPointPos;
     MCerPhotEvt    *fEvt;
     MGeomCam       *fCam;