// -------------------------------------------------------------------------
// -----                PndDchTrackFinderIdealCylHit source file              -----
// -----                       Created 27/03/08                         -----
// -------------------------------------------------------------------------

#include "PndDchTrackFinderIdealCylHit.h"

// Pnd includes
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairBaseParSet.h"
#include "PndDchPoint.h"
#include "PndDchDigi.h"
#include "PndDchCylinderHit.h"
#include "FairRootManager.h"
#include "PndDetectorList.h"
// ROOT includes
#include "TClonesArray.h"
#include "TGeoManager.h"

// C++ includes
#include <iostream>
#include <map>
#include <cmath>
using std::cout;
using std::endl;
using std::map;

// -----   Default constructor   -------------------------------------------
PndDchTrackFinderIdealCylHit::PndDchTrackFinderIdealCylHit() {
	fMCTrackArray = NULL;
	fMCPointArray = NULL;
	fDigiArray    = NULL;
	fNofEvents    = 0;
}

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

void PndDchTrackFinderIdealCylHit::Init() {

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

	// Get the pointer to the singleton FairRunAna object
	FairRunAna* ana = FairRunAna::Instance();
	if(NULL == ana) {
		cout << "-E- "<< GetName() <<"::Init :"
		<<" no FairRunAna object!" << endl;
		return;
	}
	
	// Get the pointer to run-time data base
	FairRuntimeDb* rtdb = ana->GetRuntimeDb();
	if(NULL == rtdb) {
		cout << "-E- "<< GetName() <<"::Init :"
		<<" no runtime database!" << endl;
		return;
	}

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

	// Get PndDchPoint (MCPoint) array
	fMCPointArray  = (TClonesArray*) ioman->GetObject("PndDchPoint");
	if( !fMCPointArray ) {
		cout << "-E- "<< GetName() <<"::Init: No MCPoint array!" << endl;
		return;
	}

	// Get PndDchDigi array
	fDigiArray  = (TClonesArray*) ioman->GetObject("PndDchDigi");
	if( !fDigiArray ) {
		cout << "-E- "<< GetName() <<"::Init: No Digi array!" << endl;
		return;
	}
	
	// Geometry loading
//	TFile *infile = ioman->GetInFile();
//	TGeoManager *geoMan = (TGeoManager*) infile->Get("FAIRGeom");
	fDchStructure = PndDchStructure::Instance(gGeoManager);

	std::cout << "-I- "<< GetName() <<": Intialization successfull" << std::endl;
}

