// -------------------------------------------------------------------------
// -----                  CbmSttFindTracksQa source file               -----
// -----                  Created 30/03/06  by V. Friese               -----
// -------------------------------------------------------------------------


#include "iostream.h"

#include "TObjArray.h"
#include "TVector3.h"

#include "CbmMCTrack.h"
#include "CbmRootManager.h"
#include "CbmRuntimeDb.h"
#include "CbmSttFindTracksQa.h"
#include "CbmSttTrackMatch.h"



// -----   Default constructor   -------------------------------------------
CbmSttFindTracksQa::CbmSttFindTracksQa() 
    : CbmTask("STT Track Finder QA", 1) 
{ 
  fTotalAcc = 0.;
  fTotalRec = 0.;
  fTotalGhost = 0.;
  fTotalClone = 0.;

  fMinPoints = 2;
  fQuota     = 0.7;
  CreateHistograms();
}
// -------------------------------------------------------------------------



// -----   Standard constructor   ------------------------------------------
CbmSttFindTracksQa::CbmSttFindTracksQa(Int_t minPoints, Double_t quota,
				       Int_t iVerbose) 
  : CbmTask("STTTrack Finder QA", iVerbose) 
{
  fTotalAcc = 0.;
  fTotalRec = 0.;
  fTotalGhost = 0.;
  fTotalClone = 0.;
  
  fMinPoints = minPoints;
  fQuota     = quota;
  CreateHistograms();
}
// -------------------------------------------------------------------------

  
  
// -----   Destructor   ----------------------------------------------------
CbmSttFindTracksQa::~CbmSttFindTracksQa() { }
// -------------------------------------------------------------------------



// -----   Public method SetParContainers   --------------------------------
void CbmSttFindTracksQa::SetParContainers() 
{
}
// -------------------------------------------------------------------------



// -----   Public method Init   --------------------------------------------
InitStatus CbmSttFindTracksQa::Init() {

  cout << "==========================================================="
       << endl;;
  cout << GetName() << ": Initialising..." << endl;

  // Get RootManager
  CbmRootManager* ioman = CbmRootManager::Instance();
  if (! ioman) {
    cout << "-E- " << GetName() << "::Init: "
	 << "RootManager not instantised!" << endl;
    return kFATAL;
  }

  // Get MCTrack array
  fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
  if ( ! fMCTracks ) {
    cout << "-E- " << GetName() << "::Init: No MCTrack array!" << endl;
    return kFATAL;
  }
   
  // Get SttTrack array
  fSttTracks = (TClonesArray*) ioman->GetObject("STTTrack");
  if ( ! fSttTracks ) {
    cout << "-E- " << GetName() << "::Init: No SttTrack array!" << endl;
    return kERROR;
  }

  // Get SttTrackMatch array
  fMatches = (TClonesArray*) ioman->GetObject("STTTrackMatch");
  if ( ! fMatches ) {
    cout << "-E- " << GetName() << "::Init: No SttTrackMatch array!" 
	 << endl;
    return kERROR;
  }

  fTargetPos.SetXYZ(0., 0., 0.);

  // Output
  cout << "   Minimum number of STT points : " << fMinPoints << endl;
  cout << "   Matching quota               : " << fQuota << endl;
  cout << "   Target position ( " << fTargetPos.X() << ", " 
       << fTargetPos.Y() << ", " << fTargetPos.Z() << ") " << endl;
  if (fActive) cout << "   *****   Task is ACTIVE   *****" << endl;
  cout << "==========================================================="
       << endl << endl;

  return kSUCCESS;

}
// -------------------------------------------------------------------------



// -----   Public method ReInit   ------------------------------------------
InitStatus CbmSttFindTracksQa::ReInit() {

  cout << "==========================================================="
       << endl;;
  cout << GetName() << ": Reinitialising..." << endl;

  fTargetPos.SetXYZ(0., 0., 0.);

  // Output
  cout << "   Target position ( " << fTargetPos.X() << ", " 
       << fTargetPos.Y() << ", " << fTargetPos.Z() << ") " << endl;
  if (fActive) cout << "   *****   Task is ACTIVE   *****" << endl;
  cout << "==========================================================="
       << endl << endl;

  return kSUCCESS;

}
// -------------------------------------------------------------------------


