SiCalo/SiliconSeedAna/SiliconCaloMatching.cc

0001 #include "SiliconCaloMatching.h"
0002 #include "SiliconCaloTrack_v1.h"
0003 #include "SiliconCaloTrackMap_v1.h"
0004
0005 #include <trackbase/TrkrDefs.h>
0006 #include <trackbase/ActsGeometry.h>
0007 #include <trackbase/TrkrCluster.h>
0008 #include <trackbase/TrkrClusterContainer.h>
0009
0010 //track map include
0011 #include <trackbase_historic/SvtxTrack.h>
0012 #include <trackbase_historic/SvtxTrackMap.h>
0013
0014 //vertex include
0015 #include <globalvertex/SvtxVertex.h>
0016 #include <globalvertex/SvtxVertexMap.h>
0017
0018 //calo include
0019 #include <calobase/RawCluster.h>
0020 #include <calobase/RawClusterContainer.h>
0021
0022 #include <phool/PHCompositeNode.h>
0023 #include <phool/getClass.h>
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025
0026 #include <TVector3.h>
0027
0028 // for sort algorithm
0029 bool compLayer(TrkrDefs::cluskey a, TrkrDefs::cluskey b)
0030 {
0031   int layer_a = TrkrDefs::getLayer(a);
0032   int layer_b = TrkrDefs::getLayer(b);
0033   return (layer_a < layer_b);
0034 }
0035
0036 //____________________________________________________________________________..
0037 SiliconCaloMatching::SiliconCaloMatching(const std::string &name)
0038     : SubsysReco(name)
0039 {
0040 }
0041
0042 #define LOG(msg) std::cout << "[SiliconCaloMatching] " << msg << std::endl;
0043
0044
0045 //____________________________________________________________________________..
0046 int SiliconCaloMatching::InitRun(PHCompositeNode * topNode)
0047 {
0048
0049   PHNodeIterator iter(topNode);
0050
0051   //Looking for the DST node
0052   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0053   if (!dstNode)
0054   {
0055     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0056     return Fun4AllReturnCodes::ABORTRUN;
0057   }
0058   PHNodeIterator iter_dst(dstNode);
0059
0060   // Create the SiliconCalo and association nodes if required
0061   _sicalomap = findNode::getClass<SiliconCaloTrackMap>(dstNode, "SiliconCaloTrack");
0062   if (!_sicalomap)
0063   {
0064     PHNodeIterator dstiter(dstNode);
0065     PHCompositeNode* DetNode =
0066       dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "SICALO"));
0067     if (!DetNode)
0068     {
0069       DetNode = new PHCompositeNode("SICALO");
0070       dstNode->addNode(DetNode);
0071     }
0072
0073     _sicalomap = new SiliconCaloTrackMap_v1;
0074     PHIODataNode<PHObject>* SiliconCaloTrackMapNode =
0075       new PHIODataNode<PHObject>(_sicalomap, "SiliconCaloTrack", "PHObject");
0076     DetNode->addNode(SiliconCaloTrackMapNode);
0077   }
0078
0079   return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 //____________________________________________________________________________..
0083 int SiliconCaloMatching::process_event(PHCompositeNode* topNode)
0084 {
0085   if(false)
0086     std::cout<<"SiliconCaloMatching::process_event(PHCompositeNode* topNode)"<<std::endl;
0087
0088   if(!getNodes(topNode))
0089     return Fun4AllReturnCodes::EVENT_OK;
0090
0091   //if(false)
0092     //if((evt%1000)==0) std::cout << "start track map  EVENT " << evt << " is OK" << std::endl;
0093
0094   // vertex position from crossing=0
0095   TVector3 vtxPos(0,0,0);
0096   //if (!b_skipvtx)
0097   {
0098     if (!_vertexmap->empty())
0099     {
0100       // crossing==0 has high priority
0101       for (const auto &[k, v] : *_vertexmap)
0102       {
0103         if (v && v->get_beam_crossing() == 0)
0104         {
0105           vtxPos = TVector3(v->get_x(), v->get_y(), v->get_z()) ;
0106           break;
0107         }
0108       }
0109     }
0110   }
0111
0112 ///////////////////////////////////////////////
0113   for (auto &iter : *_trackmap)
0114   {
0115     auto track = iter.second;
0116     if (!track)
0117       continue;
0118
0119     int trkcrossing = track->get_crossing();
0120     if (trkcrossing != 0 && !isMC)
0121       continue;
0122
0123     //--int   t_id      = track->get_id();
0124     float t_eta     = track->get_eta();
0125     float t_phi     = track->get_phi();
0126     float t_pt      = track->get_pt();
0127     //--int   t_charge  = track->get_charge();
0128     //--float t_chi2ndf = track->get_quality();
0129     float t_x       = track->get_x();
0130     float t_y       = track->get_y();
0131     float t_z       = track->get_z();
0132     //--float t_px      = track->get_px();
0133     //--float t_py      = track->get_py();
0134     //--float t_pz      = track->get_pz();
0135     //--int t_crossing  = trkcrossing;
0136     //--if(t_crossing==0) evt_nsiseed0++;
0137
0138 //    int t_nmaps = 0, t_nintt = 0, t_inner = 0, t_outer = 0;
0139
0140
0141     if (false)
0142       std::cout << "track_x : " << t_x << ", track_y: " << t_y << ", track_z: " << t_z
0143                 << ", track_eta: " << t_eta << ", track_phi: " << t_phi << ", track_pt: " << t_pt << std::endl;
0144
0145     // track projection to EMCAL
0146     SvtxTrackState *emcalState    = track->get_state(_caloRadiusEMCal);
0147  //   SvtxTrackState *emcalOutState = track->get_state(_caloRadiusEMCal+_caloThicknessEMCal);
0148
0149
0150     if(!emcalState){
0151       std::cout<<"No TrackState (EMC projection object)"<<std::endl;
0152       continue;
0153     }
0154
0155     // track projection to EMC surface
0156     TVector3 projPos(emcalState->get_x(), emcalState->get_y(), emcalState->get_z());
0157     //float phi_proj = projPos.phi();
0158     //float z_proj =   projPos.z();
0159
0160     // --- EMCal cluster matching ---
0161     // Find the closest EMCal cluster in eta/phi to this state
0162     float best_dR  = 9999.0;
0163     int   best_idx = -1;
0164
0165     float best_dphi=9999., best_dz=9999.;
0166
0167     const float dphi_sigma = 3.; // dphi_sigma is factor 3 wider than dz_sigma
0168
0169
0170     //search closest EMCAL cluster
0171     /*******************
0172      need to speed up
0173     *******************/
0174     //--std::cout<<"start EMCalClusmap " << std::endl;
0175     //--std::cout << "EMCalClusmap size: " << EMCalClusmap->size() << std::endl;
0176     for (unsigned int i = 0; i < _emcClusmap->size(); i++)
0177     {
0178       auto *emc = _emcClusmap->getCluster(i);
0179       if (!emc)
0180         continue;
0181       if (emc->get_energy() < _emcal_low_cut)
0182         continue;
0183
0184       //TVector3 emcPos(emc->get_x(), emc->get_y(), emc->get_z());
0185
0186       float dz   = emc->get_z()   - projPos.z(); // cm unit
0187
0188       float dphi = emc->get_phi() - projPos.Phi(); // radian
0189       if (dphi >  M_PI) dphi -= 2 * M_PI;
0190       if (dphi < -M_PI) dphi += 2 * M_PI;
0191
0192
0193       float dphi_cm = dphi * _caloRadiusEMCal; // (Radius=93.5cm) cm
0194
0195       float dR = sqrt(pow( (dphi_cm/dphi_sigma), 2) + pow(dz, 2));
0196
0197       if(dR<best_dR){
0198         best_dR  = dR;
0199         best_idx =  i;
0200
0201         best_dphi = dphi;
0202         best_dz   = dz;
0203       }
0204     }
0205
0206     //
0207     if(best_idx==-1) {
0208       std::cout<<"associated cluseter not found"<<std::endl;
0209       continue;
0210     }
0211
0212     // associated cluster
0213     RawCluster *emcClus = _emcClusmap->getCluster(best_idx);
0214
0215     // momentum update
0216     //
0217     //
0218     float pt_calo = calculatePt(track, emcClus);
0219     //--std::cout << "caloPt : " << pt_calo<<std::endl;
0220
0221
0222     // Fill updated value to SvtxTrack and SiliconCaloTrack object
0223     //track->set_px( pt_calo*cos(track->get_phi()) );
0224     //track->set_py( pt_calo*sin(track->get_phi()) );
0225     //
0226     auto sicalo = std::make_unique<SiliconCaloTrack_v1>();
0227
0228     sicalo->set_id(track->get_id());
0229     sicalo->set_calo_id(best_idx);
0230     sicalo->set_pt(pt_calo);
0231     sicalo->set_phi(track->get_phi());
0232     sicalo->set_eta(track->get_eta());
0233     sicalo->set_cal_dphi(best_dphi);
0234     sicalo->set_cal_dz(best_dz);
0235
0236     //_sicalomap->insert(track->get_id(), sicalo.release());
0237     _sicalomap->insertWithKey(sicalo.release(), track->get_id());
0238
0239   } // track loop
0240
0241   //--std::cout << "EVENT " << evt << " is OK" << std::endl;
0242
0243 ///////////////////////////////////////////////
0244
0245   return Fun4AllReturnCodes::EVENT_OK;
0246 }
0247
0248 float SiliconCaloMatching::calculatePt(SvtxTrack* track, RawCluster* emc)
0249 {
0250   if(track==nullptr||emc==nullptr) {
0251     return -9999;
0252   }
0253
0254   TrackSeed* si_seed = track->get_silicon_seed();
0255
0256   // Si position
0257   if(si_seed!=nullptr&& si_seed->size_cluster_keys()<=2) return -9999.;
0258
0259   // checking cluster list to get clusters from 2 outer layers
0260   // loop over SiClusters with the track
0261   std::vector<TrkrDefs::cluskey> vClusKey;
0262   for (auto key_iter = si_seed->begin_cluster_keys(); key_iter != si_seed->end_cluster_keys(); ++key_iter)
0263   {
0264     const auto &cluster_key  = *key_iter;
0265     TrkrCluster* trkrCluster = _clustermap->findCluster(cluster_key);
0266     if (!trkrCluster)
0267     {
0268       continue;
0269     }
0270     vClusKey.push_back(cluster_key);
0271   }
0272   std::sort(vClusKey.begin(), vClusKey.end(), compLayer); // sort by layer
0273   //std::cout<<"Nclus : "<<si_seed->size_cluster_keys()<<" "<<vClusKey.size()<<std::endl;
0274   //
0275
0276   std::vector<TrkrDefs::cluskey>::reverse_iterator ritr = vClusKey.rbegin();
0277   TrkrDefs::cluskey ckey_outer = *ritr; ritr++;
0278   TrkrDefs::cluskey ckey_inner = *ritr;
0279   if(false)
0280     std::cout<<"c layer "<<(int)TrkrDefs::getLayer(ckey_inner)<<" "<<(int)TrkrDefs::getLayer(ckey_outer)<<std::endl;
0281
0282   TrkrCluster* oClus = _clustermap->findCluster(ckey_outer);
0283   TrkrCluster* iClus = _clustermap->findCluster(ckey_inner);
0284
0285   Acts::Vector3 oCpos = _geometry->getGlobalPosition(ckey_outer, oClus);
0286   Acts::Vector3 iCpos = _geometry->getGlobalPosition(ckey_inner, iClus);
0287
0288   //float phi_intt = (oCpos-iCpos).phi();
0289
0290
0291   // EMC position
0292   float emc_x, emc_y, emc_z;
0293   getCaloPosition(emc, emc_x, emc_y, emc_z);
0294
0295
0296   // pT calculation
0297   float phi_intt = atan2(oCpos.y()-iCpos.y(), oCpos.x()-iCpos.x());
0298   float phi_calo = atan2(emc_y - oCpos.y(), emc_x - oCpos.x());
0299
0300   int charge = track->get_charge();
0301
0302   float dphi = phi_calo - phi_intt;
0303
0304   //float pt_calo = 0.21*pow((double)(-charge*dphi), -0.986);//cal_CaloPt(dphi);
0305   float pt_calo = 0.21*pow(fabs(dphi), -0.986);//cal_CaloPt(dphi);
0306   //float pt_inv = charge*0.02+4.9*(-charge*dphi)-0.6*pow(-charge*dphi, 2);
0307   //float pt_calo = 1./pt_inv;
0308
0309   // cout
0310   float pt = track->get_pt();
0311
0312   if(false)
0313     std::cout<<"phi : "<<pt_calo<<" "<<phi_intt<<" "<<phi_calo<<" "<<dphi<<", "<<charge<<", "<<pt<<" "<<emc_x<<" "<<emc_y<<" "<<emc_z<<std::endl;
0314
0315   return pt_calo;
0316 }
0317
0318
0319 //____________________________________________________________________________..
0320 int SiliconCaloMatching::EndRun(const int /*runnumber*/)
0321 {
0322
0323   return Fun4AllReturnCodes::EVENT_OK;
0324 }
0325
0326 //____________________________________________________________________________..
0327 int SiliconCaloMatching::End(PHCompositeNode * /*unused*/)
0328 {
0329   return Fun4AllReturnCodes::EVENT_OK;
0330 }
0331
0332 bool SiliconCaloMatching::getNodes(PHCompositeNode* topNode)
0333 {
0334
0335   _trackmap   = findNode::getClass<SvtxTrackMap>(        topNode, _trackMapName);
0336   _vertexmap  = findNode::getClass<SvtxVertexMap>(       topNode, _vertexMapName);
0337   _geometry   = findNode::getClass<ActsGeometry>(        topNode, _actsgeometryName);
0338   _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, _clusterContainerName);
0339   _emcClusmap = findNode::getClass<RawClusterContainer>( topNode, _emcalClusName);
0340
0341
0342   if (!_trackmap)
0343   {
0344     std::cout << PHWHERE << "Missing trackmap, can't continue" << std::endl;
0345     return false;
0346   }
0347   if (!_vertexmap)
0348   {
0349     std::cout << PHWHERE << "Missing vertexmap, can't continue" << std::endl;
0350     return false;
0351   }
0352   if (!_geometry)
0353   {
0354     std::cout << PHWHERE << "Missing geometry, can't continue" << std::endl;
0355     return false;
0356   }
0357   if (!_clustermap)
0358   {
0359     std::cout << PHWHERE << "Missing clustermap, can't continue" << std::endl;
0360     return false;
0361   }
0362   if (!_emcClusmap)
0363   {
0364     std::cout << PHWHERE << "Missing EMC clustermap, can't continue" << std::endl;
0365     return false;
0366   }
0367
0368   return true;
0369 }
0370
0371 void SiliconCaloMatching::getCaloPosition(RawCluster *calo, float &x, float &y, float &z)
0372 {
0373   if(calo==nullptr) {
0374     x = y = z = std::numeric_limits<float>::quiet_NaN();
0375     return;
0376   }
0377
0378   x   = calo->get_x();
0379   y   = calo->get_y();
0380   z   = calo->get_z();
0381   float r   = calo->get_r();
0382   float phi = calo->get_phi();
0383   float rr = sqrt(x*x + y*y);
0384   //std::cout<<"emc_x, y: "<<emc_x<<" "<<emc_y<<std::endl;
0385   if(fabs(rr - r)>1)
0386   {
0387     //std::cout<<"no emc_x, y: "<<emc_x<<" "<<emc_y<<" "<<rr
0388     //         <<", replaced by r + phi"<<emc_r<<" "<<emc_phi<<" : ";
0389     x = r * cos(phi);
0390     y = r * sin(phi);
0391     //std::cout<<emc_x<<" "<<emc_y<<std::endl;
0392   }
0393 }
0394