#define V0Analysis_cxx
#include "V0Analysis.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "iostream.h"
#include "TVector3.h"
#include "TPolyLine.h"
#include "TEllipse.h"
#include "TLorentzVector.h"


bool _2D_ = true; 


// Interaction Point definition --------
float IP_x = -0.024;
float IP_y = -0.0489;
//-------------------------------------

// Average mass
float K0MASS = 0.497672;
float LAMBDAMASS = 1.115683;
float JPSIMASS = 1.115683;
//-------------------------------------

//Tracking layers
float SMT_1 = 2.75;
float SMT_2 = 10.43;
float CFT_1 = 20.08;
float CFT_2 = 51.5;
//-------------------------------------

// Selection cuts
float R    = SMT_1; // maximum decay length
float CA   = 0.9999;
float CHI2 = 15.0;

bool vetoLambda = true; // for K0 search
bool vetoK0 = true; // for Lambda search
//-------------------------------------


void V0Analysis::readJets(TObjArray& jets)
{
  jets.Clear();
  for(int i=0; i<JCCB_JCCBnjet; i++) {
    bool hot = false;
    float deta = JCCB_JCCBdEta[i]; float dphi = JCCB_JCCBdPhi[i];
    if(dphi>=30.5 && dphi<=32 && deta>=-10 && deta<=-8) hot = true;
    if(dphi>=31.2 && dphi<=31.7 && deta>=15.8 && deta<=17.8) hot = true;
    if(dphi>=38.4 && dphi<=39.2 && deta>=-8.4 && deta<=-6) hot = true;
    if(dphi>=38 && dphi<=44 && deta>=6 && deta<=10) hot = true;
    if(dphi>=49.4 && dphi<=50.4 && deta>=-6 && deta<=-5.5) hot = true;
    if(dphi>=60.8 && dphi<=61.2 && deta>=-7 && deta<=-6) hot = true;

    if(!hot) {
      bool cEMF = JCCB_JCCBEMF[i]>=0.05 && JCCB_JCCBEMF[i]<=0.9;
      bool cCHF = JCCB_JCCBCHF[i] <= 0.4;
      bool cN90 = JCCB_JCCBn90[i] >= 2;
      bool cHF = JCCB_JCCBHotF[i] <= 10;
      if(cEMF && cCHF && cN90 && cHF & JCCB_JCCBPT[i]>=10.0) {
	TLorentzVector* jet = new TLorentzVector(JCCB_JCCBPx[i],JCCB_JCCBPy[i],
						 JCCB_JCCBPz[i],JCCB_JCCBE[i]);
	
	float phi = jet->Phi(); float eta = jet->Eta();
	if(phi>=0.1 && phi<=0.2 && eta>=-2 && eta<=-1.5) hot = true;
	if(phi>=0.25 && phi<=0.35 && eta>=-1 && eta<=-0.5) hot = true;
	if(phi>=1 && phi<=1.1 && eta>=-3.1 && eta<=-2.8) hot = true;
	if(phi>=1.25 && phi<=1.35 && eta>=-0.8 && eta<=0) hot = true;

	if(!hot) jets.Add(jet);
      }
    }
  }
}




//constructor
V0Analysis::V0Analysis(TTree *tree) : Global(tree) {}

//destructor
V0Analysis::~V0Analysis() {}

// Decay length
float V0Analysis::decayLength(float xPV, float yPV, float zPV, int j)
{
  float lx = V0_V0_x[j]-xPV;
  float ly = V0_V0_y[j]-yPV;
  float lz = V0_V0_z[j]-zPV;

  float px = V0_V0_px[j];
  float py = V0_V0_py[j];
  float pz = V0_V0_pz[j];

  if(_2D_) {
    lz = 0;
    pz = 0;
  }

  float sign = lx*px+ly*py+lz*pz;
  float l = TMath::Sqrt(lx*lx+ly*ly+lz*lz);

  if(sign<0.0) l = l*(-1);
  return l;
}
// Decay length
float V0Analysis::decayLength_sksel(float xPV, float yPV, float zPV, int j)
{
  float lx = SKSEL333_kssx[j]-xPV;
  float ly = SKSEL333_kssy[j]-yPV;
  float lz = SKSEL333_kssz[j]-zPV;

  float px = SKSEL333_ksspx[j];
  float py = SKSEL333_ksspy[j];
  float pz = SKSEL333_ksspz[j];

  if(_2D_) {
    lz = 0;
    pz = 0;
  }

  float sign = lx*px+ly*py+lz*pz;
  float l = TMath::Sqrt(lx*lx+ly*ly+lz*lz);

  if(sign<0.0) l = l*(-1);
  return l;
}

// Lifetime
float V0Analysis::properLifeTime(float xPV, float yPV, float zPV, int j)
{
  float lx = V0_V0_x[j]-xPV;
  float ly = V0_V0_y[j]-yPV;
  float lz = V0_V0_z[j]-zPV;
  
  float px = V0_V0_px[j];
  float py = V0_V0_py[j];
  float pz = V0_V0_pz[j];
  
  if(_2D_) {
    lz = 0;
    pz = 0;
  }
  
  float p2 = px*px+py*py+pz*pz;
  float L = (lx*px+ly*py+lz*pz)*K0MASS/p2;

  return L;
}

// colinearity angle
float V0Analysis::colinearityAngle(float xPV, float yPV, float zPV, int j)
{
  float lx = V0_V0_x[j]-xPV;
  float ly = V0_V0_y[j]-yPV;
  float lz = V0_V0_z[j]-zPV;

  float px = V0_V0_px[j];
  float py = V0_V0_py[j];
  float pz = V0_V0_pz[j];

  if(_2D_) {
    lz = 0;
    pz = 0;
  }

  float l = TMath::Sqrt(lx*lx+ly*ly+lz*lz);
  float p = TMath::Sqrt(px*px+py*py+pz*pz);

  float oa = (lx*px+ly*py+lz*pz)/(p*l);
  return oa;
}


void V0Analysis::signedIP(int track, float jetPhi, float& dca, float& sdca)
{
  dca = -999; sdca = -999;
  float pT = TMath::Abs(1./GTRACK333_qoverpt[track]);
  float phi = GTRACK333_phidir[track]; if(phi<0) phi += 2*TMath::Pi();
  double px = pT*TMath::Cos(phi);
  double py = pT*TMath::Sin(phi);
  double pz = pT*GTRACK333_tanl[track];
  double E = TMath::Sqrt(px*px+py*py+pz*pz+0.13956995*0.13956995);  
  TLorentzVector trk(px,py,pz,E);
  float deltaPhi = trk.Phi() - jetPhi;
  float dca0 = GTRACK333_dca0[track];
  double sign = dca0*TMath::Sin(deltaPhi);
  if(GTRACK333_dcaerr0[track]>0) { 
    dca = TMath::Abs(dca0);
    if(sign<0) dca = -dca;
    sdca = dca/GTRACK333_dcaerr0[track];
  }
}



// Tracks associated to primary vertex
void V0Analysis::getPVtracks(Int_t& n, Int_t* indices)
{
  n = 0;
  float max = 2; int j = -1;
  for(int i=0; i<PVSEL333_psnvtx333; i++) {
    float dx = PVSEL333_psx[i]-IP_x; float dy = PVSEL333_psy[i]-IP_y;
    float D = TMath::Sqrt(dx*dx+dy*dy);
    if(D<0.5 && PVSEL333_psntrk[i]>max) {
      max = PVSEL333_psntrk[i]; j = i;
    }
  }
  if(max>=2) trackIds(j,n,indices);
}
    
// V0 tracks kinematics (original tracks at DCA)
void V0Analysis::get_kinematics(int& v0, float& q1, float& px1, float& py1, 
			   float& pz1, float& p1, float& dedx1, int& nsmt1,
			   int& ncft1, float& dca1, int& id1,
			   float& q2, float& px2, float& py2, float& pz2,  
			   float& p2, float& dedx2, int& nsmt2, int& ncft2,
			   float& dca2, int& id2)
{
  float pt1, pt2;
  for(int i=0; i<GTRACK333_ngtrk333;i++) {
    if(GTRACK333_gtindex[i]==V0_V0_id1[v0] && GTRACK333_qoverpt[i]!=0) {
      q1 = GTRACK333_qoverpt[i]/TMath::Abs(GTRACK333_qoverpt[i]);
      pt1 =TMath::Abs(1./GTRACK333_qoverpt[i]);
      px1 = pt1*TMath::Cos(GTRACK333_phidir[i]);
      py1 = pt1*TMath::Sin(GTRACK333_phidir[i]);
      pz1 = pt1*GTRACK333_tanl[i];
      p1 = TMath::Sqrt(pt1*pt1+pz1*pz1);
      dedx1 = GTRACK333_dedx[i];
      nsmt1 = GTRACK333_nsmt[i];
      ncft1 = GTRACK333_ncft[i];
      dca1 = TMath::Abs(GTRACK333_rsig[i]);
      id1 = i;
    }
    if(GTRACK333_gtindex[i]==V0_V0_id2[v0] && GTRACK333_qoverpt[i]!=0) {
      q2 = GTRACK333_qoverpt[i]/TMath::Abs(GTRACK333_qoverpt[i]);
      pt2 =TMath::Abs(1./GTRACK333_qoverpt[i]);
      px2 = pt2*TMath::Cos(GTRACK333_phidir[i]);
      py2 = pt2*TMath::Sin(GTRACK333_phidir[i]);
      pz2 = pt2*GTRACK333_tanl[i];
      p2 = TMath::Sqrt(pt2*pt2+pz2*pz2);
      dedx2 = GTRACK333_dedx[i];
      nsmt2 = GTRACK333_nsmt[i];
      ncft2 = GTRACK333_ncft[i];
      dca2 = TMath::Abs(GTRACK333_rsig[i]);
      id2 = i;
    }
  }
}

// Checks the number of V0 tracks belonging to the primary vertex.
bool V0Analysis::belongsToPV(int& v0, int n, Int_t* pvIndices)
{
  bool answer = false;
  int id1 = V0_V0_id1[v0]; int id2 = V0_V0_id2[v0];

  for(int g=0; g<n; g++) {
    if(GTRACK333_gtindex[pvIndices[g]]==id1 || GTRACK333_gtindex[pvIndices[g]]==id2)
      answer = true;
  }
  return answer;
}

