Changeset 4584 for trunk/MagicSoft/Mars

Timestamp:

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

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• macros/ONOFFAnalysis.C (modified) (51 diffs)
• manalysis/MPad.cc (modified) (87 diffs)
• manalysis/MPad.h (modified) (5 diffs)
• manalysis/MSigmabar.cc (modified) (7 diffs)
• manalysis/MSigmabar.h (modified) (1 diff)
• manalysis/MSigmabarCalc.cc (modified) (5 diffs)
• manalysis/MSigmabarCalc.h (modified) (2 diffs)
• mfilter/MFSelBasic.cc (modified) (3 diffs)
• mfilter/MFSelBasic.h (modified) (2 diffs)
• mfilter/Makefile (modified) (1 diff)
• mhist/MHSigmaTheta.cc (modified) (19 diffs)
• mhist/MHSigmaTheta.h (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -19 +19 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+
+
+  2004/08/12 : Wolfgang Wittek
+
+    * manalysis/MSigmabarCalc.[h,cc]
+      - replace MMcEvt by MPointingPos
+
+    * manalysis/MSigmabar.[h,cc]
+      - in member function Calc() return fSigmabarInner,
+        not fSigmabar
+      - update comments
+      - sigmabar is the sqrt of the average (pedRMS^2/area)
+
+    * manalysis/MPad.[h,cc]
+      - replace MMcEvt by MPointingPos
+      - remove bugs
+
+    * mfilter/MFSelBasic.[h,cc]
+      - replace MMcEvt by MPointingPos
+
+    * mfilter/Makefile
+      - add -I../mpointing
+
+
+    * mhist/MHSigmaTheta.[h,cc]
+      - replace MMcEvt by MPointingPos     
+      - replace 'MCerPhotPix cerpix' by 'MCerPhotPix &cerpix'
+      - add plot "Sigmabar(Outer) versus Theta"
+
+    * macros/ONOFFAnalysis.C
+      - Job A : got the padding working, work in progress
 
 

trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C

@@ -171 +171 @@
       //-----------------------------------------------
       //TString tag = "040126";
-      TString tag = "040127";
-      //TString tag = "040215";
+      //TString tag = "040127";
+      TString tag = "040215";
 
       //const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; 
@@ -183 +183 @@
       //const char *offfile  = "*.OFF";
       // 15 Feb 04
-      const char *offfile  = "*.OFF";
+      //const char *offfile  = "*.OFF";
+      const char *offfile  = "174*.OFF";
 
 
@@ -191 +192 @@
       //const char *onfile  = "MCerPhot_output";
       //const char *onfile  = "*.ON";
-      const char *onfile  = "1216*.ON";
+      //const char *onfile  = "1216*.ON";
       // 27 Jan 04
       //const char *onfile  = "12*.ON";
@@ -198 +199 @@
       // 15 Feb 04
       //const char *onfile  = "*.ON";
-
-
-     const char *mcfile  = "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root";
+      const char *onfile  = "174*.ON";
+
+
+      //const char *mcfile  = "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root";
       //const char *mcfile  = "/data/MAGIC/mc_eth/magLQE_4/gh/0/0/G_M0_00_0_550*.root";
       //const char *mcfile  = "/data/MAGIC/mc_eth/magLQE_5/gh/0/0/G_M0_00_0_550*.root";
+      //const char *mcfile = "calibrated_gammas";
+      const char *mcfile = "calibrated_data_david*";
+
       //-----------------------------------------------
 
@@ -224 +229 @@
 
      // 15 Feb 04, David
-     if (tag == "040215")
-       TString inPath = "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15/";
-   
+     //if (tag == "040215")
+     //  TString inPath = "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15/";
+     if (tag == "040215")   
+       TString inPathON = "/.magic/magicserv01/scratch/Daniel/CalibratedData/MispointingTest/2004_02_15ForgottenData/finesam20/";
+       TString inPathOFF = "/.magic/magicserv01/scratch/David/CalibratedData/Mkn421/2004_02_15/";
+       TString inPathMC = "/.magic/magicserv01/scratch/David/MCData/MCApril2004/Data/gammas_highnoise/";
 
       // path for output from Mars
       //TString outPath = "~wittek/datacrab_feb04/";
      //TString outPath = "~wittek/datacrab_01march04/";
-      TString outPath = "~wittek/datacrab_19mar04/";
+      TString outPath = "~wittek/datacrab_27april04/";
 
       //-----------------------------------------------
@@ -248 +256 @@
 
     Bool_t JobA    = kTRUE;  
-    Bool_t GPad    = kFALSE;    // generate padding histograms?
-    Bool_t WPad    = kFALSE;   // write out padding histograms ?
-    Bool_t RPad    = kFALSE;   // read in padding histograms ?
-    Bool_t Wout    = kTRUE;   // write out root file ON1.root 
-                               // (or OFF1.root or MC1.root)?
+    Bool_t GPadON  = kFALSE;    // \  generate padding histograms 
+    Bool_t GPadOFF = kFALSE;    //  | and write them onto disk
+    Bool_t GPadMC  = kFALSE;    // /
+    Bool_t Merge   = kFALSE;   // merge padding histograms
+                               // and write them onto disk
+    Bool_t Wout    = kTRUE;   // read in merged padding histograms and
+                              // write out root file of padded data
+                              // (ON1.root, OFF1.root or MC1.root) 
 
 
@@ -310 +321 @@
     // Job E_XX : extended version of E_XX (including flux plots)  
     //  - select g/h separation method XX
-    //  - read MC root file 
-    //  - calculate eff. collection area
+    //  - read MC root file
+    //  - calculate eff. collection area
     //  - optimize energy estimator
     //  - read ON root file 
@@ -346 +357 @@
     gLog << "Macro ONOFFAnalysis : Start of Job A" << endl;
     gLog << "" << endl;
-    gLog << "Macro ONOFFAnalysis : JobA, WPad, RPad, Wout = " 
-         << (JobA ? "kTRUE" : "kFALSE")  << ",  " 
-         << (WPad ? "kTRUE" : "kFALSE")  << ",  " 
-         << (RPad ? "kTRUE" : "kFALSE")  << ",  " 
-         << (Wout ? "kTRUE" : "kFALSE")  << endl;
+    gLog << "Macro ONOFFAnalysis : JobA, GPadON, GPadOFF, GPadMC, Merge, Wout = " 
+         << (JobA    ? "kTRUE" : "kFALSE")  << ",  " 
+         << (GPadON  ? "kTRUE" : "kFALSE")  << ",  " 
+         << (GPadOFF ? "kTRUE" : "kFALSE")  << ",  " 
+         << (GPadMC  ? "kTRUE" : "kFALSE")  << ",  " 
+         << (Merge   ? "kTRUE" : "kFALSE")  << ",  " 
+         << (Wout    ? "kTRUE" : "kFALSE")  << endl;
     
 
@@ -358 +371 @@
 
 
-    TString fileON  = inPath;
+    TString fileON  = inPathON;
     fileON += onfile;
     fileON += ".root";
 
-    TString fileOFF = inPath;
+    TString fileOFF = inPathOFF;
     fileOFF += offfile;
     fileOFF += ".root";
 
-    TString fileMC  = mcfile;
+    TString fileMC = inPathMC;
     fileMC += mcfile;
     fileMC += ".root";
@@ -373 +386 @@
          << fileOFF << ",  " << fileMC   << endl;
 
-    // name of file to conatin the histograms for the padding
+    // name of file to conatin the merged histograms for the padding
     TString outNameSigTh = outPath;
     outNameSigTh += "SigmaTheta";
@@ -379 +392 @@
 
     //--------------------------------------------------
+    // name of files to contain the paddding histograms of ON, OFF and MC data
+      TString NamePadON(outPath);
+      NamePadON += "PadON";
+      NamePadON += ".root";
+
+      TString NamePadOFF(outPath);
+      NamePadOFF += "PadOFF";
+      NamePadOFF += ".root";
+
+      TString NamePadMC(outPath);
+      NamePadMC += "PadMC";
+      NamePadMC += ".root";
+
+    //--------------------------------------------------
     // type of data to be padded 
-    TString typeInput = "ON";
-    //TString typeInput = "OFF";
-    //TString typeInput = "MC";
+      //TString typeInput = "ON";
+      //TString typeInput = "OFF";
+      TString typeInput = "MC";
     gLog << "typeInput = " << typeInput << endl;
 
@@ -416 +443 @@
     // 
     // read ON, OFF and MC data files
-    // generate (or read in) the padding histograms for ON and OFF data
-    //                       and merge these histograms
+    // generate (or read in) the padding histograms for ON, OFF and MC data
+    //
 
     MPad pad;
@@ -423 +450 @@
     pad.SetDataType(typeInput);
 
+    if (GPadON || GPadOFF || GPadMC)
+    {
     // generate the padding histograms
-    if (GPad)
-    {
       gLog << "=====================================================" << endl;
       gLog << "Start generating the padding histograms" << endl;
+    }
+
       //-----------------------------------------
       // ON events
+
+      if (GPadON)
+      {
 
       gLog << "-----------" << endl;
@@ -438 +470 @@
       MParList pliston;
 
+    MObservatory observon;
+
     MReadMarsFile  readON("Events", fileON);
     readON.DisableAutoScheme();
-    //MCT1ReadPreProc readON(fileON);
-
-      //MFSelBasic selthetaon;
-      //selthetaon.SetCuts(-100.0, 29.5, 35.5);
-      //MContinue contthetaon(&selthetaon);
-
-      MBlindPixelCalc blindon;
-      blindon.SetUseBlindPixels();
+
+    MGeomApply        apply;
+
+    MSourcePosfromStarPos sourcefromstaron;
+    sourcefromstaron.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0);
+
+      MBlindPixelsCalc2 blindon;
+      //blindon.SetUseBlindPixels();
+      blindon.SetCheckPedestalRms();
 
       MFSelBasic selbasicon;
@@ -459 +494 @@
 
       MHSigmaTheta sigthON("SigmaThetaON");
-      MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt");
+      MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MPointingPos");
       fillsigthetaON.SetName("FillSigTheta");    
  
@@ -467 +502 @@
       pliston.AddToList(&tliston);
       InitBinnings(&pliston);
+      pliston.AddToList(&observon);
       pliston.AddToList(&blindON);
       pliston.AddToList(&sigthON);
@@ -475 +511 @@
     
       tliston.AddToList(&readON);
-      //tliston.AddToList(&contthetaon);
+      tliston.AddToList(&apply);
+      tliston.AddToList(&sourcefromstaron);
 
       tliston.AddToList(&blindon);
@@ -489 +526 @@
       evtloopon.SetProgressBar(&baron);
 
-      Int_t maxeventson = -1;
-      //Int_t maxeventson = 10000;
+      //Int_t maxeventson = -1;
+      Int_t maxeventson = 10000;
       if ( !evtloopon.Eventloop(maxeventson) )
           return;
@@ -507 +544 @@
       TH2D *blindnthon  = blindON.GetBlindN();
 
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : write padding plots for ON data from file "
+           << NamePadON << endl;
+
+      TFile filewriteon(NamePadON, "RECREATE", "");
+      sigthon->Write();
+      sigpixthon->Write();
+      diffpixthon->Write();
+      blindidthon->Write();
+      blindnthon->Write();
+      filewriteon.Close();
+
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : padding plots for ON data written onto file "
+           << NamePadON << endl;
+      }
+      // end of GPadON
+
+
       //-----------------------------------------
       // OFF events
+
+      if (GPadOFF)
+      {
 
       gLog << "------------" << endl;
@@ -517 +576 @@
       MParList plistoff;
 
+    MObservatory observoff;
+
     MReadMarsFile  readOFF("Events", fileOFF);
     readOFF.DisableAutoScheme();
-    //      MCT1ReadPreProc readOFF(fileOFF);
-
-      MFSelBasic selthetaoff;
-      selthetaoff.SetCuts(-100.0, 29.5, 35.5);
-      MContinue contthetaoff(&selthetaoff);
-
-      MBlindPixelCalc blindoff;
-      blindoff.SetUseBlindPixels();
+
+    MSourcePosfromStarPos sourcefromstaroff;
+    sourcefromstaroff.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0);
+
+      MBlindPixelsCalc2 blindoff;
+      //blindoff.SetUseBlindPixels();
+      blindoff.SetCheckPedestalRms();
 
       MFSelBasic selbasicoff;
@@ -538 +598 @@
 
       MHSigmaTheta sigthOFF("SigmaThetaOFF");
-      MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt");
+      MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MPointingPos");
       fillsigthetaOFF.SetName("FillSigThetaOFF");     
 
@@ -546 +606 @@
       plistoff.AddToList(&tlistoff);
       InitBinnings(&plistoff);
+      plistoff.AddToList(&observoff);
       plistoff.AddToList(&blindOFF);
       plistoff.AddToList(&sigthOFF);
@@ -554 +615 @@
     
       tlistoff.AddToList(&readOFF);
-      //tlistoff.AddToList(&contthetaoff);
+      tlistoff.AddToList(&sourcefromstaroff);
 
       tlistoff.AddToList(&blindoff);
@@ -568 +629 @@
       evtloopoff.SetProgressBar(&baroff);
 
-      Int_t maxeventsoff = -1;
-      //Int_t maxeventsoff = 20000;
+      //Int_t maxeventsoff = -1;
+      Int_t maxeventsoff = 10000;
       if ( !evtloopoff.Eventloop(maxeventsoff) )
           return;
@@ -586 +647 @@
       TH2D *blindnthoff  = blindOFF.GetBlindN();
 
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : write padding plots for OFF data from file "
+           << NamePadOFF << endl;
+
+      TFile filewriteoff(NamePadOFF, "RECREATE", "");
+      sigthoff->Write();
+      sigpixthoff->Write();
+      diffpixthoff->Write();
+      blindidthoff->Write();
+      blindnthoff->Write();
+      filewriteoff.Close();
+
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : padding plots for OFF data written onto file "
+           << NamePadOFF << endl;
+      }
+      // end of GPadOFF
+
+
 
       //-----------------------------------------
       // MC events
