Changeset 2301 for trunk

Timestamp:

08/19/03 10:59:28 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/macros/CT1Analysis.C (modified) (27 diffs)

trunk/MagicSoft/Mars/macros/CT1Analysis.C

@@ -1 +1 @@
 
 #include "CT1EgyEst.C"
+
+#include "MBinning.h"
+#include "MBlindPixelCalc.h"
+#include "MContinue.h"
+#include "MCT1PointingCorrCalc.h"
+
+#include "MFCT1SelBasic.h"
+#include "MFCT1SelStandard.h"
+#include "MFCT1SelFinal.h"
+#include "MFillH.h"
+
+#include "MHillasCalc.h"
+#include "MHillasSrcCalc.h"
+#include "MImgCleanStd.h"
+
+#include "MParList.h"
+#include "MSigmabarCalc.h"
+#include "MSrcPosCam.h"
+#include "MTaskList.h"
+#include "MWriteRootFile.h"
+
+
 
 void InitBinnings(MParList *plist)
@@ -6 +28 @@
         gLog << "InitBinnings" << endl;
 
+        //--------------------------------------------
         MBinning *binse = new MBinning("BinningE");
         //binse->SetEdgesLog(30, 1.0e2, 1.0e5);
-        binse->SetEdges(30, 2, 5);
+
+        //This is Daniel's binning in energy:
+        binse->SetEdgesLog(14, 296.296, 86497.6);
         plist->AddToList(binse);
+
+        //--------------------------------------------
 
         MBinning *binssize = new MBinning("BinningSize");
@@ -51 +78 @@
         MBinning *binth = new MBinning("BinningTheta");
         Double_t yedge[9] = 
-                       {7.5, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
+                       {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
         TArrayD yed(9,yedge);
         binth->SetEdges(yed);
@@ -69 +96 @@
         binsdiff->SetEdges(100, -5.0, 20.0);
         plist->AddToList(binsdiff);
+
+        // robert ----------------------------------------------
+        MBinning *binsalphaf = new MBinning("BinningAlphaFlux");
+        binsalphaf->SetEdges(100, -100, 100);
+        plist->AddToList(binsalphaf);
+
+        MBinning *binsdifftime = new MBinning("BinningTimeDiff");
+        binsdifftime->SetEdges(50, 0, 10);
+        plist->AddToList(binsdifftime);
+
+        MBinning *binstime = new MBinning("BinningTime");
+        binstime->SetEdges(50, 44500, 61000);
+        plist->AddToList(binstime);
 }
 
@@ -112 +152 @@
 
         bin = plist->FindObject("BinningDiffsigma2");
+        if (bin) delete bin;
+
+        //robert
+        bin = plist->FindObject("BinningAlphaFlux");
+        if (bin) delete bin;
+
+        bin = plist->FindObject("BinningTimeDiff");
+        if (bin) delete bin;
+
+        bin = plist->FindObject("BinningTime");
         if (bin) delete bin;
 }
@@ -173 +223 @@
     Bool_t RLookNN     = kFALSE;  // read in look-alike events
     Bool_t WNN         = kFALSE;  // update input root file ?
+
+
+    // Job B_SC_UP : read ON1.root (or MC1.root) file 
+    //  - depending on WParSC : create (or read in) supercuts parameter values
+    //  - calculate hadroness for the SUPERCUTS
+    //  - update the input files with the hadroness (ON1.root or MC1.root)
+
+    Bool_t JobB_SC_UP  = kTRUE;
+    Bool_t RMatrix     = kFALSE;  // read matrices from file  
+    Bool_t WParSC      = kTRUE;  // do optimization and write supercuts 
+                                  // parameter values onto a file
+    Bool_t WSC         = kFALSE;  // update input root file ?
 
 
@@ -209 +271 @@
 
 
-
     // Job E_XX :  
     //  - select g/h separation method XX
-    //  - read MC_XX2.root file 
+    //  - read MC root file 
     //  - calculate eff. collection area
-    //  - read ON_XX2.root file 
+    //  - optimize energy estimator
+    //  - read ON root file 
+    //  - apply final cuts
+    //  - write root file for ON data after final cuts 
+
+    Bool_t JobE_XX  = kFALSE;  
+    Bool_t WXX      = kFALSE;  // write out root file ?
+
+
+    // Job F_XX : extended version of E_XX (including flux plots)  
+    //  - select g/h separation method XX
+    //  - read MC root file 
+    //  - calculate eff. collection area
+    //  - optimize energy estimator
+    //  - read ON root file 
     //  - apply final cuts
     //  - calculate flux
