void nndedx()
{
  gROOT->Reset(); 
  //  createinput("trackhisto_proton500_5000mev_mode0.root","trackhisto_pi500_5000mev_mode0.root","proton_vs_pi.root");
  //  nn("proton_vs_pi",3,4,0,500);
  // algorithm, neurons in 1st layer, neurons in 2nd layer, training cycles
  plotsb("proton_vs_pi",3,4,0,2800,0,1.5,0,0.8);
  //  plotsb("proton_vs_pi",3,4,0,2800,0,1.5,0.8,2.5);
  //  plotsb("proton_vs_pi",3,4,0,2800,1.5,10);

  //  createinput("trackhisto_e500_5000mev_mode0.root","trackhisto_pi500_5000mev_mode0.root","e_vs_pi.root");
  //  nn("e_vs_pi",3,4,0,500);
  // algorithm, neurons in 1st layer, neurons in 2nd layer, training cycles
  //  plotsb("e_vs_pi",3,4,0,3000,0,1.5,0,0.8);
  //  plotsb("e_vs_pi",3,4,0,3000,0,1.5,0.8,2.5);
  //  plotsb("e_vs_pi",3,4,0,3000,1.5,10);


  /* Algorithms:
  1) Stochastic
  2) Batch
  3) SteepestDescent
  4) RibierePolak
  5) FletcherReeves
  6) BFGS
  */
}


void createinput(Char_t input1[40],Char_t input2[40],Char_t output[40])
{
  cout << "I'm merging " << input1 << " and " << input2 << " into " << output << "." << endl;
  Int_t n;
  Double_t eta;
  Double_t p;
  Double_t dp;
  Double_t dr;
  Double_t dp2;
  Double_t harm2;

  Double_t n_copy;//note int->double conversion!
  Double_t eta_copy;
  Double_t p_copy;
  Double_t dpdr;
  Double_t dp2_copy;
  Double_t dedx;
  Double_t truth;

  TFile f1(input1);
  TTree *tree1 = (TTree*)gDirectory->Get("T");
  tree1->SetBranchAddress("n",&n);
  tree1->SetBranchAddress("eta",&eta);
  tree1->SetBranchAddress("p",&p);
  tree1->SetBranchAddress("dp",&dp);
  tree1->SetBranchAddress("dr",&dr);
  tree1->SetBranchAddress("dp2",&dp2);
  tree1->SetBranchAddress("harm2",&harm2);
  TFile f2(input2);
  TTree *tree2 = (TTree*)gDirectory->Get("T");
  tree2->SetBranchAddress("n",&n);
  tree2->SetBranchAddress("eta",&eta);
  tree2->SetBranchAddress("p",&p);
  tree2->SetBranchAddress("dp",&dp);
  tree2->SetBranchAddress("dr",&dr);
  tree2->SetBranchAddress("dp2",&dp2);
  tree2->SetBranchAddress("harm2",&harm2);

  TFile *f3 = new TFile(output,"RECREATE");
  TTree *nntree;
  nntree = new TTree("T","nn inputs");
  nntree->Branch("n",&n_copy,"n_copy/D");
  nntree->Branch("eta",&eta_copy,"eta_copy/D");
  nntree->Branch("p",&p_copy,"p_copy/D");
  nntree->Branch("dpdr",&dpdr,"dpdr/D");
  nntree->Branch("dp2",&dp2_copy,"dp2_copy/D");
  nntree->Branch("dedx",&dedx,"dedx/D");
  nntree->Branch("truth",&truth,"truth/D");

  Int_t minNumEntries=min(tree1->GetEntries(),tree2->GetEntries());

  // loop on file 1:
  truth=1.;
  cout << "Signal: " << tree1->GetEntries() << " entries." << endl;
  for (int i = 0; i < minNumEntries; i++) {
    tree1->GetEntry(i);
    vector<double> v;
    v.clear();
    v.push_back(1.*n);
    v.push_back(eta);
    v.push_back(p);
    v.push_back(dp/dr);
    v.push_back(dp2);
    v.push_back(harm2);
    vector<double> w=regularize(v);
    n_copy=w.at(0);
    eta_copy=w.at(1);
    p_copy=w.at(2);
    dpdr=w.at(3);
    dp2_copy=w.at(4);
    dedx=w.at(5);
    bool ok=true;
    for (int k=0;k<w.size();k++) {if (fabs(w.at(k))>10.) ok=false;} // exclude outliers!
    if (ok) nntree->Fill();
  } 
  // loop on file 2:
  truth=0.;
  cout << "Bkg: " << tree2->GetEntries() << " entries." << endl;
  for (int i = 0; i < minNumEntries; i++) {
    tree2->GetEntry(i);
    vector<double> v;
    v.clear();
    v.push_back(n);
    v.push_back(eta);
    v.push_back(p);
    v.push_back(dp/dr);
    v.push_back(dp2);
    v.push_back(harm2);
    vector<double> w=regularize(v);
    n_copy=w.at(0);
    eta_copy=w.at(1);
    p_copy=w.at(2);
    dpdr=w.at(3);
    dp2_copy=w.at(4);
    dedx=w.at(5);
    bool ok=true;
    for (int k=0;k<w.size();k++) {if (fabs(w.at(k))>10.) ok=false;} // exclude outliers!
    if (ok) nntree->Fill();
  } 

  nntree->Write();
  f3->Write();
  f3->Close();
  cout << "Done." << endl;
}