+
+      if (GPadMC)
+      {
 
       gLog << "------------" << endl;
@@ -597 +680 @@
       MParList plistmc;
 
+    MObservatory observmc;
+
     MReadMarsFile  readMC("Events", fileMC);
     readMC.DisableAutoScheme();
-    //      MCT1ReadPreProc readMC(fileMC);
-
-      MFSelBasic selthetamc;
-      selthetamc.SetCuts(-100.0, 29.5, 35.5);
-      MContinue contthetamc(&selthetamc);
-
-      MBlindPixelCalc blindmc;
-      blindmc.SetUseBlindPixels();
+
+    MSourcePosfromStarPos sourcefromstarmc;
+
+      MBlindPixelsCalc2 blindmc;
+      //blindmc.SetUseBlindPixels();
+      blindmc.SetCheckPedestalRms();
 
       MFSelBasic selbasicmc;
@@ -618 +701 @@
 
       MHSigmaTheta sigthMC("SigmaThetaMC");
-      MFillH fillsigthetaMC ("SigmaThetaMC[MHSigmaTheta]", "MMcEvt");
+      MFillH fillsigthetaMC ("SigmaThetaMC[MHSigmaTheta]", "MPointingPos");
       fillsigthetaMC.SetName("FillSigThetaMC");     
 
@@ -626 +709 @@
       plistmc.AddToList(&tlistmc);
       InitBinnings(&plistmc);
+      plistmc.AddToList(&observmc);
       plistmc.AddToList(&blindMC);
       plistmc.AddToList(&sigthMC);
@@ -634 +718 @@
     
       tlistmc.AddToList(&readMC);
-      //tlistmc.AddToList(&contthetamc);
+      tlistmc.AddToList(&sourcefromstarmc);
 
       tlistmc.AddToList(&blindmc);
@@ -648 +732 @@
       evtloopmc.SetProgressBar(&barmc);
 
-      Int_t maxeventsmc = -1;
-      //Int_t maxeventsmc = 20000;
+      //Int_t maxeventsmc = -1;
+      Int_t maxeventsmc = 10000;
       if ( !evtloopmc.Eventloop(maxeventsmc) )
           return;
@@ -666 +750 @@
       TH2D *blindnthmc  = blindMC.GetBlindN();
 
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : write padding plots for MC data from file "
+           << NamePadMC << endl;
+
+      TFile filewritemc(NamePadMC, "RECREATE", "");
+      sigthmc->Write();
+      sigpixthmc->Write();
+      diffpixthmc->Write();
+      blindidthmc->Write();
+      blindnthmc->Write();
+      filewritemc.Close();
+
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : padding plots for MC data written onto file "
+           << NamePadMC << endl;
+      }
+      // end of GPadMC
 
       //-----------------------------------------
 
