source: trunk/MagicSoft/Mars/mtemp/mucm/macros/OptimizeSuperCuts.C@ 6315

/*############################################################################
  #
  #  How to use:
  #
  #   - 1. From a hillas file select a training subsample, for a given size
  #     range, with will be write into a MHMatrix. The the matrix will be read
  #     and the optimization will start.
  #     Once you've created the matrix you can set to false GenerateTrainMatrix
  #     and start directly the optimization from reading the matrix.
  #
  #   - 2. It is recomended to start optimazing the static cuts and then, the 
  #     dinanical cuts, this can be done:
  #      -seting in 'TString initParams' the name of the previous optimization
  #       results.
  #      -and using the the funtion MMyFindSuperCuts.SetOptimizationMode(Mode)
  #
  #     The supercuts function used is:
  #   
  #      Length_low = par[0]*S1 + par[1]*S2 +par[2]*ZA +par[3]*dist*dist
  #      Length_up  = par[4]*S1 + par[5]*S2 +par[6]*ZA +par[7]*dist*dist
  #                 ....
  #      and simiar with: Width, Dist,Asym, Conc, Leakage
  #
  #         where  S1 = log(size)-log(LowerSizeCut)
  #                S2 = S1*S1
  #                ZA = cos(theta/kRad2Deg) - 1;
  #
  #  NOTE:  Currently with neglect ZA depence, so we put ZA = 0;
  #  -----
  #
  #     The optimization modes are:
  #
  #        Mode   0 = optimize only supercuts params,i.e. 
  #                   only optimaze par0, par4, par8,..
  #               1 = optimize only dynamical params, i.e. 
  #                   all but par0, par4, par8
  #           other = optimize all the params at the same time.  
  #
  #
  #  Warning: Some times Minuit fails to converge, or its take too long. Then,
  #           it's better to try again with diferent nuumber of parameter,
  #           initial values or stepsizes.
  #
  #                                                 author: Marcos Lopez
  #                                                   email: marcos@gae.ucm.es
  ###########################################################################*/


Bool_t GenerateTrainMatrix = kTRUE;

void OptimizeSuperCuts()
{

    // --- Input -----------------------------------------

    Int_t  Mode = 0;
    Double_t SizeCutLow = 2000;
    Double_t SizeCutUp = 1e7;

    // Hillas file
    TString datafile = "~/MAGIC/data/Period022/2004_10_10/CrabNebula_20041010_6_30-25.hillas.root";
    
    // Intermediate file to speed up otpimization (from the original hillas 
    //file, a given number of events in a given size cut range will be stored 
    //in a MHatrix)
    TString matrixfile = "SC_matrix.150k.root";  // matrix with the 

    // Output file where the supercuts parameters wiil be saved
    TString outfile = "SC.150k.Above2000.root"; 


    // File with the starting supercuts parameters (if any)
    //
    //TString initParams="~/MAGIC/results/SuperCutsParams/AllDynCutsParams_Above2000ph.150k.CrabNebula_20041010_6_30-25.root";
    TString initParams="";

    // Number of train events to use from hillas file
    Int_t howmany = 150000;  

    MH3 h3("");

   
    //
    // fStepsizes 
    // if == 0.0    the parameter will be fixed in the minimization
    //    != 0.0    initial step sizes for the parameters
    //
    TArrayD fStepsizes(48);

    // LengthUp
    fStepsizes[0] = 0.1;
    fStepsizes[1] = 0.1;
    fStepsizes[2] = 0.1;
    fStepsizes[3] = 0.1;
    // LengthLo
    fStepsizes[4]  = 0.1;
    fStepsizes[5]  = 0.1;
    fStepsizes[6]  = 0.1;
    fStepsizes[7] = 0.1;
    // WidthUp
    fStepsizes[8] = 0.1;
    fStepsizes[9] = 0.1;
    fStepsizes[10] = 0.1;
    fStepsizes[11] = 0.1;
    // WidthLo
    fStepsizes[12] = 0.1;
    fStepsizes[13] = 0.1;
    fStepsizes[14] = 0.1;
    fStepsizes[15] = 0.1;
    // DistUp
    fStepsizes[16] = 0.1;
    fStepsizes[17] = 0.1;
    fStepsizes[18] = 0.1;
    fStepsizes[19] = 0.;
    // DistLo
    fStepsizes[20] = 0.1;
    fStepsizes[21] = 0.1;
    fStepsizes[22] = 0.1;
    fStepsizes[23] = 0.;

    for(int i=24;i<48;i++)
        fStepsizes[i] = .0;
 
    //------------------------------------------------------


    MMyFindSuperCuts a;
    a.SetHadronnessName("MHadronness");
    a.SetFilenameParam(outfile);
    a.SetSizeCuts(SizeCutLow, SizeCutUp);

    if(GenerateTrainMatrix)
        a.DefineTrainMatrix(datafile, h3 ,howmany, matrixfile); 


    a.ReadMatrix(matrixfile,"");
    a.SetOptimizationMode(Mode);


    //
    // Start optimization
    // 
    a.FindParams(initParams,0,fStepsizes);
}