BmnGemRaw2Digit.cxx

Go to the documentation of this file.

#include "BmnGemRaw2Digit.h"
 
#include "TH1F.h"
// Boost
#include <boost/program_options.hpp>
// BmnRoot
#include <BmnADCDigit.h>
#include <BmnGemStripDigit.h>
#include <BmnRawDataDecoder.h>
#include <UniDetectorParameter.h>
 
namespace po = boost::program_options;
 
BmnGemRaw2Digit::BmnGemRaw2Digit(Int_t period, Int_t run, TString mapFileName, BmnSetup bmnSetup)
    : BmnAdcProcessor(period, run, "GEM", ADC_N_CHANNELS, ADC32_N_SAMPLES, bmnSetup)
    , fGemStationSet(BmnGemStripStationSet::Create(fPeriod, fSetup))
    , fMapFileName(mapFileName)
{
    InitStripNoiseArrays(fGemStationSet);
    InitStripRanges(fGemStationSet);
 
    Byte_t nMods = 4;
    fBigHot.resize(nMods);
    fBig.resize(nMods);
 
    if ((fPeriod >= 8) && (fSetup != kSRCSETUP)) {
        ReadLocalMap("X0_Bottom_Left.txt", fBigHot[2], 2, 2);
        ReadLocalMap("X0_Bottom_Right.txt", fBigHot[3], 2, 3);
        ReadLocalMap("X_Bottom_Left.txt", fBig[2], 0, 2);
        ReadLocalMap("X_Bottom_Right.txt", fBig[3], 0, 3);
        ReadLocalMap("Y0_Bottom_Left.txt", fBigHot[2], 3, 2);
        ReadLocalMap("Y0_Bottom_Right.txt", fBigHot[3], 3, 3);   // 506
        ReadLocalMap("Y_Bottom_Left.txt", fBig[2], 1, 2);
        ReadLocalMap("Y_Bottom_Right.txt", fBig[3], 1, 3);
 
        ReadLocalMap("X0_Top_Left.txt", fBigHot[0], 2, 0);
        ReadLocalMap("X0_Top_Right.txt", fBigHot[1], 2, 1);
        ReadLocalMap("X_Top_Left.txt", fBig[0], 0, 0);
        ReadLocalMap("X_Top_Right.txt", fBig[1], 0, 1);
        ReadLocalMap("Y0_Top_Left.txt", fBigHot[0], 3, 0);
        ReadLocalMap("Y0_Top_Right.txt", fBigHot[1], 3, 1);   // 506
        ReadLocalMap("Y_Top_Left.txt", fBig[0], 1, 0);
        ReadLocalMap("Y_Top_Right.txt", fBig[1], 1, 1);
    } else {
        ReadMap("GEM_X0_middle", fMid, 2, 0);
        ReadMap("GEM_Y0_middle", fMid, 3, 0);
        ReadMap("GEM_X1_middle", fMid, 0, 0);
        ReadMap("GEM_Y1_middle", fMid, 1, 0);
 
        ReadMap("GEM_X0_Big_Left", fBigHot[1], 2, 1);
        ReadMap("GEM_Y0_Big_Left", fBigHot[1], 3, 1);
        ReadMap("GEM_X1_Big_Left", fBig[1], 0, 1);
        ReadMap("GEM_Y1_Big_Left", fBig[1], 1, 1);
        ReadMap("GEM_X0_Big_Right", fBigHot[0], 2, 0);
        ReadMap("GEM_Y0_Big_Right", fBigHot[0], 3, 0);
        ReadMap("GEM_X1_Big_Right", fBig[0], 0, 0);
        ReadMap("GEM_Y1_Big_Right", fBig[0], 1, 0);
    }
    ReadLocalMap("GEM_X_small.txt", fSmall, 0, 0);
    ReadLocalMap("GEM_Y_small.txt", fSmall, 1, 0);
    ReadGlobalMap(fMapFileName);
    InitChipNoiseArrays();
}
 
void BmnGemRaw2Digit::InitStripRanges(unique_ptr<BmnGemStripStationSet>& ss)
{
    for (int32_t iSt = 0; iSt < ss->GetNStations(); ++iSt) {
        auto* st = ss->GetStation(iSt);
        for (int32_t iMod = 0; iMod < st->GetNModules(); ++iMod) {
            auto* mod = st->GetModule(iMod);
            MapLayer2MaxStrip layer2MaxStrip;
            for (int32_t iLay = 0; iLay < mod->GetNStripLayers(); ++iLay) {
                auto& lay = mod->GetStripLayer(iLay);
                layer2MaxStrip[iLay] = lay.GetNStrips() - 1;
            }
            fSM2MapMaxStrip[{iSt, iMod}] = move(layer2MaxStrip);
        }
    }
}
 
