Belle II Software development
HitLevelInfoWriter.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <cdc/modules/HitLevelInfoWriter/HitLevelInfoWriter.h>
10
11#include <analysis/dataobjects/ParticleList.h>
12#include <analysis/dataobjects/Particle.h>
13#include <analysis/utility/ReferenceFrame.h>
14#include <mdst/dataobjects/HitPatternCDC.h>
15#include <mdst/dataobjects/PIDLikelihood.h>
16#include <mdst/dbobjects/BeamSpot.h>
17
18#include <TMath.h>
19
20using namespace Belle2;
21using namespace std;
22
23REG_MODULE(HitLevelInfoWriter);
24
26{
27
28 setDescription("Extract dE/dx information for calibration development.");
29
30 addParam("outputBaseName", m_strOutputBaseName, "Suffix for output file name", string("HLInfo.root"));
31 addParam("particleLists", m_strParticleList, "Vector of ParticleLists to save", vector<string>());
32 addParam("enableHitLevel", m_isHitLevel, "True or False for Hit level variables", false);
33 addParam("m_isExtraVar", m_isExtraVar, "True or False for extra track/hit level variables", false);
34 addParam("nodeadwire", m_isDeadwire, "True or False for deadwire hit variables", false);
35 addParam("relativeCorrections", m_isRelative, "If true, apply corrections relative to those used in production", false);
36 addParam("corrections", m_isCorrection, "If true, apply corrections", true);
37
38}
39
41
43{
44
45 B2INFO("Creating a ROOT file for the hit level information...");
46
47 // required inputs
48 m_dedxTracks.isRequired();
49 m_tracks.isRequired();
50 m_trackFitResults.isRequired();
51 m_eclClusters.isOptional();
52 m_klmClusters.isOptional();
53
54 // build a map to relate input strings to the right particle type
55 map<string, string> pdgMap = {{"pi+", "Const::pion.getPDGCode()"}, {"K+", "Const::kaon.getPDGCode()"}, {"mu+", "Const::muon.getPDGCode()"}, {"e+", "Const::electron.getPDGCode()"}, {"p+", "Const::proton.getPDGCode()"}, {"deuteron", "Const::deuteron.getPDGCode()"}};
56
57 // if no particle lists are given, write out all tracks
59
60 // create a new output file for each particle list specified
61 for (unsigned int i = 0; i < m_strParticleList.size(); i++) {
62 // strip the name of the particle lists to make this work
63 string pdg = pdgMap[m_strParticleList[i].substr(0, m_strParticleList[i].find(":"))];
64 string filename = string(m_strOutputBaseName + "_PID" + pdg + ".root");
65 bookOutput(filename);
66 }
67}
68
70{
71
72 int nParticleList = m_strParticleList.size();
73
74 // **************************************************
75 // LOOP OVER dE/dx measurements for all tracks if
76 // no particle list is specified.
77 // **************************************************
78 if (nParticleList == 0) {
79 for (int idedx = 0; idedx < m_dedxTracks.getEntries(); idedx++) {
80
81 CDCDedxTrack* dedxTrack = m_dedxTracks[idedx];
82 if (!dedxTrack) {
83 B2WARNING("No dedx related track...");
84 continue;
85 }
86
87 const Track* track = dedxTrack->getRelatedFrom<Track>();
88 if (!track) {
89 B2WARNING("No related track...");
90 continue;
91 }
92
93 const TrackFitResult* fitResult = track->getTrackFitResultWithClosestMass(Const::pion);
94 if (!fitResult) {
95 B2WARNING("No related fit for this track...");
96 continue;
97 }
98
99 if (dedxTrack->size() == 0 || dedxTrack->size() > 200) continue;
100 if (dedxTrack->getCosTheta() < TMath::Cos(150.0 * TMath::DegToRad()))continue; //-0.866
101 if (dedxTrack->getCosTheta() > TMath::Cos(17.0 * TMath::DegToRad())) continue; //0.95
102
103 // fill the event meta data
104 StoreObjPtr<EventMetaData> evtMetaData;
105 m_expID = evtMetaData->getExperiment();
106 m_runID = evtMetaData->getRun();
107 m_eventID = evtMetaData->getEvent();
108
109 m_injring = dedxTrack->getInjectionRing();
110 m_injtime = dedxTrack->getInjectionTime();
111
112 //--------REMOVEBAL--------
113 //when CDST are reproduced with injection time
114 if (m_injtime == -1 || m_injring == -1) {
115 if (m_TTDInfo.isValid() && m_TTDInfo->hasInjection()) {
116 m_injring = m_TTDInfo->isHER();
117 m_injtime = m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds();
118 }
119 }
120 //--------REMOVEBAL--------
121
122 // fill the E/P
123 const ECLCluster* eclCluster = track->getRelated<ECLCluster>();
124 if (eclCluster and eclCluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
125 m_eop = (eclCluster->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons)) / (fitResult->getMomentum().R());
127 if (m_isExtraVar) {
128 m_e1_9 = eclCluster->getE1oE9();
129 m_e9_21 = eclCluster->getE9oE21();
130 m_eclsnHits = eclCluster->getNumberOfCrystals();
131 }
132 // fill the muon depth
133 const KLMCluster* klmCluster = eclCluster->getRelated<KLMCluster>();
134 if (klmCluster) m_klmLayers = klmCluster->getLayers();
135 }
136
137 // fill the TTree with the Track information
138 fillTrack(fitResult);
139
140 // fill the TTree with the CDCDedxTrack information
141 fillDedx(dedxTrack);
142
143 // fill the TTree
144 m_tree[0]->Fill();
145 }
146 }
147
148 // **************************************************
149 //
150 // LOOP OVER particles in the given particle lists
151 //
152 // **************************************************
153 for (int iList = 0; iList < nParticleList; iList++) {
154 // make sure the list exists and is not empty
155 StoreObjPtr<ParticleList> particlelist(m_strParticleList[iList]);
156 if (!particlelist or particlelist->getListSize(true) == 0) {
157 //B2WARNING("ParticleList " << m_strParticleList[iList] << " not found or empty, skipping");
158 continue;
159 }
160
161 // loop over the particles in the list and follow the links to the
162 // dE/dx information (Particle -> PIDLikelihood -> Track -> CDCDedxTrack)
163 for (unsigned int iPart = 0; iPart < particlelist->getListSize(true); iPart++) {
164 Particle* part = particlelist->getParticle(iPart, true);
165 if (!part) {
166 B2WARNING("No particles...");
167 continue;
168 }
170 if (!pid) {
171 B2WARNING("No related PID likelihood...");
172 continue;
173 }
174 Track* track = pid->getRelatedFrom<Track>();
175 if (!track) {
176 B2WARNING("No related track...");
177 continue;
178 }
179 CDCDedxTrack* dedxTrack = track->getRelatedTo<CDCDedxTrack>();
180 if (!dedxTrack) {
181 B2WARNING("No related CDCDedxTrack...");
182 continue;
183 }
184 string ptype = m_strParticleList[iList].substr(0, m_strParticleList[iList].find(":"));
185 const TrackFitResult* fitResult = track->getTrackFitResult(Const::pion);
186 if (ptype != "pi+") {
187 if (ptype == "K+") fitResult = track->getTrackFitResultWithClosestMass(Const::kaon);
188 else if (ptype == "p+") fitResult = track->getTrackFitResultWithClosestMass(Const::proton);
189 else if (ptype == "deuteron") fitResult = track->getTrackFitResultWithClosestMass(Const::deuteron);
190 else if (ptype == "mu+") fitResult = track->getTrackFitResultWithClosestMass(Const::muon);
191 else if (ptype == "e+") fitResult = track->getTrackFitResultWithClosestMass(Const::electron);
192 }
193 if (!fitResult) {
194 B2WARNING("No related fit for this track...");
195 continue;
196 }
197
198 if (dedxTrack->size() == 0 || dedxTrack->size() > 200) continue;
199 //if out CDC (dont add as we dont correct via correctionmodules)
200 if (dedxTrack->getCosTheta() < TMath::Cos(150.0 * TMath::DegToRad()))continue; //-0.866
201 if (dedxTrack->getCosTheta() > TMath::Cos(17.0 * TMath::DegToRad())) continue; //0.95
202
203 // fill the event meta data
204 StoreObjPtr<EventMetaData> evtMetaData;
205 m_expID = evtMetaData->getExperiment();
206 m_runID = evtMetaData->getRun();
207 m_eventID = evtMetaData->getEvent();
208
209 m_injring = dedxTrack->getInjectionRing();
210 m_injtime = dedxTrack->getInjectionTime();
211
212 //--------REMOVEBAL--------
213 //when CDST are reproduced with injection time
214 if (m_injtime == -1 || m_injring == -1) {
215 if (m_TTDInfo.isValid() && m_TTDInfo->hasInjection()) {
216 m_injring = m_TTDInfo->isHER();
217 m_injtime = m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds();
218 }
219 }
220 //--------REMOVEBAL--------
221
222 // fill the E/P
223 const ECLCluster* eclCluster = track->getRelated<ECLCluster>();
224 if (eclCluster and eclCluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
225 m_eop = (eclCluster->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons)) / (fitResult->getMomentum().R());
227 if (m_isExtraVar) {
228 m_e1_9 = eclCluster->getE1oE9();
229 m_e9_21 = eclCluster->getE9oE21();
230 m_eclsnHits = eclCluster->getNumberOfCrystals();
231 }
232 // fill the muon depth
233 const KLMCluster* klmCluster = eclCluster->getRelated<KLMCluster>();
234 if (klmCluster) m_klmLayers = klmCluster->getLayers();
235 }
236
237 // fill the TTree with the Track information
238 fillTrack(fitResult);
239
240 // fill the TTree with the CDCDedxTrack information
241 fillDedx(dedxTrack);
242
243 // fill the TTree
244 m_tree[iList]->Fill();
245 }
246 }
247}
248
249//---------------------------------------------------------------------------------------
251{
252
253 for (unsigned int i = 0; i < m_file.size(); i++) {
254 B2INFO("Done writing out the hit level information...\t" << m_tree[i]->GetEntries() << " tracks");
255
256 // write the ttree to a root file
257 m_file[i]->cd();
258 m_tree[i]->Write();
259 m_file[i]->Close();
260 }
261}
262
263//---------------------------------------------------------------------------------------
265{
266 ROOT::Math::XYZVector trackMom = fitResult->getMomentum();
267 m_p = trackMom.R();
268 m_pt = trackMom.Rho();
269 m_phi = trackMom.Phi();
270
271 m_theta = trackMom.Theta() * 180. / TMath::Pi(); //in degree
272 if (m_theta > 17. && m_theta < 150.)m_inCDC = 1;
273 else m_inCDC = 0;
274
275 if (fitResult->getChargeSign() < 0) {
276 m_p *= -1;
277 m_pt *= -1;
278 }
279
280 ROOT::Math::XYZVector trackPos = fitResult->getPosition();
281 m_vx0 = trackPos.X();
282 m_vy0 = trackPos.Y();
283 m_vz0 = trackPos.Z();
284
285 m_d0 = fitResult->getD0();
286 m_z0 = fitResult->getZ0();
287 m_chi2 = fitResult->getPValue();
288 m_tanlambda = fitResult->getTanLambda();
289 m_phi0 = fitResult->getPhi0();
290 m_nCDChits = fitResult->getHitPatternCDC().getNHits();
291
292 static DBObjPtr<BeamSpot> beamSpotDB;
293 const auto& frame = ReferenceFrame::GetCurrent();
294 UncertainHelix helix = fitResult->getUncertainHelix();
295 helix.passiveMoveBy(ROOT::Math::XYZVector(beamSpotDB->getIPPosition()));
296 m_dr = frame.getVertex(ROOT::Math::XYZVector(helix.getPerigee())).Rho();
297 m_dphi = frame.getVertex(ROOT::Math::XYZVector(helix.getPerigee())).Phi();
298 m_dz = frame.getVertex(ROOT::Math::XYZVector(helix.getPerigee())).Z();
299
300}
301
302//---------------------------------------------------------------------------------------
304{
305 // clear the containers first
306 clearEntries();
307
308 m_trackID = dedxTrack->trackID();
309 m_length = dedxTrack->getLength();
310 m_charge = dedxTrack->getCharge();
311 m_cosTheta = dedxTrack->getCosTheta();
312 m_PDG = dedxTrack->getPDG();
313
314 m_pCDC = dedxTrack->getMomentum();
315 if (m_charge < 0) {
316 m_pCDC *= -1;
317 }
318
319 h_nhits = dedxTrack->size();
320
321 // Get the calibration constants
322 m_scale = m_DBScaleFactor->getScaleFactor();
323 m_runGain = m_DBRunGain->getRunGain();
325 m_timeGain = m_DBInjectTime->getCorrection("mean", m_injring, m_injtime);
326 m_timeReso = m_DBInjectTime->getCorrection("reso", m_injring, m_injtime);
327
330 } else {
331 m_cosEdgeCor = 1.0;
332 }
333 m_hadronpars = m_DBHadronCor->getHadronPars();
334
335 //variable to save layer variables
336 map<int, vector<double>>l_var;
337 double cdcChi[Const::ChargedStable::c_SetSize];
338
339 //Modify the hit level dedx
340 if (m_isCorrection) {
341 recalculateDedx(dedxTrack, l_var, cdcChi);
342 m_chie = cdcChi[0];
343 m_chimu = cdcChi[1];
344 m_chipi = cdcChi[2];
345 m_chik = cdcChi[3];
346 m_chip = cdcChi[4];
347 m_chid = cdcChi[5];
348 } else {
349 m_mean = dedxTrack->getDedxMean();
350 m_trunc = dedxTrack->getDedx();
351 m_truncNoSat = dedxTrack->getDedxNoSat();
352 m_error = dedxTrack->getDedxError();
353 m_chie = dedxTrack->getChi(0);
354 m_chimu = dedxTrack->getChi(1);
355 m_chipi = dedxTrack->getChi(2);
356 m_chik = dedxTrack->getChi(3);
357 m_chip = dedxTrack->getChi(4);
358 m_chid = dedxTrack->getChi(5);
359 }
360
361 m_pmeane = dedxTrack->getPmean(0);
362 m_pmeanmu = dedxTrack->getPmean(1);
363 m_pmeanpi = dedxTrack->getPmean(2);
364 m_pmeank = dedxTrack->getPmean(3);
365 m_pmeanp = dedxTrack->getPmean(4);
366 m_pmeand = dedxTrack->getPmean(5);
367
368 m_prese = dedxTrack->getPreso(0);
369 m_presmu = dedxTrack->getPreso(1);
370 m_prespi = dedxTrack->getPreso(2);
371 m_presk = dedxTrack->getPreso(3);
372 m_presp = dedxTrack->getPreso(4);
373 m_presd = dedxTrack->getPreso(5);
374
375 // Get the vector of dE/dx values for all layers
376 double lout = 0, lin = 0, increment = 0;
377 int lastlayer = 0;
378
379 if (m_isCorrection) {
380 l_nhits = l_var[0].size();
381 const int lEdgeTrunc = int(l_nhits * 0.05 + 0.51);
382 const int hEdgeTrunc = int(l_nhits * (1 - 0.25) + 0.51);
383 l_nhitsused = hEdgeTrunc - lEdgeTrunc;
384 } else {
385 l_nhits = dedxTrack->getNLayerHits();
386 l_nhitsused = dedxTrack->getNLayerHitsUsed();
387 }
388
389 for (int il = 0; il < l_nhits; ++il) {
390 if (m_isCorrection) {
391 l_nhitscombined[il] = l_var[0][il];
392 l_wirelongesthit[il] = l_var[1][il];
393 l_layer[il] = l_var[2][il];
394 l_path[il] = l_var[3][il];
395 l_dedx[il] = l_var[4][il];
396 } else {
397 l_nhitscombined[il] = dedxTrack->getNHitsCombined(il);
398 l_wirelongesthit[il] = dedxTrack->getWireLongestHit(il);
399 l_layer[il] = dedxTrack->getLayer(il);
400 l_path[il] = dedxTrack->getLayerPath(il);
401 l_dedx[il] = dedxTrack->getLayerDedx(il);
402 }
403 if (l_layer[il] > lastlayer) lout++;
404 else if (l_layer[il] < lastlayer) lin++;
405 else continue;
406
407 lastlayer = l_layer[il];
408 increment++;
409 }
410 m_ioasym = (lout - lin) / increment;
411
412 // Get the vector of dE/dx values for all hits
413 if (m_isHitLevel) {
414 for (int ihit = 0; ihit < h_nhits; ++ihit) {
415
416 if (m_isDeadwire) continue;
417
418 h_lwire[ihit] = dedxTrack->getWireInLayer(ihit);
419 h_wire[ihit] = dedxTrack->getWire(ihit);
420 h_layer[ihit] = dedxTrack->getHitLayer(ihit);
421 h_path[ihit] = dedxTrack->getPath(ihit);
422 h_dedx[ihit] = dedxTrack->getDedx(ihit);
423 h_adcraw[ihit] = dedxTrack->getADCBaseCount(ihit);
424 h_adccorr[ihit] = dedxTrack->getADCCount(ihit);
425 h_doca[ihit] = dedxTrack->getDoca(ihit);
426 h_ndoca[ihit] = h_doca[ihit] / dedxTrack->getCellHalfWidth(ihit);
427 h_ndocaRS[ihit] = dedxTrack->getDocaRS(ihit) / dedxTrack->getCellHalfWidth(ihit);
428 h_enta[ihit] = dedxTrack->getEnta(ihit);
429 h_entaRS[ihit] = dedxTrack->getEntaRS(ihit);
430 h_driftT[ihit] = dedxTrack->getDriftT(ihit);
431 h_driftD[ihit] = dedxTrack->getDriftD(ihit);
432
433 // Get extra variable from tracking
434 if (m_isExtraVar) {
435 h_WeightPionHypo[ihit] = dedxTrack->getWeightPionHypo(ihit);
436 h_WeightKaonHypo[ihit] = dedxTrack->getWeightKaonHypo(ihit);
437 h_WeightProtonHypo[ihit] = dedxTrack->getWeightProtonHypo(ihit);
438 h_foundByTrackFinder[ihit] = dedxTrack->getFoundByTrackFinder(ihit);
439 }
440 // Get the calibration constants
441 h_facnladc[ihit] = dedxTrack->getNonLADCCorrection(ihit);
442 h_wireGain[ihit] = m_DBWireGains->getWireGain(h_wire[ihit]);
443 h_twodCor[ihit] = m_DB2DCell->getMean(h_layer[ihit], h_ndocaRS[ihit], h_entaRS[ihit]);
444 h_onedCor[ihit] = m_DB1DCell->getMean(h_layer[ihit], h_entaRS[ihit]);
445 }
446 }
447
448}
449
450//---------------------------------------------------------------------------------------
451void HitLevelInfoWriterModule::recalculateDedx(CDCDedxTrack* dedxTrack, map<int, vector<double>>& l_var,
452 double (&cdcChi)[Const::ChargedStable::c_SetSize])
453{
454 vector<double> newLayerHits;
455 double newLayerDe = 0, newLayerDx = 0;
456 int nhitscombined = 0; // number of hits combined per layer
457 int wirelongesthit = 0; // wire number of longest hit
458 double longesthit = 0; // path length of longest hit
459
460 for (int ihit = 0; ihit < h_nhits; ++ihit) {
461 int jadcbase = dedxTrack->getADCBaseCount(ihit);
462 int jLayer = dedxTrack->getHitLayer(ihit);
463 double jWire = dedxTrack->getWire(ihit);
464 double jNDocaRS = dedxTrack->getDocaRS(ihit) / dedxTrack->getCellHalfWidth(ihit);
465 double jEntaRS = dedxTrack->getEntaRS(ihit);
466 double jPath = dedxTrack->getPath(ihit);
467
468 double correction = dedxTrack->getScaleFactor() * dedxTrack->getRunGain() * dedxTrack->getTimeMean() *
469 dedxTrack->getCosineCorrection() * dedxTrack->getCosEdgeCorrection() *
470 dedxTrack->getTwoDCorrection(ihit) * dedxTrack->getOneDCorrection(ihit) *
471 dedxTrack->getNonLADCCorrection(ihit);
472 if (dedxTrack->getWireGain(ihit) > 0) correction *= dedxTrack->getWireGain(ihit); //also keep dead wire
473
474 if (m_isRelative) {
475 //get same base adc + rel correction factor
476 correction *= GetCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, m_cosTheta, m_injring, m_injtime);
477 if (!m_DBWireGains && dedxTrack->getWireGain(ihit) == 0) correction = 0;
478 } else {
479 //get modified adc + abs correction factor
480 correction = GetCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, m_cosTheta, m_injring, m_injtime);
481 }
482
483 double newhitdedx = 1.0;
484
485 newhitdedx *= jadcbase * sqrt(1 - m_cosTheta * m_cosTheta) / jPath;
486
487 if (correction != 0) {
488 newhitdedx /= correction;
489 if (m_DBWireGains->getWireGain(jWire) != 0) {
490 newLayerDe += jadcbase / correction;
491 newLayerDx += jPath;
492 if (jPath > longesthit) {
493 longesthit = jPath;
494 wirelongesthit = jWire;
495 }
496 nhitscombined++;
497 }
498 } else newhitdedx = 0;
499
500 dedxTrack->setDedx(ihit, newhitdedx);
501
502 if (ihit + 1 < h_nhits && dedxTrack->getHitLayer(ihit + 1) == jLayer) {
503 continue;
504 } else {
505 if (newLayerDx != 0) {
506 double totalDistance = newLayerDx / sqrt(1 - m_cosTheta * m_cosTheta);
507 double newLayerDedx = newLayerDe / totalDistance ;
508 newLayerHits.push_back(newLayerDedx);
509 l_var[0].push_back(nhitscombined);
510 l_var[1].push_back(wirelongesthit);
511 l_var[2].push_back(jLayer);
512 l_var[3].push_back(totalDistance);
513 l_var[4].push_back(newLayerDedx);
514 }
515
516 newLayerDe = 0;
517 newLayerDx = 0;
518 nhitscombined = 0;
519 wirelongesthit = 0;
520 longesthit = 0;
521 }
522 }
523
524 // recalculate the truncated means
525 double dedxmean, dedxtrunc, dedxtruncNoSat, dedxerror;
526
527 calculateMeans(&dedxmean, &dedxtruncNoSat, &dedxerror, newLayerHits);
528
529 dedxtrunc = dedxtruncNoSat;
530
531 HadronCorrection(m_cosTheta, dedxtrunc);
532 m_mean = dedxmean;
533 m_trunc = dedxtrunc;
534 m_truncNoSat = dedxtruncNoSat;
535 m_error = dedxerror;
536
537 // save the PID information
538 saveChiValue(cdcChi, dedxTrack, m_trunc);
539}
540
541//---------------------------------------------------------------------------------------
542double HitLevelInfoWriterModule::GetCorrection(int& adc, int layer, int wireID, double doca, double enta,
543 double costheta, double ring, double time) const
544{
545 double correction = 1.0;
546 correction *= m_DBScaleFactor->getScaleFactor();
547 correction *= m_DBRunGain->getRunGain();
548 correction *= m_DB2DCell->getMean(layer, doca, enta);
549 correction *= m_DB1DCell->getMean(layer, enta);
550 correction *= m_DBCosineCor->getMean(costheta);
551 correction *= m_DBInjectTime->getCorrection("mean", ring, time);
552 if (costheta <= -0.850 || costheta >= 0.950) correction *= m_DBCosEdgeCor->getMean(costheta);
553 if (m_DBWireGains->getWireGain(wireID) > 0) correction *= m_DBWireGains->getWireGain(wireID);
554
555 //last is only for abs constant
556 if (!m_isRelative) adc = m_DBNonlADC->getCorrectedADC(adc, layer);
557
558 return correction;
559}
560
561//---------------------------------------------------------------------------------------
562void HitLevelInfoWriterModule::HadronCorrection(double costheta, double& dedx) const
563{
564 dedx = D2I(costheta, I2D(costheta, 1.00) / 1.00 * dedx);
565}
566
567//---------------------------------------------------------------------------------------
568double HitLevelInfoWriterModule::D2I(const double cosTheta, const double D) const
569{
570 double absCosTheta = fabs(cosTheta);
571 double projection = pow(absCosTheta, m_hadronpars[3]) + m_hadronpars[2];
572 if (projection == 0) {
573 B2WARNING("Something wrong with dE/dx hadron constants!");
574 return D;
575 }
576
577 double chargeDensity = D / projection;
578 double numerator = 1 + m_hadronpars[0] * chargeDensity;
579 double denominator = 1 + m_hadronpars[1] * chargeDensity;
580
581 if (denominator == 0) {
582 B2WARNING("Something wrong with dE/dx hadron constants!");
583 return D;
584 }
585
586 double I = D * m_hadronpars[4] * numerator / denominator;
587 return I;
588}
589
590//---------------------------------------------------------------------------------------
591double HitLevelInfoWriterModule::I2D(const double cosTheta, const double I) const
592{
593 double absCosTheta = fabs(cosTheta);
594 double projection = pow(absCosTheta, m_hadronpars[3]) + m_hadronpars[2];
595
596 if (projection == 0 || m_hadronpars[4] == 0) {
597 B2WARNING("Something wrong with dE/dx hadron constants!");
598 return I;
599 }
600
601 double a = m_hadronpars[0] / projection;
602 double b = 1 - m_hadronpars[1] / projection * (I / m_hadronpars[4]);
603 double c = -1.0 * I / m_hadronpars[4];
604
605 if (b == 0 && a == 0) {
606 B2WARNING("both a and b coefficiants for hadron correction are 0");
607 return I;
608 }
609
610 double discr = b * b - 4.0 * a * c;
611 if (discr < 0) {
612 B2WARNING("negative discriminant; return uncorrectecd value");
613 return I;
614 }
615
616 double D = (a != 0) ? (-b + sqrt(discr)) / (2.0 * a) : -c / b;
617 if (D < 0) {
618 B2WARNING("D is less 0! will try another solution");
619 D = (a != 0) ? (-b - sqrt(discr)) / (2.0 * a) : -c / b;
620 if (D < 0) {
621 B2WARNING("D is still less 0! just return uncorrectecd value");
622 return I;
623 }
624 }
625
626 return D;
627}
628
629//---------------------------------------------------------------------------------------
630void HitLevelInfoWriterModule::calculateMeans(double* mean, double* truncMean, double* truncMeanErr,
631 const vector<double>& dedx) const
632{
633 // Calculate the truncated average by skipping the lowest & highest
634 // events in the array of dE/dx values
635 vector<double> sortedDedx = dedx;
636 sort(sortedDedx.begin(), sortedDedx.end());
637 sortedDedx.erase(remove(sortedDedx.begin(), sortedDedx.end(), 0), sortedDedx.end());
638 sortedDedx.shrink_to_fit();
639
640 double truncMeanTmp = 0.0, meanTmp = 0.0, sumOfSqs = 0.0;
641 int nValTrunc = 0;
642 const int numDedx = sortedDedx.size();
643
644 // add a factor of 0.51 here to make sure we are rounding appropriately...
645 const int lEdgeTrunc = int(numDedx * 0.05 + 0.51);
646 const int hEdgeTrunc = int(numDedx * (1 - 0.25) + 0.51);
647 for (int i = 0; i < numDedx; i++) {
648 meanTmp += sortedDedx[i];
649 if (i >= lEdgeTrunc and i < hEdgeTrunc) {
650 truncMeanTmp += sortedDedx[i];
651 sumOfSqs += sortedDedx[i] * sortedDedx[i];
652 nValTrunc++;
653 }
654 }
655
656 if (numDedx != 0) meanTmp /= numDedx;
657
658 if (nValTrunc != 0) truncMeanTmp /= nValTrunc;
659 else truncMeanTmp = meanTmp;
660
661 *mean = meanTmp;
662 *truncMean = truncMeanTmp;
663
664 if (nValTrunc > 1)
665 *truncMeanErr = sqrt(sumOfSqs / double(nValTrunc) - truncMeanTmp * truncMeanTmp) / double(nValTrunc - 1);
666 else *truncMeanErr = 0;
667
668}
669
670//---------------------------------------------------------------------------------------
672 double dedx) const
673{
674 // determine a chi value for each particle type
676 for (const Const::ChargedStable pdgIter : set) {
677
678 // determine the predicted mean and resolution
679 double mean = dedxTrack->getPmean(pdgIter.getIndex());
680 double sigma = dedxTrack->getPreso(pdgIter.getIndex());
681
682 // fill the chi value for this particle type
683 if (sigma != 0) chi[pdgIter.getIndex()] = ((dedx - mean) / (sigma));
684 }
685}
686
687//---------------------------------------------------------------------------------------
689{
690 for (int il = 0; il < 200; ++il) {
691 l_nhitscombined[il] = 0;
692 l_wirelongesthit[il] = 0;
693 l_layer[il] = 0;
694 l_path[il] = 0;
695 l_dedx[il] = 0;
696 }
697
698 if (m_isHitLevel) {
699 for (int ihit = 0; ihit < 200; ++ihit) {
700 h_lwire[ihit] = 0;
701 h_wire[ihit] = 0;
702 h_layer[ihit] = 0;
703 h_path[ihit] = 0;
704 h_dedx[ihit] = 0;
705 h_adcraw[ihit] = 0;
706 h_adccorr[ihit] = 0;
707 h_doca[ihit] = 0;
708 h_ndoca[ihit] = 0;
709 h_ndocaRS[ihit] = 0;
710 h_enta[ihit] = 0;
711 h_entaRS[ihit] = 0;
712 h_driftT[ihit] = 0;
713 // h_driftD[ihit] = 0;
714 h_facnladc[ihit] = 0;
715 h_wireGain[ihit] = 0;
716 h_twodCor[ihit] = 0;
717 h_onedCor[ihit] = 0;
718 h_WeightPionHypo[ihit] = 0;
719 h_WeightKaonHypo[ihit] = 0;
720 h_WeightProtonHypo[ihit] = 0;
721 h_foundByTrackFinder[ihit] = 0;
722 }
723 }
724}
725
726//---------------------------------------------------------------------------------------
728{
729 // register output root file
730 m_file.push_back(new TFile(filename.c_str(), "RECREATE"));
731 m_tree.push_back(new TTree("track", "dE/dx information"));
732
733 int i = m_tree.size() - 1;
734 m_tree[i]->SetDirectory(0);
735
736 // event level information (from emd)
737 m_tree[i]->Branch("exp", &m_expID, "exp/I");
738 m_tree[i]->Branch("run", &m_runID, "run/I");
739 m_tree[i]->Branch("event", &m_eventID, "event/I");
740
741 // track level information (from tfr)
742 m_tree[i]->Branch("dz", &m_dz, "dz/D");
743 m_tree[i]->Branch("dr", &m_dr, "dr/D");
744 m_tree[i]->Branch("chi2", &m_chi2, "chi2/D");
745 m_tree[i]->Branch("injtime", &m_injtime, "injtime/D");
746 m_tree[i]->Branch("isher", &m_injring, "isher/D");
747
748 if (m_isExtraVar) {
749 m_tree[i]->Branch("d0", &m_d0, "d0/D");
750 m_tree[i]->Branch("z0", &m_z0, "z0/D");
751 m_tree[i]->Branch("dphi", &m_dphi, "dphi/D");
752 m_tree[i]->Branch("vx0", &m_vx0, "vx0/D");
753 m_tree[i]->Branch("vy0", &m_vy0, "vy0/D");
754 m_tree[i]->Branch("vz0", &m_vz0, "vz0/D");
755 m_tree[i]->Branch("tanlambda", &m_tanlambda, "tanlambda/D");
756 m_tree[i]->Branch("phi0", &m_phi0, "phi0/D");
757 m_tree[i]->Branch("e1_9", &m_e1_9, "e1_9/D");
758 m_tree[i]->Branch("e9_21", &m_e9_21, "e9_21/D");
759 m_tree[i]->Branch("eclsnhits", &m_eclsnHits, "eclsnhits/D");
760 m_tree[i]->Branch("klmLayers", &m_klmLayers, "klmLayers/I");
761 }
762
763 // track level information (from cdt)
764 m_tree[i]->Branch("nCDChits", &m_nCDChits, "nCDChits/D");
765 m_tree[i]->Branch("inCDC", &m_inCDC, "inCDC/I");
766 m_tree[i]->Branch("track", &m_trackID, "track/I");
767 m_tree[i]->Branch("length", &m_length, "length/D");
768 m_tree[i]->Branch("charge", &m_charge, "charge/I");
769 m_tree[i]->Branch("costh", &m_cosTheta, "costh/D");
770 m_tree[i]->Branch("pF", &m_pCDC, "pF/D");
771 m_tree[i]->Branch("p", &m_p, "p/D");
772 m_tree[i]->Branch("pt", &m_pt, "pt/D");
773 m_tree[i]->Branch("ioasym", &m_ioasym, "ioasym/D");
774 m_tree[i]->Branch("phi", &m_phi, "phi/D");
775 m_tree[i]->Branch("pdg", &m_PDG, "pdg/D");
776 // track level dE/dx measurements
777 m_tree[i]->Branch("mean", &m_mean, "mean/D");
778 m_tree[i]->Branch("dedx", &m_trunc, "dedx/D");
779 m_tree[i]->Branch("dedxnosat", &m_truncNoSat, "dedxnosat/D");
780 m_tree[i]->Branch("dedxerr", &m_error, "dedxerr/D");
781
782 // other track level information
783 m_tree[i]->Branch("eop", &m_eop, "eop/D");
784 m_tree[i]->Branch("e", &m_e, "e/D");
785
786 // calibration constants
787 m_tree[i]->Branch("scale", &m_scale, "scale/D");
788 m_tree[i]->Branch("coscor", &m_cosCor, "coscor/D");
789 m_tree[i]->Branch("cosedgecor", &m_cosEdgeCor, "cosedgecor/D");
790 m_tree[i]->Branch("rungain", &m_runGain, "rungain/D");
791 m_tree[i]->Branch("timegain", &m_timeGain, "timegain/D");
792 m_tree[i]->Branch("timereso", &m_timeReso, "timereso/D");
793
794 // PID values
795 m_tree[i]->Branch("chiE", &m_chie, "chiE/D");
796 m_tree[i]->Branch("chiMu", &m_chimu, "chiMu/D");
797 m_tree[i]->Branch("chiPi", &m_chipi, "chiPi/D");
798 m_tree[i]->Branch("chiK", &m_chik, "chiK/D");
799 m_tree[i]->Branch("chiP", &m_chip, "chiP/D");
800 m_tree[i]->Branch("chiD", &m_chid, "chiD/D");
801
802 m_tree[i]->Branch("chiEOld", &m_chieOld, "chiEOld/D");
803 m_tree[i]->Branch("chiMuOld", &m_chimuOld, "chiMuOld/D");
804 m_tree[i]->Branch("chiPiOld", &m_chipiOld, "chiPiOld/D");
805 m_tree[i]->Branch("chiKOld", &m_chikOld, "chiKOld/D");
806 m_tree[i]->Branch("chiPOld", &m_chipOld, "chiPOld/D");
807 m_tree[i]->Branch("chiDOld", &m_chidOld, "chiDOld/D");
808
809 m_tree[i]->Branch("pmeanE", &m_pmeane, "pmeanE/D");
810 m_tree[i]->Branch("pmeanMu", &m_pmeanmu, "pmeanMu/D");
811 m_tree[i]->Branch("pmeanPi", &m_pmeanpi, "pmeanPi/D");
812 m_tree[i]->Branch("pmeanK", &m_pmeank, "pmeanK/D");
813 m_tree[i]->Branch("pmeanP", &m_pmeanp, "pmeanP/D");
814 m_tree[i]->Branch("pmeanD", &m_pmeand, "pmeanD/D");
815
816 m_tree[i]->Branch("presE", &m_prese, "presE/D");
817 m_tree[i]->Branch("presMu", &m_presmu, "presMu/D");
818 m_tree[i]->Branch("presPi", &m_prespi, "presPi/D");
819 m_tree[i]->Branch("presK", &m_presk, "presK/D");
820 m_tree[i]->Branch("presP", &m_presp, "presP/D");
821 m_tree[i]->Branch("presD", &m_presd, "presD/D");
822
823 // layer level information
824 m_tree[i]->Branch("lNHits", &l_nhits, "lNHits/I");
825 m_tree[i]->Branch("lNHitsUsed", &l_nhitsused, "lNHitsUsed/I");
826 m_tree[i]->Branch("lNHitsCombined", l_nhitscombined, "lNHitsCombined[lNHits]/I");
827 m_tree[i]->Branch("lWireLongestHit", l_wirelongesthit, "lWireLongestHit[lNHits]/I");
828 m_tree[i]->Branch("lLayer", l_layer, "lLayer[lNHits]/I");
829 m_tree[i]->Branch("lPath", l_path, "lPath[lNHits]/D");
830 m_tree[i]->Branch("lDedx", l_dedx, "lDedx[lNHits]/D");
831
832 // hit level information
833 if (m_isHitLevel) {
834 m_tree[i]->Branch("hNHits", &h_nhits, "hNHits/I");
835 m_tree[i]->Branch("hLWire", h_lwire, "hLWire[hNHits]/I");
836 m_tree[i]->Branch("hWire", h_wire, "hWire[hNHits]/I");
837 m_tree[i]->Branch("hLayer", h_layer, "hLayer[hNHits]/I");
838 m_tree[i]->Branch("hPath", h_path, "hPath[hNHits]/D");
839 m_tree[i]->Branch("hDedx", h_dedx, "hDedx[hNHits]/D");
840 m_tree[i]->Branch("hADCRaw", h_adcraw, "hADCRaw[hNHits]/D");
841 m_tree[i]->Branch("hADCCorr", h_adccorr, "hADCCorr[hNHits]/D");
842 m_tree[i]->Branch("hDoca", h_doca, "hDoca[hNHits]/D");
843 m_tree[i]->Branch("hNDoca", h_ndoca, "hNDoca[hNHits]/D");
844 m_tree[i]->Branch("hNDocaRS", h_ndocaRS, "hNDocaRS[hNHits]/D");
845 m_tree[i]->Branch("hEnta", h_enta, "hEnta[hNHits]/D");
846 m_tree[i]->Branch("hEntaRS", h_entaRS, "hEntaRS[hNHits]/D");
847 m_tree[i]->Branch("hDriftT", h_driftT, "hDriftT[hNHits]/D");
848 m_tree[i]->Branch("hDriftD", h_driftD, "hDriftD[hNHits]/D");
849 m_tree[i]->Branch("hFacnlADC", h_facnladc, "hFacnlADC[hNHits]/D");
850 m_tree[i]->Branch("hWireGain", h_wireGain, "hWireGain[hNHits]/D");
851 m_tree[i]->Branch("hTwodcor", h_twodCor, "hTwodcor[hNHits]/D");
852 m_tree[i]->Branch("hOnedcor", h_onedCor, "hOnedcor[hNHits]/D");
853
854 if (m_isExtraVar) {
855 m_tree[i]->Branch("hWeightPionHypo", h_WeightPionHypo, "hWeightPionHypo[hNHits]/D");
856 m_tree[i]->Branch("hWeightKaonHypo", h_WeightKaonHypo, "hWeightKaonHypo[hNHits]/D");
857 m_tree[i]->Branch("hWeightProtonHypo", h_WeightProtonHypo, "hWeightProtonHypo[hNHits]/D");
858 m_tree[i]->Branch("hFoundByTrackFinder", h_foundByTrackFinder, "hFoundByTrackFinder[hNHits]/I");
859 }
860 }
861}
double R
typedef autogenerated by FFTW
Debug output for CDCDedxPID module.
double getTimeMean() const
Return the injection gain for this track.
int getHitLayer(int i) const
Return the (global) layer number for a hit.
int getADCBaseCount(int i) const
Return the base adcCount (no non-linearity) for this hit.
int getADCCount(int i) const
Return the adcCount for this hit.
double getDedxMean() const
Get the dE/dx mean for this track.
double getDriftD(int i) const
Return the drift distance for this hit.
int getLayer(int i) const
Return the (global) layer number for a layer hit.
double getDoca(int i) const
Return the distance of closest approach to the sense wire for this hit.
double getDedx() const
Get dE/dx truncated mean for this track.
double getPath(int i) const
Return the path length through the cell for this hit.
int getWireInLayer(int i) const
Return the sensor ID for this hit: wire number in the layer.
double getPmean(int i) const
Return the PID (predicted mean) value.
double getOneDCorrection(int i) const
Return the 1D correction for this hit.
int getWireLongestHit(int i) const
Return the wire number of the longest hit per layer.
double getWeightKaonHypo(int i) const
Return the max weights from KalmanFitterInfo using kaon hypothesis.
void setDedx(double mean)
Set the dE/dx truncated average for this track.
double getCosineCorrection() const
Return the cosine correction for this track.
double getEntaRS(int i) const
Return rescaled enta value for cell height=width assumption.
double getLayerDedx(int i) const
Return the total dE/dx for this layer.
double getCellHalfWidth(int i) const
Return the half-width of the CDC cell.
double getPDG() const
Get the identity of the particle.
double getTwoDCorrection(int i) const
Return the 2D correction for this hit.
double getPreso(int i) const
Return the PID (predicted reso) value.
int getNLayerHits() const
Return the number of layer hits for this track.
int getFoundByTrackFinder(int i) const
Return the TrackFinder which added the given hit to track.
int getDriftT(int i) const
Return the drift time for this hit.
double getCosTheta() const
Return cos(theta) for this track.
double getInjectionRing() const
Return cos(theta) for this track.
int getCharge() const
Return the charge for this track.
double getWireGain(int i) const
Return the wire gain for this hit.
double getWeightPionHypo(int i) const
Return the max weights from KalmanFitterInfo using pion hypothesis.
int getWire(int i) const
Return the sensor ID for this hit: wire number for CDC (0-14336)
double getDocaRS(int i) const
Return rescaled doca value for cell height=width assumption.
double getCosEdgeCorrection() const
Return the cosine correction for this track.
double getRunGain() const
Return the run gain for this track.
double getNonLADCCorrection(int i) const
Return the factor introduce for adcCount (non-linearity) correction.
double getDedxError() const
Get the error on the dE/dx truncated mean for this track.
double getWeightProtonHypo(int i) const
Return the max weights from KalmanFitterInfo using proton hypothesis.
double getEnta(int i) const
Return the entrance angle in the CDC cell for this hit.
double getInjectionTime() const
Return cos(theta) for this track.
double trackID() const
Return the track ID.
double getLayerPath(int i) const
Return the distance travelled in this layer.
double getChi(int i) const
Return the PID (chi) value.
double getScaleFactor() const
Return the scale factor for this track.
double getLength() const
Return the total path length for this track.
int getNHitsCombined(int i) const
Return the number of hits combined per layer.
double getDedxNoSat() const
Get dE/dx truncated mean without the saturation correction for this track.
int size() const
Return the number of hits for this track.
double getNLayerHitsUsed() const
Return the number of hits used to determine the truncated mean.
double getMomentum() const
Return the track momentum valid in the CDC.
Provides a type-safe way to pass members of the chargedStableSet set.
Definition Const.h:589
static const unsigned int c_SetSize
Number of elements (for use in array bounds etc.)
Definition Const.h:615
A set of ParticleType objects, with defined order.
Definition Const.h:517
static const ChargedStable muon
muon particle
Definition Const.h:660
static const ParticleSet chargedStableSet
set of charged stable particles
Definition Const.h:618
static const ChargedStable pion
charged pion particle
Definition Const.h:661
static const ChargedStable proton
proton particle
Definition Const.h:663
static const ChargedStable kaon
charged kaon particle
Definition Const.h:662
static const ChargedStable electron
electron particle
Definition Const.h:659
static const ChargedStable deuteron
deuteron particle
Definition Const.h:664
Class for accessing objects in the database.
Definition DBObjPtr.h:21
ECL cluster data.
Definition ECLCluster.h:27
bool hasHypothesis(EHypothesisBit bitmask) const
Return if specific hypothesis bit is set.
Definition ECLCluster.h:351
double getE1oE9() const
Return E1/E9 (shower shape variable).
Definition ECLCluster.h:277
double getE9oE21() const
Return E9/E21 (shower shape variable).
Definition ECLCluster.h:280
double getEnergy(EHypothesisBit hypothesis) const
Return Energy (GeV).
Definition ECLCluster.cc:23
double getNumberOfCrystals() const
Return number of a crystals in a shower (sum of weights).
Definition ECLCluster.h:295
@ c_nPhotons
CR is split into n photons (N1)
Definition ECLCluster.h:41
double h_ndoca[kMaxHits]
normalized distance of closest approach
double h_driftD[kMaxHits]
drift distance
double h_onedCor[kMaxHits]
calibration 1D cleanup correction
int l_nhits
the total number of layer hits for this Track
double m_cosEdgeCor
calibration cosine edge correction
double m_presk
pred reso for kaon hypothesis
double m_phi0
Angle of the transverse momentum in the r-phi plane.
double m_d0
Signed distance to the POCA in the r-phi plane.
double m_presp
pred reso for proton hypothesis
double D2I(const double cosTheta, const double D) const
hadron saturation parameterization part 2
int l_wirelongesthit[kMaxHits]
the wire number of longest hit in this layer
double m_ioasym
asymmetry in increasing vs decreasing layer numbers per track
StoreObjPtr< EventLevelTriggerTimeInfo > m_TTDInfo
Store Object Ptr: EventLevelTriggerTimeInfo.
std::vector< double > m_hadronpars
hadron saturation parameters
double m_p
momentum from tracking
double m_e9_21
ratio of energies of the central 3x3 crystal vs 5x5 crystals
int m_trackID
ID number of the Track.
DBObjPtr< CDCDedxRunGain > m_DBRunGain
Run gain DB object.
double m_vx0
X coordinate of track POCA to origin.
virtual ~HitLevelInfoWriterModule()
Destructor.
StoreArray< KLMCluster > m_klmClusters
Required array of input KLMClusters.
double m_pmeand
pred mean for deuteron hypothesis
virtual void initialize() override
Initialize the module.
DBObjPtr< CDCDedxHadronCor > m_DBHadronCor
hadron saturation parameters
double h_enta[kMaxHits]
entrance angle
double m_chipOld
chi value for proton hypothesis
std::string m_strOutputBaseName
Base name for the output ROOT files.
virtual void event() override
This method is called for each event.
double h_WeightKaonHypo[kMaxHits]
weight for kaon hypothesis from KalmanFitterInfo
DBObjPtr< CDCDedxCosineEdge > m_DBCosEdgeCor
cosine edge calibration
HitLevelInfoWriterModule()
Default constructor.
double m_pmeanp
pred mean for proton hypothesis
double I2D(const double cosTheta, const double I) const
hadron saturation parameterization part 1
double h_WeightProtonHypo[kMaxHits]
weight for proton hypothesis from KalmanFitterInfo
double m_chid
modified chi value for deuteron hypothesis
double m_chikOld
chi value for kaon hypothesis
double m_dz
vertex or POCA in case of tracks z in respect to IPs
bool m_isHitLevel
Flag to switch on/off hit level info.
DBObjPtr< CDCDedxADCNonLinearity > m_DBNonlADC
hadron saturation non linearity
double m_pmeank
pred mean for kaon hypothesis
double h_facnladc[kMaxHits]
calibration hit gain
DBObjPtr< CDCDedx1DCell > m_DB1DCell
1D correction DB object
void saveChiValue(double(&chi)[Const::ChargedStable::c_SetSize], CDCDedxTrack *dedxTrack, double dedx) const
for all particles, save chi values into 'chi' chi array of chi values to be modified
StoreArray< TrackFitResult > m_trackFitResults
Required array of input TrackFitResults.
int h_nhits
the number of good hits for this Track
virtual void terminate() override
End of the event processing.
double h_adccorr[kMaxHits]
charge per hit corr by nonlinear ADC
std::vector< TTree * > m_tree
output ROOT trees
double m_presmu
pred reso for muon hypothesis
double m_cosCor
calibration cosine correction
double m_tanlambda
Slope of the track in the r-z plane.
int h_layer[kMaxHits]
layer number
double h_adcraw[kMaxHits]
charge per hit
int m_expID
experiment in which this Track was found
void fillTrack(const TrackFitResult *fitResult)
Fill the TTree with the information from the track fit.
double h_path[kMaxHits]
path length in cell
double m_chimu
modified chi value for muon hypothesis
double m_vy0
Y coordinate of track POCA to origin.
double m_chik
modified chi value for kaon hypothesis
void recalculateDedx(CDCDedxTrack *dedxTrack, std::map< int, std::vector< double > > &l_var, double(&cdcChi)[Const::ChargedStable::c_SetSize])
Function to recalculate the dedx with latest constants.
double m_chi2
chi^2 from track fit
void HadronCorrection(double costheta, double &dedx) const
Function to apply the hadron correction.
void calculateMeans(double *mean, double *truncatedMean, double *truncatedMeanErr, const std::vector< double > &dedx) const
Save arithmetic and truncated mean for the 'dedx' values.
double h_wireGain[kMaxHits]
calibration hit gain
double m_theta
cos(theta) for the track
double h_twodCor[kMaxHits]
calibration 2D correction
double m_chie
modified chi value for electron hypothesis
DBObjPtr< CDCDedx2DCell > m_DB2DCell
2D correction DB object
double m_runGain
calibration run gain
int m_eventID
event in which this Track was found
double m_timeGain
calibration injection time gain
void bookOutput(std::string filename)
Create the output TFiles and TTrees.
double m_chidOld
chi value for deuteron hypothesis
bool m_isExtraVar
Flag to switch on/off extra level info and some available w/ release/5 only.
void clearEntries()
Clear the arrays before filling an event.
std::vector< TFile * > m_file
output ROOT files
double m_dr
track d0 relative to IP
double m_prese
pred reso for electron hypothesis
DBObjPtr< CDCDedxCosineCor > m_DBCosineCor
Electron saturation correction DB object.
int h_foundByTrackFinder[kMaxHits]
the 'found by track finder' flag for the given hit
double m_timeReso
calibration injection time reso
StoreArray< Track > m_tracks
Required array of input Tracks.
double l_dedx[kMaxHits]
dE/dx for this layer
double l_path[kMaxHits]
distance travelled in this layer
int l_layer[kMaxHits]
layer number
double m_chipiOld
chi value for pion hypothesis
double m_chip
modified chi value for proton hypothesis
bool m_isCorrection
Flag to switch on/off corrections.
double m_pmeane
pred mean for electron hypothesis
double m_cosTheta
cos(theta) for the track
int l_nhitsused
the total number of layer hits used for this Track
StoreArray< ECLCluster > m_eclClusters
Required array of input ECLClusters.
double h_dedx[kMaxHits]
charge per path length
double m_chieOld
chi value for electron hypothesis
double h_driftT[kMaxHits]
drift time
double m_dphi
POCA in degrees in respect to IP.
DBObjPtr< CDCDedxWireGain > m_DBWireGains
Wire gain DB object.
double m_eop
energy over momentum in the calorimeter
bool m_isDeadwire
write only active wires
double m_presd
pred reso for deuteron hypothesis
double m_scale
calibration scale factor
double m_klmLayers
number of klm layers with hits
int l_nhitscombined[kMaxHits]
the number of hits combined this layer
double m_error
standard deviation of the truncated mean
double m_pCDC
momentum valid in CDC
double m_truncNoSat
dE/dx averaged, truncated mean, with corrections (not hadron)
double m_length
total path length of the Track
double m_injring
HER injection status.
double m_nCDChits
Number of CDC hits associated to the track.
double h_ndocaRS[kMaxHits]
normalized +RS distance of closest approach
double m_vz0
Z coordinate of track POCA to origin.
double m_trunc
dE/dx averaged, truncated mean, with corrections
double GetCorrection(int &adc, int layer, int wireID, double doca, double enta, double costheta, double ring, double time) const
Function to get the correction factor.
int m_inCDC
frack is CDC acceptance or not
double h_doca[kMaxHits]
distance of closest approach
double m_chimuOld
chi value for muon hypothesis
double m_z0
z coordinate of the POCA
int m_charge
the charge for this Track
double h_WeightPionHypo[kMaxHits]
weight for pion hypothesis from KalmanFitterInfo
int m_runID
run in which this Track was found
double m_e
energy in the calorimeter
StoreArray< CDCDedxTrack > m_dedxTracks
Required array of CDCDedxTracks.
double m_injtime
time since last injection in micro seconds
double m_pmeanmu
pred mean for muon hypothesis
double h_entaRS[kMaxHits]
normalized + RS distance of entrance angle
bool m_isRelative
Flag to switch on/off relative constants.
double m_prespi
pred reso for pion hypothesis
DBObjPtr< CDCDedxInjectionTime > m_DBInjectTime
time gain/reso DB object
int h_wire[kMaxHits]
sense wire ID
DBObjPtr< CDCDedxScaleFactor > m_DBScaleFactor
Scale factor to make electrons ~1.
int h_lwire[kMaxHits]
sense wire within layer
double m_chipi
modified chi value for pion hypothesis
std::vector< std::string > m_strParticleList
Vector of ParticleLists to write out.
double m_pt
transverse momentum from tracking
double m_pmeanpi
pred mean for pion hypothesis
void fillDedx(CDCDedxTrack *dedxTrack)
Fill the TTree with the information from a CDCDedxTrack object.
double m_e1_9
ratio of energies of the central 1 crystal vs 3x3 crystals
unsigned short getNHits() const
Get the total Number of CDC hits in the fit.
KLM cluster data.
Definition KLMCluster.h:29
int getLayers() const
Get number of layers with hits.
Definition KLMCluster.h:67
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
Module()
Constructor.
Definition Module.cc:30
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Class to store reconstructed particles.
Definition Particle.h:76
static const ReferenceFrame & GetCurrent()
Get current rest frame.
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
FROM * getRelatedFrom(const std::string &name="", const std::string &namedRelation="") const
Get the object from which this object has a relation.
Type-safe access to single objects in the data store.
Definition StoreObjPtr.h:96
Values of the result of a track fit with a given particle hypothesis.
short getChargeSign() const
Return track charge (1 or -1).
double getPValue() const
Getter for Chi2 Probability of the track fit.
double getD0() const
Getter for d0.
double getTanLambda() const
Getter for tanLambda.
double getZ0() const
Getter for z0.
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
ROOT::Math::XYZVector getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
HitPatternCDC getHitPatternCDC() const
Getter for the hit pattern in the CDC;.
UncertainHelix getUncertainHelix() const
Conversion to framework Uncertain Helix (i.e., with covariance).
double getPhi0() const
Getter for phi0.
Class that bundles various TrackFitResults.
Definition Track.h:25
This class represents an ideal helix in perigee parameterization including the covariance matrix of t...
double passiveMoveBy(const ROOT::Math::XYZVector &by)
Moves origin of the coordinate system (passive transformation) by the given vector.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.
STL namespace.