Changeset 1809

Timestamp:

03/08/03 14:00:30 (22 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft

Files:

• Mars/Changelog (modified) (2 diffs)
• Mars/macros/AnalyseCT1.C (modified) (33 diffs)
• Mars/manalysis/MMultiDimDistCalc.cc (modified) (1 diff)
• Mars/manalysis/MPadSchweizer.cc (modified) (17 diffs)
• Mars/manalysis/MPadSchweizer.h (modified) (1 diff)
• Mars/manalysis/MSelBasic.cc (modified) (5 diffs)
• Mars/manalysis/MSelBasic.h (modified) (1 diff)
• Mars/manalysis/MSelFinal.cc (modified) (7 diffs)
• Mars/manalysis/MSelFinal.h (modified) (2 diffs)
• Mars/manalysis/MSelStandard.cc (modified) (11 diffs)
• Mars/manalysis/MSelStandard.h (modified) (1 diff)
• Mars/mfileio/MCT1ReadPreProc.cc (modified) (16 diffs)
• Mars/mfileio/MCT1ReadPreProc.h (modified) (1 diff)
• Mars/mfilter/MFEventSelector.cc (modified) (5 diffs)
• Mars/mfilter/MFEventSelector.h (modified) (1 diff)
• Mars/mhist/MHHadronness.cc (modified) (5 diffs)
• Mars/mhist/MHMatrix.cc (modified) (8 diffs)
• Mars/mhist/MHSigmaTheta.cc (modified) (5 diffs)
• include-Classes/MMcFormat/MMcEvt.hxx (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+
+
+ 2003/03/08: Wolfgang Wittek
+
+    * macros/AnalyseCT1.C
+      - for the CT1 data analysis
+
+    * mhist/MHMatrix.[h,cc]
+      - let refcolumn start at 1 (not at 0)
+
+    * mhist/MHSigmaTheta.[h,cc]
+      - Draw replaced by DrawCopy
+      - add SetDirectory(NULL)
+
+    * manalysis/ MSelBasic.[h,cc]
+                 MSelStandard.[h,cc]
+                 MSelFinal.[h,cc]
+      - more detailed output for errors
+      - bugs removed
+      
+    * manalysis/ MPadSchweizer.[h,cc]
+      - add SetDirectory(NULL)
+      - add fErrors
+
+    * mfilter/ MFEventSelector.[h,cc]
+      - add fErrors
+
+    * manalysis/MMultiDimDistCalc.[h,cc]
+      - check division by zero
+
+    * mhist/MHHadronness.[h,cc]
+      - check division by zero
+      - normalize distributions of hadronness
+
+    * mfileio/MCT1ReadPreProc.[h,cc]
+      - add event number (event.isecs_since_midday)
+      - change definition of "fIsMcFile", 
+        because outpars.bmontecarlo is set wrongly sometimes
+      - copy pedestalRMS for each event from the header information
+      - check for the presence of a footer record even after reading 
+        a run header
+
+    * mmc/ MMcEvt.[hxx,cxx]]
+      - add GetEvtNumber()
+
+
  2003/02/27: Abelardo Moralejo
 
@@ -63 +109 @@
     * mmc/MMcTrigHeader.hxx
       - Added GetMultiplicity() and GetMeanThreshold()
-
 
 

trunk/MagicSoft/Mars/macros/AnalyseCT1.C

@@ -1 +1 @@
 void AnalyseCT1()
 {
-    const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; 
-    const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; 
+  //-----------------------------------------------
+    //const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; 
+    //const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; 
     //const char *mcfile = "/hd43/rwagner/CT1/PreprocData/mc_Gammas.preproc.0.6";
     //const char *mcfile = "/hd43/rwagner/mkn421_mc.preproc";
-    const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc";
+    //const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc";
+
+  //-----------------------------------------------
+  const char *offfile = "~/Mars200203/CT1data/offdata.preproc"; 
+
+  const char *onfile  = "~/Mars200203/CT1data/mkn421_on.preproc"; 
+  //const char *onfile  = "~/Mars200203/CT1data/Crab_preproc_daniel_0.6"; 
+
+  //const char *mcfile  = "~/Mars200203/CT1data/mc_gammas.preproc"; //Version 0.5
+  //const char *mcfile  = "~/Mars200203/CT1data/mc_Gammas.preproc.0.6";
+  const char *mcfile = "~/Mars200203/CT1data/mkn421_mc_pedrms_0.001.preproc";
+
+  //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6";
+  //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6_040303";
+
+  //-----------------------------------------------
+  //const char *offfile = "~magican/home/ct1test/PreprocData/offdata.preproc"; 
+  //const char *onfile  = "~magican/home/ct1test/PreprocData/mkn421_on.preproc"; 
+  //const char *mcfile  = "~magican/home/ct1test/PreprocData/mc_gammas.preproc";
+
+
     const char *filename;
 
-    Bool_t Data  = kFALSE;   // loop over data ?
-    Bool_t WLook = kFALSE;   // write out look-alike events for hadrons ?
-    Bool_t WHist = kFALSE;   // write out histograms for padding ?
-
-    Bool_t MC      = kTRUE;  // loop over MC gamma data ?
-    Bool_t RHist   = kTRUE;  // read in histograms for padding ?
-    Bool_t WLookMC = kTRUE;  // write out look-alike events for gammas ?
-
-
-    cout << "" << endl;
-    cout << "Macro AnalyseCT1 : Data, WHist = " << Data << ",  "
-         << WHist << endl;
-
-    cout << "                         RHist, MC   = " << RHist << ",  "
-         << MC << endl;
-    cout << "" << endl;
+    //************************************************************************
+    // Job 1 : read OFF data  
+    // (generate sigmabar vs. Theta plot; write root file for OFF data (OFF1);
+    //  generate hadron matrix for g/h separation)
+
+    Bool_t Job1     = kFALSE;   
+    Bool_t WHistOFF = kTRUE;   // write out histogram sigmabar vs. Theta ?
+    Bool_t WLookOFF = kTRUE;   // write out look-alike events for hadrons ?
+    Bool_t WOFF1    = kTRUE;   // write out root file OFF1 ?
+
+
+    // Job 2 : read ON data
+    // (generate sigmabar vs. Theta plot; write root file for ON data (ON1))
+
+    Bool_t Job2    = kFALSE;  
+    Bool_t WHistON = kFALSE;    // write out histogram sigmabar vs. Theta ?
+    Bool_t WLookON = kFALSE;   // write out look-alike events for hadrons ?
+    Bool_t WON1    = kFALSE;   // write out root file ON1 ?
+
+    // Job 3 : read MC gamma data, read sigmabar vs. Theta plot from OFF data  
+    // (do padding; write root file for MC gammas (MC1);
+    //  generate gamma matrix for g/h separation)
+
+    Bool_t Job3     = kFALSE;  
+    Bool_t WLookMC  = kTRUE;  // write out look-alike events for gammas ?
+    Bool_t WMC1     = kTRUE;  // write out root file MC1 ?
+
+
+
+    // Job 4 : read OFF1 and MC1 files, read hadron and gamma matrix
+    // (produce the Neyman-Pearson plot)
+    Bool_t Job4  = kTRUE;  
+
+
+    // Job 5 : read MC1 file, read hadron and gamma matrix 
+    // (do g/h separation; write root file for MC gammas after final cuts (MC2))
+    Bool_t Job5   = kFALSE;  
+    Bool_t WMC2   = kFALSE;  // write out root file MC2 ?
+
+
+    // Job 6 : read ON data, read hadron and gamma matrix
+    // (do g/h separation; write root file for ON data after final cuts (ON2))
+    Bool_t Job6  = kFALSE;  
+    Bool_t WON2  = kTRUE;  // write out root file ON2 ?
+    //************************************************************************
+
+
+
 
   //---------------------------------------------------------------------
-  // start of section for ON data
-  // (in the CT1 analysis the ON data are also used as a sample representing
-  //  the hadrons for the g/h separation)  
-
-    // read ON data file 
-
-    //  - to produce the 2D-histogram "sigmabar versus Theta"
+    // Job 1
+    //======
+
+    // read OFF data file 
+
+    //  - produce the 2D-histogram "sigmabar versus Theta" for OFF data
     //    (to be used for the padding of the MC gamma data)
 
-    //  - to write a file of look-alike events (for hadrons)
+    //  - write a file of look-alike events (hadron matrix)
     //    (to be used later together with the corresponding MC gamma file
     //     for the g/h separation in the analysis of the ON data)
 
-    //  - to write a file of ON events (after the standard cuts, 
-    //                                  before the g/h separation)
+    //  - write a file of OFF events (OFF1) 
+    //    (after the standard cuts, before the g/h separation)
     //    (to be used together with the corresponding MC gamma file
     //     for the optimization of the g/h separation)
 
- if (Data)
+ if (Job1)
  {
     cout << "=====================================================" << endl;
-    cout << "Macro AnalyseCT1 : Start of section for ON data" << endl;
+    cout << "Macro AnalyseCT1 : Start of Job 1" << endl;
+
+    cout << "" << endl;
+    cout << "Macro AnalyseCT1 : Job1, WHistOFF, WLookOFF, WOFF1 = " 
+         << Job1  << ",  " << WHistOFF << ",  " << WLookOFF << ",  " 
+         << WOFF1 << endl;
+
+
+    TString typeData = "OFF";
+    cout << "typeData = " << typeData << endl;
 
     MTaskList tliston;
@@ -112 +173 @@
     //
 
-    filename = onfile;
-    readname = "ReadCT1data";
-    MCT1ReadPreProc read(filename, readname);
+    filename = offfile;
+    RName = "ReadCT1data";
+    MCT1ReadPreProc read(filename, RName);
 
     MSelBasic selbasic;
@@ -150 +211 @@
     MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
 
-    MSelStandard selstand(fHillas, fHillasSrc);
+    MSelStandard selstand(fHilName, fHilSrcName);
 
     MFillH hfill1("MHHillas",    fHilName);
@@ -168 +229 @@
 
     const char* mtxName = "MatrixHadrons";
-    Int_t howMany = 30000;
-    char* outName = "matrix_hadrons.root";
+    Int_t howMany = 50000;
+    TString outName = "matrix_hadrons_";
+    outName += typeData;
+    outName += ".root";
 
     cout << "" << endl;
@@ -182 +245 @@
 
     // matrix limitation for look-alike events (approximate number)
-    MFEventSelector selector("", "", readname);
+    MFEventSelector selector("", "", RName);
     selector.SetNumSelectEvts(howMany);
     MFilterList flist;
@@ -192 +255 @@
     MHMatrix *pmatrix = new MHMatrix(mtxName);
     MHMatrix& matrix = *pmatrix;
+
+    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
+    matrix.AddColumn("MSigmabar.fSigmabar");
+    matrix.AddColumn("log10(Hillas.fSize)");
+    matrix.AddColumn("HillasSrc.fDist");
     matrix.AddColumn("Hillas.fWidth");
     matrix.AddColumn("Hillas.fLength");
-    matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize");
-    matrix.AddColumn("abs(Hillas.fAsym)");
+    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
     matrix.AddColumn("abs(Hillas.fM3Long)");
-    matrix.AddColumn("abs(Hillas.fM3Trans)");
-    matrix.AddColumn("abs(HillasSrc.fHeadTail)");
     matrix.AddColumn("Hillas.fConc");
     matrix.AddColumn("Hillas.fConc1");
-    matrix.AddColumn("HillasSrc.fDist");
-    matrix.AddColumn("log10(Hillas.fSize)");
-    matrix.AddColumn("HillasSrc.fAlpha");
-    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
-    //matrix.AddColumn("MMcEvt.fTheta");
+
     pliston.AddToList(pmatrix);
 
@@ -211 +272 @@
     fillmatg.SetFilter(&flist);
 
+
+    if (WOFF1)
+    {    
+      TString outNameImage = typeData;
+      outNameImage += "1.root";
+      MWriteRootFile write(outNameImage);
+      write.AddContainer("MSigmabar",     "Events");
+      write.AddContainer("Hillas",        "Events");
+      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("HillasSrc",     "Events");
+      write.AddContainer("MRawRunHeader", "RunHeaders");
+      //write.AddContainer("MMcRunHeader","RunHeaders");
+      write.AddContainer("MSrcPosCam",    "RunHeaders");
+    }
+
     //+++++++++++++++++++++++++++++++++++++++++++++++++++
 
-    MSelFinal selfinal(fHillas, fHillasSrc);
 
     //*****************************
@@ -229 +304 @@
 
     tliston.AddToList(&selstand);
+    if (WOFF1)
+      tliston.AddToList(&write);
+
+    tliston.AddToList(&flist);
+    tliston.AddToList(&fillmatg);    
+
     tliston.AddToList(&hfill1);
     tliston.AddToList(&hfill2);
@@ -234 +315 @@
     tliston.AddToList(&hfill4);
 
-    tliston.AddToList(&flist);
-    tliston.AddToList(&fillmatg);
-
-    //tliston.AddToList(&selfinal);
     //*****************************
 
@@ -243 +320 @@
     // Execute event loop
     //
+    MProgressBar bar;
     MEvtLoop evtloop;
     evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
 
     Int_t maxevents = 1000000000;
+    //Int_t maxevents = 1000;
     if ( !evtloop.Eventloop(maxevents) )
         return;
@@ -255 +335 @@
     // Display the histograms
     //
-    pliston.FindObject("MHHillas")->DrawClone();
-    pliston.FindObject("MHHillasSrc")->DrawClone();
+    TObject *c;
+    c = pliston.FindObject("MHHillas")->DrawClone();
+
+    c = pliston.FindObject("MHHillasSrc")->DrawClone();
 
     //pliston.FindObject("MHHillasExt")->DrawClone();
-    pliston.FindObject(nxt)->DrawClone();
-
-    pliston.FindObject("MHStarMap")->DrawClone();
+    c = pliston.FindObject(nxt)->DrawClone();
+
+    c = pliston.FindObject("MHStarMap")->DrawClone();
 
 
@@ -277 +359 @@
     Int_t nmaxevts  = 2000;
 
-    // target distribution for the variable in column refcolumn;    
+    // target distribution for the variable in column refcolumn (start at 1);
     // set refcolumn negative if distribution is not to be changed
-    Int_t   refcolumn = 12;
+    Int_t   refcolumn = 1;
     Int_t   nbins   =  10;
     Float_t frombin = 0.5;
@@ -311 +393 @@
     // write out look-alike events (for hadrons)
     //
-    if (WLook)
+    if (WLookOFF)
     {
       TFile *writejens = new TFile(outName, "RECREATE", "");
@@ -324 +406 @@
     //-------------------------------------------
     // Write histograms onto a file
-  if (WHist)
+  if (WHistOFF)
   {
       MHSigmaTheta *sigtheta = 
@@ -330 +412 @@
       if (!sigtheta)
         {
-          *fLog << "Object 'SigmaTheta' not found" << endl;
+          cout << "Object 'SigmaTheta' not found" << endl;
           return;
         }
@@ -337 +419 @@
       TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
 
-      TFile outfile("SigmaTheta.root", "RECREATE");
+      TString outNameSigTh = "SigmaTheta_";
+      outNameSigTh += typeData;
+      outNameSigTh += ".root";
+
+      TFile outfile(outNameSigTh, "RECREATE");
       fHSigTh->Write();
       fHSigPixTh->Write();
@@ -343 +429 @@
       outfile.Close();
      
-      cout << "File 'SigmaTheta.root' was written out" << endl;
+      cout << "File " << outNameSigTh << " was written out" << endl;
   }
 
 
-    cout << "Macro AnalyseCT1 : End of section for ON data" << endl;
+    cout << "Macro AnalyseCT1 : End of Job 1" << endl;
     cout << "===================================================" << endl;
  }
@@ -353 +439 @@
 
   //---------------------------------------------------------------------
-  // start of section for MC gamma data
+  // Job 2
+  //======
+    // read ON data file 
+
+    //  - produce the 2D-histogram "sigmabar versus Theta" for ON data
+    //    (to be used for the padding of the MC gamma data)
+
+    //  - write a file of look-alike events (hadron matrix)
+    //    (to be used later together with the corresponding MC gamma file
+    //     for the g/h separation in the analysis of the ON data)
+
+    //  - write a file of ON events (ON1) 
+    //    (after the standard cuts, before the g/h separation)
+    //    (to be used together with the corresponding MC gamma file
+    //     for the optimization of the g/h separation)
+
+ if (Job2)
+ {
+    cout << "=====================================================" << endl;
+    cout << "Macro AnalyseCT1 : Start of Job 2" << endl;
+
+    cout << "" << endl;
+    cout << "Macro AnalyseCT1 : Job2, WHistON, WLookON, WON1 = " 
+         << Job2  << ",  " << WHistON << ",  " << WLookON << ",  " 
+         << WON1 << endl;
+
+    TString typeData = "ON";
+    cout << "typeData = " << typeData << endl;
+
+    MTaskList tliston;
+    MParList pliston;
+    pliston.AddToList(&tliston);
+
+
+    //::::::::::::::::::::::::::::::::::::::::::
+
+    MBinning binssize("BinningSize");
+    binssize.SetEdgesLog(50, 10, 1.0e5);
+    pliston.AddToList(&binssize);
+
+    MBinning binsdistc("BinningDist");
+    binsdistc.SetEdges(50, 0, 1.4);
+    pliston.AddToList(&binsdistc);
+
+    MBinning binswidth("BinningWidth");
+    binswidth.SetEdges(50, 0, 1.0);
+    pliston.AddToList(&binswidth);
+
+    MBinning binslength("BinningLength");
+    binslength.SetEdges(50, 0, 1.0);
+    pliston.AddToList(&binslength);
+
+    MBinning binsalpha("BinningAlpha");
+    binsalpha.SetEdges(100, -100, 100);
+    pliston.AddToList(&binsalpha);
+
+    MBinning binsasym("BinningAsym");
+    binsasym.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsasym);
+
+    MBinning binsht("BinningHeadTail");
+    binsht.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsht);
+
+    MBinning binsm3l("BinningM3Long");
+    binsm3l.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsm3l);
+
+    MBinning binsm3t("BinningM3Trans");
+    binsm3t.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsm3t);
+
+   
+    //.....
+    MBinning binsb("BinningSigmabar");
+    binsb.SetEdges( 100,  0.0,  5.0);
+    pliston.AddToList(&binsb);
+
+    MBinning binth("BinningTheta");
+    binth.SetEdges( 70, -0.5, 69.5);    
+    pliston.AddToList(&binth);
+
+    MBinning binsdiff("BinningDiffsigma2");
+    binsdiff.SetEdges(100, -5.0, 20.0);
+    pliston.AddToList(&binsdiff);
+    //::::::::::::::::::::::::::::::::::::::::::
+
+
+    //-------------------------------------------
+    // create the tasks which should be executed 
+    //
+
+    filename = onfile;
+    RName = "ReadCT1data";
+    MCT1ReadPreProc read(filename, RName);
+
+    MSelBasic selbasic;
+
+    MBlindPixelCalc blind;
+    blind.SetUseInterpolation();
+    blind.SetUseBlindPixels();
+    // blind.SetUseCetralPixel();
+
+    // create an object of class MSigmabar, 
+    // because class MHSigmaTheta will look for it
+    MSigmabar sigbar;
+    pliston.AddToList(&sigbar);
+    MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
+    fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
+
+    MImgCleanStd    clean; 
+
+    // calculate also extended image parameters
+    TString fHilName = "Hillas";
+    MHillasExt hext(fHilName);
+    pliston.AddToList(&hext);
+    MHillasExt *fHillas = &hext;
+
+    // name of output container for MHillasCalc
+    //      = name of MHillas object
+    MHillasCalc    hcalc(fHilName);
+
+    // name of output container for MHillasSrcCalc
+    //      = name of MHillasSrc object
+    TString fHilSrcName = "HillasSrc";
+    MHillasSrc hilsrc(fHilSrcName);
+    MHillasSrc *fHillasSrc = &hilsrc;
+    pliston.AddToList(&hilsrc);
+    MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
+
+    MSelStandard selstand(fHilName, fHilSrcName);
+
+    MFillH hfill1("MHHillas",    fHilName);
+    MFillH hfill2("MHStarMap",   fHilName);
+
+    TString nxt = "HillasExt";
+    MHHillasExt hhilext(nxt, "", fHilName);
+    pliston.AddToList(&hhilext); 
+    TString namext = nxt;
+    namext += "[MHHillasExt]";
+    MFillH hfill3(namext, fHilSrcName);
+
+    MFillH hfill4("MHHillasSrc", fHilSrcName);
+
+    //+++++++++++++++++++++++++++++++++++++++++++++++++++
+    // prepare writing of look-alike events (for hadrons)
+
+    const char* mtxName = "MatrixHadrons";
+    Int_t howMany = 50000;
+    TString outName = "matrix_hadrons_";
+    outName += typeData;
+    outName += ".root";
+
+
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    cout << "Matrix for (hadrons)" << endl;
+    cout << "input file                    = " << filename<< endl;
+    cout << "matrix name                   = " << mtxName << endl;
+    cout << "max. no.of look-alike events  = " << howMany << endl;
+    cout << "name of output root file      = " << outName << endl;
+    cout << "" << endl;
+
+
+    // matrix limitation for look-alike events (approximate number)
+    MFEventSelector selector("", "", RName);
+    selector.SetNumSelectEvts(howMany);
+    MFilterList flist;
+    flist.AddToList(&selector);
+
+    //
+    // --- setup the matrix and add it to the parlist
+    //
+    MHMatrix *pmatrix = new MHMatrix(mtxName);
+    MHMatrix& matrix = *pmatrix;
+
+    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
+    matrix.AddColumn("MSigmabar.fSigmabar");
+    matrix.AddColumn("log10(Hillas.fSize)");
+    matrix.AddColumn("HillasSrc.fDist");
+    matrix.AddColumn("Hillas.fWidth");
+    matrix.AddColumn("Hillas.fLength");
+    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
+    matrix.AddColumn("abs(Hillas.fM3Long)");
+    matrix.AddColumn("Hillas.fConc");
+    matrix.AddColumn("Hillas.fConc1");
+
+    pliston.AddToList(pmatrix);
+
+    MFillH fillmatg(mtxName);
+    fillmatg.SetFilter(&flist);
+
+
+    if (WON1)
+    {    
+      TString outNameImage = typeData;
+      outNameImage += "1.root";
+      MWriteRootFile write(outNameImage);
+      write.AddContainer("MSigmabar",     "Events");
+      write.AddContainer("Hillas",        "Events");
+      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("HillasSrc",     "Events");
+      write.AddContainer("MRawRunHeader", "RunHeaders");
+      //write.AddContainer("MMcRunHeader","RunHeaders");
+      write.AddContainer("MSrcPosCam",    "RunHeaders");
+    }
+
+    //+++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+    //*****************************
+    // tasks to be executed
+    
+    tliston.AddToList(&read);
+
+    tliston.AddToList(&selbasic);
+    tliston.AddToList(&blind);
+    tliston.AddToList(&fillsigtheta);
+    tliston.AddToList(&clean);
+
+    tliston.AddToList(&hcalc);
+    tliston.AddToList(&csrc1);
+
+    tliston.AddToList(&selstand);
+    if (WON1)
+      tliston.AddToList(&write);
+
+    tliston.AddToList(&flist);
+    tliston.AddToList(&fillmatg);    
+
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+
+
+    //*****************************
+
+    //-------------------------------------------
+    // Execute event loop
+    //
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxevents = 1000000000;
+    //Int_t maxevents = 1000;
+    if ( !evtloop.Eventloop(maxevents) )
+        return;
+
+    tliston.PrintStatistics(0, kTRUE);
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+    TObject *c;
+    c = pliston.FindObject("MHHillas")->DrawClone();
+
+    c = pliston.FindObject("MHHillasSrc")->DrawClone();
+
+    //pliston.FindObject("MHHillasExt")->DrawClone();
+    c = pliston.FindObject(nxt)->DrawClone();
+
+    c = pliston.FindObject("MHStarMap")->DrawClone();
+
+
+    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+    //
+    // ----------------------------------------------------------
+    //  Definition of the reference sample of look-alike events
+    //  (this is a subsample of the original sample)
+    // ----------------------------------------------------------
+    //
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    // Select a maximum of nmaxevts events from the sample of look-alike 
+    // events. They will form the reference sample.
+    Int_t nmaxevts  = 2000;
+
+    // target distribution for the variable in column refcolumn (start at 1);
+    // set refcolumn negative if distribution is not to be changed
+    Int_t   refcolumn = 1;
+    Int_t   nbins   =  10;
+    Float_t frombin = 0.5;
+    Float_t tobin   = 1.0;
+    TH1F *thsh = new TH1F("thsh","target distribution", 
+                           nbins, frombin, tobin);
+    thsh->SetXTitle("cos( \\Theta)");
+    thsh->SetYTitle("Counts");
+    Float_t dbin = (tobin-frombin)/nbins;
+    Float_t lbin = frombin +dbin*0.5;
+    for (Int_t j=1; j<=nbins; j++) 
+    {
+      thsh->Fill(lbin,1.0);
+      lbin += dbin;
+    }
+
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    cout << "Macro AnalyseCT1 : define reference sample for hadrons" << endl; 
+    cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
+         << refcolumn << ",  " << nmaxevts << endl;    
+
+    if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) 
+    {
+      cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl;
+      return;
+    }
+    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+    //-------------------------------------------
+    // write out look-alike events (for hadrons)
+    //
+    if (WLookON)
+    {
+      TFile *writejens = new TFile(outName, "RECREATE", "");
+      matrix.Write();
+      cout << "Macro AnalyseCT1 : matrix of look-alike events for hadrons written onto file "
+           << outName << endl;
+
+      writejens->Close();
+      delete writejens;
+    }
+
+    //-------------------------------------------
+    // Write histograms onto a file
+  if (WHistON)
+  {
+      MHSigmaTheta *sigtheta = 
+            (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
+      if (!sigtheta)
+        {
+          cout << "Object 'SigmaTheta' not found" << endl;
+          return;
+        }
+      TH2D *fHSigTh    = sigtheta->GetSigmaTheta();
+      TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta();
+      TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
+
+      TString outNameSigTh = "SigmaTheta_";
+      outNameSigTh += typeData;
+      outNameSigTh += ".root";
+
+      TFile outfile(outNameSigTh, "RECREATE");
+      fHSigTh->Write();
+      fHSigPixTh->Write();
+      fHDifPixTh->Write();
+      outfile.Close();
+     
+      cout << "File " << outNameSigTh << " was written out" << endl;
+  }
+
+
+    cout << "Macro AnalyseCT1 : End of Job 2" << endl;
+    cout << "===================================================" << endl;
+ }
+
+
+  //---------------------------------------------------------------------
+   // Job 3
+   //======
 
     // read MC gamma data  
