#include <iostream>
#include <iomanip>
#include <algorithm>

#include <TVector3.h>
#include <TFile.h>

#include <GFRaveVertex.h>
#include <GFRaveTrackParameters.h>

#include "TpcVertex.h"
#include "LambdaChain.h"

#include "FopiLambdaCand.h"

LambdaChain::LambdaChain(bool MC)
  :tracksBranchName("Track"),
  simpleVertexBranchName("TpcSimpleVertex"),
  vertexBranchName("TpcVertex"),
  vertexRaveBranchName("GFVertex"),
  avrVertexBranchName("avr_tpc_vertex"),
  avrVertexRaveBranchName("avr_rave_vertex"),
  lambdaVertexBranchName("LambdaVertex"),
  lambdaVertexRaveBranchName("LambdaRaveVertex"),
  lambdaCandBranchName("LambdaCandRave"),
  mc(MC),xVtxCut(3),yVtxCut(3),zVtxCut(5)
{
  chain=new TChain("cbmsim");

  tracksArray           = new TClonesArray("GFTrack");
  arrays.push_back(tracksArray);
  simpleVertexArray     = new TClonesArray("TpcVertex");
  arrays.push_back(simpleVertexArray);

  vertexArray           = new TClonesArray("TpcVertex");
  arrays.push_back(vertexArray);
  avrVertexArray        = new TClonesArray("TpcVertex");
  arrays.push_back(avrVertexArray);

  vertexRaveArray       = new TClonesArray("GFRaveVertex");
  arrays.push_back(vertexRaveArray);
  avrVertexRaveArray    = new TClonesArray("GFRaveVertex");
  arrays.push_back(avrVertexRaveArray);

  lambdaVertexArray     = new TClonesArray("TpcVertex");
  arrays.push_back(lambdaVertexArray);
  lambdaVertexRaveArray = new TClonesArray("GFRaveVertex");
  arrays.push_back(lambdaVertexRaveArray);

  lambdaCandArray       = new TClonesArray("FopiLambdaCand");
  arrays.push_back(lambdaCandArray);


  simpleCount=0;
}
LambdaChain::LambdaChain(TChain *chain_, bool MC)
  :tracksBranchName("Track"),
  simpleVertexBranchName("TpcSimpleVertex"),
  vertexBranchName("TpcVertex"),
  vertexRaveBranchName("GFVertex"),
  avrVertexBranchName("avr_tpc_vertex"),
  avrVertexRaveBranchName("avr_rave_vertex"),
  lambdaVertexBranchName("LambdaVertexNoSmoothing"),
  lambdaVertexRaveBranchName("LambdaRaveVertexNoSmoothing"),
  lambdaCandBranchName("LambdaCandRave"),
  mc(MC),xVtxCut(3),yVtxCut(3),zVtxCut(5){
  chain=chain_;

  tracksArray           = new TClonesArray("GFTrack");
  arrays.push_back(tracksArray);
  simpleVertexArray     = new TClonesArray("TpcVertex");
  arrays.push_back(simpleVertexArray);

  vertexArray           = new TClonesArray("TpcVertex");
  arrays.push_back(vertexArray);
  avrVertexArray        = new TClonesArray("TpcVertex");
  arrays.push_back(avrVertexArray);

  vertexRaveArray       = new TClonesArray("GFRaveVertex");
  arrays.push_back(vertexRaveArray);
  avrVertexRaveArray    = new TClonesArray("GFRaveVertex");
  arrays.push_back(avrVertexRaveArray);

  lambdaVertexArray     = new TClonesArray("TpcVertex");
  arrays.push_back(lambdaVertexArray);
  lambdaVertexRaveArray = new TClonesArray("GFRaveVertex");
  arrays.push_back(lambdaVertexRaveArray);

  lambdaCandArray       = new TClonesArray("FopiLambdaCand");
  arrays.push_back(lambdaCandArray);

}

