Skip to content

Commit c3677cd

Browse files
matthias-kleinermatthias-kleiner
authored
Merge pull request #9078 from matthias-kleiner/dedxdev
TPC toy cluster sim: small fixes
1 parent 1926689 commit c3677cd

1 file changed

Lines changed: 66 additions & 52 deletions

File tree

‎Detectors/TPC/simulation/macro/toyCluster.C‎

Lines changed: 66 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,6 @@
4141
#include "TRandom.h"
4242

4343
// O2 includes
44-
#include "DataFormatsTPC/TrackTPC.h"
4544
#include "TPCBase/Mapper.h"
4645
#include "TPCBase/ParameterDetector.h"
4746
#include "TPCBase/ParameterElectronics.h"
@@ -55,7 +54,6 @@
5554
#include "CommonUtils/TreeStreamRedirector.h"
5655
#include "CommonUtils/ConfigurableParam.h"
5756
#include "TPCBase/ParameterGas.h"
58-
#include "TPCSpaceCharge/PoissonSolverHelpers.h"
5957
#include "TPCReconstruction/HwClusterer.h"
6058
#include "TPCSimulation/GEMAmplification.h"
6159
#endif
@@ -66,7 +64,8 @@ GlobalPosition3D getPointFromPhi(float phi, float theta, GlobalPosition3D global
6664
GlobalPosition3D getPosBTrack(const GlobalPosition3D& posA, const GlobalPosition3D& posB);
6765
GlobalPosition3D getGlobalPositionTrk(float lambda, const GlobalPosition3D& refA, const GlobalPosition3D& refB);
6866

69-
const int mSector = 4; ///< consider only this mSector
67+
const int mSector = 4; ///< consider only this mSector
68+
const float mMaxDrift = 270; ///< maximum drift length in cm
7069

7170
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
7271
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -81,7 +80,7 @@ const int mSector = 4; ///< consider only this mSector
8180
/// \param dedxMax maximum dE/dx of the tracks (primary ionization: number of electrons per cm)
8281
/// \param maxSinPhi maximum sin(phi) of the tracks
8382
/// \param maxTanTheta maximum tan(theta) of the tracks
84-
void simulateTracks(const int maxEvents = 1, const char* outputFolder = "./", const float dedxMin = 14, const float dedxMax = 50, const float maxSinPhi = 1.f, const float maxTanTheta = 1.6)
83+
void simulateTracks(const int maxEvents = 1, const char* outputFolder = "./", const int dedxMin = 14, const int dedxMax = 50, const float maxSinPhi = 1.f, const float maxTanTheta = 1.6f)
8584
{
8685
// output file / tree
8786
TFile fOut(fmt::format("{}/o2sim_HitsTPC.root", outputFolder).data(), "RECREATE");
@@ -101,14 +100,17 @@ void simulateTracks(const int maxEvents = 1, const char* outputFolder = "./", co
101100
// set random seed
102101
gRandom->SetSeed(0);
103102

103+
// bias dedx to lower values
104+
TF1 fdEdx("fdEdx", "-log(x) + 10", 1, dedxMax);
105+
104106
// loop over events (tracks)
105107
for (int iEvent = 0; iEvent < maxEvents; ++iEvent) {
106108
hitGroupSector.clear();
107109

108110
// draw random track parameters
109111
phi = std::asin(gRandom->Uniform(0, maxSinPhi));
110112
theta = std::atan(gRandom->Uniform(0, maxTanTheta));
111-
dedx = gRandom->Uniform(dedxMin, dedxMax);
113+
dedx = (dedxMin == dedxMax) ? dedxMin : fdEdx.GetRandom();
112114

113115
// fill the tree (simulate the primary electrons along the track)
114116
fillTPCHits(theta, phi, dedx, hitGroupSector, trackInfo);
@@ -132,12 +134,16 @@ void fillTPCHits(const float theta, const float phi, const float dedx, std::vect
132134
const float xPos = gRandom->Uniform(-2, 2); // x starting position of track
133135
const int ireg = gRandom->Integer(Mapper::NREGIONS); // region where the track starts
134136
static auto& detParam = ParameterDetector::Instance();
135-
const float zCoordinate = gRandom->Uniform(0, detParam.TPClength - 1); // draw flat z position
137+
const float minZ = detParam.TPClength - mMaxDrift;
138+
static TF1 fz("fz", "tanh(x/60 - 3.5) * 0.6 + 1.6", minZ, detParam.TPClength);
139+
const float zCoordinate = fz.GetRandom(); // draw z position
140+
136141
const float flightTime = 0; // flight time of the particle. This value is not needed and set therefor to 0.
137142
const float radiusStart = mapper.getPadRegionInfo(ireg).getRadiusFirstRow(); // region start
138143
const auto padLength = mapper.getPadRegionInfo(ireg).getPadHeight();
139144
const float radius = radiusStart + Mapper::ROWSPERREGION[ireg] / 2 * padLength; // start in the center of the region
140145
HitGroup hitGroup{}; // create HitGroup and push_back to TPCHitsShiftedSector
146+
const float radiusEnd = radiusStart + Mapper::ROWSPERREGION[ireg] * padLength;
141147

142148
// first point of the track
143149
const GlobalPosition3D posATrack(xPos, radius, zCoordinate);
@@ -160,10 +166,10 @@ void fillTPCHits(const float theta, const float phi, const float dedx, std::vect
160166
const float radiusCurr = std::sqrt(xTmp * xTmp + yTmp * yTmp);
161167

162168
// check if the track is in the TPC
163-
const float rMinTPC = TPCParameters<float>::IFCRADIUS;
164-
const float rMaxTPC = TPCParameters<float>::OFCRADIUS;
165-
static auto& detParam = ParameterDetector::Instance();
166-
if (std::abs(zTmp) > detParam.TPClength || radiusCurr > rMaxTPC || yTmp < rMinTPC || zTmp < 0 || mapper.isOutOfSector(posTrk, Sector(mSector)) || radiusCurr < rMinTPC) {
169+
const float rMinTPC = radiusStart;
170+
const float rMaxTPC = radiusEnd;
171+
172+
if (std::abs(zTmp) > detParam.TPClength || radiusCurr > rMaxTPC || yTmp < rMinTPC || zTmp < minZ || mapper.isOutOfSector(posTrk, Sector(mSector)) || radiusCurr < rMinTPC) {
167173
continue;
168174
}
169175

@@ -259,6 +265,26 @@ void setZBinWidth(const int facZWidth = 1)
259265
LOG(info) << "zbinWidth: " << zbinWidth;
260266
}
261267

268+
GlobalPosition3D getElectronDrift(const GlobalPosition3D& posEle)
269+
{
270+
static o2::math_utils::RandomRing<> randomGaus;
271+
272+
const auto& detParam = ParameterDetector::Instance();
273+
const auto& gasParam = ParameterGas::Instance();
274+
float driftl = detParam.TPClength - posEle.Z();
275+
if (driftl < 0.01) {
276+
driftl = 0.01;
277+
}
278+
driftl = std::sqrt(driftl);
279+
const float sigT = driftl * gasParam.DiffT;
280+
const float sigL = driftl * gasParam.DiffL;
281+
282+
/// The position is smeared by a Gaussian with mean around the actual position and a width according to the diffusion
283+
/// coefficient times sqrt(drift length)
284+
GlobalPosition3D posEleDiffusion((randomGaus.getNextValue() * sigT) + posEle.X(), (randomGaus.getNextValue() * sigT) + posEle.Y(), (randomGaus.getNextValue() * sigL) + posEle.Z());
285+
return posEleDiffusion;
286+
}
287+
262288
/// creating the digits from the simulated hits (similar to the digitizer in O2)
263289
/// \param inpFileSim input sim file
264290
/// \param outName output path
@@ -269,6 +295,9 @@ void setZBinWidth(const int facZWidth = 1)
269295
/// \param disableNoise do not simulate any noise, commonMode and pedestal
270296
void createDigitsFromSim(const char* inpFileSim = "o2sim_HitsTPC.root", const std::string outName = "digits.root", const float eleAttachmentFac = 1, const int facZWidth = 1, const int maxSecondaries = 3, const int nEleGEM = -1, const bool disableNoise = false)
271297
{
298+
// set random seed
299+
gRandom->SetSeed(0);
300+
272301
auto& cdb = CDBInterface::instance();
273302
cdb.setUseDefaults();
274303

@@ -288,8 +317,8 @@ void createDigitsFromSim(const char* inpFileSim = "o2sim_HitsTPC.root", const st
288317
auto& detParam = ParameterDetector::Instance();
289318
auto& eleParam = ParameterElectronics::Instance();
290319
const float zbinWidth = eleParam.ZbinWidth;
291-
auto& gemParam = ParameterGEM::Instance();
292320
static GEMAmplification& gemAmplification = GEMAmplification::instance();
321+
const auto& gasParam = ParameterGas::Instance();
293322

294323
ElectronTransport& electronTransport = ElectronTransport::instance();
295324
electronTransport.updateParameters();
@@ -326,10 +355,9 @@ void createDigitsFromSim(const char* inpFileSim = "o2sim_HitsTPC.root", const st
326355
// loop over hits
327356
const int nEvents = hitTree->GetEntries(); // number of simulated events per hit file
328357
for (int iev = 0; iev < nEvents; ++iev) {
329-
const double eventTime = 100 + iev * 10; // some random time
358+
const double eventTime = 0;
330359
if (iev % 50 == false) {
331360
LOG(info) << "event: " << iev + 1 << " from " << nEvents << " events";
332-
LOG(info) << "TPC: Event time " << eventTime << " us";
333361
}
334362
hitTree->GetEntry(iev);
335363

@@ -347,17 +375,16 @@ void createDigitsFromSim(const char* inpFileSim = "o2sim_HitsTPC.root", const st
347375
GlobalPosition3D posEle(eh.GetX(), eh.GetY(), eh.GetZ());
348376

349377
const int nPrimaryElectrons = static_cast<int>(eh.GetEnergyLoss());
350-
const float hitTime = eh.GetTime() * 0.001;
378+
const float hitTime = eh.GetTime() * 0.001f;
351379
if (nPrimaryElectrons <= 0) {
352380
continue;
353381
}
354382

355-
float driftTime = 0.f;
356383
for (int iele = 0; iele < nPrimaryElectrons; iele++) {
357-
const GlobalPosition3D posEleDiff = electronTransport.getElectronDrift(posEle, driftTime);
384+
const GlobalPosition3D posEleDiff = getElectronDrift(posEle);
358385

359386
// add secondaries
360-
int nSecondaries = maxSecondaries == 0 ? 0 : getNsec();
387+
int nSecondaries = (maxSecondaries == 0) ? 0 : getNsec();
361388

362389
// restrict secondaries to avoid clusters with high charge (tail in the qMax qTot distributions)
363390
if (nSecondaries > maxSecondaries) {
@@ -374,24 +401,20 @@ void createDigitsFromSim(const char* inpFileSim = "o2sim_HitsTPC.root", const st
374401
// electrons are created randomly
375402
const double x = gRandom->Gaus(0, std::abs(posEle.X() - posEleDiff.X())) + posEle.X();
376403
const double y = gRandom->Gaus(0, std::abs(posEle.Y() - posEleDiff.Y())) + posEle.Y();
377-
double z = gRandom->Gaus(0, std::abs(posEle.Z() - posEleDiff.Z())) + posEle.Z();
378-
if (z < 0) {
379-
z = 0.01;
380-
}
404+
const double z = gRandom->Gaus(0, std::abs(posEle.Z() - posEleDiff.Z())) + posEle.Z();
381405
posTotElectrons.emplace_back(GlobalPosition3D(x, y, z));
382406
}
383407

384408
// loop over all electrons
385409
for (unsigned int j = 0; j < posTotElectrons.size(); ++j) {
386410
auto posEleTmp = posTotElectrons[j];
387-
driftTime = electronTransport.getDriftTime(posEleTmp.Z());
411+
const float driftTime = (detParam.TPClength - posEleTmp.Z()) / gasParam.DriftV;
388412

389413
const float eleTime = driftTime + hitTime; /// in us
390414
if (eleTime > maxEleTime) {
391415
LOG(warning) << "Skipping electron with driftTime " << driftTime << " from hit at time " << hitTime;
392416
continue;
393417
}
394-
const float absoluteTime = eleTime + eventTime; /// in us
395418

396419
// Attachment
397420
if (electronTransport.isElectronAttachment(driftTime)) {
@@ -404,11 +427,14 @@ void createDigitsFromSim(const char* inpFileSim = "o2sim_HitsTPC.root", const st
404427
}
405428

406429
// When the electron is not in the mSector we're processing, abandon
407-
if (mapper.isOutOfSector(posEleTmp, mSector)) {
430+
// create dummy pos at A-Side
431+
auto posEleTmpTmp = posEleTmp;
432+
posEleTmpTmp.SetZ(1);
433+
if (mapper.isOutOfSector(posEleTmpTmp, mSector)) {
408434
continue;
409435
}
410436

411-
const DigitPos digiPadPos = mapper.findDigitPosFromGlobalPosition(posEleTmp, mSector);
437+
const DigitPos digiPadPos = mapper.findDigitPosFromGlobalPosition(posEleTmpTmp, mSector);
412438
if (!digiPadPos.isValid()) {
413439
continue;
414440
}
@@ -432,11 +458,9 @@ void createDigitsFromSim(const char* inpFileSim = "o2sim_HitsTPC.root", const st
432458
const int sourceID = 0; // TPC
433459
const o2::MCCompLabel label(MCTrackID, eventID, sourceID, false);
434460

435-
for (float i = 0; i < nShapedPoints; ++i) {
436-
const float time = absoluteTime + i * zbinWidth;
437-
const auto timebin = sampaProcessing.getTimeBinFromTime(time);
438-
float fillCharge = signalArray[i];
439-
digitContainer.addDigit(label, cru, timebin, globalPad, fillCharge);
461+
for (int i = 0; i < nShapedPoints; ++i) {
462+
const float timebin = driftTime / eleParam.ZbinWidth + i;
463+
digitContainer.addDigit(label, cru, timebin, globalPad, signalArray[i]);
440464
}
441465
}
442466
} // electron loop
@@ -501,10 +525,13 @@ void createCluster(const char* inpHits = "o2sim_HitsTPC.root", const char* outFi
501525
{
502526
gRandom->SetSeed(0);
503527

504-
static auto& detParam = ParameterDetector::Instance();
505528
auto& eleParam = ParameterElectronics::Instance();
506529
auto& gasParam = ParameterGas::Instance();
530+
auto& cdb = CDBInterface::instance();
531+
cdb.setUseDefaults();
532+
507533
const static Mapper& mapper = Mapper::instance();
534+
SAMPAProcessing& sampaProcessing = SAMPAProcessing::instance();
508535

509536
// load the theta, phi and dE/dx angles from the simulation
510537
TFile fTrackInf(inpHits, "READ");
@@ -539,13 +566,6 @@ void createCluster(const char* inpHits = "o2sim_HitsTPC.root", const char* outFi
539566
std::vector<o2::tpc::Digit>* digits = new std::vector<o2::tpc::Digit>;
540567
tDigi->SetBranchAddress(fmt::format("TPCDigit_{}", mSector).data(), &digits);
541568

542-
// cut on the sigmaTime (extracted from plotting sigmaTime vs tan(theta))
543-
TF1 fSigmaTimeCut("sigmaTimeCut", "pol3", 0, 3);
544-
fSigmaTimeCut.SetParameter(0, 0.194862);
545-
fSigmaTimeCut.SetParameter(1, 0.818718);
546-
fSigmaTimeCut.SetParameter(2, -0.213279);
547-
fSigmaTimeCut.SetParameter(3, 0.0181574);
548-
549569
std::vector<float> vdedx;
550570
std::vector<float> vrelTime;
551571
std::vector<float> vpad;
@@ -619,27 +639,25 @@ void createCluster(const char* inpHits = "o2sim_HitsTPC.root", const char* outFi
619639
for (auto cont : clusterOutput) {
620640
auto container = cont.getContainer();
621641
const CRU cru(container->CRU);
622-
const PadRegionInfo& region = mapper.getPadRegionInfo(cru.region());
623-
const int rowOffset = region.getGlobalRowOffset();
642+
const int rowOffset = mapper.getPadRegionInfo(cru.region()).getGlobalRowOffset();
624643

625644
for (int clusterCount = 0; clusterCount < container->numberOfClusters; ++clusterCount) {
626645
const auto timeBinOffset = container->timeBinOffset;
627646
auto& cluster = container->clusters[clusterCount];
628647
const int qTot = cluster.getQTot();
629648
const int qMax = cluster.getQMax();
630-
const float pad = cluster.getPad() + 0.5;
649+
const float pad = cluster.getPad() + 0.5f;
631650
const float time = cluster.getTimeLocal() + timeBinOffset;
632651
const int padrow = rowOffset + cluster.getRow();
633652
const int region = cru.region();
634653
const float sigmaPad = std::sqrt(cluster.getSigmaPad2());
635654
const float sigmaTime = std::sqrt(cluster.getSigmaTime2());
636-
const float zPos = std::abs(time * eleParam.ZbinWidth * gasParam.DriftV - 0.6 - detParam.TPClength); // factor 0.6? see distToTrackZ
637-
const float relPad = pad - static_cast<int>(pad);
638-
const float relTime = time - static_cast<int>(time + 0.5);
639-
const auto tanTheta = std::tan(theta);
655+
const float zPos = sampaProcessing.getZfromTimeBin(time, Side::A) + eleParam.ZbinWidth * gasParam.DriftV;
656+
const float relPad = cluster.getPad() - static_cast<int>(pad);
657+
const float relTime = time - static_cast<int>(time + 0.5f);
640658

641659
// check for mSector edge pad
642-
const int off = 1;
660+
const int off = 2;
643661
const int offPad = 2;
644662
const int localPadRow = Mapper::getLocalRowFromGlobalRow(padrow);
645663
bool isEdge = false;
@@ -687,8 +705,7 @@ void createCluster(const char* inpHits = "o2sim_HitsTPC.root", const char* outFi
687705
const float distXT = (globX * xTrkF + globY * yTrkF) / (xTrkF * xTrkF + yTrkF * yTrkF);
688706
const float distanceToTrackXY = std::sqrt(std::pow(globX - distXT * xTrkF, 2) + std::pow(globY - distXT * yTrkF, 2));
689707

690-
const float distZ = (globZ * zTrkF) / (zTrkF * zTrkF);
691-
const float distanceToTrackZ = std::sqrt(std::pow(globZ - dist * zTrkF, 2));
708+
const float distanceToTrackZ = globZ - dist * zTrkF;
692709

693710
vLargestqTotinRow.emplace_back(largestqTot);
694711
vdedx.emplace_back(dedx);
@@ -711,8 +728,6 @@ void createCluster(const char* inpHits = "o2sim_HitsTPC.root", const char* outFi
711728
vsinglePadOrTime.emplace_back(singlePadOrTime);
712729
vLocalRow.emplace_back(localPadRow);
713730
zeroSuppOut.emplace_back(zeroSuppression);
714-
const bool lowTimeCut = sigmaTime < fSigmaTimeCut.Eval(tanTheta);
715-
vSigmaTimeCut.emplace_back(lowTimeCut);
716731
}
717732
}
718733

@@ -735,7 +750,6 @@ void createCluster(const char* inpHits = "o2sim_HitsTPC.root", const char* outFi
735750
<< "z=" << vzPos // z of the cluster
736751
<< "isEdge=" << visEdge // true if the cluster is edge pad
737752
<< "singlePadOrTime=" << vsinglePadOrTime // true if the cluster is single pad or single time cluster
738-
<< "sigmaTimeCut=" << vSigmaTimeCut // cut on the sigma time
739753
<< "eleAttFac=" << eleAttFac // electron attachement factor
740754
<< "zeroSupp=" << zeroSuppOut // absolute zero supression value in ADC counts
741755
<< "isLargestqTot=" << vLargestqTotinRow // is true cluster has the highest qTot from all clusters in the same pad row (cluster with lower charge are noise)
@@ -746,7 +760,7 @@ void createCluster(const char* inpHits = "o2sim_HitsTPC.root", const char* outFi
746760
}
747761

748762
TTree* tree = (TTree*)(pcstream.GetFile()->Get("cl"));
749-
tree->SetAlias("cut", "isLargestqTot==1 && singlePadOrTime==0 && isEdge==0 && sigmaTimeCut==0 && sigmaTime>0.5"); // cut to filter clusters
763+
tree->SetAlias("cut", "isLargestqTot==1 && singlePadOrTime==0 && isEdge==0"); // cut to filter clusters
750764
pcstream.Close();
751765
fTrackInf.Close();
752766
}

0 commit comments

Comments
 (0)