void nn(Char_t name[40], Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain)
{
  cout << "I'm training the NN on " << name << ".root" << endl;

  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
    cout << endl << "ROOT neural nets library has just been loaded." << endl << endl;
  }

  Double_t n;
  Double_t p;
  Double_t eta;
  Double_t dpdr;
  Double_t dp2;
  Double_t dedx;
  Double_t truth;

  Char_t rootfile[40];
  sprintf(rootfile,"%s.root",name);
  TFile f(rootfile);
  TTree *tree = (TTree*)gDirectory->Get("T");
  tree->SetBranchAddress("n",&n);
  tree->SetBranchAddress("p",&p);
  tree->SetBranchAddress("dpdr",&dpdr);
  tree->SetBranchAddress("dp2",&dp2);
  tree->SetBranchAddress("eta",&eta);
  tree->SetBranchAddress("dedx",&dedx);
  tree->SetBranchAddress("truth",&truth);

  // istructions: http://root.cern.ch/root/html310/TMultiLayerPerceptron.html
  TString structure("n,p,eta,dpdr,dp2,dedx"); // inputs
  structure+=":";
  structure+=neuro1;
  if (neuro2>0) { // two layers
    structure+=":";
    structure+=neuro2;
  }
  structure+=":";
  structure+="truth"; //target output; @ to normalize

  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(structure,tree);

  TEventList train;
  TEventList test;
  cout << tree->GetEntries() << " entries." << endl;
  for (int i = 0; i < tree->GetEntries(); i++) {
    tree->GetEntry(i);
    if (i%10 == 0) // only 1/10 of the available statistics
      train.Enter(i);
    else if (i%10 == 1) // only 1/10 of the available statistics
      test.Enter(i);
  }
  train.Print();
  test.Print();

  mlp->SetTrainingDataSet(&train);
  mlp->SetTestDataSet(&test);


  if (nAlgo == 1) mlp->SetLearningMethod(TMultiLayerPerceptron::kStochastic);
  if (nAlgo == 2) mlp->SetLearningMethod(TMultiLayerPerceptron::kBatch);
  if (nAlgo == 3) mlp->SetLearningMethod(TMultiLayerPerceptron::kSteepestDescent);
  if (nAlgo == 4) mlp->SetLearningMethod(TMultiLayerPerceptron::kRibierePolak);
  if (nAlgo == 5) mlp->SetLearningMethod(TMultiLayerPerceptron::kFletcherReeves);
  if (nAlgo == 6) mlp->SetLearningMethod(TMultiLayerPerceptron::kBFGS);
  if (nAlgo<1 || nAlgo>6) cout << "Algorithm unknown." << endl;
  TString algo;
  if (nAlgo == 1) algo="Stochastic";
  if (nAlgo == 2) algo="Batch";
  if (nAlgo == 3) algo="SteepestDescent";
  if (nAlgo == 4) algo="RibierePolak";
  if (nAlgo == 5) algo="FletcherReeves";
  if (nAlgo == 6) algo="BFGS";
  cout << "Algorithm: " << algo << endl;

  // training:
  mlp->Randomize();
  mlp->Train(ntrain, "text,graph,update=10");
  // loop:
  Int_t scelta=-1;
  while (scelta!=0) {
    cout << "How many other cycles?" << endl;
    cin >> scelta;
    if (scelta!=0) mlp->Train(scelta, "text,graph,update=10,+");
    ntrain = ntrain+scelta;
  }

  // Use TMLPAnalyzer to see what it looks for
  TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas",
				     "Network analysis");
  mlpa_canvas->SetWindowPosition(160,80);
  TMLPAnalyzer ana(mlp);
  // Initialisation
  ana.GatherInformations();
  // output to the console
  ana.CheckNetwork();
  mlp->Draw();

  // save weights:
  TString nomepesi("txt/");
  nomepesi+=name;
  nomepesi+="_";
  nomepesi+=algo;
  nomepesi+="_";
  nomepesi+=neuro1;
  nomepesi+="_";
  nomepesi+=neuro2;
  nomepesi+="_";
  nomepesi+=ntrain;
  nomepesi+=".txt";
  mlp->DumpWeights(nomepesi);
  cout << "File with weights: " << nomepesi << endl;


  cout << "Done." << endl;
  // http://root.cern.ch/root/html310/TMultiLayerPerceptron.html
}