+    if (GPadON || GPadOFF || GPadMC)
+    {
+      gLog << "" << endl;
       gLog << "End of generating the padding histograms" << endl;
       gLog << "=====================================================" << endl;
-
-      pad.MergeONOFFMC(sigthmc,      diffpixthmc,  sigpixthmc, 
-                       blindidthmc,  blindnthmc,
-                       sigthon,      diffpixthon,  sigpixthon, 
-                       blindidthon,  blindnthon,
-                       sigthoff,     diffpixthoff, sigpixthoff, 
-                       blindidthoff, blindnthoff);
-
-
-      if (WPad)
-      {
-        // write the padding histograms onto a file  ---------
-        pad.WritePaddingDist(outNameSigTh);     
-      }
     }
 
-    // read the padding histograms ---------------------------
-    if (RPad)
+    if (Merge)
     {
-      pad.ReadPaddingDist(outNameSigTh);
+
+      //-----------------------------------
+      TH2D sigon;
+      TH3D sigpixon;
+      TH3D diffpixon;
+      TH2D blindidon;
+      TH2D blindnon;
+
+      TFile filereadon(NamePadON);
+      filereadon.ls();
+      sigon.Read("2D-ThetaSigmabar(Inner)");
+      sigpixon.Read("3D-ThetaPixSigma");
+      diffpixon.Read("3D-ThetaPixDiff");
+      blindidon.Read("2D-IdBlindPixels");
+      blindnon.Read("2D-NBlindPixels");
+
+      TH2D *sigthon     = new TH2D( (const TH2D&)sigon     );
+      TH3D *sigpixthon  = new TH3D( (const TH3D&)sigpixon  );
+      TH3D *diffpixthon = new TH3D( (const TH3D&)diffpixon );
+      TH2D *blindidthon = new TH2D( (const TH2D&)blindidon );
+      TH2D *blindnthon  = new TH2D( (const TH2D&)blindnon  );
+
+      /*
+      TH2D *sigthon     = (TH2D*)sigon.Clone();
+      TH3D *sigpixthon  = (TH3D*)sigpixon.Clone();
+      TH3D *diffpixthon = (TH3D*)diffpixon.Clone();
+      TH2D *blindidthon = (TH2D*)blindidon.Clone();
+      TH2D *blindnthon  = (TH2D*)blindnon.Clone();
+      */
+
+      //filereadon.Close();
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : padding plots for ON data read from file "
+           << NamePadON << endl;
+
+
+      //-----------------------------------
+      TH2D sigoff;
+      TH3D sigpixoff;
+      TH3D diffpixoff;
+      TH2D blindidoff;
+      TH2D blindnoff;
+
+      TFile filereadoff(NamePadOFF);
+      filereadoff.ls();
+      sigoff.Read("2D-ThetaSigmabar(Inner)");
+      sigpixoff.Read("3D-ThetaPixSigma");
+      diffpixoff.Read("3D-ThetaPixDiff");
+      blindidoff.Read("2D-IdBlindPixels");
+      blindnoff.Read("2D-NBlindPixels");
+ 
+      TH2D *sigthoff     = new TH2D( (const TH2D&)sigoff     );
+      TH3D *sigpixthoff  = new TH3D( (const TH3D&)sigpixoff  );
+      TH3D *diffpixthoff = new TH3D( (const TH3D&)diffpixoff );
+      TH2D *blindidthoff = new TH2D( (const TH2D&)blindidoff );
+      TH2D *blindnthoff  = new TH2D( (const TH2D&)blindnoff  );
+
+      /*
+      TH2D *sigthoff     = (TH2D*)sigoff.Clone();
+      TH3D *sigpixthoff  = (TH3D*)sigpixoff.Clone();
+      TH3D *diffpixthoff = (TH3D*)diffpixoff.Clone();
+      TH2D *blindidthoff = (TH2D*)blindidoff.Clone();
+      TH2D *blindnthoff  = (TH2D*)blindnoff.Clone();
+      */
+
+      //filereadoff.Close();
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : padding plots for OFF data read from file "
+           << NamePadOFF << endl;
+
+
+      //-----------------------------------
+      TH2D sigmc;
+      TH3D sigpixmc;
+      TH3D diffpixmc;
+      TH2D blindidmc;
+      TH2D blindnmc;
+
+      TFile filereadmc(NamePadMC);
+      filereadmc.ls();
+      sigmc.Read("2D-ThetaSigmabar(Inner)");
+      sigpixmc.Read("3D-ThetaPixSigma");
+      diffpixmc.Read("3D-ThetaPixDiff");
+      blindidmc.Read("2D-IdBlindPixels");
+      blindnmc.Read("2D-NBlindPixels");
+ 
+      TH2D *sigthmc     = new TH2D( (const TH2D&)sigmc     );
+      TH3D *sigpixthmc  = new TH3D( (const TH3D&)sigpixmc  );
+      TH3D *diffpixthmc = new TH3D( (const TH3D&)diffpixmc );
+      TH2D *blindidthmc = new TH2D( (const TH2D&)blindidmc );
+      TH2D *blindnthmc  = new TH2D( (const TH2D&)blindnmc  );
+
+      /*
+      TH2D *sigthmc     = (TH2D*)sigmc.Clone();
+      TH3D *sigpixthmc  = (TH3D*)sigpixmc.Clone();
+      TH3D *diffpixthmc = (TH3D*)diffpixmc.Clone();
+      TH2D *blindidthmc = (TH2D*)blindidmc.Clone();
+      TH2D *blindnthmc  = (TH2D*)blindnmc.Clone();
+      */
+
+      //filereadmc.Close();
+      gLog << "" << endl;
+      gLog << "Macro ONOFFAnalysis : padding plots for MC data read from file "
+           << NamePadMC << endl;
+
+      pad.MergeONOFFMC(*sigthmc,     *diffpixthmc, *sigpixthmc, 
+                       *blindidthmc, *blindnthmc,
+                       *sigthon,     *diffpixthon, *sigpixthon, 
+                       *blindidthon, *blindnthon,
+                       *sigthoff,    *diffpixthoff,*sigpixthoff, 
+                       *blindidthoff,*blindnthoff);
+
+      // write the target padding histograms onto a file  ---------
+      pad.WritePaddingDist(outNameSigTh);     
     }
