CbmStsIdealDigitize.cxx

Go to the documentation of this file.

//* $Id: */
 
// -------------------------------------------------------------------------
// -----                    CbmStsDigitize source file                 -----
// -----                  Created 30/08/06  by V. Friese               -----
// -------------------------------------------------------------------------
 
// Includes from STS
#include "CbmStsIdealDigitize.h"
 
#include "CbmGeoStsPar.h"
#include "CbmStsDigi.h"
#include "CbmStsDigiMatch.h"
#include "CbmStsDigiPar.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsPoint.h"
#include "CbmStsSector.h"
#include "CbmStsSensor.h"
#include "CbmStsStation.h"
 
// Includes from base
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
 
// Includes from ROOT
#include "TClonesArray.h"
#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TMath.h"
#include "TObjArray.h"
 
#include <iomanip>
#include <iostream>
#include <map>
 
using std::cerr;
using std::cout;
using std::endl;
using std::fixed;
using std::left;
using std::map;
using std::pair;
using std::right;
using std::setprecision;
using std::setw;
 
// -----   Default constructor   ------------------------------------------
CbmStsIdealDigitize::CbmStsIdealDigitize()
    : FairTask("STS Digitizer", 1)
    , fGeoPar(NULL)
    , fDigiPar(NULL)
    , fDigiScheme(NULL)
    , fPoints(NULL)
    , fDigis(NULL)
    , fDigiMatches(NULL)
    , fNPoints(0)
    , fNFailed(0)
    , fNOutside(0)
    , fNMulti(0)
    , fNDigis(0)
    , fTimer()
    , fChannelMap()
{
    Reset();
}
// -------------------------------------------------------------------------
 
// -----   Standard constructor   ------------------------------------------
CbmStsIdealDigitize::CbmStsIdealDigitize(Int_t iVerbose)
    : FairTask("STSDigitize", iVerbose)
    , fGeoPar(NULL)
    , fDigiPar(NULL)
    , fDigiScheme(NULL)
    , fPoints(NULL)
    , fDigis(NULL)
    , fDigiMatches(NULL)
    , fNPoints(0)
    , fNFailed(0)
    , fNOutside(0)
    , fNMulti(0)
    , fNDigis(0)
    , fTimer()
    , fChannelMap()
{
    Reset();
}
// -------------------------------------------------------------------------
 
