Skip to content
Snippets Groups Projects
analyze.cpp 6.69 KiB
Newer Older
  • Learn to ignore specific revisions
  • A-Compton's avatar
    A-Compton committed
    #include <iostream>
    #include <iomanip>
    
    Chao's avatar
    Chao committed
    #include <fstream>
    
    Chao Peng's avatar
    Chao Peng committed
    #include <algorithm>
    
    #include <exception>
    
    Chao's avatar
    Chao committed
    #include "TTree.h"
    #include "TFile.h"
    
    #include "EvChannel.h"
    
    #include "Fadc250Decoder.h"
    
    Chao Peng's avatar
    Chao Peng committed
    #include "WfAnalyzer.h"
    
    Zhiwen Zhao's avatar
    Zhiwen Zhao committed
    #include "SSPDecoder.h"
    
    #include "ConfigArgs.h"
    #include "read_modules.h"
    
    Chao's avatar
    Chao committed
    
    #define PROGRESS_COUNT 1000
    
    
    
    Chao Peng's avatar
    Chao Peng committed
    void write_raw_data(const std::string &dpath, const std::string &opath, const std::string &mpath, int nev,
    
                        int res = 3, double thres = 20, int npeds = 5, double flat = 1.0);
    
    Chao's avatar
    Chao committed
    
    
    // event types
    
    Chao's avatar
    Chao committed
    enum EvType {
        CODA_PRST = 0xffd1,
        CODA_GO = 0xffd2,
        CODA_END = 0xffd4,
        CODA_PHY1 = 0xff50,
        CODA_PHY2 = 0xff70,
    };
    
    
    Chao Peng's avatar
    Chao Peng committed
    int main(int argc, char*argv[])
    
    Chao's avatar
    Chao committed
    {
    
    Chao Peng's avatar
    Chao Peng committed
        ConfigArgs arg_parser;
        arg_parser.AddHelps({"-h", "--help"});
    
        arg_parser.AddPositional("raw_data", "raw data in evio format");
        arg_parser.AddPositional("root_file", "output data file in root format");
    
        arg_parser.AddArg<int>("-n", "nev", "number of events to process (< 0 means all)", -1);
    
        arg_parser.AddArgs<std::string>({"-m", "--module"}, "module", "json file for module configuration",
    
    Chao Peng's avatar
    Chao Peng committed
                                        "database/modules/mapmt_module.json");
    
        arg_parser.AddArg<int>("-r", "res", "resolution for waveform analysis", 3);
        arg_parser.AddArg<double>("-t", "thres", "peak threshold for waveform analysis", 20.0);
        arg_parser.AddArg<int>("-p", "npeds", "sample window width for pedestal searching", 8);
        arg_parser.AddArg<double>("-f", "flat", "flatness requirement for pedestal searching", 1.0);
    
    Chao Peng's avatar
    Chao Peng committed
    
        auto args = arg_parser.ParseArgs(argc, argv);
    
    
    Chao Peng's avatar
    Chao Peng committed
        for(auto &it : args) {
            std::cout << it.first << ": " << it.second.String() << std::endl;
    
    Chao's avatar
    Chao committed
        }
    
        write_raw_data(args["raw_data"].String(),
                       args["root_file"].String(),
    
    Chao Peng's avatar
    Chao Peng committed
                       args["module"].String(),
                       args["nev"].Int(),
                       args["res"].Int(),
                       args["thres"].Double(),
                       args["npeds"].Int(),
    
                       args["flat"].Double());
    
        return 0;
    }
    
    Chao's avatar
    Chao committed
    
    
    Chao Peng's avatar
    Chao Peng committed
    // read raw data in evio format, and extract information
    
    Chao Peng's avatar
    Chao Peng committed
    void write_raw_data(const std::string &dpath, const std::string &opath, const std::string &mpath, int nev,
    
                        int res, double thres, int npeds, double flat)
    
    Chao's avatar
    Chao committed
    {
    
        auto modules = read_modules(mpath);
    
        if (modules.empty()) {
            std::cout << "Cannot find any modules in configuration file \"" << mpath << "\"" << std::endl;
            return;
        }
    
    
    Chao Peng's avatar
    Chao Peng committed
        // raw data
    
        evc::EvChannel evchan;
    
        if (evchan.Open(dpath) != evc::status::success) {
            std::cout << "Cannot open evchannel at " << dpath << std::endl;
            return;
        }
    
    Chao's avatar
    Chao committed
    
    
    Chao Peng's avatar
    Chao Peng committed
        // output
    
    Chao's avatar
    Chao committed
        auto *hfile = new TFile(opath.c_str(), "RECREATE", "MAPMT test results");
        auto tree = new TTree("EvTree", "Cherenkov Test Events");
    
    Chao Peng's avatar
    Chao Peng committed
    
    
        // loop over the modules and set branch address for channels
    
    Chao Peng's avatar
    Chao Peng committed
        std::vector<uint32_t> dbanks;
    
        for (auto &m : modules) {
    
    Chao Peng's avatar
    Chao Peng committed
            if (std::find(dbanks.begin(), dbanks.end(), m.bank) == dbanks.end()) {
                dbanks.push_back(m.bank);
            }
    
            switch (m.type) {
            case kFADC250:
                {
                    // 16 channels
                    auto event = new fdec::Fadc250Event(0, 16);
                    m.event = static_cast<void*>(event);
                    for (auto &ch : m.channels) {
                        auto n = ch.name;
                        tree->Branch(n.c_str(), &event->channels[ch.id], 32000, 0);
                    }
                }
                break;
            case kSSP:
                {
                    // assign a large enough buffer
                    auto event = new ssp::SSPEvent(10000);
                    m.event = static_cast<void*>(event);
    
    Zhiwen Zhao's avatar
    Zhiwen Zhao committed
                    tree->Branch("SSP_evt"),        &event->tTrigNum,   "evt/I");
                    tree->Branch("SSP_trigtime"),   &event->tTrigTime,  "trigtime/D");
                    tree->Branch("SSP_nedge"),      &event->Nedges,     "nedge/I");
                    tree->Branch("SSP_ch"),         &event->channel[0], "ch[nedge]/I");
                    tree->Branch("SSP_pol"),        &event->edge[0],    "pol[nedge]/I");
                    tree->Branch("SSP_time"),       &event->time[0],    "time[nedge]/I");
                    tree->Branch("SSP_fiber"),      &event->fiber[0],   "fiber[nedge]/I");
                    tree->Branch("SSP_slot"),       &event->slot[0],    "slot[nedge]/I");
    
                }
                break;
            default:
                std::cout << "Unsupported module type " << m.type << std::endl;
                break;
    
        // waveform analyzer
    
    Chao Peng's avatar
    Chao Peng committed
        fdec::Analyzer analyzer(res, thres, npeds, flat);
    
        // decoders
        fdec::Fadc250Decoder fdecoder;
        ssp::SSPDecoder sdecoder;
    
    Chao Peng's avatar
    Chao Peng committed
        int count = 0;
    
    Chao Peng's avatar
    Chao Peng committed
        auto &ref = modules.front();
    
        while ((evchan.Read() == evc::status::success) && (nev-- != 0)) {
    
            if((count % PROGRESS_COUNT) == 0) {
    
    Chao's avatar
    Chao committed
                std::cout << "Processed events - " << count << "\r" << std::flush;
            }
    
    
            switch(evchan.GetEvHeader().tag) {
            // only want physics events
            case CODA_PHY1:
            case CODA_PHY2:
                break;
            case CODA_PRST:
            case CODA_GO:
            case CODA_END:
            default:
    
    A-Compton's avatar
    A-Compton committed
                continue;
            }
    
    
    Chao Peng's avatar
    Chao Peng committed
            evchan.ScanBanks(dbanks);
    
            // get block level
    
    Chao Peng's avatar
    Chao Peng committed
            int blvl = evchan.GetEvBuffer(ref.crate, ref.bank, ref.slot).size();
    
            const uint32_t *dbuf;
    
            for (int ii = 0; ii < blvl; ++ii) {
                // parse module data
                for (auto &mod : modules) {
    
                    // get data buffer
                    try {
                        dbuf = evchan.GetEvBuffer(mod.crate, mod.bank, mod.slot, ii, buflen);
                    } catch (std::exception &e) {
    
    Chao Peng's avatar
    Chao Peng committed
                        std::cout << "warning: " << e.what() << "\n";
    
                        continue;
                    }
                    // decode by module type
                    switch (mod.type) {
                    case kFADC250:
                        {
                            auto event = static_cast<fdec::Fadc250Event*>(mod.event);
                            fdecoder.DecodeEvent(*event, dbuf, buflen);
                            for (auto &ch : event->channels) {
                                analyzer.Analyze(ch);
                            }
                        }
                        break;
                    case kSSP:
                        sdecoder.DecodeEvent(*static_cast<ssp::SSPEvent*>(mod.event), dbuf, buflen);
                        break;
                    default:
                        std::cout << "Unsupported module type " << mod.type << std::endl;
                        break;
    
    Chao's avatar
    Chao committed
                }
    
                tree->Fill();
    
    Chao's avatar
    Chao committed
        }
        std::cout << "Processed events - " << count << std::endl;
    
    
    Chao Peng's avatar
    Chao Peng committed
        evchan.Close();
    
    Chao's avatar
    Chao committed
        hfile->Write();
        hfile->Close();
    }