Skip to content
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

Add Geant track info (g2t) to StFcsHit #379

Merged
merged 8 commits into from
Sep 6, 2022
Merged

Add Geant track info (g2t) to StFcsHit #379

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
868e811
Adding g2t track info to StFcsHit, and adding code to trace back to p…
akioogawa Jul 27, 2022
bcdf46f
Update based on feedbacks on PR#379
akioogawa Jul 29, 2022
8b45edf
Change StFcsHit::addGeantTrack to do reverse search for speed
akioogawa Jul 29, 2022
d853257
StFcsHit::addGeantTrack()
akioogawa Aug 26, 2022
dc2c660
Making a private getG2tTrack() to avoid duplicate logic/code.
akioogawa Aug 27, 2022
b9f423d
Remove sorting of geant tracks in StFcsFastSimulatorMaker.cxx
akioogawa Aug 29, 2022
824b3b6
Merge branch 'main' into fcsMc
akioogawa Aug 31, 2022
33b6115
Changes in white space only
plexoos Sep 6, 2022
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
3 changes: 3 additions & 0 deletions StRoot/StEvent/StEventLinkDef.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -238,5 +238,8 @@
#pragma link C++ function operator<<(ostream&, const StFmsHit&);

#pragma link C++ class vector<StFmsPointPair*>+;
#pragma link C++ class pair<unsigned int,float>+;
#pragma link C++ class vector<pair<unsigned int,float>>+;

#endif

19 changes: 19 additions & 0 deletions StRoot/StEvent/StFcsHit.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,10 @@
*
**************************************************************************/
#include "StFcsHit.h"
#include <algorithm>

ClassImp(StFcsHit)
using namespace std;

StFcsHit::StFcsHit() { /* no operation */}

Expand Down Expand Up @@ -139,3 +141,20 @@ void StFcsHit::print(Option_t *option) const {
}
cout << endl;
}

void StFcsHit::addGeantTrack(unsigned int id, float e){
unsigned int n=mGeantTracks.size();
for(unsigned int i=0; i<n; i++){
if(mGeantTracks[i].first == id) {mGeantTracks[i].second += e; return;}
}
plexoos marked this conversation as resolved.
Show resolved Hide resolved
mGeantTracks.push_back(make_pair(id,e));
return;
}

void StFcsHit::sortGeantTracks(){ //Sort Geant Track by descending order of energy
std::sort(mGeantTracks.begin(), mGeantTracks.end(),
[](const pair<unsigned int,float>&a, const pair<unsigned int,float>&b){
return b.second < a.second;
});
}

