PhotonConversion/truthconversion/TruthConversionEval.C

0001 #include "TruthConversionEval.h"
0002 #include "Conversion.h"
0003 #include "SVReco.h"
0004 #include "VtxRegressor.h"
0005 //#include "../PhotonConversion/RaveVertexingAux.h"
0006
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/getClass.h>
0009
0010 #include <calotrigger/CaloTriggerInfo.h>
0011 #include <calobase/RawCluster.h>
0012 #include <calobase/RawClusterv1.h>
0013
0014 #include <g4main/PHG4TruthInfoContainer.h>
0015 #include <g4main/PHG4Particlev1.h>
0016 #include <g4main/PHG4Particlev2.h>
0017 #include <g4main/PHG4VtxPoint.h>
0018
0019 /*#include <trackbase_historic/SvtxHitMap.h>
0020 #include <trackbase_historic/SvtxHit.h>
0021 #include <trackbase_historic/SvtxClusterMap.h>
0022 #include <trackbase_historic/SvtxCluster.h>*/
0023 #include <trackbase_historic/SvtxVertex.h>
0024 #include <trackbase_historic/SvtxVertexMap.h>
0025 #include <trackbase_historic/SvtxTrackMap.h>
0026 #include <trackbase/TrkrClusterContainer.h>
0027 #include <trackbase/TrkrCluster.h>
0028
0029 #include <g4eval/SvtxEvalStack.h>
0030 #include <g4eval/SvtxTrackEval.h>
0031
0032 //#include <GenFit/GFRaveConverters.h>
0033 #include <GenFit/FieldManager.h>
0034 #include <GenFit/GFRaveVertex.h>
0035 #include <GenFit/GFRaveVertexFactory.h>
0036 #include <GenFit/MeasuredStateOnPlane.h>
0037 #include <GenFit/RKTrackRep.h>
0038 #include <GenFit/StateOnPlane.h>
0039 #include <GenFit/Track.h>
0040
0041 #include <phgenfit/Track.h>
0042
0043 #include <fun4all/Fun4AllReturnCodes.h>
0044
0045 #include <TFile.h>
0046 #include <TTree.h>
0047
0048 #include <math.h>
0049 #include <utility>
0050 #include <list>
0051 #include <iostream>
0052
0053 TruthConversionEval::TruthConversionEval(const std::string &name, unsigned int runnumber,
0054         int particleEmbed,  int pythiaEmbed,bool makeTTree=true,string TMVAName="",string TMVAPath="") : SubsysReco("TruthConversionEval"),
0055     _kRunNumber(runnumber),_kParticleEmbed(particleEmbed), _kPythiaEmbed(pythiaEmbed), _kMakeTTree(makeTTree)
0056 {
0057     _foutname = name;
0058     _regressor = new VtxRegressor(TMVAName,TMVAPath);
0059 }
0060
0061 TruthConversionEval::~TruthConversionEval(){
0062     cout<<"destruct"<<endl;
0063     if (_f) delete _f;
0064     if (_vertexer) delete _vertexer;
0065     if(_regressor) delete _regressor;
0066 }
0067
0068 int TruthConversionEval::InitRun(PHCompositeNode *topNode)
0069 {
0070     _vertexer = new SVReco();
0071     _vertexer->InitRun(topNode);
0072     if(_kMakeTTree){
0073         _runNumber=_kRunNumber;
0074         _f = new TFile( _foutname.c_str(), "RECREATE");
0075         _observTree = new TTree("observTree","per event observables");
0076         _observTree->SetAutoSave(300);
0077         _observTree->Branch("nMatched", &_b_nMatched);
0078         _observTree->Branch("nUnmatched", &_b_nUnmatched);
0079         _observTree->Branch("truth_pT", &_b_truth_pT);
0080         _observTree->Branch("reco_pT", &_b_reco_pT);
0081         _observTree->Branch("track_pT",&_b_alltrack_pT) ;
0082         _observTree->Branch("photon_pT",&_b_allphoton_pT) ;
0083
0084         _vtxingTree = new TTree("vtxingTree","data predicting vtx from track pair");
0085         _vtxingTree->SetAutoSave(300);
0086         _vtxingTree->Branch("vtx_radius", &_b_vtx_radius);
0087         _vtxingTree->Branch("tvtx_radius", &_b_tvtx_radius);
0088         _vtxingTree->Branch("vtx_phi", &_b_vtx_phi);
0089         _vtxingTree->Branch("vtx_eta", &_b_vtx_eta);
0090         _vtxingTree->Branch("vtx_x", &_b_vtx_x);
0091         _vtxingTree->Branch("vtx_y", &_b_vtx_y);
0092         _vtxingTree->Branch("vtx_z", &_b_vtx_z);
0093         _vtxingTree->Branch("tvtx_eta", &_b_tvtx_eta);
0094         _vtxingTree->Branch("tvtx_x", &_b_tvtx_x);
0095         _vtxingTree->Branch("tvtx_y", &_b_tvtx_y);
0096         _vtxingTree->Branch("tvtx_z", &_b_tvtx_z);
0097         _vtxingTree->Branch("tvtx_phi", &_b_tvtx_phi);
0098         _vtxingTree->Branch("vtx_chi2", &_b_vtx_chi2);
0099         _vtxingTree->Branch("track1_pt", &_b_track1_pt);
0100         _vtxingTree->Branch("track1_eta",& _b_track1_eta);
0101         _vtxingTree->Branch("track1_phi",& _b_track1_phi);
0102         _vtxingTree->Branch("track2_pt", &_b_track2_pt);
0103         _vtxingTree->Branch("track2_eta",& _b_track2_eta);
0104         _vtxingTree->Branch("track2_phi",& _b_track2_phi);
0105         _vtxingTree->Branch("track_layer",& _b_track_layer);
0106
0107         _trackBackTree = new TTree("trackBackTree","track background all single tracks");
0108         _trackBackTree->SetAutoSave(300);
0109         _trackBackTree->Branch("track_dca", &_bb_track_dca);
0110         _trackBackTree->Branch("track_pT",  &_bb_track_pT);
0111         _trackBackTree->Branch("track_layer", &_bb_track_layer);
0112         _trackBackTree->Branch("cluster_prob", &_bb_cluster_prob);
0113
0114         _pairBackTree = new TTree("pairBackTree","pair background all possible combinations");
0115         _pairBackTree->SetAutoSave(300);
0116         _pairBackTree->Branch("track_deta", &_bb_track_deta);
0117         _pairBackTree->Branch("track2_pid", &_bb_track2_pid);
0118         _pairBackTree->Branch("track1_pid", &_bb_track1_pid);
0119         _pairBackTree->Branch("parent_pid", &_bb_parent_pid);
0120         _pairBackTree->Branch("track_dphi", &_bb_track_dphi);
0121         _pairBackTree->Branch("track_dlayer",&_bb_track_dlayer);
0122         _pairBackTree->Branch("approach_dist", &_bb_approach);
0123         _pairBackTree->Branch("track_dca", &_bb_track_dca);
0124         _pairBackTree->Branch("track_pT",  &_bb_track_pT);
0125         _pairBackTree->Branch("track_layer", &_bb_track_layer);
0126         _pairBackTree->Branch("cluster_prob", &_bb_cluster_prob);
0127         //_pairBackTree->Branch("nCluster", &_bb_nCluster);
0128         _pairBackTree->Branch("cluster_dphi", &_bb_cluster_dphi);
0129         _pairBackTree->Branch("cluster_deta", &_bb_cluster_deta);
0130
0131         _vtxBackTree = new TTree("vtxBackTree","events that pass track pair cuts");
0132         _vtxBackTree->SetAutoSave(300);
0133         _vtxBackTree->Branch("track_deta", &_bb_track_deta);
0134         _vtxBackTree->Branch("track1_pid", &_bb_track1_pid);
0135         _vtxBackTree->Branch("track2_pid", &_bb_track2_pid);
0136         _vtxBackTree->Branch("track_dphi", &_bb_track_dphi);
0137         _vtxBackTree->Branch("track_dlayer",&_bb_track_dlayer);
0138         _vtxBackTree->Branch("approach_dist", &_bb_approach);
0139         _vtxBackTree->Branch("track_dca", &_bb_track_dca);
0140         _vtxBackTree->Branch("track_pT",  &_bb_track_pT);
0141         _vtxBackTree->Branch("track_layer", &_bb_track_layer);
0142         _vtxBackTree->Branch("cluster_prob", &_bb_cluster_prob);
0143         _vtxBackTree->Branch("vtx_radius",&_bb_vtx_radius);
0144         _vtxBackTree->Branch("photon_m",&_bb_photon_m);
0145         _vtxBackTree->Branch("photon_pT",&_bb_photon_pT);
0146         //_vtxBackTree->Branch("nCluster", &_bb_nCluster);
0147         _vtxBackTree->Branch("cluster_dphi", &_bb_cluster_dphi);
0148         _vtxBackTree->Branch("cluster_deta", &_bb_cluster_deta);
0149
0150         _signalCutTree = new TTree("cutTreeSignal","signal data for making track pair cuts");
0151         _signalCutTree->SetAutoSave(100);
0152         _signalCutTree->Branch("track_deta", &_b_track_deta);
0153         _signalCutTree->Branch("track_dca", &_b_track_dca);
0154         _signalCutTree->Branch("track_dphi", &_b_track_dphi);
0155         _signalCutTree->Branch("track_dlayer",&_b_track_dlayer);
0156         _signalCutTree->Branch("track_layer", &_b_track_layer);
0157         _signalCutTree->Branch("track_pT", &_b_track_pT);
0158         //_signalCutTree->Branch("ttrack_pT", &_b_ttrack_pT);
0159         _signalCutTree->Branch("approach_dist", &_b_approach);
0160         _signalCutTree->Branch("vtx_radius", &_b_vtx_radius);
0161         _signalCutTree->Branch("vtx_phi", &_b_vtx_phi);
0162         _signalCutTree->Branch("tvtx_phi", &_b_tvtx_phi);
0163         _signalCutTree->Branch("tvtx_radius", &_b_tvtx_radius);
0164         _signalCutTree->Branch("vtx_chi2", &_b_vtx_chi2);
0165         //_signalCutTree->Branch("vtxTrackRZ_dist", &_b_vtxTrackRZ_dist);
0166         //_signalCutTree->Branch("vtxTrackRPhi_dist", &_b_vtxTrackRPhi_dist);
0167         _signalCutTree->Branch("photon_m", &_b_photon_m);
0168         _signalCutTree->Branch("tphoton_m", &_b_tphoton_m);
0169         _signalCutTree->Branch("photon_pT", &_b_photon_pT);
0170         _signalCutTree->Branch("tphoton_pT", &_b_tphoton_pT);
0171         _signalCutTree->Branch("cluster_prob", &_b_cluster_prob);
0172         //_signalCutTree->Branch("nCluster", &_b_nCluster);
0173         _signalCutTree->Branch("cluster_dphi", &_b_cluster_dphi);
0174         _signalCutTree->Branch("cluster_deta", &_b_cluster_deta);
0175     }
0176     return 0;
0177 }
0178
0179 bool TruthConversionEval::doNodePointers(PHCompositeNode* topNode){
0180     bool goodPointers=true;
0181     _mainClusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0182     _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0183     _clusterMap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0184     _allTracks = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0185     //  _hitMap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");
0186     //if(!_mainClusterContainer||!_truthinfo||!_clusterMap||!_hitMap){
0187     if(!_mainClusterContainer||!_truthinfo||!_clusterMap||!_allTracks){
0188         cerr<<Name()<<": critical error-bad nodes \n";
0189         if(!_mainClusterContainer){
0190             cerr<<"\t RawClusterContainer is bad";
0191         }
0192         if(!_truthinfo){
0193             cerr<<"\t TruthInfoContainer is bad";
0194         }
0195         if(!_clusterMap){
0196             cerr<<"\t TrkrClusterMap is bad";
0197         }
0198         if(!_allTracks){
0199             cerr<<"\t SvtxTrackMap is bad";
0200         }
0201         cerr<<endl;
0202         goodPointers=false;
0203     }
0204     return goodPointers;
0205 }
0206
0207 SvtxVertex* TruthConversionEval::get_primary_vertex(PHCompositeNode *topNode)const{
0208     SvtxVertexMap *vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0209     return vertexMap->get(0);
0210 }
0211
0212 int TruthConversionEval::process_event(PHCompositeNode *topNode)
0213 {
0214     if(!doNodePointers(topNode)) return Fun4AllReturnCodes::ABORTEVENT;
0215     cout<<"init vertexer event"<<endl;
0216     _vertexer->InitEvent(topNode);
0217     _conversionClusters.Reset(); //clear the list of conversion clusters
0218     PHG4TruthInfoContainer::Range range = _truthinfo->GetParticleRange(); //look at all truth particles
0219     SvtxEvalStack *stack = new SvtxEvalStack(topNode); //truth tracking info
0220     SvtxTrackEval* trackeval = stack->get_track_eval();
0221     if (!trackeval)
0222     {
0223         cerr<<"NULL track eval in " <<Name()<<" fatal error"<<endl;
0224         return 1;
0225     }
0226     //make a map of the conversions
0227     std::map<int,Conversion> mapConversions;
0228     //h is for hadronic e is for EM
0229     std::vector<SvtxTrack*> backgroundTracks;
0230     std::vector<std::pair<SvtxTrack*,SvtxTrack*>> tightbackgroundTrackPairs; //used to find the pair for unmatched conversion tracks
0231     std::vector<PHG4Particle*> truthPhotons;
0232     std::list<int> signalTracks;
0233     //reset obervation variables
0234     _b_nMatched=0;
0235     _b_nUnmatched=0;
0236     _b_truth_pT.clear();
0237     _b_reco_pT.clear();
0238     _b_alltrack_pT.clear();
0239     _b_allphoton_pT.clear();
0240     cout<<"init truth loop"<<endl;
0241
0242     for ( PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter ) {
0243         PHG4Particle* g4particle = iter->second;
0244         if(g4particle->get_pid()==22&&g4particle->get_track_id()==g4particle->get_primary_id()){
0245             _b_allphoton_pT.push_back(sqrt(g4particle->get_px()*g4particle->get_px()+g4particle->get_py()*g4particle->get_py()));
0246             truthPhotons.push_back(g4particle);
0247         }
0248         PHG4Particle* parent =_truthinfo->GetParticle(g4particle->get_parent_id());
0249         PHG4VtxPoint* vtx=_truthinfo->GetVtx(g4particle->get_vtx_id()); //get the vertex
0250         if(!vtx){
0251             cerr<<"null vtx primaryid="<<g4particle->get_primary_id()<<'\n';
0252             g4particle->identify();
0253             cout<<endl;
0254             continue;
0255         }
0256         float radius=sqrt(vtx->get_x()*vtx->get_x()+vtx->get_y()*vtx->get_y());
0257         //if outside the tracker skip this
0258         if(radius>s_kTPCRADIUS) continue;
0259         int embedID;
0260         if (parent)//if the particle is not primary
0261         {
0262             embedID=get_embed(parent,_truthinfo);
0263             float truthpT = sqrt(g4particle->get_px()*g4particle->get_px()+g4particle->get_py()*g4particle->get_py());
0264             if(parent->get_pid()==22&&TMath::Abs(g4particle->get_pid())==11&&truthpT>2.5){ //conversion check
0265                 if (Verbosity()==10)
0266                 {
0267                     std::cout<<"Conversion with radius [cm]:"<<radius<<'\n';
0268                 }
0269                 //initialize the conversion object -don't use constructor b/c setters have error handling
0270                 (mapConversions[vtx->get_id()]).setElectron(g4particle);
0271                 (mapConversions[vtx->get_id()]).setVtx(vtx);
0272                 (mapConversions[vtx->get_id()]).setPrimaryPhoton(parent,_truthinfo);
0273                 (mapConversions[vtx->get_id()]).setEmbed(embedID);
0274                 PHG4Particle* grand =_truthinfo->GetParticle(parent->get_parent_id()); //grandparent
0275                 if (grand) (mapConversions[vtx->get_id()]).setSourceId(grand->get_pid());//record pid of the photon's source
0276                 else (mapConversions[vtx->get_id()]).setSourceId(0);//or it is from the G4 generator
0277                 _b_truth_pT.push_back(truthpT);
0278                 //build a list of the ids
0279                 SvtxTrack* recoTrack = trackeval->best_track_from(g4particle);
0280                 if(recoTrack){
0281                     signalTracks.push_back(recoTrack->get_id());
0282                     _b_reco_pT.push_back(recoTrack->get_pt());
0283                     cout<<"matched truth track"<<endl;
0284                     _b_nMatched++;
0285                 }
0286                 else{
0287                     _b_nUnmatched++;
0288                     cerr<<"WARNING no matching track for conversion"<<endl;
0289                 }
0290             }
0291         }// not primary
0292     }//truth particle loop
0293     signalTracks.sort();
0294     cout<<"intit track loop"<<endl;
0295     //build the current backgroundSample
0296     for ( SvtxTrackMap::Iter iter = _allTracks->begin(); iter != _allTracks->end(); ++iter) {
0297         _b_alltrack_pT.push_back(iter->second->get_pt());
0298         auto inCheck = std::find(signalTracks.begin(),signalTracks.end(),iter->first);
0299         //if the track is not in the list of signal tracks
0300         if (inCheck==signalTracks.end())
0301         {
0302             TLorentzVector track_tlv = tracktoTLV(iter->second);
0303             bool addTrack=true;
0304             //make sure the track is not too close to a photon
0305             for (std::vector<PHG4Particle*>::iterator iPhoton = truthPhotons.begin(); iPhoton != truthPhotons.end()&&addTrack; ++iPhoton)
0306             {
0307                 TLorentzVector photon_tlv = particletoTLV(*iPhoton);
0308                 if (photon_tlv.DeltaR(track_tlv)<.2)
0309                 {
0310                     addTrack=false;
0311                 }
0312             }
0313             if(addTrack){
0314                 backgroundTracks.push_back(iter->second);
0315                 //see if it is a good tight background cannidate
0316                 for ( SvtxTrackMap::Iter jter = std::next(iter,1); jter != _allTracks->end()&&iter->second->get_pt()>_kTightPtMin; ++jter){
0317                     if(jter->second->get_pt()>_kTightPtMin&&TMath::Abs(jter->second->get_eta()-iter->second->get_eta())<_kTightDetaMax&&
0318                             iter->second->get_charge()==jter->second->get_charge()*-1){
0319                         tightbackgroundTrackPairs.push_back(std::pair<SvtxTrack*,SvtxTrack*>(iter->second,jter->second));
0320                     }
0321                 }
0322
0323             }
0324         }
0325     }
0326     //pass the map to this helper method which fills the fields for the signal ttree
0327     numUnique(&mapConversions,trackeval,_mainClusterContainer,&tightbackgroundTrackPairs);
0328
0329     if (_kMakeTTree)
0330     {
0331         cout<<"intit background process"<<endl;
0332         if(_b_truth_pT.size()!=0) _observTree->Fill();
0333         processTrackBackground(&backgroundTracks,trackeval);
0334         cout<<"finish back"<<endl;
0335     }
0336     delete stack;
0337     cout<<"finish process"<<endl;
0338     return 0;
0339 }
0340
0341 void TruthConversionEval::numUnique(std::map<int,Conversion> *mymap=NULL,SvtxTrackEval* trackeval=NULL,RawClusterContainer *mainClusterContainer=NULL,
0342         std::vector<std::pair<SvtxTrack*,SvtxTrack*>>* backgroundTracks=NULL){
0343     cout<<"start conversion analysis loop"<<endl;
0344     for (std::map<int,Conversion>::iterator i = mymap->begin(); i != mymap->end(); ++i) {
0345         TLorentzVector tlv_photon,tlv_electron,tlv_positron; //make tlv for each particle
0346         PHG4Particle *photon = i->second.getPrimaryPhoton(); //set the photon
0347
0348         if(photon)tlv_photon.SetPxPyPzE(photon->get_px(),photon->get_py(),photon->get_pz(),photon->get_e()); //intialize
0349         else cerr<<"No truth photon for conversion"<<endl;
0350         PHG4Particle *e1=i->second.getElectron(); //set the first child
0351         if(e1){
0352             tlv_electron.SetPxPyPzE(e1->get_px(),e1->get_py(),e1->get_pz(),e1->get_e());
0353         }
0354         else cerr<<"No truth electron for conversion"<<endl;
0355         PHG4Particle *e2=i->second.getPositron();
0356         if(e2){ //this will be false for conversions with 1 truth track
0357             tlv_positron.SetPxPyPzE(e2->get_px(),e2->get_py(),e2->get_pz(),e2->get_e()); //init the tlv
0358             if (TMath::Abs(tlv_electron.Eta())<_kRAPIDITYACCEPT&&TMath::Abs(tlv_positron.Eta())<_kRAPIDITYACCEPT)
0359             {
0360                 unsigned int nRecoTracks = i->second.setRecoTracks(trackeval); //find the reco tracks for this conversion
0361                 switch(nRecoTracks)
0362                 {
0363                     case 2: //there are 2 reco tracks
0364                         {
0365                             recordConversion(&(i->second),&tlv_photon,&tlv_electron,&tlv_positron);
0366                             break;
0367                         }
0368                     case 1: //there's one reco track try to find the other
0369                         {
0370                             PHG4Particle *truthe;
0371                             //get the truth particle that is missing a reco track
0372                             if(!i->second.getRecoTrack(e1->get_track_id()))truthe = e1;
0373                             else truthe=e2;
0374                             if(!truthe) cerr<<"critical error"<<endl;
0375                             //look through the background for the missing pair
0376                             for(auto pair : *backgroundTracks){
0377                                 if(!pair.first||!pair.second) cerr<<"critical error2"<<endl;
0378                                 if ((pair.first->get_charge()>0&&truthe->get_pid()<0)||(pair.first->get_charge()<0&&truthe->get_pid()>0))
0379                                 {
0380                                     TVector3 truth_tlv(truthe->get_px(),truthe->get_py(),truthe->get_pz());
0381                                     TVector3 reco_tlv(pair.first->get_px(),pair.first->get_py(),pair.first->get_pz());
0382                                     if (reco_tlv.DeltaR(truth_tlv)<.1)
0383                                     {
0384                                         i->second.setRecoTrack(truthe->get_track_id(),pair.first);
0385                                     }
0386                                 }
0387                                 else if ((pair.second->get_charge()>0&&truthe->get_pid()<0)||(pair.second->get_charge()<0&&truthe->get_pid()>0))
0388                                 {
0389                                     TVector3 truth_tlv(truthe->get_px(),truthe->get_py(),truthe->get_pz());
0390                                     TVector3 reco_tlv(pair.second->get_px(),pair.second->get_py(),pair.second->get_pz());
0391                                     if (reco_tlv.DeltaR(truth_tlv)<.1)
0392                                     {
0393                                         i->second.setRecoTrack(truthe->get_track_id(),pair.second);
0394                                     }
0395                                 }
0396                             }
0397                             recordConversion(&(i->second),&tlv_photon,&tlv_electron,&tlv_positron);
0398                             break;
0399                         }
0400                     case 0: //no reco tracks
0401                         //TODO maybe try to find the reco tracks? for now just record the bad info
0402                         recordConversion(&(i->second),&tlv_photon,&tlv_electron,&tlv_positron);
0403                         break;
0404                     default:
0405                         if (Verbosity()>1)
0406                         {
0407                             cerr<<Name()<<" error setting reco tracks"<<endl;
0408                         }
0409                         break;
0410                 }//switch
0411             }//rapidity cut
0412         }// has 2 truth tracks
0413         cout<<"map loop"<<endl;
0414     }//map loop
0415     cout<<"done num"<<endl;
0416 }
0417
0418 //only call if _kMakeTTree is true
0419 void TruthConversionEval::processTrackBackground(std::vector<SvtxTrack*> *v_tracks,SvtxTrackEval* trackeval){
0420     Conversion pairMath;
0421     float lastpT=-1.;
0422     unsigned nNullTrack=0;
0423     cout<<"The total possible background track count is "<<v_tracks->size()<<'\n';
0424     for (std::vector<SvtxTrack*>::iterator iter = v_tracks->begin(); iter != v_tracks->end(); ++iter) {
0425         cout<<"looping"<<endl;
0426         SvtxTrack* iTrack = *iter;
0427         //get the SvtxTrack it must not be NULL
0428         if(!iTrack){
0429             nNullTrack++;
0430             continue;
0431         }
0432
0433         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?
0434         lastpT=iTrack->get_pt();
0435         cout<<"\t pT="<<lastpT<<'\n';
0436
0437         auto temp_key_it=iTrack->begin_cluster_keys();//key iterator to first cluster
0438         if(temp_key_it!=iTrack->end_cluster_keys()){//if the track has clusters
0439             TrkrCluster* temp_cluster = _clusterMap->findCluster(*temp_key_it);//get the cluster
0440             if(temp_cluster) _bb_track_layer = TrkrDefs::getLayer(temp_cluster->getClusKey());//if there is a cluster record its layer
0441             else _bb_track_layer=-999;
0442         }
0443         //record track info
0444         _bb_track_dca = iTrack->get_dca();
0445         _bb_track_pT = iTrack->get_pt();
0446         //record the EMCal cluster info
0447         auto cluster1 = _mainClusterContainer->getCluster(iTrack->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)));
0448         if(cluster1) _bb_cluster_prob= cluster1->get_prob();
0449         else _bb_cluster_prob=-999;
0450         //pair with other tracks
0451         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)
0452             SvtxTrack* jTrack = *jter;
0453             if(!jTrack||TMath::Abs(jTrack->get_eta())>1.1)continue;
0454             //record pair geometry
0455             cout<<"calculations"<<endl;
0456             _bb_track_deta = pairMath.trackDEta(iTrack,jTrack);
0457             _bb_track_dphi = pairMath.trackDPhi(iTrack,jTrack);
0458             _bb_track_dlayer = pairMath.trackDLayer(_clusterMap,iTrack,jTrack);
0459             _bb_approach = pairMath.approachDistance(iTrack,jTrack);
0460             //record second EMCal cluster
0461             auto cluster2 = _mainClusterContainer->getCluster(iTrack->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)));
0462             //if the EMCal clusters can be found record their pair geometry
0463             if (cluster2&&cluster1)
0464             {
0465                 //if the two clusters are unique calculate the values
0466                 if(cluster2->get_id()!=cluster1->get_id())
0467                 {
0468                     _bb_nCluster = 2;
0469                     _bb_cluster_dphi=fabs(cluster1->get_phi()-cluster2->get_phi());
0470                     TVector3 etaCalc(cluster1->get_x(),cluster1->get_y(),cluster1->get_z());
0471                     float eta1 = etaCalc.PseudoRapidity();
0472                     etaCalc.SetXYZ(cluster2->get_x(),cluster2->get_y(),cluster2->get_z());
0473                     _bb_cluster_deta=fabs(eta1-etaCalc.PseudoRapidity());
0474                 }
0475                 else{ //if they are not unique then they are 0
0476                     _bb_nCluster = 1;
0477                     _bb_cluster_dphi=0;
0478                     _bb_cluster_deta=0;
0479                 }
0480             }
0481             else{ //clusters were not found
0482                 _bb_nCluster=0;
0483                 _bb_cluster_deta=-999;
0484                 _bb_cluster_dphi=-999;
0485             }
0486             /*_bb_track1_pid = (*iTruthTrack)->get_pid();
0487                 _bb_track2_pid = (*jTruthTrack)->get_pid();
0488                 PHG4Particle* parent  = _truthinfo->GetParticle((*iTruthTrack)->get_parent_id());
0489                 if(parent) _bb_parent_pid = parent->get_pid();
0490                 else _bb_parent_pid=0;*/
0491             //if the track pairs pass the pair cuts make them part of the vertex background
0492             if (_bb_track_layer>=0&&_bb_track_pT>_kTightPtMin&&_bb_track_deta<_kTightDetaMax&&TMath::Abs(_bb_track_dlayer)<9)
0493             {
0494                 //iTrack->identify();
0495                 //jTrack->identify();
0496                 cout<<"here vertex"<<endl;
0497                 //make the vertex
0498                 genfit::GFRaveVertex* recoVert = _vertexer->findSecondaryVertex(iTrack,jTrack);
0499                 cout<<"made vertex"<<endl;
0500                 //if the vertex is made record its info
0501                 if (recoVert)
0502                 {
0503                     recoVert=pairMath.correctSecondaryVertex(_regressor,recoVert,iTrack,jTrack);
0504                     cout<<"corrected vertex"<<endl;
0505                     TVector3 recoVertPos = recoVert->getPos();
0506                     cout<<"deleting"<<endl;
0507                     delete recoVert;
0508                     cout<<"deleted"<<endl;
0509                     //fill the tree values
0510                     _bb_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
0511                     _bb_vtx_chi2 = recoVert->getChi2();
0512                     _bb_vtxTrackRZ_dist = pairMath.vtxTrackRZ(recoVertPos,iTrack,jTrack);
0513                     _bb_vtxTrackRPhi_dist = pairMath.vtxTrackRPhi(recoVertPos,iTrack,jTrack);
0514                     //then make the photon
0515                     TLorentzVector* recoPhoton= pairMath.getRecoPhoton(iTrack,jTrack);
0516                     if(recoPhoton){
0517                         _bb_photon_m = recoPhoton->Dot(*recoPhoton);
0518                         _bb_photon_pT = recoPhoton->Pt();
0519                         delete recoPhoton;
0520                     }
0521                     else{
0522                         _bb_photon_m=-999;
0523                         _bb_photon_pT=-999;
0524                     }
0525                 }
0526                 else{ // no vtx could be found for thre track pair
0527                     _bb_vtx_radius = -999;
0528                     _bb_vtx_chi2 = -999;
0529                     _bb_vtxTrackRZ_dist =-999;
0530                     _bb_vtxTrackRPhi_dist =-999;
0531                 }
0532                 cout<<"fill vtx"<<endl;
0533                 _vtxBackTree->Fill();
0534             }//pair cuts
0535             _pairBackTree->Fill();
0536             cout<<"filled pair"<<endl;
0537         }//jTrack loop
0538         cout<<"filling track"<<endl;
0539         _trackBackTree->Fill();
0540         cout<<"filled track"<<endl;
0541     }//iTrack loop
0542     cout<<"Null track count ="<<nNullTrack<<endl;
0543 }
0544
0545 void TruthConversionEval::recordConversion(Conversion *conversion,TLorentzVector *tlv_photon,TLorentzVector *tlv_electron, TLorentzVector *tlv_positron){
0546     cout<<"recording"<<endl;
0547     if(tlv_photon&&tlv_electron&&tlv_positron&&conversion&&conversion->recoCount()==2){
0548         _b_track_deta = conversion->trackDEta();
0549         _b_track_dphi = conversion->trackDPhi();
0550         _b_track_dlayer = conversion->trackDLayer(_clusterMap);
0551         _b_track_layer = conversion->firstLayer(_clusterMap);
0552         _b_approach = conversion->approachDistance();
0553         _b_track_dca = conversion->minDca();
0554         //record pT info
0555         _b_track_pT = conversion->minTrackpT();
0556         if(tlv_electron->Pt()>tlv_positron->Pt()) _b_ttrack_pT = tlv_positron->Pt();
0557         else _b_ttrack_pT = tlv_electron->Pt();
0558         //record initial photon info
0559         TLorentzVector* recoPhoton = conversion->getRecoPhoton();
0560         if (recoPhoton)
0561         {
0562             _b_photon_m=recoPhoton->Dot(*recoPhoton);
0563             TLorentzVector truth_added_tlv = *tlv_electron+*tlv_positron;
0564             _b_tphoton_m= truth_added_tlv.Dot(truth_added_tlv);
0565             _b_photon_pT=recoPhoton->Pt();
0566             conversion->PrintPhotonRecoInfo(tlv_photon,tlv_electron,tlv_positron,_b_photon_m);
0567         }
0568         else{//photon was not reconstructed
0569             _b_photon_m =-999;
0570             _b_photon_pT=-999;
0571             conversion->PrintPhotonRecoInfo(tlv_photon,tlv_electron,tlv_positron,_b_photon_m);
0572         }
0573         _b_tphoton_pT=tlv_photon->Pt();
0574         cout<<"second"<<endl;
0575         //truth vertex info
0576         _b_tvtx_radius = sqrt(conversion->getVtx()->get_x()*conversion->getVtx()->get_x()+conversion->getVtx()->get_y()*conversion->getVtx()->get_y());
0577         TVector3 tVertPos(conversion->getVtx()->get_x(),conversion->getVtx()->get_y(),conversion->getVtx()->get_z());
0578         _b_tvtx_phi = tVertPos.Phi();
0579         _b_tvtx_eta = tVertPos.Eta();
0580         _b_tvtx_z = tVertPos.Z();
0581         _b_tvtx_x = tVertPos.X();
0582         _b_tvtx_y = tVertPos.Y();
0583         //TODO check Conversion operations for ownership transfer->memleak due to lack of delete
0584         cout<<"vertexing"<<endl;
0585         genfit::GFRaveVertex* originalVert, *recoVert;
0586         originalVert=recoVert = conversion->getSecondaryVertex(_vertexer);
0587         recoVert = conversion->correctSecondaryVertex(_regressor);
0588         cout<<"finding gf_tracks"<<endl;
0589         //std::pair<PHGenFit::Track*,PHGenFit::Track*> ph_gf_tracks = conversion->getPHGFTracks(_vertexer);
0590         if (recoVert)
0591         {
0592             //cout<<"finding refit_gf_tracks"<<endl;
0593             //std::pair<PHGenFit::Track*,PHGenFit::Track*> refit_phgf_tracks=conversion->refitTracks(_vertexer);
0594             //TODO check repetive refitting and revterexing this is issue #23
0595             /*if (ph_gf_tracks.first&&refit_phgf_tracks.first)
0596                 {
0597                 cout<<"Good Track refit with original:\n";ph_gf_tracks.first->get_mom().Print();cout<<"\n\t and refit:\n";
0598                 refit_phgf_tracks.first->get_mom().Print();
0599                 }
0600                 if (ph_gf_tracks.second&&refit_phgf_tracks.second)
0601                 {
0602                 cout<<"Good Track refit with original:\n";
0603                 ph_gf_tracks.second->get_mom().Print();
0604                 cout<<"\n\t and refit:\n";
0605                 refit_phgf_tracks.second->get_mom().Print();
0606                 }*/
0607             recoPhoton = conversion->getRefitRecoPhoton();
0608             if(recoPhoton) _b_photon_pT=recoPhoton->Pt();
0609             //fill the vtxing tree
0610             TVector3 recoVertPos = originalVert->getPos();
0611             _b_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
0612             _b_vtx_phi = recoVertPos.Phi();
0613             _b_vtx_eta = recoVertPos.Eta();
0614             _b_vtx_z = recoVertPos.Z();
0615             _b_vtx_x = recoVertPos.X();
0616             _b_vtx_y = recoVertPos.Y();
0617             _b_vtx_chi2 = recoVert->getChi2();
0618             //track info
0619             pair<float,float> pTstemp = conversion->getTrackpTs();
0620             _b_track1_pt = pTstemp.first;
0621             _b_track2_pt = pTstemp.second;
0622             pair<float,float> etasTemp = conversion->getTrackEtas();
0623             _b_track1_eta = etasTemp.first;
0624             _b_track2_eta = etasTemp.second;
0625             pair<float,float> phisTemp = conversion->getTrackPhis();
0626             _b_track1_phi = phisTemp.first;
0627             _b_track2_phi = phisTemp.second;
0628             _vtxingTree->Fill();
0629
0630             //populate the refit vertex values
0631             recoVertPos = recoVert->getPos();
0632             _b_vtx_radius = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
0633             _b_vtx_phi = recoVertPos.Phi();
0634             _b_vtx_eta = recoVertPos.Eta();
0635             _b_vtx_z = recoVertPos.Z();
0636             _b_vtx_x = recoVertPos.X();
0637             _b_vtx_y = recoVertPos.Y();
0638             _b_vtxTrackRZ_dist = conversion->vtxTrackRZ(recoVertPos);
0639             _b_vtxTrackRPhi_dist = conversion->vtxTrackRPhi(recoVertPos);
0640             cout<<"done full vtx values"<<endl;
0641         }
0642
0643         else{//vtx not reconstructed
0644             cout<<"no vtx "<<endl;
0645             _b_vtx_radius  =-999.;
0646             _b_vtx_phi = -999.;
0647             _b_vtx_eta = -999.;
0648             _b_vtx_z = -999.;
0649             _b_vtx_x = -999.;
0650             _b_vtx_y = -999.;
0651             _b_vtx_chi2 = -999.;
0652             _b_vtxTrackRZ_dist = -999.;
0653             _b_vtxTrackRPhi_dist = -999.;
0654             //track info
0655             pair<float,float> pTstemp = conversion->getTrackpTs();
0656             _b_track1_pt = pTstemp.first;
0657             _b_track2_pt = pTstemp.second;
0658             pair<float,float> etasTemp = conversion->getTrackEtas();
0659             _b_track1_eta = etasTemp.first;
0660             _b_track2_eta = etasTemp.second;
0661             pair<float,float> phisTemp = conversion->getTrackPhis();
0662             _b_track1_phi = phisTemp.first;
0663             _b_track2_phi = phisTemp.second;
0664             cout<<"done no vtx values"<<endl;
0665             _vtxingTree->Fill();
0666         }
0667         //reset the values
0668         _b_cluster_prob=-999;
0669         _b_cluster_deta=-999;
0670         _b_cluster_dphi=-999;
0671         _b_nCluster=-999;
0672         pair<int,int> clusterIds = conversion->get_cluster_ids();
0673         RawCluster *clustemp;
0674         if(_mainClusterContainer->getCluster(clusterIds.first)){//if there is matching cluster
0675             clustemp =   dynamic_cast<RawCluster*>(_mainClusterContainer->getCluster(clusterIds.first)->CloneMe());
0676             //this is for cluster subtraction which will not be implented soon
0677             // _conversionClusters.AddCluster(clustemp); //add the calo cluster to the container
0678             if (_kMakeTTree)
0679             {
0680                 _b_cluster_prob=clustemp->get_prob();
0681                 RawCluster *clus2 = _mainClusterContainer->getCluster(clusterIds.second);
0682                 if (clus2)
0683                 {
0684                     _b_cluster_dphi=fabs(clustemp->get_phi()-clus2->get_phi());
0685                     TVector3 etaCalc(clustemp->get_x(),clustemp->get_y(),clustemp->get_z());
0686                     //TODO check cluster_prob distribution for signal
0687                     if (clus2->get_prob()>_b_cluster_prob)
0688                     {
0689                         _b_cluster_prob=clus2->get_prob();
0690                     }
0691                     //calculate deta
0692                     float eta1 = etaCalc.PseudoRapidity();
0693                     etaCalc.SetXYZ(clus2->get_x(),clus2->get_y(),clus2->get_z());
0694                     _b_cluster_deta=fabs(eta1-etaCalc.PseudoRapidity());
0695                     if (clusterIds.first!=clusterIds.second) //if there are two district clusters
0696                     {
0697                         _b_nCluster=2;
0698                     }
0699                     else{
0700                         _b_nCluster=1;
0701                     }
0702                 }
0703             }
0704         }
0705     }
0706
0707     else{//not a "complete" conversion some of the data is missing
0708         cout<<"incomplete conversion"<<endl;
0709         _b_track_deta   = -999;
0710         _b_track_dphi   = -999;
0711         _b_track_dlayer = -999;
0712         _b_track_layer  = -999;
0713         _b_approach     = -999;
0714         _b_track_dca    = -999;
0715         //record pT info
0716         _b_track_pT = -999;
0717         if(!tlv_electron||!tlv_positron) _b_ttrack_pT = -999;
0718         else if(tlv_electron->Pt()>tlv_positron->Pt()) _b_ttrack_pT = tlv_positron->Pt();
0719         else _b_ttrack_pT = tlv_electron->Pt();
0720         //record initial photon info
0721         TLorentzVector* recoPhoton = conversion->getRecoPhoton();
0722         if (recoPhoton)
0723         {
0724             _b_photon_pT=recoPhoton->Pt();
0725             _b_photon_m=recoPhoton->Dot(*recoPhoton);
0726             if (tlv_positron&&tlv_electron)
0727             {
0728                 TLorentzVector truth_added_tlv = *tlv_electron+*tlv_positron;
0729                 _b_tphoton_m= truth_added_tlv.Dot(truth_added_tlv);
0730             }
0731             else{
0732                 _b_tphoton_m=-999;
0733             }
0734         }
0735         else{//photon was not reconstructed
0736             _b_photon_m =-999;
0737             _b_photon_pT=-999;
0738         }
0739         if(tlv_photon)_b_tphoton_pT=tlv_photon->Pt();
0740         else _b_tphoton_pT=-999;
0741         cout<<"here"<<endl;
0742         //truth vertex info
0743         _b_tvtx_radius = sqrt(conversion->getVtx()->get_x()*conversion->getVtx()->get_x()+conversion->getVtx()->get_y()*conversion->getVtx()->get_y());
0744         TVector3 tVertPos(conversion->getVtx()->get_x(),conversion->getVtx()->get_y(),conversion->getVtx()->get_z());
0745         cout<<"still here"<<endl;
0746         _b_tvtx_phi = tVertPos.Phi();
0747         _b_tvtx_eta = tVertPos.Eta();
0748         _b_tvtx_z = tVertPos.Z();
0749         _b_tvtx_x = tVertPos.X();
0750         _b_tvtx_y = tVertPos.Y();
0751         //TODO check Conversion operations for ownership transfer->memleak due to lack of delete
0752         //vtx not reconstructed because conversion is incomplete
0753         _b_vtx_radius = -999.;
0754         _b_vtx_phi = -999.;
0755         _b_vtx_eta = -999.;
0756         _b_vtx_z = -999.;
0757         _b_vtx_x = -999.;
0758         _b_vtx_y = -999.;
0759         _b_vtx_chi2 = -999.;
0760         _b_vtxTrackRZ_dist = -999.;
0761         _b_vtxTrackRPhi_dist = -999.;
0762         //track info not possible to reconstruct both tracks
0763         _b_track1_pt  = -999;
0764         _b_track2_pt  = -999;
0765         _b_track1_eta = -999;
0766         _b_track2_eta = -999;
0767         _b_track1_phi = -999;
0768         _b_track2_phi = -999;
0769
0770         //reset the values
0771         _b_cluster_prob=-1;
0772         _b_cluster_deta=-1;
0773         _b_cluster_dphi=-1;
0774         _b_nCluster=-1;
0775         pair<int,int> clusterIds = conversion->get_cluster_ids();
0776         RawCluster *clustemp;
0777         if(_mainClusterContainer->getCluster(clusterIds.first)){//if there is matching cluster
0778             clustemp =   dynamic_cast<RawCluster*>(_mainClusterContainer->getCluster(clusterIds.first)->CloneMe());
0779             //this is for cluster subtraction which will not be implented soon
0780
0781             // _conversionClusters.AddCluster(clustemp); //add the calo cluster to the container
0782             if (_kMakeTTree)
0783             {
0784                 _b_cluster_prob=clustemp->get_prob();
0785                 RawCluster *clus2 = _mainClusterContainer->getCluster(clusterIds.second);
0786                 if (clus2)
0787                 {
0788                     _b_cluster_dphi=fabs(clustemp->get_phi()-clus2->get_phi());
0789                     TVector3 etaCalc(clustemp->get_x(),clustemp->get_y(),clustemp->get_z());
0790                     //TODO check cluster_prob distribution for signal
0791                     if (clus2->get_prob()>_b_cluster_prob)
0792                     {
0793                         _b_cluster_prob=clus2->get_prob();
0794                     }
0795                     //calculate deta
0796                     float eta1 = etaCalc.PseudoRapidity();
0797                     etaCalc.SetXYZ(clus2->get_x(),clus2->get_y(),clus2->get_z());
0798                     _b_cluster_deta=fabs(eta1-etaCalc.PseudoRapidity());
0799                     if (clusterIds.first!=clusterIds.second) //if there are two district clusters
0800                     {
0801                         _b_nCluster=2;
0802                     }
0803                     else{
0804                         _b_nCluster=1;
0805                     }
0806                 }
0807             }
0808         }
0809     }
0810     cout<<"Filling"<<endl;
0811     _signalCutTree->Fill();
0812 }
0813
0814 const RawClusterContainer* TruthConversionEval::getClusters()const {return &_conversionClusters;}
0815
0816 int TruthConversionEval::get_embed(PHG4Particle* particle, PHG4TruthInfoContainer* truthinfo)const{
0817     return truthinfo->isEmbeded(particle->get_track_id());
0818 }
0819
0820 float TruthConversionEval::vtoR(PHG4VtxPoint* vtx)const{
0821     return (float) sqrt(vtx->get_x()*vtx->get_x()+vtx->get_y()*vtx->get_y());
0822 }
0823
0824 int TruthConversionEval::End(PHCompositeNode *topNode)
0825 {
0826     cout<<"ending"<<endl;
0827     if(_kMakeTTree){
0828         cout<<"closing"<<endl;
0829         //_signalCutTree->Write();
0830         _f->Write();
0831         _f->Close();
0832     }
0833     return 0;
0834 }