Ignore:
Timestamp:
09/24/03 13:29:00 (21 years ago)
Author:
wittek
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/macros
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r2317 r2357  
    2222#include "MWriteRootFile.h"
    2323
     24//#include "TH3D.h"
    2425
    2526
    … …  
    7778
    7879        MBinning *binth = new MBinning("BinningTheta");
    79         //Double_t yedge[9] =
    80         //               {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
    8180        // this is Daniel's binning in theta
    82         Double_t yedge[8] =
    83           {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99};
    84         //TArrayD yed(9,yedge);
     81        //Double_t yedge[8] =
     82        //  {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99};
     83        // this is our binning
     84        Double_t yedge[9] =
     85                       {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
    8586        TArrayD yed;
    86         yed.Set(8,yedge);
     87        yed.Set(9,yedge);
    8788        binth->SetEdges(yed);
    8889        plist->AddToList(binth);
    … …  
    9495          zedge[8-i] = cos(yedge[i]/kRad2Deg);
    9596        }
    96         TArrayD zed(9,zedge);
     97        TArrayD zed;
     98        zed.Set(9,zedge);
    9799        bincosth->SetEdges(zed);
    98100        plist->AddToList(bincosth);
    … …  
    116118}
    117119
     120
    118121void DeleteBinnings(MParList *plist)
    119122{
    … …  
    122125        TObject *bin;
    123126
     127        //--------------------------------------------
     128        bin = plist->FindObject("BinningE");
     129        if (bin) delete bin;
     130
     131        //--------------------------------------------
     132
    124133        bin = plist->FindObject("BinningSize");
    125134        if (bin) delete bin;
    … …  
    159168        if (bin) delete bin;
    160169
    161         //robert
     170
     171        // robert ----------------------------------------------
    162172        bin = plist->FindObject("BinningAlphaFlux");
    163173        if (bin) delete bin;
    … …  
    170180}
    171181
     182
    172183//************************************************************************
    173184void CT1Analysis()
    174185{
     186
     187  gLog << "Entry CT1Analysis()" << endl;
     188
    175189      gLog.SetNoColors();
    176190
    … …  
    189203
    190204      // path for input for Mars
    191       TString inPath = "~magican/ct1test/wittek/marsoutput/optionB/";
     205      TString inPath = "~magican/ct1test/wittek/marsoutput/optionC/";
    192206
    193207      // path for output from Mars
    194       TString outPath = "~magican/ct1test/wittek/marsoutput/optionB/";
     208      TString outPath = "~magican/ct1test/wittek/marsoutput/optionC/";
    195209
    196210      //-----------------------------------------------
    … …  
    217231    Bool_t JobA_MC  = kFALSE; 
    218232    Bool_t WMC1     = kFALSE;  // write out root file MC1.root ?
    219 
    220 
    221     // Job B_NN_UP : read ON1.root (or MC1.root) file
    222     //  - depending on RlookNN : create (or read in) hadron and gamma matrix
    223     //  - calculate hadroness for method of NEAREST NEIGHBORS
    224     //    and for the SUPERCUTS
    225     //  - update the input files with the hadronesses (ON1.root or MC1.root)
    226 
    227     Bool_t JobB_NN_UP  = kFALSE; 
    228     Bool_t RLookNN     = kFALSE;  // read in look-alike events
    229     Bool_t WNN         = kFALSE;  // update input root file ?
    230 
    231 
    232     // Job B_SC_UP : read ON1.root (or MC1.root) file
    233     //  - depending on WParSC : create (or read in) supercuts parameter values
    234     //  - calculate hadroness for the SUPERCUTS
    235     //  - update the input files with the hadroness (ON1.root or MC1.root)
    236 
    237     Bool_t JobB_SC_UP  = kTRUE;
    238     Bool_t RMatrix     = kFALSE;  // read matrices from file 
    239     Bool_t WParSC      = kTRUE;  // do optimization and write supercuts
    240                                   // parameter values onto a file
    241     Bool_t WSC         = kFALSE;  // update input root file ?
    242233
    243234
    … …  
    252243    Bool_t RTree       = kFALSE;  // read in trees
    253244    Bool_t WRF         = kFALSE;  // update input root file ?
     245
     246
     247
     248
     249    // Job B_SC_UP : read ON1.root (or MC1.root) file
     250    //  - depending on WParSC : create (or read in) supercuts parameter values
     251    //  - calculate hadroness for the SUPERCUTS
     252    //  - update the input files with the hadroness (ON1.root or MC1.root)
     253
     254    Bool_t JobB_SC_UP  = kFALSE;
     255    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file 
     256    Bool_t WParSC      = kFALSE;  // do optimization and write supercuts
     257                                 // parameter values onto a file
     258    Bool_t WSC         = kFALSE; // update input root file ?
    254259
    255260
    … …  
    299304    //  - write root file for ON data after final cuts
    300305
    301     Bool_t JobF_XX  = kFALSE; 
    302     Bool_t WFX      = kFALSE;  // write out root file  ?
     306    gLog << "before setting switches for JobF_XX" << endl;
     307
     308    Bool_t JobF_XX  = kTRUE; 
     309    Bool_t WFX      = kTRUE;  // write out root file  ?
     310    Bool_t WRobert  = kTRUE;  // write out Robert's file  ?
    303311
    304312
    … …  
    327335    gLog << "" << endl;
    328336    gLog << "Macro CT1Analysis : JobA_ON, WHistON, WON1 = "
    329          << JobA_ON  << ",  " << WHistON << ",  " << WON1 << endl;
     337         << (JobA_ON ? "kTRUE" : "kFALSE")  << ",  "
     338         << (WHistON ? "kTRUE" : "kFALSE")  << ",  "
     339         << (WON1    ? "kTRUE" : "kFALSE")  << endl;
    330340
    331341
    … …  
    424434    contstandard.SetName("SelStandard");
    425435
    426     if (WON1)
    427     {
    428       MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
     436      //MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
     437      MWriteRootFile write(outNameImage, "RECREATE");
    429438
    430439      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    438447      write.AddContainer("MHillasSrc",    "Events");
    439448      write.AddContainer("MNewImagePar",  "Events");
    440     }
     449
    441450
    442451    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    … …  
    500509
    501510    Int_t maxevents = -1;
    502     //Int_t maxevents = 1000;
    503511    if ( !evtloop.Eventloop(maxevents) )
    504512        return;
    … …  
    587595    gLog << "" << endl;
    588596    gLog << "Macro CT1Analysis : JobA_MC, WMC1 = "
    589          << JobA_MC  << ",  " << WMC1 << endl;
     597         << (JobA_MC ? "kTRUE" : "kFALSE")  << ",  "
     598         << (WMC1    ? "kTRUE" : "kFALSE")  << endl;
    590599
    591600
    … …  
    688697    MCT1ReadPreProc read(filenamein);
    689698
     699    MBlindPixelCalc blindbeforepad;
     700    blindbeforepad.SetUseBlindPixels();
     701    blindbeforepad.SetName("BlindBeforePadding");
     702
    690703    MBlindPixelCalc blind;
    691704    blind.SetUseBlindPixels();
     705    blind.SetName("BlindAfterPadding");
    692706
    693707    MFCT1SelBasic selbasic;
    … …  
    758772
    759773
    760     if (WMC1)
    761     {
    762       MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
     774      //MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
     775      MWriteRootFile write(outNameImage, "RECREATE");
    763776
    764777      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    772785      write.AddContainer("MHillasSrc",    "Events");
    773786      write.AddContainer("MNewImagePar",  "Events");
    774     }
     787
    775788
    776789
    … …  
    788801
    789802    tlist.AddToList(&read);
     803    tlist.AddToList(&blindbeforepad);
    790804    tlist.AddToList(&padthomas);
    791805    tlist.AddToList(&blind);
    … …  
    855869 }
    856870
    857 
    858 
    859871  //---------------------------------------------------------------------
    860   // Job B_NN_UP
     872  // Job B_RF_UP
    861873  //============
    862874
    863     //  - read in the matrices of look-alike events for gammas and hadrons
    864 
    865     // then read ON1.root (or MC1.root) file
    866     //
    867     //  - calculate the hadroness for the method of nearest neighbors
    868     //
    869     //  - update input root file, including the hadroness
    870 
    871 
    872  if (JobB_NN_UP)
     875
     876    //  - create (or read in) the matrices of look-alike events for gammas
     877    //    and hadrons
     878    //  - create (or read in) the trees
     879    //  - then read ON1.root (or MC1.root) file
     880    //  - calculate the hadroness for the method of RANDOM FOREST
     881    //  - update input root file with the hadroness
     882
     883
     884 if (JobB_RF_UP)
    873885 {
    874886    gLog << "=====================================================" << endl;
    875     gLog << "Macro CT1Analysis : Start of Job B_NN_UP" << endl;
     887    gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
    876888
    877889    gLog << "" << endl;
    878     gLog << "Macro CT1Analysis : JobB_NN_UP, RLookNN, WNN = "
    879          << JobB_NN_UP  << ",  " << RLookNN << ",  " << WNN << endl;
    880 
     890    gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = "
     891         << (JobB_RF_UP ? "kTRUE" : "kFALSE")  << ",  "
     892         << (RLookRF    ? "kTRUE" : "kFALSE")  << ",  "
     893         << (RTree      ? "kTRUE" : "kFALSE")  << ",  "
     894         << (WRF        ? "kTRUE" : "kFALSE")  << endl;
     895
     896
     897    //--------------------------------------------
     898    // parameters for the random forest
     899    Int_t NumTrees = 100;
     900    Int_t NumTry   =   3;
     901    Int_t NdSize   =   1;
    881902
    882903
    … …  
    898919    outNameImage += typeInput;
    899920    outNameImage += "2.root";
    900 
    901921    //TString outNameImage = filenameData;
    902922
    … …  
    909929    filenameON += "ON";
    910930    filenameON += "1.root";
    911     Int_t howManyHadrons = 500000;
     931    Int_t howManyHadrons = 35000;
    912932    gLog << "filenameON = " << filenameON << ",   howManyHadrons = "
    913933         << howManyHadrons  << endl;
    … …  
    918938    filenameMC += "MC";
    919939    filenameMC += "1.root";
    920     Int_t howManyGammas = 50000;
     940    Int_t howManyGammas = 35000;
    921941    gLog << "filenameMC = " << filenameMC << ",   howManyGammas = "
    922942         << howManyGammas   << endl;
    … …  
    926946
    927947    TString outNameGammas = outPath;
    928     outNameGammas += "matrix_gammas_";
     948    outNameGammas += "RFmatrix_gammas_";
    929949    outNameGammas += "MC";
    930950    outNameGammas += ".root";
    … …  
    934954
    935955    TString outNameHadrons = outPath;
    936     outNameHadrons += "matrix_hadrons_";
     956    outNameHadrons += "RFmatrix_hadrons_";
    937957    outNameHadrons += typeMatrixHadrons;
    938958    outNameHadrons += ".root";
    … …  
    940960
    941961    MHMatrix matrixg("MatrixGammas");
     962    matrixg.EnableGraphicalOutput();
     963
    942964    MHMatrix matrixh("MatrixHadrons");
     965    matrixh.EnableGraphicalOutput();
     966
     967    //--------------------------------------------
     968    // file of trees of the random forest
     969
     970    TString outRF = outPath;
     971    outRF += "RF.root";
     972
    943973
    944974   //*************************************************************************
    945975   // read in matrices of look-alike events
    946 if (RLookNN)
     976if (RLookRF)
    947977  {
    948978    const char* mtxName = "MatrixGammas";
    … …  
    961991
    962992    matrixg.Read(mtxName);
    963     matrixg.Print("cols");
     993    matrixg.Print("SizeCols");
    964994
    965995
    … …  
    9811011
    9821012    matrixh.Read(mtxName);
    983     matrixh.Print("cols");
     1013    matrixh.Print("SizeCols");
    9841014  }
    9851015
    … …  
    9871017   //*************************************************************************
    9881018   // create matrices of look-alike events
    989 if (!RLookNN)
     1019if (!RLookRF)
    9901020  {
    9911021    gLog << "" << endl;
    … …  
    10111041    MFParticleId fgamma("MMcEvt", '=', kGAMMA);
    10121042    fgamma.SetName("gammaID");
     1043
    10131044    MFParticleId fhadrons("MMcEvt", '!', kGAMMA);
    10141045    fhadrons.SetName("hadronID)");
    10151046
    1016     MFEventSelector selectorg;
    1017     selectorg.SetNumSelectEvts(howManyGammas);
     1047    TString mgname("costhg");
     1048    MBinning bing("Binning"+mgname);
     1049    bing.SetEdges(10, 0., 1.0);
     1050
     1051    MH3 gref("cos(MMcEvt.fTelescopeTheta)");
     1052    gref.SetName(mgname);
     1053    MH::SetBinning(&gref.GetHist(), &bing);
     1054    for (Int_t i=1; i<=gref.GetNbins(); i++)
     1055      gref.GetHist().SetBinContent(i, 1.0);
     1056
     1057
     1058    MFEventSelector2 selectorg(gref);
     1059    selectorg.SetNumMax(howManyGammas);
    10181060    selectorg.SetName("selectGammas");
    1019     MFEventSelector selectorh;
    1020     selectorh.SetNumSelectEvts(howManyHadrons);
    1021     selectorh.SetName("selectHadrons");
     1061
    10221062
    10231063    MFillH fillmatg("MatrixGammas");
    … …  
    10291069    fillmath.SetName("fillHadrons");
    10301070
    1031    
     1071
    10321072    //*****************************   fill gammas   *** 
    10331073    // entries in MFilterList
    … …  
    10651105    tlistg.PrintStatistics(0, kTRUE);
    10661106
    1067 
     1107   
    10681108    //*****************************   fill hadrons   *** 
    10691109
    10701110    gLog << " Hadrons :" << endl;
     1111
     1112    TString mhname("costhh");
     1113    MBinning binh("Binning"+mhname);
     1114    binh.SetEdges(10, 0., 1.0);
     1115
     1116    MH3 href("cos(MMcEvt.fTelescopeTheta)");
     1117    href.SetName(mhname);
     1118    MH::SetBinning(&href.GetHist(), &binh);
     1119    for (Int_t i=1; i<=href.GetNbins(); i++)
     1120      href.GetHist().SetBinContent(i, 1.0);
     1121
     1122    MFEventSelector2 selectorh(href);
     1123    //selectorh.SetNumMax(howManyHadrons);
     1124    // select as many hadrons as gammas
     1125    selectorh.SetNumMax(matrixg.GetM().GetNrows());
     1126    selectorh.SetName("selectHadrons");
     1127
     1128
    10711129    // entries in MFilterList
    10721130
    … …  
    11021160
    11031161    tlisth.PrintStatistics(0, kTRUE);
     1162
     1163
     1164
    11041165    //***************************************************** 
    11051166
    1106     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    1107     //
    1108     // ----------------------------------------------------------
    1109     //  Definition of the reference sample of look-alike events
    1110     //  (this is a subsample of the original sample)
    1111     // ----------------------------------------------------------
    1112     //
    1113     gLog << "" << endl;
    1114     gLog << "========================================================" << endl;
    1115     // Select a maximum of nmaxevts events from the sample of look-alike
    1116     // events. They will form the reference sample.
    1117     Int_t nmaxevts  = 2000;
    1118 
    1119     // target distribution for the variable in column refcolumn (start at 0);
    1120     //Int_t   refcolumn = 0;
    1121     //Int_t   nbins   =   5;
    1122     //Float_t frombin = 0.5;
    1123     //Float_t tobin   = 1.0;
    1124     //TH1F *thsh = new TH1F("thsh","target distribution",
    1125     //                       nbins, frombin, tobin);
    1126     //thsh->SetDirectory(NULL);
    1127     //thsh->SetXTitle("cos( \\Theta)");
    1128     //thsh->SetYTitle("Counts");
    1129     //Float_t dbin = (tobin-frombin)/nbins;
    1130     //Float_t lbin = frombin +dbin*0.5;
    1131     //for (Int_t j=1; j<=nbins; j++)
    1132     //{
    1133     //  thsh->Fill(lbin,1.0);
    1134     //  lbin += dbin;
    1135     //}
    1136 
    1137     Int_t   refcolumn = 0;
    1138     MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning");
    1139     TH1F *thsh = new TH1F();
    1140     thsh->SetNameTitle("thsh","target distribution");
    1141     MH::SetBinning(thsh, binscostheta);
    1142     Int_t nbins = thsh->GetNbinsX();
    1143     Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900};
    1144     for (Int_t j=1; j<=nbins; j++)
    1145     {
    1146       thsh->Fill(thsh->GetBinCenter(j), cont[j-1]);
    1147     }
    1148 
    1149     gLog << "" << endl;
    1150     gLog << "========================================================" << endl;
    1151     gLog << "Macro CT1Analysis : define reference sample for gammas" << endl;
    1152     gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
    1153          << refcolumn << ",  " << nmaxevts << endl;   
    1154 
    1155     matrixg.EnableGraphicalOutput();
    1156     Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    1157 
    1158     if ( !matrixok )
    1159     {
    1160       gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
    1161       return;
    1162     }
    1163 
    1164     gLog << "" << endl;
    1165     gLog << "========================================================" << endl;
    1166     gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl;
    1167     gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
    1168          << refcolumn << ",  " << nmaxevts << endl;   
    1169 
    1170     matrixh.EnableGraphicalOutput();
    1171     matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    1172     delete thsh;
    1173 
    1174     if ( !matrixok )
    1175     {
    1176       gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
    1177       return;
    1178     }
    1179     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    11801167
    11811168    // write out look-alike events
    … …  
    11891176      // "gammas"
    11901177      gLog << "Gammas :" << endl;   
    1191       matrixg.Print("cols");
     1178      matrixg.Print("SizeCols");
    11921179
    11931180      TFile writeg(outNameGammas, "RECREATE", "");
    … …  
    12011188      // "hadrons"
    12021189      gLog << "Hadrons :" << endl;   
    1203       matrixh.Print("cols");
     1190      matrixh.Print("SizeCols");
    12041191
    12051192      TFile writeh(outNameHadrons, "RECREATE", "");
    … …  
    12141201
    12151202
     1203    MRanForest *fRanForest;
     1204    MRanTree *fRanTree;
    12161205    //-----------------------------------------------------------------
    1217     // Update the input files with the NN and SC hadronness
    1218     //
    1219 
    1220  if (WNN)
     1206    // read in the trees of the random forest
     1207    if (RTree)
     1208    {
     1209      MParList plisttr;
     1210      MTaskList tlisttr;
     1211      plisttr.AddToList(&tlisttr);
     1212
     1213      MReadTree readtr("TREE", outRF);
     1214      readtr.DisableAutoScheme();
     1215
     1216      MRanForestFill rffill;
     1217      rffill.SetNumTrees(NumTrees);
     1218
     1219      // list of tasks for the loop over the trees
     1220
     1221      tlisttr.AddToList(&readtr);
     1222      tlisttr.AddToList(&rffill);
     1223
     1224      //-------------------
     1225      // Execute tree loop
     1226      //
     1227      MEvtLoop evtlooptr;
     1228      evtlooptr.SetParList(&plisttr);
     1229      if (!evtlooptr.Eventloop())
     1230        return;
     1231
     1232      tlisttr.PrintStatistics(0, kTRUE);
     1233
     1234
     1235    // get adresses of objects which are used in the next eventloop
     1236    fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
     1237    if (!fRanForest)
     1238    {
     1239        *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
     1240        return kFALSE;
     1241    }
     1242
     1243    fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
     1244    if (!fRanTree)                                 
     1245    {                                                                         
     1246        *fLog << err << dbginf << "MRanTree not found... aborting." << endl;   
     1247        return kFALSE;
     1248    }
     1249
     1250    }
     1251
     1252    //-----------------------------------------------------------------
     1253    // grow the trees of the random forest (event loop = tree loop)
     1254
     1255    if (!RTree)
     1256    {
     1257
     1258    gLog << "" << endl;
     1259    gLog << "========================================================" << endl;
     1260    gLog << "Macro CT1Analysis : start growing trees" << endl;
     1261
     1262    MTaskList tlist2;
     1263    MParList plist2;
     1264    plist2.AddToList(&tlist2);
     1265
     1266    plist2.AddToList(&matrixg);
     1267    plist2.AddToList(&matrixh);
     1268
     1269    MRanForestGrow rfgrow2;
     1270    rfgrow2.SetNumTrees(NumTrees);
     1271    rfgrow2.SetNumTry(NumTry);
     1272    rfgrow2.SetNdSize(NdSize);
     1273
     1274    MWriteRootFile rfwrite2(outRF);
     1275    rfwrite2.AddContainer("MRanTree", "TREE");
     1276
     1277    // list of tasks for the loop over the trees
     1278   
     1279    tlist2.AddToList(&rfgrow2);
     1280    tlist2.AddToList(&rfwrite2);
     1281
     1282    //-------------------
     1283    // Execute tree loop
     1284    //
     1285    MEvtLoop treeloop;
     1286    treeloop.SetParList(&plist2);
     1287
     1288    if ( !treeloop.Eventloop() )
     1289        return;
     1290
     1291    tlist2.PrintStatistics(0, kTRUE);
     1292
     1293    // get adresses of objects which are used in the next eventloop
     1294    fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
     1295    if (!fRanForest)
     1296    {
     1297        *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
     1298        return kFALSE;
     1299    }
     1300
     1301    fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
     1302    if (!fRanTree)                                 
     1303    {                                                                         
     1304        *fLog << err << dbginf << "MRanTree not found... aborting." << endl;   
     1305        return kFALSE;
     1306    }
     1307
     1308    }
     1309    // end of growing the trees of the random forest
     1310    //-----------------------------------------------------------------
     1311
     1312
     1313
     1314
     1315    //-----------------------------------------------------------------
     1316    // Update the input files with the RF hadronness
     1317    //
     1318
     1319 if (WRF)
    12211320  {
    12221321    gLog << "" << endl;
    12231322    gLog << "========================================================" << endl;
    12241323    gLog << "Update input file '" <<  filenameData
    1225          << "' with the NN and SC hadronness" << endl;
     1324         << "' with the RF hadronness" << endl;
    12261325
    12271326    MTaskList tliston;
    … …  
    12461345
    12471346    //.......................................................................
    1248     // calculation of hadroness for method of Nearest Neighbors
    1249 
    1250     TString hadNNName = "HadNN";
    1251     MMultiDimDistCalc nncalc;
    1252     nncalc.SetUseNumRows(25);
    1253     nncalc.SetUseKernelMethod(kFALSE);
    1254     nncalc.SetHadronnessName(hadNNName);
    1255 
    1256     //.......................................................................
    1257     // calculation of hadroness for the supercuts
    1258     // (=0.25 if fullfilled, =0.75 otherwise)
    1259 
    1260     TString hadSCName = "HadSC";
    1261     MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc);
    1262     sccalc.SetHadronnessName(hadSCName);
     1347    // calculate hadronnes for method of RANDOM FOREST
     1348
     1349    TString hadRFName = "HadRF";
     1350    MRanForestCalc rfcalc;
     1351    rfcalc.SetHadronnessName(hadRFName);
    12631352
    12641353    //.......................................................................
    … …  
    12781367      write.AddContainer("MNewImagePar",  "Events");
    12791368
    1280       write.AddContainer("HadNN",         "Events");
    1281       write.AddContainer("HadSC",         "Events");
     1369      write.AddContainer(hadRFName,       "Events");
    12821370
    12831371
    … …  
    12871375             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    12881376
     1377
    12891378    Float_t maxhadronness =  0.40;
    12901379    Float_t maxalpha      =  20.0;
    … …  
    12931382    MFCT1SelFinal selfinalgh(fHilNameSrc);
    12941383    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    1295     selfinalgh.SetHadronnessName(hadNNName);
     1384    selfinalgh.SetHadronnessName(hadRFName);
    12961385    selfinalgh.SetName("SelFinalgh");
    12971386    MContinue contfinalgh(&selfinalgh);
    12981387    contfinalgh.SetName("ContSelFinalgh");
    12991388
    1300     MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
    1301     fillhadnn.SetName("HhadNN");
    1302     MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
    1303     fillhadsc.SetName("HhadSC");
     1389    MFillH fillranfor("MHRanForest");
     1390    fillranfor.SetName("HRanForest");
     1391
     1392    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
     1393    fillhadrf.SetName("HhadRF");
    13041394
    13051395    MFCT1SelFinal selfinal(fHilNameSrc);
    13061396    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    1307     selfinal.SetHadronnessName(hadNNName);
     1397    selfinal.SetHadronnessName(hadRFName);
    13081398    selfinal.SetName("SelFinal");
    13091399    MContinue contfinal(&selfinal);
    13101400    contfinal.SetName("ContSelFinal");
    13111401
     1402    TString mh3name = "abs(Alpha)";
     1403    MBinning binsalphaabs("Binning"+mh3name);
     1404    binsalphaabs.SetEdges(50, -2.0, 98.0);
     1405
     1406    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
     1407    alphaabs.SetName(mh3name);
     1408    MFillH alpha(&alphaabs);
     1409    alpha.SetName("FillAlphaAbs");
    13121410
    13131411    MFillH hfill1("MHHillas",    fHilName);
    … …  
    13321430    InitBinnings(&pliston);
    13331431
    1334     pliston.AddToList(&matrixg);
    1335     pliston.AddToList(&matrixh);
    1336 
     1432    pliston.AddToList(fRanForest);
     1433    pliston.AddToList(fRanTree);
     1434
     1435    pliston.AddToList(&binsalphaabs);
     1436    pliston.AddToList(&alphaabs);
    13371437
    13381438    //*****************************
    … …  
    13411441    tliston.AddToList(&read);
    13421442
    1343     tliston.AddToList(&nncalc);
    1344     tliston.AddToList(&sccalc);
    1345     tliston.AddToList(&fillhadnn);
    1346     tliston.AddToList(&fillhadsc);
    1347 
     1443    tliston.AddToList(&rfcalc);
     1444    tliston.AddToList(&fillranfor);
     1445    tliston.AddToList(&fillhadrf);
     1446
     1447    tliston.AddToList(&write);
     1448    tliston.AddToList(&contfinalgh);
     1449
     1450    tliston.AddToList(&alpha);
    13481451    tliston.AddToList(&hfill1);
    13491452    tliston.AddToList(&hfill2);
    … …  
    13521455    tliston.AddToList(&hfill5);
    13531456
    1354     tliston.AddToList(&write);
    1355 
    1356     tliston.AddToList(&contfinalgh);
    13571457    tliston.AddToList(&contfinal);
    13581458
    … …  
    13731473    tliston.PrintStatistics(0, kTRUE);
    13741474
     1475
    13751476    //-------------------------------------------
    13761477    // Display the histograms
    13771478    //
    1378     pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
    1379     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
     1479    pliston.FindObject("MHRanForest")->DrawClone();
     1480    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    13801481
    13811482    pliston.FindObject("MHHillas")->DrawClone();
    … …  
    13851486    pliston.FindObject("MHStarMap")->DrawClone();
    13861487
     1488     //-------------------------------------------
     1489    // fit alpha distribution to get the number of excess events and
     1490    // calculate significance of gamma signal in the alpha plot
     1491 
     1492    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
     1493    TH1  &alphaHist = absalpha->GetHist();
     1494    alphaHist.SetXTitle("|alpha|  [\\circ]");
     1495
     1496    Double_t alphasig = 13.1;
     1497    Double_t alphamin = 30.0;
     1498    Double_t alphamax = 90.0;
     1499    Int_t    degree   =    4;
     1500    Double_t significance = -99.0;
     1501    Bool_t   drawpoly  = kTRUE;
     1502    Bool_t   fitgauss  = kTRUE;
     1503    Bool_t   print     = kTRUE;
     1504
     1505    MHFindSignificance findsig;
     1506    findsig.SetRebin(kTRUE);
     1507    findsig.SetReduceDegree(kFALSE);
     1508
     1509    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
     1510                        alphasig, drawpoly, fitgauss, print);
     1511    significance = findsig.GetSignificance();
     1512    Float_t alphasi = findsig.GetAlphasi();
     1513
     1514    gLog << "For file '" << filenameData << "' : " << endl;
     1515    gLog << "Significance of gamma signal after supercuts : "
     1516         << significance << " (for |alpha| < " << alphasi << " degrees)"
     1517         << endl;
     1518
     1519    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
     1520
     1521    //-------------------------------------------
     1522
     1523
    13871524    DeleteBinnings(&pliston);
    13881525  }
    13891526
    1390 
    1391 
    1392     gLog << "Macro CT1Analysis : End of Job B_NN_UP" << endl;
     1527    gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
    13931528    gLog << "=======================================================" << endl;
    13941529 }
    13951530  //---------------------------------------------------------------------
     1531
     1532
    13961533
    13971534
    … …  
    14131550
    14141551    gLog << "" << endl;
    1415     //gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = "
    1416     //     << JobB_SC_UP  << ",  " << WParSC << ",  " << WSC << endl;
    1417 
     1552    gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = "
     1553         << (JobB_SC_UP ? "kTRUE" : "kFALSE")  << ",  "
     1554         << (WParSC     ? "kTRUE" : "kFALSE")  << ",  "
     1555         << (WSC        ? "kTRUE" : "kFALSE")  << endl;
     1556
     1557
     1558
     1559    //--------------------------------------------
     1560    // file which contains the initial parameter values for the supercuts
     1561    // if parSCinit ="" the initial values are taken from the constructor of
     1562    //                  MCT1Supercuts
     1563
     1564    TString parSCinit = outPath;
     1565    //parSCinit += "parSC_1709d";
     1566    parSCinit = "";
     1567
     1568    gLog << "parSCinit = " << parSCinit << endl;
     1569
     1570    //---------------
     1571    // file onto which the optimal parameter values for the supercuts
     1572    // are written
     1573
     1574    TString parSCfile = outPath;
     1575    parSCfile += "parSC_2209a";
     1576
     1577    gLog << "parSCfile = " << parSCfile << endl;
    14181578
    14191579    //--------------------------------------------
    14201580    // file to be updated (either ON or MC)
    14211581
    1422     TString typeInput = "ON";
    1423     //TString typeInput = "MC";
     1582    //TString typeInput = "ON";
     1583    TString typeInput = "MC";
    14241584    gLog << "typeInput = " << typeInput << endl;
    14251585
    … …  
    14341594    outNameImage += typeInput;
    14351595    outNameImage += "3.root";
    1436     outNameImage += "xxxxx";
    14371596   
    14381597
    … …  
    14471606    TString filenameTrain = outPath;
    14481607    filenameTrain += "ON";
    1449     filenameTrain += "2.root";
     1608    filenameTrain += "1.root";
    14501609    Int_t howManyTrain = 800000;
    14511610    gLog << "filenameTrain = " << filenameTrain << ",   howManyTrain = "
    … …  
    14551614    TString filenameTest = outPath;
    14561615    filenameTest += "ON";
    1457     filenameTest += "2.root";
    1458     Int_t howManyTest = 100000;
     1616    filenameTest += "1.root";
     1617    Int_t howManyTest = 800000;
    14591618
    14601619    gLog << "filenameTest = " << filenameTest << ",   howManyTest = "
    … …  
    14791638
    14801639   
    1481     //--------------------------------------------
    1482     // file which contains the optimum parameter values for the supercuts
    1483 
    1484     TString parSCfile = outPath;
    1485     parSCfile += "parSC";
    1486 
    1487     gLog << "parSCfile = " << parSCfile << endl;
    1488 
    14891640
    14901641    //---------------------------------------------------------------------
    1491     // optimize supercuts using the training sample
    1492     // and test them using the test sample
    1493 
    1494 if (WParSC)
    1495   {
     1642    // Training and test matrices :
     1643    // - either create them and write them onto a file
     1644    // - or read them from a file
     1645
     1646
    14961647    MCT1FindSupercuts findsuper;
    14971648    findsuper.SetFilenameParam(parSCfile);
    … …  
    15001651
    15011652    //--------------------------
    1502     // create matrices
     1653    // create matrices and write them onto files
    15031654    if (!RMatrix)
    15041655    {
    … …  
    15091660      {
    15101661        if ( !findsuper.DefineTrainTestMatrix(filenameTrain,
    1511                                               howManyTrain, mh3,
    1512                                               howManyTest,  mh3) )
     1662                              howManyTrain, mh3, howManyTest, mh3,
     1663                              fileMatrixTrain, fileMatrixTest)    )
    15131664        {
    15141665          *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl;
    … …  
    15191670      else
    15201671      {
    1521         if ( !findsuper.DefineTrainMatrix(filenameTrain, howManyTrain, mh3) )
     1672        if ( !findsuper.DefineTrainMatrix(filenameTrain,
     1673                              howManyTrain, mh3, fileMatrixTrain) )
    15221674        {
    15231675          *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl;
    … …  
    15251677        }
    15261678
    1527         if ( !findsuper.DefineTestMatrix( filenameTest,  howManyTest,  mh3) )
     1679        if ( !findsuper.DefineTestMatrix( filenameTest, 
     1680                              howManyTest,  mh3, fileMatrixTest)  )
    15281681        {
    15291682          *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl;
    … …  
    15311684        }
    15321685      }
    1533  
    1534       // writing doesn't work ???
    1535       //findsuper.WriteMatrix(fileMatrixTrain, fileMatrixTest);
    1536     }
     1686     }
    15371687
    15381688    //--------------------------
    15391689    // read matrices from files
    1540     //                              NOT YET TESTED
    1541     //if (RMatrix)
    1542     //  findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
     1690    //                             
     1691
     1692    if (RMatrix)
     1693      findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
    15431694    //--------------------------
    15441695
    15451696
    1546     //--------------------------
    1547     // optimize the supercuts
    1548     Bool_t rf = findsuper.FindParams();
     1697
     1698    //---------------------------------------------------------------------
     1699    // optimize supercuts using the training sample
     1700    //
     1701    // the initial values are taken
     1702    //     - from the file parSCinit (if != "")
     1703    //     - or from the arrays params and steps (if their sizes are != 0)
     1704    //     - or from the MCT1Supercuts constructor
     1705
     1706
     1707if (WParSC)
     1708  {
     1709    TArrayD params(0);
     1710    TArrayD steps(0);
     1711 
     1712    if (parSCinit == "")
     1713    {
     1714      Double_t vparams[104] = {
     1715      // LengthUp
     1716        0.315585,  0.001455, 0.203198, 0.005532, -0.001670, -0.020362,
     1717        0.007388, -0.013463,
     1718      // LengthLo
     1719        0.151530,  0.028323, 0.510707, 0.053089,  0.013708,  2.357993,
     1720        0.000080, -0.007157,
     1721      // WidthUp
     1722        0.145412, -0.001771, 0.054462, 0.022280, -0.009893,  0.056353,
     1723        0.020711, -0.016703,
     1724      // WidthLo
     1725        0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101,
     1726        0.006126, -0.002849,
     1727      // DistUp
     1728        1.787943,  0.0,      2.942310, 0.199815,  0.0,       0.249909,
     1729        0.189697,  0.0,
     1730      // DistLo
     1731        0.589406,  0.0,     -0.083964,-0.007975,  0.0,       0.045374,
     1732       -0.001750,  0.0,
     1733      // AsymUp
     1734        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1735        0.0,       0.0,
     1736      // AsymLo
     1737       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1738        0.0,       0.0,
     1739      // ConcUp
     1740        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1741        0.0,       0.0,
     1742      // ConcLo
     1743       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1744        0.0,       0.0,
     1745      // Leakage1Up
     1746        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1747        0.0,       0.0,
     1748      // Leakage1Lo
     1749       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1750        0.0,       0.0,
     1751      // AlphaUp
     1752        13.12344,  0.0,      0.0,      0.0,       0.0,       0.0,
     1753        0.0,       0.0                                                 };
     1754
     1755      Double_t vsteps[104] = {
     1756      // LengthUp
     1757        0.03,      0.0002,   0.02,     0.0006,    0.0002,    0.002,
     1758        0.0008,    0.002,
     1759      // LengthLo
     1760        0.02,      0.003,    0.05,     0.006,     0.002,     0.3,
     1761        0.0001,    0.0008,
     1762      // WidthUp
     1763        0.02,      0.0002,   0.006,    0.003,     0.002,     0.006,
     1764        0.002,     0.002,
     1765      // WidthLo
     1766        0.009,     0.0007,   0.008,    0.0004,    0.0005,    0.002,
     1767        0.0007,    0.003,
     1768      // DistUp
     1769        0.2,       0.0,      0.3,      0.02,      0.0,       0.03,
     1770        0.02,      0.0
     1771      // DistLo
     1772        0.06,      0.0,      0.009,    0.0008,    0.0,       0.005,
     1773        0.0002,    0.0
     1774      // AsymUp 
     1775        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1776        0.0,       0.0,
     1777      // AsymLo 
     1778        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1779        0.0,       0.0,
     1780      // ConcUp 
     1781        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1782        0.0,       0.0,
     1783      // ConcLo 
     1784        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1785        0.0,       0.0,
     1786      // Leakage1Up 
     1787        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1788        0.0,       0.0,
     1789      // Leakage1Lo 
     1790        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1791        0.0,       0.0,
     1792      // AlphaUp 
     1793        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1794        0.0,       0.0                                                 };
     1795
     1796      params.Set(104, vparams);
     1797      steps.Set (104, vsteps );
     1798    }
     1799
     1800    Bool_t rf;
     1801    rf = findsuper.FindParams(parSCinit, params, steps);
    15491802
    15501803    if (!rf)
    15511804    {
    1552       *fLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
     1805       gLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
    15531806       return;
    15541807    }
    … …  
    15651818    //
    15661819
    1567 
     1820 if (WSC)
     1821 {
    15681822    gLog << "" << endl;
    15691823    gLog << "========================================================" << endl;
    … …  
    15801834    inparam.Close();
    15811835
    1582     gLog << "Optimum parameter values for supercuts were read in from file '"
     1836    gLog << "Parameter values for supercuts were read in from file '"
    15831837         << parSCfile << "'" << endl;
    15841838
    … …  
    15861840    supercutsPar =  scin.GetParameters();
    15871841
    1588     gLog << "Optimum parameter values for supercuts : " << endl;
    1589     for (Int_t i=0; i<72; i++)
     1842    TArrayD supercutsStep;
     1843    supercutsStep =  scin.GetStepsizes();
     1844
     1845    gLog << "Parameter values for supercuts : " << endl;
     1846    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
    15901847    {
    15911848      gLog << supercutsPar[i] << ",  ";
     1849    }
     1850    gLog << endl;
     1851
     1852    gLog << "Step values for supercuts : " << endl;
     1853    for (Int_t i=0; i<supercutsStep.GetSize(); i++)
     1854    {
     1855      gLog << supercutsStep[i] << ",  ";
    15921856    }
    15931857    gLog << endl;
    … …  
    16041868         << filenameData << "'" << endl;
    16051869    supercutsPar = supercuts.GetParameters();
    1606     for (Int_t i=0; i<72; i++)
     1870    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
    16071871    {
    16081872      gLog << supercutsPar[i] << ",  ";
    … …  
    16401904
    16411905
    1642  if (WSC)
    1643   {
    16441906      //MWriteRootFile write(outNameImage, "UPDATE");
    1645       MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
     1907      //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
     1908      MWriteRootFile write(outNameImage, "RECREATE");
    16461909
    16471910      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    16561919      write.AddContainer("MNewImagePar",  "Events");
    16571920
    1658       write.AddContainer("HadNN",         "Events");
    1659       write.AddContainer("HadSC",         "Events");
    1660   }
     1921      write.AddContainer("HadRF",         "Events");
     1922      write.AddContainer(hadSCName,       "Events");
    16611923
    16621924
    … …  
    16931955    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
    16941956    alphaabs.SetName(mh3name);
     1957
     1958    TH1  &alphahist = alphaabs->GetHist();
     1959
    16951960    MFillH alpha(&alphaabs);
    16961961    alpha.SetName("FillAlphaAbs");
    … …  
    17161981    pliston.AddToList(&tliston);
    17171982    InitBinnings(&pliston);
     1983
     1984    pliston.AddToList(&supercuts);
     1985
    17181986    pliston.AddToList(&binsalphaabs);
    17191987    pliston.AddToList(&alphaabs);
    1720     pliston.AddToList(&supercuts);
    17211988
    17221989    //*****************************
    … …  
    17271994    tliston.AddToList(&sccalc);
    17281995    tliston.AddToList(&fillhadsc);
     1996
     1997    tliston.AddToList(&write);
    17291998    tliston.AddToList(&contfinalgh);
    17301999
    … …  
    17352004    tliston.AddToList(&hfill4);
    17362005    tliston.AddToList(&hfill5);
    1737 
    1738     if (WSC)
    1739       tliston.AddToList(&write);
    17402006
    17412007    tliston.AddToList(&contfinal);
    … …  
    17762042    TH1  &alphaHist = absalpha->GetHist();
    17772043    alphaHist.SetXTitle("|alpha|  [\\circ]");
     2044    alphaHist.SetName("alpha-macro");
    17782045
    17792046    Double_t alphasig = 13.1;
    17802047    Double_t alphamin = 30.0;
    17812048    Double_t alphamax = 90.0;
    1782     Int_t    degree   =    4;
     2049    Int_t    degree   =    2;
    17832050    Double_t significance = -99.0;
    17842051    Bool_t   drawpoly  = kTRUE;
    … …  
    17882055    MHFindSignificance findsig;
    17892056    findsig.SetRebin(kTRUE);
    1790     findsig.SetReduceDegree(kTRUE);
     2057    findsig.SetReduceDegree(kFALSE);
    17912058
    17922059    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
    … …  
    18052072
    18062073    DeleteBinnings(&pliston);
     2074 }
     2075
    18072076
    18082077    gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl;
    … …  
    18102079 }
    18112080  //---------------------------------------------------------------------
    1812 
    1813 
    1814   //---------------------------------------------------------------------
    1815   // Job B_RF_UP
    1816   //============
    1817 
    1818 
    1819     //  - create (or read in) the matrices of look-alike events for gammas
    1820     //    and hadrons
    1821     //  - create (or read in) the trees
    1822     //  - then read ON1.root (or MC1.root) file
    1823     //  - calculate the hadroness for the method of RANDOM FOREST
    1824     //  - update input root file with the hadroness
    1825 
    1826 
    1827  if (JobB_RF_UP)
    1828  {
    1829     gLog << "=====================================================" << endl;
    1830     gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
    1831 
    1832     gLog << "" << endl;
    1833     gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = "
    1834          << JobB_RF_UP  << ",  " << RLookRF << ",  " << RTree << ",  "
    1835          << WRF << endl;
    1836 
    1837 
    1838 
    1839     //--------------------------------------------
    1840     // parameters for the random forest
    1841     Int_t NumTrees = 100;
    1842     Int_t NumTry   =   3;
    1843     Int_t NdSize   =   3;
    1844 
    1845 
    1846     //--------------------------------------------
    1847     // file to be updated (either ON or MC)
    1848 
    1849     TString typeInput = "ON";
    1850     //TString typeInput = "MC";
    1851     gLog << "typeInput = " << typeInput << endl;
    1852 
    1853     // name of input root file
    1854     TString filenameData = outPath;
    1855     filenameData += typeInput;
    1856     filenameData += "2.root";
    1857     gLog << "filenameData = " << filenameData << endl;
    1858 
    1859     // name of output root file
    1860     TString outNameImage = outPath;
    1861     outNameImage += typeInput;
    1862     outNameImage += "3.root";
    1863     //TString outNameImage = filenameData;
    1864 
    1865     gLog << "outNameImage = " << outNameImage << endl;
    1866 
    1867     //--------------------------------------------
    1868     // files to be read for generating the look-alike events
    1869     // "hadrons" :
    1870     TString filenameON = outPath;
    1871     filenameON += "ON";
    1872     filenameON += "1.root";
    1873     Int_t howManyHadrons = 1000000;
    1874     gLog << "filenameON = " << filenameON << ",   howManyHadrons = "
    1875          << howManyHadrons  << endl;
    1876    
    1877 
    1878     // "gammas" :
    1879     TString filenameMC = outPath;
    1880     filenameMC += "MC";
    1881     filenameMC += "1.root";
    1882     Int_t howManyGammas = 50000;
    1883     gLog << "filenameMC = " << filenameMC << ",   howManyGammas = "
    1884          << howManyGammas   << endl;
    1885    
    1886     //--------------------------------------------
    1887     // files of look-alike events
    1888 
    1889     TString outNameGammas = outPath;
    1890     outNameGammas += "RFmatrix_gammas_";
    1891     outNameGammas += "MC";
    1892     outNameGammas += ".root";
    1893 
    1894     TString typeMatrixHadrons = "ON";
    1895     gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
    1896 
    1897     TString outNameHadrons = outPath;
    1898     outNameHadrons += "RFmatrix_hadrons_";
    1899     outNameHadrons += typeMatrixHadrons;
    1900     outNameHadrons += ".root";
    1901 
    1902 
    1903     MHMatrix matrixg("MatrixGammas");
    1904     MHMatrix matrixh("MatrixHadrons");
    1905 
    1906     //--------------------------------------------
    1907     // file of trees of the random forest
    1908 
    1909     TString outRF = outPath;
    1910     outRF += "RF.root";
    1911 
    1912 
    1913    //*************************************************************************
    1914    // read in matrices of look-alike events
    1915 if (RLookRF)
    1916   {
    1917     const char* mtxName = "MatrixGammas";
    1918 
    1919     gLog << "" << endl;
    1920     gLog << "========================================================" << endl;
    1921     gLog << "Get matrix for (gammas)" << endl;
    1922     gLog << "matrix name        = " << mtxName << endl;
    1923     gLog << "name of root file  = " << outNameGammas << endl;
    1924     gLog << "" << endl;
    1925 
    1926 
    1927     // read in the object with the name 'mtxName' from file 'outNameGammas'
    1928     //
    1929     TFile fileg(outNameGammas);
    1930 
    1931     matrixg.Read(mtxName);
    1932     matrixg.Print("cols");
    1933 
    1934 
    1935     //*****************************************************************
    1936 
    1937     const char* mtxName = "MatrixHadrons";
    1938 
    1939     gLog << "" << endl;
    1940     gLog << "========================================================" << endl;
    1941     gLog << " Get matrix for (hadrons)" << endl;
    1942     gLog << "matrix name        = " << mtxName << endl;
    1943     gLog << "name of root file  = " << outNameHadrons << endl;
    1944     gLog << "" << endl;
    1945 
    1946 
    1947     // read in the object with the name 'mtxName' from file 'outNameHadrons'
    1948     //
    1949     TFile fileh(outNameHadrons);
    1950 
    1951     matrixh.Read(mtxName);
    1952     matrixh.Print("cols");
    1953   }
    1954 
    1955 
    1956    //*************************************************************************
    1957    // create matrices of look-alike events
    1958 if (!RLookRF)
    1959   {
    1960     gLog << "" << endl;
    1961     gLog << "========================================================" << endl;
    1962     gLog << " Create matrices of look-alike events" << endl;
    1963     gLog << " Gammas :" << endl;
    1964 
    1965 
    1966     MParList  plistg;
    1967     MTaskList tlistg;
    1968     MFilterList flistg;
    1969 
    1970     MParList  plisth;
    1971     MTaskList tlisth;
    1972     MFilterList flisth;
    1973 
    1974     MReadMarsFile  readg("Events", filenameMC);
    1975     readg.DisableAutoScheme();
    1976 
    1977     MReadMarsFile  readh("Events", filenameON);
    1978     readh.DisableAutoScheme();
    1979 
    1980     MFParticleId fgamma("MMcEvt", '=', kGAMMA);
    1981     fgamma.SetName("gammaID");
    1982     MFParticleId fhadrons("MMcEvt", '!', kGAMMA);
    1983     fhadrons.SetName("hadronID)");
    1984 
    1985     MFEventSelector selectorg;
    1986     selectorg.SetNumSelectEvts(howManyGammas);
    1987     selectorg.SetName("selectGammas");
    1988     MFEventSelector selectorh;
    1989     selectorh.SetNumSelectEvts(howManyHadrons);
    1990     selectorh.SetName("selectHadrons");
    1991 
    1992     MFillH fillmatg("MatrixGammas");
    1993     fillmatg.SetFilter(&flistg);
    1994     fillmatg.SetName("fillGammas");
    1995 
    1996     MFillH fillmath("MatrixHadrons");
    1997     fillmath.SetFilter(&flisth);
    1998     fillmath.SetName("fillHadrons");
    1999 
    2000    
    2001     //*****************************   fill gammas   *** 
    2002     // entries in MFilterList
    2003 
    2004     flistg.AddToList(&fgamma);
    2005     flistg.AddToList(&selectorg);
    2006 
    2007     //***************************** 
    2008     // entries in MParList
    2009    
    2010     plistg.AddToList(&tlistg);
    2011     InitBinnings(&plistg);
    2012 
    2013     plistg.AddToList(&matrixg);
    2014 
    2015     //*****************************
    2016     // entries in MTaskList
    2017    
    2018     tlistg.AddToList(&readg);
    2019     tlistg.AddToList(&flistg);
    2020     tlistg.AddToList(&fillmatg);
    2021 
    2022     //*****************************
    2023 
    2024     MProgressBar matrixbar;
    2025     MEvtLoop evtloopg;
    2026     evtloopg.SetParList(&plistg);
    2027     evtloopg.ReadEnv(env, "", printEnv);
    2028     evtloopg.SetProgressBar(&matrixbar);
    2029 
    2030     Int_t maxevents = -1;
    2031     if (!evtloopg.Eventloop(maxevents))
    2032         return;
    2033 
    2034     tlistg.PrintStatistics(0, kTRUE);
    2035 
    2036 
    2037     //*****************************   fill hadrons   *** 
    2038 
    2039     gLog << " Hadrons :" << endl;
    2040     // entries in MFilterList
    2041 
    2042     flisth.AddToList(&fhadrons);
    2043     flisth.AddToList(&selectorh);
    2044 
    2045     //***************************** 
    2046     // entries in MParList
    2047    
    2048     plisth.AddToList(&tlisth);
    2049     InitBinnings(&plisth);
    2050 
    2051     plisth.AddToList(&matrixh);
    2052 
    2053     //*****************************
    2054     // entries in MTaskList
    2055    
    2056     tlisth.AddToList(&readh);
    2057     tlisth.AddToList(&flisth);
    2058     tlisth.AddToList(&fillmath);
    2059 
    2060     //*****************************
    2061 
    2062     MProgressBar matrixbar;
    2063     MEvtLoop evtlooph;
    2064     evtlooph.SetParList(&plisth);
    2065     evtlooph.ReadEnv(env, "", printEnv);
    2066     evtlooph.SetProgressBar(&matrixbar);
    2067 
    2068     Int_t maxevents = -1;
    2069     if (!evtlooph.Eventloop(maxevents))
    2070         return;
    2071 
    2072     tlisth.PrintStatistics(0, kTRUE);
    2073     //***************************************************** 
    2074 
    2075     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    2076     //
    2077     // ----------------------------------------------------------
    2078     //  Definition of the reference sample of look-alike events
    2079     //  (this is a subsample of the original sample)
    2080     // ----------------------------------------------------------
    2081     //
    2082     gLog << "" << endl;
    2083     gLog << "========================================================" << endl;
    2084     // Select a maximum of nmaxevts events from the sample of look-alike
    2085     // events. They will form the reference sample.
    2086     Int_t nmaxevts  = 10000;
    2087 
    2088     // target distribution for the variable in column refcolumn (start at 0);
    2089     //Int_t   refcolumn = 0;
    2090     //Int_t   nbins   =   5;
    2091     //Float_t frombin = 0.5;
    2092     //Float_t tobin   = 1.0;
    2093     //TH1F *thsh = new TH1F("thsh","target distribution",
    2094     //                       nbins, frombin, tobin);
    2095     //thsh->SetDirectory(NULL);
    2096     //thsh->SetXTitle("cos( \\Theta)");
    2097     //thsh->SetYTitle("Counts");
    2098     //Float_t dbin = (tobin-frombin)/nbins;
    2099     //Float_t lbin = frombin +dbin*0.5;
    2100     //for (Int_t j=1; j<=nbins; j++)
    2101     //{
    2102     //  thsh->Fill(lbin,1.0);
    2103     //  lbin += dbin;
    2104     //}
    2105 
    2106     Int_t   refcolumn = 0;
    2107     MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning");
    2108     TH1F *thsh = new TH1F();
    2109     thsh->SetNameTitle("thsh","target distribution");
    2110     MH::SetBinning(thsh, binscostheta);
    2111     Int_t nbins = thsh->GetNbinsX();
    2112     Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900};
    2113     for (Int_t j=1; j<=nbins; j++)
    2114     {
    2115       thsh->Fill(thsh->GetBinCenter(j), cont[j-1]);
    2116     }
    2117 
    2118     gLog << "" << endl;
    2119     gLog << "========================================================" << endl;
    2120     gLog << "Macro CT1Analysis : define reference sample for gammas" << endl;
    2121     gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
    2122          << refcolumn << ",  " << nmaxevts << endl;   
    2123 
    2124     matrixg.EnableGraphicalOutput();
    2125     Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    2126 
    2127     if ( !matrixok )
    2128     {
    2129       gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
    2130       return;
    2131     }
    2132 
    2133     gLog << "" << endl;
    2134     gLog << "========================================================" << endl;
    2135     gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl;
    2136     gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
    2137          << refcolumn << ",  " << nmaxevts << endl;   
    2138 
    2139     matrixh.EnableGraphicalOutput();
    2140     matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    2141     delete thsh;
    2142 
    2143     if ( !matrixok )
    2144     {
    2145       gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
    2146       return;
    2147     }
    2148     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    2149 
    2150     // write out look-alike events
    2151 
    2152     gLog << "" << endl;
    2153     gLog << "========================================================" << endl;
    2154     gLog << "Write out look=alike events" << endl;
    2155 
    2156 
    2157       //-------------------------------------------
    2158       // "gammas"
    2159       gLog << "Gammas :" << endl;   
    2160       matrixg.Print("cols");
    2161 
    2162       TFile writeg(outNameGammas, "RECREATE", "");
    2163       matrixg.Write();
    2164 
    2165       gLog << "" << endl;
    2166       gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
    2167            << outNameGammas << endl;
    2168 
    2169       //-------------------------------------------
    2170       // "hadrons"
    2171       gLog << "Hadrons :" << endl;   
    2172       matrixh.Print("cols");
    2173 
    2174       TFile writeh(outNameHadrons, "RECREATE", "");
    2175       matrixh.Write();
    2176 
    2177       gLog << "" << endl;
    2178       gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
    2179            << outNameHadrons << endl;
    2180 
    2181   }
    2182    //**********   end of creating matrices of look-alike events   ***********
    2183 
    2184 
    2185     MRanForest *fRanForest;
    2186     MRanTree *fRanTree;
    2187     //-----------------------------------------------------------------
    2188     // read in the trees of the random forest
    2189     if (RTree)
    2190     {
    2191       MParList plisttr;
    2192       MTaskList tlisttr;
    2193       plisttr.AddToList(&tlisttr);
    2194 
    2195       MReadTree readtr("TREE", outRF);
    2196       readtr.DisableAutoScheme();
    2197 
    2198       MRanForestFill rffill;
    2199       rffill.SetNumTrees(NumTrees);
    2200 
    2201       // list of tasks for the loop over the trees
    2202 
    2203       tlisttr.AddToList(&readtr);
    2204       tlisttr.AddToList(&rffill);
    2205 
    2206       //-------------------
    2207       // Execute tree loop
    2208       //
    2209       MEvtLoop evtlooptr;
    2210       evtlooptr.SetParList(&plisttr);
    2211       if (!evtlooptr.Eventloop())
    2212         return;
    2213 
    2214       tlisttr.PrintStatistics(0, kTRUE);
    2215 
    2216 
    2217     // get adresses of objects which are used in the next eventloop
    2218     fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
    2219     if (!fRanForest)
    2220     {
    2221         *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
    2222         return kFALSE;
    2223     }
    2224 
    2225     fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
    2226     if (!fRanTree)                                 
    2227     {                                                                         
    2228         *fLog << err << dbginf << "MRanTree not found... aborting." << endl;   
    2229         return kFALSE;
    2230     }
    2231 
    2232     }
    2233 
    2234     //-----------------------------------------------------------------
    2235     // grow the trees of the random forest (event loop = tree loop)
    2236 
    2237     if (!RTree)
    2238     {
    2239 
    2240     gLog << "" << endl;
    2241     gLog << "========================================================" << endl;
    2242     gLog << "Macro CT1Analysis : start growing trees" << endl;
    2243 
    2244     MTaskList tlist2;
    2245     MParList plist2;
    2246     plist2.AddToList(&tlist2);
    2247 
    2248     plist2.AddToList(&matrixg);
    2249     plist2.AddToList(&matrixh);
    2250 
    2251     MRanForestGrow rfgrow2;
    2252     rfgrow2.SetNumTrees(NumTrees);
    2253     rfgrow2.SetNumTry(NumTry);
    2254     rfgrow2.SetNdSize(NdSize);
    2255 
    2256     MWriteRootFile rfwrite2(outRF);
    2257     rfwrite2.AddContainer("MRanTree", "TREE");
    2258 
    2259     // list of tasks for the loop over the trees
    2260    
    2261     tlist2.AddToList(&rfgrow2);
    2262     tlist2.AddToList(&rfwrite2);
    2263 
    2264     //-------------------
    2265     // Execute tree loop
    2266     //
    2267     MEvtLoop treeloop;
    2268     treeloop.SetParList(&plist2);
    2269 
    2270     if ( !treeloop.Eventloop() )
    2271         return;
    2272 
    2273     tlist2.PrintStatistics(0, kTRUE);
    2274 
    2275     // get adresses of objects which are used in the next eventloop
    2276     fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
    2277     if (!fRanForest)
    2278     {
    2279         *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
    2280         return kFALSE;
    2281     }
    2282 
    2283     fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
    2284     if (!fRanTree)                                 
    2285     {                                                                         
    2286         *fLog << err << dbginf << "MRanTree not found... aborting." << endl;   
    2287         return kFALSE;
    2288     }
    2289 
    2290     }
    2291     // end of growing the trees of the random forest
    2292     //-----------------------------------------------------------------
    2293 
    2294 
    2295 
    2296 
    2297     //-----------------------------------------------------------------
    2298     // Update the input files with the RF hadronness
    2299     //
    2300 
    2301  if (WRF)
    2302   {
    2303     gLog << "" << endl;
    2304     gLog << "========================================================" << endl;
    2305     gLog << "Update input file '" <<  filenameData
    2306          << "' with the RF hadronness" << endl;
    2307 
    2308     MTaskList tliston;
    2309     MParList pliston;
    2310 
    2311 
    2312     // geometry is needed in  MHHillas... classes
    2313     MGeomCam *fGeom =
    2314              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    2315 
    2316     //-------------------------------------------
    2317     // create the tasks which should be executed
    2318     //
    2319 
    2320     MReadMarsFile read("Events", filenameData);
    2321     read.DisableAutoScheme();
    2322 
    2323     TString fHilName    = "MHillas";
    2324     TString fHilNameExt = "MHillasExt";
    2325     TString fHilNameSrc = "MHillasSrc";
    2326     TString fImgParName = "MNewImagePar";
    2327 
    2328     //.......................................................................
    2329     // calculate hadronnes for method of RANDOM FOREST
    2330 
    2331     TString hadRFName = "HadRF";
    2332     MRanForestCalc rfcalc;
    2333     rfcalc.SetHadronnessName(hadRFName);
    2334 
    2335     //.......................................................................
    2336 
    2337       //MWriteRootFile write(outNameImage, "UPDATE");
    2338       MWriteRootFile write(outNameImage, "RECREATE");
    2339 
    2340       write.AddContainer("MRawRunHeader", "RunHeaders");
    2341       write.AddContainer("MTime",         "Events");
    2342       write.AddContainer("MMcEvt",        "Events");
    2343       write.AddContainer("ThetaOrig",     "Events");
    2344       write.AddContainer("MSrcPosCam",    "Events");
    2345       write.AddContainer("MSigmabar",     "Events");
    2346       write.AddContainer("MHillas",       "Events");
    2347       write.AddContainer("MHillasExt",    "Events");
    2348       write.AddContainer("MHillasSrc",    "Events");
    2349       write.AddContainer("MNewImagePar",  "Events");
    2350 
    2351       write.AddContainer("HadNN",         "Events");
    2352       write.AddContainer("HadSC",         "Events");
    2353 
    2354       write.AddContainer(hadRFName,       "Events");
    2355 
    2356 
    2357     //-----------------------------------------------------------------
    2358     // geometry is needed in  MHHillas... classes
    2359     MGeomCam *fGeom =
    2360              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    2361 
    2362 
    2363     Float_t maxhadronness =  0.40;
    2364     Float_t maxalpha      =  20.0;
    2365     Float_t maxdist       =  10.0;
    2366 
    2367     MFCT1SelFinal selfinalgh(fHilNameSrc);
    2368     selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    2369     selfinalgh.SetHadronnessName(hadRFName);
    2370     selfinalgh.SetName("SelFinalgh");
    2371     MContinue contfinalgh(&selfinalgh);
    2372     contfinalgh.SetName("ContSelFinalgh");
    2373 
    2374     MFillH fillranfor("MHRanForest");
    2375     fillranfor.SetName("HRanForest");
    2376 
    2377     MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
    2378     fillhadrf.SetName("HhadRF");
    2379 
    2380     MFCT1SelFinal selfinal(fHilNameSrc);
    2381     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    2382     selfinal.SetHadronnessName(hadRFName);
    2383     selfinal.SetName("SelFinal");
    2384     MContinue contfinal(&selfinal);
    2385     contfinal.SetName("ContSelFinal");
    2386 
    2387 
    2388     MFillH hfill1("MHHillas",    fHilName);
    2389     hfill1.SetName("HHillas");
    2390 
    2391     MFillH hfill2("MHStarMap",   fHilName);
    2392     hfill2.SetName("HStarMap");
    2393 
    2394     MFillH hfill3("MHHillasExt",    fHilNameSrc);
    2395     hfill3.SetName("HHillasExt");
    2396    
    2397     MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    2398     hfill4.SetName("HHillasSrc");   
    2399 
    2400     MFillH hfill5("MHNewImagePar", fImgParName);
    2401     hfill5.SetName("HNewImagePar");
    2402 
    2403     //*****************************
    2404     // entries in MParList
    2405 
    2406     pliston.AddToList(&tliston);
    2407     InitBinnings(&pliston);
    2408 
    2409     pliston.AddToList(fRanForest);
    2410     pliston.AddToList(fRanTree);
    2411 
    2412     //*****************************
    2413     // entries in MTaskList
    2414    
    2415     tliston.AddToList(&read);
    2416 
    2417     tliston.AddToList(&rfcalc);
    2418     tliston.AddToList(&fillranfor);
    2419     tliston.AddToList(&fillhadrf);
    2420 
    2421     tliston.AddToList(&hfill1);
    2422     tliston.AddToList(&hfill2);
    2423     tliston.AddToList(&hfill3);
    2424     tliston.AddToList(&hfill4);
    2425     tliston.AddToList(&hfill5);
    2426 
    2427     tliston.AddToList(&write);
    2428 
    2429     tliston.AddToList(&contfinalgh);
    2430     tliston.AddToList(&contfinal);
    2431 
    2432     //*****************************
    2433 
    2434     //-------------------------------------------
    2435     // Execute event loop
    2436     //
    2437     MProgressBar bar;
    2438     MEvtLoop evtloop;
    2439     evtloop.SetParList(&pliston);
    2440     evtloop.SetProgressBar(&bar);
    2441 
    2442     Int_t maxevents = -1;
    2443     if ( !evtloop.Eventloop(maxevents) )
    2444         return;
    2445 
    2446     tliston.PrintStatistics(0, kTRUE);
    2447 
    2448 
    2449     //-------------------------------------------
    2450     // Display the histograms
    2451     //
    2452     pliston.FindObject("MHRanForest")->DrawClone();
    2453     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    2454 
    2455     pliston.FindObject("MHHillas")->DrawClone();
    2456     pliston.FindObject("MHHillasExt")->DrawClone();
    2457     pliston.FindObject("MHHillasSrc")->DrawClone();
    2458     pliston.FindObject("MHNewImagePar")->DrawClone();
    2459     pliston.FindObject("MHStarMap")->DrawClone();
    2460 
    2461     DeleteBinnings(&pliston);
    2462   }
    2463 
    2464     gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
    2465     gLog << "=======================================================" << endl;
    2466  }
    2467   //---------------------------------------------------------------------
    2468 
    24692081
    24702082
    … …  
    24842096
    24852097    gLog << "" << endl;
    2486     gLog << "Macro CT1Analysis : JobC = " << JobC  << endl;
     2098    gLog << "Macro CT1Analysis : JobC = "
     2099         << (JobC       ? "kTRUE" : "kFALSE")  << endl;
    24872100
    24882101
    … …  
    25222135    // names of hadronness containers
    25232136
    2524     TString hadNNName = "HadNN";
     2137    //TString hadNNName = "HadNN";
    25252138    TString hadSCName = "HadSC";
    25262139    TString hadRFName = "HadRF";
    … …  
    25402153    MFCT1SelFinal selfinalgh(fHilNameSrc);
    25412154    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    2542     selfinalgh.SetHadronnessName(hadNNName);
     2155    selfinalgh.SetHadronnessName(hadSCName);
    25432156    selfinalgh.SetName("SelFinalgh");
    25442157    MContinue contfinalgh(&selfinalgh);
    25452158    contfinalgh.SetName("ContSelFinalgh");
    25462159
    2547     MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
    2548     fillhadnn.SetName("HhadNN");
     2160    //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
     2161    //fillhadnn.SetName("HhadNN");
    25492162    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
    25502163    fillhadsc.SetName("HhadSC");
    … …  
    25542167    MFCT1SelFinal selfinal(fHilNameSrc);
    25552168    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    2556     selfinal.SetHadronnessName(hadNNName);
     2169    selfinal.SetHadronnessName(hadSCName);
    25572170    selfinal.SetName("SelFinal");
    25582171    MContinue contfinal(&selfinal);
    … …  
    25882201    tliston.AddToList(&read);
    25892202
    2590     tliston.AddToList(&fillhadnn);
     2203    //tliston.AddToList(&fillhadnn);
    25912204    tliston.AddToList(&fillhadsc);
    25922205    tliston.AddToList(&fillhadrf);
    … …  
    26232236    //
    26242237
    2625     pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     2238    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
    26262239    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    26272240    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    … …  
    26572270    gLog << "" << endl;
    26582271    gLog << "Macro CT1Analysis : JobD = "
    2659          << JobD  << endl;
     2272         << (JobD        ? "kTRUE" : "kFALSE")  << endl;
     2273
    26602274
    26612275    // type of data to be analysed
    … …  
    28682482    gLog << "" << endl;
    28692483    gLog << "Macro CT1Analysis : JobE_XX, WXX = "
    2870          << JobE_XX  << ",  " << WXX << endl;
     2484         << (JobE_XX ? "kTRUE" : "kFALSE")  << ",  "
     2485         << (WXX     ? "kTRUE" : "kFALSE")  << endl;
     2486
    28712487
    28722488    // type of data to be analysed
    … …  
    30602676    gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
    30612677
     2678
     2679  if (WXX)
     2680  {
    30622681    //-----------------------------------------------------------
    30632682    //
    … …  
    31192738    //.......................................................................
    31202739
    3121     if (WXX)
    3122     {
    3123       MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
     2740
     2741      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
     2742      MWriteRootFile write(filenameDataout, "RECREATE");
    31242743
    31252744      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    31392758
    31402759      write.AddContainer("MEnergyEst",    "Events");
    3141     }
     2760
    31422761
    31432762    //-----------------------------------------------------------------
    … …  
    32652884
    32662885    DeleteBinnings(&pliston);
     2886
     2887  }
    32672888
    32682889    gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
    … …  
    32922913    gLog << "" << endl;
    32932914    gLog << "Macro CT1Analysis : JobF_XX, WFX = "
    3294          << JobF_XX  << ",  " << WFX << endl;
     2915         << (JobF_XX ? "kTRUE" : "kFALSE")  << ",  "
     2916         << (WFX     ? "kTRUE" : "kFALSE")  << endl;
     2917
    32952918
    32962919    // type of data to be analysed
    … …  
    34933116    gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
    34943117
     3118
     3119  if (WFX)
     3120  {
    34953121    //-----------------------------------------------------------
    34963122    //
    … …  
    35523178    //.......................................................................
    35533179
    3554     //if (WFX)
    3555     //{
    35563180      gLog << "CT1Analysis.C : write root file '" << filenameDataout
    35573181           << "'" << endl;
    35583182   
    35593183      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
     3184      /*
    35603185      MWriteRootFile write(filenameDataout, "RECREATE");
    35613186
    … …  
    35713196      write.AddContainer("MNewImagePar",  "Events");
    35723197
    3573       write.AddContainer("HadNN",         "Events");
     3198      //write.AddContainer("HadNN",         "Events");
    35743199      write.AddContainer("HadSC",         "Events");
    35753200      write.AddContainer("HadRF",         "Events");
    35763201
    35773202      write.AddContainer("MEnergyEst",    "Events");
    3578       //}
     3203      */
    35793204
    35803205    //-----------------------------------------------------------------
    … …  
    35903215    contfinalgh.SetName("ContSelFinalgh");
    35913216
    3592     MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
    3593     fillhadnn.SetName("HhadNN");
     3217    //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
     3218    //fillhadnn.SetName("HhadNN");
    35943219    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
    35953220    fillhadsc.SetName("HhadSC");
    … …  
    36003225    // calculate estimated energy
    36013226
    3602     //MEnergyEstParam eeston(fHilName);
    3603     //eeston.Add(fHilNameSrc);
    3604 
    3605     //eeston.SetCoeffA(parA);
    3606     //eeston.SetCoeffB(parB);
     3227    MEnergyEstParam eeston(fHilName);
     3228    eeston.Add(fHilNameSrc);
     3229
     3230    eeston.SetCoeffA(parA);
     3231    eeston.SetCoeffB(parB);
    36073232
    36083233    //---------------------------
    36093234    // calculate estimated energy using Daniel's parameters
    36103235
    3611     MEnergyEstParamDanielMkn421 eeston(fHilName);
    3612     eeston.Add(fHilNameSrc);
     3236    //MEnergyEstParamDanielMkn421 eeston(fHilName);
     3237    //eeston.Add(fHilNameSrc);
    36133238    //eeston.SetCoeffA(parA);
    36143239    //eeston.SetCoeffB(parB);
    … …  
    36843309    tliston.AddToList(&eeston);
    36853310
    3686     if (WFX)
    3687       tliston.AddToList(&write);
    3688 
    3689     tliston.AddToList(&fillhadnn);
     3311    //tliston.AddToList(&write);
     3312
     3313    //tliston.AddToList(&fillhadnn);
    36903314    tliston.AddToList(&fillhadsc);
    36913315    tliston.AddToList(&fillhadrf);
    … …  
    37273351    //
    37283352
    3729     pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     3353    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     3354
     3355    gLog << "before hadRF" << endl;
    37303356    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
     3357
     3358    gLog << "before hadSC" << endl;
    37313359    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    37323360
     3361    gLog << "before MHHillas" << endl;
    37333362    pliston.FindObject("MHHillas")->DrawClone();
     3363
     3364    gLog << "before MHHillasExt" << endl;
    37343365    pliston.FindObject("MHHillasExt")->DrawClone();
     3366
     3367    gLog << "before MHHillasSrc" << endl;
    37353368    pliston.FindObject("MHHillasSrc")->DrawClone();
     3369
     3370    gLog << "before MHNewImagePar" << endl;
    37363371    pliston.FindObject("MHNewImagePar")->DrawClone();
     3372
     3373    gLog << "before MHStarMap" << endl;
    37373374    pliston.FindObject("MHStarMap")->DrawClone();
    37383375
     3376    gLog << "before DeleteBinnings" << endl;
     3377
    37393378    //DeleteBinnings(&pliston);
    37403379
     3380    gLog << "before Robert's code" << endl;
     3381
    37413382
    37423383//rwagner write all relevant histograms onto a file
    37433384
     3385  if (WRobert)
     3386  {
    37443387    gLog << "=======================================================" << endl;
    37453388    gLog << "Write results onto file '" << filenameResults << "'" << endl;
    … …  
    37513394        TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
    37523395    alphaHist->Write();
    3753    
     3396    gLog << "Alpha plot has been written out" << endl;   
     3397
     3398
    37543399    MHAlphaEnergyTheta* aetH =
    37553400      (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
    … …  
    37573402    aetHist->SetName("aetHist");
    37583403    aetHist->Write();
     3404    gLog << "AlphaEnergyTheta plot has been written out" << endl;   
    37593405
    37603406    MHAlphaEnergyTime* aetH2 =
    … …  
    37643410//     aetHist2->DrawClone();
    37653411    aetHist2->Write();
     3412    gLog << "AlphaEnergyTime plot has been written out" << endl;   
    37663413
    37673414    MHThetabarTime* tbt =
    … …  
    37703417    tbtHist->SetName("tbtHist");
    37713418    tbtHist->Write();
     3419    gLog << "ThetabarTime plot has been written out" << endl;   
    37723420
    37733421    MHEnergyTime* ent =
    … …  
    37763424    entHist->SetName("entHist");
    37773425    entHist->Write();
     3426    gLog << "EnergyTime plot has been written out" << endl;   
    37783427   
    37793428    MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
    … …  
    37823431    timeHist->SetTitle("Time diffs");
    37833432    timeHist->Write();
     3433    gLog << "TimeDiffTheta plot has been written out" << endl;   
     3434
    37843435
    37853436    MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
    … …  
    37883439    timeHist2->SetTitle("Time diffs");
    37893440    timeHist2->Write();
     3441    gLog << "TimeDiffTime plot has been written out" << endl;   
    37903442
    37913443//rwagner write also collareas to same file
    … …  
    37933445    collarea->GetHAll()->Write();
    37943446    collarea->GetHSel()->Write();
    3795    
     3447    gLog << "Effective collection areas have been written out" << endl;       
     3448
    37963449//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data
    37973450//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file
    37983451
     3452    gLog << "before closing outfile" << endl;
     3453
    37993454    //outfile.Close();
    3800 
     3455    gLog << "Results were written onto file '" << filenameResults
     3456         << "'" << endl;
    38013457    gLog << "=======================================================" << endl;
    3802 
     3458  }
     3459
     3460  }
    38033461
    38043462    gLog << "Macro CT1Analysis : End of Job F_XX" << endl;

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

    r2272 r2357  
    3333#include "MHMatrix.h"
    3434#include "MEnergyEstParam.h"
    35 #include "MEnergyEstParamDanielMkn421.h"
     35//#include "MEnergyEstParamDanielMkn421.h"
    3636#include "MMatrixLoop.h"
    3737#include "MChisqEval.h"
    … …  
    206206  //
    207207
    208   //MEnergyEstParam eest2(hilName);
     208  MEnergyEstParam eest2(hilName);
     209  eest2.Add(hilSrcName);
     210
     211  eest2.SetCoeffA(parA);
     212  eest2.SetCoeffB(parB);
     213
     214  // estimate energy using Daniel's parameters
     215  //MEnergyEstParamDanielMkn421 eest2(hilName);
    209216  //eest2.Add(hilSrcName);
    210 
    211   //eest2.SetCoeffA(parA);
    212   //eest2.SetCoeffB(parB);
    213 
    214   // estimate energy using Daniel's parameters
    215   MEnergyEstParamDanielMkn421 eest2(hilName);
    216   eest2.Add(hilSrcName);
    217217
    218218
Note: See TracChangeset for help on using the changeset viewer.