@@ -360 +807 @@
     //      (using the 2D-histogram "sigmabar versus Theta" of the OFF data)
 
-    //    - to write a file of look-alike events (for gammas)
+    //    - to write a file of look-alike events (gammas matrix)
     //      (to be used later together with the corresponding hadron file
     //       for the g/h separation in the analysis of the ON data)
 
-    //    - to write a file of padded MC gamma events 
+    //    - to write a file of padded MC gamma events (MC1)
     //      (after the standard cuts, before the g/h separation)
     //      (to be used together with the corresponding hadron file
     //       for the optimization of the g/h separation)
 
-    //    - to write a file of padded MC gamma events (after the final cuts) 
-    //      (to be used for the optimization of the energy estimator
-    //       and for calculating the effective collection areas)
-
- if (MC)
+
+ if (Job3)
  {
-    cout << "=============================================================" 
-         << endl;
-    cout << "Macro : AnalyseCT1 : Start of section for MC gamma data" 
-         << endl;
-
-    MTaskList tlist;
-    MParList plist;
-
-    plist.AddToList(&tlist);
+    cout << "=====================================================" << endl;
+    cout << "Macro AnalyseCT1 : Start of Job 3" << endl;
+
+    cout << "" << endl;
+    cout << "Macro AnalyseCT1 : Job3, WLookMC, WMC1 = " 
+         << Job3  << ",  " << WLookMC << ",  " << WMC1 << endl;
+
+
+    TString typeMC = "MC";
+    cout << "typeMC = " << typeMC << endl;
+    
+    // use for padding sigmabar vs. Theta from ON or OFF data
+    TString typeData = "ON";
+    //TString typeData = "OFF";
+    cout << "typeData = " << typeData << endl;
 
     //------------------------------------
@@ -389 +839 @@
     // and                "3D-ThetaPixSigma"
     // and                "3D-ThetaPixDiff"
-  if (RHist)
-  {
-      cout << "Reading in file 'SigmaTheta.root' " << endl;
-
-      TFile *infile = new TFile("SigmaTheta.root");
+
+
+      TString outNameSigTh = "SigmaTheta_";
+      outNameSigTh += typeData;
+      outNameSigTh += ".root";
+
+
+      cout << "Reading in file " << outNameSigTh << endl;
+
+      TFile *infile = new TFile(outNameSigTh);
       infile->ls();
 
@@ -400 +855 @@
       if (!fHSigmaTheta)
         {
-          *fLog << "Object '2D-ThetaSigmabar' not found on root file" << endl;
+          cout << "Object '2D-ThetaSigmabar' not found on root file" << endl;
           return;
         }
@@ -409 +864 @@
       if (!fHSigmaPixTheta)
         {
-          *fLog << "Object '3D-ThetaPixSigma' not found on root file" << endl;
+          cout << "Object '3D-ThetaPixSigma' not found on root file" << endl;
           return;
         }
@@ -418 +873 @@
       if (!fHSigmaPixTheta)
         {
-          *fLog << "Object '3D-ThetaPixDiff' not found on root file" << endl;
+          cout << "Object '3D-ThetaPixDiff' not found on root file" << endl;
           return;
         }
       cout << "Object '3D-ThetaPixDiff' was read in" << endl;
-  }
-  else
-  {
-      MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
-      if (!sigtheta)
-      {
-        cout << "Object with name 'SigmaTheta' not found" << endl;
-        return;
-      }
-
-      TH2D *fHSigmaTheta    = sigtheta->GetSigmaTheta();
-      TH3D *fHSigmaPixTheta = sigtheta->GetSigmaPixTheta();
-      TH3D *fHDiffPixTheta  = sigtheta->GetDiffPixTheta();
-  }
+
     //------------------------------------
 
+    MTaskList tlist;
+    MParList plist;
+
+    plist.AddToList(&tlist);
 
 
@@ -539 +985 @@
     MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
 
-    MSelStandard selstand(fHillas, fHillasSrc);
+    MSelStandard selstand(fHilName, fHilSrcName);
 
     MFillH hfill1("MHHillas",    fHilName);
@@ -557 +1003 @@
 
     const char* mtxName = "MatrixGammas";
-    Int_t howMany = 30000;
-    char* outName = "matrix_gammas.root";
+    Int_t howMany = 50000;
+
+    TString outName = "matrix_gammas_";
+    outName += typeMC;
+    outName += "_";
+    outName += typeData;
+    outName += ".root";
+
 
     cout << "" << endl;
@@ -581 +1033 @@
     MHMatrix *pmatrix = new MHMatrix(mtxName);
     MHMatrix& matrix = *pmatrix;
+
+    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
+    matrix.AddColumn("MSigmabar.fSigmabar");
+    matrix.AddColumn("log10(Hillas.fSize)");
+    matrix.AddColumn("HillasSrc.fDist");
     matrix.AddColumn("Hillas.fWidth");
     matrix.AddColumn("Hillas.fLength");
-    matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize");
-    matrix.AddColumn("abs(Hillas.fAsym)");
+    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
     matrix.AddColumn("abs(Hillas.fM3Long)");
-    matrix.AddColumn("abs(Hillas.fM3Trans)");
-    matrix.AddColumn("abs(HillasSrc.fHeadTail)");
     matrix.AddColumn("Hillas.fConc");
     matrix.AddColumn("Hillas.fConc1");
-    matrix.AddColumn("HillasSrc.fDist");
-    matrix.AddColumn("log10(Hillas.fSize)");
-    matrix.AddColumn("HillasSrc.fAlpha");
-    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
-    //matrix.AddColumn("MMcEvt.fTheta");
+
     plist.AddToList(pmatrix);
 
@@ -600 +1050 @@
     fillmatg.SetFilter(&flist);
 
+
+    if (WMC1)
+    {    
+      TString outNameImage = typeMC;
+      outNameImage += "_";
+      outNameImage += typeData;
+      outNameImage += "1.root";
+      MWriteRootFile write(outNameImage);
+      write.AddContainer("MSigmabar",     "Events");
+      write.AddContainer("Hillas",        "Events");
+      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("HillasSrc",     "Events");
+      write.AddContainer("MRawRunHeader", "RunHeaders");
+      //write.AddContainer("MMcRunHeader","RunHeaders");
+      write.AddContainer("MSrcPosCam",    "RunHeaders");
+    }
+
+
     //+++++++++++++++++++++++++++++++++++++++++++++++++++
 
-
-    MSelFinal selfinal(fHillas, fHillasSrc);
 
     //*****************************
@@ -619 +1085 @@
 
     tlist.AddToList(&selstand);
+    if (WMC1)
+      tlist.AddToList(&write);
+
+    tlist.AddToList(&flist);
+    tlist.AddToList(&fillmatg);
+
     tlist.AddToList(&hfill1);
     tlist.AddToList(&hfill2);
@@ -624 +1096 @@
     tlist.AddToList(&hfill4);
 
-    tlist.AddToList(&flist);
-    tlist.AddToList(&fillmatg);
-
-    //tlist.AddToList(&selfinal);
+
     //*****************************
 
@@ -634 +1103 @@
     // Execute event loop
     //
+    MProgressBar bar;
     MEvtLoop evtloop;
     evtloop.SetParList(&plist);
+    evtloop.SetProgressBar(&bar);
 
     Int_t maxevents = 1000000000;
+    //Int_t maxevents = 100;
     if ( !evtloop.Eventloop(maxevents) )
         return;
@@ -670 +1142 @@
     // target distribution for the variable in column refcolumn;    
     // set refcolumn negative if distribution is not to be changed
-    Int_t   refcolumn = 12;
+    Int_t   refcolumn = 1;
     Int_t   nbins   =  10;
     Float_t frombin = 0.5;
@@ -711 +1183 @@
     }
 