BmnStatus BmnGemRaw2Digit::ReadGlobalMap(TString FileName)
{
    string dummy;
    Int_t id = 0;
    UInt_t ser = 0;
    UInt_t ch_lo = 0;
    UInt_t ch_hi = 0;
    UInt_t iStation = 0;
    UInt_t iMod = 0;
    UInt_t zone = 0;
    UInt_t side = 0;
    UInt_t type = 0;
 
    vector<string> channel_map;
    vector<string> spec_thr;
    Double_t final_thr = 0;
    Double_t thr_dif = 0;
    Int_t n_iters = 0;
    Double_t cmod_cut = 0;
 
    // Setup options.
    po::options_description desc("Options");
    desc.add_options()("CHANNEL_MAP.channels", po::value<vector<string>>(&channel_map)->multitoken(),
                       "ADC channel -> strip map")("SIGNAL_CONFIGURATION.spec_thr",
                                                   po::value<vector<string>>(&spec_thr)->multitoken(),
                                                   "ADC channel -> strip map")(
        "SIGNAL_CONFIGURATION.Threshold", po::value<Double_t>(&final_thr)->default_value(25), "ADC final threshold")(
        "SIGNAL_CONFIGURATION.ThresholdDif", po::value<Double_t>(&thr_dif)->default_value(10), "ADC final threshold")(
        "SIGNAL_CONFIGURATION.NIterations", po::value<Int_t>(&n_iters)->default_value(5), "ADC final threshold")(
        "SIGNAL_CONFIGURATION.CModCut", po::value<Double_t>(&cmod_cut)->default_value(0), "ADC final threshold");
 
    // Load config file.
    po::variables_map vm;
    TString name = TString(getenv("VMCWORKDIR")) + TString("/input/") + fMapFileName;
    LOG(info) << "Loading GEM Map: Period " << fPeriod << ", Run " << fRun << " " << name;
    ifstream config_file(name.Data(), ifstream::in);
    if (!config_file.is_open()) {
        LOG(error) << "Error opening map-file (" << name << ")!";
        return kBMNERROR;
    }
    po::store(po::parse_config_file(config_file, desc), vm);
    config_file.close();
    po::notify(vm);
 
    // set thresholds
    SetThreshold(final_thr, thr_dif, n_iters, cmod_cut);
 
    // set specific thresholds
    Int_t ch = 0;
    Double_t thr = 0;
    for (auto& str : spec_thr) {
        istringstream ss(str);
        ss >> std::hex >> ser >> std::dec >> ch >> thr;
        fSpecThreshMap.insert(make_pair(make_pair(ser, ch), thr));
        LOGF(debug, "Serial  0x%08X Channel %2d  Threshold %3.3f", ser, ch, thr);
    }
 
    // set channel map
    for (auto& str : channel_map) {
        istringstream ss(str);
        ss >> std::hex >> ser >> std::dec >> ch_lo >> ch_hi >> id >> side >> type >> iStation >> iMod >> zone;
        GemMapLine* record = new GemMapLine();
        record->Ch_hi = ch_hi / GetNSamples();
        record->Serial = ser;
        record->Ch_lo = ch_lo / GetNSamples();
        record->Station = iStation;
        record->Zone = zone;
        record->GEM_id = id;
        record->Side = side;
        record->Type = type;
        record->Module = iMod;
        // adjust channel shifts
        if (record->Module > 1)
            record->Side += 2;
        if (record->GEM_id < 5) {
            if (record->GEM_id < 0) {   // small GEM
                record->Channel_shift = -512;
                record->Map_arr_ptr = &fSmall;
            } else {
                record->Map_arr_ptr = &fMid;
                if ((record->Ch_hi - record->Ch_lo) < 128 / GetNSamples())
                    record->Channel_shift = 2048 - record->Ch_lo * GetNSamples();
            }
        } else {
            if (record->Zone == 0) {   // hot zone
                record->Channel_shift = (record->Ch_lo == 0) ? 0 : -1024;
                if (fPeriod < 8)
                    record->Channel_shift += 1024;
                // For 6-7 run
                // Int_t chShift = (gemM->Ch_lo == 0) ? 1024 : 0;
                record->Map_arr_ptr = &fBigHot[record->Side];
            } else {   // big zone
                record->Map_arr_ptr = &fBig[record->Side];
                if ((record->Ch_hi - record->Ch_lo) < 128 / GetNSamples())
                    record->Channel_shift = 2048 - record->Ch_lo * GetNSamples();
            }
        }
        GrabNewSerial(ser);
        auto it = fOuterMap.find(ser);
        if (it == fOuterMap.end()) {   // create inner channel map for the serial
            InChanMap inner;
            inner.insert(make_pair(record->Ch_lo - 1, nullptr));
            inner.insert(make_pair(record->Ch_hi, record));
            fOuterMap.insert(make_pair(ser, move(inner)));
        } else {   // add range to the existing inner channel map
            InChanMap& inner = it->second;
            auto innerItHi = inner.find(record->Ch_hi);
            auto innerItLo = inner.find(record->Ch_lo - 1);
            if (innerItHi == inner.end())
                inner.insert(make_pair(record->Ch_hi, record));
            else if (innerItHi->second == nullptr) {
                inner.erase(innerItHi);
                inner.insert(make_pair(record->Ch_hi, record));
            } else {
                fprintf(stderr, "Wrong %s map! Overlapping intervals for %08X!\n", fDetName.Data(), ser);
                return kBMNERROR;
            }
            if (innerItLo == inner.end())
                inner.insert(make_pair(record->Ch_lo - 1, nullptr));
        }
    }
    return kBMNSUCCESS;
}
 
