#include "StHbtMaker/Infrastructure/StHbtTrack.hh"
#include "StHbtMaker/Cut/fxtTrackCutMonitor.h"
#include <cstdio>
#include <string>
#include "StLorentzVector.hh"

#ifdef __ROOT__ 
ClassImp(fxtTrackCutMonitor)
#endif

//_________________
fxtTrackCutMonitor::fxtTrackCutMonitor(const char* name,
				       const double ParticleMass){ // default constructor
  monMass = ParticleMass;
  string s("fxtTrackCutMonitor");
  string n(name);
  mDCAGlobal= new StHbt1FHisto( (s+n+"mDCAGlobal").c_str(),
				"DCA Global; DCA Global (cm)",
				100,0.,5.); 
  mNhits = new StHbt1FHisto( (s+n+"mNhits").c_str(),"nHits;nHits",
			     50, 0., 50.);
  mP = new StHbt1FHisto( (s+n+"mP").c_str(),"Momentum;P (GeV/c)",
			  180,0.,1.8);
  mPt = new StHbt1FHisto( (s+n+"mPt").c_str(),
			  "Transverse momentum; Pt (GeV/c)",
			  180,0.,1.8);
  mPtVsNsigmaPion = new StHbt2FHisto( (s+n+"mPtVsNsigmaPion").c_str(),
				     "p{T} vs n#sigma(#pi);p_{T} (GeV/c);n#sigma(#pi)",
				      36,0.1,1.9,100,-10.,10.);
  mPtVsNsigmaKaon = new StHbt2FHisto( (s+n+"mPtVsNsigmaKaon").c_str(),
				     "p_{T} vs n#sigma(K);p_{T} (GeV/c);n#sigma(K)",
				     36,0.1,1.9,100,-10.,10.);
  mPtVsNsigmaProton = new StHbt2FHisto( (s+n+"mPtVsNsigmaProton").c_str(),
				    "p_{T} vs n#sigma(p);p_{T} (GeV/c);n#sigma(p)",
					36,0.1,1.9,100,-10.,10.);
  mPvsDedx= new StHbt2FHisto( (s+n+"mPvsDedx").c_str(),
			      "p vs dE/dx (a.u.);p (GeV/c); dE/dx (a.u.)",
			      180, 0.1, 1.9, 100, 0., 10.); 
  mRapidity = new StHbt1FHisto( (s+n+"mRapidity").c_str(),
				"Rapidity;y",
				70, -1.2, 1.2);
  mPseudoRapidity = new StHbt1FHisto( (s+n+"mPseudoRapidity").c_str(),
				      "PseudoRapidity;#eta",
				      120, -1.2, 1.2);
  mPvsMassSqr = new StHbt2FHisto( (s+n+"mPvsMassSqr").c_str(),
				  "p_{T} vs m^{2};p_{T} (GeV/c);m^{2} (GeV/c^{2})^{2}",
				   180,0.1,1.9,120, -0.1, 1.1);
  mPvsInvBeta = new StHbt2FHisto( (s+n+"mPvsInvBeta").c_str(),
				  "1/#beta vs P;p (GeV/c);1/#beta",
				  180, 0.1, 1.9, 60, 0.8, 2.);
#ifdef TPC_DNDX
  mPtVsDndxNsigmaPion = new StHbt2FHisto( (s+n+"mPtVsDndxNsigmaPion").c_str(),
                                          "p_{T} vs n#sigma^{dN/dx}(#pi);p_{T} (GeV/c);n#sigma^{dN/dx}(#pi)",
                                          36,0.1,1.9, 100, -10., 10.);
  mPtVsDndxNsigmaKaon = new StHbt2FHisto( (s+n+"mPtVsDndxNsigmaKaon").c_str(),
                                          "p_{T} vs n#sigma^{dN/dx}(K);p_{T} (GeV/c);n#sigma^{dN/dx}(K)",
                                          36,0.1,1.9, 100, -10., 10.);
  mPtVsDndxNsigmaProton = new StHbt2FHisto( (s+n+"mPtVsDndxNsigmaProton").c_str(),
                                          "p_{T} vs n#sigma^{dN/dx}(p);p_{T} (GeV/c);n#sigma^{dN/dx}(p)",
                                          36,0.1,1.9, 100, -10., 10.);
  mPvsDndx= new StHbt2FHisto( (s+n+"mPvsDndx").c_str(),
	                  		      "p vs dN/dx (a.u.);p (GeV/c); dN/dx (a.u.)",
			                        180, 0.1, 1.9, 400, 20., 100.);
#endif
}

