CbmStsRealFindHits.cxx

Go to the documentation of this file.

//* $Id: */
 
// -------------------------------------------------------------------------
// -----                    CbmStsRealFindHits source file             -----
// -----                  Created 26/06/2008 by R. Karabowicz          -----
// -------------------------------------------------------------------------
 
#include "CbmStsRealFindHits.h"
 
#include "CbmGeoStsPar.h"
#include "CbmStsCluster.h"
#include "CbmStsDigiPar.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsHit.h"
#include "CbmStsSector.h"
#include "CbmStsStation.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "TClonesArray.h"
#include "TMath.h"
 
#include <iomanip>
 
using std::cerr;
using std::cout;
using std::endl;
using std::fixed;
using std::flush;
using std::left;
using std::map;
using std::right;
using std::set;
using std::setprecision;
using std::setw;
 
// -----   Default constructor   ------------------------------------------
CbmStsRealFindHits::CbmStsRealFindHits()
    : FairTask("STS Hit Finder", 1)
{
    fGeoPar = NULL;
    fDigiPar = NULL;
    fClusters = NULL;
    fHits = NULL;
    fNHits = 0;
    fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------
 
// -----   Standard constructor   ------------------------------------------
CbmStsRealFindHits::CbmStsRealFindHits(Int_t iVerbose)
    : FairTask("STSRealFindHits", iVerbose)
{
    fGeoPar = NULL;
    fDigiPar = NULL;
    fClusters = NULL;
    fHits = NULL;
    fNHits = 0;
    fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------
 
// -----   Constructor with name   -----------------------------------------
CbmStsRealFindHits::CbmStsRealFindHits(const char* name, Int_t iVerbose)
    : FairTask(name, iVerbose)
{
    fGeoPar = NULL;
    fDigiPar = NULL;
    fClusters = NULL;
    fHits = NULL;
    fNHits = 0;
    fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
CbmStsRealFindHits::~CbmStsRealFindHits()
{
    if (fHits) {
        fHits->Delete();
        delete fHits;
    }
}
// -------------------------------------------------------------------------
 
// -----   Public method Exec   --------------------------------------------
void CbmStsRealFindHits::Exec(Option_t* opt)
{
 
    fTimer.Start();
    Bool_t warn = kFALSE;
 
    // Check for digi scheme
    if (!fDigiScheme) {
        cerr << "-E- " << fName << "::Exec: No digi scheme!" << endl;
        return;
    }
 
    // Clear output array
    //  fHits->Clear();
    fHits->Delete();
 
    // Sort STS digis with respect to sectors
    SortClusters();
 
    // Find hits in sectors
    Int_t nClustersF = 0;
    Int_t nClustersB = 0;
    Int_t nHits = 0;
    Int_t nStations = fDigiScheme->GetNStations();
    CbmStsStation* station = NULL;
    for (Int_t iStation = 0; iStation < nStations; iStation++) {
        station = fDigiScheme->GetStation(iStation);
 
        Int_t nClustersFInStation = 0;
        Int_t nClustersBInStation = 0;
        Int_t nHitsInStation = 0;
        Int_t nSectors = station->GetNSectors();
        for (Int_t iSector = 0; iSector < nSectors; iSector++) {
            CbmStsSector* sector = station->GetSector(iSector);
            set<Int_t> fSet, bSet;
            if (fClusterMapF.find(sector) == fClusterMapF.end()) {
                cout << "-E- " << fName << "::Exec: sector " << sector->GetSectorNr() << " of station "
                     << station->GetStationNr() << "not found in front map!" << endl;
                warn = kTRUE;
                continue;
            }
            fSet = fClusterMapF[sector];
            if (sector->GetType() == 2 || sector->GetType() == 3) {
                if (fClusterMapB.find(sector) == fClusterMapB.end()) {
                    cout << "-E- " << fName << "::Exec: sector " << sector->GetSectorNr() << " of station "
                         << station->GetStationNr() << "not found in back map!" << endl;
                    warn = kTRUE;
                    continue;
                }
            }
            bSet = fClusterMapB[sector];
            Int_t nClustersFInSector = fSet.size();
            Int_t nClustersBInSector = bSet.size();
            Int_t nHitsInSector = FindHits(station, sector, fSet, bSet);
            if (fVerbose > 2)
                cout << "Sector " << sector->GetSectorNr() << ", Clusters front " << nClustersFInSector
                     << ", Clusters Back " << nClustersBInSector << ", Hits " << nHitsInSector << endl;
            nHitsInStation += nHitsInSector;
            nClustersFInStation += nClustersFInSector;
            nClustersBInStation += nClustersBInSector;
        }   // Sector loop
 
        if (fVerbose > 1)
            cout << "Total for station " << station->GetStationNr() << ": Clusters front " << nClustersFInStation
                 << ", clusters back " << nClustersBInStation << ", hits " << nHitsInStation << endl;
        nClustersB += nClustersBInStation;
        nClustersF += nClustersFInStation;
        nHits += nHitsInStation;
 
    }   // Station loop
 
    fTimer.Stop();
    if (fVerbose) {
        cout << endl;
        cout << "-I- " << fName << ":Event summary" << endl;
        cout << "    Clusters front side       : " << nClustersF << endl;
        cout << "    Clusters back side        : " << nClustersB << endl;
        cout << "    Hits created              : " << nHits << endl;
        cout << "    Real time                 : " << fTimer.RealTime() << endl;
    } else {
        if (warn)
            cout << "- ";
        else
            cout << "+ ";
        cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8) << fixed << right << fTimer.RealTime()
             << " s, clusters " << nClustersF << " / " << nClustersB << ", hits: " << nHits << endl;
    }
 
    fNHits += nHits;
}
// -------------------------------------------------------------------------
 
// -----   Private method SetParContainers   -------------------------------
void CbmStsRealFindHits::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 parameter container
    fGeoPar = (CbmGeoStsPar*)db->getContainer("CbmGeoStsPar");
 
    // Get STS digitisation parameter container
    fDigiPar = (CbmStsDigiPar*)db->getContainer("CbmStsDigiPar");
}
// -------------------------------------------------------------------------
 
// -----   Private method Init   -------------------------------------------
InitStatus CbmStsRealFindHits::Init()
{
 
    // Get input array
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman)
        Fatal("Init", "No FairRootManager");
    fClusters = (TClonesArray*)ioman->GetObject("StsCluster");
 
    // Register output array
    fHits = new TClonesArray("CbmStsHit", 1000);
    ioman->Register("StsHit", "Hit in STS", fHits, kTRUE);
 
    // Build digitisation scheme
    Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
    if (!success)
        return kERROR;
 
    // Create sectorwise cluster sets
    MakeSets();
 
    // Control output
    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;
}
// -------------------------------------------------------------------------
 
