CbmStsRealDigitize.cxx

Go to the documentation of this file.

//* $Id: */
 
// -------------------------------------------------------------------------
// -----                    CbmStsRealDigitize source file             -----
// -----                  Created 08/07/2008  by R. Karabowicz         -----
// -------------------------------------------------------------------------
 
// Includes from ROOT
#include "TClonesArray.h"
#include "TF1.h"
#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TMath.h"
#include "TObjArray.h"
#include "TRandom3.h"
 
// Includes from base
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
 
// Includes from STS
#include "CbmGeoStsPar.h"
#include "CbmStsDigi.h"
#include "CbmStsDigiMatch.h"
#include "CbmStsDigiPar.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsPoint.h"
#include "CbmStsRealDigitize.h"
#include "CbmStsSector.h"
#include "CbmStsSensor.h"
#include "CbmStsStation.h"
 
#include <iomanip>
#include <iostream>
#include <map>
 
using std::cerr;
using std::cout;
using std::endl;
using std::fixed;
using std::flush;
using std::ios_base;
using std::left;
using std::map;
using std::pair;
using std::right;
using std::set;
using std::setprecision;
using std::setw;
 
// -----   Default constructor   ------------------------------------------
CbmStsRealDigitize::CbmStsRealDigitize()
    : FairTask("STS Digitizer", 1)
{
    fGeoPar = NULL;
    fDigiPar = NULL;
    fPoints = NULL;
    fDigis = NULL;
    fDigiMatches = NULL;
    fRealistic = kFALSE;
    fDigiScheme = new CbmStsDigiScheme();
    Reset();
 
    fStep = 0.001;
 
    fFThreshold = 0.001;
    fBThreshold = 0.001;
    fFNoiseWidth = 0.0;
    fBNoiseWidth = 0.0;
 
    fFNofBits = 20;
    fBNofBits = 20;
    fFMinStep = 0.01;
    fBMinStep = 0.01;
 
    fNEvents = 0.;
}
// -------------------------------------------------------------------------
 
// -----   Standard constructor   ------------------------------------------
CbmStsRealDigitize::CbmStsRealDigitize(Int_t iVerbose)
    : FairTask("STSDigitize", iVerbose)
{
    fGeoPar = NULL;
    fDigiPar = NULL;
    fPoints = NULL;
    fDigis = NULL;
    fDigiMatches = NULL;
    fRealistic = kFALSE;
    fDigiScheme = new CbmStsDigiScheme();
    Reset();
 
    fStep = 0.001;
 
    fFThreshold = 0.001;
    fBThreshold = 0.001;
    fFNoiseWidth = 0.0;
    fBNoiseWidth = 0.0;
 
    fFNofBits = 20;
    fBNofBits = 20;
    fFMinStep = 0.01;
    fBMinStep = 0.01;
 
    fNEvents = 0.;
}
// -------------------------------------------------------------------------
 
