Belle II Software  release-05-01-25
NoKickRTSel.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Valerio Bertacchi *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <tracking/trackFindingVXD/sectorMapTools/NoKickRTSel.h>
12 #include <tracking/trackFindingVXD/sectorMapTools/NoKickCuts.h>
13 
14 #include <tracking/dataobjects/hitXPDerivate.h>
15 
16 #include <mdst/dataobjects/MCParticle.h>
17 #include <pxd/dataobjects/PXDTrueHit.h>
18 #include <svd/dataobjects/SVDTrueHit.h>
19 
20 using namespace Belle2;
21 
22 
23 void NoKickRTSel::hitXPBuilder(const RecoTrack& track)
24 {
25  StoreArray<SVDCluster> SVDClusters;
26  StoreArray<SVDTrueHit> SVDTrueHits;
27  StoreArray<MCParticle> MCParticles;
28  StoreArray<RecoTrack> recoTracks;
29 
30  if (track.getRelationsTo<MCParticle>().size() > 0) {
31 
32  const MCParticle* particle = track.getRelationsTo<MCParticle>()[0];
33 
34  std::vector<Belle2::RecoHitInformation::UsedSVDHit*> clusterListSVD = track.getSVDHitList();
35  for (const SVDCluster* cluster : clusterListSVD) {
36  for (const SVDTrueHit& hit : cluster->getRelationsTo<SVDTrueHit>()) {
37  RelationVector<MCParticle> hitFromParticle = hit.getRelationsFrom<MCParticle>();
38  if (hitFromParticle.size() > 0) {
39  if (hitFromParticle[0]->getIndex() != particle->getIndex()) {
40  continue;
41  }
42  } else continue;
43  VxdID trueHitSensorID = hit.getSensorID();
44  const VXD::SensorInfoBase& sensorInfo = VXD::GeoCache::getInstance().getSensorInfo(trueHitSensorID);
45  hitXPDerivate entry(hit, *cluster, *particle, sensorInfo);
46  int NClusterU = 0;
47  int NClusterV = 0;
48  for (SVDCluster Ncluster : hit.getRelationsFrom<SVDCluster>()) {
49  if (Ncluster.isUCluster()) NClusterU++;
50  else NClusterV++;
51  }
52  entry.setClusterU(NClusterU);
53  entry.setClusterV(NClusterV);
54 
55  bool isReconstructed(false);
56  for (const RecoTrack& aRecoTrack : particle->getRelationsFrom<RecoTrack>())
57  isReconstructed |= aRecoTrack.hasSVDHits();
58  entry.setReconstructed(isReconstructed);
59 
60  m_setHitXP.insert(entry);
61  }
62  }
63 
64 
65  StoreArray<PXDCluster> PXDClusters;
66  StoreArray<PXDTrueHit> PXDTrueHits;
67 
68  std::vector<Belle2::RecoHitInformation::UsedPXDHit*> clusterListPXD = track.getPXDHitList();
69  for (const PXDCluster* cluster : clusterListPXD) {
70  for (const PXDTrueHit& hit : cluster->getRelationsTo<PXDTrueHit>()) {
71  RelationVector<MCParticle> hitFromParticle = hit.getRelationsFrom<MCParticle>();
72  if (hitFromParticle.size() > 0) {
73  if (hitFromParticle[0]->getIndex() != particle->getIndex()) {
74  continue;
75  }
76  } else continue;
77  VxdID trueHitSensorID = hit.getSensorID();
78  const VXD::SensorInfoBase& sensorInfo = VXD::GeoCache::getInstance().getSensorInfo(trueHitSensorID);
79  hitXPDerivate entry(hit, *particle, sensorInfo);
80 
81  bool isReconstructed(false);
82  for (const RecoTrack& aRecoTrack : particle->getRelationsFrom<RecoTrack>())
83  isReconstructed |= aRecoTrack.hasPXDHits();
84  entry.setReconstructed(isReconstructed);
85 
86  m_setHitXP.insert(entry);
87  }
88  }
89 
90  for (auto element : m_setHitXP) {
91  m_hitXP.push_back(element);
92  }
93 
94  }
95 
96 }
97 
98 void NoKickRTSel::hit8TrackBuilder(const RecoTrack& track)
99 {
100  hitXPBuilder(track);
101  for (const hitXP XP : m_hitXP) {
102  if (m_8hitTrack.size() < 1) {
103  m_8hitTrack.push_back(XP);
104  } else if (XP.m_sensorLayer != m_8hitTrack.back().m_sensorLayer ||
105  XP.m_sensorLadder != m_8hitTrack.back().m_sensorLadder ||
106  XP.m_sensorSensor != m_8hitTrack.back().m_sensorSensor) {
107  m_8hitTrack.push_back(XP);
108  }
109  }
110 
111 }
112 
113 bool NoKickRTSel::globalCut(const std::vector<hitXP>& track8)
114 {
115  int flagd0 = 1;
116  int flagz0 = 1;
117  int lay3 = 0;
118  int lay4 = 0;
119  int lay5 = 0;
120  int lay6 = 0;
121  for (hitXP XP : track8) {
122  if (XP.getSensorLayer() == 3) lay3 = 1;
123  if (XP.getSensorLayer() == 4) lay4 = 1;
124  if (XP.getSensorLayer() == 5) lay5 = 1;
125  if (XP.getSensorLayer() == 6) lay6 = 1;
126  // if (fabs(XP.getD0Entry()) > 1.) flagd0 = 0;
127  // if (fabs(XP.getZ0Entry()) > 1.) flagz0 = 0;
128  }
129  int N_lay = lay3 + lay4 + lay5 + lay6;
130  if (N_lay >= 3) N_lay = 1;
131  else N_lay = 0;
132  int flagTot = flagd0 * flagz0 * N_lay;
133  if (flagTot == 1) return true;
134  else return false;
135 }
136 
137 bool NoKickRTSel::segmentSelector(hitXP hit1, hitXP hit2, std::vector<double> selCut, NoKickCuts::EParameters par, bool is0)
138 {
139  if (hit2.m_sensorLayer - hit1.m_sensorLayer > 1) return true;
140  else {
141  double deltaPar = 0;
142  switch (par) {
143  case NoKickCuts::c_Omega:
144  deltaPar = fabs(hit1.getOmegaEntry() - hit2.getOmegaEntry());
145  if (is0) deltaPar = fabs(hit1.getOmega0() - hit2.getOmegaEntry());
146  // selCutPXD =0.4;
147  break;
148 
149  case NoKickCuts::c_D0:
150  deltaPar = hit1.getD0Entry() - hit2.getD0Entry();
151  if (is0) deltaPar = hit1.getD00() - hit2.getD0Entry();
152  // selCutPXD =1;
153  break;
154 
155  case NoKickCuts::c_Phi0:
156  deltaPar = asin(sin(hit1.getPhi0Entry())) - asin(sin(hit2.getPhi0Entry()));
157  if (is0) deltaPar = asin(sin(hit1.getPhi00())) - asin(sin(hit2.getPhi0Entry()));
158  // selCutPXD =0.3;
159  break;
160 
161  case NoKickCuts::c_Z0:
162  deltaPar = hit1.getZ0Entry() - hit2.getZ0Entry();
163  if (is0) deltaPar = hit1.getZ00() - hit2.getZ0Entry();
164  // selCutPXD =1;
165  break;
166 
167  case NoKickCuts::c_Tanlambda:
168  deltaPar = hit1.getTanLambdaEntry() - hit2.getTanLambdaEntry();
169  if (is0) deltaPar = hit1.getTanLambda0() - hit2.getTanLambdaEntry();
170  // selCutPXD =0.3;
171  break;
172  }
173 
174  double usedCut = 0;
175  if (fabs(selCut.at(0)) > fabs(selCut.at(1))) {
176  usedCut = fabs(selCut.at(0));
177  } else usedCut = fabs(selCut.at(1));
178 
179  if (deltaPar > -usedCut && deltaPar < usedCut) return true;
180  else {
181  B2DEBUG(20, "--------------------------");
182  B2DEBUG(20, "lay1=" << hit1.m_sensorLayer);
183  B2DEBUG(20, "lay2=" << hit2.m_sensorLayer);
184  B2DEBUG(20, "parameter=" << par);
185  B2DEBUG(20, "Min=" << selCut.at(0));
186  B2DEBUG(20, "Max=" << selCut.at(1));
187  B2DEBUG(20, "deltaPar=" << deltaPar);
188  B2DEBUG(20, "momentum=" << hit1.m_momentum0.Mag());
189  return false;
190  }
191  }
192 }
193 
194 
195 
196 bool NoKickRTSel::trackSelector(const RecoTrack& track)
197 {
198  initNoKickRTSel();
199  hit8TrackBuilder(track);
200 
201  if (m_outputFlag) {
202  m_pdgID = m_8hitTrack[0].getPDGID();
203  m_pMag = track.getMomentumSeed().Mag();
204  m_pt = sqrt(track.getMomentumSeed().X() * track.getMomentumSeed().X() + track.getMomentumSeed().Y() * track.getMomentumSeed().Y());
205  }
206 
207  bool good = globalCut(m_8hitTrack);
208  if (good == false) {
209  if (m_outputFlag) {
210  m_numberOfCuts = -1; //it means "no specific cuts applied, but rejected for global cuts"
211  m_Ncuts = m_numberOfCuts;
212  m_isCutted = true;
213  m_nCutHit->Fill(m_numberOfCuts);
214  m_momCut->Fill(track.getMomentumSeed().Mag());
215  m_PDGIDCut->Fill(track.getRelationsTo<MCParticle>()[0]->getPDG());
216  m_noKickTree->Fill();
217  }
218  return false;
219  }
220 
221  if (track.getMomentumSeed().Mag() > m_pmax) {
222  if (m_outputFlag) {
223  m_Ncuts = m_numberOfCuts;
224  m_isCutted = false;
225  m_nCutHit->Fill(m_numberOfCuts);
226  m_momSel->Fill(track.getMomentumSeed().Mag());
227  m_PDGIDSel->Fill(track.getRelationsTo<MCParticle>()[0]->getPDG());
228  m_noKickTree->Fill();
229  }
230  return good;
231  }
232  for (int i = 0; i < (int)(m_8hitTrack.size() - 2); i++) {
233  double sinTheta = fabs(m_8hitTrack.at(i).m_momentum0.Y()) /
234  sqrt(pow(m_8hitTrack.at(i).m_momentum0.Y(), 2) +
235  pow(m_8hitTrack.at(i).m_momentum0.Z(), 2));
236 
237  double momentum =
238  sqrt(pow(m_8hitTrack.at(i).m_momentum0.X(), 2) +
239  pow(m_8hitTrack.at(i).m_momentum0.Y(), 2) +
240  pow(m_8hitTrack.at(i).m_momentum0.Z(), 2));
241 
242  for (int j = NoKickCuts::c_Omega; j <= NoKickCuts::c_Tanlambda; j++) { //track parameters loop
243  std::vector<double> selCut = m_trackCuts.cutSelector(sinTheta, momentum, m_8hitTrack.at(i).m_sensorLayer,
244  m_8hitTrack.at(i + 1).m_sensorLayer, (NoKickCuts::EParameters) j);
245  bool goodSeg = segmentSelector(m_8hitTrack.at(i), m_8hitTrack.at(i + 1), selCut, (NoKickCuts::EParameters) j);
246 
247  if (!goodSeg) {
248  good = false;
249  m_numberOfCuts++;
250  }
251  if (i == 0) { //beampipe crossing
252  bool goodSeg0 = segmentSelector(m_8hitTrack.at(i), m_8hitTrack.at(i), selCut, (NoKickCuts::EParameters) j, true);
253  if (!goodSeg0) {
254  good = false;
255  m_numberOfCuts++;
256  }
257  }
258 
259  }
260  }
261  if (m_outputFlag) {
262  m_Ncuts = m_numberOfCuts;
263  m_nCutHit->Fill(m_numberOfCuts);
264  if (good) {
265  m_isCutted = false;
266  m_momSel->Fill(track.getMomentumSeed().Mag());
267  m_PDGIDSel->Fill(track.getRelationsTo<MCParticle>()[0]->getPDG());
268  } else {
269  m_isCutted = true;
270  m_momCut->Fill(track.getMomentumSeed().Mag());
271  m_PDGIDCut->Fill(track.getRelationsTo<MCParticle>()[0]->getPDG());
272  }
273  m_noKickTree->Fill();
274  }
275  return good;
276 }
277 
278 
279 void NoKickRTSel::produceHistoNoKick()
280 {
281  if (m_outputFlag) {
282  m_noKickOutputTFile->cd();
283  m_momSel->Write();
284  m_momCut->Write();
285 
286  m_momEff->Add(m_momSel, 1);
287  m_momEff->Add(m_momCut, 1);
288  m_momEff->Divide(m_momSel, m_momEff, 1, 1);
289  m_momEff->Write();
290 
291  m_PDGIDSel->Write();
292  m_PDGIDCut->Write();
293 
294  m_PDGIDEff->Add(m_PDGIDSel, 1);
295  m_PDGIDEff->Add(m_PDGIDCut, 1);
296  m_PDGIDEff->Divide(m_PDGIDSel, m_PDGIDEff, 1, 1);
297  m_PDGIDEff->Write();
298 
299  m_nCutHit->Write();
300 
301  m_noKickTree->Write();
302 
303  delete m_noKickOutputTFile;
304  }
305 }
Belle2::RelationVector::size
size_t size() const
Get number of relations.
Definition: RelationVector.h:98
Belle2::hitXP::getPhi00
double getPhi00() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:377
Belle2::NoKickRTSel::m_momEff
TH1F * m_momEff
histogram for efficiency
Definition: NoKickRTSel.h:57
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::NoKickRTSel::m_outputFlag
bool m_outputFlag
true=produce validation output
Definition: NoKickRTSel.h:52
Belle2::NoKickRTSel::m_noKickTree
TTree * m_noKickTree
TTree to which the information is written.
Definition: NoKickRTSel.h:67
Belle2::NoKickRTSel::trackSelector
bool trackSelector(const RecoTrack &track)
This method return true if every segment (see segmentSelector) of the input track respects the cuts c...
Definition: NoKickRTSel.cc:196
Belle2::NoKickRTSel::m_trackCuts
NoKickCuts m_trackCuts
auxiliary member to apply the cuts
Definition: NoKickRTSel.h:49
Belle2::PXDTrueHit
Class PXDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition: PXDTrueHit.h:35
Belle2::NoKickRTSel::m_setHitXP
std::set< hitXP, hitXP::timeCompare > m_setHitXP
set of hit to order the hit in time
Definition: NoKickRTSel.h:47
Belle2::NoKickRTSel::m_isCutted
bool m_isCutted
Indicator if cut is applied.
Definition: NoKickRTSel.h:62
Belle2::NoKickRTSel::m_pdgID
double m_pdgID
pdg Code
Definition: NoKickRTSel.h:65
Belle2::NoKickRTSel::m_momSel
TH1F * m_momSel
histogram of selected tracks
Definition: NoKickRTSel.h:55
Belle2::SVDTrueHit
Class SVDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition: SVDTrueHit.h:35
Belle2::NoKickCuts::EParameters
EParameters
enum for parameters name
Definition: NoKickCuts.h:54
Belle2::VXD::SensorInfoBase
Base class to provide Sensor Information for PXD and SVD.
Definition: SensorInfoBase.h:40
Belle2::hitXP::getOmegaEntry
double getOmegaEntry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:335
Belle2::hitXP::m_sensorLayer
int m_sensorLayer
layer of the hit
Definition: hitXP.h:62
Belle2::NoKickCuts::cutSelector
std::vector< double > cutSelector(double sintheta, double momentum, int layer1, int layer2, EParameters par)
This methods selects 2 cuts (minimum and maximum inside a vector) from the information of theta,...
Definition: NoKickCuts.cc:20
Belle2::NoKickRTSel::m_Ncuts
int m_Ncuts
number of times the cut is applied on a particle
Definition: NoKickRTSel.h:66
Belle2::NoKickRTSel::m_noKickOutputTFile
TFile * m_noKickOutputTFile
validartion output TFile
Definition: NoKickRTSel.h:54
Belle2::NoKickRTSel::hit8TrackBuilder
void hit8TrackBuilder(const RecoTrack &track)
this metod build a vector of hitXP from a track selecting the first hit on each layer of VXD (8 hit f...
Definition: NoKickRTSel.cc:98
Belle2::NoKickRTSel::m_hitXP
std::vector< hitXP > m_hitXP
vector of hit, to convert the track
Definition: NoKickRTSel.h:46
Belle2::NoKickRTSel::m_pt
double m_pt
transverse momentum
Definition: NoKickRTSel.h:64
Belle2::hitXP::setClusterV
void setClusterV(int cluster)
get the relative memeber
Definition: hitXP.h:280
Belle2::hitXP
This class collects some information of a TrueHit, using SVDCLuster and MCParticle information too.
Definition: hitXP.h:42
Belle2::NoKickRTSel::hitXPBuilder
void hitXPBuilder(const RecoTrack &track)
this method build a vector of hitXP from a track.
Definition: NoKickRTSel.cc:23
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2::hitXP::getZ0Entry
double getZ0Entry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:383
Belle2::NoKickRTSel::initNoKickRTSel
void initNoKickRTSel()
Inizialize the class cleaning the member vectors.
Definition: NoKickRTSel.h:86
Belle2::NoKickRTSel::m_PDGIDEff
TH1F * m_PDGIDEff
histogram for efficiency for each PDGID
Definition: NoKickRTSel.h:60
Belle2::hitXP::getOmega0
double getOmega0() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:341
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:215
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticle::getPDG
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:123
Belle2::NoKickRTSel::m_momCut
TH1F * m_momCut
histrogram of cutted tracks
Definition: NoKickRTSel.h:56
Belle2::hitXP::m_momentum0
TVector3 m_momentum0
momentum at IP
Definition: hitXP.h:59
Belle2::NoKickRTSel::m_pmax
double m_pmax
range analyzed with cuts
Definition: NoKickRTSel.h:50
Belle2::hitXP::getD0Entry
double getD0Entry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:359
Belle2::NoKickRTSel::produceHistoNoKick
void produceHistoNoKick()
This method produce the validation histograms (to be used the endrun combined with the filling in tra...
Definition: NoKickRTSel.cc:279
Belle2::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:41
Belle2::NoKickRTSel::m_pMag
double m_pMag
momentum magnitut
Definition: NoKickRTSel.h:63
Belle2::hitXP::setReconstructed
void setReconstructed(bool isReconstructed)
get the relative memeber
Definition: hitXP.h:286
Belle2::hitXP::getZ00
double getZ00() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:389
Belle2::NoKickRTSel::segmentSelector
bool segmentSelector(hitXP hit1, hitXP hit2, std::vector< double > selCut, NoKickCuts::EParameters par, bool is0=false)
This method return true if a couple of hits resects the cuts constraints.
Definition: NoKickRTSel.cc:137
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:38
Belle2::NoKickRTSel::m_8hitTrack
std::vector< hitXP > m_8hitTrack
vector of selected hit
Definition: NoKickRTSel.h:48
Belle2::hitXP::setClusterU
void setClusterU(int cluster)
get the relative memeber
Definition: hitXP.h:274
Belle2::hitXP::getD00
double getD00() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:365
Belle2::NoKickRTSel::m_PDGIDSel
TH1F * m_PDGIDSel
histogram for PDGID of selected track
Definition: NoKickRTSel.h:59
Belle2::VXD::GeoCache::getSensorInfo
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a referecne to the SensorInfo of a given SensorID.
Definition: GeoCache.cc:68
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::hitXP::getTanLambdaEntry
double getTanLambdaEntry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:347
Belle2::hitXPDerivate
This class is the derivate of HitXP, and complete it with a constructor that use all other complex ty...
Definition: hitXPDerivate.h:37
Belle2::StoreArray< SVDCluster >
Belle2::NoKickRTSel::globalCut
bool globalCut(const std::vector< hitXP > &track8)
This method make some global cuts on the tracks (layer 3 and 6 required, d0 and z0 inside beam pipe).
Definition: NoKickRTSel.cc:113
Belle2::NoKickRTSel::m_PDGIDCut
TH1F * m_PDGIDCut
histogram for PDGID of cutted track
Definition: NoKickRTSel.h:58
Belle2::hitXP::getPhi0Entry
double getPhi0Entry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:371
Belle2::hitXP::getTanLambda0
double getTanLambda0() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:353
Belle2::NoKickRTSel::m_nCutHit
TH1F * m_nCutHit
histogram for number of cutted hist per track
Definition: NoKickRTSel.h:61
Belle2::NoKickRTSel::m_numberOfCuts
int m_numberOfCuts
number of catastrophic interaction for each track
Definition: NoKickRTSel.h:51