Ignore:
Timestamp:
06/03/05 21:37:54 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp
Files:
85 deleted
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestDISP.C

    r6028 r7137  
    1 /* ======================================================================== *\
    2    !
    3    ! *
    4    ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
    5    ! * Software. It is distributed to you in the hope that it can be a useful
    6    ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
    7    ! * It is distributed WITHOUT ANY WARRANTY.
    8    ! *
    9    ! * Permission to use, copy, modify and distribute this software and its
    10    ! * documentation for any purpose is hereby granted without fee,
    11    ! * provided that the above copyright notice appear in all copies and
    12    ! * that both that copyright notice and this permission notice appear
    13    ! * in supporting documentation. It is provided "as is" without express
    14    ! * or implied warranty.
    15    ! *
    16    !
    17    !
    18    !   Author(s): Abelardo Moralejo 1/2005 <mailto:moralejo@pd.infn.it>!
    19    !
    20    !   Copyright: MAGIC Software Development, 2000-2003
    21    !
    22    !
    23    \* ======================================================================== */
    24 
    25 //////////////////////////////////////////////////////////////////////////////////
    26 //
    27 // macro RanForestDISP.C
    28 //
    29 // Version of Random Forest macro, intended to run on real data without assuming
    30 // any particular source position in the camera FOV. So we do not use the DIST
    31 // parameter, but rather (one version of) DISP method. We do not use either the
    32 // sign of the third longitudinal moment with respect to the center of the camera
    33 // (since the source can be somewhere else like in the case of wobble mode).
    34 //
    35 //////////////////////////////////////////////////////////////////////////////////
    36 
    37 
    38 #include "MMcEvt.hxx"
    39 
    40 void RanForestDISP(Int_t minsize=1000)
    41 {
    42     Int_t mincorepixels = 6; // minimum number of core pixels.
    43 //  Int_t mincorepixels = 0;
    44 
    45     TString skipevents;
    46     //  skipevents = "(MMcEvt.fImpact>15000.)";
    47     //  skipevents = "(MNewImagePar.fLeakage1>0.001)";
    48     //  skipevents = "(MHillas.fSize>1000)";
    49 
    50     // Skip calibration events, leaking events and flash events:
    51     skipevents = "(MNewImagePar.fLeakage1>0.1) || ({1.5+log10(MNewImagePar.fConc)*2.33333}-{6.5-log10(MHillas.fSize)}>0.)";
    52 
    53     //
    54     // Create a empty Parameter List and an empty Task List
    55     // The tasklist is identified in the eventloop by its name
    56     //
    57     MParList  plist;
    58 
    59     MTaskList tlist;
    60     plist.AddToList(&tlist);
    61 
    62 
    63     MReadMarsFile  read("Events", "/data1/magic/mc_data/root/Period021_0.73_mirror/FixedWindowPeakSearch/WOBBLE_Calib_6slices/tc_3.84_2.56/star_gamma_train_new.root");
    64     Float_t numgammas = read.GetEntries();
    65 
    66     read.AddFile("star_train_20050110_CrabNebulaW1.root");
    67     Float_t numhadrons = read.GetEntries() - numgammas;
    68 
    69     // Fraction of gammas to be used in training:
    70     //    Float_t gamma_frac = 1.5*(numhadrons / numgammas);
    71     Float_t gamma_frac = 1.;
    72 
    73     // Total hadron number to be processed in training:
    74     // If you run RF over different data runs, using always the same # of training gammas,
    75     // you should also use always a fixed number of hadrons for training (although the data
    76     // runs may be of different duration). Otherwise the meaning of the output hadronness
    77     // would be different for the different subsamples, since its absolute values depend on
    78     // the statistics of the training sample. The number set here is the number BEFORE the
    79     // "a priori cuts" set above through "skipevents".
    80 
    81     Float_t select_n_hadrons = 37000.; 
    82     Float_t hadron_frac = select_n_hadrons / numhadrons;
    83 
    84     read.DisableAutoScheme();
    85 
    86     tlist.AddToList(&read);
    87 
    88     MFParticleId fgamma("MMcEvt", '=', MMcEvt::kGAMMA);
    89     tlist.AddToList(&fgamma);
    90 
    91     MFParticleId fhadrons("MMcEvt", '!', MMcEvt::kGAMMA);
    92     tlist.AddToList(&fhadrons);
    93 
    94 
    95     MHMatrix matrixg("MatrixGammas");
    96 
    97     // TEST TEST TEST
    98 //     matrixg.AddColumn("MMcEvt.fImpact");
    99 //     matrixg.AddColumn("MMcEvt.fMuonCphFraction");
    100 
    101     //    matrixg.AddColumn("MNewImagePar.fLeakage1");
    102     //    matrixg.AddColumn("floor(0.5+(cos(MMcEvt.fTelescopeTheta)/0.01))");
    103 
    104     //    matrixg.AddColumn("MMcEvt.fTelescopePhi");
    105     //    matrixg.AddColumn("MSigmabar.fSigmabar");
    106 
    107 
    108     matrixg.AddColumn("MHillas.fSize");
    109 
    110     // We do not use DIST, since it depends on knowing the source position on camera.
    111     // In wobble mode the position changes with time. Besides, we intend to do a 2-d map
    112     // with the DISP method, and so we cannot use DIST in the cuts (since it would amde the
    113     // camera center a "provoleged" position).
    114 
    115     matrixg.AddColumn("(1.-(MHillas.fWidth/MHillas.fLength))*(0.9776+(0.101062*log10(MHillas.fSize)))");
    116 
    117 //     matrixg.AddColumn("MHillas.fWidth");
    118 //     matrixg.AddColumn("MHillas.fLength");
    119 
    120 
    121     // Scaling from real data?
    122     //    matrixg.AddColumn("MHillas.fWidth/(-8.37847+(9.58826*log10(MHillas.fSize)))");   
    123     //    matrixg.AddColumn("MHillas.fLength/(-40.1409+(29.3044*log10(MHillas.fSize)))");
    124 
    125     // Scaling from MC:
    126     matrixg.AddColumn("(MHillas.fWidth*3.37e-3)/(-0.028235+(0.03231*log10(MHillas.fSize)))");   
    127     matrixg.AddColumn("(MHillas.fLength*3.37e-3)/(-0.13527+(0.09876*log10(MHillas.fSize)))");
    128 
    129     matrixg.AddColumn("MNewImagePar.fConc");
    130 
    131     // TEST TEST TEST TEST
    132     //    matrixg.AddColumn("MConcentration.Conc7");
    133 
    134     //    matrixg.AddColumn("MNewImagePar.fConc1");
    135     //    matrixg.AddColumn("MNewImagePar.fAsym");
    136     //    matrixg.AddColumn("MHillasExt.fM3Trans");
    137     //    matrixg.AddColumn("abs(MHillasSrc.fHeadTail)");
    138 
    139 
    140     // We cannot use the 3rd moment with a sign referred to the nominal source position, if we
    141     // assume we do not know it a priori. We cannot use either the sign as referred to the
    142     // reconstructed source position through the DISP method, since that method uses precisely
    143     // fM3Long to solve the "head-tail" asimetry and hence all showers would be "gamma-like" from
    144     // that point of view. We add however the 3rd moment with no sign to the separation parameters.
    145 
    146     matrixg.AddColumn("abs(MHillasExt.fM3Long)");
    147 
    148     //    matrixg.AddColumn("MHillasSrc.fAlpha");
    149 
    150     TString dummy("(MHillas.fSize<");
    151     dummy += minsize;
    152 
    153     // AM TEST
    154     dummy += ") || (MNewImagePar.fNumCorePixels<";
    155     dummy += mincorepixels;
    156 
    157     // AM, useful FOR High Energy only:
    158     //    dummy += ") || (MImagePar.fNumSatPixelsLG>0";
    159     dummy += ")";
    160 
    161     MContinue SizeCut(dummy);
    162     tlist.AddToList(&SizeCut);
    163 
    164 
    165     if (skipevents != "")
    166     {
    167         MContinue SCut(skipevents);
    168         tlist.AddToList(&SCut);
    169     }
    170 
    171     plist.AddToList(&matrixg);
    172 
    173     MHMatrix matrixh("MatrixHadrons");
    174     matrixh.AddColumns(matrixg.GetColumns());
    175     plist.AddToList(&matrixh);
    176 
    177     MFEventSelector reduce_training_hadrons;
    178     reduce_training_hadrons.SetSelectionRatio(hadron_frac);
    179     MFilterList hadfilter;
    180     hadfilter.AddToList(&reduce_training_hadrons);
    181     hadfilter.AddToList(&fhadrons);
    182     tlist.AddToList(&hadfilter);
    183 
    184     MFillH fillmath("MatrixHadrons");
    185     fillmath.SetFilter(&hadfilter);
    186     tlist.AddToList(&fillmath);
    187 
    188     MContinue skiphad(&fhadrons);
    189     tlist.AddToList(&skiphad);
    190 
    191     MFEventSelector reduce_training_gammas;
    192     reduce_training_gammas.SetSelectionRatio(gamma_frac);
    193     tlist.AddToList(&reduce_training_gammas);
    194 
    195     MFillH fillmatg("MatrixGammas");
    196     fillmatg.SetFilter(&reduce_training_gammas);
    197     tlist.AddToList(&fillmatg);
    198 
    199     //
    200     // Create and setup the eventloop
    201     //
    202     MEvtLoop evtloop;
    203     evtloop.SetParList(&plist);
    204 
    205     MProgressBar *bar = new MProgressBar;
    206     evtloop.SetProgressBar(bar);
    207     //
    208     // Execute your analysis
    209     //
    210     if (!evtloop.Eventloop())
    211         return;
    212 
    213     tlist.PrintStatistics();
    214 
    215 
    216     // ---------------------------------------------------------------
    217     // ---------------------------------------------------------------
    218     // second event loop: the trees of the random forest are grown,
    219     // the event loop is now actually a tree loop (loop of training
    220     // process)
    221     // ---------------------------------------------------------------
    222     // ---------------------------------------------------------------
    223 
    224     MTaskList tlist2;
    225     plist.Replace(&tlist2);
    226 
    227     //
    228     //    matrixh.ShuffleRows(9999);
    229     //    matrixg.ShuffleRows(1111);
    230     //
    231 
    232     MRanForestGrow rfgrow2;
    233     rfgrow2.SetNumTrees(100);
    234     rfgrow2.SetNumTry(3);
    235     rfgrow2.SetNdSize(10);
    236 
    237     MWriteRootFile rfwrite2("RF.root");
    238     rfwrite2.AddContainer("MRanTree","Tree");       //save all trees
    239     MFillH fillh2("MHRanForestGini");
    240 
    241     tlist2.AddToList(&rfgrow2);
    242     tlist2.AddToList(&rfwrite2);
    243     tlist2.AddToList(&fillh2);
    244 
    245     // gRandom is accessed from MRanForest (-> bootstrap aggregating)
    246     // and MRanTree (-> random split selection) and should be initialized
    247     // here if you want to set a certain random number generator
    248     if(gRandom)
    249         delete gRandom;
    250     //    gRandom = new TRandom3(0); // Takes seed from the computer clock
    251     gRandom = new TRandom3(); // Uses always  same seed
    252     //
    253     // Execute tree growing (now the eventloop is actually a treeloop)
    254     //
    255 
    256     evtloop.SetProgressBar(bar);
    257     if (!evtloop.Eventloop())
    258         return;
    259 
    260     tlist2.PrintStatistics();
    261 
    262     plist.FindObject("MHRanForestGini")->DrawClone();
    263 
    264     // ---------------------------------------------------------------
    265     // ---------------------------------------------------------------
    266     // third event loop: the control sample (star2.root) is processed
    267     // through the previously grown random forest,
    268     //
    269     // the histograms MHHadronness (quality of g/h-separation) and
    270     // MHRanForest are created and displayed.
    271     // MHRanForest shows the convergence of the classification error
    272     // as function of the number of grown (and combined) trees
    273     // and tells the user how many trees are actually needed in later
    274     // classification tasks.
    275     // ---------------------------------------------------------------
    276     // ---------------------------------------------------------------
    277 
    278     MTaskList tlist3;
    279 
    280     plist.Replace(&tlist3);
    281 
    282 
    283     MReadMarsFile  read3("Events", "star_20050110_CrabNebulaW1.root");
    284     //    MReadMarsFile  read3("Events", "/data1/magic/mc_data/root/Period021_0.73_mirror/FixedWindowPeakSearch/WOBBLE_Calib_6slices/tc_3.84_2.56/star_gamma_test.root");
    285 
    286     read3.DisableAutoScheme();
    287     tlist3.AddToList(&read3);
    288     tlist3.AddToList(&SizeCut);
    289 
    290     if (skipevents != "")
    291         tlist3.AddToList(&SCut);
    292 
    293     MRanForestCalc calc;
    294     tlist3.AddToList(&calc);
    295 
    296 
    297     TString outfname = "star_S";
    298     outfname += minsize;
    299     outfname += "_20050110_CrabNebulaW1.root";
    300 
    301     //    outfname += "_mcgammas.root";
    302 
    303     // parameter containers you want to save
    304     //
    305 
    306     MWriteRootFile write(outfname, "recreate");
    307 
    308     //    write.AddContainer("MSigmabar", "Events");
    309 
    310     write.AddContainer("MRawEvtHeader",  "Events");
    311     write.AddContainer("MMcEvt",         "Events", kFALSE);
    312     write.AddContainer("MHillas",        "Events");
    313     write.AddContainer("MHillasExt",     "Events");
    314     write.AddContainer("MImagePar",      "Events");
    315     write.AddContainer("MNewImagePar",   "Events");
    316     write.AddContainer("MHillasSrc",     "Events");
    317     write.AddContainer("MHadronness",    "Events");
    318     write.AddContainer("MPointingPos",   "Events");
    319     write.AddContainer("MTime",          "Events", kFALSE);
    320     write.AddContainer("MConcentration", "Events", kFALSE);
    321 
    322     tlist3.AddToList(&write);
    323 
    324     evtloop.SetProgressBar(bar);
    325 
    326     //
    327     // Execute your analysis
    328     //
    329     if (!evtloop.Eventloop())
    330         return;
    331 
    332     tlist3.PrintStatistics();
    333 
    334     //    bar->DestroyWindow();
    335 
    336     return;
    337 }
     1matrixg.AddColumn("(1.-(MHillas.fWidth/MHillas.fLength))*(0.9776+(0.101062*log10(MHillas.fSize)))");
Note: See TracChangeset for help on using the changeset viewer.