ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtTrackTransformer.cxx
Go to the documentation of this file.
1 #include "AtTrackTransformer.h"
2 // IWYU pragma: no_include <ext/alloc_traits.h>
3 
4 #include "AtHit.h" // for AtHit, AtHit::XYZPoint
5 #include "AtHitCluster.h" // for AtHitCluster
6 #include "AtTrack.h" // for XYZPoint, AtTrack
7 
8 #include <Math/Point3D.h> // for PositionVector3D, Cart...
9 #include <Math/Point3Dfwd.h>
10 #include <Math/Vector3D.h> // for DisplacementVector3D
11 #include <TMath.h> // for Power, Sqrt, ATan2, Pi
12 #include <TMatrixDSymfwd.h> // for TMatrixDSym
13 #include <TMatrixTSym.h> // for TMatrixTSym
14 #include <TVector3.h> // for TVector3
15 
16 #include <algorithm> // for max, for_each, copy_if
17 #include <iterator> // for back_insert_iterator
18 #include <memory> // for shared_ptr, __shared_p...
19 #include <vector> // for vector
20 
21 AtTools::AtTrackTransformer::AtTrackTransformer() = default;
22 AtTools::AtTrackTransformer::~AtTrackTransformer() = default;
23 using XYZPoint = ROOT::Math::XYZPoint;
24 
25 void AtTools::AtTrackTransformer::ClusterizeSmooth3D(AtTrack &track, Float_t distance, Float_t radius)
26 {
27  std::vector<AtHit> hitArray = track.GetHitArrayObject();
28  std::vector<AtHit> hitTBArray;
29  int clusterID = 0;
30 
31  // std::cout<<" ================================================================= "<<"\n";
32  // std::cout<<" Clusterizing track : "<<track.GetTrackID()<<"\n";
33 
34  /*for(auto iHits=0;iHits<hitArray->size();++iHits)
35  {
36  TVector3 pos = hitArray->at(iHits).GetPosition();
37  double Q = hitArray->at(iHits).GetCharge();
38  int TB = hitArray->at(iHits).GetTimeStamp();
39  //std::cout<<" Pos : "<<pos.X()<<" - "<<pos.Y()<<" - "<<pos.Z()<<" - TB : "<<TB<<" - Charge : "<<Q<<"\n";
40  }*/
41 
42  // Diffusion coefficients (TODO: Get them from the parameter file)
43  Double_t driftVel = 1.0; // cm/us
44  Double_t samplingRate = 0.320; // us
45  Double_t d_t = 0.0009; // cm^2/us
46  Double_t d_l = 0.0009; // cm^2/us
47  Double_t D_T = TMath::Sqrt((2.0 * d_t) / driftVel);
48  Double_t D_L = TMath::Sqrt((2.0 * d_l) / driftVel);
49 
50  if (hitArray.size() > 0) {
51 
52  auto refPos = hitArray.at(0).GetPosition(); // First hit
53  // TODO: Create a clustered hit from the very first hit (test)
54 
55  for (auto iHit = 0; iHit < hitArray.size(); ++iHit) {
56 
57  auto hit = hitArray.at(iHit);
58 
59  // Check distance with respect to reference Hit
60  Double_t distRef = TMath::Sqrt((hit.GetPosition() - refPos).Mag2());
61 
62  if (distRef < distance) {
63 
64  continue;
65 
66  } else {
67 
68  // std::cout<<" Clustering "<<iHit<<" of "<<hitArray->size()<<"\n";
69  // std::cout<<" Distance to reference : "<<distRef<<"\n";
70  // std::cout<<" Reference position : "<<refPos.X()<<" - "<<refPos.Y()<<" - "<<refPos.Z()<<" -
71  // "<<refPos.Mag()<<"\n";
72 
73  Double_t clusterQ = 0.0;
74  hitTBArray.clear();
75  std::copy_if(
76  hitArray.begin(), hitArray.end(), std::back_inserter(hitTBArray),
77  [&refPos, radius](AtHit &hitIn) { return TMath::Sqrt((hitIn.GetPosition() - refPos).Mag2()) < radius; });
78 
79  // std::cout<<" Clustered "<<hitTBArray.size()<<" Hits "<<"\n";
80 
81  if (hitTBArray.size() > 0) {
82  double x = 0, y = 0, z = 0;
83  double sigma_x = 0, sigma_y = 0, sigma_z = 0;
84 
85  int timeStamp;
86  std::shared_ptr<AtHitCluster> hitCluster = std::make_shared<AtHitCluster>();
87  hitCluster->SetClusterID(clusterID);
88  Double_t hitQ = 0.0;
89  std::for_each(hitTBArray.begin(), hitTBArray.end(),
90  [&x, &y, &z, &hitQ, &timeStamp, &sigma_x, &sigma_y, &sigma_z, &D_T, &D_L, &driftVel,
91  &samplingRate](AtHit &hitInQ) {
92  XYZPoint pos = hitInQ.GetPosition();
93  x += pos.X() * hitInQ.GetCharge();
94  y += pos.Y() * hitInQ.GetCharge();
95  z += pos.Z();
96  hitQ += hitInQ.GetCharge();
97  timeStamp += hitInQ.GetTimeStamp();
98 
99  // Calculation of variance (DOI: 10.1051/,00010 (2017)715001EPJ Web of
100  // Conferences50epjconf/2010010)
101  sigma_x += hitInQ.GetCharge() *
102  TMath::Sqrt(TMath::Power(0.2, 2) +
103  pos.Z() * TMath::Power(D_T, 2)); // 0.2 mm of position resolution
104  sigma_y += sigma_x;
105  sigma_z += TMath::Sqrt((1.0 / 6.0) * TMath::Power(driftVel * samplingRate, 2) +
106  pos.Z() * TMath::Power(D_L, 2));
107  });
108  x /= hitQ;
109  y /= hitQ;
110  z /= hitTBArray.size();
111  timeStamp /= hitTBArray.size();
112 
113  sigma_x /= hitQ;
114  sigma_y /= hitQ;
115  sigma_z /= hitTBArray.size();
116 
117  XYZPoint clustPos(x, y, z);
118  Bool_t checkDistance = kTRUE;
119 
120  // Check distance with respect to existing clusters
121  for (auto iClusterHit : *track.GetHitClusterArray()) {
122  if (TMath::Sqrt((iClusterHit.GetPosition() - clustPos).Mag2()) < distance) {
123  // std::cout<<" Cluster with less than : "<<distance<<" found "<<"\n";
124  checkDistance = kFALSE;
125  continue;
126  }
127  }
128 
129  if (checkDistance) {
130  hitCluster->SetCharge(hitQ);
131  hitCluster->SetPosition({x, y, z});
132  hitCluster->SetTimeStamp(timeStamp);
133  TMatrixDSym cov(3); // TODO: Setting covariant matrix based on pad size and drift time resolution.
134  // Using estimations for the moment.
135  cov(0, 1) = 0;
136  cov(1, 2) = 0;
137  cov(2, 0) = 0;
138  cov(0, 0) = TMath::Power(sigma_x, 2); // 0.04;
139  cov(1, 1) = TMath::Power(sigma_y, 2); // 0.04;
140  cov(2, 2) = TMath::Power(sigma_z, 2); // 0.01;
141  hitCluster->SetCovMatrix(cov);
142  ++clusterID;
143  track.AddClusterHit(hitCluster);
144  }
145  }
146  }
147 
148  // Sanity check
149  /*std::cout<<" Hits for cluster "<<iHit<<" centered in "<<refPos.X()<<" - "<<refPos.Y()<<"-"<<refPos.Z()<<"\n";
150  for(auto iHits=0;iHits<hitTBArray.size();++iHits)
151  {
152  TVector3 pos = hitTBArray.at(iHits).GetPosition();
153  double Q = hitTBArray.at(iHits).GetCharge();
154  int TB = hitTBArray.at(iHits).GetTimeStamp();
155  std::cout<<" Pos : "<<pos.X()<<" - "<<pos.Y()<<" - "<<pos.Z()<<" - TB : "<<TB<<" - Charge : "<<Q<<"\n";
156  std::cout<<" Distance to cluster center "<<TMath::Abs((track.GetHitClusterArray()->back().GetPosition() -
157  pos).Mag())<<"\n";
158  }
159  std::cout<<"=================================================="<<"\n";*/
160 
161  refPos = hitArray.at(iHit).GetPosition();
162 
163  //} // if distance
164 
165  } // for
166 
167  // Smoothing track
168  std::vector<AtHitCluster> *hitClusterArray = track.GetHitClusterArray();
169  radius /= 2.0;
170  std::vector<std::shared_ptr<AtHitCluster>> hitClusterBuffer;
171 
172  // std::cout<<" Hit cluster array size "<<hitClusterArray->size()<<"\n";
173 
174  if (hitClusterArray->size() > 2) {
175 
176  for (auto iHitCluster = 0; iHitCluster < hitClusterArray->size() - 1;
177  ++iHitCluster) // Calculating distances between pairs of clusters
178  {
179 
180  XYZPoint clusBack = hitClusterArray->at(iHitCluster).GetPosition();
181  XYZPoint clusForw = hitClusterArray->at(iHitCluster + 1).GetPosition();
182  XYZPoint clusMidPos = clusBack + (clusForw - clusBack) * 0.5;
183  std::vector<XYZPoint> renormClus{clusBack, clusMidPos};
184 
185  if (iHitCluster == (hitClusterArray->size() - 2))
186  renormClus.push_back(clusForw);
187 
188  // Create a new cluster and renormalize the charge of the other with half the radius.
189  for (auto iClus : renormClus) {
190  hitTBArray.clear();
191  std::copy_if(hitArray.begin(), hitArray.end(), std::back_inserter(hitTBArray),
192  [&iClus, radius](AtHit &hitIn) {
193  return TMath::Sqrt((hitIn.GetPosition() - iClus).Mag2()) < radius;
194  });
195 
196  if (hitTBArray.size() > 0) {
197  double x = 0, y = 0, z = 0;
198  double sigma_x = 0, sigma_y = 0, sigma_z = 0;
199 
200  int timeStamp = 0;
201  std::shared_ptr<AtHitCluster> hitCluster = std::make_shared<AtHitCluster>();
202  hitCluster->SetClusterID(clusterID);
203  Double_t hitQ = 0.0;
204  std::for_each(hitTBArray.begin(), hitTBArray.end(),
205  [&x, &y, &z, &hitQ, &timeStamp, &sigma_x, &sigma_y, &sigma_z, &D_T, &D_L, &driftVel,
206  &samplingRate](AtHit &hitInQ) {
207  auto pos = hitInQ.GetPosition();
208  x += pos.X() * hitInQ.GetCharge();
209  y += pos.Y() * hitInQ.GetCharge();
210  z += pos.Z();
211  hitQ += hitInQ.GetCharge();
212  timeStamp += hitInQ.GetTimeStamp();
213 
214  // Calculation of variance (DOI: 10.1051/,00010 (2017)715001EPJ Web of
215  // Conferences50epjconf/2010010)
216  sigma_x +=
217  hitInQ.GetCharge() *
218  TMath::Sqrt(TMath::Power(0.2, 2) +
219  pos.Z() * TMath::Power(D_T, 2)); // 0.2 mm of position resolution
220  sigma_y += sigma_x;
221  sigma_z += TMath::Sqrt((1.0 / 6.0) * TMath::Power(driftVel * samplingRate, 2) +
222  pos.Z() * TMath::Power(D_L, 2));
223  });
224  x /= hitQ;
225  y /= hitQ;
226  z /= hitTBArray.size();
227  timeStamp /= hitTBArray.size();
228 
229  sigma_x /= hitQ;
230  sigma_y /= hitQ;
231  sigma_z /= hitTBArray.size();
232 
233  TVector3 clustPos(x, y, z);
234  hitCluster->SetCharge(hitQ);
235  hitCluster->SetPosition({x, y, z});
236  hitCluster->SetTimeStamp(timeStamp);
237  TMatrixDSym cov(3); // TODO: Setting covariant matrix based on pad size and drift time resolution.
238  // Using estimations for the moment.
239  cov(0, 1) = 0;
240  cov(1, 2) = 0;
241  cov(2, 0) = 0;
242  cov(0, 0) = TMath::Power(sigma_x, 2); // 0.04;
243  cov(1, 1) = TMath::Power(sigma_y, 2); // 0.04;
244  cov(2, 2) = TMath::Power(sigma_z, 2); // 0.01;
245  hitCluster->SetCovMatrix(cov);
246  ++clusterID;
247  hitClusterBuffer.push_back(hitCluster);
248 
249  } // hitTBArray size
250 
251  } // for iClus
252 
253  } // for HitArray
254 
255  // Remove previous clusters
256  track.ResetHitClusterArray();
257 
258  // Adding new clusters
259  for (auto iHitClusterRe : hitClusterBuffer) {
260 
261  track.AddClusterHit(iHitClusterRe);
262  }
263 
264  } // Cluster array size
265 
266  } // if array size
267 }
AtHit::SetTimeStamp
void SetTimeStamp(Int_t Time)
Definition: AtHit.h:72
AtTools::AtTrackTransformer::ClusterizeSmooth3D
void ClusterizeSmooth3D(AtTrack &track, Float_t distance, Float_t radius)
Definition: AtTrackTransformer.cxx:25
AtTrack::ResetHitClusterArray
void ResetHitClusterArray()
Definition: AtTrack.h:111
AtHitCluster.h
XYZPoint
ROOT::Math::XYZPoint XYZPoint
Definition: AtTrackTransformer.cxx:23
AtTrack::GetHitArrayObject
std::vector< AtHit > GetHitArrayObject()
Definition: AtTrack.h:82
AtHitCluster::SetCovMatrix
void SetCovMatrix(TMatrixDSym matrix)
Definition: AtHitCluster.h:65
XYZPoint
ROOT::Math::XYZPoint XYZPoint
Definition: AtPatternCircle2D.h:18
AtHit.h
AtTrack::GetHitClusterArray
std::vector< AtHitCluster > * GetHitClusterArray()
Definition: AtTrack.h:90
AtHit::SetPosition
void SetPosition(const XYZPoint &pos)
Definition: AtHit.h:65
y
const double * y
Definition: lmcurve.cxx:20
AtTrack
Definition: AtTrack.h:25
AtTrack.h
AtTrack::AddClusterHit
void AddClusterHit(std::shared_ptr< AtHitCluster > hitCluster)
Definition: AtTrack.cxx:38
AtTools::AtTrackTransformer::AtTrackTransformer
AtTrackTransformer()
AtHit::SetCharge
void SetCharge(Double_t charge)
Definition: AtHit.h:63
AtTrackTransformer.h
AtTools::AtTrackTransformer::~AtTrackTransformer
~AtTrackTransformer()
AtHitCluster::SetClusterID
void SetClusterID(Int_t id)
Definition: AtHitCluster.h:69
AtHit
Point in space with charge.
Definition: AtHit.h:27
Generated on Sat Nov 11 2023 15:30:39 for ATTPCROOT by   doxygen 1.8.18