TrackAndClusterMatchingQA/src/TrackAndClusterMatchingQA.cc

0001 #include "TrackAndClusterMatchingQA.h"
0002
0003 /// Cluster/Calorimeter includes
0004 #include <calobase/RawCluster.h>
0005 #include <calobase/RawClusterContainer.h>
0006 #include <calobase/RawClusterUtility.h>
0007 #include <calobase/RawTower.h>
0008 #include <calobase/RawTowerContainer.h>
0009 #include <calobase/RawTowerDefs.h>
0010 #include <calobase/RawTowerGeom.h>
0011 #include <calobase/RawTowerGeomContainer.h>
0012 #include <calobase/TowerInfo.h>
0013 #include <calobase/TowerInfoContainer.h>
0014 #include <calotrigger/CaloTriggerInfo.h>
0015
0016 #include <phool/phool.h>
0017
0018 /// Tracking includes
0019 #include <globalvertex/GlobalVertex.h>
0020 #include <globalvertex/GlobalVertexMap.h>
0021 #include <globalvertex/SvtxVertexMap.h>
0022 #include <globalvertex/SvtxVertex.h>
0023 #include <trackbase_historic/SvtxTrack.h>
0024 #include <trackbase_historic/SvtxTrackMap.h>
0025 #include <trackbase_historic/SvtxTrackState.h>
0026 #include <trackbase_historic/TrackAnalysisUtils.h>
0027 #include <trackbase_historic/TrackSeed.h>
0028
0029
0030 #include <Acts/Definitions/Algebra.hpp>
0031
0032 #include <CLHEP/Vector/ThreeVector.h>
0033
0034 /// Fun4All includes
0035 #include <fun4all/Fun4AllHistoManager.h>
0036 #include <fun4all/Fun4AllReturnCodes.h>
0037 #include <g4main/PHG4Hit.h>
0038 #include <g4main/PHG4Particle.h>
0039 #include <g4main/PHG4TruthInfoContainer.h>
0040 #include <phool/PHCompositeNode.h>
0041 #include <phool/getClass.h>
0042
0043
0044 #include <g4main/PHG4Particle.h>            // for PHG4Particle
0045 #include <g4main/PHG4TruthInfoContainer.h>  // for PHG4TruthInfoContainer
0046 #include <g4main/PHG4VtxPoint.h>            // for PHG4VtxPoint
0047 #include <trackbase_historic/SvtxPHG4ParticleMap_v1.h>
0048 #include <kfparticle_sphenix/KFParticle_truthAndDetTools.h>
0049
0050 /// ROOT includes
0051 #include <TFile.h>
0052 #include <TH1F.h>
0053 #include <TH2F.h>
0054 #include <TMath.h>
0055 #include <TNtuple.h>
0056 #include <TTree.h>
0057 #include <TDatabasePDG.h>
0058 #include <TParticlePDG.h>
0059
0060 /// C++ includes
0061 #include <cassert>
0062 #include <sstream>
0063 #include <string>
0064
0065 using namespace std;
0066
0067 /**
0068  * TrackAndClusterMatchingQA is a class developed for jet QA
0069  * Author: Antonio Silva (antonio.sphenix@gmail.com)
0070  */
0071
0072 /**
0073  * Constructor of module
0074  */
0075 TrackAndClusterMatchingQA::TrackAndClusterMatchingQA(const std::string &name, const std::string &filename)
0076   : SubsysReco(name)
0077   , _outfilename(filename)
0078 {
0079   /// Initialize variables and trees so we don't accidentally access
0080   /// memory that was never allocated
0081 }
0082
0083 /**
0084  * Destructor of module
0085  */
0086 TrackAndClusterMatchingQA::~TrackAndClusterMatchingQA()
0087 {
0088
0089 }
0090
0091 /**
0092  * Initialize the module and prepare looping over events
0093  */
0094 int TrackAndClusterMatchingQA::Init(PHCompositeNode *topNode)
0095 {
0096   if (Verbosity() > 5)
0097   {
0098     cout << "Beginning Init in TrackAndClusterMatchingQA" << endl;
0099   }
0100
0101   _outfile = new TFile(_outfilename.c_str(), "RECREATE");
0102
0103   _h2trackPt_vs_clusterEt = new TH2F("h2trackPt_vs_clusterEt", ";track #it{p}_{T} (GeV/#it{c});EMCal cluster #it{E}_{T} (GeV)", 40, 0., 20., 40, 0., 20.);
0104   _h1EMCal_TowerEnergy = new TH1F("_h1EMCal_TowerEnergy", ";EMCal tower #it{E} (GeV);Entries", 40, 0., 20.);
0105   _h1EMCal_TowerEnergy_Retowered = new TH1F("_h1EMCal_TowerEnergy_Retowered", ";EMCal retorewed tower #it{E} (GeV);Entries", 40, 0., 20.);
0106   _h1HCalIN_TowerEnergy = new TH1F("_h1HCalIN_TowerEnergy", ";HCalIN tower #it{E} (GeV);Entries", 40, 0., 20.);
0107   _h1HCalOUT_TowerEnergy = new TH1F("_h1HCalOUT_TowerEnergy", ";HCalOUT tower #it{E} (GeV);Entries", 40, 0., 20.);
0108   _h1EMCal_TowerEnergy_afterSelection = new TH1F("_h1EMCal_TowerEnergy_afterSelection", ";EMCal tower #it{E} (GeV);Entries", 40, 0., 20.);
0109   _h1EMCal_TowerEnergy_Retowered_afterSelection = new TH1F("_h1EMCal_TowerEnergy_Retowered_afterSelection", ";EMCal retorewed tower #it{E} (GeV);Entries", 40, 0., 20.);
0110   _h1HCalIN_TowerEnergy_afterSelection = new TH1F("_h1HCalIN_TowerEnergy_afterSelection", ";HCalIN tower #it{E} (GeV);Entries", 40, 0., 20.);
0111   _h1HCalOUT_TowerEnergy_afterSelection = new TH1F("_h1HCalOUT_TowerEnergy_afterSelection", ";HCalOUT tower #it{E} (GeV);Entries", 40, 0., 20.);
0112   _h2TowerEnergy_EMCal_vs_HCalIN = new TH2F("_h2TowerEnergy_EMCal_vs_HCalIN", ";EMCal Retowered Towers #it{E} (GeV);Inner HCal Towers #it{E} (GeV)", 40, 0., 20., 40, 0., 20.);
0113   _h2TowerEnergy_EMCal_vs_HCalOUT = new TH2F("_h2TowerEnergy_EMCal_vs_HCalOUT", ";EMCal Retowered Towers #it{E} (GeV);Outer HCal Towers #it{E} (GeV)", 40, 0., 20., 40, 0., 20.);
0114   _h2TowerEnergy_HCalIN_vs_HCalOUT = new TH2F("_h2TowerEnergy_HCalIN_vs_HCalOUT", ";EMCal Retowered Towers #it{E} (GeV);Outer HCal Towers #it{E} (GeV)", 40, 0., 20., 40, 0., 20.);
0115   _h1EMCal_RetowerEnergyFraction = new TH1F("_h1EMCal_RetowerEnergyFraction", ";#it{E}^{max}_{tower}/#it{E}_{retowered};Entries", 20, 0., 1.);
0116
0117   _h1Track_Quality = new TH1F("_h1Track_Quality", ";Quality;Entries", 100, 0., 50);
0118   _h1Track_DCAxy = new TH1F("_h1Track_DCAxy", ";DCA xy;Entries", 200, -0.02, 0.02);
0119   _h1Track_DCAz = new TH1F("_h1Track_DCAz", ";DCA z;Entries", 200, -0.02, 0.02);
0120   _h1Track_TPC_Hits = new TH1F("_h1Track_TPC_Hits", ";Number of TPC hits;Entries", 50, -0.5, 49.5);
0121   _h1Track_Silicon_Hits = new TH1F("_h1Track_Silicon_Hits", ";Number of TPC hits;Entries", 50, -0.5, 49.5);
0122   _h1Track_Pt_beforeSelections = new TH1F("_h1Track_Pt_beforeSelections", ";track #it{p}_{T};Entries", 100, 0., 50);
0123   _h1Track_Pt_afterSelections = new TH1F("_h1Track_Pt_afterSelections", ";track #it{p}_{T};Entries", 100, 0., 50);
0124
0125   _h1Track_TPC_Hits_Selected = new TH1F("_h1Track_TPC_Hits_Selected", ";Number of TPC hits;Entries", 50, -0.5, 49.5);
0126   _h2Track_TPC_Hits_vs_Phi = new TH2F("_h2Track_TPC_Hits_vs_Phi", ";Number of TPC hits;track #phi", 50, -0.5, 49.5, 63, -M_PI, M_PI);
0127   _h2Track_TPC_Hits_vs_Eta = new TH2F("_h2Track_TPC_Hits_vs_Eta", ";Number of TPC hits;track #eta", 50, -0.5, 49.5, 22, -1.1, 1.1);
0128   _h2Track_TPC_Hits_vs_Pt = new TH2F("_h2Track_TPC_Hits_vs_Pt", ";Number of TPC hits;track #it{p}_{T} (GeV/#it{c})", 50, -0.5, 49.5, 40, 0., 20);
0129
0130   _h1deta = new TH1F("hdeta", "Cluster #deta; #deta; Entries", 20, -0.2, 0.2);    //deta distribution
0131   _h1dphi  = new TH1F("hdphi", "Cluster #dphi; #dphi; Entries", 50, -0.15, 0.15);   //dphi distribution
0132   _h1min_dR  = new TH1F("hdR", "Cluster #dR; #dR; Entries", 100,0.,5.);   //jet delta radius
0133   _h2phi_vs_deta = new TH2F("_h2phi_vs_deta", ";#dphi; #deta",  50, -0.15, 0.15, 20, -0.2, 0.2); // deta Vs. dphi distribution
0134
0135   return Fun4AllReturnCodes::EVENT_OK;
0136 }
0137
0138 /**
0139  * Main workhorse function where each event is looped over and
0140  * data from each event is collected from the node tree for analysis
0141  */
0142 int TrackAndClusterMatchingQA::process_event(PHCompositeNode *topNode)
0143 {
0144   if (Verbosity() > 5)
0145   {
0146     cout << "Beginning process_event in TrackAndClusterMatchingQA" << endl;
0147   }
0148
0149   GlobalVertex *vtx = nullptr;
0150   CLHEP::Hep3Vector vertex(0., 0., 0.);
0151
0152   bool has_vertex = false;
0153
0154   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0155   if (vertexmap)
0156   {
0157     if (vertexmap->empty())
0158     {
0159       if(Verbosity() > 5) std::cout << "GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0160     }
0161     else
0162     {
0163       vtx = vertexmap->begin()->second;
0164       vertex.setX(vtx->get_x());
0165       vertex.setY(vtx->get_y());
0166       vertex.setZ(vtx->get_z());
0167       has_vertex = true;
0168     }
0169   }
0170   else
0171   {
0172     if(Verbosity() > 5) std::cout << "GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0173   }
0174
0175   SvtxVertex *svtx_vtx = nullptr;
0176
0177   SvtxVertexMap *vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0178
0179   if(vertexMap)
0180   {
0181     if(!vertexMap->empty())
0182     {
0183       svtx_vtx = vertexMap->begin()->second;
0184       vertex.setX(svtx_vtx->get_x());
0185       vertex.setY(svtx_vtx->get_y());
0186       vertex.setZ(svtx_vtx->get_z());
0187       has_vertex = true;
0188     }
0189   }
0190
0191   if(has_vertex == false && (!_use_origin_when_no_vertex))
0192   {
0193     return Fun4AllReturnCodes::ABORTEVENT;
0194   }
0195
0196   _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0197
0198   if (!_track_map)
0199   {
0200     cout << PHWHERE
0201          << "SvtxTrackMap node is missing, can't collect tracks"
0202          << endl;
0203     return Fun4AllReturnCodes::ABORTEVENT;
0204   }
0205
0206   _em_clusters = findNode::getClass<RawClusterContainer>(topNode, _em_clusters_name);
0207
0208   // read in tower geometries
0209   RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0210   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0211   RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0212
0213   if(!geomEM || !geomIH || !geomOH)
0214   {
0215     std::cout << "FATAL ERROR, cannot find tower geometry containers" << std::endl;
0216     return Fun4AllReturnCodes::ABORTEVENT;
0217   }
0218
0219   //Radius of the front face of the EMCal
0220   float cemcradius = geomEM->get_radius();
0221
0222   if (!_em_clusters)
0223   {
0224     std::cout << PHWHERE << _em_clusters_name << " node is missing, can't collect clusters" << std::endl;
0225     return Fun4AllReturnCodes::ABORTEVENT;
0226   }
0227
0228   TowerInfoContainer *towerinfosEM_retower = findNode::getClass<TowerInfoContainer>(topNode, _em_retowered_towers_name);
0229   TowerInfoContainer *towerinfosEM = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0230   TowerInfoContainer *towerinfosIH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0231   TowerInfoContainer *towerinfosOH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0232
0233   if (!towerinfosEM)
0234   {
0235     std::cout << " TrackAndClusterMatchingQA::process_event : container with EM towers does not exist, aborting " << std::endl;
0236     return Fun4AllReturnCodes::ABORTEVENT;
0237   }
0238   if (!towerinfosIH)
0239   {
0240     std::cout << " TrackAndClusterMatchingQA::process_event : container with HCalIN towers does not exist, aborting " << std::endl;
0241     return Fun4AllReturnCodes::ABORTEVENT;
0242   }
0243   if (!towerinfosOH)
0244   {
0245     std::cout << " TrackAndClusterMatchingQA::process_event : container with HCalOUT towers does not exist, aborting " << std::endl;
0246     return Fun4AllReturnCodes::ABORTEVENT;
0247   }
0248
0249   SvtxTrack *track = nullptr;
0250
0251   for (SvtxTrackMap::Iter iter = _track_map->begin(); iter != _track_map->end();++iter)
0252   {
0253     track = iter->second;
0254
0255     if(!track) continue;
0256
0257     if(track->get_pt() < _track_min_pt)
0258     {
0259       continue;
0260     }
0261
0262     _h1Track_Pt_beforeSelections->Fill(track->get_pt());
0263
0264     if(!IsAcceptableTrack(track, vtx))
0265     {
0266       continue;
0267     }
0268
0269     _h1Track_Pt_afterSelections->Fill(track->get_pt());
0270
0271     //Get track projection to the front face of the EMCal
0272     SvtxTrackState* cemcstate = track->get_state(cemcradius);
0273
0274     float track_phi = track->get_phi();
0275     float track_eta = track->get_eta();
0276
0277     //Calculate eta and phi of the track projection to the EMCal
0278     if(cemcstate)
0279     {
0280       track_phi = atan2(cemcstate->get_y(), cemcstate->get_x());
0281       track_eta = asinh(cemcstate->get_z()/sqrt(cemcstate->get_x()*cemcstate->get_x() + cemcstate->get_y()*cemcstate->get_y()));
0282     }
0283
0284     float min_dR = 999.;
0285     RawCluster *cluster_match = nullptr;
0286
0287     RawClusterContainer::ConstRange begin_end_EMC = _em_clusters->getClusters();
0288     RawClusterContainer::ConstIterator clusIter_EMC;
0289
0290     /// Loop over the EMCal clusters
0291     for (clusIter_EMC = begin_end_EMC.first; clusIter_EMC != begin_end_EMC.second; ++clusIter_EMC)
0292     {
0293       RawCluster *cluster = clusIter_EMC->second;
0294
0295       float cluster_phi = RawClusterUtility::GetAzimuthAngle(*cluster, vertex);
0296       float cluster_eta = RawClusterUtility::GetPseudorapidity(*cluster, vertex);
0297
0298       float dR = GetDeltaR(track_phi, cluster_phi, track_eta, cluster_eta);
0299
0300       if(dR > 0.2)
0301       {
0302         continue;
0303       }
0304
0305       if(dR < min_dR)
0306       {
0307         min_dR = dR;
0308         cluster_match = cluster;
0309       }
0310     }
0311
0312     if(cluster_match)
0313     {
0314       _h2trackPt_vs_clusterEt->Fill(track->get_pt(), RawClusterUtility::GetET(*cluster_match, vertex));
0315
0316       float cluster_phi = RawClusterUtility::GetAzimuthAngle(*cluster_match, vertex);
0317       float cluster_eta = RawClusterUtility::GetPseudorapidity(*cluster_match, vertex);
0318       float deta = track_eta - cluster_eta;   //added for eta analysis
0319       float dphi = track_phi - cluster_phi;   //added for phi analysis
0320
0321       _h1deta->Fill(deta);
0322       _h1dphi->Fill(dphi);
0323       _h2phi_vs_deta->Fill(dphi,deta);
0324
0325     }
0326
0327       _h1min_dR->Fill(min_dR);
0328
0329   }
0330
0331   TowerInfo *EMCal_towerInfo = nullptr;
0332   TowerInfo *HCalIN_towerInfo = nullptr;
0333   TowerInfo *HCalOUT_towerInfo = nullptr;
0334   for(unsigned int i = 0; i < towerinfosEM_retower->size(); i++)
0335   {
0336     EMCal_towerInfo = towerinfosEM_retower->get_tower_at_channel(i);
0337     HCalIN_towerInfo = towerinfosIH->get_tower_at_channel(i);
0338     HCalOUT_towerInfo = towerinfosOH->get_tower_at_channel(i);
0339
0340     _h1EMCal_TowerEnergy_Retowered->Fill(EMCal_towerInfo->get_energy());
0341     _h1HCalIN_TowerEnergy->Fill(HCalIN_towerInfo->get_energy());
0342     _h1HCalOUT_TowerEnergy->Fill(HCalOUT_towerInfo->get_energy());
0343
0344     if(_apply_tower_selection)
0345     {
0346       if((!EMCal_towerInfo->get_isGood()) || (!HCalIN_towerInfo->get_isGood()) || (!HCalOUT_towerInfo->get_isGood()))
0347       {
0348         continue;
0349       }
0350     }
0351
0352     _h1EMCal_TowerEnergy_Retowered_afterSelection->Fill(EMCal_towerInfo->get_energy());
0353     _h1HCalIN_TowerEnergy_afterSelection->Fill(HCalIN_towerInfo->get_energy());
0354     _h1HCalOUT_TowerEnergy_afterSelection->Fill(HCalOUT_towerInfo->get_energy());
0355
0356     _h2TowerEnergy_EMCal_vs_HCalIN->Fill(EMCal_towerInfo->get_energy(), HCalIN_towerInfo->get_energy());
0357     _h2TowerEnergy_EMCal_vs_HCalOUT->Fill(EMCal_towerInfo->get_energy(), HCalOUT_towerInfo->get_energy());
0358     _h2TowerEnergy_HCalIN_vs_HCalOUT->Fill(HCalIN_towerInfo->get_energy(), HCalOUT_towerInfo->get_energy());
0359   }
0360
0361   const int NETA = geomIH->get_etabins();
0362   const int NPHI = geomIH->get_phibins();
0363
0364   float EMCAL_RETOWER_E[NETA][NPHI];
0365   float EMCAL_RETOWER_HIGHEST_E[NETA][NPHI];
0366
0367   for (int eta = 0; eta < NETA; eta++)
0368   {
0369     for (int phi = 0; phi < NPHI; phi++)
0370     {
0371       EMCAL_RETOWER_E[eta][phi] = 0;
0372       EMCAL_RETOWER_HIGHEST_E[eta][phi] = 0;
0373     }
0374   }
0375
0376   for(unsigned int i = 0; i < towerinfosEM->size(); i++)
0377   {
0378     EMCal_towerInfo = towerinfosEM->get_tower_at_channel(i);
0379     _h1EMCal_TowerEnergy->Fill(EMCal_towerInfo->get_energy());
0380
0381     if(_apply_tower_selection && (!EMCal_towerInfo->get_isGood()))
0382     {
0383       continue;
0384     }
0385
0386     _h1EMCal_TowerEnergy_afterSelection->Fill(EMCal_towerInfo->get_energy());
0387
0388     //Reproducing what is done in the RetowerCEMC module
0389     unsigned int channelkey = towerinfosEM->encode_key(i);
0390       int ieta = towerinfosEM->getTowerEtaBin(channelkey);
0391       int iphi = towerinfosEM->getTowerPhiBin(channelkey);
0392       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0393       RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
0394       int this_IHetabin = geomIH->get_etabin(tower_geom->get_eta());
0395
0396     double fractionalcontribution[3] = {0};
0397
0398     std::pair<double, double> range_embin= geomEM->get_etabounds(tower_geom->get_bineta());
0399     for(int etabin_iter = -1;etabin_iter <= 1;etabin_iter++)
0400     {
0401       if(this_IHetabin+etabin_iter < 0 || this_IHetabin+etabin_iter >= NETA){continue;}
0402           std::pair<double, double> range_ihbin= geomIH->get_etabounds(this_IHetabin + etabin_iter);
0403           if(range_ihbin.first <= range_embin.first && range_ihbin.second >= range_embin.second)
0404       {
0405         fractionalcontribution[etabin_iter+1] = 1;
0406       }
0407           else if(range_ihbin.first <= range_embin.first && range_ihbin.second < range_embin.second  && range_embin.first < range_ihbin.second)
0408       {
0409         fractionalcontribution[etabin_iter+1] = (range_ihbin.second - range_embin.first) / (range_embin.second- range_embin.first);
0410       }
0411           else if(range_ihbin.first > range_embin.first && range_ihbin.second >= range_embin.second && range_embin.second > range_ihbin.first)
0412       {
0413         fractionalcontribution[etabin_iter+1] = (range_embin.second - range_ihbin.first) / (range_embin.second- range_embin.first);
0414       }
0415           else
0416       {
0417         fractionalcontribution[etabin_iter+1] = 0;
0418       }
0419     }
0420
0421     int this_IHphibin = geomIH->get_phibin(tower_geom->get_phi());
0422       float this_E = EMCal_towerInfo->get_energy();
0423       for (int etabin_iter = -1 ; etabin_iter <= 1;etabin_iter++)
0424     {
0425       if (this_IHetabin+etabin_iter < 0 || this_IHetabin+etabin_iter >= NETA)
0426       {
0427         continue;
0428       }
0429       EMCAL_RETOWER_E[this_IHetabin+etabin_iter][this_IHphibin] += this_E * fractionalcontribution[etabin_iter+1];
0430       if(this_E * fractionalcontribution[etabin_iter+1] > EMCAL_RETOWER_HIGHEST_E[this_IHetabin+etabin_iter][this_IHphibin])
0431       {
0432         EMCAL_RETOWER_HIGHEST_E[this_IHetabin+etabin_iter][this_IHphibin] = this_E * fractionalcontribution[etabin_iter+1];
0433       }
0434     }
0435
0436   }
0437
0438   //Calculate the ratio between the highest energy tower in the retowered tower
0439   for (int eta = 0; eta < NETA; eta++)
0440   {
0441     for (int phi = 0; phi < NPHI; phi++)
0442     {
0443       if(EMCAL_RETOWER_E[eta][phi] > 0)
0444       {
0445         _h1EMCal_RetowerEnergyFraction->Fill(EMCAL_RETOWER_HIGHEST_E[eta][phi]/EMCAL_RETOWER_E[eta][phi]);
0446       }
0447     }
0448   }
0449
0450
0451   return Fun4AllReturnCodes::EVENT_OK;
0452 }
0453
0454 /**
0455  * End the module and finish any data collection. Clean up any remaining
0456  * loose ends
0457  */
0458 int TrackAndClusterMatchingQA::End(PHCompositeNode *topNode)
0459 {
0460   if (Verbosity() > 1)
0461   {
0462     cout << "Ending TrackAndClusterMatchingQA analysis package" << endl;
0463   }
0464
0465   /// Change to the outfile
0466   _outfile->cd();
0467
0468   /// Write and close the outfile
0469   _outfile->Write();
0470   _outfile->Close();
0471
0472   /// Clean up pointers and associated histos/trees in TFile
0473   delete _outfile;
0474
0475   if (Verbosity() > 1)
0476   {
0477     cout << "Finished TrackAndClusterMatchingQA analysis package" << endl;
0478   }
0479
0480   return 0;
0481 }
0482
0483 float TrackAndClusterMatchingQA::GetDeltaR(float track_phi, float cluster_phi, float track_eta, float cluster_eta)
0484 {
0485   float deta = track_eta - cluster_eta;
0486   float dphi = track_phi - cluster_phi;
0487   if(dphi > M_PI) dphi -= 2*M_PI;
0488   if(dphi < -M_PI) dphi += 2*M_PI;
0489
0490   return sqrt(dphi*dphi + deta*deta);
0491 }
0492
0493 bool TrackAndClusterMatchingQA::IsAcceptableTrack(SvtxTrack *track, GlobalVertex *vtx)
0494 {
0495   Acts::Vector3 vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0496
0497   auto dcapair = TrackAnalysisUtils::get_dca(track, vertex);
0498
0499   float dcaxy = dcapair.first.first;
0500   //float dcaxy_err = dcapair.first.second;
0501   _h1Track_DCAxy->Fill(dcaxy);
0502   float dcaz = dcapair.second.first;
0503   //float dcaz_err = dcapair.second.second;
0504   _h1Track_DCAz->Fill(dcaz);
0505
0506   float quality = track->get_quality();
0507   _h1Track_Quality->Fill(quality);
0508
0509
0510   auto silicon_seed = track->get_silicon_seed();
0511   auto tpc_seed = track->get_tpc_seed();
0512   int nsiliconhits = 0;
0513   int nTPChits = 0;
0514
0515   // Getting number of TPC hits as done in SvtxEvaluator.h
0516   if(tpc_seed)
0517   {
0518     for (TrackSeed::ConstClusterKeyIter iter = tpc_seed->begin_cluster_keys(); iter != tpc_seed->end_cluster_keys(); ++iter)
0519     {
0520       TrkrDefs::cluskey cluster_key = *iter;
0521       unsigned int layer = TrkrDefs::getLayer(cluster_key);
0522       if (_nlayers_tpc > 0 && layer >= (_nlayers_maps + _nlayers_intt) && layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
0523       {
0524         nTPChits++;
0525       }
0526     }
0527   }
0528
0529   if (silicon_seed)
0530   {
0531     for (TrackSeed::ConstClusterKeyIter iter = silicon_seed->begin_cluster_keys(); iter != silicon_seed->end_cluster_keys(); ++iter)
0532     {
0533       TrkrDefs::cluskey cluster_key = *iter;
0534       unsigned int layer = TrkrDefs::getLayer(cluster_key);
0535       if (_nlayers_tpc > 0 && layer >= (_nlayers_maps + _nlayers_intt) && layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
0536       {
0537         nTPChits++;
0538       }
0539     }
0540     nsiliconhits = silicon_seed->size_cluster_keys(); // Get number of hits on silicon detectors (INTT+MVTX)
0541   }
0542
0543   _h1Track_TPC_Hits->Fill(nTPChits);
0544   _h2Track_TPC_Hits_vs_Phi->Fill(nTPChits, track->get_phi());
0545   _h2Track_TPC_Hits_vs_Eta->Fill(nTPChits, track->get_eta());
0546   _h2Track_TPC_Hits_vs_Pt->Fill(nTPChits, track->get_pt());
0547   _h1Track_Silicon_Hits->Fill(nsiliconhits);
0548
0549   if(quality > _track_quality) return false;
0550
0551   // DCA units and values have to be cross-checked
0552   //std::cout << "DCAxy: " << dcaxy << std::endl;
0553   if(std::fabs(dcaxy) > _track_max_dcaxy) return false; // will accept everything, need to check units (microns?)
0554   //std::cout << "DCAz: " << dcaz << std::endl;
0555   if(std::fabs(dcaz) > _track_max_dcaz) return false; // will accept everything, need to check units (microns?)
0556
0557   if (nsiliconhits < _track_min_silicon_hits) return false;
0558
0559   _h1Track_TPC_Hits_Selected->Fill(nTPChits);
0560
0561   if(nTPChits < _track_min_tpc_hits) return false;
0562
0563
0564   return true;
0565 }