Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TruthConversionEval.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TruthConversionEval.C
1 #include "TruthConversionEval.h"
2 #include "Conversion.h"
3 #include "SVReco.h"
4 #include "VtxRegressor.h"
5 //#include "../PhotonConversion/RaveVertexingAux.h"
6 
7 #include <phool/PHCompositeNode.h>
8 #include <phool/getClass.h>
9 
10 #include <calotrigger/CaloTriggerInfo.h>
11 #include <calobase/RawCluster.h>
12 #include <calobase/RawClusterv1.h>
13 
14 #include <g4main/PHG4TruthInfoContainer.h>
15 #include <g4main/PHG4Particlev1.h>
16 #include <g4main/PHG4Particlev2.h>
17 #include <g4main/PHG4VtxPoint.h>
18 
19 /*#include <trackbase_historic/SvtxHitMap.h>
20 #include <trackbase_historic/SvtxHit.h>
21 #include <trackbase_historic/SvtxClusterMap.h>
22 #include <trackbase_historic/SvtxCluster.h>*/
23 #include <trackbase_historic/SvtxVertex.h>
24 #include <trackbase_historic/SvtxVertexMap.h>
25 #include <trackbase_historic/SvtxTrackMap.h>
26 #include <trackbase/TrkrClusterContainer.h>
27 #include <trackbase/TrkrCluster.h>
28 
29 #include <g4eval/SvtxEvalStack.h>
30 #include <g4eval/SvtxTrackEval.h>
31 
32 //#include <GenFit/GFRaveConverters.h>
33 #include <GenFit/FieldManager.h>
34 #include <GenFit/GFRaveVertex.h>
35 #include <GenFit/GFRaveVertexFactory.h>
36 #include <GenFit/MeasuredStateOnPlane.h>
37 #include <GenFit/RKTrackRep.h>
38 #include <GenFit/StateOnPlane.h>
39 #include <GenFit/Track.h>
40 
41 #include <phgenfit/Track.h>
42 
43 #include <fun4all/Fun4AllReturnCodes.h>
44 
45 #include <TFile.h>
46 #include <TTree.h>
47 
48 #include <math.h>
49 #include <utility>
50 #include <list>
51 #include <iostream>
52 
53 TruthConversionEval::TruthConversionEval(const std::string &name, unsigned int runnumber,
54  int particleEmbed, int pythiaEmbed,bool makeTTree=true,string TMVAName="",string TMVAPath="") : SubsysReco("TruthConversionEval"),
55  _kRunNumber(runnumber),_kParticleEmbed(particleEmbed), _kPythiaEmbed(pythiaEmbed), _kMakeTTree(makeTTree)
56 {
57  _foutname = name;
58  _regressor = new VtxRegressor(TMVAName,TMVAPath);
59 }
60 
61 TruthConversionEval::~TruthConversionEval(){
62  cout<<"destruct"<<endl;
63  if (_f) delete _f;
64  if (_vertexer) delete _vertexer;
65  if(_regressor) delete _regressor;
66 }
67 
68 int TruthConversionEval::InitRun(PHCompositeNode *topNode)
69 {
70  _vertexer = new SVReco();
71  _vertexer->InitRun(topNode);
72  if(_kMakeTTree){
73  _runNumber=_kRunNumber;
74  _f = new TFile( _foutname.c_str(), "RECREATE");
75  _observTree = new TTree("observTree","per event observables");
76  _observTree->SetAutoSave(300);
77  _observTree->Branch("nMatched", &_b_nMatched);
78  _observTree->Branch("nUnmatched", &_b_nUnmatched);
79  _observTree->Branch("truth_pT", &_b_truth_pT);
80  _observTree->Branch("reco_pT", &_b_reco_pT);
81  _observTree->Branch("track_pT",&_b_alltrack_pT) ;
82  _observTree->Branch("photon_pT",&_b_allphoton_pT) ;
83 
84  _vtxingTree = new TTree("vtxingTree","data predicting vtx from track pair");
85  _vtxingTree->SetAutoSave(300);
86  _vtxingTree->Branch("vtx_radius", &_b_vtx_radius);
87  _vtxingTree->Branch("tvtx_radius", &_b_tvtx_radius);
88  _vtxingTree->Branch("vtx_phi", &_b_vtx_phi);
89  _vtxingTree->Branch("vtx_eta", &_b_vtx_eta);
90  _vtxingTree->Branch("vtx_x", &_b_vtx_x);
91  _vtxingTree->Branch("vtx_y", &_b_vtx_y);
92  _vtxingTree->Branch("vtx_z", &_b_vtx_z);
93  _vtxingTree->Branch("tvtx_eta", &_b_tvtx_eta);
94  _vtxingTree->Branch("tvtx_x", &_b_tvtx_x);
95  _vtxingTree->Branch("tvtx_y", &_b_tvtx_y);
96  _vtxingTree->Branch("tvtx_z", &_b_tvtx_z);
97  _vtxingTree->Branch("tvtx_phi", &_b_tvtx_phi);
98  _vtxingTree->Branch("vtx_chi2", &_b_vtx_chi2);
99  _vtxingTree->Branch("track1_pt", &_b_track1_pt);
100  _vtxingTree->Branch("track1_eta",& _b_track1_eta);
101  _vtxingTree->Branch("track1_phi",& _b_track1_phi);
102  _vtxingTree->Branch("track2_pt", &_b_track2_pt);
103  _vtxingTree->Branch("track2_eta",& _b_track2_eta);
104  _vtxingTree->Branch("track2_phi",& _b_track2_phi);
105  _vtxingTree->Branch("track_layer",& _b_track_layer);
106 
107  _trackBackTree = new TTree("trackBackTree","track background all single tracks");
108  _trackBackTree->SetAutoSave(300);
109  _trackBackTree->Branch("track_dca", &_bb_track_dca);
110  _trackBackTree->Branch("track_pT", &_bb_track_pT);
111  _trackBackTree->Branch("track_layer", &_bb_track_layer);
112  _trackBackTree->Branch("cluster_prob", &_bb_cluster_prob);
113 
114  _pairBackTree = new TTree("pairBackTree","pair background all possible combinations");
115  _pairBackTree->SetAutoSave(300);
116  _pairBackTree->Branch("track_deta", &_bb_track_deta);
117  _pairBackTree->Branch("track2_pid", &_bb_track2_pid);
118  _pairBackTree->Branch("track1_pid", &_bb_track1_pid);
119  _pairBackTree->Branch("parent_pid", &_bb_parent_pid);
120  _pairBackTree->Branch("track_dphi", &_bb_track_dphi);
121  _pairBackTree->Branch("track_dlayer",&_bb_track_dlayer);
122  _pairBackTree->Branch("approach_dist", &_bb_approach);
123  _pairBackTree->Branch("track_dca", &_bb_track_dca);
124  _pairBackTree->Branch("track_pT", &_bb_track_pT);
125  _pairBackTree->Branch("track_layer", &_bb_track_layer);
126  _pairBackTree->Branch("cluster_prob", &_bb_cluster_prob);
127  //_pairBackTree->Branch("nCluster", &_bb_nCluster);
128  _pairBackTree->Branch("cluster_dphi", &_bb_cluster_dphi);
129  _pairBackTree->Branch("cluster_deta", &_bb_cluster_deta);
130 
131  _vtxBackTree = new TTree("vtxBackTree","events that pass track pair cuts");
132  _vtxBackTree->SetAutoSave(300);
133  _vtxBackTree->Branch("track_deta", &_bb_track_deta);
134  _vtxBackTree->Branch("track1_pid", &_bb_track1_pid);
135  _vtxBackTree->Branch("track2_pid", &_bb_track2_pid);
136  _vtxBackTree->Branch("track_dphi", &_bb_track_dphi);
137  _vtxBackTree->Branch("track_dlayer",&_bb_track_dlayer);
138  _vtxBackTree->Branch("approach_dist", &_bb_approach);
139  _vtxBackTree->Branch("track_dca", &_bb_track_dca);
140  _vtxBackTree->Branch("track_pT", &_bb_track_pT);
141  _vtxBackTree->Branch("track_layer", &_bb_track_layer);
142  _vtxBackTree->Branch("cluster_prob", &_bb_cluster_prob);
143  _vtxBackTree->Branch("vtx_radius",&_bb_vtx_radius);
144  _vtxBackTree->Branch("photon_m",&_bb_photon_m);
145  _vtxBackTree->Branch("photon_pT",&_bb_photon_pT);
146  //_vtxBackTree->Branch("nCluster", &_bb_nCluster);
147  _vtxBackTree->Branch("cluster_dphi", &_bb_cluster_dphi);
148  _vtxBackTree->Branch("cluster_deta", &_bb_cluster_deta);
149 
150  _signalCutTree = new TTree("cutTreeSignal","signal data for making track pair cuts");
151  _signalCutTree->SetAutoSave(100);
152  _signalCutTree->Branch("track_deta", &_b_track_deta);
153  _signalCutTree->Branch("track_dca", &_b_track_dca);
154  _signalCutTree->Branch("track_dphi", &_b_track_dphi);
155  _signalCutTree->Branch("track_dlayer",&_b_track_dlayer);
156  _signalCutTree->Branch("track_layer", &_b_track_layer);
157  _signalCutTree->Branch("track_pT", &_b_track_pT);
158  //_signalCutTree->Branch("ttrack_pT", &_b_ttrack_pT);
159  _signalCutTree->Branch("approach_dist", &_b_approach);
160  _signalCutTree->Branch("vtx_radius", &_b_vtx_radius);
161  _signalCutTree->Branch("vtx_phi", &_b_vtx_phi);
162  _signalCutTree->Branch("tvtx_phi", &_b_tvtx_phi);
163  _signalCutTree->Branch("tvtx_radius", &_b_tvtx_radius);
164  _signalCutTree->Branch("vtx_chi2", &_b_vtx_chi2);
165  //_signalCutTree->Branch("vtxTrackRZ_dist", &_b_vtxTrackRZ_dist);
166  //_signalCutTree->Branch("vtxTrackRPhi_dist", &_b_vtxTrackRPhi_dist);
167  _signalCutTree->Branch("photon_m", &_b_photon_m);
168  _signalCutTree->Branch("tphoton_m", &_b_tphoton_m);
169  _signalCutTree->Branch("photon_pT", &_b_photon_pT);
170  _signalCutTree->Branch("tphoton_pT", &_b_tphoton_pT);
171  _signalCutTree->Branch("cluster_prob", &_b_cluster_prob);
172  //_signalCutTree->Branch("nCluster", &_b_nCluster);
173  _signalCutTree->Branch("cluster_dphi", &_b_cluster_dphi);
174  _signalCutTree->Branch("cluster_deta", &_b_cluster_deta);
175  }
176  return 0;
177 }
178 
179 bool TruthConversionEval::doNodePointers(PHCompositeNode* topNode){
180  bool goodPointers=true;
181  _mainClusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
182  _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
183  _clusterMap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
184  _allTracks = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
185  // _hitMap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");
186  //if(!_mainClusterContainer||!_truthinfo||!_clusterMap||!_hitMap){
187  if(!_mainClusterContainer||!_truthinfo||!_clusterMap||!_allTracks){
188  cerr<<Name()<<": critical error-bad nodes \n";
189  if(!_mainClusterContainer){
190  cerr<<"\t RawClusterContainer is bad";
191  }
192  if(!_truthinfo){
193  cerr<<"\t TruthInfoContainer is bad";
194  }
195  if(!_clusterMap){
196  cerr<<"\t TrkrClusterMap is bad";
197  }
198  if(!_allTracks){
199  cerr<<"\t SvtxTrackMap is bad";
200  }
201  cerr<<endl;
202  goodPointers=false;
203  }
204  return goodPointers;
205 }
206 
207 SvtxVertex* TruthConversionEval::get_primary_vertex(PHCompositeNode *topNode)const{
208  SvtxVertexMap *vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
209  return vertexMap->get(0);
210 }
211 
212 int TruthConversionEval::process_event(PHCompositeNode *topNode)
213 {
214  if(!doNodePointers(topNode)) return Fun4AllReturnCodes::ABORTEVENT;
215  cout<<"init vertexer event"<<endl;
216  _vertexer->InitEvent(topNode);
217  _conversionClusters.Reset(); //clear the list of conversion clusters
218  PHG4TruthInfoContainer::Range range = _truthinfo->GetParticleRange(); //look at all truth particles
219  SvtxEvalStack *stack = new SvtxEvalStack(topNode); //truth tracking info
220  SvtxTrackEval* trackeval = stack->get_track_eval();
221  if (!trackeval)
222  {
223  cerr<<"NULL track eval in " <<Name()<<" fatal error"<<endl;
224  return 1;
225  }
226  //make a map of the conversions
227  std::map<int,Conversion> mapConversions;
228  //h is for hadronic e is for EM
229  std::vector<SvtxTrack*> backgroundTracks;
230  std::vector<std::pair<SvtxTrack*,SvtxTrack*>> tightbackgroundTrackPairs; //used to find the pair for unmatched conversion tracks
231  std::vector<PHG4Particle*> truthPhotons;
232  std::list<int> signalTracks;
233  //reset obervation variables
234  _b_nMatched=0;
235  _b_nUnmatched=0;
236  _b_truth_pT.clear();
237  _b_reco_pT.clear();
238  _b_alltrack_pT.clear();
239  _b_allphoton_pT.clear();
240  cout<<"init truth loop"<<endl;
241 
242  for ( PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter ) {
243  PHG4Particle* g4particle = iter->second;
244  if(g4particle->get_pid()==22&&g4particle->get_track_id()==g4particle->get_primary_id()){
245  _b_allphoton_pT.push_back(sqrt(g4particle->get_px()*g4particle->get_px()+g4particle->get_py()*g4particle->get_py()));
246  truthPhotons.push_back(g4particle);
247  }
248  PHG4Particle* parent =_truthinfo->GetParticle(g4particle->get_parent_id());
249  PHG4VtxPoint* vtx=_truthinfo->GetVtx(g4particle->get_vtx_id()); //get the vertex
250  if(!vtx){
251  cerr<<"null vtx primaryid="<<g4particle->get_primary_id()<<'\n';
252  g4particle->identify();
253  cout<<endl;
254  continue;
255  }
256  float radius=sqrt(vtx->get_x()*vtx->get_x()+vtx->get_y()*vtx->get_y());
257  //if outside the tracker skip this
258  if(radius>s_kTPCRADIUS) continue;
259  int embedID;
260  if (parent)//if the particle is not primary
261  {
262  embedID=get_embed(parent,_truthinfo);
263  float truthpT = sqrt(g4particle->get_px()*g4particle->get_px()+g4particle->get_py()*g4particle->get_py());
264  if(parent->get_pid()==22&&TMath::Abs(g4particle->get_pid())==11&&truthpT>2.5){ //conversion check
265  if (Verbosity()==10)
266  {
267  std::cout<<"Conversion with radius [cm]:"<<radius<<'\n';
268  }
269  //initialize the conversion object -don't use constructor b/c setters have error handling
270  (mapConversions[vtx->get_id()]).setElectron(g4particle);
271  (mapConversions[vtx->get_id()]).setVtx(vtx);
272  (mapConversions[vtx->get_id()]).setPrimaryPhoton(parent,_truthinfo);
273  (mapConversions[vtx->get_id()]).setEmbed(embedID);
274  PHG4Particle* grand =_truthinfo->GetParticle(parent->get_parent_id()); //grandparent
275  if (grand) (mapConversions[vtx->get_id()]).setSourceId(grand->get_pid());//record pid of the photon's source
276  else (mapConversions[vtx->get_id()]).setSourceId(0);//or it is from the G4 generator
277  _b_truth_pT.push_back(truthpT);
278  //build a list of the ids
279  SvtxTrack* recoTrack = trackeval->best_track_from(g4particle);
280  if(recoTrack){
281  signalTracks.push_back(recoTrack->get_id());
282  _b_reco_pT.push_back(recoTrack->get_pt());
283  cout<<"matched truth track"<<endl;
284  _b_nMatched++;
285  }
286  else{
287  _b_nUnmatched++;
288  cerr<<"WARNING no matching track for conversion"<<endl;
289  }
290  }
291  }// not primary
292  }//truth particle loop
293  signalTracks.sort();
294  cout<<"intit track loop"<<endl;
295  //build the current backgroundSample
296  for ( SvtxTrackMap::Iter iter = _allTracks->begin(); iter != _allTracks->end(); ++iter) {
297  _b_alltrack_pT.push_back(iter->second->get_pt());
298  auto inCheck = std::find(signalTracks.begin(),signalTracks.end(),iter->first);
299  //if the track is not in the list of signal tracks
300  if (inCheck==signalTracks.end())
301  {
302  TLorentzVector track_tlv = tracktoTLV(iter->second);
303  bool addTrack=true;
304  //make sure the track is not too close to a photon
305  for (std::vector<PHG4Particle*>::iterator iPhoton = truthPhotons.begin(); iPhoton != truthPhotons.end()&&addTrack; ++iPhoton)
306  {
307  TLorentzVector photon_tlv = particletoTLV(*iPhoton);
308  if (photon_tlv.DeltaR(track_tlv)<.2)
309  {
310  addTrack=false;
311  }
312  }
313  if(addTrack){
314  backgroundTracks.push_back(iter->second);
315  //see if it is a good tight background cannidate
316  for ( SvtxTrackMap::Iter jter = std::next(iter,1); jter != _allTracks->end()&&iter->second->get_pt()>_kTightPtMin; ++jter){
317  if(jter->second->get_pt()>_kTightPtMin&&TMath::Abs(jter->second->get_eta()-iter->second->get_eta())<_kTightDetaMax&&
318  iter->second->get_charge()==jter->second->get_charge()*-1){
319  tightbackgroundTrackPairs.push_back(std::pair<SvtxTrack*,SvtxTrack*>(iter->second,jter->second));
320  }
321  }
322 
323  }
324  }
325  }
326  //pass the map to this helper method which fills the fields for the signal ttree
327  numUnique(&mapConversions,trackeval,_mainClusterContainer,&tightbackgroundTrackPairs);
328 
329  if (_kMakeTTree)
330  {
331  cout<<"intit background process"<<endl;
332  if(_b_truth_pT.size()!=0) _observTree->Fill();
333  processTrackBackground(&backgroundTracks,trackeval);
334  cout<<"finish back"<<endl;
335  }
336  delete stack;
337  cout<<"finish process"<<endl;
338  return 0;
339 }
340 
341 void TruthConversionEval::numUnique(std::map<int,Conversion> *mymap=NULL,SvtxTrackEval* trackeval=NULL,RawClusterContainer *mainClusterContainer=NULL,
342  std::vector<std::pair<SvtxTrack*,SvtxTrack*>>* backgroundTracks=NULL){
343  cout<<"start conversion analysis loop"<<endl;
344  for (std::map<int,Conversion>::iterator i = mymap->begin(); i != mymap->end(); ++i) {
345  TLorentzVector tlv_photon,tlv_electron,tlv_positron; //make tlv for each particle
346  PHG4Particle *photon = i->second.getPrimaryPhoton(); //set the photon
347 
348  if(photon)tlv_photon.SetPxPyPzE(photon->get_px(),photon->get_py(),photon->get_pz(),photon->get_e()); //intialize
349  else cerr<<"No truth photon for conversion"<<endl;
350  PHG4Particle *e1=i->second.getElectron(); //set the first child
351  if(e1){
352  tlv_electron.SetPxPyPzE(e1->get_px(),e1->get_py(),e1->get_pz(),e1->get_e());
353  }
354  else cerr<<"No truth electron for conversion"<<endl;
355  PHG4Particle *e2=i->second.getPositron();
356  if(e2){ //this will be false for conversions with 1 truth track
357  tlv_positron.SetPxPyPzE(e2->get_px(),e2->get_py(),e2->get_pz(),e2->get_e()); //init the tlv
358  if (TMath::Abs(tlv_electron.Eta())<_kRAPIDITYACCEPT&&TMath::Abs(tlv_positron.Eta())<_kRAPIDITYACCEPT)
359  {
360  unsigned int nRecoTracks = i->second.setRecoTracks(trackeval); //find the reco tracks for this conversion
361  switch(nRecoTracks)
362  {
363  case 2: //there are 2 reco tracks
364  {
365  recordConversion(&(i->second),&tlv_photon,&tlv_electron,&tlv_positron);
366  break;
367  }
368  case 1: //there's one reco track try to find the other
369  {
370  PHG4Particle *truthe;
371  //get the truth particle that is missing a reco track
372  if(!i->second.getRecoTrack(e1->get_track_id()))truthe = e1;
373  else truthe=e2;
374  if(!truthe) cerr<<"critical error"<<endl;
375  //look through the background for the missing pair
376  for(auto pair : *backgroundTracks){
377  if(!pair.first||!pair.second) cerr<<"critical error2"<<endl;
378  if ((pair.first->get_charge()>0&&truthe->get_pid()<0)||(pair.first->get_charge()<0&&truthe->get_pid()>0))
379  {
380  TVector3 truth_tlv(truthe->get_px(),truthe->get_py(),truthe->get_pz());
381  TVector3 reco_tlv(pair.first->get_px(),pair.first->get_py(),pair.first->get_pz());
382  if (reco_tlv.DeltaR(truth_tlv)<.1)
383  {
384  i->second.setRecoTrack(truthe->get_track_id(),pair.first);
385  }
386  }
387  else if ((pair.second->get_charge()>0&&truthe->get_pid()<0)||(pair.second->get_charge()<0&&truthe->get_pid()>0))
388  {
389  TVector3 truth_tlv(truthe->get_px(),truthe->get_py(),truthe->get_pz());
390  TVector3 reco_tlv(pair.second->get_px(),pair.second->get_py(),pair.second->get_pz());
391  if (reco_tlv.DeltaR(truth_tlv)<.1)
392  {
393  i->second.setRecoTrack(truthe->get_track_id(),pair.second);
394  }
395  }
396  }
397  recordConversion(&(i->second),&tlv_photon,&tlv_electron,&tlv_positron);
398  break;
399  }
400  case 0: //no reco tracks
401  //TODO maybe try to find the reco tracks? for now just record the bad info
402  recordConversion(&(i->second),&tlv_photon,&tlv_electron,&tlv_positron);
403  break;
404  default:
405  if (Verbosity()>1)
406  {
407  cerr<<Name()<<" error setting reco tracks"<<endl;
408  }
409  break;
410  }//switch
411  }//rapidity cut
412  }// has 2 truth tracks
413  cout<<"map loop"<<endl;
414  }//map loop
415  cout<<"done num"<<endl;
416 }
417 
418 //only call if _kMakeTTree is true
419 void TruthConversionEval::processTrackBackground(std::vector<SvtxTrack*> *v_tracks,SvtxTrackEval* trackeval){
420  Conversion pairMath;
421  float lastpT=-1.;
422  unsigned nNullTrack=0;
423  cout<<"The total possible background track count is "<<v_tracks->size()<<'\n';
424  for (std::vector<SvtxTrack*>::iterator iter = v_tracks->begin(); iter != v_tracks->end(); ++iter) {
425  cout<<"looping"<<endl;
426  SvtxTrack* iTrack = *iter;
427  //get the SvtxTrack it must not be NULL
428  if(!iTrack){
429  nNullTrack++;
430  continue;
431  }
432 
433  if(TMath::Abs(iTrack->get_eta())>1.1||iTrack->get_pt()==lastpT)continue; //TODO this skips duplicate tracks but why are they there in the first place?
434  lastpT=iTrack->get_pt();
435  cout<<"\t pT="<<lastpT<<'\n';
436 
437  auto temp_key_it=iTrack->begin_cluster_keys();//key iterator to first cluster
438  if(temp_key_it!=iTrack->end_cluster_keys()){//if the track has clusters
439  TrkrCluster* temp_cluster = _clusterMap->findCluster(*temp_key_it);//get the cluster
440  if(temp_cluster) _bb_track_layer = TrkrDefs::getLayer(temp_cluster->getClusKey());//if there is a cluster record its layer
441  else _bb_track_layer=-999;
442  }
443  //record track info
444  _bb_track_dca = iTrack->get_dca();
445  _bb_track_pT = iTrack->get_pt();
446  //record the EMCal cluster info
447  auto cluster1 = _mainClusterContainer->getCluster(iTrack->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)));
448  if(cluster1) _bb_cluster_prob= cluster1->get_prob();
449  else _bb_cluster_prob=-999;
450  //pair with other tracks
451  for(std::vector<SvtxTrack*>::iterator jter =std::next(iter,1);jter!=v_tracks->end(); ++jter){//posible bias by filling the track level variables with iTrack instead of min(iTrack,jTrack)
452  SvtxTrack* jTrack = *jter;
453  if(!jTrack||TMath::Abs(jTrack->get_eta())>1.1)continue;
454  //record pair geometry
455  cout<<"calculations"<<endl;
456  _bb_track_deta = pairMath.trackDEta(iTrack,jTrack);
457  _bb_track_dphi = pairMath.trackDPhi(iTrack,jTrack);
458  _bb_track_dlayer = pairMath.trackDLayer(_clusterMap,iTrack,jTrack);
459  _bb_approach = pairMath.approachDistance(iTrack,jTrack);
460  //record second EMCal cluster
461  auto cluster2 = _mainClusterContainer->getCluster(iTrack->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)));
462  //if the EMCal clusters can be found record their pair geometry
463  if (cluster2&&cluster1)
464  {
465  //if the two clusters are unique calculate the values
466  if(cluster2->get_id()!=cluster1->get_id())
467  {
468  _bb_nCluster = 2;
469  _bb_cluster_dphi=fabs(cluster1->get_phi()-cluster2->get_phi());
470  TVector3 etaCalc(cluster1->get_x(),cluster1->get_y(),cluster1->get_z());
471  float eta1 = etaCalc.PseudoRapidity();
472  etaCalc.SetXYZ(cluster2->get_x(),cluster2->get_y(),cluster2->get_z());
473  _bb_cluster_deta=fabs(eta1-etaCalc.PseudoRapidity());
474  }
475  else{ //if they are not unique then they are 0
476  _bb_nCluster = 1;
477  _bb_cluster_dphi=0;
478  _bb_cluster_deta=0;
479  }
480  }
481  else{ //clusters were not found
482  _bb_nCluster=0;
483  _bb_cluster_deta=-999;
484  _bb_cluster_dphi=-999;
485  }
486  /*_bb_track1_pid = (*iTruthTrack)->get_pid();
487  _bb_track2_pid = (*jTruthTrack)->get_pid();
488  PHG4Particle* parent = _truthinfo->GetParticle((*iTruthTrack)->get_parent_id());
489  if(parent) _bb_parent_pid = parent->get_pid();
490  else _bb_parent_pid=0;*/
491  //if the track pairs pass the pair cuts make them part of the vertex background
492  if (_bb_track_layer>=0&&_bb_track_pT>_kTightPtMin&&_bb_track_deta<_kTightDetaMax&&TMath::Abs(_bb_track_dlayer)<9)
493  {
494  //iTrack->identify();
495  //jTrack->identify();
496  cout<<"here vertex"<<endl;
497  //make the vertex
498  genfit::GFRaveVertex* recoVert = _vertexer->findSecondaryVertex(iTrack,jTrack);
499  cout<<"made vertex"<<endl;
500  //if the vertex is made record its info
501  if (recoVert)
502  {
503  recoVert=pairMath.correctSecondaryVertex(_regressor,recoVert,iTrack,jTrack);
504  cout<<"corrected vertex"<<endl;
505  TVector3 recoVertPos = recoVert->getPos();
506  cout<<"deleting"<<endl;
507  delete recoVert;
508  cout<<"deleted"<<endl;
509  //fill the tree values
510  _bb_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
511  _bb_vtx_chi2 = recoVert->getChi2();
512  _bb_vtxTrackRZ_dist = pairMath.vtxTrackRZ(recoVertPos,iTrack,jTrack);
513  _bb_vtxTrackRPhi_dist = pairMath.vtxTrackRPhi(recoVertPos,iTrack,jTrack);
514  //then make the photon
515  TLorentzVector* recoPhoton= pairMath.getRecoPhoton(iTrack,jTrack);
516  if(recoPhoton){
517  _bb_photon_m = recoPhoton->Dot(*recoPhoton);
518  _bb_photon_pT = recoPhoton->Pt();
519  delete recoPhoton;
520  }
521  else{
522  _bb_photon_m=-999;
523  _bb_photon_pT=-999;
524  }
525  }
526  else{ // no vtx could be found for thre track pair
527  _bb_vtx_radius = -999;
528  _bb_vtx_chi2 = -999;
529  _bb_vtxTrackRZ_dist =-999;
530  _bb_vtxTrackRPhi_dist =-999;
531  }
532  cout<<"fill vtx"<<endl;
533  _vtxBackTree->Fill();
534  }//pair cuts
535  _pairBackTree->Fill();
536  cout<<"filled pair"<<endl;
537  }//jTrack loop
538  cout<<"filling track"<<endl;
539  _trackBackTree->Fill();
540  cout<<"filled track"<<endl;
541  }//iTrack loop
542  cout<<"Null track count ="<<nNullTrack<<endl;
543 }
544 
545 void TruthConversionEval::recordConversion(Conversion *conversion,TLorentzVector *tlv_photon,TLorentzVector *tlv_electron, TLorentzVector *tlv_positron){
546  cout<<"recording"<<endl;
547  if(tlv_photon&&tlv_electron&&tlv_positron&&conversion&&conversion->recoCount()==2){
548  _b_track_deta = conversion->trackDEta();
549  _b_track_dphi = conversion->trackDPhi();
550  _b_track_dlayer = conversion->trackDLayer(_clusterMap);
551  _b_track_layer = conversion->firstLayer(_clusterMap);
552  _b_approach = conversion->approachDistance();
553  _b_track_dca = conversion->minDca();
554  //record pT info
555  _b_track_pT = conversion->minTrackpT();
556  if(tlv_electron->Pt()>tlv_positron->Pt()) _b_ttrack_pT = tlv_positron->Pt();
557  else _b_ttrack_pT = tlv_electron->Pt();
558  //record initial photon info
559  TLorentzVector* recoPhoton = conversion->getRecoPhoton();
560  if (recoPhoton)
561  {
562  _b_photon_m=recoPhoton->Dot(*recoPhoton);
563  TLorentzVector truth_added_tlv = *tlv_electron+*tlv_positron;
564  _b_tphoton_m= truth_added_tlv.Dot(truth_added_tlv);
565  _b_photon_pT=recoPhoton->Pt();
566  conversion->PrintPhotonRecoInfo(tlv_photon,tlv_electron,tlv_positron,_b_photon_m);
567  }
568  else{//photon was not reconstructed
569  _b_photon_m =-999;
570  _b_photon_pT=-999;
571  conversion->PrintPhotonRecoInfo(tlv_photon,tlv_electron,tlv_positron,_b_photon_m);
572  }
573  _b_tphoton_pT=tlv_photon->Pt();
574  cout<<"second"<<endl;
575  //truth vertex info
576  _b_tvtx_radius = sqrt(conversion->getVtx()->get_x()*conversion->getVtx()->get_x()+conversion->getVtx()->get_y()*conversion->getVtx()->get_y());
577  TVector3 tVertPos(conversion->getVtx()->get_x(),conversion->getVtx()->get_y(),conversion->getVtx()->get_z());
578  _b_tvtx_phi = tVertPos.Phi();
579  _b_tvtx_eta = tVertPos.Eta();
580  _b_tvtx_z = tVertPos.Z();
581  _b_tvtx_x = tVertPos.X();
582  _b_tvtx_y = tVertPos.Y();
583  //TODO check Conversion operations for ownership transfer->memleak due to lack of delete
584  cout<<"vertexing"<<endl;
585  genfit::GFRaveVertex* originalVert, *recoVert;
586  originalVert=recoVert = conversion->getSecondaryVertex(_vertexer);
587  recoVert = conversion->correctSecondaryVertex(_regressor);
588  cout<<"finding gf_tracks"<<endl;
589  //std::pair<PHGenFit::Track*,PHGenFit::Track*> ph_gf_tracks = conversion->getPHGFTracks(_vertexer);
590  if (recoVert)
591  {
592  //cout<<"finding refit_gf_tracks"<<endl;
593  //std::pair<PHGenFit::Track*,PHGenFit::Track*> refit_phgf_tracks=conversion->refitTracks(_vertexer);
594  //TODO check repetive refitting and revterexing this is issue #23
595  /*if (ph_gf_tracks.first&&refit_phgf_tracks.first)
596  {
597  cout<<"Good Track refit with original:\n";ph_gf_tracks.first->get_mom().Print();cout<<"\n\t and refit:\n";
598  refit_phgf_tracks.first->get_mom().Print();
599  }
600  if (ph_gf_tracks.second&&refit_phgf_tracks.second)
601  {
602  cout<<"Good Track refit with original:\n";
603  ph_gf_tracks.second->get_mom().Print();
604  cout<<"\n\t and refit:\n";
605  refit_phgf_tracks.second->get_mom().Print();
606  }*/
607  recoPhoton = conversion->getRefitRecoPhoton();
608  if(recoPhoton) _b_photon_pT=recoPhoton->Pt();
609  //fill the vtxing tree
610  TVector3 recoVertPos = originalVert->getPos();
611  _b_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
612  _b_vtx_phi = recoVertPos.Phi();
613  _b_vtx_eta = recoVertPos.Eta();
614  _b_vtx_z = recoVertPos.Z();
615  _b_vtx_x = recoVertPos.X();
616  _b_vtx_y = recoVertPos.Y();
617  _b_vtx_chi2 = recoVert->getChi2();
618  //track info
619  pair<float,float> pTstemp = conversion->getTrackpTs();
620  _b_track1_pt = pTstemp.first;
621  _b_track2_pt = pTstemp.second;
622  pair<float,float> etasTemp = conversion->getTrackEtas();
623  _b_track1_eta = etasTemp.first;
624  _b_track2_eta = etasTemp.second;
625  pair<float,float> phisTemp = conversion->getTrackPhis();
626  _b_track1_phi = phisTemp.first;
627  _b_track2_phi = phisTemp.second;
628  _vtxingTree->Fill();
629 
630  //populate the refit vertex values
631  recoVertPos = recoVert->getPos();
632  _b_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
633  _b_vtx_phi = recoVertPos.Phi();
634  _b_vtx_eta = recoVertPos.Eta();
635  _b_vtx_z = recoVertPos.Z();
636  _b_vtx_x = recoVertPos.X();
637  _b_vtx_y = recoVertPos.Y();
638  _b_vtxTrackRZ_dist = conversion->vtxTrackRZ(recoVertPos);
639  _b_vtxTrackRPhi_dist = conversion->vtxTrackRPhi(recoVertPos);
640  cout<<"done full vtx values"<<endl;
641  }
642 
643  else{//vtx not reconstructed
644  cout<<"no vtx "<<endl;
645  _b_vtx_radius =-999.;
646  _b_vtx_phi = -999.;
647  _b_vtx_eta = -999.;
648  _b_vtx_z = -999.;
649  _b_vtx_x = -999.;
650  _b_vtx_y = -999.;
651  _b_vtx_chi2 = -999.;
652  _b_vtxTrackRZ_dist = -999.;
653  _b_vtxTrackRPhi_dist = -999.;
654  //track info
655  pair<float,float> pTstemp = conversion->getTrackpTs();
656  _b_track1_pt = pTstemp.first;
657  _b_track2_pt = pTstemp.second;
658  pair<float,float> etasTemp = conversion->getTrackEtas();
659  _b_track1_eta = etasTemp.first;
660  _b_track2_eta = etasTemp.second;
661  pair<float,float> phisTemp = conversion->getTrackPhis();
662  _b_track1_phi = phisTemp.first;
663  _b_track2_phi = phisTemp.second;
664  cout<<"done no vtx values"<<endl;
665  _vtxingTree->Fill();
666  }
667  //reset the values
668  _b_cluster_prob=-999;
669  _b_cluster_deta=-999;
670  _b_cluster_dphi=-999;
671  _b_nCluster=-999;
672  pair<int,int> clusterIds = conversion->get_cluster_ids();
673  RawCluster *clustemp;
674  if(_mainClusterContainer->getCluster(clusterIds.first)){//if there is matching cluster
675  clustemp = dynamic_cast<RawCluster*>(_mainClusterContainer->getCluster(clusterIds.first)->CloneMe());
676  //this is for cluster subtraction which will not be implented soon
677  // _conversionClusters.AddCluster(clustemp); //add the calo cluster to the container
678  if (_kMakeTTree)
679  {
680  _b_cluster_prob=clustemp->get_prob();
681  RawCluster *clus2 = _mainClusterContainer->getCluster(clusterIds.second);
682  if (clus2)
683  {
684  _b_cluster_dphi=fabs(clustemp->get_phi()-clus2->get_phi());
685  TVector3 etaCalc(clustemp->get_x(),clustemp->get_y(),clustemp->get_z());
686  //TODO check cluster_prob distribution for signal
687  if (clus2->get_prob()>_b_cluster_prob)
688  {
689  _b_cluster_prob=clus2->get_prob();
690  }
691  //calculate deta
692  float eta1 = etaCalc.PseudoRapidity();
693  etaCalc.SetXYZ(clus2->get_x(),clus2->get_y(),clus2->get_z());
694  _b_cluster_deta=fabs(eta1-etaCalc.PseudoRapidity());
695  if (clusterIds.first!=clusterIds.second) //if there are two district clusters
696  {
697  _b_nCluster=2;
698  }
699  else{
700  _b_nCluster=1;
701  }
702  }
703  }
704  }
705  }
706 
707  else{//not a "complete" conversion some of the data is missing
708  cout<<"incomplete conversion"<<endl;
709  _b_track_deta = -999;
710  _b_track_dphi = -999;
711  _b_track_dlayer = -999;
712  _b_track_layer = -999;
713  _b_approach = -999;
714  _b_track_dca = -999;
715  //record pT info
716  _b_track_pT = -999;
717  if(!tlv_electron||!tlv_positron) _b_ttrack_pT = -999;
718  else if(tlv_electron->Pt()>tlv_positron->Pt()) _b_ttrack_pT = tlv_positron->Pt();
719  else _b_ttrack_pT = tlv_electron->Pt();
720  //record initial photon info
721  TLorentzVector* recoPhoton = conversion->getRecoPhoton();
722  if (recoPhoton)
723  {
724  _b_photon_pT=recoPhoton->Pt();
725  _b_photon_m=recoPhoton->Dot(*recoPhoton);
726  if (tlv_positron&&tlv_electron)
727  {
728  TLorentzVector truth_added_tlv = *tlv_electron+*tlv_positron;
729  _b_tphoton_m= truth_added_tlv.Dot(truth_added_tlv);
730  }
731  else{
732  _b_tphoton_m=-999;
733  }
734  }
735  else{//photon was not reconstructed
736  _b_photon_m =-999;
737  _b_photon_pT=-999;
738  }
739  if(tlv_photon)_b_tphoton_pT=tlv_photon->Pt();
740  else _b_tphoton_pT=-999;
741  cout<<"here"<<endl;
742  //truth vertex info
743  _b_tvtx_radius = sqrt(conversion->getVtx()->get_x()*conversion->getVtx()->get_x()+conversion->getVtx()->get_y()*conversion->getVtx()->get_y());
744  TVector3 tVertPos(conversion->getVtx()->get_x(),conversion->getVtx()->get_y(),conversion->getVtx()->get_z());
745  cout<<"still here"<<endl;
746  _b_tvtx_phi = tVertPos.Phi();
747  _b_tvtx_eta = tVertPos.Eta();
748  _b_tvtx_z = tVertPos.Z();
749  _b_tvtx_x = tVertPos.X();
750  _b_tvtx_y = tVertPos.Y();
751  //TODO check Conversion operations for ownership transfer->memleak due to lack of delete
752  //vtx not reconstructed because conversion is incomplete
753  _b_vtx_radius = -999.;
754  _b_vtx_phi = -999.;
755  _b_vtx_eta = -999.;
756  _b_vtx_z = -999.;
757  _b_vtx_x = -999.;
758  _b_vtx_y = -999.;
759  _b_vtx_chi2 = -999.;
760  _b_vtxTrackRZ_dist = -999.;
761  _b_vtxTrackRPhi_dist = -999.;
762  //track info not possible to reconstruct both tracks
763  _b_track1_pt = -999;
764  _b_track2_pt = -999;
765  _b_track1_eta = -999;
766  _b_track2_eta = -999;
767  _b_track1_phi = -999;
768  _b_track2_phi = -999;
769 
770  //reset the values
771  _b_cluster_prob=-1;
772  _b_cluster_deta=-1;
773  _b_cluster_dphi=-1;
774  _b_nCluster=-1;
775  pair<int,int> clusterIds = conversion->get_cluster_ids();
776  RawCluster *clustemp;
777  if(_mainClusterContainer->getCluster(clusterIds.first)){//if there is matching cluster
778  clustemp = dynamic_cast<RawCluster*>(_mainClusterContainer->getCluster(clusterIds.first)->CloneMe());
779  //this is for cluster subtraction which will not be implented soon
780 
781  // _conversionClusters.AddCluster(clustemp); //add the calo cluster to the container
782  if (_kMakeTTree)
783  {
784  _b_cluster_prob=clustemp->get_prob();
785  RawCluster *clus2 = _mainClusterContainer->getCluster(clusterIds.second);
786  if (clus2)
787  {
788  _b_cluster_dphi=fabs(clustemp->get_phi()-clus2->get_phi());
789  TVector3 etaCalc(clustemp->get_x(),clustemp->get_y(),clustemp->get_z());
790  //TODO check cluster_prob distribution for signal
791  if (clus2->get_prob()>_b_cluster_prob)
792  {
793  _b_cluster_prob=clus2->get_prob();
794  }
795  //calculate deta
796  float eta1 = etaCalc.PseudoRapidity();
797  etaCalc.SetXYZ(clus2->get_x(),clus2->get_y(),clus2->get_z());
798  _b_cluster_deta=fabs(eta1-etaCalc.PseudoRapidity());
799  if (clusterIds.first!=clusterIds.second) //if there are two district clusters
800  {
801  _b_nCluster=2;
802  }
803  else{
804  _b_nCluster=1;
805  }
806  }
807  }
808  }
809  }
810  cout<<"Filling"<<endl;
811  _signalCutTree->Fill();
812 }
813 
814 const RawClusterContainer* TruthConversionEval::getClusters()const {return &_conversionClusters;}
815 
816 int TruthConversionEval::get_embed(PHG4Particle* particle, PHG4TruthInfoContainer* truthinfo)const{
817  return truthinfo->isEmbeded(particle->get_track_id());
818 }
819 
820 float TruthConversionEval::vtoR(PHG4VtxPoint* vtx)const{
821  return (float) sqrt(vtx->get_x()*vtx->get_x()+vtx->get_y()*vtx->get_y());
822 }
823 
824 int TruthConversionEval::End(PHCompositeNode *topNode)
825 {
826  cout<<"ending"<<endl;
827  if(_kMakeTTree){
828  cout<<"closing"<<endl;
829  //_signalCutTree->Write();
830  _f->Write();
831  _f->Close();
832  }
833  return 0;
834 }