Skip to content

Regain possibility to do MC truth handling for TPC digits manually #385

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sawenzel merged 1 commit into from
May 19, 2017
Merged

Regain possibility to do MC truth handling for TPC digits manually #385

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 53 additions & 7 deletions Detectors/TPC/simulation/include/TPCSimulation/DigitMC.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#endif

#include "TPCBase/Digit.h"
#include "TPCSimulation/DigitPad.h"

#include "FairTimeStamp.h"

Expand All @@ -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
Expand All @@ -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<long> 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<int>(FairTimeStamp::GetTimeStamp()); };
int getTimeStamp() const final { return static_cast<int>(DigitBase::GetTimeStamp()); };

/// Get the common mode signal of the DigitMC
/// \return common mode signal of the DigitMC
Expand All @@ -60,25 +91,40 @@ class DigitMC : public FairTimeStamp, public Digit {
#ifndef __CINT__
friend class boost::serialization::access;
#endif

#ifndef TPC_DIGIT_USEFAIRLINKS
std::vector<long> 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);
};

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<long> 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

}
}
Expand Down
30 changes: 28 additions & 2 deletions Detectors/TPC/simulation/include/TPCSimulation/DigitPad.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<long> &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<long, int> 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<MCID> mMCID;
#endif
};

Expand All @@ -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];
Copy link
Copy Markdown
Member

@ktf ktf May 19, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not simply using a struct as the key?

struct MCID{
unsigned int eventId : 15; // 32k event Id possible
unsigned int trackId: 17; // 128K tracks possible
};

Also why using a map rather than a sorted vector with key and value in the same bitfield?

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<MCID> mMCID;

this has the added value you can use less bits for the occurrences and should be almost two orders of magnitude more compact.

Copy link
Copy Markdown
Collaborator Author

@sawenzel sawenzel May 19, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What you say is correct and it has been discussed in similar form among the TPC developers already. However, this commit merely "reverts" to what was already available before ... just to debug and compare to the FairLinks. I will add your comments as future ideas.

Copy link
Copy Markdown
Member

@ktf ktf May 19, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok.

#endif
mChargePad += charge;
}
Expand All @@ -101,6 +125,8 @@ void DigitPad::reset()
mChargePad = 0;
#ifdef TPC_DIGIT_USEFAIRLINKS
mMCLinks.Reset();
#else
mMCID.clear();
#endif
}

Expand Down
26 changes: 25 additions & 1 deletion Detectors/TPC/simulation/src/DigitPad.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<long> 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<long> &sortedMCLabels) const
{
/// Dump the map into a vector of pairs
std::vector<std::pair<long, int> > pairMClabels(mMCID.begin(), mMCID.end());
/// Sort by the number of occurrences
std::sort(pairMClabels.begin(), pairMClabels.end(), boost::bind(&std::pair<long, int>::second, _1) < boost::bind(&std::pair<long, int>::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
2 changes: 2 additions & 0 deletions Detectors/TPC/simulation/src/DigitizerTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand All @@ -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]);
}
}
Expand Down
1 change: 1 addition & 0 deletions Detectors/TPC/simulation/src/TPCSimulationLinkDef.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#pragma link C++ class std::vector<o2::TPC::ElementalHit>+;
#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<o2::TPC::Cluster>+;

Expand Down
2 changes: 2 additions & 0 deletions Detectors/TPC/simulation/test/testTPCSimulation.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,13 +33,15 @@ 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);
BOOST_CHECK(testdigit.getRow() == 3);
BOOST_CHECK(testdigit.getPad() == 4);
BOOST_CHECK_CLOSE(testdigit.getTimeStamp(), 5.f, 1E-12);
BOOST_CHECK_CLOSE(testdigit.getCommonMode(),6.f,1E-12);
*/
}
}
}
13 changes: 11 additions & 2 deletions macro/run_sim_tpc.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down Expand Up @@ -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);
Expand All @@ -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();
Expand Down