#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
#include <time.h>
#include "TROOT.h"
#include "TStyle.h"
#include "TChain.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TLorentzVector.h"
#include "TLorentzRotation.h"
#include "TVector3.h"
#include "TMath.h"
#include "TH1I.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TTree.h"
//#include "defineHistos.h"
#include "TH1F.h"

float FidFunc_pos(float *, float *);

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



int main()
{
TLorentzVector MyProton, MyPip, MyPim, MyProton_Rest, My_Initial, MyPi0, MyPhoton, MM,MyPPip;
         TVector3 pro_mag, mm_mag;
         Float_t m_p = 0.938;
  float omega_p;
 float sum = 0;
  
  int Ebinmax =1;
  int beampol;
  char name[500],name1[500], name2[500],name9[500],name3[500];
  float temp[50],temp8[1500],temp9[1500],temp2[1500];
  float runnum, eventnum, top,  realphi, cosrealtheta, cos_theta_p_cm, beamE, p_px, p_py, p_pz, p_E, pip_px, pip_py, pip_pz, pip_E, pim_px, pim_py, pim_pz, pim_E;
  float cl_nop, cl_pi0, z, kin_beamE, kin_p_px, kin_p_py, kin_p_pz, kin_p_E, kin_pip_px, kin_pip_py, kin_pip_pz, kin_pip_E, kin_pim_px, kin_pim_py, kin_pim_pz, kin_pim_E, q_value, q_error;
  float kin_mm_px, kin_mm_py, kin_mm_pz, kin_mm_E, kin_mm_m, kin_o_px, kin_o_py, kin_o_pz, kin_o_E, o_m, kin_p_m, beta_cm, gamma_cm, o_pz_cm, o_pt, theta_cm, cos_theta_cm;
  float filenum;
  float kin_pip_mm_m, kin_pim_mm_m, pip_KE, pim_KE, mm_KE,X,Y, pip_m, pim_m, Q_o;
  float protonmass = .938272013;
  
  bool firstevent = true;
  int type = 0;
Float_t costetastar, degreepol, phistar, cosprocm, p_sec, pip_sec,pim_sec, p_sc_id_part, pip_sc_id_part,pim_sc_id_part;
double phi_rad_p, phi_rad_pip, phi_rad_pim;



//TFile *f = new TFile("/d/grid13/akbar/newanalysis/triggermap4.root");
TFile *f = new TFile("/d/grid13/akbar/newanalysis/triggermap_may2017v2.root");
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] = (TH2D*)f->Get(name1);
 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] = (TH2D*)f->Get(name1);
 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] = (TH2D*)f->Get(name1);
 pim_hit[i]->Sumw2();
}

TFile fff ("allcut2.root", "RECREATE", "Histograms from ntuples" );    
TH1D *rnd = new TH1D("random_test","random_test",1000,0,1);
TH1D *three_sec = new TH1D("thresec","threesec",86,1.1,5.4);
TH1D *two_sec = new TH1D("twosec","twosec",86,1.1,5.4);
TH1D *all_sec = new TH1D("allsec","allsec",86,1.1,5.4);
TH1D *stat_sec = new TH1D("statsec","statsec",12,0,6);

         qin.open("/d/grid13/akbar/g12_analysis_new/ksigma_mc_final_oct_new_trigger_sectorcut.txt",ios::out);
       //qin.open("/d/grid13/akbar/g12_analysis_new/ksigma_mc_allcut_oct_new_trigger_sectorcut.txt",ios::out);
       //qin.open("/d/grid13/akbar/g12_analysis_new/omega_mc_allmc_allcut_sept22_newtrigger_onlyfortest_usebelow.txt",ios::out);
      //qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_allcut_sept18_test_newtrigger.txt",ios::out);
      //qin.open("/d/grid13/akbar/g12_analysis_new/ksigma_mc_allcut_sept28_new_trigger_sectorcut.txt",ios::out); 
   //qin.open("/d/grid13/akbar/g12_analysis_new/ksigma_mc_allcut_sept28_new_trigger.txt",ios::out);
    //qin.open("/d/grid13/akbar/g12_analysis_new/omega_mc_allmc_allcut_sept22.txt",ios::out);
  //qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_sept19_allcut_trigger.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/ksigma_mc_allmc_allcut_sept19_trigger.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_allcut_sept18_test.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/omega_test_notrigger_sept18.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_allcut_august31_trigger.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/ksigma_mc_august17.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_allcut_august15_trigger.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_allcut_august15.txt",ios::out);
    //qin.open("/d/grid13/akbar/g12_analysis_new/omega_mc_allmc_allcut_august12.txt",ios::out);
     // qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_allcut_august7.txt",ios::out);
  // qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_additional_allcut_august7.txt",ios::out);
 // qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_june13_allcut_onlynew.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_june13_allcut_old.txt",ios::out);
  //qin.open("/d/grid13/akbar/g12_analysis_new/eta_mc_june13_allcut.txt",ios::out);