+    // end of Merge
+
 
 
@@ -698 +902 @@
   if (Wout)
   {
+    // read the target padding histograms ---------------------------
+    pad.ReadPaddingDist(outNameSigTh);
+
+
     gLog << "=====================================================" << endl;
     gLog << "Start the padding" << endl;
@@ -723 +931 @@
     MGeomApply        apply;
 
-    //MPedestalWorkaround waround;
-
     // a way to find out whether one is dealing with MC :
-    MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5);  // MC
-    f1.SetName("Select MC");
-    MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5);  // data
-    f2.SetName("Select Data");
-
-    //if (typeInput ==  "ON")
-    //{
-    //   MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", 
-    //                                             "MCT1PointingCorrCalc");
-    //}
+    //MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5);  // MC
+    //f1.SetName("Select MC");
+    //MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5);  // data
+    //f2.SetName("Select Data");
+
     MSourcePosfromStarPos sourcefromstar;
     if (typeInput == "ON")
     {
-      sourcefromstar.SetSourceAndStarPosition(
-                               "Crab", 22,  0, 52, 5, 34, 32.0,
-                           "Zeta-Tau", 21,  8, 33, 5, 37, 38.7);
-      sourcefromstar.AddStar("Tau114", 21, 56, 13, 5, 27, 38.1);
-      //sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsOn.4.txt",0);
-
-      if(inPath == "/.magic/magicserv01/scratch/calibrated/")
-      {
-        sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions2stars.txt", 0);
-        //sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsNoStar.txt", 0);
-      }
-      else if(inPath == "/.magic/magicserv01/scratch/calibrated26/")
-        sourcefromstar.AddFile("~wittek/datacrab_26feb04/stars_2004_01_26", 0);
-
-      else if(inPath == "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15/")
-        sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0);
+      //sourcefromstar.SetSourceAndStarPosition(
+      //                         "Crab", 22,  0, 52, 5, 34, 32.0,
+      //                           "Zeta-Tau", 21,  8, 33, 5, 37, 38.7);
+      //sourcefromstar.AddStar("Tau114", 21, 56, 13, 5, 27, 38.1);
+      ////sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsOn.4.txt",0);
+
+      //if(inPath == "/.magic/magicserv01/scratch/calibrated/")
+      //{
+        //sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions2stars.txt", 0);
+        ////sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsNoStar.txt", 0);
+      //}
+      //else if(inPath == "/.magic/magicserv01/scratch/calibrated26/")
+      //  sourcefromstar.AddFile("~wittek/datacrab_26feb04/stars_2004_01_26", 0);
+
+      //else if(inPath == "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15/")
+      sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0);
     }
     else if (typeInput == "OFF")
     {
-      if(inPath == "/.magic/magicserv01/scratch/calibrated/")
-        sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsOff.txt", 0);
-      else if(inPath == "/.magic/magicserv01/scratch/calibrated26/")
-        sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions_2004_01_26", 0);
-      else if(inPath == "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15")
-        sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0);
+      //if(inPath == "/.magic/magicserv01/scratch/calibrated/")
+      //  sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsOff.txt", 0);
+      //else if(inPath == "/.magic/magicserv01/scratch/calibrated26/")
+      //  sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions_2004_01_26", 0);
+      //else if(inPath == "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15")
+      sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0);
     }