BmnStatus BmnGemRaw2Digit::ReadLocalMap(TString FileName, StripMapArray& m, Int_t lay, Int_t mod)
{
    TString name = TString(getenv("VMCWORKDIR")) + TString("/input/") + FileName;
    // printf("%s\n", name.Data());
    ifstream inFile(name.Data());
    if (!inFile.is_open()) {
        LOG(error) << "Error opening map-file (" << name << ")!";
        return kBMNERROR;
    }
 
    Int_t iStrip = 0;
    Int_t chan = 0;
    while (inFile >> chan) {
        m[chan] = BmnGemMap(iStrip++, lay, mod);
        // printf("\t\t strip %4d  lay %4d mod %4d  chan %4d\n", iStrip, lay, mod, chan);
    }
 
    return kBMNSUCCESS;
}
 
BmnStatus BmnGemRaw2Digit::ReadMap(TString parName, StripMapArray& m, Int_t lay, Int_t mod)
{
    // Int_t size = 0;
    UniDetectorParameter* par = UniDetectorParameter::GetDetectorParameter("GEM", parName, GetPeriod(), GetRun());
    vector<UniValue*> iiArr;
    if (par != NULL)
        par->GetValue(iiArr);
    delete par;
    //    printf("%20s  mod %d lay %d\n", parName.Data(), mod, lay);
    for (size_t i = 0; i < iiArr.size(); ++i) {
        IIValue* pValue = (IIValue*)iiArr[i];
        m[pValue->value2] = BmnGemMap(pValue->value1 - 1, lay, mod);   // Strip begins from 0
        //        printf("\t\t strip %4d  lay %4d mod %4d  chan %4d %s\n", pValue->value1, lay, mod, pValue->value2,
        //        (pValue->value1 == 0) ? "WFT!!!" : "");
    }
 
    if (!iiArr.empty())
        for (size_t i = 0; i < iiArr.size(); i++)
            delete iiArr[i];
    return kBMNSUCCESS;
}
 
BmnGemRaw2Digit::~BmnGemRaw2Digit()
{
    for (auto& it : fOuterMap)
        for (auto& inner : it.second)
            delete inner.second;
    DeleteStripNoiseArrays(fGemStationSet);
}
 
BmnStatus BmnGemRaw2Digit::FillProfiles(TClonesArray* adc)
{
    (this->*PrecalcEventModsImp)(adc);
#ifdef BUILD_DEBUG
    CalcEventMods();
#else
    CalcEventMods_simd();
#endif
    ProcessAdc(nullptr, kTRUE);
    return kBMNSUCCESS;
}
 
