//* $Id: */

// -------------------------------------------------------------------------
// -----                    CbmStsMatchHits source file                -----
// -----                  Created 27/11/06  by V. Friese               -----
// -------------------------------------------------------------------------
#include <iostream>
#include <iomanip>

// --- Includes from ROOT
#include "TClonesArray.h"

// --- Includes from base
#include "CbmRootManager.h"
#include "CbmRunAna.h"
#include "CbmRuntimeDb.h"

// --- Includes from STS
#include "CbmGeoStsPar.h"
#include "CbmStsDigiPar.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsDigiMatch.h"
#include "CbmStsHit.h"
#include "CbmStsMatchHits.h"
#include "CbmStsPoint.h"
#include "CbmStsSector.h"
#include "CbmStsStation.h"

#include "TMath.h"

using std::cout;
using std::endl;
using std::setw;
using std::left;
using std::right;
using std::fixed;
using std::setprecision;


// -----   Default constructor   -------------------------------------------
CbmStsMatchHits::CbmStsMatchHits() : CbmTask("STSMatchHits", 1) {
  fGeoPar      = NULL;
  fDigiPar     = NULL;
  fPoints      = NULL;
  fDigis       = NULL;
  fDigiMatches = NULL;
  fHits        = NULL;
  fDigiScheme  = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------



// -----   Standard constructor   ------------------------------------------
CbmStsMatchHits::CbmStsMatchHits(Int_t iVerbose)
  : CbmTask("STSMatchHits", iVerbose) {
  fGeoPar      = NULL;
  fDigiPar     = NULL;
  fPoints      = NULL;
  fDigis       = NULL;
  fDigiMatches = NULL;
  fHits        = NULL;
  fDigiScheme  = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------



// -----   Constructor with name   -----------------------------------------
CbmStsMatchHits::CbmStsMatchHits(const char* name, Int_t iVerbose)
  : CbmTask(name, iVerbose) {
  fGeoPar      = NULL;
  fDigiPar     = NULL;
  fPoints      = NULL;
  fDigis       = NULL;
  fDigiMatches = NULL;
  fHits        = NULL;
  fDigiScheme  = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
CbmStsMatchHits::~CbmStsMatchHits() {
  if ( fGeoPar )     delete fGeoPar;
  if ( fDigiPar )    delete fDigiPar;
  if ( fDigiScheme ) delete fDigiScheme;
}
// -------------------------------------------------------------------------



// -----   Public method Exec   --------------------------------------------
void CbmStsMatchHits::Exec(Option_t* opt) {

  // Timer
  fTimer.Start();
  Bool_t warn = kFALSE;

  // Counters
  Int_t nHits    = fHits->GetEntriesFast();
  Int_t nNoDigi  = 0;
  Int_t nBackgrd = 0;
  Int_t nDistant = 0;
  Int_t nMatched = 0;

  // Loop over all StsHits
  for (Int_t iHit=0; iHit<nHits; iHit++) {
    CbmStsHit* hit = (CbmStsHit*) fHits->At(iHit);
    if ( ! hit ) {
      cout << "-W- " << GetName() << "::Exec: Empty hit at index "
	   << iHit << endl;
      warn = kTRUE;
      continue;
    }

    // Determine sector type and channel numbers
    Int_t iStation = hit->GetStationNr();
    Int_t iSector  = hit->GetSectorNr();
    CbmStsStation* station = fDigiScheme->GetStationByNr(iStation);
    CbmStsSector*  sector = fDigiScheme->GetSector(iStation, iSector);
    Int_t iType = sector->GetType();
    CbmStsDigiMatch* dMatchF = NULL;
    CbmStsDigiMatch* dMatchB = NULL;
    

    // Hit coordinates and errors
    Double_t xH = hit->GetX();
    Double_t yH = hit->GetY();
    Double_t dX = hit->GetDx();
    Double_t dY = hit->GetDy();

    // Get front side DigiMatch corresponding to hit
    Int_t  iMatchF  = hit->GetDigi(0);
    if ( iMatchF >= 0 ) 
      dMatchF = (CbmStsDigiMatch*) fDigiMatches->At(iMatchF);
    if ( ! dMatchF ) {
      cout << "-E- " << GetName() << "::Exec: "
	   << "No DigiMatchF for hit " << iHit << endl;
      hit->SetRefIndex(-1);
      nNoDigi++;
      warn = kTRUE;
      continue;
    }

    // Get back side DigiMatch of hit (for strip sensors only)
    if ( iType != 1 ) {
      Int_t  iMatchB  = hit->GetDigi(1);    
      if ( iMatchB >= 0 ) 
	dMatchB = (CbmStsDigiMatch*) fDigiMatches->At(iMatchB);
      if ( ! dMatchB ) {
	cout << "-E- " << GetName() << "::Exec: "
	     << "No DigiMatchB for hit " << iHit << endl;
	hit->SetRefIndex(-1);
	nNoDigi++;
	warn = kTRUE;
	continue;
      }
    }

    // Map candidate points to distance to hit
    fCandMap.clear();
    Int_t nPointsF = 0;
    Int_t nPointsB = 0;

    if ( iType == 1 ) {
      // Case pixel: Candidates are all points of the digi.
      for (Int_t iMatchF=0; iMatchF<3; iMatchF++) {
	Int_t iPointF = dMatchF->GetRefIndex(iMatchF);
	if ( iPointF < 0 ) continue;
	nPointsF++;
	// Calculate distance to hit
	CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPointF);
	if ( ! point ) {
	  cout << "-E- " << GetName() << "::Exec: "
	       << "No StsPoint (" << iPointF << ") for pixel hit " 
	       << iHit << endl;
	  warn = kTRUE;
	  continue;
	}
	// point coordinates in station midplane
	Double_t xP = point->GetX(station->GetZ());
	Double_t yP = point->GetY(station->GetZ());
	Double_t dist = TMath::Sqrt( (xP-xH)*(xP-xH) + (yP-yH)*(yP-yH) );
	fCandMap[dist] = iPointF;
      }     // front digi loop
    }       // pixel sensor

    else if ( iType == 2 || iType == 3 ) {
      // Case strip: Candidates are points corresponding to front
      // and back side digis
      for ( Int_t iMatchF=0; iMatchF<3; iMatchF++) {
	Int_t iPointF = dMatchF->GetRefIndex(iMatchF);
	if ( iPointF < 0 ) continue; 
	nPointsF++;
	for ( Int_t iMatchB=0; iMatchB<3; iMatchB++) {
	  Int_t iPointB = dMatchB->GetRefIndex(iMatchB);
	  if ( iPointB < 0 ) continue; 
	  if ( iMatchF == 0 ) nPointsB++;
	  if ( iPointB != iPointF ) continue; // chance combination
	  // Calculate distance to hit
	  CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPointF);
	  if ( ! point ) {
	    cout << "-E- " << GetName() << "::Exec: "
		 << "No StsPoint (" << iPointF << ") for strip hit " 
		 << iHit << endl;
	    warn = kTRUE;
	    continue;
	  }
	  // point coordinates in station midplane
	  Double_t xP = point->GetX(station->GetZ());
	  Double_t yP = point->GetY(station->GetZ());
	  Double_t dist = TMath::Sqrt( (xP-xH)*(xP-xH) + 
				       (yP-yH)*(yP-yH) );
	  fCandMap[dist] = iPointF;
	}   // back digi loop
      }     // front digi loop
    }       // strip sensor

    else {
      // Unknown sensor type
      cout << "-E- " << GetName() << "::Exec: Unknwon sensor type "
	   << iType << endl;
      Fatal("Exec", "Unknwon sensor type");
    }

    if ( fVerbose > 1 ) cout << "-I- " << GetName() << ": Hit " 
			     << iHit << ", type " << iType 
			     << ", points " << nPointsF << " / " 
			     << nPointsB << ", candidates " 
			     << fCandMap.size();


    // Check for at least one candidate point. Else background hit.
    // This can happen for noise digis or for fake combinations
    // of strip digis.
    if ( fCandMap.empty() ) {
      hit->SetRefIndex(-1);
      if (fVerbose>1) cout << ", background " << endl;
      nBackgrd++;
      continue;
    }
      
    // Select closest point from candidates.
    Double_t distMin = 999999.;
    Int_t    iPoint  = -1;
    for (fIter=fCandMap.begin(); fIter!=fCandMap.end(); fIter++) {
      if ( (*fIter).first < distMin ) {
	distMin = (*fIter).first;
	iPoint  = (*fIter).second;
      }
    }
    if ( iPoint < 0 ) {
      cout << "-E- " << GetName() << "::Exec: "
	   << "No closest point found in candidate map!" << endl;
      Fatal("Exec", "No closest point");
    }
    CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPoint);
    if (fVerbose>1) cout << ", matched to " << iPoint << ", distance " 
			 << distMin << " cm";
 	    
    // Check whether closest point is more than 5 sigma away from hit.
    // This should not happen in case of pixel or strip OSU sensors,
    // but can happen for strip GSI sensors due to multiple
    // crossings of front and back side strips.
    Double_t xP = point->GetX(station->GetZ());
    Double_t yP = point->GetY(station->GetZ());
    if ( TMath::Abs(xP-xH) > 5. * dX || 
	 TMath::Abs(yP-yH) > 5. * dY ) {
      hit->SetRefIndex(-2);
      nDistant++;
      if (fVerbose>1) cout << ", distant" << endl;
      if ( iType == 1 || iType == 2) {
	cout << "-W- " << fName << "::Exec: "
	     << "Distant hit in pixel / strip MSU station" << endl;
	cout << "Hit " << iHit << " coordinates " << xH
	     << ", " << yH << endl;
	cout << "Matched point " << iPoint << ", distance "
	     << distMin << endl;
	cout << "Distance y " << TMath::Abs(xP-xH) << " Error " << dX << endl;
	cout << "Distance y " << TMath::Abs(yP-yH) << " Error " << dY << endl;
      }

      
      continue;
    }

    // Match closest StsPoint to hit
    hit->SetRefIndex(iPoint);
    nMatched++;
    if (fVerbose>1) cout << ", good match" << endl;

    // Check whether hit and point are in the same station.
    // Else something has really gone wrong.
    if ( TMath::Abs(hit->GetZ() - point->GetZ()) > 1. ) {
      cout << "-E- " << GetName() << "::Exec: Hit " << iHit
	   << " is at z = " << hit->GetZ() << " cm, but matched point "
	   << iPoint << " is at z = " << point->GetZ() << "cm " << endl;
      Fatal("Exec", "Different stations for hit and point");
    }

  }   // Loop over StsHits

  // Event statistics
  fTimer.Stop();
  if ( fVerbose ) {
    if ( warn ) cout << "- ";
    else        cout << "+ ";
    cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8) 
	 << fixed << right << fTimer.RealTime()
	 << " s, hits " << nHits << ", matched " << nMatched << ", distant " 
	 << nDistant << ", background " << setw(6) << nBackgrd << endl;
  }

  // Run statistics
  if ( warn) fNEventsFailed++;
  else {
    fNEvents++;
    fTime     += fTimer.RealTime();
    fNHits    += Double_t(nHits);
    fNMatched += Double_t(nMatched);
    fNDistant += Double_t(nDistant);
    fNBackgrd += Double_t(nBackgrd);
  }

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