void LambdaChain::init(){
  std::cout<<"LambdaChain::init()"<<std::endl;
  chain->SetBranchAddress(tracksBranchName,&tracksArray);
  chain->SetBranchAddress(simpleVertexBranchName,&simpleVertexArray);
  chain->SetBranchAddress(vertexBranchName,&vertexArray);
  chain->SetBranchAddress(vertexRaveBranchName,&vertexRaveArray);

  chain->SetBranchAddress(avrVertexBranchName,&avrVertexArray);
  chain->SetBranchAddress(avrVertexRaveBranchName,&avrVertexRaveArray);

  chain->SetBranchAddress(lambdaVertexBranchName,&lambdaVertexArray);
  chain->SetBranchAddress(lambdaVertexRaveBranchName,&lambdaVertexRaveArray);
  chain->SetBranchAddress(lambdaCandBranchName,&lambdaCandArray);
  //vertexing hists



  //simple
  histograms1D["vtx_simple_z"]=new TH1D("vtx_simple_z","vtx_simple_z;z (cm)",1000,-100,100);
  histograms1D["vtx_simple_x"]=new TH1D("vtx_simple_x","vtx_simple_x;x (cm)",1000,-10,10);
  histograms1D["vtx_simple_y"]=new TH1D("vtx_simple_y","vtx_simple_y;y (cm)",1000,-10,10);
  histograms2D["vtx_simple_xy"]=new TH2D("vtx_simple_xy","vtx_simple_xy;x (cm);y (cm)",500,-10,10,500,-10,10);
  histograms2D["vtx_simple_zx"]=new TH2D("vtx_simple_zx","vtx_simple_zx;z (cm);x (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_simple_zy"]=new TH2D("vtx_simple_zy","vtx_simple_zy;x (cm);y (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_simple_xy_z_m10"]=new TH2D("vtx_simple_xy_z_m10","vtx_simple_xy z<-10;x (cm);y (cm)",500,-10,10,500,-10,10);
//RAVE
  histograms1D["vtx_rave_4tr_z"]=new TH1D("vtx_rave_4tr_z","vtx_rave_4tr_z;z (cm)",1000,-100,100);
  histograms1D["vtx_rave_4tr_x"]=new TH1D("vtx_rave_4tr_x","vtx_rave_4tr_x;x (cm)",1000,-10,10);
  histograms1D["vtx_rave_4tr_y"]=new TH1D("vtx_rave_4tr_y","vtx_rave_4tr_y;y (cm)",1000,-10,10);
  histograms2D["vtx_rave_4tr_xy"]=new TH2D("vtx_rave_4tr_xy","vtx_rave_4tr_xy;x (cm);y (cm)",500,-10,10,500,-10,10);
  histograms2D["vtx_rave_4tr_zx"]=new TH2D("vtx_rave_4tr_zx","vtx_rave_4tr_zx;z (cm);x (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_rave_4tr_zy"]=new TH2D("vtx_rave_4tr_zy","vtx_rave_4tr_zy;x (cm);y (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_rave_4tr_xy_z_m10"]=new TH2D("vtx_rave4tr_xy_z_m10","vtx_rave_4tr_xy z<-10;x (cm);y (cm)",500,-10,10,500,-10,10);

  histograms1D["vtx_rave_z"]=new TH1D("vtx_rave_z","vtx_rave_z;z (cm)",1000,-100,100);
  histograms1D["vtx_rave_x"]=new TH1D("vtx_rave_x","vtx_rave_x;x (cm)",1000,-10,10);
  histograms1D["vtx_rave_y"]=new TH1D("vtx_rave_y","vtx_rave_y;y (cm)",1000,-10,10);
  histograms2D["vtx_rave_xy"]=new TH2D("vtx_rave_xy","vtx_rave_xy;x (cm);y (cm)",500,-10,10,500,-10,10);
  histograms2D["vtx_rave_zx"]=new TH2D("vtx_rave_zx","vtx_rave_zx;z (cm);x (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_rave_zy"]=new TH2D("vtx_rave_zy","vtx_rave_zy;x (cm);y (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_rave_xy_z_m10"]=new TH2D("vtx_rave_xy_z_m10","vtx_rave_xy z<-10;x (cm);y (cm)",500,-10,10,500,-10,10);

  histograms1D["chi_ndf_gfvtx"]=new TH1D("RaveChi2_ndf","Rave chi2/ndf;chi2/ndf",1000,0,100);
  histograms1D["ndf_gfvtx"]=new TH1D("Rave ndf","Rave ndf;ndf",100,0,100);;
  histograms1D["nTracks_rave"]=new TH1D("nTracks_rave","# Tracks pr vertex; # Tracks",35,0,35);
//RAVE avr
  histograms1D["vtx_rave_avr_4tr_z"]=new TH1D("vtx_rave_avr_4tr_z","vtx_rave_avr_4tr_z;z (cm)",1000,-100,100);
  histograms1D["vtx_rave_avr_4tr_x"]=new TH1D("vtx_rave_avr_4tr_x","vtx_rave_avr_4tr_x;x (cm)",1000,-10,10);
  histograms1D["vtx_rave_avr_4tr_y"]=new TH1D("vtx_rave_avr_4tr_y","vtx_rave_avr_4tr_y;y (cm)",1000,-10,10);
  histograms2D["vtx_rave_avr_4tr_xy"]=new TH2D("vtx_rave_avr_4tr_xy","vtx_rave_avr_4tr_xy;x (cm);y (cm)",500,-10,10,500,-10,10);
  histograms2D["vtx_rave_avr_4tr_zx"]=new TH2D("vtx_rave_avr_4tr_zx","vtx_rave_avr_4tr_zx;z (cm);x (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_rave_avr_4tr_zy"]=new TH2D("vtx_rave_avr_4tr_zy","vtx_rave_avr_4tr_zy;x (cm);y (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_rave_avr_4tr_xy_z_m10"]=new TH2D("vtx_rave_avr4tr_xy_z_m10","vtx_rave_avr_4tr_xy z<-10;x (cm);y (cm)",500,-10,10,500,-10,10);

  histograms1D["vtx_rave_avr_z"]=new TH1D("vtx_rave_avr_z","vtx_rave_avr_z;z (cm)",1000,-100,100);
  histograms1D["vtx_rave_avr_x"]=new TH1D("vtx_rave_avr_x","vtx_rave_avr_x;x (cm)",1000,-10,10);
  histograms1D["vtx_rave_avr_y"]=new TH1D("vtx_rave_avr_y","vtx_rave_avr_y;y (cm)",1000,-10,10);
  histograms2D["vtx_rave_avr_xy"]=new TH2D("vtx_rave_avr_xy","vtx_rave_avr_xy;x (cm);y (cm)",500,-10,10,500,-10,10);
  histograms2D["vtx_rave_avr_zx"]=new TH2D("vtx_rave_avr_zx","vtx_rave_avr_zx;z (cm);x (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_rave_avr_zy"]=new TH2D("vtx_rave_avr_zy","vtx_rave_avr_zy;x (cm);y (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_rave_avr_xy_z_m10"]=new TH2D("vtx_rave_avr_xy_z_m10","vtx_rave_avr_xy z<-10;x (cm);y (cm)",500,-10,10,500,-10,10);

  histograms1D["chi_ndf_gfvtx_avr"]=new TH1D("RaveChi2_ndf_avr","Rave avr chi2/ndf;chi2/ndf",1000,0,100);
  histograms1D["ndf_gfvtx_avr"]=new TH1D("Rave_ndf_avr","Rave avr ndf;ndf",100,0,100);;
  histograms1D["nTracks_rave_avr"]=new TH1D("nTracks_rave_avr","# Tracks pr vertex avr; # Tracks",35,0,35);
  histograms1D["nvtx_rave_avr"]=new TH1D("nvtx_rave_avr","# Vtx/Event avr; # vtx",35,0,35);
  histograms2D["nVtxAvrVsnLambdaCand"]=new TH2D("nVtxAvrVsnLambdaCand","# VtxAvr vs # LambdaCands;# vtx;# #Lambda-cands",20,0,20,20,0,20);
  histograms2D["nVtxAvrVsnLambdaCandGood"]=new TH2D("nVtxAvrVsnLambdaCandGood","# VtxAvr vs # LambdaCands;# vtx;# #Lambda-cands(good)",20,0,20,20,0,20);

  //Lambda
  histograms1D["vtx_lambda_z"]=new TH1D("vtx_lambda_z","vtx_lambda_z;z (cm)",1000,-100,100);
  histograms1D["vtx_lambda_x"]=new TH1D("vtx_lambda_x","vtx_lambda_x;x (cm)",1000,-10,10);
  histograms1D["vtx_lambda_y"]=new TH1D("vtx_lambda_y","vtx_lambda_y;y (cm)",1000,-10,10);
  histograms2D["vtx_lambda_xy"]=new TH2D("vtx_lambda_xy","vtx_lambda_xy;x (cm);y (cm)",500,-10,10,500,-10,10);
  histograms2D["vtx_lambda_zx"]=new TH2D("vtx_lambda_zx","vtx_lambda_zx;z (cm);x (cm)",500,-100,100,500,-10,10);
  histograms2D["vtx_lambda_zy"]=new TH2D("vtx_lambda_zy","vtx_lambda_zy;x (cm);y (cm)",500,-100,100,500,-10,10);

  histograms1D["vtx_lambda_dPIP"]=new TH1D("vtx_lambda_dPIP","vtx_lambda_dPIP;cm",1000,0,100);
  histograms1D["chi_ndf_lambdaRave_vtx"]=new TH1D("LambdaRaveChi2_ndf","Lambda Rave chi2/ndf;chi2/ndf",1000,0,100);;
  histograms1D["ndf_lambdaRave_vtx"]=new TH1D("Lambda_Rave_ndf","Lambda Rave ndf;ndf",100,0,100);;

  //mass hists
  histograms1D["L_mass"]=new TH1D("L_mass","#Lambda mass;GeV/c^{2}",1000,1,2);
  histograms2D["L_mass_v_dist"]=new TH2D("L_mass_v_dist","#Lambda mass vs dist;GeV/c^{2}",1000,1,2,1000,0,10);

  histograms2D[TString::Format("L_mass_v_dist_z_%f",zVtxCut).Data()]=new TH2D(TString::Format("L_mass_v_dist_z%f",zVtxCut),TString::Format("#Lambda mass vs dist, z_{vtx}>%f;GeV/c^{2}",zVtxCut),1000,1,2,1000,0,10);
  histograms2D[TString::Format("L_mass_v_dist_z_X_or_y_Xor_x_X").Data()]=new TH2D(TString::Format("L_mass_v_dist_z%f_x%f_y%f3",zVtxCut,yVtxCut,xVtxCut),
                                                                                                       TString::Format("#Lambda mass vs dist,z_{vtx}>%f,x>%f,y>%f;GeV/c^{2}",zVtxCut,yVtxCut,xVtxCut),1000,1,2,1000,0,10);

  histograms1D[TString::Format("L_mass_z_%f",zVtxCut).Data()]=new TH1D(TString::Format("L_mass_z_%f",zVtxCut),TString::Format("#Lambda mass, z_{vtx} > %f, dist #pi-p<5;GeV/c^{2}",zVtxCut),1000,1,2);
  histograms1D[TString::Format("L_mass_z_X_or_y_X_or_x_X")]=new TH1D(TString::Format("L_mass_z_%f_or_y_%for_x_%f",zVtxCut,yVtxCut,xVtxCut),
                                                                                               TString::Format("#Lambda mass, z_{vtx} > %f | |y_{vtx}| > %f | |x_{vtx}| > %f, dist #pi-p<5;GeV/c^{2}",zVtxCut,yVtxCut,xVtxCut),1000,1,2);
  histograms1D[TString::Format("L_mass_d2_z_%f",zVtxCut).Data()]=new TH1D(TString::Format("L_mass_d2_z_%f",zVtxCut),
                                                                   TString::Format("#Lambda mass, z_{vtx} > %f, dist #pi-p<2;GeV/c^{2}",zVtxCut),1000,1,2);
  histograms1D[TString::Format("L_mass_d2_z_X_or_y_Xor_x_X").Data()]=new TH1D(TString::Format("L_mass_d2_z_%f_or_y_%for_x_%f",zVtxCut,yVtxCut,xVtxCut),
                                                                                                  TString::Format("#Lambda mass, z_{vtx} > %f | |y_{vtx}| > %f | |x_{vtx}| > %f, dist #pi-p<2;GeV/c^{2}",zVtxCut,yVtxCut,xVtxCut),1000,1,2);

  histograms1D[TString::Format("L_mass_d3_z_%f",zVtxCut).Data()]=new TH1D(TString::Format("L_mass_d3_z_%f",zVtxCut),TString::Format("#Lambda mass, z_{vtx} > %f, dist #pi-p<3;GeV/c^{2}",zVtxCut),1000,1,2);
  histograms1D[TString::Format("L_mass_d3_z_X_or_y_Xor_x_X").Data()]=new TH1D(TString::Format("L_mass_d3_z_%f_or_y_%for_x_%f",zVtxCut,yVtxCut,xVtxCut),
                                                                                                  TString::Format("#Lambda mass, z_{vtx} > 5 | |y_{vtx}| > 3 | |x_{vtx}| > 3, dist #pi-p<3;GeV/c^{2}",zVtxCut,yVtxCut,xVtxCut),1000,1,2);

  //track hists
  histograms1D["nTracks"]=new TH1D("nTracks","# Tracks;# Tracks",35,0,35);
  std::cout<<"LambdaChain::init() end"<<std::endl;

}