//_________________
fxtTrackCutMonitor::fxtTrackCutMonitor( const  fxtTrackCutMonitor& cutMoni)  {
  mDCAGlobal =new StHbt1FHisto(*(cutMoni.mDCAGlobal));
  mNhits = new StHbt1FHisto(*(cutMoni.mNhits));
  mP = new StHbt1FHisto(*(cutMoni.mP));
  mPt = new StHbt1FHisto(*(cutMoni.mPt));
  mPtVsNsigmaPion = new StHbt2FHisto(*(cutMoni.mPtVsNsigmaPion));
  mPtVsNsigmaKaon = new StHbt2FHisto(*(cutMoni.mPtVsNsigmaKaon));
  mPtVsNsigmaProton = new StHbt2FHisto(*(cutMoni.mPtVsNsigmaProton));
  mPvsDedx =new StHbt2FHisto(*(cutMoni.mPvsDedx));
  mRapidity= new StHbt1FHisto(*(cutMoni.mRapidity));
  mPseudoRapidity= new StHbt1FHisto(*(cutMoni.mPseudoRapidity));
  mPvsMassSqr = new StHbt2FHisto(*(cutMoni.mPvsMassSqr));
  mPvsInvBeta = new StHbt2FHisto(*(cutMoni.mPvsInvBeta));
#ifdef TPC_DNDX
  mPtVsDndxNsigmaPion = new StHbt2FHisto(*(cutMoni.mPtVsDndxNsigmaPion));
  mPtVsDndxNsigmaKaon = new StHbt2FHisto(*(cutMoni.mPtVsDndxNsigmaKaon));
  mPtVsDndxNsigmaProton = new StHbt2FHisto(*(cutMoni.mPtVsDndxNsigmaProton));
  mPvsDndx =new StHbt2FHisto(*(cutMoni.mPvsDndx));
#endif
}

//_________________
fxtTrackCutMonitor::~fxtTrackCutMonitor(){
  delete mDCAGlobal;
  delete mNhits;
  delete mP;
  delete mPt;
  delete mPtVsNsigmaPion;
  delete mPtVsNsigmaKaon;
  delete mPtVsNsigmaProton;
  delete mPvsDedx ;
  delete mRapidity;
  delete mPseudoRapidity;
  delete mPvsMassSqr;
  delete mPvsInvBeta;
#ifdef TPC_DNDX
  delete mPtVsDndxNsigmaPion;
  delete mPtVsDndxNsigmaKaon;
  delete mPtVsDndxNsigmaProton;
  delete mPvsDndx;
#endif
}

//_________________
void  fxtTrackCutMonitor::Fill(const StHbtTrack* track){

  float TEnergy = ::sqrt(track->P().mag2()+monMass*monMass);
  float TRapidity = 0.5*::log((TEnergy+track->P().z())/
			      (TEnergy-track->P().z()));
  /*
    float TPhi = atan2(track->P().y(),track->P().x());
    if(TPhi<0.0) { TPhi += 2.*TMath::Pi(); }
    TPhi *= 180./TMath::Pi();
  */
  
  mDCAGlobal->Fill( track->DCAxyGlobal(), 1.); //using for DCAxyz global
  mNhits->Fill(track->NHits(), 1.);
  mP->Fill(track->P().mag(), 1.);
  mPt->Fill(track->Pt(), 1.);
  mPtVsNsigmaPion->Fill(track->Pt(), track->NSigmaPion(), 1.);
  mPtVsNsigmaKaon->Fill(track->Pt(), track->NSigmaKaon(), 1.);
  mPtVsNsigmaProton->Fill(track->Pt(), track->NSigmaProton(), 1.);
  mPvsDedx->Fill( fabs(track->P().mag()), track->dEdx()*1000000, 1.);
  mRapidity->Fill(TRapidity, 1.);
  mPseudoRapidity->Fill(track->Eta(), 1.);
  mPvsMassSqr->Fill(track->P().mag(), track->TofMassSqr(), 1.);
  mPvsInvBeta->Fill( fabs(track->P().mag()), (1./track->TofBeta()), 1.);
#ifdef TPC_DNDX
  mPvsDndx->Fill(track->P().mag(), track->dNdx()*1000000, 1.);
  mPtVsDndxNsigmaPion->Fill(track->Pt(), track->DndxNSigmaPion(), 1.);
  mPtVsDndxNsigmaKaon->Fill(track->Pt(), track->DndxNSigmaKaon(), 1.);
  mPtVsDndxNsigmaProton->Fill(track->Pt(), track->DndxNSigmaProton(), 1.);
#endif
}