-    cout << "Macro AnalyseCT1 : End of section for MC gamma data" 
+    cout << "Macro AnalyseCT1 : End of Job 3" 
          << endl;
     cout << "=========================================================" 
          << endl;
  }
+
+
+
+  //---------------------------------------------------------------------
+  // Job 4
+  //======
+
+    // read OFF1 and MC1 data files  
+    //
+    //  - produce Neyman-Pearson plot
+ 
+ if (Job4)
+ {
+    cout << "=====================================================" << endl;
+    cout << "Macro AnalyseCT1 : Start of Job 4" << endl;
+
+    cout << "" << endl;
+    cout << "Macro AnalyseCT1 : Job4 = " << Job4  << endl;
+
+
+    TString typeMC = "MC";
+    cout << "typeMC = " << typeMC << endl;
+
+    // use hadron matrix from ON or OFF data
+    TString typeData = "ON";
+    //TString typeData = "OFF";
+    cout << "typeData = " << typeData << endl;
+
+   //*************************************************************************
+   // read in matrices of look-alike events
+
+    const char* mtxName = "MatrixGammas";
+
+
+    TString outName = "matrix_gammas_";
+    outName += typeMC;
+    outName += "_";
+    outName += typeData;
+    outName += ".root";
+
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    cout << "Get matrix for (gammas)" << endl;
+    cout << "matrix name        = " << mtxName << endl;
+    cout << "name of root file  = " << outName << endl;
+    cout << "" << endl;
+
+
+    // read in the object with the name 'mtxName' from file 'outName'
+    //
+    TFile *fileg = new TFile(outName); 
+
+    MHMatrix matrixg;
+    matrixg.Read(mtxName);
+    pmatrixg = &matrixg;
+
+    delete fileg;
+
+
+    //***************************************************************** 
+
+    const char* mtxName = "MatrixHadrons";
+
+    TString outName = "matrix_hadrons_";
+    outName += typeData;
+    outName += ".root";
+
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    cout << " Get matrix for (hadrons)" << endl;
+    cout << "matrix name        = " << mtxName << endl;
+    cout << "name of root file  = " << outName << endl;
+    cout << "" << endl;
+
+
+    // read in the object with the name mtxName
+    //
+    TFile *fileh = new TFile(outName); 
+
+    MHMatrix matrixh;
+    matrixh.Read(mtxName);
+    pmatrixh = &matrixh;
+
+    delete fileh;
+
+    //***************************************************************** 
+
+    MTaskList tliston;
+    MParList pliston;
+    pliston.AddToList(&tliston);
+
+    pliston.AddToList(pmatrixg);
+    pliston.AddToList(pmatrixh);
+
+    //::::::::::::::::::::::::::::::::::::::::::
+
+    MBinning binssize("BinningSize");
+    binssize.SetEdgesLog(50, 10, 1.0e5);
+    pliston.AddToList(&binssize);
+
+    MBinning binsdistc("BinningDist");
+    binsdistc.SetEdges(50, 0, 1.4);
+    pliston.AddToList(&binsdistc);
+
+    MBinning binswidth("BinningWidth");
+    binswidth.SetEdges(50, 0, 1.0);
+    pliston.AddToList(&binswidth);
+
+    MBinning binslength("BinningLength");
+    binslength.SetEdges(50, 0, 1.0);
+    pliston.AddToList(&binslength);
+
+    MBinning binsalpha("BinningAlpha");
+    binsalpha.SetEdges(100, -100, 100);
+    pliston.AddToList(&binsalpha);
+
+    MBinning binsasym("BinningAsym");
+    binsasym.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsasym);
+
+    MBinning binsht("BinningHeadTail");
+    binsht.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsht);
+
+    MBinning binsm3l("BinningM3Long");
+    binsm3l.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsm3l);
+
+    MBinning binsm3t("BinningM3Trans");
+    binsm3t.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsm3t);
+
+   
+    //.....
+    MBinning binsb("BinningSigmabar");
+    binsb.SetEdges( 100,  0.0,  5.0);
+    pliston.AddToList(&binsb);
+
+    MBinning binth("BinningTheta");
+    binth.SetEdges( 70, -0.5, 69.5);    
+    pliston.AddToList(&binth);
+
+    MBinning binsdiff("BinningDiffsigma2");
+    binsdiff.SetEdges(100, -5.0, 20.0);
+    pliston.AddToList(&binsdiff);
+    //::::::::::::::::::::::::::::::::::::::::::
+
+
+    //-------------------------------------------
+    // create the tasks which should be executed 
+    //
+
+    TString filenameData = typeData;
+    filenameData += "1.root";
+
+    TString filenameMC = typeMC;
+    filenameMC += "_";
+    filenameMC += typeData;
+    filenameMC += "1.root";
+   
+
+    cout << "filenameData = " << filenameData << endl;
+    cout << "filenameMC   = " << filenameMC   << endl;
+
+    MReadMarsFile read("Events", filenameMC);
+    read.AddFile(filenameData);
+    read.DisableAutoScheme();
+
+    //.......................................................................
+    // g/h separation
+
+    MMultiDimDistCalc ghsep;
+    ghsep.SetUseNumRows(25);
+    ghsep.SetUseKernelMethod(kFALSE);
+
+    //.......................................................................
+
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    TString fHilName    = "Hillas"; 
+    TString fHilSrcName = "HillasSrc"; 
+
+    Float_t hadcut = 0.2;
+    MSelFinal selfinalgh(fHilName, fHilSrcName);
+    selfinalgh.SetHadronnessCut(hadcut);
+    selfinalgh.SetAlphaCut(100.0);
+
+    MFillH fillh("MHHadronness");
+
+    MSelFinal selfinal(fHilName, fHilSrcName);
+    selfinal.SetHadronnessCut(hadcut);
+    selfinal.SetAlphaCut(20.0);
+
+    MFillH hfill1("MHHillas",    fHilName);
+    MFillH hfill2("MHStarMap",   fHilName);
+
+    TString nxt = "HillasExt";
+    MHHillasExt hhilext(nxt, "", fHilName);
+    pliston.AddToList(&hhilext); 
+    TString namext = nxt;
+    namext += "[MHHillasExt]";
+    MFillH hfill3(namext, fHilSrcName);
+
+    MFillH hfill4("MHHillasSrc", fHilSrcName);
+
+
+    //*****************************
+    // tasks to be executed
+    
+    tliston.AddToList(&read);
+
+    tliston.AddToList(&ghsep);
+    tliston.AddToList(&fillh);
+   
+    tliston.AddToList(&selfinalgh);
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+
+    tliston.AddToList(&selfinal);
+
+    //*****************************
+
+    //-------------------------------------------
+    // Execute event loop
+    //
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxevents = 1000000000;
+    //Int_t maxevents = 1000;
+    if ( !evtloop.Eventloop(maxevents) )
+        return;
+
+    tliston.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+    TObject *c;
+
+    c = pliston.FindObject("MHHadronness")->DrawClone();
+    c->Print();
+
+    //c = pliston.FindObject("MHHillas")->DrawClone();
+
+    c = pliston.FindObject("MHHillasSrc")->DrawClone();
+
+    //pliston.FindObject("MHHillasExt")->DrawClone();
+    //c = pliston.FindObject(nxt)->DrawClone();
+
+    c = pliston.FindObject("MHStarMap")->DrawClone();
+
+
+
+    cout << "Macro AnalyseCT1 : End of Job 4" << endl;
+    cout << "===================================================" << endl;
+ }
+
+
+  //---------------------------------------------------------------------
+  // Job 5
+  //======
+
+    //  - read in the matrices of look-alike events for gammas and hadrons
+
+    // then read MC1 data file 
+    //
+    //  - perform the g/h separation
+    //  - apply the final cuts
+    //
+    //  - write a file of MC gamma events (MC2)
+    //    (after the g/h separation and after the final cuts)
+
+ if (Job5)
+ {
+    cout << "=====================================================" << endl;
+    cout << "Macro AnalyseCT1 : Start of Job 5" << endl;
+
+    cout << "" << endl;
+    cout << "Macro AnalyseCT1 : Job5, WMC2 = " 
+         << Job5  << ",  " << WMC2 << endl;
+
+    TString typeInput = "MC";
+    cout << "typeInput = " << typeInput << endl;
+
+    TString typeMC   = "MC";
+    cout << "typeMC = " << typeMC << endl;
+
+    // use hadron matrix from ON or OFF data
+    TString typeData = "ON";
+    //TString typeData = "OFF";
+    cout << "typeData = " << typeData << endl;
+
+   //*************************************************************************
+   // read in matrices of look-alike events
+
+    const char* mtxName = "MatrixGammas";
+
+    TString outName = "matrix_gammas_";
+    outName += typeMC;
+    outName += "_";
+    outName += typeData;
+    outName += ".root";
+
+
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    cout << "Get matrix for (gammas)" << endl;
+    cout << "matrix name        = " << mtxName << endl;
+    cout << "name of root file  = " << outName << endl;
+    cout << "" << endl;
+
+
+    // read in the object with the name 'mtxName' from file 'outName'
+    //
+    TFile *fileg = new TFile(outName); 
+
+    MHMatrix matrixg;
+    matrixg.Read(mtxName);
+    pmatrixg = &matrixg;
+
+    delete fileg;
+
+    //***************************************************************** 
+
+    const char* mtxName = "MatrixHadrons";
+
+    TString outName = "matrix_hadrons_";
+    outName += typeData;
+    outName += ".root";
+
+
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    cout << " Get matrix for (hadrons)" << endl;
+    cout << "matrix name        = " << mtxName << endl;
+    cout << "name of root file  = " << outName << endl;
+    cout << "" << endl;
+
+
+    // read in the object with the name mtxName
+    //
+    TFile *fileh = new TFile(outName); 
+
+    MHMatrix matrixh;
+    matrixh.Read(mtxName);
+    pmatrixh = &matrixh;
+
+    delete fileh;
+
+   //*************************************************************************
+
+
+    MTaskList tliston;
+    MParList pliston;
+    pliston.AddToList(&tliston);
+
+    pliston.AddToList(pmatrixg);
+    pliston.AddToList(pmatrixh);
+
+    //::::::::::::::::::::::::::::::::::::::::::
+
+    MBinning binssize("BinningSize");
+    binssize.SetEdgesLog(50, 10, 1.0e5);
+    pliston.AddToList(&binssize);
+
+    MBinning binsdistc("BinningDist");
+    binsdistc.SetEdges(50, 0, 1.4);
+    pliston.AddToList(&binsdistc);
+
+    MBinning binswidth("BinningWidth");
+    binswidth.SetEdges(50, 0, 1.0);
+    pliston.AddToList(&binswidth);
+
+    MBinning binslength("BinningLength");
+    binslength.SetEdges(50, 0, 1.0);
+    pliston.AddToList(&binslength);
+
+    MBinning binsalpha("BinningAlpha");
+    binsalpha.SetEdges(100, -100, 100);
+    pliston.AddToList(&binsalpha);
+
+    MBinning binsasym("BinningAsym");
+    binsasym.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsasym);
+
+    MBinning binsht("BinningHeadTail");
+    binsht.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsht);
+
+    MBinning binsm3l("BinningM3Long");
+    binsm3l.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsm3l);
+
+    MBinning binsm3t("BinningM3Trans");
+    binsm3t.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsm3t);
+
+   
+    //.....
+    MBinning binsb("BinningSigmabar");
+    binsb.SetEdges( 100,  0.0,  5.0);
+    pliston.AddToList(&binsb);
+
+    MBinning binth("BinningTheta");
+    binth.SetEdges( 70, -0.5, 69.5);    
+    pliston.AddToList(&binth);
+
+    MBinning binsdiff("BinningDiffsigma2");
+    binsdiff.SetEdges(100, -5.0, 20.0);
+    pliston.AddToList(&binsdiff);
+    //::::::::::::::::::::::::::::::::::::::::::
+
+
+    //-------------------------------------------
+    // create the tasks which should be executed 
+    //
+
+    TString filenameMC = typeInput;
+    filenameMC += "_";
+    filenameMC += typeData;
+    filenameMC += "1.root";
+
+    readname = "ReadCT1MCdata";
+    MReadMarsFile read("Events", filenameMC);
+    read.DisableAutoScheme();
+
+
+    //.......................................................................
+    // g/h separation
+
+    MMultiDimDistCalc ghsep;
+    ghsep.SetUseNumRows(25);
+    ghsep.SetUseKernelMethod(kFALSE);
+
+
+
+    //.......................................................................
+
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    TString fHilName    = "Hillas"; 
+    TString fHilSrcName = "HillasSrc"; 
+
+    Float_t hadcut = 0.2;
+    MSelFinal selfinalgh(fHilName, fHilSrcName);
+    selfinalgh.SetHadronnessCut(hadcut);
+    selfinalgh.SetAlphaCut(100.0);
+
+    MFillH fillh("MHHadronness");
+
+    MSelFinal selfinal(fHilName, fHilSrcName);
+    selfinal.SetHadronnessCut(hadcut);
+    selfinal.SetAlphaCut(20.0);
+
+    if (WMC2)
+    {    
+      TString outNameImage = typeInput;
+      outNameImage += "_";
+      outNameImage += typeData;
+      outNameImage += "2.root";
+      MWriteRootFile write(outNameImage);
+      write.AddContainer("MHadronness",   "Events");
+      write.AddContainer("MSigmabar",     "Events");
+      write.AddContainer("Hillas",        "Events");
+      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("HillasSrc",     "Events");
+      write.AddContainer("MRawRunHeader", "RunHeaders");
+      //write.AddContainer("MMcRunHeader","RunHeaders");
+      write.AddContainer("MSrcPosCam",    "RunHeaders");
+    }
+
+
+    MFillH hfill1("MHHillas",    fHilName);
+    MFillH hfill2("MHStarMap",   fHilName);
+
+    TString nxt = "HillasExt";
+    MHHillasExt hhilext(nxt, "", fHilName);
+    pliston.AddToList(&hhilext); 
+    TString namext = nxt;
+    namext += "[MHHillasExt]";
+    MFillH hfill3(namext, fHilSrcName);
+
+    MFillH hfill4("MHHillasSrc", fHilSrcName);
+
+
+
+    //*****************************
+    // tasks to be executed
+    
+    tliston.AddToList(&read);
+
+    tliston.AddToList(&ghsep);
+    tliston.AddToList(&fillh);
+
+    tliston.AddToList(&selfinalgh);
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+
+    tliston.AddToList(&selfinal);
+    if (WMC2)
+      tliston.AddToList(&write);
+
+    //*****************************
+
+    //-------------------------------------------
+    // Execute event loop
+    //
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxevents = 1000000000;
+    //Int_t maxevents = 1000;
+    if ( !evtloop.Eventloop(maxevents) )
+        return;
+
+    tliston.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+    TObject *c;
+
+    c = pliston.FindObject("MHHadronness")->DrawClone();
+    c->Print();
+
+    c = pliston.FindObject("MHHillas")->DrawClone();
+
+    c = pliston.FindObject("MHHillasSrc")->DrawClone();
+
+    //pliston.FindObject("MHHillasExt")->DrawClone();
+    c = pliston.FindObject(nxt)->DrawClone();
+
+    c = pliston.FindObject("MHStarMap")->DrawClone();
+
+
+
+    cout << "Macro AnalyseCT1 : End of Job 5" << endl;
+    cout << "===================================================" << endl;
+ }
+
+
+  //---------------------------------------------------------------------
+  // Job 6
+  //======
+
+    //  - read in the matrices of look-alike events for gammas and hadrons
+
+    // then read ON1 data file 
+    //
+    //  - perform the g/h separation
+    //  - apply the final cuts
+    //
+    //  - write a file of ON events (ON2)
+    //    (after the g/h separation and after the final cuts)
+    //    (to be used for the flux calculation)
+    //
+    //  - do the flux calculation
+
+ if (Job6)
+ {
+    cout << "=====================================================" << endl;
+    cout << "Macro AnalyseCT1 : Start of Job 6" << endl;
+
+    cout << "" << endl;
+    cout << "Macro AnalyseCT1 : Job6, WON2 = " 
+         << Job6  << ",  " << WON2 << endl;
+
+    TString typeInput = "ON";
+    cout << "typeInput = " << typeInput << endl;
+
+    TString typeMC   = "MC";
+    cout << "typeMC = " << typeMC << endl;
+
+    // use hadron matrix from ON or OFF data
+    TString typeData = "ON";
+    //TString typeData = "OFF";
+    cout << "typeData = " << typeData << endl;
+
+
+   //*************************************************************************
+   // read in matrices of look-alike events
+
+    const char* mtxName = "MatrixGammas";
+
+    TString outName = "matrix_gammas_";
+    outName += typeMC;
+    outName += "_";
+    outName += typeData;
+    outName += ".root";
+
+
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    cout << "Get matrix for (gammas)" << endl;
+    cout << "matrix name        = " << mtxName << endl;
+    cout << "name of root file  = " << outName << endl;
+    cout << "" << endl;
+
+
+    // read in the object with the name 'mtxName' from file 'outName'
+    //
+    TFile *fileg = new TFile(outName); 
+
+    MHMatrix matrixg;
+    matrixg.Read(mtxName);
+    pmatrixg = &matrixg;
+
+    delete fileg;
+
+    //***************************************************************** 
+
+    const char* mtxName = "MatrixHadrons";
+
+    TString outName = "matrix_hadrons_";
+    outName += typeData;
+    outName += ".root";
+
+
+    cout << "" << endl;
+    cout << "========================================================" << endl;
+    cout << " Get matrix for (hadrons)" << endl;
+    cout << "matrix name        = " << mtxName << endl;
+    cout << "name of root file  = " << outName << endl;
+    cout << "" << endl;
+
+
+    // read in the object with the name mtxName
+    //
+    TFile *fileh = new TFile(outName); 
+
+    MHMatrix matrixh;
+    matrixh.Read(mtxName);
+    pmatrixh = &matrixh;
+
+    delete fileh;
+
+   //*************************************************************************
+
+
+    MTaskList tliston;
+    MParList pliston;
+    pliston.AddToList(&tliston);
+
+    pliston.AddToList(pmatrixg);
+    pliston.AddToList(pmatrixh);
+
+    //::::::::::::::::::::::::::::::::::::::::::
+
+    MBinning binssize("BinningSize");
+    binssize.SetEdgesLog(50, 10, 1.0e5);
+    pliston.AddToList(&binssize);
+
+    MBinning binsdistc("BinningDist");
+    binsdistc.SetEdges(50, 0, 1.4);
+    pliston.AddToList(&binsdistc);
+
+    MBinning binswidth("BinningWidth");
+    binswidth.SetEdges(50, 0, 1.0);
+    pliston.AddToList(&binswidth);
+
+    MBinning binslength("BinningLength");
+    binslength.SetEdges(50, 0, 1.0);
+    pliston.AddToList(&binslength);
+
+    MBinning binsalpha("BinningAlpha");
+    binsalpha.SetEdges(100, -100, 100);
+    pliston.AddToList(&binsalpha);
+
+    MBinning binsasym("BinningAsym");
+    binsasym.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsasym);
+
+    MBinning binsht("BinningHeadTail");
+    binsht.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsht);
+
+    MBinning binsm3l("BinningM3Long");
+    binsm3l.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsm3l);
+
+    MBinning binsm3t("BinningM3Trans");
+    binsm3t.SetEdges(50, -1.5, 1.5);
+    pliston.AddToList(&binsm3t);
+
+   
+    //.....
+    MBinning binsb("BinningSigmabar");
+    binsb.SetEdges( 100,  0.0,  5.0);
+    pliston.AddToList(&binsb);
+
+    MBinning binth("BinningTheta");
+    binth.SetEdges( 70, -0.5, 69.5);    
+    pliston.AddToList(&binth);
+
+    MBinning binsdiff("BinningDiffsigma2");
+    binsdiff.SetEdges(100, -5.0, 20.0);
+    pliston.AddToList(&binsdiff);
+    //::::::::::::::::::::::::::::::::::::::::::
+
+
+    //-------------------------------------------
+    // create the tasks which should be executed 
+    //
+
+    TString filenameData = typeInput;
+    filenameData += "1.root";
+
+    readname = "ReadCT1ONdata";
+    MReadMarsFile read("Events", filenameData);
+    read.DisableAutoScheme();
+
+
+    //.......................................................................
+    // g/h separation
+
+    MMultiDimDistCalc ghsep;
+    ghsep.SetUseNumRows(25);
+    ghsep.SetUseKernelMethod(kFALSE);
+
+
+
+    //.......................................................................
+
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    TString fHilName    = "Hillas"; 
+    TString fHilSrcName = "HillasSrc"; 
+
+    Float_t hadcut = 0.2;
+    MSelFinal selfinalgh(fHilName, fHilSrcName);
+    selfinalgh.SetHadronnessCut(hadcut);
+    selfinalgh.SetAlphaCut(100.0);
+
+    MFillH fillh("MHHadronness");
+
+    MSelFinal selfinal(fHilName, fHilSrcName);
+    selfinal.SetHadronnessCut(hadcut);
+    selfinal.SetAlphaCut(20.0);
+
+    if (WON2)
+    {    
+      TString outNameImage = typeInput;
+      outNameImage += "_";
+      outNameImage += typeData;
+      outNameImage += "2.root";
+      MWriteRootFile write(outNameImage);
+      write.AddContainer("MHadronness",   "Events");
+      write.AddContainer("MSigmabar",     "Events");
+      write.AddContainer("Hillas",        "Events");
+      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("HillasSrc",     "Events");
+      write.AddContainer("MRawRunHeader", "RunHeaders");
+      //write.AddContainer("MMcRunHeader","RunHeaders");
+      write.AddContainer("MSrcPosCam",    "RunHeaders");
+    }
+
+
+    MFillH hfill1("MHHillas",    fHilName);
+    MFillH hfill2("MHStarMap",   fHilName);
+
+    TString nxt = "HillasExt";
+    MHHillasExt hhilext(nxt, "", fHilName);
+    pliston.AddToList(&hhilext); 
+    TString namext = nxt;
+    namext += "[MHHillasExt]";
+    MFillH hfill3(namext, fHilSrcName);
+
+    MFillH hfill4("MHHillasSrc", fHilSrcName);
+
+
+
+    //*****************************
+    // tasks to be executed
+    
+    tliston.AddToList(&read);
+
+    tliston.AddToList(&ghsep);
+    tliston.AddToList(&fillh);
+
+    tliston.AddToList(&selfinalgh);
+
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+
+    tliston.AddToList(&selfinal);
+    if (WON2)
+      tliston.AddToList(&write);
+
+
+
+    //*****************************
+
+    //-------------------------------------------
+    // Execute event loop
+    //
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxevents = 1000000000;
+    //Int_t maxevents = 1000;
+    if ( !evtloop.Eventloop(maxevents) )
+        return;
+
+    tliston.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+    TObject *c;
+
+    c = pliston.FindObject("MHHadronness")->DrawClone();
+    c->Print();
+
+    //c = pliston.FindObject("MHHillas")->DrawClone();
+
+    c = pliston.FindObject("MHHillasSrc")->DrawClone();
+
+    ////pliston.FindObject("MHHillasExt")->DrawClone();
+    //c = pliston.FindObject(nxt)->DrawClone();
+
+    c = pliston.FindObject("MHStarMap")->DrawClone();
+
+
+
+    cout << "Macro AnalyseCT1 : End of Job 6" << endl;
+    cout << "=======================================================" << endl;
+ }
+  //---------------------------------------------------------------------
 }
 