void LambdaChain::loop(unsigned int nEv){
  std::cout<<"LambdaChain::loop()"<<std::endl;
  unsigned int nEntries = chain->GetEntries();
  std::cout<<"chain has : "<<nEntries<<" entries, looping over ";

  if(nEv!=0)
    nEntries=nEv;
  std::cout<<nEntries<<std::endl;
  int timeStart=clock();
  for(unsigned int iEv=0;iEv<nEntries;++iEv){
    for(unsigned int i=0;i<arrays.size();++i){
      arrays[i]->Delete();
    }
    chain->GetEvent(iEv);
    if((iEv%250)==0){
      printEv(iEv,nEntries,timeStart);
    }
    unsigned int nLambda=lambdaCandArray->GetEntriesFast();
    unsigned int nSimpleVtx=simpleVertexArray->GetEntriesFast();
    unsigned int nVtx=vertexArray->GetEntriesFast();
    unsigned int nAvrVtx=avrVertexArray->GetEntriesFast();
    histograms1D["nvtx_rave_avr"]->Fill(nAvrVtx);
    histograms2D["nVtxAvrVsnLambdaCand"]->Fill(nAvrVtx,nLambda);

    unsigned int nTracks=tracksArray->GetEntriesFast();
    histograms1D["nTracks"]->Fill(nTracks);
    unsigned int maxLoop=std::max(nLambda,std::max(nSimpleVtx,std::max(nVtx,nAvrVtx)));

    //    std::cout<<"nLambda "<<nLambda<<std::endl;
    //    std::cout<<"nSimpleVtx "<<nSimpleVtx<<std::endl;
    //    std::cout<<"nVtx "<<nVtx<<std::endl;
    //    std::cout<<"maxLoop "<<maxLoop<<std::endl;
    int goodLambda=0;
    for(unsigned int i=0;i<maxLoop;++i){
      if(i<nLambda){
        FopiLambdaCand* cand= (FopiLambdaCand*)lambdaCandArray->At(i);
        TpcVertex* lVtx=(TpcVertex*)lambdaVertexArray->At(i);
        GFRaveVertex* lRaveVtx=(GFRaveVertex*)lambdaVertexRaveArray->At(i);
        goodLambda+=fillLambdaHists(cand,lVtx,lRaveVtx);
      }
      if(i<nSimpleVtx){
        TpcVertex* simpleVtx=(TpcVertex*)simpleVertexArray->At(i);
        fillSimpleVtxHists(simpleVtx);
      }
      if(i<nVtx){
        TpcVertex* theVtx=(TpcVertex*)vertexArray->At(i);
        GFRaveVertex* theRaveVtx=(GFRaveVertex*)vertexRaveArray->At(i);
        fillRaveVtxHists(theVtx,theRaveVtx);
      }
      if(i<nAvrVtx){
        TpcVertex* theVtx=(TpcVertex*)avrVertexArray->At(i);
        GFRaveVertex* theRaveVtx=(GFRaveVertex*)avrVertexRaveArray->At(i);
        fillRaveAvrVtxHists(theVtx,theRaveVtx);
      }
    }
    histograms2D["nVtxAvrVsnLambdaCandGood"]->Fill(nAvrVtx,goodLambda);

  }
  double diffms=((double)(clock()-timeStart))/(CLOCKS_PER_SEC);
  std::cout << std::endl<< "time from start to end " << diffms << std::endl;
  std::cout<<"LambdaChain::loop() end"<<std::endl;
}

