match_kinvtx(int nevts=1000, int whichfit=0, bool fVerbose=false)
{
    //int nevts=1000; //10k
	enum kFitType {kPoca, kFastPrg, kFullPrg, kKinVtx, kChiVtx};
	const char* kFitNames[5] = {"Poca", "FastPrg", "FullPrg", "KinVtx", "ChiVtx"};
	//int whichfit = kFullPrg;
    //Bool_t fVerbose = false;

    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");
    //theAnalysis->SetVerbose(false);

    ///vertex fit
    TH1F *x_km = new TH1F("x_km","x_km",200,-0.25,0.25);
    TH1F *y_km = new TH1F("y_km","y_km",200,-0.25,0.25);
    TH1F *z_km = new TH1F("z_km","z_km",500,-1.,1.);
    //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 ,charged,all;
    //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;
    RhoCandidate* truth=0;
    int i=0,j=0;
    int ok1,ok2,ok3,ok4,okd0,ok5,ok6;
    int flag;
    if (nevts==0) nevts= theAnalysis->GetEntries();

    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%50)==0)
          std::cout<<"evt " << i << std::endl;

/*        theAnalysis->FillList(all,"All");
       // std::cout<<"ALL:      "<<all<<std::endl;
        for (j=0;j<all.GetLength();++j)
        {
          RhoCandidate* onecand = all[j];
          std::cout<<"pinter test: [] vs. Get: "<<onecand<<" vs. "<<all.Get(j)<<std::endl;
          std::cout<<"all:     "<<*onecand<<std::endl;
          RhoCandidate* onetruth=onecand->GetMcTruth();
          std::cout<<"truth pointer test: "<<onetruth<<std::endl;
          std::cout<<"truth:   "<<*onetruth<<std::endl;
          
        }*/
        