// Obtain the track indices of a vertex
void V0Analysis::trackIds(int index, Int_t& n, Int_t* indices)
{
  int bits[20] = {0};

  bits[0] = PVSEL333_psattgt0[index];
  bits[1] = PVSEL333_psattgt1[index];
  bits[2] = PVSEL333_psattgt2[index];
  bits[3] = PVSEL333_psattgt3[index];
  bits[4] = PVSEL333_psattgt4[index];
  bits[5] = PVSEL333_psattgt5[index];
  bits[6] = PVSEL333_psattgt6[index];
  bits[7] = PVSEL333_psattgt7[index];
  bits[8] = PVSEL333_psattgt8[index];
  bits[9] = PVSEL333_psattgt9[index];

  n = 0;
  for(int itk=0; itk<GTRACK333_ngtrk333; itk++) {
    int ibit = itk%32;
    int iset = itk/32;
    if( 1 & (bits[iset] >> ibit) ) {
      indices[n] = itk;
      n++;
    }
  }
}

//display tracks
void V0Analysis::PlotTrack(Float_t p1, Float_t p2, Float_t p3, Float_t p4, 
			   Float_t p5, Int_t color, Int_t style, 
			   Float_t dim)
{  
  Float_t CLIGHT = 2.99792458e10;
  Float_t BFAC = 1.0e-13 * CLIGHT;

  Double_t phi_dir = p3;
  Double_t xp =  p1 * sin(phi_dir);
  Double_t yp = -p1 * cos(phi_dir);
  Double_t BMAG = -2.0;            // magnetic field is reversed...
  Double_t CONSB = BFAC * BMAG;    // BFAC = 1.0e-13*CLIGHT (trfutil/TRFMath.h)
  Double_t  r =  1./p5/CONSB;
  Double_t x0 = xp + r*sin(phi_dir);
  Double_t y0 = yp - r*cos(phi_dir);
  Double_t R = -r;
 
  Float_t sign = -1.;

  Float_t range;
  //if(TMath::Abs(1/p5)<4) {
  //range = TMath::Abs(1/p5)/4*dim/TMath::Abs(R);}
  //else { range = dim/TMath::Abs(R);}
  range = dim/TMath::Abs(R);

  Double_t x[100],y[100]; 
  for (Int_t i=0; i<100; i++) {
    x[i] = x0 + R*TMath::Sin((i-50)*sign*range/50.0+phi_dir);
    y[i] = y0 - R*TMath::Cos((i-50)*sign*range/50.0+phi_dir);
  }
  
  TPolyLine* track = new TPolyLine(100, x, y);
  track->SetLineWidth(1);
  track->SetLineColor(color);
  track->SetLineStyle(style);
  track->Draw();
}


///---------------------------------------------
/// Event Display
///---------------------------------------------
void V0Analysis::Display(int eventNumber, float dim)
{
  gStyle->SetOptStat(000000);
  
  if (fChain == 0) return; 
  Int_t nentries = Int_t(fChain->GetEntries());
  Int_t nbytes = 0, nb = 0;
  if (LoadTree(eventNumber) < 0) return;
  nb = fChain->GetEntry(eventNumber);   nbytes += nb;
  cout << endl << "Event #" << eventNumber << endl;
    
  //1- Create Canvas
  TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
  if (c1) {c1->Clear();} 
  else {
    c1 = new TCanvas("c1","c1",700,600);
  }  
  c1->SetGrid();
  c1->SetLeftMargin(0.15);
  c1->SetBottomMargin(0.15);
  c1->SetLogy(0);
  gStyle->SetCanvasColor(10);
  gStyle->SetStatColor(10);
  gStyle->SetTitleColor(10);
  gStyle->SetPadColor(10);
  gStyle->SetPaperSize(20,24);
  gStyle->SetStatFont(62);
  gStyle->SetOptStat(000000);

  //2- Create frame
  TH2F *hr = new TH2F("hr","",2,-dim,dim,2,-dim,dim);
  hr->SetXTitle("x [cm]");
  hr->SetYTitle("y [cm]");
  hr->GetXaxis()->SetTitleSize(0.05);
  hr->GetYaxis()->SetTitleSize(0.05);
  hr->SetNdivisions(6,"X");
  hr->SetNdivisions(6,"Y");
  hr->Draw();
  //3- Draw Silicon layers
  TEllipse* cir1 = new TEllipse(0,0,SMT_1);
  cir1->SetLineWidth(2);
  cir1->SetLineColor(1);
  cir1->Draw();
  TEllipse* cir2 = new TEllipse(0,0,SMT_2);
  cir2->SetLineWidth(2);
  cir2->SetLineColor(1);
  cir2->Draw();
  TEllipse* cir3 = new TEllipse(0,0,CFT_1);
  cir3->SetLineWidth(2);
  cir3->SetLineColor(1);
  cir3->Draw();
  TEllipse* cir4 = new TEllipse(0,0,CFT_2);
  cir4->SetLineWidth(2);
  cir4->SetLineColor(1);
  cir4->Draw();
  //4- Draw Tracks
  cout << endl << GTRACK333_ngtrk333 << " tracks" << endl;
  for(int i=0; i<GTRACK333_ngtrk333; i++) {
    PlotTrack(GTRACK333_rsig[i],GTRACK333_zca[i],GTRACK333_phidir[i],GTRACK333_tanl[i],
	      GTRACK333_qoverpt[i],4,1,dim);
    cout << "gtrack id = " << GTRACK333_gtindex[i] << ", pT = ";
    cout << 1./GTRACK333_qoverpt[i] << endl;
  }
  cout << "-------------------------------" << endl;
  //5- Draw primary vertex
  if(PVSEL333_psnvtx333>0) {
    int n; int indices[100];
    TEllipse* vertex = new TEllipse(PVSEL333_psx[0],PVSEL333_psy[0],dim/35);
    vertex->SetLineColor(1); vertex->SetFillColor(1);
    vertex->Draw();
    cout << "Primary vertex: ";
    cout << "(x,y) = " << PVSEL333_psx[0] << ", " << PVSEL333_psy[0] << endl;
    cout << "chi2 = " << PVSEL333_pschisq[0] << ", ntrk = ";
    cout << PVSEL333_psntrk[0] << endl;
    trackIds(0,n,indices);
    cout <<"indices: ";
    for(int g=0; g<n; g++) {
      cout << GTRACK333_gtindex[indices[g]] << " ";
      PlotTrack(GTRACK333_rsig[indices[g]],GTRACK333_zca[indices[g]],
		GTRACK333_phidir[indices[g]],GTRACK333_tanl[indices[g]],
		GTRACK333_qoverpt[indices[g]],3,1,dim);
    }
    cout << endl;
  }    
  cout << "-------------------------------" << endl;
  //6- Draw V0s
  for(int i=0; i<V0_V0_nv0;i++) {
    if(V0_V0_ch[i]==-1) {
      float dl = decayLength(IP_x, IP_y, 0.0, i);
      float oa = colinearityAngle(IP_x, IP_y, 0.0, i);
      if(dl<10.0 && dl>0.2 && oa>0.9999) {
	cout << "V0 vertex (x,y,decayLength) = ";
	cout << V0_V0_x[i] << ", " << V0_V0_y[i] << ";  " << dl << endl;
	cout << "  colinearity = " << oa << ";  track ids = ";
	cout << V0_V0_id1[i] << ", " << V0_V0_id2[i] << endl;
	TEllipse* vertex = new TEllipse(V0_V0_x[i],V0_V0_y[i],dim/35);
	vertex->SetLineColor(2); vertex->SetFillColor(2);
	vertex->Draw();
      }
    }
  }
  cout << "-------------------------------" << endl;
}



///---------------------------------------------------
/// Tracking studies
///
///---------------------------------------------------
void V0Analysis::AnalyzeTracking(char* fileName)
{
  TH1F* h_ntrk_all = new TH1F("h_ntrk_all","Track multiplicity",31,-0.5,30.5);

  TH1F* h_ntrk_pv = new TH1F("h_ntrk_pv","Primary Vertex  multiplicity",
			     29,1.5,30.5);

  TH1F* h_ntrk_kpv = new TH1F("h_ntrk_kpv","Primary Vertex  multiplicity [Kalman Filter]",29,1.5,30.5);

  TH2F* h_ntrk_cor = new TH2F("h_ntrk_cor","",31,-0.5,30.5,29,1.5,30.5);
  h_ntrk_cor->SetXTitle("GTracks >=11 hits");
  h_ntrk_cor->SetYTitle("Primary Vertex multiplicity");

  TH1F* h_dca = new TH1F("h_dca","Impact Parameter [p_{T} > 1.5 GeV/c]",
			 200,-0.2,0.2);
  h_dca->SetXTitle("cm");

  TH1F* h_dcaSig_p = new TH1F("h_dcaSig_p","Signed Impact Parameter Significance [p_{T} > 1.5 GeV/c]",200,0,30);
  TH1F* h_dcaSig_n = new TH1F("h_dcaSig_n","Signed Impact Parameter Significance [p_{T} > 1.5 GeV/c]",200,0,30);

  TH1F* h_dcaSig = new TH1F("h_dcaSig","Signed Impact Parameter Significance [p_{T} > 1.5 GeV/c]",200,-30,30);


  int aa = 1000; int nn = 0;
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());  
  Int_t nbytes = 0, nb = 0;
  for (Int_t jentry=0; jentry<nentries;jentry++) {
    Int_t ientry = LoadTree(jentry); 
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    if(jentry/aa>nn) cout << "Event #" << jentry << endl;
    nn = jentry/aa;

    //count # tracks with 11 hits
    int ntrk = 0;
    for(int i=0; i<GTRACK333_ngtrk333; i++) 
      if(GTRACK333_nsmt[i]+GTRACK333_ncft[i]>=11) ntrk++;

    h_ntrk_all->Fill(ntrk);

    if(PVSEL333_psntrk>0) {
      h_ntrk_pv->Fill(PVSEL333_psntrk[0]);
      h_ntrk_cor->Fill(ntrk,PVSEL333_psntrk[0]);
    }
    else {
      h_ntrk_pv->Fill(0);
      h_ntrk_cor->Fill(ntrk,0);
    }

    if(PKSEL333_kpsntrk>0)  
      h_ntrk_kpv->Fill(PKSEL333_kpsntrk[0]);
    else
      h_ntrk_kpv->Fill(0);


    if(PVSEL333_psnvtx333>0 && PVSEL333_psntrk[0]>=3) {//3 or more tracks in PV
      //impact parameter (all tracks)
      for(int i=0; i<GTRACK333_ngtrk333; i++) {
	if(TMath::Abs(1./GTRACK333_qoverpt[i])>1.5 && 
	   (GTRACK333_nsmt[i]+GTRACK333_ncft[i]>=10)) {
	  h_dca->Fill(GTRACK333_dca0[i]);
	}


	TObjArray jets; readJets(jets);
	for(int y=0; y<jets.GetEntries(); y++) {
	  TLorentzVector* jet = (TLorentzVector*) jets[y];
 
	  float jETA = jet->Eta();
	  float jPHI = jet->Phi();
	  if(jPHI<-TMath::Pi()) jPHI += 2*TMath::Pi();
	  if(jPHI>TMath::Pi()) jPHI -= 2*TMath::Pi();
	  //match tracks
	  for(int j=0; j<GTRACK333_ngtrk333; j++) {
	    if(TMath::Abs(1./GTRACK333_qoverpt[i])>1.5 && 
	       (GTRACK333_nsmt[j]+GTRACK333_ncft[j]>=10)) {
	      float phi = GTRACK333_phidir[j];
	      if(phi<-TMath::Pi()) phi += 2*TMath::Pi();
	      if(phi>TMath::Pi()) phi -= 2*TMath::Pi();
	      //if(phi<0) phi += 2*TMath::Pi();
	      float theta = TMath::ATan(GTRACK333_tanth[j]); 
	      if(theta<0.0) theta += TMath::Pi();
	      float eta = -log( tan( theta/2.0 ) );
	      double deltaPhi = phi-jPHI;
	      double DP = deltaPhi;
	      if(DP>TMath::Pi()) DP -= 2*TMath::Pi();
	      if(DP<-TMath::Pi()) DP += 2*TMath::Pi();
	      double deltaEta = eta-jETA;
	      double dr = TMath::Sqrt(DP*DP+deltaEta*deltaEta);
	      if(dr<=0.5) {
		float dca = GTRACK333_dca0[j];
		double sign = dca*TMath::Sin(deltaPhi);
		if(GTRACK333_dcaerr0[j]>0 && TMath::Abs(dca)<=0.1) {
		  float dcaSig = TMath::Abs(dca)/GTRACK333_dcaerr0[j];
		  if(sign>=0) {
		    h_dcaSig_p->Fill(dcaSig);
		    h_dcaSig->Fill(dcaSig);
		  }
		  else {
		    h_dcaSig_n->Fill(dcaSig);
		    h_dcaSig->Fill(-dcaSig);
		  }
		}
	      }
	    }
	  }
	}
      }
    } // loop over jets

  } //entries

  TFile file(fileName,"recreate");

  h_ntrk_all->Write();
  h_ntrk_pv->Write();
  h_ntrk_kpv->Write();
  h_ntrk_cor->Write();

  h_dca->Write();
  h_dcaSig->Write(); h_dcaSig_p->Write(); h_dcaSig_n->Write();

  file.Close();
}




