Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:15:03 

0001 #include "Conversion.h"
0002 #include "SVReco.h"
0003 #include "VtxRegressor.h"
0004 #include <phool/PHCompositeNode.h>
0005 #include <phool/getClass.h>
0006 #include <phgenfit/Track.h>
0007 #include <g4main/PHG4TruthInfoContainer.h>
0008 #include <trackbase_historic/SvtxCluster.h>
0009 #include <trackbase_historic/SvtxHitMap.h>
0010 #include <trackbase_historic/SvtxTrackMap.h>
0011 #include <trackbase_historic/SvtxVertex_v1.h>
0012 #include <trackbase/TrkrClusterContainer.h>
0013 #include <trackbase/TrkrClusterv1.h>
0014 #include <TRandom3.h>
0015 #include <assert.h>
0016 
0017 
0018 const float Conversion::_kElectronRestM=.0005109989461;
0019 
0020 Conversion::Conversion(SvtxTrackEval* trackeval,int verbosity){
0021   this->trackeval=trackeval;
0022   this->verbosity=verbosity;
0023   _refit_phgf_tracks.first=NULL;
0024   _refit_phgf_tracks.second=NULL;
0025   pairTruthReco1.second=0; 
0026   pairTruthReco2.second=0; 
0027   pairTruthReco1.first=0; 
0028   pairTruthReco2.first=0; 
0029 }
0030 Conversion::Conversion(PHG4VtxPoint* vtx,int verbosity){
0031   this->vtx=vtx;
0032   this->verbosity=verbosity;
0033   _refit_phgf_tracks.first=NULL;
0034   _refit_phgf_tracks.second=NULL;
0035   pairTruthReco1.second=0; 
0036   pairTruthReco2.second=0; 
0037   pairTruthReco1.first=0; 
0038   pairTruthReco2.first=0; 
0039 }
0040 
0041 Conversion::Conversion(PHG4VtxPoint* vtx,SvtxTrackEval *trackeval,int verbosity){
0042   this->trackeval=trackeval;
0043   this->vtx=vtx;
0044   this->verbosity=verbosity;
0045   _refit_phgf_tracks.first=NULL;
0046   _refit_phgf_tracks.second=NULL;
0047   pairTruthReco1.second=0; 
0048   pairTruthReco2.second=0; 
0049   pairTruthReco1.first=0; 
0050   pairTruthReco2.first=0; 
0051 }
0052 
0053 Conversion::~Conversion(){
0054   if(recoVertex) delete recoVertex;
0055   if(recoPhoton) delete recoPhoton;
0056   if(truthSvtxVtx) delete truthSvtxVtx;
0057   if (_refit_phgf_tracks.first) delete _refit_phgf_tracks.first;
0058   if (_refit_phgf_tracks.second) delete _refit_phgf_tracks.second;
0059   _refit_phgf_tracks.first=NULL;
0060   _refit_phgf_tracks.second=NULL;
0061   recoVertex=NULL;
0062   recoPhoton=NULL;
0063   truthSvtxVtx=NULL;
0064   //dont delete the points as you are not the owner and did not make your own copies
0065 }
0066 
0067 void Conversion::setElectron(PHG4Particle* e){
0068   if (e1)
0069   {
0070     if (e2&&verbosity>0)
0071     {
0072       std::cout<<"WARNING in Conversion oversetting conversion electrons"<<std::endl;
0073     }
0074     else{
0075       e2=e;
0076       pairTruthReco2.second=e->get_track_id();
0077     }
0078   }
0079   else{
0080     e1=e;
0081     pairTruthReco1.second=e->get_track_id();
0082   }
0083 }
0084 
0085 
0086 bool Conversion::setElectron(){
0087   if (hasPair())
0088   {
0089     if (e1->get_pid()<0)
0090     {
0091       PHG4Particle *temp = e1;
0092       e1=e2;
0093       e2=temp;
0094     }
0095     if (e1->get_pid()<0)
0096     {
0097       if (verbosity>0)
0098       {
0099         std::cout<<"Warning in Conversion bad charges in conversion"<<std::endl;
0100       }
0101       return false;
0102     }
0103   }
0104   else{
0105     if (!e1)
0106     {
0107       if (verbosity>0)
0108       {
0109         std::cout<<"Warning in Conversion no truth tracks set"<<std::endl;
0110       }
0111       return false;
0112     }
0113     else if (e1->get_pid()<0) return false;
0114   }
0115   return true;
0116 }
0117 
0118 PHG4Particle* Conversion::getElectron(){
0119   setElectron();
0120   return e1;
0121 }
0122 
0123 PHG4Particle* Conversion::getPositron(){
0124   if(setElectron()){
0125     return e2;
0126   }
0127   else{
0128     return e1;
0129   }
0130 }
0131 
0132 bool Conversion::setParent(PHG4Particle* parent){
0133   bool r =true;
0134   if(!photon) photon=parent;
0135   else{
0136     if(!(*photon==*parent)) cerr<<"Bad photon matching!"<<endl;
0137     r=false;
0138   }
0139   return r;
0140 }
0141 
0142 void Conversion::setPrimaryPhoton(PHG4Particle* parent,PHG4TruthInfoContainer* truthinfo){
0143   if(!setParent(parent)) cerr<<"Bad photon matching during primary photon set"<<endl;
0144   if(photon->get_track_id()==parent->get_primary_id()){
0145     primaryPhoton=photon;
0146   }
0147   else{
0148     primaryPhoton=truthinfo->GetParticle(parent->get_primary_id());
0149   }
0150 }
0151 
0152 void Conversion::setRecoTrack(int truthID, SvtxTrack* recoTrack){
0153   if(!recoTrack)return;
0154   setElectron();
0155   setRecoTracks();
0156   if (e1&&e1->get_track_id()==truthID&&!reco1)
0157   {
0158     reco1=recoTrack;
0159     pairTruthReco1.first=e1->get_track_id();
0160     pairTruthReco1.second=recoTrack->get_id();
0161   }
0162   else if (e2&&e2->get_track_id()==truthID&&!reco2)
0163   {
0164     reco2=recoTrack;
0165     pairTruthReco2.first=e2->get_track_id();
0166     pairTruthReco2.second=recoTrack->get_id();
0167   }
0168 }
0169 
0170 int Conversion::setRecoTracks(SvtxTrackEval* trackeval){    
0171   this->trackeval=trackeval;
0172   setElectron();
0173   if (e1)
0174   {
0175     reco1=trackeval->best_track_from(e1);  
0176   }
0177   if (e2)
0178   {
0179     reco2=trackeval->best_track_from(e2);
0180   }
0181   if(reco2==reco1){
0182     reco2=NULL;
0183   }
0184   int r=0;
0185   if (reco1)
0186   {
0187     r++;
0188     pairTruthReco1.first = e1->get_track_id();
0189     pairTruthReco1.second = reco1->get_id();
0190   }
0191   if (reco2)
0192   {
0193     r++;
0194     pairTruthReco2.first = e2->get_track_id();
0195     pairTruthReco2.second = reco2->get_id();
0196   }
0197   setRecoPhoton();
0198   return r;
0199 }
0200 
0201 int Conversion::setRecoTracks(){  
0202   assert(trackeval);
0203   if (e1)
0204   {
0205     reco1=trackeval->best_track_from(e1); // have not checked that these are in range 
0206   }
0207   if (e2)
0208   {
0209     reco2=trackeval->best_track_from(e2);
0210   }
0211   if(reco2==reco1){
0212     reco2=NULL;
0213   }
0214   int r=0;
0215   if (reco1)
0216   {
0217     r++;
0218     pairTruthReco1.first = e1->get_track_id();
0219     pairTruthReco1.second = reco1->get_id();
0220   }
0221   if (reco2)
0222   {
0223     r++;
0224     pairTruthReco2.first = e2->get_track_id();
0225     pairTruthReco2.second = reco2->get_id();
0226   }
0227   setRecoPhoton();
0228   return r;
0229 }
0230 
0231 SvtxTrack* Conversion::getRecoTrack(unsigned truthID) const{
0232   if (pairTruthReco1.second==truthID)
0233   {
0234     if(reco1&&reco1->get_id()==pairTruthReco1.first)
0235       return reco1;
0236     else if(reco2&&reco2->get_id()==pairTruthReco1.first)
0237       return reco2;
0238   }
0239   else if (pairTruthReco2.second==truthID)
0240   {
0241     if(reco1&&reco1->get_id()==pairTruthReco2.first)
0242       return reco1;
0243     else if(reco2&&reco2->get_id()==pairTruthReco2.first)
0244       return reco2;
0245   }
0246   return NULL;
0247 }
0248 
0249 TLorentzVector* Conversion::setRecoPhoton(){
0250   if (reco1&&reco2)
0251   {
0252     TLorentzVector tlv1(reco1->get_px(),reco1->get_py(),reco1->get_pz(),
0253         sqrt(_kElectronRestM*_kElectronRestM+reco1->get_p()*reco1->get_p()));
0254     TLorentzVector tlv2(reco2->get_px(),reco2->get_py(),reco2->get_pz(),
0255         sqrt(_kElectronRestM*_kElectronRestM+reco2->get_p()*reco2->get_p()));
0256     if (recoPhoton) delete recoPhoton;
0257     recoPhoton= new TLorentzVector(tlv1+tlv2);
0258   }
0259   return recoPhoton;
0260 }
0261 
0262 TLorentzVector* Conversion::getRecoPhoton(){
0263   return setRecoPhoton();
0264 }
0265 
0266 TLorentzVector* Conversion::getRecoPhoton(SvtxTrack* reco1, SvtxTrack* reco2){
0267   if(!(reco1&&reco2)) return NULL;
0268   TLorentzVector tlv1(reco1->get_px(),reco1->get_py(),reco1->get_pz(),
0269       sqrt(_kElectronRestM*_kElectronRestM+reco1->get_p()*reco1->get_p()));
0270   TLorentzVector tlv2(reco2->get_px(),reco2->get_py(),reco2->get_pz(),
0271       sqrt(_kElectronRestM*_kElectronRestM+reco2->get_p()*reco2->get_p()));
0272   return new TLorentzVector(tlv1+tlv2);
0273 }
0274 
0275 TLorentzVector* Conversion::getRefitRecoPhoton(){
0276   std::pair<TLorentzVector*,TLorentzVector*> refit_tlvs =getRefitRecoTlvs();
0277   if (refit_tlvs.first&&refit_tlvs.second)
0278   {
0279     if(recoPhoton) delete recoPhoton;
0280     recoPhoton= new TLorentzVector(*refit_tlvs.first + *refit_tlvs.second);
0281   }
0282   return recoPhoton;
0283 }
0284 
0285 std::pair<TLorentzVector*,TLorentzVector*> Conversion::getRecoTlvs(){
0286   std::pair<TLorentzVector*,TLorentzVector*> r;
0287   switch(recoCount()){
0288     case 2:
0289       r.first = new TLorentzVector();
0290       r.first->SetPxPyPzE (reco1->get_px(),reco1->get_py(),reco1->get_pz(),
0291           sqrt(_kElectronRestM*_kElectronRestM+reco1->get_p()*reco1->get_p()));
0292       r.second =   new TLorentzVector();
0293       r.second->SetPxPyPzE (reco2->get_px(),reco2->get_py(),reco2->get_pz(),
0294           sqrt(_kElectronRestM*_kElectronRestM+reco2->get_p()*reco2->get_p()));
0295       break;
0296     case 1:
0297       if(reco1){
0298         r.first = new TLorentzVector();
0299         r.first->SetPxPyPzE (reco1->get_px(),reco1->get_py(),reco1->get_pz(),
0300             sqrt(_kElectronRestM*_kElectronRestM+reco1->get_p()*reco1->get_p()));
0301         r.second = NULL;
0302       }
0303       else{
0304         r.first = NULL;
0305         r.second = new TLorentzVector ();
0306         r.second->SetPxPyPzE(reco2->get_px(),reco2->get_py(),reco2->get_pz(),
0307             sqrt(_kElectronRestM*_kElectronRestM+reco2->get_p()*reco2->get_p()));
0308       }
0309       break;
0310     default:
0311       r.first=NULL;
0312       r.second=NULL;
0313       break;
0314   }
0315   return r;
0316 }
0317 
0318 std::pair<TLorentzVector*,TLorentzVector*> Conversion::getRefitRecoTlvs(){
0319   std::pair<TLorentzVector*,TLorentzVector*> r;
0320   if (_refit_phgf_tracks.first&&_refit_phgf_tracks.second)
0321   {
0322     r.first = new TLorentzVector();
0323     r.first->SetVectM (_refit_phgf_tracks.first->get_mom(),_kElectronRestM);
0324     r.second =   new TLorentzVector();
0325     r.second->SetVectM (_refit_phgf_tracks.second->get_mom(),_kElectronRestM);
0326   }
0327   else{
0328     r.first=NULL;
0329     r.second=NULL;
0330   }  
0331   return r;
0332 }
0333 
0334 PHG4Particle* Conversion::getTruthPhoton(PHG4TruthInfoContainer* truthinfo){
0335   return photon;
0336 }
0337 
0338 
0339 int Conversion::get_cluster_id(){
0340   if(!reco1){
0341     assert(trackeval);
0342     reco1=trackeval->best_track_from(e1);
0343     if(!reco1){
0344       return -1;
0345     }
0346   }
0347   return reco1->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1));//id of the emcal
0348 }
0349 
0350 int Conversion::get_cluster_id()const{
0351   if(!reco1){
0352     return -1;
0353   }
0354   else return reco1->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1));//id of the emcal
0355 }
0356 
0357 int Conversion::get_cluster_id(SvtxTrackEval *trackeval){
0358   this->trackeval=trackeval;
0359   if (!reco1)
0360   {
0361     reco1=trackeval->best_track_from(e1);
0362     if(!reco1){
0363       cout<<"bad reco"<<endl;
0364       return -1;
0365     }
0366   }
0367   return reco1->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1));//id of the emcal
0368 }
0369 
0370 std::pair<int,int> Conversion::get_cluster_ids(){
0371   switch(recoCount()){
0372     case 2:
0373       return std::pair<int,int>(reco1->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)),reco2->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)));
0374       break;
0375     case 1:
0376       if (reco1)
0377       {
0378         return std::pair<int,int>(reco1->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)),-1);
0379       }
0380       else return std::pair<int,int>(reco2->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)),-1);
0381       break;
0382     default:
0383       return std::pair<int,int>(-1,-1);
0384       break;
0385   }
0386 }
0387 
0388 
0389 
0390 bool Conversion::hasSilicon(SvtxClusterMap* svtxClusterMap){
0391   switch(recoCount()){
0392     case 2:
0393       {
0394         SvtxCluster *c1 = svtxClusterMap->get(*(reco1->begin_clusters()));
0395         SvtxCluster *c2 = svtxClusterMap->get(*(reco2->begin_clusters()));
0396         return c1->get_layer()<=_kNSiliconLayer||c2->get_layer()<=_kNSiliconLayer;
0397       }
0398       break;
0399     case 1:
0400       {
0401         SvtxCluster *c1;
0402         if (reco1)
0403         {
0404           c1 = svtxClusterMap->get(*(reco1->begin_clusters()));
0405         }
0406         else{
0407           c1 = svtxClusterMap->get(*(reco2->begin_clusters()));
0408         }
0409         return c1->get_layer()<=_kNSiliconLayer;
0410       }
0411       break;
0412     default:
0413       return false;
0414       break;
0415   }
0416 }
0417 
0418 float Conversion::vtxTrackRZ(TVector3 vertpos){
0419   float d1=vtxTrackRZ(vertpos,reco1);
0420   float d2=vtxTrackRZ(vertpos,reco2);
0421   return d1>d2?d1:d2;
0422 }
0423 float Conversion::vtxTrackRZ(TVector3 vertpos,SvtxTrack* reco1,SvtxTrack* reco2){
0424   float d1=vtxTrackRZ(vertpos,reco1);
0425   float d2=vtxTrackRZ(vertpos,reco2);
0426   return d1>d2?d1:d2;
0427 }
0428 
0429 float Conversion::vtxTrackRPhi(TVector3 vertpos){
0430   float d1=vtxTrackRPhi(vertpos,reco1);
0431   float d2=vtxTrackRPhi(vertpos,reco2);
0432   return d1>d2?d1:d2;
0433 }
0434 
0435 float Conversion::vtxTrackRPhi(TVector3 vertpos,SvtxTrack* reco1,SvtxTrack* reco2){
0436   float d1=vtxTrackRPhi(vertpos,reco1);
0437   float d2=vtxTrackRPhi(vertpos,reco2);
0438   return d1>d2?d1:d2;
0439 }
0440 float Conversion::vtxTrackRZ(TVector3 recoVertPos,SvtxTrack *track){
0441   float dR = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y())-sqrt(track->get_x()*track->get_x()+track->get_y()*track->get_y());
0442   float dZ = recoVertPos.z()-track->get_z();
0443   return sqrt(dR*dR+dZ*dZ);
0444 }
0445 
0446 float Conversion::vtxTrackRPhi(TVector3 recoVertPos,SvtxTrack *track){
0447   float vtxR=sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
0448   float trackR=sqrt(track->get_x()*track->get_x()+track->get_y()*track->get_y());
0449   return sqrt(vtxR*vtxR+trackR*trackR-2*vtxR*trackR*cos(recoVertPos.Phi()-track->get_phi()));
0450 }
0451 
0452 int Conversion::trackDLayer(SvtxClusterMap* svtxClusterMap,SvtxHitMap *hitmap){
0453   if (recoCount()==2){
0454     SvtxCluster *c1 = svtxClusterMap->get(*(reco1->begin_clusters()));
0455     SvtxCluster *c2 = svtxClusterMap->get(*(reco2->begin_clusters()));
0456     SvtxHit *h1 = hitmap->get(*(c1->begin_hits()));
0457     SvtxHit *h2 = hitmap->get(*(c2->begin_hits()));
0458     int l1 = h1->get_layer();
0459     int l2 = h2->get_layer();
0460     return abs(l1-l2);
0461   }
0462   else return -1;
0463 }
0464 
0465 float Conversion::minDca(){
0466   if(reco1->get_dca()>reco2->get_dca()) return reco2->get_dca();
0467   else return reco1->get_dca();
0468 }
0469 
0470 int Conversion::trackDLayer(TrkrClusterContainer* clusterMap){
0471   if (recoCount()==2){
0472     TrkrCluster *c1 = clusterMap->findCluster(*(reco1->begin_cluster_keys()));
0473     TrkrCluster *c2 = clusterMap->findCluster(*(reco2->begin_cluster_keys()));
0474     unsigned l1 = TrkrDefs::getLayer(c1->getClusKey());
0475     unsigned l2 = TrkrDefs::getLayer(c2->getClusKey());
0476     return abs(l1-l2);
0477   }
0478   else return -1;
0479 }
0480 int Conversion::trackDLayer(TrkrClusterContainer* clusterMap,SvtxTrack* reco1, SvtxTrack* reco2){
0481   TrkrCluster *c1 = clusterMap->findCluster(*(reco1->begin_cluster_keys()));
0482   TrkrCluster *c2 = clusterMap->findCluster(*(reco2->begin_cluster_keys()));
0483   unsigned l1 = TrkrDefs::getLayer(c1->getClusKey());
0484   unsigned l2 = TrkrDefs::getLayer(c2->getClusKey());
0485   return abs(l1-l2);
0486 }
0487 
0488 int Conversion::firstLayer(SvtxClusterMap* svtxClusterMap,SvtxHitMap *hitmap){
0489   switch(recoCount()){
0490     case 2:
0491       {
0492         SvtxCluster *c1 = svtxClusterMap->get(*(reco1->begin_clusters()));
0493         SvtxCluster *c2 = svtxClusterMap->get(*(reco2->begin_clusters()));
0494         SvtxHit *h1 = hitmap->get(*(c1->begin_hits()));
0495         SvtxHit *h2 = hitmap->get(*(c2->begin_hits()));
0496         if(h1->get_layer()>h2->get_layer()){
0497           return h2->get_layer();
0498         }
0499         else return h1->get_layer();
0500       }
0501     case 1:
0502       {
0503         if (reco1)return svtxClusterMap->get(*(reco1->begin_clusters()))->get_layer();
0504         else return svtxClusterMap->get(*(reco2->begin_clusters()))->get_layer();
0505       }
0506     default:
0507       return -1;
0508   }
0509 }
0510 
0511 int Conversion::firstLayer(TrkrClusterContainer* clusterMap){
0512   switch(recoCount()){
0513     case 2:
0514       {
0515         TrkrCluster *c1 = clusterMap->findCluster(*(reco1->begin_cluster_keys()));
0516         TrkrCluster *c2 = clusterMap->findCluster(*(reco2->begin_cluster_keys()));
0517         unsigned l1 = TrkrDefs::getLayer(c1->getClusKey());
0518         unsigned l2 = TrkrDefs::getLayer(c2->getClusKey());
0519         //maybe i need this? TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluskey);
0520         if(l1>l2){
0521           return l1;
0522         }
0523         else return l2;
0524       }
0525     case 1:
0526       {
0527         if (reco1)return TrkrDefs::getLayer(*(reco1->begin_cluster_keys()));
0528         else return TrkrDefs::getLayer(*(reco2->begin_cluster_keys()));
0529       }
0530     default:
0531       return -1;
0532   }
0533 }
0534 
0535 double Conversion::dist(PHG4VtxPoint *recovtx,SvtxClusterMap* svtxClusterMap){
0536   SvtxCluster *c1 = svtxClusterMap->get(*(reco1->begin_clusters()));
0537   SvtxCluster *c2 = svtxClusterMap->get(*(reco2->begin_clusters()));
0538   double r1 = sqrt(abs(c1->get_x()-recovtx->get_x())+abs(c1->get_y()-recovtx->get_y())+abs(c1->get_z()-recovtx->get_z()));
0539   double r2 = sqrt(abs(c2->get_x()-recovtx->get_x())+abs(c2->get_y()-recovtx->get_y())+abs(c2->get_z()-recovtx->get_z()));
0540   if (r1>r2)
0541   {
0542     return r1;
0543   }
0544   else return r2;
0545 }
0546 
0547 double Conversion::dist(TVector3 *recovtx,SvtxClusterMap* svtxClusterMap){
0548   SvtxCluster *c1 = svtxClusterMap->get(*(reco1->begin_clusters()));
0549   SvtxCluster *c2 = svtxClusterMap->get(*(reco2->begin_clusters()));
0550   double r1 = sqrt(abs(c1->get_x()-recovtx->x())+abs(c1->get_y()-recovtx->y())+abs(c1->get_z()-recovtx->z()));
0551   double r2 = sqrt(abs(c2->get_x()-recovtx->x())+abs(c2->get_y()-recovtx->y())+abs(c2->get_z()-recovtx->z()));
0552   if (r1>r2)
0553   {
0554     return r1;
0555   }
0556   else return r2;
0557 }
0558 
0559 double Conversion::dist(TVector3 *recovtx,TrkrClusterContainer* clusterMap){
0560   TrkrCluster *c1 = clusterMap->findCluster(*(reco1->begin_cluster_keys()));
0561   TrkrCluster *c2 = clusterMap->findCluster(*(reco2->begin_cluster_keys()));
0562   double r1 = sqrt(abs(c1->getX()-recovtx->x())+abs(c1->getY()-recovtx->y())+abs(c1->getZ()-recovtx->z()));
0563   double r2 = sqrt(abs(c2->getX()-recovtx->x())+abs(c2->getY()-recovtx->y())+abs(c2->getZ()-recovtx->z()));
0564   if (r1>r2)
0565   {
0566     return r1;
0567   }
0568   else return r2;
0569 }
0570 
0571 void Conversion::printTruth(){
0572   cout<<"Conversion with truth info:\n";
0573   if (e1) e1->identify();
0574   if (e2) e2->identify();
0575   if (vtx) recoVertex->Print();
0576   if (recoPhoton) recoPhoton->Print();
0577 }
0578 void Conversion::printReco(){
0579   cout<<"Conversion with reco info:\n";
0580   if (reco1) reco1->identify();
0581   if (reco2) reco2->identify();
0582   if (vtx) vtx->identify();
0583   if (photon) photon->identify();
0584 }
0585 
0586 void Conversion::PrintPhotonRecoInfo(){
0587   if(!recoPhoton) cerr<<"No photon reconstructed"<<endl;
0588   else{
0589     cout<<"Truth Track: ";e1->identify();
0590     cout<<"Reco Track: ";getRecoTrack(e1->get_track_id())->identify(); 
0591     cout<<"Truth Track: ";e2->identify();
0592     cout<<"Reco Track: ";getRecoTrack(e2->get_track_id())->identify(); 
0593     cout<<"Truth Photon: ";photon->identify();
0594     cout<<"Reco Photon: ";recoPhoton->Print();
0595   }
0596 }
0597 void Conversion::PrintPhotonRecoInfo(TLorentzVector *tlv_photon,TLorentzVector *tlv_electron, TLorentzVector *tlv_positron,float mass){
0598   if(!recoPhoton) cerr<<"No photon reconstructed"<<endl;
0599   else{
0600     cout<<"Truth Track: ";tlv_electron->Print();
0601     cout<<"Reco Track: ";reco1->identify(); 
0602     cout<<"Truth Track: ";tlv_positron->Print();
0603     cout<<"Reco Track: ";reco2->identify(); 
0604     cout<<"Truth Photon: ";tlv_photon->Print();
0605     cout<<"Reco Photon with mass: "<<mass<<": ";recoPhoton->Print();
0606   }
0607 }
0608 
0609 /*This is deprecated
0610  * float Conversion::setRecoVtx(SvtxVertex *recovtx,SvtxClusterMap* svtxClusterMap){
0611  recoVertex=recovtx;
0612  SvtxCluster *c1 = svtxClusterMap->get(*(reco1->begin_clusters()));
0613  SvtxCluster *c2 = svtxClusterMap->get(*(reco2->begin_clusters()));
0614  float r1 = sqrt(abs(c1->get_x()-recovtx->get_x())+abs(c1->get_y()-recovtx->get_y())+abs(c1->get_z()-recovtx->get_z()));
0615  float r2 = sqrt(abs(c2->get_x()-recovtx->get_x())+abs(c2->get_y()-recovtx->get_y())+abs(c2->get_z()-recovtx->get_z()));
0616  if (r1>r2)
0617  {
0618  return r1;
0619  }
0620  else return r2;
0621  }*/
0622 
0623 double Conversion::approachDistance()const{
0624   if (recoCount()==2)
0625   {
0626     static const double eps = 0.000001;
0627     TVector3 u(reco1->get_px(),reco1->get_py(),reco1->get_pz());
0628     TVector3 v(reco2->get_px(),reco2->get_py(),reco2->get_pz());
0629     TVector3 w(reco1->get_x()-reco2->get_x(),reco1->get_x()-reco2->get_y(),reco1->get_x()-reco2->get_z());
0630 
0631     double a = u.Dot(u);
0632     double b = u.Dot(v);
0633     double c = v.Dot(v);
0634     double d = u.Dot(w);
0635     double e = v.Dot(w);
0636 
0637     double D = a*c - b*b;
0638     double sc, tc;
0639     // compute the line parameters of the two closest points
0640     if (D < eps) {         // the lines are almost parallel
0641       sc = 0.0;
0642       tc = (b>c ? d/b : e/c);   // use the largest denominator
0643     }
0644     else {
0645       sc = (b*e - c*d) / D;
0646       tc = (a*e - b*d) / D;
0647     }
0648     // get the difference of the two closest points
0649     u*=sc;
0650     v*=tc;
0651     w+=u;
0652     w-=v;
0653     return w.Mag();   // return the closest distance 
0654   }
0655   else return -1; 
0656 }
0657 
0658 float Conversion::trackDEta()const{
0659   if (recoCount()==2)
0660   {
0661     return fabs(reco1->get_eta()-reco2->get_eta());
0662   }
0663   else return -1.;
0664 }
0665 
0666 double Conversion::approachDistance(SvtxTrack* reco1, SvtxTrack* reco2){
0667   static const double eps = 0.000001;
0668   TVector3 u(reco1->get_px(),reco1->get_py(),reco1->get_pz());
0669   TVector3 v(reco2->get_px(),reco2->get_py(),reco2->get_pz());
0670   TVector3 w(reco1->get_x()-reco2->get_x(),reco1->get_x()-reco2->get_y(),reco1->get_x()-reco2->get_z());
0671 
0672   double a = u.Dot(u);
0673   double b = u.Dot(v);
0674   double c = v.Dot(v);
0675   double d = u.Dot(w);
0676   double e = v.Dot(w);
0677 
0678   double D = a*c - b*b;
0679   double sc, tc;
0680   // compute the line parameters of the two closest points
0681   if (D < eps) {         // the lines are almost parallel
0682     sc = 0.0;
0683     tc = (b>c ? d/b : e/c);   // use the largest denominator
0684   }
0685   else {
0686     sc = (b*e - c*d) / D;
0687     tc = (a*e - b*d) / D;
0688   }
0689   // get the difference of the two closest points
0690   u*=sc;
0691   v*=tc;
0692   w+=u;
0693   w-=v;
0694   return w.Mag();   // return the closest distance 
0695 }
0696 
0697 float Conversion::trackDEta(SvtxTrack* reco1, SvtxTrack* reco2){
0698   return fabs(reco1->get_eta()-reco2->get_eta());
0699 }
0700 
0701 float Conversion::minTrackpT(){
0702   switch(recoCount()){
0703     case 2:
0704       if (reco1->get_pt()<reco2->get_pt())
0705       {
0706         return reco1->get_pt();
0707       }
0708       else return reco2->get_pt();
0709       break;
0710     case 1:
0711       if (reco1)
0712       {
0713         return reco1->get_pt();
0714       }
0715       else return reco2->get_pt();
0716       break;
0717     default:
0718       return -1;
0719       break;
0720   }
0721 }
0722 
0723 std::pair<float,float> Conversion::getTrackpTs(){
0724   switch(recoCount()){
0725     case 2:
0726       return std::pair<float,float>(reco1->get_pt(),reco2->get_pt());
0727     case 1:
0728       if (reco1) return std::pair<float,float>(reco1->get_pt(),-1);
0729       else return std::pair<float,float>(-1,reco2->get_pt());
0730       break;
0731     default:
0732       return std::pair<float,float>(-1,-1);
0733       break;
0734   }
0735 }
0736 
0737 std::pair<float,float> Conversion::getTrackEtas(){
0738   switch(recoCount()){
0739     case 2:
0740       return std::pair<float,float>(reco1->get_eta(),reco2->get_eta());
0741     case 1:
0742       if (reco1) return std::pair<float,float>(reco1->get_eta(),-1);
0743       else return std::pair<float,float>(-1,reco2->get_eta());
0744       break;
0745     default:
0746       return std::pair<float,float>(-1,-1);
0747       break;
0748   }
0749 }
0750 
0751 std::pair<float,float> Conversion::getTrackPhis(){
0752   switch(recoCount()){
0753     case 2:
0754       return std::pair<float,float>(reco1->get_phi(),reco2->get_phi());
0755     case 1:
0756       if (reco1) return std::pair<float,float>(reco1->get_phi(),-1);
0757       else return std::pair<float,float>(-1,reco2->get_phi());
0758       break;
0759     default:
0760       return std::pair<float,float>(-1,-1);
0761       break;
0762   }
0763 }
0764 /*This is the ideal case but I do not have RaveVtx to SvtxVertex matching yet
0765  * genfit::GFRaveVertex* Conversion::getSecondaryVertex(SVReco* vertexer){
0766  if(recoCount()==2){
0767  if (recoVertex) delete recoVertex;
0768  recoVertex= vertexer->findSecondaryVertex(reco1,reco2);
0769  }
0770  return recoVertex;
0771  }*/
0772 genfit::GFRaveVertex* Conversion::getSecondaryVertex(SVReco* vertexer){
0773   if(recoCount()==2){
0774     //this might seg fault
0775     if(recoVertex) delete recoVertex;
0776     recoVertex= vertexer->findSecondaryVertex(reco1,reco2);
0777   }
0778   return recoVertex;
0779 }
0780 
0781 genfit::GFRaveVertex* Conversion::correctSecondaryVertex(VtxRegressor* regressor){
0782   if(!recoVertex) {
0783     cerr<<"WARNING: no vertex to correct"<<endl;
0784     return NULL;
0785   }
0786   if (recoCount()!=2)
0787   {
0788     cerr<<"WARNING: no reco tracks to do vertex correction"<<endl;
0789     return NULL;
0790   }
0791 
0792   TVector3 nextPos = recoVertex->getPos();
0793   nextPos.SetMagThetaPhi(regressor->regress(reco1,reco2,recoVertex),nextPos.Theta(),nextPos.Phi());
0794 
0795   using namespace genfit;
0796   // GFRaveVertex* temp = recoVertex;
0797   std::vector<GFRaveTrackParameters*> tracks;
0798   for(unsigned i =0; i<recoVertex->getNTracks();i++){
0799     tracks.push_back(recoVertex->getParameters(i));
0800   }
0801   recoVertex = new GFRaveVertex(nextPos,recoVertex->getCov(),tracks,recoVertex->getNdf(),recoVertex->getChi2(),recoVertex->getId());
0802   //  delete temp; //this caused outside references to seg fault //TODO shared_ptr is better 
0803   return recoVertex;
0804 }
0805 
0806 genfit::GFRaveVertex* Conversion::correctSecondaryVertex(VtxRegressor* regressor,genfit::GFRaveVertex* recoVertex, SvtxTrack* reco1, SvtxTrack* reco2){
0807   if(!recoVertex) {
0808     cerr<<"WARNING: no vertex to correct"<<endl;
0809     return NULL;
0810   }
0811   if (!(reco1&&reco2))
0812   {
0813     cerr<<"WARNING: no reco tracks to do vertex correction"<<endl;
0814     return NULL;
0815   }
0816 
0817   TVector3 nextPos = recoVertex->getPos();
0818   nextPos.SetMagThetaPhi(regressor->regress(reco1,reco2,recoVertex),nextPos.Theta(),nextPos.Phi());
0819 
0820   using namespace genfit;
0821   // GFRaveVertex* temp = recoVertex;
0822   std::vector<GFRaveTrackParameters*> tracks;
0823   for(unsigned i =0; i<recoVertex->getNTracks();i++){
0824     tracks.push_back(recoVertex->getParameters(i));
0825   }
0826   recoVertex = new GFRaveVertex(nextPos,recoVertex->getCov(),tracks,recoVertex->getNdf(),recoVertex->getChi2(),recoVertex->getId());
0827   //  delete temp; //this caused outside references to seg fault //TODO shared_ptr is better 
0828   return recoVertex;
0829 }
0830 
0831 std::pair<PHGenFit::Track*,PHGenFit::Track*> Conversion::getPHGFTracks(SVReco* vertexer){
0832   std::pair<PHGenFit::Track*,PHGenFit::Track*> r;
0833   r.first = vertexer->getPHGFTrack(reco1);
0834   r.second = vertexer->getPHGFTrack(reco2);
0835   return r;
0836 }
0837 
0838 std::pair<PHGenFit::Track*,PHGenFit::Track*> Conversion::refitTracksTruthVtx(SVReco* vertexer,SvtxVertex* seedVtx){
0839   PHG4VtxPointToSvtxVertex(seedVtx);
0840   if(_refit_phgf_tracks.first) delete _refit_phgf_tracks.first;
0841   if(_refit_phgf_tracks.second) delete _refit_phgf_tracks.second;
0842   _refit_phgf_tracks.first=vertexer->refitTrack(truthSvtxVtx,reco1);
0843   _refit_phgf_tracks.second=vertexer->refitTrack(truthSvtxVtx,reco2);
0844   return _refit_phgf_tracks;
0845 }
0846 
0847 std::pair<PHGenFit::Track*,PHGenFit::Track*> Conversion::refitTracksTruthVtx(SVReco* vertexer){
0848   PHG4VtxPointToSvtxVertex();
0849   if(_refit_phgf_tracks.first) delete _refit_phgf_tracks.first;
0850   if(_refit_phgf_tracks.second) delete _refit_phgf_tracks.second;
0851   _refit_phgf_tracks.first=vertexer->refitTrack(truthSvtxVtx,reco1);
0852   _refit_phgf_tracks.second=vertexer->refitTrack(truthSvtxVtx,reco2);
0853   return _refit_phgf_tracks;
0854 }
0855 std::pair<PHGenFit::Track*,PHGenFit::Track*> Conversion::refitTracks(SVReco* vertexer){
0856   if (!recoVertex) getSecondaryVertex(vertexer);
0857   if(!recoVertex)
0858   {
0859     cerr<<"WARNING: No vertex to refit tracks"<<endl;
0860   }
0861   else{
0862     if(_refit_phgf_tracks.first) delete _refit_phgf_tracks.first;
0863     if(_refit_phgf_tracks.second) delete _refit_phgf_tracks.second;
0864     _refit_phgf_tracks.first=vertexer->refitTrack(recoVertex,reco1);
0865     _refit_phgf_tracks.second=vertexer->refitTrack(recoVertex,reco2);
0866   }
0867   return _refit_phgf_tracks;
0868 }
0869 
0870 std::pair<PHGenFit::Track*,PHGenFit::Track*> Conversion::refitTracks(SVReco* vertexer, SvtxVertex* recoVtx){
0871   if (!recoVtx)
0872   {
0873     cerr<<"WARNING: No vertex to refit tracks"<<endl;
0874   }
0875   else{
0876     if(_refit_phgf_tracks.first) delete _refit_phgf_tracks.first;
0877     if(_refit_phgf_tracks.second) delete _refit_phgf_tracks.second;
0878     _refit_phgf_tracks.first=vertexer->refitTrack(recoVtx,reco1);
0879     _refit_phgf_tracks.second=vertexer->refitTrack(recoVtx,reco2);
0880   }
0881   return _refit_phgf_tracks;
0882 }
0883 
0884 void Conversion::PHG4VtxPointToSvtxVertex(){
0885   truthSvtxVtx = new SvtxVertex_v1();
0886   truthSvtxVtx->set_x(vtx->get_x());
0887   truthSvtxVtx->set_y(vtx->get_y());
0888   truthSvtxVtx->set_z(vtx->get_z());
0889   truthSvtxVtx->set_t0(vtx->get_t());
0890   truthSvtxVtx->set_chisq(1.);
0891   truthSvtxVtx->set_ndof(1);
0892 
0893   TRandom3 rand(0);
0894   double ae = 0.0007;  //7 um
0895   double i = 0.003;//30 um
0896   double d = rand.Gaus(0, ae);
0897   double g = rand.Gaus(0, ae);
0898   double h = rand.Gaus(0, i);
0899   truthSvtxVtx->set_error(0,0,ae*ae);
0900   truthSvtxVtx->set_error(1,1,d*d+ae*ae);
0901   truthSvtxVtx->set_error(2,2,g*g+h*h+i*i);
0902   truthSvtxVtx->set_error(0,1,ae*d);
0903   truthSvtxVtx->set_error(1,0,ae*d);
0904   truthSvtxVtx->set_error(2,0,ae*g);
0905   truthSvtxVtx->set_error(0,2,ae*g);
0906   truthSvtxVtx->set_error(1,2,d*g+ae*h);
0907   truthSvtxVtx->set_error(2,1,d*g+ae*h);
0908   switch(recoCount()){
0909     case 2:
0910       truthSvtxVtx->insert_track(reco1->get_id());
0911       truthSvtxVtx->insert_track(reco2->get_id());
0912     case 1:
0913       if (reco1) return truthSvtxVtx->insert_track(reco1->get_id());
0914       else return truthSvtxVtx->insert_track(reco2->get_id());
0915       break;
0916     default:
0917       break;
0918   }
0919 }
0920 
0921 void Conversion::PHG4VtxPointToSvtxVertex(SvtxVertex* seedVtx){
0922   truthSvtxVtx = new SvtxVertex_v1();
0923   truthSvtxVtx->set_x(vtx->get_x());
0924   truthSvtxVtx->set_y(vtx->get_y());
0925   truthSvtxVtx->set_z(vtx->get_z());
0926   truthSvtxVtx->set_t0(vtx->get_t());
0927   truthSvtxVtx->set_chisq(1.);
0928   truthSvtxVtx->set_ndof(1);
0929 
0930   for (int i = 0; i < 3; ++i)
0931   {
0932     for (int j = 0; j < 3; ++j)
0933     {
0934       truthSvtxVtx->set_error(i,j,seedVtx->get_error(i,j));
0935     }
0936   }
0937   switch(recoCount()){
0938     case 2:
0939       truthSvtxVtx->insert_track(reco1->get_id());
0940       truthSvtxVtx->insert_track(reco2->get_id());
0941     case 1:
0942       if (reco1) return truthSvtxVtx->insert_track(reco1->get_id());
0943       else return truthSvtxVtx->insert_track(reco2->get_id());
0944       break;
0945     default:
0946       break;
0947   }
0948 }
0949
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!