/*        theAnalysis->FillList(charged,"Charged");
        for (j=0;j<charged.GetLength();++j)
        {
          RhoCandidate* onecharged = charged[j];
          std::cout<<"charged: "<<*onecharged<<std::endl;
          std::cout<<"a#"<<std::endl;
          RhoCandidate* truth=onecharged->GetMcTruth();
          std::cout<<"b#"<<std::endl;
          std::cout<<"truth:   "<<*truth<<std::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(fVerbose) printf("List size check: kp:%i \tkm:%i \tpip1:%i \tpip2:%i \tpim1:%i \tpim2:%i\n",kp.GetLength(),km.GetLength(),pip1.GetLength(),pip2.GetLength(),pim1.GetLength(),pim2.GetLength());
        if (kp.GetLength()!=0)
        {
            for (j=0;j<kp.GetLength();++j)
            {
              //std::cout<<"kp:  "<<kp[j]<<std::endl;
              //std::cout<<"a"<<std::endl;
                //PndPidCandidate mic =(PndPidCandidate&) kp[j]->GetRecoCandidate();
                truth=kp[j]->GetMcTruth();
              //std::cout<<"b"<<std::endl;
                if(!truth) continue;
              //std::cout<<"c"<<std::endl;
                if(truth->PdgCode()==321)
                    //if (mic.GetMcIndex()==4)
                {
                    //std::cout<<"4\t";
                    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)
                truth=pim1[j]->GetMcTruth();
                if(!truth) continue;
                if(truth->PdgCode()==211)
                {
                    //std::cout<<"5\t";
                    recopim1.Put(pim1[j]);
                    ok2=1;
                }
                //std::cout<<"pim1[j]->GetMcTruth()->PdgCode()\t="<<pim1[j]->GetMcTruth()->PdgCode()<<std::endl;
            }
        }//end of if(pim1.GetLength()!=0)

        if (km.GetLength()!=0)
        {
            for (j=0;j<km.GetLength();++j)
            {
                //PndPidCandidate mic =(PndPidCandidate&) km[j]->GetRecoCandidate();
                truth=km[j]->GetMcTruth();
                if(!truth) continue;
                if (truth->PdgCode()==-321)
                    //if (mic.GetMcIndex()==1)
                {
                   // std::cout<<"1\t";
                    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 (pip1[j]->GetMcTruth()->PdgCode()==211)
                    //if (mic2.GetMcIndex()==2)
                truth=pip1[j]->GetMcTruth();
                if(!truth) continue;
                if(truth->PdgCode()==211)
                {
                   // std::cout<<"2\t";
                    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 (pip2[j]->GetMcTruth()->PdgCode()==211)
                    //if (mic2.GetMcIndex()==0)
                truth=pip2[j]->GetMcTruth();
                if(!truth) continue;
                if(truth->PdgCode()==211)
                {
                   // std::cout<<"0\t";
                    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 (pim2[j]->GetMcTruth()->PdgCode()==-211)
                    //if (mic2.GetMcIndex()==3)
                truth=pim2[j]->GetMcTruth();
                if(!truth) continue;
                if(truth->PdgCode()==211)
                {
                   // std::cout<<"3\t";
                    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;
				switch(whichfit){
                case kChiVtx:
                    PndChiVtxFitter fitb(d0[j]);
                    //fitb.SetVerbose(kTRUE);
                    fitb.Fit();
                    chi2 = fitb.GetChi2();
                    RhoCandidate *d0fit=d0[j]->GetFit();
                    RhoCandidate *pionfit=d0fit->Daughter(0);
                    vtx=pionfit->GetPosition();
                case kKinVtx:
                    PndKinVtxFitter fitt(d0[j]);
                    //fitt.SetVerbose(kTRUE);
                    fitt.Fit();
                    chi2 = fitt.GetChi2();
                    RhoCandidate *d0fit=d0[j]->GetFit();
                    RhoCandidate *pionfit=d0fit->Daughter(0);
                    vtx=pionfit->GetPosition();
					break;
                case kPoca:
                    PndVtxPoca poca;
                    chi2 = poca.GetPocaVtx(vtx,d0[j]);
					break;
                case kFastPrg:
                    PndVtxPRG prg(d0[j]);
                    TMatrixD vtxcov(3,3);
                    chi2 = prg.FitVertexFast(vtx,vtxcov,true);//skipping COV
                    prg.SetExpansionPoint(vtx);
                    chi2 = prg.FitVertexFast(vtx,vtxcov);
					break;
                case kFullPrg:
                    PndVtxPRG prg(d0[j]);
                    //prg.SetDebug(true);
                    //prg.SetVerbose();
                    prg.SetNIterations(3);
                    prg.Fit();
 					chi2 = prg.GetChi2();
                    RhoCandidate *d0fit=d0[j]->GetFit();
                    RhoCandidate *pionfit=d0fit->Daughter(0);
                    vtx=pionfit->GetPosition();
					break;
				}
	                RhoCandidate *kaone=d0[j]->Daughter(0);
                // Int_t mcIdx = kaone->GetMcIdx();
                //if (mcIdx>=0)// && mcIdx<verita.GetLength())
                RhoCandidate *truth = kaone->GetMcTruth();
                //std::cout<<"truth pointer: "<<truth<<std::endl;
                if (truth)
                {
                    //RhoCandidate truth = verita[mcIdx];
                    if ( truth->PdgCode() == kaone->PdgCode() )
                    {
                        TVector3 mc_ver= truth->GetPosition();
                        if (fVerbose) std::cout << "reco: x " << vtx.X() << " y " << vtx.Y() << " z " << vtx.Z() << std::endl;
                        if (fVerbose) std::cout << "mc:   x " << mc_ver.X() << " y " << mc_ver.Y() << " z " << mc_ver.Z() << std::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->Fit("gaus");
    x_km->Draw();
    c4->cd(2);
    //y_km->Fit("gaus");
    y_km->Draw();
    c4->cd(3);
    //z_km->Fit("gaus");
    z_km->Draw();
    c4->cd(4);
    chi2_vertex_d0->Draw();
	c4->Print(Form("vtxtest_%s.pdf",kFitNames[whichfit]));

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

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

}