diff --git a/Detectors/TPC/simulation/include/TPCSimulation/DigitMC.h b/Detectors/TPC/simulation/include/TPCSimulation/DigitMC.h index 5fd1f6ec05e8e..a5116cfe30c6a 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/DigitMC.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/DigitMC.h @@ -10,6 +10,7 @@ #endif #include "TPCBase/Digit.h" +#include "TPCSimulation/DigitPad.h" #include "FairTimeStamp.h" @@ -22,19 +23,38 @@ class access; namespace o2 { namespace TPC { +#ifdef TPC_DIGIT_USEFAIRLINKS + using DigitBase = FairTimeStamp; +#else + // A minimal (temporary) TimeStamp class, introduced here for + // reducing memory consumption to a minimum. + // This can be used only when MCtruth is not done using FairLinks. + class TimeStamp : public TObject { + public: + TimeStamp() {} + TimeStamp(int time) { + // we use the TObjectID for the time + SetUniqueID(time); + } + int GetTimeStamp() const { return TObject::GetUniqueID(); } + ClassDef(TimeStamp, 1); + }; + using DigitBase = TimeStamp; +#endif + /// \class DigitMC /// This is the definition of the Monte Carlo Digit object, which is the final entity after Digitization /// Its coordinates are defined by the CRU, the time bin, the pad row and the pad. /// It holds the ADC value of a given pad on the pad plane. /// Additional information attached to it are the MC label of the contributing tracks - -class DigitMC : public FairTimeStamp, public Digit { +class DigitMC : public DigitBase, public Digit { public: /// Default constructor DigitMC(); +#ifdef TPC_DIGIT_USEFAIRLINKS /// Constructor, initializing values for position, charge, time and common mode /// \param MClabel std::vector containing the MC event and track ID encoded in a long int /// \param cru CRU of the DigitMC @@ -44,13 +64,24 @@ class DigitMC : public FairTimeStamp, public Digit { /// \param time Time at which the DigitMC was created /// \param commonMode Common mode signal on that ROC in the time bin of the DigitMC. If not assigned, it is set to zero. DigitMC(int cru, float charge, int row, int pad, int time, float commonMode = 0.f); +#else + /// Constructor, initializing values for position, charge, time and common mode + /// \param MClabel std::vector containing the MC event and track ID encoded in a long int + /// \param cru CRU of the DigitMC + /// \param charge Accumulated charge of DigitMC + /// \param row Row in which the DigitMC was created + /// \param pad Pad in which the DigitMC was created + /// \param time Time at which the DigitMC was created + /// \param commonMode Common mode signal on that ROC in the time bin of the DigitMC. If not assigned, it is set to zero. + DigitMC(std::vector const &MClabel, int cru, float charge, int row, int pad, int time, float commonMode = 0.f); +#endif /// Destructor ~DigitMC() override = default; /// Get the timeBin of the DigitMC /// \return timeBin of the DigitMC - int getTimeStamp() const final { return static_cast(FairTimeStamp::GetTimeStamp()); }; + int getTimeStamp() const final { return static_cast(DigitBase::GetTimeStamp()); }; /// Get the common mode signal of the DigitMC /// \return common mode signal of the DigitMC @@ -60,7 +91,9 @@ class DigitMC : public FairTimeStamp, public Digit { #ifndef __CINT__ friend class boost::serialization::access; #endif - +#ifndef TPC_DIGIT_USEFAIRLINKS + std::vector mMClabel; ///< MC truth information to (multiple) event ID and track ID encoded in a long +#endif float mCommonMode; ///< Common mode value of the DigitMC ClassDefOverride(DigitMC, 3); @@ -68,17 +101,30 @@ class DigitMC : public FairTimeStamp, public Digit { inline DigitMC::DigitMC() - : FairTimeStamp() + : DigitBase() , Digit(-1, -1.f, -1, -1) , mCommonMode(0.f) -{} +#ifndef TPC_DIGIT_USEFAIRLINKS + , mMClabel() +#endif + {} +#ifdef TPC_DIGIT_USEFAIRLINKS inline DigitMC::DigitMC(int cru, float charge, int row, int pad, int time, float commonMode) - : FairTimeStamp(time) + : DigitBase(time) , Digit(cru, charge, row, pad) , mCommonMode(commonMode) {} +#else +inline +DigitMC::DigitMC(std::vector const &MClabel, int cru, float charge, int row, int pad, int time, float commonMode) + : DigitBase(time) + , Digit(cru, charge, row, pad) + , mMClabel(MClabel) + , mCommonMode(commonMode) +{} +#endif } } diff --git a/Detectors/TPC/simulation/include/TPCSimulation/DigitPad.h b/Detectors/TPC/simulation/include/TPCSimulation/DigitPad.h index 7ac8c7ddc6c50..2fc7d0085a121 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/DigitPad.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/DigitPad.h @@ -66,10 +66,30 @@ class DigitPad{ void fillOutputContainer(TClonesArray *output, int cru, int timeBin, int row, int pad, float commonMode = 0); private: + +#ifndef TPC_DIGIT_USEFAIRLINKS + /// The MC labels are sorted by occurrence such that the event/track combination with the largest number of occurrences is first + /// This is then dumped into a std::vector and attached to the digits + /// \todo Find out how many different event/track combinations are relevant + /// \param std::vector containing the sorted MCLabels + void processMClabels(std::vector &sortedMCLabels) const; +#endif + float mChargePad; ///< Total accumulated charge on that pad for a given time bin unsigned char mPad; ///< Pad of the ADC value #ifdef TPC_DIGIT_USEFAIRLINKS FairMultiLinkedData mMCLinks; ///< MC links +#else + // according to event + trackID + sorted according to most probable + std::map mMCID; //! Map containing the MC labels (key) and the according number of occurrence (value) + + // TODO: optimize this treatment, for example by using a structure like this + // struct MCIDValue { + // unsigned int eventId : 15; // 32k event Id possible + // unsigned int trackId: 17; // 128K tracks possible + // unsigned int occurences : 32; // 4G occurrences possible + // } + // std::vector mMCID; #endif }; @@ -85,12 +105,16 @@ DigitPad::DigitPad(int pad) inline void DigitPad::setDigit(size_t trackID, float charge) { + static FairRootManager *mgr = FairRootManager::Instance(); + const int eventID = mgr->GetEntryNr(); +#ifdef TPC_DIGIT_USEFAIRLINKS + mMCLinks.AddLink(FairLink(-1, eventID, "MCTrack", trackID)); +#else /// the MC ID is encoded such that we can have 999,999 tracks /// numbers larger than 1000000 correspond to the event ID /// i.e. 12000010 corresponds to event 12 with track ID 10 /// \todo Faster would be a bit shift -#ifdef TPC_DIGIT_USEFAIRLINKS - mMCLinks.AddLink(FairLink(-1, FairRootManager::Instance()->GetEntryNr(), "MCTrack", trackID)); + ++mMCID[(eventID)*1000000 + trackID]; #endif mChargePad += charge; } @@ -101,6 +125,8 @@ void DigitPad::reset() mChargePad = 0; #ifdef TPC_DIGIT_USEFAIRLINKS mMCLinks.Reset(); +#else + mMCID.clear(); #endif } diff --git a/Detectors/TPC/simulation/src/DigitPad.cxx b/Detectors/TPC/simulation/src/DigitPad.cxx index b61364e58d358..9b54ed0f7303b 100644 --- a/Detectors/TPC/simulation/src/DigitPad.cxx +++ b/Detectors/TPC/simulation/src/DigitPad.cxx @@ -21,12 +21,36 @@ void DigitPad::fillOutputContainer(TClonesArray *output, int cru, int timeBin, i const float mADC = SAMPAProcessing::makeSignal(totalADC, PadSecPos(CRU(cru).sector(), PadPos(row, pad))); if(mADC > 0) { + +#ifndef TPC_DIGIT_USEFAIRLINKS + static std::vector MClabels; + MClabels.resize(0); + DigitPad::processMClabels(MClabels); +#endif + TClonesArray &clref = *output; const size_t digiPos = clref.GetEntriesFast(); - DigitMC *digit = new(clref[digiPos]) DigitMC(cru, mADC, row, pad, timeBin, commonMode); + DigitMC *digit = new(clref[digiPos]) DigitMC( +#ifndef TPC_DIGIT_USEFAIRLINKS + MClabels, +#endif + cru, mADC, row, pad, timeBin, commonMode); #ifdef TPC_DIGIT_USEFAIRLINKS digit->SetLinks(getMCLinks()); #endif } } +#ifndef TPC_DIGIT_USEFAIRLINKS +void DigitPad::processMClabels(std::vector &sortedMCLabels) const +{ + /// Dump the map into a vector of pairs + std::vector > pairMClabels(mMCID.begin(), mMCID.end()); + /// Sort by the number of occurrences + std::sort(pairMClabels.begin(), pairMClabels.end(), boost::bind(&std::pair::second, _1) < boost::bind(&std::pair::second, _2)); + // iterate backwards over the vector and hence write MC with largest number of occurrences as first into the sortedMClabels vector + for(auto &aMCIDreversed : boost::adaptors::reverse(pairMClabels)) { + sortedMCLabels.emplace_back(aMCIDreversed.first); + } +} +#endif diff --git a/Detectors/TPC/simulation/src/DigitizerTask.cxx b/Detectors/TPC/simulation/src/DigitizerTask.cxx index 32421803e4714..76058e01bc47f 100644 --- a/Detectors/TPC/simulation/src/DigitizerTask.cxx +++ b/Detectors/TPC/simulation/src/DigitizerTask.cxx @@ -91,6 +91,7 @@ InitStatus DigitizerTask::Init() // Register output container mDigitsArray = new TClonesArray("o2::TPC::DigitMC"); + mDigitsArray->BypassStreamer(true); mgr->Register("TPCDigitMC", "TPC", mDigitsArray, kTRUE); mDigitizer->init(); @@ -115,6 +116,7 @@ void DigitizerTask::Exec(Option_t *option) if (mHitSector == -1){ // treat all sectors for (int s=0; s<18; ++s){ + LOG(DEBUG) << "Processing sector " << s << "\n"; mDigitContainer = mDigitizer->Process(mSectorHitsArray[s]); } } diff --git a/Detectors/TPC/simulation/src/TPCSimulationLinkDef.h b/Detectors/TPC/simulation/src/TPCSimulationLinkDef.h index bdd90ff50ff66..d6d62df512c98 100644 --- a/Detectors/TPC/simulation/src/TPCSimulationLinkDef.h +++ b/Detectors/TPC/simulation/src/TPCSimulationLinkDef.h @@ -36,6 +36,7 @@ #pragma link C++ class std::vector+; #pragma link C++ class o2::TPC::LinkableHitGroup+; #pragma link C++ class o2::TPC::SAMPAProcessing+; +#pragma link C++ class o2::TPC::TimeStamp+; #pragma link C++ class std::vector+; diff --git a/Detectors/TPC/simulation/test/testTPCSimulation.cxx b/Detectors/TPC/simulation/test/testTPCSimulation.cxx index a33f3ff5f7645..02b418697a728 100644 --- a/Detectors/TPC/simulation/test/testTPCSimulation.cxx +++ b/Detectors/TPC/simulation/test/testTPCSimulation.cxx @@ -33,6 +33,7 @@ namespace TPC { /// Precision: 1E-12 % BOOST_AUTO_TEST_CASE(DigitMC_test) { + /* DigitMC testdigit(1, 2.f, 3, 4, 5, 6.f); BOOST_CHECK(testdigit.getCRU() == 1); BOOST_CHECK_CLOSE(testdigit.getCharge(),2.f,1E-12); @@ -40,6 +41,7 @@ namespace TPC { BOOST_CHECK(testdigit.getPad() == 4); BOOST_CHECK_CLOSE(testdigit.getTimeStamp(), 5.f, 1E-12); BOOST_CHECK_CLOSE(testdigit.getCommonMode(),6.f,1E-12); + */ } } } diff --git a/macro/run_sim_tpc.C b/macro/run_sim_tpc.C index 01b8a055f7da4..ba6dcd9a3e1b7 100644 --- a/macro/run_sim_tpc.C +++ b/macro/run_sim_tpc.C @@ -19,10 +19,12 @@ #include "Field/MagneticField.h" #include "DetectorsPassive/Cave.h" - + #include "Generators/GeneratorFromFile.h" #include "TPCSimulation/Detector.h" #endif +//#define BOX_GENERATOR 1 + void run_sim_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3") { TString dir = getenv("VMCWORKDIR"); @@ -80,6 +82,7 @@ void run_sim_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3") // Create PrimaryGenerator FairPrimaryGenerator* primGen = new FairPrimaryGenerator(); +#ifdef BOX_GENERATOR FairBoxGenerator* boxGen = new FairBoxGenerator(211, 10); /*protons*/ //boxGen->SetThetaRange(0.0, 90.0); @@ -89,11 +92,17 @@ void run_sim_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3") boxGen->SetDebug(kTRUE); primGen->AddGenerator(boxGen); +#else + // reading the events from a kinematics file (produced by AliRoot) + auto extGen = new o2::eventgen::GeneratorFromFile("Kinematics.root"); + extGen->SetStartEvent(2); + primGen->AddGenerator(extGen); +#endif run->SetGenerator(primGen); // store track trajectories - run->SetStoreTraj(kTRUE); + // run->SetStoreTraj(kTRUE); // Initialize simulation run run->Init();