Changeset 2475 for trunk/MagicSoft


Ignore:
Timestamp:
11/05/03 15:11:25 (21 years ago)
Author:
wittek
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
8 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r2474 r2475  
    11                                                 -*-*- END OF LINE -*-*-
     2
     3  2003/11/05: Wolfgang Wittek
     4
     5  * macros/CT1Analysis.C
     6          /ONOFFCT1Analysis.C
     7    - current versions of the macros for the analysis of CT1 data
     8
     9  * manalysis/MMarquardt.[h,cc]
     10    - very pleliminary version of a class performing a minimization
     11      using the Marquardt method
     12
     13  * mimage/M2dimFunctionFit.[h,cc]
     14    - very pleliminary version of a class which fits a 2-dim function to the
     15      shower image using the maximum likelihood method
     16
     17  * mimage/M2dimFunction.[h,cc]
     18    - very pleliminary version of a container which contains the parameters
     19      of the 2-dim function describing the shower image
     20
     21  * mimage/MH2dimFunction.[h,cc]
     22    - very pleliminary version of a container holding the histograms
     23      for the parameters of the 2-dim function describing the shower image
     24
     25
     26
    227 2003/11/05: Marcos Lopez
    328 
    … …  
    205230   * macros/weights.C
    206231      - Added macro showing how to transform the spectrum of the MC showers.
     232
     233
    207234
    208235

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

    r2437 r2475  
    241241
    242242    Bool_t JobB_RF_UP  = kFALSE; 
     243    Bool_t CTrainRF    = kFALSE;  // create matrices of training events
    243244    Bool_t RTrainRF    = kFALSE;  // read in matrices of training events
    244     Bool_t CTrainRF    = kFALSE;  // create matrices of training events
    245     Bool_t RTree       = kFALSE;  // read in trees
     245    Bool_t RTree       = kFALSE;  // read in trees (otherwise grow them)
    246246    Bool_t WRF         = kFALSE;  // update input root file ?
    247247
    … …  
    255255
    256256    Bool_t JobB_SC_UP  = kFALSE;
    257     Bool_t RMatrix     = kFALSE;  // read training and test matrices from file 
     257    Bool_t CMatrix     = kFALSE;  // create training and test matrices
     258    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file
    258259    Bool_t WOptimize   = kFALSE;  // do optimization using the training sample
    259260                                  // and write supercuts parameter values
    … …  
    12931294    rfwrite2.AddContainer("MRanTree", "TREE");
    12941295
     1296    MFillH fillh2("MHRanForestGini");
     1297
    12951298    // list of tasks for the loop over the trees
    12961299   
    12971300    tlist2.AddToList(&rfgrow2);
    12981301    tlist2.AddToList(&rfwrite2);
     1302    tlist2.AddToList(&fillh2);
    12991303
    13001304    //-------------------
    … …  
    13091313
    13101314    tlist2.PrintStatistics(0, kTRUE);
     1315
     1316    plist2.FindObject("MHRanForestGini")->DrawClone();
     1317
    13111318
    13121319    // get adresses of objects which are used in the next eventloop
    … …  
    15571564
    15581565    gLog << "" << endl;
    1559     gLog << "Macro CT1Analysis : JobB_SC_UP, RMatrix, WOptimize, RTest, WSC = "
     1566    gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = "
    15601567         << (JobB_SC_UP ? "kTRUE" : "kFALSE")  << ",  "
     1568         << (CMatrix    ? "kTRUE" : "kFALSE")  << ",  "
    15611569         << (RMatrix    ? "kTRUE" : "kFALSE")  << ",  "
    15621570         << (WOptimize  ? "kTRUE" : "kFALSE")  << ",  "
    … …  
    16601668    //--------------------------
    16611669    // create matrices and write them onto files
    1662     if (!RMatrix)
     1670    if (CMatrix)
    16631671    {
    16641672      MH3 &mh3 = *(new MH3("MHillas.fSize"));
    … …  
    21112119    //  - read ON1 and MC1 data files 
    21122120    //    which should have been updated to contain the hadronnesses
    2113     //    for the method of NEAREST NEIGHBORS and for the SUOERCUTS
     2121    //    for the method of Random Forest and for the SUPERCUTS
    21142122    //  - produce Neyman-Pearson plots
    21152123 

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

    r2437 r2475  
     1
    12
    23#include "CT1EgyEst.C"
    … …  
    78        gLog << "InitBinnings" << endl;
    89
     10        //--------------------------------------------
    911        MBinning *binse = new MBinning("BinningE");
    1012        //binse->SetEdgesLog(30, 1.0e2, 1.0e5);
    11         binse->SetEdges(30, 2, 5);
     13
     14        //This is Daniel's binning in energy:
     15        binse->SetEdgesLog(14, 296.296, 86497.6);
    1216        plist->AddToList(binse);
     17
     18        //--------------------------------------------
    1319
    1420        MBinning *binssize = new MBinning("BinningSize");
    … …  
    5157
    5258        MBinning *binth = new MBinning("BinningTheta");
     59        // this is Daniel's binning in theta
     60        //Double_t yedge[8] =
     61        //  {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99};
     62        // this is our binning
    5363        Double_t yedge[9] =
    5464                       {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
    55         TArrayD yed(9,yedge);
     65        TArrayD yed;
     66        yed.Set(9,yedge);
    5667        binth->SetEdges(yed);
    5768        plist->AddToList(binth);
    … …  
    6374          zedge[8-i] = cos(yedge[i]/kRad2Deg);
    6475        }
    65         TArrayD zed(9,zedge);
     76        TArrayD zed;
     77        zed.Set(9,zedge);
    6678        bincosth->SetEdges(zed);
    6779        plist->AddToList(bincosth);
    … …  
    7082        binsdiff->SetEdges(100, -5.0, 20.0);
    7183        plist->AddToList(binsdiff);
     84
     85        // robert ----------------------------------------------
     86        MBinning *binsalphaf = new MBinning("BinningAlphaFlux");
     87        binsalphaf->SetEdges(100, -100, 100);
     88        plist->AddToList(binsalphaf);
     89
     90        MBinning *binsdifftime = new MBinning("BinningTimeDiff");
     91        binsdifftime->SetEdges(50, 0, 10);
     92        plist->AddToList(binsdifftime);
     93
     94        MBinning *binstime = new MBinning("BinningTime");
     95        binstime->SetEdges(50, 44500, 61000);
     96        plist->AddToList(binstime);
    7297}
     98
    7399
    74100void DeleteBinnings(MParList *plist)
    … …  
    76102        gLog << "DeleteBinnings" << endl;
    77103
    78         if (!plist)
    79         {
    80           gLog << "Deletebinnins : MParlist no longer existing" << endl;
    81           return;
    82         }
    83 
    84104        TObject *bin;
     105
     106        //--------------------------------------------
     107        bin = plist->FindObject("BinningE");
     108        if (bin) delete bin;
     109
     110        //--------------------------------------------
    85111
    86112        bin = plist->FindObject("BinningSize");
    … …  
    121147        if (bin) delete bin;
    122148
    123         gLog << "exit DeleteBinnings" << endl;
     149
     150        // robert ----------------------------------------------
     151        bin = plist->FindObject("BinningAlphaFlux");
     152        if (bin) delete bin;
     153
     154        bin = plist->FindObject("BinningTimeDiff");
     155        if (bin) delete bin;
     156
     157        bin = plist->FindObject("BinningTime");
     158        if (bin) delete bin;
    124159}
     160
     161
    125162
    126163//************************************************************************
    … …  
    162199    //  - write root file for ON (or OFF or MC) data (ON1.root, ...);
    163200
    164     Bool_t JobA    = kTRUE; 
    165     Bool_t WPad    = kTRUE;   // write out padding histograms ?
     201    Bool_t JobA    = kFALSE; 
     202    Bool_t WPad    = kFALSE;   // write out padding histograms ?
    166203    Bool_t RPad    = kFALSE;   // read in padding histograms ?
    167204    Bool_t Wout    = kFALSE;   // write out root file ON1.root
    … …  
    178215
    179216    Bool_t JobB_RF_UP  = kFALSE; 
    180     Bool_t RTrainRF     = kFALSE;  // read in matrices of training events
    181     Bool_t CTrainRF     = kFALSE;  // create matrices of training events
    182     Bool_t RTree       = kFALSE;  // read in trees
     217    Bool_t CTrainRF    = kFALSE;  // create matrices of training events
     218    Bool_t RTrainRF    = kFALSE;  // read in matrices of training events
     219    Bool_t RTree       = kFALSE;  // read in trees (otherwise grow them)
    183220    Bool_t WRF         = kFALSE;  // update input root file ?
    184221
    185222
     223    // Job B_SC_UP : read ON2.root (or MC2.root) file
     224    //  - depending on WParSC : create (or read in) supercuts parameter values
     225    //  - calculate hadroness for the SUPERCUTS
     226    //  - update the input files with the hadroness (==>ON3.root or MC3.root)
     227
     228    Bool_t JobB_SC_UP  = kFALSE;
     229    Bool_t CMatrix     = kFALSE;  // create training and test matrices
     230    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file
     231    Bool_t WOptimize   = kFALSE;  // do optimization using the training sample
     232                                  // and write supercuts parameter values
     233                                  // onto the file parSCfile
     234    Bool_t RTest       = kFALSE;  // test the supercuts using the test matrix
     235    Bool_t WSC         = kFALSE;  // update input root file ?
     236
    186237
    187238    // Job C:
    188     //  - read ON1.root and MC1.root files
     239    //  - read ON3.root and MC3.root files
    189240    //    which should have been updated to contain the hadronnesses 
    190241    //    for the method of
    191     //              NEAREST NEIGHBORS 
     242    //              RF
    192243    //              SUPERCUTS and
    193     //              RF
    194244    //  - produce Neyman-Pearson plots
    195245
    … …  
    199249    // Job D : 
    200250    //  - select g/h separation method XX
    201     //  - read ON2 (or MC2) root file
     251    //  - read ON3 (or MC3) root file
    202252    //  - apply cuts in hadronness
    203253    //  - make plots
    … …  
    207257
    208258
    209     // Job E_EST_UP :
    210     //  - read MC1.root file
     259    // Job E_XX : extended version of E_XX (including flux plots) 
    211260    //  - select g/h separation method XX
    212     //  - optimize energy estimation for events passing the final cuts
    213     //  - write parameters of energy estimator onto file
    214     //  - update ON1.root, OFF1.root and MC1.root files with estimated energy
    215     //    (ON_XX1.root, OFF_XX1.root and MC_XX1.root)
    216 
    217     Bool_t JobE_EST_UP  = kFALSE; 
    218     Bool_t WESTUP       = kFALSE;  // update root files ?
    219 
    220 
    221 
    222     // Job F_XX : 
    223     //  - select g/h separation method XX
    224     //  - read MC_XX2.root file
     261    //  - read MC root file
    225262    //  - calculate eff. collection area
    226     //  - read ON_XX2.root file
     263    //  - optimize energy estimator
     264    //  - read ON root file
    227265    //  - apply final cuts
    228266    //  - calculate flux
    229     //  - write root file for ON data after final cuts (ON3.root))
    230 
    231     Bool_t JobF_XX  = kFALSE; 
    232     Bool_t WXX      = kFALSE;  // write out root file ON3.root ?
     267    //  - write root file for ON data after final cuts
     268
     269
     270    Bool_t JobE_XX  = kTRUE; 
     271    Bool_t OEEst    = kFALSE;  // optimize energy estimator
     272    Bool_t WEX      = kTRUE;  // update root file  ?
     273    Bool_t WRobert  = kTRUE;  // write out Robert's file  ?
     274
    233275
    234276
    … …  
    277319    //--------------------------------------------------
    278320    // type of data to be padded
    279     TString typeInput = "ON";
     321    //TString typeInput = "ON";
    280322    //TString typeInput = "OFF";
    281     //TString typeInput = "MC";
     323    TString typeInput = "MC";
    282324    gLog << "typeInput = " << typeInput << endl;
    283325
    … …  
    740782    // file to be updated (ON, OFF or MC)
    741783
    742     TString typeInput = "ON";
     784    //TString typeInput = "ON";
    743785    //TString typeInput = "OFF";
    744     //TString typeInput = "MC";
     786    TString typeInput = "MC";
    745787    gLog << "typeInput = " << typeInput << endl;
    746788
    … …  
    937979    MProgressBar matrixbar;
    938980    MEvtLoop evtloopg;
     981    evtloopg.SetName("FillGammaMatrix");
    939982    evtloopg.SetParList(&plistg);
    940983    //evtloopg.ReadEnv(env, "", printEnv);
    … …  
    9971040    MProgressBar matrixbar;
    9981041    MEvtLoop evtlooph;
     1042    evtlooph.SetName("FillHadronMatrix");
    9991043    evtlooph.SetParList(&plisth);
    10001044    //evtlooph.ReadEnv(env, "", printEnv);
    … …  
    10091053
    10101054
    1011     // write out matrices of training events events
     1055    // write out matrices of training events
    10121056
    10131057    gLog << "" << endl;
    … …  
    11191163    rfwrite2.AddContainer("MRanTree", "TREE");
    11201164
     1165    MFillH fillh2("MHRanForestGini");
     1166
    11211167    // list of tasks for the loop over the trees
    11221168   
    11231169    tlist2.AddToList(&rfgrow2);
    11241170    tlist2.AddToList(&rfwrite2);
     1171    tlist2.AddToList(&fillh2);
    11251172
    11261173    //-------------------
    … …  
    11351182
    11361183    tlist2.PrintStatistics(0, kTRUE);
     1184
     1185    plist2.FindObject("MHRanForestGini")->DrawClone();
     1186
    11371187
    11381188    // get adresses of objects which are used in the next eventloop
    … …  
    12121262
    12131263    //-----------------------------------------------------------------
    1214     // geometry is needed in  MHHillas... classes
    1215     MGeomCam *fGeom =
    1216              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    1217 
    12181264
    12191265
    … …  
    13041350    MProgressBar bar;
    13051351    MEvtLoop evtloop;
     1352    evtloop.SetName("UpdateRootFile");
    13061353    evtloop.SetParList(&pliston);
    13071354    evtloop.SetProgressBar(&bar);
    … …  
    13731420
    13741421
    1375 
    13761422  //---------------------------------------------------------------------
    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)
     1423  // Job B_SC_UP
     1424  //============
     1425
     1426    //  - create (or read in) optimum supercuts parameter values
     1427    //
     1428    //  - calculate the hadroness for the supercuts
     1429    //
     1430    //  - update input root file, including the hadroness
     1431
     1432
     1433 if (JobB_SC_UP)
    13861434 {
    13871435    gLog << "=====================================================" << endl;
    1388     gLog << "Macro CT1Analysis : Start of Job C" << endl;
     1436    gLog << "Macro CT1Analysis : Start of Job B_SC_UP" << endl;
    13891437
    13901438    gLog << "" << endl;
    1391     gLog << "Macro CT1Analysis : JobC = " << JobC  << endl;
    1392 
    1393 
    1394     // name of input data file
     1439    gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = "
     1440         << (JobB_SC_UP ? "kTRUE" : "kFALSE")  << ",  "
     1441         << (CMatrix    ? "kTRUE" : "kFALSE")  << ",  "
     1442         << (RMatrix    ? "kTRUE" : "kFALSE")  << ",  "
     1443         << (WOptimize  ? "kTRUE" : "kFALSE")  << ",  "
     1444         << (RTest      ? "kTRUE" : "kFALSE")  << ",  "
     1445         << (WSC        ? "kTRUE" : "kFALSE")  << endl;
     1446
     1447
     1448    //--------------------------------------------
     1449    // file which contains the initial parameter values for the supercuts
     1450    // if parSCinit ="" the initial values are taken from the constructor of
     1451    //                  MCT1Supercuts
     1452
     1453    TString parSCinit = outPath;
     1454    //parSCinit += "parSC_1709d";
     1455    parSCinit = "";
     1456
     1457    gLog << "parSCinit = " << parSCinit << endl;
     1458
     1459    //---------------
     1460    // file onto which the optimal parameter values for the supercuts
     1461    // are written
     1462
     1463    TString parSCfile = outPath;
     1464    parSCfile += "parSC_2310a";
     1465
     1466    gLog << "parSCfile = " << parSCfile << endl;
     1467
     1468    //--------------------------------------------
     1469    // file to be updated (either ON or MC)
     1470
     1471    //TString typeInput = "ON";
     1472    //TString typeInput = "OFF";
     1473    TString typeInput = "MC";
     1474    gLog << "typeInput = " << typeInput << endl;
     1475
     1476    // name of input root file
    13951477    TString filenameData = outPath;
    1396     filenameData += "OFF";
     1478    filenameData += typeInput;
    13971479    filenameData += "2.root";
    1398     gLog << "filenameData = " << filenameData << endl;
    1399 
    1400     // name of input MC file
    1401     TString filenameMC = outPath;
    1402     filenameMC += "MC";
    1403     filenameMC += "2.root";
    1404     gLog << "filenameMC   = " << filenameMC   << endl;
     1480    gLog << "filenameData = " << filenameData << endl;
     1481
     1482    // name of output root file
     1483    TString outNameImage = outPath;
     1484    outNameImage += typeInput;
     1485    outNameImage += "3.root";
     1486   
     1487
     1488    //TString outNameImage = filenameData;
     1489
     1490    gLog << "outNameImage = " << outNameImage << endl;
     1491
     1492    //--------------------------------------------
     1493    // files to be read for optimizing the supercuts
     1494    //
     1495    // for the training
     1496    TString filenameTrain = outPath;
     1497    filenameTrain += "ON";
     1498    filenameTrain += "1.root";
     1499    Int_t howManyTrain = 800000;
     1500    gLog << "filenameTrain = " << filenameTrain << ",   howManyTrain = "
     1501         << howManyTrain  << endl;
     1502
     1503    // for testing
     1504    TString filenameTest = outPath;
     1505    filenameTest += "ON";
     1506    filenameTest += "1.root";
     1507    Int_t howManyTest = 800000;
     1508
     1509    gLog << "filenameTest = " << filenameTest << ",   howManyTest = "
     1510         << howManyTest  << endl;
     1511   
     1512
     1513    //--------------------------------------------
     1514    // files to contain the matrices (generated from filenameTrain and
     1515    //                                               filenameTest)
     1516    //
     1517    // for the training
     1518    TString fileMatrixTrain = outPath;
     1519    fileMatrixTrain += "MatrixTrainSC";
     1520    fileMatrixTrain += ".root";
     1521    gLog << "fileMatrixTrain = " << fileMatrixTrain << endl;
     1522
     1523    // for testing
     1524    TString fileMatrixTest = outPath;
     1525    fileMatrixTest += "MatrixTestSC";
     1526    fileMatrixTest += ".root";
     1527    gLog << "fileMatrixTest = " << fileMatrixTest << endl;
     1528
     1529   
     1530
     1531    //---------------------------------------------------------------------
     1532    // Training and test matrices :
     1533    // - either create them and write them onto a file
     1534    // - or read them from a file
     1535
     1536
     1537    MCT1FindSupercuts findsuper;
     1538    findsuper.SetFilenameParam(parSCfile);
     1539    findsuper.SetHadronnessName("HadSC");
     1540
     1541
     1542    //--------------------------
     1543    // create matrices and write them onto files
     1544    if (CMatrix)
     1545    {
     1546      MH3 &mh3 = *(new MH3("MHillas.fSize"));
     1547      mh3.SetName("Target distribution for SIZE");
     1548
     1549      if (filenameTrain == filenameTest)
     1550      {
     1551        if ( !findsuper.DefineTrainTestMatrix(filenameTrain,
     1552                              howManyTrain, mh3, howManyTest,  mh3,
     1553                              fileMatrixTrain, fileMatrixTest)     )
     1554        {
     1555          *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl;
     1556          return;
     1557        }
     1558
     1559      }
     1560      else
     1561      {
     1562        if ( !findsuper.DefineTrainMatrix(filenameTrain,
     1563                              howManyTrain, mh3, fileMatrixTrain) )
     1564        {
     1565          *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl;
     1566          return;
     1567        }
     1568
     1569        if ( !findsuper.DefineTestMatrix( filenameTest, 
     1570                              howManyTest,  mh3, fileMatrixTest)  )
     1571        {
     1572          *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl;
     1573          return;
     1574        }
     1575      }
     1576     }
     1577
     1578    //--------------------------
     1579    // read matrices from files
     1580    //                             
     1581
     1582    if (RMatrix)
     1583      findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
     1584    //--------------------------
     1585
     1586
     1587
     1588    //---------------------------------------------------------------------
     1589    // optimize supercuts using the training sample
     1590    //
     1591    // the initial values are taken
     1592    //     - from the file parSCinit (if != "")
     1593    //     - or from the arrays params and steps (if their sizes are != 0)
     1594    //     - or from the MCT1Supercuts constructor
     1595
     1596
     1597if (WOptimize)
     1598  {
     1599    gLog << "CT1Analysis.C : optimize the supercuts using the training matrix"
     1600         << endl;
     1601
     1602    TArrayD params(0);
     1603    TArrayD steps(0);
     1604 
     1605    if (parSCinit == "")
     1606    {
     1607      Double_t vparams[104] = {
     1608      // LengthUp
     1609        0.315585,  0.001455, 0.203198, 0.005532, -0.001670, -0.020362,
     1610        0.007388, -0.013463,
     1611      // LengthLo
     1612        0.151530,  0.028323, 0.510707, 0.053089,  0.013708,  2.357993,
     1613        0.000080, -0.007157,
     1614      // WidthUp
     1615        0.145412, -0.001771, 0.054462, 0.022280, -0.009893,  0.056353,
     1616        0.020711, -0.016703,
     1617      // WidthLo
     1618        0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101,
     1619        0.006126, -0.002849,
     1620      // DistUp
     1621        1.787943,  0.0,      2.942310, 0.199815,  0.0,       0.249909,
     1622        0.189697,  0.0,
     1623      // DistLo
     1624        0.589406,  0.0,     -0.083964,-0.007975,  0.0,       0.045374,
     1625       -0.001750,  0.0,
     1626      // AsymUp
     1627        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1628        0.0,       0.0,
     1629      // AsymLo
     1630       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1631        0.0,       0.0,
     1632      // ConcUp
     1633        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1634        0.0,       0.0,
     1635      // ConcLo
     1636       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1637        0.0,       0.0,
     1638      // Leakage1Up
     1639        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1640        0.0,       0.0,
     1641      // Leakage1Lo
     1642       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
     1643        0.0,       0.0,
     1644      // AlphaUp
     1645        13.12344,  0.0,      0.0,      0.0,       0.0,       0.0,
     1646        0.0,       0.0                                                 };
     1647
     1648      Double_t vsteps[104] = {
     1649      // LengthUp
     1650        0.03,      0.0002,   0.02,     0.0006,    0.0002,    0.002,
     1651        0.0008,    0.002,
     1652      // LengthLo
     1653        0.02,      0.003,    0.05,     0.006,     0.002,     0.3,
     1654        0.0001,    0.0008,
     1655      // WidthUp
     1656        0.02,      0.0002,   0.006,    0.003,     0.002,     0.006,
     1657        0.002,     0.002,
     1658      // WidthLo
     1659        0.009,     0.0007,   0.008,    0.0004,    0.0005,    0.002,
     1660        0.0007,    0.003,
     1661      // DistUp
     1662        0.2,       0.0,      0.3,      0.02,      0.0,       0.03,
     1663        0.02,      0.0
     1664      // DistLo
     1665        0.06,      0.0,      0.009,    0.0008,    0.0,       0.005,
     1666        0.0002,    0.0
     1667      // AsymUp 
     1668        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1669        0.0,       0.0,
     1670      // AsymLo 
     1671        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1672        0.0,       0.0,
     1673      // ConcUp 
     1674        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1675        0.0,       0.0,
     1676      // ConcLo 
     1677        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1678        0.0,       0.0,
     1679      // Leakage1Up 
     1680        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1681        0.0,       0.0,
     1682      // Leakage1Lo 
     1683        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1684        0.0,       0.0,
     1685      // AlphaUp 
     1686        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
     1687        0.0,       0.0                                                 };
     1688
     1689      params.Set(104, vparams);
     1690      steps.Set (104, vsteps );
     1691    }
     1692
     1693    Bool_t rf;
     1694    rf = findsuper.FindParams(parSCinit, params, steps);
     1695
     1696    if (!rf)
     1697    {
     1698       gLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
     1699       return;
     1700    }
     1701  }
     1702
     1703    //--------------------------------------
     1704    // test the supercuts on the test sample
     1705    //   
     1706
     1707 if (RTest)
     1708 {
     1709    gLog << "CT1Analysis.C : test the supercuts on the test matrix" << endl;
     1710    Bool_t rt = findsuper.TestParams();
     1711    if (!rt)
     1712    {
     1713       gLog << "CT1Analysis.C : test of supercuts on the test matrix failed"
     1714            << endl;
     1715    }
     1716
     1717 }
    14051718
    14061719
    14071720    //-----------------------------------------------------------------
    1408 
     1721    // Update the input files with the SC hadronness
     1722    //
     1723
     1724 if (WSC)
     1725 {
     1726    gLog << "" << endl;
     1727    gLog << "========================================================" << endl;
     1728    gLog << "Update input file '" <<  filenameData
     1729         << "' with the SC hadronness" << endl;
     1730
     1731
     1732    //----------------------------------------------------
     1733    // read in optimum parameter values for the supercuts
     1734
     1735    TFile inparam(parSCfile);
     1736    MCT1Supercuts scin;
     1737    scin.Read("MCT1Supercuts");
     1738    inparam.Close();
     1739
     1740    gLog << "Parameter values for supercuts were read in from file '"
     1741         << parSCfile << "'" << endl;
     1742
     1743    TArrayD supercutsPar;
     1744    supercutsPar =  scin.GetParameters();
     1745
     1746    TArrayD supercutsStep;
     1747    supercutsStep =  scin.GetStepsizes();
     1748
     1749    gLog << "Parameter values for supercuts : " << endl;
     1750    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
     1751    {
     1752      gLog << supercutsPar[i] << ",  ";
     1753    }
     1754    gLog << endl;
     1755
     1756    gLog << "Step values for supercuts : " << endl;
     1757    for (Int_t i=0; i<supercutsStep.GetSize(); i++)
     1758    {
     1759      gLog << supercutsStep[i] << ",  ";
     1760    }
     1761    gLog << endl;
     1762
     1763
     1764    //----------------------------------------------------
    14091765    MTaskList tliston;
    14101766    MParList pliston;
     1767
     1768    // set the parameters of the supercuts
     1769    MCT1Supercuts supercuts;
     1770    supercuts.SetParameters(supercutsPar);
     1771    gLog << "parameter values for the supercuts used for updating the input file ' "
     1772         << filenameData << "'" << endl;
     1773    supercutsPar = supercuts.GetParameters();
     1774    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
     1775    {
     1776      gLog << supercutsPar[i] << ",  ";
     1777    }
     1778    gLog << endl;
    14111779
    14121780
    … …  
    14191787    //
    14201788
    1421     MReadMarsFile read("Events", filenameMC);
    1422     read.AddFile(filenameData);
     1789    MReadMarsFile read("Events", filenameData);
    14231790    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 
    14351791
    14361792    TString fHilName    = "MHillas";
    … …  
    14391795    TString fImgParName = "MNewImagePar";
    14401796
    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");
     1797
     1798    //.......................................................................
     1799    // calculation of hadroness for the supercuts
     1800    // (=0.25 if fullfilled, =0.75 otherwise)
     1801
     1802    TString hadSCName = "HadSC";
     1803    MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc);
     1804    sccalc.SetHadronnessName(hadSCName);
     1805
     1806
     1807    //.......................................................................
     1808
     1809
     1810      //MWriteRootFile write(outNameImage, "UPDATE");
     1811      //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
     1812
    14751813   
    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 
    1938     //.......................................................................
    1939 
    1940       MWriteRootFile write(filenameOFFup);
     1814      MWriteRootFile write(outNameImage, "RECREATE");
    19411815
    19421816      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    19511825      write.AddContainer("MNewImagePar",  "Events");
    19521826
    1953       //write.AddContainer("HadNN",         "Events");
    1954       write.AddContainer("HadSC",         "Events");
    19551827      write.AddContainer("HadRF",         "Events");
    1956 
    1957       write.AddContainer("MEnergyEst",    "Events");
     1828      write.AddContainer(hadSCName,       "Events");
     1829   
    19581830
    19591831    //-----------------------------------------------------------------
     1832    // geometry is needed in  MHHillas... classes
     1833    MGeomCam *fGeom =
     1834             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     1835
     1836    Float_t maxhadronness =  0.40;
     1837    Float_t maxalpha      =  20.0;
     1838    Float_t maxdist       =  10.0;
     1839
     1840    MFCT1SelFinal selfinalgh(fHilNameSrc);
     1841    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
     1842    selfinalgh.SetHadronnessName(hadSCName);
     1843    selfinalgh.SetName("SelFinalgh");
     1844    MContinue contfinalgh(&selfinalgh);
     1845    contfinalgh.SetName("ContSelFinalgh");
     1846
     1847    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
     1848    fillhadsc.SetName("HhadSC");
     1849
     1850    MFCT1SelFinal selfinal(fHilNameSrc);
     1851    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
     1852    selfinal.SetHadronnessName(hadSCName);
     1853    selfinal.SetName("SelFinal");
     1854    MContinue contfinal(&selfinal);
     1855    contfinal.SetName("ContSelFinal");
     1856
     1857    TString mh3name = "abs(Alpha)";
     1858    MBinning binsalphaabs("Binning"+mh3name);
     1859    binsalphaabs.SetEdges(50, -2.0, 98.0);
     1860
     1861    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
     1862    alphaabs.SetName(mh3name);
     1863
     1864    TH1  &alphahist = alphaabs->GetHist();
     1865
     1866    MFillH alpha(&alphaabs);
     1867    alpha.SetName("FillAlphaAbs");
     1868
     1869    MFillH hfill1("MHHillas",    fHilName);
     1870    hfill1.SetName("HHillas");
     1871
     1872    MFillH hfill2("MHStarMap",   fHilName);
     1873    hfill2.SetName("HStarMap");
     1874
     1875    MFillH hfill3("MHHillasExt",    fHilNameSrc);
     1876    hfill3.SetName("HHillasExt");
     1877   
     1878    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
     1879    hfill4.SetName("HHillasSrc");   
     1880
     1881    MFillH hfill5("MHNewImagePar", fImgParName);
     1882    hfill5.SetName("HNewImagePar");
     1883
     1884    //*****************************
     1885    // entries in MParList
     1886
     1887    pliston.AddToList(&tliston);
     1888    InitBinnings(&pliston);
     1889
     1890    pliston.AddToList(&supercuts);
     1891
     1892    pliston.AddToList(&binsalphaabs);
     1893    pliston.AddToList(&alphaabs);
     1894
     1895    //*****************************
     1896    // entries in MTaskList
     1897   
     1898    tliston.AddToList(&read);
     1899
     1900    tliston.AddToList(&sccalc);
     1901    tliston.AddToList(&fillhadsc);
     1902
     1903    tliston.AddToList(&write);
     1904    tliston.AddToList(&contfinalgh);
     1905
     1906    tliston.AddToList(&alpha);
     1907    tliston.AddToList(&hfill1);
     1908    tliston.AddToList(&hfill2);
     1909    tliston.AddToList(&hfill3);
     1910    tliston.AddToList(&hfill4);
     1911    tliston.AddToList(&hfill5);
     1912
     1913    tliston.AddToList(&contfinal);
     1914
     1915    //*****************************
     1916
     1917    //-------------------------------------------
     1918    // Execute event loop
     1919    //
     1920    MProgressBar bar;
     1921    MEvtLoop evtloop;
     1922    evtloop.SetParList(&pliston);
     1923    evtloop.SetProgressBar(&bar);
     1924
     1925    Int_t maxevents = -1;
     1926    if ( !evtloop.Eventloop(maxevents) )
     1927        return;
     1928
     1929    tliston.PrintStatistics(0, kTRUE);
     1930
     1931
     1932    //-------------------------------------------
     1933    // Display the histograms
     1934    //
     1935    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
     1936
     1937    pliston.FindObject("MHHillas")->DrawClone();
     1938    pliston.FindObject("MHHillasExt")->DrawClone();
     1939    pliston.FindObject("MHHillasSrc")->DrawClone();
     1940    pliston.FindObject("MHNewImagePar")->DrawClone();
     1941    pliston.FindObject("MHStarMap")->DrawClone();
     1942
     1943     //-------------------------------------------
     1944    // fit alpha distribution to get the number of excess events and
     1945    // calculate significance of gamma signal in the alpha plot
     1946 
     1947    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
     1948    TH1  &alphaHist = absalpha->GetHist();
     1949    alphaHist.SetXTitle("|alpha|  [\\circ]");
     1950    alphaHist.SetName("alpha-macro");
     1951
     1952    Double_t alphasig = 13.1;
     1953    Double_t alphamin = 30.0;
     1954    Double_t alphamax = 90.0;
     1955    Int_t    degree   =    2;
     1956    Double_t significance = -99.0;
     1957    Bool_t   drawpoly  = kTRUE;
     1958    Bool_t   fitgauss  = kTRUE;
     1959    Bool_t   print     = kTRUE;
     1960
     1961    MHFindSignificance findsig;
     1962    findsig.SetRebin(kTRUE);
     1963    findsig.SetReduceDegree(kFALSE);
     1964
     1965    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
     1966                        alphasig, drawpoly, fitgauss, print);
     1967    significance = findsig.GetSignificance();
     1968    Float_t alphasi = findsig.GetAlphasi();
     1969
     1970    gLog << "For file '" << filenameData << "' : " << endl;
     1971    gLog << "Significance of gamma signal after supercuts : "
     1972         << significance << " (for |alpha| < " << alphasi << " degrees)"
     1973         << endl;
     1974
     1975    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
     1976
     1977    //-------------------------------------------
     1978
     1979    DeleteBinnings(&pliston);
     1980 }
     1981
     1982
     1983    gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl;
     1984    gLog << "=======================================================" << endl;
     1985 }
     1986  //---------------------------------------------------------------------
     1987
     1988
     1989
     1990  //---------------------------------------------------------------------
     1991  // Job C 
     1992  //======
     1993
     1994    //  - read ON1 and MC1 data files 
     1995    //    which should have been updated to contain the hadronnesses
     1996    //    for the method of Random Forest and for the SUPERCUTS
     1997    //  - produce Neyman-Pearson plots
     1998 
     1999 if (JobC)
     2000 {
     2001    gLog << "=====================================================" << endl;
     2002    gLog << "Macro CT1Analysis : Start of Job C" << endl;
     2003
     2004    gLog << "" << endl;
     2005    gLog << "Macro CT1Analysis : JobC = "
     2006         << (JobC ? "kTRUE" : "kFALSE")  << endl;
     2007
     2008
     2009    // name of input data file
     2010    TString filenameData = outPath;
     2011    filenameData += "OFF";
     2012    filenameData += "3.root";
     2013    gLog << "filenameData = " << filenameData << endl;
     2014
     2015    // name of input MC file
     2016    TString filenameMC = outPath;
     2017    filenameMC += "MC";
     2018    filenameMC += "3.root";
     2019    gLog << "filenameMC   = " << filenameMC   << endl;
     2020
     2021
     2022    //-----------------------------------------------------------------
     2023
     2024    MTaskList tliston;
     2025    MParList pliston;
     2026
     2027
     2028    // geometry is needed in  MHHillas... classes
     2029    MGeomCam *fGeom =
     2030             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     2031
     2032    //-------------------------------------------
     2033    // create the tasks which should be executed
     2034    //
     2035
     2036    MReadMarsFile read("Events", filenameMC);
     2037    read.AddFile(filenameData);
     2038    read.DisableAutoScheme();
     2039
     2040
     2041    //.......................................................................
     2042    // names of hadronness containers
     2043
     2044    TString hadSCName = "HadSC";
     2045    TString hadRFName = "HadRF";
     2046
     2047    //.......................................................................
     2048
     2049
     2050    TString fHilName    = "MHillas";
     2051    TString fHilNameExt = "MHillasExt";
     2052    TString fHilNameSrc = "MHillasSrc";
     2053    TString fImgParName = "MNewImagePar";
     2054
     2055    Float_t maxhadronness =  0.40;
     2056    Float_t maxalpha      =  20.0;
     2057    Float_t maxdist       =  10.0;
     2058
     2059    MFCT1SelFinal selfinalgh(fHilNameSrc);
     2060    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
     2061    selfinalgh.SetHadronnessName(hadRFName);
     2062    selfinalgh.SetName("SelFinalgh");
     2063    MContinue contfinalgh(&selfinalgh);
     2064    contfinalgh.SetName("ContSelFinalgh");
     2065
     2066
     2067    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
     2068    fillhadsc.SetName("HhadSC");
     2069    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
     2070    fillhadrf.SetName("HhadRF");
     2071
     2072    MFCT1SelFinal selfinal(fHilNameSrc);
     2073    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
     2074    selfinal.SetHadronnessName(hadRFName);
     2075    selfinal.SetName("SelFinal");
     2076    MContinue contfinal(&selfinal);
     2077    contfinal.SetName("ContSelFinal");
     2078
     2079
     2080    MFillH hfill1("MHHillas",    fHilName);
     2081    hfill1.SetName("HHillas");
     2082
     2083    MFillH hfill2("MHStarMap",   fHilName);
     2084    hfill2.SetName("HStarMap");
     2085
     2086    MFillH hfill3("MHHillasExt",    fHilNameSrc);
     2087    hfill3.SetName("HHillasExt");
     2088   
     2089    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
     2090    hfill4.SetName("HHillasSrc");   
     2091
     2092    MFillH hfill5("MHNewImagePar", fImgParName);
     2093    hfill5.SetName("HNewImagePar");
     2094
     2095
     2096    //*****************************
     2097    // entries in MParList
     2098
     2099    pliston.AddToList(&tliston);
     2100    InitBinnings(&pliston);
     2101
     2102
     2103    //*****************************
     2104    // entries in MTaskList
     2105   
     2106    tliston.AddToList(&read);
     2107
     2108    tliston.AddToList(&fillhadsc);
     2109    tliston.AddToList(&fillhadrf);
     2110   
     2111    tliston.AddToList(&contfinalgh);
     2112    tliston.AddToList(&hfill1);
     2113    tliston.AddToList(&hfill2);
     2114    tliston.AddToList(&hfill3);
     2115    tliston.AddToList(&hfill4);
     2116    tliston.AddToList(&hfill5);
     2117
     2118    tliston.AddToList(&contfinal);
     2119
     2120    //*****************************
     2121
     2122    //-------------------------------------------
     2123    // Execute event loop
     2124    //
     2125    MProgressBar bar;
     2126    MEvtLoop evtloop;
     2127    evtloop.SetParList(&pliston);
     2128    evtloop.SetProgressBar(&bar);
     2129
     2130    Int_t maxevents = -1;
     2131    //Int_t maxevents = 35000;
     2132    if ( !evtloop.Eventloop(maxevents) )
     2133        return;
     2134
     2135    tliston.PrintStatistics(0, kTRUE);
     2136
     2137
     2138    //-------------------------------------------
     2139    // Display the histograms
     2140    //
     2141
     2142    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
     2143    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
     2144
     2145    pliston.FindObject("MHHillas")->DrawClone();
     2146    pliston.FindObject("MHHillasExt")->DrawClone();
     2147    pliston.FindObject("MHHillasSrc")->DrawClone();
     2148    pliston.FindObject("MHNewImagePar")->DrawClone();
     2149    pliston.FindObject("MHStarMap")->DrawClone();
     2150
     2151    DeleteBinnings(&pliston);
     2152
     2153    gLog << "Macro CT1Analysis : End of Job C" << endl;
     2154    gLog << "===================================================" << endl;
     2155 }
     2156
     2157
     2158
     2159  //---------------------------------------------------------------------
     2160  // Job D
     2161  //======
     2162
     2163    //  - select g/h separation method XX
     2164    //  - read ON2 (or MC2) root file
     2165    //  - apply cuts in hadronness
     2166    //  - make plots
     2167
     2168
     2169 if (JobD)
     2170 {
     2171    gLog << "=====================================================" << endl;
     2172    gLog << "Macro CT1Analysis : Start of Job D" << endl;
     2173
     2174    gLog << "" << endl;
     2175    gLog << "Macro CT1Analysis : JobD = "
     2176         << (JobD        ? "kTRUE" : "kFALSE")  << endl;
     2177
     2178
     2179    // type of data to be analysed
     2180    TString typeData = "ON";
     2181    //TString typeData = "OFF";
     2182    //TString typeData = "MC";
     2183    gLog << "typeData = " << typeData << endl;
     2184
     2185    TString ext      = "3.root";
     2186
     2187
     2188    //------------------------------
     2189    // selection of g/h separation method
     2190    // and definition of final selections
     2191
     2192    //TString XX("SC");
     2193    TString XX("RF");
     2194    TString fhadronnessName("Had");
     2195    fhadronnessName += XX;
     2196    gLog << "fhadronnessName = " << fhadronnessName << endl;
     2197
     2198    // maximum values of the hadronness, |ALPHA| and DIST
     2199    Float_t maxhadronness   = 0.233;
     2200    Float_t maxalpha        = 20.0;
     2201    Float_t maxdist         = 10.0;
     2202    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
     2203         << maxhadronness << ",  " << maxalpha << ",  "
     2204         << maxdist << endl;
     2205
     2206
     2207    //------------------------------
     2208    // name of data file to be analysed
     2209    TString filenameData(outPath);
     2210    filenameData += typeData;
     2211    filenameData += ext;
     2212    gLog << "filenameData = " << filenameData << endl;
     2213
     2214
     2215
     2216    //*************************************************************************
     2217    //
     2218    // Analyse the data
     2219    //
     2220
     2221    MTaskList tliston;
     2222    MParList pliston;
     2223
     2224    // geometry is needed in  MHHillas... classes
     2225    MGeomCam *fGeom =
     2226             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     2227
     2228
     2229    TString fHilName    = "MHillas";
     2230    TString fHilNameExt = "MHillasExt";
     2231    TString fHilNameSrc = "MHillasSrc";
     2232    TString fImgParName = "MNewImagePar";
     2233
     2234    //-------------------------------------------
     2235    // create the tasks which should be executed
     2236    //
     2237
     2238    MReadMarsFile read("Events", filenameData);
     2239    read.DisableAutoScheme();
     2240
     2241
     2242    //-----------------------------------------------------------------
     2243    // geometry is needed in  MHHillas... classes
     2244    MGeomCam *fGeom =
     2245             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     2246
     2247    MFCT1SelFinal selfinalgh(fHilNameSrc);
     2248    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
     2249    selfinalgh.SetHadronnessName(fhadronnessName);
     2250    selfinalgh.SetName("SelFinalgh");
     2251    MContinue contfinalgh(&selfinalgh);
     2252    contfinalgh.SetName("ContSelFinalgh");
     2253
     2254    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
     2255    fillhadsc.SetName("HhadSC");
     2256    MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
     2257    fillhadrf.SetName("HhadRF");
     2258
     2259    TString mh3name = "abs(Alpha)";
     2260    MBinning binsalphaabs("Binning"+mh3name);
     2261    binsalphaabs.SetEdges(50, -2.0, 98.0);
     2262
     2263    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
     2264    alphaabs.SetName(mh3name);
     2265
     2266    TH1  &alphahist = alphaabs->GetHist();
     2267
     2268    MFillH alpha(&alphaabs);
     2269    alpha.SetName("FillAlphaAbs");
     2270
     2271    MFillH hfill1("MHHillas",    fHilName);
     2272    hfill1.SetName("HHillas");
     2273
     2274    MFillH hfill2("MHStarMap",   fHilName);
     2275    hfill2.SetName("HStarMap");
     2276
     2277    MFillH hfill3("MHHillasExt",   fHilNameSrc);
     2278    hfill3.SetName("HHillasExt");   
     2279
     2280    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
     2281    hfill4.SetName("HHillasSrc");   
     2282
     2283    MFillH hfill5("MHNewImagePar", fImgParName);
     2284    hfill5.SetName("HNewImagePar");   
    19602285
    19612286    MFCT1SelFinal selfinal(fHilNameSrc);
    19622287    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    19632288    selfinal.SetHadronnessName(fhadronnessName);
     2289    selfinal.SetName("SelFinal");
    19642290    MContinue contfinal(&selfinal);
     2291    contfinal.SetName("ContSelFinal");
    19652292
    19662293
    … …  
    19682295    // entries in MParList
    19692296
    1970     plistoff.AddToList(&tlistoff);
    1971     InitBinnings(&plistoff);
    1972 
     2297    pliston.AddToList(&tliston);
     2298    InitBinnings(&pliston);
     2299    pliston.AddToList(&binsalphaabs);
     2300    pliston.AddToList(&alphaabs);
    19732301
    19742302    //*****************************
    19752303    // entries in MTaskList
    19762304   
    1977     tlistoff.AddToList(&read);
    1978     tlistoff.AddToList(&eest2);
    1979     tlistoff.AddToList(&write);
    1980     tlistoff.AddToList(&contfinal);
     2305    tliston.AddToList(&read);
     2306
     2307    tliston.AddToList(&contfinalgh);
     2308
     2309    tliston.AddToList(&fillhadsc);
     2310    tliston.AddToList(&fillhadrf);
     2311
     2312    tliston.AddToList(&alpha);
     2313    tliston.AddToList(&hfill1);
     2314    tliston.AddToList(&hfill2);
     2315    tliston.AddToList(&hfill3);
     2316    tliston.AddToList(&hfill4);
     2317    tliston.AddToList(&hfill5);
     2318
     2319    tliston.AddToList(&contfinal);
    19812320
    19822321    //*****************************
    … …  
    19872326    MProgressBar bar;
    19882327    MEvtLoop evtloop;
    1989     evtloop.SetParList(&plistoff);
     2328    evtloop.SetParList(&pliston);
    19902329    evtloop.SetProgressBar(&bar);
    19912330
    19922331    Int_t maxevents = -1;
    1993     //Int_t maxevents = 1000;
     2332    //Int_t maxevents = 10000;
    19942333    if ( !evtloop.Eventloop(maxevents) )
    19952334        return;
    19962335
    1997     tlistoff.PrintStatistics(0, kTRUE);
    1998     DeleteBinnings(&plistoff);
    1999 
    2000     //---------------------------------------------------
    2001 
    2002     //---------------------------------------------------
    2003     // Update ON data
     2336    tliston.PrintStatistics(0, kTRUE);
     2337
     2338
     2339    //-------------------------------------------
     2340    // Display the histograms
     2341    //
     2342
     2343    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
     2344    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
     2345
     2346    pliston.FindObject("MHHillas")->DrawClone();
     2347    pliston.FindObject("MHHillasExt")->DrawClone();
     2348    pliston.FindObject("MHHillasSrc")->DrawClone();
     2349    pliston.FindObject("MHNewImagePar")->DrawClone();
     2350    pliston.FindObject("MHStarMap")->DrawClone();
     2351
     2352
     2353    //-------------------------------------------
     2354
     2355    // fit alpha distribution to get the number of excess events and
     2356    // calculate significance of gamma signal in the alpha plot
     2357 
     2358    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
     2359    TH1  &alphaHist = absalpha->GetHist();
     2360    alphaHist.SetXTitle("|alpha|  [\\circ]");
     2361    alphaHist.SetName("alpha-JobD");
     2362
     2363    Double_t alphasig = 13.1;
     2364    Double_t alphamin = 30.0;
     2365    Double_t alphamax = 90.0;
     2366    Int_t    degree   =    2;
     2367    Double_t significance = -99.0;
     2368    Bool_t   drawpoly  = kTRUE;
     2369    Bool_t   fitgauss  = kTRUE;
     2370    Bool_t   print     = kTRUE;
     2371
     2372    MHFindSignificance findsig;
     2373    findsig.SetRebin(kTRUE);
     2374    findsig.SetReduceDegree(kFALSE);
     2375
     2376    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
     2377                        alphasig, drawpoly, fitgauss, print);
     2378    significance = findsig.GetSignificance();
     2379    Float_t alphasi = findsig.GetAlphasi();
     2380
     2381    gLog << "For file '" << filenameData << "' : " << endl;
     2382    gLog << "Significance of gamma signal after supercuts : "
     2383         << significance << " (for |alpha| < " << alphasi << " degrees)"
     2384         << endl;
     2385
     2386    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
     2387
     2388    //-------------------------------------------
     2389
     2390
     2391    DeleteBinnings(&pliston);
     2392
     2393    gLog << "Macro CT1Analysis : End of Job D" << endl;
     2394    gLog << "=======================================================" << endl;
     2395 }
     2396  //---------------------------------------------------------------------
     2397
     2398
     2399
     2400
     2401
     2402  //---------------------------------------------------------------------
     2403  // Job E_XX
     2404  //=========
     2405
     2406    //  - select g/h separation method XX
     2407    //  - read MC_XX2.root file
     2408    //  - calculate eff. collection area
     2409    //  - read ON_XX2.root file
     2410    //  - apply final cuts
     2411    //  - calculate flux
     2412    //  - write root file for ON data after final cuts (ON_XX3.root))
     2413
     2414
     2415 if (JobE_XX)
     2416 {
     2417    gLog << "=====================================================" << endl;
     2418    gLog << "Macro CT1Analysis : Start of Job E_XX" << endl;
     2419
     2420    gLog << "" << endl;
     2421    gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = "
     2422         << (JobE_XX ? "kTRUE" : "kFALSE")  << ",  "
     2423         << (OEEst ?   "kTRUE" : "kFALSE")  << ",  "
     2424         << (WEX     ? "kTRUE" : "kFALSE")  << endl;
     2425
     2426
     2427    // type of data to be analysed
     2428    //TString typeData = "ON";
     2429    //TString typeData = "OFF";
     2430    TString typeData = "MC";
     2431    gLog << "typeData = " << typeData << endl;
     2432
     2433    TString typeMC   = "MC";
     2434    TString ext      = "3.root";
     2435    TString extout   = "4.root";
     2436
     2437    //------------------------------
     2438    // selection of g/h separation method
     2439    // and definition of final selections
     2440
     2441    //TString XX("SC");
     2442    TString XX("RF");
     2443    TString fhadronnessName("Had");
     2444    fhadronnessName += XX;
     2445    gLog << "fhadronnessName = " << fhadronnessName << endl;
     2446
     2447    // maximum values of the hadronness, |ALPHA| and DIST
     2448    Float_t maxhadronness   = 0.23;
     2449    Float_t maxalpha        = 20.0;
     2450    Float_t maxdist         = 10.0;
     2451    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
     2452         << maxhadronness << ",  " << maxalpha << ",  "
     2453         << maxdist << endl;
     2454
     2455    //------------------------------
     2456    // name of MC file to be used for optimizing the energy estimator
     2457    TString filenameOpt(outPath);
     2458    filenameOpt += typeMC;
     2459    filenameOpt += ext;
     2460    gLog << "filenameOpt = " << filenameOpt << endl;
     2461
     2462    //------------------------------
     2463    // name of file containing the parameters of the energy estimator
     2464    TString energyParName(outPath);
     2465    energyParName += "energyest_";
     2466    energyParName += XX;
     2467    energyParName += ".root";
     2468    gLog << "energyParName = " << energyParName << endl;
     2469
     2470    //------------------------------
     2471    // name of MC file to be used for calculating the eff. collection areas
     2472    TString filenameArea(outPath);
     2473    filenameArea += typeMC;
     2474    filenameArea += ext;
     2475    gLog << "filenameArea = " << filenameArea << endl;
     2476
     2477    //------------------------------
     2478    // name of file containing the eff. collection areas
     2479    TString collareaName(outPath);
     2480    collareaName += "area_";
     2481    collareaName += XX;
     2482    collareaName += ".root";
     2483    gLog << "collareaName = " << collareaName << endl;
     2484
     2485    //------------------------------
     2486    // name of data file to be analysed
     2487    TString filenameData(outPath);
     2488    filenameData += typeData;
     2489    filenameData += ext;
     2490    gLog << "filenameData = " << filenameData << endl;
     2491
     2492    //------------------------------
     2493    // name of output data file (after the final cuts)
     2494    TString filenameDataout(outPath);
     2495    filenameDataout += typeData;
     2496    filenameDataout += "_";
     2497    filenameDataout += XX;
     2498    filenameDataout += extout;
     2499    gLog << "filenameDataout = " << filenameDataout << endl;
     2500
     2501    //------------------------------
     2502    // name of file containing histograms for flux calculastion
     2503    TString filenameResults(outPath);
     2504    filenameResults += typeData;
     2505    filenameResults += "Results_";
     2506    filenameResults += XX;
     2507    filenameResults += extout;
     2508    gLog << "filenameResults = " << filenameResults << endl;
     2509
     2510
     2511    //====================================================================
     2512    gLog << "-----------------------------------------------" << endl;
     2513    gLog << "Start calculation of effective collection areas" << endl;
     2514    MParList  parlist;
     2515    MTaskList tasklist;
     2516
     2517    //---------------------------------------
     2518    // Setup the tasks to be executed
     2519    //
     2520    MReadMarsFile reader("Events", filenameArea);
     2521    reader.DisableAutoScheme();
     2522
     2523    MFCT1SelFinal cuthadrons;
     2524    cuthadrons.SetHadronnessName(fhadronnessName);
     2525    cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
     2526
     2527    MContinue conthadrons(&cuthadrons);
     2528
     2529    MHMcCT1CollectionArea collarea;
     2530    collarea.SetEaxis(MHMcCT1CollectionArea::kLinear);
     2531
     2532    MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
     2533    filler.SetName("CollectionArea");
     2534
     2535    //********************************
     2536    // entries in MParList
     2537
     2538    parlist.AddToList(&tasklist);
     2539    InitBinnings(&parlist);
     2540    parlist.AddToList(&collarea);
     2541
     2542    //********************************
     2543    // entries in MTaskList
     2544
     2545    tasklist.AddToList(&reader);   
     2546    tasklist.AddToList(&conthadrons);
     2547    tasklist.AddToList(&filler);
     2548
     2549    //********************************
     2550
     2551    //-----------------------------------------
     2552    // Execute event loop
     2553    //
     2554    MEvtLoop magic;
     2555    magic.SetParList(&parlist);
     2556
     2557    MProgressBar bar;
     2558    magic.SetProgressBar(&bar);
     2559    if (!magic.Eventloop())
     2560        return;
     2561
     2562    tasklist.PrintStatistics(0, kTRUE);
     2563
     2564    // Calculate effective collection areas
     2565    // and display the histograms
     2566    //
     2567    //MHMcCT1CollectionArea *collarea =
     2568    //     (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
     2569    collarea.CalcEfficiency();
     2570    collarea.DrawClone();
     2571
     2572    // save binnings for call to CT1EEst
     2573    MBinning *binsE     = (MBinning*)parlist.FindObject("BinningE");
     2574    if (!binsE)
     2575        {
     2576          gLog << "Object 'BinningE' not found in MParList" << endl;
     2577          return;
     2578        }
     2579    MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
     2580    if (!binsTheta)
     2581        {
     2582          gLog << "Object 'BinningTheta' not found in MParList" << endl;
     2583          return;
     2584        }
     2585
     2586
     2587    //---------------------------------------------
     2588    // Write histograms to a file
     2589    //
     2590
     2591    TFile f(collareaName, "RECREATE");
     2592    collarea.GetHist()->Write();
     2593    collarea.GetHAll()->Write();
     2594    collarea.GetHSel()->Write();
     2595    f.Close();
     2596
     2597    gLog << "Collection area plots written onto file " << collareaName << endl;
     2598
     2599    gLog << "Calculation of effective collection areas done" << endl;
     2600    gLog << "-----------------------------------------------" << endl;   
     2601    //------------------------------------------------------------------
     2602
     2603
     2604    TString fHilName    = "MHillas";
     2605    TString fHilNameExt = "MHillasExt";
     2606    TString fHilNameSrc = "MHillasSrc";
     2607    TString fImgParName = "MNewImagePar";
     2608
     2609
     2610 if (OEEst)
     2611 {
     2612   //===========================================================
     2613    //
     2614    // Optimization of energy estimator
     2615    //
     2616    gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl;
     2617
     2618    TString inpath("");
     2619    TString outpath("");
     2620    Int_t howMany = 2000;
     2621    CT1EEst(inpath,   filenameOpt,   outpath, energyParName,
     2622            fHilName, fHilNameSrc,   fhadronnessName,
     2623            howMany,  maxhadronness, maxalpha, maxdist,
     2624            binsE, binsTheta);
     2625    gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
     2626 }
     2627
     2628 if (WEX)
     2629 {
     2630    //-----------------------------------------------------------
     2631    //
     2632    // Read in parameters of energy estimator ("MMcEnergyEst")
     2633    //                   and migration matrix ("MHMcEnergyMigration")
     2634    //
     2635    gLog << "================================================================"
     2636         << endl;
     2637    gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
     2638         << energyParName << "'" << endl;
     2639    TFile enparam(energyParName);
     2640    enparam.ls();
     2641    MMcEnergyEst mcest("MMcEnergyEst");
     2642    mcest.Read("MMcEnergyEst");
     2643
     2644    //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
     2645    gLog << "Parameters of energy estimator were read in" << endl;
     2646
     2647
     2648    gLog << "Read in Migration matrix" << endl;   
     2649
     2650    MHMcEnergyMigration mighiston("MHMcEnergyMigration");
     2651    mighiston.Read("MHMcEnergyMigration");
     2652    //MHMcEnergyMigration &mighiston =
     2653    //      *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
     2654
     2655    gLog << "Migration matrix was read in" << endl;
     2656
     2657
     2658    TArrayD parA(mcest.GetNumCoeffA());
     2659    TArrayD parB(mcest.GetNumCoeffB());
     2660    for (Int_t i=0; i<parA.GetSize(); i++)
     2661      parA[i] = mcest.GetCoeff(i);
     2662    for (Int_t i=0; i<parB.GetSize(); i++)
     2663      parB[i] = mcest.GetCoeff( i+parA.GetSize() );
     2664
     2665    //*************************************************************************
     2666    //
     2667    // Analyse the data
    20042668    //
    20052669    gLog << "============================================================"
    20062670         << endl;
    2007     gLog << "Macro CT1Analysis.C : update file '" << filenameON
    2008          << "'" << endl;
     2671    gLog << "Analyse the data" << endl;
    20092672
    20102673    MTaskList tliston;
    20112674    MParList pliston;
    2012 
    20132675
    20142676    // geometry is needed in  MHHillas... classes
    … …  
    20162678             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    20172679
     2680
    20182681    //-------------------------------------------
    20192682    // create the tasks which should be executed
    20202683    //
    20212684
    2022     MReadMarsFile read("Events", filenameON);
     2685    MReadMarsFile read("Events", filenameData);
    20232686    read.DisableAutoScheme();
    20242687
    2025     //---------------------------
    2026     // calculate estimated energy
    2027 
    2028     MEnergyEstParam eest2(fHilName);
    2029     eest2.Add(fHilNameSrc);
    2030 
    2031     eest2.SetCoeffA(parA);
    2032     eest2.SetCoeffB(parB);
    2033 
    20342688    //.......................................................................
    20352689
    2036       MWriteRootFile write(filenameONup);
     2690      gLog << "CT1Analysis.C : write root file '" << filenameDataout
     2691           << "'" << endl;
     2692   
     2693      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
     2694
     2695
     2696      MWriteRootFile write(filenameDataout, "RECREATE");
    20372697
    20382698      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    20532713      write.AddContainer("MEnergyEst",    "Events");
    20542714
    2055     //-----------------------------------------------------------------
    2056 
    2057     MFCT1SelFinal selfinal(fHilNameSrc);
    2058     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    2059     selfinal.SetHadronnessName(fhadronnessName);
    2060     MContinue contfinal(&selfinal);
    2061 
    2062 
    2063     //*****************************
    2064     // entries in MParList
    2065 
    2066     pliston.AddToList(&tliston);
    2067     InitBinnings(&pliston);
    2068 
    2069 
    2070     //*****************************
    2071     // entries in MTaskList
    2072    
    2073     tliston.AddToList(&read);
    2074     tliston.AddToList(&eest2);
    2075     tliston.AddToList(&write);
    2076     tliston.AddToList(&contfinal);
    2077 
    2078     //*****************************
    2079 
    2080     //-------------------------------------------
    2081     // Execute event loop
    2082     //
    2083     MProgressBar bar;
    2084     MEvtLoop evtloop;
    2085     evtloop.SetParList(&pliston);
    2086     evtloop.SetProgressBar(&bar);
    2087 
    2088     Int_t maxevents = -1;
    2089     //Int_t maxevents = 1000;
    2090     if ( !evtloop.Eventloop(maxevents) )
    2091         return;
    2092 
    2093     tliston.PrintStatistics(0, kTRUE);
    2094     DeleteBinnings(&pliston);
    2095 
    2096     //---------------------------------------------------
    2097 
    2098     //---------------------------------------------------
    2099     // Update MC data
    2100     //
    2101     gLog << "============================================================"
    2102          << endl;
    2103     gLog << "Macro CT1Analysis.C : update file '" << filenameMC
    2104          << "'" << endl;
    2105 
    2106     MTaskList tlistmc;
    2107     MParList plistmc;
    2108 
    2109     //-------------------------------------------
    2110     // create the tasks which should be executed
    2111     //
    2112 
    2113     MReadMarsFile read("Events", filenameMC);
    2114     read.DisableAutoScheme();
    2115 
    2116     //---------------------------
    2117     // calculate estimated energy
    2118 
    2119     MEnergyEstParam eest2(fHilName);
    2120     eest2.Add(fHilNameSrc);
    2121 
    2122     eest2.SetCoeffA(parA);
    2123     eest2.SetCoeffB(parB);
    2124 
    2125     //.......................................................................
    2126 
    2127       MWriteRootFile write(filenameMCup);
    2128 
    2129       write.AddContainer("MRawRunHeader", "RunHeaders");
    2130       write.AddContainer("MTime",         "Events");
    2131       write.AddContainer("MMcEvt",        "Events");
    2132       write.AddContainer("ThetaOrig",     "Events");
    2133       write.AddContainer("MSrcPosCam",    "Events");
    2134       write.AddContainer("MSigmabar",     "Events");
    2135       write.AddContainer("MHillas",       "Events");
    2136       write.AddContainer("MHillasExt",    "Events");
    2137       write.AddContainer("MHillasSrc",    "Events");
    2138       write.AddContainer("MNewImagePar",  "Events");
    2139 
    2140       //write.AddContainer("HadNN",         "Events");
    2141       write.AddContainer("HadSC",         "Events");
    2142       write.AddContainer("HadRF",         "Events");
    2143 
    2144       write.AddContainer("MEnergyEst",    "Events");
    2145 
    2146     //-----------------------------------------------------------------
    2147 
    2148     MFCT1SelFinal selfinal(fHilNameSrc);
    2149     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    2150     selfinal.SetHadronnessName(fhadronnessName);
    2151     MContinue contfinal(&selfinal);
    2152 
    2153 
    2154     //*****************************
    2155     // entries in MParList
    2156 
    2157     plistmc.AddToList(&tlistmc);
    2158     InitBinnings(&plistmc);
    2159 
    2160 
    2161     //*****************************
    2162     // entries in MTaskList
    2163    
    2164     tlistmc.AddToList(&read);
    2165     tlistmc.AddToList(&eest2);
    2166     tlistmc.AddToList(&write);
    2167     tlistmc.AddToList(&contfinal);
    2168 
    2169     //*****************************
    2170 
    2171     //-------------------------------------------
    2172     // Execute event loop
    2173     //
    2174     MProgressBar bar;
    2175     MEvtLoop evtloop;
    2176     evtloop.SetParList(&plistmc);
    2177     evtloop.SetProgressBar(&bar);
    2178 
    2179     Int_t maxevents = -1;
    2180     //Int_t maxevents = 1000;
    2181     if ( !evtloop.Eventloop(maxevents) )
    2182         return;
    2183 
    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";
    2383 
    2384     //-------------------------------------------
    2385     // create the tasks which should be executed
    2386     //
    2387 
    2388     MReadMarsFile read("Events", filenameData);
    2389     read.DisableAutoScheme();
    2390 
    2391     //.......................................................................
    2392 
    2393 
    2394       MWriteRootFile write(filenameDataout);
    2395 
    2396       write.AddContainer("MRawRunHeader", "RunHeaders");
    2397       write.AddContainer("MTime",         "Events");
    2398       write.AddContainer("MMcEvt",        "Events");
    2399       write.AddContainer("ThetaOrig",     "Events");
    2400       write.AddContainer("MSrcPosCam",    "Events");
    2401       write.AddContainer("MSigmabar",     "Events");
    2402       write.AddContainer("MHillas",       "Events");
    2403       write.AddContainer("MHillasExt",    "Events");
    2404       write.AddContainer("MHillasSrc",    "Events");
    2405       write.AddContainer("MNewImagePar",  "Events");
    2406 
    2407       //write.AddContainer("HadNN",         "Events");
    2408       write.AddContainer("HadSC",         "Events");
    2409       write.AddContainer("HadRF",         "Events");
    2410 
    2411       write.AddContainer("MEnergyEst",    "Events");
    2412 
    24132715
    24142716    //-----------------------------------------------------------------
    … …  
    24312733    fillhadrf.SetName("HhadRF");
    24322734
     2735    //---------------------------
     2736    // calculate estimated energy
     2737
     2738    MEnergyEstParam eeston(fHilName);
     2739    eeston.Add(fHilNameSrc);
     2740
     2741    eeston.SetCoeffA(parA);
     2742    eeston.SetCoeffB(parB);
     2743
     2744    //---------------------------
     2745    // calculate estimated energy using Daniel's parameters
     2746
     2747    //MEnergyEstParamDanielMkn421 eeston(fHilName);
     2748    //eeston.Add(fHilNameSrc);
     2749    //eeston.SetCoeffA(parA);
     2750    //eeston.SetCoeffB(parB);
     2751
     2752
     2753    //---------------------------
     2754
    24332755
    24342756    MFillH hfill1("MHHillas",    fHilName);
    … …  
    24462768    MFillH hfill5("MHNewImagePar", fImgParName);
    24472769    hfill5.SetName("HNewImagePar");   
     2770
     2771    //---------------------------
     2772    // new from Robert
     2773
     2774    MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
     2775    hfill6.SetName("HTimeDiffTheta");
     2776
     2777    MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
     2778    hfill6a.SetName("HTimeDiffTime");
     2779
     2780    MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
     2781    hfill7.SetName("HAlphaEnergyTheta");
     2782
     2783    MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
     2784    hfill7a.SetName("HAlphaEnergyTime");
     2785
     2786    MFillH hfill7b("MHThetabarTime", fHilNameSrc);
     2787    hfill7b.SetName("HThetabarTime");
     2788
     2789    MFillH hfill7c("MHEnergyTime", "MMcEvt");
     2790    hfill7c.SetName("HEnergyTime");
     2791
     2792
     2793    //---------------------------
    24482794
    24492795    MFCT1SelFinal selfinal(fHilNameSrc);
    … …  
    24672813    tliston.AddToList(&read);
    24682814
     2815    // robert     
     2816    tliston.AddToList(&hfill6);   //timediff
     2817    tliston.AddToList(&hfill6a);   //timediff
     2818
    24692819    tliston.AddToList(&contfinalgh);
    2470 
    2471     if (WXX)
    2472       tliston.AddToList(&write);
     2820    tliston.AddToList(&eeston);
     2821
     2822    tliston.AddToList(&write);
    24732823
    24742824    //tliston.AddToList(&fillhadnn);
    … …  
    24822832    tliston.AddToList(&hfill5);
    24832833
     2834    //robert
     2835    tliston.AddToList(&hfill7);
     2836    tliston.AddToList(&hfill7a);
     2837    tliston.AddToList(&hfill7b);
     2838    tliston.AddToList(&hfill7c);
     2839
    24842840    tliston.AddToList(&contfinal);
    24852841
    … …  
    24952851
    24962852    Int_t maxevents = -1;
    2497     //Int_t maxevents = 10000;
    24982853    if ( !evtloop.Eventloop(maxevents) )
    24992854        return;
    … …  
    25072862
    25082863    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     2864
     2865    gLog << "before hadRF" << endl;
    25092866    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
     2867
     2868    gLog << "before hadSC" << endl;
    25102869    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    25112870
     2871    gLog << "before MHHillas" << endl;
    25122872    pliston.FindObject("MHHillas")->DrawClone();
     2873
     2874    gLog << "before MHHillasExt" << endl;
    25132875    pliston.FindObject("MHHillasExt")->DrawClone();
     2876
     2877    gLog << "before MHHillasSrc" << endl;
    25142878    pliston.FindObject("MHHillasSrc")->DrawClone();
     2879
     2880    gLog << "before MHNewImagePar" << endl;
    25152881    pliston.FindObject("MHNewImagePar")->DrawClone();
     2882
     2883    gLog << "before MHStarMap" << endl;
    25162884    pliston.FindObject("MHStarMap")->DrawClone();
    25172885
     2886    gLog << "before DeleteBinnings" << endl;
     2887
    25182888    DeleteBinnings(&pliston);
    25192889
    2520     gLog << "Macro CT1Analysis : End of Job F_XX" << endl;
     2890    gLog << "before Robert's code" << endl;
     2891
     2892
     2893//rwagner write all relevant histograms onto a file
     2894
     2895  if (WRobert)
     2896  {
     2897    gLog << "=======================================================" << endl;
     2898    gLog << "Write results onto file '" << filenameResults << "'" << endl;
     2899
     2900    TFile outfile(filenameResults,"recreate");
     2901
     2902    MHHillasSrc* hillasSrc =
     2903      (MHHillasSrc*)(pliston->FindObject("MHHillasSrc"));
     2904        TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
     2905    alphaHist->Write();
     2906    gLog << "Alpha plot has been written out" << endl;   
     2907
     2908
     2909    MHAlphaEnergyTheta* aetH =
     2910      (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
     2911    TH3D* aetHist = (TH3D*)(aetH->GetHist());
     2912    aetHist->SetName("aetHist");
     2913    aetHist->Write();
     2914    gLog << "AlphaEnergyTheta plot has been written out" << endl;   
     2915
     2916    MHAlphaEnergyTime* aetH2 =
     2917      (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime"));
     2918    TH3D* aetHist2 = (TH3D*)(aetH2->GetHist());
     2919    aetHist2->SetName("aetimeHist");
     2920//     aetHist2->DrawClone();
     2921    aetHist2->Write();
     2922    gLog << "AlphaEnergyTime plot has been written out" << endl;   
     2923
     2924    MHThetabarTime* tbt =
     2925      (MHThetabarTime*)(pliston->FindObject("MHThetabarTime"));
     2926    TProfile* tbtHist = (TProfile*)(tbt->GetHist());
     2927    tbtHist->SetName("tbtHist");
     2928    tbtHist->Write();
     2929    gLog << "ThetabarTime plot has been written out" << endl;   
     2930
     2931    MHEnergyTime* ent =
     2932      (MHEnergyTime*)(pliston->FindObject("MHEnergyTime"));
     2933    TH2D* entHist = (TH2D*)(ent->GetHist());
     2934    entHist->SetName("entHist");
     2935    entHist->Write();
     2936    gLog << "EnergyTime plot has been written out" << endl;   
     2937   
     2938    MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
     2939    TH2D* timeHist = (TH2D*)(time->GetHist());
     2940    timeHist->SetName("MHTimeDiffTheta");
     2941    timeHist->SetTitle("Time diffs");
     2942    timeHist->Write();
     2943    gLog << "TimeDiffTheta plot has been written out" << endl;   
     2944
     2945
     2946    MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
     2947    TH2D* timeHist2 = (TH2D*)(time2->GetHist());
     2948    timeHist2->SetName("MHTimeDiffTime");
     2949    timeHist2->SetTitle("Time diffs");
     2950    timeHist2->Write();
     2951    gLog << "TimeDiffTime plot has been written out" << endl;   
     2952
     2953//rwagner write also collareas to same file
     2954    collarea->GetHist()->Write();
     2955    collarea->GetHAll()->Write();
     2956    collarea->GetHSel()->Write();
     2957    gLog << "Effective collection areas have been written out" << endl;       
     2958
     2959//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data
     2960//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file
     2961
     2962    gLog << "before closing outfile" << endl;
     2963
     2964    //outfile.Close();
     2965    gLog << "Results were written onto file '" << filenameResults
     2966         << "'" << endl;
     2967    gLog << "=======================================================" << endl;
     2968  }
     2969
     2970  }
     2971
     2972    gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
    25212973    gLog << "=======================================================" << endl;
    25222974 }
    … …  
    25272979
    25282980
     2981
     2982
     2983
     2984
     2985
     2986
     2987
     2988
     2989
     2990
     2991
     2992
     2993
     2994
     2995
     2996
     2997
     2998
     2999
     3000
     3001
     3002
     3003
     3004
     3005
     3006
Note: See TracChangeset for help on using the changeset viewer.