#include <iostream>
using std::cout;
using std::endl;
using namespace TMath;
#include <cstdlib>
#include <iostream>
#include <map>
#include <string>

#include "TChain.h"
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TObjString.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TLorentzVector.h"

#include <fstream>
#include <sstream>
#if not defined(__CINT__) || defined(__MAKECINT__)
// needs to be included when makecint runs (ACLIC)
#include "TMVA/Factory.h"
#include "TMVA/Tools.h"
#endif
double htriton_mass;

int snn = 3 ;
const int ncentcut=3;

//Significance = 30.8, Sign = 19.7, piDca > 5.0 cm, proDca > 0.4 cm, deDca < 1.3 cm, decayL > 10.0 cm, ht_chi2topo < 8.0, ht_chi2ndf < 2.0, mean = 2.9924, sigma = 0.0013; 3,0,3,4,3,1

//Significance = 26.1, Sign = 24.8, piDca > 3.0 cm, proDca > 0.7 cm, deDca < 0.8 cm, decayL > 10.0 cm, ht_chi2topo < 6.0, ht_chi2ndf < 3.0, mean = 2.9930, sigma = 0.0015; 1,1,2,4,2,2

const int ncutdcapionbins = 5;
float dcaPionCuts[ncutdcapionbins] = {1, 2, 3, 5, 7};

const int ncutdcaprobins = 4;
float dcaProCuts[ncutdcaprobins] = {0.3, 0.6, 0.7, 1.0};

const int ncutdcadebins = 4;
float dcaDeCuts[ncutdcadebins] = {0.3, 0.6, 1, 2};

const int ncutdecaylengthbins = 4;
float decayLCuts[ncutdecaylengthbins] = {2, 4, 10, 12};

const int ncutchi2topo = 5;
float chi2topoCuts[ncutchi2topo] = {2, 4, 6, 8, 10};
//const int ncutchi2topo = 3;
//float chi2topoCuts[ncutchi2topo] = { 6, 8, 10};

const int ncutchi2ndf = 3;
float chi2ndfCuts[ncutchi2ndf] = {1, 2, 4};

const int ndP=11;
float dPCuts[ndP]={0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9,2.1,2.5, 3.0};

const int npP=8;
float pPCuts[npP]={1., 1.2, 1.4, 1.6, 1.8, 2., 2.5, 3.0};
TH1F* sHtmass_scanP[npP][ndP];

TH1F* sHtmass[ncutdcapionbins][ncutdcaprobins][ncutdcadebins][ncutdecaylengthbins][ncutchi2topo][ncutchi2ndf];
TH1F* sHtmassPRot[ncutdcapionbins][ncutdcaprobins][ncutdcadebins][ncutdecaylengthbins][ncutchi2topo][ncutchi2ndf];

Int_t           brunid;
Int_t           beventid;
Float_t         bVz;
Float_t         bVx;
Float_t         bVy;
Int_t           brefmult;
Int_t           btofmult;
Double_t        refmultcor;
Double_t        reweight;
Int_t           cent9;
Int_t           bparticleid;
Float_t         bparticlemass;
Float_t         bpx;
Float_t         bpy;
Float_t         bpz;
Float_t         ht_chi2;
Float_t         ht_NDF;
Float_t         chi2primary_proton;
Float_t         chi2primary_pi;
Float_t         chi2primary_d;
Float_t         ht_ldl;
Float_t         ht_dl;
Float_t         ht_l;
Float_t         ht_chi2topo;
Float_t         ht_chi2ndf;
Float_t         bpionpx;
Float_t         bpionpy;
Float_t         bpionpz;
Float_t         bprotonpx;
Float_t         bprotonpy;
Float_t         bprotonpz;
Float_t         bdpx;
Float_t         bdpy;
Float_t         bdpz;
Float_t         bzdeuteron;
Float_t         bpionnsigma;
Float_t         bprotonsigma;
Float_t         dca_pi;
Float_t         dedx_pi;
Float_t         dca_proton;
Float_t         dedx_proton;
Float_t         dca_deuteron;
Float_t         dedx_deuteron;
Int_t           nhits_pi;
Int_t           nhits_proton;
Int_t           nhits_deuteron;
Float_t         bdedx;
Float_t         bmcpx;
Float_t         bmcpy;
Float_t         bmcpz;
Int_t           bismc;