// -----   Public method Exec   --------------------------------------------
void CbmSttFindTracksQa::Exec(Option_t* opt) 
{
  // Clear matching map and reset eventwise counters
  fMatchMap.clear();
  Int_t nAll     = 0;
  Int_t nAllAcc  = 0;
  Int_t nAllRec  = 0;
  Int_t nPrim    = 0;
  Int_t nPrimAcc = 0;
  Int_t nPrimRec = 0;
  Int_t nSec     = 0;
  Int_t nSecAcc  = 0;
  Int_t nSecRec  = 0;
  Int_t nRef     = 0;
  Int_t nRefRec  = 0;
  Int_t nGhost   = 0;
  Int_t nClone  = 0;

  cout << "1" << endl;

  // Fill map from MCTrack index to matched SttTrack Index
  Int_t nRec = fSttTracks->GetEntriesFast();
  Int_t nMtc = fMatches->GetEntriesFast();
  if ( nMtc != nRec ) 
    {
      cout << "-E- " << GetName() << ":Exec Number of SttMatches ("
	   << nMtc << ") does not equal number of SttTracks ("
	   << nRec << ")" << endl;
      cout << "                               Skipping event " << endl;
      return;
    }

  cout << "2" << endl;

  for (Int_t iRec=0; iRec<nRec; iRec++) 
    {
      CbmSttTrackMatch
	  *match = (CbmSttTrackMatch*) fMatches->At(iRec);

      CbmSttTrack
	  *track = (CbmSttTrack*) fSttTracks->At(iRec);
      
      if (! match) 
      {
	  cout << "-W- " << GetName() << "::Exec: Empty SttTrackMatch at "
	       << iRec << endl;
	  continue;
      }

      Int_t iMC = match->GetMCTrackID();

      if (iMC == -1 ) 
      {
	  cout << "-I- " << GetName() << ": SttTrack " << iRec 
	       << " has no corresponding MCTrack " << endl;

	  FillGhostHistos(track);
	  
	  nGhost++;
	  continue;
      }
      
      if ( fMatchMap.find(iMC) == fMatchMap.end() ) 
	  fMatchMap[iMC] = iRec;
      else 
      {
	  if ( fVerbose>1 ) 
	      cout << "-W- " << GetName() << "::Exec: "
		   << "MCTrack " << iMC << " doubly matched."
		   << "Current match " << iRec 
		   << ", previous match " << fMatchMap[iMC] 
		   << endl;

	  FillCloneHistos(track);
	  nClone++;
      }
    }  // Loop over SttTrackMatches
  

  cout << "3" << endl;

  // Some variables outside loop
  TVector3 vertex;
  CbmMCTrack*       mcTrack  = NULL;
  CbmSttTrack*      sttTrack = NULL;
  CbmSttTrackMatch* match    = NULL;
  Bool_t isAcc  = kFALSE;
  Bool_t isPrim = kFALSE;
  Bool_t isRef  = kFALSE;
  Bool_t isRec  = kFALSE;
  Double_t mom  =  0.;
  Int_t iRec    = -1;
  Int_t nHits   =  0;
  Int_t nTrue   =  0;
  Int_t nPoints =  0;

  cout << "4" << endl;

  // Loop over MCTracks
  Int_t nMC = fMCTracks->GetEntriesFast();
  for (Int_t iMC=0; iMC<nMC; iMC++) 
    {
      mcTrack = (CbmMCTrack*) fMCTracks->At(iMC);
      if ( ! mcTrack ) 
      {
	  cout << "-W- " << GetName() << "::Exec: Empty MCTrack at " << iMC
	       << endl;
	  continue;
      }
      
      // Check geometrical acceptance
      isAcc = kFALSE;

      nPoints = mcTrack->GetSttPoints();

      if ( nPoints>= fMinPoints ) 
	{
	  isAcc = kTRUE;
	}

      FillNormHistos(mcTrack, isAcc);

      // Check origin of MCTrack
      vertex = mcTrack->GetStartVertex();
      isPrim = kFALSE;
      if ( (vertex-fTargetPos).Mag() < 1. ) 
	isPrim = kTRUE;
      
      // Counters
      nAll++;
      if ( isPrim ) 
	nPrim++;
      else          
	nSec++;

      if ( fVerbose > 2 ) 
	cout << "-I- " << GetName() << ": MCTrack " << iMC
	     << ", points " << nPoints << ", not accepted"
	     << endl;

      if ( ! isAcc) 
	continue;
      
      // Only accepted MCTracks from now on
      nAllAcc++;
      if ( isPrim ) 
	nPrimAcc++;
      else          
	nSecAcc++;

      // Check momentum
      mom = mcTrack->GetMomentum().Mag();
      isRef = kFALSE;
      if ( mom > 1. && isPrim) 
	{
	  isRef = kTRUE;
	  nRef++;
	}
      
      // Get attached SttTrack
      if (fMatchMap.find(iMC) == fMatchMap.end() ) 
	{
	  isRec = kFALSE;
	  iRec = -1;
	  nHits = nTrue = 0;
	}
      else 
	{
	  iRec = fMatchMap[iMC];
	  sttTrack = (CbmSttTrack*) fSttTracks->At(iRec);
	  if ( ! sttTrack ) 
	    {
	      cout << "-E- " << GetName() << "::Exec: No SttTrack for MCTrack "
		   << iMC << endl;
	      continue;
	    }
	  match = (CbmSttTrackMatch*) fMatches->At(iRec);
	  if ( ! match ) 
	    {
	      cout << "-E- " << GetName() << "::Exec: No SttTrackMatch for MCTrack "
		   << iMC << endl;
	      continue;
	    }

	  // Check for reconstruction criterion
	  nHits = sttTrack->GetNHits();
	  nTrue = match->GetNofTrueHits();
	  isRec = kFALSE;

	  if ( (Double_t(nTrue) / Double_t(nHits) ) >= fQuota ) 
	  {
	      isRec = kTRUE;
	      FillEffiencyHistos(mcTrack, isAcc);
	  }
	  else 
	  {
	      FillGhostHistos(sttTrack);
	      nGhost++;
	  }
	}

      // Verbose output
      if ( fVerbose > 1 ) 
	{
	  cout << "-I- " << GetName() << ": MCTrack " << iMC << ", points "
	       << nPoints << ", SttTrack " << iRec << ", hits " << nHits 
	       << ", true hits " << nTrue;
	  if ( isRec ) 
	    cout << ", reconstructed" << endl;
	  else   
	    cout << ", not reconstructed" << endl;
	}
      if ( isRec ) 
	{
	  nAllRec++;
	  if ( isPrim ) 
	    nPrimRec++;
	  else          
	    nSecRec++;
	  if (isRef)    
	    nRefRec++;
	}
    }

  // Calculate efficiencies
  Double_t 
    allEff = 0.,
    primEff = 0.,
    refEff = 0.;

  if (nAllAcc > 0.)
    allEff  = Double_t(nAllRec)  / Double_t(nAllAcc);

  if (nPrimAcc > 0.)
    primEff = Double_t(nPrimRec) / Double_t(nPrimAcc);

  if (nRef > 0.)
    refEff  = Double_t(nRefRec)  / Double_t(nRef);

  cout << "6" << endl;
  // Event summary
  if ( fVerbose > 0 ) 
    {
      cout << "-------   SttFindTracksQa : Event summary   ---------" << endl;
      cout << "MCTracks   : " << nAll << ", accepted: " << nAllAcc
	   << ", reconstructed: " << nAllRec << endl;
      cout << "Primaries  : " << nPrim << ", accepted: " << nPrimAcc
	   << ", reconstructed: " << nPrimRec << endl;
      cout << "Secondaries: " << nSec << ", accepted: " << nSecAcc
	   << ", reconstructed: " << nSecRec << endl;
      cout << "Efficiencies: all " << allEff*100. << " %, primaries "
	   << primEff*100. << "%, reference " << refEff*100. << "%" << endl;
      cout << "STTTracks " << nRec << ", ghosts " << nGhost 
	   << ", clones " << nClone << endl;
      cout << "-----------------------------------------------------" << endl;
      cout << endl;
    }
  else cout << "-I- " << GetName() << ": all " << allEff*100. << " %, prim. "
	    << primEff*100. << " %, ref. " << refEff*100. << " %" << endl;
  
  if (nAllRec > 0)
    {
	fTotalAcc += nAllAcc;
	fTotalRec += nAllRec;
	fTotalGhost += nGhost;
	fTotalClone += nClone;
	
	GoodTracks->Fill(nAllAcc, nAllRec);
	GhostTracks->Fill(nAllAcc, nGhost);
	CloneTracks->Fill(nAllAcc, nClone);
    }

  cout << "7" << endl;
  cout << "total eff: " << fTotalAcc   << " tracks accepted" << endl;
  cout << "           " << fTotalRec   << " tracks reconstructed (" << (fTotalRec / fTotalAcc) * 100. << "%)" << endl;
  cout << "           " << fTotalGhost << " ghost tracks (" << (fTotalGhost / fTotalAcc) * 100. << "%)" << endl;
  cout << "           " << fTotalClone - fTotalGhost << " clone tracks (" << ((fTotalClone - fTotalGhost) / fTotalAcc) * 100. << "%)" << endl;
}

