void calpid(Int_t ntrain=40) {
   
  gROOT->SetStyle("Plain");
  gStyle->SetFillColor(10);
  if(!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
    gSystem->Load("libTreePlayer");
  }
  
  TFile input("../data/plus.root");
  TTree *tpp = (TTree *) input.Get("pp");
  TTree *tpip = (TTree *) input.Get("pip");

  TFile output("result.root", "RECREATE");
  TTree *simu = new TTree("simu", "P Pi Events");
  
  float pmom, eoverp, emcfac, emcvz, hacvz;
  int emccells, haccells;
  Int_t type;

  tpp->SetBranchAddress("pmom", &pmom);
  tpp->SetBranchAddress("eoverp", &eoverp);
  tpp->SetBranchAddress("emcfac", &emcfac);
  tpp->SetBranchAddress("emcvz", &emcvz);
  tpp->SetBranchAddress("hacvz", &hacvz);
  tpp->SetBranchAddress("emccells", &emccells);
  tpp->SetBranchAddress("haccells", &haccells);
  tpip->SetBranchAddress("pmom", &pmom);
  tpip->SetBranchAddress("eoverp", &eoverp);
  tpip->SetBranchAddress("emcfac", &emcfac);
  tpip->SetBranchAddress("emcvz", &emcvz);
  tpip->SetBranchAddress("hacvz", &hacvz);
  tpip->SetBranchAddress("emccells", &emccells);
  tpip->SetBranchAddress("haccells", &haccells);
 
  simu->Branch("pmom", &pmom, "pmom/F");
  simu->Branch("eoverp", &eoverp, "eoverp/F");
  simu->Branch("emcfac", &emcfac, "emcfac/F");
  simu->Branch("emcvz", &emcvz, "emcvz/F");
  simu->Branch("hacvz", &hacvz, "hacvz/F");
  simu->Branch("emccells", &emccells, "emccells/I");
  simu->Branch("haccells", &haccells, "haccells/I");
  simu->Branch("type", &type, "type/I");

  type = 1;
  Int_t i;
  for(i = 0; i < tpp->GetEntries(); i++) 
    {
      tpp->GetEntry(i);
      simu->Fill();
    }
  type = -1;
  for(i = 0; i < tpip->GetEntries(); i++) 
    {
      tpip->GetEntry(i);
      simu->Fill();
   }

  //save the tree for checking
  simu->Write();

   // Prepare event lists
   TEventList train;
   TEventList test;
   for (i = 0; i < simu->GetEntries(); i++) {
      if (i % 2)
         train.Enter(i);
      else
         test.Enter(i);
   }
   train.Print();
   test.Print();
   // Build and train the NN 
   TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron("pmom,eoverp,emcfac,emcvz,hacvz,emccells,haccells:10:type", simu, &train, &test);
   //   mlp->SetLearningMethod(TMultiLayerPerceptron::kStochastic);
   //   mlp->SetLearningMethod(TMultiLayerPerceptron::kBatch);
   //   mlp->SetLearningMethod(TMultiLayerPerceptron::kSteepestDescent);
   //   mlp->SetLearningMethod(TMultiLayerPerceptron::kRibierePolak); //crash
   //   mlp->SetLearningMethod(TMultiLayerPerceptron::kFletcherReeves); //crash
   //   mlp->SetLearningMethod(TMultiLayerPerceptron::kBFGS);
   mlp->Train(ntrain, "text,graph,update=10");
   // Use the NN to plot the results for each sample
   TH1F *hpp = new TH1F("hpp", "NN output", 500, -1.5, 1.5);
   TH1F *hpip = new TH1F("hpip", "NN output", 500, -1.5, 1.5);
   hpp->SetDirectory(0);
   hpip->SetDirectory(0);
   Double_t params[7];
   for (i = 0; i < tpip->GetEntries(); i++) 
     {
       tpip->GetEntry(i);
       params[0] = pmom;
       params[1] = eoverp;
       params[2] = emcfac;
       params[3] = emcvz;
       params[4] = hacvz;
       params[5] = emccells;
       params[6] = haccells;
       hpip->Fill(mlp->Evaluate(0, params));
     }
   for (i = 0; i < tpp->GetEntries(); i++) 
     {
       tpp->GetEntry(i);
       params[0] = pmom;
       params[1] = eoverp;
       params[2] = emcfac;
       params[3] = emcvz;
       params[4] = hacvz;
       params[5] = emccells;
       params[6] = haccells;
       hpp->Fill(mlp->Evaluate(0,params));
     }

   TCanvas *cv = new TCanvas("NNout_cv", "Neural net output");
   hpip->SetLineColor(kBlue);
   hpip->SetFillColor(kBlue);
   hpip->SetStats(0);

   hpp->SetLineColor(kRed);
   hpp->SetFillColor(kRed);
   hpp->SetFillStyle(3002);   
   hpp->SetStats(0);

   float resulthistomaximum = hpip->GetMaximum();
   if(resulthistomaximum < hpp->GetMaximum()) 
     resulthistomaximum = hpp->GetMaximum();
   
   hpip->SetMaximum(1.1 * resulthistomaximum);
   hpip->Draw();
   hpp->Draw("same");
   TLegend *legend = new TLegend(0.728,0.724576,0.876437,0.875,NULL,"brNDC");
   legend->SetBorderSize(0);
   legend->SetFillColor(10);
   legend->SetTextFont(62);
   legend->AddEntry(hpip, "#pi^{+}", "f");
   legend->AddEntry(hpp, "p", "f");
   legend->Draw();
   
   output.Close();
}