// -----   Private method ReInit   -----------------------------------------
InitStatus CbmStsRealFindHits::ReInit()
{
 
    // Clear digitisation scheme
    fDigiScheme->Clear();
 
    // Build new digitisation scheme
    Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
    if (!success)
        return kERROR;
 
    // Create sectorwise digi sets
    MakeSets();
 
    // Control output
    if (fVerbose == 1 || fVerbose == 2)
        fDigiScheme->Print(kFALSE);
    else if (fVerbose > 2)
        fDigiScheme->Print(kTRUE);
    cout << "-I- " << fName << "::Init: "
         << "STS digitisation scheme succesfully reinitialised" << endl;
    cout << "    Stations: " << fDigiScheme->GetNStations() << ", Sectors: " << fDigiScheme->GetNSectors()
         << ", Channels: " << fDigiScheme->GetNChannels() << endl;
 
    return kSUCCESS;
}
// -------------------------------------------------------------------------
 
// -----   Private method MakeSets   ---------------------------------------
void CbmStsRealFindHits::MakeSets()
{
 
    fClusterMapF.clear();
    fClusterMapB.clear();
    Int_t nStations = fDigiScheme->GetNStations();
    for (Int_t iStation = 0; iStation < nStations; iStation++) {
        CbmStsStation* station = fDigiScheme->GetStation(iStation);
        Int_t nSectors = station->GetNSectors();
        for (Int_t iSector = 0; iSector < nSectors; iSector++) {
            CbmStsSector* sector = station->GetSector(iSector);
            set<Int_t> a;
            fClusterMapF[sector] = a;
            if (sector->GetType() == 2 || sector->GetType() == 3) {
                set<Int_t> b;
                fClusterMapB[sector] = b;
            }
        }
    }
}
// -------------------------------------------------------------------------
 