-
-
-    //MBlindPixelCalc blindbeforepad;
+    else if (typeInput == "MC")
+    {
+    }
+
+    MBlindPixelsCalc2 blindbeforepad;
     //blindbeforepad.SetUseBlindPixels();
-    //blindbeforepad.SetName("BlindBeforePadding");
-
-    //MBlindPixelCalc blind;
-    //blind.SetUseBlindPixels();
-    //blind.SetUseInterpolation();
-    //blind.SetName("BlindAfterPadding");
-
-    MSigmabarCalc sigbar;
+    blindbeforepad.SetCheckPedestalRms();
+    blindbeforepad.SetName("BlindBeforePadding");
 
     //MBadPixelCalcRms blind;
@@ -794 +991 @@
     MSigmabarCalc sigbarcalc;
 
-    MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
+    MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MPointingPos");
     fillsigtheta.SetName("HSigmaTheta");
 
@@ -845 +1042 @@
       //write.AddContainer("MMcRunHeader",  "RunHeaders", kFALSE);
       //write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      //write.AddContainer("ThetaOrig",     "Events");
+      write.AddContainer("MPointingPos", "Events");
       write.AddContainer("MSrcPosCam",    "Events");
       write.AddContainer("MSigmabar",     "Events");
@@ -868 +1064 @@
     
     tliston.AddToList(&read);
-
     //tliston.AddToList(&f1);
     //tliston.AddToList(&f2);
     tliston.AddToList(&apply);
-    //tliston.AddToList(&waround);
     tliston.AddToList(&sourcefromstar);
 
-    //tliston.AddToList(&blindbeforepad);
-    //  tliston.AddToList(&pad); 
-    //  if (typeInput ==  "ON")
-    //  tliston.AddToList(&pointcorr);
-
-    tliston.AddToList(&sigbar);
+    tliston.AddToList(&blindbeforepad);
+
+    tliston.AddToList(&contbasic);
+    tliston.AddToList(&pad); 
+
     tliston.AddToList(&blind);
-    tliston.AddToList(&contbasic);
-
     tliston.AddToList(&fillblind);
-    //tliston.AddToList(&sigbarcalc);
+    tliston.AddToList(&sigbarcalc);
     tliston.AddToList(&fillsigtheta);
     tliston.AddToList(&clean);
