BmnSiliconDigitizer.cxx

Go to the documentation of this file.

#include "BmnSiliconDigitizer.h"
 
#include "BmnSiliconPoint.h"
#include "CbmMCTrack.h"
#include "FairRootManager.h"
#include "TSystem.h"
 
#include <TRandom.h>   //AZ-240823
 
static Float_t workTime = 0.0;
static int entries = 0;
 
BmnSiliconDigitizer::BmnSiliconDigitizer()
    : fOnlyPrimary(kFALSE)
    , fStripMatching(kTRUE)
    , fUseRealEffects(kFALSE)
{
 
    fInputBranchName = "SiliconPoint";
    fOutputDigitsBranchName = "BmnSiliconDigit";
    fOutputDigitMatchesBranchName = "BmnSiliconDigitMatch";
 
    fVerbose = 1;
 
    fCurrentConfig = BmnSiliconConfiguration::None;
    StationSet = nullptr;
    fField = nullptr;
    GlobalModTransfSet = nullptr;
    fBmnSiliconDigitsArray = nullptr;
    fBmnSiliconDigitMatchesArray = nullptr;
 
    fSigma = -1;    // AZ-240823
    fThresh = 12;   // AZ-240823 - initial value
}
 
BmnSiliconDigitizer::BmnSiliconDigitizer(UInt_t run_period)
    : BmnSiliconDigitizer()
{
    switch (run_period) {
        case 6:   // BM@N RUN-6
            fCurrentConfig = BmnSiliconConfiguration::RunSpring2017;
            break;
        case 7:   // BM@N RUN-7
            fCurrentConfig = BmnSiliconConfiguration::RunSpring2018;
            break;
        case 8:   // BM@N RUN-8
            fCurrentConfig = BmnSiliconConfiguration::Run8;
            break;
        case 9:   // BM@N RUN-9
            fCurrentConfig = BmnSiliconConfiguration::Run9;
            break;
 
        default:
            fCurrentConfig = BmnSiliconConfiguration::None;
    }
}
 
BmnSiliconDigitizer::~BmnSiliconDigitizer()
{
    if (StationSet) {
        delete StationSet;
    }
 
    if (GlobalModTransfSet) {
        delete GlobalModTransfSet;
    }
    if (fBmnSiliconDigitsArray) {
        fBmnSiliconDigitsArray->Delete();
        delete fBmnSiliconDigitsArray;
    }
    if (fBmnSiliconDigitMatchesArray) {
        fBmnSiliconDigitMatchesArray->Delete();
        delete fBmnSiliconDigitMatchesArray;
    }
}
 