///---------------------------------------------------
/// Kshort/Lambda studies
///
///---------------------------------------------------
void V0Analysis::AnalyzeV0(char* fileName)
{
  cout << endl << endl;
  cout << "K0/Lambda analysis" << endl;
  cout << endl;
  cout << "Interaction Point = " << IP_x << ", " << IP_y << endl;
  cout << "R cut = " << R << endl;
  cout << "Colinearity angle = " << CA << endl;
  cout << "Chi-square = " << CHI2 << endl;
  cout << "Veto Lambda in K0 = " << vetoLambda << endl;
  cout << "Veto K0 in Lambda = " << vetoK0 << endl;
  cout << endl;
  cout << fChain->GetEntries() << " Events to process " << endl;
  cout << endl;


  // 1- Kshorts
  // -----------------------------------------------------

  float minMass = 0.2; float maxMass = 1.2;
  float div = 0.01;
  int nBins = (maxMass - minMass)/div; 

  //decay length cuts
  TH1F* hm_opos_p_035 = new TH1F("hm_opos_p_035","",nBins,minMass,maxMass);
  TH1F* hm_opos_p_075 = new TH1F("hm_opos_p_075","",nBins,minMass,maxMass);
  TH1F* hpull_opos_p_035=new TH1F("hpull_opos_p_035","",nBins,-5,5);
  TH1F* hpull_opos_p_075=new TH1F("hpull_opos_p_075","",nBins,-5,5);

  TH1F* hm_opos_n_035 = new TH1F("hm_opos_n_035","",nBins,minMass,maxMass);
  TH1F* hm_opos_n_075 = new TH1F("hm_opos_n_075","",nBins,minMass,maxMass);
  TH1F* hpull_opos_n_035=new TH1F("hpull_opos_n_035","",nBins,-5,5);
  TH1F* hpull_opos_n_075=new TH1F("hpull_opos_n_075","",nBins,-5,5);

  TH1F* hm_same_p_035 = new TH1F("hm_same_p_035","",nBins,minMass,maxMass);
  TH1F* hm_same_p_075 = new TH1F("hm_same_p_075","",nBins,minMass,maxMass);
  TH1F* hm_same_n_035 = new TH1F("hm_same_n_035","",nBins,minMass,maxMass);
  TH1F* hm_same_n_075 = new TH1F("hm_same_n_075","",nBins,minMass,maxMass);

  //Properties
  TH1F*  hlife035_s = new TH1F("hlife035_s","",200,-2.5,2.5);
  TH1F*  hlife075_s = new TH1F("hlife075_s","",200,-2.5,2.5);
  TH1F*  hlife035_b = new TH1F("hlife035_b","",200,-2.5,2.5);
  TH1F*  hlife075_b = new TH1F("hlife075_b","",200,-2.5,2.5);

  TH1F* hlifetime075_s = new TH1F("hlifetime075_s","",200,-2.5,2.5);
  TH1F* hlifetime075_b = new TH1F("hlifetime075_b","",200,-2.5,2.5);

  TH1F*  hpt_035_s= new TH1F("hpt_035_s","",80,0.5,3.0);
  TH1F*  hpt_075_s= new TH1F("hpt_075_s","",80,0.5,3.0);
  TH1F*  hpt_035_b= new TH1F("hpt_035_b","",80,0.5,3.0);
  TH1F*  hpt_075_b= new TH1F("hpt_075_b","",80,0.5,3.0);

  TH1F*  hpipt_035_s = new TH1F("hpipt_035_s","",100,0.3,3.0);
  TH1F*  hpipt_035_b = new TH1F("hpipt_035_b","",100,0.3,3.0);
  TH1F*  hpipt_075_s = new TH1F("hpipt_075_s","",100,0.3,3.0);
  TH1F*  hpipt_075_b = new TH1F("hpipt_075_b","",100,0.3,3.0);

  // Individual cut analysis. Fix mass and loose 1 variable at the time
  TH1F*  hm_noCA = new TH1F("hm_noCA","",nBins,minMass,maxMass);
  TH1F*  h_noCA_s= new TH1F("h_noCA_s","",200,0.990,1.0);
  TH1F*  h_noCA_b= new TH1F("h_noCA_b","",200,0.990,1.0);

  TH1F*  hm_noIP= new TH1F("hm_noIP","",nBins,minMass,maxMass);
  TH1F*  h_noIP_s= new TH1F("h_noIP_s","",100,0,1.0);
  TH1F*  h_noIP_b= new TH1F("h_noIP_b","",100,0,1.0);
  TH2F*  h_noIP2_s= new TH2F("h_noIP2_s","",100,0,1.0,100,0,1.0);
  TH2F*  h_noIP2_b= new TH2F("h_noIP2_b","",100,0,1.0,100,0,1.0);

  TH1F*  hm_noDL = new TH1F("hm_noDL","",nBins,minMass,maxMass);
  TH1F*  h_noDL_s= new TH1F("h_noDL_s","",75,-R,R);
  TH1F*  h_noDL_b= new TH1F("h_noDL_b","",75,-R,R);
  TH1F*  hlife_s = new TH1F("hlife_s","",200,-2.5,2.5);
  TH1F*  hlife_b = new TH1F("hlife_b","",200,-2.5,2.5);

  //Effect of 1 cut at the time
  TH1F* hm1 = new TH1F("hm1","Opposite charged tracks",nBins,minMass,maxMass);
  TH1F* hm2 = new TH1F("hm2","Colinearity > 0.9999",nBins,minMass,maxMass);
  TH1F* hm3 = new TH1F("hm3","Decay lenght> 0.75",nBins,minMass,maxMass);
  TH1F* hm4 = new TH1F("hm4","|dca|>0.1cm",nBins,minMass,maxMass);
  TH1F* hm5 = new TH1F("hm5","Colinearity Angle > 0.9999",nBins,minMass,
		       maxMass);
  TH1F* hm6 = new TH1F("hm6","Decay lenght> 0.35",nBins,minMass,maxMass);
  TH1F* hm7 = new TH1F("hm7","#Delta#phi> 0.2",nBins,minMass,maxMass);



  // 2- Lambdas
  // -----------------------------------------------------
  minMass = 1.0; maxMass = 1.3;
  div = 0.004; //4 MeV
  nBins = (maxMass - minMass)/div;

  TH1F*  hl_50= new TH1F("hl_50","",nBins,minMass,maxMass);
  TH1F*  hall_50= new TH1F("hall_50","",nBins,minMass,maxMass);
  TH1F*  hlpull_50= new TH1F("hlpull_50","",nBins,minMass,maxMass);
  TH1F*  hlwrong_50= new TH1F("hlwrong_50","",nBins,minMass,maxMass);
  TH1F*  ha_50= new TH1F("ha_50","",nBins,minMass,maxMass);
  TH1F*  hapull_50= new TH1F("hapull_50","",nBins,minMass,maxMass);
  TH1F*  hawrong_50= new TH1F("hawrong_50","",nBins,minMass,maxMass);

  TH1F*  hl_100= new TH1F("hl_100","",nBins,minMass,maxMass);
  TH1F*  hall_100= new TH1F("hall_100","",nBins,minMass,maxMass);
  TH1F*  hlpull_100= new TH1F("hlpull_100","",nBins,minMass,maxMass);
  TH1F*  hlwrong_100= new TH1F("hlwrong_100","",nBins,minMass,maxMass);
  TH1F*  ha_100= new TH1F("ha_100","",nBins,minMass,maxMass);
  TH1F*  hapull_100= new TH1F("hapull_100","",nBins,minMass,maxMass);
  TH1F*  hawrong_100= new TH1F("hawrong_100","",nBins,minMass,maxMass);

  TH2F* h_dalitz = new TH2F("h_dalitz","",100,0.2,1.3,105,0.2,1.3);


  float px1,px2,py1,py2,pz1,pz2,p1,p2,q1,q2,dedx1,dedx2,dca1,dca2;  
  float pt1, pt2;
  int nsmt1,ncft1,nsmt2,ncft2;
  int id1, id2;
  bool select1, select2;


  int aa = 1000; int nn = 0;
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());  

  cout << nentries << " Events to process " << endl;

  Int_t nbytes = 0, nb = 0;
  for (Int_t jentry=1672; jentry<1673;jentry++) {
    //for (Int_t jentry=0; jentry<nentries;jentry++) {
    Int_t ientry = LoadTree(jentry); 
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    if(jentry/aa>nn) cout << "Event #" << jentry << endl;
    nn = jentry/aa;
    

       
    //loop over V0s
    for(int i=0; i<V0_V0_nv0; i++) {
      
      float dl    = decayLength(IP_x, IP_y, 0.0, i);
      float lifeTime = properLifeTime(IP_x, IP_y, 0.0, i);
      float r     = TMath::Abs(decayLength(0.0,0.0,0.0,i));
      float oa    = colinearityAngle(IP_x,IP_y,0.0,i);

      float K0mass  = V0_V0_mpi[i];
      float K0massError  = V0_V0_smpi[i];
      float K0massPull = (K0mass-K0MASS)/K0massError;

      float lambdaMass = V0_V0_ml[i];
      float lambdaMassError = V0_V0_sml[i];
      float lambdaMassPull =(lambdaMass-LAMBDAMASS)/lambdaMassError;

      float antiLambdaMass = V0_V0_mlb[i];
      float antiLambdaMassError= V0_V0_smlb[i];
      float antiLambdaMassPull=(antiLambdaMass-LAMBDAMASS)/antiLambdaMassError;

      float chi2  = V0_V0_chisq[i];
      float ch    = V0_V0_ch[i];
      float pt    = TMath::Sqrt(V0_V0_px[i]*V0_V0_px[i]+
				V0_V0_py[i]*V0_V0_py[i]);
      
      get_kinematics(i,q1,px1,py1,pz1,p1,dedx1,nsmt1,ncft1,dca1,id1,
		     q2,px2,py2,pz2,p2,dedx2,nsmt2,ncft2,dca2,id2);


      TVector3 vv(V0_V0_x[i],V0_V0_y[i],0.0);
      float phidir = vv.Phi();
      TVector3 pp1(px1,py1,pz1);
      TVector3 pp2(px2,py2,pz2);
      float deltaphi = pp1.DeltaPhi(pp2);

 
      pt1 = TMath::Sqrt(p1*p1-pz1*pz1);
      pt2 = TMath::Sqrt(p2*p2-pz2*pz2);

      select1 = pt1>=0.3 && (nsmt1+ncft1)>=11;
      select2 = pt2>=0.3 && (nsmt2+ncft2)>=11;

      //basic cuts: good tracks within silicon detector
      if(select1 && select2 && r<R) {
	
	//MASS HYPOTHESIS for Lambda/AntiLambda
	bool isLambda = false;
	if(pt1>pt2 && q1>0) isLambda = true;
	if(pt2>pt1 && q2>0) isLambda = true;

	bool isL0 = false;
	if( (isLambda && lambdaMass>1.113 && lambdaMass<1.121) || 
	    (!isLambda && antiLambdaMass>1.113 && antiLambdaMass<1.121) ) 
	  isL0 = true; 
	
	bool isKs = false; bool isSideband = false;
	if(K0mass>=0.475 && K0mass<=0.525) isKs = true;
	if(K0mass>=0.7 && K0mass<=1.2) isSideband = true;
	
 

	// ----------------------------------------
	// K0 analysis
	// ----------------------------------------
	if(!vetoLambda || (vetoLambda && !isL0)) {

	  // 1- Study effect of different cuts on K0 mass peak
	  // ------------------------------------------------
	  if(ch==-1 && chi2<=CHI2) {
	    hm1->Fill(K0mass);
	    if(oa>0.9999) {
	      hm2->Fill(K0mass);
	      if(dl>0.75) {
		hm3->Fill(K0mass);
		if(deltaphi>0.2) hm7->Fill(K0mass);
	      }
	    }
	    if(dca1>0.1 && dca2>0.1) {
	      hm4->Fill(K0mass);
	      if(oa>0.9999) {
		hm5->Fill(K0mass);
		if(dl>0.75) {
		  hm6->Fill(K0mass);
		}
	      }
	    }
	  }
	  
	  // 2- Decay length > 0.35cm
	  // ------------------------------------------------
	  if(TMath::Abs(dl)>=0.35 && TMath::Abs(oa)>=CA && chi2<=CHI2 &&
	     dca1>0.1 && dca2>0.1) {

	    cout << jentry << "  " << id1 << " " << id2 << " - " << K0mass<<endl;
	    cout << px1 << " " << py1 << "    " << px2 << " " << py2 << endl;
	    cout << dca1 << " " << dca2 << endl;

	    if(ch==-1 && dl>0 && oa>0) {  // ct > 0
	      hm_opos_p_035->Fill(K0mass);
	      hpull_opos_p_035->Fill(K0massPull);
	      if(isKs) {
		hpt_035_s->Fill(pt);
		hlife035_s->Fill(lifeTime);
		hpipt_035_s->Fill(pt1); hpipt_035_s->Fill(pt2);
		cout << "K0 candidate event #" << jentry << endl;
	      } 
	      if(isSideband) {
		hpt_035_b->Fill(pt);
		hpipt_035_b->Fill(pt1); hpipt_035_b->Fill(pt2);
		hlife035_b->Fill(lifeTime);
	      }
	    }
	    if(ch==-1 && dl<0 && oa<0) {  //ct < 0
	      hm_opos_n_035->Fill(K0mass);
	      hpull_opos_n_035->Fill(K0massPull);
	    }
	    if(ch==1 && dl>0 && oa>0) {  // ct > 0 wrong sign
	      hm_same_p_035->Fill(K0mass);
	    }
	    if(ch==1 && dl<0 && oa<0) {  //ct < 0  wrong sign
	      hm_same_n_035->Fill(K0mass);
	    }
	  }
	  
	  // 3- Decay length > 0.75cm
	  // ------------------------------------------------
	  if(TMath::Abs(dl)>=0.75 && TMath::Abs(oa)>=CA && chi2<=CHI2 &&
	     dca1>0.1 && dca2>0.1) {
	    if(ch==-1 && isKs) hlifetime075_s->Fill(lifeTime);
	    if(ch==-1 && isSideband) hlifetime075_b->Fill(lifeTime);
	    
	    if(ch==-1 && dl>0 && oa>0) {  // ct > 0
	      hm_opos_p_075->Fill(K0mass);
	      hpull_opos_p_075->Fill(K0massPull);
	      if(isKs) {
	      hpt_075_s->Fill(pt);
	      hlife075_s->Fill(lifeTime);
	      hpipt_075_s->Fill(pt1); hpipt_075_s->Fill(pt2);
	      } 
	      if(isSideband) {
		hpt_075_b->Fill(pt);
		hlife075_b->Fill(lifeTime);
		hpipt_075_b->Fill(pt1); hpipt_075_b->Fill(pt2);
	      }
	    }
	    if(ch==-1 && dl<0 && oa<0) {  //ct < 0
	      hm_opos_n_075->Fill(K0mass);
	      hpull_opos_n_075->Fill(K0massPull);
	    }
	    if(ch==1 && dl>0 && oa>0) {  // ct > 0 wrong sign
	      hm_same_p_075->Fill(K0mass);
	    }
	    if(ch==1 && dl<0 && oa<0) {  //ct < 0  wrong sign
	      hm_same_n_075->Fill(K0mass);
	    }
	  }
	  
	  // 4- Study colinearity angle cut
	  // ------------------------------------------------
	  if(dl>=0.75 && chi2<=CHI2 && ch==-1 && dca1>0.1 && dca2>0.1) {  
	    hm_noCA->Fill(K0mass);
	    if(isKs) h_noCA_s->Fill(oa);
	    if(isSideband) h_noCA_b->Fill(oa);	    
	  }
	  
	  // 5- Study decay length cut
	  // ------------------------------------------------
	  if(TMath::Abs(oa)>=0.9999 && chi2<=CHI2 && ch==-1 && dca1>0.1 
	     && dca2>0.1) {  
	    hm_noDL->Fill(K0mass);
	    if(isKs) {
	      h_noDL_s->Fill(dl);
	      hlife_s->Fill(lifeTime);
	    } 
	    if(isSideband) {
	      h_noDL_b->Fill(dl);
	      hlife_b->Fill(lifeTime);
	    }
	  }
	  
	  // 6- Study Impact Parameter cut
	  // ------------------------------------------------
	  if(dl>=0.75 && oa>=CA && ch==-1 && chi2<=CHI2) {  
	    hm_noIP->Fill(K0mass);
	    if(isKs) {
	      h_noIP_s->Fill(dca1);
	      h_noIP_s->Fill(dca2);
	      h_noIP2_s->Fill(dca1,dca2);
	    } 
	    if(isSideband) {
	      h_noIP_b->Fill(dca1);
	      h_noIP_b->Fill(dca2);
	      h_noIP2_b->Fill(dca1,dca2);
	    }
	  }
	} //veto lambda
	 
	  
	  
	// ----------------------------------------
	// Lambda analysis
	// ----------------------------------------
	if(!vetoK0 || (vetoK0 && !isKs)) {
	  
	  // Decay length > 0.5cm
	  if(dl>=0.5 && oa>=0.9999 && chi2<=15 && ch==-1) {
	    if(isLambda) {
	      hl_50->Fill(lambdaMass);
	      hall_50->Fill(lambdaMass);
	      hlpull_50->Fill(lambdaMassPull);
	      hlwrong_50->Fill(antiLambdaMass);
	    } else {
	      ha_50->Fill(antiLambdaMass);
	      hall_50->Fill(antiLambdaMass);
	      hapull_50->Fill(antiLambdaMassPull);
	      hawrong_50->Fill(lambdaMass);
	    }
	  }

	  // Decay length > 1.0cm
	  if(dl>=1.0 && oa>=0.9999 && chi2<=15 && ch==-1) {
	    
	    if(isLambda) {
	      hl_100->Fill(lambdaMass);
	      hall_100->Fill(lambdaMass);
	      hlpull_100->Fill(lambdaMassPull);
	      hlwrong_100->Fill(antiLambdaMass);
	    } else {
	      ha_100->Fill(antiLambdaMass);
	      hall_100->Fill(antiLambdaMass);
	      hapull_100->Fill(antiLambdaMassPull);
	      hawrong_100->Fill(lambdaMass);
	    }
	  }
	} //veto K0
	
	
	if(dl>=0.75 && oa>=0.9999 && chi2<=15 && ch==-1 && isLambda)
	  h_dalitz->Fill(K0mass,lambdaMass);
       

      } // basic selection 
    }  // V0 loop  
  }//nentries
  
  
  //save histogramas
  TFile file(fileName,"recreate");

  hm1->Write(); hm2->Write(); hm3->Write(); 
  hm4->Write(); hm5->Write(); hm6->Write();
  hm7->Write(); 

  hm_opos_p_035->Write(); hpull_opos_p_035->Write();
  hm_opos_p_075->Write(); hpull_opos_p_075->Write();
  hm_opos_n_035->Write(); hpull_opos_n_035->Write();
  hm_opos_n_075->Write(); hpull_opos_n_075->Write();
  hm_same_p_035->Write(); hm_same_n_035->Write();
  hm_same_p_075->Write(); hm_same_n_075->Write();

  hlife035_s->Write();  hlife035_b->Write();
  hlife075_s->Write();  hlife075_b->Write();

  hpt_035_s->Write(); hpt_035_b->Write();
  hpt_075_s->Write(); hpt_075_b->Write();

  hpipt_035_s->Write(); hpipt_035_b->Write();
  hpipt_075_s->Write(); hpipt_075_b->Write();

  hm_noCA->Write();
  h_noCA_s->Write(); h_noCA_b->Write();

  hm_noDL->Write(); 
  h_noDL_s->Write(); h_noDL_b->Write(); 
  hlife_s->Write();hlife_b->Write();

  hm_noIP->Write();
  h_noIP_s->Write(); h_noIP_b->Write();
  h_noIP2_s->Write(); h_noIP2_b->Write();



  hl_50->Write();
  hall_50->Write();
  hlpull_50->Write();
  hlwrong_50->Write();
  ha_50->Write();
  hapull_50->Write();
  hawrong_50->Write();
   
  hl_100->Write();
  hall_100->Write();
  hlpull_100->Write();
  hlwrong_100->Write();
  ha_100->Write();
  hapull_100->Write();
  hawrong_100->Write();
  
  h_dalitz->Write();

  file.Close();
}





