HF-Particle/HFTrigger/HFTrigger.cc

0001 #include "HFTrigger.h"
0002
0003 /*
0004  * Basic heavy flavor software trigger
0005  * Used in study of hardware trigger
0006  * Cameron Dean
0007  * 01/02/2021
0008  */
0009 std::map<std::string, bool> triggerDecisions =
0010 {
0011   {"oneTrack", false}
0012 , {"twoTrack", false}
0013 , {"lowMultiplicity", false}
0014 , {"highMultiplicity", false}
0015 };
0016
0017 namespace HFTriggerRequirement
0018 {
0019   float meanHighMult = 100; //Set min num. INTT hits
0020   float asymmHighMult = 0.1; //Set highMult asymm between INTT layers (fraction)
0021   float meanLowMult = 20;  //Set max num. INTT hits
0022   float asymmLowMult = 0.1; //Set lowMult asymm between INTT layers (fraction)
0023   float trackPT = 0.5; //Min track pT in GeV
0024   float trackVertexDCA = 0.01; //Min. DCA of a track with any vertex in cm
0025   float trackTrackDCA = 0.03; //Max. DCA of a track with other tracks in cm
0026 }
0027
0028 HFTrigger::HFTrigger()
0029   : SubsysReco("HFTRIGGER")
0030   , m_useOneTrackTrigger(false)
0031   , m_useTwoTrackTrigger(false)
0032   , m_useLowMultiplicityTrigger(false)
0033   , m_useHighMultiplicityTrigger(false)
0034 {
0035 }
0036
0037 HFTrigger::HFTrigger(const std::string &name)
0038   : SubsysReco(name)
0039   , m_useOneTrackTrigger(false)
0040   , m_useTwoTrackTrigger(false)
0041   , m_useLowMultiplicityTrigger(false)
0042   , m_useHighMultiplicityTrigger(false)
0043 {
0044 }
0045
0046 int HFTrigger::Init(PHCompositeNode *topNode)
0047 {
0048   return 0;
0049 }
0050
0051 int HFTrigger::process_event(PHCompositeNode *topNode)
0052 {
0053   bool successfulTrigger = runTrigger(topNode);
0054
0055   if (Verbosity() >= VERBOSITY_MORE)
0056   {
0057     if (successfulTrigger) std::cout << "One of the heavy flavor triggers fired" << std::endl;
0058     else std::cout << "No heavy flavor triggers fired" << std::endl;
0059   }
0060
0061   int failedTriggerDecisions = 0;
0062   if (m_useOneTrackTrigger && !triggerDecisions.find("oneTrack")->second)
0063   {
0064     if (Verbosity() >= VERBOSITY_SOME) std::cout << "The oneTrackTrigger did not fire for this event, skipping." << std::endl;
0065     failedTriggerDecisions += 1;
0066   }
0067   if (m_useTwoTrackTrigger && !triggerDecisions.find("twoTrack")->second)
0068   {
0069     if (Verbosity() >= VERBOSITY_SOME) std::cout << "The twoTrackTrigger did not fire for this event, skipping." << std::endl;
0070     failedTriggerDecisions += 1;
0071   }
0072   if (m_useLowMultiplicityTrigger && !triggerDecisions.find("lowMultiplicity")->second)
0073   {
0074     if (Verbosity() >= VERBOSITY_SOME) std::cout << "The lowMultiplicityTrigger did not fire for this event, skipping." << std::endl;
0075     failedTriggerDecisions += 1;
0076   }
0077   if (m_useHighMultiplicityTrigger && !triggerDecisions.find("highMultiplicity")->second)
0078   {
0079     if (Verbosity() >= VERBOSITY_SOME) std::cout << "The highMultiplicityTrigger did not fire for this event, skipping." << std::endl;
0080     failedTriggerDecisions += 1;
0081   }
0082
0083   if (failedTriggerDecisions > 0) return Fun4AllReturnCodes::ABORTEVENT;
0084   else return Fun4AllReturnCodes::EVENT_OK;
0085 }
0086
0087 int HFTrigger::End(PHCompositeNode *topNode)
0088 {
0089   return 0;
0090 }
0091
0092 bool HFTrigger::runTrigger(PHCompositeNode *topNode)
0093 {
0094   bool anyTriggerFired = false;
0095
0096   std::vector<Track> allTracks = makeAllTracks(topNode);
0097   std::vector<Vertex> allVertices = makeAllPrimaryVertices(topNode);
0098
0099   float meanMult = 0;
0100   float asymmMult = 0;
0101   calculateMultiplicity(topNode, meanMult, asymmMult);
0102
0103   if (m_useOneTrackTrigger) triggerDecisions.find("oneTrack")->second = runOneTrackTrigger(allTracks, allVertices);
0104   if (m_useTwoTrackTrigger) triggerDecisions.find("twoTrack")->second = runTwoTrackTrigger(allTracks, allVertices);
0105   if (m_useLowMultiplicityTrigger) triggerDecisions.find("lowMultiplicity")->second = runLowMultiplicityTrigger(meanMult, asymmMult);
0106   if (m_useHighMultiplicityTrigger) triggerDecisions.find("highMultiplicity")->second = runHighMultiplicityTrigger(meanMult, asymmMult);
0107
0108   if (Verbosity() >= VERBOSITY_SOME) printTrigger();
0109
0110   const int numberOfFiredTriggers = std::accumulate( begin(triggerDecisions), end(triggerDecisions), 0,
0111                                                    [](const int previous, const std::pair<const std::string, bool>& element)
0112                                                    { return previous + element.second; });
0113
0114   if (numberOfFiredTriggers != 0) anyTriggerFired = true;
0115
0116   return anyTriggerFired;
0117 }
0118
0119 bool HFTrigger::runOneTrackTrigger(std::vector<Track> Tracks, std::vector<Vertex> Vertices)
0120 { //Can maybe use return when we get a positive trigger decision to break the loops
0121   bool oneTrackTriggerDecision = false;
0122
0123   for (Track track : Tracks)
0124   {
0125     float pT = sqrt(pow(track(3,0), 2) + pow(track(4,0), 2));
0126     float minIP = FLT_MAX;
0127     for (Vertex vertex : Vertices)
0128     {
0129       float thisIP = calcualteTrackVertexDCA(track, vertex);
0130       minIP = std::min(minIP, thisIP);
0131     }
0132     if (pT > HFTriggerRequirement::trackPT && minIP > HFTriggerRequirement::trackVertexDCA) oneTrackTriggerDecision = true;
0133   }
0134
0135   return oneTrackTriggerDecision;
0136 }
0137
0138 bool HFTrigger::runTwoTrackTrigger(std::vector<Track> Tracks, std::vector<Vertex> Vertices)
0139 {
0140   bool twoTrackTriggerDecision = false;
0141   std::vector<Track> goodTracks;
0142
0143   for (Track track : Tracks)
0144   {
0145     std::vector<Track> singleTrack = {track};
0146     bool oneTrackTriggerDecision = runOneTrackTrigger(singleTrack, Vertices);
0147     if (oneTrackTriggerDecision) goodTracks.push_back(track);
0148   }
0149
0150   if (goodTracks.size() > 1) //We need at least two good track to start a two track trigger
0151   {
0152     for (unsigned int i = 0; i < goodTracks.size(); i++)
0153     {
0154       for (unsigned int j = i + 1; j < goodTracks.size(); j++)
0155       {
0156         float trackDCA = calcualteTrackTrackDCA(goodTracks[i], goodTracks[j]);
0157         if (trackDCA < HFTriggerRequirement::trackTrackDCA) twoTrackTriggerDecision = true;
0158       }
0159     }
0160   }
0161
0162   return twoTrackTriggerDecision;
0163 }
0164
0165 void HFTrigger::calculateMultiplicity(PHCompositeNode *topNode, float& meanMultiplicity, float& asymmetryMultiplicity)
0166 {
0167   TrkrHitSetContainer* hitContainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0168
0169   TrkrHitSetContainer::ConstRange inttHitSetRange[2];
0170   for (int i = 0; i < 2; i++) inttHitSetRange[i] = hitContainer->getHitSets(TrkrDefs::TrkrId::inttId, i + 4);
0171
0172   int inttHits[2] = {0};
0173
0174   for (int i = 0; i < 2; i++)
0175   {
0176     for (TrkrHitSetContainer::ConstIterator hitsetitr = inttHitSetRange[i].first;
0177          hitsetitr != inttHitSetRange[i].second;
0178          ++hitsetitr)
0179     {
0180       TrkrHitSet *hitset = hitsetitr->second;
0181       inttHits[i] += hitset->size();
0182     }
0183   }
0184
0185   meanMultiplicity = (inttHits[0] + inttHits[1])/2;
0186   asymmetryMultiplicity = (inttHits[0] - inttHits[1])/(inttHits[0] + inttHits[1]);
0187 }
0188
0189 bool HFTrigger::runHighMultiplicityTrigger(float meanMultiplicity, float asymmetryMultiplicity)
0190 {
0191   bool multiplicityTriggerDecision = false;
0192   float meanRequirement = HFTriggerRequirement::meanHighMult;
0193   float asymmRequirement = HFTriggerRequirement::asymmHighMult;
0194
0195   if ( meanMultiplicity > meanRequirement && asymmetryMultiplicity < asymmRequirement)
0196     multiplicityTriggerDecision = true;
0197
0198   return multiplicityTriggerDecision;
0199 }
0200
0201 bool HFTrigger::runLowMultiplicityTrigger(float meanMultiplicity, float asymmetryMultiplicity)
0202 {
0203   bool multiplicityTriggerDecision = false;
0204   float meanRequirement = HFTriggerRequirement::meanLowMult;
0205   float asymmRequirement = HFTriggerRequirement::asymmLowMult;
0206
0207   if ( meanMultiplicity < meanRequirement && asymmetryMultiplicity < asymmRequirement)
0208     multiplicityTriggerDecision = true;
0209
0210   return multiplicityTriggerDecision;
0211 }
0212
0213 Vertex HFTrigger::makeVertex(PHCompositeNode *topNode)
0214 {
0215   Vertex vertex;
0216
0217   vertex(0,0) = m_dst_vertex->get_x();
0218   vertex(1,0) = m_dst_vertex->get_y();
0219   vertex(2,0) = m_dst_vertex->get_z();
0220
0221   return vertex;
0222 }
0223
0224 std::vector<Vertex> HFTrigger::makeAllPrimaryVertices(PHCompositeNode *topNode)
0225 {
0226   std::vector<Vertex> primaryVertices;
0227   m_dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0228
0229   for (SvtxVertexMap::ConstIter iter = m_dst_vertexmap->begin(); iter != m_dst_vertexmap->end(); ++iter)
0230   {
0231     m_dst_vertex = iter->second;
0232     primaryVertices.push_back(makeVertex(topNode));
0233   }
0234
0235   return primaryVertices;
0236 }
0237
0238 Track HFTrigger::makeTrack(PHCompositeNode *topNode)
0239 {
0240   Track track;
0241
0242   track(0,0) = m_dst_track->get_x();
0243   track(1,0) = m_dst_track->get_y();
0244   track(2,0) = m_dst_track->get_z();
0245   track(3,0) = m_dst_track->get_px();
0246   track(4,0) = m_dst_track->get_py();
0247   track(5,0) = m_dst_track->get_pz();
0248
0249   return track;
0250 }
0251
0252 std::vector<Track> HFTrigger::makeAllTracks(PHCompositeNode *topNode)
0253 {
0254   std::vector<Track> tracks;
0255   m_dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0256
0257   for (SvtxTrackMap::Iter iter = m_dst_trackmap->begin(); iter != m_dst_trackmap->end(); ++iter)
0258   {
0259     m_dst_track = iter->second;
0260     tracks.push_back(makeTrack(topNode));
0261   }
0262
0263   return tracks;
0264 }
0265
0266 int HFTrigger::decomposeTrack(Track track, TrackX& trackPosition, TrackP& trackMomentum)
0267 {
0268   for (unsigned int i = 0; i < 3; i++)
0269   {
0270     trackPosition(i,0) = track(i,0);
0271     trackMomentum(i,0) = track(i+3,0);
0272   }
0273
0274   return 0;
0275 }
0276
0277 float HFTrigger::calcualteTrackVertex2DDCA(Track track, Vertex vertex)
0278 {
0279   TrackX pos;
0280   TrackP mom;
0281   DCA dcaVertex;
0282
0283   decomposeTrack(track, pos, mom);
0284
0285   dcaVertex = pos - vertex - mom.dot(pos - vertex)*mom/(mom.dot(mom));
0286
0287   float twoD_DCA = std::sqrt(std::pow(dcaVertex(0,0), 2) + std::pow(dcaVertex(1,0), 2));
0288
0289   return twoD_DCA;
0290 }
0291
0292 float HFTrigger::calcualteTrackVertexDCA(Track track, Vertex vertex)
0293 {
0294   TrackX pos;
0295   TrackP mom;
0296   DCA dcaVertex;
0297
0298   decomposeTrack(track, pos, mom);
0299
0300   dcaVertex = pos - vertex - mom.dot(pos - vertex)*mom/(mom.dot(mom));
0301
0302   return std::abs(dcaVertex.norm());
0303 }
0304
0305 float HFTrigger::calcualteTrackTrackDCA(Track trackOne, Track trackTwo)
0306 {
0307   TrackX posOne;
0308   TrackP momOne;
0309   TrackX posTwo;
0310   TrackP momTwo;
0311   float dcaTrackSize;
0312
0313   decomposeTrack(trackOne, posOne, momOne);
0314   decomposeTrack(trackTwo, posTwo, momTwo);
0315
0316   Eigen::Vector3f momOneCrossmomTwo = momOne.cross(momTwo);
0317   Eigen::Vector3f posOneMinusposTwo = posOne - posTwo;
0318   dcaTrackSize = std::abs(momOneCrossmomTwo.dot(posOneMinusposTwo)/(momOneCrossmomTwo.norm()));
0319
0320   return dcaTrackSize;
0321 }
0322
0323 void HFTrigger::printTrigger()
0324 {
0325   std::cout << "\n---------------HFTrigger information---------------" << std::endl;
0326   if (m_useOneTrackTrigger) std::cout << "The oneTrack trigger decision is " << triggerDecisions.find("oneTrack")->second << std::endl;
0327   if (m_useTwoTrackTrigger) std::cout << "The twoTrack trigger decision is " << triggerDecisions.find("twoTrack")->second << std::endl;
0328   if (m_useLowMultiplicityTrigger) std::cout << "The lowMultiplicity trigger decision is " << triggerDecisions.find("lowMultiplicity")->second << std::endl;
0329   if (m_useHighMultiplicityTrigger) std::cout << "The highMultiplicity trigger decision is " << triggerDecisions.find("highMultiplicity")->second << std::endl;
0330   std::cout << "---------------------------------------------------\n" << std::endl;
0331 }
0332