// -----   Constructor with name   -----------------------------------------
CbmStsIdealDigitize::CbmStsIdealDigitize(const char* name, Int_t iVerbose)
    : FairTask(name, iVerbose)
    , fGeoPar(NULL)
    , fDigiPar(NULL)
    , fDigiScheme(NULL)
    , fPoints(NULL)
    , fDigis(NULL)
    , fDigiMatches(NULL)
    , fNPoints(0)
    , fNFailed(0)
    , fNOutside(0)
    , fNMulti(0)
    , fNDigis(0)
    , fTimer()
    , fChannelMap()
{
    fDigiScheme = new CbmStsDigiScheme();
    Reset();
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
CbmStsIdealDigitize::~CbmStsIdealDigitize()
{
    if (fGeoPar)
        delete fGeoPar;
    if (fDigiPar)
        delete fDigiPar;
    if (fDigis) {
        fDigis->Delete();
        delete fDigis;
    }
    if (fDigiMatches) {
        fDigiMatches->Delete();
        delete fDigiMatches;
    }
    if (fDigiScheme)
        delete fDigiScheme;
    Reset();
}
// -------------------------------------------------------------------------
 
// -----   Public method Exec   --------------------------------------------
void CbmStsIdealDigitize::Exec(Option_t* opt)
{
 
    // Reset all eventwise counters
    fTimer.Start();
    Reset();
 
    // Verbose screen output
    if (fVerbose > 2) {
        cout << endl << "-I- " << fName << ": executing event" << endl;
        cout << "----------------------------------------------" << endl;
    }
 
    // Loop over all StsPoints
    if (!fPoints) {
        cerr << "-W- " << fName << "::Exec: No input array (STSPoint) " << endl;
        cout << "- " << fName << endl;
        return;
    }
 
    for (Int_t iPoint = 0; iPoint < fPoints->GetEntriesFast(); iPoint++) {
        CbmStsPoint* point = (CbmStsPoint*)fPoints->At(iPoint);
        fNPoints++;
 
        Double_t xIn = point->GetXIn();
        Double_t yIn = point->GetYIn();
        Double_t zIn = point->GetZIn();
 
        gGeoManager->FindNode(xIn, yIn, zIn);
        TGeoNode* curNode = gGeoManager->GetCurrentNode();
 
        // AZ CbmStsSensor* sensor = fDigiScheme->GetSensorByName(curNode->GetName());
        // AZ
        CbmStsSensor* sensor = NULL;
        if (fDigiScheme->IsNewGeometry()) {
            TString curPath = fDigiScheme->GetCurrentPath();
            sensor = fDigiScheme->GetSensorByName(curPath);
        } else
            sensor = fDigiScheme->GetSensorByName(curNode->GetName());
        // AZ
 
        if (!sensor)
            LOG(fatal) << "No sensor in DigiScheme with name: " << curNode->GetName();
 
        Int_t stationNr = sensor->GetStationNr();
        Int_t sectorNr = sensor->GetSectorNr();
 
        // Take point coordinates in the midplane of the sensor
        if (!point->IsUsable())
            continue;
        Double_t xpt = point->GetX(sensor->GetZ0());
        Double_t ypt = point->GetY(sensor->GetZ0());
 
        Int_t iChan = -1;
        Int_t channelF = -1;
        Int_t channelB = -1;
 
        iChan = sensor->GetChannel(xpt, ypt, 0);
 
        // point outside active sensor area
        if (iChan < 0) {
            if (fVerbose) {
                Float_t tempX = sensor->GetX0();
                Float_t tempY = sensor->GetY0();
                Float_t rot = sensor->GetRotation();
                Float_t tempLX = sensor->GetLx() / 2. * TMath::Cos(rot) + sensor->GetLy() / 2. * TMath::Sin(rot);
                Float_t tempLY = sensor->GetLy() / 2. * TMath::Cos(rot) + sensor->GetLx() / 2. * TMath::Sin(rot);
 
                cout.precision(8);
                cout << fName << ". Warning: Point " << xpt << ", " << ypt << ", " << sensor->GetZ0()
                     << " not inside sector." << endl;
            }
            continue;
        }
 
        channelF = iChan;
 
        Int_t sectorDetId = 2 << 24 | stationNr << 16 | sectorNr << 4;
 
        // Treat front side
        pair<Int_t, Int_t> a(sectorDetId, channelF);
        if (fChannelMap.find(a) == fChannelMap.end()) {
            // Channel not yet active. Create new Digi and Match.
            new ((*fDigis)[fNDigis]) CbmStsDigi(stationNr, sectorNr, 0, channelF, 0, 0);
            new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(iPoint);
            fChannelMap[a] = fNDigis;
            fNDigis++;
        } else {
            // Channel already active. Update DigiMatch.
            Int_t iDigi = fChannelMap[a];
            CbmStsDigiMatch* match = dynamic_cast<CbmStsDigiMatch*>(fDigiMatches->At(iDigi));
            if (match) {
                match->AddPoint(iPoint);
                fNMulti++;
            }
        }
 
        // Treat back side (strip sensors only)
        if (sensor->GetType() < 2 || sensor->GetType() > 3)
            continue;
 
        iChan = sensor->GetChannel(xpt, ypt, 1);
        if (iChan < 0) {
            cerr << "-W- " << fName << "::Exec: No back side channel "
                 << " for StsPoint " << iPoint << ", station " << stationNr << ", sector " << sectorNr
                 << ", front side channel " << channelF << endl;
            continue;
        }
        sectorDetId = sectorDetId | (1 << 0);   // for back side channel
        channelB = iChan;
        pair<Int_t, Int_t> b(sectorDetId, channelB);
        if (fChannelMap.find(b) == fChannelMap.end()) {
            // Channel not yet active. Create new Digi and Match.
            new ((*fDigis)[fNDigis]) CbmStsDigi(stationNr, sectorNr, 1, channelB, 0, 0);
            new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(iPoint);
            fChannelMap[b] = fNDigis;
            fNDigis++;
        } else {
            // Channel already active. Update DigiMatch.
            Int_t iDigi = fChannelMap[b];
            CbmStsDigiMatch* match = dynamic_cast<CbmStsDigiMatch*>(fDigiMatches->At(iDigi));
            if (match) {
                match->AddPoint(iPoint);
                fNMulti++;
            }
        }
 
        // Verbose screen output
        if (fVerbose > 2) {
            cout.precision(6);
            cout << "StsPoint " << iPoint << ", station " << stationNr << ", (" << xpt << ", " << ypt << ") cm, sector "
                 << sectorNr << ", front " << channelF << ", back " << channelB << endl;
        }
 
        // Not found in any sector?
        if (sectorNr == -1)
            fNOutside++;
 
    }   // StsPoint loop
 
    // Screen output
    fTimer.Stop();
    if (fVerbose) {
        cout << "+ ";
        cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8) << fixed << right << fTimer.RealTime()
             << " s, points " << fNPoints << ", failed " << fNFailed << ", outside " << fNOutside << ", multihits "
             << fNMulti << ", digis " << fNDigis << endl;
    }
}
// -------------------------------------------------------------------------
 
// -----   Private method SetParContainers   -------------------------------
void CbmStsIdealDigitize::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 CbmStsIdealDigitize::Init()
{
 
    // Get input array
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman)
        Fatal("Init", "No FairRootManager");
    fPoints = (TClonesArray*)ioman->GetObject("StsPoint");
 
    // Register output array StsDigi
    fDigis = new TClonesArray("CbmStsDigi", 1000);
    ioman->Register("StsDigi", "Digital response in STS", fDigis, kTRUE);
 
    // Register output array StsDigiMatches
    fDigiMatches = new TClonesArray("CbmStsDigiMatch", 1000);
    ioman->Register("StsDigiMatch", "Digi Match in STS", fDigiMatches, kTRUE);
 
    // 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 CbmStsIdealDigitize::ReInit()
{
 
    // Clear digitisation scheme
    fDigiScheme->Clear();
 
    // Build new digitisation scheme
    if (fDigiScheme->Init(fGeoPar, fDigiPar))
        return kSUCCESS;
 
    return kERROR;
}
// -------------------------------------------------------------------------
 
// -----   Private method Reset   ------------------------------------------
void CbmStsIdealDigitize::Reset()
{
    fNPoints = fNFailed = fNOutside = fNMulti = fNDigis = 0;
    fChannelMap.clear();
    //  if ( fDigis ) fDigis->Clear();
    //  if ( fDigiMatches ) fDigiMatches->Clear();
    if (fDigis)
        fDigis->Delete();
    if (fDigiMatches)
        fDigiMatches->Delete();
}
// -------------------------------------------------------------------------