InitStatus BmnSiliconDigitizer::Init()
{
    if (fVerbose)
        cout << "\nBmnSiliconDigitizer::Init()\n ";
 
    if (fVerbose && fOnlyPrimary)
        cout << "  Only primary particles are processed!!! " << endl;
 
    if (fVerbose && fStripMatching)
        cout << "  Strip Matching is activated!!! " << endl;
    else
        cout << "  Strip matching is deactivated!!! " << endl;
 
    FairRootManager* ioman = FairRootManager::Instance();
 
    fBmnSiliconPointsArray = (TClonesArray*)ioman->GetObject(fInputBranchName);
    if (!fBmnSiliconPointsArray) {
        cout << "BmnSiliconDigitizer::Init(): branch SiliconPoint not found! Task "
                "will be deactivated"
             << endl;
        SetActive(kFALSE);
        return kERROR;
    }
 
    fMCTracksArray = (TClonesArray*)ioman->GetObject("MCTrack");
 
    fBmnSiliconDigitsArray = new TClonesArray(fOutputDigitsBranchName);
    ioman->Register(fOutputDigitsBranchName, "SILICON", fBmnSiliconDigitsArray, kTRUE);
 
    if (fStripMatching) {
        fBmnSiliconDigitMatchesArray = new TClonesArray("BmnMatch");
        ioman->Register(fOutputDigitMatchesBranchName, "SILICON", fBmnSiliconDigitMatchesArray, kTRUE);
    }
 
    TString gPathSiliconConfig = gSystem->Getenv("VMCWORKDIR");
    gPathSiliconConfig += "/parameters/silicon/XMLConfigs/";
 
    // Create SILICON detector -------------------------------------------------
    switch (fCurrentConfig) {
        case BmnSiliconConfiguration::RunSpring2017:
            XMLConfigFile = gPathSiliconConfig + "SiliconRunSpring2017.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : RunSpring2017"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::RunSpring2018:
            XMLConfigFile = gPathSiliconConfig + "SiliconRunSpring2018.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : RunSpring2018"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::RunSRCSpring2018:
            XMLConfigFile = gPathSiliconConfig + "SiliconRunSRCSpring2018.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : RunSRCSpring2018"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::Run8_3stations:
            XMLConfigFile = gPathSiliconConfig + "SiliconRun8_3stations.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : SiliconRun8_3stations"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::Run8_4stations:
            XMLConfigFile = gPathSiliconConfig + "SiliconRun8_4stations.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : SiliconRun8_4stations"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::Run8_5stations:
            XMLConfigFile = gPathSiliconConfig + "SiliconRun8_5stations.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : SiliconRun8_5stations"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::Run8_mods_6_10_14_18:
            XMLConfigFile = gPathSiliconConfig + "SiliconRun8_mods_6_10_14_18.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : SiliconRun8_mods_6_10_14_18"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::Run8_mods_10_14rot_18:
            XMLConfigFile = gPathSiliconConfig + "SiliconRun8_mods_10_14rot_18.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : SiliconRun8_mods_10_14rot_18"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::Run8:
            XMLConfigFile = gPathSiliconConfig + "SiliconRun8.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : Run8"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::Run9:
            XMLConfigFile = gPathSiliconConfig + "SiliconRun9.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : Run9"
                     << "\n";
            break;
 
        case BmnSiliconConfiguration::Run_5stations_2026:
            XMLConfigFile = gPathSiliconConfig + "Silicon_5stations_2026.xml";
            if (fVerbose)
                cout << "   Current SILICON Configuration : 5stations_2026"
                     << "\n";
            break;
 
        default:;
    }
 
    if (!gSystem->AccessPathName(XMLConfigFile)) {   // returns FALSE if the file exists
        if (fVerbose)
            cout << "   XMLConfigFile : " << XMLConfigFile << "\n";
        StationSet = new BmnSiliconStationSet(XMLConfigFile);
        GlobalModTransfSet = new BmnSiliconTransform();
        GlobalModTransfSet->LoadFromXMLFile(XMLConfigFile);
    } else {
        Fatal("BmnSiliconDigitizer::Init()", " !!! Current configuration is not set !!! ");
    }
    //--------------------------------------------------------------------------
 
    if (fVerbose)
        cout << "BmnSiliconDigitizer::Init() finished\n\n";
    return kSUCCESS;
}
 
void BmnSiliconDigitizer::Exec(Option_t* opt)
{
    clock_t tStart = clock();
    fBmnSiliconDigitsArray->Delete();
 
    if (fStripMatching) {
        fBmnSiliconDigitMatchesArray->Delete();
    }
 
    if (!fBmnSiliconPointsArray) {
        Error("BmnSiliconDigitizer::Exec()", " !!! Unknown branch name !!! ");
        return;
    }
 
    if (fVerbose) {
        // cout << " BmnSiliconDigitizer::Exec(), Number of BmnSiliconPoints = " <<
        // fBmnSiliconPointsArray->GetEntriesFast() << "\n";
    }
 
    ProcessMCPoints();
 
    // if (fVerbose) cout << " BmnSiliconDigitizer::Exec() finished\n\n";
    entries++;
    clock_t tFinish = clock();
    workTime += ((Float_t)(tFinish - tStart)) / CLOCKS_PER_SEC;
}
 