BmnStatus BmnGemRaw2Digit::FillNoisyChannels()
{
    const Int_t kNStripsInBunch = BmnRawDataDecoder::samplesGem;
    // cout << "GEM: " << kNStripsInBunch << endl;
    //     const Int_t kNBunches = kNStrips / kNStripsInBunch;
    const double kNThresh = BmnRawDataDecoder::threshGem;
 
    for (size_t iCr = 0; iCr < GetNSerials(); ++iCr) {
        auto it = fOuterMap.find(GetSerials()[iCr]);
        if (it == fOuterMap.end())
            continue;
        for (uint16_t iCh = 0; iCh < GetNChannels(); ++iCh) {
            InChanMap& inner = it->second;
            auto innerIt = inner.lower_bound(iCh);
            if (innerIt == inner.end())
                continue;
            GemMapLine* rec = innerIt->second;
            if (!rec)
                continue;
            Int_t iStation = rec->Station;
            Int_t iMod = rec->Module;
            MapLayer2MaxStrip& layer2MaxStrip = fSM2MapMaxStrip[{iStation, iMod}];
            for (Int_t iSmpl = 0; iSmpl < GetNSamples(); ++iSmpl)
                if (GetNoisyChipChannels()[iCr][iCh][iSmpl] == kTRUE) {
                    Int_t strip = -1;
                    Int_t iLayer = -1;
                    MapStrip(rec, iCh, iSmpl, iLayer, strip);
                    // check range violation
                    if (strip < 0)
                        continue;
                    if (strip > layer2MaxStrip[iLayer]) {
                        LOGF(debug2,
                             "fNoisyChannels strip outside range!\n\t"
                             " Station: %d,  Module: %2d, Layer: %2d iStrip: %d set true",
                             iStation, iMod, iLayer, strip);
                        LOGF(debug2, "lay  : [%3d, %3d]", 0, layer2MaxStrip[iLayer]);
                        continue;
                    }
                    fNoisyChannels[iStation][iMod][iLayer][strip] = kTRUE;
                }
            //                }
        }
    }
    LOGF(debug, "Chip noise channels to strip finished");
    for (Int_t iSt = 0; iSt < fGemStationSet->GetNStations(); ++iSt) {
        if ((iSt == 7) && (fPeriod == 9))   // crutch to avoid marking strips on the beam noisy
            continue;
        auto* st = fGemStationSet->GetStation(iSt);
        for (Int_t iMod = 0; iMod < st->GetNModules(); ++iMod) {
            auto* mod = st->GetModule(iMod);
            for (Int_t iLay = 0; iLay < mod->GetNStripLayers(); ++iLay) {
                TH1F* prof = fSigProf[iSt][iMod][iLay];
                auto& lay = mod->GetStripLayer(iLay);
                Int_t kNBunches = lay.GetNStrips() / kNStripsInBunch;
                for (Int_t iBunch = 0; iBunch < kNBunches; ++iBunch) {
                    Double_t mean = 0.0;
                    for (Int_t iStrip = 0; iStrip < kNStripsInBunch; ++iStrip) {
                        Int_t strip = iStrip + iBunch * kNStripsInBunch;
                        if (fNoisyChannels[iSt][iMod][iLay][strip] == kTRUE)
                            continue;
                        Double_t curr = prof->GetBinContent(strip + 1);
                        // Double_t prev = prof->GetBinContent(strip);
                        mean += curr;
                    }
                    mean /= kNStripsInBunch;
                    for (Int_t iStrip = 0; iStrip < kNStripsInBunch; ++iStrip) {
                        Int_t strip = iStrip + iBunch * kNStripsInBunch;
                        if (fNoisyChannels[iSt][iMod][iLay][strip] == kTRUE)
                            continue;
                        Double_t curr = prof->GetBinContent(strip + 1);
                        // Double_t prev = prof->GetBinContent(strip);
                        //                         if (kNThresh * meanDiff < curr - prev)
                        if (kNThresh * mean < Abs(curr - mean))
                            fNoisyChannels[iSt][iMod][iLay][strip] = kTRUE;
                    }
                }
            }
        }
    }
    LOGF(debug, "Marking noisy strips finished");
    for (size_t iCr = 0; iCr < GetNSerials(); ++iCr) {
        auto it = fOuterMap.find(GetSerials()[iCr]);
        if (it == fOuterMap.end())
            continue;
        for (uint16_t iCh = 0; iCh < GetNChannels(); ++iCh) {
            InChanMap& inner = it->second;
            auto innerIt = inner.lower_bound(iCh);
            if (innerIt == inner.end())
                continue;
            GemMapLine* rec = innerIt->second;
            if (!rec)
                continue;
            Int_t iStation = rec->Station;
            Int_t iMod = rec->Module;
            MapLayer2MaxStrip& layer2MaxStrip = fSM2MapMaxStrip[{iStation, iMod}];
            for (Int_t iSmpl = 0; iSmpl < GetNSamples(); ++iSmpl) {
                Int_t strip = -1;
                Int_t iLayer = -1;
                MapStrip(rec, iCh, iSmpl, iLayer, strip);
                // check range violation
                if (strip < 0)
                    continue;
                if (strip > layer2MaxStrip[iLayer]) {
                    LOGF(debug2,
                         "fNoisyChannels strip outside range!\n\t"
                         " Station: %d,  Module: %2d, Layer: %2d iStrip: %d set true back",
                         iStation, iMod, iLayer, strip);
                    LOGF(debug2, "lay  : [%3d, %3d]", 0, layer2MaxStrip[iLayer]);
                    continue;
                }
                if (fNoisyChannels[iStation][iMod][iLayer][strip] == kTRUE)
                    GetNoisyChipChannels()[iCr][iCh][iSmpl] = kTRUE;
            }
        }
    }
    LOGF(debug, "noisy strips to chips finished");
    return kBMNSUCCESS;
}
 