@@ -721 +2043 @@
 
 
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc

@@ -208 +208 @@
     }
 
-    //fHadronness->SetHadronness(dg/(dg+dh));
-    fHadronness->SetHadronness(exp(-dh/dg));
+    Double_t arg;
+
+    if (dg+dh != 0.0)
+      arg = dg / (dg+dh);
+    else
+      arg = 1.e10;
+    //fHadronness->SetHadronness(arg);
+
+    if (dg != 0.0)
+      arg = exp(-dh/dg);
+    else
+      arg = 0.0;
+    fHadronness->SetHadronness(arg);
+      
 
     return kTRUE;

trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc

@@ -1 @@
-
 /* ======================================================================== *\
 !
@@ -102 +101 @@
   fHDiffPixTheta  = fHist3Diff;
 
+  fHSigmaTheta->SetDirectory(NULL);
+  fHSigmaPixTheta->SetDirectory(NULL);
+  fHDiffPixTheta->SetDirectory(NULL);
+
   Print();
 }
@@ -239 +242 @@
    fHDiffPixTh->SetZTitle("Sigma^2-Sigmabar^2");
 
+   //--------------------------------------------------------------------
+
+   memset(fErrors, 0, sizeof(fErrors));
+
    return kTRUE;
 }
@@ -251 +258 @@
   //*fLog << "Entry MPadSchweizer::Process();" << endl;
 
+  Int_t rc;
+
   const UInt_t npix = fEvt->GetNumPixels();
+
+  Double_t SigbarOld;
+
+  //*fLog << "before padding : " << endl;
+  //SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);
+  //fSigmabar->Print("");
+
 
   //$$$$$$$$$$$$$$$$$$$$$$$$$$
   // to simulate the situation that before the padding the NSB and 
   // electronic noise are zero : set Sigma = 0 for all pixels
-  for (UInt_t i=0; i<npix; i++) 
-  {   
-    MCerPhotPix &pix = fEvt->operator[](i);      
-    Int_t j = pix.GetPixId();
-
-    MPedestalPix &ppix = fPed->operator[](j);
-    ppix.SetMeanRms(0.0);
-  }
+  //for (UInt_t i=0; i<npix; i++) 
+  //{   
+  //  MCerPhotPix &pix = fEvt->operator[](i);      
+  //  Int_t j = pix.GetPixId();
+
+  //  MPedestalPix &ppix = fPed->operator[](j);
+  //  ppix.SetMeanRms(0.0);
+  //}
   //$$$$$$$$$$$$$$$$$$$$$$$$$$
 
@@ -269 +285 @@
   // Calculate average sigma of the event
   //
-  Double_t SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);
+  SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);
   //fSigmabar->Print("");
 
-  //if (SigbarOld > 0.0)
-  //{
+  if (SigbarOld > 0.0)
+  {
     //*fLog << "MPadSchweizer::Process(); Sigmabar of event to be padded is > 0 : "
     //      << SigbarOld << ". Stop event loop " << endl;
     // input data should have Sigmabar = 0; stop event loop
-    // return kFALSE; 
-  //}
+  
+    rc = 1;
+    fErrors[rc]++;
+    return kCONTINUE; 
+  }
 
   Double_t Theta  = kRad2Deg*fMcEvt->GetTelescopeTheta();
@@ -296 +315 @@
           << Theta << ",  " << binNumber << ";  Skip event " << endl;
     // event cannot be padded; skip event
+
+    rc = 2;
+    fErrors[rc]++;
     return kCONTINUE;
   }
@@ -308 +330 @@
       // event cannot be padded; skip event
       delete fHSigma;
+
+      rc = 3;
+      fErrors[rc]++;
       return kCONTINUE;
     }
@@ -331 +356 @@
     *fLog << "MPadSchweizer::Process(); target Sigmabar is less than SigbarOld : "
           << Sigmabar << ",  " << SigbarOld << ",   Skip event" << endl;
+
+    rc = 4;
+    fErrors[rc]++;
     return kCONTINUE;     
   }
@@ -426 +454 @@
               << binTheta << "  and pixel bin " << binPixel  
               << " has no entries;  aborting " << endl;
-        return kERROR;
+
+        rc = 5;
+        fErrors[rc]++;
+        return kCONTINUE;     
       }
 
@@ -465 +496 @@
               << binTheta << "  and pixel bin " << binPixel  
               << " has no entries;  aborting " << endl;
-        return kERROR;
+
+        rc = 6;
+        fErrors[rc]++;
+        return kCONTINUE;     
       }
 
@@ -553 +587 @@
 
   // Calculate Sigmabar again and crosscheck
+
   Double_t SigbarNew = fSigmabar->Calc(*fCam, *fPed, *fEvt);
+  //*fLog << "after padding : " << endl;
   //fSigmabar->Print("");
 
@@ -585 +621 @@
   //*fLog << "Exit MPadSchweizer::Process();" << endl;
 
+  rc = 0;
+  fErrors[rc]++;
+
   return kTRUE;
 
@@ -594 +633 @@
 Bool_t MPadSchweizer::PostProcess()
 {
+    if (GetNumExecutions() != 0)
+    {
+
+    *fLog << inf << endl;
+    *fLog << GetDescriptor() << " execution statistics:" << endl;
+    *fLog << dec << setfill(' ');
+    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
+          << (int)(fErrors[1]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
+
+    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) 
+          << (int)(fErrors[2]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: Zenith angle out of range" << endl;
+
+    *fLog << " " << setw(7) << fErrors[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 << " " << fErrors[0] << " (" 
+          << (int)(fErrors[0]*100/GetNumExecutions()) 
+          << "%) Evts survived the padding!" << endl;
+    *fLog << endl;
+
+    }
+
+    //---------------------------------------------------------------
     TCanvas &c = *(MH::MakeDefCanvas("PadSchweizer", "", 900, 1200)); 
     c.Divide(3, 4);
@@ -600 +677 @@
 
     c.cd(1);
+    fHSigmaTheta->SetDirectory(NULL);
     fHSigmaTheta->SetTitle("2D : Sigmabar, \\Theta (reference sample)");
-    fHSigmaTheta->DrawClone();     
+    fHSigmaTheta->DrawCopy();     
     fHSigmaTheta->SetBit(kCanDelete);     
 
       //c.cd(1);
-      //fHSigmaPixTheta->DrawClone();     
+      //fHSigmaPixTheta->DrawCopy();     
       //fHSigmaPixTheta->SetBit(kCanDelete);     
 
     c.cd(4);
-    fHSigbarTheta->DrawClone();     
+    fHSigbarTheta->SetDirectory(NULL);
+    fHSigbarTheta->DrawCopy();     
     fHSigbarTheta->SetBit(kCanDelete);     
 
 
     c.cd(7);
-    fHNSB->DrawClone();
+    fHNSB->SetDirectory(NULL);
+    fHNSB->DrawCopy();
     fHNSB->SetBit(kCanDelete);    
 
@@ -625 +705 @@
     c.cd(2);
     l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zx");
-
+    l->SetDirectory(NULL);
     l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
     l->SetXTitle("\\Theta [\\circ]");
     l->SetYTitle("Sigma^2-Sigmabar^2");
 
-    *fLog << "before box" << endl;
-
-    l->Draw("box");
+    l->DrawCopy("box");
     l->SetBit(kCanDelete);;
 
     c.cd(5);
     l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zy");
+    l->SetDirectory(NULL);
     l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
     l->SetXTitle("pixel");
     l->SetYTitle("Sigma^2-Sigmabar^2");
 
-    l->Draw("box");
+    l->DrawCopy("box");
     l->SetBit(kCanDelete);;
 
@@ -655 +734 @@
     c.cd(3);
     k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zx");
+    k->SetDirectory(NULL);
     k->SetTitle("Sigma vs. \\Theta (all pixels)");
     k->SetXTitle("\\Theta [\\circ]");
     k->SetYTitle("Sigma");
 
-    k->Draw("box");
+    k->DrawCopy("box");
     k->SetBit(kCanDelete);
 
     c.cd(6);
     k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zy");
+    k->SetDirectory(NULL);
     k->SetTitle("Sigma vs. pixel number (all \\Theta)");
     k->SetXTitle("pixel");
     k->SetYTitle("Sigma");
 
-    k->Draw("box");
+    k->DrawCopy("box");
     k->SetBit(kCanDelete);;
 
@@ -677 +758 @@
 
     c.cd(10);
-    fHSigmaPedestal->DrawClone();
+    fHSigmaPedestal->SetDirectory(NULL);
+    fHSigmaPedestal->DrawCopy();
     fHSigmaPedestal->SetBit(kCanDelete);    
 
     c.cd(11);
-    fHPhotons->DrawClone();
+    fHPhotons->SetDirectory(NULL);
+    fHPhotons->DrawCopy();
     fHPhotons->SetBit(kCanDelete);    
 

trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h

@@ -39 +39 @@
     Int_t          fRunType;
     Int_t          fGroup;
+
+    Int_t          fErrors[7];
 
     // plots used for the padding

trunk/MagicSoft/Mars/manalysis/MSelBasic.cc

@@ -44 +44 @@
 #include "MGeomCam.h"
 #include "MPedestalCam.h"
+#include "MPedestalPix.h"
 #include "MGeomPix.h"
 
@@ -119 +120 @@
 Bool_t MSelBasic::Process()
 {
+  /*
+  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+  *fLog << "=========================================================" << endl;
+  *fLog << "" << endl;
+  *fLog << "MmcEvt data : " << endl;
+  *fLog << "" << endl;
+  fMcEvt->Print();
+  *fLog << "" << endl;
+  *fLog << "PartId() = " << fMcEvt->GetPartId() << endl;
+  *fLog << "Energy() = " << fMcEvt->GetEnergy() << endl;
+  *fLog << "Theta()  = " << fMcEvt->GetTheta() << endl;
+
+  *fLog << "Phi()    = " << fMcEvt->GetPhi() << endl;
+  //*fLog << "CoreD()  = " << fMcEvt->GetCoreD() << endl;
+  //*fLog << "CoreX()  = " << fMcEvt->GetCoreX() << endl;
+
+  //*fLog << "CoreY()  = " << fMcEvt->GetCoreY() << endl;
+  *fLog << "Impact() = " << fMcEvt->GetImpact() << endl;
+  //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl;
+
+  //*fLog << "PassPhotAtm() = " << fMcEvt->GetPassPhotAtm() << endl;
+  //*fLog << "PassPhotRef() = " << fMcEvt->GetPassPhotRef() << endl;
+  //*fLog << "PassPhotCone()  = " << fMcEvt->GetPassPhotCone() << endl;
+
+  //*fLog << "PhotElfromShower()    = " << fMcEvt->GetPhotElfromShower() << endl;
+  //*fLog << "PhotElinCamera()  = " << fMcEvt->GetPhotElinCamera() << endl;
+  *fLog << "TelescopePhi()  = " << fMcEvt->GetTelescopePhi() << endl;
+
+  *fLog << "TelescopeTheta()  = " << fMcEvt->GetTelescopeTheta() << endl;
+  *fLog << "Impact() = " << fMcEvt->GetImpact() << endl;
+  //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl;
+  *fLog << "" << endl;
+  *fLog << "=========================================================" << endl;
+
+  *fLog << "=========================================================" << endl;
+  *fLog << "" << endl;
+  *fLog << "MPedestalPix data : " << endl;
+  *fLog << "" << endl;
+
+  Int_t ntot;
+  ntot = fPed->GetSize();
+    *fLog << "MeanRms() :" << endl;    
+  for (Int_t i=0; i<ntot; i++)
+  {
+    MPedestalPix &pix = (*fPed)[i];
+    //*fLog << "Mean()    = " << i << ",  " << pix.GetMean() << endl;    
+    //*fLog << "Sigma()   = " << i << ",  " << pix.GetSigma() << endl;    
+    *fLog << pix.GetMeanRms() << " ";    
+    //*fLog << "SigmaRms()= " << i << ",  " << pix.GetSigmaRms() << endl;    
+  }
+  *fLog << "" << endl;
+
+  *fLog << "" << endl;
+  ntot = fEvt->GetNumPixels();
+    *fLog << "MCerPhotPix :  pix.GetNumPhotons()" << endl;    
+  for (Int_t i=0; i<ntot; i++)
+  {
+    MCerPhotPix &pix = (*fEvt)[i];
+    *fLog << pix.GetNumPhotons() << " ";    
+  }
+  *fLog << "" << endl;
+
+  *fLog << "" << endl;
+  ntot = fEvt->GetNumPixels();
+    *fLog << "MCerPhotPix :  pix.GetErrorPhot()" << endl;    
+  for (Int_t i=0; i<ntot; i++)
+  {
+    MCerPhotPix &pix = (*fEvt)[i];
+    *fLog << pix.GetErrorPhot() << " ";    
+  }
+  *fLog << "" << endl;
+
+  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+  */
+
+
     Int_t rc = 0;
 
