{
  TString inPidFile  = "pid_complete.root";    // this file contains the PndPidCandidates
  TString inRecoFile  = "reco_complete.root";
  TString inSimFile = "sim_complete.root";  // this file contains the MC truth
  TString inParFile = "simparams.root";
  std::string nomeout=Form("mcmatch_kinvtx_after.root");
  
  TFile *output = new TFile(nomeout.c_str(),"RECREATE");
  
  TString OutFile="output.root";  
  FairRunAna* fRun = new FairRunAna();
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  fRun->SetInputFile(inSimFile);
  fRun->AddFriend(inPidFile);
  fRun->AddFriend(inRecoFile);
  FairParRootFileIo* parIO = new FairParRootFileIo();
  parIO->open(inParFile);
  rtdb->setFirstInput(parIO);
  rtdb->setOutput(parIO);   
  fRun->SetOutputFile(OutFile);
  fRun->Init(); 
  PndMcTruthMatch mcm;
  PndAnalysis *theAnalysis = new PndAnalysis("SttMvdGemGenTrack","FtsIdealGenTrack");


  ///vertex fit
  TH1F *x_km = new TH1F("x_km","x_km",200,-0.005,0.005);
  TH1F *y_km = new TH1F("y_km","y_km",200,-0.005,0.005);
  TH1F *z_km = new TH1F("z_km","z_km",500,-0.5,0.5);
  //TH1F *x_km = new TH1F("x_km","x_km",200,-0.5,0.5);
  //TH1F *y_km = new TH1F("y_km","y_km",200,-0.5,0.5);
  //TH1F *z_km = new TH1F("z_km","z_km",500,-5,5);
  //TH1F *x_km = new TH1F("x_km","x_km",200,-0.1,0.1);
  //TH1F *y_km = new TH1F("y_km","y_km",200,-0.1,0.1);
  //TH1F *z_km = new TH1F("z_km","z_km",500,-0.2,0.2);

  //isto for pocaquality
  TH1F *chi2_vertex_d0  = new TH1F("chi2_vertex_d0","chi vertex",200,0,5);

  TNtuple *d0ntu = new TNtuple("d0ntu","ntu","indicemc:xk:yk:zk");

 
  RhoCandList verita ;
  //Starting cand list
  RhoCandList  pip1,pim1,pip2,pim2,kp,km;
  //cand list after mc match
  RhoCandList recokp,recopim1,recokm,recopip1,recopim2,recopip2;
  //cand combine
  RhoCandList d0;
  RhoCandidate reco=0;
  int nevts=0;
  int i=0,j=0;
  int ok1,ok2,ok3,ok4,okd0,ok5,ok6;
  int flag;
  if (nevts==0) nevts= theAnalysis->GetEntries();
  Bool_t fVerbose = kFALSE;
  
  int indicemc;

  while (theAnalysis->GetEvent() && i++<nevts)
    {
      flag=9000;
      //flag out of range
      //reco ad0
      ok1=5000;	ok2=5000;
      //reco d0
      ok3=5000;	ok4=5000;	okd0=5000;
      ok5=5000; ok6=5000;
      //clean the start list
      kp.Cleanup();
      km.Cleanup();
      pip1.Cleanup();
      pim1.Cleanup();
      pip2.Cleanup();
      pim2.Cleanup();

      //clean the matched list
      recokp.Cleanup();
      recopim1.Cleanup();
      recokm.Cleanup();
      recopip1.Cleanup();
      recopip2.Cleanup();
      recopim2.Cleanup();

    
      d0.Cleanup();
      verita.Cleanup();

      indicemc=-5000;
      if ((i%100)==0) cout<<"evt " << i << endl;

	  
      theAnalysis->FillList(kp, "KaonAllPlus");
      theAnalysis->FillList(km, "KaonAllMinus");
      theAnalysis->FillList(pip1, "PionAllPlus");
      theAnalysis->FillList(pim1, "PionAllMinus");
      theAnalysis->FillList(pip2, "PionAllPlus");
      theAnalysis->FillList(pim2, "PionAllMinus");
      theAnalysis->FillList(verita, "McTruth");

      if(kp.GetLength()!=0)
	{
	  for (j=0;j<kp.GetLength();++j) 
	    {
	      PndPidCandidate mic =(PndPidCandidate&) kp[j].GetRecoCandidate();
	      if (mic.GetMcIndex()==4) 
		{		         
		  recokp.Put(kp[j]);ok1=1;
		}
	    }
	}//end of if(kp.GetLength()!=0)

      if(pim1.GetLength()!=0)
	{
	  for (j=0;j<pim1.GetLength();++j) 
	    {
	      PndPidCandidate mic2 =(PndPidCandidate&) pim1[j].GetRecoCandidate();
	      if (mic2.GetMcIndex()==5) 
		{		         
		  recopim1.Put(pim1[j]);ok2=1;
		}
	    }
	}//end of if(pim1.GetLength()!=0)

      if(km.GetLength()!=0)
	{
	  for (j=0;j<km.GetLength();++j) 
	    {
	      PndPidCandidate mic =(PndPidCandidate&) km[j].GetRecoCandidate();
	      if (mic.GetMcIndex()==1) 
		{		         
		  recokm.Put(km[j]);ok3=1;
		}
	    }
	}//end of  if(km.GetLength()!=0)

      if(pip1.GetLength()!=0)
	{
	  for (j=0;j<pip1.GetLength();++j) 
	    {
	      PndPidCandidate mic2 =(PndPidCandidate&) pip1[j].GetRecoCandidate();
	      if (mic2.GetMcIndex()==2) 
		{		         
		  recopip1.Put(pip1[j]);ok4=1;
		}
	    }
	}//end of if(pip1.GetLength()!=0)

      if(pip2.GetLength()!=0)
	{
	  for (j=0;j<pip2.GetLength();++j) 
	    {
	      PndPidCandidate mic2 =(PndPidCandidate&) pip2[j].GetRecoCandidate();
	      if (mic2.GetMcIndex()==0) 
		{		         
		  recopip2.Put(pip2[j]);ok5=1;
		}
	    }
	}//end of if(pip2.GetLength()!=0)

      if(pim2.GetLength()!=0)
	{
	  for (j=0;j<pim2.GetLength();++j) 
	  {
	    PndPidCandidate mic2 =(PndPidCandidate&) pim2[j].GetRecoCandidate();
	    if (mic2.GetMcIndex()==3) 
	      {		         
		recopim2.Put(pim2[j]);ok6=1;
	      }
	  }
	}//end of if(pim2.GetLength()!=0)
      
      
    if(ok3==ok4==1)
	{
	  d0.Combine(recokm,recopip1);
	  okd0=1;
	}

      if(d0.GetLength()!=0 && okd0==1)
	{

	  for(j=0;j<d0.GetLength();++j) 
	    {
	      TVector3 vtx;
	      // determine pseudo vertex
	      Float_t chi2 = -1;
	      if (kTRUE)
		{
		  PndKinVtxFitter fitt(d0[j]);
		  fitt.Fit();
		  chi2 = fitt->GetChi2();
		  RhoCandidate *d0fit=d0[j].GetFit();
		  vtx=d0fit->GetPosition();
		}
	      if (kFALSE)
		{
		  PndVtxPoca poca;
		  chi2 = poca.GetPocaVtx(vtx,d0[j]);
		}
	      if (kFALSE)
		{
		  PndVtxPRG prg(d0[j]);
		  TMatrixD vtxcov(3,3);
		  prg.SetNIterations(10);
		  chi2 = prg.FitVertexFast(vtx,vtxcov);
		}
	      RhoCandidate *kaone=d0[j].Daughter(0);
	      Int_t mcIdx = kaone->GetMcIdx();
 	      if (mcIdx>=0)// && mcIdx<verita.GetLength())     
		{
		  RhoCandidate truth = verita[mcIdx]; 
		  //RhoCandidate *truth = kaone->GetMcTruth();   // DOES NOT WORK!
		  if ( truth.PdgCode() == kaone.PdgCode() )  
		    {
		      TVector3 mc_ver= truth.GetPosition(); 
		      if (fVerbose) cout << "reco: x " << vtx.X() << " y " << vtx.Y() << " z " << vtx.Z() << endl;
		      if (fVerbose) cout << "mc:   x " << mc_ver.X() << " y " << mc_ver.Y() << " z " << mc_ver.Z() << endl;
		      x_km->Fill(vtx.X()-mc_ver.X());
		      y_km->Fill(vtx.Y()-mc_ver.Y());
		      z_km->Fill(vtx.Z()-mc_ver.Z());
		      chi2_vertex_d0->Fill(chi2);
		      indicemc=i-1;
		      d0ntu->Fill(indicemc,vtx.X(),vtx.Y(),vtx.Z());
		    }
		}
	    }
	}//end of if(d0.GetLength()!=0 && okd0==1)

    }
   
 


  TCanvas *c4 = new TCanvas("c4","",1200,1000);
  c4->Divide(2,2);
  c4->cd(1);x_km->Draw();
  c4->cd(2);y_km->Draw();
  c4->cd(3);z_km->Draw();
  c4->cd(4);chi2_vertex_d0->Draw();


  output->cd();
  x_km->Write();
  y_km->Write();
  z_km->Write();
  chi2_vertex_d0->Write();
  d0ntu->Write();

  output->Write();
  output->Close();


    
}