// -----   Private method SortClusters   --------------------------------------
void CbmStsRealFindHits::SortClusters()
{
 
    // Check input array
    if (!fClusters) {
        cout << "-E- " << fName << "::SortClusters: No input array!" << endl;
        return;
    }
 
    // Clear sector cluster sets
    map<CbmStsSector*, set<Int_t>>::iterator mapIt;
    for (mapIt = fClusterMapF.begin(); mapIt != fClusterMapF.end(); mapIt++)
        ((*mapIt).second).clear();
    for (mapIt = fClusterMapB.begin(); mapIt != fClusterMapB.end(); mapIt++)
        ((*mapIt).second).clear();
 
    // Fill clusters into sets
    CbmStsCluster* cluster = NULL;
    CbmStsSector* sector = NULL;
    Int_t stationNr = -1;
    Int_t sectorNr = -1;
    Int_t iSide = -1;
    Int_t nClusters = fClusters->GetEntriesFast();
    for (Int_t iClus = 0; iClus < nClusters; iClus++) {
        cluster = (CbmStsCluster*)fClusters->At(iClus);
        stationNr = cluster->GetStationNr();
        sectorNr = cluster->GetSectorNr();
        iSide = cluster->GetSide();
        sector = fDigiScheme->GetSector(stationNr, sectorNr);
        if (iSide == 0) {
            if (fClusterMapF.find(sector) == fClusterMapF.end()) {
                cerr << "-E- " << fName << "::SortClusters:: sector " << sectorNr << " of station " << stationNr
                     << " not found in digi scheme (F)!" << endl;
                continue;
            }
            fClusterMapF[sector].insert(iClus);
        } else if (iSide == 1) {
            if (fClusterMapB.find(sector) == fClusterMapB.end()) {
                cerr << "-E- " << fName << "::SortClusters:: sector " << sectorNr << " of station " << stationNr
                     << " not found in digi scheme (B)!" << endl;
                continue;
            }
            fClusterMapB[sector].insert(iClus);
        }
    }
}
// -------------------------------------------------------------------------
 