-    //if ( fRawRun->GetRunNumber() < 1025 ) 
+    //if ( fRawRun->GetRunNumber() < 16279 ) 
     //{
     //   rc = 1;
@@ -128 +205 @@
 
     Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
-    if (Theta > 45.0  ||  !SwTrigger() )
-    {
-      //*fLog << "MSelBasic::Process; Theta = " << Theta << endl;
+    if ( Theta < 0.0 )
+    {
+      *fLog << "MSelBasic::Process; Run, Event, Theta = " 
+            << fRawRun->GetRunNumber()<< ",  " 
+            << fMcEvt->GetEvtNumber() << ",  " << Theta << endl;
       rc = 1;
+    }    
+
+    else if ( Theta > 45.0 )
+    {
+      rc = 2;
+    }    
+
+    else if ( !SwTrigger() )
+    {
+      //*fLog << "MSelBasic::Process; SwTrigger = " << SwTrigger << endl;
+      rc = 3;
     }    
 
@@ -223 +313 @@
     *fLog << GetDescriptor() << " execution statistics:" << endl;
     *fLog << dec << setfill(' ');
-    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Basic selections are not fullfilled" << endl;
+    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Zenith angle < 0" << endl;
+
+    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) << (int)(fErrors[2]*100/GetNumExecutions()) << "%) Evts skipped due to: Zenith angle too high" << endl;
+
+    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) << (int)(fErrors[3]*100/GetNumExecutions()) << "%) Evts skipped due to: Software trigger not fullfilled" << endl;
 
     *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Basic selections!" << endl;
