#define MyClass_cxx
#include "MyClass.h"
#include "TH2.h"
#include "TH1.h"
#include "TStyle.h"
#include "TCanvas.h"

void MyClass::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L MyClass.C
//      Root > MyClass t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch

 //  gROOT->Reset();

   if (fChain == 0) return;

   Int_t nentries = Int_t(fChain->GetEntriesFast());

   Int_t nbytes = 0, nb = 0;

   TH1F Neutrino_Energy = TH1F("Neutrino_Energy", "Neutrino Energy", 100,0,1);
   TH1F Neutrino_Mass = TH1F("Neutrino_Mass", "Neutrino Mass", 100,0,1);

   for (Int_t jentry=0; jentry<nentries;jentry++) {
   
      Int_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
      
      // NEUTRINO ENERGY
      Float_t Neutrino_E;
      Neutrino_E = 939.56563 - 938.272 - 0.510999 - ElectronEnergy - ProtonEnergy;

      // NEUTRINO MOMENTUM
      
      // Momentum Factor sqrt(T^2+ 2Tmass)
      Float_t Mom_Factor_e;
      Mom_Factor_e = sqrt( (ElectronEnergy)**2 + 2*ElectronEnergy* 0.510999 );

      Float_t Mom_Factor_p;
      Mom_Factor_p = sqrt( (ProtonEnergy)**2 + 2*ProtonEnergy*938.272 );
      
      // Electron and Proton position vector sizes
      Float_t mod_pos_e;
      mod_pos_e = (ElectronPosition_fX)**2 + (ElectronPosition_fY)**2 + (ElectronPosition_fZ)**2;
      mod_pos_e = sqrt(mod_pos_e);

      Float_t mod_pos_p;
      mod_pos_p = (ProtonPosition_fX)**2 + (ProtonPosition_fY)**2 + (ProtonPosition_fZ)**2;
      mod_pos_p = sqrt(mod_pos_p);

      // Electron Momentum Components
      Float_t e_Mom_x;
      e_Mom_x = (Mom_Factor_e *ElectronPosition_fX) / (mod_pos_e);

      Float_t e_Mom_y;
      e_Mom_y = (Mom_Factor_e *ElectronPosition_fY) / (mod_pos_e);

      Float_t e_Mom_z;
      e_Mom_z = (Mom_Factor_e *ElectronPosition_fZ) / (mod_pos_e);

      // Proton Momentum Components
      Float_t p_Mom_x;
      p_Mom_x = (Mom_Factor_p *ProtonPosition_fX) / (mod_pos_p);

      Float_t p_Mom_y;
      p_Mom_y = (Mom_Factor_p *ProtonPosition_fY) / (mod_pos_p);

      Float_t p_Mom_z;
      p_Mom_z = (Mom_Factor_p *ProtonPosition_fZ) / (mod_pos_p);

      // Neutrino Momentum Components
      Float_t nu_Mom_x;
      nu_Mom_x = -e_Mom_x - p_Mom_x;

      Float_t nu_Mom_y;
      nu_Mom_y = -e_Mom_y - p_Mom_y;

      Float_t nu_Mom_z;
      nu_Mom_z = -e_Mom_z - p_Mom_z;

      // Neutrino Momentum vector Size
      Float_t mod_mom_nu;
      mod_mom_nu = (nu_Mom_x)**2 + (nu_Mom_y)**2 + (nu_Mom_z)**2;
      
      // Neutrino Masss
      Float_t nu_Mass;
      nu_Mass = sqrt( (Neutrino_E)**2 - (mod_mom_nu) );

        //if (del) //GLOBAL CUTS FIRST       ????????????????

      Neutrino_Energy->Fill(Neutrino_E);
      Neutrino_Mass->Fill( nu_Mass);
  }

//++++++++++++++++++++++++++++

    Canvas = new TCanvas("Canvas","Plotting Canvas");
    Canvas->Divide(2);
    Canvas->cd(1);
    Neutrino_Energy->Draw();
    Canvas->cd(2);
    Neutrino_Mass->Draw();
    Canvas->Update();

    Int_t a;
   cin>>a;
//  if (Canvas) delete Canvas;

//----------------------------


}