#define read_ntuple_cxx
#include "read_ntuple.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TRandom.h>

void read_ntuple::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L read_ntuple.C
//      Root > read_ntuple 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
   if (fChain == 0) return;
   ofstream  of;
   of.open("ccbar_d_mu.oscar");

   Long64_t nentries = fChain->GetEntriesFast();

   Int_t nbytes = 0, nb = 0;
   of<<"# OSC1999A"<<endl;
   of<<"# final_id_p_x"<<endl;
   of<<"# PYTHIA event generator ccbar->D->mu"<<endl;
   of<<"#"<<endl;

   float pt,theta;
   float deg_rad = 3.1415926535/180.;
   float min_ang = 6.*deg_rad;
   float max_ang = 42.*deg_rad;
 
   TRandom* tr = new TRandom();

   cout <<"Number of entries: "<<nentries<<endl;

   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      
      pt = sqrt(px_mu*px_mu+py_mu*py_mu);
      theta = atan(pt/fabs(pz_mu));

      if (theta>max_ang || theta<min_ang) continue;
      //cout <<"ientry - theta: "<<ientry<<" "<<theta<<endl;

      pt = sqrt(px_mu*px_mu+py_mu*py_mu);
      theta = atan(pt/fabs(pz_mu));
      if (theta>max_ang || theta<min_ang) continue;
      
      if (pz_parent*pz_mu <0 ) continue;

                        // take care of the mu vertex:
      vx_mu *= 1.0e13;    // vx_mu is in mm. However, the oscar output definition
      vy_mu *= 1.0e13;    // requires that it be written im fm
      vz_mu *= 1.0e13;    // 

      zvtx = -1000.;             // beam diamond in mm:
      while(fabs(zvtx)>400.) {
        zvtx = tr->Gaus(0.,220.);
      }
      //vz_mu = vz_mu + zvtx*1e13;    // add diamond to z

      of<<"0 1"<<endl;
      of<<"0 13 0 "<<px_mu<<" "<<py_mu<<" "<<pz_mu<<" "<<e_mu<<" 0 "<<vx_mu<<" "<<vy_mu<<" "<<vz_mu<<" 0"<<endl;
      of<<"0 0"<<endl;
      // if (Cut(ientry) < 0) continue;
   }
   of.close();
}