// -------------------------------------------------------------------------
void CbmSttFindTracksQa::FillGhostHistos(CbmSttTrack *recTrack)
{
    Double_t
	recTrackMomentumTx = recTrack->GetParamLast()->GetTx(),
	recTrackMomentumTy = recTrack->GetParamLast()->GetTy(),
	recTrackTheta = acos(1. / sqrt(recTrackMomentumTx * recTrackMomentumTx + 
				       recTrackMomentumTy * recTrackMomentumTy + 1.)),
	recTrackPhi = atan(recTrackMomentumTy / recTrackMomentumTx);
    
    fPhiGhosts->Fill(recTrackPhi);
    fThetaGhosts->Fill(recTrackTheta);
}

void CbmSttFindTracksQa::FillCloneHistos(CbmSttTrack *recTrack)
{
    Double_t
	recTrackMomentumTx = recTrack->GetParamLast()->GetTx(),
	recTrackMomentumTy = recTrack->GetParamLast()->GetTy(),
	recTrackTheta = acos(1. / sqrt(recTrackMomentumTx * recTrackMomentumTx + 
				       recTrackMomentumTy * recTrackMomentumTy + 1.)),
	recTrackPhi = atan(recTrackMomentumTy / recTrackMomentumTx);
    
    fPhiClones->Fill(recTrackPhi);
    fThetaClones->Fill(recTrackTheta);
}

