Changeset 2437 for trunk/MagicSoft


Ignore:
Timestamp:
10/28/03 07:33:49 (21 years ago)
Author:
wittek
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r2435 r2437  
    11                                                 -*-*- END OF LINE -*-*-
     2
     3  2003/10/28: Wolfgang Wittek
     4
     5  * manalysis/MCT1PadONOFF.cc
     6    - replace GetMeanRms() by GetPedestalRms()
     7    - replace SetMeanRms() by SetPedestalRms()
     8    - reactivate code which was commented out by tgb
     9                     (no compilation errors on Alpha OSF)
     10
     11  * manalysis/AnalysisLinkDef.h
     12             /Makefile
     13    - put back MCT1PadONOFF
     14
     15  * macros/CT1Analysis.C
     16          /ONOFFCT1Analysis.C
     17    - current versions of macros for the analysis of CT1 data
     18
     19
    220
    321  2003/10/26: Oscar Blanch Bigas
    … …  
    5169     - removed FillLevels - obsolete
    5270     - added SetLevels instead
    53 
    5471
    5572

  • trunk/MagicSoft/Mars/macros/CT1Analysis.C

    r2436 r2437  
    10811081    fillmatg.SetName("fillGammas");
    10821082
    1083     MFillH fillmath("MatrixHadrons");
    1084     fillmath.SetFilter(&flisth);
    1085     fillmath.SetName("fillHadrons");
    10861083
    10871084
    … …  
    11421139    selectorh.SetName("selectHadrons");
    11431140
     1141    MFillH fillmath("MatrixHadrons");
     1142    fillmath.SetFilter(&flisth);
     1143    fillmath.SetName("fillHadrons");
     1144
    11441145
    11451146    // entries in MFilterList
    … …  
    11821183
    11831184
    1184     // write out training events
     1185    // write out matrices of training events
    11851186
    11861187    gLog << "" << endl;
    11871188    gLog << "========================================================" << endl;
    1188     gLog << "Write out training events" << endl;
     1189    gLog << "Write out matrices of training events" << endl;
    11891190
    11901191
    … …  
    12421243      //
    12431244      MEvtLoop evtlooptr;
     1245      evtlooptr.SetName("ReadRFTrees");
    12441246      evtlooptr.SetParList(&plisttr);
    12451247      if (!evtlooptr.Eventloop())
    … …  
    13001302    //
    13011303    MEvtLoop treeloop;
     1304    treeloop.SetName("GrowRFTrees"); 
    13021305    treeloop.SetParList(&plist2);
    13031306

  • trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C

    r2255 r2437  
    11
    22#include "CT1EgyEst.C"
     3
    34
    45void InitBinnings(MParList *plist)
    … …  
    142143
    143144      // path for input for Mars
    144       TString inPath = "~magican/ct1test/wittek/marsoutput/optionC/";
     145      TString inPath = "~magican/ct1test/wittek/marsoutput/optionD/";
    145146
    146147      // path for output from Mars
    147       TString outPath = "~magican/ct1test/wittek/marsoutput/optionC/";
     148      TString outPath = "~magican/ct1test/wittek/marsoutput/optionD/";
    148149
    149150      //-----------------------------------------------
    150151
    151       TEnv env("macros/CT1env.rc");
    152       Bool_t printEnv = kFALSE;
     152      //TEnv env("macros/CT1env.rc");
     153      //Bool_t printEnv = kFALSE;
    153154
    154155    //************************************************************************
    … …  
    167168                               // (or OFF1.root or MC1.root)?
    168169
    169     // Job B_NN_UP : read ON1.root (OFF1.root or MC1.root) file
    170     //  - depending on RlookNN : create (or read in) hadron and gamma matrix
    171     //  - calculate hadroness for method of NEAREST NEIGHBORS
    172     //  - update the input files with the hadroness (ON1.root, OFF1.root
    173     //     or MC1.root)
    174 
    175     Bool_t JobB_NN_UP  = kFALSE; 
    176     Bool_t RLookNN     = kFALSE;  // read in look-alike events
    177     Bool_t WNN         = kFALSE;  // update input root file ?
    178 
    179170
    180171    // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file
    181     //  - depending on RLook : create (or read in) hadron and gamma matrix
    182     //  - depending on RTree : create (or read in) trees
     172    //  - if RTrainRF = TRUE : read in training matrices for hadrons and gammas
     173    //  - if CTrainRF = TRUE : create matrices of training events
     174    //  - if RTree    = TRUE : read in trees, otherwise create trees
    183175    //  - calculate hadroness for method of RANDOM FOREST
    184     //    and for the SUPERCUTS
    185     //  - update the input files with the hadronesses (ON1.root, OFF1.root
    186     //     or MC1.root)
     176    //  - update the input files with the hadronesses (ON2.root, OFF2.root
     177    //     or MC2.root)
    187178
    188179    Bool_t JobB_RF_UP  = kFALSE; 
    189     Bool_t RLookRF     = kFALSE;  // read in look-alike events
     180    Bool_t RTrainRF     = kFALSE;  // read in matrices of training events
     181    Bool_t CTrainRF     = kFALSE;  // create  matrices of training events
    190182    Bool_t RTree       = kFALSE;  // read in trees
    191183    Bool_t WRF         = kFALSE;  // update input root file ?
    … …  
    265257    gLog << "" << endl;
    266258    gLog << "Macro CT1Analysis : JobA, WPad, RPad, Wout = "
    267          << JobA  << ",  " << WPad << ",  " << RPad << ",  " << Wout << endl;
    268 
     259         << (JobA ? "kTRUE" : "kFALSE")  << ",  "
     260         << (WPad ? "kTRUE" : "kFALSE")  << ",  "
     261         << (RPad ? "kTRUE" : "kFALSE")  << ",  "
     262         << (Wout ? "kTRUE" : "kFALSE")  << endl;
     263   
    269264
    270265    //--------------------------------------------------
    … …  
    334329      MCT1ReadPreProc readON(fileON);
    335330
    336       MFCT1SelBasic selthetaon;
    337       selthetaon.SetCuts(-100.0, 29.5, 35.5);
    338       MContinue contthetaon(&selthetaon);
     331      //MFCT1SelBasic selthetaon;
     332      //selthetaon.SetCuts(-100.0, 29.5, 35.5);
     333      //MContinue contthetaon(&selthetaon);
    339334
    340335      MBlindPixelCalc blindon;
    … …  
    346341      MHBlindPixels blindON("BlindPixelsON");
    347342      MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels");
     343      fillblindON.SetName("FillBlind");
    348344
    349345      MSigmabarCalc sigbarcalcon;
    … …  
    351347      MHSigmaTheta sigthON("SigmaThetaON");
    352348      MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt");
    353      
     349      fillsigthetaON.SetName("FillSigTheta");   
     350 
    354351      //*****************************
    355352      // entries in MParList
    … …  
    380377
    381378      Int_t maxeventson = -1;
    382       //Int_t maxeventson = 20000;
     379      //Int_t maxeventson = 10000;
    383380      if ( !evtloopon.Eventloop(maxeventson) )
    384381          return;
    … …  
    420417
    421418      MHBlindPixels blindOFF("BlindPixelsOFF");
    422       MHBlindPixels *pblindOFF = &blindOFF;
    423       gLog << "pblindOFF = " << pblindOFF << endl;
    424 
    425419      MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels");
     420      fillblindOFF.SetName("FillBlindOFF");
    426421
    427422      MSigmabarCalc sigbarcalcoff;
    … …  
    429424      MHSigmaTheta sigthOFF("SigmaThetaOFF");
    430425      MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt");
    431      
     426      fillsigthetaOFF.SetName("FillSigThetaOFF");     
     427
    432428      //*****************************
    433429      // entries in MParList
    … …  
    657653    MEvtLoop evtloop;
    658654    evtloop.SetParList(&pliston);
    659     evtloop.ReadEnv(env, "", printEnv);
     655    //evtloop.ReadEnv(env, "", printEnv);
    660656    evtloop.SetProgressBar(&bar);
    661657    //  evtloop.Write();
    … …  
    694690
    695691
     692
     693
    696694  //---------------------------------------------------------------------
    697   // Job B_NN_UP
     695  // Job B_RF_UP
    698696  //============
    699697
    700     //  - read in the matrices of look-alike events for gammas and hadrons
    701 
    702     // then read ON1.root (or MC1.root) file
    703     //
    704     //  - calculate the hadroness for the method of nearest neighbors
    705     //
    706     //  - update input root file, including the hadroness
    707 
    708 
    709  if (JobB_NN_UP)
     698
     699    //  - create (or read in) the matrices of training events for gammas
     700    //    and hadrons
     701    //  - create (or read in) the trees
     702    //  - then read ON1.root (or MC1.root) file
     703    //  - calculate the hadroness for the method of RANDOM FOREST
     704    //  - update input root file with the hadroness
     705
     706
     707 if (JobB_RF_UP)
    710708 {
    711709    gLog << "=====================================================" << endl;
    712     gLog << "Macro CT1Analysis : Start of Job B_NN_UP" << endl;
     710    gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
    713711
    714712    gLog << "" << endl;
    715     gLog << "Macro CT1Analysis : JobB_NN_UP, RLookNN, WNN = "
    716          << JobB_NN_UP  << ",  " << RLookNN << ",  " << WNN << endl;
    717 
     713    gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = "
     714         << (JobB_RF_UP ? "kTRUE" : "kFALSE")  << ",  "
     715         << (RTrainRF ?   "kTRUE" : "kFALSE")  << ",  "
     716         << (CTrainRF ?   "kTRUE" : "kFALSE")  << ",  "
     717         << (RTree ?      "kTRUE" : "kFALSE")  << ",  "
     718         << (WRF ?        "kTRUE" : "kFALSE")  << endl;
     719
     720
     721    //--------------------------------------------
     722    // parameters for the random forest
     723    Int_t NumTrees = 100;
     724    Int_t NumTry   =   3;
     725    Int_t NdSize   =   1;
     726
     727
     728    TString hadRFName = "HadRF";
     729    Float_t maxhadronness =  0.23;
     730    Float_t maxalpha      =  20.0;
     731    Float_t maxdist       =  10.0;
     732
     733    TString fHilName    = "MHillas";
     734    TString fHilNameExt = "MHillasExt";
     735    TString fHilNameSrc = "MHillasSrc";
     736    TString fImgParName = "MNewImagePar";
    718737
    719738
    … …  
    721740    // file to be updated (ON, OFF or MC)
    722741
    723     //TString typeInput = "ON";
     742    TString typeInput = "ON";
    724743    //TString typeInput = "OFF";
    725     TString typeInput = "MC";
     744    //TString typeInput = "MC";
    726745    gLog << "typeInput = " << typeInput << endl;
    727746
    … …  
    736755    outNameImage += typeInput;
    737756    outNameImage += "2.root";
    738 
    739757    //TString outNameImage = filenameData;
    740758
    … …  
    742760
    743761    //--------------------------------------------
    744     // files to be read for generating the look-alike events
     762    // files to be read for generating the matrices of training events
    745763    // "hadrons" :
    746764    TString filenameHad = outPath;
    747765    filenameHad += "OFF";
    748766    filenameHad += "1.root";
    749     Int_t howManyHadrons = 500000;
     767    Int_t howManyHadrons = 1000000;
    750768    gLog << "filenameHad = " << filenameHad << ",   howManyHadrons = "
    751769         << howManyHadrons  << endl;
    … …  
    761779   
    762780    //--------------------------------------------
    763     // files of look-alike events
     781    // files of training events
    764782
    765783    TString outNameGammas = outPath;
    766     outNameGammas += "matrix_gammas_";
     784    outNameGammas += "RFmatrix_gammas_";
    767785    outNameGammas += "MC";
    768786    outNameGammas += ".root";
    … …  
    772790
    773791    TString outNameHadrons = outPath;
    774     outNameHadrons += "matrix_hadrons_";
     792    outNameHadrons += "RFmatrix_hadrons_";
    775793    outNameHadrons += typeMatrixHadrons;
    776794    outNameHadrons += ".root";
    … …  
    778796
    779797    MHMatrix matrixg("MatrixGammas");
     798    matrixg.EnableGraphicalOutput();
     799
     800    matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)");
     801    matrixg.AddColumn("MSigmabar.fSigmabar");
     802    matrixg.AddColumn("log10(MHillas.fSize)");
     803    matrixg.AddColumn("MHillasSrc.fDist");
     804    matrixg.AddColumn("MHillas.fWidth");
     805    matrixg.AddColumn("MHillas.fLength");
     806    matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
     807    matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
     808    matrixg.AddColumn("MNewImagePar.fConc");
     809    matrixg.AddColumn("MNewImagePar.fLeakage1");
     810
    780811    MHMatrix matrixh("MatrixHadrons");
     812    matrixh.EnableGraphicalOutput();
     813
     814    matrixh.AddColumns(matrixg.GetColumns());
     815
     816    //--------------------------------------------
     817    // file of trees of the random forest
     818
     819    TString outRF = outPath;
     820    outRF += "RF.root";
     821
    781822
    782823   //*************************************************************************
    783    // read in matrices of look-alike events
    784 if (RLookNN)
     824   // read in matrices of training events
     825if (RTrainRF)
    785826  {
    786827    const char* mtxName = "MatrixGammas";
    … …  
    799840
    800841    matrixg.Read(mtxName);
    801     matrixg.Print("cols");
     842    matrixg.Print("SizeCols");
    802843
    803844
    … …  
    819860
    820861    matrixh.Read(mtxName);
    821     matrixh.Print("cols");
     862    matrixh.Print("SizeCols");
    822863  }
    823864
    824865
    825866   //*************************************************************************
    826    // create matrices of look-alike events
    827 if (!RLookNN)
     867   // create matrices of training events
     868if (CTrainRF)
    828869  {
    829870    gLog << "" << endl;
    830871    gLog << "========================================================" << endl;
    831     gLog << " Create matrices of look-alike events" << endl;
     872    gLog << " Create matrices of training events" << endl;
    832873    gLog << " Gammas :" << endl;
    833874
    … …  
    852893    fhadrons.SetName("hadronID)");
    853894
    854     MFEventSelector selectorg;
    855     selectorg.SetNumSelectEvts(howManyGammas);
     895    TString mgname("costhg");
     896    MBinning bing("Binning"+mgname);
     897    bing.SetEdges(10, 0., 1.0);
     898
     899    MH3 gref("cos(MMcEvt.fTelescopeTheta)");
     900    gref.SetName(mgname);
     901    MH::SetBinning(&gref.GetHist(), &bing);
     902    for (Int_t i=1; i<=gref.GetNbins(); i++)
     903      gref.GetHist().SetBinContent(i, 1.0);
     904
     905    MFEventSelector2 selectorg(gref);
     906    selectorg.SetNumMax(howManyGammas);
    856907    selectorg.SetName("selectGammas");
    857     MFEventSelector selectorh;
    858     selectorh.SetNumSelectEvts(howManyHadrons);
    859     selectorh.SetName("selectHadrons");
    860908
    861909    MFillH fillmatg("MatrixGammas");
    862910    fillmatg.SetFilter(&flistg);
    863911    fillmatg.SetName("fillGammas");
    864 
    865     MFillH fillmath("MatrixHadrons");
    866     fillmath.SetFilter(&flisth);
    867     fillmath.SetName("fillHadrons");
    868912
    869913   
    … …  
    894938    MEvtLoop evtloopg;
    895939    evtloopg.SetParList(&plistg);
    896     evtloopg.ReadEnv(env, "", printEnv);
     940    //evtloopg.ReadEnv(env, "", printEnv);
    897941    evtloopg.SetProgressBar(&matrixbar);
    898942
    … …  
    907951
    908952    gLog << " Hadrons :" << endl;
     953
     954    TString mhname("costhh");
     955    MBinning binh("Binning"+mhname);
     956    binh.SetEdges(10, 0., 1.0);
     957
     958    MH3 href("cos(MMcEvt.fTelescopeTheta)");
     959    href.SetName(mhname);
     960    MH::SetBinning(&href.GetHist(), &binh);
     961    for (Int_t i=1; i<=href.GetNbins(); i++)
     962      href.GetHist().SetBinContent(i, 1.0);
     963
     964    MFEventSelector2 selectorh(href);
     965    //selectorh.SetNumMax(howManyHadrons);
     966    // select as many hadrons as gammas
     967    selectorh.SetNumMax(matrixg.GetM().GetNrows());
     968    selectorh.SetName("selectHadrons");
     969
     970    MFillH fillmath("MatrixHadrons");
     971    fillmath.SetFilter(&flisth);
     972    fillmath.SetName("fillHadrons");
     973
     974
    909975    // entries in MFilterList
    910976
    … …  
    932998    MEvtLoop evtlooph;
    933999    evtlooph.SetParList(&plisth);
    934     evtlooph.ReadEnv(env, "", printEnv);
     1000    //evtlooph.ReadEnv(env, "", printEnv);
    9351001    evtlooph.SetProgressBar(&matrixbar);
    9361002
    … …  
    9421008    //***************************************************** 
    9431009
    944     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    945     //
    946     // ----------------------------------------------------------
    947     //  Definition of the reference sample of look-alike events
    948     //  (this is a subsample of the original sample)
    949     // ----------------------------------------------------------
    950     //
     1010
     1011    // write out matrices of training events events
     1012
    9511013    gLog << "" << endl;
    9521014    gLog << "========================================================" << endl;
    953     // Select a maximum of nmaxevts events from the sample of look-alike
    954     // events. They will form the reference sample.
    955     Int_t nmaxevts  = 2000;
    956 
    957     // target distribution for the variable in column refcolumn (start at 0);
    958     //Int_t   refcolumn = 0;
    959     //Int_t   nbins   =   5;
    960     //Float_t frombin = 0.5;
    961     //Float_t tobin   = 1.0;
    962     //TH1F *thsh = new TH1F("thsh","target distribution",
    963     //                       nbins, frombin, tobin);
    964     //thsh->SetDirectory(NULL);
    965     //thsh->SetXTitle("cos( \\Theta)");
    966     //thsh->SetYTitle("Counts");
    967     //Float_t dbin = (tobin-frombin)/nbins;
    968     //Float_t lbin = frombin +dbin*0.5;
    969     //for (Int_t j=1; j<=nbins; j++)
    970     //{
    971     //  thsh->Fill(lbin,1.0);
    972     //  lbin += dbin;
    973     //}
    974 
    975     Int_t   refcolumn = 0;
    976     MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning");
    977     TH1F *thsh = new TH1F();
    978     thsh->SetNameTitle("thsh","target distribution");
    979     MH::SetBinning(thsh, binscostheta);
    980     Int_t nbins = thsh->GetNbinsX();
    981     Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900};
    982     for (Int_t j=1; j<=nbins; j++)
    983     {
    984       thsh->Fill(thsh->GetBinCenter(j), cont[j-1]);
    985     }
    986 
    987     gLog << "" << endl;
    988     gLog << "========================================================" << endl;
    989     gLog << "Macro CT1Analysis : define reference sample for gammas" << endl;
    990     gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
    991          << refcolumn << ",  " << nmaxevts << endl;   
    992 
    993     matrixg.EnableGraphicalOutput();
    994     Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    995 
    996     if ( !matrixok )
    997     {
    998       gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
    999       return;
    1000     }
    1001 
    1002     gLog << "" << endl;
    1003     gLog << "========================================================" << endl;
    1004     gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl;
    1005     gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
    1006          << refcolumn << ",  " << nmaxevts << endl;   
    1007 
    1008     matrixh.EnableGraphicalOutput();
    1009     matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    1010     delete thsh;
    1011 
    1012     if ( !matrixok )
    1013     {
    1014       gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
    1015       return;
    1016     }
    1017     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    1018 
    1019     // write out look-alike events
    1020 
    1021     gLog << "" << endl;
    1022     gLog << "========================================================" << endl;
    1023     gLog << "Write out look=alike events" << endl;
     1015    gLog << "Write out matrices of training events" << endl;
    10241016
    10251017
    … …  
    10271019      // "gammas"
    10281020      gLog << "Gammas :" << endl;   
    1029       matrixg.Print("cols");
     1021      matrixg.Print("SizeCols");
    10301022
    10311023      TFile writeg(outNameGammas, "RECREATE", "");
    … …  
    10331025
    10341026      gLog << "" << endl;
    1035       gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
     1027      gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file "
    10361028           << outNameGammas << endl;
    10371029
    … …  
    10391031      // "hadrons"
    10401032      gLog << "Hadrons :" << endl;   
    1041       matrixh.Print("cols");
     1033      matrixh.Print("SizeCols");
    10421034
    10431035      TFile writeh(outNameHadrons, "RECREATE", "");
    … …  
    10451037
    10461038      gLog << "" << endl;
    1047       gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
     1039      gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file "
    10481040           << outNameHadrons << endl;
    10491041
    10501042  }
    1051    //**********   end of creating matrices of look-alike events   ***********
    1052 
    1053 
     1043   //**********   end of creating matrices of training events   ***********
     1044
     1045
     1046    MRanForest *fRanForest;
     1047    MRanTree *fRanTree;
    10541048    //-----------------------------------------------------------------
    1055     // Update the input files with the NN and SC hadronness
    1056     //
    1057 
    1058  if (WNN)
     1049    // read in the trees of the random forest
     1050    if (RTree)
     1051    {
     1052      MParList plisttr;
     1053      MTaskList tlisttr;
     1054      plisttr.AddToList(&tlisttr);
     1055
     1056      MReadTree readtr("TREE", outRF);
     1057      readtr.DisableAutoScheme();
     1058
     1059      MRanForestFill rffill;
     1060      rffill.SetNumTrees(NumTrees);
     1061
     1062      // list of tasks for the loop over the trees
     1063
     1064      tlisttr.AddToList(&readtr);
     1065      tlisttr.AddToList(&rffill);
     1066
     1067      //-------------------
     1068      // Execute tree loop
     1069      //
     1070      MEvtLoop evtlooptr;
     1071      evtlooptr.SetName("ReadRFTrees");
     1072      evtlooptr.SetParList(&plisttr);
     1073      if (!evtlooptr.Eventloop())
     1074        return;
     1075
     1076      tlisttr.PrintStatistics(0, kTRUE);
     1077
     1078
     1079    // get adresses of objects which are used in the next eventloop
     1080    fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
     1081    if (!fRanForest)
     1082    {
     1083        gLog << err << dbginf << "MRanForest not found... aborting." << endl;
     1084        return kFALSE;
     1085    }
     1086
     1087    fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
     1088    if (!fRanTree)                                 
     1089    {                                                                         
     1090        gLog << err << dbginf << "MRanTree not found... aborting." << endl;   
     1091        return kFALSE;
     1092    }
     1093
     1094    }
     1095
     1096    //-----------------------------------------------------------------
     1097    // grow the trees of the random forest (event loop = tree loop)
     1098
     1099    if (!RTree)
     1100    {
     1101
     1102    gLog << "" << endl;
     1103    gLog << "========================================================" << endl;
     1104    gLog << "Macro CT1Analysis : start growing trees" << endl;
     1105
     1106    MTaskList tlist2;
     1107    MParList plist2;
     1108    plist2.AddToList(&tlist2);
     1109
     1110    plist2.AddToList(&matrixg);
     1111    plist2.AddToList(&matrixh);
     1112
     1113    MRanForestGrow rfgrow2;
     1114    rfgrow2.SetNumTrees(NumTrees);
     1115    rfgrow2.SetNumTry(NumTry);
     1116    rfgrow2.SetNdSize(NdSize);
     1117
     1118    MWriteRootFile rfwrite2(outRF);
     1119    rfwrite2.AddContainer("MRanTree", "TREE");
     1120
     1121    // list of tasks for the loop over the trees
     1122   
     1123    tlist2.AddToList(&rfgrow2);
     1124    tlist2.AddToList(&rfwrite2);
     1125
     1126    //-------------------
     1127    // Execute tree loop
     1128    //
     1129    MEvtLoop treeloop;
     1130    treeloop.SetName("GrowRFTrees");
     1131    treeloop.SetParList(&plist2);
     1132
     1133    if ( !treeloop.Eventloop() )
     1134        return;
     1135
     1136    tlist2.PrintStatistics(0, kTRUE);
     1137
     1138    // get adresses of objects which are used in the next eventloop
     1139    fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
     1140    if (!fRanForest)
     1141    {
     1142        gLog << err << dbginf << "MRanForest not found... aborting." << endl;
     1143        return kFALSE;
     1144    }
     1145
     1146    fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
     1147    if (!fRanTree)                                 
     1148    {                                                                         
     1149        gLog << err << dbginf << "MRanTree not found... aborting." << endl;   
     1150        return kFALSE;
     1151    }
     1152
     1153    }
     1154    // end of growing the trees of the random forest
     1155    //-----------------------------------------------------------------
     1156
     1157
     1158
     1159
     1160    //-----------------------------------------------------------------
     1161    // Update the input files with the RF hadronness
     1162    //
     1163
     1164 if (WRF)
    10591165  {
    10601166    gLog << "" << endl;
    10611167    gLog << "========================================================" << endl;
    10621168    gLog << "Update input file '" <<  filenameData
    1063          << "' with the NN and SC hadronness" << endl;
     1169         << "' with the RF hadronness" << endl;
    10641170
    10651171    MTaskList tliston;
    … …  
    10781184    read.DisableAutoScheme();
    10791185
    1080     TString fHilName    = "MHillas";
    1081     TString fHilNameExt = "MHillasExt";
    1082     TString fHilNameSrc = "MHillasSrc";
    1083     TString fImgParName = "MNewImagePar";
    10841186
    10851187    //.......................................................................
    1086     // calculation of hadroness for method of Nearest Neighbors
    1087 
    1088     TString hadNNName = "HadNN";
    1089     MMultiDimDistCalc nncalc;
    1090     nncalc.SetUseNumRows(25);
    1091     nncalc.SetUseKernelMethod(kFALSE);
    1092     nncalc.SetHadronnessName(hadNNName);
     1188    // calculate hadronnes for method of RANDOM FOREST
     1189
     1190
     1191    MRanForestCalc rfcalc;
     1192    rfcalc.SetHadronnessName(hadRFName);
     1193
    10931194
    10941195    //.......................................................................
    … …  
    11081209      write.AddContainer("MNewImagePar",  "Events");
    11091210
    1110       write.AddContainer("HadRF",         "Events");
    1111       write.AddContainer("HadSC",         "Events");
    1112       write.AddContainer("HadNN",         "Events");
    1113 
     1211      write.AddContainer(hadRFName,       "Events");
    11141212
    11151213    //-----------------------------------------------------------------
    … …  
    11181216             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    11191217
    1120     Float_t maxhadronness =  0.40;
    1121     Float_t maxalpha      =  20.0;
    1122     Float_t maxdist       =  10.0;
     1218
    11231219
    11241220    MFCT1SelFinal selfinalgh(fHilNameSrc);
    11251221    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    1126     selfinalgh.SetHadronnessName(hadNNName);
     1222    selfinalgh.SetHadronnessName(hadRFName);
    11271223    selfinalgh.SetName("SelFinalgh");
    11281224    MContinue contfinalgh(&selfinalgh);
    11291225    contfinalgh.SetName("ContSelFinalgh");
    11301226
    1131     MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
    1132     fillhadnn.SetName("HhadNN");
     1227    MFillH fillranfor("MHRanForest");
     1228    fillranfor.SetName("HRanForest");
     1229
     1230    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
     1231    fillhadrf.SetName("HhadRF");
    11331232
    11341233    MFCT1SelFinal selfinal(fHilNameSrc);
    11351234    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    1136     selfinal.SetHadronnessName(hadNNName);
     1235    selfinal.SetHadronnessName(hadRFName);
    11371236    selfinal.SetName("SelFinal");
    11381237    MContinue contfinal(&selfinal);
    11391238    contfinal.SetName("ContSelFinal");
     1239
     1240    TString mh3name = "abs(Alpha)";
     1241    MBinning binsalphaabs("Binning"+mh3name);
     1242    binsalphaabs.SetEdges(50, -2.0, 98.0);
     1243
     1244    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
     1245    alphaabs.SetName(mh3name);
     1246    MFillH alpha(&alphaabs);
     1247    alpha.SetName("FillAlphaAbs");
    11401248
    11411249
    … …  
    11611269    InitBinnings(&pliston);
    11621270
    1163     pliston.AddToList(&matrixg);
    1164     pliston.AddToList(&matrixh);
     1271    pliston.AddToList(fRanForest);
     1272    pliston.AddToList(fRanTree);
     1273
     1274    pliston.AddToList(&binsalphaabs);
     1275    pliston.AddToList(&alphaabs);
    11651276
    11661277
    … …  
    11701281    tliston.AddToList(&read);
    11711282
    1172     tliston.AddToList(&nncalc);
    1173     tliston.AddToList(&fillhadnn);
    1174 
     1283    tliston.AddToList(&rfcalc);
     1284    tliston.AddToList(&fillranfor);
     1285    tliston.AddToList(&fillhadrf);
     1286
     1287    tliston.AddToList(&write);
     1288    tliston.AddToList(&contfinalgh);
     1289
     1290    tliston.AddToList(&alpha);
    11751291    tliston.AddToList(&hfill1);
    11761292    tliston.AddToList(&hfill2);
    … …  
    11791295    tliston.AddToList(&hfill5);
    11801296
    1181     tliston.AddToList(&write);
    1182 
    1183     tliston.AddToList(&contfinalgh);
    11841297    tliston.AddToList(&contfinal);
    11851298
    … …  
    12001313    tliston.PrintStatistics(0, kTRUE);
    12011314
     1315
    12021316    //-------------------------------------------
    12031317    // Display the histograms
    12041318    //
    1205     pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     1319    pliston.FindObject("MHRanForest")->DrawClone();
     1320    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    12061321
    12071322    pliston.FindObject("MHHillas")->DrawClone();
    … …  
    12111326    pliston.FindObject("MHStarMap")->DrawClone();
    12121327
     1328
     1329     //-------------------------------------------
     1330    // fit alpha distribution to get the number of excess events and
     1331    // calculate significance of gamma signal in the alpha plot
     1332 
     1333    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
     1334    TH1  &alphaHist = absalpha->GetHist();
     1335    alphaHist.SetXTitle("|alpha|  [\\circ]");
     1336    alphaHist.SetName("alpha-macro");
     1337
     1338    Double_t alphasig = 13.1;
     1339    Double_t alphamin = 30.0;
     1340    Double_t alphamax = 90.0;
     1341    Int_t    degree   =    2;
     1342    Double_t significance = -99.0;
     1343    Bool_t   drawpoly  = kTRUE;
     1344    Bool_t   fitgauss  = kTRUE;
     1345    Bool_t   print     = kTRUE;
     1346
     1347    MHFindSignificance findsig;
     1348    findsig.SetRebin(kTRUE);
     1349    findsig.SetReduceDegree(kFALSE);
     1350
     1351    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
     1352                        alphasig, drawpoly, fitgauss, print);
     1353    significance = findsig.GetSignificance();
     1354    Float_t alphasi = findsig.GetAlphasi();
     1355
     1356    gLog << "For file '" << filenameData << "' : " << endl;
     1357    gLog << "Significance of gamma signal after supercuts : "
     1358         << significance << " (for |alpha| < " << alphasi << " degrees)"
     1359         << endl;
     1360
     1361    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
     1362
     1363    //-------------------------------------------
     1364
     1365
    12131366    DeleteBinnings(&pliston);
    12141367  }
    12151368
    1216 
    1217 
    1218     gLog << "Macro CT1Analysis : End of Job B_NN_UP" << endl;
     1369    gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
    12191370    gLog << "=======================================================" << endl;
    12201371 }
    … …  
    12221373
    12231374
     1375
    12241376  //---------------------------------------------------------------------
    1225   // Job B_RF_UP
    1226   //============
    1227 
    1228 
    1229     //  - create (or read in) the matrices of look-alike events for gammas
    1230     //    and hadrons
    1231     //  - create (or read in) the trees
    1232     //  - then read ON1.root (or MC1.root) file
    1233     //  - calculate the hadroness for the method of RANDOM FOREST
    1234     //  - update input root file with the hadroness
    1235 
    1236 
    1237  if (JobB_RF_UP)
     1377  // Job C 
     1378  //======
     1379
     1380    //  - read ON1 and MC1 data files 
     1381    //    which should have been updated to contain the hadronnesses
     1382    //    for the method of NEAREST NEIGHBORS and for the SUOERCUTS
     1383    //  - produce Neyman-Pearson plots
     1384 
     1385 if (JobC)
    12381386 {
    12391387    gLog << "=====================================================" << endl;
    1240     gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
     1388    gLog << "Macro CT1Analysis : Start of Job C" << endl;
    12411389
    12421390    gLog << "" << endl;
    1243     gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = "
    1244          << JobB_RF_UP  << ",  " << RLookRF << ",  " << RTree << ",  "
    1245          << WRF << endl;
    1246 
    1247 
    1248 
    1249     //--------------------------------------------
    1250     // parameters for the random forest
    1251     Int_t NumTrees = 100;
    1252     Int_t NumTry   =   3;
    1253     Int_t NdSize   =   3;
    1254 
    1255 
    1256     //--------------------------------------------
    1257     // file to be updated (ON, OFF or MC)
    1258 
    1259     TString typeInput = "ON";
    1260     //TString typeInput = "OFF";
    1261     //TString typeInput = "MC";
    1262     gLog << "typeInput = " << typeInput << endl;
    1263 
    1264     // name of input root file
     1391    gLog << "Macro CT1Analysis : JobC = " << JobC  << endl;
     1392
     1393
     1394    // name of input data file
    12651395    TString filenameData = outPath;
    1266     filenameData += typeInput;
    1267     filenameData += "1.root";
    1268     gLog << "filenameData = " << filenameData << endl;
    1269 
    1270     // name of output root file
    1271     TString outNameImage = outPath;
    1272     outNameImage += typeInput;
    1273     outNameImage += "2.root";
    1274     //TString outNameImage = filenameData;
    1275 
    1276     gLog << "outNameImage = " << outNameImage << endl;
    1277 
    1278     //--------------------------------------------
    1279     // files to be read for generating the look-alike events
    1280     // "hadrons" :
    1281     TString filenameHad = outPath;
    1282     filenameHad += "OFF";
    1283     filenameHad += "1.root";
    1284     Int_t howManyHadrons = 1000000;
    1285     gLog << "filenameHad = " << filenameHad << ",   howManyHadrons = "
    1286          << howManyHadrons  << endl;
    1287    
    1288 
    1289     // "gammas" :
     1396    filenameData += "OFF";
     1397    filenameData += "2.root";
     1398    gLog << "filenameData = " << filenameData << endl;
     1399
     1400    // name of input MC file
    12901401    TString filenameMC = outPath;
    12911402    filenameMC += "MC";
    1292     filenameMC += "1.root";
    1293     Int_t howManyGammas = 50000;
    1294     gLog << "filenameMC = " << filenameMC << ",   howManyGammas = "
    1295          << howManyGammas   << endl;
    1296    
    1297     //--------------------------------------------
    1298     // files of look-alike events
    1299 
    1300     TString outNameGammas = outPath;
    1301     outNameGammas += "RFmatrix_gammas_";
    1302     outNameGammas += "MC";
    1303     outNameGammas += ".root";
    1304 
    1305     TString typeMatrixHadrons = "OFF";
    1306     gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
    1307 
    1308     TString outNameHadrons = outPath;
    1309     outNameHadrons += "RFmatrix_hadrons_";
    1310     outNameHadrons += typeMatrixHadrons;
    1311     outNameHadrons += ".root";
    1312 
    1313 
    1314     MHMatrix matrixg("MatrixGammas");
    1315     MHMatrix matrixh("MatrixHadrons");
    1316 
    1317     //--------------------------------------------
    1318     // file of trees of the random forest
    1319 
    1320     TString outRF = outPath;
    1321     outRF += "RF.root";
    1322 
    1323 
    1324    //*************************************************************************
    1325    // read in matrices of look-alike events
    1326 if (RLookRF)
    1327   {
    1328     const char* mtxName = "MatrixGammas";
    1329 
    1330     gLog << "" << endl;
    1331     gLog << "========================================================" << endl;
    1332     gLog << "Get matrix for (gammas)" << endl;
    1333     gLog << "matrix name        = " << mtxName << endl;
    1334     gLog << "name of root file  = " << outNameGammas << endl;
    1335     gLog << "" << endl;
    1336 
    1337 
    1338     // read in the object with the name 'mtxName' from file 'outNameGammas'
    1339     //
    1340     TFile fileg(outNameGammas);
    1341 
    1342     matrixg.Read(mtxName);
    1343     matrixg.Print("cols");
    1344 
    1345 
    1346     //*****************************************************************
    1347 
    1348     const char* mtxName = "MatrixHadrons";
    1349 
    1350     gLog << "" << endl;
    1351     gLog << "========================================================" << endl;
    1352     gLog << " Get matrix for (hadrons)" << endl;
    1353     gLog << "matrix name        = " << mtxName << endl;
    1354     gLog << "name of root file  = " << outNameHadrons << endl;
    1355     gLog << "" << endl;
    1356 
    1357 
    1358     // read in the object with the name 'mtxName' from file 'outNameHadrons'
    1359     //
    1360     TFile fileh(outNameHadrons);
    1361 
    1362     matrixh.Read(mtxName);
    1363     matrixh.Print("cols");
    1364   }
    1365 
    1366 
    1367    //*************************************************************************
    1368    // create matrices of look-alike events
    1369 if (!RLookRF)
    1370   {
    1371     gLog << "" << endl;
    1372     gLog << "========================================================" << endl;
    1373     gLog << " Create matrices of look-alike events" << endl;
    1374     gLog << " Gammas :" << endl;
    1375 
    1376 
    1377     MParList  plistg;
    1378     MTaskList tlistg;
    1379     MFilterList flistg;
    1380 
    1381     MParList  plisth;
    1382     MTaskList tlisth;
    1383     MFilterList flisth;
    1384 
    1385     MReadMarsFile  readg("Events", filenameMC);
    1386     readg.DisableAutoScheme();
    1387 
    1388     MReadMarsFile  readh("Events", filenameHad);
    1389     readh.DisableAutoScheme();
    1390 
    1391     MFParticleId fgamma("MMcEvt", '=', kGAMMA);
    1392     fgamma.SetName("gammaID");
    1393     MFParticleId fhadrons("MMcEvt", '!', kGAMMA);
    1394     fhadrons.SetName("hadronID)");
    1395 
    1396     MFEventSelector selectorg;
    1397     selectorg.SetNumSelectEvts(howManyGammas);
    1398     selectorg.SetName("selectGammas");
    1399     MFEventSelector selectorh;
    1400     selectorh.SetNumSelectEvts(howManyHadrons);
    1401     selectorh.SetName("selectHadrons");
    1402 
    1403     MFillH fillmatg("MatrixGammas");
    1404     fillmatg.SetFilter(&flistg);
    1405     fillmatg.SetName("fillGammas");
    1406 
    1407     MFillH fillmath("MatrixHadrons");
    1408     fillmath.SetFilter(&flisth);
    1409     fillmath.SetName("fillHadrons");
    1410 
    1411    
    1412     //*****************************   fill gammas   *** 
    1413     // entries in MFilterList
    1414 
    1415     flistg.AddToList(&fgamma);
    1416     flistg.AddToList(&selectorg);
    1417 
    1418     //***************************** 
    1419     // entries in MParList
    1420    
    1421     plistg.AddToList(&tlistg);
    1422     InitBinnings(&plistg);
    1423 
    1424     plistg.AddToList(&matrixg);
    1425 
    1426     //*****************************
    1427     // entries in MTaskList
    1428    
    1429     tlistg.AddToList(&readg);
    1430     tlistg.AddToList(&flistg);
    1431     tlistg.AddToList(&fillmatg);
    1432 
    1433     //*****************************
    1434 
    1435     MProgressBar matrixbar;
    1436     MEvtLoop evtloopg;
    1437     evtloopg.SetParList(&plistg);
    1438     evtloopg.ReadEnv(env, "", printEnv);
    1439     evtloopg.SetProgressBar(&matrixbar);
    1440 
    1441     Int_t maxevents = -1;
    1442     if (!evtloopg.Eventloop(maxevents))
    1443         return;
    1444 
    1445     tlistg.PrintStatistics(0, kTRUE);
    1446 
    1447 
    1448     //*****************************   fill hadrons   *** 
    1449 
    1450     gLog << " Hadrons :" << endl;
    1451     // entries in MFilterList
    1452 
    1453     flisth.AddToList(&fhadrons);
    1454     flisth.AddToList(&selectorh);
    1455 
    1456     //***************************** 
    1457     // entries in MParList
    1458    
    1459     plisth.AddToList(&tlisth);
    1460     InitBinnings(&plisth);
    1461 
    1462     plisth.AddToList(&matrixh);
    1463 
    1464     //*****************************
    1465     // entries in MTaskList
    1466    
    1467     tlisth.AddToList(&readh);
    1468     tlisth.AddToList(&flisth);
    1469     tlisth.AddToList(&fillmath);
    1470 
    1471     //*****************************
    1472 
    1473     MProgressBar matrixbar;
    1474     MEvtLoop evtlooph;
    1475     evtlooph.SetParList(&plisth);
    1476     evtlooph.ReadEnv(env, "", printEnv);
    1477     evtlooph.SetProgressBar(&matrixbar);
    1478 
    1479     Int_t maxevents = -1;
    1480     if (!evtlooph.Eventloop(maxevents))
    1481         return;
    1482 
    1483     tlisth.PrintStatistics(0, kTRUE);
    1484     //***************************************************** 
    1485 
    1486     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    1487     //
    1488     // ----------------------------------------------------------
    1489     //  Definition of the reference sample of look-alike events
    1490     //  (this is a subsample of the original sample)
    1491     // ----------------------------------------------------------
    1492     //
    1493     gLog << "" << endl;
    1494     gLog << "========================================================" << endl;
    1495     // Select a maximum of nmaxevts events from the sample of look-alike
    1496     // events. They will form the reference sample.
    1497     Int_t nmaxevts  = 10000;
    1498 
    1499     // target distribution for the variable in column refcolumn (start at 0);
    1500     //Int_t   refcolumn = 0;
    1501     //Int_t   nbins   =   5;
    1502     //Float_t frombin = 0.5;
    1503     //Float_t tobin   = 1.0;
    1504     //TH1F *thsh = new TH1F("thsh","target distribution",
    1505     //                       nbins, frombin, tobin);
    1506     //thsh->SetDirectory(NULL);
    1507     //thsh->SetXTitle("cos( \\Theta)");
    1508     //thsh->SetYTitle("Counts");
    1509     //Float_t dbin = (tobin-frombin)/nbins;
    1510     //Float_t lbin = frombin +dbin*0.5;
    1511     //for (Int_t j=1; j<=nbins; j++)
    1512     //{
    1513     //  thsh->Fill(lbin,1.0);
    1514     //  lbin += dbin;
    1515     //}
    1516 
    1517     Int_t   refcolumn = 0;
    1518     MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning");
    1519     TH1F *thsh = new TH1F();
    1520     thsh->SetNameTitle("thsh","target distribution");
    1521     MH::SetBinning(thsh, binscostheta);
    1522     Int_t nbins = thsh->GetNbinsX();
    1523     Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900};
    1524     for (Int_t j=1; j<=nbins; j++)
    1525     {
    1526       thsh->Fill(thsh->GetBinCenter(j), cont[j-1]);
    1527     }
    1528 
    1529     gLog << "" << endl;
    1530     gLog << "========================================================" << endl;
    1531     gLog << "Macro CT1Analysis : define reference sample for gammas" << endl;
    1532     gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
    1533          << refcolumn << ",  " << nmaxevts << endl;   
    1534 
    1535     matrixg.EnableGraphicalOutput();
    1536     Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    1537 
    1538     if ( !matrixok )
    1539     {
    1540       gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
    1541       return;
    1542     }
    1543 
    1544     gLog << "" << endl;
    1545     gLog << "========================================================" << endl;
    1546     gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl;
    1547     gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
    1548          << refcolumn << ",  " << nmaxevts << endl;   
    1549 
    1550     matrixh.EnableGraphicalOutput();
    1551     matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    1552     delete thsh;
    1553 
    1554     if ( !matrixok )
    1555     {
    1556       gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
    1557       return;
    1558     }
    1559     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    1560 
    1561     // write out look-alike events
    1562 
    1563     gLog << "" << endl;
    1564     gLog << "========================================================" << endl;
    1565     gLog << "Write out look=alike events" << endl;
    1566 
    1567 
    1568       //-------------------------------------------
    1569       // "gammas"
    1570       gLog << "Gammas :" << endl;   
    1571       matrixg.Print("cols");
    1572 
    1573       TFile writeg(outNameGammas, "RECREATE", "");
    1574       matrixg.Write();
    1575 
    1576       gLog << "" << endl;
    1577       gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
    1578            << outNameGammas << endl;
    1579 
    1580       //-------------------------------------------
    1581       // "hadrons"
    1582       gLog << "Hadrons :" << endl;   
    1583       matrixh.Print("cols");
    1584 
    1585       TFile writeh(outNameHadrons, "RECREATE", "");
    1586       matrixh.Write();
    1587 
    1588       gLog << "" << endl;
    1589       gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
    1590            << outNameHadrons << endl;
    1591 
    1592   }
    1593    //**********   end of creating matrices of look-alike events   ***********
    1594 
    1595 
    1596     MRanForest *fRanForest;
    1597     MRanTree *fRanTree;
     1403    filenameMC += "2.root";
     1404    gLog << "filenameMC   = " << filenameMC   << endl;
     1405
     1406
    15981407    //-----------------------------------------------------------------
    1599     // read in the trees of the random forest
    1600     if (RTree)
    1601     {
    1602       MParList plisttr;
    1603       MTaskList tlisttr;
    1604       plisttr.AddToList(&tlisttr);
    1605 
    1606       MReadTree readtr("TREE", outRF);
    1607       readtr.DisableAutoScheme();
    1608 
    1609       MRanForestFill rffill;
    1610       rffill.SetNumTrees(NumTrees);
    1611 
    1612       // list of tasks for the loop over the trees
    1613 
    1614       tlisttr.AddToList(&readtr);
    1615       tlisttr.AddToList(&rffill);
    1616 
    1617       //-------------------
    1618       // Execute tree loop
    1619       //
    1620       MEvtLoop evtlooptr;
    1621       evtlooptr.SetParList(&plisttr);
    1622       if (!evtlooptr.Eventloop())
    1623         return;
    1624 
    1625       tlisttr.PrintStatistics(0, kTRUE);
    1626 
    1627 
    1628     // get adresses of objects which are used in the next eventloop
    1629     fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
    1630     if (!fRanForest)
    1631     {
    1632         gLog << err << dbginf << "MRanForest not found... aborting." << endl;
    1633         return kFALSE;
    1634     }
    1635 
    1636     fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
    1637     if (!fRanTree)                                 
    1638     {                                                                         
    1639         gLog << err << dbginf << "MRanTree not found... aborting." << endl;   
    1640         return kFALSE;
    1641     }
    1642 
    1643     }
    1644 
    1645     //-----------------------------------------------------------------
    1646     // grow the trees of the random forest (event loop = tree loop)
    1647 
    1648     if (!RTree)
    1649     {
    1650 
    1651     gLog << "" << endl;
    1652     gLog << "========================================================" << endl;
    1653     gLog << "Macro CT1Analysis : start growing trees" << endl;
    1654 
    1655     MTaskList tlist2;
    1656     MParList plist2;
    1657     plist2.AddToList(&tlist2);
    1658 
    1659     plist2.AddToList(&matrixg);
    1660     plist2.AddToList(&matrixh);
    1661 
    1662     MRanForestGrow rfgrow2;
    1663     rfgrow2.SetNumTrees(NumTrees);
    1664     rfgrow2.SetNumTry(NumTry);
    1665     rfgrow2.SetNdSize(NdSize);
    1666 
    1667     MWriteRootFile rfwrite2(outRF);
    1668     rfwrite2.AddContainer("MRanTree", "TREE");
    1669 
    1670     // list of tasks for the loop over the trees
    1671    
    1672     tlist2.AddToList(&rfgrow2);
    1673     tlist2.AddToList(&rfwrite2);
    1674 
    1675     //-------------------
    1676     // Execute tree loop
    1677     //
    1678     MEvtLoop treeloop;
    1679     treeloop.SetParList(&plist2);
    1680 
    1681     if ( !treeloop.Eventloop() )
    1682         return;
    1683 
    1684     tlist2.PrintStatistics(0, kTRUE);
    1685 
    1686     // get adresses of objects which are used in the next eventloop
    1687     fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
    1688     if (!fRanForest)
    1689     {
    1690         gLog << err << dbginf << "MRanForest not found... aborting." << endl;
    1691         return kFALSE;
    1692     }
    1693 
    1694     fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
    1695     if (!fRanTree)                                 
    1696     {                                                                         
    1697         gLog << err << dbginf << "MRanTree not found... aborting." << endl;   
    1698         return kFALSE;
    1699     }
    1700 
    1701     }
    1702     // end of growing the trees of the random forest
    1703     //-----------------------------------------------------------------
    1704 
    1705 
    1706 
    1707 
    1708     //-----------------------------------------------------------------
    1709     // Update the input files with the RF hadronness
    1710     //
    1711 
    1712  if (WRF)
    1713   {
    1714     gLog << "" << endl;
    1715     gLog << "========================================================" << endl;
    1716     gLog << "Update input file '" <<  filenameData
    1717          << "' with the RF hadronness" << endl;
    17181408
    17191409    MTaskList tliston;
    … …  
    17291419    //
    17301420
    1731     MReadMarsFile read("Events", filenameData);
     1421    MReadMarsFile read("Events", filenameMC);
     1422    read.AddFile(filenameData);
    17321423    read.DisableAutoScheme();
     1424
     1425
     1426    //.......................................................................
     1427    // names of hadronness containers
     1428
     1429    TString hadNNName = "HadNN";
     1430    TString hadSCName = "HadSC";
     1431    TString hadRFName = "HadRF";
     1432
     1433    //.......................................................................
     1434
    17331435
    17341436    TString fHilName    = "MHillas";
    … …  
    17371439    TString fImgParName = "MNewImagePar";
    17381440
     1441    Float_t maxhadronness =  0.40;
     1442    Float_t maxalpha      =  20.0;
     1443    Float_t maxdist       =  10.0;
     1444
     1445    MFCT1SelFinal selfinalgh(fHilNameSrc);
     1446    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
     1447    selfinalgh.SetHadronnessName(hadRFName);
     1448    selfinalgh.SetName("SelFinalgh");
     1449    MContinue contfinalgh(&selfinalgh);
     1450    contfinalgh.SetName("ContSelFinalgh");
     1451
     1452    //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
     1453    //fillhadnn.SetName("HhadNN");
     1454    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
     1455    fillhadsc.SetName("HhadSC");
     1456    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
     1457    fillhadrf.SetName("HhadRF");
     1458
     1459    MFCT1SelFinal selfinal(fHilNameSrc);
     1460    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
     1461    selfinal.SetHadronnessName(hadRFName);
     1462    selfinal.SetName("SelFinal");
     1463    MContinue contfinal(&selfinal);
     1464    contfinal.SetName("ContSelFinal");
     1465
     1466
     1467    MFillH hfill1("MHHillas",    fHilName);
     1468    hfill1.SetName("HHillas");
     1469
     1470    MFillH hfill2("MHStarMap",   fHilName);
     1471    hfill2.SetName("HStarMap");
     1472
     1473    MFillH hfill3("MHHillasExt",    fHilNameSrc);
     1474    hfill3.SetName("HHillasExt");
     1475   
     1476    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
     1477    hfill4.SetName("HHillasSrc");   
     1478
     1479    MFillH hfill5("MHNewImagePar", fImgParName);
     1480    hfill5.SetName("HNewImagePar");
     1481
     1482
     1483    //*****************************
     1484    // entries in MParList
     1485
     1486    pliston.AddToList(&tliston);
     1487    InitBinnings(&pliston);
     1488
     1489
     1490    //*****************************
     1491    // entries in MTaskList
     1492   
     1493    tliston.AddToList(&read);
     1494
     1495    //tliston.AddToList(&fillhadnn);
     1496    tliston.AddToList(&fillhadsc);
     1497    tliston.AddToList(&fillhadrf);
     1498   
     1499    tliston.AddToList(&contfinalgh);
     1500    tliston.AddToList(&hfill1);
     1501    tliston.AddToList(&hfill2);
     1502    tliston.AddToList(&hfill3);
     1503    tliston.AddToList(&hfill4);
     1504    tliston.AddToList(&hfill5);
     1505
     1506    tliston.AddToList(&contfinal);
     1507
     1508    //*****************************
     1509
     1510    //-------------------------------------------
     1511    // Execute event loop
     1512    //
     1513    MProgressBar bar;
     1514    MEvtLoop evtloop;
     1515    evtloop.SetParList(&pliston);
     1516    evtloop.SetProgressBar(&bar);
     1517
     1518    Int_t maxevents = -1;
     1519    //Int_t maxevents = 35000;
     1520    if ( !evtloop.Eventloop(maxevents) )
     1521        return;
     1522
     1523    tliston.PrintStatistics(0, kTRUE);
     1524
     1525
     1526    //-------------------------------------------
     1527    // Display the histograms
     1528    //
     1529
     1530    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     1531    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
     1532    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
     1533
     1534    pliston.FindObject("MHHillas")->DrawClone();
     1535    pliston.FindObject("MHHillasExt")->DrawClone();
     1536    pliston.FindObject("MHHillasSrc")->DrawClone();
     1537    pliston.FindObject("MHNewImagePar")->DrawClone();
     1538    pliston.FindObject("MHStarMap")->DrawClone();
     1539
     1540    DeleteBinnings(&pliston);
     1541
     1542    gLog << "Macro CT1Analysis : End of Job C" << endl;
     1543    gLog << "===================================================" << endl;
     1544 }
     1545
     1546
     1547  //---------------------------------------------------------------------
     1548  // Job D
     1549  //======
     1550
     1551    //  - select g/h separation method XX
     1552    //  - read ON2 (or MC2) root file
     1553    //  - apply cuts in hadronness
     1554    //  - make plots
     1555
     1556
     1557 if (JobD)
     1558 {
     1559    gLog << "=====================================================" << endl;
     1560    gLog << "Macro CT1Analysis : Start of Job D" << endl;
     1561
     1562    gLog << "" << endl;
     1563    gLog << "Macro CT1Analysis : JobD = "
     1564         << JobD  << endl;
     1565
     1566    // type of data to be analysed
     1567    TString typeData = "ON";
     1568    //TString typeData = "OFF";
     1569    //TString typeData = "MC";
     1570    gLog << "typeData = " << typeData << endl;
     1571
     1572    TString ext      = "2.root";
     1573
     1574
     1575    //------------------------------
     1576    // selection of g/h separation method
     1577    // and definition of final selections
     1578
     1579    //TString XX("NN");
     1580    //TString XX("SC");
     1581    TString XX("RF");
     1582    TString fhadronnessName("Had");
     1583    fhadronnessName += XX;
     1584    gLog << "fhadronnessName = " << fhadronnessName << endl;
     1585
     1586    // maximum values of the hadronness, |ALPHA| and DIST
     1587    Float_t maxhadronness   = 0.30;
     1588    Float_t maxalpha        = 20.0;
     1589    Float_t maxdist         = 10.0;
     1590    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
     1591         << maxhadronness << ",  " << maxalpha << ",  "
     1592         << maxdist << endl;
     1593
     1594
     1595    //------------------------------
     1596    // name of data file to be analysed
     1597    TString filenameData(outPath);
     1598    filenameData += typeData;
     1599    filenameData += ext;
     1600    gLog << "filenameData = " << filenameData << endl;
     1601
     1602
     1603
     1604    //*************************************************************************
     1605    //
     1606    // Analyse the data
     1607    //
     1608
     1609    MTaskList tliston;
     1610    MParList pliston;
     1611
     1612    // geometry is needed in  MHHillas... classes
     1613    MGeomCam *fGeom =
     1614             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     1615
     1616
     1617    TString fHilName    = "MHillas";
     1618    TString fHilNameExt = "MHillasExt";
     1619    TString fHilNameSrc = "MHillasSrc";
     1620    TString fImgParName = "MNewImagePar";
     1621
     1622    //-------------------------------------------
     1623    // create the tasks which should be executed
     1624    //
     1625
     1626    MReadMarsFile read("Events", filenameData);
     1627    read.DisableAutoScheme();
     1628
     1629
     1630    //-----------------------------------------------------------------
     1631    // geometry is needed in  MHHillas... classes
     1632    MGeomCam *fGeom =
     1633             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     1634
     1635    MFCT1SelFinal selfinalgh(fHilNameSrc);
     1636    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
     1637    selfinalgh.SetHadronnessName(fhadronnessName);
     1638    selfinalgh.SetName("SelFinalgh");
     1639    MContinue contfinalgh(&selfinalgh);
     1640    contfinalgh.SetName("ContSelFinalgh");
     1641
     1642    //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
     1643    //fillhadnn.SetName("HhadNN");
     1644    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
     1645    fillhadsc.SetName("HhadSC");
     1646    MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
     1647    fillhadrf.SetName("HhadRF");
     1648
     1649
     1650    MFillH hfill1("MHHillas",    fHilName);
     1651    hfill1.SetName("HHillas");
     1652
     1653    MFillH hfill2("MHStarMap",   fHilName);
     1654    hfill2.SetName("HStarMap");
     1655
     1656    MFillH hfill3("MHHillasExt",   fHilNameSrc);
     1657    hfill3.SetName("HHillasExt");   
     1658
     1659    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
     1660    hfill4.SetName("HHillasSrc");   
     1661
     1662    MFillH hfill5("MHNewImagePar", fImgParName);
     1663    hfill5.SetName("HNewImagePar");   
     1664
     1665    MFCT1SelFinal selfinal(fHilNameSrc);
     1666    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
     1667    selfinal.SetHadronnessName(fhadronnessName);
     1668    selfinal.SetName("SelFinal");
     1669    MContinue contfinal(&selfinal);
     1670    contfinal.SetName("ContSelFinal");
     1671
     1672
     1673    //*****************************
     1674    // entries in MParList
     1675
     1676    pliston.AddToList(&tliston);
     1677    InitBinnings(&pliston);
     1678
     1679
     1680    //*****************************
     1681    // entries in MTaskList
     1682   
     1683    tliston.AddToList(&read);
     1684
     1685    tliston.AddToList(&contfinalgh);
     1686
     1687    //tliston.AddToList(&fillhadnn);
     1688    tliston.AddToList(&fillhadsc);
     1689    tliston.AddToList(&fillhadrf);
     1690
     1691    tliston.AddToList(&hfill1);
     1692    tliston.AddToList(&hfill2);
     1693    tliston.AddToList(&hfill3);
     1694    tliston.AddToList(&hfill4);
     1695    tliston.AddToList(&hfill5);
     1696
     1697    tliston.AddToList(&contfinal);
     1698
     1699    //*****************************
     1700
     1701    //-------------------------------------------
     1702    // Execute event loop
     1703    //
     1704    MProgressBar bar;
     1705    MEvtLoop evtloop;
     1706    evtloop.SetParList(&pliston);
     1707    evtloop.SetProgressBar(&bar);
     1708
     1709    Int_t maxevents = -1;
     1710    //Int_t maxevents = 10000;
     1711    if ( !evtloop.Eventloop(maxevents) )
     1712        return;
     1713
     1714    tliston.PrintStatistics(0, kTRUE);
     1715
     1716
     1717
     1718    //-------------------------------------------
     1719    // Display the histograms
     1720    //
     1721
     1722    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     1723    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
     1724    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
     1725
     1726    pliston.FindObject("MHHillas")->DrawClone();
     1727    pliston.FindObject("MHHillasExt")->DrawClone();
     1728    pliston.FindObject("MHHillasSrc")->DrawClone();
     1729    pliston.FindObject("MHNewImagePar")->DrawClone();
     1730    pliston.FindObject("MHStarMap")->DrawClone();
     1731
     1732    //-------------------------------------------
     1733    // fit alpha distribution to get the number of excess events
     1734    //
     1735
     1736    MHHillasSrc* hillasSrc =
     1737      (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
     1738    TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
     1739 
     1740    MHOnSubtraction onsub;
     1741    onsub.Calc(alphaHist, &pliston, kTRUE, 13.1);
     1742    //-------------------------------------------
     1743
     1744    DeleteBinnings(&pliston);
     1745
     1746    gLog << "Macro CT1Analysis : End of Job D" << endl;
     1747    gLog << "=======================================================" << endl;
     1748 }
     1749  //---------------------------------------------------------------------
     1750
     1751
     1752  //---------------------------------------------------------------------
     1753  // Job E_EST_UP
     1754  //============
     1755
     1756    //  - read MC2.root file
     1757    //  - select g/h separation method XX
     1758    //  - optimize energy estimator for events passing final cuts
     1759    //  - write parameters of energy estimator onto file "energyest_XX.root"
     1760    //
     1761    //  - read ON2.root and MC2.root files
     1762    //  - update input root file with the estimated energy
     1763    //    (ON_XX2.root, MC_XX2.root)
     1764
     1765
     1766 if (JobE_EST_UP)
     1767 {
     1768    gLog << "=====================================================" << endl;
     1769    gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl;
     1770
     1771    gLog << "" << endl;
     1772    gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = "
     1773         << JobE_EST_UP  << ",  " << WESTUP << endl;
     1774
     1775
     1776    TString typeON  = "ON";
     1777    TString typeOFF = "OFF";
     1778    TString typeMC  = "MC";
     1779    TString ext    = "2.root";
     1780    TString extout = "3.root";
     1781
     1782    //------------------------------
     1783    // name of MC file to be used for optimizing the energy estimator
     1784    TString filenameOpt(outPath);
     1785    filenameOpt += typeMC;
     1786    filenameOpt += ext;
     1787    gLog << "filenameOpt = " << filenameOpt << endl;
     1788
     1789    //------------------------------
     1790    // selection of g/h separation method
     1791    // and definition of final selections
     1792
     1793    //TString XX("NN");
     1794    //TString XX("SC");
     1795    TString XX("RF");
     1796    TString fhadronnessName("Had");
     1797    fhadronnessName += XX;
     1798    gLog << "fhadronnessName = " << fhadronnessName << endl;
     1799
     1800    // maximum values of the hadronness, |alpha| and dist
     1801    Float_t maxhadronness   = 0.40;
     1802    Float_t maxalpha        = 20.0;
     1803    Float_t maxdist         = 10.0;
     1804    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
     1805         << maxhadronness << ",  " << maxalpha << ",  "
     1806         << maxdist << endl;
     1807
     1808    // name of file containing the parameters of the energy estimator
     1809    TString energyParName(outPath);
     1810    energyParName += "energyest_";
     1811    energyParName += XX;
     1812    energyParName += ".root";
     1813    gLog << "energyParName = " << energyParName << endl;
     1814
     1815
     1816    //------------------------------
     1817    // name of ON file to be updated
     1818    TString filenameON(outPath);
     1819    filenameON += typeON;
     1820    filenameON += ext;
     1821    gLog << "filenameON = " << filenameON << endl;
     1822
     1823    // name of OFF file to be updated
     1824    TString filenameOFF(outPath);
     1825    filenameOFF += typeOFF;
     1826    filenameOFF += ext;
     1827    gLog << "filenameOFF = " << filenameOFF << endl;
     1828
     1829    // name of MC file to be updated
     1830    TString filenameMC(outPath);
     1831    filenameMC += typeMC;
     1832    filenameMC += ext;
     1833    gLog << "filenameMC = " << filenameMC << endl;
     1834
     1835    //------------------------------
     1836    // name of updated ON file
     1837    TString filenameONup(outPath);
     1838    filenameONup += typeON;
     1839    filenameONup += "_";
     1840    filenameONup += XX;
     1841    filenameONup += extout;
     1842    gLog << "filenameONup = " << filenameONup << endl;
     1843
     1844    // name of updated OFF file
     1845    TString filenameOFFup(outPath);
     1846    filenameOFFup += typeOFF;
     1847    filenameOFFup += "_";
     1848    filenameOFFup += XX;
     1849    filenameOFFup += extout;
     1850    gLog << "filenameOFFup = " << filenameOFFup << endl;
     1851
     1852    // name of updated MC file
     1853    TString filenameMCup(outPath);
     1854    filenameMCup += typeMC;
     1855    filenameMCup += "_";
     1856    filenameMCup += XX;
     1857    filenameMCup += extout;
     1858    gLog << "filenameMCup = " << filenameMCup << endl;
     1859
     1860    //-----------------------------------------------------------
     1861
     1862    TString fHilName    = "MHillas";
     1863    TString fHilNameExt = "MHillasExt";
     1864    TString fHilNameSrc = "MHillasSrc";
     1865    TString fImgParName = "MNewImagePar";
     1866
     1867    //===========================================================
     1868    //
     1869    // Optimization of energy estimator
     1870    //
     1871
     1872    TString inpath("");
     1873    TString outpath("");
     1874    Int_t howMany = 2000;
     1875    CT1EEst(inpath,   filenameOpt,   outpath, energyParName,
     1876            fHilName, fHilNameSrc,   fhadronnessName,
     1877            howMany,  maxhadronness, maxalpha, maxdist);
     1878
     1879    //-----------------------------------------------------------
     1880    //
     1881    // Read in parameters of energy estimator
     1882    //
     1883    gLog << "================================================================"
     1884         << endl;
     1885    gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '"
     1886         << energyParName << "'" << endl;
     1887    TFile enparam(energyParName);
     1888    MMcEnergyEst mcest("MMcEnergyEst");
     1889    mcest.Read("MMcEnergyEst");
     1890    enparam.Close();
     1891
     1892    TArrayD parA(5);
     1893    TArrayD parB(7);
     1894    for (Int_t i=0; i<parA.GetSize(); i++)
     1895      parA[i] = mcest.GetCoeff(i);
     1896    for (Int_t i=0; i<parB.GetSize(); i++)
     1897      parB[i] = mcest.GetCoeff( i+parA.GetSize() );
     1898
     1899
     1900   if (WESTUP)
     1901   {
     1902    //==========   start update   ============================================
     1903    //
     1904    // Update ON, OFF and MC root files with the estimated energy
     1905
     1906    //---------------------------------------------------
     1907    // Update OFF data
     1908    //
     1909    gLog << "============================================================"
     1910         << endl;
     1911    gLog << "Macro CT1Analysis.C : update file '" << filenameOFF
     1912         << "'" << endl;
     1913
     1914    MTaskList tlistoff;
     1915    MParList plistoff;
     1916
     1917
     1918    // geometry is needed in  MHHillas... classes
     1919    MGeomCam *fGeom =
     1920             (MGeomCam*)plistoff->FindCreateObj("MGeomCamCT1", "MGeomCam");
     1921
     1922    //-------------------------------------------
     1923    // create the tasks which should be executed
     1924    //
     1925
     1926    MReadMarsFile read("Events", filenameOFF);
     1927    read.DisableAutoScheme();
     1928
     1929    //---------------------------
     1930    // calculate estimated energy
     1931
     1932    MEnergyEstParam eest2(fHilName);
     1933    eest2.Add(fHilNameSrc);
     1934
     1935    eest2.SetCoeffA(parA);
     1936    eest2.SetCoeffB(parB);
     1937
    17391938    //.......................................................................
    1740     // calculate hadronnes for method of RANDOM FOREST
    1741 
    1742     TString hadRFName = "HadRF";
    1743     MRanForestCalc rfcalc;
    1744     rfcalc.SetHadronnessName(hadRFName);
    1745 
    1746     //.......................................................................
    1747     // calculation of hadroness for the supercuts
    1748     // (=0.25 if fullfilled, =0.75 otherwise)
    1749 
    1750     TString hadSCName = "HadSC";
    1751     MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc);
    1752     sccalc.SetHadronnessName(hadSCName);
    1753 
    1754 
    1755     //.......................................................................
    1756 
    1757       //MWriteRootFile write(outNameImage, "UPDATE");
    1758       MWriteRootFile write(outNameImage, "RECREATE");
     1939
     1940      MWriteRootFile write(filenameOFFup);
    17591941
    17601942      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    17691951      write.AddContainer("MNewImagePar",  "Events");
    17701952
     1953      //write.AddContainer("HadNN",         "Events");
     1954      write.AddContainer("HadSC",         "Events");
    17711955      write.AddContainer("HadRF",         "Events");
    1772       write.AddContainer("HadSC",         "Events");
     1956
     1957      write.AddContainer("MEnergyEst",    "Events");
    17731958
    17741959    //-----------------------------------------------------------------
     1960
     1961    MFCT1SelFinal selfinal(fHilNameSrc);
     1962    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
     1963    selfinal.SetHadronnessName(fhadronnessName);
     1964    MContinue contfinal(&selfinal);
     1965
     1966
     1967    //*****************************
     1968    // entries in MParList
     1969
     1970    plistoff.AddToList(&tlistoff);
     1971    InitBinnings(&plistoff);
     1972
     1973
     1974    //*****************************
     1975    // entries in MTaskList
     1976   
     1977    tlistoff.AddToList(&read);
     1978    tlistoff.AddToList(&eest2);
     1979    tlistoff.AddToList(&write);
     1980    tlistoff.AddToList(&contfinal);
     1981
     1982    //*****************************
     1983
     1984    //-------------------------------------------
     1985    // Execute event loop
     1986    //
     1987    MProgressBar bar;
     1988    MEvtLoop evtloop;
     1989    evtloop.SetParList(&plistoff);
     1990    evtloop.SetProgressBar(&bar);
     1991
     1992    Int_t maxevents = -1;
     1993    //Int_t maxevents = 1000;
     1994    if ( !evtloop.Eventloop(maxevents) )
     1995        return;
     1996
     1997    tlistoff.PrintStatistics(0, kTRUE);
     1998    DeleteBinnings(&plistoff);
     1999
     2000    //---------------------------------------------------
     2001
     2002    //---------------------------------------------------
     2003    // Update ON data
     2004    //
     2005    gLog << "============================================================"
     2006         << endl;
     2007    gLog << "Macro CT1Analysis.C : update file '" << filenameON
     2008         << "'" << endl;
     2009
     2010    MTaskList tliston;
     2011    MParList pliston;
     2012
     2013
    17752014    // geometry is needed in  MHHillas... classes
    17762015    MGeomCam *fGeom =
    17772016             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    17782017
    1779 
    1780     Float_t maxhadronness =  0.40;
    1781     Float_t maxalpha      =  20.0;
    1782     Float_t maxdist       =  10.0;
    1783 
    1784     MFCT1SelFinal selfinalgh(fHilNameSrc);
    1785     selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    1786     selfinalgh.SetHadronnessName(hadRFName);
    1787     selfinalgh.SetName("SelFinalgh");
    1788     MContinue contfinalgh(&selfinalgh);
    1789     contfinalgh.SetName("ContSelFinalgh");
    1790 
    1791     MFillH fillranfor("MHRanForest");
    1792     fillranfor.SetName("HRanForest");
    1793 
    1794     MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
    1795     fillhadrf.SetName("HhadRF");
    1796     MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
    1797     fillhadsc.SetName("HhadSC");
    1798 
    1799 
    1800     MFCT1SelFinal selfinal(fHilNameSrc);
    1801     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    1802     selfinal.SetHadronnessName(hadRFName);
    1803     selfinal.SetName("SelFinal");
    1804     MContinue contfinal(&selfinal);
    1805     contfinal.SetName("ContSelFinal");
    1806 
    1807 
    1808     MFillH hfill1("MHHillas",    fHilName);
    1809     hfill1.SetName("HHillas");
    1810 
    1811     MFillH hfill2("MHStarMap",   fHilName);
    1812     hfill2.SetName("HStarMap");
    1813 
    1814     MFillH hfill3("MHHillasExt",    fHilNameSrc);
    1815     hfill3.SetName("HHillasExt");
    1816    
    1817     MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    1818     hfill4.SetName("HHillasSrc");   
    1819 
    1820     MFillH hfill5("MHNewImagePar", fImgParName);
    1821     hfill5.SetName("HNewImagePar");
    1822 
    1823     //*****************************
    1824     // entries in MParList
    1825 
    1826     pliston.AddToList(&tliston);
    1827     InitBinnings(&pliston);
    1828 
    1829     pliston.AddToList(fRanForest);
    1830     pliston.AddToList(fRanTree);
    1831 
    1832     //*****************************
    1833     // entries in MTaskList
    1834    
    1835     tliston.AddToList(&read);
    1836 
    1837     tliston.AddToList(&rfcalc);
    1838     tliston.AddToList(&sccalc);
    1839     tliston.AddToList(&fillranfor);
    1840     tliston.AddToList(&fillhadrf);
    1841     tliston.AddToList(&fillhadsc);
    1842 
    1843     tliston.AddToList(&hfill1);
    1844     tliston.AddToList(&hfill2);
    1845     tliston.AddToList(&hfill3);
    1846     tliston.AddToList(&hfill4);
    1847     tliston.AddToList(&hfill5);
    1848 
    1849     tliston.AddToList(&write);
    1850 
    1851     tliston.AddToList(&contfinalgh);
    1852     tliston.AddToList(&contfinal);
    1853 
    1854     //*****************************
    1855 
    1856     //-------------------------------------------
    1857     // Execute event loop
    1858     //
    1859     MProgressBar bar;
    1860     MEvtLoop evtloop;
    1861     evtloop.SetParList(&pliston);
    1862     evtloop.SetProgressBar(&bar);
    1863 
    1864     Int_t maxevents = -1;
    1865     if ( !evtloop.Eventloop(maxevents) )
    1866         return;
    1867 
    1868     tliston.PrintStatistics(0, kTRUE);
    1869 
    1870 
    1871     //-------------------------------------------
    1872     // Display the histograms
    1873     //
    1874     pliston.FindObject("MHRanForest")->DrawClone();
    1875     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    1876     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    1877 
    1878     pliston.FindObject("MHHillas")->DrawClone();
    1879     pliston.FindObject("MHHillasExt")->DrawClone();
    1880     pliston.FindObject("MHHillasSrc")->DrawClone();
    1881     pliston.FindObject("MHNewImagePar")->DrawClone();
    1882     pliston.FindObject("MHStarMap")->DrawClone();
    1883 
    1884     DeleteBinnings(&pliston);
    1885   }
    1886 
    1887     gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
    1888     gLog << "=======================================================" << endl;
    1889  }
    1890   //---------------------------------------------------------------------
    1891 
    1892 
    1893 
    1894   //---------------------------------------------------------------------
    1895   // Job C 
    1896   //======
    1897 
    1898     //  - read ON1 and MC1 data files 
    1899     //    which should have been updated to contain the hadronnesses
    1900     //    for the method of NEAREST NEIGHBORS and for the SUOERCUTS
    1901     //  - produce Neyman-Pearson plots
    1902  
    1903  if (JobC)
    1904  {
    1905     gLog << "=====================================================" << endl;
    1906     gLog << "Macro CT1Analysis : Start of Job C" << endl;
    1907 
    1908     gLog << "" << endl;
    1909     gLog << "Macro CT1Analysis : JobC = " << JobC  << endl;
    1910 
    1911 
    1912     // name of input data file
    1913     TString filenameData = outPath;
    1914     filenameData += "OFF";
    1915     filenameData += "2.root";
    1916     gLog << "filenameData = " << filenameData << endl;
    1917 
    1918     // name of input MC file
    1919     TString filenameMC = outPath;
    1920     filenameMC += "MC";
    1921     filenameMC += "2.root";
    1922     gLog << "filenameMC   = " << filenameMC   << endl;
    1923 
    1924 
    1925     //-----------------------------------------------------------------
    1926 
    1927     MTaskList tliston;
    1928     MParList pliston;
    1929 
    1930 
    1931     // geometry is needed in  MHHillas... classes
    1932     MGeomCam *fGeom =
    1933              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    1934 
    19352018    //-------------------------------------------
    19362019    // create the tasks which should be executed
    19372020    //
    19382021
    1939     MReadMarsFile read("Events", filenameMC);
    1940     read.AddFile(filenameData);
    1941     read.DisableAutoScheme();
    1942 
    1943 
    1944     //.......................................................................
    1945     // names of hadronness containers
    1946 
    1947     TString hadNNName = "HadNN";
    1948     TString hadSCName = "HadSC";
    1949     TString hadRFName = "HadRF";
    1950 
    1951     //.......................................................................
    1952 
    1953 
    1954     TString fHilName    = "MHillas";
    1955     TString fHilNameExt = "MHillasExt";
    1956     TString fHilNameSrc = "MHillasSrc";
    1957     TString fImgParName = "MNewImagePar";
    1958 
    1959     Float_t maxhadronness =  0.40;
    1960     Float_t maxalpha      =  20.0;
    1961     Float_t maxdist       =  10.0;
    1962 
    1963     MFCT1SelFinal selfinalgh(fHilNameSrc);
    1964     selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    1965     selfinalgh.SetHadronnessName(hadRFName);
    1966     selfinalgh.SetName("SelFinalgh");
    1967     MContinue contfinalgh(&selfinalgh);
    1968     contfinalgh.SetName("ContSelFinalgh");
    1969 
    1970     //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
    1971     //fillhadnn.SetName("HhadNN");
    1972     MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
    1973     fillhadsc.SetName("HhadSC");
    1974     MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
    1975     fillhadrf.SetName("HhadRF");
    1976 
    1977     MFCT1SelFinal selfinal(fHilNameSrc);
    1978     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    1979     selfinal.SetHadronnessName(hadRFName);
    1980     selfinal.SetName("SelFinal");
    1981     MContinue contfinal(&selfinal);
    1982     contfinal.SetName("ContSelFinal");
    1983 
    1984 
    1985     MFillH hfill1("MHHillas",    fHilName);
    1986     hfill1.SetName("HHillas");
    1987 
    1988     MFillH hfill2("MHStarMap",   fHilName);
    1989     hfill2.SetName("HStarMap");
    1990 
    1991     MFillH hfill3("MHHillasExt",    fHilNameSrc);
    1992     hfill3.SetName("HHillasExt");
    1993    
    1994     MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    1995     hfill4.SetName("HHillasSrc");   
    1996 
    1997     MFillH hfill5("MHNewImagePar", fImgParName);
    1998     hfill5.SetName("HNewImagePar");
    1999 
    2000 
    2001     //*****************************
    2002     // entries in MParList
    2003 
    2004     pliston.AddToList(&tliston);
    2005     InitBinnings(&pliston);
    2006 
    2007 
    2008     //*****************************
    2009     // entries in MTaskList
    2010    
    2011     tliston.AddToList(&read);
    2012 
    2013     //tliston.AddToList(&fillhadnn);
    2014     tliston.AddToList(&fillhadsc);
    2015     tliston.AddToList(&fillhadrf);
    2016    
    2017     tliston.AddToList(&contfinalgh);
    2018     tliston.AddToList(&hfill1);
    2019     tliston.AddToList(&hfill2);
    2020     tliston.AddToList(&hfill3);
    2021     tliston.AddToList(&hfill4);
    2022     tliston.AddToList(&hfill5);
    2023 
    2024     tliston.AddToList(&contfinal);
    2025 
    2026     //*****************************
    2027 
    2028     //-------------------------------------------
    2029     // Execute event loop
    2030     //
    2031     MProgressBar bar;
    2032     MEvtLoop evtloop;
    2033     evtloop.SetParList(&pliston);
    2034     evtloop.SetProgressBar(&bar);
    2035 
    2036     Int_t maxevents = -1;
    2037     //Int_t maxevents = 35000;
    2038     if ( !evtloop.Eventloop(maxevents) )
    2039         return;
    2040 
    2041     tliston.PrintStatistics(0, kTRUE);
    2042 
    2043 
    2044     //-------------------------------------------
    2045     // Display the histograms
    2046     //
    2047 
    2048     //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
    2049     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    2050     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    2051 
    2052     pliston.FindObject("MHHillas")->DrawClone();
    2053     pliston.FindObject("MHHillasExt")->DrawClone();
    2054     pliston.FindObject("MHHillasSrc")->DrawClone();
    2055     pliston.FindObject("MHNewImagePar")->DrawClone();
    2056     pliston.FindObject("MHStarMap")->DrawClone();
    2057 
    2058     DeleteBinnings(&pliston);
    2059 
    2060     gLog << "Macro CT1Analysis : End of Job C" << endl;
    2061     gLog << "===================================================" << endl;
    2062  }
    2063 
    2064 
    2065   //---------------------------------------------------------------------
    2066   // Job D
    2067   //======
    2068 
    2069     //  - select g/h separation method XX
    2070     //  - read ON2 (or MC2) root file
    2071     //  - apply cuts in hadronness
    2072     //  - make plots
    2073 
    2074 
    2075  if (JobD)
    2076  {
    2077     gLog << "=====================================================" << endl;
    2078     gLog << "Macro CT1Analysis : Start of Job D" << endl;
    2079 
    2080     gLog << "" << endl;
    2081     gLog << "Macro CT1Analysis : JobD = "
    2082          << JobD  << endl;
    2083 
    2084     // type of data to be analysed
    2085     TString typeData = "ON";
    2086     //TString typeData = "OFF";
    2087     //TString typeData = "MC";
    2088     gLog << "typeData = " << typeData << endl;
    2089 
    2090     TString ext      = "2.root";
    2091 
    2092 
    2093     //------------------------------
    2094     // selection of g/h separation method
    2095     // and definition of final selections
    2096 
    2097     //TString XX("NN");
    2098     //TString XX("SC");
    2099     TString XX("RF");
    2100     TString fhadronnessName("Had");
    2101     fhadronnessName += XX;
    2102     gLog << "fhadronnessName = " << fhadronnessName << endl;
    2103 
    2104     // maximum values of the hadronness, |ALPHA| and DIST
    2105     Float_t maxhadronness   = 0.30;
    2106     Float_t maxalpha        = 20.0;
    2107     Float_t maxdist         = 10.0;
    2108     gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
    2109          << maxhadronness << ",  " << maxalpha << ",  "
    2110          << maxdist << endl;
    2111 
    2112 
    2113     //------------------------------
    2114     // name of data file to be analysed
    2115     TString filenameData(outPath);
    2116     filenameData += typeData;
    2117     filenameData += ext;
    2118     gLog << "filenameData = " << filenameData << endl;
    2119 
    2120 
    2121 
    2122     //*************************************************************************
    2123     //
    2124     // Analyse the data
    2125     //
    2126 
    2127     MTaskList tliston;
    2128     MParList pliston;
    2129 
    2130     // geometry is needed in  MHHillas... classes
    2131     MGeomCam *fGeom =
    2132              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    2133 
    2134 
    2135     TString fHilName    = "MHillas";
    2136     TString fHilNameExt = "MHillasExt";
    2137     TString fHilNameSrc = "MHillasSrc";
    2138     TString fImgParName = "MNewImagePar";
    2139 
    2140     //-------------------------------------------
    2141     // create the tasks which should be executed
    2142     //
    2143 
    2144     MReadMarsFile read("Events", filenameData);
    2145     read.DisableAutoScheme();
    2146 
    2147 
    2148     //-----------------------------------------------------------------
    2149     // geometry is needed in  MHHillas... classes
    2150     MGeomCam *fGeom =
    2151              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    2152 
    2153     MFCT1SelFinal selfinalgh(fHilNameSrc);
    2154     selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    2155     selfinalgh.SetHadronnessName(fhadronnessName);
    2156     selfinalgh.SetName("SelFinalgh");
    2157     MContinue contfinalgh(&selfinalgh);
    2158     contfinalgh.SetName("ContSelFinalgh");
    2159 
    2160     //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
    2161     //fillhadnn.SetName("HhadNN");
    2162     MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
    2163     fillhadsc.SetName("HhadSC");
    2164     MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
    2165     fillhadrf.SetName("HhadRF");
    2166 
    2167 
    2168     MFillH hfill1("MHHillas",    fHilName);
    2169     hfill1.SetName("HHillas");
    2170 
    2171     MFillH hfill2("MHStarMap",   fHilName);
    2172     hfill2.SetName("HStarMap");
    2173 
    2174     MFillH hfill3("MHHillasExt",   fHilNameSrc);
    2175     hfill3.SetName("HHillasExt");   
    2176 
    2177     MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    2178     hfill4.SetName("HHillasSrc");   
    2179 
    2180     MFillH hfill5("MHNewImagePar", fImgParName);
    2181     hfill5.SetName("HNewImagePar");   
    2182 
    2183     MFCT1SelFinal selfinal(fHilNameSrc);
    2184     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    2185     selfinal.SetHadronnessName(fhadronnessName);
    2186     selfinal.SetName("SelFinal");
    2187     MContinue contfinal(&selfinal);
    2188     contfinal.SetName("ContSelFinal");
    2189 
    2190 
    2191     //*****************************
    2192     // entries in MParList
    2193 
    2194     pliston.AddToList(&tliston);
    2195     InitBinnings(&pliston);
    2196 
    2197 
    2198     //*****************************
    2199     // entries in MTaskList
    2200    
    2201     tliston.AddToList(&read);
    2202 
    2203     tliston.AddToList(&contfinalgh);
    2204 
    2205     //tliston.AddToList(&fillhadnn);
    2206     tliston.AddToList(&fillhadsc);
    2207     tliston.AddToList(&fillhadrf);
    2208 
    2209     tliston.AddToList(&hfill1);
    2210     tliston.AddToList(&hfill2);
    2211     tliston.AddToList(&hfill3);
    2212     tliston.AddToList(&hfill4);
    2213     tliston.AddToList(&hfill5);
    2214 
    2215     tliston.AddToList(&contfinal);
    2216 
    2217     //*****************************
    2218 
    2219     //-------------------------------------------
    2220     // Execute event loop
    2221     //
    2222     MProgressBar bar;
    2223     MEvtLoop evtloop;
    2224     evtloop.SetParList(&pliston);
    2225     evtloop.SetProgressBar(&bar);
    2226 
    2227     Int_t maxevents = -1;
    2228     //Int_t maxevents = 10000;
    2229     if ( !evtloop.Eventloop(maxevents) )
    2230         return;
    2231 
    2232     tliston.PrintStatistics(0, kTRUE);
    2233 
    2234 
    2235 
    2236     //-------------------------------------------
    2237     // Display the histograms
    2238     //
    2239 
    2240     //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
    2241     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    2242     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    2243 
    2244     pliston.FindObject("MHHillas")->DrawClone();
    2245     pliston.FindObject("MHHillasExt")->DrawClone();
    2246     pliston.FindObject("MHHillasSrc")->DrawClone();
    2247     pliston.FindObject("MHNewImagePar")->DrawClone();
    2248     pliston.FindObject("MHStarMap")->DrawClone();
    2249 
    2250     //-------------------------------------------
    2251     // fit alpha distribution to get the number of excess events
    2252     //
    2253 
    2254     MHHillasSrc* hillasSrc =
    2255       (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
    2256     TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
    2257  
    2258     MHOnSubtraction onsub;
    2259     onsub.Calc(alphaHist, &pliston, kTRUE, 13.1);
    2260     //-------------------------------------------
    2261 
    2262     DeleteBinnings(&pliston);
    2263 
    2264     gLog << "Macro CT1Analysis : End of Job D" << endl;
    2265     gLog << "=======================================================" << endl;
    2266  }
    2267   //---------------------------------------------------------------------
    2268 
    2269 
    2270   //---------------------------------------------------------------------
    2271   // Job E_EST_UP
    2272   //============
    2273 
    2274     //  - read MC2.root file
    2275     //  - select g/h separation method XX
    2276     //  - optimize energy estimator for events passing final cuts
    2277     //  - write parameters of energy estimator onto file "energyest_XX.root"
    2278     //
    2279     //  - read ON2.root and MC2.root files
    2280     //  - update input root file with the estimated energy
    2281     //    (ON_XX2.root, MC_XX2.root)
    2282 
    2283 
    2284  if (JobE_EST_UP)
    2285  {
    2286     gLog << "=====================================================" << endl;
    2287     gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl;
    2288 
    2289     gLog << "" << endl;
    2290     gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = "
    2291          << JobE_EST_UP  << ",  " << WESTUP << endl;
    2292 
    2293 
    2294     TString typeON  = "ON";
    2295     TString typeOFF = "OFF";
    2296     TString typeMC  = "MC";
    2297     TString ext    = "2.root";
    2298     TString extout = "3.root";
    2299 
    2300     //------------------------------
    2301     // name of MC file to be used for optimizing the energy estimator
    2302     TString filenameOpt(outPath);
    2303     filenameOpt += typeMC;
    2304     filenameOpt += ext;
    2305     gLog << "filenameOpt = " << filenameOpt << endl;
    2306 
    2307     //------------------------------
    2308     // selection of g/h separation method
    2309     // and definition of final selections
    2310 
    2311     //TString XX("NN");
    2312     //TString XX("SC");
    2313     TString XX("RF");
    2314     TString fhadronnessName("Had");
    2315     fhadronnessName += XX;
    2316     gLog << "fhadronnessName = " << fhadronnessName << endl;
    2317 
    2318     // maximum values of the hadronness, |alpha| and dist
    2319     Float_t maxhadronness   = 0.40;
    2320     Float_t maxalpha        = 20.0;
    2321     Float_t maxdist         = 10.0;
    2322     gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
    2323          << maxhadronness << ",  " << maxalpha << ",  "
    2324          << maxdist << endl;
    2325 
    2326     // name of file containing the parameters of the energy estimator
    2327     TString energyParName(outPath);
    2328     energyParName += "energyest_";
    2329     energyParName += XX;
    2330     energyParName += ".root";
    2331     gLog << "energyParName = " << energyParName << endl;
    2332 
    2333 
    2334     //------------------------------
    2335     // name of ON file to be updated
    2336     TString filenameON(outPath);
    2337     filenameON += typeON;
    2338     filenameON += ext;
    2339     gLog << "filenameON = " << filenameON << endl;
    2340 
    2341     // name of OFF file to be updated
    2342     TString filenameOFF(outPath);
    2343     filenameOFF += typeOFF;
    2344     filenameOFF += ext;
    2345     gLog << "filenameOFF = " << filenameOFF << endl;
    2346 
    2347     // name of MC file to be updated
    2348     TString filenameMC(outPath);
    2349     filenameMC += typeMC;
    2350     filenameMC += ext;
    2351     gLog << "filenameMC = " << filenameMC << endl;
    2352 
    2353     //------------------------------
    2354     // name of updated ON file
    2355     TString filenameONup(outPath);
    2356     filenameONup += typeON;
    2357     filenameONup += "_";
    2358     filenameONup += XX;
    2359     filenameONup += extout;
    2360     gLog << "filenameONup = " << filenameONup << endl;
    2361 
    2362     // name of updated OFF file
    2363     TString filenameOFFup(outPath);
    2364     filenameOFFup += typeOFF;
    2365     filenameOFFup += "_";
    2366     filenameOFFup += XX;
    2367     filenameOFFup += extout;
    2368     gLog << "filenameOFFup = " << filenameOFFup << endl;
    2369 
    2370     // name of updated MC file
    2371     TString filenameMCup(outPath);
    2372     filenameMCup += typeMC;
    2373     filenameMCup += "_";
    2374     filenameMCup += XX;
    2375     filenameMCup += extout;
    2376     gLog << "filenameMCup = " << filenameMCup << endl;
    2377 
    2378     //-----------------------------------------------------------
    2379 
    2380     TString fHilName    = "MHillas";
    2381     TString fHilNameExt = "MHillasExt";
    2382     TString fHilNameSrc = "MHillasSrc";
    2383     TString fImgParName = "MNewImagePar";
    2384 
    2385     //===========================================================
    2386     //
    2387     // Optimization of energy estimator
    2388     //
    2389 
    2390     TString inpath("");
    2391     TString outpath("");
    2392     Int_t howMany = 2000;
    2393     CT1EEst(inpath,   filenameOpt,   outpath, energyParName,
    2394             fHilName, fHilNameSrc,   fhadronnessName,
    2395             howMany,  maxhadronness, maxalpha, maxdist);
    2396 
    2397     //-----------------------------------------------------------
    2398     //
    2399     // Read in parameters of energy estimator
    2400     //
    2401     gLog << "================================================================"
    2402          << endl;
    2403     gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '"
    2404          << energyParName << "'" << endl;
    2405     TFile enparam(energyParName);
    2406     MMcEnergyEst mcest("MMcEnergyEst");
    2407     mcest.Read("MMcEnergyEst");
    2408     enparam.Close();
    2409 
    2410     TArrayD parA(5);
    2411     TArrayD parB(7);
    2412     for (Int_t i=0; i<parA.GetSize(); i++)
    2413       parA[i] = mcest.GetCoeff(i);
    2414     for (Int_t i=0; i<parB.GetSize(); i++)
    2415       parB[i] = mcest.GetCoeff( i+parA.GetSize() );
    2416 
    2417 
    2418    if (WESTUP)
    2419    {
    2420     //==========   start update   ============================================
    2421     //
    2422     // Update ON, OFF and MC root files with the estimated energy
    2423 
    2424     //---------------------------------------------------
    2425     // Update OFF data
    2426     //
    2427     gLog << "============================================================"
    2428          << endl;
    2429     gLog << "Macro CT1Analysis.C : update file '" << filenameOFF
    2430          << "'" << endl;
    2431 
    2432     MTaskList tlistoff;
    2433     MParList plistoff;
    2434 
    2435 
    2436     // geometry is needed in  MHHillas... classes
    2437     MGeomCam *fGeom =
    2438              (MGeomCam*)plistoff->FindCreateObj("MGeomCamCT1", "MGeomCam");
    2439 
    2440     //-------------------------------------------
    2441     // create the tasks which should be executed
    2442     //
    2443 
    2444     MReadMarsFile read("Events", filenameOFF);
     2022    MReadMarsFile read("Events", filenameON);
    24452023    read.DisableAutoScheme();
    24462024
    … …  
    24562034    //.......................................................................
    24572035
    2458       MWriteRootFile write(filenameOFFup);
     2036      MWriteRootFile write(filenameONup);
    24592037
    24602038      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    24862064    // entries in MParList
    24872065
    2488     plistoff.AddToList(&tlistoff);
    2489     InitBinnings(&plistoff);
     2066    pliston.AddToList(&tliston);
     2067    InitBinnings(&pliston);
    24902068
    24912069
    … …  
    24932071    // entries in MTaskList
    24942072   
    2495     tlistoff.AddToList(&read);
    2496     tlistoff.AddToList(&eest2);
    2497     tlistoff.AddToList(&write);
    2498     tlistoff.AddToList(&contfinal);
     2073    tliston.AddToList(&read);
     2074    tliston.AddToList(&eest2);
     2075    tliston.AddToList(&write);
     2076    tliston.AddToList(&contfinal);
    24992077
    25002078    //*****************************
    … …  
    25052083    MProgressBar bar;
    25062084    MEvtLoop evtloop;
    2507     evtloop.SetParList(&plistoff);
     2085    evtloop.SetParList(&pliston);
    25082086    evtloop.SetProgressBar(&bar);
    25092087
    … …  
    25132091        return;
    25142092
    2515     tlistoff.PrintStatistics(0, kTRUE);
    2516     DeleteBinnings(&plistoff);
     2093    tliston.PrintStatistics(0, kTRUE);
     2094    DeleteBinnings(&pliston);
    25172095
    25182096    //---------------------------------------------------
    25192097
    25202098    //---------------------------------------------------
    2521     // Update ON data
     2099    // Update MC data
    25222100    //
    25232101    gLog << "============================================================"
    25242102         << endl;
    2525     gLog << "Macro CT1Analysis.C : update file '" << filenameON
     2103    gLog << "Macro CT1Analysis.C : update file '" << filenameMC
    25262104         << "'" << endl;
    25272105
    2528     MTaskList tliston;
    2529     MParList pliston;
    2530 
    2531 
    2532     // geometry is needed in  MHHillas... classes
    2533     MGeomCam *fGeom =
    2534              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     2106    MTaskList tlistmc;
     2107    MParList plistmc;
    25352108
    25362109    //-------------------------------------------
    … …  
    25382111    //
    25392112
    2540     MReadMarsFile read("Events", filenameON);
     2113    MReadMarsFile read("Events", filenameMC);
    25412114    read.DisableAutoScheme();
    25422115
    … …  
    25522125    //.......................................................................
    25532126
    2554       MWriteRootFile write(filenameONup);
     2127      MWriteRootFile write(filenameMCup);
    25552128
    25562129      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    25822155    // entries in MParList
    25832156
    2584     pliston.AddToList(&tliston);
    2585     InitBinnings(&pliston);
     2157    plistmc.AddToList(&tlistmc);
     2158    InitBinnings(&plistmc);
    25862159
    25872160
    … …  
    25892162    // entries in MTaskList
    25902163   
    2591     tliston.AddToList(&read);
    2592     tliston.AddToList(&eest2);
    2593     tliston.AddToList(&write);
    2594     tliston.AddToList(&contfinal);
     2164    tlistmc.AddToList(&read);
     2165    tlistmc.AddToList(&eest2);
     2166    tlistmc.AddToList(&write);
     2167    tlistmc.AddToList(&contfinal);
    25952168
    25962169    //*****************************
    … …  
    26012174    MProgressBar bar;
    26022175    MEvtLoop evtloop;
    2603     evtloop.SetParList(&pliston);
     2176    evtloop.SetParList(&plistmc);
    26042177    evtloop.SetProgressBar(&bar);
    26052178
    … …  
    26092182        return;
    26102183
    2611     tliston.PrintStatistics(0, kTRUE);
    2612     DeleteBinnings(&pliston);
    2613 
    2614     //---------------------------------------------------
    2615 
    2616     //---------------------------------------------------
    2617     // Update MC data
    2618     //
    2619     gLog << "============================================================"
    2620          << endl;
    2621     gLog << "Macro CT1Analysis.C : update file '" << filenameMC
    2622          << "'" << endl;
    2623 
    2624     MTaskList tlistmc;
    2625     MParList plistmc;
     2184    tlistmc.PrintStatistics(0, kTRUE);
     2185    DeleteBinnings(&plistmc);
     2186
     2187
     2188    //==========   end update   ============================================
     2189   }
     2190   
     2191    enparam.Close();
     2192
     2193    gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl;
     2194    gLog << "=======================================================" << endl;
     2195 }
     2196  //---------------------------------------------------------------------
     2197
     2198
     2199  //---------------------------------------------------------------------
     2200  // Job F_XX
     2201  //=========
     2202
     2203    //  - select g/h separation method XX
     2204    //  - read MC_XX2.root file
     2205    //  - calculate eff. collection area
     2206    //  - read ON_XX2.root file
     2207    //  - apply final cuts
     2208    //  - calculate flux
     2209    //  - write root file for ON data after final cuts (ON_XX3.root))
     2210
     2211
     2212 if (JobF_XX)
     2213 {
     2214    gLog << "=====================================================" << endl;
     2215    gLog << "Macro CT1Analysis : Start of Job F_XX" << endl;
     2216
     2217    gLog << "" << endl;
     2218    gLog << "Macro CT1Analysis : JobF_XX, WXX = "
     2219         << JobF_XX  << ",  " << WXX << endl;
     2220
     2221    // type of data to be analysed
     2222    //TString typeData = "ON";
     2223    //TString typeData = "OFF";
     2224    TString typeData = "MC";
     2225    gLog << "typeData = " << typeData << endl;
     2226
     2227    TString typeMC   = "MC";
     2228    TString ext      = "3.root";
     2229    TString extout   = "4.root";
     2230
     2231    //------------------------------
     2232    // selection of g/h separation method
     2233    // and definition of final selections
     2234
     2235    //TString XX("NN");
     2236    //TString XX("SC");
     2237    TString XX("RF");
     2238    TString fhadronnessName("Had");
     2239    fhadronnessName += XX;
     2240    gLog << "fhadronnessName = " << fhadronnessName << endl;
     2241
     2242    // maximum values of the hadronness, |ALPHA| and DIST
     2243    Float_t maxhadronness   = 0.40;
     2244    Float_t maxalpha        = 20.0;
     2245    Float_t maxdist         = 10.0;
     2246    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
     2247         << maxhadronness << ",  " << maxalpha << ",  "
     2248         << maxdist << endl;
     2249
     2250
     2251    //------------------------------
     2252    // name of MC file to be used for calculating the eff. collection areas
     2253    TString filenameArea(outPath);
     2254    filenameArea += typeMC;
     2255    filenameArea += "_";
     2256    filenameArea += XX;
     2257    filenameArea += ext;
     2258    gLog << "filenameArea = " << filenameArea << endl;
     2259
     2260    //------------------------------
     2261    // name of file containing the eff. collection areas
     2262    TString collareaName(outPath);
     2263    collareaName += "area_";
     2264    collareaName += XX;
     2265    collareaName += ".root";
     2266    gLog << "collareaName = " << collareaName << endl;
     2267
     2268    //------------------------------
     2269    // name of data file to be analysed
     2270    TString filenameData(outPath);
     2271    filenameData += typeData;
     2272    filenameData += "_";
     2273    filenameData += XX;
     2274    filenameData += ext;
     2275    gLog << "filenameData = " << filenameData << endl;
     2276
     2277    //------------------------------
     2278    // name of output data file (after the final cuts)
     2279    TString filenameDataout(outPath);
     2280    filenameDataout += typeData;
     2281    filenameDataout += "_";
     2282    filenameDataout += XX;
     2283    filenameDataout += extout;
     2284    gLog << "filenameDataout = " << filenameDataout << endl;
     2285
     2286
     2287    //====================================================================
     2288    gLog << "-----------------------------------------------" << endl;
     2289    gLog << "Start calculation of effective collection areas" << endl;
     2290    MParList  parlist;
     2291    MTaskList tasklist;
     2292
     2293    //---------------------------------------
     2294    // Setup the tasks to be executed
     2295    //
     2296    MReadMarsFile reader("Events", filenameArea);
     2297    reader.DisableAutoScheme();
     2298
     2299    MFCT1SelFinal cuthadrons;
     2300    cuthadrons.SetHadronnessName(fhadronnessName);
     2301    cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
     2302
     2303    MContinue conthadrons(&cuthadrons);
     2304
     2305    //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
     2306    //MHMcCT1CollectionArea* collarea;
     2307
     2308    MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
     2309    filler.SetName("CollectionArea");
     2310
     2311    //********************************
     2312    // entries in MParList
     2313
     2314    parlist.AddToList(&tasklist);
     2315    InitBinnings(&parlist);
     2316    //parlist.AddToList(collarea);
     2317
     2318    //********************************
     2319    // entries in MTaskList
     2320
     2321    tasklist.AddToList(&reader);   
     2322    tasklist.AddToList(&conthadrons);
     2323    tasklist.AddToList(&filler);
     2324
     2325    //********************************
     2326
     2327    //-----------------------------------------
     2328    // Execute event loop
     2329    //
     2330    MEvtLoop magic;
     2331    magic.SetParList(&parlist);
     2332
     2333    MProgressBar bar;
     2334    magic.SetProgressBar(&bar);
     2335    if (!magic.Eventloop())
     2336        return;
     2337
     2338    tasklist.PrintStatistics(0, kTRUE);
     2339
     2340    // Calculate effective collection areas
     2341    // and display the histograms
     2342    //
     2343
     2344    MHMcCT1CollectionArea *collarea =
     2345        (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
     2346
     2347    collarea->CalcEfficiency();
     2348    collarea->DrawClone("lego");
     2349
     2350    //---------------------------------------------
     2351    // Write histograms to a file
     2352    //
     2353
     2354    TFile f(collareaName, "RECREATE");
     2355    collarea->GetHist()->Write();
     2356    collarea->GetHAll()->Write();
     2357    collarea->GetHSel()->Write();
     2358    f.Close();
     2359
     2360    //delete collarea;
     2361
     2362    gLog << "Calculation of effective collection areas done" << endl;
     2363    gLog << "-----------------------------------------------" << endl;   
     2364    //------------------------------------------------------------------
     2365
     2366
     2367    //*************************************************************************
     2368    //
     2369    // Analyse the data
     2370    //
     2371
     2372    MTaskList tliston;
     2373    MParList pliston;
     2374
     2375    // geometry is needed in  MHHillas... classes
     2376    MGeomCam *fGeom =
     2377             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     2378
     2379    TString fHilName    = "MHillas";
     2380    TString fHilNameExt = "MHillasExt";
     2381    TString fHilNameSrc = "MHillasSrc";
     2382    TString fImgParName = "MNewImagePar";
    26262383
    26272384    //-------------------------------------------
    … …  
    26292386    //
    26302387
    2631     MReadMarsFile read("Events", filenameMC);
     2388    MReadMarsFile read("Events", filenameData);
    26322389    read.DisableAutoScheme();
    26332390
    2634     //---------------------------
    2635     // calculate estimated energy
    2636 
    2637     MEnergyEstParam eest2(fHilName);
    2638     eest2.Add(fHilNameSrc);
    2639 
    2640     eest2.SetCoeffA(parA);
    2641     eest2.SetCoeffB(parB);
    2642 
    26432391    //.......................................................................
    26442392
    2645       MWriteRootFile write(filenameMCup);
     2393
     2394      MWriteRootFile write(filenameDataout);
    26462395
    26472396      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    26622411      write.AddContainer("MEnergyEst",    "Events");
    26632412
    2664     //-----------------------------------------------------------------
    2665 
    2666     MFCT1SelFinal selfinal(fHilNameSrc);
    2667     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    2668     selfinal.SetHadronnessName(fhadronnessName);
    2669     MContinue contfinal(&selfinal);
    2670 
    2671 
    2672     //*****************************
    2673     // entries in MParList
    2674 
    2675     plistmc.AddToList(&tlistmc);
    2676     InitBinnings(&plistmc);
    2677 
    2678 
    2679     //*****************************
    2680     // entries in MTaskList
    2681    
    2682     tlistmc.AddToList(&read);
    2683     tlistmc.AddToList(&eest2);
    2684     tlistmc.AddToList(&write);
    2685     tlistmc.AddToList(&contfinal);
    2686 
    2687     //*****************************
    2688 
    2689     //-------------------------------------------
    2690     // Execute event loop
    2691     //
    2692     MProgressBar bar;
    2693     MEvtLoop evtloop;
    2694     evtloop.SetParList(&plistmc);
    2695     evtloop.SetProgressBar(&bar);
    2696 
    2697     Int_t maxevents = -1;
    2698     //Int_t maxevents = 1000;
    2699     if ( !evtloop.Eventloop(maxevents) )
    2700         return;
    2701 
    2702     tlistmc.PrintStatistics(0, kTRUE);
    2703     DeleteBinnings(&plistmc);
    2704 
    2705 
    2706     //==========   end update   ============================================
    2707    }
    2708    
    2709     enparam.Close();
    2710 
    2711     gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl;
    2712     gLog << "=======================================================" << endl;
    2713  }
    2714   //---------------------------------------------------------------------
    2715 
    2716 
    2717   //---------------------------------------------------------------------
    2718   // Job F_XX
    2719   //=========
    2720 
    2721     //  - select g/h separation method XX
    2722     //  - read MC_XX2.root file
    2723     //  - calculate eff. collection area
    2724     //  - read ON_XX2.root file
    2725     //  - apply final cuts
    2726     //  - calculate flux
    2727     //  - write root file for ON data after final cuts (ON_XX3.root))
    2728 
    2729 
    2730  if (JobF_XX)
    2731  {
    2732     gLog << "=====================================================" << endl;
    2733     gLog << "Macro CT1Analysis : Start of Job F_XX" << endl;
    2734 
    2735     gLog << "" << endl;
    2736     gLog << "Macro CT1Analysis : JobF_XX, WXX = "
    2737          << JobF_XX  << ",  " << WXX << endl;
    2738 
    2739     // type of data to be analysed
    2740     //TString typeData = "ON";
    2741     //TString typeData = "OFF";
    2742     TString typeData = "MC";
    2743     gLog << "typeData = " << typeData << endl;
    2744 
    2745     TString typeMC   = "MC";
    2746     TString ext      = "3.root";
    2747     TString extout   = "4.root";
    2748 
    2749     //------------------------------
    2750     // selection of g/h separation method
    2751     // and definition of final selections
    2752 
    2753     //TString XX("NN");
    2754     //TString XX("SC");
    2755     TString XX("RF");
    2756     TString fhadronnessName("Had");
    2757     fhadronnessName += XX;
    2758     gLog << "fhadronnessName = " << fhadronnessName << endl;
    2759 
    2760     // maximum values of the hadronness, |ALPHA| and DIST
    2761     Float_t maxhadronness   = 0.40;
    2762     Float_t maxalpha        = 20.0;
    2763     Float_t maxdist         = 10.0;
    2764     gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
    2765          << maxhadronness << ",  " << maxalpha << ",  "
    2766          << maxdist << endl;
    2767 
    2768 
    2769     //------------------------------
    2770     // name of MC file to be used for calculating the eff. collection areas
    2771     TString filenameArea(outPath);
    2772     filenameArea += typeMC;
    2773     filenameArea += "_";
    2774     filenameArea += XX;
    2775     filenameArea += ext;
    2776     gLog << "filenameArea = " << filenameArea << endl;
    2777 
    2778     //------------------------------
    2779     // name of file containing the eff. collection areas
    2780     TString collareaName(outPath);
    2781     collareaName += "area_";
    2782     collareaName += XX;
    2783     collareaName += ".root";
    2784     gLog << "collareaName = " << collareaName << endl;
    2785 
    2786     //------------------------------
    2787     // name of data file to be analysed
    2788     TString filenameData(outPath);
    2789     filenameData += typeData;
    2790     filenameData += "_";
    2791     filenameData += XX;
    2792     filenameData += ext;
    2793     gLog << "filenameData = " << filenameData << endl;
    2794 
    2795     //------------------------------
    2796     // name of output data file (after the final cuts)
    2797     TString filenameDataout(outPath);
    2798     filenameDataout += typeData;
    2799     filenameDataout += "_";
    2800     filenameDataout += XX;
    2801     filenameDataout += extout;
    2802     gLog << "filenameDataout = " << filenameDataout << endl;
    2803 
    2804 
    2805     //====================================================================
    2806     gLog << "-----------------------------------------------" << endl;
    2807     gLog << "Start calculation of effective collection areas" << endl;
    2808     MParList  parlist;
    2809     MTaskList tasklist;
    2810 
    2811     //---------------------------------------
    2812     // Setup the tasks to be executed
    2813     //
    2814     MReadMarsFile reader("Events", filenameArea);
    2815     reader.DisableAutoScheme();
    2816 
    2817     MFCT1SelFinal cuthadrons;
    2818     cuthadrons.SetHadronnessName(fhadronnessName);
    2819     cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
    2820 
    2821     MContinue conthadrons(&cuthadrons);
    2822 
    2823     //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
    2824     //MHMcCT1CollectionArea* collarea;
    2825 
    2826     MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
    2827     filler.SetName("CollectionArea");
    2828 
    2829     //********************************
    2830     // entries in MParList
    2831 
    2832     parlist.AddToList(&tasklist);
    2833     InitBinnings(&parlist);
    2834     //parlist.AddToList(collarea);
    2835 
    2836     //********************************
    2837     // entries in MTaskList
    2838 
    2839     tasklist.AddToList(&reader);   
    2840     tasklist.AddToList(&conthadrons);
    2841     tasklist.AddToList(&filler);
    2842 
    2843     //********************************
    2844 
    2845     //-----------------------------------------
    2846     // Execute event loop
    2847     //
    2848     MEvtLoop magic;
    2849     magic.SetParList(&parlist);
    2850 
    2851     MProgressBar bar;
    2852     magic.SetProgressBar(&bar);
    2853     if (!magic.Eventloop())
    2854         return;
    2855 
    2856     tasklist.PrintStatistics(0, kTRUE);
    2857 
    2858     // Calculate effective collection areas
    2859     // and display the histograms
    2860     //
    2861 
    2862     MHMcCT1CollectionArea *collarea =
    2863         (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
    2864 
    2865     collarea->CalcEfficiency();
    2866     collarea->DrawClone("lego");
    2867 
    2868     //---------------------------------------------
    2869     // Write histograms to a file
    2870     //
    2871 
    2872     TFile f(collareaName, "RECREATE");
    2873     collarea->GetHist()->Write();
    2874     collarea->GetHAll()->Write();
    2875     collarea->GetHSel()->Write();
    2876     f.Close();
    2877 
    2878     //delete collarea;
    2879 
    2880     gLog << "Calculation of effective collection areas done" << endl;
    2881     gLog << "-----------------------------------------------" << endl;   
    2882     //------------------------------------------------------------------
    2883 
    2884 
    2885     //*************************************************************************
    2886     //
    2887     // Analyse the data
    2888     //
    2889 
    2890     MTaskList tliston;
    2891     MParList pliston;
    2892 
    2893     // geometry is needed in  MHHillas... classes
    2894     MGeomCam *fGeom =
    2895              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    2896 
    2897     TString fHilName    = "MHillas";
    2898     TString fHilNameExt = "MHillasExt";
    2899     TString fHilNameSrc = "MHillasSrc";
    2900     TString fImgParName = "MNewImagePar";
    2901 
    2902     //-------------------------------------------
    2903     // create the tasks which should be executed
    2904     //
    2905 
    2906     MReadMarsFile read("Events", filenameData);
    2907     read.DisableAutoScheme();
    2908 
    2909     //.......................................................................
    2910 
    2911 
    2912       MWriteRootFile write(filenameDataout);
    2913 
    2914       write.AddContainer("MRawRunHeader", "RunHeaders");
    2915       write.AddContainer("MTime",         "Events");
    2916       write.AddContainer("MMcEvt",        "Events");
    2917       write.AddContainer("ThetaOrig",     "Events");
    2918       write.AddContainer("MSrcPosCam",    "Events");
    2919       write.AddContainer("MSigmabar",     "Events");
    2920       write.AddContainer("MHillas",       "Events");
    2921       write.AddContainer("MHillasExt",    "Events");
    2922       write.AddContainer("MHillasSrc",    "Events");
    2923       write.AddContainer("MNewImagePar",  "Events");
    2924 
    2925       //write.AddContainer("HadNN",         "Events");
    2926       write.AddContainer("HadSC",         "Events");
    2927       write.AddContainer("HadRF",         "Events");
    2928 
    2929       write.AddContainer("MEnergyEst",    "Events");
    2930 
    29312413
    29322414    //-----------------------------------------------------------------

  • trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h

    r2426 r2437  
    4242
    4343#pragma link C++ class MCT1PadSchweizer+;
    44 //#pragma link C++ class MCT1PadONOFF+;
     44#pragma link C++ class MCT1PadONOFF+;
    4545
    4646#pragma link C++ class MCT1PointingCorrCalc+;
    … …  
    5959#pragma link C++ class MMcPedestalRead+;
    6060
     61
    6162#endif
    6263
    … …  
    7071
    7172
     73
     74

  • trunk/MagicSoft/Mars/manalysis/MCT1PadONOFF.cc

    r2255 r2437  
    200200  {
    201201    edgesx[i] = ax->GetBinLowEdge(i+1);
    202     //*fLog << "i, theta low edge = " << i << ",  " << edgesx[i] << endl;
     202    // *fLog << "i, theta low edge = " << i << ",  " << edgesx[i] << endl;
    203203  }
    204204  MBinning binth;
    … …  
    213213  {
    214214    edgesy[i] = ay->GetBinLowEdge(i+1);
    215     //*fLog << "i, sigmabar low edge = " << i << ",  " << edgesy[i] << endl;
     215    // *fLog << "i, sigmabar low edge = " << i << ",  " << edgesy[i] << endl;
    216216  }
    217217  MBinning binsg;
    … …  
    989989Int_t MCT1PadONOFF::Process()
    990990{
    991   //*fLog << "Entry MCT1PadONOFF::Process();" << endl;
     991  // *fLog << "Entry MCT1PadONOFF::Process();" << endl;
    992992
    993993  //------------------------------------------------
    … …  
    10141014      continue;
    10151015
    1016     fEvt->AddPixel(i, 0.0, (*fPed)[i].GetMeanRms());
     1016    fEvt->AddPixel(i, 0.0, (*fPed)[i].GetPedestalRms());
    10171017  }
    10181018
    … …  
    10261026
    10271027  //fSigmabar->Calc(*fCam, *fPed, *fEvt);
    1028   //*fLog << "before padding : " << endl;
     1028  // *fLog << "before padding : " << endl;
    10291029  //fSigmabar->Print("");
    10301030
    … …  
    10551055    if (sigbarold > 0)
    10561056    {
    1057       //*fLog << "MCT1PadONOFF::Process(); sigmabar of event to be padded is > 0 : "
     1057      // *fLog << "MCT1PadONOFF::Process(); sigmabar of event to be padded is > 0 : "
    10581058      //      << sigbarold << ". Stop event loop " << endl;
    10591059      // input data should have sigmabar = 0; stop event loop
    … …  
    10661066
    10671067  const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
    1068   //*fLog << "theta = " << theta << endl;
     1068  // *fLog << "theta = " << theta << endl;
    10691069
    10701070  Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta);
    10711071  if ( binTheta < 1  ||  binTheta > fHBlindPixNTheta->GetNbinsX() )
    10721072  {
    1073     //*fLog << "MCT1PadONOFF::Process(); binNumber out of range : theta, binTheta = "
     1073    // *fLog << "MCT1PadONOFF::Process(); binNumber out of range : theta, binTheta = "
    10741074    //      << theta << ",  " << binTheta << ";  Skip event " << endl;
    10751075    // event cannot be padded; skip event
    … …  
    10951095    nSigma = hn->GetBinContent(binSigma);
    10961096
    1097     //*fLog << "Theta, sigbarold, binTheta, binSigma, nTheta, nSigma = "
     1097    // *fLog << "Theta, sigbarold, binTheta, binSigma, nTheta, nSigma = "
    10981098    //      << theta << ",  " << sigbarold << ",  " << binTheta << ",  "
    10991099    //      << binSigma << ",  " << nTheta << ",  " << nSigma   << endl;     
    … …  
    11111111    nSigma = hn->GetBinContent(binSigma);
    11121112
    1113     //*fLog << "Theta, sigbarold, binTheta, binSigma, nTheta, nSigma = "
     1113    // *fLog << "Theta, sigbarold, binTheta, binSigma, nTheta, nSigma = "
    11141114    //      << theta << ",  " << sigbarold << ",  " << binTheta << ",  "
    11151115    //      << binSigma << ",  " << nTheta << ",  " << nSigma   << endl;     
    … …  
    11411141  if ( nblind->Integral() == 0.0 )
    11421142  {
    1143     //*fLog << "MCT1PadONOFF::Process(); projection for Theta bin "
     1143    // *fLog << "MCT1PadONOFF::Process(); projection for Theta bin "
    11441144    //      << binTheta << " has no entries; Skip event " << endl;
    11451145    // event cannot be padded; skip event
    … …  
    11561156  delete nblind;
    11571157
    1158 #warn Code commented out due no compilation errors on Alpha OSF (tgb)
    1159 /*
     1158  //warn Code commented out due no compilation errors on Alpha OSF (tgb)
     1159
    11601160  // throw the Id of numBlind different pixels in this event
    11611161  TH1D *hblind;
    … …  
    11851185
    11861186      fBlinds->SetPixelBlind(idBlind);
    1187       //*fLog << "idBlind = " << idBlind << endl;
     1187      // *fLog << "idBlind = " << idBlind << endl;
    11881188    }
    11891189  fBlinds->SetReadyToSave();
    … …  
    11921192
    11931193  }
    1194 */
     1194
    11951195  //******************************************************************
    11961196  // has the event to be padded ?
    … …  
    11991199
    12001200  Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold);
    1201   //*fLog << "binSig, sigbarold = " << binSig << ",  " << sigbarold << endl;
     1201  // *fLog << "binSig, sigbarold = " << binSig << ",  " << sigbarold << endl;
    12021202
    12031203  Double_t prob;
    … …  
    12191219      prob = 0.0;
    12201220
    1221     //*fLog << "nTheta, nSigma, prob = " << nTheta << ",  " << nSigma
     1221    // *fLog << "nTheta, nSigma, prob = " << nTheta << ",  " << nSigma
    12221222    //      << ",  " << prob << endl;
    12231223  }
    … …  
    12381238    delete hpad;
    12391239    // event should not be padded
    1240     //*fLog << "event is not padded" << endl;
     1240    // *fLog << "event is not padded" << endl;
    12411241
    12421242    rc = 8;
    … …  
    12451245  }
    12461246  // event should be padded
    1247   //*fLog << "event is padded" << endl; 
     1247  // *fLog << "event is padded" << endl; 
    12481248
    12491249
    … …  
    12661266    sigmabar = hpad->GetRandom();
    12671267
    1268     //*fLog << "sigmabar = " << sigmabar << endl;
     1268    // *fLog << "sigmabar = " << sigmabar << endl;
    12691269
    12701270    delete hpad;
    … …  
    12991299    {
    13001300      sigmabar = hsigma->GetRandom();
    1301        //*fLog << "Theta, bin number = " << theta << ",  " << binTheta 
     1301       // *fLog << "Theta, bin number = " << theta << ",  " << binTheta 
    13021302       //      << ",  sigmabar = " << sigmabar << endl
    13031303    }
    … …  
    13091309  //-------------------------------------------
    13101310
    1311   //*fLog << "MCT1PadONOFF::Process(); sigbarold, sigmabar = "
     1311  // *fLog << "MCT1PadONOFF::Process(); sigbarold, sigmabar = "
    13121312  //      << sigbarold << ",  "<< sigmabar << endl;
    13131313
    … …  
    13671367
    13681368  elNoise2 = TMath::Min(RMS,  sigmabar2 - sigbarold2);
    1369   //*fLog << "elNoise2 = " << elNoise2 << endl;
     1369  // *fLog << "elNoise2 = " << elNoise2 << endl;
    13701370
    13711371  lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess;     // [photons]
    … …  
    14051405
    14061406    MPedestalPix &ppix = (*fPed)[j];
    1407     Double_t oldsigma = ppix.GetMeanRms();
     1407    Double_t oldsigma = ppix.GetPedestalRms();
    14081408    Double_t oldsigma2 = oldsigma*oldsigma;
    14091409
    … …  
    14631463      if (!ok)
    14641464      {
    1465         //*fLog << "theta, j, count, sigmabar, diff = " << theta << ",  "
     1465        // *fLog << "theta, j, count, sigmabar, diff = " << theta << ",  "
    14661466        //      << j << ",  " << count << ",  " << sigmabar << ",  "
    14671467        //      << diff << endl;
    … …  
    15241524      if (!ok)
    15251525      {
    1526         //*fLog << "theta, j, count, sigmabar, sig = " << theta << ",  "
     1526        // *fLog << "theta, j, count, sigmabar, sig = " << theta << ",  "
    15271527        //      << j << ",  " << count << ",  " << sigmabar << ",  "
    15281528        //      << sig << endl;
    … …  
    16021602
    16031603    Double_t newsigma = sqrt( oldsigma2 + addSig2 );
    1604     ppix.SetMeanRms( newsigma );
     1604    ppix.SetPedestalRms( newsigma );
    16051605
    16061606    fHSigmaPedestal->Fill( oldsigma, newsigma );
    … …  
    16111611
    16121612  //fSigmabar->Calc(*fCam, *fPed, *fEvt);
    1613   //*fLog << "after padding : " << endl;
     1613  // *fLog << "after padding : " << endl;
    16141614  //fSigmabar->Print("");
    16151615

  • trunk/MagicSoft/Mars/manalysis/Makefile

    r2426 r2437  
    6666           MMinuitInterface.cc \
    6767           MFiltercutsCalc.cc \
     68           MCT1PadONOFF.cc \
    6869           MMcPedestalRead.cc
    69 #          MCT1PadONOFF.cc \
     70
    7071
    7172SRCS    = $(SRCFILES)
Note: See TracChangeset for help on using the changeset viewer.