PhotonConversion/truthconversion/VtxTest.C

0001 #include "VtxTest.h"
0002 #include "Conversion.h"
0003 #include "SVReco.h"
0004 //#include "../PhotonConversion/RaveVertexingAux.h"
0005
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/getClass.h>
0008
0009 #include <calotrigger/CaloTriggerInfo.h>
0010 #include <calobase/RawCluster.h>
0011 #include <calobase/RawClusterv1.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015 #include <g4main/PHG4TruthInfoContainer.h>
0016 #include <g4main/PHG4Particle.h>
0017 #include <g4main/PHG4Particlev1.h>
0018 #include <g4main/PHG4Particlev2.h>
0019 #include <g4main/PHG4VtxPoint.h>
0020
0021 /*#include <trackbase_historic/SvtxHitMap.h>
0022 #include <trackbase_historic/SvtxHit.h>
0023 #include <trackbase_historic/SvtxClusterMap.h>
0024 #include <trackbase_historic/SvtxCluster.h>*/
0025 #include <trackbase/TrkrClusterContainer.h>
0026
0027 #include <g4eval/SvtxEvalStack.h>
0028 #include <g4eval/SvtxTrackEval.h>
0029
0030 //#include <GenFit/GFRaveConverters.h>
0031 #include <GenFit/FieldManager.h>
0032 #include <GenFit/GFRaveVertex.h>
0033 #include <GenFit/GFRaveVertexFactory.h>
0034 #include <GenFit/MeasuredStateOnPlane.h>
0035 #include <GenFit/RKTrackRep.h>
0036 #include <GenFit/StateOnPlane.h>
0037 #include <GenFit/Track.h>
0038
0039 #include <TFile.h>
0040 #include <TTree.h>
0041 #include <TLorentzVector.h>
0042
0043 #include <utility>
0044 #include <iostream>
0045 #include <math.h>
0046
0047 VtxTest::VtxTest(const std::string &name, unsigned int runnumber,
0048     int particleEmbed,  int pythiaEmbed,bool makeTTree=true) : SubsysReco("VtxTest"),
0049   _kRunNumber(runnumber),_kParticleEmbed(particleEmbed), _kPythiaEmbed(pythiaEmbed), _kMakeTTree(makeTTree)
0050 {
0051   _foutname = name;
0052 }
0053
0054 VtxTest::~VtxTest(){
0055   if (_f) delete _f;
0056   if (_vertexer) delete _vertexer;
0057 }
0058
0059 int VtxTest::InitRun(PHCompositeNode *topNode)
0060 {
0061   _vertexer = new SVReco();
0062   _vertexer->InitRun(topNode);
0063   if(_kMakeTTree){
0064     _b_event=0;
0065     _runNumber=_kRunNumber;
0066     _f = new TFile( _foutname.c_str(), "RECREATE");
0067     _tree = new TTree("ttree","conversion data");
0068     _tree->SetAutoSave(300);
0069     _tree->Branch("runNumber",&_runNumber);
0070     _tree->Branch("event",&_b_event);
0071     _tree->Branch("nVtx", &_b_nVtx);
0072     _tree->Branch("nTpair", &_b_Tpair);
0073     _tree->Branch("nRpair", &_b_Rpair);
0074     _tree->Branch("rVtx", _b_rVtx,"rVtx[nVtx]/D");
0075     _tree->Branch("pythia", _b_pythia,"pythia[nVtx]/B");
0076     _tree->Branch("electron_pt", _b_electron_pt,"electron_pt[nVtx]/F");
0077     _tree->Branch("positron_reco_pt", _b_positron_reco_pt,"positron_reco_pt[nTpair]/F");
0078     _tree->Branch("electron_reco_pt", _b_electron_reco_pt,"electron_reco_pt[nTpair]/F");
0079     _tree->Branch("positron_pt", _b_positron_pt,"positron_pt[nVtx]/F");
0080     _tree->Branch("photon_pt",   _b_parent_pt    ,"photon_pt[nVtx]/F");
0081     _tree->Branch("photon_eta",  _b_parent_eta  ,"photon_eta[nVtx]/F");
0082     _tree->Branch("photon_phi",  _b_parent_phi  ,"photon_phi[nVtx]/F");
0083     _tree->Branch("e_deta",  _b_e_deta  ,"e_deta[nTpair]/F");
0084     _tree->Branch("e_dphi",  _b_e_dphi  ,"e_dphi[nTpair]/F");
0085     _tree->Branch("fLayer",_b_fLayer,"fLayer[nTpair]/I");
0086     _tree->Branch("photon_source_id",  _b_grandparent_id  ,"photon_source_id[nVtx]/I");
0087     _tree->Branch("nCluster",_b_nCluster,"nCluster[nRpair]/I");
0088     _tree->Branch("clus_dphi",_b_cluster_dphi,"clus_dphi[nRpair]/F");
0089     _tree->Branch("clus_deta",_b_cluster_deta,"clus_deta[nRpair]/F");
0090     _tree->Branch("Scluster_prob", &_b_Scluster_prob,"Scluster_prob[nRpair]/F");
0091     _tree->Branch("Mcluster_prob", &_b_Mcluster_prob,"Mcluster_prob[nRpair]/F");
0092
0093     _signalCutTree = new TTree("cutTreeSignal","signal data for making track pair cuts");
0094     _signalCutTree->SetAutoSave(300);
0095     _signalCutTree->Branch("track_deta", &_b_track_deta);
0096     _signalCutTree->Branch("track_dphi", &_b_track_dphi);
0097     _signalCutTree->Branch("track_dlayer",&_b_track_dlayer);
0098     _signalCutTree->Branch("track_layer", &_b_track_layer);
0099     _signalCutTree->Branch("track_pT", &_b_track_pT);
0100     _signalCutTree->Branch("ttrack_pT", &_b_ttrack_pT);
0101     _signalCutTree->Branch("approach_dist", &_b_approach);
0102     _signalCutTree->Branch("vtx_radius", &_b_vtx_radius);
0103     _signalCutTree->Branch("tvtx_radius", &_b_tvtx_radius);
0104     _signalCutTree->Branch("vtx_phi", &_b_vtx_phi);
0105     _signalCutTree->Branch("vtx_eta", &_b_vtx_eta);
0106     _signalCutTree->Branch("vtx_x", &_b_vtx_x);
0107     _signalCutTree->Branch("vtx_y", &_b_vtx_y);
0108     _signalCutTree->Branch("vtx_z", &_b_vtx_z);
0109     _signalCutTree->Branch("tvtx_eta", &_b_tvtx_eta);
0110     _signalCutTree->Branch("tvtx_x", &_b_tvtx_x);
0111     _signalCutTree->Branch("tvtx_y", &_b_tvtx_y);
0112     _signalCutTree->Branch("tvtx_z", &_b_tvtx_z);
0113     _signalCutTree->Branch("tvtx_phi", &_b_tvtx_phi);
0114     _signalCutTree->Branch("vtx_chi2", &_b_vtx_chi2);
0115     _signalCutTree->Branch("vtxTrack_dist", &_b_vtxTrack_dist);
0116     _signalCutTree->Branch("photon_m", &_b_photon_m);
0117     _signalCutTree->Branch("photon_pT", &_b_photon_pT);
0118     _signalCutTree->Branch("cluster_prob", &_b_cluster_prob);
0119
0120     _h_backgroundCutTree = new TTree("cutTreeBackh","background data for making track pair cuts");
0121     _h_backgroundCutTree->SetAutoSave(300);
0122     _h_backgroundCutTree->Branch("track_deta", &_bb_track_deta);
0123     _h_backgroundCutTree->Branch("track_dphi", &_bb_track_dphi);
0124     _h_backgroundCutTree->Branch("track_dlayer", &_bb_track_dlayer);
0125     _h_backgroundCutTree->Branch("track_layer", &_bb_track_layer);
0126     _h_backgroundCutTree->Branch("track_pT", &_bb_track_pT);
0127     _h_backgroundCutTree->Branch("vtx_radius", &_bb_vtx_radius);
0128     _h_backgroundCutTree->Branch("vtx_chi2", &_bb_vtx_chi2);
0129     _h_backgroundCutTree->Branch("approach_dist", &_bb_approach);
0130     _h_backgroundCutTree->Branch("vtxTrack_dist", &_bb_vtxTrack_dist);
0131     _h_backgroundCutTree->Branch("photon_m", &_bb_photon_m);
0132     _h_backgroundCutTree->Branch("photon_pT", &_bb_photon_pT);
0133     _h_backgroundCutTree->Branch("cluster_prob", &_bb_cluster_prob);
0134     _h_backgroundCutTree->Branch("pid", &_bb_pid);
0135
0136     _e_backgroundCutTree = new TTree("cutTreeBacke","background data for making track pair cuts");
0137     _e_backgroundCutTree->SetAutoSave(300);
0138     _e_backgroundCutTree->Branch("track_deta", &_bb_track_deta);
0139     _e_backgroundCutTree->Branch("track_dphi", &_bb_track_dphi);
0140     _e_backgroundCutTree->Branch("track_dlayer", &_bb_track_dlayer);
0141     _e_backgroundCutTree->Branch("track_layer", &_bb_track_layer);
0142     _e_backgroundCutTree->Branch("track_pT", &_bb_track_pT);
0143     _e_backgroundCutTree->Branch("vtx_radius", &_bb_vtx_radius);
0144     _e_backgroundCutTree->Branch("vtx_chi2", &_bb_vtx_chi2);
0145     _e_backgroundCutTree->Branch("approach_dist", &_bb_approach);
0146     _e_backgroundCutTree->Branch("vtxTrack_dist", &_bb_vtxTrack_dist);
0147     _e_backgroundCutTree->Branch("photon_m", &_bb_photon_m);
0148     _e_backgroundCutTree->Branch("photon_pT", &_bb_photon_pT);
0149     _e_backgroundCutTree->Branch("cluster_prob", &_bb_cluster_prob);
0150   }
0151   return 0;
0152 }
0153
0154 void VtxTest::doNodePointers(PHCompositeNode* topNode){
0155   _mainClusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0156   _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0157   _clusterMap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0158   //  _hitMap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");
0159   //if(!_mainClusterContainer||!_truthinfo||!_clusterMap||!_hitMap){
0160   if(!_mainClusterContainer||!_truthinfo||!_clusterMap){
0161     cerr<<Name()<<": critical error-bad nodes \n";
0162     if(!_mainClusterContainer){
0163       cerr<<"\t RawClusterContainer is bad";
0164     }
0165     if(!_truthinfo){
0166       cerr<<"\t TruthInfoContainer is bad";
0167     }
0168     if(!_clusterMap){
0169       cerr<<"\t TrkrClusterMap is bad";
0170     }
0171     /*if(!_hitMap){
0172       cerr<<"\t SvtxHitMap is bad";
0173       }*/
0174     cerr<<endl;
0175   }
0176 }
0177
0178 int VtxTest::process_event(PHCompositeNode *topNode)
0179 {
0180   doNodePointers(topNode);
0181   _vertexer->InitEvent(topNode);
0182   PHG4TruthInfoContainer::Range range = _truthinfo->GetParticleRange(); //look at all truth particles
0183   SvtxEvalStack *stack = new SvtxEvalStack(topNode); //truth tracking info
0184   SvtxTrackEval* trackeval = stack->get_track_eval();
0185   SvtxTrack *svtxtrack1=NULL;
0186   SvtxTrack *svtxtrack2=NULL;
0187   PHG4VtxPoint *truth_vtx=NULL;
0188   for ( PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter ) {
0189     PHG4Particle* g4particle = iter->second;
0190     if(g4particle&&g4particle->get_pid()==211&&!_truthinfo->GetParticle(g4particle->get_parent_id())){
0191       cout<<"found pion\n";
0192       TLorentzVector truth;
0193       truth.SetPxPyPzE(g4particle->get_px(),g4particle->get_py(),g4particle->get_pz(),g4particle->get_e());
0194       if(!svtxtrack1)svtxtrack1=trackeval->best_track_from(g4particle);
0195       else svtxtrack2=trackeval->best_track_from(g4particle);
0196       if(truth_vtx==NULL) truth_vtx=_truthinfo->GetVtx(g4particle->get_vtx_id()); //get the vertex
0197       else{
0198         PHG4VtxPoint *vtemp = truth_vtx;
0199         truth_vtx=_truthinfo->GetVtx(g4particle->get_vtx_id());
0200         if(!(*truth_vtx==*vtemp)) cout<<"vtx does not agree"<<endl;
0201       }
0202     }
0203   }
0204   if(!truth_vtx){
0205     cout<<"did not find truth_vtx"<<endl;
0206     return Fun4AllReturnCodes::ABORTEVENT;
0207   }
0208   if(!svtxtrack1||!svtxtrack2){
0209     cout<<"did not find tracks"<<endl;
0210     return Fun4AllReturnCodes::ABORTEVENT;
0211   }
0212   genfit::GFRaveVertex* recoVert=_vertexer->findSecondaryVertex(svtxtrack1,svtxtrack2);
0213   cout<<"here"<<endl;
0214   if (recoVert)
0215   {
0216     TVector3 recoVertPos = recoVert->getPos();
0217     _b_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
0218     _b_tvtx_radius = sqrt(truth_vtx->get_x()*truth_vtx->get_x()+truth_vtx->get_y()*truth_vtx->get_y());
0219     _b_vtx_phi = recoVertPos.Phi();
0220     _b_vtx_eta = recoVertPos.Eta();
0221     _b_vtx_z = recoVertPos.Z();
0222     _b_vtx_x = recoVertPos.X();
0223     _b_vtx_y = recoVertPos.Y();
0224     TVector3 tVertPos(truth_vtx->get_x(),truth_vtx->get_y(),truth_vtx->get_z());
0225     _b_tvtx_phi = tVertPos.Phi();
0226     _b_tvtx_eta = tVertPos.Eta();
0227     _b_tvtx_z = tVertPos.Z();
0228     _b_tvtx_x = tVertPos.X();
0229     _b_tvtx_y = tVertPos.Y();
0230     _b_vtx_chi2 = recoVert->getChi2();
0231   }
0232   else{
0233     _b_vtx_radius = 0;
0234     _b_vtx_phi = 0;
0235     _b_vtx_eta = 0;
0236     _b_vtx_z = 0;
0237     _b_vtx_x = 0;
0238     _b_vtx_y = 0;
0239     _b_tvtx_radius = sqrt(truth_vtx->get_x()*truth_vtx->get_x()+truth_vtx->get_y()*truth_vtx->get_y());
0240     TVector3 tVertPos(truth_vtx->get_x(),truth_vtx->get_y(),truth_vtx->get_z());
0241     _b_tvtx_phi = tVertPos.Phi();
0242     _b_tvtx_eta = tVertPos.Eta();
0243     _b_tvtx_z = tVertPos.Z();
0244     _b_tvtx_x = tVertPos.X();
0245     _b_tvtx_y = tVertPos.Y();
0246   }
0247   delete stack;
0248   _signalCutTree->Fill();
0249   return 0;
0250 }
0251
0252 std::queue<std::pair<int,int>> VtxTest::numUnique(std::map<int,Conversion> *mymap=NULL,SvtxTrackEval* trackeval=NULL,RawClusterContainer *mainClusterContainer=NULL){
0253   _b_nVtx = 0;
0254   _b_Tpair=0;
0255   _b_Rpair=0;
0256   std::queue<std::pair<int,int>> missingChildren;//this is deprecated
0257   for (std::map<int,Conversion>::iterator i = mymap->begin(); i != mymap->end(); ++i) {
0258     PHG4VtxPoint *vtx =i->second.getVtx(); //get the vtx
0259     PHG4Particle *temp = i->second.getPhoton(); //set the photon
0260     TLorentzVector tlv_photon,tlv_electron,tlv_positron; //make tlv for each particle
0261     tlv_photon.SetPxPyPzE(temp->get_px(),temp->get_py(),temp->get_pz(),temp->get_e()); //intialize
0262     temp=i->second.getElectron(); //set the first child
0263     tlv_electron.SetPxPyPzE(temp->get_px(),temp->get_py(),temp->get_pz(),temp->get_e());
0264     if(_kMakeTTree){//fill tree values
0265       cout<<"numUnique filling tree"<<endl;
0266       _b_rVtx[_b_nVtx] = sqrt(vtx->get_x()*vtx->get_x()+vtx->get_y()*vtx->get_y()); //find the radius
0267       _b_parent_pt[_b_nVtx] =tlv_photon.Pt();
0268       _b_parent_phi[_b_nVtx]=tlv_photon.Phi();
0269       _b_parent_eta[_b_nVtx]=tlv_photon.Eta();
0270       _b_grandparent_id[_b_nVtx]=i->second.getSourceId();
0271       _b_e_deta[_b_nVtx]=-1.;
0272       _b_e_dphi[_b_nVtx]=-1.;
0273       _b_positron_pt[_b_nVtx]=-1.;
0274       _b_electron_pt[_b_nVtx]=tlv_electron.Pt();
0275     }
0276     temp=i->second.getPositron();
0277     if(temp){ //this will be false for conversions with 1 truth track
0278       tlv_positron.SetPxPyPzE(temp->get_px(),temp->get_py(),temp->get_pz(),temp->get_e()); //init the tlv
0279       if(_kMakeTTree) _b_positron_pt[_b_nVtx]=tlv_positron.Pt(); //fill tree
0280       if (TMath::Abs(tlv_electron.Eta())<_kRAPIDITYACCEPT&&TMath::Abs(tlv_positron.Eta())<_kRAPIDITYACCEPT)
0281       {
0282         unsigned int nRecoTracks = i->second.setRecoTracks(trackeval); //find the reco tracks for this conversion
0283         if(_kMakeTTree){
0284           _b_e_deta[_b_Tpair]=TMath::Abs(tlv_electron.Eta()-tlv_positron.Eta());
0285           _b_e_dphi[_b_Tpair]=TMath::Abs(tlv_electron.Phi()-tlv_positron.Phi());
0286           _b_fLayer[_b_Tpair]=_b_track_layer = i->second.firstLayer(_clusterMap);
0287           _b_fLayer[_b_Tpair]=-1;
0288           _b_Tpair++;
0289         }//tree
0290         switch(nRecoTracks)
0291         {
0292           case 2: //there are 2 reco tracks
0293             {
0294               if(_kMakeTTree){
0295                 _b_track_deta = i->second.trackDEta();
0296                 _b_track_dphi = i->second.trackDPhi();
0297                 _b_track_dlayer = i->second.trackDLayer(_clusterMap);
0298                 _b_track_dlayer=-1;
0299                 _b_track_pT = i->second.minTrackpT();
0300                 if(tlv_electron.Pt()>tlv_positron.Pt()) _b_ttrack_pT = tlv_positron.Pt();
0301                 else _b_ttrack_pT = tlv_electron.Pt();
0302                 _b_approach = i->second.approachDistance();
0303                 pair<SvtxTrack*, SvtxTrack*> reco_tracks=i->second.getRecoTracks();
0304                 genfit::GFRaveVertex* recoVert = _vertexer->findSecondaryVertex(reco_tracks.first,reco_tracks.second);
0305                 if (recoVert)
0306                 {
0307                   TVector3 recoVertPos = recoVert->getPos();
0308                   _b_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
0309                   _b_tvtx_radius = _b_rVtx[_b_nVtx];
0310                   _b_vtx_phi = recoVertPos.Phi();
0311                   _b_vtx_eta = recoVertPos.Eta();
0312                   _b_vtx_z = recoVertPos.Z();
0313                   _b_vtx_x = recoVertPos.X();
0314                   _b_vtx_y = recoVertPos.Y();
0315                   TVector3 tVertPos(vtx->get_x(),vtx->get_y(),vtx->get_z());
0316                   _b_tvtx_phi = tVertPos.Phi();
0317                   _b_tvtx_eta = tVertPos.Eta();
0318                   _b_tvtx_z = tVertPos.Z();
0319                   _b_tvtx_x = tVertPos.X();
0320                   _b_tvtx_y = tVertPos.Y();
0321                   _b_vtx_chi2 = recoVert->getChi2();
0322                   _b_vtxTrack_dist = i->second.dist(&recoVertPos,_clusterMap);
0323                 }
0324
0325                 TLorentzVector* recoPhoton = i->second.setRecoPhoton();
0326                 if (recoPhoton)
0327                 {
0328                   _b_photon_m=recoPhoton->Dot(*recoPhoton);
0329                   _b_photon_pT=recoPhoton->Pt();
0330                 }
0331                 else{ //photon was not reconstructed
0332                   _b_photon_m =-1;
0333                   _b_photon_pT=-1;
0334                 }
0335                 //reset the values
0336                 _b_cluster_prob=-1;
0337                 _b_cluster_deta[_b_Rpair]=-1.;
0338                 _b_cluster_dphi[_b_Rpair]=-1.;
0339                 _b_nCluster[_b_Rpair]=0;
0340                 _b_Scluster_prob[_b_Rpair]=-1;
0341                 _b_Mcluster_prob[_b_Rpair]=-1;
0342
0343               }
0344               pair<int,int> clusterIds = i->second.get_cluster_ids();
0345               RawCluster *clustemp;
0346               if(mainClusterContainer->getCluster(clusterIds.first)){//if there is matching cluster
0347                 clustemp =   dynamic_cast<RawCluster*>(mainClusterContainer->getCluster(clusterIds.first)->Clone());
0348                 _conversionClusters.AddCluster(clustemp); //add the calo cluster to the container
0349                 if (_kMakeTTree)
0350                 {
0351                   _b_cluster_prob=clustemp->get_prob();
0352                   RawCluster *clus2 = mainClusterContainer->getCluster(clusterIds.second);
0353                   if (clus2)
0354                   {
0355                     _b_cluster_dphi[_b_Rpair]=fabs(clustemp->get_phi()-clus2->get_phi());
0356                     TVector3 etaCalc(clustemp->get_x(),clustemp->get_y(),clustemp->get_z());
0357                     if (clus2->get_prob()>_b_cluster_prob)
0358                     {
0359                       _b_cluster_prob=clus2->get_prob();
0360                     }
0361                     //calculate deta
0362                     float eta1 = etaCalc.PseudoRapidity();
0363                     etaCalc.SetXYZ(clus2->get_x(),clus2->get_y(),clus2->get_z());
0364                     _b_cluster_deta[_b_Rpair]=fabs(eta1-etaCalc.PseudoRapidity());
0365                     /* this may be a logical bug. What are S and M? If this is false the identical cluster prob is recorded twice.*/
0366                     if (clusterIds.first!=clusterIds.second) //if there are two district clusters
0367                     {
0368                       _b_nCluster[_b_Rpair]=2;
0369                       _b_Scluster_prob[_b_Rpair]=clustemp->get_prob();
0370                     }
0371                     else{
0372                       _b_nCluster[_b_Rpair]=1;
0373                       _b_Mcluster_prob[_b_Rpair]=clustemp->get_prob();
0374                     }
0375                   }
0376                 }
0377               }
0378               _signalCutTree->Fill();
0379               _b_Rpair++;
0380               break;
0381             }
0382           case 1: //there's one reco track
0383             {
0384               int clustidtemp=i->second.get_cluster_id(); //get the cluster id of the current conversion
0385               if(mainClusterContainer->getCluster(clustidtemp)){//if thre is matching cluster
0386                 RawCluster *clustemp =   dynamic_cast<RawCluster*>(mainClusterContainer->getCluster(clustidtemp)->Clone());
0387                 _conversionClusters.AddCluster(clustemp); //add the calo cluster to the container
0388               }
0389               //cout<<"Matched 1 reco with layer="<<i->second.firstLayer(_clusterMap)<<"pTs:"<<tlv_electron.Pt()<<"-"<<tlv_positron.Pt()<<'\n';
0390               break;
0391             }
0392           case 0: //no reco tracks
0393             break;
0394           default:
0395             if (Verbosity()>1)
0396             {
0397               cerr<<Name()<<" error setting reco tracks"<<endl;
0398             }
0399             break;
0400         }//switch
0401       }//rapidity cut
0402       if (_kMakeTTree)
0403       {
0404         _b_pythia[_b_nVtx]=i->second.getEmbed()==_kPythiaEmbed;
0405         _b_nVtx++;
0406       }
0407     }// has 2 truth tracks
0408   }//map loop
0409   return missingChildren;
0410 }
0411
0412 //only call if _kMakeTTree is true
0413 void VtxTest::processBackground(std::map<int,Conversion> *mymap,SvtxTrackEval* trackeval,TTree* tree){
0414   for (std::map<int,Conversion>::iterator i = mymap->begin(); i != mymap->end(); ++i) {
0415     int nReco=i->second.setRecoTracks(trackeval);
0416     if (nReco==2)
0417     {
0418       _bb_track_deta = i->second.trackDEta();
0419       _bb_track_dphi = i->second.trackDPhi();
0420       _bb_track_dlayer = i->second.trackDLayer(_clusterMap);
0421       _bb_track_layer = i->second.firstLayer(_clusterMap);
0422       _bb_track_layer=-1;
0423       _bb_track_dlayer=-1;
0424       _bb_track_pT = i->second.minTrackpT();
0425       _bb_approach = i->second.approachDistance();
0426       _bb_pid = i->second.getElectron()->get_pid();
0427       pair<SvtxTrack*, SvtxTrack*> reco_tracks=i->second.getRecoTracks();
0428       genfit::GFRaveVertex* recoVert = _vertexer->findSecondaryVertex(reco_tracks.first,reco_tracks.second);
0429       if (recoVert)
0430       {
0431         TVector3 recoVertPos = recoVert->getPos();
0432         _bb_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
0433         _bb_vtx_chi2 = recoVert->getChi2();
0434         _bb_vtxTrack_dist = i->second.dist(&recoVertPos,_clusterMap);
0435       }
0436       TLorentzVector* recoPhoton = i->second.setRecoPhoton();
0437       if (recoPhoton)
0438       {
0439         _bb_photon_m=recoPhoton->Dot(*recoPhoton);
0440         _bb_photon_pT=recoPhoton->Pt();
0441       }
0442       else{
0443         _bb_photon_m =0;
0444         _bb_photon_pT=0;
0445       }
0446       _bb_cluster_prob= _mainClusterContainer->getCluster(i->second.get_cluster_id())->get_prob();
0447       tree->Fill();
0448     }
0449   }
0450 }
0451
0452 void VtxTest::findChildren(std::queue<std::pair<int,int>> missingChildren,PHG4TruthInfoContainer* _truthinfo){
0453   if (Verbosity()>=6)
0454   {
0455     while(!missingChildren.empty()){
0456       for (PHG4TruthInfoContainer::ConstIterator iter = _truthinfo->GetParticleRange().first; iter != _truthinfo->GetParticleRange().second; ++iter)
0457       {
0458         if(Verbosity()>1&& iter->second->get_parent_id()==missingChildren.front().first&&iter->second->get_track_id()!=missingChildren.front().second){
0459           cout<<"Found Child:\n";
0460           iter->second->identify();
0461           cout<<"With mother:\n";
0462
0463         }
0464       }
0465       missingChildren.pop();
0466     }
0467   }
0468 }
0469
0470 const RawClusterContainer* VtxTest::getClusters()const {return &_conversionClusters;}
0471
0472 int VtxTest::get_embed(PHG4Particle* particle, PHG4TruthInfoContainer* truthinfo)const{
0473   return truthinfo->isEmbeded(particle->get_track_id());
0474 }
0475
0476 float VtxTest::vtoR(PHG4VtxPoint* vtx)const{
0477   return (float) sqrt(vtx->get_x()*vtx->get_x()+vtx->get_y()*vtx->get_y());
0478 }
0479
0480
0481 int VtxTest::End(PHCompositeNode *topNode)
0482 {
0483   if(_kMakeTTree){
0484     _f->Write();
0485     _f->Close();
0486   }
0487   return 0;
0488 }