Changeset 49234af in git for external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc
- Timestamp:
- Dec 9, 2014, 1:27:13 PM (10 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 37deb3b, 9e991f8
- Parents:
- f6b6ee7 (diff), e7e90df (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the(diff)
links above to see all the changes relative to each parent. - File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc
rf6b6ee7 r49234af 5 5 // Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason 6 6 // 7 // $Id: MeasureFunction.cc 670 2014-06-06 01:24:42Z jthaler $ 7 8 //---------------------------------------------------------------------- 8 9 // This file is part of FastJet contrib. … … 36 37 37 38 // Return all of the necessary TauComponents for specific input particles and axes 38 TauComponents MeasureFunction::result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) { 39 TauComponents MeasureFunction::result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const { 40 41 // first find partition 42 // this sets jetPartitionStorage and beamPartitionStorage 43 PseudoJet beamPartitionStorage; 44 std::vector<fastjet::PseudoJet> jetPartitionStorage = get_partition(particles,axes,&beamPartitionStorage); 45 46 // then return result calculated from partition 47 return result_from_partition(jetPartitionStorage,axes,&beamPartitionStorage); 48 } 39 49 40 std::vector<double> jetPieces(axes.size(), 0.0); 41 double beamPiece = 0.0; 50 std::vector<fastjet::PseudoJet> MeasureFunction::get_partition(const std::vector<fastjet::PseudoJet>& particles, 51 const std::vector<fastjet::PseudoJet>& axes, 52 PseudoJet * beamPartitionStorage) const { 42 53 43 // Calculates the unnormalized sub-tau values, i.e. a std::vector of the contributions to tau_N of each Voronoi region (or region within R_0) 54 std::vector<std::vector<PseudoJet> > jetPartition(axes.size()); 55 std::vector<PseudoJet> beamPartition; 56 57 // Figures out the partiting of the input particles into the various jet pieces 58 // Based on which axis the parition is closest to 44 59 for (unsigned i = 0; i < particles.size(); i++) { 45 60 … … 60 75 61 76 if (j_min == -1) { 62 if (_has_beam) beamP iece += beam_numerator(particles[i]);77 if (_has_beam) beamPartition.push_back(particles[i]); 63 78 else assert(_has_beam); // this should never happen. 64 79 } else { 65 jetP ieces[j_min] += jet_numerator(particles[i],axes[j_min]);80 jetPartition[j_min].push_back(particles[i]); 66 81 } 67 82 } 68 83 69 // Calculates normalization for tau and subTau if _has_denominator is true, otherwise returns 1.0 (i.e. no normalization) 70 double tauDen = 0.0; 71 if (_has_denominator) { 72 for (unsigned i = 0; i < particles.size(); i++) { 73 tauDen += denominator(particles[i]); 74 } 75 } else { 76 tauDen = 1.0; // if no denominator, then 1.0 for no normalization factor 84 // Store beam partition 85 if (beamPartitionStorage) { 86 *beamPartitionStorage = join(beamPartition); 87 } 88 89 // Store jet partitions 90 std::vector<PseudoJet> jetPartitionStorage(axes.size(),PseudoJet(0,0,0,0)); 91 for (unsigned j = 0; j < axes.size(); j++) { 92 jetPartitionStorage[j] = join(jetPartition[j]); 77 93 } 78 94 95 return jetPartitionStorage; 96 } 97 98 // does partition, but only stores index of PseudoJets 99 std::vector<std::list<int> > MeasureFunction::get_partition_list(const std::vector<fastjet::PseudoJet>& particles, 100 const std::vector<fastjet::PseudoJet>& axes) const { 101 102 std::vector<std::list<int> > jetPartition(axes.size()); 103 104 // Figures out the partiting of the input particles into the various jet pieces 105 // Based on which axis the parition is closest to 106 for (unsigned i = 0; i < particles.size(); i++) { 107 108 // find minimum distance; start with beam (-1) for reference 109 int j_min = -1; 110 double minRsq; 111 if (_has_beam) minRsq = beam_distance_squared(particles[i]); 112 else minRsq = std::numeric_limits<double>::max(); // make it large value 113 114 // check to see which axis the particle is closest to 115 for (unsigned j = 0; j < axes.size(); j++) { 116 double tempRsq = jet_distance_squared(particles[i],axes[j]); // delta R distance 117 if (tempRsq < minRsq) { 118 minRsq = tempRsq; 119 j_min = j; 120 } 121 } 122 123 if (j_min == -1) { 124 assert(_has_beam); // consistency check 125 } else { 126 jetPartition[j_min].push_back(i); 127 } 128 } 129 130 return jetPartition; 131 } 132 133 134 // Uses existing partition and calculates result 135 // TODO: Can we cache this for speed up when doing area subtraction? 136 TauComponents MeasureFunction::result_from_partition(const std::vector<fastjet::PseudoJet>& jet_partition, 137 const std::vector<fastjet::PseudoJet>& axes, 138 PseudoJet * beamPartitionStorage) const { 139 140 std::vector<double> jetPieces(axes.size(), 0.0); 141 double beamPiece = 0.0; 142 143 double tauDen = 0.0; 144 if (!_has_denominator) tauDen = 1.0; // if no denominator, then 1.0 for no normalization factor 145 146 // first find jet pieces 147 for (unsigned j = 0; j < axes.size(); j++) { 148 std::vector<PseudoJet> thisPartition = jet_partition[j].constituents(); 149 for (unsigned i = 0; i < thisPartition.size(); i++) { 150 jetPieces[j] += jet_numerator(thisPartition[i],axes[j]); //numerator jet piece 151 if (_has_denominator) tauDen += denominator(thisPartition[i]); // denominator 152 } 153 } 154 155 // then find beam piece 156 if (_has_beam) { 157 assert(beamPartitionStorage); // make sure I have beam information 158 std::vector<PseudoJet> beamPartition = beamPartitionStorage->constituents(); 159 160 for (unsigned i = 0; i < beamPartition.size(); i++) { 161 beamPiece += beam_numerator(beamPartition[i]); //numerator beam piece 162 if (_has_denominator) tauDen += denominator(beamPartition[i]); // denominator 163 } 164 } 79 165 return TauComponents(jetPieces, beamPiece, tauDen, _has_denominator, _has_beam); 80 166 } 167 168 169 170 81 171 82 172 } //namespace contrib
Note:
See TracChangeset
for help on using the changeset viewer.