@@ -230 +324 @@
     return kTRUE;
 }
+

trunk/MagicSoft/Mars/manalysis/MSelBasic.h

@@ -30 +30 @@
     const MRawRunHeader *fRawRun;       
 
-          Int_t        fErrors[2];
+          Int_t        fErrors[4];
 
 public:

trunk/MagicSoft/Mars/manalysis/MSelFinal.cc

@@ -41 +41 @@
 
 #include "MHillas.h"
+#include "MHillasExt.h"
 #include "MHillasSrc.h"
 #include "MCerPhotEvt.h"
@@ -57 +58 @@
 // Default constructor.
 //
-MSelFinal::MSelFinal(MHillas *parhil, MHillasSrc *parhilsrc,
+MSelFinal::MSelFinal(const char *HilName, const char *HilSrcName,
                      const char *name, const char *title)
 {
@@ -63 +64 @@
     fTitle = title ? title : "Task to evaluate the Final Cuts";
 
-    fHil    = parhil;
-    fHilsrc = parhilsrc;
+    fHilName    = HilName;
+    fHilSrcName = HilSrcName;
+
+    fHadronnessCut =   0.2;
+    fAlphaCut      = 100.0; //degrees
 }
 
@@ -75 +79 @@
 Bool_t MSelFinal::PreProcess(MParList *pList)
 {
+    fHil    = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt");
+    if (!fHil)
+    {
+      *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl;
+      return kFALSE;
+    }
+
+    fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
+    if (!fHilSrc)
+    {
+      *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl;
+      return kFALSE;
+    }
+
+
     fHadronness = (MHadronness*)pList->FindObject("MHadronness");
     if (!fHadronness)
@@ -90 +109 @@
     }
 
-    fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
-    if (!fEvt)
-    {
-        *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
-        return kFALSE;
-    }
-
-
-    fCam = (MGeomCam*)pList->FindObject("MGeomCam");
-    if (!fCam)
-    {
-        *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
-        return kFALSE;
-    }
-    fMm2Deg = fCam->GetConvMm2Deg();
-
-    *fLog << "fMm2Deg = " << fMm2Deg << endl;
-
     memset(fErrors, 0, sizeof(fErrors));
 
@@ -122 +123 @@
 Bool_t MSelFinal::Process()
 {
+  //*fLog << "Entry MSelFinal; fHilSrc = " << fHilSrc << endl;
+
+    
+
     Int_t rc = 0;
 
-    Double_t alphacut = 20.0;
+    Double_t modalpha = fabs( fHilSrc->GetAlpha() ); 
 
-    Double_t modalpha = fabs( fHilsrc->GetAlpha() ); 
     Double_t h = fHadronness->GetHadronness();
 
-    if ( h>0.5  ||  modalpha > alphacut )
+    if ( h>fHadronnessCut )
     {
       //*fLog << "MSelFinal::Process; h, alpha = " << h << ",  " 
-      //      << fHilsrc->GetAlpha() << endl;
+      //      << fHilSrc->GetAlpha() << endl;
       rc = 1;
+    }    
+
+    else if ( modalpha > fAlphaCut )
+    {
+      //*fLog << "MSelFinal::Process; h, alpha = " << h << ",  " 
+      //      << fHilSrc->GetAlpha() << endl;
+      rc = 2;
     }    
 
@@ -153 +164 @@
     *fLog << GetDescriptor() << " execution statistics:" << endl;
     *fLog << dec << setfill(' ');
-    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Final selections are not fullfilled" << endl;
+    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
+          << (int)(fErrors[1]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: g/h separation cut (" << fHadronnessCut
+          << ")" << endl;
 
-    *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Final selections!" << endl;
+    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) 
+          << (int)(fErrors[2]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: cut in ALPHA (" << fAlphaCut
+          << " degrees)" << endl;
+
+    *fLog << " " << fErrors[0] << " (" 
+          << (int)(fErrors[0]*100/GetNumExecutions()) 
+          << "%) Evts survived Final selections!" << endl;
     *fLog << endl;
 
     return kTRUE;
 }
+
+
+

trunk/MagicSoft/Mars/manalysis/MSelFinal.h

@@ -28 +28 @@
     MMcEvt      *fMcEvt;       
     MHillas     *fHil;       
-    MHillasSrc  *fHilsrc;       
+    MHillasSrc  *fHilSrc;       
     MHadronness *fHadronness;       
 
     Double_t     fMm2Deg;   // conversion mm to degrees in camera
-    Int_t        fErrors[2];
+    Int_t        fErrors[3];
+    TString      fHilName;
+    TString      fHilSrcName;
+ 
+    Float_t      fHadronnessCut;
+    Float_t      fAlphaCut;
 
 public:
-    MSelFinal(MHillas *fHil, MHillasSrc *fHilsrc,
+    MSelFinal(const char *HilName, const char *HilSrcName,
               const char *name=NULL, const char *title=NULL);
 
@@ -42 +47 @@
     Bool_t PostProcess();
 
+    void SetHadronnessCut(Float_t hadcut) { fHadronnessCut = hadcut; }
+    void SetAlphaCut(Float_t alpha)       { fAlphaCut      = alpha;  }
 
     ClassDef(MSelFinal, 0)   // Task to evaluate final cuts

trunk/MagicSoft/Mars/manalysis/MSelStandard.cc

@@ +1 @@
+
 /* ======================================================================== *\
 !
@@ -39 +40 @@
 
 #include "MHillas.h"
+#include "MHillasExt.h"
 #include "MHillasSrc.h"
 #include "MCerPhotEvt.h"
@@ -54 +56 @@
 // Default constructor.
 //
-MSelStandard::MSelStandard(const MHillas *parhil, const MHillasSrc *parhilsrc,
+MSelStandard::MSelStandard(const char *HilName, const char *HilSrcName,
                            const char *name, const char *title)
 {
@@ -60 +62 @@
     fTitle = title ? title : "Task to evaluate the Standard Cuts";
 
-    fHil    = parhil;
-    fHilsrc = parhilsrc;
+    fHilName    = HilName;
+    fHilSrcName = HilSrcName;
 }
 
@@ -72 +74 @@
 Bool_t MSelStandard::PreProcess(MParList *pList)
 {
+    fHil    = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt");
+    if (!fHil)
+    {
+      *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl;
+      return kFALSE;
+    }
+
+    fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
+    if (!fHilSrc)
+    {
+      *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl;
+      return kFALSE;
+    }
+
+
     fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
     if (!fMcEvt)
@@ -95 +112 @@
     fMm2Deg = fCam->GetConvMm2Deg();
 
-    *fLog << "fMm2Deg = " << fMm2Deg << endl;
+    //*fLog << "fMm2Deg = " << fMm2Deg << endl;
 
     memset(fErrors, 0, sizeof(fErrors));
@@ -114 +131 @@
     Int_t rc = 0;
 
-    //Double_t fLength       = fHil->GetLength() * fMm2Deg;
-    //Double_t fWidth        = fHil->GetWidth()  * fMm2Deg;
-    Double_t fDist         = fHilsrc->GetDist()* fMm2Deg;
+    Double_t fLength       = fHil->GetLength() * fMm2Deg;
+    Double_t fWidth        = fHil->GetWidth()  * fMm2Deg;
+    Double_t fDist         = fHilSrc->GetDist()* fMm2Deg;
     //Double_t fDelta        = fHil->GetDelta()  * kRad2Deg;
     Double_t fSize         = fHil->GetSize();
@@ -122 +139 @@
     Int_t fNumCorePixels   = fHil->GetNumCorePixels();
 
-    if (     fSize <= 60.0         ||  fDist< 0.4           ||  fDist > 1.1  
-         ||  fNumUsedPixels >= 92  ||  fNumCorePixels < 4)
+    if ( fNumUsedPixels >= 92  ||  fNumCorePixels < 4 )
     {
       //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = "
@@ -129 +145 @@
       //      << fNumCorePixels << endl;
       rc = 1;
+    }    
+
+    else if ( fSize <= 60.0         ||  fDist< 0.4           ||  fDist > 1.1 )
+    {
+      //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = "
+      //      << fSize << ",  " << fDist << ",  " << fNumUsedPixels << ",  "
+      //      << fNumCorePixels << endl;
+      rc = 2;
+    }    
+
+    else if ( fLength <= 0.0         ||  fWidth <= 0.0 )
+    {
+      //*fLog << "MSelStandard::Process; fLength, fWidth = "
+      //      << fLength << ",  " << fWidth << endl;
+      rc = 3;
     }    
 
@@ -149 +180 @@
     *fLog << GetDescriptor() << " execution statistics:" << endl;
     *fLog << dec << setfill(' ');
-    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Standard selections are not fullfilled" << endl;
+    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Requirements on no.of used or core pxels not fullfilled" << endl;
+
+    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) << (int)(fErrors[2]*100/GetNumExecutions()) << "%) Evts skipped due to: Requirements on SIZE or DIST not fullfilled" << endl;
+
+    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) << (int)(fErrors[3]*100/GetNumExecutions()) << "%) Evts skipped due to: Length or Width is <= 0" << endl;
 
     *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Standard selections!" << endl;
@@ -156 +191 @@
     return kTRUE;
 }
+
+

trunk/MagicSoft/Mars/manalysis/MSelStandard.h

@@ -23 +23 @@
 {
 private:
-    const MGeomCam    *fCam;      // Camera Geometry 
-    const MCerPhotEvt *fEvt;      // Cerenkov Photon Event 
-    const MMcEvt      *fMcEvt;       
-    const MHillas     *fHil;       
-    const MHillasSrc  *fHilsrc;       
+    MGeomCam    *fCam;      // Camera Geometry 
+    MCerPhotEvt *fEvt;      // Cerenkov Photon Event 
+    MMcEvt      *fMcEvt;       
+    MHillas     *fHil;       
+    MHillasSrc  *fHilSrc;       
 
-          Double_t     fMm2Deg;   // conversion mm to degrees in camera
-          Int_t        fErrors[2];
+    Double_t     fMm2Deg;   // conversion mm to degrees in camera
+    Int_t        fErrors[4];
+    TString      fHilName;
+    TString      fHilSrcName; 
 
 public:
-    MSelStandard(const MHillas *fHil, const MHillasSrc *fHilsrc,
+    MSelStandard(const char *HilName, const char *HilSrcName,
                  const char *name=NULL, const char *title=NULL);
 

trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc

@@ -122 +122 @@
 {
     const char *name = gSystem->ExpandPathName(txt);
+
     TString fname(name);
     delete [] name;
@@ -336 +337 @@
     // float   frms_pedsig_phot[iMAXNUMPIX];      // standard deviation of the calibrated signals from the pedestal run */
     fPedest->InitSize(iMAXNUMPIX);
+    *fLog << "PedestalRMS : ";
     for (Int_t i=0; i<iMAXNUMPIX; i++)
+    {
         (*fPedest)[i].SetMeanRms(outpars.frms_pedsig_phot[i]);
+        *fLog << outpars.frms_pedsig_phot[i] << " "; 
+ 
+        //$$$$$$$$$$$$$$$$$$$$$$$$$
+        savePedRMS[i] =  outpars.frms_pedsig_phot[i];
+        //$$$$$$$$$$$$$$$$$$$$$$$$$
+    }
+    *fLog << endl;
 
     fPedest->SetReadyToSave();
@@ -373 +383 @@
      *         be treated further. */
 
+    *fLog << "outpars.bmontecarlo = " << outpars.bmontecarlo << endl;
+    *fLog << "outpars.imcparticle = " << outpars.imcparticle << endl;
+    *fLog << "outpars.dsourcera_hours = " << outpars.dsourcera_hours << endl;
+    *fLog << "outpars.dsourcedec_deg = " << outpars.dsourcedec_deg << endl;
+
+    //*fLog << "File is a ";
+    //*fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data");
+    //*fLog << " file." << endl;
+
+
+    // Next statement commented out because bmontecarlo was set wrongly
+    //fIsMcFile = outpars.bmontecarlo==TRUE;
+    fIsMcFile = (outpars.dsourcera_hours == 0.0 &&
+                 outpars.dsourcedec_deg  == 0.0 &&
+                 outpars.imcparticle != 0          );
+
     *fLog << "File is a ";
-    *fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data");
+    *fLog << (fIsMcFile ? "Monte Carlo" : "Real Data");
     *fLog << " file." << endl;
-
-    fIsMcFile = outpars.bmontecarlo==TRUE;
-
+    *fLog << " " << endl;
 }
 
@@ -433 +457 @@
     ProcessRunHeader(outpars);
 
+
     //rwagner: ReInit whenever new run commences
     // rc==-1 means: ReInit didn't work out
+
     if (!fTaskList->ReInit(fParList))
         return -1;
@@ -453 +479 @@
     TString m = cEND_EVENTS_TEMPLATE;
     Int_t p = m.First('%');
+
 
     if (!s.BeginsWith(m(0,p)))
@@ -515 +542 @@
 
     //
-    // Check for the existance of a next file to read
+    // Check for the existence of a next file to read
     //
     TNamed *file = (TNamed*)fFileNames->First();
@@ -538 +565 @@
 
     if (!CheckHeader(fname))
-        return kFALSE;
+    {
+        *fLog << "OpenNextFile : CheckHeader(fname) is FALSE" << endl;
+        return kFALSE;
+    }
 
     fIn = new ifstream(fname);
 
     *fLog << inf << "-----------------------------------------------------------------------" << endl;
+
 
     switch (ReadRunHeader())
@@ -553 +584 @@
         return kFALSE;
       default:
+        *fLog << "OpenNextFile : after ReadRunHeader, FnumEventsInRun = " 
+              << fNumEventsInRun << endl;
         return kTRUE;
       }
+
+
+
 }
 
@@ -752 +788 @@
 void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event)
 {
+  if (fRawRunHeader->GetRunNumber() == 1)
+  {
+  *fLog << "eventrecord" << endl;
+  *fLog << "isecs_since_midday = " << event.isecs_since_midday << endl;
+  *fLog << "isecfrac_200ns     = " << event.isecfrac_200ns << endl;
+  *fLog << "snot_ok_flags      = " << event.snot_ok_flags << endl;
+  *fLog << "ialt_arcs          = " << event.ialt_arcs << endl;
+  *fLog << "iaz_arcs           = " << event.iaz_arcs << endl;
+  *fLog << "ipreproc_alt_arcs  = " << event.ipreproc_alt_arcs << endl;
+  *fLog << "ipreproc_az_arcs   = " << event.ipreproc_az_arcs << endl;
+  *fLog << "ifieldrot_arcs     = " << event.ifieldrot_arcs << endl;
+
+  *fLog << "srate_millihz      = " << event.srate_millihz << endl;
+  *fLog << "fhourangle         = " << event.fhourangle << endl;
+  *fLog << "fmcenergy_tev      = " << event.fmcenergy_tev << endl;
+  *fLog << "fmcsize_phel       = " << event.fmcsize_phel << endl;
+  *fLog << "imcimpact_m        = " << event.imcimpact_m << endl;
+  *fLog << "imcparticle        = " << event.imcparticle << endl;
+  *fLog << "imctriggerflag     = " << event.imctriggerflag << endl;
+  }  
+
+
+    for (Int_t i=0; i<iMAXNUMPIX; i++)
+    {
+        (*fPedest)[i].SetMeanRms(savePedRMS[i]);
+    }
+
+
     //  int   isecs_since_midday; // seconds passed since midday before sunset (JD of run start)
     //  int   isecfrac_200ns;     // fractional part of isecs_since_midday
@@ -788 +852 @@
     // actual number of pixels (outputpars.inumpixels) is written out
     // short spixsig_10thphot[iMAXNUMPIX];
+    //*fLog << "spixsig_10thphot : " << endl;
     for (Int_t i=0; i<iMAXNUMPIX; i++)
     {
+      //*fLog << event.spixsig_10thphot[i] << " ";
+
         if (event.spixsig_10thphot[i]==0)
             continue;
@@ -796 +863 @@
                          (*fPedest)[i].GetMeanRms());
     }
+    //*fLog << "" << endl;
+
     fNphot->SetReadyToSave();
 
@@ -801 +870 @@
     // int ipreproc_az_arcs;  // "should be" az according to preproc (arcseconds)
 
-    fMcEvt->Fill(0, /*fEvtNum*/
+    fMcEvt->Fill(event.isecs_since_midday,     //0, /*fEvtNum*/
                  fIsMcFile ? event.imcparticle : 0, /*corsika particle type*/
                  fIsMcFile ? event.fmcenergy_tev*1000 : 0,
@@ -867 +936 @@
         // an event
         //
-        if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0])
+
+      // "goto TONI" was introduced in order to check for a footer record
+      // after reading a run header; this is necessary for runs with
+      // zero events after the filter
+      TONI:
+
+      if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0])
             return kTRUE;
 