void BmnSiliconDigitizer::ProcessMCPoints()
{
 
    FairMCPoint* SiliconPoint;
    Int_t NNotPrimaries = 0;
 
    for (Int_t ipoint = 0; ipoint < fBmnSiliconPointsArray->GetEntriesFast(); ipoint++) {
        SiliconPoint = (FairMCPoint*)fBmnSiliconPointsArray->At(ipoint);
 
        if (fOnlyPrimary) {
            CbmMCTrack* track = (CbmMCTrack*)fMCTracksArray->UncheckedAt(SiliconPoint->GetTrackID());
            if (!track)
                continue;
            if (track->GetMotherId() != -1) {
                NNotPrimaries++;
                continue;
            }
        }
 
        Double_t x = -SiliconPoint->GetX();   // invert because in current geometry +x is left, -x is right
        Double_t y = SiliconPoint->GetY();
        Double_t z = SiliconPoint->GetZ();
 
        Double_t px = -SiliconPoint->GetPx();   // invert because in current geometry +x is left, -x is right
        Double_t py = SiliconPoint->GetPy();
        Double_t pz = SiliconPoint->GetPz();
 
        Double_t dEloss = SiliconPoint->GetEnergyLoss() * 1e6;   // in keV
        Int_t refId = ipoint;
 
        // Information from MC-points
        Int_t mc_station_num = ((BmnSiliconPoint*)SiliconPoint)->GetStation();
        Int_t mc_module_num = ((BmnSiliconPoint*)SiliconPoint)->GetModule();
 
        // test output
        // cout << "mc_station_num = " << mc_station_num << "\n";
        // cout << "mc_module_num = " << mc_module_num << "\n";
 
        // Apply inverse global module transforms to MC-point (in right-handed system)
        if (GlobalModTransfSet && mc_station_num < StationSet->GetNStations()) {
            if (mc_module_num < StationSet->GetStation(mc_station_num)->GetNModules()) {
                Plane3D::Point point_transformed =
                    GlobalModTransfSet->ApplyInverseTransforms(Plane3D::Point(-x, y, z), mc_station_num, mc_module_num);
                Plane3D::Point direction_transformed = GlobalModTransfSet->ApplyInverseTransforms(
                    Plane3D::Point(-(px + x), (py + y), (pz + z)), mc_station_num, mc_module_num);
                x = -point_transformed.X();   // convert to detector system (x-inverted)
                y = point_transformed.Y();
                z = point_transformed.Z();
 
                px = -(direction_transformed.X() - point_transformed.X());   // convert to detector system (x-inverted)
                py = direction_transformed.Y() - point_transformed.Y();
                pz = direction_transformed.Z() - point_transformed.Z();
            }
        }
 
        // AZ-191223 - gain fluctuations
        Double_t gain = 1;
        // gain += gRandom->Gaus(0, 0.5);
        if (gain < 0)
            gain = 0;
        dEloss *= gain;
        //
        StationSet->AddPointToDetector(x, y, z, px, py, pz, dEloss, refId);
    }
 
    // Int_t NAddedPoints = StationSet->CountNAddedToDetectorPoints();
    // if (fVerbose && fOnlyPrimary) cout << "   Number of not primaries points : " << NNotPrimaries << "\n";
    // if (fVerbose) cout << "   Processed MC points  : " << NAddedPoints << "\n";
 
    for (Int_t iStation = 0; iStation < StationSet->GetNStations(); ++iStation) {
        BmnSiliconStation* station = StationSet->GetSiliconStation(iStation);
 
        for (Int_t iModule = 0; iModule < station->GetNModules(); ++iModule) {
            BmnSiliconModule* module = station->GetModule(iModule);
            vector<Int_t> processed_zones;
 
            for (Int_t iLayer = 0; iLayer < module->GetNStripLayers(); ++iLayer) {
                BmnSiliconLayer layer = module->GetStripLayer(iLayer);
 
                Int_t zone_id = layer.GetZoneID();
 
                Bool_t is_processed_zone = false;
                for (size_t iz = 0; iz < processed_zones.size(); ++iz)
                    if (zone_id == processed_zones[iz])
                        is_processed_zone = true;
 
                if (!is_processed_zone) {
                    Int_t first_strip_in_zone = module->GetFirstStripInZone(zone_id);
                    Int_t last_strip_in_zone = module->GetLastStripInZone(zone_id);
                    map<int, int> noiseMap;   // AZ-270823 - add noise to strips (including left neighbour)
 
                    for (Int_t iStrip = first_strip_in_zone; iStrip < last_strip_in_zone + 1; ++iStrip) {
                        Double_t signal = module->GetStripSignalInZone(zone_id, iStrip);
 
                        if (signal > 0.0) {
                            signal = signal / (120. / 500.0);   // el to ADC
                            signal *= 1.17;                     // AZ-110923 - for beam data
                            // AZ-130923 - gain fluctuations
                            //*
                            Double_t gain = 1;
                            // AZ-191223 gain += gRandom->Gaus(0, 0.5);
                            if (gain < 0)
                                gain = 0;
                            signal *= gain;
                            // AZ-240823 if (signal <= 12) continue;
                            // AZ-240823--->
                            if (fSigma > 0) {
                                // Add noise
                                signal += gRandom->Gaus(0, fSigma);
                            }
                            // AZ-130923 - gain fluctuations
                            /*
                            Double_t gain = 1;
                                        //AZ-161223 gain += gRandom->Gaus(0, 0.2); //20%
                                        gain += gRandom->Gaus(0, 0.5); //50%
                            if (gain < 0) gain = 0;
                                        signal *= gain;
                            */
                            noiseMap[iStrip] = 1;   // noise flag
                            if (signal <= fThresh)
                                continue;
                            //<---AZ
                            if (signal >= 2048)
                                signal = 2047;
                            new ((*fBmnSiliconDigitsArray)[fBmnSiliconDigitsArray->GetEntriesFast()]) BmnSiliconDigit(
                                iStation, iModule, zone_id, iStrip, signal);   // zone_id == layer_num !!!!!!
 
                            if (fStripMatching) {
                                new ((*fBmnSiliconDigitMatchesArray)[fBmnSiliconDigitMatchesArray->GetEntriesFast()])
                                    BmnMatch(module->GetStripMatchInZone(zone_id, iStrip));
                            }
                            // AZ-270823 - add noise to empty channel to the leht from the active one
                            // if (fSigma > 0 && iStrip > first_strip_in_zone && noiseMap[iStrip-1] < 1) {
                            if (0) {
                                noiseMap[iStrip - 1] = 1;
                                signal = gRandom->Gaus(0, fSigma);
                                if (signal <= fThresh)
                                    continue;
 
                                if (signal >= 2048)
                                    signal = 2047;
                                new ((*fBmnSiliconDigitsArray)[fBmnSiliconDigitsArray->GetEntriesFast()])
                                    BmnSiliconDigit(iStation, iModule, zone_id, iStrip - 1,
                                                    signal);   // zone_id == layer_num !!!!!!
 
                                if (fStripMatching) {
                                    new (
                                        (*fBmnSiliconDigitMatchesArray)[fBmnSiliconDigitMatchesArray->GetEntriesFast()])
                                        BmnMatch(module->GetStripMatchInZone(zone_id, iStrip - 1));
                                }
                            }
 
                        } else if (fSigma > 0) {
                            //} else if (0) {
                            // AZ-270823 - add noise to empty channel to the right from the active one
                            // if (iStrip < last_strip_in_zone && noiseMap[iStrip-1] > 0) {
                            if (1) {
                                signal += gRandom->Gaus(0, fSigma);
                                // AZ-130923 - gain fluctuations
                                /*
                                Double_t gain = 1;
                                //AZ-161223 gain += gRandom->Gaus(0, 0.2);
                                gain += gRandom->Gaus(0, 0.5); //50%
                                if (gain < 0) gain = 0;
                                signal *= gain;
                                */
                                if (signal <= fThresh)
                                    continue;
                                if (signal >= 2048)
                                    signal = 2047;
                                new ((*fBmnSiliconDigitsArray)[fBmnSiliconDigitsArray->GetEntriesFast()])
                                    BmnSiliconDigit(iStation, iModule, zone_id, iStrip,
                                                    signal);   // zone_id == layer_num !!!!!!
 
                                if (fStripMatching) {
                                    new (
                                        (*fBmnSiliconDigitMatchesArray)[fBmnSiliconDigitMatchesArray->GetEntriesFast()])
                                        BmnMatch(module->GetStripMatchInZone(zone_id, iStrip));
                                }
                            }
                        }
                    }
 
                    processed_zones.push_back(zone_id);
                }
            }
        }
    }
 
    StationSet->Reset();
}
 
void BmnSiliconDigitizer::Finish()
{
    if (StationSet) {
        delete StationSet;
        StationSet = nullptr;
    }
 
    if (GlobalModTransfSet) {
        delete GlobalModTransfSet;
        GlobalModTransfSet = nullptr;
    }
 
    cout << "Work time of the Silicon digitizer: " << workTime << endl;
}