///---------------------------------------------------
/// J/Psi studies
///
///---------------------------------------------------
void V0Analysis::AnalyzeJpsi(char* fileName)
{
  CHI2 = 5;

  cout << endl << endl;
  cout << "JPsi analysis" << endl;
  cout << endl;
  cout << "Interaction Point = " << "-> selected PVSEL333" << endl;
  cout << "R cut = " << R << endl;
  cout << "Colinearity angle = " << CA << endl;
  cout << "Chi-square = " << CHI2 << endl;
  cout << endl;
  cout << fChain->GetEntries() << " Events to process " << endl;
  cout << endl;

  //nice peak can be obtained with masses from 2.0 to 4.0 with 0.05 interval 
  float minMass = 2.0; float maxMass = 4;
  float div = 0.05;
  int nBins = (maxMass - minMass)/div; 

  //pT  cuts
  TH1F* hm_pt1 = new TH1F("hm_pt1","",nBins,minMass,maxMass);
  TH1F* hm_pt2 = new TH1F("hm_pt2","",nBins,minMass,maxMass);
  TH1F* hm_pt3 = new TH1F("hm_pt3","",nBins,minMass,maxMass);

  div = 0.04;
  nBins = (maxMass - minMass)/div; 
  TH1F* hm_pt1_i = new TH1F("hm_pt1_i","",nBins,minMass,maxMass);
  TH1F* hm_pt2_i = new TH1F("hm_pt2_i","",nBins,minMass,maxMass);
  TH1F* hm_pt3_i = new TH1F("hm_pt3_i","",nBins,minMass,maxMass);

  div = 0.03;
  nBins = (maxMass - minMass)/div; 
  TH1F* hm_pt1_ii = new TH1F("hm_pt1_ii","",nBins,minMass,maxMass);
  TH1F* hm_pt2_ii = new TH1F("hm_pt2_ii","",nBins,minMass,maxMass);
  TH1F* hm_pt3_ii = new TH1F("hm_pt3_ii","",nBins,minMass,maxMass);



  //Properties
  // 1- Decay length
  TH1F*  hdl_s = new TH1F("hdl_s","",200,-0.5,0.5);
  TH1F*  hdl_b = new TH1F("hdl_b","",200,-0.5,0.5);

  // 2- Impact parameter
  TH1F*  hdca_s = new TH1F("hddca_s","",200,-0.1,0.1);
  TH1F*  hdca_b = new TH1F("hddca_b","",200,-0.1,0.1);


  // 2- Lifetime
  TH1F*  hlife_s = new TH1F("hlife_s","",200,-1.0,1.0);
  TH1F*  hlife_b = new TH1F("hlife_b","",200,-1.0,1.0);

  // 3- Transverse momentum
  TH1F*  hpt_s= new TH1F("hpt_s","",80,0.5,3.0);
  TH1F*  hpt_b= new TH1F("hpt_b","",80,0.5,3.0);

  // 4- Muon pT
  TH1F*  hmupt_s = new TH1F("hmupt_s","",100,0.3,3.0);
  TH1F*  hmupt_b = new TH1F("hmupt_b","",100,0.3,3.0);


  TH1F*  h_dphi = new TH1F("h_dphi","#Delta#phi di-muon - vertex",100,0,3.15);
  

  float px1,px2,py1,py2,pz1,pz2,p1,p2,q1,q2,dedx1,dedx2,dca1,dca2;  
  float pt1, pt2;
  int nsmt1,ncft1,nsmt2,ncft2;
  int id1, id2;


  int aa = 1000; int nn = 0;
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());  
  Int_t nbytes = 0, nb = 0;
  for (Int_t jentry=0; jentry<nentries;jentry++) {
    Int_t ientry = LoadTree(jentry); 
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    if(jentry/aa>nn) cout << "Event #" << jentry << endl;
    nn = jentry/aa;
    
    //Primary vertex
    bool pv = false;
    if(PVSEL333_psnvtx333>0 && PVSEL333_psntrk[0]>=2) pv = true;
    float PV_x = PVSEL333_psx[0];
    float PV_y = PVSEL333_psy[0];
    float PV_z = PVSEL333_psz[0];
    
    //Muons
    if(Muonid_Nmuonid<2) continue;
    TVector3 mu1(Muonid_Idpx[0],Muonid_Idpy[0],Muonid_Idpz[0]);
    TVector3 mu2(Muonid_Idpx[1],Muonid_Idpy[1],Muonid_Idpz[1]);
    TVector3 diMuon = mu1+mu2;
    
    //loop over V0s
    for(int i=0; i<V0_V0_nv0; i++) {
      
      float dl    = decayLength(PV_x, PV_y, 0.0, i);
      float lifeTime = properLifeTime(PV_x, PV_y, 0.0, i);
      float r     = TMath::Abs(decayLength(0.0,0.0,0.0,i));
      float oa    = colinearityAngle(PV_x,PV_y,0.0,i);

      float MUmass  = V0_V0_mmu[i];
      float MUmassError  = V0_V0_smmu[i];
      float MUmassPull = (MUmass-JPSIMASS)/MUmassError;

      float chi2  = V0_V0_chisq[i];
      float ch    = V0_V0_ch[i];
      float pt    = TMath::Sqrt(V0_V0_px[i]*V0_V0_px[i]+
				V0_V0_py[i]*V0_V0_py[i]);
      
      get_kinematics(i,q1,px1,py1,pz1,p1,dedx1,nsmt1,ncft1,dca1,id1,
		     q2,px2,py2,pz2,p2,dedx2,nsmt2,ncft2,dca2,id2);
      
      
      pt1 = TMath::Sqrt(p1*p1-pz1*pz1);
      pt2 = TMath::Sqrt(p2*p2-pz2*pz2);
      
      TVector3 vtx(V0_V0_px[i],V0_V0_py[i],V0_V0_pz[i]);
      float deltaPhi = diMuon.DeltaPhi(vtx);
 
      //basic cuts: good opposite charge tracks within silicon detector
      if(ch==-1 && chi2<=CHI2 && nsmt1+ncft1>=11 && nsmt2+ncft2>=11) {

	
	if((pt1>=3 && pt2>=2) || (pt1>=2 && pt2>=3)) h_dphi->Fill(deltaPhi);

	if(deltaPhi<0.3) {
	
	  bool isJpsi = false; bool isSideband = false;
	  if((pt1>=3 && pt2>=2) || (pt1>=2 && pt2>=3)) {
	    if(MUmass>=3.06-0.1 && MUmass<=3.06+0.1) isJpsi = true;
	    if((MUmass>=3.5 && MUmass<=3.7) || (MUmass>=2.4 && MUmass<=2.6)) 
	      isSideband = true;
	  }
	  
	  if( (pt1>=2 && pt2>=2.5) ||(pt1>=2.5 && pt2>=2.0) ) {
	    hm_pt1->Fill(MUmass);
	    hm_pt1_i->Fill(MUmass);
	    hm_pt1_ii->Fill(MUmass);
	  }
	  
	  if(pt1>=2 && pt2>=2) {
	    hm_pt2->Fill(MUmass);
	    hm_pt2_i->Fill(MUmass);
	    hm_pt2_ii->Fill(MUmass);
	  }
	  
	  if((pt1>=3 && pt2>=2) || (pt1>=2 && pt2>=3)) {
	    hm_pt3->Fill(MUmass);
	    hm_pt3_i->Fill(MUmass);
	    hm_pt3_ii->Fill(MUmass);
	  }
	  
	  
	  //properties
	  if(isJpsi && pv) {
	    hdl_s->Fill(dl);
	    hlife_s->Fill(lifeTime);
	    hpt_s->Fill(pt);
	    hdca_s->Fill(dca1); hdca_s->Fill(dca2); 
	    hmupt_s->Fill(pt1); hmupt_s->Fill(pt2);
	  }
	  if(isSideband) {
	    hdl_b->Fill(dl);
	    hlife_b->Fill(lifeTime);
	    hpt_b->Fill(pt);
	    hdca_b->Fill(dca1); hdca_b->Fill(dca2); 
	    hmupt_b->Fill(pt1); hmupt_b->Fill(pt2);
	  }
	  
	} //delta-phi
      } // basic selection
    } // V0 loop
  } // entries
  
  //save histogramas
  TFile file(fileName,"recreate");

  h_dphi->Write(); 

  hm_pt1->Write(); hm_pt2->Write(); hm_pt3->Write(); 
  hm_pt1_i->Write(); hm_pt2_i->Write(); hm_pt3_i->Write(); 
  hm_pt1_ii->Write(); hm_pt2_ii->Write(); hm_pt3_ii->Write(); 

  hdl_s->Write(); hdl_b->Write();
  hlife_s->Write(); hlife_b->Write(); 
  hpt_s->Write(); hpt_b->Write(); 
  hdca_s->Write(); hdca_b->Write(); 
  hmupt_s->Write(); hmupt_b->Write(); 
}