@@ -913 +1104 @@
     //  evtloop.Write();
 
-    Int_t maxevents = -1;
-    //Int_t maxevents = 1000;
+    //Int_t maxevents = -1;
+    Int_t maxevents = 2000;
     if ( !evtloop.Eventloop(maxevents) )
         return;
@@ -1109 +1300 @@
     matrixg.EnableGraphicalOutput();
 
-    matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)");
+    matrixg.AddColumn("cos((MPointingPos.fZd)/kRad2Deg)");
     matrixg.AddColumn("MSigmabar.fSigmabar");
     matrixg.AddColumn("log10(MHillas.fSize)");
@@ -1198 +1389 @@
     bing.SetEdges(10, 0., 1.0);
 
-    MH3 gref("cos(MMcEvt.fTelescopeTheta)");
+    MH3 gref("cos((MPointingPos.fZd)/kRad2Deg)");
     gref.SetName(mgname);
     MH::SetBinning(&gref.GetHist(), &bing);
@@ -1231 +1422 @@
       writetraing.AddContainer("MRawRunHeader", "RunHeaders");
       writetraing.AddContainer("MTime",         "Events");
-      writetraing.AddContainer("MMcEvt",        "Events");
+      writetraing.AddContainer("MPointingPos",  "Events");
       writetraing.AddContainer("ThetaOrig",     "Events");
       writetraing.AddContainer("MSrcPosCam",    "Events");
@@ -1248 +1439 @@
       writetestg.AddContainer("MRawRunHeader", "RunHeaders");
       writetestg.AddContainer("MTime",         "Events");
-      writetestg.AddContainer("MMcEvt",        "Events");
+      writetestg.AddContainer("MPointingPos",  "Events");
       writetestg.AddContainer("ThetaOrig",     "Events");
       writetestg.AddContainer("MSrcPosCam",    "Events");
@@ -1331 +1522 @@
     binh.SetEdges(10, 0., 1.0);
 
-    //MH3 href("cos(MMcEvt.fTelescopeTheta)");
+    //MH3 href("cos((MPointingPos.fZd)/kRad2Deg)");
     //href.SetName(mhname);
     //MH::SetBinning(&href.GetHist(), &binh);
@@ -1366 +1557 @@
       writetrainh.AddContainer("MRawRunHeader", "RunHeaders");
       writetrainh.AddContainer("MTime",         "Events");
-      writetrainh.AddContainer("MMcEvt",        "Events");
+      writetrainh.AddContainer("MPointingPos",  "Events");
       writetrainh.AddContainer("ThetaOrig",     "Events");
       writetrainh.AddContainer("MSrcPosCam",    "Events");
@@ -1382 +1573 @@
       writetesth.AddContainer("MRawRunHeader", "RunHeaders");
       writetesth.AddContainer("MTime",         "Events");
-      writetesth.AddContainer("MMcEvt",        "Events");
+      writetesth.AddContainer("MPointingPos",  "Events");
       writetesth.AddContainer("ThetaOrig",     "Events");
       writetesth.AddContainer("MSrcPosCam",    "Events");
@@ -1667 +1858 @@
       write.AddContainer("MRawRunHeader", "RunHeaders");
       write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("MPointingPos",  "Events");
       write.AddContainer("ThetaOrig",     "Events");
       write.AddContainer("MSrcPosCam",    "Events");
@@ -1992 +2183 @@
       bin.SetEdges(10, 0., 1.0);
 
-      MH3 mh3("cos(MMcEvt.fTelescopeTheta)");
+      MH3 mh3("cos((MPointingPos.fZd)/kRad2Deg)");
       mh3.SetName(mname);
       MH::SetBinning(&mh3.GetHist(), &bin);
@@ -2355 +2546 @@
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
-      //write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      //write.AddContainer("ThetaOrig",     "Events");
+      write.AddContainer("MTime",         "Events");
+      write.AddContainer("MPointingPos",  "Events");
       write.AddContainer("MSrcPosCam",    "Events");
       write.AddContainer("MSigmabar",     "Events");
@@ -3137 +3327 @@
 
 
-    MFillH filler("MHMcCollectionArea", "MMcEvt");
+    MFillH filler("MHMcCollectionArea", "MPointingPos");
     filler.SetName("CollectionArea");
 
@@ -3326 +3516 @@
       write.AddContainer("MRawRunHeader", "RunHeaders");
       write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("MPointingPos",  "Events");
       write.AddContainer("ThetaOrig",     "Events");
       write.AddContainer("MSrcPosCam",    "Events");
@@ -3400 +3590 @@
     // new from Robert
 
-    MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
+    MFillH hfill6("MHTimeDiffTheta", "MPointingPos");
     hfill6.SetName("HTimeDiffTheta");
 
-    MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
+    MFillH hfill6a("MHTimeDiffTime", "MPointingPos");
     hfill6a.SetName("HTimeDiffTime");
 
@@ -3611 +3801 @@
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+   
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

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;

trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc

@@ -46 +46 @@
 #include "MParList.h"
 
-#include "MMcEvt.hxx"
+#include "MPointingPos.h"
 
 #include "MCerPhotEvt.h"
@@ -104 +104 @@
     }
 