// -----   Constructor with name   -----------------------------------------
CbmStsRealDigitize::CbmStsRealDigitize(const char* name, Int_t iVerbose)
    : FairTask(name, iVerbose)
{
    fGeoPar = NULL;
    fDigiPar = NULL;
    fPoints = NULL;
    fDigis = NULL;
    fDigiMatches = NULL;
    fRealistic = kFALSE;
    fDigiScheme = new CbmStsDigiScheme();
    Reset();
 
    fStep = 0.001;
 
    fFThreshold = 0.001;
    fBThreshold = 0.001;
    fFNoiseWidth = 0.0;
    fBNoiseWidth = 0.0;
 
    fFNofBits = 20;
    fBNofBits = 20;
    fFMinStep = 0.01;
    fBMinStep = 0.01;
 
    fNEvents = 0.;
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
CbmStsRealDigitize::~CbmStsRealDigitize()
{
    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 CbmStsRealDigitize::Exec(Option_t* opt)
{
 
    //   cout << "threshold = " << fFThreshold << endl;
 
    // Reset all eventwise counters
    fTimer.Start();
    Reset();
    TF1* digiGausDist = new TF1("digiGausDist", "gaus", -5., 5.);
 
    // Verbose screen output
    if (fVerbose > 2) {
        cout << endl << "-I- " << fName << ": executing event" << endl;
        cout << "----------------------------------------------" << endl;
    }
 
    Double_t nPoints = 0.;
    Double_t nOutside = 0.;
    Double_t nDigisF = 0.;
    Double_t nDigisB = 0.;
 
    // Loop over all StsPoints
    if (!fPoints) {
        cerr << "-W- " << fName << "::Exec: No input array (STSPoint) " << endl;
        cout << "- " << fName << endl;
        return;
    }
 
    map<CbmStsSensor*, set<Int_t>>::iterator mapIt;
    for (mapIt = fPointMap.begin(); mapIt != fPointMap.end(); mapIt++)
        ((*mapIt).second).clear();
 
    for (Int_t iPoint = 0; iPoint < fPoints->GetEntriesFast(); iPoint++) {
        CbmStsPoint* point = (CbmStsPoint*)fPoints->At(iPoint);
 
        Double_t xin = point->GetXIn();
        Double_t yin = point->GetYIn();
        Double_t zin = point->GetZIn();
        gGeoManager->FindNode(xin, yin, zin);
        TGeoNode* curNode = gGeoManager->GetCurrentNode();
 
        CbmStsSensor* sensor = fDigiScheme->GetSensorByName(curNode->GetName());
 
        if (fPointMap.find(sensor) == fPointMap.end()) {
            cerr << "-E- " << fName << "::Exec:: sensor " << curNode->GetName() << " not found in digi scheme!" << endl;
            continue;
        }
        fPointMap[sensor].insert(iPoint);
    }
 
    for (Int_t iStation = fDigiScheme->GetNStations(); iStation > 0;) {
        CbmStsStation* station = fDigiScheme->GetStation(--iStation);
        for (Int_t iSector = station->GetNSectors(); iSector > 0;) {
            CbmStsSector* sector = station->GetSector(--iSector);
 
            map<Int_t, set<Int_t>>::iterator mapCh;
 
            for (mapCh = fFChannelPointsMap.begin(); mapCh != fFChannelPointsMap.end(); mapCh++)
                ((*mapCh).second).clear();
            for (mapCh = fBChannelPointsMap.begin(); mapCh != fBChannelPointsMap.end(); mapCh++)
                ((*mapCh).second).clear();
 
            // simulating detector+cables+electronics noise
            // should be more sophisticated...
            // the question is: sectorwise or sensorwise???
            Int_t nChannels = sector->GetNChannelsFront();
            for (Int_t iChannel = nChannels; iChannel > 0;) {
                //  fStripSignalF[--iChannel] = fGen->Landau(.1,.02);
                //  fStripSignalB[  iChannel] = fGen->Landau(.1,.02);
                //  fStripSignalF[--iChannel] = 0.;
                //  fStripSignalB[  iChannel] = 0.;
                fStripSignalF[--iChannel] = TMath::Abs(gRandom->Gaus(0., fFNoiseWidth));
                fStripSignalB[iChannel] = TMath::Abs(gRandom->Gaus(0., fBNoiseWidth));
            }
 
            for (Int_t iSensor = sector->GetNSensors(); iSensor > 0;) {
                CbmStsSensor* sensor = sector->GetSensor(--iSensor);
 
                ProduceHitResponse(sensor);
            }
 
            Int_t stationNr = sector->GetStationNr();
            Int_t sectorNr = sector->GetSectorNr();
            Int_t sectorDetId = sector->GetDetectorId();
 
            for (Int_t ifstr = 0; ifstr < nChannels; ifstr++) {
                if (fStripSignalF[ifstr] < fFThreshold)
                    continue;
                //  cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr
                //       << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;
                Int_t digiFSignal = 1 + (Int_t)((fStripSignalF[ifstr] - fFThreshold) / fFMinStep);
                if (digiFSignal >= fFNofSteps)
                    digiFSignal = fFNofSteps - 1;
                new ((*fDigis)[fNDigis]) CbmStsDigi(0, stationNr, sectorNr, 0, ifstr, digiFSignal, 0);
                //                        fStripSignalF[ifstr], 0);
                set<Int_t>::iterator it1;
                set<Int_t> chPnt = fFChannelPointsMap[ifstr];
                Int_t pnt;
                CbmStsDigiMatch* match;
                if (chPnt.size() == 0) {
                    //    cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr
                    //         << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;
                    new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(-666);
                } else {
                    for (it1 = chPnt.begin(); it1 != chPnt.end(); it1++) {
                        pnt = (*it1);
                        if (it1 == chPnt.begin())
                            match = new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(pnt);
                        else {
                            match->AddPoint(pnt);
                            fNMulti++;
                        }
                    }
                }
                fNDigis++;
            }
 
            for (Int_t ibstr = 0; ibstr < nChannels; ibstr++) {
                if (fStripSignalB[ibstr] < fBThreshold)
                    continue;
                Int_t digiBSignal = 1 + (Int_t)((fStripSignalB[ibstr] - fBThreshold) / fBMinStep);
                if (digiBSignal >= fBNofSteps)
                    digiBSignal = fBNofSteps - 1;
                new ((*fDigis)[fNDigis]) CbmStsDigi(0, stationNr, sectorNr, 1, ibstr, digiBSignal, 0);
                //                        fStripSignalB[ibstr], 0);
                set<Int_t>::iterator it1;
                set<Int_t> chPnt = fBChannelPointsMap[ibstr];
                Int_t pnt;
                CbmStsDigiMatch* match;
                if (chPnt.size() == 0) {
                    //    cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr
                    //         << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;
                    new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(-666);
                } else {
                    for (it1 = chPnt.begin(); it1 != chPnt.end(); it1++) {
                        pnt = (*it1);
                        if (it1 == chPnt.begin())
                            match = new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(pnt);
                        else {
                            match->AddPoint(pnt);
                            fNMulti++;
                        }
                    }
                }
                fNDigis++;
            }
        }
    }
 
    fTimer.Stop();
    cout << "+ " << flush;
    cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8) << fixed << right << fTimer.RealTime()
         << " s, digis " << nDigisF << " / " << nDigisB << endl;
 
    fNEvents += 1.;
    fNPoints += Double_t(nPoints);
    fNOutside += Double_t(nOutside);
    fNDigisFront += Double_t(nDigisF);
    fNDigisBack += Double_t(nDigisB);
    fTime += fTimer.RealTime();
}
// -------------------------------------------------------------------------
 
// -----   Private method ProduceHitResponse   --------------------------------
void CbmStsRealDigitize::ProduceHitResponse(CbmStsSensor* sensor)
{
    set<Int_t> pSet;
    if (fPointMap.find(sensor) == fPointMap.end()) {
        cerr << "-E- " << fName << "::Exec:: sensor"
             << " not found in digi scheme!" << endl;
        return;
    }
    pSet = fPointMap[sensor];
 
    Int_t iPoint = -1;
    CbmStsPoint* point = NULL;
 
    set<Int_t>::iterator it1;
 
    for (it1 = pSet.begin(); it1 != pSet.end(); it1++) {
        iPoint = (*it1);
        point = (CbmStsPoint*)fPoints->At(iPoint);
 
        Double_t xin = point->GetXIn();
        Double_t yin = point->GetYIn();
        Double_t zin = point->GetZIn();
 
        Double_t xvec = point->GetXOut() - xin;
        Double_t yvec = point->GetYOut() - yin;
        Double_t zvec = point->GetZOut() - zin;
 
        Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec * xvec + yvec * yvec + zvec * zvec) / fStep + 1);
 
        Double_t stepEL = fEnergyLossToSignal * point->GetEnergyLoss() / (nofSteps + 1);
 
        xvec = xvec / ((Double_t)nofSteps);
        yvec = yvec / ((Double_t)nofSteps);
        zvec = zvec / ((Double_t)nofSteps);
 
        for (Int_t istep = 0; istep <= nofSteps; istep++) {
            //      Float_t iFChan = sensor->GetChannelPlus(xin, yin, 0);
            Int_t iIChan = sensor->GetFrontChannel(xin, yin, zin);
            //(Int_t)iFChan;
 
            if (iIChan != -1) {
                fStripSignalF[iIChan] += stepEL;
                fFChannelPointsMap[iIChan].insert(iPoint);
            }
 
            //      iFChan = sensor->GetChannelPlus(xin, yin, 1);
            //      iIChan = (Int_t)iFChan;
            iIChan = sensor->GetBackChannel(xin, yin, zin);
 
            if (iIChan != -1) {
                fStripSignalB[iIChan] += stepEL;
                fBChannelPointsMap[iIChan].insert(iPoint);
            }
 
            xin += xvec;
            yin += yvec;
            zin += zvec;
        }
    }
}
// -------------------------------------------------------------------------
 