11 changes: 8 additions & 3 deletions StRoot/StEvent/StFcsHit.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,9 +87,13 @@ class StFcsHit : public StObject {
void setCluster(StFcsCluster* clu) {mCluster = clu;}
StFcsCluster *cluster() {return mCluster;}

const vector<pair<unsigned int, float>>& getGeantTracks() const {return mGeantTracks;}
Copy link
Member

Choose a reason for hiding this comment

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

Do we really need this information? It's not complete enough to learn something about the shower, but too complicated to learn about what caused the hit. I don't see why the analyzers couldn't get by with just unsigned int idTruth() and float qaTruth() (aka fraction in your current PR). This should be enough for two-particle studies (energy fraction that is not associated with the leading particle must be coming from the other one). I suppose, unlike for TPC, you also wanted to store the particle that arrived to the detector (getParentG2tTrack) in the current PR, but again, this is not equivalent to storing energy and id of every hit.

I'm mostly asking because I'm pretty sure the interface from StEvent is going to migrate to MuDst without a change.

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes, it is needed. A single ID truth is insufficient to understand the relative contributions of electromagnetic and hadronic components to the cluster energy response (You don't know what the "rest" of the energy is...)

With the list of track IDs hitting each cluster, their contribution to the energy deposition in the cluster, and the MC track's kinematics from the MC event record, you have the complete information needed to understand the clustering performance.

Copy link
Member

Choose a reason for hiding this comment

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

Yes, it is needed. A single ID truth is insufficient to understand the relative contributions of electromagnetic and hadronic components to the cluster energy response (You don't know what the "rest" of the energy is...)

Okay, but if that is the case, I need to also save PDG id. If not, I need to look into geant.root, and at this point I might as well not save anything to event.root. Right?

With the list of track IDs hitting each cluster, their contribution to the energy deposition in the cluster, and the MC track's kinematics from the MC event record, you have the complete information needed to understand the clustering performance.

Again, for any foreseeable study, there are just two particles that could merge into the same cluster, one fraction per cell describes that case it well enough.

Copy link
Contributor

Choose a reason for hiding this comment

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

Okay, but if that is the case, I need to also save PDG id. If not, I need to look into geant.root, and at this point I might as well not save anything to event.root. Right?

PDG ID is available in the MuDst MuMcTrack branch.

Again, for any foreseeable study, there are just two particles that could merge into the same cluster, one fraction per cell describes that case it well enough.

Relative contribution of hadronic and electromagnetic particles to the energy of a cluster is one such study which requires this information.

Copy link
Member

Choose a reason for hiding this comment

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

PDG ID is available in the MuDst MuMcTrack branch.

I'm not very familiar with that facility. And this answer also implies that this structure will be propagated to MuDst.

Again, for any foreseeable study, there are just two particles that could merge into the same cluster, one fraction per cell describes that case it well enough.

Relative contribution of hadronic and electromagnetic particles to the energy of a cluster is one such study which requires this information.

Okay, so that would be a study that looks at cluster merging and hadronic/electromagnetic fractions at the same time. Thank you for the explanations. I hope other experts/analyzers can chime in.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

  • Yes, idTruth is not enough, especially when many towers contributing to a cluster, and hadron showers are fat and big.
  • Yes, when we get to mudst, since Mudst has StMcTrack, but not g2t_hits.
  • No, when reading event.root file and geant.root file in sync. But this means we essentially re-run StFcsFastSimulationMake and creating new StFcsHits. (for historical reason, we don't have StMcTrack in StEvent).

void addGeantTrack(unsigned int id, float e);
void sortGeantTracks();
plexoos marked this conversation as resolved.
Show resolved Hide resolved

void print(Option_t *option="") const;

protected:
protected:
UShort_t mDetId=0; // 1 bit ZS, 3 bits DetectorId, 12 bits id
UShort_t mDepCh=0; // 1 bit for NS, 2 bits for EHP, 5 bits for DEP, 8 bits for channal
UInt_t mAdcSum=0; // ADC sum
Expand All @@ -100,9 +104,10 @@ class StFcsHit : public StObject {
Float_t mEnergy=0.0; // corrected energy
StFcsCluster* mCluster=0; // pointer to cluster this hit belongs
TArrayS* mData=0; // 12bit ADC values + flag at highest 4 bits, array of timebin

ClassDef(StFcsHit,5)

vector<pair<unsigned int, float>> mGeantTracks; // parent G2T track id and dE
akioogawa marked this conversation as resolved.
Show resolved Hide resolved

ClassDef(StFcsHit,6)
};

#endif
43 changes: 43 additions & 0 deletions StRoot/StFcsClusterMaker/StFcsClusterMaker.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,7 @@
#include "StEvent/StFcsHit.h"
#include "StEvent/StFcsCluster.h"
#include "StFcsDbMaker/StFcsDb.h"
#include "tables/St_g2t_track_Table.h"

#include "StMuDSTMaker/COMMON/StMuTypes.hh"
#include "StMuDSTMaker/COMMON/StMuDst.h"
Expand Down Expand Up @@ -217,6 +218,48 @@ int StFcsClusterMaker::makeCluster(int det) {
if(ret==kStErr) ret=clusterMomentAnalysis(clu,0.0); //Redo with 0 threshold
categorization(clu);
}

//debug MC info
if(GetDebug()>=5){
g2t_track_st* g2ttrk=0;
St_g2t_track* trackTable = static_cast<St_g2t_track*>(GetDataSet("g2t_track"));
if(!trackTable) {
LOG_INFO << "g2t_track Table not found" << endm;
}else{
const int nTrk = trackTable->GetNRows();
LOG_INFO << Form("g2t_track table has %d tracks",nTrk) << endm;
if(nTrk>0){
g2ttrk = trackTable->GetTable();
if(!g2ttrk){
LOG_INFO << "g2t_track GetTable failed" << endm;
}
}
}
if(g2ttrk){
float frc=0;
int nh = hits.size();
for(int i=0; i<nh; i++){
StFcsHit* hit=hits[i];
g2t_track_st* trk = mDb->getParentG2tTrack(hit,g2ttrk,frc);
printf("Det=%1d Id=%3d E=%8.3f Parent Id=%4d Pid=%4d E=%8.3f Frc=%6.3f\n",
klendathu2k marked this conversation as resolved.
Show resolved Hide resolved
det,hit->id(),hit->energy(),trk->id,trk->ge_pid,trk->e,frc);
g2t_track_st* ptrk = mDb->getPrimaryG2tTrack(hit,g2ttrk,frc);
printf("Det=%1d Id=%3d E=%8.3f Primary Id=%4d Pid=%4d E=%8.3f Frc=%6.3f\n",
det,hit->id(),hit->energy(),ptrk->id,ptrk->ge_pid,ptrk->e,frc);
}
int nc = clusters.size();
for(int j=0; j<nc; j++){
StFcsCluster* clu=clusters[j];
g2t_track_st* trk = mDb->getParentG2tTrack(clu,g2ttrk,frc);
printf("Det=%1d C#=%3d E=%8.3f Parent Id=%4d Pid=%4d E=%8.3f Frc=%6.3f\n",
det,j,clu->energy(),trk->id,trk->ge_pid,trk->e,frc);
g2t_track_st* ptrk = mDb->getPrimaryG2tTrack(clu,g2ttrk,frc);
printf("Det=%1d C#=%3d E=%8.3f Primary Id=%4d Pid=%4d E=%8.3f Frc=%6.3f\n",
det,j,clu->energy(),ptrk->id,ptrk->ge_pid,ptrk->e,frc);
}
}
}

return kStOk;
}

Expand Down
116 changes: 116 additions & 0 deletions StRoot/StFcsDbMaker/StFcsDb.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1979,3 +1979,119 @@ void StFcsDb::readGainCorrFromText(){
}
fclose(F);
}

//g2t track info
unsigned int StFcsDb::backTraceG2tTrack(unsigned int id, g2t_track_st* g2ttrk){
int i = id - 1;
while(g2ttrk[i].next_parent_p !=0){
if(mDebug>3)
LOG_INFO<<Form(" BackTrace from=%3d id=%3d Epid=%4d Gpid=%3d Vtx=%3d Parent=%3d E=%6.2f",
id,g2ttrk[i].id,g2ttrk[i].eg_pid,g2ttrk[i].ge_pid,g2ttrk[i].start_vertex_p,
g2ttrk[i].next_parent_p,g2ttrk[i].e)<<endm;
i = g2ttrk[i].next_parent_p - 1;
}
if(mDebug>3) LOG_INFO<<Form(" BackTrace from=%3d id=%3d Epid=%4d Gpid=%3d Vtx=%3d Parent=%3d E=%6.2f Primary!!",
id,g2ttrk[i].id,g2ttrk[i].eg_pid,g2ttrk[i].ge_pid,g2ttrk[i].start_vertex_p,
g2ttrk[i].next_parent_p,g2ttrk[i].e)<<endm;
return i + 1;
}

g2t_track_st* StFcsDb::getParentG2tTrack(StFcsHit* h, g2t_track_st* g2ttrk, float& fraction, unsigned int order){
const vector<pair<unsigned int, float>>& gt = h->getGeantTracks();
plexoos marked this conversation as resolved.
Show resolved Hide resolved
int ntrk=gt.size();
float detot=0;
if(order > gt.size()) {fraction=0; return 0;}
for(int itrk=0; itrk<ntrk; itrk++) detot += gt[itrk].second;
fraction = gt[order].second / detot;
return &g2ttrk[gt[order].first-1];
}

g2t_track_st* StFcsDb::getPrimaryG2tTrack(StFcsHit* h, g2t_track_st* g2ttrk, float& fraction, unsigned int order){
const vector<pair<unsigned int, float>>& gt = h->getGeantTracks();
vector<pair<unsigned int, float>> primary;
int ntrk=gt.size();
float detot=0;
for(int itrk=0; itrk<ntrk; itrk++){
float de=gt[itrk].second;
unsigned int id=backTraceG2tTrack(gt[itrk].first, g2ttrk);
int found=0;
unsigned int np=primary.size();
for(unsigned int jtrk=0; jtrk<np; jtrk++){
if(primary[jtrk].first == id) {primary[jtrk].second += de; found=1; break;}
}
if(found==0) primary.push_back(make_pair(id,de));
plexoos marked this conversation as resolved.
Show resolved Hide resolved
detot+=de;
}
if(order > primary.size()) {fraction=0; return 0;}
std::sort(primary.begin(), primary.end(),
[](const pair<unsigned int,float>&a, const pair<unsigned int,float>&b){
return b.second < a.second;
});
plexoos marked this conversation as resolved.
Show resolved Hide resolved
if(order > primary.size()) {fraction=0; return 0;}
fraction = primary[order].second / detot;
return &g2ttrk[primary[order].first-1];
}

g2t_track_st* StFcsDb::getParentG2tTrack(StFcsCluster* c, g2t_track_st* g2ttrk, float& fraction, unsigned int order){
float detot=0;
vector<pair<unsigned int, float>> parents;
StPtrVecFcsHit& hits = c->hits();
int nhit = hits.size();
for(int ihit=0; ihit<nhit; ihit++){
const vector<pair<unsigned int, float>>& gt = hits[ihit]->getGeantTracks();
plexoos marked this conversation as resolved.
Show resolved Hide resolved
int ntrk=gt.size();
for(int itrk=0; itrk<ntrk; itrk++){
unsigned int id=gt[itrk].first;
float de=gt[itrk].second;
int found=0;
unsigned int np=parents.size();
for(unsigned int jtrk=0; jtrk<np; jtrk++){
if(parents[jtrk].first == id) {parents[jtrk].second += de; found=1; break;}
}
if(found==0) parents.push_back(make_pair(id,de));
detot+=de;
}
}
if(order > parents.size()) {fraction=0; return 0;}
std::sort(parents.begin(), parents.end(),
[](const pair<unsigned int,float>&a, const pair<unsigned int,float>&b){
return b.second < a.second;
});
fraction = parents[order].second / detot;
return &g2ttrk[parents[order].first-1];
}

g2t_track_st* StFcsDb::getPrimaryG2tTrack(StFcsCluster* c, g2t_track_st* g2ttrk, float& fraction, unsigned int order){
float detot=0;
vector<pair<unsigned int, float>> primary;
StPtrVecFcsHit& hits = c->hits();
int nhit = hits.size();
for(int ihit=0; ihit<nhit; ihit++){
const vector<pair<unsigned int, float>>& gt = hits[ihit]->getGeantTracks();
int ntrk=gt.size();
for(int itrk=0; itrk<ntrk; itrk++){
unsigned int id=backTraceG2tTrack(gt[itrk].first,g2ttrk);
float de=gt[itrk].second;
int found=0;
unsigned int np=primary.size();
for(unsigned int jtrk=0; jtrk<np; jtrk++){
if(primary[jtrk].first == id) {primary[jtrk].second += de; found=1; break;}
}
if(found==0) primary.push_back(make_pair(id,de));
detot+=de;
}
}
if(order > primary.size()) {fraction=0; return 0;}
std::sort(primary.begin(), primary.end(),
[](const pair<unsigned int,float>&a, const pair<unsigned int,float>&b){
return b.second < a.second;
});
if(mDebug>3){
unsigned int np=primary.size();
for(unsigned int jtrk=0; jtrk<np; jtrk++){
LOG_INFO << Form("Sorted Primary G2T Track %3d id=%3d dE=%f",jtrk,primary[jtrk].first,primary[jtrk].second)<<endm;
}
}
fraction = primary[order].second / detot;
return &g2ttrk[primary[order].first-1];
}
8 changes: 8 additions & 0 deletions StRoot/StFcsDbMaker/StFcsDb.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,7 @@
#include "tables/St_fcsHcalGainCorr_Table.h"
#include "tables/St_fcsPresValley_Table.h"
#include "tables/St_vertexSeed_Table.h"
#include "tables/St_g2t_track_Table.h"
class StFcsHit;
class StFcsCluster;
class StFcsPoint;
Expand Down Expand Up @@ -256,6 +257,13 @@ class StFcsDb : public TDataSet {
void setPedestal(int ehp, int ns, int dep, int ch, float ped); //! setting pedestal
void readPedFromText(const char* file="fcsped.txt"); //! reading pedestal from text

//g2t track info
unsigned int backTraceG2tTrack(unsigned int id, g2t_track_st* g2ttrk);
g2t_track_st* getParentG2tTrack(StFcsHit* h, g2t_track_st* g2ttrk, float& fraction, unsigned int order=0);
g2t_track_st* getParentG2tTrack(StFcsCluster* c, g2t_track_st* g2ttrk, float& fraction, unsigned int order=0);
g2t_track_st* getPrimaryG2tTrack(StFcsHit* h, g2t_track_st* g2ttrk, float& fraction, unsigned int order=0);
g2t_track_st* getPrimaryG2tTrack(StFcsCluster* c, g2t_track_st* g2ttrk, float& fraction, unsigned int order=0);
akioogawa marked this conversation as resolved.
Show resolved Hide resolved

private:
int mDbAccess=1; //! enable(1) or disabe(0) DB access
int mRun=0; //! run#
Expand Down
5 changes: 5 additions & 0 deletions StRoot/StFcsFastSimulatorMaker/StFcsFastSimulatorMaker.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,7 @@ void StFcsFastSimulatorMaker::fillStEvent(StEvent* event) {
fcshit = mEcalMap[ns][id];
fcshit->setEnergy(fcshit->energy() + de);
}
fcshit->addGeantTrack(hit->track_p,de);
hit++;
}
}
Expand Down Expand Up @@ -200,6 +201,7 @@ void StFcsFastSimulatorMaker::fillStEvent(StEvent* event) {
fcshit = mHcalMap[ns][id];
fcshit->setEnergy(fcshit->energy() + de);
}
fcshit->addGeantTrack(hit->track_p,de);
hit++;
}else{ //leaky hcal with up to 4 WLSP getting lights from a tower
float de[4];
Expand Down Expand Up @@ -234,6 +236,7 @@ void StFcsFastSimulatorMaker::fillStEvent(StEvent* event) {
fcshit = mHcalMap[ns][id2];
fcshit->setEnergy(fcshit->energy() + de[j]);
}
fcshit->addGeantTrack(hit->track_p,de[j]);
}
hit++;
}
Expand Down Expand Up @@ -291,6 +294,7 @@ void StFcsFastSimulatorMaker::fillStEvent(StEvent* event) {
fcshit->setEnergy(fcshit->energy() + de);
}
hit++;
fcshit->addGeantTrack(hit->track_p,de);
}
}
}
Expand Down Expand Up @@ -318,6 +322,7 @@ void StFcsFastSimulatorMaker::fillStEvent(StEvent* event) {
int ehp = mFcsDb->ecalHcalPres(det);
hits[i]->setAdc(0,adc);
hits[i]->setEnergy(digi_energy);
hits[i]->sortGeantTracks();
plexoos marked this conversation as resolved.
Show resolved Hide resolved
fcscollection->addHit(det,hits[i]);
etot[ehp] += digi_energy;
nhit[ehp]++;
Expand Down