-    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
-    if (!fMcEvt)
-    {
-        *fLog << dbginf << "MMcEvt not found... aborting." << endl;
+    fPointPos = (MPointingPos*)pList->FindObject("MPointingPos");
+    if (!fPointPos)
+    {
+        *fLog << dbginf << "MPointingPos not found... aborting." << endl;
         return kFALSE;
     }
@@ -146 +146 @@
 Int_t MFSelBasic::Process()
 {
-    const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
+    const Double_t theta = fPointPos->GetZd();
 
     fResult  = kFALSE;

trunk/MagicSoft/Mars/mfilter/MFSelBasic.h

@@ -14 +14 @@
 #endif
 
-class MMcEvt;
+class MPointingPos;
 class MGeomCam;
 class MCerPhotEvt;
@@ -23 +23 @@
 {
 private:
-    const MMcEvt        *fMcEvt;       
+    const MPointingPos  *fPointPos;       
     const MGeomCam      *fCam;      // Camera Geometry 
     const MCerPhotEvt   *fEvt;      // Cerenkov Photon Event 

trunk/MagicSoft/Mars/mfilter/Makefile

@@ -13 +13 @@
 INCLUDES = -I. -I../mbase -I../mfbase -I../mraw -I../mmc -I../mdata \
            -I../manalysis -I../mfileio  -I../mgeom -I../mimage      \
-           -I../mhbase -I../mmain -I../mgui -I../msignal            \
+           -I../mhbase -I../mmain -I../mgui -I../msignal -I../mpointing  \
            -I../mpedestal
 

trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc

@@ -26 +26 @@
 //////////////////////////////////////////////////////////////////////////////
 //                                                                          //
-//  MHSigmaTheta (extension of Robert's MHSigmabarTheta)                    //
+//  MHSigmaTheta                                                            //
 //                                                                          //
-//  calculates - the 2D-histogram   sigmabar vs. Theta, and                 //
+//  calculates - the 2D-histogram   sigmabar vs. Theta (Inner), and         //
+//             - the 2D-histogram   sigmabar vs. Theta (Outer)              //
 //             - the 3D-histogram   sigma, pixel no., Theta                 //
 //             - the 3D-histogram   (sigma^2-sigmabar^2), pixel no., Theta  //
@@ -39 +40 @@
 
 #include "MTime.h"
-#include "MMcEvt.hxx"
+#include "MPointingPos.h"
 
 #include "MBinning.h"
@@ -46 +47 @@
 
 #include "MGeomCam.h"
+#include "MBlindPixels.h"
 
 #include "MPedPhotCam.h"
@@ -70 +72 @@
 
     fSigmaTheta.SetDirectory(NULL);
-    fSigmaTheta.SetName("2D-ThetaSigmabar");
+    fSigmaTheta.SetName("2D-ThetaSigmabar(Inner)");
     fSigmaTheta.SetTitle("2D: \\bar{\\sigma}, \\Theta");
     fSigmaTheta.SetXTitle("\\Theta [\\circ]");
-    fSigmaTheta.SetYTitle("Sigmabar");
+    fSigmaTheta.SetYTitle("Sigmabar(Inner) / SQRT(Area)");
+
+    fSigmaThetaOuter.SetDirectory(NULL);
+    fSigmaThetaOuter.SetName("2D-ThetaSigmabar(Outer)");
+    fSigmaThetaOuter.SetTitle("2D: \\bar{\\sigma}, \\Theta");
+    fSigmaThetaOuter.SetXTitle("\\Theta [\\circ]");
+    fSigmaThetaOuter.SetYTitle("Sigmabar(Outer) / SQRT(Area)");
 
     fSigmaPixTheta.SetDirectory(NULL);
@@ -80 +88 @@
     fSigmaPixTheta.SetXTitle("\\Theta [\\circ]");
     fSigmaPixTheta.SetYTitle("Pixel Id");
-    fSigmaPixTheta.SetZTitle("\\sigma");
+    fSigmaPixTheta.SetZTitle("Sigma");
 
     fDiffPixTheta.SetDirectory(NULL);
@@ -87 +95 @@
     fDiffPixTheta.SetXTitle("\\Theta [\\circ]");
     fDiffPixTheta.SetYTitle("Pixel Id");
-    fDiffPixTheta.SetZTitle("{\\sigma}^{2}-\\bar{\\sigma}^{2}");
+    fDiffPixTheta.SetZTitle("(Sigma2 - Sigmabar2)/Area");
 
     // Set default binning
@@ -102 +110 @@
     binspix.SetEdges(578, -0.5, 577.5);
 
-    SetBinning(&fSigmaTheta,    &binst, &binsb);
-    SetBinning(&fSigmaPixTheta, &binst, &binspix, &binsb);
-    SetBinning(&fDiffPixTheta,  &binst, &binspix, &binsd);
+    SetBinning(&fSigmaTheta,      &binst, &binsb);
+    SetBinning(&fSigmaThetaOuter, &binst, &binsb);
+    SetBinning(&fSigmaPixTheta,   &binst, &binspix, &binsb);
+    SetBinning(&fDiffPixTheta,    &binst, &binspix, &binsd);
 }
 
@@ -120 +129 @@
     }
 
-    fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
-    if (!fMcEvt)
-        *fLog << warn << "MMcEvt not found... aborting." << endl;
+    fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
+    if (!fPointPos)
+        *fLog << warn << "MPointingPos not found... aborting." << endl;
 
 
@@ -132 +141 @@
     }
     fPed->InitSize(fCam->GetNumPixels());
+
+    
+    fBlindPix = (MBlindPixels*)plist->FindObject("MBlindPixels");
+    if (!fBlindPix)
+    {  
+       *fLog << err << "MBlindPixels not found... continue. " << endl; 
+    }
 
 
@@ -177 +193 @@
     if (binssigma)
     {
-        SetBinning(&fSigmaTheta, binstheta, binssigma);
-        SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma);
+        SetBinning(&fSigmaTheta,      binstheta, binssigma);
+        SetBinning(&fSigmaThetaOuter, binstheta, binssigma);
+        SetBinning(&fSigmaPixTheta,   binstheta, &binspix, binssigma);
     }
 
@@ -191 +208 @@
 //  Fill the histograms
 //