Char_t path_ch_b[10000];
Char_t path_ch_s[10000];
Char_t ofile_b[10000];
Char_t ofile_s[10000];
float mmin;
float mmax;
Float_t he3M2_max = 2.4;
Float_t he3M2_min = 1.7;
void ana_tree_htriton_Invmass(int pid=3004, bool Scan_dca_L_chi2=true, bool Scan_P=false){

    htriton_mass = 2.991;

    TChain bhtriton_tree("htriton_tree");
    TChain shtriton_tree("htriton3_tree");

    int nChains1 = -999999;
    int nChains2 = -999999;
    ifstream fin1;
    ifstream fin2;
    TFile* out_file_b;
    TFile* out_file_s;

    double massmin;
    double massmax;
    double nbin;

    if(pid==3004){massmin=2.95; massmax=3.05; nbin=100; mmin=2.9885; mmax=2.9975;}
    if(pid==3005){massmin=3.85, massmax=4.0; nbin=75; mmin=3.9192; mmax=3.9288;}
    cout<<"Now the particle is "<<pid<<endl;    

    //    if(snn==3){ycm = -1.045;}

    fin2.open("/home/fengchu/star_workdir/33.standard_reco/taxi_ht3b/fullfilelist/filelist_s.txt");
    if(gSystem->AccessPathName(Form("Invmass_sig_%d_PID2_dE1.5_he_ptcut.root",pid))){
        cout<<"WARNING! File already Exit! Do you want to recreate it?(y/n)"<<endl;
        char choice=cin.get();
        if(choice!="y") return;
    }
    out_file_s = new TFile(Form("Invmass_sig_%d_PID2_dE1.5_he_ptcut.root",pid),"recreate");
    //fin2.open("/home/fengchu/star_workdir/33.standard_reco/taxi_ht3b_pRot/filelist_PRot.txt");
    //out_file_s = new TFile(Form("Invmass_PRot_sig_%d_PID2_dE1.5_he_ptcut.root",pid),"recreate");
    nChains2 = 1;

    for(int i=0; i<nChains2; i++){
        fin2 >> path_ch_s;
        shtriton_tree.Add(path_ch_s);
        cout<<path_ch_s<<"  "<<shtriton_tree.GetEntries() <<endl;
    }//ichain
    cout << " done TChain-ing.. " << endl;



    shtriton_tree.SetBranchAddress("brunid", &brunid);
    shtriton_tree.SetBranchAddress("beventid", &beventid);
    shtriton_tree.SetBranchAddress("bVz", &bVz);
    shtriton_tree.SetBranchAddress("bVx", &bVx);
    shtriton_tree.SetBranchAddress("bVy", &bVy);
    shtriton_tree.SetBranchAddress("brefmult", &brefmult);
    shtriton_tree.SetBranchAddress("btofmult", &btofmult);
    shtriton_tree.SetBranchAddress("refmultcor", &refmultcor);
    shtriton_tree.SetBranchAddress("reweight", &reweight);
    shtriton_tree.SetBranchAddress("cent9", &cent9);
    shtriton_tree.SetBranchAddress("bparticleid", &bparticleid);
    shtriton_tree.SetBranchAddress("bparticlemass", &bparticlemass);
    shtriton_tree.SetBranchAddress("bpx", &bpx);
    shtriton_tree.SetBranchAddress("bpy", &bpy);
    shtriton_tree.SetBranchAddress("bpz", &bpz);
    shtriton_tree.SetBranchAddress("ht_chi2", &ht_chi2);
    shtriton_tree.SetBranchAddress("ht_NDF", &ht_NDF);
    shtriton_tree.SetBranchAddress("chi2primary_proton", &chi2primary_proton);
    shtriton_tree.SetBranchAddress("chi2primary_pi", &chi2primary_pi);
    shtriton_tree.SetBranchAddress("chi2primary_d", &chi2primary_d);
    shtriton_tree.SetBranchAddress("ht_ldl", &ht_ldl);
    shtriton_tree.SetBranchAddress("ht_dl", &ht_dl);
    shtriton_tree.SetBranchAddress("ht_l", &ht_l);
    shtriton_tree.SetBranchAddress("ht_chi2topo", &ht_chi2topo);
    shtriton_tree.SetBranchAddress("ht_chi2ndf", &ht_chi2ndf);
    shtriton_tree.SetBranchAddress("bpionpx", &bpionpx);
    shtriton_tree.SetBranchAddress("bpionpy", &bpionpy);
    shtriton_tree.SetBranchAddress("bpionpz", &bpionpz);
    shtriton_tree.SetBranchAddress("bprotonpx", &bprotonpx);
    shtriton_tree.SetBranchAddress("bprotonpy", &bprotonpy);
    shtriton_tree.SetBranchAddress("bprotonpz", &bprotonpz);
    shtriton_tree.SetBranchAddress("bdpx", &bdpx);
    shtriton_tree.SetBranchAddress("bdpy", &bdpy);
    shtriton_tree.SetBranchAddress("bdpz", &bdpz);
    shtriton_tree.SetBranchAddress("bzdeuteron", &bzdeuteron);
    shtriton_tree.SetBranchAddress("bpionnsigma", &bpionnsigma);
    shtriton_tree.SetBranchAddress("bprotonsigma", &bprotonsigma);
    shtriton_tree.SetBranchAddress("dca_pi", &dca_pi);
    shtriton_tree.SetBranchAddress("dedx_pi", &dedx_pi);
    shtriton_tree.SetBranchAddress("dca_proton", &dca_proton);
    shtriton_tree.SetBranchAddress("dedx_proton", &dedx_proton);
    shtriton_tree.SetBranchAddress("dca_deuteron", &dca_deuteron);
    shtriton_tree.SetBranchAddress("dedx_deuteron", &dedx_deuteron);
    shtriton_tree.SetBranchAddress("nhits_pi", &nhits_pi);
    shtriton_tree.SetBranchAddress("nhits_proton", &nhits_proton);
    shtriton_tree.SetBranchAddress("nhits_deuteron", &nhits_deuteron);
    shtriton_tree.SetBranchAddress("bdedx", &bdedx);
    shtriton_tree.SetBranchAddress("bmcpx", &bmcpx);
    shtriton_tree.SetBranchAddress("bmcpy", &bmcpy);
    shtriton_tree.SetBranchAddress("bmcpz", &bmcpz);
    shtriton_tree.SetBranchAddress("bismc", &bismc);

    //create th1f
    //if(Scan_dca_L_chi2){
    cout<<"Scan_dca_L_chi2: "<<Scan_dca_L_chi2 << endl;
    int cutIndex[6]={0,0,3,2,4,2};
    if(Scan_dca_L_chi2){
        for(int ipiondca=0; ipiondca<ncutdcapionbins; ipiondca++){
            for(int iprodca=0; iprodca<ncutdcaprobins; iprodca++){
                for(int idedca=0; idedca<ncutdcadebins; idedca++){
                    for(int idecayL=0; idecayL<ncutdecaylengthbins; idecayL++){
                        for(int ichi2topo=0; ichi2topo<ncutchi2topo; ichi2topo++){
                            for(int ichi2ndf=0; ichi2ndf<ncutchi2ndf; ichi2ndf++){

                                    /*
                                    if(ipiondca!=cutIndex[0]) continue;
                                    if(iprodca!=cutIndex[1]) continue;
                                    if(idedca!=cutIndex[2]) continue;
                                    if(idecayL!=cutIndex[3]) continue;
                                    if(ichi2topo!=cutIndex[4]) continue;
                                    if(ichi2ndf!=cutIndex[5]) continue;
                                    */
                                sHtmass[ipiondca][iprodca][idedca][idecayL][ichi2topo][ichi2ndf] = new TH1F(Form("sHtmass_pidca%d_prodca%d_dedca%d_decayL%d_chi2topo%d_chi2ndf%d",ipiondca,iprodca,idedca,idecayL,ichi2topo,ichi2ndf),Form("sHtmass_pidca%d_prodca%d_dedca%d_decayL%d_chi2topo%d_chi2ndf%d",ipiondca,iprodca,idedca,idecayL,ichi2topo,ichi2ndf),nbin,massmin,massmax);
                                //sHtmassPRot[ipiondca][iprodca][idedca][idecayL][ichi2topo][ichi2ndf] = new TH1F(Form("sHtmassPRot_pidca%d_prodca%d_dedca%d_decayL%d_chi2topo%d_chi2ndf%d",ipiondca,iprodca,idedca,idecayL,ichi2topo,ichi2ndf),Form("sHtmassPRot_pidca%d_prodca%d_dedca%d_decayL%d_chi2topo%d_chi2ndf%d",ipiondca,iprodca,idedca,idecayL,ichi2topo,ichi2ndf),nbin,massmin,massmax);
                            }
                        }
                    }
                }
            }
        }

        Long64_t n_signal_Entries = shtriton_tree.GetEntries();
        cout<<"--------------------------------------------"<<endl;
        cout<<"Total entries for signal "<<n_signal_Entries<< endl;

        for (Long64_t iEntry = 0; iEntry<=n_signal_Entries; iEntry++){
//        for (Long64_t iEntry = 0; iEntry<=1000000; iEntry++){
            shtriton_tree.GetEntry(iEntry);

            const double Mpion=0.1396;
            const double Mp=0.938;
            TVector3 Ppion(bpionpx,bpionpy,bpionpz);
            TVector3 Pp(bprotonpx,bprotonpy,bprotonpz);
            TVector3 Pd(bdpx,bdpy,bdpz);
            float pion_pt = sqrt(bpionpx*bpionpx+bpionpy*bpionpy);
            float pion_p  = Ppion.Mag();
            float p_pt = sqrt(bprotonpx*bprotonpx+bprotonpy*bprotonpy);
            float p_p  = Pp.Mag();
            float d_pt = sqrt(bdpx*bdpx+bdpy*bdpy);
            float d_p  = Pd.Mag();
            float Epi  = sqrt(Mpion*Mpion+pion_p*pion_p); 
            float Ep  = sqrt(Mp*Mp+p_p*p_p); 
            float Mppi=sqrt((Ep+Epi)*(Ep+Epi)-(Pp+Ppion).Mag2());

            //fixed cuts
            int c2=chi2primary_proton>3; //satisfied by default
            int c3=chi2primary_pi>3; //satisfied by default
            int c4=chi2primary_d>3; //satisfied by defaut
            int c5=dca_proton<5;
            //    int c5="ht_chi2topo<2 && ht_chi2topo>0";
            //    int c6="ht_chi2ndf<5";
            //    int c7="ht_l>4";
            int c8=ht_ldl>3; //satisfied by default
            int c9=bparticlemass>2.95 && bparticlemass<3.05;
            //    int c10="dca_pi>3&&dca_proton>1.&&dca_proton<5&&dca_deuteron<2";
            int c11=pion_pt>0.1; //not used now
            int c12=p_p<2.0;
            int c13=d_p<3.0;
            int c14=Mppi>(1.115+1*0.0015) || Mppi<(1.115-1*0.0015);

            if(iEntry%100000==0){
                cout<<"Processing entry:" <<iEntry<<endl;
            }
            //----- choose particle to reconstruct;
            //        if(bparticleid!=pid) continue;

            if(ht_chi2topo<0) continue;
            if(!(c5*c9*c12*c13*c14)) continue;

            //----- Fill histograms

            for(int ipiondca=0; ipiondca<ncutdcapionbins; ipiondca++){
                for(int iprodca=0; iprodca<ncutdcaprobins; iprodca++){
                    for(int idedca=0; idedca<ncutdcadebins; idedca++){
                        for(int idecayL=0; idecayL<ncutdecaylengthbins; idecayL++){
                            for(int ichi2topo=0; ichi2topo<ncutchi2topo; ichi2topo++){
                                for(int ichi2ndf=0; ichi2ndf<ncutchi2ndf; ichi2ndf++){
                                    if(dca_pi < dcaPionCuts[ipiondca]) continue;
                                    if(dca_proton< dcaProCuts[iprodca]) continue;
                                    if(dca_deuteron > dcaDeCuts[idedca]) continue;
                                    if(ht_l  < decayLCuts[idecayL]) continue;
                                    if(ht_chi2topo > chi2topoCuts[ichi2topo]) continue;
                                    if(ht_chi2ndf> chi2ndfCuts[ichi2ndf]) continue;
                                    sHtmass[ipiondca][iprodca][idedca][idecayL][ichi2topo][ichi2ndf]->Fill(bparticlemass);
                                    /*
                                    if(ipiondca!=cutIndex[0]) continue;
                                    if(iprodca!=cutIndex[1]) continue;
                                    if(idedca!=cutIndex[2]) continue;
                                    if(idecayL!=cutIndex[3]) continue;
                                    if(ichi2topo!=cutIndex[4]) continue;
                                    if(ichi2ndf!=cutIndex[5]) continue;
                                    sHtmassPRot[ipiondca][iprodca][idedca][idecayL][ichi2topo][ichi2ndf]->Fill(bparticlemass);
                                    */
                                }
                            }
                        }
                    }
                }
            }
        }
        //----- Write histograms

        for(int ipiondca=0; ipiondca<ncutdcapionbins; ipiondca++){
            for(int iprodca=0; iprodca<ncutdcaprobins; iprodca++){
                for(int idedca=0; idedca<ncutdcadebins; idedca++){
                    for(int idecayL=0; idecayL<ncutdecaylengthbins; idecayL++){
                        for(int ichi2topo=0; ichi2topo<ncutchi2topo; ichi2topo++){
                            for(int ichi2ndf=0; ichi2ndf<ncutchi2ndf; ichi2ndf++){
                                sHtmass[ipiondca][iprodca][idedca][idecayL][ichi2topo][ichi2ndf]-> Write();
                                    /*
                                    if(ipiondca!=cutIndex[0]) continue;
                                    if(iprodca!=cutIndex[1]) continue;
                                    if(idedca!=cutIndex[2]) continue;
                                    if(idecayL!=cutIndex[3]) continue;
                                    if(ichi2topo!=cutIndex[4]) continue;
                                    if(ichi2ndf!=cutIndex[5]) continue;
                                    sHtmassPRot[ipiondca][iprodca][idedca][idecayL][ichi2topo][ichi2ndf]->Fill(bparticlemass);
                                    */
                            }
                        }
                    } 
                }
            }
        }
    }

    //timer.Stop();
    //double rtime = timer.RealTime();
    //double ctime = timer.CpuTime();
    //cout<<"Real time is: "<<rtime<<" s, CPU time is: "<<ctime
    //    <<"s"<<endl;
    cout<<"Scan_P: "<<Scan_P << endl;
    if(Scan_P){
        for(int ipP=0; ipP<npP; ipP++){
            for(int idP=0; idP<ndP; idP++){
                sHtmass_scanP[ipP][idP]= new TH1F(Form("sHtmass_pP%d_dP%d",ipP,idP),Form("sHtmass_pP%d_dP%d",ipP,idP),nbin,massmin,massmax);
            }
        }

        Long64_t n_signal_Entries = shtriton_tree.GetEntries();
        cout<<"--------------------------------------------"<<endl;
        cout<<"Total entries for signal "<<n_signal_Entries<< endl;
        for (Long64_t iEntry = 0; iEntry<=n_signal_Entries; iEntry++){
            shtriton_tree.GetEntry(iEntry);
            const double Mpion=0.1396;
            const double Mp=0.938;
            TVector3 Ppion(bpionpx,bpionpy,bpionpz);
            TVector3 Pp(bprotonpx,bprotonpy,bprotonpz);
            TVector3 Pd(bdpx,bdpy,bdpz);
            float pion_pt = sqrt(bpionpx*bpionpx+bpionpy*bpionpy);
            float pion_p  = Ppion.Mag();
            float p_pt = sqrt(bprotonpx*bprotonpx+bprotonpy*bprotonpy);
            float p_p  = Pp.Mag();
            float d_pt = sqrt(bdpx*bdpx+bdpy*bdpy);
            float d_p  = Pd.Mag();
            float Epi  = sqrt(Mpion*Mpion+pion_p*pion_p); 
            float Ep  = sqrt(Mp*Mp+p_p*p_p); 
            float Mppi=sqrt((Ep+Epi)*(Ep+Epi)-(Pp+Ppion).Mag2());

            //fixed cuts
            int c1=bzdeuteron<0.2 && bzdeuteron >0;
            int c2=chi2primary_proton>3;
            int c3=chi2primary_pi>3;
            int c4=chi2primary_d>3;
            int c5=dca_proton<5&&dca_proton>0.5;
            int c6=ht_chi2topo<7;
            int c7=ht_chi2ndf<4;
            int c8=ht_l>6;
            int c9=ht_ldl>3;
            int c10=bparticlemass>2.95 && bparticlemass<3.05;
            int c11=dca_pi>4&&dca_deuteron<1;
            int c14=Mppi>(1.115+1*0.0015) || Mppi<(1.115-1*0.0015);
            //int c14=sqrt(0.1396*0.1396+0.938*0.938+2*sqrt((bprotonpx*bprotonpx+bprotonpy*bprotonpy+bprotonpz*bprotonpz+0.938*0.938)*(bpionpx*bpionpx+bpionpy*bpionpy+bpionpz*bpionpz+0.1396*0.1396))-2*bprotonpx*bpionpx-2*bprotonpy*bpionpy-2*bprotonpz*bpionpz)>(1.115+1*0.0015) || sqrt(0.1396*0.1396+0.938*0.938+2*sqrt((bprotonpx*bprotonpx+bprotonpy*bprotonpy+bprotonpz*bprotonpz+0.938*0.938)*(bpionpx*bpionpx+bpionpy*bpionpy+bpionpz*bpionpz+0.1396*0.1396))-2*bprotonpx*bpionpx-2*bprotonpy*bpionpy-2*bprotonpz*bpionpz)<(1.115-1*0.0015);

            if(iEntry%100000==0){
                cout<<"Processing entry:" <<iEntry<<endl;
            }
            //----- choose particle to reconstruct;
            //        if(bparticleid!=pid) continue;

            if(ht_chi2topo<0) continue;
            if(!(c2*c3*c4*c5*c6*c7*c8*c9*c10*c11*c14)) continue;

            //----- Fill histograms
            for(int ipP=0; ipP<npP; ipP++){
                for(int idP=0; idP<ndP; idP++){
                    if(p_p>pPCuts[ipP]) continue;
                    if(d_p>dPCuts[idP]) continue;
                    sHtmass_scanP[ipP][idP]->Fill(bparticlemass);
                }
            }

        }
        //----- Write histograms
        for(int ipP=0; ipP<npP; ipP++){
            for(int idP=0; idP<ndP; idP++){
                cout<<"for nP= "<<ipP<<", nd= "<<idP<<"; entries is: "<<sHtmass_scanP[ipP][idP]->GetEntries()<<endl;
                sHtmass_scanP[ipP][idP]->Write();
            }
        }
    }

    out_file_s->Close();


}