@@ -874 +949 @@
         // must be an event
         //
+
         switch (ReadRunFooter())
         {
@@ -906 +982 @@
         *fLog << "-----------------------------------------------------------------------" << endl;
 
+
         if (ReadRunHeader() < 0)
         {
@@ -911 +988 @@
             return kFALSE;
         }
+
+        goto TONI;
+
         return kTRUE;
     }

trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h

@@ -48 +48 @@
     UInt_t fNumFilterEvts;  // number of events mentioned in the runs footers
 
+    Float_t savePedRMS[127];
+
+
     Bool_t OpenNextFile();
 

trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc

@@ -118 +118 @@
 Bool_t MFEventSelector::PreProcess(MParList *plist)
 {
+    memset(fErrors, 0, sizeof(fErrors));
+
     fNumSelectedEvts = 0;
     if (fSelRatio>0)
@@ -153 +155 @@
 Bool_t MFEventSelector::Process()
 {
+    Int_t rc;
+
     const Float_t evt = gRandom->Uniform();
 
@@ -161 +165 @@
 
     if (fResult)
+    {
         fNumSelectedEvts++;
 
+        rc = 0;
+        fErrors[rc]++;
+        return kTRUE;
+    }
+    
+    rc = 1;
+    fErrors[rc]++;
+
     return kTRUE;
 }
@@ -172 +185 @@
 Bool_t MFEventSelector::PostProcess()
 {
+    //---------------------------------
+    if (GetNumExecutions() != 0)
+    {
+      *fLog << inf << endl;
+      *fLog << GetDescriptor() << " execution statistics:" << endl;
+      *fLog << dec << setfill(' ');
+      *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
+            << (int)(fErrors[1]*100/GetNumExecutions()) 
+            << "%) Events not selected" << endl;
+
+      *fLog << " " << fErrors[0] << " (" 
+            << (int)(fErrors[0]*100/GetNumExecutions()) 
+            << "%) Events selected!" << endl;
+      *fLog << endl;
+    }
+
+    //---------------------------------
     if (TestBit(kNumTotalFromFile))
         fNumTotalEvts = -1;
@@ -199 +229 @@
 
 }
-

trunk/MagicSoft/Mars/mfilter/MFEventSelector.h

@@ -37 +37 @@
      */
 
+    Int_t fErrors[2];
+
 public:
     // MFEventSelector();

trunk/MagicSoft/Mars/mhist/MHHadronness.cc

@@ -235 +235 @@
     }
 
