Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VtxTest.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file VtxTest.C
1 #include "VtxTest.h"
2 #include "Conversion.h"
3 #include "SVReco.h"
4 //#include "../PhotonConversion/RaveVertexingAux.h"
5 
6 #include <phool/PHCompositeNode.h>
7 #include <phool/getClass.h>
8 
9 #include <calotrigger/CaloTriggerInfo.h>
10 #include <calobase/RawCluster.h>
11 #include <calobase/RawClusterv1.h>
12 
13 #include <fun4all/Fun4AllReturnCodes.h>
14 
15 #include <g4main/PHG4TruthInfoContainer.h>
16 #include <g4main/PHG4Particle.h>
17 #include <g4main/PHG4Particlev1.h>
18 #include <g4main/PHG4Particlev2.h>
19 #include <g4main/PHG4VtxPoint.h>
20 
21 /*#include <trackbase_historic/SvtxHitMap.h>
22 #include <trackbase_historic/SvtxHit.h>
23 #include <trackbase_historic/SvtxClusterMap.h>
24 #include <trackbase_historic/SvtxCluster.h>*/
25 #include <trackbase/TrkrClusterContainer.h>
26 
27 #include <g4eval/SvtxEvalStack.h>
28 #include <g4eval/SvtxTrackEval.h>
29 
30 //#include <GenFit/GFRaveConverters.h>
31 #include <GenFit/FieldManager.h>
32 #include <GenFit/GFRaveVertex.h>
33 #include <GenFit/GFRaveVertexFactory.h>
34 #include <GenFit/MeasuredStateOnPlane.h>
35 #include <GenFit/RKTrackRep.h>
36 #include <GenFit/StateOnPlane.h>
37 #include <GenFit/Track.h>
38 
39 #include <TFile.h>
40 #include <TTree.h>
41 #include <TLorentzVector.h>
42 
43 #include <utility>
44 #include <iostream>
45 #include <math.h>
46 
47 VtxTest::VtxTest(const std::string &name, unsigned int runnumber,
48  int particleEmbed, int pythiaEmbed,bool makeTTree=true) : SubsysReco("VtxTest"),
49  _kRunNumber(runnumber),_kParticleEmbed(particleEmbed), _kPythiaEmbed(pythiaEmbed), _kMakeTTree(makeTTree)
50 {
51  _foutname = name;
52 }
53 
54 VtxTest::~VtxTest(){
55  if (_f) delete _f;
56  if (_vertexer) delete _vertexer;
57 }
58 
59 int VtxTest::InitRun(PHCompositeNode *topNode)
60 {
61  _vertexer = new SVReco();
62  _vertexer->InitRun(topNode);
63  if(_kMakeTTree){
64  _b_event=0;
65  _runNumber=_kRunNumber;
66  _f = new TFile( _foutname.c_str(), "RECREATE");
67  _tree = new TTree("ttree","conversion data");
68  _tree->SetAutoSave(300);
69  _tree->Branch("runNumber",&_runNumber);
70  _tree->Branch("event",&_b_event);
71  _tree->Branch("nVtx", &_b_nVtx);
72  _tree->Branch("nTpair", &_b_Tpair);
73  _tree->Branch("nRpair", &_b_Rpair);
74  _tree->Branch("rVtx", _b_rVtx,"rVtx[nVtx]/D");
75  _tree->Branch("pythia", _b_pythia,"pythia[nVtx]/B");
76  _tree->Branch("electron_pt", _b_electron_pt,"electron_pt[nVtx]/F");
77  _tree->Branch("positron_reco_pt", _b_positron_reco_pt,"positron_reco_pt[nTpair]/F");
78  _tree->Branch("electron_reco_pt", _b_electron_reco_pt,"electron_reco_pt[nTpair]/F");
79  _tree->Branch("positron_pt", _b_positron_pt,"positron_pt[nVtx]/F");
80  _tree->Branch("photon_pt", _b_parent_pt ,"photon_pt[nVtx]/F");
81  _tree->Branch("photon_eta", _b_parent_eta ,"photon_eta[nVtx]/F");
82  _tree->Branch("photon_phi", _b_parent_phi ,"photon_phi[nVtx]/F");
83  _tree->Branch("e_deta", _b_e_deta ,"e_deta[nTpair]/F");
84  _tree->Branch("e_dphi", _b_e_dphi ,"e_dphi[nTpair]/F");
85  _tree->Branch("fLayer",_b_fLayer,"fLayer[nTpair]/I");
86  _tree->Branch("photon_source_id", _b_grandparent_id ,"photon_source_id[nVtx]/I");
87  _tree->Branch("nCluster",_b_nCluster,"nCluster[nRpair]/I");
88  _tree->Branch("clus_dphi",_b_cluster_dphi,"clus_dphi[nRpair]/F");
89  _tree->Branch("clus_deta",_b_cluster_deta,"clus_deta[nRpair]/F");
90  _tree->Branch("Scluster_prob", &_b_Scluster_prob,"Scluster_prob[nRpair]/F");
91  _tree->Branch("Mcluster_prob", &_b_Mcluster_prob,"Mcluster_prob[nRpair]/F");
92 
93  _signalCutTree = new TTree("cutTreeSignal","signal data for making track pair cuts");
94  _signalCutTree->SetAutoSave(300);
95  _signalCutTree->Branch("track_deta", &_b_track_deta);
96  _signalCutTree->Branch("track_dphi", &_b_track_dphi);
97  _signalCutTree->Branch("track_dlayer",&_b_track_dlayer);
98  _signalCutTree->Branch("track_layer", &_b_track_layer);
99  _signalCutTree->Branch("track_pT", &_b_track_pT);
100  _signalCutTree->Branch("ttrack_pT", &_b_ttrack_pT);
101  _signalCutTree->Branch("approach_dist", &_b_approach);
102  _signalCutTree->Branch("vtx_radius", &_b_vtx_radius);
103  _signalCutTree->Branch("tvtx_radius", &_b_tvtx_radius);
104  _signalCutTree->Branch("vtx_phi", &_b_vtx_phi);
105  _signalCutTree->Branch("vtx_eta", &_b_vtx_eta);
106  _signalCutTree->Branch("vtx_x", &_b_vtx_x);
107  _signalCutTree->Branch("vtx_y", &_b_vtx_y);
108  _signalCutTree->Branch("vtx_z", &_b_vtx_z);
109  _signalCutTree->Branch("tvtx_eta", &_b_tvtx_eta);
110  _signalCutTree->Branch("tvtx_x", &_b_tvtx_x);
111  _signalCutTree->Branch("tvtx_y", &_b_tvtx_y);
112  _signalCutTree->Branch("tvtx_z", &_b_tvtx_z);
113  _signalCutTree->Branch("tvtx_phi", &_b_tvtx_phi);
114  _signalCutTree->Branch("vtx_chi2", &_b_vtx_chi2);
115  _signalCutTree->Branch("vtxTrack_dist", &_b_vtxTrack_dist);
116  _signalCutTree->Branch("photon_m", &_b_photon_m);
117  _signalCutTree->Branch("photon_pT", &_b_photon_pT);
118  _signalCutTree->Branch("cluster_prob", &_b_cluster_prob);
119 
120  _h_backgroundCutTree = new TTree("cutTreeBackh","background data for making track pair cuts");
121  _h_backgroundCutTree->SetAutoSave(300);
122  _h_backgroundCutTree->Branch("track_deta", &_bb_track_deta);
123  _h_backgroundCutTree->Branch("track_dphi", &_bb_track_dphi);
124  _h_backgroundCutTree->Branch("track_dlayer", &_bb_track_dlayer);
125  _h_backgroundCutTree->Branch("track_layer", &_bb_track_layer);
126  _h_backgroundCutTree->Branch("track_pT", &_bb_track_pT);
127  _h_backgroundCutTree->Branch("vtx_radius", &_bb_vtx_radius);
128  _h_backgroundCutTree->Branch("vtx_chi2", &_bb_vtx_chi2);
129  _h_backgroundCutTree->Branch("approach_dist", &_bb_approach);
130  _h_backgroundCutTree->Branch("vtxTrack_dist", &_bb_vtxTrack_dist);
131  _h_backgroundCutTree->Branch("photon_m", &_bb_photon_m);
132  _h_backgroundCutTree->Branch("photon_pT", &_bb_photon_pT);
133  _h_backgroundCutTree->Branch("cluster_prob", &_bb_cluster_prob);
134  _h_backgroundCutTree->Branch("pid", &_bb_pid);
135 
136  _e_backgroundCutTree = new TTree("cutTreeBacke","background data for making track pair cuts");
137  _e_backgroundCutTree->SetAutoSave(300);
138  _e_backgroundCutTree->Branch("track_deta", &_bb_track_deta);
139  _e_backgroundCutTree->Branch("track_dphi", &_bb_track_dphi);
140  _e_backgroundCutTree->Branch("track_dlayer", &_bb_track_dlayer);
141  _e_backgroundCutTree->Branch("track_layer", &_bb_track_layer);
142  _e_backgroundCutTree->Branch("track_pT", &_bb_track_pT);
143  _e_backgroundCutTree->Branch("vtx_radius", &_bb_vtx_radius);
144  _e_backgroundCutTree->Branch("vtx_chi2", &_bb_vtx_chi2);
145  _e_backgroundCutTree->Branch("approach_dist", &_bb_approach);
146  _e_backgroundCutTree->Branch("vtxTrack_dist", &_bb_vtxTrack_dist);
147  _e_backgroundCutTree->Branch("photon_m", &_bb_photon_m);
148  _e_backgroundCutTree->Branch("photon_pT", &_bb_photon_pT);
149  _e_backgroundCutTree->Branch("cluster_prob", &_bb_cluster_prob);
150  }
151  return 0;
152 }
153 
154 void VtxTest::doNodePointers(PHCompositeNode* topNode){
155  _mainClusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
156  _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
157  _clusterMap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
158  // _hitMap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");
159  //if(!_mainClusterContainer||!_truthinfo||!_clusterMap||!_hitMap){
160  if(!_mainClusterContainer||!_truthinfo||!_clusterMap){
161  cerr<<Name()<<": critical error-bad nodes \n";
162  if(!_mainClusterContainer){
163  cerr<<"\t RawClusterContainer is bad";
164  }
165  if(!_truthinfo){
166  cerr<<"\t TruthInfoContainer is bad";
167  }
168  if(!_clusterMap){
169  cerr<<"\t TrkrClusterMap is bad";
170  }
171  /*if(!_hitMap){
172  cerr<<"\t SvtxHitMap is bad";
173  }*/
174  cerr<<endl;
175  }
176 }
177 
178 int VtxTest::process_event(PHCompositeNode *topNode)
179 {
180  doNodePointers(topNode);
181  _vertexer->InitEvent(topNode);
182  PHG4TruthInfoContainer::Range range = _truthinfo->GetParticleRange(); //look at all truth particles
183  SvtxEvalStack *stack = new SvtxEvalStack(topNode); //truth tracking info
184  SvtxTrackEval* trackeval = stack->get_track_eval();
185  SvtxTrack *svtxtrack1=NULL;
186  SvtxTrack *svtxtrack2=NULL;
187  PHG4VtxPoint *truth_vtx=NULL;
188  for ( PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter ) {
189  PHG4Particle* g4particle = iter->second;
190  if(g4particle&&g4particle->get_pid()==211&&!_truthinfo->GetParticle(g4particle->get_parent_id())){
191  cout<<"found pion\n";
192  TLorentzVector truth;
193  truth.SetPxPyPzE(g4particle->get_px(),g4particle->get_py(),g4particle->get_pz(),g4particle->get_e());
194  if(!svtxtrack1)svtxtrack1=trackeval->best_track_from(g4particle);
195  else svtxtrack2=trackeval->best_track_from(g4particle);
196  if(truth_vtx==NULL) truth_vtx=_truthinfo->GetVtx(g4particle->get_vtx_id()); //get the vertex
197  else{
198  PHG4VtxPoint *vtemp = truth_vtx;
199  truth_vtx=_truthinfo->GetVtx(g4particle->get_vtx_id());
200  if(!(*truth_vtx==*vtemp)) cout<<"vtx does not agree"<<endl;
201  }
202  }
203  }
204  if(!truth_vtx){
205  cout<<"did not find truth_vtx"<<endl;
206  return Fun4AllReturnCodes::ABORTEVENT;
207  }
208  if(!svtxtrack1||!svtxtrack2){
209  cout<<"did not find tracks"<<endl;
210  return Fun4AllReturnCodes::ABORTEVENT;
211  }
212  genfit::GFRaveVertex* recoVert=_vertexer->findSecondaryVertex(svtxtrack1,svtxtrack2);
213  cout<<"here"<<endl;
214  if (recoVert)
215  {
216  TVector3 recoVertPos = recoVert->getPos();
217  _b_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
218  _b_tvtx_radius = sqrt(truth_vtx->get_x()*truth_vtx->get_x()+truth_vtx->get_y()*truth_vtx->get_y());
219  _b_vtx_phi = recoVertPos.Phi();
220  _b_vtx_eta = recoVertPos.Eta();
221  _b_vtx_z = recoVertPos.Z();
222  _b_vtx_x = recoVertPos.X();
223  _b_vtx_y = recoVertPos.Y();
224  TVector3 tVertPos(truth_vtx->get_x(),truth_vtx->get_y(),truth_vtx->get_z());
225  _b_tvtx_phi = tVertPos.Phi();
226  _b_tvtx_eta = tVertPos.Eta();
227  _b_tvtx_z = tVertPos.Z();
228  _b_tvtx_x = tVertPos.X();
229  _b_tvtx_y = tVertPos.Y();
230  _b_vtx_chi2 = recoVert->getChi2();
231  }
232  else{
233  _b_vtx_radius = 0;
234  _b_vtx_phi = 0;
235  _b_vtx_eta = 0;
236  _b_vtx_z = 0;
237  _b_vtx_x = 0;
238  _b_vtx_y = 0;
239  _b_tvtx_radius = sqrt(truth_vtx->get_x()*truth_vtx->get_x()+truth_vtx->get_y()*truth_vtx->get_y());
240  TVector3 tVertPos(truth_vtx->get_x(),truth_vtx->get_y(),truth_vtx->get_z());
241  _b_tvtx_phi = tVertPos.Phi();
242  _b_tvtx_eta = tVertPos.Eta();
243  _b_tvtx_z = tVertPos.Z();
244  _b_tvtx_x = tVertPos.X();
245  _b_tvtx_y = tVertPos.Y();
246  }
247  delete stack;
248  _signalCutTree->Fill();
249  return 0;
250 }
251 
252 std::queue<std::pair<int,int>> VtxTest::numUnique(std::map<int,Conversion> *mymap=NULL,SvtxTrackEval* trackeval=NULL,RawClusterContainer *mainClusterContainer=NULL){
253  _b_nVtx = 0;
254  _b_Tpair=0;
255  _b_Rpair=0;
256  std::queue<std::pair<int,int>> missingChildren;//this is deprecated
257  for (std::map<int,Conversion>::iterator i = mymap->begin(); i != mymap->end(); ++i) {
258  PHG4VtxPoint *vtx =i->second.getVtx(); //get the vtx
259  PHG4Particle *temp = i->second.getPhoton(); //set the photon
260  TLorentzVector tlv_photon,tlv_electron,tlv_positron; //make tlv for each particle
261  tlv_photon.SetPxPyPzE(temp->get_px(),temp->get_py(),temp->get_pz(),temp->get_e()); //intialize
262  temp=i->second.getElectron(); //set the first child
263  tlv_electron.SetPxPyPzE(temp->get_px(),temp->get_py(),temp->get_pz(),temp->get_e());
264  if(_kMakeTTree){//fill tree values
265  cout<<"numUnique filling tree"<<endl;
266  _b_rVtx[_b_nVtx] = sqrt(vtx->get_x()*vtx->get_x()+vtx->get_y()*vtx->get_y()); //find the radius
267  _b_parent_pt[_b_nVtx] =tlv_photon.Pt();
268  _b_parent_phi[_b_nVtx]=tlv_photon.Phi();
269  _b_parent_eta[_b_nVtx]=tlv_photon.Eta();
270  _b_grandparent_id[_b_nVtx]=i->second.getSourceId();
271  _b_e_deta[_b_nVtx]=-1.;
272  _b_e_dphi[_b_nVtx]=-1.;
273  _b_positron_pt[_b_nVtx]=-1.;
274  _b_electron_pt[_b_nVtx]=tlv_electron.Pt();
275  }
276  temp=i->second.getPositron();
277  if(temp){ //this will be false for conversions with 1 truth track
278  tlv_positron.SetPxPyPzE(temp->get_px(),temp->get_py(),temp->get_pz(),temp->get_e()); //init the tlv
279  if(_kMakeTTree) _b_positron_pt[_b_nVtx]=tlv_positron.Pt(); //fill tree
280  if (TMath::Abs(tlv_electron.Eta())<_kRAPIDITYACCEPT&&TMath::Abs(tlv_positron.Eta())<_kRAPIDITYACCEPT)
281  {
282  unsigned int nRecoTracks = i->second.setRecoTracks(trackeval); //find the reco tracks for this conversion
283  if(_kMakeTTree){
284  _b_e_deta[_b_Tpair]=TMath::Abs(tlv_electron.Eta()-tlv_positron.Eta());
285  _b_e_dphi[_b_Tpair]=TMath::Abs(tlv_electron.Phi()-tlv_positron.Phi());
286  _b_fLayer[_b_Tpair]=_b_track_layer = i->second.firstLayer(_clusterMap);
287  _b_fLayer[_b_Tpair]=-1;
288  _b_Tpair++;
289  }//tree
290  switch(nRecoTracks)
291  {
292  case 2: //there are 2 reco tracks
293  {
294  if(_kMakeTTree){
295  _b_track_deta = i->second.trackDEta();
296  _b_track_dphi = i->second.trackDPhi();
297  _b_track_dlayer = i->second.trackDLayer(_clusterMap);
298  _b_track_dlayer=-1;
299  _b_track_pT = i->second.minTrackpT();
300  if(tlv_electron.Pt()>tlv_positron.Pt()) _b_ttrack_pT = tlv_positron.Pt();
301  else _b_ttrack_pT = tlv_electron.Pt();
302  _b_approach = i->second.approachDistance();
303  pair<SvtxTrack*, SvtxTrack*> reco_tracks=i->second.getRecoTracks();
304  genfit::GFRaveVertex* recoVert = _vertexer->findSecondaryVertex(reco_tracks.first,reco_tracks.second);
305  if (recoVert)
306  {
307  TVector3 recoVertPos = recoVert->getPos();
308  _b_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
309  _b_tvtx_radius = _b_rVtx[_b_nVtx];
310  _b_vtx_phi = recoVertPos.Phi();
311  _b_vtx_eta = recoVertPos.Eta();
312  _b_vtx_z = recoVertPos.Z();
313  _b_vtx_x = recoVertPos.X();
314  _b_vtx_y = recoVertPos.Y();
315  TVector3 tVertPos(vtx->get_x(),vtx->get_y(),vtx->get_z());
316  _b_tvtx_phi = tVertPos.Phi();
317  _b_tvtx_eta = tVertPos.Eta();
318  _b_tvtx_z = tVertPos.Z();
319  _b_tvtx_x = tVertPos.X();
320  _b_tvtx_y = tVertPos.Y();
321  _b_vtx_chi2 = recoVert->getChi2();
322  _b_vtxTrack_dist = i->second.dist(&recoVertPos,_clusterMap);
323  }
324 
325  TLorentzVector* recoPhoton = i->second.setRecoPhoton();
326  if (recoPhoton)
327  {
328  _b_photon_m=recoPhoton->Dot(*recoPhoton);
329  _b_photon_pT=recoPhoton->Pt();
330  }
331  else{ //photon was not reconstructed
332  _b_photon_m =-1;
333  _b_photon_pT=-1;
334  }
335  //reset the values
336  _b_cluster_prob=-1;
337  _b_cluster_deta[_b_Rpair]=-1.;
338  _b_cluster_dphi[_b_Rpair]=-1.;
339  _b_nCluster[_b_Rpair]=0;
340  _b_Scluster_prob[_b_Rpair]=-1;
341  _b_Mcluster_prob[_b_Rpair]=-1;
342 
343  }
344  pair<int,int> clusterIds = i->second.get_cluster_ids();
345  RawCluster *clustemp;
346  if(mainClusterContainer->getCluster(clusterIds.first)){//if there is matching cluster
347  clustemp = dynamic_cast<RawCluster*>(mainClusterContainer->getCluster(clusterIds.first)->Clone());
348  _conversionClusters.AddCluster(clustemp); //add the calo cluster to the container
349  if (_kMakeTTree)
350  {
351  _b_cluster_prob=clustemp->get_prob();
352  RawCluster *clus2 = mainClusterContainer->getCluster(clusterIds.second);
353  if (clus2)
354  {
355  _b_cluster_dphi[_b_Rpair]=fabs(clustemp->get_phi()-clus2->get_phi());
356  TVector3 etaCalc(clustemp->get_x(),clustemp->get_y(),clustemp->get_z());
357  if (clus2->get_prob()>_b_cluster_prob)
358  {
359  _b_cluster_prob=clus2->get_prob();
360  }
361  //calculate deta
362  float eta1 = etaCalc.PseudoRapidity();
363  etaCalc.SetXYZ(clus2->get_x(),clus2->get_y(),clus2->get_z());
364  _b_cluster_deta[_b_Rpair]=fabs(eta1-etaCalc.PseudoRapidity());
365  /* this may be a logical bug. What are S and M? If this is false the identical cluster prob is recorded twice.*/
366  if (clusterIds.first!=clusterIds.second) //if there are two district clusters
367  {
368  _b_nCluster[_b_Rpair]=2;
369  _b_Scluster_prob[_b_Rpair]=clustemp->get_prob();
370  }
371  else{
372  _b_nCluster[_b_Rpair]=1;
373  _b_Mcluster_prob[_b_Rpair]=clustemp->get_prob();
374  }
375  }
376  }
377  }
378  _signalCutTree->Fill();
379  _b_Rpair++;
380  break;
381  }
382  case 1: //there's one reco track
383  {
384  int clustidtemp=i->second.get_cluster_id(); //get the cluster id of the current conversion
385  if(mainClusterContainer->getCluster(clustidtemp)){//if thre is matching cluster
386  RawCluster *clustemp = dynamic_cast<RawCluster*>(mainClusterContainer->getCluster(clustidtemp)->Clone());
387  _conversionClusters.AddCluster(clustemp); //add the calo cluster to the container
388  }
389  //cout<<"Matched 1 reco with layer="<<i->second.firstLayer(_clusterMap)<<"pTs:"<<tlv_electron.Pt()<<"-"<<tlv_positron.Pt()<<'\n';
390  break;
391  }
392  case 0: //no reco tracks
393  break;
394  default:
395  if (Verbosity()>1)
396  {
397  cerr<<Name()<<" error setting reco tracks"<<endl;
398  }
399  break;
400  }//switch
401  }//rapidity cut
402  if (_kMakeTTree)
403  {
404  _b_pythia[_b_nVtx]=i->second.getEmbed()==_kPythiaEmbed;
405  _b_nVtx++;
406  }
407  }// has 2 truth tracks
408  }//map loop
409  return missingChildren;
410 }
411 
412 //only call if _kMakeTTree is true
413 void VtxTest::processBackground(std::map<int,Conversion> *mymap,SvtxTrackEval* trackeval,TTree* tree){
414  for (std::map<int,Conversion>::iterator i = mymap->begin(); i != mymap->end(); ++i) {
415  int nReco=i->second.setRecoTracks(trackeval);
416  if (nReco==2)
417  {
418  _bb_track_deta = i->second.trackDEta();
419  _bb_track_dphi = i->second.trackDPhi();
420  _bb_track_dlayer = i->second.trackDLayer(_clusterMap);
421  _bb_track_layer = i->second.firstLayer(_clusterMap);
422  _bb_track_layer=-1;
423  _bb_track_dlayer=-1;
424  _bb_track_pT = i->second.minTrackpT();
425  _bb_approach = i->second.approachDistance();
426  _bb_pid = i->second.getElectron()->get_pid();
427  pair<SvtxTrack*, SvtxTrack*> reco_tracks=i->second.getRecoTracks();
428  genfit::GFRaveVertex* recoVert = _vertexer->findSecondaryVertex(reco_tracks.first,reco_tracks.second);
429  if (recoVert)
430  {
431  TVector3 recoVertPos = recoVert->getPos();
432  _bb_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
433  _bb_vtx_chi2 = recoVert->getChi2();
434  _bb_vtxTrack_dist = i->second.dist(&recoVertPos,_clusterMap);
435  }
436  TLorentzVector* recoPhoton = i->second.setRecoPhoton();
437  if (recoPhoton)
438  {
439  _bb_photon_m=recoPhoton->Dot(*recoPhoton);
440  _bb_photon_pT=recoPhoton->Pt();
441  }
442  else{
443  _bb_photon_m =0;
444  _bb_photon_pT=0;
445  }
446  _bb_cluster_prob= _mainClusterContainer->getCluster(i->second.get_cluster_id())->get_prob();
447  tree->Fill();
448  }
449  }
450 }
451 
452 void VtxTest::findChildren(std::queue<std::pair<int,int>> missingChildren,PHG4TruthInfoContainer* _truthinfo){
453  if (Verbosity()>=6)
454  {
455  while(!missingChildren.empty()){
456  for (PHG4TruthInfoContainer::ConstIterator iter = _truthinfo->GetParticleRange().first; iter != _truthinfo->GetParticleRange().second; ++iter)
457  {
458  if(Verbosity()>1&& iter->second->get_parent_id()==missingChildren.front().first&&iter->second->get_track_id()!=missingChildren.front().second){
459  cout<<"Found Child:\n";
460  iter->second->identify();
461  cout<<"With mother:\n";
462 
463  }
464  }
465  missingChildren.pop();
466  }
467  }
468 }
469 
470 const RawClusterContainer* VtxTest::getClusters()const {return &_conversionClusters;}
471 
472 int VtxTest::get_embed(PHG4Particle* particle, PHG4TruthInfoContainer* truthinfo)const{
473  return truthinfo->isEmbeded(particle->get_track_id());
474 }
475 
476 float VtxTest::vtoR(PHG4VtxPoint* vtx)const{
477  return (float) sqrt(vtx->get_x()*vtx->get_x()+vtx->get_y()*vtx->get_y());
478 }
479 
480 
481 int VtxTest::End(PHCompositeNode *topNode)
482 {
483  if(_kMakeTTree){
484  _f->Write();
485  _f->Close();
486  }
487  return 0;
488 }