void CbmSttFindTracksQa::FillEffiencyHistos(CbmMCTrack *mcTrack, Bool_t isAcc)
{
/*
    Double_t
	mcTrackMomentumTx = mcTrack->GetMomentum().X() / mcTrack->GetMomentum().Z(),
	mcTrackMomentumTy = mcTrack->GetMomentum().Y() / mcTrack->GetMomentum().Z(),
	mcTrackTheta = acos(1. / sqrt(mcTrackMomentumTx * mcTrackMomentumTx + 
				      mcTrackMomentumTy * mcTrackMomentumTy + 1.)),
	mcTrackPhi = atan(mcTrackMomentumTy / mcTrackMomentumTx);

    fPhiEfficiency->Fill(mcTrackPhi);
    fThetaEfficiency->Fill(mcTrackTheta);
    
    if (isAcc)
    {
	fPhiEfficiencyAccepted->Fill(mcTrackPhi);
	fThetaEfficiencyAccepted->Fill(mcTrackTheta);
    }
*/
}

void CbmSttFindTracksQa::FillNormHistos(CbmMCTrack *mcTrack, Bool_t isAcc)
{
    /*
    Double_t
	mcTrackMomentumTx = mcTrack->GetMomentum().X() / mcTrack->GetMomentum().Z(),
	mcTrackMomentumTy = mcTrack->GetMomentum().Y() / mcTrack->GetMomentum().Z(),
	mcTrackTheta = acos(1. / sqrt(mcTrackMomentumTx * mcTrackMomentumTx + 
				      mcTrackMomentumTy * mcTrackMomentumTy + 1.)),
	mcTrackPhi = atan(mcTrackMomentumTy / mcTrackMomentumTx);

    fPhiEfficiencyNorm->Fill(mcTrackPhi);
    fThetaEfficiencyNorm->Fill(mcTrackTheta);
    
    if (isAcc)
    {
	fPhiEfficiencyAcceptedNorm->Fill(mcTrackPhi);
	fThetaEfficiencyAcceptedNorm->Fill(mcTrackTheta);
    }
    */
}

