Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "RecoConversionEval.h"
0002 #include "SVReco.h"
0003 #include "VtxRegressor.h"
0004 
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/getClass.h>
0008 #include <phgenfit/Track.h>
0009 #include <calobase/RawClusterContainer.h>
0010 #include <calobase/RawCluster.h>
0011 #include <trackbase_historic/SvtxTrack.h>
0012 #include <trackbase_historic/SvtxTrackMap.h>
0013 #include <trackbase_historic/SvtxVertexMap.h>
0014 #include <trackbase_historic/SvtxVertex.h>
0015 #include <trackbase_historic/SvtxHitMap.h>
0016 #include <trackbase_historic/SvtxHit.h>
0017 #include <trackbase_historic/SvtxClusterMap.h>
0018 #include <trackbase_historic/SvtxCluster.h>
0019 #include <trackbase/TrkrClusterContainer.h>
0020 #include <trackbase/TrkrClusterv1.h>
0021 #include <g4eval/SvtxEvalStack.h>
0022 #include <g4main/PHG4Particle.h>
0023 #include <g4main/PHG4VtxPoint.h>
0024 #include <g4main/PHG4TruthInfoContainer.h>
0025 
0026 
0027 #include <TTree.h>
0028 #include <TFile.h>
0029 #include <TLorentzVector.h>
0030 
0031 #include <iostream>
0032 #include <cmath> 
0033 #include <algorithm>
0034 #include <sstream>
0035 
0036 using namespace std;
0037 
0038 RecoConversionEval::RecoConversionEval(const std::string &name,std::string tmvamethod,std::string tmvapath) :
0039   SubsysReco("RecoConversionEval"), _fname(name) 
0040 {
0041   _regressor = new VtxRegressor(tmvamethod,tmvapath);
0042 }
0043 
0044 RecoConversionEval::~RecoConversionEval(){
0045   cout<<"deleting RCE"<<endl;
0046   if(_vertexer) delete _vertexer;
0047   if(_regressor) delete _regressor;
0048 }
0049 
0050 int RecoConversionEval::Init(PHCompositeNode *topNode) {
0051   return Fun4AllReturnCodes::EVENT_OK;
0052 }
0053 
0054 int RecoConversionEval::InitRun(PHCompositeNode *topNode) {
0055   _vertexer = new SVReco();
0056   //TODO turn this back into a subsystem and put it on the node tree
0057   _vertexer->InitRun(topNode);
0058   _file = new TFile( _fname.c_str(), "RECREATE");
0059   _treeSignal = new TTree("recoSignal","strong saharah bush");
0060   _treeSignal->SetAutoSave(300);
0061   _treeSignal->Branch("photon_m",   &_b_photon_m);
0062   _treeSignal->Branch("photon_pT",  &_b_photon_pT);
0063   _treeSignal->Branch("photon_eta", &_b_photon_eta);
0064   _treeSignal->Branch("photon_phi", &_b_photon_phi);
0065   _treeSignal->Branch("track1_pT", &_b_track1_pT);
0066   _treeSignal->Branch("track2_pT", &_b_track2_pT);
0067   _treeSignal->Branch("refit", &_b_refit);
0068 
0069   _treeBackground = new TTree("recoBackground","friendly fern");
0070   _treeBackground->SetAutoSave(300);
0071   _treeBackground->Branch("total",   &totalTracks);
0072   _treeBackground->Branch("tracksPEQ",  &passedpTEtaQ);
0073   _treeBackground->Branch("tracks_clus", &passedCluster);
0074   _treeBackground->Branch("pairs", &passedPair);
0075   _treeBackground->Branch("vtx",      &passedVtx);
0076   _treeBackground->Branch("photon",       &passedPhoton);
0077 
0078   return Fun4AllReturnCodes::EVENT_OK;
0079 }
0080 
0081 void RecoConversionEval::doNodePointers(PHCompositeNode *topNode){
0082   _allTracks = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0083   _mainClusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0084   /*These are deprecated
0085    * _svtxClusterMap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
0086    _hitMap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");*/
0087   //new version
0088   _clusterMap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0089   _vertexer->InitEvent(topNode);
0090   //to check if the id is correct
0091   _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0092 
0093 }
0094 
0095 bool RecoConversionEval::hasNodePointers()const{
0096   return _allTracks &&_mainClusterContainer && _vertexer;
0097 }
0098 
0099 int RecoConversionEval::process_event(PHCompositeNode *topNode) {
0100   doNodePointers(topNode);
0101   /*the is not optimized but is just a nlogn process*/
0102   for ( SvtxTrackMap::Iter iter = _allTracks->begin(); iter != _allTracks->end(); ++iter) {
0103     //I want to now only check e tracks so check the clusters of the |charge|=1 tracks
0104     totalTracks++;
0105     if (abs(iter->second->get_charge())==1&&iter->second->get_pt()>_kTrackPtCut&&abs(iter->second->get_eta())<1.) //should i have the layer cut?
0106     {
0107       SvtxTrack* thisTrack = iter->second;
0108       passedpTEtaQ++;
0109       RawCluster* bestCluster= _mainClusterContainer->getCluster(thisTrack->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)));
0110       //TODO what if no cluster is found?
0111       if(bestCluster&&bestCluster->get_prob()>=_kEMProbCut){
0112         //loop over the following tracks
0113         passedCluster++;
0114         for (SvtxTrackMap::Iter jter = iter; jter != _allTracks->end(); ++jter)
0115         {
0116           //check that the next track is an opposite charge electron
0117           if (thisTrack->get_charge()*-1==jter->second->get_charge()&&jter->second->get_pt()>_kTrackPtCut&&abs(iter->second->get_eta())<1.)
0118           {
0119             RawCluster* nextCluster= _mainClusterContainer->getCluster(jter->second->get_cal_cluster_id(SvtxTrack::CAL_LAYER(1)));
0120             //what if no cluster is found?
0121             if(nextCluster&&nextCluster->get_prob()>=_kEMProbCut&&pairCuts(thisTrack,jter->second)){
0122               passedPair++;
0123               genfit::GFRaveVertex* vtxCan = _vertexer->findSecondaryVertex(thisTrack,jter->second);
0124               vtxCan=correctSecondaryVertex(vtxCan,thisTrack,jter->second);
0125               if (vtxCan&&vtxCuts(vtxCan))
0126               {
0127                 passedVtx++;
0128                 _b_refit=true;
0129                 std::pair<PHGenFit::Track*,PHGenFit::Track*> refit_tracks = refitTracks(vtxCan,thisTrack,jter->second);
0130                 //attempt to set the photon to the addition of the refit tracks may return NULL
0131                 TLorentzVector* photon= reconstructPhoton(refit_tracks);
0132                 //if photon is null attempt to set it to the addition of the original tracks may return NULL
0133                 if (!photon)
0134                 {
0135                   photon = reconstructPhoton(thisTrack,jter->second);
0136                   _b_refit=false;
0137                 }
0138                 //if the photon is reconstructed record its data 
0139                 if(photonCuts(photon)){
0140                   _b_photon_m = photon->Dot(*photon);
0141                   _b_photon_pT = photon->Pt();
0142                   _b_photon_eta = photon->Eta();
0143                   _b_photon_phi = photon->Phi();
0144                   _b_track1_pT = thisTrack->get_pt();
0145                   _b_track2_pT = jter->second->get_pt();
0146                   passedPhoton++;
0147                   delete photon;
0148                 }
0149                 else{
0150                   _b_photon_m =   -999.;
0151                   _b_photon_pT =  -999.;
0152                   _b_photon_eta = -999.;
0153                   _b_photon_phi = -999.;
0154                   cout<<"photon not reconstructed"<<endl;
0155                 }
0156                 /*FIXME currently this is not a valid way to get the truthparticle because get_truth_track_id returns UINT_MAX
0157                   PHG4Particle* truthparticle = _truthinfo->GetParticle(jter->second->get_truth_track_id());
0158                   PHG4Particle* truthparticle2 = _truthinfo->GetParticle(thisTrack->get_truth_track_id());
0159                   PHG4Particle* parent;
0160                   if(truthparticle) {
0161                   parent = _truthinfo->GetParticle(truthparticle->get_parent_id());
0162                   if(TMath::Abs(truthparticle->get_pid())!=11||(parent&&parent->get_pid()!=22)||truthparticle->get_parent_id()!=truthparticle2->get_parent_id()){
0163                   _b_fake=true;
0164                   }
0165                   else if(parent&&parent->get_pid()==22){
0166                   TLorentzVector ptlv;
0167                   ptlv.SetPxPyPzE(parent->get_px(),parent->get_py(),parent->get_pz(),parent->get_e());
0168                   parent->identify();
0169                   PHG4Particle* grandparent = _truthinfo->GetParticle(parent->get_parent_id());
0170                   if(grandparent) grandparent->identify();
0171                   _b_tphoton_phi = ptlv.Phi();
0172                   _b_tphoton_eta = ptlv.Eta();
0173                   _b_tphoton_pT =  ptlv.Pt();
0174                   }
0175                   else{
0176                   _b_tphoton_phi = -999.;
0177                   _b_tphoton_eta = -999.;
0178                   _b_tphoton_pT =  -999.;
0179                   }
0180                   }//found truth particle
0181                   else{
0182                   cout<<"no truth particle"<<endl;
0183                   _b_tphoton_phi = -999.;
0184                   _b_tphoton_eta = -999.;
0185                   _b_tphoton_pT =  -999.;
0186                   }*/
0187                 _treeSignal->Fill();
0188               }//vtx cuts
0189             }
0190           }
0191         }
0192       }
0193     }
0194   }
0195   return Fun4AllReturnCodes::EVENT_OK;
0196 }
0197 
0198 genfit::GFRaveVertex* RecoConversionEval::correctSecondaryVertex(genfit::GFRaveVertex* recoVertex,SvtxTrack* reco1,SvtxTrack* reco2){
0199   if(!(recoVertex&&reco1&&reco2)) {
0200     return recoVertex;
0201   }
0202   TVector3 nextPos = recoVertex->getPos();
0203   nextPos.SetMagThetaPhi(_regressor->regress(reco1,reco2,recoVertex),nextPos.Theta(),nextPos.Phi());
0204 
0205   using namespace genfit;
0206   // GFRaveVertex* temp = recoVertex;
0207   std::vector<GFRaveTrackParameters*> tracks;
0208   for(unsigned i =0; i<recoVertex->getNTracks();i++){
0209     tracks.push_back(recoVertex->getParameters(i));
0210   }
0211   recoVertex = new GFRaveVertex(nextPos,recoVertex->getCov(),tracks,recoVertex->getNdf(),recoVertex->getChi2(),recoVertex->getId());
0212   //  delete temp; //this caused outside references to seg fault //TODO shared_ptr is better 
0213   return recoVertex;
0214 }
0215 
0216 TLorentzVector* RecoConversionEval::reconstructPhoton(std::pair<PHGenFit::Track*,PHGenFit::Track*> recos){
0217   if (recos.first&&recos.second)
0218   {
0219     cout<<"reconstructing photon from refit tracks"<<endl;
0220     TLorentzVector tlv1;
0221     tlv1.SetVectM(recos.first->get_mom(),_kElectronRestM);
0222     TLorentzVector tlv2;
0223     tlv2.SetVectM(recos.second->get_mom(),_kElectronRestM);
0224     return new TLorentzVector(tlv1+tlv2);
0225   }
0226   else return NULL;
0227 }
0228 
0229 TLorentzVector* RecoConversionEval::reconstructPhoton(SvtxTrack* reco1,SvtxTrack* reco2){
0230   if (reco1&&reco2)
0231   {
0232     cout<<"reconstructing photon from svtx tracks"<<endl;
0233     TLorentzVector tlv1(reco1->get_px(),reco1->get_py(),reco1->get_pz(),
0234         sqrt(_kElectronRestM*_kElectronRestM+reco1->get_p()*reco1->get_p()));
0235     TLorentzVector tlv2(reco2->get_px(),reco2->get_py(),reco2->get_pz(),
0236         sqrt(_kElectronRestM*_kElectronRestM+reco2->get_p()*reco2->get_p()));
0237     return new TLorentzVector(tlv1+tlv2);
0238   }
0239   else return NULL;
0240 }
0241 
0242 std::pair<PHGenFit::Track*,PHGenFit::Track*> RecoConversionEval::refitTracks(genfit::GFRaveVertex* vtx,SvtxTrack* reco1,SvtxTrack* reco2){
0243   std::pair<PHGenFit::Track*,PHGenFit::Track*> r;
0244   if(!vtx)
0245   {
0246     cerr<<"WARNING: No vertex to refit tracks"<<endl;
0247     r.first=NULL;
0248     r.second=NULL;
0249 
0250   }
0251   else{
0252     r.first=_vertexer->refitTrack(vtx,reco1);
0253     r.second=_vertexer->refitTrack(vtx,reco2);
0254   }
0255   return r;
0256 }
0257 
0258 bool RecoConversionEval::pairCuts(SvtxTrack* t1, SvtxTrack* t2)const{
0259   return detaCut(t1->get_eta(),t2->get_eta())&&hitCuts(t1,t2);//TODO add approach distance ?
0260 }
0261 
0262 bool RecoConversionEval::photonCuts(TLorentzVector* photon)const{
0263   return photon&&photon->Dot(*photon)>_kPhotonMmin&&photon->Dot(*photon)<_kPhotonMmax&&_kPhotonPTmin<photon->Pt();
0264 }
0265 
0266 
0267 bool RecoConversionEval::hitCuts(SvtxTrack* reco1, SvtxTrack* reco2)const {
0268   TrkrCluster *c1 = _clusterMap->findCluster(*(reco1->begin_cluster_keys()));
0269   TrkrCluster *c2 = _clusterMap->findCluster(*(reco2->begin_cluster_keys()));
0270   unsigned l1 = TrkrDefs::getLayer(c1->getClusKey());
0271   unsigned l2 = TrkrDefs::getLayer(c2->getClusKey());
0272   //check that the first hits are close enough
0273   return abs(l1-l2)<=_kDLayerCut;
0274 }
0275 
0276 /*bool RecoConversionEval::vtxCuts(genfit::GFRaveVertex* vtxCan, SvtxTrack* t1, SvtxTrack *t2){
0277 //TODO program the cuts invariant mass, pT
0278 return vtxRadiusCut(vtxCan->getPos());
0279 // && vtxTrackRPhiCut(vtxCan->getPos(),t1)&&vtxTrackRPhiCut(vtxCan->getPos(),t2)&& 
0280 //vtxTrackRZCut(vtxCan->getPos(),t1)&&vtxTrackRZCut(vtxCan->getPos(),t2)&&vtxCan->getChi2()>_kVtxChi2Cut;
0281 }*/
0282 
0283 bool RecoConversionEval::vtxCuts(genfit::GFRaveVertex* vtxCan){
0284   return vtxRadiusCut(vtxCan->getPos());
0285 }
0286 
0287 bool RecoConversionEval::vtxTrackRZCut(TVector3 recoVertPos, SvtxTrack* track){
0288   float dR = sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y())-sqrt(track->get_x()*track->get_x()+track->get_y()*track->get_y());
0289   float dZ = recoVertPos.z()-track->get_z();
0290   return sqrt(dR*dR+dZ*dZ)<_kVtxRZCut;
0291 }
0292 
0293 //bool RecoConversionEval::invariantMassCut()
0294 
0295 bool RecoConversionEval::vtxTrackRPhiCut(TVector3 recoVertPos, SvtxTrack* track){
0296   float vtxR=sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y());
0297   float trackR=sqrt(track->get_x()*track->get_x()+track->get_y()*track->get_y());
0298   return sqrt(vtxR*vtxR+trackR*trackR-2*vtxR*trackR*cos(recoVertPos.Phi()-track->get_phi()))<_kVtxRPhiCut;
0299 }
0300 
0301 bool RecoConversionEval::vtxRadiusCut(TVector3 recoVertPos){
0302   return sqrt(recoVertPos.x()*recoVertPos.x()+recoVertPos.y()*recoVertPos.y()) > _kVtxRCut;
0303 }
0304 
0305 int RecoConversionEval::End(PHCompositeNode *topNode) {
0306   cout<<"Did RecoConversionEval with "<<totalTracks<<" total tracks\n\t";
0307   cout<<1-(float)passedpTEtaQ/totalTracks<<"+/-"<<sqrt((float)passedpTEtaQ)/totalTracks<<" of tracks were rejected by pTEtaQ\n\t";
0308   cout<<1-(float)passedCluster/passedpTEtaQ<<"+/-"<<sqrt((float)passedCluster)/passedpTEtaQ<<" of remaining tracks were rejected by cluster\n\t";
0309   cout<<1-(float)passedPair/passedCluster<<"+/-"<<sqrt((float)passedPair)/passedCluster<<" of pairs were rejected by pair cuts\n\t";
0310   cout<<1-(float)passedVtx/passedPair<<"+/-"<<sqrt((float)passedVtx)/passedPair<<" of vtx were rejected by vtx cuts\n\t";
0311   _treeBackground->Fill();
0312   if(_file){
0313     _file->Write();
0314     _file->Close();
0315   }
0316   cout<<"good end"<<endl;
0317   return Fun4AllReturnCodes::EVENT_OK;
0318 }
0319 
0320 bool RecoConversionEval::approachDistance(SvtxTrack *t1,SvtxTrack* t2)const{
0321   static const double eps = 0.000001;
0322   TVector3 u(t1->get_px(),t1->get_py(),t1->get_pz());
0323   TVector3 v(t2->get_px(),t2->get_py(),t2->get_pz());
0324   TVector3 w(t1->get_x()-t2->get_x(),t1->get_x()-t2->get_y(),t1->get_x()-t2->get_z());
0325 
0326   double a = u.Dot(u);
0327   double b = u.Dot(v);
0328   double c = v.Dot(v);
0329   double d = u.Dot(w);
0330   double e = v.Dot(w);
0331 
0332   double D = a*c - b*b;
0333   double sc, tc;
0334   // compute the line parameters of the two closest points
0335   if (D < eps) {         // the lines are almost parallel
0336     sc = 0.0;
0337     tc = (b>c ? d/b : e/c);   // use the largest denominator
0338   }
0339   else {
0340     sc = (b*e - c*d) / D;
0341     tc = (a*e - b*d) / D;
0342   }
0343   // get the difference of the two closest points
0344   u*=sc;
0345   v*=tc;
0346   w+=u;
0347   w-=v;
0348   return w.Mag()<=_kApprochCut;   // return the closest distance 
0349 }
0350
[ 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!