-    //  - write root file for ON data after final cuts (ON3.root))
-
-    Bool_t JobE_XX  = kTRUE;  
-    Bool_t WXX      = kFALSE;  // write out root file ON3.root ?
+    //  - write root file for ON data after final cuts 
+
+    Bool_t JobF_XX  = kFALSE;  
+    Bool_t WFX      = kFALSE;  // write out root file  ?
 
 
@@ -344 +419 @@
     contstandard.SetName("SelStandard");
 
-
-      MWriteRootFile write(outNameImage);
+    if (WON1)
+    {
+      MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -357 +433 @@
       write.AddContainer("MHillasSrc",    "Events");
       write.AddContainer("MNewImagePar",  "Events");
-
+    }
 
     //$$$$$$$$$$$$$$$$$$$$$$$$$$$$
@@ -677 +753 @@
 
 
-
-      MWriteRootFile write(outNameImage);
+    if (WMC1)
+    {
+      MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -690 +767 @@
       write.AddContainer("MHillasSrc",    "Events");
       write.AddContainer("MNewImagePar",  "Events");
-
+    }
 
 
@@ -1315 +1392 @@
 
   //---------------------------------------------------------------------
+  // Job B_SC_UP
+  //============
+
+    //  - create (or read in) optimum supercuts parameter values
+    //
+    //  - calculate the hadroness for the supercuts
+    //
+    //  - update input root file, including the hadroness
+
+
+ if (JobB_SC_UP)
+ {
+    gLog << "=====================================================" << endl;
+    gLog << "Macro CT1Analysis : Start of Job B_SC_UP" << endl;
+
+    gLog << "" << endl;
+    gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = " 
+         << JobB_SC_UP  << ",  " << WParSC << ",  " << WSC << endl;
+
+
+
+    //--------------------------------------------
+    // file to be updated (either ON or MC)
+
+    TString typeInput = "ON";
+    //TString typeInput = "MC";
+    gLog << "typeInput = " << typeInput << endl;
+
+    // name of input root file
+    TString filenameData = outPath;
+    filenameData += typeInput;
+    filenameData += "2.root";
+    gLog << "filenameData = " << filenameData << endl; 
+
+    // name of output root file
+    TString outNameImage = outPath;
+    outNameImage += typeInput;
+    outNameImage += "3.root";
+    outNameImage += "xxxxx";
+    
+
+    //TString outNameImage = filenameData;
+
+    gLog << "outNameImage = " << outNameImage << endl; 
+
+    //--------------------------------------------
+    // files to be read for optimizing the supercuts
+    // 
+    // for the training
+    TString filenameTrain = outPath;
+    filenameTrain += "ON";
+    filenameTrain += "2.root";
+    Int_t howManyTrain = 800000;
+    gLog << "filenameTrain = " << filenameTrain << ",   howManyTrain = "
+         << howManyTrain  << endl; 
+
+    // for testing
+    TString filenameTest = outPath;
+    filenameTest += "ON";
+    filenameTest += "2.root";
+    Int_t howManyTest = 100000;
+
+    gLog << "filenameTest = " << filenameTest << ",   howManyTest = "
+         << howManyTest  << endl; 
+    
+
+    //--------------------------------------------
+    // files to contain the matrices (genertated from filenameTrain and
+    //                                                filenameTest)
+    // 
+    // for the training
+    TString fileMatrixTrain = outPath;
+    fileMatrixTrain += "MatrixTrainSC";
+    fileMatrixTrain += ".root";
+    gLog << "fileMatrixTrain = " << fileMatrixTrain << endl; 
+
+    // for testing
+    TString fileMatrixTest = outPath;
+    fileMatrixTest += "MatrixTestSC";
+    fileMatrixTest += ".root";
+    gLog << "fileMatrixTest = " << fileMatrixTest << endl; 
+
+    
+    //--------------------------------------------
+    // file which contains the optimum parameter values for the supercuts 
+
+    TString parSCfile = outPath;
+    parSCfile += "parSC";
+
+    gLog << "parSCfile = " << parSCfile << endl;
+
+
+    //---------------------------------------------------------------------
+    // optimize supercuts using the training sample
+    // and test them using the test sample
+
+if (WParSC)
+  {
+    MCT1FindSupercuts findsuper;
+    findsuper.SetFilenameParam(parSCfile);
+    findsuper.SetHadronnessName("HadSC");
+
+
+    //--------------------------
+    // create matrices 
+    if (!RMatrix)
+    {
+      MH3 &mh3 = *(new MH3("MHillas.fSize"));
+      mh3.SetName("Target distribution for SIZE");
+
+      if (filenameTrain == filenameTest)
+      {
+        findsuper.DefineTrainTestMatrix(filenameTrain, 
+                                        howManyTrain, mh3,
+                                        howManyTest,  mh3);
+      }
+      else
+      {
+        findsuper.DefineTrainMatrix(filenameTrain, howManyTrain, mh3);
+        findsuper.DefineTestMatrix( filenameTest,  howManyTest,  mh3);
+      }
+ 
+      // writing doesn't work ???
+      //findsuper.WriteMatrix(fileMatrixTrain, fileMatrixTest);
+    }
+
+    //--------------------------
+    // read matrices from files
+    //                              NOT YET TESTED
+    //if (RMatrix)
+    //  findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
+    //--------------------------
+
+
+    //--------------------------
+    // optimize the supercuts
+    Bool_t rf = findsuper.FindParams();
+
+    if (!rf) 
+    {
+      *fLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
+       return;
+    }
+
+    //--------------------------------------
+    // test the supercuts on the test sample
+    //    
+    Bool_t rt = findsuper.TestParams();
+  }
+
+
+    //-----------------------------------------------------------------
+    // Update the input files with the SC hadronness
+    //
+
+
+    gLog << "" << endl;
+    gLog << "========================================================" << endl;
+    gLog << "Update input file '" <<  filenameData 
+         << "' with the SC hadronness" << endl;
+
+
+    //----------------------------------------------------
+    // read in optimum parameter values for the supercuts
+
+    TFile inparam(parSCfile);
+    MCT1Supercuts scin;
+    scin.Read("MCT1Supercuts");
+    inparam.Close();
+
+    gLog << "Optimum parameter values for supercuts were read in from file '"
+         << parSCfile << "'" << endl;
+
+    TArrayD supercutsPar(72);
+    scin.GetParams(72, supercutsPar);
+
+    gLog << "Optimum parameter values for supercuts : " << endl;
+    for (Int_t i=0; i<72; i++)
+    {
+      gLog << supercutsPar[i] << ",  ";
+    }
+    gLog << endl;
+
+
+    //----------------------------------------------------
+    MTaskList tliston;
+    MParList pliston;
+
+    // set the parameters of the supercuts
+    MCT1Supercuts supercuts;
+    supercuts.SetParams(72, supercutsPar);
+    gLog << "parameter values for the supercuts used for updating the input file ' " 
+         << filenameData << "'" << endl;
+    supercuts.GetParams(72, supercutsPar);
+    for (Int_t i=0; i<72; i++)
+    {
+      gLog << supercutsPar[i] << ",  ";
+    }
+    gLog << endl;
+
+
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    //-------------------------------------------
+    // create the tasks which should be executed 
+    //
+
+    MReadMarsFile read("Events", filenameData);
+    read.DisableAutoScheme();
+
+    TString fHilName    = "MHillas"; 
+    TString fHilNameExt = "MHillasExt"; 
+    TString fHilNameSrc = "MHillasSrc"; 
+    TString fImgParName = "MNewImagePar"; 
+
+
+    //.......................................................................
+    // calculation of hadroness for the supercuts
+    // (=0.25 if fullfilled, =0.75 otherwise)
+
+    TString hadSCName = "HadSC";
+    MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc);
+    sccalc.SetHadronnessName(hadSCName);
+
+
+    //.......................................................................
+
+
+ if (WSC)
+  {
+      //MWriteRootFile write(outNameImage, "UPDATE");
+      MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
+
+      write.AddContainer("MRawRunHeader", "RunHeaders");
+      write.AddContainer("MTime",         "Events");
+      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("ThetaOrig",     "Events");
+      write.AddContainer("MSrcPosCam",    "Events");
+      write.AddContainer("MSigmabar",     "Events");
+      write.AddContainer("MHillas",       "Events");
+      write.AddContainer("MHillasExt",    "Events");
+      write.AddContainer("MHillasSrc",    "Events");
+      write.AddContainer("MNewImagePar",  "Events");
+
+      write.AddContainer("HadNN",         "Events");
+      write.AddContainer("HadSC",         "Events");
+  }
+
+
+    //-----------------------------------------------------------------
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    Float_t maxhadronness =  0.40;
+    Float_t maxalpha      =  20.0;
+    Float_t maxdist       =  10.0;
+
+    MFCT1SelFinal selfinalgh(fHilNameSrc);
+    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
+    selfinalgh.SetHadronnessName(hadSCName);
+    selfinalgh.SetName("SelFinalgh");
+    MContinue contfinalgh(&selfinalgh);
+    contfinalgh.SetName("ContSelFinalgh");
+
+    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
+    fillhadsc.SetName("HhadSC");
+
+    MFCT1SelFinal selfinal(fHilNameSrc);
+    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
+    selfinal.SetHadronnessName(hadSCName);
+    selfinal.SetName("SelFinal");
+    MContinue contfinal(&selfinal);
+    contfinal.SetName("ContSelFinal");
+
+    TString mh3name = "abs(Alpha)";
+    MBinning binsalphaabs("Binning"+mh3name);
+    binsalphaabs.SetEdges(50, -2.0, 98.0);
+
+    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
+    alphaabs.SetName(mh3name);
+    MFillH alpha(&alphaabs);
+    alpha.SetName("FillAlphaAbs");
+
+    MFillH hfill1("MHHillas",    fHilName);
+    hfill1.SetName("HHillas");
+
+    MFillH hfill2("MHStarMap",   fHilName);
+    hfill2.SetName("HStarMap");
+
+    MFillH hfill3("MHHillasExt",    fHilNameSrc);
+    hfill3.SetName("HHillasExt");
+    
+    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
+    hfill4.SetName("HHillasSrc");    
+
+    MFillH hfill5("MHNewImagePar", fImgParName);
+    hfill5.SetName("HNewImagePar");
+
+    //*****************************
+    // entries in MParList
+
+    pliston.AddToList(&tliston);
+    InitBinnings(&pliston);
+    pliston.AddToList(&binsalphaabs);
+    pliston.AddToList(&alphaabs);
+    pliston.AddToList(&supercuts);
+
+    //*****************************
+    // entries in MTaskList
+    
+    tliston.AddToList(&read);
+
+    tliston.AddToList(&sccalc);
+    tliston.AddToList(&fillhadsc);
+    tliston.AddToList(&contfinalgh);
+
+    tliston.AddToList(&alpha);
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+    tliston.AddToList(&hfill5);
+
+    if (WSC)
+      tliston.AddToList(&write);
+
+    tliston.AddToList(&contfinal);
+
+    //*****************************
+
+    //-------------------------------------------
+    // Execute event loop
+    //
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxevents = -1;
+    if ( !evtloop.Eventloop(maxevents) )
+        return;
+
+    tliston.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
+
+    pliston.FindObject("MHHillas")->DrawClone();
+    pliston.FindObject("MHHillasExt")->DrawClone();
+    pliston.FindObject("MHHillasSrc")->DrawClone();
+    pliston.FindObject("MHNewImagePar")->DrawClone();
+    pliston.FindObject("MHStarMap")->DrawClone();
+
+     //-------------------------------------------
+    // fit alpha distribution to get the number of excess events and
+    // calculate significance of gamma signal in the alpha plot
+  
+    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
+    TH1  &alphaHist = absalpha->GetHist();
+    alphaHist.SetXTitle("|alpha|  [\\circ]");
+
+    Double_t alphasig = 13.1;
+    Double_t alphamin = 30.0;
+    Double_t alphamax = 90.0;
+    Int_t    degree   =    4;
+    Double_t significance = -99.0;
+    Bool_t   drawpoly  = kTRUE;
+    Bool_t   fitgauss  = kTRUE;
+    Bool_t   print     = kTRUE;
+
+    MHFindSignificance findsig;
+    findsig.SetRebin(kTRUE);
+    findsig.SetReduceDegree(kFALSE);
+
+    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
+                        alphasig, drawpoly, fitgauss, print);
+    significance = findsig.GetSignificance();
+    Float_t alphasi = findsig.GetAlphasi();
+
+    gLog << "For file '" << filenameData << "' : " << endl;
+    gLog << "Significance of gamma signal after supercuts : "
+         << significance << " (for |alpha| < " << alphasi << " degrees)" 
+         << endl;
+
+    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
+
+    //-------------------------------------------
+
+    DeleteBinnings(&pliston);
+
+    gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl;
+    gLog << "=======================================================" << endl;
+ }
+  //---------------------------------------------------------------------
+
+
+  //---------------------------------------------------------------------
   // Job B_RF_UP
   //============