void plotsb(Char_t name[40], Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain, float pmin=0, float pmax=1000, float etamin=0, float etamax=2.5)
{
  cout << "I'm using the NN to plot S/B for " << name << ".root" << endl;

  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
    cout << endl << "ROOT neural nets library has just been loaded." << endl << endl;
  }

  Double_t n;
  Double_t p;
  Double_t eta;
  Double_t dpdr;
  Double_t dp2;
  Double_t dedx;
  Double_t truth;

  Char_t rootfile[40];
  sprintf(rootfile,"%s.root",name);
  TFile f(rootfile);
  TTree *tree = (TTree*)gDirectory->Get("T");
  tree->SetBranchAddress("n",&n);
  tree->SetBranchAddress("p",&p);
  tree->SetBranchAddress("eta",&eta);
  tree->SetBranchAddress("dpdr",&dpdr);
  tree->SetBranchAddress("dp2",&dp2);
  tree->SetBranchAddress("dedx",&dedx);
  tree->SetBranchAddress("truth",&truth);

  // istructions: http://root.cern.ch/root/html310/TMultiLayerPerceptron.html
  TString structure("n,p,eta,dpdr,dp2,dedx"); // inputs
  structure+=":";
  structure+=neuro1;
  if (neuro2>0) { // two layers
    structure+=":";
    structure+=neuro2;
  }
  structure+=":";
  structure+="truth"; //target output; @ to normalize

  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(structure,tree);

  TString algo;
  if (nAlgo == 1) algo="Stochastic";
  if (nAlgo == 2) algo="Batch";
  if (nAlgo == 3) algo="SteepestDescent";
  if (nAlgo == 4) algo="RibierePolak";
  if (nAlgo == 5) algo="FletcherReeves";
  if (nAlgo == 6) algo="BFGS";
  cout << "Algorithm: " << algo << endl;

  TString nomepesi("txt/");
  nomepesi+=name;
  nomepesi+="_";
  nomepesi+=algo;
  nomepesi+="_";
  nomepesi+=neuro1;
  nomepesi+="_";
  nomepesi+=neuro2;
  nomepesi+="_";
  nomepesi+=ntrain;
  nomepesi+=".txt";
  mlp->LoadWeights(nomepesi);
  cout << "File with weights: " << nomepesi << endl;

  TH1F *h1, *i1;
  TH1F *h2, *i2;
  Int_t nbins=50;
  h1 = new TH1F("h1","",nbins,-.5,1.5);
  h2 = new TH1F("h2","",nbins,-.5,1.5);
  i1 = new TH1F("i1","",nbins,-1.1,1.1);
  i2 = new TH1F("i2","",nbins,-1.1,1.1);
  k1 = new TH1F("k1","",nbins,-1.1,1.1);
  k2 = new TH1F("k2","",nbins,-1.1,1.1);
  l1 = new TH1F("l1","",nbins,-1.1,1.1);
  l2 = new TH1F("l2","",nbins,-1.1,1.1);
  Double_t inputs[8];
  bool cut=true;
  for (int i = 0; i < tree->GetEntries(); i++) {
    tree->GetEntry(i);
    // apply cuts:
    vector<double> w;
    w.clear();
    w.push_back(n);
    w.push_back(eta);
    w.push_back(p);
    w.push_back(dpdr);
    w.push_back(dp2);
    w.push_back(dedx);
    vector<double> v=invert(w);
    cut = (v.at(2)>pmin && v.at(2)<pmax && fabs(v.at(1))>etamin && fabs(v.at(1))<etamax);
    // evaluate NN output:
    if (i%10 > 1 && cut) {
      inputs[0]=n;
      inputs[1]=p;
      inputs[2]=eta;
      inputs[3]=dpdr;
      inputs[4]=dp2;
      inputs[5]=dedx;
      if (truth==1) h1->Fill(mlp->Evaluate(0,inputs));
      if (truth==0) h2->Fill(mlp->Evaluate(0,inputs));
      if (truth==1) i1->Fill(dedx);
      if (truth==0) i2->Fill(dedx);
      if (truth==1) k1->Fill(dpdr);
      if (truth==0) k2->Fill(dpdr);
      if (truth==1) l1->Fill(dp2);
      if (truth==0) l2->Fill(dp2);
    }
  }

  gStyle->SetOptStat(0000);
  TCanvas* c = new TCanvas("c","Signal/Bkg comparison"); 
  c->SetWindowPosition(260,100);
  c->Divide(2,2);
  c->cd(1);
  Draw(i1,i2);
  Double_t d=fround(Distance(i1,i2),1);
  cout << "* Distance S-B (dedx):    " << d << endl;
  DrawLegend(i1,i2,"dE/dx",d);
  c->cd(2);
  Draw(k1,k2);
  Double_t d=fround(Distance(k1,k2),1);
  cout << "* Distance S-B (dp/dr):   " << d << endl;
  DrawLegend(i1,i2,"dp/dr",d);
  c->cd(3);
  Draw(l1,l2);
  Double_t d=fround(Distance(l1,l2),1);
  cout << "* Distance S-B (dp2):     " << d << endl;
  DrawLegend(i1,i2,"Pin-Pout",d);
  c->cd(4);
  Draw(h1,h2);
  Double_t d=fround(Distance(h1,h2),1);
  cout << "* Distance S-B (NN):      " << d << endl;
  DrawLegend(i1,i2,"NN",d);

  TString nomegif(name);
  nomegif+="_";
  nomegif+=algo;
  nomegif+="_";
  nomegif+=neuro1;
  nomegif+="_";
  nomegif+=neuro2;
  nomegif+="_";
  nomegif+=ntrain;
  if (pmin>0 || pmax<1000) {
    nomegif+="_p";
    nomegif+=(int)(pmin*1000);
    nomegif+="_";
    nomegif+=(int)(pmax*1000);
  }
  if (etamin>0 || etamax<2.5) {
    nomegif+="_eta";
    nomegif+=(int)(etamin*10);
    nomegif+="_";
    nomegif+=(int)(etamax*10);
  }
  nomegif+=".gif";
  c->Print(nomegif);
  cout << "File with plot: " << nomegif << endl;
 
  //  cout << "* Optimal cut on dedx: " << findOptimalCut(i1,i2,20.) << endl;
  //  cout << "* Optimal cut on NN:   " << findOptimalCut(h1,h2,20.) << endl;

  delete h1;
  delete h2;
  delete i1;
  delete i2;
  delete k1;
  delete k2;
  delete l1;
  delete l2;
  cout << "Done." << endl;
}