-    return val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
+    Float_t retValue;
+    if (val2x-val1x != 0.0)
+      retValue = val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
+    else
+      retValue = 0.0;
+
+    return retValue;
 }
 
@@ -252 +258 @@
     fQfac->Set(n);
 
-    const Stat_t sumg = fGhness->Integral(1, n+1);
-    const Stat_t sump = fPhness->Integral(1, n+1);
+    Stat_t sumg;
+    Stat_t sump;
+
+    sumg = fGhness->Integral(1, n);
+    sump = fPhness->Integral(1, n);
+
+    if (sumg == 0.0  ||  sump == 0.0)
+    {
+      *fLog << "MHHadronness::Finalize; sumg or sump is zero;   sumg, sump = "
+            << sumg << ",  " << sump << ".  Cannot calculate hadronness" 
+            << endl;
+    }
+
+
+    // Normalize photon distribution
+    Stat_t con;
+    if (sumg > 0.0)
+      for (Int_t i=1; i<=n; i++)
+      {
+        con = (fGhness->GetBinContent(i)) / sumg;
+        fGhness->SetBinContent(i, con);       
+      }
+
+    // Normalize hadron distribution
+    if (sump > 0.0)
+      for (Int_t i=1; i<=n; i++)
+      {
+        con = (fPhness->GetBinContent(i)) / sump;
+        fPhness->SetBinContent(i, con);       
+      }
+
+    // Calculate acceptances
+    sumg = fGhness->Integral(1, n);
+    sump = fPhness->Integral(1, n);
+
+    *fLog << "MHHadronness::Finalize; sumg, sump = " << sumg << ",  "
+          << sump << endl;
 
     Float_t max=0;
@@ -259 +300 @@
     for (Int_t i=1; i<=n; i++)
     {
-        const Stat_t ip = fPhness->Integral(1, i)/sump;
-        const Stat_t ig = fGhness->Integral(1, i)/sumg;
+        Stat_t ip; 
+        if (sump != 0.0)
+          ip = fPhness->Integral(1, i)/sump;
+        else
+          ip = 0;
+
+        Stat_t ig;
+        if (sumg != 0.0)
+          ig = fGhness->Integral(1, i)/sumg;
+        else
+          ig = 0;
 
         fIntPhness->SetBinContent(i, ip);
@@ -409 +459 @@
 
     c.cd(1);
-    gStyle->SetOptStat(10);
+    //gStyle->SetOptStat(10);
     Getghness()->DrawCopy();
     Getphness()->SetLineColor(kRed);
     Getphness()->DrawCopy("same");
     c.cd(2);
-    gStyle->SetOptStat(0);
+    //gStyle->SetOptStat(0);
     Getighness()->DrawCopy();
     Getiphness()->SetLineColor(kRed);
@@ -483 +533 @@
 
     gPad->cd(1);
-    gStyle->SetOptStat(10);
+    //gStyle->SetOptStat(10);
     Getghness()->Draw();
     Getphness()->SetLineColor(kRed);
     Getphness()->Draw("same");
     gPad->cd(2);
-    gStyle->SetOptStat(0);
+    //gStyle->SetOptStat(0);
     Getighness()->Draw();
     Getiphness()->SetLineColor(kRed);

trunk/MagicSoft/Mars/mhist/MHMatrix.cc

@@ -504 +504 @@
 }
 
+// --------------------------------------------------------------------------
+//
 void MHMatrix::Reassign()
 {
@@ -691 +693 @@
 //
 // Define the reference matrix
-//   refcolumn  number of the column containing the variable, for which a
-//              target distribution may be given;
+//   refcolumn  number of the column (starting at 1)containing the variable, 
+//              for which a target distribution may be given;
 //              if refcolumn is negative the target distribution will be set 
 //              equal to the real distribution; the events in the reference
@@ -748 +750 @@
    //---------------------------------------------------------
    //
-   // if refcol < 0 : select reference events randomly
-   //                 i.e. set the normaliztion factotrs equal to 1.0
+   // if refcolumn < 0 : select reference events randomly
+   //                    i.e. set the normaliztion factotrs equal to 1.0
+   // refcol is the column number starting at 0; it is >= 0 
 
    if (refcolumn<0) 
    {
-     frefcol = -refcolumn;
+     frefcol = -refcolumn - 1;
    }
    else
    {
-     frefcol = refcolumn;
+     frefcol =  refcolumn - 1;
    }
    
@@ -790 +793 @@
      //      << ",  " << fM(j-1,frefcol) << endl;
    }
-
-   // if refcolumn<0 set target distribution equal to the real distribution
-   // in order to obtain normalization factors = 1.0
-   //if (refcolumn<0) 
-   //{
-   //  for (Int_t j=1; j<=fnbins; j++) 
-   //  {
-   //    float cont = fHth-> GetBinContent(j);
-   //    fHthsh->SetBinContent(j, cont);
-   //  }
-   //}
 
    //---------------------------------------------------------
@@ -816 +808 @@
 
      // if refcolumn < 0 set the correction factors equal to 1
-     if ( refcolumn>=0.0 )
+     if ( refcolumn>=0 )
        b = fHthsh->GetBinContent(j);
      else
@@ -948 +940 @@
    {
      *fLog <<ir <<" ";
-     for (Int_t ic=0;ic<13;ic++) cout<<Mnew(ir,ic)<<" ";
+     for (Int_t ic=0; ic<Mnew.GetNcols(); ic++) cout<<Mnew(ir,ic)<<" ";
      *fLog <<endl;
    }
@@ -957 +949 @@
    {
      float a = fHthaft->GetBinContent(j);
-     if (a>0) *fLog << j << "  "<< a << "   ";
+     *fLog << j << "  "<< a << "   ";
    }
    *fLog <<endl;
@@ -967 +959 @@
 
    th1->cd(1);
-   ((TH1F*)fHthsh)->Draw();      // target
+   ((TH1F*)fHthsh)->DrawCopy();      // target
 
    th1->cd(2);
-   ((TH1F*)fHth)->Draw();        // real histogram before
+   ((TH1F*)fHth)->DrawCopy();        // real histogram before
 
    th1->cd(3);
-   ((TH1F*)fHthd) -> Draw();     // correction factors
+   ((TH1F*)fHthd)->DrawCopy();       // correction factors
 
    th1->cd(4);
-   ((TH1F*)fHthaft) -> Draw();   // histogram after
-
-   //---------------------------------------------------------
-   // --- write onto output file
-   //
-   //TFile *outfile = new TFile(fileNameout, "RECREATE", "");
-   //Mnew.Write(fileNameout);
-   //outfile->Close();
+   ((TH1F*)fHthaft)->DrawCopy();     // histogram after
+
+   //---------------------------------------------------------
 
    return kTRUE;

trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc

@@ -224 +224 @@
 TObject *MHSigmaTheta::DrawClone(Option_t *opt) 
 {
-    TCanvas &c = *MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta",
+    TCanvas &c = *MakeDefCanvas("SigmaThetaPlot", "Sigmabar vs. Theta",
                                  900, 900);
     c.Divide(3, 3);
@@ -236 +236 @@
     c.cd(1);
     h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E");
+    h->SetDirectory(NULL);
     h->SetTitle("Distribution of \\Theta");
     h->SetXTitle("\\Theta [\\circ]");
     h->SetYTitle("No.of events");
 
-    h->Draw(opt);
+    h->DrawCopy(opt);
     h->SetBit(kCanDelete);;
     gPad->SetLogy();
@@ -246 +247 @@
     c.cd(4);
     h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E");
+    h->SetDirectory(NULL);
     h->SetTitle("Distribution of Sigmabar");
     h->SetXTitle("Sigmabar");
     h->SetYTitle("No.of events");
 
-    h->Draw(opt);
+    h->DrawCopy(opt);
     h->SetBit(kCanDelete);;
 
@@ -263 +265 @@
     c.cd(2);
     l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zx");
+    l->SetDirectory(NULL);
     l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
     l->SetXTitle("\\Theta [\\circ]");
     l->SetYTitle("Sigma^2-Sigmabar^2");
 
-    l->Draw("box");
+    l->DrawCopy("box");
     l->SetBit(kCanDelete);;
 
     c.cd(5);
     l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zy");
+    l->SetDirectory(NULL);
     l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
     l->SetXTitle("pixel");
     l->SetYTitle("Sigma^2-Sigmabar^2");
 
-    l->Draw("box");
+    l->DrawCopy("box");
     l->SetBit(kCanDelete);;
 
@@ -290 +294 @@
     c.cd(3);
     k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zx");
+    k->SetDirectory(NULL);
     k->SetTitle("Sigma vs. \\Theta (all pixels)");
     k->SetXTitle("\\Theta [\\circ]");
     k->SetYTitle("Sigma");
 
-    k->Draw("box");
+    k->DrawCopy("box");
     k->SetBit(kCanDelete);;
 
     c.cd(6);
     k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zy");
+    k->SetDirectory(NULL);
     k->SetTitle("Sigma vs. pixel number (all \\Theta)");
     k->SetXTitle("pixel");
     k->SetYTitle("Sigma");
 
-    k->Draw("box");
+    k->DrawCopy("box");
     k->SetBit(kCanDelete);;
 

trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx

@@ -88 +88 @@
   void Print(Option_t *opt=NULL) const;
 
+  UInt_t GetEvtNumber() const { return fEvtNumber; }  //Get Event Number
   Short_t GetPartId() const { return fPartId; }       //Get Type of particle
   Float_t GetEnergy() const { return fEnergy; }        //Get Energy