//qin.open("/d/grid13/akbar/g12_analysis_new/ksigma_mc_june11_allcut.txt",ios::out);
//qin.open("/d/grid6/akbar/g12_analysis_new/ksigma_mc_june_allcut.txt",ios::out);
//qin.open("/d/grid6/akbar/g12_analysis_new/MC/g12_omega_sim_new_jan2017_allcut_nosectorcut_30may2017.txt",ios::out); //output text file
//qin.open("/d/grid6/akbar/g12_analysis_new/g12_eta_sim_new_allcut_nosectorcut_19may2017.txt",ios::out);
    qin.setf(ios::showpoint|ios::fixed);
    qin.precision(5);

Float_t counting=0;

ifstream inputfile;
       sprintf(name,"/d/grid13/akbar/g12_analysis_new/ksigma_mc_june11.txt");
           // sprintf(name,"/d/grid6/akbar/g12_analysis_new/MC/g12_omega_sim_new_jan2017.txt");
         // sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_sept19.txt");
        // sprintf(name,"/d/grid13/akbar/g12_analysis_new/ksigma_mc_june11.txt");
        //sprintf(name,"/d/grid6/akbar/g12_analysis_new/MC/g12_omega_sim_new_jan2017.txt");
       //sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_sept19.txt");
      //sprintf(name,"/d/grid13/akbar/g12_analysis_new/ksigma_mc_june11.txt");
    //sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_august31.txt");
       //sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_august31.txt");
     //sprintf(name,"/d/grid13/akbar/g12_analysis_new/ksigma_mc_june11.txt");
    //sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_august15.txt");
   //sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_august15.txt");
   //sprintf(name,"/d/grid6/akbar/g12_analysis_new/MC/g12_omega_sim_new_jan2017.txt");
  //sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_allmc_august7.txt");
//sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_additional_august7.txt");
//sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_june13_onlynew.txt");
//sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_june13_old.txt");
//sprintf(name,"/d/grid13/akbar/g12_analysis_new/eta_mc_june13.txt");
//sprintf(name,"/d/grid13/akbar/g12_analysis_new/ksigma_mc_june11.txt");
//sprintf(name,"/d/grid6/akbar/g12_analysis_new/ksigma_mc_june.txt");
//sprintf(name,"/d/grid6/akbar/g12_analysis_new/MC/g12_omega_sim_new_jan2017.txt");
//sprintf(name,"/d/grid6/akbar/g12_analysis_new/g12_eta_sim_new_3april2017.txt");
    inputfile.open(name,ios::in);
    int counter = 0;
    while (!inputfile.eof()){
      for (int i = 0; i < /*38*/42 ; i++) inputfile >> temp[i]; 
      runnum = temp[0];
      eventnum = temp[1];
      top = temp[2];
      beampol = temp[3];
      realphi = temp[4];
      cosrealtheta = temp[5];
      cos_theta_p_cm = temp[6];
      beamE = temp[7] /1000.0;
      p_px = temp[8] /1000.0;
      p_py = temp[9] /1000.0;
      p_pz = temp[10] /1000.0;
      p_E = temp[11] /1000.0;
      pip_px = temp[12] /1000.0;
      pip_py = temp[13] /1000.0;
      pip_pz = temp[14] /1000.0;
      pip_E = temp[15] /1000.0;
      pim_px = temp[16] /1000.0;
      pim_py = temp[17] /1000.0;
      pim_pz = temp[18] /1000.0;
      pim_E = temp[19] /1000.0;
      cl_nop = temp[20];
      cl_pi0 = temp[21];
      z = temp[22];
      kin_beamE = temp[23] /1000.0;
      kin_p_px = temp[24] /1000.0;
      kin_p_py = temp[25] /1000.0;
      kin_p_pz = temp[26] /1000.0;
      kin_p_E = temp[27] /1000.0;
      kin_pip_px = temp[28] /1000.0;
      kin_pip_py = temp[29] /1000.0;
      kin_pip_pz = temp[30] /1000.0;
      kin_pip_E = temp[31] /1000.0;
      kin_pim_px = temp[32] /1000.0;
      kin_pim_py = temp[33] /1000.0;
      kin_pim_pz = temp[34] /1000.0;
      kin_pim_E = temp[35] /1000.0;
p_sc_id_part = temp[36];
pip_sc_id_part = temp[37];
pim_sc_id_part = temp[38];
p_sec = temp[39];
pip_sec = temp[40];
pim_sec = temp[41];

if(runnum != 56520 && runnum != 56544 && runnum != 56559 && runnum != 56619 && runnum != 56637 && runnum != 56585 && runnum > 56519)
      {

 ///proton 
 TVector3 temp_p = TVector3(kin_p_px, kin_p_py, kin_p_pz);
 double teta_p = temp_p.Theta()*(180/PI);


  
   phi_rad_p = atan2(p_py,p_px)* TMath::RadToDeg(); ; //find phi in radians
double p_mk = phi_rad_p;
  if(phi_rad_p<0.0000) //convert neg angles to pos
    {phi_rad_p = phi_rad_p+360.0;}
  
  Double_t phiSector_p = phi_rad_p;
 if (phiSector_p > 330.0) phiSector_p -= 360.;
  int sect_p = p_sec;


float parL[6][9];
 float parH[6][9];
 parL[0][0] = -2295.59; parL[0][1] = 0.000700362; parL[0][2] = 2266.32;
      parL[0][3] = -198822; parL[0][4] = 0.000363836; parL[0][5] = 198706;
      parL[0][6] = -198764; parL[0][7] = 3.78601e-05; parL[0][8] = 198763;
      parL[1][0] = -2295.45; parL[1][1] = 0.000812151; parL[1][2] = 2266.18;
      parL[1][3] = -198812; parL[1][4] = 0.000280686; parL[1][5] = 198716;
      parL[1][6] = -198764; parL[1][7] = 3.58984e-05; parL[1][8] = 198764;
      parL[2][0] = -2295.01; parL[2][1] = 0.000510464; parL[2][2] = 2266.32;
      parL[2][3] = -198801; parL[2][4] = 0.000167734; parL[2][5] = 198726;
      parL[2][6] = -198762; parL[2][7] = 1.81963e-05; parL[2][8] = 198765;
      parL[3][0] = -2295.19; parL[3][1] = 0.000596059; parL[3][2] = 2265.85;
      parL[3][3] = -198819; parL[3][4] = 0.000275454; parL[3][5] = 198709;
      parL[3][6] = -198764; parL[3][7] = 2.2017e-05; parL[3][8] = 198764;
      parL[4][0] = -2294.95; parL[4][1] = 0.000629762; parL[4][2] = 2265.82;
      parL[4][3] = -198817; parL[4][4] = 0.000300642; parL[4][5] = 198710;
      parL[4][6] = -198764; parL[4][7] = 2.36482e-05; parL[4][8] = 198764;
      parL[5][0] = -2294.57; parL[5][1] = 0.000587603; parL[5][2] = 2265.88;
      parL[5][3] = -198801; parL[5][4] = 0.000264759; parL[5][5] = 198726;
      parL[5][6] = -198763; parL[5][7] = 3.49751e-05; parL[5][8] = 198764;
      
      parH[0][0] = -2266.96; parH[0][1] = -0.000136134; parH[0][2] = 2294.8;
      parH[0][3] = -198725; parH[0][4] = -4.85729e-06; parH[0][5] = 198803;
      parH[0][6] = -198764; parH[0][7] = 2.91321e-05; parH[0][8] = 198764;
      parH[1][0] = -2266.48; parH[1][1] = -0.00054583; parH[1][2] = 2294.99;
      parH[1][3] = -198715; parH[1][4] = -0.00022386; parH[1][5] = 198813;
      parH[1][6] = -198763; parH[1][7] = 2.46949e-05; parH[1][8] = 198764;
      parH[2][0] = -2266.63; parH[2][1] = -0.00041358; parH[2][2] = 2294.55;
      parH[2][3] = -198732; parH[2][4] = -0.000165719; parH[2][5] = 198796;
      parH[2][6] = -198762; parH[2][7] = 2.36424e-05; parH[2][8] = 198765;
      parH[3][0] = -2266.56; parH[3][1] = -0.000353913; parH[3][2] = 2294.34;
      parH[3][3] = -198733; parH[3][4] = -0.000160694; parH[3][5] = 198795;
      parH[3][6] = -198762; parH[3][7] = 2.0461e-05; parH[3][8] = 198766;
      parH[4][0] = -2266.22; parH[4][1] = -0.000632337; parH[4][2] = 2294.37;
      parH[4][3] = -198728; parH[4][4] = -0.000223575; parH[4][5] = 198799;
      parH[4][6] = -198763; parH[4][7] = 3.39674e-05; parH[4][8] = 198765;
      parH[5][0] = -2266.3; parH[5][1] = -0.000349999; parH[5][2] = 2294;
      parH[5][3] = -198726; parH[5][4] = -0.00012299; parH[5][5] = 198802;
      parH[5][6] = -198762; parH[5][7] = 1.54008e-05; parH[5][8] = 198765;

 float phidiff = (60.*(float)(sect_p - 1)) - phiSector_p;


double p_p_mag = temp_p.Mag();
  
    
      float fidvar[2] = {p_p_mag,teta_p};
      
      float lowReg = 0;
      lowReg = FidFunc_pos(fidvar,parL[sect_p - 1]);
      float highReg = 0;
      highReg = FidFunc_pos(fidvar,parH[sect_p - 1]);

 float pass_low, pass_high;
 int fid_p = 0;

        pass_low = lowReg;
        pass_high = highReg;
        if ((phidiff > pass_low) && (phidiff < pass_high))
        {
          fid_p = 1;
        }
         


///pip
 TVector3 temp_pip = TVector3(kin_pip_px, kin_pip_py, kin_pip_pz);
 double teta_pip = temp_pip.Theta()*(180/PI);



    phi_rad_pip = atan2(pip_py,pip_px)* TMath::RadToDeg(); ; //find phi in radians
double pip_mk = phi_rad_pip;
  if(phi_rad_pip<0.0000) //convert neg angles to pos
    {phi_rad_pip = phi_rad_pip+360.0;}
  
  Double_t phiSector_pip = phi_rad_pip;
 if (phiSector_pip > 330.0) phiSector_pip -= 360.;
 
 int sect_pip = pip_sec;
 

//float parL[6][9];
 //float parH[6][9];
 parL[0][0] = -2295.59; parL[0][1] = 0.000700362; parL[0][2] = 2266.32;
      parL[0][3] = -198822; parL[0][4] = 0.000363836; parL[0][5] = 198706;
      parL[0][6] = -198764; parL[0][7] = 3.78601e-05; parL[0][8] = 198763;
      parL[1][0] = -2295.45; parL[1][1] = 0.000812151; parL[1][2] = 2266.18;
      parL[1][3] = -198812; parL[1][4] = 0.000280686; parL[1][5] = 198716;
      parL[1][6] = -198764; parL[1][7] = 3.58984e-05; parL[1][8] = 198764;
      parL[2][0] = -2295.01; parL[2][1] = 0.000510464; parL[2][2] = 2266.32;
      parL[2][3] = -198801; parL[2][4] = 0.000167734; parL[2][5] = 198726;
      parL[2][6] = -198762; parL[2][7] = 1.81963e-05; parL[2][8] = 198765;
      parL[3][0] = -2295.19; parL[3][1] = 0.000596059; parL[3][2] = 2265.85;
      parL[3][3] = -198819; parL[3][4] = 0.000275454; parL[3][5] = 198709;
      parL[3][6] = -198764; parL[3][7] = 2.2017e-05; parL[3][8] = 198764;
      parL[4][0] = -2294.95; parL[4][1] = 0.000629762; parL[4][2] = 2265.82;
      parL[4][3] = -198817; parL[4][4] = 0.000300642; parL[4][5] = 198710;
      parL[4][6] = -198764; parL[4][7] = 2.36482e-05; parL[4][8] = 198764;
      parL[5][0] = -2294.57; parL[5][1] = 0.000587603; parL[5][2] = 2265.88;
      parL[5][3] = -198801; parL[5][4] = 0.000264759; parL[5][5] = 198726;
      parL[5][6] = -198763; parL[5][7] = 3.49751e-05; parL[5][8] = 198764;
      
      parH[0][0] = -2266.96; parH[0][1] = -0.000136134; parH[0][2] = 2294.8;
      parH[0][3] = -198725; parH[0][4] = -4.85729e-06; parH[0][5] = 198803;
      parH[0][6] = -198764; parH[0][7] = 2.91321e-05; parH[0][8] = 198764;
      parH[1][0] = -2266.48; parH[1][1] = -0.00054583; parH[1][2] = 2294.99;
      parH[1][3] = -198715; parH[1][4] = -0.00022386; parH[1][5] = 198813;
      parH[1][6] = -198763; parH[1][7] = 2.46949e-05; parH[1][8] = 198764;
      parH[2][0] = -2266.63; parH[2][1] = -0.00041358; parH[2][2] = 2294.55;
      parH[2][3] = -198732; parH[2][4] = -0.000165719; parH[2][5] = 198796;
      parH[2][6] = -198762; parH[2][7] = 2.36424e-05; parH[2][8] = 198765;
      parH[3][0] = -2266.56; parH[3][1] = -0.000353913; parH[3][2] = 2294.34;
      parH[3][3] = -198733; parH[3][4] = -0.000160694; parH[3][5] = 198795;
      parH[3][6] = -198762; parH[3][7] = 2.0461e-05; parH[3][8] = 198766;
      parH[4][0] = -2266.22; parH[4][1] = -0.000632337; parH[4][2] = 2294.37;
      parH[4][3] = -198728; parH[4][4] = -0.000223575; parH[4][5] = 198799;
      parH[4][6] = -198763; parH[4][7] = 3.39674e-05; parH[4][8] = 198765;
      parH[5][0] = -2266.3; parH[5][1] = -0.000349999; parH[5][2] = 2294;
      parH[5][3] = -198726; parH[5][4] = -0.00012299; parH[5][5] = 198802;
      parH[5][6] = -198762; parH[5][7] = 1.54008e-05; parH[5][8] = 198765;

  phidiff = (60.*(float)(sect_pip - 1)) - phiSector_pip;


double pip_p_mag = temp_pip.Mag();
 
    
      float fidvarpip[2] = {pip_p_mag,teta_pip};
      
       lowReg = 0;
      lowReg = FidFunc_pos(fidvarpip,parL[sect_pip - 1]);
       highReg = 0;
      highReg = FidFunc_pos(fidvarpip,parH[sect_pip - 1]);

 //float pass_low, pass_high;
 int fid_pip = 0;

        pass_low = lowReg;
        pass_high = highReg;
        if ((phidiff > pass_low) && (phidiff < pass_high))
        {
          fid_pip = 1;
        }
         
   
///pim
TVector3 temp_pim = TVector3(kin_pim_px, kin_pim_py, kin_pim_pz);

 double teta_pim = temp_pim.Theta()*(180/PI);



   phi_rad_pim = atan2(pim_py,pim_px)* TMath::RadToDeg(); ; //find phi in radians
double pim_mk = phi_rad_pim;
  if(phi_rad_pim<0.0000) //convert neg angles to pos
    {phi_rad_pim = phi_rad_pim+360.0;}
  
  Double_t phiSector_pim = phi_rad_pim;
 if (phiSector_pim > 330.0) phiSector_pim -= 360.;
  int sect_pim = pim_sec;
 

//float parL[6][9];
// float parH[6][9];
 parL[0][0] = -2295.59; parL[0][1] = 0.000700362; parL[0][2] = 2266.32;
       parL[0][0] = -15810.2; parL[0][1] = 9.9663e-05; parL[0][2] = 15783.2;
      parL[0][3] = -538215; parL[0][4] = 9.84788e-05; parL[0][5] = 538135;
      parL[0][6] = -678156; parL[0][7] = 1.42357e-05; parL[0][8] = 678167;
      
      parL[1][0] = -15810.7; parL[1][1] = 0.0001; parL[1][2] = 15783.7;
      parL[1][3] = -538209; parL[1][4] = 8.31014e-05; parL[1][5] = 538141;
      parL[1][6] = -678155; parL[1][7] = 1.32519e-05; parL[1][8] = 678167;
      
      parL[2][0] = -37808.6; parL[2][1] = 4.38676e-05; parL[2][2] = 37782.3;
      parL[2][3] = -678187; parL[2][4] = 3.0748e-05; parL[2][5] = 678136;
      parL[2][6] = -678155; parL[2][7] = 1.13289e-05; parL[2][8] = 678168;
      
      parL[3][0] = -37809.1; parL[3][1] = 5.52776e-05; parL[3][2] = 37781.8;
      parL[3][3] = -678197; parL[3][4] = 4.21772e-05; parL[3][5] = 678126;
      parL[3][6] = -678156; parL[3][7] = 9.59678e-06; parL[3][8] = 678167;
      
      parL[4][0] = -37809.2; parL[4][1] = 5.79073e-05; parL[4][2] = 37781.7;
      parL[4][3] = -678207; parL[4][4] = 7.72987e-05; parL[4][5] = 678115;
      parL[4][6] = -678156; parL[4][7] = 1.63781e-05; parL[4][8] = 678166;
 parL[5][0] = -37808.9; parL[5][1] = 4.70858e-05; parL[5][2] = 37781.9;
            parL[5][3] = -678198; parL[5][4] = 9.27108e-05; parL[5][5] = 678125;
            parL[5][6] = -678155; parL[5][7] = 1.63781e-05; parL[5][8] = 678167;
      
 parH[0][0] = -15783.8; parH[0][1] = -0.0001; parH[0][2] = 15809.5;
      parH[0][3] = -538148; parH[0][4] = -8.77852e-05; parH[0][5] = 538201;
      parH[0][6] = -678155; parH[0][7] = 1.67856e-05; parH[0][8] = 678168;
      parH[1][0] = -37782.2; parH[1][1] = -9.78359e-05; parH[1][2] = 37808.7;
      parH[1][3] = -678123; parH[1][4] = -9.99919e-05; parH[1][5] = 678199;
      parH[1][6] = -678155; parH[1][7] = 1.71348e-05; parH[1][8] = 678168;
      parH[2][0] = -37782.3; parH[2][1] = -6.62596e-05; parH[2][2] = 37808.6;
      parH[2][3] = -678133; parH[2][4] = -4.84207e-05; parH[2][5] = 678190;
      parH[2][6] = -678155; parH[2][7] = 1.47442e-05; parH[2][8] = 678168;
      parH[3][0] = -37782.5; parH[3][1] = -3.18908e-05; parH[3][2] = 37808.4;
      parH[3][3] = -678135; parH[3][4] = -1.21321e-05; parH[3][5] = 678188;
      parH[3][6] = -678155; parH[3][7] = 1.22068e-05; parH[3][8] = 678168;
      
      parH[4][0] = -37782.3; parH[4][1] = -6.28596e-05; parH[4][2] = 37808.5;
      parH[4][3] = -678136; parH[4][4] = -0.0001; parH[4][5] = 678187;
      parH[4][6] = -678154; parH[4][7] = 1.79547e-05; parH[4][8] = 678168;
      parH[5][0] = -37782.4; parH[5][1] = -6.08308e-05; parH[5][2] = 37808.4;
      parH[5][3] = -678131; parH[5][4] = -5.88378e-05; parH[5][5] = 678192;
      parH[5][6] = -678155; parH[5][7] = 1.44427e-05; parH[5][8] = 678168;
      
     

 phidiff = (60.*(float)(sect_pim - 1)) - phiSector_pim;


double pim_p_mag = temp_pim.Mag();
  
    
  float fidvarpim[2] = {pim_p_mag,teta_pim};
      
    lowReg = 0;
      lowReg = FidFunc_pos(fidvarpim,parL[sect_pim - 1]);
   highReg = 0;
      highReg = FidFunc_pos(fidvarpim,parH[sect_pim - 1]);

 //float pass_low, pass_high;
 int fid_pim = 0;

        pass_low = lowReg;
        pass_high = highReg;
        if ((phidiff > pass_low) && (phidiff < pass_high))
        {
          fid_pim = 1;
        }

int scid;
int passed_p =1;
scid = int (p_sc_id_part);

if( p_sec == 1 && (scid == 6 || scid == 35 || scid == 40 || scid == 41 || scid == 50 || scid == 56)) // bad paddles in sector 1
				passed_p = 0;
			if( p_sec == 2 && (scid == 2 || scid == 8 || scid == 34 || scid == 35 || scid == 41 || scid == 44 || scid == 50 || scid == 54 || scid == 56)) //bad paddles in sector 2
				passed_p =0;
			if( p_sec == 3 && (scid == 11 || scid == 35 || scid == 40 || scid == 41 || scid == 56)) // bad paddles in sector 3
				passed_p = 0;
			if( p_sec == 4 && (scid == 41 || scid == 48)) // bad paddles in sector 4
				passed_p = 0;
			if( p_sec == 5 && scid == 48) // bad paddles in sector 5
				passed_p = 0;
			if( p_sec == 6 && (scid == 1 || scid == 5 || scid == 33 || scid == 56)) // bad paddles in sector 6
				passed_p = 0;
if( p_sec == 1 && (scid == 25 || scid == 26) ) // bad paddles in sector 1
					passed_p = 0;

				if( p_sec == 2 && (scid == 18 || scid == 25 || scid ==27) ) // bad paddles in sector 2
					passed_p = 0;

				if( p_sec == 3 && (scid == 1 || scid == 18 || scid == 23)) // bad paddles in sector 3
					passed_p = 0;

				if( p_sec == 4 && (scid == 8 || scid == 19)) // bad paddles in sector 4
					passed_p = 0;

				if( p_sec == 6 && (scid == 24)) // bad paddles in sector 6
					passed_p = 0;



int passed_pip =1;
scid = int (pip_sc_id_part);

if( pip_sec == 1 && (scid == 6 || scid == 35 || scid == 40 || scid == 41 || scid == 50 || scid == 56)) // bad paddles in sector 1
				passed_pip = 0;
			if( pip_sec == 2 && (scid == 2 || scid == 8 || scid == 34 || scid == 35 || scid == 41 || scid == 44 || scid == 50 || scid == 54 || scid == 56)) 
				passed_pip =0;
			if(pip_sec == 3 && (scid == 11 || scid == 35 || scid == 40 || scid == 41 || scid == 56)) // bad paddles in sector 3
				passed_pip = 0;
			if( pip_sec == 4 && (scid == 41 || scid == 48)) // bad paddles in sector 4
				passed_pip = 0;
			if( pip_sec == 5 && scid == 48) // bad paddles in sector 5
				passed_pip = 0;
			if( pip_sec == 6 && (scid == 1 || scid == 5 || scid == 33 || scid == 56)) // bad paddles in sector 6
				passed_pip = 0;
if( pip_sec == 1 && (scid == 25 || scid == 26) ) // bad paddles in sector 1
					passed_pip = 0;

				if( pip_sec == 2 && (scid == 18 || scid == 25 || scid ==27) ) // bad paddles in sector 2
					passed_pip = 0;

				if( pip_sec == 3 && (scid == 1 || scid == 18 || scid == 23)) // bad paddles in sector 3
					passed_pip = 0;

				if( pip_sec == 4 && (scid == 8 || scid == 19)) // bad paddles in sector 4
					passed_pip = 0;

				if( pip_sec == 6 && (scid == 24)) // bad paddles in sector 6
					passed_pip = 0;



int passed_pim =1;
scid = int (pim_sc_id_part);

if( pim_sec == 1 && (scid == 6 || scid == 35 || scid == 40 || scid == 41 || scid == 50 || scid == 56)) // bad paddles in sector 1
				passed_pim = 0;
			if( pim_sec == 2 && (scid == 2 || scid == 8 || scid == 34 || scid == 35 || scid == 41 || scid == 44 || scid == 50 || scid == 54 || scid == 56)) 
				passed_pim =0;
			if(pim_sec == 3 && (scid == 11 || scid == 35 || scid == 40 || scid == 41 || scid == 56)) // bad paddles in sector 3
				passed_pim = 0;
			if( pim_sec == 4 && (scid == 41 || scid == 48)) // bad paddles in sector 4
				passed_pim = 0;
			if( pim_sec == 5 && scid == 48) // bad paddles in sector 5
				passed_pim = 0;
			if( pim_sec == 6 && (scid == 1 || scid == 5 || scid == 33 || scid == 56)) // bad paddles in sector 6
				passed_pim = 0;
if( pim_sec == 1 && (scid == 25 || scid == 26) ) // bad paddles in sector 1
					passed_pim = 0;

				if( pim_sec == 2 && (scid == 18 || scid == 25 || scid ==27) ) // bad paddles in sector 2
					passed_pim = 0;

				if( pim_sec == 3 && (scid == 1 || scid == 18 || scid == 23)) // bad paddles in sector 3
					passed_pim = 0;

				if( pim_sec == 4 && (scid == 8 || scid == 19)) // bad paddles in sector 4
					passed_pim = 0;

				if( pim_sec == 6 && (scid == 24)) // bad paddles in sector 6
					passed_pim = 0;

//cospi0 cut
MyProton.SetPxPyPzE( kin_p_px,kin_p_py, kin_p_pz, kin_p_E );
                          MyPip.SetPxPyPzE( kin_pip_px,kin_pip_py, kin_pip_pz, kin_pip_E );
                          MyPim.SetPxPyPzE( kin_pim_px,kin_pim_py, kin_pim_pz, kin_pim_E);
                          MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
                          My_Initial.SetPxPyPzE(0.0,0.0,kin_beamE, kin_beamE+0.938);
                         MyPhoton.SetPxPyPzE(0.0,0.0,kin_beamE, kin_beamE);
			 TLorentzVector pi0_init = My_Initial - MyProton - MyPip - MyPim;
                    
			     TLorentzRotation cm_r;
			     TVector3 cm_v = My_Initial.BoostVector();
			     cm_r.Boost(-cm_v);
			     TLorentzVector pi0_cm2 = cm_r * pi0_init;
                             


double costheta_pi0 = pi0_cm2.CosTheta();
//
//trigger simulation
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;}
 }

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;}
 }

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;}
 }


