Belle II Software development
CDCDedxElectronCollectorModule.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/CDCDedxElectronCollector/CDCDedxElectronCollectorModule.h>
10
11#include <mdst/dataobjects/ECLCluster.h>
12#include <TTree.h>
13#include <TH1D.h>
14#include <TH1I.h>
15#include <TMath.h>
16
17#include <cmath>
18
19using namespace Belle2;
20//-----------------------------------------------------------------
21// Register the Module
22//-----------------------------------------------------------------
23REG_MODULE(CDCDedxElectronCollector);
24
25//-----------------------------------------------------------------
26// Implementation
27//-----------------------------------------------------------------
29{
30 // Set module properties
31 setDescription("A collector module for CDC dE/dx electron calibrations");
32
33 // Parameter definitions
34
35 addParam("cleanupCuts", m_cuts, "boolean to apply cleanup cuts", true);
36 addParam("maxHits", m_maxHits, "maximum number of hits per track ", int(100));
37 addParam("setEoP", m_setEoP, "Set E over p Cut values. ", double(0.25));
38 addParam("isCosth", m_isCosth, "true for adding costh tree branch. ", false);
39 addParam("isMom", m_isMom, "true for adding momentum tree branch. ", false);
40 addParam("isPt", m_isPt, "true for adding transverse momentum tree branch. ", false);
41 addParam("isCharge", m_isCharge, "true for charge dedx tree branch. ", false);
42 addParam("isRun", m_isRun, "true for adding run number tree branch. ", false);
43 addParam("isWire", m_isWire, "true for adding wires tree branch. ", false);
44 addParam("isLayer", m_isLayer, "true for adding layers tree branch. ", false);
45 addParam("isDoca", m_isDoca, "true for adding doca tree branch. ", false);
46 addParam("isEnta", m_isEnta, "true for adding enta tree branch. ", false);
47 addParam("isInjTime", m_isInjTime, "true for adding time var tree branch. ", false);
48 addParam("isDocaRS", m_isDocaRS, "true for adding doca tree branch. ", false);
49 addParam("isEntaRS", m_isEntaRS, "true for adding enta tree branch. ", false);
50 addParam("isDedxhit", m_isDedxhit, "true for adding dedxhit tree branch. ", false);
51 addParam("isADCcorr", m_isADCcorr, "true for adding adc tree branch. ", false);
52 addParam("isBhabha", m_isBhabha, "true for bhabha events", true);
53 addParam("isRadee", m_isRadee, "true for radee events", false);
54 addParam("isTrgSel", m_isTrgSel, "true to enable trigger sel inside module", false);
55}
56
57//-----------------------------------------------------------------
58// Create ROOT objects
59//-----------------------------------------------------------------
60
62{
63 m_trgResult.isOptional();
64 m_dedxTracks.isRequired();
65 m_tracks.isRequired();
66 m_trackFitResults.isRequired();
67
68 // Data object creation
69 auto means = new TH1D("means", "CDC dE/dx truncated means", 250, 0, 2.5);
70 auto ttree = new TTree("tree", "Tree with dE/dx information");
71
72 auto hestats = new TH1I("hestats", "Event Stats", 6, -0.5, 5.5);
73 hestats->SetFillColor(kYellow);
74 hestats->GetXaxis()->SetBinLabel(1, "all");
75 hestats->GetXaxis()->SetBinLabel(2, "notrig");
76 hestats->GetXaxis()->SetBinLabel(3, "noskim");
77 hestats->GetXaxis()->SetBinLabel(4, "wrongskim");
78 hestats->GetXaxis()->SetBinLabel(5, "unclean");
79 hestats->GetXaxis()->SetBinLabel(6, "selected");
80
81 auto htstats = new TH1I("htstats", "track Stats", 7, -0.5, 6.5);
82 htstats->SetFillColor(kYellow);
83 htstats->GetXaxis()->SetBinLabel(1, "alltrk");
84 htstats->GetXaxis()->SetBinLabel(2, "vtx");
85 htstats->GetXaxis()->SetBinLabel(3, "inCDC");
86 htstats->GetXaxis()->SetBinLabel(4, "whits");
87 htstats->GetXaxis()->SetBinLabel(5, "weop");
88 htstats->GetXaxis()->SetBinLabel(6, "radee");
89 htstats->GetXaxis()->SetBinLabel(7, "selected");
90
91 if (m_isInjTime) {
92 ttree->Branch<double>("injtime", &m_injTime);
93 ttree->Branch<double>("injring", &m_injRing);
94 ttree->Branch<double>("costh", &m_costh);
95 ttree->Branch<int>("nhits", &m_nhits);
96 ttree->Branch<double>("p", &m_p);
97 }
98
99 ttree->Branch<double>("dedx", &m_dedx);
100 if (m_isCosth)ttree->Branch<double>("costh", &m_costh);
101 if (m_isMom)ttree->Branch<double>("p", &m_p);
102 if (m_isPt)ttree->Branch<double>("pt", &m_pt);
103 if (m_isCharge)ttree->Branch<int>("charge", &m_charge);
104 if (m_isRun)ttree->Branch<int>("run", &m_run);
105
106 if (m_isWire)ttree->Branch("wire", &m_wire);
107 if (m_isLayer)ttree->Branch("layer", &m_layer);
108 if (m_isDoca)ttree->Branch("doca", &m_doca);
109 if (m_isEnta)ttree->Branch("enta", &m_enta);
110 if (m_isDocaRS)ttree->Branch("docaRS", &m_docaRS);
111 if (m_isEntaRS)ttree->Branch("entaRS", &m_entaRS);
112 if (m_isDedxhit)ttree->Branch("dedxhit", &m_dedxhit);
113 if (m_isADCcorr)ttree->Branch("adccorr", &m_adccorr);
114
115 // Collector object registration
116 registerObject<TH1D>("means", means);
117 registerObject<TTree>("tree", ttree);
118 registerObject<TH1I>("hestats", hestats);
119 registerObject<TH1I>("htstats", htstats);
120}
121
122//-----------------------------------------------------------------
123// Fill ROOT objects
124//-----------------------------------------------------------------
126{
127
128 auto hestats = getObjectPtr<TH1I>("hestats");
129 hestats->Fill(0);
130
131 if (m_isTrgSel) {
132 if (!m_trgResult.isValid()) {
133 B2WARNING("SoftwareTriggerResult required to select bhabha/radee event is not found");
134 hestats->Fill(1);
135 return;
136 }
137
138 //release05: bhabha_all is grand skim = bhabha+bhabhaecl+radee
139 const std::map<std::string, int>& fresults = m_trgResult->getResults();
140 if (fresults.find("software_trigger_cut&skim&accept_bhabha") == fresults.end() and
141 fresults.find("software_trigger_cut&skim&accept_radee") == fresults.end()) {
142 B2WARNING("Can't find required bhabha/radee trigger identifiers");
143 hestats->Fill(2);
144 return;
145 }
146
147 const bool eBhabha = (m_trgResult->getResult("software_trigger_cut&skim&accept_bhabha") ==
149
150 const bool eRadBhabha = (m_trgResult->getResult("software_trigger_cut&skim&accept_radee") ==
152
153 if (!m_isBhabha && !m_isRadee) {
154 B2WARNING("requested not-supported event type: going back");
155 hestats->Fill(3);
156 return;
157 } else if (m_isBhabha && !m_isRadee && !eBhabha) {
158 B2WARNING("requested bhabha only but event not found: going back");
159 hestats->Fill(3);
160 return;
161 } else if (m_isRadee && !m_isBhabha && !eRadBhabha) {
162 B2WARNING("requested radee only but event not found: going back");
163 hestats->Fill(3);
164 return;
165 }
166 } else {
167 hestats->GetXaxis()->SetBinLabel(2, "inact1");
168 hestats->GetXaxis()->SetBinLabel(3, "inact2");
169 hestats->GetXaxis()->SetBinLabel(4, "inact3");
170 }
171
172 StoreObjPtr<EventMetaData> eventMetaDataPtr;
173 int run = eventMetaDataPtr->getRun();
174 if (m_isRun)m_run = run;
175 int nTracks = m_dedxTracks.getEntries();
176 if (nTracks >= 4) {
177 B2WARNING("too many tracks: unclean bhabha or radee event: " << nTracks);
178 hestats->Fill(4);
179 return;
180 }
181
182 hestats->Fill(5);
183
184 //Collector object access
185 auto tree = getObjectPtr<TTree>("tree");
186 auto htstats = getObjectPtr<TH1I>("htstats");
187 auto hmeans = getObjectPtr<TH1D>("means");
188
189 for (int idedx = 0; idedx < nTracks; idedx++) {
190
191 CDCDedxTrack* dedxTrack = m_dedxTracks[idedx];
192 if (!dedxTrack) {
193 B2WARNING("No dedx track: Going back: " << idedx);
194 continue;
195 }
196
197 const Track* track = dedxTrack->getRelatedFrom<Track>();
198 if (!track) {
199 B2WARNING("No track: Going back: " << idedx);
200 continue;
201 }
202
203 const TrackFitResult* fitResult = track->getTrackFitResultWithClosestMass(Const::pion);
204 if (!fitResult) {
205 B2WARNING("No related fit for this track...");
206 continue;
207 }
208
209 m_dedx = dedxTrack->getDedxNoSat();
210 m_p = dedxTrack->getMomentum();
211 m_costh = dedxTrack->getCosTheta();
212 m_charge = fitResult->getChargeSign();
213 m_pt = fitResult->getTransverseMomentum();
214 m_injTime = dedxTrack->getInjectionTime();
215 m_injRing = dedxTrack->getInjectionRing();
216 m_nhits = dedxTrack->size();
217
218 htstats->Fill(0);
219
220 if (m_cuts) {
221 // apply cleanup cuts
222 if (std::abs(fitResult->getD0()) >= 1.0)continue;
223 if (std::abs(fitResult->getZ0()) >= 1.0) continue;
224 htstats->Fill(1);
225
226 //if outside CDC
227 if (m_costh < TMath::Cos(150.0 * TMath::DegToRad()))continue; //-0.866
228 if (m_costh > TMath::Cos(17.0 * TMath::DegToRad())) continue; //0.95
229 htstats->Fill(2);
230
231 if (m_nhits > m_maxHits) continue;
232
233 //making some cuts based on acceptance
234 if (m_costh > -0.55 && m_costh < 0.820) {
235 if (dedxTrack->getNLayerHits() < 25)continue; //all CDC layer available here
236 } else {
238 if (dedxTrack->getNLayerHits() < 10)continue; //less layer available here
239 if (m_costh > 0 && dedxTrack->getNLayerHits() < 13)continue;
240 } else {
241 if (dedxTrack->getNLayerHits() < 18)continue;
242 }
243 }
244 htstats->Fill(3);
245
246 const ECLCluster* eclCluster = track->getRelated<ECLCluster>();
247 if (eclCluster and eclCluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
248 double TrkEoverP = (eclCluster->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons)) / (fitResult->getMomentum().R());
249 if (std::abs(TrkEoverP - 1.0) > m_setEoP)continue;
250 }
251 htstats->Fill(4);
252 }
253
254 //if dealing with radee here (do a safe side cleanup)
255 if (m_isRadee) {
256 if (nTracks != 2)continue; //exactly 2 tracks
257 bool goodradee = false;
258 //checking if dedx of other track is restricted
259 //will not do too much as radee is clean enough
260 for (int jdedx = 0; jdedx < nTracks; jdedx++) {
261 CDCDedxTrack* dedxOtherTrack = m_dedxTracks[std::abs(jdedx - 1)];
262 if (!dedxOtherTrack)continue;
263 if (std::abs(dedxOtherTrack->getDedxNoSat() - 1.0) > 0.25)continue; //loose for uncalibrated
264 goodradee = true;
265 break;
266 }
267 if (!goodradee)continue;
268 htstats->Fill(5);
269 }
270
271
272 // Make sure to remove all the data in vectors from the previous track
273 if (m_isWire)m_wire.clear();
274 if (m_isLayer)m_layer.clear();
275 if (m_isDoca)m_doca.clear();
276 if (m_isEnta)m_enta.clear();
277 if (m_isDocaRS)m_docaRS.clear();
278 if (m_isEntaRS)m_entaRS.clear();
279 if (m_isDedxhit)m_dedxhit.clear();
280 if (m_isADCcorr)m_adccorr.clear();
281
282 // Simple numbers don't need to be cleared
283 // make sure to use the truncated mean without the hadron saturation correction
284
285 for (int i = 0; i < m_nhits; ++i) {
286 // if (m_DBWireGains->getWireGain(dedxTrack->getWire(i)) == 0)continue;
287 if (m_isWire)m_wire.push_back(dedxTrack->getWire(i));
288 if (m_isLayer)m_layer.push_back(dedxTrack->getHitLayer(i));
289 if (m_isDoca)m_doca.push_back(dedxTrack->getDoca(i));
290 if (m_isEnta)m_enta.push_back(dedxTrack->getEnta(i));
291 if (m_isDocaRS)m_docaRS.push_back(dedxTrack->getDocaRS(i) / dedxTrack->getCellHalfWidth(i));
292 if (m_isEntaRS)m_entaRS.push_back(dedxTrack->getEntaRS(i));
293 if (m_isDedxhit)m_dedxhit.push_back(dedxTrack->getDedx(i));
294 if (m_isADCcorr)m_adccorr.push_back(dedxTrack->getADCCount(i));
295 }
296
297 // Track and/or hit information filled as per config
298 htstats->Fill(6);
299 hmeans->Fill(m_dedx);
300 tree->Fill();
301 }
302}
double R
typedef autogenerated by FFTW
bool m_isEntaRS
flag to write rescaled enta in tree
bool m_isLayer
flag to write layer number in tree
bool m_isPt
flag to write trans momentum in treet
bool m_isADCcorr
flag to write adc corrected in tree
bool m_isWire
flag to write wire number in tree
std::vector< double > m_adccorr
adc corrected for the hit
std::vector< int > m_wire
hit level information
std::vector< int > m_layer
continuous layer number for the hit
StoreArray< TrackFitResult > m_trackFitResults
Required array for TrackFitResults.
bool m_isTrgSel
flag to enable trigger skim selected in the module (off default)
StoreObjPtr< SoftwareTriggerResult > m_trgResult
required input
std::vector< double > m_dedxhit
dE/dx for the hit
std::vector< double > m_doca
distance of closest approach for the hit
int m_nhits
number of dE/dx hits on the track
virtual void collect() override
Fill ROOT objects.
StoreArray< Track > m_tracks
Required array for Tracks.
virtual void prepare() override
Create and book ROOT objects.
CDCDedxElectronCollectorModule()
Constructor: Sets the description, the properties and the parameters of the module.
std::vector< double > m_entaRS
rescaled entrance angle for the hit
bool m_isInjTime
flag to enable trigger skim (off default)
std::vector< double > m_enta
entrance angle for the hit
StoreArray< CDCDedxTrack > m_dedxTracks
Required array for CDCDedxTracks.
bool m_isDocaRS
flag to write rescaled doca in tree
std::vector< double > m_docaRS
rescaled distance of closest approach for the hit
Debug output for CDCDedxPID module.
int getHitLayer(int i) const
Return the (global) layer number for a hit.
int getADCCount(int i) const
Return the adcCount for this 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 getEntaRS(int i) const
Return rescaled enta value for cell height=width assumption.
double getCellHalfWidth(int i) const
Return the half-width of the CDC cell.
int getNLayerHits() const
Return the number of layer hits for this track.
double getCosTheta() const
Return cos(theta) for this track.
double getInjectionRing() const
Return cos(theta) for this track.
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 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 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 getMomentum() const
Return the track momentum valid in the CDC.
void registerObject(std::string name, T *obj)
Register object with a name, takes ownership, do not access the pointer beyond prepare()
CalibrationCollectorModule()
Constructor. Sets the default prefix for calibration dataobjects.
T * getObjectPtr(std::string name)
Calls the CalibObjManager to get the requested stored collector data.
static const ChargedStable pion
charged pion particle
Definition Const.h:661
ECL cluster data.
Definition ECLCluster.h:27
bool hasHypothesis(EHypothesisBit bitmask) const
Return if specific hypothesis bit is set.
Definition ECLCluster.h:351
double getEnergy(EHypothesisBit hypothesis) const
Return Energy (GeV).
Definition ECLCluster.cc:23
@ c_nPhotons
CR is split into n photons (N1)
Definition ECLCluster.h:41
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
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 getD0() const
Getter for d0.
double getTransverseMomentum() const
Getter for the absolute value of the transverse momentum at the perigee.
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.
Class that bundles various TrackFitResults.
Definition Track.h:25
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
Abstract base class for different kinds of events.