This method is called for each event.
All processing of the event takes place in this method.
92 {
93
94 if (not m_hits.isValid()) {
95 m_warnCount++;
96 if (m_warnCount < 10) {
97 B2WARNING("StoreArray 'CDCDedxHits' does not exist, returning. Probably running on old cdst.");
98 } else if (m_warnCount == 10) {
99 B2WARNING("StoreArray 'CDCDedxHits' does not exist, returning. ...message will be suppressed now.");
100 }
101 return;
102 }
103
104
105 m_likelihoods.clear();
106 m_dedxTracks.clear();
107
108
109 const auto& cdcgeo = CDCGeometryPar::Instance();
110 LinearGlobalADCCountTranslator adcTranslator;
111 RealisticTDCCountTranslator tdcTranslator;
112
113
114 bool isData = not Environment::Instance().isMC();
115
116
117 double runGain = isData ? m_DBRunGain->getRunGain() : 1.0;
118 double timeGain = 1;
119 double timeReso = 1;
120 if (isData and m_TTDInfo.isValid() and m_TTDInfo->hasInjection()) {
121 timeGain = m_DBInjectTime->getCorrection("mean", m_TTDInfo->isHER(), m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds());
122 timeReso = m_DBInjectTime->getCorrection("reso", m_TTDInfo->isHER(), m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds());
123 }
124 double scale = m_DBScaleFactor->getScaleFactor();
125 if (scale == 0) {
126 B2ERROR("Scale factor from DB is zero! Will be set to one");
127 scale = 1;
128 }
129
130
131 for (const auto& track : m_tracks) {
132
133 const auto* fitResult = track.getTrackFitResultWithClosestMass(Const::pion);
134 if (not fitResult) {
135 B2WARNING("No related fit for track, skip it.");
136 continue;
137 }
138
139
140 const auto hits = track.getRelationsTo<CDCDedxHit>();
141 if (hits.size() == 0) continue;
142
143
144 const auto& trackMom = fitResult->getMomentum();
145 double theta = trackMom.Theta();
146 double cosTheta = std::cos(theta);
147 double sinTheta = std::sin(theta);
148
149
150 double cosCor = isData ? m_DBCosineCor->getMean(cosTheta) : 1.0;
151 bool isEdge = std::abs(cosTheta + 0.860) < 0.010 or std::abs(cosTheta - 0.955) <= 0.005;
152 double cosEdgeCor = (isData and isEdge) ? m_DBCosEdgeCor->getMean(cosTheta) : 1.0;
153
154
155 const auto* mcParticle = isData ? nullptr : track.getRelated<MCParticle>();
156
157
158 CDCDedxTrack* dedxTrack = m_enableDebugOutput ? m_dedxTracks.appendNew() : nullptr;
159 if (dedxTrack) {
160 dedxTrack->m_track = track.getArrayIndex();
161 dedxTrack->m_charge = fitResult->getChargeSign();
162 dedxTrack->m_cosTheta = cosTheta;
163 dedxTrack->m_p = trackMom.R();
164 if (isData and m_TTDInfo.isValid() and m_TTDInfo->hasInjection()) {
165 dedxTrack->m_injring = m_TTDInfo->isHER();
166 dedxTrack->m_injtime = m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds();
167 }
168 if (mcParticle) {
169 dedxTrack->m_pdg = mcParticle->getPDG();
170 dedxTrack->m_mcmass = mcParticle->getMass();
171 const auto* mother = mcParticle->getMother();
172 dedxTrack->m_motherPDG = mother ? mother->getPDG() : 0;
173 const auto& trueMom = mcParticle->getMomentum();
174 dedxTrack->m_pTrue = trueMom.R();
175 dedxTrack->m_cosThetaTrue = std::cos(trueMom.Theta());
176 }
177 dedxTrack->m_scale = scale;
178 dedxTrack->m_cosCor = cosCor;
179 dedxTrack->m_cosEdgeCor = cosEdgeCor;
180 dedxTrack->m_runGain = runGain;
181 dedxTrack->m_timeGain = timeGain;
182 dedxTrack->m_timeReso = timeReso;
183 }
184
185
186 int lastLayer = -1;
187 double pCDC = 0;
188 std::map<int, DEDX> dedxWires;
189 for (const auto& hit : hits) {
190
191 const auto& wireID = hit.getWireID();
192 int layer = wireID.getILayer();
193 int superlayer = wireID.getISuperLayer();
194 int currentLayer = (superlayer == 0) ? layer : (8 + (superlayer - 1) * 6 + layer);
195 if (not m_useBackHalfCurlers and currentLayer < lastLayer) break;
196 lastLayer = currentLayer;
197
198
199 if (pCDC == 0) pCDC = hit.getPOCAMomentum().R();
200
201
202 double innerRadius = cdcgeo.innerRadiusWireLayer()[currentLayer];
203 double outerRadius = cdcgeo.outerRadiusWireLayer()[currentLayer];
204 const ROOT::Math::XYZVector& wirePosF = cdcgeo.wireForwardPosition(wireID, CDCGeometryPar::c_Aligned);
205 double wireRadius = wirePosF.Rho();
206 int nWires = cdcgeo.nWiresInLayer(currentLayer);
207 double topHeight = outerRadius - wireRadius;
208 double bottomHeight = wireRadius - innerRadius;
209 double topHalfWidth = M_PI * outerRadius / nWires;
210 double bottomHalfWidth = M_PI * innerRadius / nWires;
211 DedxDriftCell cell(DedxPoint(-topHalfWidth, topHeight),
212 DedxPoint(topHalfWidth, topHeight),
213 DedxPoint(bottomHalfWidth, -bottomHeight),
214 DedxPoint(-bottomHalfWidth, -bottomHeight));
215
216
217 double doca = hit.getSignedDOCAXY();
218 double entAng = hit.getEntranceAngle();
219 double celldx = cell.dx(doca, entAng) / sinTheta;
220 if (not cell.isValid()) continue;
221
222
223 int wire = wireID.getIWire();
224 int iwire = (superlayer == 0) ? 160 * layer + wire : m_nLayerWires[superlayer - 1] + (160 + 32 * (superlayer - 1)) * layer + wire;
225 double wiregain = isData ? m_DBWireGains->getWireGain(iwire) : 1.0;
226
227
228 double cellHalfWidth = M_PI * wireRadius / nWires;
229 double cellHeight = topHeight + bottomHeight;
230 double cellR = 2 * cellHalfWidth / cellHeight;
231 double tana = std::max(std::min(std::tan(entAng), 1e10), -1e10);
232 double docaRS = doca * std::sqrt((1 + cellR * cellR * tana * tana) / (1 + tana * tana));
233 double normDocaRS = docaRS / cellHalfWidth;
234 double entAngRS = std::atan(tana / cellR);
235
236
237 double onedcor = isData ? m_DB1DCell->getMean(currentLayer, entAngRS) : 1.0;
238 double twodcor = isData ? m_DB2DCell->getMean(currentLayer, normDocaRS, entAngRS) : 1.0;
239
240
241 double correction = runGain * cosCor * cosEdgeCor * timeGain * wiregain * twodcor * onedcor;
242
243
244 double adcCount = isData ? m_DBNonlADC->getCorrectedADC(hit.getADCCount(), currentLayer) : hit.getADCCount();
245 double adcCalibrated = correction != 0 ? adcCount / scale / correction : 0;
246
247
248 if (correction != 0) dedxWires[iwire].add(hit, iwire, currentLayer, celldx, adcCalibrated);
249
250
251 if (dedxTrack) {
252 dedxTrack->m_pCDC = pCDC;
253 const auto& pocaMom = hit.getPOCAMomentum();
254 double pocaPhi = pocaMom.Phi();
255 double pocaTheta = pocaMom.Theta();
256 double pocaZ = hit.getPOCAOnWire().Z();
257 double hitCharge = adcTranslator.getCharge(adcCount, wireID, false, pocaZ, pocaPhi);
258 double driftDRealistic = tdcTranslator.getDriftLength(hit.getTDCCount(), wireID, 0, true, pocaZ, pocaPhi, pocaTheta);
259 double driftDRealisticRes = tdcTranslator.getDriftLengthResolution(driftDRealistic, wireID, true, pocaZ, pocaPhi, pocaTheta);
260 double cellDedx = adcCalibrated / celldx;
261
262 dedxTrack->addHit(wire, iwire, currentLayer, doca, docaRS, entAng, entAngRS,
263 adcCount, hit.getADCCount(), hitCharge, celldx * sinTheta, cellDedx, cellHeight, cellHalfWidth,
264 hit.getTDCCount(), driftDRealistic, driftDRealisticRes, wiregain, twodcor, onedcor,
265 hit.getFoundByTrackFinder(), hit.getWeightPionHypo(), hit.getWeightKaonHypo(), hit.getWeightProtonHypo());
266 }
267
268 }
269
270
271 std::map<int, DEDX> dedxLayers;
272 for (const auto& dedxWire : dedxWires) {
273 const auto& dedx = dedxWire.second;
274 dedxLayers[dedx.cLayer].add(dedx);
275 }
276
277
278 std::vector<double> dedxValues;
279 for (const auto& dedxLayer : dedxLayers) {
280 const auto& dedx = dedxLayer.second;
281 if (dedx.dx > 0 and dedx.dE > 0) {
282 dedxValues.push_back(dedx.dE / dedx.dx);
283
284 if (dedxTrack) dedxTrack->addDedx(dedx.nhits, dedx.cWire, dedx.cLayer, dedx.dx, dedxValues.back());
285 }
286 }
287 if (dedxValues.empty()) continue;
288
289
290 std::sort(dedxValues.begin(), dedxValues.end());
291
292
293 double mean = 0;
294 for (auto x : dedxValues) mean += x;
295 mean /= dedxValues.size();
296
297
298 int lowEdgeTrunc = int(dedxValues.size() * m_removeLowest + 0.51);
299 int highEdgeTrunc = int(dedxValues.size() * (1 - m_removeHighest) + 0.51);
300 double truncatedMean = 0;
301 double sumOfSquares = 0;
302 int numValues = 0;
303 for (int i = lowEdgeTrunc; i < highEdgeTrunc; i++) {
304 double x = dedxValues[i];
305 truncatedMean += x;
306 sumOfSquares += x * x;
307 numValues++;
308 }
309 if (numValues > 0) {
310 truncatedMean /= numValues;
311 } else {
312 truncatedMean = mean;
313 numValues = dedxValues.size();
314 }
315 double truncatedError = numValues > 1 ? std::sqrt(sumOfSquares / numValues - truncatedMean * truncatedMean) / (numValues - 1) : 0;
316
317
318 double correctedMean = isData ? m_DBHadronCor->getCorrectedMean(truncatedMean, cosTheta) : truncatedMean;
319
320
321 if (m_trackLevel and mcParticle) {
322 double mass = mcParticle->getMass();
323 if (mass > 0) {
324 double mcMean = m_DBMeanPars->getMean(pCDC / mass);
325 double mcSigma = m_DBSigmaPars->getSigma(mcMean, numValues, cosTheta, timeReso);
326 correctedMean = gRandom->Gaus(mcMean, mcSigma);
327 while (correctedMean < 0) correctedMean = gRandom->Gaus(mcMean, mcSigma);
328
329 if (dedxTrack) dedxTrack->m_simDedx = correctedMean;
330 }
331 }
332
333
334 double cdcLogL[Const::ChargedStable::c_SetSize] = {0};
335 for (const auto& chargedStable : Const::chargedStableSet) {
336 double betagamma = pCDC / chargedStable.getMass();
337 double predictedMean = m_DBMeanPars->getMean(betagamma);
338 double predictedSigma = m_DBSigmaPars->getSigma(predictedMean, numValues, cosTheta, timeReso);
339 if (predictedSigma <= 0) B2ERROR("Predicted sigma is not positive for PDG = " << chargedStable.getPDGCode());
340 double chi = (correctedMean - predictedMean) / predictedSigma;
341 int index = chargedStable.getIndex();
342 cdcLogL[index] = -0.5 * chi * chi;
343
344 if (dedxTrack) {
345 dedxTrack->m_predmean[index] = predictedMean;
346 dedxTrack->m_predres[index] = predictedSigma;
347 dedxTrack->m_cdcChi[index] = chi;
348 dedxTrack->m_cdcLogl[index] = cdcLogL[index];
349 }
350 }
351
352
353 auto* likelihoods = m_likelihoods.appendNew(cdcLogL);
354 track.addRelationTo(likelihoods);
355
356
357 if (dedxTrack) {
358 double fullLength = 0;
359 for (const auto& dedxLayer : dedxLayers) fullLength += dedxLayer.second.dx;
360 dedxTrack->m_length = fullLength;
361 dedxTrack->m_dedxAvg = mean;
362 dedxTrack->m_dedxAvgTruncatedNoSat = truncatedMean;
363 dedxTrack->m_dedxAvgTruncatedErr = truncatedError;
364 dedxTrack->m_dedxAvgTruncated = correctedMean;
365 dedxTrack->m_lNHitsUsed = numValues;
366 track.addRelationTo(dedxTrack);
367 }
368
369 }
370
371 }