Ignore:
Timestamp:
03/08/03 14:00:30 (22 years ago)
Author:
wittek
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1767 r1809  
    11void AnalyseCT1()
    22{
    3     const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc";
    4     const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc";
     3  //-----------------------------------------------
     4    //const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc";
     5    //const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc";
    56    //const char *mcfile = "/hd43/rwagner/CT1/PreprocData/mc_Gammas.preproc.0.6";
    67    //const char *mcfile = "/hd43/rwagner/mkn421_mc.preproc";
    7     const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc";
     8    //const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc";
     9
     10  //-----------------------------------------------
     11  const char *offfile = "~/Mars200203/CT1data/offdata.preproc";
     12
     13  const char *onfile  = "~/Mars200203/CT1data/mkn421_on.preproc";
     14  //const char *onfile  = "~/Mars200203/CT1data/Crab_preproc_daniel_0.6";
     15
     16  //const char *mcfile  = "~/Mars200203/CT1data/mc_gammas.preproc"; //Version 0.5
     17  //const char *mcfile  = "~/Mars200203/CT1data/mc_Gammas.preproc.0.6";
     18  const char *mcfile = "~/Mars200203/CT1data/mkn421_mc_pedrms_0.001.preproc";
     19
     20  //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6";
     21  //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6_040303";
     22
     23  //-----------------------------------------------
     24  //const char *offfile = "~magican/home/ct1test/PreprocData/offdata.preproc";
     25  //const char *onfile  = "~magican/home/ct1test/PreprocData/mkn421_on.preproc";
     26  //const char *mcfile  = "~magican/home/ct1test/PreprocData/mc_gammas.preproc";
     27
     28
    829    const char *filename;
    930
    10     Bool_t Data  = kFALSE;   // loop over data ?
    11     Bool_t WLook = kFALSE;   // write out look-alike events for hadrons ?
    12     Bool_t WHist = kFALSE;   // write out histograms for padding ?
    13 
    14     Bool_t MC      = kTRUE;  // loop over MC gamma data ?
    15     Bool_t RHist   = kTRUE;  // read in histograms for padding ?
    16     Bool_t WLookMC = kTRUE;  // write out look-alike events for gammas ?
    17 
    18 
    19     cout << "" << endl;
    20     cout << "Macro AnalyseCT1 : Data, WHist = " << Data << ",  "
    21          << WHist << endl;
    22 
    23     cout << "                         RHist, MC   = " << RHist << ",  "
    24          << MC << endl;
    25     cout << "" << endl;
     31    //************************************************************************
     32    // Job 1 : read OFF data 
     33    // (generate sigmabar vs. Theta plot; write root file for OFF data (OFF1);
     34    //  generate hadron matrix for g/h separation)
     35
     36    Bool_t Job1     = kFALSE;   
     37    Bool_t WHistOFF = kTRUE;   // write out histogram sigmabar vs. Theta ?
     38    Bool_t WLookOFF = kTRUE;   // write out look-alike events for hadrons ?
     39    Bool_t WOFF1    = kTRUE;   // write out root file OFF1 ?
     40
     41
     42    // Job 2 : read ON data
     43    // (generate sigmabar vs. Theta plot; write root file for ON data (ON1))
     44
     45    Bool_t Job2    = kFALSE; 
     46    Bool_t WHistON = kFALSE;    // write out histogram sigmabar vs. Theta ?
     47    Bool_t WLookON = kFALSE;   // write out look-alike events for hadrons ?
     48    Bool_t WON1    = kFALSE;   // write out root file ON1 ?
     49
     50    // Job 3 : read MC gamma data, read sigmabar vs. Theta plot from OFF data 
     51    // (do padding; write root file for MC gammas (MC1);
     52    //  generate gamma matrix for g/h separation)
     53
     54    Bool_t Job3     = kFALSE; 
     55    Bool_t WLookMC  = kTRUE;  // write out look-alike events for gammas ?
     56    Bool_t WMC1     = kTRUE;  // write out root file MC1 ?
     57
     58
     59
     60    // Job 4 : read OFF1 and MC1 files, read hadron and gamma matrix
     61    // (produce the Neyman-Pearson plot)
     62    Bool_t Job4  = kTRUE; 
     63
     64
     65    // Job 5 : read MC1 file, read hadron and gamma matrix
     66    // (do g/h separation; write root file for MC gammas after final cuts (MC2))
     67    Bool_t Job5   = kFALSE; 
     68    Bool_t WMC2   = kFALSE;  // write out root file MC2 ?
     69
     70
     71    // Job 6 : read ON data, read hadron and gamma matrix
     72    // (do g/h separation; write root file for ON data after final cuts (ON2))
     73    Bool_t Job6  = kFALSE; 
     74    Bool_t WON2  = kTRUE;  // write out root file ON2 ?
     75    //************************************************************************
     76
     77
     78
    2679
    2780  //---------------------------------------------------------------------
    28   // start of section for ON data
    29   // (in the CT1 analysis the ON data are also used as a sample representing
    30   //  the hadrons for the g/h separation) 
    31 
    32     // read ON data file
    33 
    34     //  - to produce the 2D-histogram "sigmabar versus Theta"
     81    // Job 1
     82    //======
     83
     84    // read OFF data file
     85
     86    //  - produce the 2D-histogram "sigmabar versus Theta" for OFF data
    3587    //    (to be used for the padding of the MC gamma data)
    3688
    37     //  - to write a file of look-alike events (for hadrons)
     89    //  - write a file of look-alike events (hadron matrix)
    3890    //    (to be used later together with the corresponding MC gamma file
    3991    //     for the g/h separation in the analysis of the ON data)
    4092
    41     //  - to write a file of ON events (after the standard cuts,
    42     //                                  before the g/h separation)
     93    //  - write a file of OFF events (OFF1)
     94    //    (after the standard cuts, before the g/h separation)
    4395    //    (to be used together with the corresponding MC gamma file
    4496    //     for the optimization of the g/h separation)
    4597
    46  if (Data)
     98 if (Job1)
    4799 {
    48100    cout << "=====================================================" << endl;
    49     cout << "Macro AnalyseCT1 : Start of section for ON data" << endl;
     101    cout << "Macro AnalyseCT1 : Start of Job 1" << endl;
     102
     103    cout << "" << endl;
     104    cout << "Macro AnalyseCT1 : Job1, WHistOFF, WLookOFF, WOFF1 = "
     105         << Job1  << ",  " << WHistOFF << ",  " << WLookOFF << ",  "
     106         << WOFF1 << endl;
     107
     108
     109    TString typeData = "OFF";
     110    cout << "typeData = " << typeData << endl;
    50111
    51112    MTaskList tliston;
    … …  
    112173    //
    113174
    114     filename = onfile;
    115     readname = "ReadCT1data";
    116     MCT1ReadPreProc read(filename, readname);
     175    filename = offfile;
     176    RName = "ReadCT1data";
     177    MCT1ReadPreProc read(filename, RName);
    117178
    118179    MSelBasic selbasic;
    … …  
    150211    MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
    151212
    152     MSelStandard selstand(fHillas, fHillasSrc);
     213    MSelStandard selstand(fHilName, fHilSrcName);
    153214
    154215    MFillH hfill1("MHHillas",    fHilName);
    … …  
    168229
    169230    const char* mtxName = "MatrixHadrons";
    170     Int_t howMany = 30000;
    171     char* outName = "matrix_hadrons.root";
     231    Int_t howMany = 50000;
     232    TString outName = "matrix_hadrons_";
     233    outName += typeData;
     234    outName += ".root";
    172235
    173236    cout << "" << endl;
    … …  
    182245
    183246    // matrix limitation for look-alike events (approximate number)
    184     MFEventSelector selector("", "", readname);
     247    MFEventSelector selector("", "", RName);
    185248    selector.SetNumSelectEvts(howMany);
    186249    MFilterList flist;
    … …  
    192255    MHMatrix *pmatrix = new MHMatrix(mtxName);
    193256    MHMatrix& matrix = *pmatrix;
     257
     258    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
     259    matrix.AddColumn("MSigmabar.fSigmabar");
     260    matrix.AddColumn("log10(Hillas.fSize)");
     261    matrix.AddColumn("HillasSrc.fDist");
    194262    matrix.AddColumn("Hillas.fWidth");
    195263    matrix.AddColumn("Hillas.fLength");
    196     matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize");
    197     matrix.AddColumn("abs(Hillas.fAsym)");
     264    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
    198265    matrix.AddColumn("abs(Hillas.fM3Long)");
    199     matrix.AddColumn("abs(Hillas.fM3Trans)");
    200     matrix.AddColumn("abs(HillasSrc.fHeadTail)");
    201266    matrix.AddColumn("Hillas.fConc");
    202267    matrix.AddColumn("Hillas.fConc1");
    203     matrix.AddColumn("HillasSrc.fDist");
    204     matrix.AddColumn("log10(Hillas.fSize)");
    205     matrix.AddColumn("HillasSrc.fAlpha");
    206     matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
    207     //matrix.AddColumn("MMcEvt.fTheta");
     268
    208269    pliston.AddToList(pmatrix);
    209270
    … …  
    211272    fillmatg.SetFilter(&flist);
    212273
     274
     275    if (WOFF1)
     276    {   
     277      TString outNameImage = typeData;
     278      outNameImage += "1.root";
     279      MWriteRootFile write(outNameImage);
     280      write.AddContainer("MSigmabar",     "Events");
     281      write.AddContainer("Hillas",        "Events");
     282      write.AddContainer("MMcEvt",        "Events");
     283      write.AddContainer("HillasSrc",     "Events");
     284      write.AddContainer("MRawRunHeader", "RunHeaders");
     285      //write.AddContainer("MMcRunHeader","RunHeaders");
     286      write.AddContainer("MSrcPosCam",    "RunHeaders");
     287    }
     288
    213289    //+++++++++++++++++++++++++++++++++++++++++++++++++++
    214290
    215     MSelFinal selfinal(fHillas, fHillasSrc);
    216291
    217292    //*****************************
    … …  
    229304
    230305    tliston.AddToList(&selstand);
     306    if (WOFF1)
     307      tliston.AddToList(&write);
     308
     309    tliston.AddToList(&flist);
     310    tliston.AddToList(&fillmatg);   
     311
    231312    tliston.AddToList(&hfill1);
    232313    tliston.AddToList(&hfill2);
    … …  
    234315    tliston.AddToList(&hfill4);
    235316
    236     tliston.AddToList(&flist);
    237     tliston.AddToList(&fillmatg);
    238 
    239     //tliston.AddToList(&selfinal);
    240317    //*****************************
    241318
    … …  
    243320    // Execute event loop
    244321    //
     322    MProgressBar bar;
    245323    MEvtLoop evtloop;
    246324    evtloop.SetParList(&pliston);
     325    evtloop.SetProgressBar(&bar);
    247326
    248327    Int_t maxevents = 1000000000;
     328    //Int_t maxevents = 1000;
    249329    if ( !evtloop.Eventloop(maxevents) )
    250330        return;
    … …  
    255335    // Display the histograms
    256336    //
    257     pliston.FindObject("MHHillas")->DrawClone();
    258     pliston.FindObject("MHHillasSrc")->DrawClone();
     337    TObject *c;
     338    c = pliston.FindObject("MHHillas")->DrawClone();
     339
     340    c = pliston.FindObject("MHHillasSrc")->DrawClone();
    259341
    260342    //pliston.FindObject("MHHillasExt")->DrawClone();
    261     pliston.FindObject(nxt)->DrawClone();
    262 
    263     pliston.FindObject("MHStarMap")->DrawClone();
     343    c = pliston.FindObject(nxt)->DrawClone();
     344
     345    c = pliston.FindObject("MHStarMap")->DrawClone();
    264346
    265347
    … …  
    277359    Int_t nmaxevts  = 2000;
    278360
    279     // target distribution for the variable in column refcolumn;   
     361    // target distribution for the variable in column refcolumn (start at 1);
    280362    // set refcolumn negative if distribution is not to be changed
    281     Int_t   refcolumn = 12;
     363    Int_t   refcolumn = 1;
    282364    Int_t   nbins   =  10;
    283365    Float_t frombin = 0.5;
    … …  
    311393    // write out look-alike events (for hadrons)
    312394    //
    313     if (WLook)
     395    if (WLookOFF)
    314396    {
    315397      TFile *writejens = new TFile(outName, "RECREATE", "");
    … …  
    324406    //-------------------------------------------
    325407    // Write histograms onto a file
    326   if (WHist)
     408  if (WHistOFF)
    327409  {
    328410      MHSigmaTheta *sigtheta =
    … …  
    330412      if (!sigtheta)
    331413        {
    332           *fLog << "Object 'SigmaTheta' not found" << endl;
     414          cout << "Object 'SigmaTheta' not found" << endl;
    333415          return;
    334416        }
    … …  
    337419      TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
    338420
    339       TFile outfile("SigmaTheta.root", "RECREATE");
     421      TString outNameSigTh = "SigmaTheta_";
     422      outNameSigTh += typeData;
     423      outNameSigTh += ".root";
     424
     425      TFile outfile(outNameSigTh, "RECREATE");
    340426      fHSigTh->Write();
    341427      fHSigPixTh->Write();
    … …  
    343429      outfile.Close();
    344430     
    345       cout << "File 'SigmaTheta.root' was written out" << endl;
     431      cout << "File " << outNameSigTh << " was written out" << endl;
    346432  }
    347433
    348434
    349     cout << "Macro AnalyseCT1 : End of section for ON data" << endl;
     435    cout << "Macro AnalyseCT1 : End of Job 1" << endl;
    350436    cout << "===================================================" << endl;
    351437 }
    … …  
    353439
    354440  //---------------------------------------------------------------------
    355   // start of section for MC gamma data
     441  // Job 2
     442  //======
     443    // read ON data file
     444
     445    //  - produce the 2D-histogram "sigmabar versus Theta" for ON data
     446    //    (to be used for the padding of the MC gamma data)
     447
     448    //  - write a file of look-alike events (hadron matrix)
     449    //    (to be used later together with the corresponding MC gamma file
     450    //     for the g/h separation in the analysis of the ON data)
     451
     452    //  - write a file of ON events (ON1)
     453    //    (after the standard cuts, before the g/h separation)
     454    //    (to be used together with the corresponding MC gamma file
     455    //     for the optimization of the g/h separation)
     456
     457 if (Job2)
     458 {
     459    cout << "=====================================================" << endl;
     460    cout << "Macro AnalyseCT1 : Start of Job 2" << endl;
     461
     462    cout << "" << endl;
     463    cout << "Macro AnalyseCT1 : Job2, WHistON, WLookON, WON1 = "
     464         << Job2  << ",  " << WHistON << ",  " << WLookON << ",  "
     465         << WON1 << endl;
     466
     467    TString typeData = "ON";
     468    cout << "typeData = " << typeData << endl;
     469
     470    MTaskList tliston;
     471    MParList pliston;
     472    pliston.AddToList(&tliston);
     473
     474
     475    //::::::::::::::::::::::::::::::::::::::::::
     476
     477    MBinning binssize("BinningSize");
     478    binssize.SetEdgesLog(50, 10, 1.0e5);
     479    pliston.AddToList(&binssize);
     480
     481    MBinning binsdistc("BinningDist");
     482    binsdistc.SetEdges(50, 0, 1.4);
     483    pliston.AddToList(&binsdistc);
     484
     485    MBinning binswidth("BinningWidth");
     486    binswidth.SetEdges(50, 0, 1.0);
     487    pliston.AddToList(&binswidth);
     488
     489    MBinning binslength("BinningLength");
     490    binslength.SetEdges(50, 0, 1.0);
     491    pliston.AddToList(&binslength);
     492
     493    MBinning binsalpha("BinningAlpha");
     494    binsalpha.SetEdges(100, -100, 100);
     495    pliston.AddToList(&binsalpha);
     496
     497    MBinning binsasym("BinningAsym");
     498    binsasym.SetEdges(50, -1.5, 1.5);
     499    pliston.AddToList(&binsasym);
     500
     501    MBinning binsht("BinningHeadTail");
     502    binsht.SetEdges(50, -1.5, 1.5);
     503    pliston.AddToList(&binsht);
     504
     505    MBinning binsm3l("BinningM3Long");
     506    binsm3l.SetEdges(50, -1.5, 1.5);
     507    pliston.AddToList(&binsm3l);
     508
     509    MBinning binsm3t("BinningM3Trans");
     510    binsm3t.SetEdges(50, -1.5, 1.5);
     511    pliston.AddToList(&binsm3t);
     512
     513   
     514    //.....
     515    MBinning binsb("BinningSigmabar");
     516    binsb.SetEdges( 100,  0.0,  5.0);
     517    pliston.AddToList(&binsb);
     518
     519    MBinning binth("BinningTheta");
     520    binth.SetEdges( 70, -0.5, 69.5);   
     521    pliston.AddToList(&binth);
     522
     523    MBinning binsdiff("BinningDiffsigma2");
     524    binsdiff.SetEdges(100, -5.0, 20.0);
     525    pliston.AddToList(&binsdiff);
     526    //::::::::::::::::::::::::::::::::::::::::::
     527
     528
     529    //-------------------------------------------
     530    // create the tasks which should be executed
     531    //
     532
     533    filename = onfile;
     534    RName = "ReadCT1data";
     535    MCT1ReadPreProc read(filename, RName);
     536
     537    MSelBasic selbasic;
     538
     539    MBlindPixelCalc blind;
     540    blind.SetUseInterpolation();
     541    blind.SetUseBlindPixels();
     542    // blind.SetUseCetralPixel();
     543
     544    // create an object of class MSigmabar,
     545    // because class MHSigmaTheta will look for it
     546    MSigmabar sigbar;
     547    pliston.AddToList(&sigbar);
     548    MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
     549    fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
     550
     551    MImgCleanStd    clean;
     552
     553    // calculate also extended image parameters
     554    TString fHilName = "Hillas";
     555    MHillasExt hext(fHilName);
     556    pliston.AddToList(&hext);
     557    MHillasExt *fHillas = &hext;
     558
     559    // name of output container for MHillasCalc
     560    //      = name of MHillas object
     561    MHillasCalc    hcalc(fHilName);
     562
     563    // name of output container for MHillasSrcCalc
     564    //      = name of MHillasSrc object
     565    TString fHilSrcName = "HillasSrc";
     566    MHillasSrc hilsrc(fHilSrcName);
     567    MHillasSrc *fHillasSrc = &hilsrc;
     568    pliston.AddToList(&hilsrc);
     569    MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
     570
     571    MSelStandard selstand(fHilName, fHilSrcName);
     572
     573    MFillH hfill1("MHHillas",    fHilName);
     574    MFillH hfill2("MHStarMap",   fHilName);
     575
     576    TString nxt = "HillasExt";
     577    MHHillasExt hhilext(nxt, "", fHilName);
     578    pliston.AddToList(&hhilext);
     579    TString namext = nxt;
     580    namext += "[MHHillasExt]";
     581    MFillH hfill3(namext, fHilSrcName);
     582
     583    MFillH hfill4("MHHillasSrc", fHilSrcName);
     584
     585    //+++++++++++++++++++++++++++++++++++++++++++++++++++
     586    // prepare writing of look-alike events (for hadrons)
     587
     588    const char* mtxName = "MatrixHadrons";
     589    Int_t howMany = 50000;
     590    TString outName = "matrix_hadrons_";
     591    outName += typeData;
     592    outName += ".root";
     593
     594
     595    cout << "" << endl;
     596    cout << "========================================================" << endl;
     597    cout << "Matrix for (hadrons)" << endl;
     598    cout << "input file                    = " << filename<< endl;
     599    cout << "matrix name                   = " << mtxName << endl;
     600    cout << "max. no.of look-alike events  = " << howMany << endl;
     601    cout << "name of output root file      = " << outName << endl;
     602    cout << "" << endl;
     603
     604
     605    // matrix limitation for look-alike events (approximate number)
     606    MFEventSelector selector("", "", RName);
     607    selector.SetNumSelectEvts(howMany);
     608    MFilterList flist;
     609    flist.AddToList(&selector);
     610
     611    //
     612    // --- setup the matrix and add it to the parlist
     613    //
     614    MHMatrix *pmatrix = new MHMatrix(mtxName);
     615    MHMatrix& matrix = *pmatrix;
     616
     617    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
     618    matrix.AddColumn("MSigmabar.fSigmabar");
     619    matrix.AddColumn("log10(Hillas.fSize)");
     620    matrix.AddColumn("HillasSrc.fDist");
     621    matrix.AddColumn("Hillas.fWidth");
     622    matrix.AddColumn("Hillas.fLength");
     623    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
     624    matrix.AddColumn("abs(Hillas.fM3Long)");
     625    matrix.AddColumn("Hillas.fConc");
     626    matrix.AddColumn("Hillas.fConc1");
     627
     628    pliston.AddToList(pmatrix);
     629
     630    MFillH fillmatg(mtxName);
     631    fillmatg.SetFilter(&flist);
     632
     633
     634    if (WON1)
     635    {   
     636      TString outNameImage = typeData;
     637      outNameImage += "1.root";
     638      MWriteRootFile write(outNameImage);
     639      write.AddContainer("MSigmabar",     "Events");
     640      write.AddContainer("Hillas",        "Events");
     641      write.AddContainer("MMcEvt",        "Events");
     642      write.AddContainer("HillasSrc",     "Events");
     643      write.AddContainer("MRawRunHeader", "RunHeaders");
     644      //write.AddContainer("MMcRunHeader","RunHeaders");
     645      write.AddContainer("MSrcPosCam",    "RunHeaders");
     646    }
     647
     648    //+++++++++++++++++++++++++++++++++++++++++++++++++++
     649
     650
     651    //*****************************
     652    // tasks to be executed
     653   
     654    tliston.AddToList(&read);
     655
     656    tliston.AddToList(&selbasic);
     657    tliston.AddToList(&blind);
     658    tliston.AddToList(&fillsigtheta);
     659    tliston.AddToList(&clean);
     660
     661    tliston.AddToList(&hcalc);
     662    tliston.AddToList(&csrc1);
     663
     664    tliston.AddToList(&selstand);
     665    if (WON1)
     666      tliston.AddToList(&write);
     667
     668    tliston.AddToList(&flist);
     669    tliston.AddToList(&fillmatg);   
     670
     671    tliston.AddToList(&hfill1);
     672    tliston.AddToList(&hfill2);
     673    tliston.AddToList(&hfill3);
     674    tliston.AddToList(&hfill4);
     675
     676
     677    //*****************************
     678
     679    //-------------------------------------------
     680    // Execute event loop
     681    //
     682    MProgressBar bar;
     683    MEvtLoop evtloop;
     684    evtloop.SetParList(&pliston);
     685    evtloop.SetProgressBar(&bar);
     686
     687    Int_t maxevents = 1000000000;
     688    //Int_t maxevents = 1000;
     689    if ( !evtloop.Eventloop(maxevents) )
     690        return;
     691
     692    tliston.PrintStatistics(0, kTRUE);
     693
     694    //-------------------------------------------
     695    // Display the histograms
     696    //
     697    TObject *c;
     698    c = pliston.FindObject("MHHillas")->DrawClone();
     699
     700    c = pliston.FindObject("MHHillasSrc")->DrawClone();
     701
     702    //pliston.FindObject("MHHillasExt")->DrawClone();
     703    c = pliston.FindObject(nxt)->DrawClone();
     704
     705    c = pliston.FindObject("MHStarMap")->DrawClone();
     706
     707
     708    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     709    //
     710    // ----------------------------------------------------------
     711    //  Definition of the reference sample of look-alike events
     712    //  (this is a subsample of the original sample)
     713    // ----------------------------------------------------------
     714    //
     715    cout << "" << endl;
     716    cout << "========================================================" << endl;
     717    // Select a maximum of nmaxevts events from the sample of look-alike
     718    // events. They will form the reference sample.
     719    Int_t nmaxevts  = 2000;
     720
     721    // target distribution for the variable in column refcolumn (start at 1);
     722    // set refcolumn negative if distribution is not to be changed
     723    Int_t   refcolumn = 1;
     724    Int_t   nbins   =  10;
     725    Float_t frombin = 0.5;
     726    Float_t tobin   = 1.0;
     727    TH1F *thsh = new TH1F("thsh","target distribution",
     728                           nbins, frombin, tobin);
     729    thsh->SetXTitle("cos( \\Theta)");
     730    thsh->SetYTitle("Counts");
     731    Float_t dbin = (tobin-frombin)/nbins;
     732    Float_t lbin = frombin +dbin*0.5;
     733    for (Int_t j=1; j<=nbins; j++)
     734    {
     735      thsh->Fill(lbin,1.0);
     736      lbin += dbin;
     737    }
     738
     739    cout << "" << endl;
     740    cout << "========================================================" << endl;
     741    cout << "Macro AnalyseCT1 : define reference sample for hadrons" << endl;
     742    cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
     743         << refcolumn << ",  " << nmaxevts << endl;   
     744
     745    if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) )
     746    {
     747      cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl;
     748      return;
     749    }
     750    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     751
     752    //-------------------------------------------
     753    // write out look-alike events (for hadrons)
     754    //
     755    if (WLookON)
     756    {
     757      TFile *writejens = new TFile(outName, "RECREATE", "");
     758      matrix.Write();
     759      cout << "Macro AnalyseCT1 : matrix of look-alike events for hadrons written onto file "
     760           << outName << endl;
     761
     762      writejens->Close();
     763      delete writejens;
     764    }
     765
     766    //-------------------------------------------
     767    // Write histograms onto a file
     768  if (WHistON)
     769  {
     770      MHSigmaTheta *sigtheta =
     771            (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
     772      if (!sigtheta)
     773        {
     774          cout << "Object 'SigmaTheta' not found" << endl;
     775          return;
     776        }
     777      TH2D *fHSigTh    = sigtheta->GetSigmaTheta();
     778      TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta();
     779      TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
     780
     781      TString outNameSigTh = "SigmaTheta_";
     782      outNameSigTh += typeData;
     783      outNameSigTh += ".root";
     784
     785      TFile outfile(outNameSigTh, "RECREATE");
     786      fHSigTh->Write();
     787      fHSigPixTh->Write();
     788      fHDifPixTh->Write();
     789      outfile.Close();
     790     
     791      cout << "File " << outNameSigTh << " was written out" << endl;
     792  }
     793
     794
     795    cout << "Macro AnalyseCT1 : End of Job 2" << endl;
     796    cout << "===================================================" << endl;
     797 }
     798
     799
     800  //---------------------------------------------------------------------
     801   // Job 3
     802   //======
    356803
    357804    // read MC gamma data 
    … …  
    360807    //      (using the 2D-histogram "sigmabar versus Theta" of the OFF data)
    361808
    362     //    - to write a file of look-alike events (for gammas)
     809    //    - to write a file of look-alike events (gammas matrix)
    363810    //      (to be used later together with the corresponding hadron file
    364811    //       for the g/h separation in the analysis of the ON data)
    365812
    366     //    - to write a file of padded MC gamma events
     813    //    - to write a file of padded MC gamma events (MC1)
    367814    //      (after the standard cuts, before the g/h separation)
    368815    //      (to be used together with the corresponding hadron file
    369816    //       for the optimization of the g/h separation)
    370817
    371     //    - to write a file of padded MC gamma events (after the final cuts)
    372     //      (to be used for the optimization of the energy estimator
    373     //       and for calculating the effective collection areas)
    374 
    375  if (MC)
     818
     819 if (Job3)
    376820 {
    377     cout << "============================================================="
    378          << endl;
    379     cout << "Macro : AnalyseCT1 : Start of section for MC gamma data"
    380          << endl;
    381 
    382     MTaskList tlist;
    383     MParList plist;
    384 
    385     plist.AddToList(&tlist);
     821    cout << "=====================================================" << endl;
     822    cout << "Macro AnalyseCT1 : Start of Job 3" << endl;
     823
     824    cout << "" << endl;
     825    cout << "Macro AnalyseCT1 : Job3, WLookMC, WMC1 = "
     826         << Job3  << ",  " << WLookMC << ",  " << WMC1 << endl;
     827
     828
     829    TString typeMC = "MC";
     830    cout << "typeMC = " << typeMC << endl;
     831   
     832    // use for padding sigmabar vs. Theta from ON or OFF data
     833    TString typeData = "ON";
     834    //TString typeData = "OFF";
     835    cout << "typeData = " << typeData << endl;
    386836
    387837    //------------------------------------
    … …  
    389839    // and                "3D-ThetaPixSigma"
    390840    // and                "3D-ThetaPixDiff"
    391   if (RHist)
    392   {
    393       cout << "Reading in file 'SigmaTheta.root' " << endl;
    394 
    395       TFile *infile = new TFile("SigmaTheta.root");
     841
     842
     843      TString outNameSigTh = "SigmaTheta_";
     844      outNameSigTh += typeData;
     845      outNameSigTh += ".root";
     846
     847
     848      cout << "Reading in file " << outNameSigTh << endl;
     849
     850      TFile *infile = new TFile(outNameSigTh);
    396851      infile->ls();
    397852
    … …  
    400855      if (!fHSigmaTheta)
    401856        {
    402           *fLog << "Object '2D-ThetaSigmabar' not found on root file" << endl;
     857          cout << "Object '2D-ThetaSigmabar' not found on root file" << endl;
    403858          return;
    404859        }
    … …  
    409864      if (!fHSigmaPixTheta)
    410865        {
    411           *fLog << "Object '3D-ThetaPixSigma' not found on root file" << endl;
     866          cout << "Object '3D-ThetaPixSigma' not found on root file" << endl;
    412867          return;
    413868        }
    … …  
    418873      if (!fHSigmaPixTheta)
    419874        {
    420           *fLog << "Object '3D-ThetaPixDiff' not found on root file" << endl;
     875          cout << "Object '3D-ThetaPixDiff' not found on root file" << endl;
    421876          return;
    422877        }
    423878      cout << "Object '3D-ThetaPixDiff' was read in" << endl;
    424   }
    425   else
    426   {
    427       MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
    428       if (!sigtheta)
    429       {
    430         cout << "Object with name 'SigmaTheta' not found" << endl;
    431         return;
    432       }
    433 
    434       TH2D *fHSigmaTheta    = sigtheta->GetSigmaTheta();
    435       TH3D *fHSigmaPixTheta = sigtheta->GetSigmaPixTheta();
    436       TH3D *fHDiffPixTheta  = sigtheta->GetDiffPixTheta();
    437   }
     879
    438880    //------------------------------------
    439881
     882    MTaskList tlist;
     883    MParList plist;
     884
     885    plist.AddToList(&tlist);
    440886
    441887
    … …  
    539985    MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
    540986
    541     MSelStandard selstand(fHillas, fHillasSrc);
     987    MSelStandard selstand(fHilName, fHilSrcName);
    542988
    543989    MFillH hfill1("MHHillas",    fHilName);
    … …  
    5571003
    5581004    const char* mtxName = "MatrixGammas";
    559     Int_t howMany = 30000;
    560     char* outName = "matrix_gammas.root";
     1005    Int_t howMany = 50000;
     1006
     1007    TString outName = "matrix_gammas_";
     1008    outName += typeMC;
     1009    outName += "_";
     1010    outName += typeData;
     1011    outName += ".root";
     1012
    5611013
    5621014    cout << "" << endl;
    … …  
    5811033    MHMatrix *pmatrix = new MHMatrix(mtxName);
    5821034    MHMatrix& matrix = *pmatrix;
     1035
     1036    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
     1037    matrix.AddColumn("MSigmabar.fSigmabar");
     1038    matrix.AddColumn("log10(Hillas.fSize)");
     1039    matrix.AddColumn("HillasSrc.fDist");
    5831040    matrix.AddColumn("Hillas.fWidth");
    5841041    matrix.AddColumn("Hillas.fLength");
    585     matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize");
    586     matrix.AddColumn("abs(Hillas.fAsym)");
     1042    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
    5871043    matrix.AddColumn("abs(Hillas.fM3Long)");
    588     matrix.AddColumn("abs(Hillas.fM3Trans)");
    589     matrix.AddColumn("abs(HillasSrc.fHeadTail)");
    5901044    matrix.AddColumn("Hillas.fConc");
    5911045    matrix.AddColumn("Hillas.fConc1");
    592     matrix.AddColumn("HillasSrc.fDist");
    593     matrix.AddColumn("log10(Hillas.fSize)");
    594     matrix.AddColumn("HillasSrc.fAlpha");
    595     matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
    596     //matrix.AddColumn("MMcEvt.fTheta");
     1046
    5971047    plist.AddToList(pmatrix);
    5981048
    … …  
    6001050    fillmatg.SetFilter(&flist);
    6011051
     1052
     1053    if (WMC1)
     1054    {   
     1055      TString outNameImage = typeMC;
     1056      outNameImage += "_";
     1057      outNameImage += typeData;
     1058      outNameImage += "1.root";
     1059      MWriteRootFile write(outNameImage);
     1060      write.AddContainer("MSigmabar",     "Events");
     1061      write.AddContainer("Hillas",        "Events");
     1062      write.AddContainer("MMcEvt",        "Events");
     1063      write.AddContainer("HillasSrc",     "Events");
     1064      write.AddContainer("MRawRunHeader", "RunHeaders");
     1065      //write.AddContainer("MMcRunHeader","RunHeaders");
     1066      write.AddContainer("MSrcPosCam",    "RunHeaders");
     1067    }
     1068
     1069
    6021070    //+++++++++++++++++++++++++++++++++++++++++++++++++++
    6031071
    604 
    605     MSelFinal selfinal(fHillas, fHillasSrc);
    6061072
    6071073    //*****************************
    … …  
    6191085
    6201086    tlist.AddToList(&selstand);
     1087    if (WMC1)
     1088      tlist.AddToList(&write);
     1089
     1090    tlist.AddToList(&flist);
     1091    tlist.AddToList(&fillmatg);
     1092
    6211093    tlist.AddToList(&hfill1);
    6221094    tlist.AddToList(&hfill2);
    … …  
    6241096    tlist.AddToList(&hfill4);
    6251097
    626     tlist.AddToList(&flist);
    627     tlist.AddToList(&fillmatg);
    628 
    629     //tlist.AddToList(&selfinal);
     1098
    6301099    //*****************************
    6311100
    … …  
    6341103    // Execute event loop
    6351104    //
     1105    MProgressBar bar;
    6361106    MEvtLoop evtloop;
    6371107    evtloop.SetParList(&plist);
     1108    evtloop.SetProgressBar(&bar);
    6381109
    6391110    Int_t maxevents = 1000000000;
     1111    //Int_t maxevents = 100;
    6401112    if ( !evtloop.Eventloop(maxevents) )
    6411113        return;
    … …  
    6701142    // target distribution for the variable in column refcolumn;   
    6711143    // set refcolumn negative if distribution is not to be changed
    672     Int_t   refcolumn = 12;
     1144    Int_t   refcolumn = 1;
    6731145    Int_t   nbins   =  10;
    6741146    Float_t frombin = 0.5;
    … …  
    7111183    }
    7121184
    713     cout << "Macro AnalyseCT1 : End of section for MC gamma data"
     1185    cout << "Macro AnalyseCT1 : End of Job 3"
    7141186         << endl;
    7151187    cout << "========================================================="
    7161188         << endl;
    7171189 }
     1190
     1191
     1192
     1193  //---------------------------------------------------------------------
     1194  // Job 4
     1195  //======
     1196
     1197    // read OFF1 and MC1 data files 
     1198    //
     1199    //  - produce Neyman-Pearson plot
     1200 
     1201 if (Job4)
     1202 {
     1203    cout << "=====================================================" << endl;
     1204    cout << "Macro AnalyseCT1 : Start of Job 4" << endl;
     1205
     1206    cout << "" << endl;
     1207    cout << "Macro AnalyseCT1 : Job4 = " << Job4  << endl;
     1208
     1209
     1210    TString typeMC = "MC";
     1211    cout << "typeMC = " << typeMC << endl;
     1212
     1213    // use hadron matrix from ON or OFF data
     1214    TString typeData = "ON";
     1215    //TString typeData = "OFF";
     1216    cout << "typeData = " << typeData << endl;
     1217
     1218   //*************************************************************************
     1219   // read in matrices of look-alike events
     1220
     1221    const char* mtxName = "MatrixGammas";
     1222
     1223
     1224    TString outName = "matrix_gammas_";
     1225    outName += typeMC;
     1226    outName += "_";
     1227    outName += typeData;
     1228    outName += ".root";
     1229
     1230    cout << "" << endl;
     1231    cout << "========================================================" << endl;
     1232    cout << "Get matrix for (gammas)" << endl;
     1233    cout << "matrix name        = " << mtxName << endl;
     1234    cout << "name of root file  = " << outName << endl;
     1235    cout << "" << endl;
     1236
     1237
     1238    // read in the object with the name 'mtxName' from file 'outName'
     1239    //
     1240    TFile *fileg = new TFile(outName);
     1241
     1242    MHMatrix matrixg;
     1243    matrixg.Read(mtxName);
     1244    pmatrixg = &matrixg;
     1245
     1246    delete fileg;
     1247
     1248
     1249    //*****************************************************************
     1250
     1251    const char* mtxName = "MatrixHadrons";
     1252
     1253    TString outName = "matrix_hadrons_";
     1254    outName += typeData;
     1255    outName += ".root";
     1256
     1257    cout << "" << endl;
     1258    cout << "========================================================" << endl;
     1259    cout << " Get matrix for (hadrons)" << endl;
     1260    cout << "matrix name        = " << mtxName << endl;
     1261    cout << "name of root file  = " << outName << endl;
     1262    cout << "" << endl;
     1263
     1264
     1265    // read in the object with the name mtxName
     1266    //
     1267    TFile *fileh = new TFile(outName);
     1268
     1269    MHMatrix matrixh;
     1270    matrixh.Read(mtxName);
     1271    pmatrixh = &matrixh;
     1272
     1273    delete fileh;
     1274
     1275    //*****************************************************************
     1276
     1277    MTaskList tliston;
     1278    MParList pliston;
     1279    pliston.AddToList(&tliston);
     1280
     1281    pliston.AddToList(pmatrixg);
     1282    pliston.AddToList(pmatrixh);
     1283
     1284    //::::::::::::::::::::::::::::::::::::::::::
     1285
     1286    MBinning binssize("BinningSize");
     1287    binssize.SetEdgesLog(50, 10, 1.0e5);
     1288    pliston.AddToList(&binssize);
     1289
     1290    MBinning binsdistc("BinningDist");
     1291    binsdistc.SetEdges(50, 0, 1.4);
     1292    pliston.AddToList(&binsdistc);
     1293
     1294    MBinning binswidth("BinningWidth");
     1295    binswidth.SetEdges(50, 0, 1.0);
     1296    pliston.AddToList(&binswidth);
     1297
     1298    MBinning binslength("BinningLength");
     1299    binslength.SetEdges(50, 0, 1.0);
     1300    pliston.AddToList(&binslength);
     1301
     1302    MBinning binsalpha("BinningAlpha");
     1303    binsalpha.SetEdges(100, -100, 100);
     1304    pliston.AddToList(&binsalpha);
     1305
     1306    MBinning binsasym("BinningAsym");
     1307    binsasym.SetEdges(50, -1.5, 1.5);
     1308    pliston.AddToList(&binsasym);
     1309
     1310    MBinning binsht("BinningHeadTail");
     1311    binsht.SetEdges(50, -1.5, 1.5);
     1312    pliston.AddToList(&binsht);
     1313
     1314    MBinning binsm3l("BinningM3Long");
     1315    binsm3l.SetEdges(50, -1.5, 1.5);
     1316    pliston.AddToList(&binsm3l);
     1317
     1318    MBinning binsm3t("BinningM3Trans");
     1319    binsm3t.SetEdges(50, -1.5, 1.5);
     1320    pliston.AddToList(&binsm3t);
     1321
     1322   
     1323    //.....
     1324    MBinning binsb("BinningSigmabar");
     1325    binsb.SetEdges( 100,  0.0,  5.0);
     1326    pliston.AddToList(&binsb);
     1327
     1328    MBinning binth("BinningTheta");
     1329    binth.SetEdges( 70, -0.5, 69.5);   
     1330    pliston.AddToList(&binth);
     1331
     1332    MBinning binsdiff("BinningDiffsigma2");
     1333    binsdiff.SetEdges(100, -5.0, 20.0);
     1334    pliston.AddToList(&binsdiff);
     1335    //::::::::::::::::::::::::::::::::::::::::::
     1336
     1337
     1338    //-------------------------------------------
     1339    // create the tasks which should be executed
     1340    //
     1341
     1342    TString filenameData = typeData;
     1343    filenameData += "1.root";
     1344
     1345    TString filenameMC = typeMC;
     1346    filenameMC += "_";
     1347    filenameMC += typeData;
     1348    filenameMC += "1.root";
     1349   
     1350
     1351    cout << "filenameData = " << filenameData << endl;
     1352    cout << "filenameMC   = " << filenameMC   << endl;
     1353
     1354    MReadMarsFile read("Events", filenameMC);
     1355    read.AddFile(filenameData);
     1356    read.DisableAutoScheme();
     1357
     1358    //.......................................................................
     1359    // g/h separation
     1360
     1361    MMultiDimDistCalc ghsep;
     1362    ghsep.SetUseNumRows(25);
     1363    ghsep.SetUseKernelMethod(kFALSE);
     1364
     1365    //.......................................................................
     1366
     1367    // geometry is needed in  MHHillas... classes
     1368    MGeomCam *fGeom =
     1369             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     1370
     1371    TString fHilName    = "Hillas";
     1372    TString fHilSrcName = "HillasSrc";
     1373
     1374    Float_t hadcut = 0.2;
     1375    MSelFinal selfinalgh(fHilName, fHilSrcName);
     1376    selfinalgh.SetHadronnessCut(hadcut);
     1377    selfinalgh.SetAlphaCut(100.0);
     1378
     1379    MFillH fillh("MHHadronness");
     1380
     1381    MSelFinal selfinal(fHilName, fHilSrcName);
     1382    selfinal.SetHadronnessCut(hadcut);
     1383    selfinal.SetAlphaCut(20.0);
     1384
     1385    MFillH hfill1("MHHillas",    fHilName);
     1386    MFillH hfill2("MHStarMap",   fHilName);
     1387
     1388    TString nxt = "HillasExt";
     1389    MHHillasExt hhilext(nxt, "", fHilName);
     1390    pliston.AddToList(&hhilext);
     1391    TString namext = nxt;
     1392    namext += "[MHHillasExt]";
     1393    MFillH hfill3(namext, fHilSrcName);
     1394
     1395    MFillH hfill4("MHHillasSrc", fHilSrcName);
     1396
     1397
     1398    //*****************************
     1399    // tasks to be executed
     1400   
     1401    tliston.AddToList(&read);
     1402
     1403    tliston.AddToList(&ghsep);
     1404    tliston.AddToList(&fillh);
     1405   
     1406    tliston.AddToList(&selfinalgh);
     1407    tliston.AddToList(&hfill1);
     1408    tliston.AddToList(&hfill2);
     1409    tliston.AddToList(&hfill3);
     1410    tliston.AddToList(&hfill4);
     1411
     1412    tliston.AddToList(&selfinal);
     1413
     1414    //*****************************
     1415
     1416    //-------------------------------------------
     1417    // Execute event loop
     1418    //
     1419    MProgressBar bar;
     1420    MEvtLoop evtloop;
     1421    evtloop.SetParList(&pliston);
     1422    evtloop.SetProgressBar(&bar);
     1423
     1424    Int_t maxevents = 1000000000;
     1425    //Int_t maxevents = 1000;
     1426    if ( !evtloop.Eventloop(maxevents) )
     1427        return;
     1428
     1429    tliston.PrintStatistics(0, kTRUE);
     1430
     1431
     1432    //-------------------------------------------
     1433    // Display the histograms
     1434    //
     1435    TObject *c;
     1436
     1437    c = pliston.FindObject("MHHadronness")->DrawClone();
     1438    c->Print();
     1439
     1440    //c = pliston.FindObject("MHHillas")->DrawClone();
     1441
     1442    c = pliston.FindObject("MHHillasSrc")->DrawClone();
     1443
     1444    //pliston.FindObject("MHHillasExt")->DrawClone();
     1445    //c = pliston.FindObject(nxt)->DrawClone();
     1446
     1447    c = pliston.FindObject("MHStarMap")->DrawClone();
     1448
     1449
     1450
     1451    cout << "Macro AnalyseCT1 : End of Job 4" << endl;
     1452    cout << "===================================================" << endl;
     1453 }
     1454
     1455
     1456  //---------------------------------------------------------------------
     1457  // Job 5
     1458  //======
     1459
     1460    //  - read in the matrices of look-alike events for gammas and hadrons
     1461
     1462    // then read MC1 data file
     1463    //
     1464    //  - perform the g/h separation
     1465    //  - apply the final cuts
     1466    //
     1467    //  - write a file of MC gamma events (MC2)
     1468    //    (after the g/h separation and after the final cuts)
     1469
     1470 if (Job5)
     1471 {
     1472    cout << "=====================================================" << endl;
     1473    cout << "Macro AnalyseCT1 : Start of Job 5" << endl;
     1474
     1475    cout << "" << endl;
     1476    cout << "Macro AnalyseCT1 : Job5, WMC2 = "
     1477         << Job5  << ",  " << WMC2 << endl;
     1478
     1479    TString typeInput = "MC";
     1480    cout << "typeInput = " << typeInput << endl;
     1481
     1482    TString typeMC   = "MC";
     1483    cout << "typeMC = " << typeMC << endl;
     1484
     1485    // use hadron matrix from ON or OFF data
     1486    TString typeData = "ON";
     1487    //TString typeData = "OFF";
     1488    cout << "typeData = " << typeData << endl;
     1489
     1490   //*************************************************************************
     1491   // read in matrices of look-alike events
     1492
     1493    const char* mtxName = "MatrixGammas";
     1494
     1495    TString outName = "matrix_gammas_";
     1496    outName += typeMC;
     1497    outName += "_";
     1498    outName += typeData;
     1499    outName += ".root";
     1500
     1501
     1502    cout << "" << endl;
     1503    cout << "========================================================" << endl;
     1504    cout << "Get matrix for (gammas)" << endl;
     1505    cout << "matrix name        = " << mtxName << endl;
     1506    cout << "name of root file  = " << outName << endl;
     1507    cout << "" << endl;
     1508
     1509
     1510    // read in the object with the name 'mtxName' from file 'outName'
     1511    //
     1512    TFile *fileg = new TFile(outName);
     1513
     1514    MHMatrix matrixg;
     1515    matrixg.Read(mtxName);
     1516    pmatrixg = &matrixg;
     1517
     1518    delete fileg;
     1519
     1520    //*****************************************************************
     1521
     1522    const char* mtxName = "MatrixHadrons";
     1523
     1524    TString outName = "matrix_hadrons_";
     1525    outName += typeData;
     1526    outName += ".root";
     1527
     1528
     1529    cout << "" << endl;
     1530    cout << "========================================================" << endl;
     1531    cout << " Get matrix for (hadrons)" << endl;
     1532    cout << "matrix name        = " << mtxName << endl;
     1533    cout << "name of root file  = " << outName << endl;
     1534    cout << "" << endl;
     1535
     1536
     1537    // read in the object with the name mtxName
     1538    //
     1539    TFile *fileh = new TFile(outName);
     1540
     1541    MHMatrix matrixh;
     1542    matrixh.Read(mtxName);
     1543    pmatrixh = &matrixh;
     1544
     1545    delete fileh;
     1546
     1547   //*************************************************************************
     1548
     1549
     1550    MTaskList tliston;
     1551    MParList pliston;
     1552    pliston.AddToList(&tliston);
     1553
     1554    pliston.AddToList(pmatrixg);
     1555    pliston.AddToList(pmatrixh);
     1556
     1557    //::::::::::::::::::::::::::::::::::::::::::
     1558
     1559    MBinning binssize("BinningSize");
     1560    binssize.SetEdgesLog(50, 10, 1.0e5);
     1561    pliston.AddToList(&binssize);
     1562
     1563    MBinning binsdistc("BinningDist");
     1564    binsdistc.SetEdges(50, 0, 1.4);
     1565    pliston.AddToList(&binsdistc);
     1566
     1567    MBinning binswidth("BinningWidth");
     1568    binswidth.SetEdges(50, 0, 1.0);
     1569    pliston.AddToList(&binswidth);
     1570
     1571    MBinning binslength("BinningLength");
     1572    binslength.SetEdges(50, 0, 1.0);
     1573    pliston.AddToList(&binslength);
     1574
     1575    MBinning binsalpha("BinningAlpha");
     1576    binsalpha.SetEdges(100, -100, 100);
     1577    pliston.AddToList(&binsalpha);
     1578
     1579    MBinning binsasym("BinningAsym");
     1580    binsasym.SetEdges(50, -1.5, 1.5);
     1581    pliston.AddToList(&binsasym);
     1582
     1583    MBinning binsht("BinningHeadTail");
     1584    binsht.SetEdges(50, -1.5, 1.5);
     1585    pliston.AddToList(&binsht);
     1586
     1587    MBinning binsm3l("BinningM3Long");
     1588    binsm3l.SetEdges(50, -1.5, 1.5);
     1589    pliston.AddToList(&binsm3l);
     1590
     1591    MBinning binsm3t("BinningM3Trans");
     1592    binsm3t.SetEdges(50, -1.5, 1.5);
     1593    pliston.AddToList(&binsm3t);
     1594
     1595   
     1596    //.....
     1597    MBinning binsb("BinningSigmabar");
     1598    binsb.SetEdges( 100,  0.0,  5.0);
     1599    pliston.AddToList(&binsb);
     1600
     1601    MBinning binth("BinningTheta");
     1602    binth.SetEdges( 70, -0.5, 69.5);   
     1603    pliston.AddToList(&binth);
     1604
     1605    MBinning binsdiff("BinningDiffsigma2");
     1606    binsdiff.SetEdges(100, -5.0, 20.0);
     1607    pliston.AddToList(&binsdiff);
     1608    //::::::::::::::::::::::::::::::::::::::::::
     1609
     1610
     1611    //-------------------------------------------
     1612    // create the tasks which should be executed
     1613    //
     1614
     1615    TString filenameMC = typeInput;
     1616    filenameMC += "_";
     1617    filenameMC += typeData;
     1618    filenameMC += "1.root";
     1619
     1620    readname = "ReadCT1MCdata";
     1621    MReadMarsFile read("Events", filenameMC);
     1622    read.DisableAutoScheme();
     1623
     1624
     1625    //.......................................................................
     1626    // g/h separation
     1627
     1628    MMultiDimDistCalc ghsep;
     1629    ghsep.SetUseNumRows(25);
     1630    ghsep.SetUseKernelMethod(kFALSE);
     1631
     1632
     1633
     1634    //.......................................................................
     1635
     1636    // geometry is needed in  MHHillas... classes
     1637    MGeomCam *fGeom =
     1638             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     1639
     1640    TString fHilName    = "Hillas";
     1641    TString fHilSrcName = "HillasSrc";
     1642
     1643    Float_t hadcut = 0.2;
     1644    MSelFinal selfinalgh(fHilName, fHilSrcName);
     1645    selfinalgh.SetHadronnessCut(hadcut);
     1646    selfinalgh.SetAlphaCut(100.0);
     1647
     1648    MFillH fillh("MHHadronness");
     1649
     1650    MSelFinal selfinal(fHilName, fHilSrcName);
     1651    selfinal.SetHadronnessCut(hadcut);
     1652    selfinal.SetAlphaCut(20.0);
     1653
     1654    if (WMC2)
     1655    {   
     1656      TString outNameImage = typeInput;
     1657      outNameImage += "_";
     1658      outNameImage += typeData;
     1659      outNameImage += "2.root";
     1660      MWriteRootFile write(outNameImage);
     1661      write.AddContainer("MHadronness",   "Events");
     1662      write.AddContainer("MSigmabar",     "Events");
     1663      write.AddContainer("Hillas",        "Events");
     1664      write.AddContainer("MMcEvt",        "Events");
     1665      write.AddContainer("HillasSrc",     "Events");
     1666      write.AddContainer("MRawRunHeader", "RunHeaders");
     1667      //write.AddContainer("MMcRunHeader","RunHeaders");
     1668      write.AddContainer("MSrcPosCam",    "RunHeaders");
     1669    }
     1670
     1671
     1672    MFillH hfill1("MHHillas",    fHilName);
     1673    MFillH hfill2("MHStarMap",   fHilName);
     1674
     1675    TString nxt = "HillasExt";
     1676    MHHillasExt hhilext(nxt, "", fHilName);
     1677    pliston.AddToList(&hhilext);
     1678    TString namext = nxt;
     1679    namext += "[MHHillasExt]";
     1680    MFillH hfill3(namext, fHilSrcName);
     1681
     1682    MFillH hfill4("MHHillasSrc", fHilSrcName);
     1683
     1684
     1685
     1686    //*****************************
     1687    // tasks to be executed
     1688   
     1689    tliston.AddToList(&read);
     1690
     1691    tliston.AddToList(&ghsep);
     1692    tliston.AddToList(&fillh);
     1693
     1694    tliston.AddToList(&selfinalgh);
     1695    tliston.AddToList(&hfill1);
     1696    tliston.AddToList(&hfill2);
     1697    tliston.AddToList(&hfill3);
     1698    tliston.AddToList(&hfill4);
     1699
     1700    tliston.AddToList(&selfinal);
     1701    if (WMC2)
     1702      tliston.AddToList(&write);
     1703
     1704    //*****************************
     1705
     1706    //-------------------------------------------
     1707    // Execute event loop
     1708    //
     1709    MProgressBar bar;
     1710    MEvtLoop evtloop;
     1711    evtloop.SetParList(&pliston);
     1712    evtloop.SetProgressBar(&bar);
     1713
     1714    Int_t maxevents = 1000000000;
     1715    //Int_t maxevents = 1000;
     1716    if ( !evtloop.Eventloop(maxevents) )
     1717        return;
     1718
     1719    tliston.PrintStatistics(0, kTRUE);
     1720
     1721
     1722    //-------------------------------------------
     1723    // Display the histograms
     1724    //
     1725    TObject *c;
     1726
     1727    c = pliston.FindObject("MHHadronness")->DrawClone();
     1728    c->Print();
     1729
     1730    c = pliston.FindObject("MHHillas")->DrawClone();
     1731
     1732    c = pliston.FindObject("MHHillasSrc")->DrawClone();
     1733
     1734    //pliston.FindObject("MHHillasExt")->DrawClone();
     1735    c = pliston.FindObject(nxt)->DrawClone();
     1736
     1737    c = pliston.FindObject("MHStarMap")->DrawClone();
     1738
     1739
     1740
     1741    cout << "Macro AnalyseCT1 : End of Job 5" << endl;
     1742    cout << "===================================================" << endl;
     1743 }
     1744
     1745
     1746  //---------------------------------------------------------------------
     1747  // Job 6
     1748  //======
     1749
     1750    //  - read in the matrices of look-alike events for gammas and hadrons
     1751
     1752    // then read ON1 data file
     1753    //
     1754    //  - perform the g/h separation
     1755    //  - apply the final cuts
     1756    //
     1757    //  - write a file of ON events (ON2)
     1758    //    (after the g/h separation and after the final cuts)
     1759    //    (to be used for the flux calculation)
     1760    //
     1761    //  - do the flux calculation
     1762
     1763 if (Job6)
     1764 {
     1765    cout << "=====================================================" << endl;
     1766    cout << "Macro AnalyseCT1 : Start of Job 6" << endl;
     1767
     1768    cout << "" << endl;
     1769    cout << "Macro AnalyseCT1 : Job6, WON2 = "
     1770         << Job6  << ",  " << WON2 << endl;
     1771
     1772    TString typeInput = "ON";
     1773    cout << "typeInput = " << typeInput << endl;
     1774
     1775    TString typeMC   = "MC";
     1776    cout << "typeMC = " << typeMC << endl;
     1777
     1778    // use hadron matrix from ON or OFF data
     1779    TString typeData = "ON";
     1780    //TString typeData = "OFF";
     1781    cout << "typeData = " << typeData << endl;
     1782
     1783
     1784   //*************************************************************************
     1785   // read in matrices of look-alike events
     1786
     1787    const char* mtxName = "MatrixGammas";
     1788
     1789    TString outName = "matrix_gammas_";
     1790    outName += typeMC;
     1791    outName += "_";
     1792    outName += typeData;
     1793    outName += ".root";
     1794
     1795
     1796    cout << "" << endl;
     1797    cout << "========================================================" << endl;
     1798    cout << "Get matrix for (gammas)" << endl;
     1799    cout << "matrix name        = " << mtxName << endl;
     1800    cout << "name of root file  = " << outName << endl;
     1801    cout << "" << endl;
     1802
     1803
     1804    // read in the object with the name 'mtxName' from file 'outName'
     1805    //
     1806    TFile *fileg = new TFile(outName);
     1807
     1808    MHMatrix matrixg;
     1809    matrixg.Read(mtxName);
     1810    pmatrixg = &matrixg;
     1811
     1812    delete fileg;
     1813
     1814    //*****************************************************************
     1815
     1816    const char* mtxName = "MatrixHadrons";
     1817
     1818    TString outName = "matrix_hadrons_";
     1819    outName += typeData;
     1820    outName += ".root";
     1821
     1822
     1823    cout << "" << endl;
     1824    cout << "========================================================" << endl;
     1825    cout << " Get matrix for (hadrons)" << endl;
     1826    cout << "matrix name        = " << mtxName << endl;
     1827    cout << "name of root file  = " << outName << endl;
     1828    cout << "" << endl;
     1829
     1830
     1831    // read in the object with the name mtxName
     1832    //
     1833    TFile *fileh = new TFile(outName);
     1834
     1835    MHMatrix matrixh;
     1836    matrixh.Read(mtxName);
     1837    pmatrixh = &matrixh;
     1838
     1839    delete fileh;
     1840
     1841   //*************************************************************************
     1842
     1843
     1844    MTaskList tliston;
     1845    MParList pliston;
     1846    pliston.AddToList(&tliston);
     1847
     1848    pliston.AddToList(pmatrixg);
     1849    pliston.AddToList(pmatrixh);
     1850
     1851    //::::::::::::::::::::::::::::::::::::::::::
     1852
     1853    MBinning binssize("BinningSize");
     1854    binssize.SetEdgesLog(50, 10, 1.0e5);
     1855    pliston.AddToList(&binssize);
     1856
     1857    MBinning binsdistc("BinningDist");
     1858    binsdistc.SetEdges(50, 0, 1.4);
     1859    pliston.AddToList(&binsdistc);
     1860
     1861    MBinning binswidth("BinningWidth");
     1862    binswidth.SetEdges(50, 0, 1.0);
     1863    pliston.AddToList(&binswidth);
     1864
     1865    MBinning binslength("BinningLength");
     1866    binslength.SetEdges(50, 0, 1.0);
     1867    pliston.AddToList(&binslength);
     1868
     1869    MBinning binsalpha("BinningAlpha");
     1870    binsalpha.SetEdges(100, -100, 100);
     1871    pliston.AddToList(&binsalpha);
     1872
     1873    MBinning binsasym("BinningAsym");
     1874    binsasym.SetEdges(50, -1.5, 1.5);
     1875    pliston.AddToList(&binsasym);
     1876
     1877    MBinning binsht("BinningHeadTail");
     1878    binsht.SetEdges(50, -1.5, 1.5);
     1879    pliston.AddToList(&binsht);
     1880
     1881    MBinning binsm3l("BinningM3Long");
     1882    binsm3l.SetEdges(50, -1.5, 1.5);
     1883    pliston.AddToList(&binsm3l);
     1884
     1885    MBinning binsm3t("BinningM3Trans");
     1886    binsm3t.SetEdges(50, -1.5, 1.5);
     1887    pliston.AddToList(&binsm3t);
     1888
     1889   
     1890    //.....
     1891    MBinning binsb("BinningSigmabar");
     1892    binsb.SetEdges( 100,  0.0,  5.0);
     1893    pliston.AddToList(&binsb);
     1894
     1895    MBinning binth("BinningTheta");
     1896    binth.SetEdges( 70, -0.5, 69.5);   
     1897    pliston.AddToList(&binth);
     1898
     1899    MBinning binsdiff("BinningDiffsigma2");
     1900    binsdiff.SetEdges(100, -5.0, 20.0);
     1901    pliston.AddToList(&binsdiff);
     1902    //::::::::::::::::::::::::::::::::::::::::::
     1903
     1904
     1905    //-------------------------------------------
     1906    // create the tasks which should be executed
     1907    //
     1908
     1909    TString filenameData = typeInput;
     1910    filenameData += "1.root";
     1911
     1912    readname = "ReadCT1ONdata";
     1913    MReadMarsFile read("Events", filenameData);
     1914    read.DisableAutoScheme();
     1915
     1916
     1917    //.......................................................................
     1918    // g/h separation
     1919
     1920    MMultiDimDistCalc ghsep;
     1921    ghsep.SetUseNumRows(25);
     1922    ghsep.SetUseKernelMethod(kFALSE);
     1923
     1924
     1925
     1926    //.......................................................................
     1927
     1928    // geometry is needed in  MHHillas... classes
     1929    MGeomCam *fGeom =
     1930             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
     1931
     1932    TString fHilName    = "Hillas";
     1933    TString fHilSrcName = "HillasSrc";
     1934
     1935    Float_t hadcut = 0.2;
     1936    MSelFinal selfinalgh(fHilName, fHilSrcName);
     1937    selfinalgh.SetHadronnessCut(hadcut);
     1938    selfinalgh.SetAlphaCut(100.0);
     1939
     1940    MFillH fillh("MHHadronness");
     1941
     1942    MSelFinal selfinal(fHilName, fHilSrcName);
     1943    selfinal.SetHadronnessCut(hadcut);
     1944    selfinal.SetAlphaCut(20.0);
     1945
     1946    if (WON2)
     1947    {   
     1948      TString outNameImage = typeInput;
     1949      outNameImage += "_";
     1950      outNameImage += typeData;
     1951      outNameImage += "2.root";
     1952      MWriteRootFile write(outNameImage);
     1953      write.AddContainer("MHadronness",   "Events");
     1954      write.AddContainer("MSigmabar",     "Events");
     1955      write.AddContainer("Hillas",        "Events");
     1956      write.AddContainer("MMcEvt",        "Events");
     1957      write.AddContainer("HillasSrc",     "Events");
     1958      write.AddContainer("MRawRunHeader", "RunHeaders");
     1959      //write.AddContainer("MMcRunHeader","RunHeaders");
     1960      write.AddContainer("MSrcPosCam",    "RunHeaders");
     1961    }
     1962
     1963
     1964    MFillH hfill1("MHHillas",    fHilName);
     1965    MFillH hfill2("MHStarMap",   fHilName);
     1966
     1967    TString nxt = "HillasExt";
     1968    MHHillasExt hhilext(nxt, "", fHilName);
     1969    pliston.AddToList(&hhilext);
     1970    TString namext = nxt;
     1971    namext += "[MHHillasExt]";
     1972    MFillH hfill3(namext, fHilSrcName);
     1973
     1974    MFillH hfill4("MHHillasSrc", fHilSrcName);
     1975
     1976
     1977
     1978    //*****************************
     1979    // tasks to be executed
     1980   
     1981    tliston.AddToList(&read);
     1982
     1983    tliston.AddToList(&ghsep);
     1984    tliston.AddToList(&fillh);
     1985
     1986    tliston.AddToList(&selfinalgh);
     1987
     1988    tliston.AddToList(&hfill1);
     1989    tliston.AddToList(&hfill2);
     1990    tliston.AddToList(&hfill3);
     1991    tliston.AddToList(&hfill4);
     1992
     1993    tliston.AddToList(&selfinal);
     1994    if (WON2)
     1995      tliston.AddToList(&write);
     1996
     1997
     1998
     1999    //*****************************
     2000
     2001    //-------------------------------------------
     2002    // Execute event loop
     2003    //
     2004    MProgressBar bar;
     2005    MEvtLoop evtloop;
     2006    evtloop.SetParList(&pliston);
     2007    evtloop.SetProgressBar(&bar);
     2008
     2009    Int_t maxevents = 1000000000;
     2010    //Int_t maxevents = 1000;
     2011    if ( !evtloop.Eventloop(maxevents) )
     2012        return;
     2013
     2014    tliston.PrintStatistics(0, kTRUE);
     2015
     2016
     2017    //-------------------------------------------
     2018    // Display the histograms
     2019    //
     2020    TObject *c;
     2021
     2022    c = pliston.FindObject("MHHadronness")->DrawClone();
     2023    c->Print();
     2024
     2025    //c = pliston.FindObject("MHHillas")->DrawClone();
     2026
     2027    c = pliston.FindObject("MHHillasSrc")->DrawClone();
     2028
     2029    ////pliston.FindObject("MHHillasExt")->DrawClone();
     2030    //c = pliston.FindObject(nxt)->DrawClone();
     2031
     2032    c = pliston.FindObject("MHStarMap")->DrawClone();
     2033
     2034
     2035
     2036    cout << "Macro AnalyseCT1 : End of Job 6" << endl;
     2037    cout << "=======================================================" << endl;
     2038 }
     2039  //---------------------------------------------------------------------
    7182040}
    7192041
    … …  
    7212043
    7222044
     2045
     2046
     2047
     2048
     2049
     2050
Note: See TracChangeset for help on using the changeset viewer.