SimReaderSpec.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
#include "SimReaderSpec.h"
 
#include "Framework/DataProcessorSpec.h"
#include "Framework/DataRefUtils.h"
#include "Framework/ControlService.h"
#include "Framework/ConfigParamRegistry.h"
#include "Framework/Lifetime.h"
#include "Steer/HitProcessingManager.h"
#include "SimulationDataFormat/InteractionSampler.h"
#include "CommonDataFormat/InteractionRecord.h"
#include "CommonUtils/NameConf.h"
#include "DataFormatsTPC/TPCSectorHeader.h"
#include <DataFormatsParameters/GRPLHCIFData.h>
#include "DetectorsRaw/HBFUtils.h"
#include <CCDB/BasicCCDBManager.h>
#include <fairlogger/Logger.h>
#include <memory>     // std::unique_ptr
#include <cstring>    // memcpy
#include <string>     // std::string
#include <cassert>
#include <chrono>
#include <thread>
#include <algorithm>
#include <filesystem>
#include <boost/interprocess/sync/named_semaphore.hpp>
 
using namespace o2::framework;
namespace o2lhc = o2::constants::lhc;
 
using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType;
namespace o2
{
namespace steer
{
 
std::vector<o2::ctp::CTPDigit>* ctptrigger = nullptr;
float gIntRate = -1.;
 
DataProcessorSpec getSimReaderSpec(SubspecRange range, const std::vector<std::string>& simprefixes, const std::vector<int>& tpcsectors, bool withTrigger)
{
  uint64_t activeSectors = 0;
  for (const auto& tpcsector : tpcsectors) {
    activeSectors |= (uint64_t)0x1 << tpcsector;
  }
 
  auto doit = [range, tpcsectors, activeSectors, withTrigger](ProcessingContext& pc) {
    auto& mgr = steer::HitProcessingManager::instance();
    auto& context = mgr.getDigitizationContext();
    auto eventrecords = context.getEventRecords();
 
    if (withTrigger) {
      // fetch the digits and transport them as part of the context
      LOG(info) << "Setting CTP trigger object to " << ctptrigger;
      context.setCTPDigits(ctptrigger);
    }
 
    // inject the global interaction rate information
    context.setDigitizerInteractionRate(gIntRate);
 
    for (auto const& sector : tpcsectors) {
      // Note: the TPC sector header was serving the sector to lane mapping before
      // now the only remaining purpose is to propagate the mask of valid sectors
      // in principle even that is not necessary any more because every sector has
      // a dedicated route bound to the sector numbers as subspecification and the
      // merging can be done based on that. However if we in the future go over to
      // a multipart scheme again, we again will need the information, so we keep it
      // For the moment, sector member in the sector header is in sync with
      // subspecification of the route
      o2::tpc::TPCSectorHeader header{sector};
      header.activeSectors = activeSectors;
      pc.outputs().snapshot(OutputRef{"collisioncontext", static_cast<SubSpecificationType>(sector), {header}},
                            context);
    }
 
    // the first 36 channel numbers are reserved for the TPC, now publish the remaining
    // channels
    for (int subchannel = range.min; subchannel < range.max; ++subchannel) {
      LOG(info) << "SENDING SOMETHING TO OTHERS";
      pc.outputs().snapshot(
        OutputRef{"collisioncontext", static_cast<SubSpecificationType>(subchannel)},
        context);
    }
    pc.outputs().snapshot(OutputRef{"bunchFilling"}, mgr.getInteractionSampler().getBunchFilling());
    // digitizer workflow runs only once
    // send endOfData control event and mark the reader as ready to finish
    pc.services().get<ControlService>().endOfStream();
    pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
  };
 
  // init function return a lambda taking a ProcessingContext
  auto initIt = [simprefixes, doit, withTrigger](InitContext& ctx) {
    // initialize fundamental objects
    auto& mgr = steer::HitProcessingManager::instance();
 
    if (withTrigger) {
      // fetch the ctp trigger/digit information from the CTP file
      auto triggerf = TFile(ctx.options().get<std::string>("triggerfile").c_str(), "OPEN");
      auto tr = (TTree*)triggerf.Get("o2sim");
      if (!tr) {
        LOG(fatal) << "Did not find CTP TTree";
      }
      auto br = tr->GetBranch("CTPDigits");
      if (!br) {
        LOG(fatal) << "Did not find CTPDigit branch";
      }
      br->SetAddress(&ctptrigger);
      br->GetEntry(0);
      LOG(info) << " Read " << ctptrigger->size() << " CTP digits ";
    }
 
    // init gRandom
    gRandom->SetSeed(ctx.options().get<int>("seed"));
 
    if (simprefixes.size() == 0) {
      LOG(error) << "No simulation prefix available";
    } else {
      LOG(info) << "adding " << simprefixes[0] << "\n";
      mgr.addInputFile(simprefixes[0]);
      for (int part = 1; part < simprefixes.size(); ++part) {
        mgr.addInputSignalFile(simprefixes[part]);
      }
    }
 
    gIntRate = ctx.options().get<float>("interactionRate"); // is interaction rate requested?
    if (gIntRate < 1.f) {
      gIntRate = 1.f;
    }
    // do we start from an existing context
    auto incontextstring = ctx.options().get<std::string>("incontext");
    LOG(info) << "INCONTEXTSTRING " << incontextstring;
    if (incontextstring.size() > 0) {
      auto success = mgr.setupRunFromExistingContext(incontextstring.c_str());
      if (!success) {
        LOG(fatal) << "Could not read collision context from " << incontextstring;
      }
    } else {
      LOG(info) << "Imposing hadronic interaction rate " << gIntRate << "Hz";
      mgr.getInteractionSampler().setInteractionRate(gIntRate);
      o2::raw::HBFUtils::Instance().print();
      o2::raw::HBFUtils::Instance().checkConsistency();
      mgr.getInteractionSampler().setFirstIR({0, o2::raw::HBFUtils::Instance().orbitFirstSampled});
      mgr.getDigitizationContext().setFirstOrbitForSampling(o2::raw::HBFUtils::Instance().orbitFirstSampled);
 
      auto setBCFillingHelper = [](auto& sampler, auto& bcPatternString) {
        if (bcPatternString == "ccdb") {
          LOG(info) << "Fetch bcPattern information from CCDB";
          // fetch the GRP Object
          auto& ccdb = o2::ccdb::BasicCCDBManager::instance();
          auto grpLHC = ccdb.get<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF");
          LOG(info) << "Fetched injection scheme " << grpLHC->getInjectionScheme() << " from CCDB";
          sampler.setBunchFilling(grpLHC->getBunchFilling());
        } else {
          sampler.setBunchFilling(bcPatternString);
        }
      };
 
      auto bcPatternFile = ctx.options().get<std::string>("bcPatternFile");
      if (!bcPatternFile.empty()) {
        setBCFillingHelper(mgr.getInteractionSampler(), bcPatternFile);
      }
 
      mgr.getInteractionSampler().init();
      mgr.getInteractionSampler().print();
 
      // doing a random event selection/subsampling?
      mgr.setRandomEventSequence(ctx.options().get<int>("randomsample") > 0);
 
      // finalize collisions (with number of collisions asked)
      auto col = ctx.options().get<int>("ncollisions");
      if (col != 0) {
        mgr.setupRun(col);
      } else {
        mgr.setupRun();
      }
 
      // --- we add QED contributions to the digitization context
      // --- for now in between first and last real collision
      auto qedprefix = ctx.options().get<std::string>("simPrefixQED");
      if (qedprefix.size() > 0) {
        o2::steer::InteractionSampler qedInteractionSampler;
        if (!bcPatternFile.empty()) {
          setBCFillingHelper(qedInteractionSampler, bcPatternFile);
        }
 
        // get first and last "hadronic" interaction records and let
        // QED events range from the first bunch crossing to the last bunch crossing
        // in this range
        auto first = mgr.getDigitizationContext().getEventRecords().front();
        auto last = mgr.getDigitizationContext().getEventRecords().back();
        first.bc = 0;
        last.bc = o2::constants::lhc::LHCMaxBunches;
 
        const float ratio = ctx.options().get<float>("qed-x-section-ratio");
        if (ratio <= 0.) {
          throw std::runtime_error("no meaningful qed-x-section-ratio was provided");
        }
        const float hadronicrate = ctx.options().get<float>("interactionRate");
        const float qedrate = ratio * hadronicrate;
        LOG(info) << "QED RATE " << qedrate;
        qedInteractionSampler.setInteractionRate(qedrate);
        qedInteractionSampler.setFirstIR(first);
        qedInteractionSampler.init();
        qedInteractionSampler.print();
        std::vector<o2::InteractionTimeRecord> qedinteractionrecords;
        o2::InteractionTimeRecord t;
        LOG(info) << "GENERATING COL TIMES";
        t = qedInteractionSampler.generateCollisionTime();
        while ((t = qedInteractionSampler.generateCollisionTime()) < last) {
          qedinteractionrecords.push_back(t);
        }
        LOG(info) << "DONE GENERATING COL TIMES";
 
        // get digitization context and add QED stuff
        mgr.getDigitizationContext().fillQED(qedprefix, qedinteractionrecords);
        mgr.getDigitizationContext().printCollisionSummary(true, 2000); // print with QED but truncate output
      }
      // --- end addition of QED contributions
 
      LOG(info) << "Initializing Spec ... have " << mgr.getDigitizationContext().getEventRecords().size() << " times ";
      LOG(info) << "Serializing Context for later reuse";
      mgr.writeDigitizationContext(ctx.options().get<std::string>("outcontext").c_str());
    }
 
    return doit;
  };
 
  std::vector<OutputSpec> outputs;
  for (auto const& tpcsector : tpcsectors) {
    outputs.emplace_back(
      OutputSpec{{"collisioncontext"}, "SIM", "COLLISIONCONTEXT", static_cast<SubSpecificationType>(tpcsector), Lifetime::Timeframe});
  }
  for (int subchannel = range.min; subchannel < range.max; ++subchannel) {
    outputs.emplace_back(
      OutputSpec{{"collisioncontext"}, "SIM", "COLLISIONCONTEXT", static_cast<SubSpecificationType>(subchannel), Lifetime::Timeframe});
  }
 
  outputs.emplace_back(OutputSpec{{"bunchFilling"}, "SIM", "BUNCHFILLING", 0, Lifetime::Timeframe});
 
  return DataProcessorSpec{
    /*ID*/ "SimReader",
    /*INPUT CHANNELS*/ Inputs{}, outputs,
    /* ALGORITHM */
    AlgorithmSpec{initIt},
    /* OPTIONS */
    Options{
      {"interactionRate", VariantType::Float, 50000.0f, {"Total hadronic interaction rate (Hz)"}},
      {"bcPatternFile", VariantType::String, "", {"Interacting BC pattern file (e.g. from CreateBCPattern.C)"}},
      {"simPrefixQED", VariantType::String, "", {"Sim (QED) input prefix (example: path/o2qed). The prefix allows to find files like path/o2qed_Kine.root etc."}},
      {"qed-x-section-ratio", VariantType::Float, -1.f, {"Ratio of cross sections QED/hadronic events. Determines QED interaction rate from hadronic interaction rate."}},
      {"outcontext", VariantType::String, "collisioncontext.root", {"Output file for collision context"}},
      {"incontext", VariantType::String, "", {"Take collision context from this file"}},
      {"triggerfile", VariantType::String, "ctpdigits.root", {"Name of the CTP trigger/digit file to use"}},
      {"seed", VariantType::Int, 0, {"Random seed for collision context generation"}},
      {"ncollisions,n",
       VariantType::Int,
       0,
       {"number of collisions to sample (default is given by number of entries in chain"}},
      {"randomsample", VariantType::Int, 0, {"Draw collisions random instead of linear sequence. (Default no = 0)"}}}};
}
} // namespace steer
} // namespace o2