#include <cstdlib>
#include <vector>
#include <iostream>
#include <map>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TPluginManager.h"
#include "TStopwatch.h"
//
#include "TROOT.h"
#include "TRint.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TVector3.h"
#include "TMath.h"
#include "TFile.h"
#include "TSystem.h"
#include "TGStatusBar.h"
#include "TSystem.h"
#include "TXMLEngine.h"
#include "TTree.h"
#include "TLorentzVector.h"
#include "TNtuple.h"
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TGraph.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TList.h"
#include "TMatrixT.h"
#include "TMatrixD.h"
#include "TLatex.h"

#include <cmath>
#include <iostream>
#include <string>
#include <cstdlib>
#include<iomanip>
#include <fstream>
#include <map>
#include <sstream>
#include <stdexcept>
#include <vector>

//#include "TMVAGui.C"

#if not defined(__CINT__) || defined(__MAKECINT__)
#include "TMVA/Tools.h"
#include "TMVA/Reader.h"
#include "TMVA/MethodCuts.h"
#endif

using namespace std;
#define SQR(a) ((a)*(a))
#define PI 3.14159
ofstream qin;

TDirectoryFile *FitDir[500];
void triggerstudy2()
{   
Float_t  p_tof, pim_tof, pip_tof, beamE;
Float_t p_sec, pip_sec, pim_sec;
Float_t bitcon; 
//Int_t trigbit;
Float_t trigbit;


  int prefix;
  //Int_t NRun = 0;
  Float_t event;
  
  Float_t kin_beamE;
  Float_t mm_m;

  Float_t pip_pim_mm;

  //things not in gamp2root

  Float_t pip_pim, run;
  Float_t p_pip_pim;

  Float_t pip_mm;
  Float_t pim_mm;
  Float_t p_pip;
  Float_t p_pim;
  Float_t p_mm;
  Float_t p_pip_mm;
  Float_t p_pim_mm;

  Float_t z;
  Float_t x;
  Float_t y;

  Float_t mmsq;

  Float_t betapim;
  Float_t betapim_tof;
  Float_t betapip,betaproton;
  Float_t betapip_tof,betaproton_tof;
  Float_t dbetapim;
  Float_t dbetapip, dbetap;

  Float_t vtd;

  Float_t p_ngrf;
  Float_t pip_ngrf;
  Float_t pim_ngrf;
  Float_t p_tagr_id;
  Float_t pip_tagr_id;
  Float_t pim_tagr_id;

  Float_t mm_E;


  Float_t pip_px;
  Float_t pip_py;
  Float_t pip_pz;
  Float_t pim_px;
  Float_t pim_py;
  Float_t pim_pz;
  Float_t kin_pip_px;
  Float_t kin_pip_py;
  Float_t kin_pip_pz;
  Float_t kin_pim_px;
  Float_t kin_pim_py;
  Float_t kin_pim_pz;
  Float_t mm_px;
  Float_t mm_py;
  Float_t mm_pz;
  Float_t p_px;
  Float_t p_py;
  Float_t p_pz;
  Float_t kin_p_px;
  Float_t kin_p_py;
  Float_t kin_p_pz;
  Float_t pip_E;
  Float_t pim_E;
  Float_t p_E;
  Float_t kin_pip_E;
  Float_t kin_pim_E;
  Float_t kin_p_E;
  
  /*Float_t pip_p;
  Float_t pim_p;
  Float_t p_p;*/

  Float_t theta_mesons;
  Float_t t;
  Float_t p_m;

  Float_t beta_cm;
  Float_t gamma_cm,gamma_cm2;
  Float_t x_pz,x_px,x_py;
  Float_t x_pz_cm,p0_E_cm,p0_pz_cm,beamE_cm,p_pz_cm,p_E_cm;
  Float_t x_E;
  Float_t x_pt,theta_cm,cos_theta_cm;
  Int_t pi0_fail,mm_E_fail,dbeta_fail,vtd_fail;

  Float_t bec_pi0_prob;
  Float_t bec_nop_prob;

Float_t p_stvtime;
Float_t pip_stvtime;
Float_t pim_stvtime;

Float_t p_sc_id_part;
Float_t pip_sc_id_part;
Float_t pim_sc_id_part;

Float_t d_pim1;
Float_t d_pim2;


char name1[500];

float phi_p =0;
//float phi_rad_p=0;


float phi_pip=0;
//float phi_rad_pip=0;


float phi_pim=0;
//float phi_rad_pim=0;


TLorentzVector MyProton, MyPip, MyPim, MyProton_Rest, My_Initial, MyPi0, MyPhoton, MM;
         TVector3 omega_pvec;
         Float_t m_p = 0.938;
TFile f ("triggermap_may2017v2.root", "RECREATE", "Histograms from ntuples" );   
TH2D *plot = new TH2D("plot","plot",57,0,57,180,0,360);
  
TH2D *proton_hit[6];
for (int i=0; i<6; i++)
{
 int sector = 1+i;
 sprintf(name1,"hit_proton_sector_%i",sector);
 proton_hit[i] = new TH2D(name1,name1,57,0,57,30,-30,30);
 proton_hit[i]->Sumw2();
}

TH2D *pip_hit[6];
for (int i=0; i<6; i++)
{
 int sector = 1+i;
 sprintf(name1,"hit_pip_sector_%i",sector);
 pip_hit[i] = new TH2D(name1,name1,57,0,57,30,-30,30);
 pip_hit[i]->Sumw2();
}

TH2D *pim_hit[6];
for (int i=0; i<6; i++)
{
 int sector = 1+i;
 sprintf(name1,"hit_pim_sector_%i",sector);
 pim_hit[i] = new TH2D(name1,name1,57,0,57,30,-30,30);
 pim_hit[i]->Sumw2();
}

TH2D *proton_total[6];
for (int i=0; i<6; i++)
{
 int sector = 1+i;
 sprintf(name1,"total_proton_sector_%i",sector);
 proton_total[i] = new TH2D(name1,name1,57,0,57,30,-30,30);
 proton_total[i]->Sumw2();
}

TH2D *pip_total[6];
for (int i=0; i<6; i++)
{
 int sector = 1+i;
 sprintf(name1,"total_pip_sector_%i",sector);
 pip_total[i] = new TH2D(name1,name1,57,0,57,30,-30,30);
 pip_total[i]->Sumw2();
}

TH2D *pim_total[6];
for (int i=0; i<6; i++)
{
 int sector = 1+i;
 sprintf(name1,"total_pim_sector_%i",sector);
 pim_total[i] = new TH2D(name1,name1,57,0,57,30,-30,30);
 pim_total[i]->Sumw2();
}


TChain *ttree;
ttree = new TChain("ntp");

ttree->Add("/d/grid6/akbar/g12/ROOT-new/5652*.*.ppm.bos.root");
ttree->Add("/d/grid6/akbar/g12/ROOT-new/5653*.*.ppm.bos.root");
ttree->Add("/d/grid6/akbar/g12/ROOT-new/5654*.*.ppm.bos.root");
/*ttree->Add("/d/grid6/akbar/g12/ROOT-new/5655*.*.ppm.bos.root");
ttree->Add("/d/grid6/akbar/g12/ROOT-new/5656*.*.ppm.bos.root");

ttree->Add("/d/grid6/akbar/g12/ROOT-new/5657*.*.ppm.bos.root");
ttree->Add("/d/grid6/akbar/g12/ROOT-new/5658*.*.ppm.bos.root");
ttree->Add("/d/grid6/akbar/g12/ROOT-new/5659*.*.ppm.bos.root");*/



    


TTree *tree = (TTree*)gROOT->FindObject("ntp");

tree->SetBranchAddress( "p_sec", &p_sec);
tree->SetBranchAddress( "pip_sec", &pip_sec);
tree->SetBranchAddress( "pim_sec", &pim_sec);
tree->SetBranchAddress( "p_tof", &p_tof);
tree->SetBranchAddress( "pim_tof", &pim_tof);
tree->SetBranchAddress( "pip_tof", &pip_tof);
tree->SetBranchAddress( "beamE", &beamE);
tree->SetBranchAddress( "trigbit", &trigbit);
tree->SetBranchAddress("mm_m",&mm_m);
  tree->SetBranchAddress("mm_E",&mm_E);
  tree->SetBranchAddress("beamE",&beamE);


  tree->SetBranchAddress("pip_pim_mm",&pip_pim_mm);

  tree->SetBranchAddress("runNumber",&run);
  tree->SetBranchAddress("eventNumber",&event);
  tree->SetBranchAddress("pip_pim",&pip_pim);
  tree->SetBranchAddress("p_pip_pim",&p_pip_pim);

  tree->SetBranchAddress("pip_mm",&pip_mm);
  tree->SetBranchAddress("pim_mm",&pim_mm);
  tree->SetBranchAddress("p_pip_mm",&p_pip_mm);
  tree->SetBranchAddress("p_pim_mm",&p_pim_mm);
  tree->SetBranchAddress("p_pip",&p_pip);
  tree->SetBranchAddress("p_pim",&p_pim);
  tree->SetBranchAddress("p_mm",&p_mm);

  tree->SetBranchAddress("betapim",&betapim);
  tree->SetBranchAddress("betapim_tof",&betapim_tof);
  tree->SetBranchAddress("betapip",&betapip);
  tree->SetBranchAddress("betapip_tof",&betapip_tof);
  tree->SetBranchAddress("betaproton",&betaproton);
  tree->SetBranchAddress("betaproton_tof",&betaproton_tof);

  tree->SetBranchAddress("mmsqorig",&mmsq);


  tree->SetBranchAddress("z",&z);
  tree->SetBranchAddress("x",&x);
  tree->SetBranchAddress("y",&y);

  tree->SetBranchAddress("stVtimeDiff",&vtd);
  tree->SetBranchAddress("kin_pip_pz",&kin_pip_pz);
  tree->SetBranchAddress("kin_pip_px",&kin_pip_px);
  tree->SetBranchAddress("kin_pip_py",&kin_pip_py);
  tree->SetBranchAddress("kin_pim_pz",&kin_pim_pz);
  tree->SetBranchAddress("kin_pim_px",&kin_pim_px);
  tree->SetBranchAddress("kin_pim_py",&kin_pim_py);
  tree->SetBranchAddress("kin_p_pz",&kin_p_pz);
  tree->SetBranchAddress("kin_p_px",&kin_p_px);
  tree->SetBranchAddress("kin_p_py",&kin_p_py);
  tree->SetBranchAddress("kin_pip_E",&kin_pip_E);
  tree->SetBranchAddress("kin_pim_E",&kin_pim_E);
  tree->SetBranchAddress("kin_p_E",&kin_p_E);
  tree->SetBranchAddress("kin_beamE",&kin_beamE);
  tree->SetBranchAddress("p_ngrf",&p_ngrf);
  tree->SetBranchAddress("pip_ngrf",&pip_ngrf);
  tree->SetBranchAddress("pim_ngrf",&pim_ngrf);
  tree->SetBranchAddress("p_tagr_id",&p_tagr_id);
  tree->SetBranchAddress("pip_tagr_id",&pip_tagr_id);
  tree->SetBranchAddress("pim_tagr_id",&pim_tagr_id);

   tree->SetBranchAddress("p_m",&p_m);


  tree->SetBranchAddress("pip_pz",&pip_pz);
  tree->SetBranchAddress("pip_px",&pip_px);
  tree->SetBranchAddress("pip_py",&pip_py);
  tree->SetBranchAddress("pim_pz",&pim_pz);
  tree->SetBranchAddress("pim_px",&pim_px);
  tree->SetBranchAddress("pim_py",&pim_py);


  tree->SetBranchAddress("mm_pz",&mm_pz);
  tree->SetBranchAddress("mm_px",&mm_px);
  tree->SetBranchAddress("mm_py",&mm_py);
  tree->SetBranchAddress("p_pz",&p_pz);
  tree->SetBranchAddress("p_px",&p_px);
  tree->SetBranchAddress("p_py",&p_py);

  tree->SetBranchAddress("pip_E",&pip_E);
  tree->SetBranchAddress("pim_E",&pim_E);
  tree->SetBranchAddress("p_E",&p_E);

  
    tree->SetBranchAddress("bec_pi0_prob",&bec_pi0_prob);
    tree->SetBranchAddress("bec_nop_prob",&bec_nop_prob);


tree->SetBranchAddress("pim_stvtime",&pim_stvtime);
tree->SetBranchAddress("pip_stvtime",&pip_stvtime);
tree->SetBranchAddress("p_stvtime",&p_stvtime);
tree->SetBranchAddress("stVtimeDiff",&vtd);
tree->SetBranchAddress("p_sc_id_part",&p_sc_id_part);
tree->SetBranchAddress("pip_sc_id_part",&pip_sc_id_part);
tree->SetBranchAddress("pim_sc_id_part",&pim_sc_id_part);



//
Int_t nentries = (Int_t)tree->GetEntries();

vector<int> bit(13,0);
  for (int ijj=0; ijj<12; ijj++)
  {
    bit[ijj+1] = 1 << ijj;
  }

Int_t counting;
counting = 0;


for (Int_t zz=0; zz<nentries ;zz++) {
tree->GetEntry(zz);

dbetapim = betapim - betapim_tof;
      dbetapip = betapip - betapip_tof;
      dbetap = betaproton - betaproton_tof;




if ( bec_pi0_prob > 0.01  && (p_sec != pip_sec) && (p_sec != pim_sec) && (pip_sec != pim_sec) )
{
 counting = counting + 1;

cerr<<counting<<" "<<nentries<<endl;
//cerr<<p_sec<<" "<<pip_sec<<" "<<pim_sec<<endl;
//
int sect_p = p_sec;
int sect_pip = pip_sec;
int sect_pim = pim_sec;

double phi_rad_p = atan2(p_py,p_px)* TMath::RadToDeg(); ; //find phi in radians
  if(phi_rad_p<0.0000) //convert neg angles to pos
    {phi_rad_p = phi_rad_p+360.0;}

if(p_sec > 1)
 {
  phi_rad_p = phi_rad_p-(p_sec-1)*60;
 }
else
 {
   if (phi_rad_p >= 330)
   {phi_rad_p = phi_rad_p-360;}
 }

if (trigbit & bit[11])
{
 //cerr<<"pass the event"<<endl;
}
else
{
proton_total[sect_p-1]->Fill(p_sc_id_part, phi_rad_p);

if ((trigbit & bit[p_sec]) || (trigbit & bit[12]))
{proton_hit[sect_p-1]->Fill(p_sc_id_part, phi_rad_p);}
}
//
double phi_rad_pip = atan2(pip_py,pip_px)* TMath::RadToDeg(); ; //find phi in radians
  if(phi_rad_pip<0.0000) //convert neg angles to pos
    {phi_rad_pip = phi_rad_pip+360.0;}

if(pip_sec > 1)
 {
  phi_rad_pip = phi_rad_pip-(pip_sec-1)*60;
 }
else
 {
   if (phi_rad_pip >= 330)
   {phi_rad_pip = phi_rad_pip-360;}
 }

if (trigbit & bit[11])
{
 //cerr<<"pass the event"<<endl;
}
else
{
pip_total[sect_pip-1]->Fill(pip_sc_id_part, phi_rad_pip);

if ((trigbit & bit[pip_sec]) || (trigbit & bit[12]))
{pip_hit[sect_pip-1]->Fill(pip_sc_id_part, phi_rad_pip);}
}
//

double phi_rad_pim = atan2(pim_py,pim_px)* TMath::RadToDeg(); ; //find phi in radians
  if(phi_rad_pim<0.0000) //convert neg angles to pos
    {phi_rad_pim = phi_rad_pim+360.0;}

if(pim_sec > 1)
 {
  phi_rad_pim = phi_rad_pim-(pim_sec-1)*60;
 }
else
 {
   if (phi_rad_pim >= 330)
   {phi_rad_pim = phi_rad_pim-360;}
 }

if (trigbit & bit[11])
{
 //cerr<<"pass the event"<<endl;
}
else
{
pim_total[sect_pim-1]->Fill(pim_sc_id_part, phi_rad_pim);

if ((trigbit & bit[pim_sec]) || (trigbit & bit[12] ))
{pim_hit[sect_pim-1]->Fill(pim_sc_id_part, phi_rad_pim);}
}

  


}
}
//

for (int i=0; i<6; i++)
{
 proton_hit[i]->Divide(proton_total[i]);
 pip_hit[i]->Divide(pip_total[i]);
 pim_hit[i]->Divide(pim_total[i]);
 proton_hit[i]->Write();
 pip_hit[i]->Write();
 pim_hit[i]->Write();
}



cerr<<"FINISH FIRST LOOP";


}