// -----   Private method SetParContainers   -------------------------------
void CbmStsMatchHits::SetParContainers() {

  // Get run and runtime database
  CbmRunAna* run = CbmRunAna::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");
  CbmRuntimeDb* db = run->GetRuntimeDb();
  if ( ! db ) Fatal("SetParContainers", "No runtime database");

  // Get STS geometry and digitisation parameter container
  fGeoPar  = (CbmGeoStsPar*)  db->getContainer("CbmGeoStsPar");
  fDigiPar = (CbmStsDigiPar*) db->getContainer("CbmStsDigiPar");

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



// -----   Private method Init   -------------------------------------------
InitStatus CbmStsMatchHits::Init() {

  // Reset counters
  fNEvents = fNEventsFailed = 0;
  fTime = fNHits = fNMatched = fNDistant = fNBackgrd = 0.;

  // Get input arrays
  CbmRootManager* ioman = CbmRootManager::Instance();
  if ( ! ioman ) Fatal("Init", "No CbmRootManager");
  fPoints = (TClonesArray*) ioman->GetObject("STSPoint");
  if ( ! fPoints ) {
    cout << "-E- " << GetName() << "::Init: No STSPoint array!" << endl;
    return kERROR;
  }
  fDigis  = (TClonesArray*) ioman->GetObject("STSDigi");
  if ( ! fDigis ) {
    cout << "-E- " << GetName() << "::Init: No STSDigi array!" << endl;
    return kERROR;
  }
  fDigiMatches = (TClonesArray*) ioman->GetObject("STSDigiMatch");
  if ( ! fDigiMatches ) {
    cout << "-E- " << GetName() << "::Init: No STSDigiMatch array!" 
	 << endl;
    return kERROR;
  }
  fHits = (TClonesArray*) ioman->GetObject("STSHit");
  if ( ! fHits ) {
    cout << "-E- " << GetName() << "::Init: No STSHit array!" << endl;
    return kERROR;
  }

  // Build digitisation scheme
  if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) {
    if      (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
    else if (fVerbose >  2) fDigiScheme->Print(kTRUE);
    cout << "-I- " << fName << "::Init: "
	 << "STS digitisation scheme succesfully initialised" << endl;
    cout << "    Stations: " << fDigiScheme->GetNStations() 
	 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: " 
	 << fDigiScheme->GetNChannels() << endl;
    return kSUCCESS;
  }

  return kERROR;

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



// -----   Private method ReInit   -----------------------------------------
InitStatus CbmStsMatchHits::ReInit() {

  // Clear digitisation scheme
  fDigiScheme->Clear();

  // Build new digitisation scheme
  if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) return kSUCCESS;

  return kERROR;

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



// -----   Private method Finish   -----------------------------------------
void CbmStsMatchHits::Finish() {

  cout << endl;
  cout << "============================================================"
       << endl;
  cout << "=====   " << GetName() << ": Run summary " << endl;
  cout << "===== " << endl;
  cout << "===== Good events   : " << setw(6) << fNEvents << endl;
  cout << "===== Failed events : " << setw(6) << fNEventsFailed << endl;
  cout << "===== Average time  : " << setprecision(4) << setw(8) << right
       << fTime / Double_t(fNEvents)  << " s" << endl;
  cout << "===== " << endl;
  cout << "===== Hits per event  : " << fixed << setprecision(0)
       << fNHits / Double_t(fNEvents) << endl;
  cout << setprecision(2);
  cout << "===== Matched hits    : " << fixed << setw(6) << right
       << fNMatched / fNHits * 100. << " %" << endl;
  cout << "===== Distant hits    : " << fixed << setw(6) << right
       << fNDistant / fNHits * 100. << " %" << endl;
  cout << "===== Background hits : " << fixed << setw(6) << right
       << fNBackgrd / fNHits * 100. << " % " << endl;
  cout << "============================================================"
       << endl;

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



ClassImp(CbmStsMatchHits)