Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/VertexFinder.cc

    r77e9ae1 r95b4e9f  
    77 */
    88
     9
    910#include "modules/VertexFinder.h"
    1011#include "classes/DelphesClasses.h"
    … …  
    1314#include "classes/DelphesPileUpReader.h"
    1415
     16#include "ExRootAnalysis/ExRootResult.h"
     17#include "ExRootAnalysis/ExRootFilter.h"
    1518#include "ExRootAnalysis/ExRootClassifier.h"
    16 #include "ExRootAnalysis/ExRootFilter.h"
    17 #include "ExRootAnalysis/ExRootResult.h"
    18 
     19
     20#include "TMath.h"
     21#include "TString.h"
     22#include "TFormula.h"
     23#include "TRandom3.h"
     24#include "TObjArray.h"
    1925#include "TDatabasePDG.h"
    20 #include "TFormula.h"
    2126#include "TLorentzVector.h"
    22 #include "TMath.h"
    2327#include "TMatrixT.h"
    24 #include "TObjArray.h"
    25 #include "TRandom3.h"
    26 #include "TString.h"
    2728#include "TVector3.h"
    2829
     30#include <utility>
    2931#include <algorithm>
     32#include <stdexcept>
    3033#include <iostream>
     34#include <vector>
    3135#include <map>
    32 #include <stdexcept>
    3336#include <string>
    34 #include <utility>
    35 #include <vector>
    3637
    3738using namespace std;
    3839
    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;
     40static const Double_t mm  = 1.;
     41static const Double_t m = 1000.*mm;
     42static const Double_t ns  = 1.;
     43static const Double_t s = 1.e+9 *ns;
     44static const Double_t c_light   = 2.99792458e+8 * m/s;
    4445
    4546//------------------------------------------------------------------------------
    … …  
    8384//------------------------------------------------------------------------------
    8485
    85 static Bool_t secondAscending(pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
     86static Bool_t secondAscending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
    8687{
    8788  return (pair0.second < pair1.second);
    8889}
    8990
    90 static Bool_t secondDescending(pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
     91static Bool_t secondDescending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
    9192{
    9293  return (pair0.second > pair1.second);
    … …  
    100101
    101102  // Clear the track and cluster maps before starting
    102   trackIDToDouble.clear();
    103   trackIDToInt.clear();
    104   trackIDToBool.clear();
    105   clusterIDToDouble.clear();
    106   clusterIDToInt.clear();
    107   clusterIDToBool.clear();
    108   trackPT.clear();
    109   clusterSumPT2.clear();
     103  trackIDToDouble.clear ();
     104  trackIDToInt.clear ();
     105  trackIDToBool.clear ();
     106  clusterIDToDouble.clear ();
     107  clusterIDToInt.clear ();
     108  clusterIDToBool.clear ();
     109  trackPT.clear ();
     110  clusterSumPT2.clear ();
    110111
    111112  // Create the initial cluster seeds
    112   createSeeds();
     113  createSeeds ();
    113114
    114115  // In order of descending seed pt, grow each cluster. If a cluster ends up with
    115116  // fewer than MinNDF tracks, release the tracks for other clusters to claim.
    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   }
     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    }
    143144
    144145  // Add tracks to the output array after updating their ClusterIndex.
    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   }
     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    }
    153154
    154155  // Add clusters with at least MinNDF tracks to the output array in order of
    155156  // descending sum(pt**2).
    156   clusterSumPT2.clear();
    157   for(map<UInt_t, map<string, Int_t> >::const_iterator cluster = clusterIDToInt.begin(); cluster != clusterIDToInt.end(); cluster++)
     157  clusterSumPT2.clear ();
     158  for (map<UInt_t, map<string, Int_t> >::const_iterator cluster = clusterIDToInt.begin (); cluster != clusterIDToInt.end (); cluster++)
    158159  {
    159 
    160     if(cluster->second.at("ndf") < fMinNDF)
     160   
     161    if (cluster->second.at ("ndf") < fMinNDF)
    161162      continue;
    162     clusterSumPT2.push_back(make_pair(cluster->first, clusterIDToDouble.at(cluster->first).at("sumPT2")));
     163    clusterSumPT2.push_back (make_pair (cluster->first, clusterIDToDouble.at (cluster->first).at ("sumPT2")));
    163164  }
    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++)
     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++)
    167168  {
    168169    DelphesFactory *factory = GetFactory();
    … …  
    170171
    171172    candidate->ClusterIndex = cluster->first;
    172     candidate->ClusterNDF = clusterIDToInt.at(cluster->first).at("ndf");
     173    candidate->ClusterNDF = clusterIDToInt.at (cluster->first).at ("ndf");
    173174    candidate->ClusterSigma = fSigma;
    174175    candidate->SumPT2 = cluster->second;
    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);
     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);
    177178
    178179    fVertexOutputArray->Add(candidate);
    … …  
    182183//------------------------------------------------------------------------------
    183184
    184 void VertexFinder::createSeeds()
     185void VertexFinder::createSeeds ()
    185186{
    186187  Candidate *candidate;
    … …  
    189190  // Loop over all tracks, initializing some variables.
    190191  fItInputArray->Reset();
    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   }
     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    }
    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())
    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)
     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)
     270                {
     271                  nearestID = track->first;
     272                  nearestDistance = distance;
     273                }
     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)
     287                {
     288                  nearestID = *track;
     289                  nearestDistance = distance;
     290                }
     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        {
    266298          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       }
    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)
    287         {
    288           nearestID = *track;
    289           nearestDistance = distance;
    290         }
    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);
     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.