// -----   Private method FindFiredStrips   --------------------------------
void CbmStsRealDigitize::FindFiredStrips(CbmStsPoint* pnt,
                                         Int_t& nofStr,
                                         Int_t*& strips,
                                         Double_t*& signals,
                                         Int_t side)
{
 
    nofStr = 0;
 
    Double_t xin = pnt->GetXIn();
    Double_t yin = pnt->GetYIn();
    Double_t zin = pnt->GetZIn();
 
    gGeoManager->FindNode(xin, yin, zin);
    TGeoNode* curNode = gGeoManager->GetCurrentNode();
 
    CbmStsSensor* sensor = fDigiScheme->GetSensorByName(curNode->GetName());
 
    Double_t xvec = pnt->GetXOut() - xin;
    Double_t yvec = pnt->GetYOut() - yin;
    Double_t zvec = pnt->GetZOut() - zin;
 
    Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec * xvec + yvec * yvec + zvec * zvec) / fStep + 1);
 
    Double_t stepEL = fEnergyLossToSignal * pnt->GetEnergyLoss() / (nofSteps + 1);
 
    xvec = xvec / ((Double_t)nofSteps);
    yvec = yvec / ((Double_t)nofSteps);
    zvec = zvec / ((Double_t)nofSteps);
 
    for (Int_t istep = 0; istep <= nofSteps; istep++) {
        Float_t iFChan = sensor->GetChannelPlus(xin, yin, side);
        Int_t iIChan = (Int_t)iFChan;
 
        xin += xvec;
        yin += yvec;
        zin += zvec;
 
        if (iIChan == -1)
            continue;
 
        for (Int_t ifstr = 0; ifstr < nofStr; ifstr++) {
            if (strips[ifstr] == iIChan) {
                signals[ifstr] += stepEL;
                iIChan = -1;
                break;
            }
        }
        if (iIChan == -1)
            continue;
 
        strips[nofStr] = iIChan;
        signals[nofStr] = stepEL;
 
        nofStr++;
    }
}
// -------------------------------------------------------------------------
 