@@ -2336 +2816 @@
   
     MHOnSubtraction onsub;
-    onsub.Calc(alphaHist, &pliston, kTRUE, 13.1);
+    onsub.Calc(alphaHist, kTRUE, 13.1);
     //-------------------------------------------
 
@@ -2387 +2867 @@
 
     //TString XX("NN");
-    //TString XX("SC");
-    TString XX("RF");
+    TString XX("SC");
+    //TString XX("RF");
     TString fhadronnessName("Had");
     fhadronnessName += XX;
@@ -2412 +2892 @@
     TString energyParName(outPath);
     energyParName += "energyest_";
-    //energyParName += XX;
+    energyParName += XX;
     energyParName += ".root";
     gLog << "energyParName = " << energyParName << endl;
@@ -2420 +2900 @@
     TString filenameArea(outPath);
     filenameArea += typeMC;
-    //filenameArea += "_";
-    //filenameArea += XX;
     filenameArea += ext; 
     gLog << "filenameArea = " << filenameArea << endl;
@@ -2429 +2907 @@
     TString collareaName(outPath);
     collareaName += "area_";
-    //collareaName += XX;
+    collareaName += XX;
     collareaName += ".root";
     gLog << "collareaName = " << collareaName << endl;
@@ -2437 +2915 @@
     TString filenameData(outPath);
     filenameData += typeData;