Double_t max(Double_t a, Double_t b)
{
  Double_t res=a;
  if (b>a) res=b;

  return res;
}

Double_t min(Double_t a, Double_t b)
{
  Double_t res=a;
  if (b<a) res=b;

  return res;
}

void Draw(TH1F *h1, TH1F *h2)
{
  h1->SetLineWidth(2);
  h2->SetLineWidth(2);
  h1->SetLineColor(kRed);
  h2->SetLineColor(kBlack);
  normalize(h1);
  normalize(h2);  
  if (h1->GetMaximum() < h2->GetMaximum()) {
    h2->DrawCopy(); 
    h1->DrawCopy("same");
  } else {
    h1->DrawCopy(); 
    h2->DrawCopy("same");
  }

}

void DrawLegend(TH1F *hist1, TH1F *hist2, Char_t var[10], Double_t d)
{
  TLegend *legend = new TLegend(0.0,0.60,0.35,0.95,var,"NDC");
  legend->AddEntry(hist1,"signal","l");
  legend->AddEntry(hist2,"background","l");
  legend->Draw();
  TPaveText *pave = new TPaveText(0.0,0.45,0.35,0.60,"NDC");
  TString txt("");
  txt+=d;
  txt+=" sigma";
  pave->SetTextAlign(12);
  pave->AddText(txt);
  pave->Draw();
}