// -----   Private method SetParContainers   -------------------------------
void CbmStsRealDigitize::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 CbmStsRealDigitize::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);
 
    //   fGen = new TRandom3();
    //   time_t curtime;
    //   time(&curtime);
    //   fGen->SetSeed(curtime);
 
    fStripSignalF = new Double_t[2000];
    fStripSignalB = new Double_t[2000];
 
    fEnergyLossToSignal = 200000.;
 
    fFNofSteps = (Int_t)TMath::Power(2, (Double_t)fFNofBits);
    fBNofSteps = (Int_t)TMath::Power(2, (Double_t)fBNofBits);
 
    // Build digitisation scheme
    if (fDigiScheme->Init(fGeoPar, fDigiPar)) {
        MakeSets();
 
        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 CbmStsRealDigitize::ReInit()
{
 
    // Clear digitisation scheme
    fDigiScheme->Clear();
 
    // Build new digitisation scheme
    if (fDigiScheme->Init(fGeoPar, fDigiPar)) {
        MakeSets();
        return kSUCCESS;
    }
 
    return kERROR;
}
// -------------------------------------------------------------------------
 
// -----   Private method MakeSets   ---------------------------------------
void CbmStsRealDigitize::MakeSets()
{
 
    fPointMap.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);
            Int_t nSensors = sector->GetNSensors();
            for (Int_t iSensor = 0; iSensor < nSensors; iSensor++) {
                CbmStsSensor* sensor = sector->GetSensor(iSensor);
                set<Int_t> a;
                fPointMap[sensor] = a;
            }
        }
    }
    fFChannelPointsMap.clear();
    fBChannelPointsMap.clear();
    for (Int_t ichan = 2000; ichan > 0;) {
        set<Int_t> a;
        fFChannelPointsMap[--ichan] = a;
        set<Int_t> b;
        fBChannelPointsMap[ichan] = b;
    }
}
// -------------------------------------------------------------------------
 
// -----   Private method Reset   ------------------------------------------
void CbmStsRealDigitize::Reset()
{
    //  fNPoints = fNOutside = fNDigisFront = fNDigisBack = fTime = 0.;
    fNDigis = fNMulti = 0;
    fFChannelPointsMap.clear();
    fBChannelPointsMap.clear();
    //  if ( fDigis ) fDigis->Clear();
    //  if ( fDigiMatches ) fDigiMatches->Clear();
    if (fDigis)
        fDigis->Delete();
    if (fDigiMatches)
        fDigiMatches->Delete();
}
// -------------------------------------------------------------------------
 
// -----   Virtual method Finish   -----------------------------------------
void CbmStsRealDigitize::Finish()
{
    cout << endl;
    cout << "============================================================" << endl;
    cout << "===== " << fName << ": Run summary " << endl;
    cout << "===== " << endl;
    cout << "===== Events processed          : " << setw(8) << fNEvents << endl;
    cout.setf(ios_base::fixed, ios_base::floatfield);
    cout << "===== Real time per event       : " << setw(8) << setprecision(4) << fTime / fNEvents << " s" << endl;
    cout << "===== StsPoints per event       : " << setw(8) << setprecision(2) << fNPoints / fNEvents << endl;
    cout << "===== Outside hits per event    : " << setw(8) << setprecision(2) << fNOutside / fNEvents << " = "
         << setw(6) << setprecision(2) << fNOutside / fNPoints * 100. << " %" << endl;
    cout << "===== Front digis per point     : " << setw(8) << setprecision(2) << fNDigisFront / (fNPoints - fNOutside)
         << endl;
    cout << "===== Back digis per point      : " << setw(8) << setprecision(2) << fNDigisBack / (fNPoints - fNOutside)
         << endl;
    cout << "============================================================" << endl;
}
// -------------------------------------------------------------------------