void CbmSttFindTracksQa::WriteHistograms(char *filename)
{
    TFile
	*hfile = new TFile(filename,"RECREATE","histograms for CbmSttFindTracksQa");

    fPhiEfficiency->Write();
    fThetaEfficiency->Write();

    fPhiEfficiencyAccepted->Write();
    fThetaEfficiencyAccepted->Write();

    fPhiEfficiencyNorm->Write();
    fThetaEfficiencyNorm->Write();

    fPhiEfficiencyAcceptedNorm->Write();
    fThetaEfficiencyAcceptedNorm->Write();

    fPhiGhosts->Write();
    fThetaGhosts->Write();

    fPhiClones->Write();
    fThetaClones->Write();

    GoodTracks->Write();
    GhostTracks->Write();
    CloneTracks->Write();

    hfile->Close();
    delete hfile;
}

void CbmSttFindTracksQa::CreateHistograms()
{
    fPhiEfficiencyNorm = new TH1F("Number of thrown tracks (phi)", "Number of thrown tracks (phi)", 100, -3.20, 3.20);
    fThetaEfficiencyNorm = new TH1F("Number of thrown tracks (theta)", "Number of thrown tracks (theta)", 100, -0.10, 3.20);
    fPhiEfficiencyAcceptedNorm = new TH1F("Number of thrown tracks in acc. (phi)", "Number of thrown tracks in acc. (phi)", 100, -3.20, 3.20);
    fThetaEfficiencyAcceptedNorm = new TH1F("Number of thrown tracks in acc. (theta)", "Number of thrown tracks in acc. (theta)", 100, -0.10, 3.20);

    fPhiEfficiency = new TH1F("Number of rec tracks (phi)", "Number of rec tracks (phi)", 100, -3.20, 3.20);
    fThetaEfficiency = new TH1F("Number of rec tracks (theta)", "Number of rec tracks (theta)", 100, -0.10, 3.20);
    fPhiEfficiencyAccepted = new TH1F("Number of rec tracks in acc. (phi)", "Number of rec tracks in acc. (phi)", 100, -3.20, 3.20);
    fThetaEfficiencyAccepted = new TH1F("Number of rec tracks in acc. (theta)", "Number of rec tracks in acc. (theta)", 100, -0.10, 3.20);

    fPhiGhosts = new TH1F("Number of ghost tracks (phi)", "Number of ghost tracks (phi)", 100, -3.20, 3.20);
    fThetaGhosts = new TH1F("Number of ghost tracks (theta)", "Number of ghost tracks (theta)", 100, -0.10, 3.20);
    fPhiClones = new TH1F("Number of clone tracks (phi)", "Number of clone tracks (phi)", 100, -3.20, 3.20);
    fThetaClones = new TH1F("Number of clone tracks (theta)", "Number of clone tracks (theta)", 100, -0.10, 3.20);

    GoodTracks = new TH2F("Number of reconstructed vs. number of accepted tracks", 
			  "Number of reconstructed vs. number of accepted tracks", 11, -0.5, 10.5, 11, -0.5, 10.5);
    GhostTracks = new TH2F("Number of ghosts vs. number of accepted tracks", 
			   "Number of ghosts vs. number of accepted tracks", 11, -0.5, 10.5, 11, -0.5, 10.5);
    CloneTracks = new TH2F("Number of clones vs. number of accepted tracks", 
			   "Number of clones vs. number of accepted tracks", 11, -0.5, 10.5, 11, -0.5, 10.5);

}


ClassImp(CbmSttFindTracksQa)