// author: C. Delaere
// Exercise on Neural Networks.
// Note also the examples in root tutorials: root/tutorials/mlp

void NNsimple() 
{
  // we will look for 2 variables x,y
  // x,y are in [0,1]
  // the signal region is defined as (x>1/3 && y>2/3)||(x>2/3 && y>1/3)
  // the background region is defined as y < 4/3-x
  
  // Fill signal tree
  TTree* tree = new TTree("MonteCarlo", "Filtered Monte Carlo Events");
  Float_t x = -1;
  Float_t y = -1;
  Int_t signal = 1;
  tree->Branch("x",&x,"x/F");
  tree->Branch("y",&y,"y/F");
  tree->Branch("signal",&signal,"signal/I");
  // TODO: you way want to change the number of signal and background event that you generate
  for(int i=0;i<100;++i) {
    // one signal event ...
    x=-1; y=-1; signal=1;
    while(!((x>1/3. && y>2/3.)||(x>2/3. && y>1/3.))) {
      x = gRandom->Uniform();
      y = gRandom->Uniform();
    }
    tree->Fill();
  }
  for(int i=0;i<100;++i) {
    // one background event ...
    x=1; y=1; signal = 0;
    while(!(y < 4/3.-x)) {
      x = gRandom->Uniform();
      y = gRandom->Uniform();
    }
    tree->Fill();
  }

  // Build and train the network
  // TODO: here, you have to set the proper structure
  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron("x,y:XXXXXXXX:signal",tree);
  // TODO: change the number of iterations in the training (here 1000)
  mlp->Train(1000, "text,graph,update=10");

  // Result canvas
  TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas","Network analysis");
  mlpa_canvas->Divide(1,2);
  TMLPAnalyzer ana(mlp);
  ana.GatherInformations();
  ana.CheckNetwork();
  mlpa_canvas->cd(1);
  mlp->Draw();
  mlpa_canvas->cd(2);
  ana.DrawNetwork(0,"signal==1","signal==0"); // that draws the result

  // this is an example to show how to use the NN later on: we plot the NN in 2D
  Double_t vx[121];
  Double_t vy[121];
  Double_t  f[121];
  Double_t  v[2];
  // we simply run on a 2D grid and store x, y, NN(x,y).
  for (Int_t ix=0; ix<=10; ix++) {
     v[0]=ix/10.;
     for (Int_t iy=0; iy<=10; iy++) {
        v[1]=iy/10.;
        Int_t idx=ix*11+iy;
        vx[idx]=v[0];
        vy[idx]=v[1];
        f[idx]=mlp->Evaluate(0, v); // this is how we evaluate the NN at point v={x,y}...
     }
  }
  // now we can create a TGraph2D an draw it !
  TGraph2D* NN = new TGraph2D("NetworkOutput","Network Output",121, vx, vy, f);
  NN->Draw("cont4z");

  // TODO: now do the rejection/efficiency study and find the best working point !

}