int binphi_p = (phi_rad_p-(-30))/2+1;
int binphi_pip = (phi_rad_pip-(-30))/2+1;
int binphi_pim = (phi_rad_pim-(-30))/2+1;

int bintof_p = p_sc_id_part;
int bintof_pip = pip_sc_id_part;
int bintof_pim = pim_sc_id_part;

double prob_p = proton_hit[sect_p-1]->GetBinContent(bintof_p,binphi_p);
double prob_pip = pip_hit[sect_pip-1]->GetBinContent(bintof_pip,binphi_pip);
double prob_pim = pim_hit[sect_pim-1]->GetBinContent(bintof_pim,binphi_pim);

double prob_trigger_hit=1;
//

double v_p = rand() % 1000 +1;
double random_p = v_p/1000.0 ;

double v_pip = rand() % 1000 +1;
double random_pip = v_pip/1000.0 ;

double v_pim = rand() % 1000 +1;
double random_pim = v_pim/1000.0 ;

double v_tagger = rand() % 1000 +1;
double random_tagger = v_tagger/1000.0 ;

int hit_p, hit_pim, hit_pip, hit_total;
hit_p =0;
hit_pim=0;
hit_pip=0;
hit_total=0;

if (random_p <= prob_p)
{
 hit_p =1;
}

