Changeset 1809 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
03/08/03 14:00:30 (22 years ago)
Author:
wittek
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
18 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r1800 r1809  
    11                                                 -*-*- END OF LINE -*-*-
     2
     3
     4 2003/03/08: Wolfgang Wittek
     5
     6    * macros/AnalyseCT1.C
     7      - for the CT1 data analysis
     8
     9    * mhist/MHMatrix.[h,cc]
     10      - let refcolumn start at 1 (not at 0)
     11
     12    * mhist/MHSigmaTheta.[h,cc]
     13      - Draw replaced by DrawCopy
     14      - add SetDirectory(NULL)
     15
     16    * manalysis/ MSelBasic.[h,cc]
     17                 MSelStandard.[h,cc]
     18                 MSelFinal.[h,cc]
     19      - more detailed output for errors
     20      - bugs removed
     21     
     22    * manalysis/ MPadSchweizer.[h,cc]
     23      - add SetDirectory(NULL)
     24      - add fErrors
     25
     26    * mfilter/ MFEventSelector.[h,cc]
     27      - add fErrors
     28
     29    * manalysis/MMultiDimDistCalc.[h,cc]
     30      - check division by zero
     31
     32    * mhist/MHHadronness.[h,cc]
     33      - check division by zero
     34      - normalize distributions of hadronness
     35
     36    * mfileio/MCT1ReadPreProc.[h,cc]
     37      - add event number (event.isecs_since_midday)
     38      - change definition of "fIsMcFile",
     39        because outpars.bmontecarlo is set wrongly sometimes
     40      - copy pedestalRMS for each event from the header information
     41      - check for the presence of a footer record even after reading
     42        a run header
     43
     44    * mmc/ MMcEvt.[hxx,cxx]]
     45      - add GetEvtNumber()
     46
     47
    248 2003/02/27: Abelardo Moralejo
    349
     
    63109    * mmc/MMcTrigHeader.hxx
    64110      - Added GetMultiplicity() and GetMeanThreshold()
    65 
    66111
    67112
  • 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
  • trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc

    r1660 r1809  
    208208    }
    209209
    210     //fHadronness->SetHadronness(dg/(dg+dh));
    211     fHadronness->SetHadronness(exp(-dh/dg));
     210    Double_t arg;
     211
     212    if (dg+dh != 0.0)
     213      arg = dg / (dg+dh);
     214    else
     215      arg = 1.e10;
     216    //fHadronness->SetHadronness(arg);
     217
     218    if (dg != 0.0)
     219      arg = exp(-dh/dg);
     220    else
     221      arg = 0.0;
     222    fHadronness->SetHadronness(arg);
     223     
    212224
    213225    return kTRUE;
  • trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc

    r1772 r1809  
    1 
    21/* ======================================================================== *\
    32!
     
    102101  fHDiffPixTheta  = fHist3Diff;
    103102
     103  fHSigmaTheta->SetDirectory(NULL);
     104  fHSigmaPixTheta->SetDirectory(NULL);
     105  fHDiffPixTheta->SetDirectory(NULL);
     106
    104107  Print();
    105108}
     
    239242   fHDiffPixTh->SetZTitle("Sigma^2-Sigmabar^2");
    240243
     244   //--------------------------------------------------------------------
     245
     246   memset(fErrors, 0, sizeof(fErrors));
     247
    241248   return kTRUE;
    242249}
     
    251258  //*fLog << "Entry MPadSchweizer::Process();" << endl;
    252259
     260  Int_t rc;
     261
    253262  const UInt_t npix = fEvt->GetNumPixels();
     263
     264  Double_t SigbarOld;
     265
     266  //*fLog << "before padding : " << endl;
     267  //SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);
     268  //fSigmabar->Print("");
     269
    254270
    255271  //$$$$$$$$$$$$$$$$$$$$$$$$$$
    256272  // to simulate the situation that before the padding the NSB and
    257273  // electronic noise are zero : set Sigma = 0 for all pixels
    258   for (UInt_t i=0; i<npix; i++)
    259   {   
    260     MCerPhotPix &pix = fEvt->operator[](i);     
    261     Int_t j = pix.GetPixId();
    262 
    263     MPedestalPix &ppix = fPed->operator[](j);
    264     ppix.SetMeanRms(0.0);
    265   }
     274  //for (UInt_t i=0; i<npix; i++)
     275  //{   
     276  //  MCerPhotPix &pix = fEvt->operator[](i);     
     277  //  Int_t j = pix.GetPixId();
     278
     279  //  MPedestalPix &ppix = fPed->operator[](j);
     280  //  ppix.SetMeanRms(0.0);
     281  //}
    266282  //$$$$$$$$$$$$$$$$$$$$$$$$$$
    267283
     
    269285  // Calculate average sigma of the event
    270286  //
    271   Double_t SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);
     287  SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);
    272288  //fSigmabar->Print("");
    273289
    274   //if (SigbarOld > 0.0)
    275   //{
     290  if (SigbarOld > 0.0)
     291  {
    276292    //*fLog << "MPadSchweizer::Process(); Sigmabar of event to be padded is > 0 : "
    277293    //      << SigbarOld << ". Stop event loop " << endl;
    278294    // input data should have Sigmabar = 0; stop event loop
    279     // return kFALSE;
    280   //}
     295 
     296    rc = 1;
     297    fErrors[rc]++;
     298    return kCONTINUE;
     299  }
    281300
    282301  Double_t Theta  = kRad2Deg*fMcEvt->GetTelescopeTheta();
     
    296315          << Theta << ",  " << binNumber << ";  Skip event " << endl;
    297316    // event cannot be padded; skip event
     317
     318    rc = 2;
     319    fErrors[rc]++;
    298320    return kCONTINUE;
    299321  }
     
    308330      // event cannot be padded; skip event
    309331      delete fHSigma;
     332
     333      rc = 3;
     334      fErrors[rc]++;
    310335      return kCONTINUE;
    311336    }
     
    331356    *fLog << "MPadSchweizer::Process(); target Sigmabar is less than SigbarOld : "
    332357          << Sigmabar << ",  " << SigbarOld << ",   Skip event" << endl;
     358
     359    rc = 4;
     360    fErrors[rc]++;
    333361    return kCONTINUE;     
    334362  }
     
    426454              << binTheta << "  and pixel bin " << binPixel 
    427455              << " has no entries;  aborting " << endl;
    428         return kERROR;
     456
     457        rc = 5;
     458        fErrors[rc]++;
     459        return kCONTINUE;     
    429460      }
    430461
     
    465496              << binTheta << "  and pixel bin " << binPixel 
    466497              << " has no entries;  aborting " << endl;
    467         return kERROR;
     498
     499        rc = 6;
     500        fErrors[rc]++;
     501        return kCONTINUE;     
    468502      }
    469503
     
    553587
    554588  // Calculate Sigmabar again and crosscheck
     589
    555590  Double_t SigbarNew = fSigmabar->Calc(*fCam, *fPed, *fEvt);
     591  //*fLog << "after padding : " << endl;
    556592  //fSigmabar->Print("");
    557593
     
    585621  //*fLog << "Exit MPadSchweizer::Process();" << endl;
    586622
     623  rc = 0;
     624  fErrors[rc]++;
     625
    587626  return kTRUE;
    588627
     
    594633Bool_t MPadSchweizer::PostProcess()
    595634{
     635    if (GetNumExecutions() != 0)
     636    {
     637
     638    *fLog << inf << endl;
     639    *fLog << GetDescriptor() << " execution statistics:" << endl;
     640    *fLog << dec << setfill(' ');
     641    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
     642          << (int)(fErrors[1]*100/GetNumExecutions())
     643          << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
     644
     645    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
     646          << (int)(fErrors[2]*100/GetNumExecutions())
     647          << "%) Evts skipped due to: Zenith angle out of range" << endl;
     648
     649    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
     650          << (int)(fErrors[3]*100/GetNumExecutions())
     651          << "%) Evts skipped due to: No data for generating Sigmabar" << endl;
     652
     653    *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
     654          << (int)(fErrors[4]*100/GetNumExecutions())
     655          << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;
     656
     657    *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3)
     658          << (int)(fErrors[5]*100/GetNumExecutions())
     659          << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
     660
     661    *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3)
     662          << (int)(fErrors[6]*100/GetNumExecutions())
     663          << "%) Evts skipped due to: No data for generating Sigma" << endl;
     664
     665    *fLog << " " << fErrors[0] << " ("
     666          << (int)(fErrors[0]*100/GetNumExecutions())
     667          << "%) Evts survived the padding!" << endl;
     668    *fLog << endl;
     669
     670    }
     671
     672    //---------------------------------------------------------------
    596673    TCanvas &c = *(MH::MakeDefCanvas("PadSchweizer", "", 900, 1200));
    597674    c.Divide(3, 4);
     
    600677
    601678    c.cd(1);
     679    fHSigmaTheta->SetDirectory(NULL);
    602680    fHSigmaTheta->SetTitle("2D : Sigmabar, \\Theta (reference sample)");
    603     fHSigmaTheta->DrawClone();     
     681    fHSigmaTheta->DrawCopy();     
    604682    fHSigmaTheta->SetBit(kCanDelete);     
    605683
    606684      //c.cd(1);
    607       //fHSigmaPixTheta->DrawClone();     
     685      //fHSigmaPixTheta->DrawCopy();     
    608686      //fHSigmaPixTheta->SetBit(kCanDelete);     
    609687
    610688    c.cd(4);
    611     fHSigbarTheta->DrawClone();     
     689    fHSigbarTheta->SetDirectory(NULL);
     690    fHSigbarTheta->DrawCopy();     
    612691    fHSigbarTheta->SetBit(kCanDelete);     
    613692
    614693
    615694    c.cd(7);
    616     fHNSB->DrawClone();
     695    fHNSB->SetDirectory(NULL);
     696    fHNSB->DrawCopy();
    617697    fHNSB->SetBit(kCanDelete);   
    618698
     
    625705    c.cd(2);
    626706    l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zx");
    627 
     707    l->SetDirectory(NULL);
    628708    l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
    629709    l->SetXTitle("\\Theta [\\circ]");
    630710    l->SetYTitle("Sigma^2-Sigmabar^2");
    631711
    632     *fLog << "before box" << endl;
    633 
    634     l->Draw("box");
     712    l->DrawCopy("box");
    635713    l->SetBit(kCanDelete);;
    636714
    637715    c.cd(5);
    638716    l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zy");
     717    l->SetDirectory(NULL);
    639718    l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
    640719    l->SetXTitle("pixel");
    641720    l->SetYTitle("Sigma^2-Sigmabar^2");
    642721
    643     l->Draw("box");
     722    l->DrawCopy("box");
    644723    l->SetBit(kCanDelete);;
    645724
     
    655734    c.cd(3);
    656735    k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zx");
     736    k->SetDirectory(NULL);
    657737    k->SetTitle("Sigma vs. \\Theta (all pixels)");
    658738    k->SetXTitle("\\Theta [\\circ]");
    659739    k->SetYTitle("Sigma");
    660740
    661     k->Draw("box");
     741    k->DrawCopy("box");
    662742    k->SetBit(kCanDelete);
    663743
    664744    c.cd(6);
    665745    k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zy");
     746    k->SetDirectory(NULL);
    666747    k->SetTitle("Sigma vs. pixel number (all \\Theta)");
    667748    k->SetXTitle("pixel");
    668749    k->SetYTitle("Sigma");
    669750
    670     k->Draw("box");
     751    k->DrawCopy("box");
    671752    k->SetBit(kCanDelete);;
    672753
     
    677758
    678759    c.cd(10);
    679     fHSigmaPedestal->DrawClone();
     760    fHSigmaPedestal->SetDirectory(NULL);
     761    fHSigmaPedestal->DrawCopy();
    680762    fHSigmaPedestal->SetBit(kCanDelete);   
    681763
    682764    c.cd(11);
    683     fHPhotons->DrawClone();
     765    fHPhotons->SetDirectory(NULL);
     766    fHPhotons->DrawCopy();
    684767    fHPhotons->SetBit(kCanDelete);   
    685768
  • trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h

    r1771 r1809  
    3939    Int_t          fRunType;
    4040    Int_t          fGroup;
     41
     42    Int_t          fErrors[7];
    4143
    4244    // plots used for the padding
  • trunk/MagicSoft/Mars/manalysis/MSelBasic.cc

    r1762 r1809  
    4444#include "MGeomCam.h"
    4545#include "MPedestalCam.h"
     46#include "MPedestalPix.h"
    4647#include "MGeomPix.h"
    4748
     
    119120Bool_t MSelBasic::Process()
    120121{
     122  /*
     123  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
     124  *fLog << "=========================================================" << endl;
     125  *fLog << "" << endl;
     126  *fLog << "MmcEvt data : " << endl;
     127  *fLog << "" << endl;
     128  fMcEvt->Print();
     129  *fLog << "" << endl;
     130  *fLog << "PartId() = " << fMcEvt->GetPartId() << endl;
     131  *fLog << "Energy() = " << fMcEvt->GetEnergy() << endl;
     132  *fLog << "Theta()  = " << fMcEvt->GetTheta() << endl;
     133
     134  *fLog << "Phi()    = " << fMcEvt->GetPhi() << endl;
     135  //*fLog << "CoreD()  = " << fMcEvt->GetCoreD() << endl;
     136  //*fLog << "CoreX()  = " << fMcEvt->GetCoreX() << endl;
     137
     138  //*fLog << "CoreY()  = " << fMcEvt->GetCoreY() << endl;
     139  *fLog << "Impact() = " << fMcEvt->GetImpact() << endl;
     140  //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl;
     141
     142  //*fLog << "PassPhotAtm() = " << fMcEvt->GetPassPhotAtm() << endl;
     143  //*fLog << "PassPhotRef() = " << fMcEvt->GetPassPhotRef() << endl;
     144  //*fLog << "PassPhotCone()  = " << fMcEvt->GetPassPhotCone() << endl;
     145
     146  //*fLog << "PhotElfromShower()    = " << fMcEvt->GetPhotElfromShower() << endl;
     147  //*fLog << "PhotElinCamera()  = " << fMcEvt->GetPhotElinCamera() << endl;
     148  *fLog << "TelescopePhi()  = " << fMcEvt->GetTelescopePhi() << endl;
     149
     150  *fLog << "TelescopeTheta()  = " << fMcEvt->GetTelescopeTheta() << endl;
     151  *fLog << "Impact() = " << fMcEvt->GetImpact() << endl;
     152  //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl;
     153  *fLog << "" << endl;
     154  *fLog << "=========================================================" << endl;
     155
     156  *fLog << "=========================================================" << endl;
     157  *fLog << "" << endl;
     158  *fLog << "MPedestalPix data : " << endl;
     159  *fLog << "" << endl;
     160
     161  Int_t ntot;
     162  ntot = fPed->GetSize();
     163    *fLog << "MeanRms() :" << endl;   
     164  for (Int_t i=0; i<ntot; i++)
     165  {
     166    MPedestalPix &pix = (*fPed)[i];
     167    //*fLog << "Mean()    = " << i << ",  " << pix.GetMean() << endl;   
     168    //*fLog << "Sigma()   = " << i << ",  " << pix.GetSigma() << endl;   
     169    *fLog << pix.GetMeanRms() << " ";   
     170    //*fLog << "SigmaRms()= " << i << ",  " << pix.GetSigmaRms() << endl;   
     171  }
     172  *fLog << "" << endl;
     173
     174  *fLog << "" << endl;
     175  ntot = fEvt->GetNumPixels();
     176    *fLog << "MCerPhotPix :  pix.GetNumPhotons()" << endl;   
     177  for (Int_t i=0; i<ntot; i++)
     178  {
     179    MCerPhotPix &pix = (*fEvt)[i];
     180    *fLog << pix.GetNumPhotons() << " ";   
     181  }
     182  *fLog << "" << endl;
     183
     184  *fLog << "" << endl;
     185  ntot = fEvt->GetNumPixels();
     186    *fLog << "MCerPhotPix :  pix.GetErrorPhot()" << endl;   
     187  for (Int_t i=0; i<ntot; i++)
     188  {
     189    MCerPhotPix &pix = (*fEvt)[i];
     190    *fLog << pix.GetErrorPhot() << " ";   
     191  }
     192  *fLog << "" << endl;
     193
     194  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
     195  */
     196
     197
    121198    Int_t rc = 0;
    122199
    123     //if ( fRawRun->GetRunNumber() < 1025 )
     200    //if ( fRawRun->GetRunNumber() < 16279 )
    124201    //{
    125202    //   rc = 1;
     
    128205
    129206    Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
    130     if (Theta > 45.0  ||  !SwTrigger() )
    131     {
    132       //*fLog << "MSelBasic::Process; Theta = " << Theta << endl;
     207    if ( Theta < 0.0 )
     208    {
     209      *fLog << "MSelBasic::Process; Run, Event, Theta = "
     210            << fRawRun->GetRunNumber()<< ",  "
     211            << fMcEvt->GetEvtNumber() << ",  " << Theta << endl;
    133212      rc = 1;
     213    }   
     214
     215    else if ( Theta > 45.0 )
     216    {
     217      rc = 2;
     218    }   
     219
     220    else if ( !SwTrigger() )
     221    {
     222      //*fLog << "MSelBasic::Process; SwTrigger = " << SwTrigger << endl;
     223      rc = 3;
    134224    }   
    135225
     
    223313    *fLog << GetDescriptor() << " execution statistics:" << endl;
    224314    *fLog << dec << setfill(' ');
    225     *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Basic selections are not fullfilled" << endl;
     315    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Zenith angle < 0" << endl;
     316
     317    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) << (int)(fErrors[2]*100/GetNumExecutions()) << "%) Evts skipped due to: Zenith angle too high" << endl;
     318
     319    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) << (int)(fErrors[3]*100/GetNumExecutions()) << "%) Evts skipped due to: Software trigger not fullfilled" << endl;
    226320
    227321    *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Basic selections!" << endl;
     
    230324    return kTRUE;
    231325}
     326
  • trunk/MagicSoft/Mars/manalysis/MSelBasic.h

    r1762 r1809  
    3030    const MRawRunHeader *fRawRun;       
    3131
    32           Int_t        fErrors[2];
     32          Int_t        fErrors[4];
    3333
    3434public:
  • trunk/MagicSoft/Mars/manalysis/MSelFinal.cc

    r1781 r1809  
    4141
    4242#include "MHillas.h"
     43#include "MHillasExt.h"
    4344#include "MHillasSrc.h"
    4445#include "MCerPhotEvt.h"
     
    5758// Default constructor.
    5859//
    59 MSelFinal::MSelFinal(MHillas *parhil, MHillasSrc *parhilsrc,
     60MSelFinal::MSelFinal(const char *HilName, const char *HilSrcName,
    6061                     const char *name, const char *title)
    6162{
     
    6364    fTitle = title ? title : "Task to evaluate the Final Cuts";
    6465
    65     fHil    = parhil;
    66     fHilsrc = parhilsrc;
     66    fHilName    = HilName;
     67    fHilSrcName = HilSrcName;
     68
     69    fHadronnessCut =   0.2;
     70    fAlphaCut      = 100.0; //degrees
    6771}
    6872
     
    7579Bool_t MSelFinal::PreProcess(MParList *pList)
    7680{
     81    fHil    = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt");
     82    if (!fHil)
     83    {
     84      *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl;
     85      return kFALSE;
     86    }
     87
     88    fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
     89    if (!fHilSrc)
     90    {
     91      *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl;
     92      return kFALSE;
     93    }
     94
     95
    7796    fHadronness = (MHadronness*)pList->FindObject("MHadronness");
    7897    if (!fHadronness)
     
    90109    }
    91110
    92     fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
    93     if (!fEvt)
    94     {
    95         *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
    96         return kFALSE;
    97     }
    98 
    99 
    100     fCam = (MGeomCam*)pList->FindObject("MGeomCam");
    101     if (!fCam)
    102     {
    103         *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
    104         return kFALSE;
    105     }
    106     fMm2Deg = fCam->GetConvMm2Deg();
    107 
    108     *fLog << "fMm2Deg = " << fMm2Deg << endl;
    109 
    110111    memset(fErrors, 0, sizeof(fErrors));
    111112
     
    122123Bool_t MSelFinal::Process()
    123124{
     125  //*fLog << "Entry MSelFinal; fHilSrc = " << fHilSrc << endl;
     126
     127   
     128
    124129    Int_t rc = 0;
    125130
    126     Double_t alphacut = 20.0;
     131    Double_t modalpha = fabs( fHilSrc->GetAlpha() );
    127132
    128     Double_t modalpha = fabs( fHilsrc->GetAlpha() );
    129133    Double_t h = fHadronness->GetHadronness();
    130134
    131     if ( h>0.5  ||  modalpha > alphacut )
     135    if ( h>fHadronnessCut )
    132136    {
    133137      //*fLog << "MSelFinal::Process; h, alpha = " << h << ",  "
    134       //      << fHilsrc->GetAlpha() << endl;
     138      //      << fHilSrc->GetAlpha() << endl;
    135139      rc = 1;
     140    }   
     141
     142    else if ( modalpha > fAlphaCut )
     143    {
     144      //*fLog << "MSelFinal::Process; h, alpha = " << h << ",  "
     145      //      << fHilSrc->GetAlpha() << endl;
     146      rc = 2;
    136147    }   
    137148
     
    153164    *fLog << GetDescriptor() << " execution statistics:" << endl;
    154165    *fLog << dec << setfill(' ');
    155     *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Final selections are not fullfilled" << endl;
     166    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
     167          << (int)(fErrors[1]*100/GetNumExecutions())
     168          << "%) Evts skipped due to: g/h separation cut (" << fHadronnessCut
     169          << ")" << endl;
    156170
    157     *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Final selections!" << endl;
     171    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
     172          << (int)(fErrors[2]*100/GetNumExecutions())
     173          << "%) Evts skipped due to: cut in ALPHA (" << fAlphaCut
     174          << " degrees)" << endl;
     175
     176    *fLog << " " << fErrors[0] << " ("
     177          << (int)(fErrors[0]*100/GetNumExecutions())
     178          << "%) Evts survived Final selections!" << endl;
    158179    *fLog << endl;
    159180
    160181    return kTRUE;
    161182}
     183
     184
     185
  • trunk/MagicSoft/Mars/manalysis/MSelFinal.h

    r1781 r1809  
    2828    MMcEvt      *fMcEvt;       
    2929    MHillas     *fHil;       
    30     MHillasSrc  *fHilsrc;       
     30    MHillasSrc  *fHilSrc;       
    3131    MHadronness *fHadronness;       
    3232
    3333    Double_t     fMm2Deg;   // conversion mm to degrees in camera
    34     Int_t        fErrors[2];
     34    Int_t        fErrors[3];
     35    TString      fHilName;
     36    TString      fHilSrcName;
     37 
     38    Float_t      fHadronnessCut;
     39    Float_t      fAlphaCut;
    3540
    3641public:
    37     MSelFinal(MHillas *fHil, MHillasSrc *fHilsrc,
     42    MSelFinal(const char *HilName, const char *HilSrcName,
    3843              const char *name=NULL, const char *title=NULL);
    3944
     
    4247    Bool_t PostProcess();
    4348
     49    void SetHadronnessCut(Float_t hadcut) { fHadronnessCut = hadcut; }
     50    void SetAlphaCut(Float_t alpha)       { fAlphaCut      = alpha;  }
    4451
    4552    ClassDef(MSelFinal, 0)   // Task to evaluate final cuts
  • trunk/MagicSoft/Mars/manalysis/MSelStandard.cc

    r1762 r1809  
     1
    12/* ======================================================================== *\
    23!
     
    3940
    4041#include "MHillas.h"
     42#include "MHillasExt.h"
    4143#include "MHillasSrc.h"
    4244#include "MCerPhotEvt.h"
     
    5456// Default constructor.
    5557//
    56 MSelStandard::MSelStandard(const MHillas *parhil, const MHillasSrc *parhilsrc,
     58MSelStandard::MSelStandard(const char *HilName, const char *HilSrcName,
    5759                           const char *name, const char *title)
    5860{
     
    6062    fTitle = title ? title : "Task to evaluate the Standard Cuts";
    6163
    62     fHil    = parhil;
    63     fHilsrc = parhilsrc;
     64    fHilName    = HilName;
     65    fHilSrcName = HilSrcName;
    6466}
    6567
     
    7274Bool_t MSelStandard::PreProcess(MParList *pList)
    7375{
     76    fHil    = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt");
     77    if (!fHil)
     78    {
     79      *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl;
     80      return kFALSE;
     81    }
     82
     83    fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
     84    if (!fHilSrc)
     85    {
     86      *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl;
     87      return kFALSE;
     88    }
     89
     90
    7491    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    7592    if (!fMcEvt)
     
    95112    fMm2Deg = fCam->GetConvMm2Deg();
    96113
    97     *fLog << "fMm2Deg = " << fMm2Deg << endl;
     114    //*fLog << "fMm2Deg = " << fMm2Deg << endl;
    98115
    99116    memset(fErrors, 0, sizeof(fErrors));
     
    114131    Int_t rc = 0;
    115132
    116     //Double_t fLength       = fHil->GetLength() * fMm2Deg;
    117     //Double_t fWidth        = fHil->GetWidth()  * fMm2Deg;
    118     Double_t fDist         = fHilsrc->GetDist()* fMm2Deg;
     133    Double_t fLength       = fHil->GetLength() * fMm2Deg;
     134    Double_t fWidth        = fHil->GetWidth()  * fMm2Deg;
     135    Double_t fDist         = fHilSrc->GetDist()* fMm2Deg;
    119136    //Double_t fDelta        = fHil->GetDelta()  * kRad2Deg;
    120137    Double_t fSize         = fHil->GetSize();
     
    122139    Int_t fNumCorePixels   = fHil->GetNumCorePixels();
    123140
    124     if (     fSize <= 60.0         ||  fDist< 0.4           ||  fDist > 1.1 
    125          ||  fNumUsedPixels >= 92  ||  fNumCorePixels < 4)
     141    if ( fNumUsedPixels >= 92  ||  fNumCorePixels < 4 )
    126142    {
    127143      //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = "
     
    129145      //      << fNumCorePixels << endl;
    130146      rc = 1;
     147    }   
     148
     149    else if ( fSize <= 60.0         ||  fDist< 0.4           ||  fDist > 1.1 )
     150    {
     151      //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = "
     152      //      << fSize << ",  " << fDist << ",  " << fNumUsedPixels << ",  "
     153      //      << fNumCorePixels << endl;
     154      rc = 2;
     155    }   
     156
     157    else if ( fLength <= 0.0         ||  fWidth <= 0.0 )
     158    {
     159      //*fLog << "MSelStandard::Process; fLength, fWidth = "
     160      //      << fLength << ",  " << fWidth << endl;
     161      rc = 3;
    131162    }   
    132163
     
    149180    *fLog << GetDescriptor() << " execution statistics:" << endl;
    150181    *fLog << dec << setfill(' ');
    151     *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Standard selections are not fullfilled" << endl;
     182    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Requirements on no.of used or core pxels not fullfilled" << endl;
     183
     184    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) << (int)(fErrors[2]*100/GetNumExecutions()) << "%) Evts skipped due to: Requirements on SIZE or DIST not fullfilled" << endl;
     185
     186    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) << (int)(fErrors[3]*100/GetNumExecutions()) << "%) Evts skipped due to: Length or Width is <= 0" << endl;
    152187
    153188    *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Standard selections!" << endl;
     
    156191    return kTRUE;
    157192}
     193
     194
  • trunk/MagicSoft/Mars/manalysis/MSelStandard.h

    r1762 r1809  
    2323{
    2424private:
    25     const MGeomCam    *fCam;      // Camera Geometry
    26     const MCerPhotEvt *fEvt;      // Cerenkov Photon Event
    27     const MMcEvt      *fMcEvt;       
    28     const MHillas     *fHil;       
    29     const MHillasSrc  *fHilsrc;       
     25    MGeomCam    *fCam;      // Camera Geometry
     26    MCerPhotEvt *fEvt;      // Cerenkov Photon Event
     27    MMcEvt      *fMcEvt;       
     28    MHillas     *fHil;       
     29    MHillasSrc  *fHilSrc;       
    3030
    31           Double_t     fMm2Deg;   // conversion mm to degrees in camera
    32           Int_t        fErrors[2];
     31    Double_t     fMm2Deg;   // conversion mm to degrees in camera
     32    Int_t        fErrors[4];
     33    TString      fHilName;
     34    TString      fHilSrcName;
    3335
    3436public:
    35     MSelStandard(const MHillas *fHil, const MHillasSrc *fHilsrc,
     37    MSelStandard(const char *HilName, const char *HilSrcName,
    3638                 const char *name=NULL, const char *title=NULL);
    3739
  • trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc

    r1749 r1809  
    122122{
    123123    const char *name = gSystem->ExpandPathName(txt);
     124
    124125    TString fname(name);
    125126    delete [] name;
     
    336337    // float   frms_pedsig_phot[iMAXNUMPIX];      // standard deviation of the calibrated signals from the pedestal run */
    337338    fPedest->InitSize(iMAXNUMPIX);
     339    *fLog << "PedestalRMS : ";
    338340    for (Int_t i=0; i<iMAXNUMPIX; i++)
     341    {
    339342        (*fPedest)[i].SetMeanRms(outpars.frms_pedsig_phot[i]);
     343        *fLog << outpars.frms_pedsig_phot[i] << " ";
     344 
     345        //$$$$$$$$$$$$$$$$$$$$$$$$$
     346        savePedRMS[i] =  outpars.frms_pedsig_phot[i];
     347        //$$$$$$$$$$$$$$$$$$$$$$$$$
     348    }
     349    *fLog << endl;
    340350
    341351    fPedest->SetReadyToSave();
     
    373383     *         be treated further. */
    374384
     385    *fLog << "outpars.bmontecarlo = " << outpars.bmontecarlo << endl;
     386    *fLog << "outpars.imcparticle = " << outpars.imcparticle << endl;
     387    *fLog << "outpars.dsourcera_hours = " << outpars.dsourcera_hours << endl;
     388    *fLog << "outpars.dsourcedec_deg = " << outpars.dsourcedec_deg << endl;
     389
     390    //*fLog << "File is a ";
     391    //*fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data");
     392    //*fLog << " file." << endl;
     393
     394
     395    // Next statement commented out because bmontecarlo was set wrongly
     396    //fIsMcFile = outpars.bmontecarlo==TRUE;
     397    fIsMcFile = (outpars.dsourcera_hours == 0.0 &&
     398                 outpars.dsourcedec_deg  == 0.0 &&
     399                 outpars.imcparticle != 0          );
     400
    375401    *fLog << "File is a ";
    376     *fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data");
     402    *fLog << (fIsMcFile ? "Monte Carlo" : "Real Data");
    377403    *fLog << " file." << endl;
    378 
    379     fIsMcFile = outpars.bmontecarlo==TRUE;
    380 
     404    *fLog << " " << endl;
    381405}
    382406
     
    433457    ProcessRunHeader(outpars);
    434458
     459
    435460    //rwagner: ReInit whenever new run commences
    436461    // rc==-1 means: ReInit didn't work out
     462
    437463    if (!fTaskList->ReInit(fParList))
    438464        return -1;
     
    453479    TString m = cEND_EVENTS_TEMPLATE;
    454480    Int_t p = m.First('%');
     481
    455482
    456483    if (!s.BeginsWith(m(0,p)))
     
    515542
    516543    //
    517     // Check for the existance of a next file to read
     544    // Check for the existence of a next file to read
    518545    //
    519546    TNamed *file = (TNamed*)fFileNames->First();
     
    538565
    539566    if (!CheckHeader(fname))
    540         return kFALSE;
     567    {
     568        *fLog << "OpenNextFile : CheckHeader(fname) is FALSE" << endl;
     569        return kFALSE;
     570    }
    541571
    542572    fIn = new ifstream(fname);
    543573
    544574    *fLog << inf << "-----------------------------------------------------------------------" << endl;
     575
    545576
    546577    switch (ReadRunHeader())
     
    553584        return kFALSE;
    554585      default:
     586        *fLog << "OpenNextFile : after ReadRunHeader, FnumEventsInRun = "
     587              << fNumEventsInRun << endl;
    555588        return kTRUE;
    556589      }
     590
     591
     592
    557593}
    558594
     
    752788void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event)
    753789{
     790  if (fRawRunHeader->GetRunNumber() == 1)
     791  {
     792  *fLog << "eventrecord" << endl;
     793  *fLog << "isecs_since_midday = " << event.isecs_since_midday << endl;
     794  *fLog << "isecfrac_200ns     = " << event.isecfrac_200ns << endl;
     795  *fLog << "snot_ok_flags      = " << event.snot_ok_flags << endl;
     796  *fLog << "ialt_arcs          = " << event.ialt_arcs << endl;
     797  *fLog << "iaz_arcs           = " << event.iaz_arcs << endl;
     798  *fLog << "ipreproc_alt_arcs  = " << event.ipreproc_alt_arcs << endl;
     799  *fLog << "ipreproc_az_arcs   = " << event.ipreproc_az_arcs << endl;
     800  *fLog << "ifieldrot_arcs     = " << event.ifieldrot_arcs << endl;
     801
     802  *fLog << "srate_millihz      = " << event.srate_millihz << endl;
     803  *fLog << "fhourangle         = " << event.fhourangle << endl;
     804  *fLog << "fmcenergy_tev      = " << event.fmcenergy_tev << endl;
     805  *fLog << "fmcsize_phel       = " << event.fmcsize_phel << endl;
     806  *fLog << "imcimpact_m        = " << event.imcimpact_m << endl;
     807  *fLog << "imcparticle        = " << event.imcparticle << endl;
     808  *fLog << "imctriggerflag     = " << event.imctriggerflag << endl;
     809  } 
     810
     811
     812    for (Int_t i=0; i<iMAXNUMPIX; i++)
     813    {
     814        (*fPedest)[i].SetMeanRms(savePedRMS[i]);
     815    }
     816
     817
    754818    //  int   isecs_since_midday; // seconds passed since midday before sunset (JD of run start)
    755819    //  int   isecfrac_200ns;     // fractional part of isecs_since_midday
     
    788852    // actual number of pixels (outputpars.inumpixels) is written out
    789853    // short spixsig_10thphot[iMAXNUMPIX];
     854    //*fLog << "spixsig_10thphot : " << endl;
    790855    for (Int_t i=0; i<iMAXNUMPIX; i++)
    791856    {
     857      //*fLog << event.spixsig_10thphot[i] << " ";
     858
    792859        if (event.spixsig_10thphot[i]==0)
    793860            continue;
     
    796863                         (*fPedest)[i].GetMeanRms());
    797864    }
     865    //*fLog << "" << endl;
     866
    798867    fNphot->SetReadyToSave();
    799868
     
    801870    // int ipreproc_az_arcs;  // "should be" az according to preproc (arcseconds)
    802871
    803     fMcEvt->Fill(0, /*fEvtNum*/
     872    fMcEvt->Fill(event.isecs_since_midday,     //0, /*fEvtNum*/
    804873                 fIsMcFile ? event.imcparticle : 0, /*corsika particle type*/
    805874                 fIsMcFile ? event.fmcenergy_tev*1000 : 0,
     
    867936        // an event
    868937        //
    869         if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0])
     938
     939      // "goto TONI" was introduced in order to check for a footer record
     940      // after reading a run header; this is necessary for runs with
     941      // zero events after the filter
     942      TONI:
     943
     944      if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0])
    870945            return kTRUE;
    871946
     
    874949        // must be an event
    875950        //
     951
    876952        switch (ReadRunFooter())
    877953        {
     
    906982        *fLog << "-----------------------------------------------------------------------" << endl;
    907983
     984
    908985        if (ReadRunHeader() < 0)
    909986        {
     
    911988            return kFALSE;
    912989        }
     990
     991        goto TONI;
     992
    913993        return kTRUE;
    914994    }
  • trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h

    r1735 r1809  
    4848    UInt_t fNumFilterEvts;  // number of events mentioned in the runs footers
    4949
     50    Float_t savePedRMS[127];
     51
     52
    5053    Bool_t OpenNextFile();
    5154
  • trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc

    r1762 r1809  
    118118Bool_t MFEventSelector::PreProcess(MParList *plist)
    119119{
     120    memset(fErrors, 0, sizeof(fErrors));
     121
    120122    fNumSelectedEvts = 0;
    121123    if (fSelRatio>0)
     
    153155Bool_t MFEventSelector::Process()
    154156{
     157    Int_t rc;
     158
    155159    const Float_t evt = gRandom->Uniform();
    156160
     
    161165
    162166    if (fResult)
     167    {
    163168        fNumSelectedEvts++;
    164169
     170        rc = 0;
     171        fErrors[rc]++;
     172        return kTRUE;
     173    }
     174   
     175    rc = 1;
     176    fErrors[rc]++;
     177
    165178    return kTRUE;
    166179}
     
    172185Bool_t MFEventSelector::PostProcess()
    173186{
     187    //---------------------------------
     188    if (GetNumExecutions() != 0)
     189    {
     190      *fLog << inf << endl;
     191      *fLog << GetDescriptor() << " execution statistics:" << endl;
     192      *fLog << dec << setfill(' ');
     193      *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
     194            << (int)(fErrors[1]*100/GetNumExecutions())
     195            << "%) Events not selected" << endl;
     196
     197      *fLog << " " << fErrors[0] << " ("
     198            << (int)(fErrors[0]*100/GetNumExecutions())
     199            << "%) Events selected!" << endl;
     200      *fLog << endl;
     201    }
     202
     203    //---------------------------------
    174204    if (TestBit(kNumTotalFromFile))
    175205        fNumTotalEvts = -1;
     
    199229
    200230}
    201 
  • trunk/MagicSoft/Mars/mfilter/MFEventSelector.h

    r1762 r1809  
    3737     */
    3838
     39    Int_t fErrors[2];
     40
    3941public:
    4042    // MFEventSelector();
  • trunk/MagicSoft/Mars/mhist/MHHadronness.cc

    r1668 r1809  
    235235    }
    236236
    237     return val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
     237    Float_t retValue;
     238    if (val2x-val1x != 0.0)
     239      retValue = val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
     240    else
     241      retValue = 0.0;
     242
     243    return retValue;
    238244}
    239245
     
    252258    fQfac->Set(n);
    253259
    254     const Stat_t sumg = fGhness->Integral(1, n+1);
    255     const Stat_t sump = fPhness->Integral(1, n+1);
     260    Stat_t sumg;
     261    Stat_t sump;
     262
     263    sumg = fGhness->Integral(1, n);
     264    sump = fPhness->Integral(1, n);
     265
     266    if (sumg == 0.0  ||  sump == 0.0)
     267    {
     268      *fLog << "MHHadronness::Finalize; sumg or sump is zero;   sumg, sump = "
     269            << sumg << ",  " << sump << ".  Cannot calculate hadronness"
     270            << endl;
     271    }
     272
     273
     274    // Normalize photon distribution
     275    Stat_t con;
     276    if (sumg > 0.0)
     277      for (Int_t i=1; i<=n; i++)
     278      {
     279        con = (fGhness->GetBinContent(i)) / sumg;
     280        fGhness->SetBinContent(i, con);       
     281      }
     282
     283    // Normalize hadron distribution
     284    if (sump > 0.0)
     285      for (Int_t i=1; i<=n; i++)
     286      {
     287        con = (fPhness->GetBinContent(i)) / sump;
     288        fPhness->SetBinContent(i, con);       
     289      }
     290
     291    // Calculate acceptances
     292    sumg = fGhness->Integral(1, n);
     293    sump = fPhness->Integral(1, n);
     294
     295    *fLog << "MHHadronness::Finalize; sumg, sump = " << sumg << ",  "
     296          << sump << endl;
    256297
    257298    Float_t max=0;
     
    259300    for (Int_t i=1; i<=n; i++)
    260301    {
    261         const Stat_t ip = fPhness->Integral(1, i)/sump;
    262         const Stat_t ig = fGhness->Integral(1, i)/sumg;
     302        Stat_t ip;
     303        if (sump != 0.0)
     304          ip = fPhness->Integral(1, i)/sump;
     305        else
     306          ip = 0;
     307
     308        Stat_t ig;
     309        if (sumg != 0.0)
     310          ig = fGhness->Integral(1, i)/sumg;
     311        else
     312          ig = 0;
    263313
    264314        fIntPhness->SetBinContent(i, ip);
     
    409459
    410460    c.cd(1);
    411     gStyle->SetOptStat(10);
     461    //gStyle->SetOptStat(10);
    412462    Getghness()->DrawCopy();
    413463    Getphness()->SetLineColor(kRed);
    414464    Getphness()->DrawCopy("same");
    415465    c.cd(2);
    416     gStyle->SetOptStat(0);
     466    //gStyle->SetOptStat(0);
    417467    Getighness()->DrawCopy();
    418468    Getiphness()->SetLineColor(kRed);
     
    483533
    484534    gPad->cd(1);
    485     gStyle->SetOptStat(10);
     535    //gStyle->SetOptStat(10);
    486536    Getghness()->Draw();
    487537    Getphness()->SetLineColor(kRed);
    488538    Getphness()->Draw("same");
    489539    gPad->cd(2);
    490     gStyle->SetOptStat(0);
     540    //gStyle->SetOptStat(0);
    491541    Getighness()->Draw();
    492542    Getiphness()->SetLineColor(kRed);
  • trunk/MagicSoft/Mars/mhist/MHMatrix.cc

    r1772 r1809  
    504504}
    505505
     506// --------------------------------------------------------------------------
     507//
    506508void MHMatrix::Reassign()
    507509{
     
    691693//
    692694// Define the reference matrix
    693 //   refcolumn  number of the column containing the variable, for which a
    694 //              target distribution may be given;
     695//   refcolumn  number of the column (starting at 1)containing the variable,
     696//              for which a target distribution may be given;
    695697//              if refcolumn is negative the target distribution will be set
    696698//              equal to the real distribution; the events in the reference
     
    748750   //---------------------------------------------------------
    749751   //
    750    // if refcol < 0 : select reference events randomly
    751    //                 i.e. set the normaliztion factotrs equal to 1.0
     752   // if refcolumn < 0 : select reference events randomly
     753   //                    i.e. set the normaliztion factotrs equal to 1.0
     754   // refcol is the column number starting at 0; it is >= 0
    752755
    753756   if (refcolumn<0)
    754757   {
    755      frefcol = -refcolumn;
     758     frefcol = -refcolumn - 1;
    756759   }
    757760   else
    758761   {
    759      frefcol = refcolumn;
     762     frefcol =  refcolumn - 1;
    760763   }
    761764   
     
    790793     //      << ",  " << fM(j-1,frefcol) << endl;
    791794   }
    792 
    793    // if refcolumn<0 set target distribution equal to the real distribution
    794    // in order to obtain normalization factors = 1.0
    795    //if (refcolumn<0)
    796    //{
    797    //  for (Int_t j=1; j<=fnbins; j++)
    798    //  {
    799    //    float cont = fHth-> GetBinContent(j);
    800    //    fHthsh->SetBinContent(j, cont);
    801    //  }
    802    //}
    803795
    804796   //---------------------------------------------------------
     
    816808
    817809     // if refcolumn < 0 set the correction factors equal to 1
    818      if ( refcolumn>=0.0 )
     810     if ( refcolumn>=0 )
    819811       b = fHthsh->GetBinContent(j);
    820812     else
     
    948940   {
    949941     *fLog <<ir <<" ";
    950      for (Int_t ic=0;ic<13;ic++) cout<<Mnew(ir,ic)<<" ";
     942     for (Int_t ic=0; ic<Mnew.GetNcols(); ic++) cout<<Mnew(ir,ic)<<" ";
    951943     *fLog <<endl;
    952944   }
     
    957949   {
    958950     float a = fHthaft->GetBinContent(j);
    959      if (a>0) *fLog << j << "  "<< a << "   ";
     951     *fLog << j << "  "<< a << "   ";
    960952   }
    961953   *fLog <<endl;
     
    967959
    968960   th1->cd(1);
    969    ((TH1F*)fHthsh)->Draw();      // target
     961   ((TH1F*)fHthsh)->DrawCopy();      // target
    970962
    971963   th1->cd(2);
    972    ((TH1F*)fHth)->Draw();        // real histogram before
     964   ((TH1F*)fHth)->DrawCopy();        // real histogram before
    973965
    974966   th1->cd(3);
    975    ((TH1F*)fHthd) -> Draw();     // correction factors
     967   ((TH1F*)fHthd)->DrawCopy();       // correction factors
    976968
    977969   th1->cd(4);
    978    ((TH1F*)fHthaft) -> Draw();   // histogram after
    979 
    980    //---------------------------------------------------------
    981    // --- write onto output file
    982    //
    983    //TFile *outfile = new TFile(fileNameout, "RECREATE", "");
    984    //Mnew.Write(fileNameout);
    985    //outfile->Close();
     970   ((TH1F*)fHthaft)->DrawCopy();     // histogram after
     971
     972   //---------------------------------------------------------
    986973
    987974   return kTRUE;
  • trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc

    r1785 r1809  
    224224TObject *MHSigmaTheta::DrawClone(Option_t *opt)
    225225{
    226     TCanvas &c = *MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta",
     226    TCanvas &c = *MakeDefCanvas("SigmaThetaPlot", "Sigmabar vs. Theta",
    227227                                 900, 900);
    228228    c.Divide(3, 3);
     
    236236    c.cd(1);
    237237    h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E");
     238    h->SetDirectory(NULL);
    238239    h->SetTitle("Distribution of \\Theta");
    239240    h->SetXTitle("\\Theta [\\circ]");
    240241    h->SetYTitle("No.of events");
    241242
    242     h->Draw(opt);
     243    h->DrawCopy(opt);
    243244    h->SetBit(kCanDelete);;
    244245    gPad->SetLogy();
     
    246247    c.cd(4);
    247248    h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E");
     249    h->SetDirectory(NULL);
    248250    h->SetTitle("Distribution of Sigmabar");
    249251    h->SetXTitle("Sigmabar");
    250252    h->SetYTitle("No.of events");
    251253
    252     h->Draw(opt);
     254    h->DrawCopy(opt);
    253255    h->SetBit(kCanDelete);;
    254256
     
    263265    c.cd(2);
    264266    l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zx");
     267    l->SetDirectory(NULL);
    265268    l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
    266269    l->SetXTitle("\\Theta [\\circ]");
    267270    l->SetYTitle("Sigma^2-Sigmabar^2");
    268271
    269     l->Draw("box");
     272    l->DrawCopy("box");
    270273    l->SetBit(kCanDelete);;
    271274
    272275    c.cd(5);
    273276    l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zy");
     277    l->SetDirectory(NULL);
    274278    l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
    275279    l->SetXTitle("pixel");
    276280    l->SetYTitle("Sigma^2-Sigmabar^2");
    277281
    278     l->Draw("box");
     282    l->DrawCopy("box");
    279283    l->SetBit(kCanDelete);;
    280284
     
    290294    c.cd(3);
    291295    k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zx");
     296    k->SetDirectory(NULL);
    292297    k->SetTitle("Sigma vs. \\Theta (all pixels)");
    293298    k->SetXTitle("\\Theta [\\circ]");
    294299    k->SetYTitle("Sigma");
    295300
    296     k->Draw("box");
     301    k->DrawCopy("box");
    297302    k->SetBit(kCanDelete);;
    298303
    299304    c.cd(6);
    300305    k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zy");
     306    k->SetDirectory(NULL);
    301307    k->SetTitle("Sigma vs. pixel number (all \\Theta)");
    302308    k->SetXTitle("pixel");
    303309    k->SetYTitle("Sigma");
    304310
    305     k->Draw("box");
     311    k->DrawCopy("box");
    306312    k->SetBit(kCanDelete);;
    307313
Note: See TracChangeset for help on using the changeset viewer.