Changes in modules/VertexFinder.cc [95b4e9f:77e9ae1] in git
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/VertexFinder.cc
r95b4e9f r77e9ae1 7 7 */ 8 8 9 10 9 #include "modules/VertexFinder.h" 11 10 #include "classes/DelphesClasses.h" … … 14 13 #include "classes/DelphesPileUpReader.h" 15 14 15 #include "ExRootAnalysis/ExRootClassifier.h" 16 #include "ExRootAnalysis/ExRootFilter.h" 16 17 #include "ExRootAnalysis/ExRootResult.h" 17 #include "ExRootAnalysis/ExRootFilter.h" 18 #include "ExRootAnalysis/ExRootClassifier.h" 19 18 19 #include "TDatabasePDG.h" 20 #include "TFormula.h" 21 #include "TLorentzVector.h" 20 22 #include "TMath.h" 23 #include "TMatrixT.h" 24 #include "TObjArray.h" 25 #include "TRandom3.h" 21 26 #include "TString.h" 22 #include "TFormula.h"23 #include "TRandom3.h"24 #include "TObjArray.h"25 #include "TDatabasePDG.h"26 #include "TLorentzVector.h"27 #include "TMatrixT.h"28 27 #include "TVector3.h" 29 28 29 #include <algorithm> 30 #include <iostream> 31 #include <map> 32 #include <stdexcept> 33 #include <string> 30 34 #include <utility> 31 #include <algorithm>32 #include <stdexcept>33 #include <iostream>34 35 #include <vector> 35 #include <map>36 #include <string>37 36 38 37 using namespace std; 39 38 40 static const Double_t mm 41 static const Double_t m = 1000. *mm;42 static const Double_t ns 43 static const Double_t s = 1.e+9 * ns;44 static const Double_t c_light = 2.99792458e+8 * m/s;39 static const Double_t mm = 1.; 40 static const Double_t m = 1000. * mm; 41 static const Double_t ns = 1.; 42 static const Double_t s = 1.e+9 * ns; 43 static const Double_t c_light = 2.99792458e+8 * m / s; 45 44 46 45 //------------------------------------------------------------------------------ … … 84 83 //------------------------------------------------------------------------------ 85 84 86 static Bool_t secondAscending 85 static Bool_t secondAscending(pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1) 87 86 { 88 87 return (pair0.second < pair1.second); 89 88 } 90 89 91 static Bool_t secondDescending 90 static Bool_t secondDescending(pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1) 92 91 { 93 92 return (pair0.second > pair1.second); … … 101 100 102 101 // Clear the track and cluster maps before starting 103 trackIDToDouble.clear 104 trackIDToInt.clear 105 trackIDToBool.clear 106 clusterIDToDouble.clear 107 clusterIDToInt.clear 108 clusterIDToBool.clear 109 trackPT.clear 110 clusterSumPT2.clear 102 trackIDToDouble.clear(); 103 trackIDToInt.clear(); 104 trackIDToBool.clear(); 105 clusterIDToDouble.clear(); 106 clusterIDToInt.clear(); 107 clusterIDToBool.clear(); 108 trackPT.clear(); 109 clusterSumPT2.clear(); 111 110 112 111 // Create the initial cluster seeds 113 createSeeds 112 createSeeds(); 114 113 115 114 // In order of descending seed pt, grow each cluster. If a cluster ends up with 116 115 // fewer than MinNDF tracks, release the tracks for other clusters to claim. 117 sort (clusterSumPT2.begin (), clusterSumPT2.end(), secondDescending);118 for (vector<pair<UInt_t, Double_t> >::const_iterator cluster = clusterSumPT2.begin (); cluster != clusterSumPT2.end(); cluster++)119 120 121 if (!clusterIDToInt.at (cluster->first).at("ndf"))122 123 124 125 if(fGrowSeeds)126 growCluster(cluster->first);127 128 129 130 131 if ((Int_t) clusterIDToInt.at (cluster->first).at("ndf") < fMinNDF)132 133 for (map<UInt_t, map<string, Int_t> >::iterator track = trackIDToInt.begin (); track != trackIDToInt.end(); track++)134 135 if (track->second.at ("clusterIndex") != (Int_t)cluster->first)136 137 138 139 140 141 142 trackIDToBool[clusterIDToInt.at (cluster->first).at("seed")]["claimed"] = true;143 116 sort(clusterSumPT2.begin(), clusterSumPT2.end(), secondDescending); 117 for(vector<pair<UInt_t, Double_t> >::const_iterator cluster = clusterSumPT2.begin(); cluster != clusterSumPT2.end(); cluster++) 118 { 119 // Skip the cluster if it no longer has any tracks 120 if(!clusterIDToInt.at(cluster->first).at("ndf")) 121 continue; 122 123 // Grow the cluster if GrowSeeds is true 124 if(fGrowSeeds) 125 growCluster(cluster->first); 126 127 // If the cluster still has fewer than MinNDF tracks, release the tracks; 128 // otherwise, mark the seed track as claimed 129 130 if((Int_t)clusterIDToInt.at(cluster->first).at("ndf") < fMinNDF) 131 { 132 for(map<UInt_t, map<string, Int_t> >::iterator track = trackIDToInt.begin(); track != trackIDToInt.end(); track++) 133 { 134 if(track->second.at("clusterIndex") != (Int_t)cluster->first) 135 continue; 136 track->second["clusterIndex"] = -1; 137 trackIDToBool[track->first]["claimed"] = false; 138 } 139 } 140 else 141 trackIDToBool[clusterIDToInt.at(cluster->first).at("seed")]["claimed"] = true; 142 } 144 143 145 144 // Add tracks to the output array after updating their ClusterIndex. 146 fItInputArray->Reset 147 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))148 149 if (candidate->Momentum.Pt () < fMinPT || fabs (candidate->Momentum.Eta()) > fMaxEta)150 151 candidate->ClusterIndex = trackIDToInt.at (candidate->GetUniqueID ()).at("clusterIndex");152 153 145 fItInputArray->Reset(); 146 while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) 147 { 148 if(candidate->Momentum.Pt() < fMinPT || fabs(candidate->Momentum.Eta()) > fMaxEta) 149 continue; 150 candidate->ClusterIndex = trackIDToInt.at(candidate->GetUniqueID()).at("clusterIndex"); 151 fOutputArray->Add(candidate); 152 } 154 153 155 154 // Add clusters with at least MinNDF tracks to the output array in order of 156 155 // descending sum(pt**2). 157 clusterSumPT2.clear 158 for (map<UInt_t, map<string, Int_t> >::const_iterator cluster = clusterIDToInt.begin (); cluster != clusterIDToInt.end(); cluster++)159 { 160 161 if (cluster->second.at("ndf") < fMinNDF)162 continue; 163 clusterSumPT2.push_back (make_pair (cluster->first, clusterIDToDouble.at (cluster->first).at("sumPT2")));164 } 165 sort (clusterSumPT2.begin (), clusterSumPT2.end(), secondDescending);166 167 for (vector<pair<UInt_t, Double_t> >::const_iterator cluster = clusterSumPT2.begin (); cluster != clusterSumPT2.end(); cluster++)156 clusterSumPT2.clear(); 157 for(map<UInt_t, map<string, Int_t> >::const_iterator cluster = clusterIDToInt.begin(); cluster != clusterIDToInt.end(); cluster++) 158 { 159 160 if(cluster->second.at("ndf") < fMinNDF) 161 continue; 162 clusterSumPT2.push_back(make_pair(cluster->first, clusterIDToDouble.at(cluster->first).at("sumPT2"))); 163 } 164 sort(clusterSumPT2.begin(), clusterSumPT2.end(), secondDescending); 165 166 for(vector<pair<UInt_t, Double_t> >::const_iterator cluster = clusterSumPT2.begin(); cluster != clusterSumPT2.end(); cluster++) 168 167 { 169 168 DelphesFactory *factory = GetFactory(); … … 171 170 172 171 candidate->ClusterIndex = cluster->first; 173 candidate->ClusterNDF = clusterIDToInt.at (cluster->first).at("ndf");172 candidate->ClusterNDF = clusterIDToInt.at(cluster->first).at("ndf"); 174 173 candidate->ClusterSigma = fSigma; 175 174 candidate->SumPT2 = cluster->second; 176 candidate->Position.SetXYZT(0.0, 0.0, clusterIDToDouble.at (cluster->first).at("z"), 0.0);177 candidate->PositionError.SetXYZT(0.0, 0.0, clusterIDToDouble.at (cluster->first).at("ez"), 0.0);175 candidate->Position.SetXYZT(0.0, 0.0, clusterIDToDouble.at(cluster->first).at("z"), 0.0); 176 candidate->PositionError.SetXYZT(0.0, 0.0, clusterIDToDouble.at(cluster->first).at("ez"), 0.0); 178 177 179 178 fVertexOutputArray->Add(candidate); … … 183 182 //------------------------------------------------------------------------------ 184 183 185 void VertexFinder::createSeeds 184 void VertexFinder::createSeeds() 186 185 { 187 186 Candidate *candidate; … … 190 189 // Loop over all tracks, initializing some variables. 191 190 fItInputArray->Reset(); 192 while((candidate = static_cast<Candidate*>(fItInputArray->Next()))) 193 { 194 if (candidate->Momentum.Pt () < fMinPT || fabs (candidate->Momentum.Eta ()) > fMaxEta) 195 continue; 196 197 trackIDToDouble[candidate->GetUniqueID ()]["pt"] = candidate->Momentum.Pt (); 198 trackIDToDouble[candidate->GetUniqueID ()]["ept"] = candidate->ErrorPT ? candidate->ErrorPT : 1.0e-15;; 199 trackIDToDouble[candidate->GetUniqueID ()]["eta"] = candidate->Momentum.Eta (); 200 201 trackIDToDouble[candidate->GetUniqueID ()]["z"] = candidate->DZ; 202 trackIDToDouble[candidate->GetUniqueID ()]["ez"] = candidate->ErrorDZ ? candidate->ErrorDZ : 1.0e-15; 203 204 trackIDToInt[candidate->GetUniqueID ()]["clusterIndex"] = -1; 205 trackIDToInt[candidate->GetUniqueID ()]["interactionIndex"] = candidate->IsPU; 206 207 trackIDToBool[candidate->GetUniqueID ()]["claimed"] = false; 208 209 trackPT.push_back (make_pair (candidate->GetUniqueID (), candidate->Momentum.Pt ())); 210 } 191 while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) 192 { 193 if(candidate->Momentum.Pt() < fMinPT || fabs(candidate->Momentum.Eta()) > fMaxEta) 194 continue; 195 196 trackIDToDouble[candidate->GetUniqueID()]["pt"] = candidate->Momentum.Pt(); 197 trackIDToDouble[candidate->GetUniqueID()]["ept"] = candidate->ErrorPT ? candidate->ErrorPT : 1.0e-15; 198 ; 199 trackIDToDouble[candidate->GetUniqueID()]["eta"] = candidate->Momentum.Eta(); 200 201 trackIDToDouble[candidate->GetUniqueID()]["z"] = candidate->DZ; 202 trackIDToDouble[candidate->GetUniqueID()]["ez"] = candidate->ErrorDZ ? candidate->ErrorDZ : 1.0e-15; 203 204 trackIDToInt[candidate->GetUniqueID()]["clusterIndex"] = -1; 205 trackIDToInt[candidate->GetUniqueID()]["interactionIndex"] = candidate->IsPU; 206 207 trackIDToBool[candidate->GetUniqueID()]["claimed"] = false; 208 209 trackPT.push_back(make_pair(candidate->GetUniqueID(), candidate->Momentum.Pt())); 210 } 211 211 212 212 // Sort tracks by pt and leave only the SeedMinPT highest pt ones in the 213 213 // trackPT vector. 214 sort (trackPT.begin (), trackPT.end(), secondDescending);215 for (vector<pair<UInt_t, Double_t> >::const_iterator track = trackPT.begin (); track != trackPT.end(); track++, maxSeeds++)216 217 if(track->second < fSeedMinPT)218 219 214 sort(trackPT.begin(), trackPT.end(), secondDescending); 215 for(vector<pair<UInt_t, Double_t> >::const_iterator track = trackPT.begin(); track != trackPT.end(); track++, maxSeeds++) 216 { 217 if(track->second < fSeedMinPT) 218 break; 219 } 220 220 // If there are no tracks with pt above MinSeedPT, create just one seed from 221 221 // the highest pt track. 222 if 222 if(!maxSeeds) 223 223 maxSeeds++; 224 if (trackPT.size() > maxSeeds)225 226 trackPT.erase (trackPT.begin () + maxSeeds, trackPT.end());227 224 if(trackPT.size() > maxSeeds) 225 { 226 trackPT.erase(trackPT.begin() + maxSeeds, trackPT.end()); 227 } 228 228 229 229 // Create the seeds from the SeedMinPT highest pt tracks. 230 for (vector<pair<UInt_t, Double_t> >::const_iterator track = trackPT.begin (); track != trackPT.end(); track++, clusterIndex++)231 232 addTrackToCluster(track->first, clusterIndex);233 clusterSumPT2.push_back (make_pair(clusterIndex, track->second * track->second));234 235 } 236 237 //------------------------------------------------------------------------------ 238 239 void VertexFinder::growCluster 230 for(vector<pair<UInt_t, Double_t> >::const_iterator track = trackPT.begin(); track != trackPT.end(); track++, clusterIndex++) 231 { 232 addTrackToCluster(track->first, clusterIndex); 233 clusterSumPT2.push_back(make_pair(clusterIndex, track->second * track->second)); 234 } 235 } 236 237 //------------------------------------------------------------------------------ 238 239 void VertexFinder::growCluster(const UInt_t clusterIndex) 240 240 { 241 241 Bool_t done = false; … … 244 244 Double_t nearestDistance; 245 245 vector<UInt_t> nearTracks; 246 nearTracks.clear 246 nearTracks.clear(); 247 247 248 248 // Grow the cluster until there are no more tracks within Sigma standard 249 249 // deviations of the cluster. 250 while (!done) 251 { 252 done = true; 253 nearestID = 0; 254 nearestDistance = -1.0; 255 256 // These two loops are for finding the nearest track to the cluster. The 257 // first time, the ID of each track within 10*Sigma of the cluster is 258 // saved in the nearTracks vector; subsequently, to save time, only the 259 // tracks in this vector are checked. 260 if (!nearTracks.size ()) 250 while(!done) 251 { 252 done = true; 253 nearestID = 0; 254 nearestDistance = -1.0; 255 256 // These two loops are for finding the nearest track to the cluster. The 257 // first time, the ID of each track within 10*Sigma of the cluster is 258 // saved in the nearTracks vector; subsequently, to save time, only the 259 // tracks in this vector are checked. 260 if(!nearTracks.size()) 261 { 262 263 for(map<UInt_t, map<string, Double_t> >::const_iterator track = trackIDToDouble.begin(); track != trackIDToDouble.end(); track++) 264 { 265 if(trackIDToBool.at(track->first).at("claimed") || trackIDToInt.at(track->first).at("clusterIndex") == (Int_t)clusterIndex) 266 continue; 267 268 Double_t distance = fabs(clusterIDToDouble.at(clusterIndex).at("z") - track->second.at("z")) / hypot(clusterIDToDouble.at(clusterIndex).at("ez"), track->second.at("ez")); 269 if(nearestDistance < 0.0 || distance < nearestDistance) 261 270 { 262 263 for (map<UInt_t, map<string, Double_t> >::const_iterator track = trackIDToDouble.begin (); track != trackIDToDouble.end (); track++) 264 { 265 if (trackIDToBool.at (track->first).at ("claimed") || trackIDToInt.at (track->first).at ("clusterIndex") == (Int_t) clusterIndex) 266 continue; 267 268 Double_t distance = fabs (clusterIDToDouble.at (clusterIndex).at ("z") - track->second.at ("z")) / hypot (clusterIDToDouble.at (clusterIndex).at ("ez"), track->second.at ("ez")); 269 if (nearestDistance < 0.0 || distance < nearestDistance) 270 { 271 nearestID = track->first; 272 nearestDistance = distance; 273 } 274 if (distance < 10.0 * fSigma) 275 nearTracks.push_back (track->first); 276 } 271 nearestID = track->first; 272 nearestDistance = distance; 277 273 } 278 279 else 274 if(distance < 10.0 * fSigma) 275 nearTracks.push_back(track->first); 276 } 277 } 278 279 else 280 { 281 for(vector<UInt_t>::const_iterator track = nearTracks.begin(); track != nearTracks.end(); track++) 282 { 283 if(trackIDToBool.at(*track).at("claimed") || trackIDToInt.at(*track).at("clusterIndex") == (Int_t)clusterIndex) 284 continue; 285 Double_t distance = fabs(clusterIDToDouble.at(clusterIndex).at("z") - trackIDToDouble.at(*track).at("z")) / hypot(clusterIDToDouble.at(clusterIndex).at("ez"), trackIDToDouble.at(*track).at("ez")); 286 if(nearestDistance < 0.0 || distance < nearestDistance) 280 287 { 281 for (vector<UInt_t>::const_iterator track = nearTracks.begin (); track != nearTracks.end (); track++) 282 { 283 if (trackIDToBool.at (*track).at ("claimed") || trackIDToInt.at (*track).at ("clusterIndex") == (Int_t) clusterIndex) 284 continue; 285 Double_t distance = fabs (clusterIDToDouble.at (clusterIndex).at ("z") - trackIDToDouble.at (*track).at ("z")) / hypot (clusterIDToDouble.at (clusterIndex).at ("ez"), trackIDToDouble.at (*track).at ("ez")); 286 if (nearestDistance < 0.0 || distance < nearestDistance) 287 { 288 nearestID = *track; 289 nearestDistance = distance; 290 } 291 } 288 nearestID = *track; 289 nearestDistance = distance; 292 290 } 293 294 // If no tracks within Sigma of the cluster were found, stop growing. 295 done = nearestDistance > fSigma || nearestDistance < 0.0; 296 if (done) 297 { 298 continue; 299 } 300 301 // Add the nearest track within Sigma to the cluster. If it already 302 // belonged to another cluster, remove it from that cluster first. 303 if (nearestDistance < fSigma) 304 { 305 oldClusterIndex = trackIDToInt.at (nearestID).at ("clusterIndex"); 306 if (oldClusterIndex >= 0) 307 removeTrackFromCluster (nearestID, oldClusterIndex); 308 309 trackIDToBool[nearestID]["claimed"] = true; 310 addTrackToCluster (nearestID, clusterIndex); 311 } 312 } 313 } 314 315 //------------------------------------------------------------------------------ 316 317 Double_t VertexFinder::weight (const UInt_t trackID) 318 { 319 return ((trackIDToDouble.at (trackID).at ("pt") / (trackIDToDouble.at (trackID).at ("ept") * trackIDToDouble.at (trackID).at ("ez"))) * (trackIDToDouble.at (trackID).at ("pt") / (trackIDToDouble.at (trackID).at ("ept") * trackIDToDouble.at (trackID).at ("ez")))); 320 } 321 322 //------------------------------------------------------------------------------ 323 324 void VertexFinder::removeTrackFromCluster (const UInt_t trackID, const UInt_t clusterID) 325 { 326 Double_t wz = weight (trackID); 291 } 292 } 293 294 // If no tracks within Sigma of the cluster were found, stop growing. 295 done = nearestDistance > fSigma || nearestDistance < 0.0; 296 if(done) 297 { 298 continue; 299 } 300 301 // Add the nearest track within Sigma to the cluster. If it already 302 // belonged to another cluster, remove it from that cluster first. 303 if(nearestDistance < fSigma) 304 { 305 oldClusterIndex = trackIDToInt.at(nearestID).at("clusterIndex"); 306 if(oldClusterIndex >= 0) 307 removeTrackFromCluster(nearestID, oldClusterIndex); 308 309 trackIDToBool[nearestID]["claimed"] = true; 310 addTrackToCluster(nearestID, clusterIndex); 311 } 312 } 313 } 314 315 //------------------------------------------------------------------------------ 316 317 Double_t VertexFinder::weight(const UInt_t trackID) 318 { 319 return ((trackIDToDouble.at(trackID).at("pt") / (trackIDToDouble.at(trackID).at("ept") * trackIDToDouble.at(trackID).at("ez"))) * (trackIDToDouble.at(trackID).at("pt") / (trackIDToDouble.at(trackID).at("ept") * trackIDToDouble.at(trackID).at("ez")))); 320 } 321 322 //------------------------------------------------------------------------------ 323 324 void VertexFinder::removeTrackFromCluster(const UInt_t trackID, const UInt_t clusterID) 325 { 326 Double_t wz = weight(trackID); 327 327 328 328 trackIDToInt[trackID]["clusterIndex"] = -1; 329 329 clusterIDToInt[clusterID]["ndf"]--; 330 330 331 clusterIDToDouble[clusterID]["sumZ"] -= wz * trackIDToDouble.at (trackID).at("z");332 clusterIDToDouble[clusterID]["errorSumZ"] -= wz * trackIDToDouble.at (trackID).at ("ez") * trackIDToDouble.at (trackID).at("ez");331 clusterIDToDouble[clusterID]["sumZ"] -= wz * trackIDToDouble.at(trackID).at("z"); 332 clusterIDToDouble[clusterID]["errorSumZ"] -= wz * trackIDToDouble.at(trackID).at("ez") * trackIDToDouble.at(trackID).at("ez"); 333 333 clusterIDToDouble[clusterID]["sumOfWeightsZ"] -= wz; 334 clusterIDToDouble[clusterID]["z"] = clusterIDToDouble.at (clusterID).at ("sumZ") / clusterIDToDouble.at (clusterID).at("sumOfWeightsZ");335 clusterIDToDouble[clusterID]["ez"] = sqrt ((1.0 / clusterIDToInt.at (clusterID).at ("ndf")) * (clusterIDToDouble.at (clusterID).at ("errorSumZ") / clusterIDToDouble.at (clusterID).at("sumOfWeightsZ")));336 clusterIDToDouble[clusterID]["sumPT2"] -= trackIDToDouble.at (trackID).at ("pt") * trackIDToDouble.at (trackID).at("pt");337 } 338 339 //------------------------------------------------------------------------------ 340 341 void VertexFinder::addTrackToCluster 342 { 343 Double_t wz = weight 344 345 if (!clusterIDToInt.count(clusterID))346 347 348 349 350 351 352 353 334 clusterIDToDouble[clusterID]["z"] = clusterIDToDouble.at(clusterID).at("sumZ") / clusterIDToDouble.at(clusterID).at("sumOfWeightsZ"); 335 clusterIDToDouble[clusterID]["ez"] = sqrt((1.0 / clusterIDToInt.at(clusterID).at("ndf")) * (clusterIDToDouble.at(clusterID).at("errorSumZ") / clusterIDToDouble.at(clusterID).at("sumOfWeightsZ"))); 336 clusterIDToDouble[clusterID]["sumPT2"] -= trackIDToDouble.at(trackID).at("pt") * trackIDToDouble.at(trackID).at("pt"); 337 } 338 339 //------------------------------------------------------------------------------ 340 341 void VertexFinder::addTrackToCluster(const UInt_t trackID, const UInt_t clusterID) 342 { 343 Double_t wz = weight(trackID); 344 345 if(!clusterIDToInt.count(clusterID)) 346 { 347 clusterIDToInt[clusterID]["ndf"] = 0; 348 clusterIDToInt[clusterID]["seed"] = trackID; 349 clusterIDToDouble[clusterID]["sumZ"] = 0.0; 350 clusterIDToDouble[clusterID]["errorSumZ"] = 0.0; 351 clusterIDToDouble[clusterID]["sumOfWeightsZ"] = 0.0; 352 clusterIDToDouble[clusterID]["sumPT2"] = 0.0; 353 } 354 354 355 355 trackIDToInt[trackID]["clusterIndex"] = clusterID; 356 356 clusterIDToInt[clusterID]["ndf"]++; 357 357 358 clusterIDToDouble[clusterID]["sumZ"] += wz * trackIDToDouble.at (trackID).at("z");359 clusterIDToDouble[clusterID]["errorSumZ"] += wz * trackIDToDouble.at (trackID).at ("ez") * trackIDToDouble.at (trackID).at("ez");358 clusterIDToDouble[clusterID]["sumZ"] += wz * trackIDToDouble.at(trackID).at("z"); 359 clusterIDToDouble[clusterID]["errorSumZ"] += wz * trackIDToDouble.at(trackID).at("ez") * trackIDToDouble.at(trackID).at("ez"); 360 360 clusterIDToDouble[clusterID]["sumOfWeightsZ"] += wz; 361 clusterIDToDouble[clusterID]["z"] = clusterIDToDouble.at (clusterID).at ("sumZ") / clusterIDToDouble.at (clusterID).at("sumOfWeightsZ");362 clusterIDToDouble[clusterID]["ez"] = sqrt ((1.0 / clusterIDToInt.at (clusterID).at ("ndf")) * (clusterIDToDouble.at (clusterID).at ("errorSumZ") / clusterIDToDouble.at (clusterID).at("sumOfWeightsZ")));363 clusterIDToDouble[clusterID]["sumPT2"] += trackIDToDouble.at (trackID).at ("pt") * trackIDToDouble.at (trackID).at("pt");364 } 365 366 //------------------------------------------------------------------------------ 361 clusterIDToDouble[clusterID]["z"] = clusterIDToDouble.at(clusterID).at("sumZ") / clusterIDToDouble.at(clusterID).at("sumOfWeightsZ"); 362 clusterIDToDouble[clusterID]["ez"] = sqrt((1.0 / clusterIDToInt.at(clusterID).at("ndf")) * (clusterIDToDouble.at(clusterID).at("errorSumZ") / clusterIDToDouble.at(clusterID).at("sumOfWeightsZ"))); 363 clusterIDToDouble[clusterID]["sumPT2"] += trackIDToDouble.at(trackID).at("pt") * trackIDToDouble.at(trackID).at("pt"); 364 } 365 366 //------------------------------------------------------------------------------
Note:
See TracChangeset
for help on using the changeset viewer.