///---------------------------------------------------
/// J/Psi studies
///
///---------------------------------------------------
void V0Analysis::AnalyzeBJpsiKs(char* fileName)
{
  CHI2 = 10;

  cout << endl << endl;
  cout << "B->JPsi+Ks analysis" << endl;
  cout << endl;
  cout << "Interaction Point = " << "-> selected PVSEL333" << endl;
  cout << "R cut = " << R << endl;
  cout << "Colinearity angle = " << CA << endl;
  cout << "Chi-square = " << CHI2 << endl;
  cout << endl;
  cout << fChain->GetEntries() << " Events to process " << endl;
  cout << endl;


  float minMass = 2.0; float maxMass = 8;
  float div = 0.05;
  int nBins = (maxMass - minMass)/div; 

  //pT  cuts
  TH1F* hm_b = new TH1F("hm_b","",nBins,minMass,maxMass);


  float px1,px2,py1,py2,pz1,pz2,p1,p2,q1,q2,dedx1,dedx2,dca1,dca2;  
  float pt1, pt2;
  int nsmt1,ncft1,nsmt2,ncft2;
  int id1, id2;


  int aa = 1000; int nn = 0;
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());  
  Int_t nbytes = 0, nb = 0;
  for (Int_t jentry=0; jentry<nentries;jentry++) {
    Int_t ientry = LoadTree(jentry); 
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    if(jentry/aa>nn) cout << "Event #" << jentry << endl;
    nn = jentry/aa;
    
    //Primary vertex
    bool pv = false;
    if(PVSEL333_psnvtx333>0 && PVSEL333_psntrk[0]>=2) pv = true;
    float PV_x = PVSEL333_psx[0];
    float PV_y = PVSEL333_psy[0];
    float PV_z = PVSEL333_psz[0];
    
    //Muons
    if(Muonid_Nmuonid<2) continue;
    TVector3 mu1(Muonid_Idpx[0],Muonid_Idpy[0],Muonid_Idpz[0]);
    TVector3 mu2(Muonid_Idpx[1],Muonid_Idpy[1],Muonid_Idpz[1]);
    TVector3 diMuon = mu1+mu2;
    
    //1- Select J/Psi
    for(int i=0; i<V0_V0_nv0; i++) {
      
      float dl_jpsi    = decayLength(PV_x, PV_y, 0.0, i);
      float oa_jpsi    = colinearityAngle(PV_x,PV_y,0.0,i);
      
      float MUmass  = V0_V0_mmu[i];

      float chi2_jpsi  = V0_V0_chisq[i];
      float ch_jpsi    = V0_V0_ch[i];
      float pt_jpsi    = TMath::Sqrt(V0_V0_px[i]*V0_V0_px[i]+
				     V0_V0_py[i]*V0_V0_py[i]);
      
      get_kinematics(i,q1,px1,py1,pz1,p1,dedx1,nsmt1,ncft1,dca1,id1,
		     q2,px2,py2,pz2,p2,dedx2,nsmt2,ncft2,dca2,id2);
      
      
      pt1 = TMath::Sqrt(p1*p1-pz1*pz1);
      pt2 = TMath::Sqrt(p2*p2-pz2*pz2);
      
      TVector3 vtx(V0_V0_px[i],V0_V0_py[i],V0_V0_pz[i]);
      float deltaPhi = diMuon.DeltaPhi(vtx);
 
      //basic cuts: good opposite charge tracks within silicon detector
      if(ch_jpsi==-1 && chi2_jpsi<=CHI2 && nsmt1+ncft1>=11 && nsmt2+ncft2>=11 
	 && ((pt1>=2 && pt2>=2) || (pt1>=2 && pt2>=2)) && deltaPhi<0.3 &&
	 (MUmass>=3.06-0.1 && MUmass<=3.06+0.1)) { // JPsi found
	
	float JPSI_x = V0_V0_x[i];
	float JPSI_y = V0_V0_y[i];
	float JPSI_z = V0_V0_z[i];
	
	//2- Select K0 pointing to jpsi (dca cut missing)
	for(int j=0; j<V0_V0_nv0; j++) {
	  
	  float dl_k0    = decayLength(JPSI_x, JPSI_y, 0.0, j);
	  float oa_k0    = colinearityAngle(JPSI_x,JPSI_y,0.0,j);
	  
	  float K0mass  = V0_V0_mpi[j];
	  float chi2_k0  = V0_V0_chisq[j];
	  float ch_k0    = V0_V0_ch[j];
	  float pt_k0    = TMath::Sqrt(V0_V0_px[j]*V0_V0_px[j]+
				       V0_V0_py[j]*V0_V0_py[j]);
	  
	  get_kinematics(i,q1,px1,py1,pz1,p1,dedx1,nsmt1,ncft1,dca1,id1,
			 q2,px2,py2,pz2,p2,dedx2,nsmt2,ncft2,dca2,id2);
	  
	  pt1 = TMath::Sqrt(p1*p1-pz1*pz1);
	  pt2 = TMath::Sqrt(p2*p2-pz2*pz2);
	  
	  bool select1 = pt1>=0.3 && (nsmt1+ncft1)>=11;
	  bool select2 = pt2>=0.3 && (nsmt2+ncft2)>=11;
	  
	  //basic cuts: good tracks within silicon detector
	  if(select1 && select2 && dl_k0>=0.1 && oa_k0>=CA && chi2_k0<CHI2 && 
	     ch_k0==-1 && (K0mass>=0.475 && K0mass<=0.525)) {
	    
	    cout << "JPsi+Ks event #" << jentry << endl;
	    
	    TLorentzVector jpsiVtx; 
	    jpsiVtx.SetXYZM(V0_V0_px[i],V0_V0_py[i],V0_V0_pz[i],V0_V0_mmu[i]);
	    TLorentzVector k0Vtx;
	    k0Vtx.SetXYZM(V0_V0_px[j],V0_V0_py[j],V0_V0_pz[j],V0_V0_mpi[j]);
	    
	    TLorentzVector B = jpsiVtx + k0Vtx;
	    cout << "B->Jpsi+Ks vertex candidate :";
	    cout << B(0) << ", " << B(1) << ", " << B(2) << ", " << B(3)<<endl;
	    cout << "Mass = " << B.M() << endl;
	    hm_b->Fill(B.M());
	  } 
	} // K0
      } 
    } // JPsi
  } // entries
  
  //save histogramas
  TFile file(fileName,"recreate");

  hm_b->Write();

  file.Close();
}







