Belle II Software  release-05-01-25
BtubeCreatorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Sourav Dey, Abi Soffer *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <analysis/modules/BtubeCreator/BtubeCreatorModule.h>
12 
13 // framework aux
14 #include <framework/gearbox/Unit.h>
15 #include <framework/gearbox/Const.h>
16 #include <framework/logging/Logger.h>
17 
18 // dataobjects
19 #include <analysis/dataobjects/Particle.h>
20 #include <analysis/dataobjects/ParticleList.h>
21 #include <analysis/DecayDescriptor/ParticleListName.h>
22 
23 // utilities
24 #include <analysis/utility/PCmsLabTransform.h>
25 #include <analysis/utility/ParticleCopy.h>
26 
27 // Magnetic field
28 #include <framework/geometry/BFieldManager.h>
29 
30 #include <TMath.h>
31 #include <Eigen/Core>
32 
33 using namespace std;
34 using namespace Belle2;
35 
36 //-----------------------------------------------------------------
37 // Register the Module
38 //-----------------------------------------------------------------
39 REG_MODULE(BtubeCreator)
40 
41 //-----------------------------------------------------------------
42 // Implementation
43 //-----------------------------------------------------------------
44 
45 BtubeCreatorModule::BtubeCreatorModule() : Module(),
46  m_Bfield(0)
47 {
48  // Set module properties
49  setDescription("BtubeCreator : This module creates Btube object, an object geometrically similar to a very long ellipsoid along a particular direction and can be used as a constraint in vertexing B particles. Upsilon4S decays to two Bs. the user selects one B with a ^ in decaystring. The order of daughters in reconstructdecay of Upsilon4S and this decaystring should be same. Using selected B as reference, the module calculates the direction in which the other B should fly and constructs a Btube object along that direction. In semi-leptonic analyses, the selected B is usually the fully reconstructed one, while for time dependent studies, the selected B is usually signal B which can be partially reconstructed");
50 
51  // Parameter definitions
52  addParam("listName", m_listName, "name of mother particle list", string(""));
53  addParam("decayString", m_decayString,
54  "decay string of the mother particle, the selected daughter specifies which daughter will be used as reference to create Btube",
55  string(""));
56  addParam("associateBtubeToBselected", m_associateBtubeToBselected,
57  "whether to associate the Btube with the selected B",
58  false);
59  addParam("verbosity", m_verbose, "print statements", true);
60 }
61 
62 void BtubeCreatorModule::initialize()
63 {
64  // magnetic field
65  m_Bfield = BFieldManager::getField(TVector3(0, 0, 0)).Z() / Unit::T;
66 
67  m_BeamSpotCenter = m_beamSpotDB->getIPPosition();
68  m_beamSpotCov.ResizeTo(3, 3);
69  m_beamSpotCov = m_beamSpotDB->getCovVertex();
70  B2INFO("BtubeCreator : magnetic field = " << m_Bfield);
71 
72  tubeArray.registerInDataStore();
73  particles.registerRelationTo(tubeArray);
74 
75  bool valid = m_decaydescriptor.init(m_decayString);
76  if (!valid)
77  B2ERROR("BtubeCreatorModule: invalid decay string" << m_decayString);
78  int nProducts = m_decaydescriptor.getNDaughters();
79  if (nProducts != 2)
80  B2ERROR("BtubeCreatorModule: decay string should contain only two daughters");
81 }
82 
83 void BtubeCreatorModule::event()
84 {
85  StoreObjPtr<ParticleList> plist(m_listName);
86  if (!plist) {
87  B2ERROR("ParticleList " << m_listName << " not found");
88  return;
89  }
90 
91  analysis::RaveSetup::initialize(1, m_Bfield);
92 
93  std::vector<unsigned int> toRemove;
94  unsigned int n = plist->getListSize();
95 
96  for (unsigned i = 0; i < n; i++) {
97  Particle* particle = plist->getParticle(i);
98 
99  std::vector<Particle*> daughtervec = particle->getDaughters();
100  std::vector<const Particle*> selParticles = m_decaydescriptor.getSelectionParticles(particle);
101  if (selParticles.size() > 1) {
102  B2FATAL("Please select only one daughter which will be used as reference to create Btube");
103  }
104 
105  int selectindex = 1;
106  if ((daughtervec.at(0))->getArrayIndex() == (selParticles.at(0))->getArrayIndex()) selectindex = 0 ;
107 
108  Particle* tubecreatorB = const_cast<Particle*>(particle->getDaughter(selectindex));
109  Particle* otherB = const_cast<Particle*>(particle->getDaughter(1 - selectindex));
110 
111  if ((tubecreatorB->getVertexErrorMatrix()(2, 2)) == 0.0) {
112  B2FATAL("Please perform a vertex fit of the selected B before calling this module");
113  }
114 
115  //make a copy of tubecreatorB so as not to modify the original object
116 
117  Particle* tubecreatorBCopy = ParticleCopy::copyParticle(tubecreatorB);
118 
119  if (m_verbose) {
120  B2DEBUG(10, "tubecreator B decay vertex: ");
121  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << tubecreatorBCopy->getVertex()[0] << "," << std::fixed <<
122  std::setprecision(
123  20) << tubecreatorBCopy->getVertex()[1] << "," << std::fixed << std::setprecision(20) << tubecreatorBCopy->getVertex()[2] << "}");
124  }
125 
126  bool ok0 = doVertexFit(tubecreatorBCopy);
127 
128  if (ok0) {
129  particle->setVertex(tubecreatorBCopy->getVertex());
130  particle->setMomentumVertexErrorMatrix(tubecreatorBCopy->getMomentumVertexErrorMatrix());
131 
132  tubecreatorB->writeExtraInfo("prod_vtx_x", tubecreatorBCopy->getVertex()[0]);
133  tubecreatorB->writeExtraInfo("prod_vtx_y", tubecreatorBCopy->getVertex()[1]);
134  tubecreatorB->writeExtraInfo("prod_vtx_z", tubecreatorBCopy->getVertex()[2]);
135  tubecreatorB->writeExtraInfo("prod_vtx_cov00", tubecreatorBCopy->getVertexErrorMatrix()(0, 0));
136  tubecreatorB->writeExtraInfo("prod_vtx_cov01", tubecreatorBCopy->getVertexErrorMatrix()(0, 1));
137  tubecreatorB->writeExtraInfo("prod_vtx_cov02", tubecreatorBCopy->getVertexErrorMatrix()(0, 2));
138  tubecreatorB->writeExtraInfo("prod_vtx_cov10", tubecreatorBCopy->getVertexErrorMatrix()(1, 0));
139  tubecreatorB->writeExtraInfo("prod_vtx_cov11", tubecreatorBCopy->getVertexErrorMatrix()(1, 1));
140  tubecreatorB->writeExtraInfo("prod_vtx_cov12", tubecreatorBCopy->getVertexErrorMatrix()(1, 2));
141  tubecreatorB->writeExtraInfo("prod_vtx_cov20", tubecreatorBCopy->getVertexErrorMatrix()(2, 0));
142  tubecreatorB->writeExtraInfo("prod_vtx_cov21", tubecreatorBCopy->getVertexErrorMatrix()(2, 1));
143  tubecreatorB->writeExtraInfo("prod_vtx_cov22", tubecreatorBCopy->getVertexErrorMatrix()(2, 2));
144 
145  tubecreatorB->writeExtraInfo("Px_after_avf", (tubecreatorBCopy->get4Vector()).Px());
146  tubecreatorB->writeExtraInfo("Py_after_avf", (tubecreatorBCopy->get4Vector()).Py());
147  tubecreatorB->writeExtraInfo("Pz_after_avf", (tubecreatorBCopy->get4Vector()).Pz());
148  tubecreatorB->writeExtraInfo("E_after_avf", (tubecreatorBCopy->get4Vector()).E());
149 
150  Eigen::Matrix<double, 3, 1> tubecreatorBOriginpos(tubecreatorBCopy->getVertex()[0], tubecreatorBCopy->getVertex()[1],
151  tubecreatorBCopy->getVertex()[2]);
152  TLorentzVector v4Final = tubecreatorBCopy->get4Vector();
153  PCmsLabTransform T;
154  TLorentzVector vec = T.rotateLabToCms() * v4Final;
155  TLorentzVector vecNew(-1 * vec.Px(), -1 * vec.Py(), -1 * vec.Pz(), vec.E());
156  TLorentzVector v4FinalNew = T.rotateCmsToLab() * vecNew;
157 
158  if (m_verbose) {
159  B2DEBUG(10, "beamspot center :");
160  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << m_BeamSpotCenter.X() << "," << std::fixed << std::setprecision(
161  20) << m_BeamSpotCenter.Y() << "," << std::fixed << std::setprecision(20) << m_BeamSpotCenter.Z() << "}");
162  B2DEBUG(10, "beamspot cov :");
163 
164  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << m_beamSpotCov(0,
165  0) << "," << std::fixed << std::setprecision(20) << m_beamSpotCov(0,
166  1) << "," << std::fixed << std::setprecision(20) << m_beamSpotCov(0, 2) << "},");
167  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << m_beamSpotCov(1,
168  0) << "," << std::fixed << std::setprecision(20) << m_beamSpotCov(1,
169  1) << "," << std::fixed << std::setprecision(20) << m_beamSpotCov(1, 2) << "},");
170  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << m_beamSpotCov(2,
171  0) << "," << std::fixed << std::setprecision(20) << m_beamSpotCov(2,
172  1) << "," << std::fixed << std::setprecision(20) << m_beamSpotCov(2, 2) << "}");
173  }
174  TMatrixFSym pp = (tubecreatorBCopy->getMomentumErrorMatrix()).GetSub(0, 2, 0, 2, "S");
175  double pe = tubecreatorBCopy->getMomentumErrorMatrix()(2, 2);
176  TMatrixFSym pv = tubecreatorBCopy->getVertexErrorMatrix();
177 
178  // start rotation
179 
180  double theta = v4FinalNew.Theta();
181  double phi = v4FinalNew.Phi();
182 
183  double st = TMath::Sin(theta);
184  double ct = TMath::Cos(theta);
185  double sp = TMath::Sin(phi);
186  double cp = TMath::Cos(phi);
187 
188  TMatrix r2z(3, 3); r2z(2, 2) = 1;
189  r2z(0, 0) = cp; r2z(0, 1) = -1 * sp;
190  r2z(1, 0) = sp; r2z(1, 1) = cp;
191 
192  TMatrix r2y(3, 3); r2y(1, 1) = 1;
193  r2y(0, 0) = ct; r2y(0, 2) = st;
194  r2y(2, 0) = -1 * st; r2y(2, 2) = ct;
195 
196  TMatrix r2(3, 3); r2.Mult(r2z, r2y);
197  TMatrix r2t(3, 3); r2t.Transpose(r2);
198 
199  TMatrix longerror(3, 3); longerror(2, 2) = 1000;
200  TMatrix longerror_temp(3, 3); longerror_temp.Mult(r2, longerror);
201  TMatrix longerrorRotated(3, 3); longerrorRotated.Mult(longerror_temp, r2t);
202 
203  TMatrix pvNew(3, 3);
204  pvNew += pv;
205  pvNew += longerrorRotated;
206 
207  TMatrixFSym errNew(7);
208  errNew.SetSub(0, 0, pp);
209  errNew.SetSub(4, 4, pvNew);
210  errNew(3, 3) = pe;
211 
212  TMatrixFSym tubeMat(3);
213  tubeMat.SetSub(0, 0, pvNew);
214 
215  TMatrixFSym tubeMatCenterError(3);
216  tubeMatCenterError.SetSub(0, 0, pv);
217 
218  if (m_verbose) {
219  B2DEBUG(10, "B origin error matrix : ");
220  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << pv(0, 0) << "," << std::fixed << std::setprecision(20) << pv(0,
221  1) << "," << std::fixed << std::setprecision(20) << pv(0, 2) << "},");
222  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << pv(1, 0) << "," << std::fixed << std::setprecision(20) << pv(1,
223  1) << "," << std::fixed << std::setprecision(20) << pv(1, 2) << "},");
224  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << pv(2, 0) << "," << std::fixed << std::setprecision(20) << pv(2,
225  1) << "," << std::fixed << std::setprecision(20) << pv(2, 2) << "}");
226 
227  B2DEBUG(10, "B tube error matrix : ");
228  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << pvNew(0, 0) << "," << std::fixed << std::setprecision(20) << pvNew(0,
229  1) << "," << std::fixed << std::setprecision(20) << pvNew(0, 2) << "},");
230  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << pvNew(1, 0) << "," << std::fixed << std::setprecision(20) << pvNew(1,
231  1) << "," << std::fixed << std::setprecision(20) << pvNew(1, 2) << "},");
232  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << pvNew(2, 0) << "," << std::fixed << std::setprecision(20) << pvNew(2,
233  1) << "," << std::fixed << std::setprecision(20) << pvNew(2, 2) << "}");
234 
235  B2DEBUG(10, "B origin ");
236  B2DEBUG(10, "{" << std::fixed << std::setprecision(20) << tubecreatorBCopy->getVertex()[0] << "," << std::fixed <<
237  std::setprecision(
238  20) << tubecreatorBCopy->getVertex()[1] << "," << std::fixed << std::setprecision(20) << tubecreatorBCopy->getVertex()[2] << "}");
239  }
240 
241  tubecreatorBCopy->setMomentumVertexErrorMatrix(errNew);
242 
243  Btube* tubeconstraint = tubeArray.appendNew(Btube());
244  if (m_associateBtubeToBselected) {
245  tubecreatorB->addRelationTo(tubeconstraint);
246  } else {
247  otherB->addRelationTo(tubeconstraint);
248  }
249 
250  tubeconstraint->setTubeCenter(tubecreatorBOriginpos);
251  tubeconstraint->setTubeMatrix(tubeMat);
252  tubeconstraint->setTubeCenterErrorMatrix(tubeMatCenterError);
253 
254  if (m_associateBtubeToBselected) {
255  addextrainfos(tubecreatorB, tubecreatorBCopy, pvNew, v4FinalNew);
256  } else {
257  addextrainfos(otherB, tubecreatorBCopy, pvNew, v4FinalNew);
258  }
259  }
260  if (!ok0) toRemove.push_back(particle->getArrayIndex());
261  }
262  plist->removeParticles(toRemove);
263 
264  analysis::RaveSetup::getInstance()->reset();
265 }
266 
267 bool BtubeCreatorModule::doVertexFit(Particle* mother)
268 {
269  analysis::RaveSetup::getInstance()->setBeamSpot(m_BeamSpotCenter, m_beamSpotCov);
270 
271  analysis::RaveVertexFitter rsg;
272  rsg.addTrack(mother);
273  int nvert = rsg.fit("avf");
274 
275  if (nvert == 1) {
276  rsg.updateDaughters();
277  } else {return false;}
278  return true;
279 }
280 
281 void BtubeCreatorModule::addextrainfos(Particle* daughter, Particle* copy, TMatrix mat, TLorentzVector TLV)
282 {
283  daughter->writeExtraInfo("TubePosX", copy->getVertex()[0]);
284  daughter->writeExtraInfo("TubePosY", copy->getVertex()[1]);
285  daughter->writeExtraInfo("TubePosZ", copy->getVertex()[2]);
286 
287  daughter->writeExtraInfo("TubeCov00", mat(0, 0));
288  daughter->writeExtraInfo("TubeCov01", mat(0, 1));
289  daughter->writeExtraInfo("TubeCov02", mat(0, 2));
290  daughter->writeExtraInfo("TubeCov10", mat(1, 0));
291  daughter->writeExtraInfo("TubeCov11", mat(1, 1));
292  daughter->writeExtraInfo("TubeCov12", mat(1, 2));
293  daughter->writeExtraInfo("TubeCov20", mat(2, 0));
294  daughter->writeExtraInfo("TubeCov21", mat(2, 1));
295  daughter->writeExtraInfo("TubeCov22", mat(2, 2));
296 
297  daughter->writeExtraInfo("TubeDirX", TLV.Px());
298  daughter->writeExtraInfo("TubeDirY", TLV.Py());
299  daughter->writeExtraInfo("TubeDirZ", TLV.Pz());
300 }
Belle2::Particle::writeExtraInfo
void writeExtraInfo(const std::string &name, const float value)
Sets the user defined extraInfo.
Definition: Particle.cc:1196
Belle2::Particle::getVertex
TVector3 getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:529
Belle2::Btube::setTubeCenterErrorMatrix
void setTubeCenterErrorMatrix(const TMatrixFSym &tubecentererrormatrix)
Sets Btube Center Error Matrix.
Definition: Btube.h:70
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::RelationsInterface::addRelationTo
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
Definition: RelationsObject.h:144
Belle2::Btube
For each MCParticle with hits in the CDC, this class stores some summarising information on those hit...
Definition: Btube.h:38
Belle2::BtubeCreatorModule
Create a B particle from a Bbar particle.
Definition: BtubeCreatorModule.h:51
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::PCmsLabTransform::rotateCmsToLab
const TLorentzRotation rotateCmsToLab() const
Returns Lorentz transformation from CMS to Lab.
Definition: PCmsLabTransform.h:84
Belle2::Btube::setTubeCenter
void setTubeCenter(const Eigen::Matrix< double, 3, 1 > &tubecenter)
Sets Btube Center.
Definition: Btube.h:57
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::Particle::getMomentumVertexErrorMatrix
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:386
Belle2::analysis::RaveVertexFitter::addTrack
void addTrack(const Particle *const aParticlePtr)
add a track (in the format of a Belle2::Particle) to set of tracks that should be fitted to a vertex
Definition: RaveVertexFitter.cc:87
Belle2::Particle::setMomentumVertexErrorMatrix
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:373
Belle2::analysis::RaveVertexFitter::updateDaughters
void updateDaughters()
update the Daughters particles
Definition: RaveVertexFitter.cc:250
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::Particle::getVertexErrorMatrix
TMatrixFSym getVertexErrorMatrix() const
Returns the 3x3 position error sub-matrix.
Definition: Particle.cc:413
Belle2::PCmsLabTransform
Class to hold Lorentz transformations from/to CMS and boost vector.
Definition: PCmsLabTransform.h:37
Belle2::PCmsLabTransform::rotateLabToCms
const TLorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Definition: PCmsLabTransform.h:74
Belle2::analysis::RaveVertexFitter::fit
int fit(std::string options="default")
do the vertex fit with all tracks previously added with the addTrack or addMother function.
Definition: RaveVertexFitter.cc:121
Belle2::Particle::get4Vector
TLorentzVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:464
Belle2::analysis::RaveVertexFitter
The RaveVertexFitter class is part of the RaveInterface together with RaveSetup.
Definition: RaveVertexFitter.h:46
Belle2::Btube::setTubeMatrix
void setTubeMatrix(const TMatrixFSym &tubematrix)
Sets Btube Matrix.
Definition: Btube.h:65
Belle2::Particle::getMomentumErrorMatrix
TMatrixFSym getMomentumErrorMatrix() const
Returns the 4x4 momentum error matrix.
Definition: Particle.cc:401