Fork me on GitHub

Changeset 341014c in git for modules/VertexFinder.cc


Ignore:
Timestamp:
Feb 12, 2019, 9:29:17 PM (6 years ago)
Author:
Pavel Demin <pavel-demin@…>
Branches:
ImprovedOutputFile, Timing, llp, master
Children:
6455202
Parents:
45e58be
Message:

apply .clang-format to all .h, .cc and .cpp files

File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/VertexFinder.cc

    r45e58be r341014c  
    77 */
    88
    9 
    109#include "modules/VertexFinder.h"
    1110#include "classes/DelphesClasses.h"
    … …  
    1413#include "classes/DelphesPileUpReader.h"
    1514
     15#include "ExRootAnalysis/ExRootClassifier.h"
     16#include "ExRootAnalysis/ExRootFilter.h"
    1617#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"
    2022#include "TMath.h"
     23#include "TMatrixT.h"
     24#include "TObjArray.h"
     25#include "TRandom3.h"
    2126#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"
    2827#include "TVector3.h"
    2928
     29#include <algorithm>
     30#include <iostream>
     31#include <map>
     32#include <stdexcept>
     33#include <string>
    3034#include <utility>
    31 #include <algorithm>
    32 #include <stdexcept>
    33 #include <iostream>
    3435#include <vector>
    35 #include <map>
    36 #include <string>
    3736
    3837using namespace std;
    3938
    40 static const Double_t mm  = 1.;
    41 static const Double_t m = 1000.*mm;
    42 static const Double_t ns  = 1.;
    43 static const Double_t s = 1.e+9 *ns;
    44 static const Double_t c_light   = 2.99792458e+8 * m/s;
     39static const Double_t mm = 1.;
     40static const Double_t m = 1000. * mm;
     41static const Double_t ns = 1.;
     42static const Double_t s = 1.e+9 * ns;
     43static const Double_t c_light = 2.99792458e+8 * m / s;
    4544
    4645//------------------------------------------------------------------------------
    … …  
    8483//------------------------------------------------------------------------------
    8584
    86 static Bool_t secondAscending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
     85static Bool_t secondAscending(pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
    8786{
    8887  return (pair0.second < pair1.second);
    8988}
    9089
    91 static Bool_t secondDescending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
     90static Bool_t secondDescending(pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
    9291{
    9392  return (pair0.second > pair1.second);
    … …  
    101100
    102101  // 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();
    111110
    112111  // Create the initial cluster seeds
    113   createSeeds ();
     112  createSeeds();
    114113
    115114  // In order of descending seed pt, grow each cluster. If a cluster ends up with
    116115  // 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       // Skip the cluster if it no longer has any tracks
    121       if (!clusterIDToInt.at (cluster->first).at ("ndf"))
    122         continue;
    123        
    124       // Grow the cluster if GrowSeeds is true
    125       if (fGrowSeeds)
    126         growCluster (cluster->first);
    127 
    128       // If the cluster still has fewer than MinNDF tracks, release the tracks;
    129       // otherwise, mark the seed track as claimed
    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                 continue;
    137               track->second["clusterIndex"] = -1;
    138               trackIDToBool[track->first]["claimed"] = false;
    139             }
    140         }
    141       else
    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  }
    144143
    145144  // 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         continue;
    151       candidate->ClusterIndex = trackIDToInt.at (candidate->GetUniqueID ()).at ("clusterIndex");
    152       fOutputArray->Add(candidate);
    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  }
    154153
    155154  // Add clusters with at least MinNDF tracks to the output array in order of
    156155  // 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++)
    168167  {
    169168    DelphesFactory *factory = GetFactory();
    … …  
    171170
    172171    candidate->ClusterIndex = cluster->first;
    173     candidate->ClusterNDF = clusterIDToInt.at (cluster->first).at ("ndf");
     172    candidate->ClusterNDF = clusterIDToInt.at(cluster->first).at("ndf");
    174173    candidate->ClusterSigma = fSigma;
    175174    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);
    178177
    179178    fVertexOutputArray->Add(candidate);
    … …  
    183182//------------------------------------------------------------------------------
    184183
    185 void VertexFinder::createSeeds ()
     184void VertexFinder::createSeeds()
    186185{
    187186  Candidate *candidate;
    … …  
    190189  // Loop over all tracks, initializing some variables.
    191190  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  }
    211211
    212212  // Sort tracks by pt and leave only the SeedMinPT highest pt ones in the
    213213  // 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         break;
    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  }
    220220  // If there are no tracks with pt above MinSeedPT, create just one seed from
    221221  // the highest pt track.
    222   if (!maxSeeds)
     222  if(!maxSeeds)
    223223    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  }
    228228
    229229  // 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 (const UInt_t clusterIndex)
     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
     239void VertexFinder::growCluster(const UInt_t clusterIndex)
    240240{
    241241  Bool_t done = false;
    … …  
    244244  Double_t nearestDistance;
    245245  vector<UInt_t> nearTracks;
    246   nearTracks.clear ();
     246  nearTracks.clear();
    247247
    248248  // Grow the cluster until there are no more tracks within Sigma standard
    249249  // 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)
    261270        {
    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;
    277273        }
    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)
    280287        {
    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;
    292290        }
    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
     317Double_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
     324void VertexFinder::removeTrackFromCluster(const UInt_t trackID, const UInt_t clusterID)
     325{
     326  Double_t wz = weight(trackID);
    327327
    328328  trackIDToInt[trackID]["clusterIndex"] = -1;
    329329  clusterIDToInt[clusterID]["ndf"]--;
    330330
    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");
    333333  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 (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     }
     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
     341void 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  }
    354354
    355355  trackIDToInt[trackID]["clusterIndex"] = clusterID;
    356356  clusterIDToInt[clusterID]["ndf"]++;
    357357
    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");
    360360  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.