class RhoCandList;
class RhoCandidate;
class TFitParams;
class PndAnaPidSelector;
class PndAnaPidCombiner;
class PndAnalysis;

#include <sstream>



void s_ana_complet(Int_t start, Double_t pbarmom)
{
	// *** some variables
	int i=0,j=0, k=0, l=0;
	int nevts=0;
    gStyle->SetOptFit(1011);
	
	 TString inSimFile =Form("sim_complete_%d_sim.root",start);
	 TString input = Form("sim_complete_%d.root",start);
	
	    TString  parAsciiFile   = "all.par";
         // TString  prefix         = "evtday1";
           //TString  input          = "psi2s_Jpsi2pi_Jpsi_mumu.dec";
              TString  output         = "fana";
              TString  friend1        = "pid";
              TString  friend2        = "";
              TString  friend3        = "";
              TString  friend4        = "";


	 
	 PndMasterRunAna *fRun= new PndMasterRunAna();
  fRun->SetInput(input);
  fRun->SetOutput(output);
  fRun->SetFriend1(friend1);
  fRun->SetFriend2(friend2);
  fRun->SetFriend3(friend3);
  fRun->SetFriend4(friend4);
  fRun->SetParamAsciiFile(parAsciiFile);
  fRun->Setup();
	 
	    fRun->Init();
	 
	 
	 
	 
	 
	TFile *out = TFile::Open(Form("sim_complete_%d_ana.root",start),"RECREATE");
	


	// *** create ntuple
    
    ///////##########################################
    //create ntuple for the reconstructed events
	
  
    
  
    
    TNtuple *pipi2 = new
   
TNtuple("pipi2","pipi2","npair:p1:theta1:phi1:p2:theta2:phi2:emc_ene1:emc_ene2");
    
   
    
     
    TNtuple *pip_mc = new TNtuple("pip_mc","pip_mc","pdgp:rmom1:rthe1:rph1:rmom1cm:rthe1cm:rph1cm:rmom1_mc:rthe1_mc:rph1_mc:rmom1cm_mc:rthe1cm_mc:rph1cm_mc");
    TNtuple *pim_mc = new TNtuple("pim_mc","pim_mc","pdgm:rmom2:rthe2:rph2:rmom2cm:rthe2cm:rph2cm:rmom2_mc:rthe2_mc:rph2_mc:rmom2cm_mc:rthe2cm_mc:rph2cm_mc");
	
	


    TNtuple *hpair = new TNtuple("hpair","hpair","npair");
    TNtuple *hpt2 = new TNtuple("hpt2","hpt2","valuefit");
   
    

    

    //######################################################################
	
		
	// *** the data reader object
	PndAnalysis* theAnalysis = new PndAnalysis();
	if (nevts==0) nevts= theAnalysis->GetEntries();
	
	std::cout<<"nevts="<<nevts<<std::endl;
	
    
    
    
     RhoCandList  pipluse, piminuse;
    
    
	
    
  
    
	// *** Pid Selection Algorithms
//	TString pidSelection = "PidAlgoEmcBayes;PidAlgoDrc;PidAlgoStt;PidAlgoMvd";
   TString pidSelection = "PidAlgoEmcBayes;PidAlgoDrc;PidAlgoStt;PidAlgoMdtHardCuts;PidAlgoMvd";

    
    
    // *** the lorentz vector of the initial
    Double_t melectron=0.0005109989;
    Double_t mpion=0.139570;
    Double_t ppbar=pbarmom;
    Double_t mproton=0.9382720;
    Double_t Epbar=sqrt(ppbar*ppbar+mproton*mproton);
    TLorentzVector lPbarBeam;
    TLorentzVector  lPTarget;
    lPTarget.SetPxPyPzE( 0.,0.,0.,mproton);
    lPbarBeam.SetPxPyPzE(0.,0.,ppbar,Epbar);
    TLorentzVector ini =  lPbarBeam + lPTarget;
    
    
    Int_t i;
    
    
    Double_t ncountss=0;
	
	while (theAnalysis->GetEvent() && i++<nevts)
	{
		if ((i%100)==0) cout<<"evt " << i << endl;
				

        
        
        
        TLorentzVector lorpip1, lorpim1,lorpipcm1,lorpimcm1 ;
        
        
        
        //define the parameter to choose the best eplus-eminus pair per event
        
        Int_t ix=-1, iy=-1;
        
        Int_t npion_pair = 0;
        Double_t besttheta=3.14;
        
        
    
     
		
	
        theAnalysis->FillList(pipluse, "ElectronAllPlus",pidSelection);
        theAnalysis->FillList(piminuse, "ElectronAllMinus",pidSelection);
        
     
        
        //loop on the possible pairs per event and choose the best back to back in the center of mass frame.
		
        for (Int_t nc1=0;nc1<pipluse.GetLength();++nc1)
        {
            
         
	    lorpip1=pipluse[nc1]->P4();
	    lorpipcm1=pipluse[nc1]->P4();
	    lorpipcm1.Boost(-(ini).BoostVector());
	    

            for (Int_t nc2=0;nc2<piminuse.GetLength();++nc2)
            {
                
		 lorpim1=piminuse[nc2]->P4();
	        
		
		 lorpimcm1=piminuse[nc2]->P4();
	         lorpimcm1.Boost(-(ini).BoostVector());
		
		
                Double_t th1= lorpipcm1.Theta();
                Double_t th2= lorpimcm1.Theta();
                
                Double_t dth=fabs(3.14-th1-th2);
                
                
                npion_pair++;
                
                if(dth<besttheta)
                {
                    besttheta=dth; 
                    ix=nc1; 
                    iy=nc2;
                }
            }
        }
        
      
       
        // the best positron is at [ix], and the best electron is at [iy]
        
        hpair->Fill((double) npion_pair);
        hpt2->Fill((double) besttheta);
        
        
       // if(pipluse2.GetLength()==0||piminuse2.GetLength()==0) continue;        
        
	if(pipluse.GetLength()!=1||piminuse.GetLength()!=1) continue; 
        
        
        TLorentzVector lorpip=pipluse[ix]->P4();
        TLorentzVector lorpipcm= lorpip;     
        lorpipcm.Boost(-(ini).BoostVector());     
		 
	 
   
   
    TLorentzVector lorpim=piminuse[iy]->P4();
    TLorentzVector lorpimcm= lorpim;     
    lorpimcm.Boost(-(ini).BoostVector());     
		 	 



PndPidCandidate *micp=(PndPidCandidate*)pipluse[ix]->GetRecoCandidate();
PndPidCandidate *micm=(PndPidCandidate*)piminuse[iy]->GetRecoCandidate();

 //std::cout<<" EMC  "<<micp->GetEmcRawEnergy()<<std::endl;

        //Fill reconstructed variables in pipi2

Float_t pipi2_ntuple[] = {
                          npion_pair,
                            lorpip.P(), lorpip.Theta(), lorpip.Phi(),
                            lorpim.P(), lorpim.Theta(), lorpim.Phi(),micp->GetEmcRawEnergy(), micm->GetEmcRawEnergy()} ;
			    

                       pipi2->Fill(pipi2_ntuple);
			    
// get truth information 

bool mctp = theAnalysis->McTruthMatch(pipluse[ix]);
RhoCandidate *true_piplus = pipluse[ix]->GetMcTruth();
			
bool mctm = theAnalysis->McTruthMatch(piminuse[iy]);
RhoCandidate *true_piminus = piminuse[iy]->GetMcTruth();

      

// Reconstruction -MC Truth events


       
       if (true_piplus && true_piplus->PdgCode()==pipluse[ix]->PdgCode())
       {
    
            int pdgmcp=pipluse[ix]->PdgCode();

       TLorentzVector lorpip_mc=true_piplus->P4(); 
       TLorentzVector lorpipcm_mc= lorpip_mc; 
       lorpipcm_mc.Boost(-(ini).BoostVector());
       
       
       Float_t pip_mc_ntuple[] = {pdgmcp, lorpip.P(), lorpip.Theta(), lorpip.Phi(),lorpipcm.P(), lorpipcm.Theta(),lorpipcm.Phi(),lorpip_mc.P(),    lorpip_mc.Theta(), lorpip_mc.Phi(), lorpipcm_mc.P(),lorpipcm_mc.Theta(),lorpipcm_mc.Phi()};
	 pip_mc->Fill(pip_mc_ntuple);
       
       }

        
               
        
	if (true_piminus && true_piminus->PdgCode()==piminuse[iy]->PdgCode())
	{
	
	

            int pdgmcm=piminuse[iy]->PdgCode();
	    
            
            TLorentzVector lorpim_mc=true_piminus->P4();
	    TLorentzVector lorpimcm_mc= lorpim_mc;
            lorpimcm_mc.Boost(-(ini).BoostVector());
            
            
            Float_t pim_mc_ntuple[] = {pdgmcm, lorpim.P(), lorpim.Theta(), lorpim.Phi(),lorpimcm.P(), lorpimcm.Theta(),lorpimcm.Phi(),lorpim_mc.P(),    lorpim_mc.Theta(), lorpim_mc.Phi(), lorpimcm_mc.P(),lorpimcm_mc.Theta(),lorpimcm_mc.Phi()};
            pim_mc->Fill(pim_mc_ntuple);
        }
        
        
	
	

	   }
    
     

    

	// *** write out all the histos
	out->cd();
        pipi2->Write();
	
	pip_mc->Write();
	pim_mc->Write();
        hpair->Write();
        hpt2->Write();
      
	out->Save();
    
    
    
	cout << " Test passed" << endl;
        cout << " All ok " << endl;
    
       exit(0);

	
}