void testParticles ( int nevts=0 )
{ 
  // *** some variables
  int i=0,j=0, k=0, l=0;
  TString plotfile="datakin/testParticles.root";
  TString OutFile="datakin/dummyOut.root";

  // *** the files coming from the simulation
  TString inPidFile  = "datakin/pid_complete.root";    // this file contains the PndPidCandidates
  TString inRecoFile = "datakin/reco_complete.root";
  TString inDigiFile = "datakin/digi_complete.root";
  TString inSimFile  = "datakin/sim_complete.root";  // this file contains the MC truth
  TString inParFile  = "datakin/simparams.root";

  gStyle->SetOptFit ( 1011 );
  gROOT->LoadMacro ( "$VMCWORKDIR/macro/run/Tools.C" );
  ImproveDefaultStyle();

  FairLogger::GetLogger()->SetLogToFile ( kFALSE );

  // *** initialization
  FairRunAna* fRun = new FairRunAna();
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  fRun->SetInputFile ( inSimFile );
  fRun->AddFriend ( inDigiFile );
  fRun->AddFriend ( inRecoFile );
  fRun->AddFriend ( inPidFile );

  FairParRootFileIo* parIO = new FairParRootFileIo();
  parIO->open ( inParFile );
  rtdb->setFirstInput ( parIO );
  rtdb->setOutput ( parIO );

  fRun->SetOutputFile ( OutFile );
  FairGeane *Geane = new FairGeane();
  fRun->AddTask(Geane);
  fRun->Init();

  int nbins=100;
  double max=10.;
  double delta=max/20.;
  double dpull=15.;
  double dvmax=0.5;
  double epmax=0.1;
  double evmax=0.05;
  // *** create some histograms
  TList* histolist = new TList();
  
  // Px
  TH1F* h_dpx = new TH1F ( "h_dpx","Momentum difference x;#Delta p_{x} / GeV/c",nbins,-delta,delta);
  histolist->Add(h_dpx);h_dpx->SetLineColor(kRed+2);
  TH1F* h_ppx = new TH1F ( "h_ppx","Momentum pull x;#Delta p_{x} / #sigma(p_{x})",nbins,-dpull,dpull);
  histolist->Add(h_ppx);h_ppx->SetLineColor(kGreen+2);
  TH1F* h_epx = new TH1F ( "h_epx","Momentum errors x;#sigma(p_{x}) / GeV/c",nbins,0,epmax);
  histolist->Add(h_epx);h_epx->SetLineColor(kBlue+2);
  
  // Py
  TH1F* h_dpy = new TH1F ( "h_dpy","Momentum difference y;#Delta p_{y} / GeV/c",nbins,-delta,delta);
  histolist->Add(h_dpy);h_dpy->SetLineColor(kRed+2);  
  TH1F* h_ppy = new TH1F ( "h_ppy","Momentum pull y;#Delta p_{y} / #sigma(p_{y})",nbins,-dpull,dpull);
  histolist->Add(h_ppy);h_ppy->SetLineColor(kGreen+2);  
  TH1F* h_epy = new TH1F ( "h_epy","Momentum errors y;#sigma(p_{y}) / GeV/c",nbins,0,epmax);
  histolist->Add(h_epy);h_epy->SetLineColor(kBlue+2);
  
  // Pz
  TH1F* h_dpz = new TH1F ( "h_dpz","Momentum difference z;#Delta p_{z} / GeV/c",nbins,-delta,delta);
  histolist->Add(h_dpz);h_dpz->SetLineColor(kRed+2);
  TH1F* h_ppz = new TH1F ( "h_ppz","Momentum pull z;#Delta p_{z} / #sigma(p_{z})",nbins,-dpull,dpull);
  histolist->Add(h_ppz);h_ppz->SetLineColor(kGreen+2);
  TH1F* h_epz = new TH1F ( "h_epz","Momentum errors z;#sigma(p_{z}) / GeV/c",nbins,0,epmax);
  histolist->Add(h_epz);h_epz->SetLineColor(kBlue+2);
  
  // E
  TH1F* h_dE = new TH1F ( "h_dE","Energy difference;#Delta E / GeV",nbins,-delta,delta);
  histolist->Add(h_dE);h_dE->SetLineColor(kRed+2);
  TH1F* h_pE = new TH1F ( "h_pE","Energy pull;#Delta E / #sigma(E)",nbins,-dpull,dpull);
  histolist->Add(h_pE);h_pE->SetLineColor(kGreen+2);
  TH1F* h_eE = new TH1F ( "h_eE","Energy errors;#sigma(E) / GeV",nbins,0,epmax);
  histolist->Add(h_eE);h_eE->SetLineColor(kBlue+2);
  
  // X
  TH1F* h_dvx = new TH1F ( "h_dvx","Position difference x;#Delta x / cm",nbins,-dvmax,dvmax);
  histolist->Add(h_dvx);h_dvx->SetLineColor(kRed+2);
  TH1F* h_pvx = new TH1F ( "h_pvx","Position pull x;#Delta x / #sigma(x)",nbins,-dpull,dpull);
  histolist->Add(h_pvx);h_pvx->SetLineColor(kGreen+2);
  TH1F* h_evx = new TH1F ( "h_evx","Position errors x;#sigma x / cm",nbins,0,evmax);
  histolist->Add(h_evx);    h_evx->SetLineColor(kBlue+2);
  
  // Y
  TH1F* h_dvy = new TH1F ( "h_dvy","Position difference y;#Delta y / cm",nbins,-dvmax,dvmax);
  histolist->Add(h_dvy);h_dvy->SetLineColor(kRed+2);
  TH1F* h_pvy = new TH1F ( "h_pvy","Position pull y;#Delta y / #sigma(y)",nbins,-dpull,dpull);
  histolist->Add(h_pvy);h_pvy->SetLineColor(kGreen+2);
  TH1F* h_evy = new TH1F ( "h_evy","Position errors y;#sigma y / cm",nbins,0,evmax);
  histolist->Add(h_evy);h_evy->SetLineColor(kBlue+2);
  
  // Z
  TH1F* h_dvz = new TH1F ( "h_dvz","Position difference z;#Delta z / cm",nbins,-dvmax,dvmax);
  histolist->Add(h_dvz);h_dvz->SetLineColor(kRed+2);
  TH1F* h_pvz = new TH1F ( "h_pvz","Position pull z;#Delta z / #sigma(z)",nbins,-dpull,dpull);
  histolist->Add(h_pvz);h_pvz->SetLineColor(kGreen+2);
  TH1F* h_evz = new TH1F ( "h_evz","Position errors z;#sigma z / cm",nbins,0,evmax);
  histolist->Add(h_evz);h_evz->SetLineColor(kBlue+2);
  
  //Params
  //Qp
  double maxQp=12;
  TH1F* h_dQp = new TH1F ( "h_dQp","SD Paramter Qp difference;#Delta Qp * GeV/c",nbins,-maxQp,maxQp);
  histolist->Add(h_dQp);h_dQp->SetLineColor(kRed+2);
  TH1F* h_pQp = new TH1F ( "h_pQp","SD Paramter Qp pull;#Delta Qp / #sigma_(Qp)",nbins,-dpull,dpull);
  histolist->Add(h_pQp);h_pQp->SetLineColor(kGreen+2);
  TH1F* h_eQp = new TH1F ( "h_eQp","SD Paramter Qp errors;#sigma_(Qp) * GeV/c",nbins,0,maxQp);
  histolist->Add(h_eQp);h_eQp->SetLineColor(kBlue+2);h_eQp->SetOption("log");
  
  //V
  double maxV=5.5;
  TH1F* h_dV = new TH1F ( "h_dV","SD Paramter V difference;#Delta V / ?",nbins,-maxV,maxV);
  histolist->Add(h_dV);h_dV->SetLineColor(kRed+2);
  TH1F* h_pV = new TH1F ( "h_pV","SD Paramter V pull;#Delta V / #sigma_(V)",nbins,-dpull,dpull);
  histolist->Add(h_pV);h_pV->SetLineColor(kGreen+2);
  TH1F* h_eV = new TH1F ( "h_eV","SD Paramter V errors;#sigma_(V) / ?",nbins,0,0.01*maxV);
  histolist->Add(h_eV);h_eV->SetLineColor(kBlue+2);h_eV->SetOption("log");
  
  //W
  TH1F* h_dW = new TH1F ( "h_dW","SD Paramter W difference;#Delta W / ?",nbins,-maxV,maxV);
  histolist->Add(h_dW);h_dW->SetLineColor(kRed+2);
  TH1F* h_pW = new TH1F ( "h_pW","SD Paramter W pull;#Delta W / #sigma_(W)",nbins,-dpull,dpull);
  histolist->Add(h_pW);h_pW->SetLineColor(kGreen+2);
  TH1F* h_eW = new TH1F ( "h_eW","SD Paramter W errors;#sigma_(W) / ?",nbins,0,0.01*maxV);
  histolist->Add(h_eW);h_eW->SetLineColor(kBlue+2);h_eW->SetOption("log");
  
  //TV
  double maxT=2.1;
  TH1F* h_dTV = new TH1F ( "h_dTV","SD Paramter TV difference;#Delta TV",nbins,-maxT,maxT);
  histolist->Add(h_dTV);h_dTV->SetLineColor(kRed+2);
  TH1F* h_pTV = new TH1F ( "h_pTV","SD Paramter TV pull;#Delta TV / #sigma_(TV)",nbins,-dpull,dpull);
  histolist->Add(h_pTV);h_pTV->SetLineColor(kGreen+2);
  TH1F* h_eTV = new TH1F ( "h_eTV","SD Paramter TV errors;#sigma_(TV)",nbins,0,maxT);
  histolist->Add(h_eTV);h_eTV->SetLineColor(kBlue+2);h_eTV->SetOption("log");
  
  //TW
  TH1F* h_dTW = new TH1F ( "h_dTW","SD Paramter TW difference;#Delta TW",nbins,-maxT,maxT);
  histolist->Add(h_dTW);h_dTW->SetLineColor(kRed+2);
  TH1F* h_pTW = new TH1F ( "h_pTW","SD Paramter TW pull;#Delta TW / #sigma_(TW)",nbins,-dpull,dpull);
  histolist->Add(h_pTW);h_pTW->SetLineColor(kGreen+2);
  TH1F* h_eTW = new TH1F ( "h_eTW","SD Paramter TW errors;#sigma_(TW)",nbins,0,maxT);
  histolist->Add(h_eTW);h_eTW->SetLineColor(kBlue+2);h_eTW->SetOption("log");
  
  TH1F* h_chi2   = new TH1F ( "h_chi2","Track fit #chi^{2};#chi^{2}",nbins,0,50 );
  histolist->Add(h_chi2);//h_chi2->SetOption("log");

  TH1F* h_prob   = new TH1F ( "h_prob","Track fit p(#chi^{2});p",nbins,0,1 );
  histolist->Add(h_prob);h_prob->SetOption("log");

  TH1I* h_stat   = new TH1I( "h_stat","Track fit status;",30,-5,10 );
  //histolist->Add(h_stat);//h_stat->SetOption("log");

  TH1F* h_m = new TH1F ( "h_m","Mass;m / GeV/c^{2}",nbins,0.,max);
  //histolist->Add(h_m);
  TH1F* h_p = new TH1F ( "h_p","Momentum;p / GeV/c",nbins,0.,max);
  histolist->Add(h_p);

  // *** the data reader object
  TDatabasePDG* pdg = RhoPdtLoader::ReadPDGTable(Form("%s/pgenerators/EvtGen/evt.pdl",gSystem->Getenv("VMCWORKDIR")));                   
  PndAnalysis* theAnalysis = new PndAnalysis ( "SttMvdGemGenTrack" , "FtsIdealGenTrack" );
  
  // adding particles from our decay which are not present in the particle table
  if ( nevts==0 ) nevts= theAnalysis->GetEntries();

  // *** RhoCandLists for the analysis
  //RhoCandList all, chrg, el;
  RhoCandList mctrk;
  RhoCandList all;

  // *** Mass selector for the jpsi cands
  //RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector ( "jpsi",3.096,0.5 );
  //RhoGoodTrackSelector *selector = new RhoGoodTrackSelector("good tracks");
  //selector->SetFitCut(1e5,0.001); //(cut chisquare above, cut prob below)


  // *** the lorentz vector of the initial psi(2S)
  TLorentzVector ini ( 0, 0, 6.231552, 7.240065 );
  TVector3 dummy(0,0,0);

  // ***
  // the event loop
  // ***
  //TMemStat mm("gnubuiltin");
  while ( theAnalysis->GetEvent() && i++<nevts ) {
    //if(i<238)continue;
    if ( ( i%100 ) ==0 )
      cout<<"evt " << i << endl;

    // *** the MC Truth objects
    //cout<<"fill mc list"<<endl;
    //theAnalysis->FillList ( mctrk,"McTruth" );
    //cout<<"fill candidate list"<<endl;
    //theAnalysis->FillList ( all, "KaonTight" );
    theAnalysis->FillList ( all, "PionTight" );
    //theAnalysis->FillList ( all, "Charged" );
    //all.Select(selector);
    for ( j=0;j<all.GetLength();++j ) {
      //cout<<"particle j="<<j<<endl;
      RhoCandidate* cand = all[j];
      if(!cand) {
        cout<<"no cand?"<<endl;
        continue;
      }
      //truthMatcher.MctMatch(*cand,mctrk);
      FairRecoCandidate* recocand=cand->GetRecoCandidate();
      int fitstat=recocand->GetFitStatus();

      if(0>fitstat) {
        //cout<<"bad fitstat: "<<fitstat<<endl;
        continue;
      }

      double chiq=recocand->GetChiSquared();
      double ndf=recocand->GetDegreesOfFreedom();
      double prob=TMath::Prob(chiq,ndf);

    //  if(prob < 0.001) continue;
      h_chi2->Fill(chiq);
      h_prob->Fill(prob);
      h_stat->Fill(fitstat);

      theAnalysis->PropagateToIp(cand);
      //cout<<"cand"<<cand<<endl;
      TLorentzVector tlvc=cand->P4();
      TVector3 vtxc=cand->Pos();
      RhoError covp4=cand->P4Cov();
      RhoError covvt=cand->PosCov();
      //cout<<"mcitd="<<cand->GetMcIdx()<<" mccd "<<mccd<<endl;
      h_m->Fill(tlvc.M());
      h_p->Fill(tlvc.P());
      
      RhoCandidate* mccd = cand->GetMcTruth();
      if(!mccd) {
        cout<<"no mccand?"<<endl;
        continue;
      }
      TLorentzVector tlvm=mccd->P4();
      TVector3 vtxm=mccd->Pos();
      TVector3 momm=mccd->P3();

      FairTrackParP parabola = theAnalysis->GetFirstPar(cand);
//       if(false){
//         //theAnalysis->SetVerbose(5);
//         //cout<<"************ can pdgcode is "<<cand->PdgCode()<<endl;
//         //theAnalysis->SetVerbose(0);
//         tlvc.SetXYZM(parabola.GetPx(),parabola.GetPy(),parabola.GetPz(),tlvc.M());
//         vtxc.SetXYZ(parabola.GetX(),parabola.GetY(),parabola.GetZ());
//         for(int ie=0;ie<3;ie++){
//           for (int je=0;je<3;je++){
//             covvt[ie][je]=0.;
//             covp4[ie][je]=0.;
//           }
//           covp4[ie][3]=0.;
//           covp4[3][ie]=0.;
//         }
//         covvt[0][0]=parabola.GetDX()*parabola.GetDX();
//         covvt[1][1]=parabola.GetDY()*parabola.GetDY();
//         covvt[2][2]=parabola.GetDZ()*parabola.GetDZ();
//         covp4[0][0]=parabola.GetDPx()*parabola.GetDPx();
//         covp4[1][1]=parabola.GetDPy()*parabola.GetDPy();
//         covp4[2][2]=parabola.GetDPz()*parabola.GetDPz();
//       }

      TLorentzVector tlvd(tlvc-tlvm);
      TVector3 vtxd(vtxc-vtxm);
      h_dpx->Fill(tlvd.Px());
      h_dpy->Fill(tlvd.Py());
      h_dpz->Fill(tlvd.Pz());
      h_dE->Fill(tlvd.E());

      h_ppx->Fill(tlvd.Px()/sqrt(covp4[0][0]));
      h_ppy->Fill(tlvd.Py()/sqrt(covp4[1][1]));
      h_ppz->Fill(tlvd.Pz()/sqrt(covp4[2][2]));
      h_pE->Fill(tlvd.E()/sqrt(covp4[3][3]));
      
      h_epx->Fill(sqrt(covp4[0][0]));
      h_epy->Fill(sqrt(covp4[1][1]));
      h_epz->Fill(sqrt(covp4[2][2]));
      h_eE->Fill(sqrt(covp4[3][3]));
      
      h_dvx->Fill(vtxd.X());
      h_dvy->Fill(vtxd.Y());
      h_dvz->Fill(vtxd.Z());
      
      h_pvx->Fill(vtxd.X()/sqrt(covvt[0][0]));
      h_pvy->Fill(vtxd.Y()/sqrt(covvt[1][1]));
      h_pvz->Fill(vtxd.Z()/sqrt(covvt[2][2]));
      
      h_evx->Fill(sqrt(covvt[0][0]));
      h_evy->Fill(sqrt(covvt[1][1]));
      h_evz->Fill(sqrt(covvt[2][2]));

      // set up mc parabola parameters
      TVector3 og(0,0,0);
      TVector3 dj(1,0,0);
      TVector3 dk(0,1,0);
      FairTrackParP parabolamc(vtxm,momm,dummy,dummy,mccd->GetCharge(),og,dj,dk);

      // get destination detector plane
      TVector3 ogc=parabola.GetOrigin();
      TVector3 djc=parabola.GetJVer();
      TVector3 dkc=parabola.GetKVer();
      
      //propagate MC truth to that plane, skipping cov matices and overwriting the parabolamc 
      theAnalysis->Propagator(3,parabolamc,mccd,ogc,true,true,djc,dkc);
      h_dQp->Fill(parabola.GetQp()-parabolamc.GetQp());
      h_dV->Fill(parabola.GetV()-parabolamc.GetV());
      h_dW->Fill(parabola.GetW()-parabolamc.GetW());
      h_dTV->Fill(parabola.GetTV()-parabolamc.GetTV());
      h_dTW->Fill(parabola.GetTW()-parabolamc.GetTW());

      h_pQp->Fill((parabola.GetQp()-parabolamc.GetQp()) / parabola.GetDQp());
      h_pV->Fill((parabola.GetV()-parabolamc.GetV()) / parabola.GetDV());
      h_pW->Fill((parabola.GetW()-parabolamc.GetW()) / parabola.GetDW());
      h_pTV->Fill((parabola.GetTV()-parabolamc.GetTV()) / parabola.GetDTV());
      h_pTW->Fill((parabola.GetTW()-parabolamc.GetTW()) / parabola.GetDTW());
      
      h_eQp->Fill(parabola.GetDQp());
      h_eV->Fill(parabola.GetDV());
      h_eW->Fill(parabola.GetDW());
      h_eTV->Fill(parabola.GetDTV());
      h_eTW->Fill(parabola.GetDTW());
      
    }
  }
  //mm.Close();

  // *** write out all the histos
  // *** create an output file for all histograms
  TFile *out = TFile::Open ( plotfile.Data(),"RECREATE" );
  out->cd();

  h_ppx->Fit("gaus","R","",-2.5,2.5);//h_ppx->GetFit()->SetLineColor(h_ppx->GetLineColor());
  h_ppy->Fit("gaus","R","",-2.5,2.5);//h_ppy->GetFit()->SetLineColor(h_ppy->GetLineColor());
  h_ppz->Fit("gaus","R","",-2.5,2.5);//h_ppz->GetFit()->SetLineColor(h_ppz->GetLineColor());
  h_pE->Fit("gaus","R","",-2.5,2.5);//h_ppz->GetFit()->SetLineColor(h_ppz->GetLineColor());
  h_pvx->Fit("gaus","R","",-2.5,2.5);
  h_pvy->Fit("gaus","R","",-2.5,2.5);
  h_pvz->Fit("gaus","R","",-2.5,2.5);
  h_pQp->Fit("gaus","R","",-2.5,2.5);
  h_pV->Fit("gaus","R","",-2.5,2.5);
  h_pW->Fit("gaus","R","",-2.5,2.5);
  h_pTV->Fit("gaus","R","",-2.5,2.5);
  h_pTW->Fit("gaus","R","",-2.5,2.5);
  
  TListIter iter ( histolist );

  TH1* tmph=NULL;

  while ( tmph= ( TH1* ) iter() ) tmph->Write();

  out->Close();

  //call a macro from Tools.C
  plothistosfromfile ( plotfile.Data(),".pdf",3,4,900 );
  //mm.Show();
}


// FairTrackParP* GetMcPar(RhoCandidate* mctruth)
// {
//   TVector3 pos = mctruth->GetPosition();
//   TVector3 mom = mctruth->GetMomentum();
//   Int_t q=mctruth->GetCharge();
//   TVector3 perr(0.,0.,0.);
//   TVector3 merr(0.,0.,0.);
//   TVector3 o(0.,0.,0.);
//   TVector3 dj(1.,0.,0.);
//   TVector3 dk(0.,1.,0.);
//   FairTrackParP* mcparabola = new FairTrackParP(pos,mom,perr,merr,q,o,dj,dk);
//   return mcparabola;
// }