-    //filenameData += "_";
-    //filenameData += XX;
     filenameData += ext;
     gLog << "filenameData = " << filenameData << endl;
@@ -2446 +2922 @@
     TString filenameDataout(outPath);
     filenameDataout += typeData;
-    //filenameDataout += "_";
-    //filenameDataout += XX;
+    filenameDataout += "_";
+    filenameDataout += XX;
     filenameDataout += extout;
     gLog << "filenameDataout = " << filenameDataout << endl;
@@ -2538 +3014 @@
     f.Close();
 
+    gLog << "Collection area plots written onto file " << collareaName << endl;
 
     gLog << "Calculation of effective collection areas done" << endl;
@@ -2554 +3031 @@
     // Optimization of energy estimator
     //
+    gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl;
 
     TString inpath("");
@@ -2562 +3040 @@
             howMany,  maxhadronness, maxalpha, maxdist,
             binsE, binsTheta);
+    gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
 
     //-----------------------------------------------------------
@@ -2622 +3101 @@
     //.......................................................................
 
-
-      MWriteRootFile write(filenameDataout);
+    if (WXX)
+    {
+      MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -2641 +3121 @@
 
       write.AddContainer("MEnergyEst",    "Events");
-
+    }
 
     //-----------------------------------------------------------------
@@ -2670 +3150 @@
     eeston.SetCoeffA(parA);
     eeston.SetCoeffB(parB);