void BmnGemRaw2Digit::ProcessAdc(TClonesArray* gem, Bool_t doFill)
{
    for (size_t iCr = 0; iCr < GetNSerials(); ++iCr) {
        auto it = fOuterMap.find(GetSerials()[iCr]);
        if (it == fOuterMap.end())
            continue;
        //        LOGF(info, "cr[%2lu] =  0x08X found", iCr, GetSerials()[iCr]);
        for (uint16_t iCh = 0; iCh < fNChannels; ++iCh) {
            auto innerIt = it->second.lower_bound(iCh);
            if (innerIt == it->second.end())
                continue;
            GemMapLine* rec = innerIt->second;
            if (!rec)
                continue;
            Int_t iStation = rec->Station;
            Int_t iMod = rec->Module;
            MapLayer2MaxStrip& layer2MaxStrip = fSM2MapMaxStrip[{iStation, iMod}];
            for (Int_t iSmpl = 0; iSmpl < GetNSamples(); ++iSmpl) {
                if (GetNoisyChipChannels()[iCr][iCh][iSmpl] == kTRUE)
                    continue;
                Int_t strip = -1;
                Int_t iLayer = -1;
                MapStrip(rec, iCh, iSmpl, iLayer, strip);
                // cout << "DECODER INFO: " <<  station << " " << mod << " " << layer << " " << strip << endl;
                if ((strip < 0) || (strip > layer2MaxStrip[iLayer])) {
                    LOGF(debug2,
                         "fNoisyChannels strip outside range!\n\t"
                         " Station: %d,  Module: %2d, Layer: %2d iStrip: %d  digi add",
                         iStation, iMod, iLayer, strip);
                    LOGF(debug2, "info  : [%3d, %3d]", 0, layer2MaxStrip[iLayer]);
                    continue;
                }
 
                Double_t sig = fAdc[iCr][iCh][iSmpl] - fPedVal[iCr][iCh][iSmpl] + fCMode[iCr][iCh] - fSMode[iCr][iCh];
                // Double_t Asig = TMath::Abs(sig);
                Double_t SpecThr(0);
                auto itThr = fSpecThreshMap.find(make_pair(GetSerials()[iCr], iCh));
                if (itThr != fSpecThreshMap.end())
                    SpecThr = itThr->second;
 
                Double_t thr = Max(FinalThr, 3.5 * GetPedestalsRMS()[iCr][iCh][iSmpl]);
                thr = Max(thr, SpecThr);
                if ((!GetApplyThreshold()) || (sig > thr)) {
                    if (doFill /* && (Abs(fCMode[iCr][iCh] - fSMode[iCr][iCh]) < cmodcut)*/) {
                        fSigProf[iStation][iMod][iLayer]->Fill(strip);
                    } else {
                        BmnGemStripDigit* resDig = new ((*gem)[gem->GetEntriesFast()])
                            BmnGemStripDigit(iStation, iMod, iLayer, strip, sig, GetPedestalsRMS()[iCr][iCh][iSmpl]);
                        //                                if ((Abs(fCMode[iCr][iCh] - fSMode[iCr][iCh]) >
                        //                                cmodcut))
                        //                                    resDig->SetIsGoodDigit(kFALSE);
                        //                                else
                        resDig->SetIsGoodDigit(kTRUE);
                    }
                }
            }
        }
    }
}
 
BmnStatus BmnGemRaw2Digit::FillEvent(TClonesArray* adc, TClonesArray* gem)
{
    (this->*PrecalcEventModsImp)(adc);
#ifdef BUILD_DEBUG
    CalcEventMods();
#else
    CalcEventMods_simd();
#endif
    ProcessAdc(gem, kFALSE);
    return kBMNSUCCESS;
}
 
inline void BmnGemRaw2Digit::MapStrip(GemMapLine* gemM, UInt_t ch, Int_t iSmpl, Int_t& lay, Int_t& strip)
{
    size_t realChannel = ch * GetNSamples() + iSmpl + gemM->Channel_shift;
    lay = (*gemM->Map_arr_ptr)[realChannel].lay;
    strip = (*gemM->Map_arr_ptr)[realChannel].strip;
    return;
}