CbmRoot
CbmStsAlgoAnaCluster.cxx
Go to the documentation of this file.
1 
6 #include "CbmStsAlgoAnaCluster.h"
7 
8 #include "CbmDigiManager.h"
9 #include "CbmStsCluster.h"
10 #include "CbmStsDigi.h"
11 #include "CbmStsParModule.h"
12 #include "CbmStsPhysics.h"
13 #include <TMath.h>
14 
15 using std::unique_ptr;
16 
17 
19 
20 
21  // ----- Constructor ----------------------------------------------------
23  : fDigiMan(CbmDigiManager::Instance()), fPhysics(CbmStsPhysics::Instance()) {}
24 // --------------------------------------------------------------------------
25 
26 
27 // ----- Algorithm for one-digi clusters --------------------------------
29  const CbmStsParModule* module) {
30 
31  assert(module);
32  Int_t index = cluster.GetDigi(0);
33  const CbmStsDigi* digi = fDigiMan->Get<CbmStsDigi>(index);
34  assert(digi);
35  auto& asic = module->GetParAsic(digi->GetChannel());
36  Double_t x = Double_t(digi->GetChannel());
37  Double_t time = digi->GetTime();
38  Double_t timeError = asic.GetTimeResol();
39  Double_t charge = asic.AdcToCharge(digi->GetCharge());
40  Double_t xError = 1. / sqrt(24.);
41 
42  cluster.SetProperties(charge, x, xError, time, timeError);
43  cluster.SetSize(1);
44 }
45 // --------------------------------------------------------------------------
46 
47 
48 // ----- Algorithm for two-digi clusters --------------------------------
50  const CbmStsParModule* modPar) {
51 
52  Int_t index1 = cluster.GetDigi(0);
53  Int_t index2 = cluster.GetDigi(1);
54  const CbmStsDigi* digi1 = fDigiMan->Get<CbmStsDigi>(index1);
55  const CbmStsDigi* digi2 = fDigiMan->Get<CbmStsDigi>(index2);
56  assert(digi1);
57  assert(digi2);
58 
59  auto& asic1 = modPar->GetParAsic(digi1->GetChannel());
60  auto& asic2 = modPar->GetParAsic(digi2->GetChannel());
61 
62  // --- Uncertainties of the charge measurements
63  Double_t eNoiseSq = 0.5
64  * (asic1.GetNoise() * asic1.GetNoise()
65  + asic2.GetNoise() * asic2.GetNoise());
66  Double_t chargePerAdc =
67  0.5
68  * (asic1.GetDynRange() / Double_t(asic1.GetNofAdc())
69  + asic2.GetDynRange() / Double_t(asic2.GetNofAdc()));
70  Double_t eDigitSq = chargePerAdc * chargePerAdc / 12.;
71 
72  UInt_t chan1 = digi1->GetChannel();
73  UInt_t chan2 = digi2->GetChannel();
74  assert(chan2 == chan1 + 1
75  || chan2 == chan1 - modPar->GetNofChannels() / 2 + 1);
76 
77  // Channel positions and charge
78  Double_t x1 = Double_t(chan1);
79  Double_t q1 = asic1.AdcToCharge(digi1->GetCharge());
80  Double_t q2 = asic2.AdcToCharge(digi2->GetCharge());
81 
82  // Periodic position for clusters round the edge
83  if (chan1 > chan2) x1 -= Double_t(modPar->GetNofChannels() / 2);
84 
85  // Uncertainties of the charge measurements
86  Double_t width1 = fPhysics->LandauWidth(q1);
87  Double_t eq1sq = width1 * width1 + eNoiseSq + eDigitSq;
88  Double_t width2 = fPhysics->LandauWidth(q2);
89  Double_t eq2sq = width2 * width2 + eNoiseSq + eDigitSq;
90 
91  // Cluster time
92  Double_t time = 0.5 * (digi1->GetTime() + digi2->GetTime());
93  Double_t timeError = 0.5 * (asic1.GetTimeResol() + asic2.GetTimeResol())
94  * 0.70710678; // 1/sqrt(2)
95 
96  // Cluster position
97  // See corresponding software note.
98  Double_t x = x1 + 0.5 + (q2 - q1) / 3. / TMath::Max(q1, q2);
99 
100  // Correct negative position for clusters around the edge
101  if (x < -0.5) x += Double_t(modPar->GetNofChannels() / 2);
102 
103  // Uncertainty on cluster position. See software note.
104  Double_t ex0sq = 0.; // error for ideal charge measurements
105  Double_t ex1sq = 0.; // error from first charge
106  Double_t ex2sq = 0.; // error from second charge
107  if (q1 < q2) {
108  ex0sq = (q2 - q1) * (q2 - q1) / q2 / q2 / 72.;
109  ex1sq = eq1sq / q2 / q2 / 9.;
110  ex2sq = eq2sq * q1 * q1 / q2 / q2 / q2 / q2 / 9.;
111  } else {
112  ex0sq = (q2 - q1) * (q2 - q1) / q1 / q1 / 72.;
113  ex1sq = eq1sq * q2 * q2 / q1 / q1 / q1 / q1 / 9.;
114  ex2sq = eq2sq / q1 / q1 / 9.;
115  }
116  Double_t xError = TMath::Sqrt(ex0sq + ex1sq + ex2sq);
117 
118  // Cluster charge
119  Double_t charge = q1 + q2;
120 
121  cluster.SetProperties(charge, x, xError, time, timeError);
122  cluster.SetSize(2);
123 }
124 // --------------------------------------------------------------------------
125 
126 
127 // ----- Algorithm for clusters with more than two digis ----------------
129  const CbmStsParModule* modPar) {
130 
131  Double_t tSum = 0.; // sum of digi times
132  Int_t chanF = 9999999; // first channel in cluster
133  Int_t chanL = -1; // last channel in cluster
134  Double_t qF = 0.; // charge in first channel
135  Double_t qM = 0.; // sum of charges in middle channels
136  Double_t qL = 0.; // charge in last cluster
137  Double_t eqFsq = 0.; // uncertainty of qF
138  Double_t eqMsq = 0.; // uncertainty of qMid
139  Double_t eqLsq = 0.; // uncertainty of qL
140  Double_t prevChannel = 0;
141  Double_t tResolSum = 0.;
142 
143  for (Int_t iDigi = 0; iDigi < cluster.GetNofDigis(); iDigi++) {
144 
145  Int_t index = cluster.GetDigi(iDigi);
146  const CbmStsDigi* digi = fDigiMan->Get<CbmStsDigi>(index);
147  assert(digi);
148  Int_t channel = digi->GetChannel();
149  auto& asic = modPar->GetParAsic(channel);
150 
151  // --- Uncertainties of the charge measurements
152  Double_t eNoiseSq = asic.GetNoise() * asic.GetNoise();
153  Double_t chargePerAdc = asic.GetDynRange() / Double_t(asic.GetNofAdc());
154  Double_t eDigitSq = chargePerAdc * chargePerAdc / 12.;
155  tResolSum += asic.GetTimeResol();
156 
157  tSum += digi->GetTime();
158  Double_t charge = asic.AdcToCharge(digi->GetCharge());
159  Double_t lWidth = fPhysics->LandauWidth(charge);
160  Double_t eChargeSq = lWidth * lWidth + eNoiseSq + eDigitSq;
161 
162  // Check ascending order of channel number
163  if (iDigi > 0)
164  assert(channel == prevChannel + 1
165  || channel == prevChannel - modPar->GetNofChannels() / 2 + 1);
166  prevChannel = channel;
167 
168  if (iDigi == 0) { // first channel
169  chanF = channel;
170  qF = charge;
171  eqFsq = eChargeSq;
172  } else if (iDigi == cluster.GetNofDigis() - 1) { // last channel
173  chanL = channel;
174  qL = charge;
175  eqLsq = eChargeSq;
176  } else { // one of the middle channels
177  qM += charge;
178  eqMsq += eChargeSq;
179  }
180 
181  } //# digis in cluster
182 
183  // Periodic channel position for clusters round the edge
184  if (chanF > chanL) chanF -= modPar->GetNofChannels() / 2;
185 
186  // Cluster time and total charge
187  tSum = tSum / Double_t(cluster.GetNofDigis());
188  Double_t tError = (tResolSum / Double_t(cluster.GetNofDigis()))
189  / TMath::Sqrt(Double_t(cluster.GetNofDigis()));
190  Double_t qSum = qF + qM + qL;
191 
192  // Average charge in middle strips
193  qM /= Double_t(cluster.GetNofDigis() - 2);
194  eqMsq /= Double_t(cluster.GetNofDigis() - 2);
195 
196  // Cluster position
197  Double_t x = 0.5 * (Double_t(chanF + chanL) + (qL - qF) / qM);
198 
199  // Correct negative cluster position for clusters round the edge
200  if (x < -0.5) x += Double_t(modPar->GetNofChannels() / 2);
201 
202  // Cluster position error
203  Double_t exFsq = eqFsq / qM / qM / 4.; // error from first charge
204  Double_t exMsq = eqMsq * (qL - qF) * (qL - qF) / qM / qM / qM / qM / 4.;
205  Double_t exLsq = eqLsq / qM / qM / 4.;
206  Double_t xError = TMath::Sqrt(exFsq + exMsq + exLsq);
207 
208  // Correction for corrupt clusters
209  if (x < chanF || x > chanL) x = WeightedMean(cluster, modPar);
210 
211  assert(x >= chanF && x <= chanL);
212 
213  cluster.SetProperties(qSum, x, xError, tSum, tError);
214  cluster.SetSize(chanL - chanF + 1);
215 }
216 // --------------------------------------------------------------------------
217 
218 
219 // ----- Algorithm execution --------------------------------------------
221  const CbmStsParModule* modPar) {
222 
223  Int_t nDigis = cluster.GetNofDigis();
224  assert(nDigis >= 0);
225 
226  switch (nDigis) {
227 
228  case 0: break;
229  case 1: AnaSize1(cluster, modPar); break;
230  case 2: AnaSize2(cluster, modPar); break;
231  default: AnaSizeN(cluster, modPar); break;
232 
233  } //? number of digis
234 }
235 // --------------------------------------------------------------------------
236 
237 
238 // ----- Weighted mean calculation --------------------------------------
240  const CbmStsParModule* modPar) {
241 
242  Double_t qSum = 0.;
243  Double_t xSum = 0.;
244  for (Int_t iDigi = 0; iDigi < cluster.GetNofDigis(); iDigi++) {
245  Int_t index = cluster.GetDigi(iDigi);
246  const CbmStsDigi* digi = fDigiMan->Get<CbmStsDigi>(index);
247  assert(digi);
248  Int_t channel = digi->GetChannel();
249  auto& asic = modPar->GetParAsic(channel);
250  Double_t charge = asic.AdcToCharge(digi->GetCharge());
251  qSum += charge;
252  xSum += charge * Double_t(channel);
253  }
254 
255  return xSum / qSum;
256 }
257 // --------------------------------------------------------------------------
CbmStsParAsic::GetNoise
Double_t GetNoise() const
Electronic noise RMS.
Definition: CbmStsParAsic.h:109
CbmStsPhysics
Auxiliary class for physics processes in Silicon.
Definition: CbmStsPhysics.h:24
CbmStsAlgoAnaCluster::fDigiMan
CbmDigiManager * fDigiMan
Definition: CbmStsAlgoAnaCluster.h:93
CbmStsAlgoAnaCluster.h
CbmStsCluster
Data class for STS clusters.
Definition: CbmStsCluster.h:31
CbmStsParModule::GetNofChannels
UInt_t GetNofChannels() const
Number of channels.
Definition: CbmStsParModule.h:74
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmStsPhysics::LandauWidth
Double_t LandauWidth(Double_t mostProbableCharge)
Half width at half max of Landau distribution in ultra-relativistic case.
Definition: CbmStsPhysics.cxx:187
CbmStsCluster::SetSize
void SetSize(Int_t size)
Set size of the cluster (number of channels)
Definition: CbmStsCluster.h:138
CbmStsAlgoAnaCluster::CbmStsAlgoAnaCluster
CbmStsAlgoAnaCluster()
Constructor.
ClassImp
ClassImp(CbmStsAlgoAnaCluster) CbmStsAlgoAnaCluster
Definition: CbmStsAlgoAnaCluster.cxx:18
CbmStsCluster::SetProperties
void SetProperties(Double_t charge, Double_t position, Double_t positionError, Double_t time=0., Double_t timeError=0.)
Definition: CbmStsCluster.h:119
CbmStsAlgoAnaCluster::AnaSizeN
void AnaSizeN(CbmStsCluster &cluster, const CbmStsParModule *modPar)
Analyse cluster with more than two digis.
Definition: CbmStsAlgoAnaCluster.cxx:128
CbmStsAlgoAnaCluster::AnaSize2
void AnaSize2(CbmStsCluster &cluster, const CbmStsParModule *modPar)
Analyse two-digi cluster.
Definition: CbmStsAlgoAnaCluster.cxx:49
CbmCluster::GetNofDigis
Int_t GetNofDigis() const
Number of digis in cluster.
Definition: CbmCluster.h:69
CbmStsDigi.h
CbmStsAlgoAnaCluster::WeightedMean
Double_t WeightedMean(CbmStsCluster &cluster, const CbmStsParModule *modPar)
Weighted mean cluster position.
Definition: CbmStsAlgoAnaCluster.cxx:239
CbmStsParModule::GetParAsic
const CbmStsParAsic & GetParAsic(UInt_t channel) const
ASIC parameters for a given channel.
Definition: CbmStsParModule.cxx:29
CbmDigiManager::Get
const Digi * Get(Int_t index) const
Get a digi object.
Definition: CbmDigiManager.h:52
CbmStsParModule
Parameters for one STS module.
Definition: CbmStsParModule.h:28
CbmStsParModule.h
CbmStsDigi::GetTime
Double_t GetTime() const
Definition: CbmStsDigi.h:83
CbmStsPhysics.h
CbmStsAlgoAnaCluster::AnaSize1
void AnaSize1(CbmStsCluster &cluster, const CbmStsParModule *modPar)
Analyse single-digi cluster.
Definition: CbmStsAlgoAnaCluster.cxx:28
CbmDigiManager
CbmDigiManager.
Definition: CbmDigiManager.h:37
CbmStsDigi
Data class for a single-channel message in the STS.
Definition: CbmStsDigi.h:29
fDigiMan
CbmDigiManager * fDigiMan
Definition: CbmTofAnaTestbeam.cxx:88
CbmStsDigi::GetCharge
Double_t GetCharge() const
Definition: CbmStsDigi.h:71
CbmStsParAsic::AdcToCharge
Double_t AdcToCharge(UShort_t adc) const
Charge from ADC channel (mean)
Definition: CbmStsParAsic.h:74
CbmDigiManager.h
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmStsAlgoAnaCluster::fPhysics
CbmStsPhysics * fPhysics
Interface to digi data.
Definition: CbmStsAlgoAnaCluster.h:94
CbmStsAlgoAnaCluster::Exec
void Exec(CbmStsCluster &cluster, const CbmStsParModule *module)
Algorithm execution.
Definition: CbmStsAlgoAnaCluster.cxx:220
CbmStsAlgoAnaCluster
Determination of cluster parameters.
Definition: CbmStsAlgoAnaCluster.h:31
CbmStsCluster.h
Data class for STS clusters.
CbmCluster::GetDigi
Int_t GetDigi(Int_t index) const
Get digi at position index.
Definition: CbmCluster.h:76
CbmStsDigi::GetChannel
UShort_t GetChannel() const
Channel number in module @value Channel number.
Definition: CbmStsDigi.h:59