void LambdaChain::printEv(unsigned int iEv, unsigned int nEntries, int timeStart){
  double percent = iEv*100./nEntries;
  std::string bla="";
  switch (simpleCount){
  case 0:
    bla="-";
    simpleCount++;
    break;
  case 1:
    bla="\\";
    simpleCount++;
    break;
  case 2:
    bla="|";
    simpleCount++;
    break;
  case 3:
    bla="/";
    simpleCount=0;
    break;
  }
  std::cout<<"\r"<<iEv<<": "<<std::setiosflags(std::ios::fixed)<<std::setprecision(3)<<percent<<"\% "<<bla<<"Time Elapsed: "<<((double)(clock()-timeStart))/(CLOCKS_PER_SEC)<<std::flush;


}
int LambdaChain::fillLambdaHists(FopiLambdaCand *cand, TpcVertex *lVtx, GFRaveVertex *lRaveVtx){
  double x=lVtx->getPos().X();
  double y=lVtx->getPos().Y();
  double z=lVtx->getPos().Z();
  histograms1D["vtx_lambda_z"]->Fill(z);
  histograms1D["vtx_lambda_x"]->Fill(x);
  histograms1D["vtx_lambda_y"]->Fill(y);
  histograms2D["vtx_lambda_xy"]->Fill(x,y);
  histograms2D["vtx_lambda_zx"]->Fill(z,x);
  histograms2D["vtx_lambda_zy"]->Fill(z,y);
  TVector3 piPos=lRaveVtx->getParameters(0)->getPos();
  TVector3 pPos=lRaveVtx->getParameters(1)->getPos();
  TVector3 dist=piPos-pPos;
  histograms1D["vtx_lambda_dPIP"]->Fill(dist.Mag());

  histograms1D["chi_ndf_lambdaRave_vtx"]->Fill(lRaveVtx->getChi2()/lRaveVtx->getNdf());
  histograms1D["ndf_lambdaRave_vtx"]->Fill(lRaveVtx->getNdf());
  histograms1D["L_mass"]->Fill(cand->getState().Mag());
  histograms2D["L_mass_v_dist"]->Fill(cand->getState().Mag(),dist.Mag());

  if((z>zVtxCut)&&dist.Mag()<5){
    histograms1D[TString::Format("L_mass_z_%f",zVtxCut).Data()]->Fill(cand->getState().Mag());
    histograms2D[TString::Format("L_mass_v_dist_z_%f",zVtxCut).Data()]->Fill(cand->getState().Mag(),dist.Mag());
    if(dist.Mag()<2){
      histograms1D[TString::Format("L_mass_d2_z_%f",zVtxCut).Data()]->Fill(cand->getState().Mag());
    }
    if(dist.Mag()<3){
      histograms1D[TString::Format("L_mass_d3_z_%f",zVtxCut).Data()]->Fill(cand->getState().Mag());
    }
  }
  if((z>zVtxCut||fabs(x)>xVtxCut||fabs(y)>yVtxCut)&&dist.Mag()<5){
    histograms1D[TString::Format("L_mass_z_X_or_y_X_or_x_X")]->Fill(cand->getState().Mag());
    histograms2D[TString::Format("L_mass_v_dist_z_X_or_y_Xor_x_X")]->Fill(cand->getState().Mag(),dist.Mag());
    if(dist.Mag()<2){
      histograms1D[TString::Format("L_mass_d2_z_X_or_y_Xor_x_X",zVtxCut,yVtxCut,xVtxCut)]->Fill(cand->getState().Mag());
    }
    if(dist.Mag()<3){
      histograms1D[TString::Format("L_mass_d3_z_X_or_y_Xor_x_X",zVtxCut,yVtxCut,xVtxCut)]->Fill(cand->getState().Mag());
    }
  }
  if(dist.Mag()<5){
    return 1;
  }else {
    return 0;
  }
}
void LambdaChain::fillRaveVtxHists(TpcVertex *theVtx, GFRaveVertex *theRaveVtx){
  double x=theVtx->getPos().X();
  double y=theVtx->getPos().Y();
  double z=theVtx->getPos().Z();
  unsigned int nTracks=theVtx->getNTracks();
  histograms1D["nTracks_rave"]->Fill(nTracks);
  histograms1D["vtx_rave_z"]->Fill(z);
  histograms1D["vtx_rave_x"]->Fill(x);
  histograms1D["vtx_rave_y"]->Fill(y);
  histograms2D["vtx_rave_xy"]->Fill(x,y);
  histograms2D["vtx_rave_zx"]->Fill(z,x);
  histograms2D["vtx_rave_zy"]->Fill(z,y);
  if(z<-10){
    histograms2D["vtx_rave_xy_z_m10"]->Fill(x,y);
  }
  if(theVtx->getNTracks()>3){
    histograms1D["vtx_rave_4tr_z"]->Fill(z);
    histograms1D["vtx_rave_4tr_x"]->Fill(x);
    histograms1D["vtx_rave_4tr_y"]->Fill(y);
    histograms2D["vtx_rave_4tr_xy"]->Fill(x,y);
    histograms2D["vtx_rave_4tr_zx"]->Fill(z,x);
    histograms2D["vtx_rave_4tr_zy"]->Fill(z,y);
    if(z<-10){
      histograms2D["vtx_rave_4tr_xy_z_m10"]->Fill(x,y);
    }
  }
  histograms1D["chi_ndf_gfvtx"]->Fill(theRaveVtx->getChi2()/theRaveVtx->getNdf());
  histograms1D["ndf_gfvtx"]->Fill(theRaveVtx->getNdf());
}
void LambdaChain::fillRaveAvrVtxHists(TpcVertex *theVtx, GFRaveVertex *theRaveVtx){
  double x=theVtx->getPos().X();
  double y=theVtx->getPos().Y();
  double z=theVtx->getPos().Z();
  unsigned int nTracks=theVtx->getNTracks();
  histograms1D["nTracks_rave_avr"]->Fill(nTracks);
  histograms1D["vtx_rave_avr_z"]->Fill(z);
  histograms1D["vtx_rave_avr_x"]->Fill(x);
  histograms1D["vtx_rave_avr_y"]->Fill(y);
  histograms2D["vtx_rave_avr_xy"]->Fill(x,y);
  histograms2D["vtx_rave_avr_zx"]->Fill(z,x);
  histograms2D["vtx_rave_avr_zy"]->Fill(z,y);
  if(z<-10){
    histograms2D["vtx_rave_avr_xy_z_m10"]->Fill(x,y);
  }
  if(theVtx->getNTracks()>3){
    histograms1D["vtx_rave_avr_4tr_z"]->Fill(z);
    histograms1D["vtx_rave_avr_4tr_x"]->Fill(x);
    histograms1D["vtx_rave_avr_4tr_y"]->Fill(y);
    histograms2D["vtx_rave_avr_4tr_xy"]->Fill(x,y);
    histograms2D["vtx_rave_avr_4tr_zx"]->Fill(z,x);
    histograms2D["vtx_rave_avr_4tr_zy"]->Fill(z,y);
    if(z<-10){
      histograms2D["vtx_rave_avr_4tr_xy_z_m10"]->Fill(x,y);
    }
  }
  histograms1D["chi_ndf_gfvtx_avr"]->Fill(theRaveVtx->getChi2()/theRaveVtx->getNdf());
  histograms1D["ndf_gfvtx_avr"]->Fill(theRaveVtx->getNdf());
}
void LambdaChain::fillSimpleVtxHists(TpcVertex *theVtx){
  double x=theVtx->getPos().X();
  double y=theVtx->getPos().Y();
  double z=theVtx->getPos().Z();
  histograms1D["vtx_simple_z"]->Fill(z);
  histograms1D["vtx_simple_x"]->Fill(x);
  histograms1D["vtx_simple_y"]->Fill(y);
  histograms2D["vtx_simple_xy"]->Fill(x,y);
  histograms2D["vtx_simple_zx"]->Fill(z,x);
  histograms2D["vtx_simple_zy"]->Fill(z,y);
  if(z<-10){
    histograms2D["vtx_simple_xy_z_m10"]->Fill(x,y);
  }
}
void LambdaChain::saveHistograms(TString file){
  std::cout<<"LambdaChain::saveHistograms("<<file<<")"<<std::endl;
  TFile* outFile= new TFile(file,"RECREATE");
  std::map<TString,TH1D*>::iterator it;
  for(it=histograms1D.begin();it!=histograms1D.end();++it){
    (*it).second->Write();
  }

  std::map<TString,TH2D*>::iterator it2;
  for(it2=histograms2D.begin();it2!=histograms2D.end();++it2){
    (*it2).second->Write();
  }
  std::cout<<"LambdaChain::saveHistograms("<<file<<") end"<<std::endl;
}