// -----   Private method FindHits   ---------------------------------------
Int_t CbmStsRealFindHits::FindHits(CbmStsStation* station, CbmStsSector* sector, set<Int_t>& fSet, set<Int_t>& bSet)
{
 
    // Counter
    Int_t nNew = 0;
 
    // Get sector parameters
    Int_t detId = sector->GetDetectorId();
    Int_t iType = sector->GetType();
 
    Double_t rot = sector->GetRotation();
    Double_t dx = sector->GetDx();
    Double_t dy = sector->GetDy();
    Double_t stereoB = sector->GetStereoB();
 
    //  Double_t z      = station->GetZ();
    Int_t stationNr = station->GetStationNr();
    Int_t sectorNr = sector->GetSectorNr();
 
    // Some auxiliary values
    Double_t sinrot = TMath::Sin(rot);
    Double_t cosrot = TMath::Cos(rot);
    Double_t tanstr = TMath::Tan(stereoB);
 
    // Calculate error matrix in sector system
    Double_t vX, vY, vXY;
    if (iType == 1) {
        vX = dx / TMath::Sqrt(12.);
        vY = dy / TMath::Sqrt(12.);
        vXY = 0.;
    } else if (iType == 2 || iType == 3) {
        vX = dx / TMath::Sqrt(12.);
        vY = dx / TMath::Sqrt(6.) / TMath::Abs(tanstr);
        vXY = -1. * dx * dx / 12. / tanstr;
    } else {
        cerr << "-E- " << fName << "::FindHits: Illegal sector type " << iType << endl;
        return 0;
    }
 
    // Transform variances into global c.s.
    Double_t wX = vX * vX * cosrot * cosrot - 2. * vXY * cosrot * sinrot + vY * vY * sinrot * sinrot;
    Double_t wY = vX * vX * sinrot * sinrot + 2. * vXY * cosrot * sinrot + vY * vY * cosrot * cosrot;
    Double_t wXY = (vX * vX - vY * vY) * cosrot * sinrot + vXY * (cosrot * cosrot - sinrot * sinrot);
    Double_t sigmaX = TMath::Sqrt(wX);
    Double_t sigmaY = TMath::Sqrt(wY);
 
    // Now perform the loop over active channels
    set<Int_t>::iterator it1;
    set<Int_t>::iterator it2;
 
    // ----- Type 1 : Pixel sector   ---------------------------------------
    if (iType == 1) {
        Fatal("FindHits", "Sorry, not implemented yet");
    }   // Pixel sensor
    // ---------------------------------------------------------------------
 
    // -----  Type 2: Strip sector OSU   -----------------------------------
    else if (iType == 2)
    {
        Fatal("FindHits", "Sorry, not implemented yet");
    }   // Strip OSU
    // ---------------------------------------------------------------------
 
    // -----  Type 3: Strip sector GSI   -----------------------------------
    else if (iType == 3)
    {
        Int_t iClusF = -1;
        Int_t iClusB = -1;
        Double_t chanF = -1;
        Double_t chanB = -1;
        Int_t nHits = fHits->GetEntriesFast();
        Double_t xHit;
        Double_t yHit;
        Double_t zHit;
        TVector3 pos, dpos;
        CbmStsCluster* clusterF = NULL;
        CbmStsCluster* clusterB = NULL;
        for (it1 = fSet.begin(); it1 != fSet.end(); it1++) {
            iClusF = (*it1);
            clusterF = (CbmStsCluster*)fClusters->At(iClusF);
            if (!clusterF) {
                cout << "-W- " << GetName() << "::FindHits: Invalid cluster index " << iClusF
                     << " in front set of sector " << sector->GetSectorNr() << ", station " << station->GetStationNr()
                     << endl;
                continue;
            }
            chanF = clusterF->GetMean();
            for (it2 = bSet.begin(); it2 != bSet.end(); it2++) {
                iClusB = (*it2);
                clusterB = (CbmStsCluster*)fClusters->At(iClusB);
                if (!clusterB) {
                    cout << "-W- " << GetName() << "::FindHits: Invalid cluster index " << iClusB
                         << " in back set of sector " << sector->GetSectorNr() << ", station "
                         << station->GetStationNr() << endl;
                    continue;
                }
                chanB = clusterB->GetMean();
 
                Int_t sensorDetId = sector->IntersectClusters(chanF, chanB, xHit, yHit, zHit);
 
                if (sensorDetId == -1)
                    continue;
 
                pos.SetXYZ(xHit, yHit, zHit);
                dpos.SetXYZ(sigmaX * clusterF->GetMeanError(), sigmaY * clusterB->GetMeanError(), 0.);
 
                Int_t statLayer = -1;
                for (Int_t istatL = station->GetNofZ(); istatL > 0; istatL--)
                    if (TMath::Abs(zHit - station->GetZ(istatL - 1)) < 0.00001) {
                        statLayer = istatL - 1;
                        break;
                    }
 
                if (statLayer == -1)
                    cout << "unknown layer for hit at z = " << zHit << endl;
 
                new ((*fHits)[nHits++])
                    CbmStsHit(sensorDetId, pos, dpos, wXY, iClusF, iClusB, (Int_t)chanF, (Int_t)chanB, statLayer);
 
                nNew++;
                if (fVerbose > 3)
                    cout << "New StsHit at (" << xHit << ", " << yHit << ", " << zHit << "), station " << stationNr
                         << ", sector " << sectorNr << ", channel " << chanF << " / " << chanB << endl;
 
            }   // back side strip loop
        }   // front side strip loop
 
    }   // strip GSI
    // ---------------------------------------------------------------------
 
    return nNew;
}
// -------------------------------------------------------------------------
 
// -----   Virtual method Finish   -----------------------------------------
void CbmStsRealFindHits::Finish()
{
    cout << endl;
    cout << "============================================================" << endl;
    cout << "===== " << fName << ": Run summary " << endl;
    cout << "===== " << endl;
    cout << "===== Number of hits                 : " << setw(8) << setprecision(2) << fNHits << endl;
    cout << "============================================================" << endl;
}
// -------------------------------------------------------------------------