///---------------------------------------------------
/// b-tagging studies
///
///---------------------------------------------------
void V0Analysis::AnalyzeBTagging_jpsi(char* fileName, bool debug)
{
  CHI2 = 5;


  // JPSI
  TH1F* h_njets = new TH1F("h_njets","Jet multiplicity",5,-0.5,4.5);
  TH2F* h_jetaphi = new TH2F("h_jetaphi","Cal-Jet #eta vs. #phi",
			     100,-3,3,100,-3.15,3.15);
  TH1F* h_jpsiJet = new TH1F("h_jpsiJet","#DeltaR J/#psi - Jet",100,0,1);
  TH1F* h_aJet = new TH1F("h_aJet","#Delta#phi J/#psi - away Jet",100,0,1);
  TH1F* h_njpsijets = new TH1F("h_njpsijets","Inclusive J/#psi + Jet multiplicity",3,-0.5,2.5);
  TH1F* h_njpsijets_matched = new TH1F("h_njpsijets_matched","Inclusive J/#psi + Jet multiplicity",3,-0.5,2.5);




  // SECONDARY VERTEX away-jet
  TH1F* hvtx_dphi = new TH1F("hvtx_dphi","#Delta#phi(SVX,J/#psi)",100,0,3.15);
  TH1F* hvtx_dl = new TH1F("hvtx_dl","Secondary Vertex decay length",
			   100,-1,1);
  hvtx_dl->SetXTitle("cm");
  TH1F* hvtx_mass = new TH1F("hvtx_mass","Secondary Vertex mass",100,0,8);
  hvtx_mass->SetXTitle("GeV/c^{2}");
  TH1F* hvtx_mult = new TH1F("hvtx_mult","Secondary Vertex multiplicity",
			     5,1.5,6.5);

  TH1F* hvtx_dl_away = new TH1F("hvtx_dl_away","Secondary Vertex decay length [away jet]",100,-1,1);
  hvtx_dl_away->SetXTitle("cm");


  // IMPACT-PARAMETER jpsi
  TH1F* h_dca_jpsi = new TH1F("h_dca_jpsi","J/#psi Signed Impact Parameter",
			      100,-0.2,0.2);
  h_dca_jpsi->SetXTitle("cm");
  TH1F* h_sdca_jpsi = new TH1F("h_sdca_jpsi","J/#psi Signed Impact Parameter Significance",100,-30,30);


  // IMPACT-PARAMETER away-jet
  TH1F* h_dcaSig_p = new TH1F("h_dcaSig_p","Signed Impact Parameter Significance [p_{T} > 1.5 GeV/c]",200,0,30);
  TH1F* h_dcaSig_n = new TH1F("h_dcaSig_n","Signed Impact Parameter Significance [p_{T} > 1.5 GeV/c]",200,0,30);




  float px1,px2,py1,py2,pz1,pz2,p1,p2,q1,q2,dedx1,dedx2,dca1,dca2;  
  float pt1, pt2;
  int nsmt1,ncft1,nsmt2,ncft2;
  int id1, id2;

  int njpsi0 = 0; int njpsi1 = 0; int njpsi2 = 0;
  int njpsijet = 0;
  int najet = 0;
  int njpsiANDajet = 0;
  int njpsiORajet = 0;


  int aa = 1000; int nn = 0;
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());  
  Int_t nbytes = 0, nb = 0;
  for (Int_t jentry=0; jentry<nentries;jentry++) {
    Int_t ientry = LoadTree(jentry); 
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    if(jentry/aa>nn) cout << "Event #" << jentry << endl;
    nn = jentry/aa;


    // 1-Primary vertex
    bool pv = false;
    if(PVSEL333_psnvtx333>0 && PVSEL333_psntrk[0]>=2) pv = true;
    if(!pv) continue;
    float PV_x = PVSEL333_psx[0];
    float PV_y = PVSEL333_psy[0];
    float PV_z = PVSEL333_psz[0];

    // 2-Jets
    TObjArray jets; readJets(jets);
    h_njets->Fill(jets.GetEntries());
    for(int y=0; y<jets.GetEntries(); y++) {
      TLorentzVector* jet = (TLorentzVector*) jets[y];
      h_jetaphi->Fill(jet->Phi(),jet->Eta());
     }
    


    // 3- Muons
    if(Muonid_Nmuonid<2) continue;
    TLorentzVector mu1(Muonid_Idpx[0],Muonid_Idpy[0],Muonid_Idpz[0],
		       Muonid_IdE[0]);
    TLorentzVector mu2(Muonid_Idpx[1],Muonid_Idpy[1],Muonid_Idpz[1],
		       Muonid_IdE[1]);
    TLorentzVector diMuon = mu1+mu2;
    

    //3- J/Psi
    int jpsi = -1; float ptmin = 0; TLorentzVector jpsiVTX;
    for(int i=0; i<V0_V0_nv0; i++) {
      float dl_jpsi    = decayLength(PV_x, PV_y, 0.0, i);
      float oa_jpsi    = colinearityAngle(PV_x,PV_y,0.0,i);
      float MUmass     = V0_V0_mmu[i];
      float chi2_jpsi  = V0_V0_chisq[i];
      float ch_jpsi    = V0_V0_ch[i];
      float pt_jpsi    = TMath::Sqrt(V0_V0_px[i]*V0_V0_px[i]+
				     V0_V0_py[i]*V0_V0_py[i]);
      
      get_kinematics(i,q1,px1,py1,pz1,p1,dedx1,nsmt1,ncft1,dca1,id1,
		     q2,px2,py2,pz2,p2,dedx2,nsmt2,ncft2,dca2,id2);
      pt1 = TMath::Sqrt(p1*p1-pz1*pz1); pt2 = TMath::Sqrt(p2*p2-pz2*pz2);
      
      TLorentzVector vtx; vtx.SetXYZM(V0_V0_px[i],V0_V0_py[i],V0_V0_pz[i],
				       V0_V0_mmu[i]);
      float deltaPhi = diMuon.DeltaPhi(vtx);
      
      //basic cuts
      if(ch_jpsi==-1 && chi2_jpsi<=CHI2 && nsmt1+ncft1>=11 && nsmt2+ncft2>=11 
	 && ((pt1>=3 && pt2>=2) || (pt1>=2 && pt2>=3)) && deltaPhi<0.3 &&
	 (MUmass>=3.06-0.1 && MUmass<=3.06+0.1)) { // JPsi found
	
	if(pt_jpsi>=ptmin) {
	  ptmin = pt_jpsi;
	  jpsi = i;
	}
      }
    }


    if(jpsi>-1) {
      njpsi0++;
      if(jets.GetEntries()>=1) njpsi1++;
      if(jets.GetEntries()>=2) njpsi2++;


      jpsiVTX.SetXYZM(V0_V0_px[jpsi],V0_V0_py[jpsi],
		      V0_V0_pz[jpsi],V0_V0_mmu[jpsi]);
      
      
      //define Jpsi-Jet 
      TLorentzVector* jpsiJET = 0; float min = 0.5;
      for(int y=0; y<jets.GetEntries(); y++) {
	TLorentzVector* jet = (TLorentzVector*) jets[y];
	float dr = jet->DeltaR(jpsiVTX);
	if(dr<min) {
	  min = dr;
	  jpsiJET = jet;
	}
      }
      if(jpsiJET) {
	h_jpsiJet->Fill(min);
	njpsijet++;
      }
      
      
      //define away-jet
      TLorentzVector* aJET = 0; min = 0.5;
      for(int y=0; y<jets.GetEntries(); y++) {
	TLorentzVector* jet = (TLorentzVector*) jets[y];
	float dphi = TMath::Abs(jet->DeltaPhi(jpsiVTX));
	if(TMath::Pi()-dphi<min) {
	  min = TMath::Pi()-dphi;
	  aJET = jet;
	}
      }
      if(aJET) {
	h_aJet->Fill(min);
	najet++;
      }

      if(jpsiJET || aJET) njpsiORajet++;
      if(jpsiJET && aJET) njpsiANDajet++;


      if(debug) {
	cout << "Jpsi    :" << jpsiVTX.Eta() << " " << jpsiVTX.Phi();
	cout << "  mass = " << jpsiVTX.M() << " " << V0_V0_mmu[jpsi] <<endl;
	if(jpsiJET) 
	  cout << "jpsi-jet:" << jpsiJET->Eta() << " " << jpsiJET->Phi()<<endl;
	if(aJET) 
	  cout << "a-jet   :" << aJET->Eta() << " " << aJET->Phi() << endl;
	cout << endl;
      }
    
    
      ///------------------------------------
      /// Secondary vertex analysis
      ///------------------------------------
      for(int u=0; u<SKSEL333_kssnvtx333; u++) {
	TLorentzVector vtx; 
	vtx.SetXYZM(SKSEL333_ksspx[u],SKSEL333_ksspy[u],
		    SKSEL333_ksspz[u],SKSEL333_kssmass[u]);
	float dl_svx  = decayLength_sksel(PV_x, PV_y, 0.0, u);
	float dphi = TMath::Abs(jpsiVTX.DeltaPhi(vtx));
	hvtx_dphi->Fill(dphi);
	hvtx_dl->Fill(dl_svx);
	hvtx_mass->Fill(vtx.M());
 	hvtx_mult->Fill(SKSEL333_kssntrk[u]);
	
	//look at vertices in away jet for different j/psi pTs
	if(dphi>=TMath::Pi()-0.5) {
	  hvtx_dl_away->Fill(dl_svx);
	}
      }
      
      
      ///------------------------------------
      /// Impact-Parameter analysis
      ///------------------------------------
      if(jpsiJET) {
 	float dca,sdca;
	get_kinematics(jpsi,q1,px1,py1,pz1,p1,dedx1,nsmt1,ncft1,dca1,id1,
		       q2,px2,py2,pz2,p2,dedx2,nsmt2,ncft2,dca2,id2);
	//track 1
	signedIP(id1,jpsiJET->Phi(),dca,sdca);
	h_dca_jpsi->Fill(dca); h_sdca_jpsi->Fill(sdca);

	//track 2
	signedIP(id2,jpsiJET->Phi(),dca,sdca);
	h_dca_jpsi->Fill(dca); h_sdca_jpsi->Fill(sdca);
      }


      if(aJET) {
	//match tracks to away-jet
	for(int j=0; j<GTRACK333_ngtrk333; j++) {
	  float pT = 1./GTRACK333_qoverpt[j];
	  int hits = GTRACK333_nsmt[j]+GTRACK333_ncft[j];
	  if(pT>1.5 && hits>=11) {
	    float phi = GTRACK333_phidir[j]; if(phi<0) phi += 2*TMath::Pi();
	    double px = pT*TMath::Cos(phi);
	    double py = pT*TMath::Sin(phi);
	    double pz = pT*GTRACK333_tanl[j];
	    double E = TMath::Sqrt(px*px+py*py+pz*pz+0.13956995*0.13956995);  
	    TLorentzVector trk(px,py,pz,E);
	    float deltaPhi = trk.Phi() - aJET->Phi();
	    float dr = aJET->DeltaR(trk);
	    if(dr<=0.5) {
	      float dca = GTRACK333_dca0[j];
	      double sign = dca*TMath::Sin(deltaPhi);
	      if(GTRACK333_dcaerr0[j]>0 && TMath::Abs(dca)<=0.1) {
		float dcaSig = TMath::Abs(dca)/GTRACK333_dcaerr0[j];
		if(sign>=0) h_dcaSig_p->Fill(dcaSig);
		else h_dcaSig_n->Fill(dcaSig);
	      }
	    }
	  }
	}
      }

      
    } // jpsi
    jets.Clear();
  } //entries

  /*
  h_njpsijets->Fill(0,njpsi0);
  h_njpsijets->Fill(1,njpsi1);
  h_njpsijets->Fill(2,njpsi2);

  h_njpsijets_matched->Fill(1,njpsiORajet);
  h_njpsijets_matched->Fill(2,njpsiANDajet);
  */

  cout << endl;
  cout << "------------------------------------" << endl;
  cout << "N jpsi = " << njpsi0 << endl;
  cout << "N jpsi-jets = " << njpsijet << endl;
  cout << "N away-jets = " << najet << endl;
  cout << "N jpsi + ajet = " << njpsiANDajet << endl;
  cout << "------------------------------------" << endl;
  cout << endl;

  
  TFile file(fileName,"recreate");

  //JPSI
  h_njets->Write();
  h_jetaphi->Write();
  h_jpsiJet->Write(); h_aJet->Write();
  h_njpsijets->Write(); h_njpsijets_matched->Write();


  //SVX
  hvtx_dphi->Write();
  hvtx_dl->Write();
  hvtx_mass->Write();
  hvtx_mult->Write();

  hvtx_dl_away->Write();



  // IP
  h_dca_jpsi->Write(); h_sdca_jpsi->Write();

  h_dcaSig_n->Write();
  h_dcaSig_p->Write();


  file.Close();
}