if (random_pip <= prob_pip)
{
 hit_pip =1;
}

if (random_pim <= prob_pim)
{
 hit_pim =1;
}
hit_total = hit_p + hit_pip + hit_pim;

int status_trigger =0;

if (hit_total ==3)
{
 status_trigger=1;
}

if (hit_total ==2 && random_tagger <= 0.51 && beamE < 3.6)
{
 status_trigger=1;
 if(fid_p ==1 && fid_pip==1 && fid_pim==1 && passed_p==1 && passed_pip==1 && passed_pim ==1 && costheta_pi0 < 0.99)
 {two_sec->Fill(beamE);}
}

if (hit_total ==2 && beamE > 3.6)
{
 status_trigger=1;
 if(fid_p ==1 && fid_pip==1 && fid_pim==1 && passed_p==1 && passed_pip==1 && passed_pim ==1 && costheta_pi0 < 0.99)
 {two_sec->Fill(beamE);}
}

//
cerr<<counting<<"  "<<hit_total<<" "<<status_trigger<<endl;
if (fid_p ==1 && fid_pip==1 && fid_pim==1 && passed_p==1 && passed_pip==1 && passed_pim ==1 && costheta_pi0 < 0.99 && /*random < prob_trigger_hit*/ (status_trigger==1) && ((sect_p != sect_pip) && (sect_p != sect_pim) && (sect_pip != sect_pim)))
{
counting=counting+1;
all_sec->Fill(beamE);
 stat_sec->Fill(5);
if ((sect_p != sect_pip) && (sect_p != sect_pim) && (sect_pip != sect_pim))
{
  three_sec->Fill(beamE);
  if (beamE <3.6)
  {stat_sec->Fill(3);}
}
else
{
 //two_sec->Fill(beamE);
 if (beamE <3.6)
 {stat_sec->Fill(2);}
}

qin << runnum << "  " << eventnum << "  4" << "  "<<beampol << endl
          << realphi<< "  "<<cosrealtheta <<"  " << cos_theta_p_cm << endl
          << beamE*1000 << endl
          << p_px*1000 << "  " << p_py*1000 << "  " << p_pz*1000 << "  " << p_E*1000 << endl
          << pip_px*1000 << "  " << pip_py*1000 << "  " << pip_pz*1000 << "  " << pip_E*1000 << endl
          << pim_px*1000 << "  " << pim_py*1000 << "  " << pim_pz*1000 << "  " << pim_E*1000 << endl
          << cl_nop << "  " << cl_pi0 << "  " << z << endl
          << kin_beamE*1000 << endl
          << kin_p_px*1000 << "  " << kin_p_py*1000 << "  " << kin_p_pz*1000 << "  " << kin_p_E*1000 << endl
          << kin_pip_px*1000 << "  " << kin_pip_py*1000 << "  " << kin_pip_pz*1000 << "  " << kin_pip_E*1000 << endl
          << kin_pim_px*1000 << "  " << kin_pim_py*1000 << "  " << kin_pim_pz*1000 << "  " << kin_pim_E*1000 << endl;
	  

}


}
}
cerr<<"FINISHHH"<<endl;
rnd->Write();
all_sec->Write();
two_sec->Write();
three_sec->Write();
three_sec->Divide(two_sec);
three_sec->Write();
stat_sec->Write();
}

float FidFunc_pos(float *x, float *par)
{
  float xx = x[0]; //momentum
  float yy = x[1]; //theta
  float a = (par[0]*pow(xx,par[1]) + par[2]);
  float b = (par[3]*pow(xx,par[4]) + par[5]);
  float c = (par[6]*pow(xx,par[7]) + par[8]);
  float f =  a -  b/(yy - c);
  return 	f;
}