+
+    //---------------------------
+    // calculate estimated energy using Daniel's parameters
+
+    //MEnergyEstParam eeston(fHilName);
+    //eeston.Add(fHilNameSrc);
+
     //---------------------------
 
@@ -2766 +3253 @@
   //---------------------------------------------------------------------
 
+
+  //---------------------------------------------------------------------
+  // Job F_XX
+  //=========
+
+    //  - select g/h separation method XX
+    //  - read MC_XX2.root file 
+    //  - calculate eff. collection area
+    //  - read ON_XX2.root file 
+    //  - apply final cuts
+    //  - calculate flux
+    //  - write root file for ON data after final cuts (ON_XX3.root))
+
+
+ if (JobF_XX)
+ {
+    gLog << "=====================================================" << endl;
+    gLog << "Macro CT1Analysis : Start of Job F_XX" << endl;
+
+    gLog << "" << endl;
+    gLog << "Macro CT1Analysis : JobF_XX, WXX = " 
+         << JobF_XX  << ",  " << WFX << endl;
+
+    // type of data to be analysed
+    TString typeData = "ON";
+    //TString typeData = "OFF";
+    //TString typeData = "MC";
+    gLog << "typeData = " << typeData << endl;
+
+    TString typeMC   = "MC";
+    TString ext      = "3.root";
+    TString extout   = "4.root";
+
+    //------------------------------
+    // selection of g/h separation method
+    // and definition of final selections
+
+    //TString XX("NN");
+    TString XX("SC");
+    //TString XX("RF");
+    TString fhadronnessName("Had");
+    fhadronnessName += XX;
+    gLog << "fhadronnessName = " << fhadronnessName << endl;
+
+    // maximum values of the hadronness, |ALPHA| and DIST
+    Float_t maxhadronness   = 0.40;
+    Float_t maxalpha        = 20.0;
+    Float_t maxdist         = 10.0;
+    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
+         << maxhadronness << ",  " << maxalpha << ",  " 
+         << maxdist << endl;
+
+    //------------------------------
+    // name of MC file to be used for optimizing the energy estimator
+    TString filenameOpt(outPath);
+    filenameOpt += typeMC;
+    filenameOpt += ext; 
+    gLog << "filenameOpt = " << filenameOpt << endl;
+
+    //------------------------------
+    // name of file containing the parameters of the energy estimator
+    TString energyParName(outPath);
+    energyParName += "energyest_";
+    energyParName += XX;
+    energyParName += ".root";
+    gLog << "energyParName = " << energyParName << endl;
+
+    //------------------------------
+    // name of MC file to be used for calculating the eff. collection areas
+    TString filenameArea(outPath);
+    filenameArea += typeMC;
+    filenameArea += ext; 
+    gLog << "filenameArea = " << filenameArea << endl;
+
+    //------------------------------
+    // name of file containing the eff. collection areas
+    TString collareaName(outPath);
+    collareaName += "area_";
+    collareaName += XX;
+    collareaName += ".root";
+    gLog << "collareaName = " << collareaName << endl;
+
+    //------------------------------
+    // name of data file to be analysed
+    TString filenameData(outPath);
+    filenameData += typeData;
+    filenameData += ext;
+    gLog << "filenameData = " << filenameData << endl;
+
+    //------------------------------
+    // name of output data file (after the final cuts)
+    TString filenameDataout(outPath);
+    filenameDataout += typeData;
+    filenameDataout += "_";
+    filenameDataout += XX;
+    filenameDataout += extout;
+    gLog << "filenameDataout = " << filenameDataout << endl;
+
+    //------------------------------
+    // name of file containing histograms for flux calculastion
+    TString filenameResults(outPath);
+    filenameResults += typeData;
+    filenameResults += "Results_";
+    filenameResults += XX;
+    filenameResults += extout;
+    gLog << "filenameResults = " << filenameResults << endl;
+
+
+    //====================================================================
+    gLog << "-----------------------------------------------" << endl;
+    gLog << "Start calculation of effective collection areas" << endl;
+    MParList  parlist;
+    MTaskList tasklist;
+
+    //---------------------------------------
+    // Setup the tasks to be executed
+    //
+    MReadMarsFile reader("Events", filenameArea);
+    reader.DisableAutoScheme();
+
+    MFCT1SelFinal cuthadrons;
+    cuthadrons.SetHadronnessName(fhadronnessName);
+    cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
+
+    MContinue conthadrons(&cuthadrons);
+
+    MHMcCT1CollectionArea collarea;
+    collarea.SetEaxis(MHMcCT1CollectionArea::kLinear);
+
+    MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
+    filler.SetName("CollectionArea");
+
+    //********************************
+    // entries in MParList
+
+    parlist.AddToList(&tasklist);
+    InitBinnings(&parlist);
+    parlist.AddToList(&collarea);
+
+    //********************************
+    // entries in MTaskList
+
+    tasklist.AddToList(&reader);   
+    tasklist.AddToList(&conthadrons);
+    tasklist.AddToList(&filler);
+
+    //********************************
+
+    //-----------------------------------------
+    // Execute event loop
+    //
+    MEvtLoop magic;
+    magic.SetParList(&parlist);
+
+    MProgressBar bar;
+    magic.SetProgressBar(&bar);
+    if (!magic.Eventloop())
+        return;
+
+    tasklist.PrintStatistics(0, kTRUE);
+
+    // Calculate effective collection areas 
+    // and display the histograms
+    //
+    //MHMcCT1CollectionArea *collarea = 
+    //     (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
+    collarea.CalcEfficiency();
+    collarea.DrawClone();
+
+    // save binnings for call to CT1EEst
+    MBinning *binsE     = (MBinning*)parlist.FindObject("BinningE");
+    if (!binsE)
+        {
+          gLog << "Object 'BinningE' not found in MParList" << endl;
+          return;
+        }
+    MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
+    if (!binsTheta)
+        {
+          gLog << "Object 'BinningTheta' not found in MParList" << endl;
+          return;
+        }
+
+
+    //---------------------------------------------
+    // Write histograms to a file 
+    //
+
+    TFile f(collareaName, "RECREATE");
+    collarea.GetHist()->Write();
+    collarea.GetHAll()->Write();
+    collarea.GetHSel()->Write();
+    f.Close();
+
+    gLog << "Collection area plots written onto file " << collareaName << endl;
+
+    gLog << "Calculation of effective collection areas done" << endl;
+    gLog << "-----------------------------------------------" << endl;    
+    //------------------------------------------------------------------
+
+
+    TString fHilName    = "MHillas"; 
+    TString fHilNameExt = "MHillasExt"; 
+    TString fHilNameSrc = "MHillasSrc"; 
+    TString fImgParName = "MNewImagePar"; 
+
+  
+    //===========================================================
+    //
+    // Optimization of energy estimator
+    //
+    gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl;
+
+    TString inpath("");
+    TString outpath("");
+    Int_t howMany = 2000;
+    CT1EEst(inpath,   filenameOpt,   outpath, energyParName, 
+            fHilName, fHilNameSrc,   fhadronnessName,
+            howMany,  maxhadronness, maxalpha, maxdist,
+            binsE, binsTheta);
+    gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
+
+    //-----------------------------------------------------------
+    //
+    // Read in parameters of energy estimator ("MMcEnergyEst")
+    //                   and migration matrix ("MHMcEnergyMigration")
+    //
+    gLog << "================================================================"
+         << endl;
+    gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
+         << energyParName << "'" << endl;
+    TFile enparam(energyParName);
+    enparam.ls();
+
+    //MMcEnergyEst mcest("MMcEnergyEst"); 
+    //mcest.Read("MMcEnergyEst");
+    MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
+
+    gLog << "Parameters of energy estimator were read in" << endl;
+
+    TArrayD parA(mcest.GetNumCoeffA());
+    TArrayD parB(mcest.GetNumCoeffB());
+    for (Int_t i=0; i<parA.GetSize(); i++)
+      parA[i] = mcest.GetCoeff(i);
+    for (Int_t i=0; i<parB.GetSize(); i++)
+      parB[i] = mcest.GetCoeff( i+parA.GetSize() );
+
+    gLog << "Read in Migration matrix" << endl;   
+
+    //MHMcEnergyMigration mighiston("MHMcEnergyMigration");
+    //mighiston.Read("MHMcEnergyMigration");
+    MHMcEnergyMigration &mighiston = 
+          *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
+
+    gLog << "Migration matrix was read in" << endl;
+
+    //*************************************************************************
+    //
+    // Analyse the data
+    //
+    gLog << "============================================================"
+         << endl;
+    gLog << "Analyse the data" << endl;
+
+    MTaskList tliston;
+    MParList pliston;
+
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+
+    //-------------------------------------------
+    // create the tasks which should be executed 
+    //
+
+    MReadMarsFile read("Events", filenameData);
+    read.DisableAutoScheme();
+
+    //.......................................................................
+
+    if (WFX)
+    {
+      MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+
+      write.AddContainer("MRawRunHeader", "RunHeaders");
+      write.AddContainer("MTime",         "Events");
+      write.AddContainer("MMcEvt",        "Events");
+      write.AddContainer("ThetaOrig",     "Events");
+      write.AddContainer("MSrcPosCam",    "Events");
+      write.AddContainer("MSigmabar",     "Events");
+      write.AddContainer("MHillas",       "Events");
+      write.AddContainer("MHillasExt",    "Events");
+      write.AddContainer("MHillasSrc",    "Events");
+      write.AddContainer("MNewImagePar",  "Events");
+
+      write.AddContainer("HadNN",         "Events");
+      write.AddContainer("HadSC",         "Events");
+      write.AddContainer("HadRF",         "Events");
+
+      write.AddContainer("MEnergyEst",    "Events");
+    }
+
+    //-----------------------------------------------------------------
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    MFCT1SelFinal selfinalgh(fHilNameSrc);
+    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
+    selfinalgh.SetHadronnessName(fhadronnessName);
+    selfinalgh.SetName("SelFinalgh");
+    MContinue contfinalgh(&selfinalgh);
+    contfinalgh.SetName("ContSelFinalgh");
+
+    MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
+    fillhadnn.SetName("HhadNN");
+    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
+    fillhadsc.SetName("HhadSC");
+    MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
+    fillhadrf.SetName("HhadRF");
+
+    //---------------------------
+    // calculate estimated energy
+
+    //MEnergyEstParam eeston(fHilName);
+    //eeston.Add(fHilNameSrc);
+
+    //eeston.SetCoeffA(parA);
+    //eeston.SetCoeffB(parB);
+
+    //---------------------------
+    // calculate estimated energy using Daniel's parameters
+
+    MEnergyEstParamDanielMkn421 eeston(fHilName);
+    eeston.Add(fHilNameSrc);
+    //eeston.SetCoeffA(parA);
+    //eeston.SetCoeffB(parB);
+
+
+    //---------------------------
+
+
+    MFillH hfill1("MHHillas",    fHilName);
+    hfill1.SetName("HHillas");
+
+    MFillH hfill2("MHStarMap",   fHilName);
+    hfill2.SetName("HStarMap");
+
+    MFillH hfill3("MHHillasExt",   fHilNameSrc);
+    hfill3.SetName("HHillasExt");    
+
+    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
+    hfill4.SetName("HHillasSrc");    
+
+    MFillH hfill5("MHNewImagePar", fImgParName);
+    hfill5.SetName("HNewImagePar");    
+
+    //---------------------------
+    // new from Robert
+
+    MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
+    hfill6.SetName("HTimeDiffTheta");
+
+    MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
+    hfill6a.SetName("HTimeDiffTime");
+
+    MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
+    hfill7.SetName("HAlphaEnergyTheta");
+
+    MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
+    hfill7a.SetName("HAlphaEnergyTime");
+
+    MFillH hfill7b("MHThetabarTime", fHilNameSrc);
+    hfill7b.SetName("HThetabarTime");
+
+    MFillH hfill7c("MHEnergyTime", "MMcEvt");
+    hfill7c.SetName("HEnergyTime");
+
+
+    //---------------------------
+
+    MFCT1SelFinal selfinal(fHilNameSrc);
+    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
+    selfinal.SetHadronnessName(fhadronnessName);
+    selfinal.SetName("SelFinal");
+    MContinue contfinal(&selfinal);
+    contfinal.SetName("ContSelFinal");
+
+
+    //*****************************
+    // entries in MParList
+
+    pliston.AddToList(&tliston);
+    InitBinnings(&pliston);
+
+
+    //*****************************
+    // entries in MTaskList
+    
+    tliston.AddToList(&read);
+
+    // robert      
+    tliston.AddToList(&hfill6);   //timediff
+    tliston.AddToList(&hfill6a);   //timediff
+
+    tliston.AddToList(&contfinalgh);
+    tliston.AddToList(&eeston);
+
+    if (WFX)
+      tliston.AddToList(&write);
+
+    tliston.AddToList(&fillhadnn);
+    tliston.AddToList(&fillhadsc);
+    tliston.AddToList(&fillhadrf);
+
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+    tliston.AddToList(&hfill5);
+
+    //robert
+    tliston.AddToList(&hfill7);
+    tliston.AddToList(&hfill7a);
+    tliston.AddToList(&hfill7b);
+    tliston.AddToList(&hfill7c);
+
+    tliston.AddToList(&contfinal);
+
+    //*****************************
+
+    //-------------------------------------------
+    // Execute event loop
+    //
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxevents = -1;
+    //Int_t maxevents = 10000;
+    if ( !evtloop.Eventloop(maxevents) )
+        return;
+
+    tliston.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+
+    pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
+    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
+    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
+
+    pliston.FindObject("MHHillas")->DrawClone();
+    pliston.FindObject("MHHillasExt")->DrawClone();
+    pliston.FindObject("MHHillasSrc")->DrawClone();
+    pliston.FindObject("MHNewImagePar")->DrawClone();
+    pliston.FindObject("MHStarMap")->DrawClone();
+
+    DeleteBinnings(&pliston);
+
+
+//rwagner write all relevant histograms onto a file
+
+    gLog << "=======================================================" << endl;
+    gLog << "Write results onto file '" << filenameResults << "'" << endl;
+
+    TFile outfile(filenameResults,"recreate");
+
+    MHHillasSrc* hillasSrc = 
+      (MHHillasSrc*)(pliston->FindObject("MHHillasSrc"));
+        TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
+    alphaHist->Write();
+    
+    MHAlphaEnergyTheta* aetH = 
+      (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
+    TH3D* aetHist = (TH3D*)(aetH->GetHist());
+    aetHist->SetName("aetHist");
+    aetHist->Write();
+
+    MHAlphaEnergyTime* aetH2 = 
+      (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime"));
+    TH3D* aetHist2 = (TH3D*)(aetH2->GetHist());
+    aetHist2->SetName("aetimeHist");
+//     aetHist2->DrawClone();
+    aetHist2->Write();
+
+    MHThetabarTime* tbt = 
+      (MHThetabarTime*)(pliston->FindObject("MHThetabarTime"));
+    TProfile* tbtHist = (TProfile*)(tbt->GetHist());
+    tbtHist->SetName("tbtHist");
+    tbtHist->Write();
+
+    MHEnergyTime* ent = 
+      (MHEnergyTime*)(pliston->FindObject("MHEnergyTime"));
+    TH2D* entHist = (TH2D*)(ent->GetHist());
+    entHist->SetName("entHist");
+    entHist->Write();
+    
+    MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
+    TH2D* timeHist = (TH2D*)(time->GetHist());
+    timeHist->SetName("MHTimeDiffTheta");
+    timeHist->SetTitle("Time diffs");
+    timeHist->Write();
+
+    MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
+    TH2D* timeHist2 = (TH2D*)(time2->GetHist());
+    timeHist2->SetName("MHTimeDiffTime");
+    timeHist2->SetTitle("Time diffs");
+    timeHist2->Write();
+
+//rwagner write also collareas to same file
+    collarea->GetHist()->Write();
+    collarea->GetHAll()->Write();
+    collarea->GetHSel()->Write();
+    
+//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data
+//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file
+
+    outfile.Close();
+
+    gLog << "=======================================================" << endl;
+
+
+    gLog << "Macro CT1Analysis : End of Job F_XX" << endl;
+    gLog << "=======================================================" << endl;
+ }
+  //---------------------------------------------------------------------
+
 }
 
@@ -2781 +3800 @@
 
 
+
+
+
+
+
+
+
+
+
+
+
+
+