CbmStsIdealMatchHits.cxx

Go to the documentation of this file.

//* $Id: */
 
// -------------------------------------------------------------------------
// -----                    CbmStsIdealMatchHits source file                -----
// -----                  Created 27/11/06  by V. Friese               -----
// -------------------------------------------------------------------------
#include <iomanip>
#include <iostream>
 
// --- 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 "CbmStsDigiMatch.h"
#include "CbmStsDigiPar.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsHit.h"
#include "CbmStsIdealMatchHits.h"
#include "CbmStsPoint.h"
#include "CbmStsSector.h"
#include "CbmStsStation.h"
#include "TMath.h"
 
using std::cout;
using std::endl;
using std::fixed;
using std::left;
using std::right;
using std::setprecision;
using std::setw;
 
// -----   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 = 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 jMatchF = 0; jMatchF < 3; jMatchF++) {
                Int_t iPointF = dMatchF->GetRefIndex(jMatchF);
                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 < 3; jMatchF++) {
                Int_t iPointF = dMatchF->GetRefIndex(jMatchF);
                if (iPointF < 0)
                    continue;
                nPointsF++;
                for (Int_t jMatchB = 0; jMatchB < 3; jMatchB++) {
                    Int_t iPointB = dMatchB->GetRefIndex(jMatchB);
                    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->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 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;
}
// -------------------------------------------------------------------------