/* round number n to d decimal points */
double fround(double n, unsigned d)
{
      return floor(n * pow(10., d) + .5) / pow(10., d);
}

void normalize(TH1F* h)
{
  Double_t scale;
  if (h->Integral()!=0) scale = 1/h->Integral();
  h->Scale(scale);
}

Double_t Distance(TH1F *h1, TH1F *h2)
{
  Double_t m1,m2,s1,s2;
  m1=h1->GetMean();
  m2=h2->GetMean();
  s1=h1->GetRMS();
  s2=h2->GetRMS();
  return (m1-m2)/sqrt(s1**2 + s2**2);
}

Double_t findOptimalCut(TH1F *hs, TH1F *hb, float B_over_S)
{
  Double_t cut=-1.;
  Double_t max=0.;
  bool quit=false;
  Int_t nbins = hs->GetNbinsX();
  for (int i=0; i<nbins; i++) {
    if (!quit) {
      Double_t value=hs->GetBinCenter(i);
      Double_t s=hs->Integral(i,nbins);
      Double_t b=B_over_S*(hb->Integral(i,nbins));
      Double_t q=0.;
      if (b>0.) {q=s/sqrt(b);} else {cut=value; quit=true;}
      if (q>max) {cut=value; max=q; }
    }
  }
  return cut;
}

vector<double> regularize(vector<double> v)
{
  vector<double> w;
  w.clear();
  double n=v.at(0);
  double eta=v.at(1);
  double p=v.at(2);
  double dpdr=v.at(3);
  double dp2=v.at(4);
  double dedx=v.at(5);
  w.push_back((1.*n-15.)/8.);
  w.push_back(eta/2.4);
  w.push_back((atan(p)-1.2)/.4);
  w.push_back((dpdr-0.007)/0.005);
  w.push_back(dp2/0.1);
  w.push_back(log(dedx)-1.);
  return w;
}

vector<double> invert(vector<double> v)
{
  vector<double> w;
  w.clear();
  double n=v.at(0);
  double eta=v.at(1);
  double p=v.at(2);
  double dpdr=v.at(3);
  double dp2=v.at(4);
  double dedx=v.at(5);
  w.push_back(8.*n+15.);
  w.push_back(eta*2.4);
  w.push_back(tan(0.4*p+1.2));
  w.push_back(dpdr*0.005+0.007);
  w.push_back(dp2*0.1);
  w.push_back(exp(dedx+1.));
  return w;
}