///---------------------------------------------------
/// b-tagging studies
///
///---------------------------------------------------
void V0Analysis::AnalyzeBTagging_muon(char* fileName, bool debug)
{
  CHI2 = 5;


  // JPSI
  TH1F* h_njets = new TH1F("h_njets","Jet multiplicity",5,-0.5,4.5);
  TH2F* h_jetaphi = new TH2F("h_jetaphi","Cal-Jet #eta vs. #phi",
			     100,-3,3,100,-3.15,3.15);
  TH1F* h_jpsiJet = new TH1F("h_jpsiJet","#DeltaR J/#psi - Jet",100,0,1);
  TH1F* h_aJet = new TH1F("h_aJet","#Delta#phi J/#psi - away Jet",100,0,1);
  TH1F* h_njpsijets = new TH1F("h_njpsijets","Inclusive J/#psi + Jet multiplicity",3,-0.5,2.5);
  TH1F* h_njpsijets_matched = new TH1F("h_njpsijets_matched","Inclusive J/#psi + Jet multiplicity",3,-0.5,2.5);




  // SECONDARY VERTEX away-jet
  TH1F* h_vtx_dphi = new TH1F("h_vtx_dphi","#Delta#phi(SVX,#mu)",100,0,3.15);
  TH1F* h_vtx_dphi_i = new TH1F("h_vtx_dphi_i","#Delta#phi(SVX,#mu)",100,0,3.15);
  TH1F* h_vtx_dl = new TH1F("h_vtx_dl","Secondary Vertex decay length",
			    100,-0.5,0.5);
  h_vtx_dl->SetXTitle("cm");
  TH1F* h_vtx_dl_i = new TH1F("h_vtx_dl_i","Secondary Vertex decay length",
			      100,-0.5,0.5);
  h_vtx_dl_i->SetXTitle("cm");
  TH1F* h_vtx_dl_ii = new TH1F("h_vtx_dl_ii","Secondary Vertex decay length [away-jet]",
			       100,-0.5,0.5);
  h_vtx_dl_ii->SetXTitle("cm");
  TH1F* h_vtx_dl_iii = new TH1F("h_vtx_dl_iii","Secondary Vertex decay length [#mu jet]",
			       100,-0.5,0.5);
  h_vtx_dl_iii->SetXTitle("cm");
  
  TH1F* h_vtx_mass = new TH1F("h_vtx_mass","Secondary Vertex mass",100,0,1);
  h_vtx_mass->SetXTitle("GeV/c^{2}");
  TH1F* h_vtx_mult = new TH1F("h_vtx_mult","Secondary Vertex multiplicity",
			      5,1.5,6.5);

  TH1F* hvtx_dl_away = new TH1F("hvtx_dl_away","Secondary Vertex decay length [away jet]",100,-1,1);
  hvtx_dl_away->SetXTitle("cm");


  // IMPACT-PARAMETER jpsi
  TH1F* h_dca_jpsi = new TH1F("h_dca_jpsi","J/#psi Signed Impact Parameter",
			      100,-0.2,0.2);
  h_dca_jpsi->SetXTitle("cm");
  TH1F* h_sdca_jpsi = new TH1F("h_sdca_jpsi","J/#psi Signed Impact Parameter Significance",100,-30,30);


  // IMPACT-PARAMETER away-jet
  TH1F* h_dcaSig_p = new TH1F("h_dcaSig_p","Signed Impact Parameter Significance [p_{T} > 1.5 GeV/c]",200,0,30);
  TH1F* h_dcaSig_n = new TH1F("h_dcaSig_n","Signed Impact Parameter Significance [p_{T} > 1.5 GeV/c]",200,0,30);




  float px1,px2,py1,py2,pz1,pz2,p1,p2,q1,q2,dedx1,dedx2,dca1,dca2;  
  float pt1, pt2;
  int nsmt1,ncft1,nsmt2,ncft2;
  int id1, id2;

  int njpsi0 = 0; int njpsi1 = 0; int njpsi2 = 0;
  int njpsijet = 0;
  int najet = 0;
  int njpsiANDajet = 0;
  int njpsiORajet = 0;


  bool ks;

  int aa = 1000; int nn = 0;
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());  
  Int_t nbytes = 0, nb = 0;
  for (Int_t jentry=0; jentry<nentries;jentry++) {
    Int_t ientry = LoadTree(jentry); 
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    if(jentry/aa>nn) cout << "Event #" << jentry << endl;
    nn = jentry/aa;


    // 1-Primary vertex
    bool pv = false;
    if(PVSEL333_psnvtx333>0 && PVSEL333_psntrk[0]>=3) pv = true;
    if(!pv) continue;
    float PV_x = PVSEL333_psx[0];
    float PV_y = PVSEL333_psy[0];
    float PV_z = PVSEL333_psz[0];

    // 2-Jets
    TObjArray jets; readJets(jets);
    h_njets->Fill(jets.GetEntries());
    for(int y=0; y<jets.GetEntries(); y++) {
      TLorentzVector* jet = (TLorentzVector*) jets[y];
      h_jetaphi->Fill(jet->Phi(),jet->Eta());
     }
    


    // 3- Muons
    int muon = 0;
    if(Muonid_Nmuonid<1) continue;
    for(int u=0; u<Muonid_Nmuonid && muon==0; u++) {
      if(Muonid_Idpt[u]>3 && TMath::Abs(Muonid_Idnseg[u])==3)
	muon = u;
    }
    if(muon==0) continue;
    TLorentzVector mu1(Muonid_Idpx[muon],Muonid_Idpy[muon],Muonid_Idpz[muon],
		       Muonid_IdE[muon]);
    

    //3- Secondary vertices
    int jpsi = -1; float ptmin = 0; TLorentzVector jpsiVTX;
    for(int i=0; i<SKSEL333_kssnvtx333; i++) {
      float dl    = decayLength_sksel(PV_x, PV_y, 0.0, i);
      if(dl>0.5) continue;
      
      TLorentzVector vtx; vtx.SetXYZM(SKSEL333_ksspx[i],SKSEL333_ksspy[i],
				      SKSEL333_ksspz[i],SKSEL333_kssmass[i]);
      float deltaPhi = TMath::Abs(mu1.DeltaPhi(vtx));
      if(deltaPhi>TMath::Pi())  2*TMath::Pi()-deltaPhi ;
      h_vtx_dphi->Fill(deltaPhi);

      h_vtx_mass->Fill(vtx.M());
      ks = false; if(vtx.M()>=0.486 && vtx.M()<0.50) ks = true;

      h_vtx_dl->Fill(dl);

      if(!ks) {
	h_vtx_dphi_i->Fill(deltaPhi);
	h_vtx_dl_i->Fill(dl);
	
	if(deltaPhi>3) {
	  h_vtx_dl_ii->Fill(dl);
	}

	if(deltaPhi<0.3) {
	  h_vtx_dl_iii->Fill(dl);
	}



      }



    }
     

 
    /*
      //define Jpsi-Jet 
      TLorentzVector* jpsiJET = 0; float min = 0.5;
      for(int y=0; y<jets.GetEntries(); y++) {
	TLorentzVector* jet = (TLorentzVector*) jets[y];
	float dr = jet->DeltaR(jpsiVTX);
	if(dr<min) {
	  min = dr;
	  jpsiJET = jet;
	}
      }
      if(jpsiJET) {
	h_jpsiJet->Fill(min);
	njpsijet++;
      }
      
      
      //define away-jet
      TLorentzVector* aJET = 0; min = 0.5;
      for(int y=0; y<jets.GetEntries(); y++) {
	TLorentzVector* jet = (TLorentzVector*) jets[y];
	float dphi = TMath::Abs(jet->DeltaPhi(jpsiVTX));
	if(TMath::Pi()-dphi<min) {
	  min = TMath::Pi()-dphi;
	  aJET = jet;
	}
      }
      if(aJET) {
	h_aJet->Fill(min);
	najet++;
      }

      if(jpsiJET || aJET) njpsiORajet++;
      if(jpsiJET && aJET) njpsiANDajet++;


      if(debug) {
	cout << "Jpsi    :" << jpsiVTX.Eta() << " " << jpsiVTX.Phi();
	cout << "  mass = " << jpsiVTX.M() << " " << V0_V0_mmu[jpsi] <<endl;
	if(jpsiJET) 
	  cout << "jpsi-jet:" << jpsiJET->Eta() << " " << jpsiJET->Phi()<<endl;
	if(aJET) 
	  cout << "a-jet   :" << aJET->Eta() << " " << aJET->Phi() << endl;
	cout << endl;
      }
    
    */

    /*
      ///------------------------------------
      /// Secondary vertex analysis
      ///------------------------------------
      for(int u=0; u<SKSEL333_kssnvtx333; u++) {
	TLorentzVector vtx; 
	vtx.SetXYZM(SKSEL333_ksspx[u],SKSEL333_ksspy[u],
		    SKSEL333_ksspz[u],SKSEL333_kssmass[u]);
	float dl_svx  = decayLength(PV_x, PV_y, 0.0, u);
	float dphi = TMath::Abs(jpsiVTX.DeltaPhi(vtx));
	hvtx_dphi->Fill(dphi);
	hvtx_dl->Fill(dl_svx);
	hvtx_mass->Fill(vtx.M());
 	hvtx_mult->Fill(SKSEL333_kssntrk[u]);
	
	//look at vertices in away jet for different j/psi pTs
	if(dphi>=TMath::Pi()-0.5) {
	  hvtx_dl_away->Fill(dl_svx);
	}
      }
    */  
      
    /*
      ///------------------------------------
      /// Impact-Parameter analysis
      ///------------------------------------
      if(jpsiJET) {
 	float dca,sdca;
	get_kinematics(jpsi,q1,px1,py1,pz1,p1,dedx1,nsmt1,ncft1,dca1,id1,
		       q2,px2,py2,pz2,p2,dedx2,nsmt2,ncft2,dca2,id2);
	//track 1
	signedIP(id1,jpsiJET->Phi(),dca,sdca);
	h_dca_jpsi->Fill(dca); h_sdca_jpsi->Fill(sdca);

	//track 2
	signedIP(id2,jpsiJET->Phi(),dca,sdca);
	h_dca_jpsi->Fill(dca); h_sdca_jpsi->Fill(sdca);
      }
    */
    /*
      if(aJET) {
	//match tracks to away-jet
	for(int j=0; j<GTRACK333_ngtrk333; j++) {
	  float pT = 1./GTRACK333_qoverpt[j];
	  int hits = GTRACK333_nsmt[j]+GTRACK333_ncft[j];
	  if(pT>1.5 && hits>=11) {
	    float phi = GTRACK333_phidir[j]; if(phi<0) phi += 2*TMath::Pi();
	    double px = pT*TMath::Cos(phi);
	    double py = pT*TMath::Sin(phi);
	    double pz = pT*GTRACK333_tanl[j];
	    double E = TMath::Sqrt(px*px+py*py+pz*pz+0.13956995*0.13956995);  
	    TLorentzVector trk(px,py,pz,E);
	    float deltaPhi = trk.Phi() - aJET->Phi();
	    float dr = aJET->DeltaR(trk);
	    if(dr<=0.5) {
	      float dca = GTRACK333_dca0[j];
	      double sign = dca*TMath::Sin(deltaPhi);
	      if(GTRACK333_dcaerr0[j]>0 && TMath::Abs(dca)<=0.1) {
		float dcaSig = TMath::Abs(dca)/GTRACK333_dcaerr0[j];
		if(sign>=0) h_dcaSig_p->Fill(dcaSig);
		else h_dcaSig_n->Fill(dcaSig);
	      }
	    }
	  }
	}
      }
    */
    jets.Clear();
  } //entries

  /*
  h_njpsijets->Fill(0,njpsi0);
  h_njpsijets->Fill(1,njpsi1);
  h_njpsijets->Fill(2,njpsi2);

  h_njpsijets_matched->Fill(1,njpsiORajet);
  h_njpsijets_matched->Fill(2,njpsiANDajet);
  */

  cout << endl;
  cout << "------------------------------------" << endl;
  cout << "N jpsi = " << njpsi0 << endl;
  cout << "N jpsi-jets = " << njpsijet << endl;
  cout << "N away-jets = " << najet << endl;
  cout << "N jpsi + ajet = " << njpsiANDajet << endl;
  cout << "------------------------------------" << endl;
  cout << endl;

  
  TFile file(fileName,"recreate");

  //JPSI
  h_njets->Write();
  h_jetaphi->Write();
  h_jpsiJet->Write(); h_aJet->Write();
  h_njpsijets->Write(); h_njpsijets_matched->Write();


  //SVX
  h_vtx_dphi->Write(); h_vtx_dphi_i->Write();
  h_vtx_dl->Write(); h_vtx_dl_i->Write(); h_vtx_dl_ii->Write(); 
  h_vtx_dl_iii->Write();

  h_vtx_mass->Write();
  h_vtx_mult->Write();

  hvtx_dl_away->Write();



  // IP
  h_dca_jpsi->Write(); h_sdca_jpsi->Write();

  h_dcaSig_n->Write();
  h_dcaSig_p->Write();


  file.Close();
}