//* $Id: */

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

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

// --- Includes from base
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

// --- Includes from STS
#include "CbmGeoStsPar.h"
#include "CbmStsDigiPar.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsHit.h"
#include "CbmStsIdealMatchHits.h"
#include "CbmStsPoint.h"
#include "CbmStsSector.h"
#include "CbmStsStation.h"
#include "CbmStsAddress.h"
#include "CbmMatch.h"
#include "CbmLink.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   -------------------------------------------
CbmStsIdealMatchHits::CbmStsIdealMatchHits()
    : FairTask("STSMatchHits", 1)
    , fGeoPar(NULL)
    , fDigiPar(NULL)
    , fPoints(NULL)
    , fDigis(NULL)
    , fDigiMatches(NULL)
    , fHits(NULL)
    , fDigiScheme(NULL)
    , fTimer()
    , fNEvents(0)
    , fNEventsFailed(0)
    , fTime(0.)
    , fNHits(0.)
    , fNMatched(0.)
    , fNDistant(0.)
    , fNBackgrd(0.)
    , fCandMap()
    , fIter()
{
    fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------

// -----   Standard constructor   ------------------------------------------
CbmStsIdealMatchHits::CbmStsIdealMatchHits(Int_t iVerbose)
    : FairTask("STSMatchHits", iVerbose)
    , fGeoPar(NULL)
    , fDigiPar(NULL)
    , fPoints(NULL)
    , fDigis(NULL)
    , fDigiMatches(NULL)
    , fHits(NULL)
    , fDigiScheme(NULL)
    , fTimer()
    , fNEvents(0)
    , fNEventsFailed(0)
    , fTime(0.)
    , fNHits(0.)
    , fNMatched(0.)
    , fNDistant(0.)
    , fNBackgrd(0.)
    , fCandMap()
    , fIter()
{
    fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------

// -----   Constructor with name   -----------------------------------------
CbmStsIdealMatchHits::CbmStsIdealMatchHits(const char* name, Int_t iVerbose)
    : FairTask(name, iVerbose)
    , fGeoPar(NULL)
    , fDigiPar(NULL)
    , fPoints(NULL)
    , fDigis(NULL)
    , fDigiMatches(NULL)
    , fHits(NULL)
    , fDigiScheme(NULL)
    , fTimer()
    , fNEvents(0)
    , fNEventsFailed(0)
    , fTime(0.)
    , fNHits(0.)
    , fNMatched(0.)
    , fNDistant(0.)
    , fNBackgrd(0.)
    , fCandMap()
    , fIter()
{
    fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------

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

// -----   Public method Exec   --------------------------------------------
void CbmStsIdealMatchHits::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 = CbmStsAddress::GetElementId(hit->GetAddress(), kStsStation) + 1; // hit->GetStationNr();
        Int_t iSector = hit->GetSectorNr();
        CbmStsStation* station = fDigiScheme->GetStationByNr(iStation);
        CbmStsSector* sector = fDigiScheme->GetSector(iStation, iSector);
        Int_t iType = sector->GetType();
        CbmMatch* dMatchF = NULL;
        CbmMatch* 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->GetFrontDigiId();
        if (iMatchF >= 0)
            dMatchF = (CbmMatch*)fDigiMatches->At(iMatchF);
        if (!dMatchF)
        {
            cout << "-E- " << GetName() << "::Exec: "
                 << "No DigiMatchF for hit " << iHit << endl;
            hit->SetRefId(-1);
            nNoDigi++;
            warn = kTRUE;
            continue;
        }

        // Get back side DigiMatch of hit (for strip sensors only)
        if (iType != 1)
        {
            Int_t iMatchB = hit->GetBackDigiId();
            if (iMatchB >= 0)
                dMatchB = (CbmMatch*)fDigiMatches->At(iMatchB);
            if (!dMatchB)
            {
                cout << "-E- " << GetName() << "::Exec: "
                     << "No DigiMatchB for hit " << iHit << endl;
                hit->SetRefId(-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 jMatchF = 0; jMatchF < dMatchF->GetNofLinks(); jMatchF++)
            {
                Int_t iPointF = dMatchF->GetLink(jMatchF).GetEntry();
                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 jMatchF = 0; jMatchF < dMatchF->GetNofLinks(); jMatchF++)
            {
                Int_t iPointF = dMatchF->GetLink(jMatchF).GetEntry();
                if (iPointF < 0)
                    continue;
                nPointsF++;
                for (Int_t jMatchB = 0; jMatchB < dMatchB->GetNofLinks(); jMatchB++)
                {
                    Int_t iPointB = dMatchB->GetLink(jMatchB).GetEntry();
                    if (iPointB < 0)
                        continue;
                    if (jMatchF == 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->SetRefId(-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->SetRefId(-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->SetRefId(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 CbmStsIdealMatchHits::SetParContainers()
{

    // Get run and runtime database
    FairRunAna* run = FairRunAna::Instance();
    if (!run)
        Fatal("SetParContainers", "No analysis run");
    FairRuntimeDb* 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 CbmStsIdealMatchHits::Init()
{

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

    // Get input arrays
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman)
        Fatal("Init", "No FairRootManager");
    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 CbmStsIdealMatchHits::ReInit()
{

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

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

    return kERROR;
}
// -------------------------------------------------------------------------

// -----   Private method Finish   -----------------------------------------
void CbmStsIdealMatchHits::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(CbmStsIdealMatchHits)