+//  ignore pixels if they are unused or blind
+//
 Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w)
 {
-    Double_t theta = fMcEvt ? fMcEvt->GetTelescopeTheta()*kRad2Deg : 0;
+    Double_t theta = fPointPos->GetZd();
     fSigmabar->Calc(*fCam, *fPed, *fEvt);
-    Double_t mysig = fSigmabar->GetSigmabarInner();
-
-    //*fLog << "theta, mysig = " << theta << ",  " << mysig << endl;
+    Double_t mysig      = fSigmabar->GetSigmabarInner();
+    Double_t mysigouter = fSigmabar->GetSigmabarOuter();
+
+    //*fLog << "theta, mysig, mysigouter = " << theta << ",  " << mysig 
+    //      << ",  " << mysigouter << endl;
 
     fSigmaTheta.Fill(theta, mysig);
+    fSigmaThetaOuter.Fill(theta, mysigouter);
 
     const UInt_t npix = fEvt->GetNumPixels();
@@ -205 +227 @@
     for (UInt_t i=0; i<npix; i++)
     {
-        MCerPhotPix cerpix = (*fEvt)[i];
+        MCerPhotPix &cerpix = (*fEvt)[i];
+        const Int_t id = cerpix.GetPixId();
+
         if (!cerpix.IsPixelUsed())
-            continue;
-
-        const Int_t id = cerpix.GetPixId();
+        {
+          //*fLog << all << "MHSigmaTheta::Fill; unused pixel found, id = "
+          //      << id << endl;
+          continue;
+        }
+
         const MPedPhotPix &pix = (*fPed)[id];
+
+        if ( fBlindPix != NULL  &&  fBlindPix->IsBlind(id) )
+        {
+          // this should never occur, because blind pixels should have
+          // been set unused by MBlindPixelsCalc2::UnMap()
+          //*fLog << all << "MHSigmaTheta::Fill; blind pixel found which is used, id = "
+          //      << id << "... go to next pixel." << endl;
+          continue;
+        }
 
         // ratio is the area of pixel 0 
@@ -219 +255 @@
         fSigmaPixTheta.Fill(theta, (Double_t)id, sigma*sqrt(ratio));
 
-        const Double_t diff = sigma*sigma*ratio - mysig*mysig;
+        Double_t diff;
+        if (ratio > 0.5)
+        {
+          // inner pixel
+          diff = sigma*sigma*ratio - mysig*mysig;
+        }
+        else
+        {
+          // outer pixel
+          diff = sigma*sigma*ratio - mysigouter*mysigouter;
+        }
+
         fDiffPixTheta.Fill(theta, (Double_t)id, diff);
     }
@@ -262 +309 @@
     AppendPad("");
 
-    pad->Divide(3, 2);
+    pad->Divide(3, 3);
 
     // draw the 2D histogram Sigmabar versus Theta
@@ -283 +330 @@
     h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)");
     h->SetXTitle("\\Theta [\\circ]");
-    h->SetYTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2}");
+    h->SetYTitle("(Sigma2 - Sigmabar2) / Area");
     h->Draw("box");
     h->SetBit(kCanDelete);
@@ -293 +340 @@
     h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)");
     h->SetXTitle("\\Theta [\\circ]");
-    h->SetYTitle("\\sigma_{ped}");
+    h->SetYTitle("Sigma");
     h->Draw("box");
     h->SetBit(kCanDelete);
@@ -313 +360 @@
     h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all  \\Theta)");
     h->SetXTitle("Pixel Id");
-    h->SetYTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2}");
+    h->SetYTitle("Sigma2 - sigmabar2) / Area");
     h->Draw("box");
     h->SetBit(kCanDelete);
@@ -323 +370 @@
     h->SetTitle("\\sigma_{ped} vs. pixel Id (all  \\Theta)");
     h->SetXTitle("Pixel Id");
-    h->SetYTitle("\\sigma_{ped}");
+    h->SetYTitle("Sigma");
     h->Draw("box");
     h->SetBit(kCanDelete);
@@ -329 +376 @@
     pad->cd(4);
     fSigmaTheta.Draw(opt);
+
+    pad->cd(7);
+    fSigmaThetaOuter.Draw(opt);
 
     //pad->cd(8);

trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h

@@ -16 +16 @@
 class MGeomCam;
 class MCerPhotEvt;
-class MMcEvt;
+class MPointingPos;
 class MPedPhotCam;
 class MSigmabar;
 class MParList;
+class MBlindPixels;
 
 
@@ -29 +30 @@
     MCerPhotEvt    *fEvt;        //!
     MSigmabar      *fSigmabar;   //!
-    MMcEvt         *fMcEvt;      //!
- 
-    TH2D fSigmaTheta;    // 2D-distribution sigmabar versus Theta; 
-                         // sigmabar is the average pedestasl sigma in an event
+    MPointingPos   *fPointPos;   //!
+    MBlindPixels   *fBlindPix;   //!
+
+                           // sigmabar is the average pedestal sigma  
+    TH2D fSigmaTheta;      // 2D-distribution sigmabar versus Theta (Inner) 
+    TH2D fSigmaThetaOuter; // 2D-distribution sigmabar versus Theta (Outer) 
+
     TH3D fSigmaPixTheta; // 3D-distr.:Theta, pixel, pedestal sigma
     TH3D fDiffPixTheta;  // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2
@@ -46 +50 @@
     const TH2D *GetSigmaTheta() { return &fSigmaTheta; }
     const TH2D *GetSigmaTheta() const { return &fSigmaTheta; }
+
+    const TH2D *GetSigmaThetaOuter() { return &fSigmaThetaOuter; }
+    const TH2D *GetSigmaThetaOuter() const { return &fSigmaThetaOuter; }
 
     const TH3D *GetSigmaPixTheta() { return &fSigmaPixTheta; }