// -----   Public method DoFind   ------------------------------------------
Int_t PndDchTrackFinderIdealCylHit::DoFind(TClonesArray* cylHitArray,
					   TClonesArray* trackArray) {
  
  // Count events
  fNofEvents++;
  if ( fVerbose ) {
    cout << endl << endl<< endl << endl;
    cout << "=======================================================" << endl;
    cout << "-I-        Event No: " <<  fNofEvents << endl;
    cout << "=======================================================" << endl;
    
    
    cout <<"-I- "<< GetName() <<"::DoFind "<< endl;
    cout << "-------------------------------------------------------" << endl;
    cout << "     ### Start DoFind" << endl;
    cout << "-------------------------------------------------------" << endl;
  }

  // Check pointers
  if( !fMCTrackArray ) {
    cout << "-E- "<< GetName() <<"::DoFind: "
	 << "MCTrack array missing! " << endl;
    return -1;
  } else {
    // Size of fMCTrackArray
    if ( fVerbose)    cout <<"# MCTracks (fMCTrackArray): "<< fMCTrackArray->GetEntriesFast() << endl;
  }
  
  if( !fMCPointArray ) {
    cout << "-E- "<< GetName() <<"::DoFind: "
	 << "MCPoint array missing! " << endl;
    return -1;
  } else {
    // Size of fMCPointArray
    if ( fVerbose)    cout <<"# MCPoints (fMCPointArray): "<< fMCPointArray->GetEntriesFast() << endl;
  }
  
  if( !fDigiArray ) {
    cout << "-E- "<< GetName() <<"::DoFind: "
	 << "Digi arrays missing! "<< endl;
    return -1;
  } else {
    // Size of fDigiArray
    if ( fVerbose)    cout <<"# Digis (fDigiArray): "<< fDigiArray->GetEntriesFast() << endl;
  }
  
  if( !cylHitArray ) {
    cout << "-E- "<< GetName() <<"::DoFind: "
	 << "CylHit arrays missing! "<< endl;
    return -1;
  } else {
    // Size of cylHitArray
    if ( fVerbose)    cout <<"# DchCylHits (cylHitArray): "<< cylHitArray->GetEntriesFast() << endl;
  }
  
  // Initialise control counters
  Int_t nNoMCTrack   = 0;
  Int_t nNoTrack     = 0;
  Int_t nNoDchPoint  = 0;
  Int_t nNoDchDigi  = 0;
  Int_t nNoDchCylHit   = 0;
  Int_t nDchCylHits6CH = 0;
  
  // Create pointers to DchCylinderHit/DchPoint/DchDigi
  PndDchCylinderHit*  dchCylHit  = NULL;
  PndDchDigi*  dchDigi  = NULL;
  FairMCPoint*  mcPoint  = NULL;
  PndMCTrack*  mcTrack  = NULL;
  PndTrackCand* dchTrack = NULL;
  
  // Declare variables outside the loop
  Int_t ptIndex = 0;       // MC point index
  Int_t mcTrackIndex = 0;  // MC track index
  Int_t trackIndex = 0;    // Dch track index
  
  
  // Create STL map from MCTrack index to number of digis
  std::map<Int_t, Int_t> cylHitMap;
  std::map<Int_t, Int_t>::iterator itCylHitMap;
  
  // Number of Dch cylHits
  Int_t nDchCylHits = cylHitArray->GetEntriesFast();
  if(fVerbose > 2) cout <<"# DchCylHits: "<< nDchCylHits << endl;
  
  
  for(Int_t iCylHit = 0; iCylHit < nDchCylHits; iCylHit++){
    // Get the pointer to Dch 
    dchCylHit = (PndDchCylinderHit*) cylHitArray->At(iCylHit);
    dchDigi = (PndDchDigi*) fDigiArray->At(dchCylHit->GetDigiIndex());
    
    if(dchDigi->GetChamber() > 2) {
      if(NULL == dchCylHit ) continue;
      nDchCylHits6CH++;
      
      // Get point index
      ptIndex = dchDigi->GetRefIndex();

      // Get pointer to MC point
      mcPoint = (FairMCPoint*) fMCPointArray->At(ptIndex);
      if(NULL == mcPoint) continue;
      
      // Get MC track index
      mcTrackIndex = mcPoint->GetTrackID();
      
      // Increment the number of digis
      cylHitMap[mcTrackIndex] += 1;
    }
  }
  
  
  // Create STL map from MCTrack index to DchTrack index
  map<Int_t, Int_t> trackMap;
  
  // Loop over reconstructable MCTracks and create corresponding DchTrack
  Int_t nMCacc = 0;  // Accepted MC tracks
  Int_t nTracks = 0; // Reconstructable MC tracks
  Int_t nMCTracks = fMCTrackArray->GetEntriesFast();
  
  if(fVerbose > 2) {
    cout <<"# of MC Tracks (nMCTracks): " << nMCTracks << endl;
  }
  
  for(Int_t iMCTrack = 0; iMCTrack < nMCTracks; iMCTrack++) {
    mcTrack = (PndMCTrack*) fMCTrackArray->At(iMCTrack);
    if( !mcTrack ) continue;
    
    
    // Either Primary or All MC tracks are considered
    if(fPrimary) {
      if(mcTrack->GetMotherID() != -1 && fabs((mcTrack->GetStartVertex()).Z())>1.) continue;
      
      if(fVerbose > 2) {
	cout <<"MC track #: "<< iMCTrack
	     << "\t MotherID: "<< mcTrack->GetMotherID()
	     << "\t StartVertex Z: "<< mcTrack->GetStartVertex().Z() << endl;
      }
    }
    //----------------------------------------------------
    // trivial check for number of DCH hits
    if ( !cylHitMap[iMCTrack] ) continue;
    //----------------------------------------------------
    nMCacc++;
    
    if(fVerbose > 2) {
      cout << iMCTrack << ": #DchPoints in MCTrack: "<< mcTrack->GetNPoints(kDCH) << endl;
    }
    
    new((*trackArray)[nTracks]) PndTrackCand();
    
    dchTrack = (PndTrackCand*) trackArray->At(nTracks);
    /*		TLorentzVector tlVec = mcTrack->Get4Momentum();
		Double_t mom = tlVec.Mag();
		if(mom != 0) {
		dchTrack->GetParamLast()->SetQp(1./mom);
		} else {
		dchTrack->GetParamLast()->SetQp(0);
		}*/
    
    dchTrack->setMcTrackId(iMCTrack);
    
    trackMap[iMCTrack] = nTracks;
    
    nTracks++;
  }

  // Loop over cylHits. Get corresponding MCPoint and MCTrack index
  for(Int_t iCylHit = 0; iCylHit < nDchCylHits; iCylHit++) {
    dchCylHit = (PndDchCylinderHit*) cylHitArray->At(iCylHit);
    if( !dchCylHit ) {
      nNoDchCylHit++;
      continue;
    }
    
    dchDigi = (PndDchDigi*) fDigiArray->At(dchCylHit->GetDigiIndex());
    ptIndex = dchDigi->GetRefIndex();
    
    if(ptIndex < 0) continue;           // fake or background digi

    mcPoint = (FairMCPoint*) fMCPointArray->At(ptIndex);

    if( !mcPoint ) {
      nNoDchPoint++;
      continue;
    }
    
    mcTrackIndex = mcPoint->GetTrackID();

    if(mcTrackIndex < 0 || mcTrackIndex > nMCTracks) {
      cout << "-E- "<< GetName() <<"::DoFind: "
	   << "MCTrack index out of range. " << mcTrackIndex << " "
	   << nMCTracks << endl;
      nNoMCTrack++;
      continue;
    }
    
    if(trackMap.find(mcTrackIndex) == trackMap.end()) continue;
    trackIndex = trackMap[mcTrackIndex];

    dchTrack = (PndTrackCand*) trackArray->At(trackIndex);
    
    if( !dchTrack ) {
      cout << "-E- "<< GetName() <<"::DoFind: "
	   << "No DchTrack pointer. " << iCylHit << " " << ptIndex
	   << " " << mcTrackIndex << " " << trackIndex << endl;
      nNoTrack++;
      continue;
    }
    
    TVector3 pos;
    mcPoint->Position(pos);

    dchTrack->AddHit(kDchHit,iCylHit,pos.Mag());
    
    if(fVerbose > 3) {
      cout << "DCH cylHit " << iCylHit << " from DCH point "
	   << ptIndex << " (" << mcTrackIndex << ") "
	   << "added to DCH track " << trackIndex << endl;
    }
  }
   
  // Sorting hits
  for(Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
    dchTrack = (PndTrackCand*) trackArray->At(iTrack);
    dchTrack->Sort();
    
    PndTrackCandHit tch = dchTrack->GetSortedHit(0);

    dchCylHit = (PndDchCylinderHit*) cylHitArray->At(tch.GetHitId());
    dchDigi = (PndDchDigi*) fDigiArray->At(dchCylHit->GetDigiIndex());
    ptIndex = dchDigi->GetRefIndex();

    TVector3 pos;
    TVector3 mom;
    
    mcPoint = (FairMCPoint*) fMCPointArray->At(ptIndex);

    mcPoint->Position(pos);
    mcPoint->Momentum(mom);

    dchTrack->setTrackSeed(pos,mom.Unit(),1./mom.Mag());
  }

  if(fVerbose) {
    cout << endl;
    cout << "-------------------------------------------------------" << endl;
    cout << "-I-    "<< GetName() <<": Event summary    -I-" << endl;
    cout << "-------------------------------------------------------" << endl;
    cout << "Total Dch cylHits: " << nDchCylHits <<"   In 6CH: "<< nDchCylHits6CH << endl;
    cout << "MC tracks total: " << nMCTracks << ", accepted: "
	 << nMCacc << ", reconstructable: " << nTracks << endl;
    if(nNoDchCylHit)  cout << "DchCylHits not found  : " << nNoDchCylHit  << endl;
    if(nNoDchPoint) cout << "DchPoints not found : " << nNoDchPoint << endl;
    if(nNoMCTrack)  cout << "MCTracks not found  : " << nNoMCTrack  << endl;
    if(nNoTrack)    cout << "DchTracks not found : " << nNoTrack    << endl;
    cout << "------------------------------------------------------" << endl;
    cout << endl;
  } else {
//     cout << "All: "        << nMCTracks
// 	 << ", Accepted: "      << nMCacc
// 	 << ", Reconstructed: " << nTracks << endl;
  }
  
  return nTracks;
}


ClassImp(PndDchTrackFinderIdealCylHit)