Fork me on GitHub

Changeset 84a097e in git


Ignore:
Timestamp:
Dec 2, 2019, 6:55:22 PM (5 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
Timing
Children:
2ad823e, 6fc566b
Parents:
c614dd7
Message:

first iteration of adding timing in Delphes

Files:
8 edited

Legend:

Unmodified
Added
Removed

  • classes/DelphesClasses.cc

    rc614dd7 r84a097e  
    147147  CtgTheta(0), ErrorCtgTheta(0),
    148148  Phi(0), ErrorPhi(0),
    149   Xd(0), Yd(0), Zd(0),
     149  Xd(0), Yd(0), Zd(0), Td(0),
     150  VertexingWeight(0),
    150151  TrackResolution(0),
    151152  NCharged(0),
    … …  
    318319  object.Phi = Phi;
    319320  object.ErrorPhi = ErrorPhi;
     321  object.Td = Td;
    320322  object.Xd = Xd;
    321323  object.Yd = Yd;
    322324  object.Zd = Zd;
     325  object.VertexingWeight = Zd;
    323326  object.TrackResolution = TrackResolution;
    324327  object.NCharged = NCharged;
    … …  
    447450  Phi = 0.0;
    448451  ErrorPhi = 0.0;
     452  Td = 0.0;
    449453  Xd = 0.0;
    450454  Yd = 0.0;
    451455  Zd = 0.0;
     456  VertexingWeight = 0.0;
    452457  TrackResolution = 0.0;
    453458  NCharged = 0;

  • classes/DelphesClasses.h

    rc614dd7 r84a097e  
    456456  Float_t ZOuter; // track position (z component) at the tracker edge
    457457
     458  Float_t Td; // T coordinate of point of closest approach to vertex
    458459  Float_t Xd; // X coordinate of point of closest approach to vertex
    459460  Float_t Yd; // Y coordinate of point of closest approach to vertex
    460461  Float_t Zd; // Z coordinate of point of closest approach to vertex
     462
     463  Float_t TOFreco; // reconstructed time of flight
     464  Float_t TOFgen; // gen time of flight
    461465
    462466  Float_t L; // track path length
    … …  
    482486  TLorentzVector P4() const;
    483487
    484   ClassDef(Track, 3)
     488  ClassDef(Track, 4)
    485489};
    486490
    … …  
    544548  Float_t ZOuter; // track position (z component) at the tracker edge
    545549
     550  Float_t Td; // X coordinate of point of closest approach to vertex
    546551  Float_t Xd; // X coordinate of point of closest approach to vertex
    547552  Float_t Yd; // Y coordinate of point of closest approach to vertex
    … …  
    664669  Float_t ErrorPhi;
    665670
     671  Float_t Td;
    666672  Float_t Xd;
    667673  Float_t Yd;
    668674  Float_t Zd;
    669675
     676  // Vertexing variables
     677  Float_t VertexingWeight;
     678
    670679  // tracking resolution
    671 
    672680  Float_t TrackResolution;
    673681

  • modules/MomentumSmearing.cc

    rc614dd7 r84a097e  
    9595{
    9696  Candidate *candidate, *mother;
    97   Double_t pt, eta, phi, e, res;
     97  Double_t pt, eta, phi, e, m, res;
    9898
    9999  fItInputArray->Reset();
    … …  
    106106    pt = candidateMomentum.Pt();
    107107    e = candidateMomentum.E();
     108    m = candidateMomentum.M();
     109
    108110    res = fFormula->Eval(pt, eta, phi, e);
    109 
    110     // apply smearing formula
    111     //pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
    112 
    113111    res = (res > 1.0) ? 1.0 : res;
    114112
    115113    pt = LogNormal(pt, res * pt);
    116 
    117     //if(pt <= 0.0) continue;
    118114
    119115    mother = candidate;
    … …  
    121117    eta = candidateMomentum.Eta();
    122118    phi = candidateMomentum.Phi();
    123     candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
    124     //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
     119    candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
     120    candidate->ErrorPT = res*pt;
     121    candidate->PT = pt;
     122   
    125123    candidate->TrackResolution = res;
    126124    candidate->AddCandidate(mother);

  • modules/TrackSmearing.cc

    rc614dd7 r84a097e  
    158158  TLorentzVector beamSpotPosition;
    159159  Candidate *candidate, *mother;
    160   Double_t pt, eta, d0, d0Error, trueD0, dz, dzError, trueDZ, p, pError, trueP, ctgTheta, ctgThetaError, trueCtgTheta, phi, phiError, truePhi;
     160  Double_t pt, eta, d0, d0Error, trueD0, dz, dzError, trueDZ, p, pError, trueP, ctgTheta, ctgThetaError, trueCtgTheta, phi, phiError, truePhi, m;
    161161  Double_t x, y, z, t, px, py, pz, theta;
    162162  Double_t q, r;
    … …  
    328328
    329329    theta = TMath::ACos(ctgTheta / TMath::Sqrt(1.0 + ctgTheta * ctgTheta));
     330    m = candidate->Momentum.M();
     331   
    330332    candidate->Momentum.SetPx(p * TMath::Cos(phi) * TMath::Sin(theta));
    331333    candidate->Momentum.SetPy(p * TMath::Sin(phi) * TMath::Sin(theta));
    332334    candidate->Momentum.SetPz(p * TMath::Cos(theta));
    333     candidate->Momentum.SetE(candidate->Momentum.Pt() * TMath::CosH(eta));
     335    //candidate->Momentum.SetE(candidate->Momentum.Pt() * TMath::CosH(eta));
     336    candidate->Momentum.SetE( TMath::Sqrt(m*m + p*p) );
    334337    candidate->PT = candidate->Momentum.Pt();
    335338
    … …  
    382385    candidate->Yd = yd * 1.0E3;
    383386    candidate->Zd = zd * 1.0E3;
    384 
     387    candidate->Td = -9999*c_light*1E-3;
     388   
    385389    if(fApplyToPileUp || !candidate->IsPU)
    386390    {

  • modules/TreeWriter.cc

    rc614dd7 r84a097e  
    360360    entry->Yd = candidate->Yd;
    361361    entry->Zd = candidate->Zd;
     362    entry->Td = candidate->Td*1.0E-3/c_light;
     363
     364    if(candidate->ClusterIndex != -1)
     365    {
     366      entry->TOFreco = 1E-3*(candidate->Position.T() - candidate->InitialPosition.T())/c_light;
     367    }
     368    else
     369    {
     370      entry->TOFreco =-1E6;
     371    }
    362372
    363373    const TLorentzVector &momentum = candidate->Momentum;
    … …  
    380390
    381391    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
     392    entry->TOFgen  = 1E-3*particle->L/(c_light*particle->Momentum.P()/particle->Momentum.E());
     393   
    382394    const TLorentzVector &initialPosition = particle->Position;
    383395

  • modules/VertexFinderDA4D.cc

    rc614dd7 r84a097e  
    33 *  Cluster vertices from tracks using deterministic annealing and timing information
    44 *
    5  *  \authors M. Selvaggi, L. Gray
     5 *  \authors O. Cerri
    66 *
    77 */
     8
    89
    910#include "modules/VertexFinderDA4D.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"
     28#include "TLatex.h"
     29#include "TVector3.h"
     30
     31#include "TAxis.h"
     32#include "TGraphErrors.h"
     33#include "TCanvas.h"
    2634#include "TString.h"
    27 #include "TVector3.h"
    28 
     35#include "TLegend.h"
     36#include "TFile.h"
     37#include "TColor.h"
     38#include "TLegend.h"
     39
     40#include <utility>
    2941#include <algorithm>
     42#include <stdexcept>
    3043#include <iostream>
    31 #include <stdexcept>
    32 #include <utility>
    3344#include <vector>
    3445
    3546using namespace std;
    3647
    37 static const Double_t mm = 1.;
    38 static const Double_t m = 1000. * mm;
    39 static const Double_t ns = 1.;
    40 static const Double_t s = 1.e+9 * ns;
    41 static const Double_t c_light = 2.99792458e+8 * m / s;
    42 
    43 struct track_t
    44 {
    45   double z; // z-coordinate at point of closest approach to the beamline
    46   double t; // t-coordinate at point of closest approach to the beamline
    47   double dz2; // square of the error of z(pca)
    48   double dtz; // covariance of z-t
    49   double dt2; // square of the error of t(pca)
    50   Candidate *tt; // a pointer to the Candidate Track
    51   double Z; // Z[i]   for DA clustering
    52   double pi; // track weight
    53   double pt;
    54   double eta;
    55   double phi;
    56 };
    57 
    58 struct vertex_t
    59 {
    60   double z;
    61   double t;
    62   double pk; // vertex weight for "constrained" clustering
    63   // --- temporary numbers, used during update
    64   double ei;
    65   double sw;
    66   double swz;
    67   double swt;
    68   double se;
    69   // ---for Tc
    70   double swE;
    71   double Tc;
    72 };
    73 
    74 static bool split(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y);
    75 static double update1(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y);
    76 static double update2(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y, double &rho0, const double dzCutOff);
    77 static void dump(const double beta, const std::vector<vertex_t> &y, const std::vector<track_t> &tks);
    78 static bool merge(std::vector<vertex_t> &);
    79 static bool merge(std::vector<vertex_t> &, double &);
    80 static bool purge(std::vector<vertex_t> &, std::vector<track_t> &, double &, const double, const double);
    81 static void splitAll(std::vector<vertex_t> &y);
    82 static double beta0(const double betamax, std::vector<track_t> &tks, std::vector<vertex_t> &y, const double coolingFactor);
    83 static double Eik(const track_t &t, const vertex_t &k);
    84 
    85 static bool recTrackLessZ1(const track_t &tk1, const track_t &tk2)
    86 {
    87   return tk1.z < tk2.z;
    88 }
    89 
    90 using namespace std;
     48namespace vtx_DAZT
     49{
     50  static const Double_t c_light = 2.99792458e+8; // [m/s]
     51}
     52using namespace vtx_DAZT;
    9153
    9254//------------------------------------------------------------------------------
    9355
    94 VertexFinderDA4D::VertexFinderDA4D() :
    95   fVerbose(0), fMinPT(0), fVertexSpaceSize(0), fVertexTimeSize(0),
    96   fUseTc(0), fBetaMax(0), fBetaStop(0), fCoolingFactor(0),
    97   fMaxIterations(0), fDzCutOff(0), fD0CutOff(0), fDtCutOff(0)
    98 {
     56VertexFinderDA4D::VertexFinderDA4D()
     57{
     58  fVerbose = 0;
     59  fMaxIterations = 0;
     60  fBetaMax = 0;
     61  fBetaStop = 0;
     62  fBetaPurge = 0;
     63  fVertexZSize = 0;
     64  fVertexTSize = 0;
     65  fCoolingFactor = 0;
     66  fDzCutOff = 0;
     67  fD0CutOff = 0;
     68  fDtCutOff = 0;
     69  fPtMin = 0;
     70  fPtMax = 0;
     71  fD2Merge = 0;
     72  fMuOutlayer = 0;
     73  fMinTrackProb = 0;
    9974}
    10075
    … …  
    10984void VertexFinderDA4D::Init()
    11085{
    111 
    112   fVerbose = GetBool("Verbose", 1);
    113   fMinPT = GetDouble("MinPT", 0.1);
    114   fVertexSpaceSize = GetDouble("VertexSpaceSize", 0.5); //in mm
    115   fVertexTimeSize = GetDouble("VertexTimeSize", 10E-12); //in s
    116   fUseTc = GetBool("UseTc", 1);
    117   fBetaMax = GetDouble("BetaMax ", 0.1);
    118   fBetaStop = GetDouble("BetaStop", 1.0);
    119   fCoolingFactor = GetDouble("CoolingFactor", 0.8);
    120   fMaxIterations = GetInt("MaxIterations", 100);
    121   fDzCutOff = GetDouble("DzCutOff", 40); // Adaptive Fitter uses 30 mm but that appears to be a bit tight here sometimes
    122   fD0CutOff = GetDouble("D0CutOff", 30);
    123   fDtCutOff = GetDouble("DtCutOff", 100E-12); // dummy
    124 
    125   // convert stuff in cm, ns
    126   fVertexSpaceSize /= 10.0;
    127   fVertexTimeSize *= 1E9;
    128   fDzCutOff /= 10.0; // Adaptive Fitter uses 3.0 but that appears to be a bit tight here sometimes
    129   fD0CutOff /= 10.0;
     86  fVerbose         = GetInt("Verbose", 0);
     87
     88  fMaxIterations   = GetInt("MaxIterations", 100);
     89  fMaxVertexNumber = GetInt("MaxVertexNumber", 500);
     90
     91  fBetaMax         = GetDouble("BetaMax", 1.5);
     92  fBetaPurge       = GetDouble("BetaPurge", 1.);
     93  fBetaStop        = GetDouble("BetaStop", 0.2);
     94
     95  fVertexZSize     = GetDouble("VertexZSize", 0.1); //in mm
     96  fVertexTSize     = 1E12*GetDouble("VertexTimeSize", 15E-12); //Convert from [s] to [ps]
     97
     98  fCoolingFactor   = GetDouble("CoolingFactor", 0.8); // Multiply T so to cooldown must be <1
     99
     100  fDzCutOff        = GetDouble("DzCutOff", 40);      // For the moment 3*DzCutOff is hard cut off for the considered tracks
     101  fD0CutOff        = GetDouble("D0CutOff", .5);       // d0/sigma_d0, used to compute the pi (weight) of the track
     102  fDtCutOff        = GetDouble("DtCutOff", 160);     // [ps], 3*DtCutOff is hard cut off for tracks
     103  fPtMin           = GetDouble("PtMin", 0.5);        // Minimum pt accepted for tracks
     104  fPtMax           = GetDouble("PtMax", 50);        // Maximum pt accepted for tracks
     105
     106
     107  fD2UpdateLim     = GetDouble("D2UpdateLim", .5);   // ((dz/ZSize)^2+(dt/TSize)^2)/nv limit for merging vertices
     108  fD2Merge         = GetDouble("D2Merge", 4.0);      // (dz/ZSize)^2+(dt/TSize)^2 limit for merging vertices
     109  fMuOutlayer      = GetDouble("MuOutlayer", 4);     // Outlayer rejection exponent
     110  fMinTrackProb    = GetDouble("MinTrackProb", 0.6); // Minimum probability to be assigned at a vertex
     111  fMinNTrack       = GetInt("MinNTrack", 10);        // Minimum number of tracks per vertex
     112
     113  fFigFolderPath   = GetString("DebugFigPath", ".");
    130114
    131115  fInputArray = ImportArray(GetString("InputArray", "TrackSmearing/tracks"));
    132116  fItInputArray = fInputArray->MakeIterator();
    133117
    134   fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
     118  fTrackOutputArray = ExportArray(GetString("TrackOutputArray", "tracks"));
    135119  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
     120
     121  fInputGenVtx = ImportArray(GetString("InputGenVtx", "PileUpMerger/vertices"));
     122  fItInputGenVtx = fInputGenVtx->MakeIterator();
     123
     124  if (fBetaMax < fBetaPurge)
     125  {
     126    fBetaPurge = fBetaMax;
     127    if (fVerbose)
     128    {
     129      cout << "BetaPurge set to " << fBetaPurge << endl;
     130    }
     131  }
     132
     133  if (fBetaPurge < fBetaStop)
     134  {
     135    fBetaStop = fBetaPurge;
     136    if (fVerbose)
     137    {
     138      cout << "BetaPurge set to " << fBetaPurge << endl;
     139    }
     140  }
    136141}
    137142
    … …  
    147152void VertexFinderDA4D::Process()
    148153{
    149   Candidate *candidate, *track;
    150   TObjArray *ClusterArray;
    151   ClusterArray = new TObjArray;
    152   TIterator *ItClusterArray;
    153   Int_t ivtx = 0;
    154 
    155154  fInputArray->Sort();
    156155
    157   TLorentzVector pos, mom;
     156  if (fVerbose)
     157  {
     158     cout<< endl << "      Start processing vertices with VertexFinderDA4D" << endl;
     159     cout<<" Found "<<fInputArray->GetEntriesFast()<<" input tracks"<<endl;
     160  }
     161
     162  // clusterize tracks
     163  TObjArray *ClusterArray = new TObjArray;
     164  clusterize(*ClusterArray);
     165
     166  if(fVerbose>10)
     167  {
     168    unsigned int N = fEnergy_rec.size();
     169    TGraph* gr1 = new TGraph(N, &fBeta_rec[0], &fNvtx_rec[0]);
     170    gr1->SetName("gr1");
     171    gr1->GetXaxis()->SetTitle("beta");
     172    gr1->GetYaxis()->SetTitle("# Vtx");
     173    TGraph* gr2 = new TGraph(N, &fBeta_rec[0], &fEnergy_rec[0]);
     174    gr2->SetName("gr2");
     175    gr2->GetXaxis()->SetTitle("beta");
     176    gr2->GetYaxis()->SetTitle("Total Energy");
     177    TGraph* gr3 = new TGraph(N, &fNvtx_rec[0], &fEnergy_rec[0]);
     178    gr3->SetName("gr3");
     179    gr3->GetXaxis()->SetTitle("# Vtx");
     180    gr3->GetYaxis()->SetTitle("Total Energy");
     181
     182    auto f = new TFile("~/Desktop/debug/EnergyStat.root", "recreate");
     183    gr1->Write("gr1");
     184    gr2->Write("gr2");
     185    gr3->Write("gr3");
     186
     187    f->Close();
     188  }
     189
     190  if (fVerbose){std::cout <<  " clustering returned  "<< ClusterArray->GetEntriesFast() << " clusters  from " << fInputArray->GetEntriesFast() << " input tracks" <<std::endl;}
     191
     192  // //loop over vertex candidates
     193  TIterator * ItClusterArray = ClusterArray->MakeIterator();
     194  ItClusterArray->Reset();
     195  Candidate *candidate;
     196  unsigned int k = 0;
     197  while((candidate = static_cast<Candidate*>(ItClusterArray->Next())))
     198  {
     199    if(fVerbose)
     200    {
     201     cout << Form("Cluster %d has %d tracks ", k, candidate->GetCandidates()->GetEntriesFast()) << endl;
     202    }
     203    if(candidate->ClusterNDF>0)
     204    {
     205      // Estimate the vertex resolution
     206      // loop over tracks belonging to this vertex
     207      TIter it1(candidate->GetCandidates());
     208      it1.Reset();
     209
     210      Candidate *track;
     211      double sum_Dt_2 = 0;
     212      double sum_Dz_2 = 0;
     213      double sum_wt = 0;
     214      double sum_wz = 0;
     215      while((track = static_cast<Candidate*>(it1.Next())))
     216      {
     217        double dz = candidate->Position.Z() - track->Zd;
     218        double dt = candidate->Position.T() - track->Td;
     219
     220        double wz = track->VertexingWeight/(track->ErrorDZ*track->ErrorDZ);
     221        double wt = track->VertexingWeight/(track->ErrorT*track->ErrorT);
     222
     223        sum_Dt_2 += wt*dt*dt;
     224        sum_Dz_2 += wz*dz*dz;
     225        sum_wt += wt;
     226        sum_wz += wz;
     227      }
     228
     229      double sigma_z = sqrt(sum_Dz_2/sum_wz);
     230      double sigma_t = sqrt(sum_Dt_2/sum_wt);
     231      candidate->PositionError.SetXYZT(0.0, 0.0, sigma_z , sigma_t);
     232      if(fVerbose > 3)
     233      {
     234        cout << "k: " << k << endl;
     235        cout << "Sigma z: " << sigma_z*1E3 << " um" << endl;
     236        cout << "Sigma t: " << sigma_t*1E9/c_light << " ps" << endl;
     237      }
     238
     239      fVertexOutputArray->Add(candidate);
     240      k++;
     241    }
     242   }// end of cluster loop
     243
     244  delete ClusterArray;
     245}
     246
     247//------------------------------------------------------------------------------
     248
     249void VertexFinderDA4D::clusterize(TObjArray &clusters)
     250{
     251  tracks_t tks;
     252  fill(tks);
     253  unsigned int nt=tks.getSize();
    158254  if(fVerbose)
    159255  {
    160     cout << " start processing vertices ..." << endl;
    161     cout << " Found " << fInputArray->GetEntriesFast() << " input tracks" << endl;
    162     //loop over input tracks
    163     fItInputArray->Reset();
    164     while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    165     {
    166       pos = candidate->InitialPosition;
    167       mom = candidate->Momentum;
    168 
    169       cout << "pt: " << mom.Pt() << ", eta: " << mom.Eta() << ", phi: " << mom.Phi() << ", z: " << candidate->DZ / 10 << endl;
    170     }
    171   }
    172 
    173   // clusterize tracks in Z
    174   clusterize(*fInputArray, *ClusterArray);
    175 
    176   if(fVerbose)
    177   {
    178     std::cout << " clustering returned  " << ClusterArray->GetEntriesFast() << " clusters  from " << fInputArray->GetEntriesFast() << " selected tracks" << std::endl;
    179   }
    180 
    181   //loop over vertex candidates
    182   ItClusterArray = ClusterArray->MakeIterator();
    183   ItClusterArray->Reset();
    184   while((candidate = static_cast<Candidate *>(ItClusterArray->Next())))
    185   {
    186 
    187     double meantime = 0.;
    188     double expv_x2 = 0.;
    189     double normw = 0.;
    190     double errtime = 0;
    191 
    192     double meanpos = 0.;
    193     double meanerr2 = 0.;
    194     double normpos = 0.;
    195     double errpos = 0.;
    196 
    197     double sumpt2 = 0.;
    198 
    199     int itr = 0;
    200 
    201     if(fVerbose) cout << "this vertex has: " << candidate->GetCandidates()->GetEntriesFast() << " tracks" << endl;
    202 
    203     // loop over tracks belonging to this vertex
    204     TIter it1(candidate->GetCandidates());
    205     it1.Reset();
    206 
    207     while((track = static_cast<Candidate *>(it1.Next())))
    208     {
    209 
    210       itr++;
    211       // TBC: the time is in ns for now TBC
    212       double t = track->InitialPosition.T() / c_light;
    213       double dt = track->ErrorT / c_light;
    214       const double time = t;
    215       const double inverr = 1.0 / dt;
    216       meantime += time * inverr;
    217       expv_x2 += time * time * inverr;
    218       normw += inverr;
    219 
    220       // compute error position TBC
    221       const double pt = track->Momentum.Pt();
    222       const double z = track->DZ / 10.0;
    223       const double err_pt = track->ErrorPT;
    224       const double err_z = track->ErrorDZ;
    225 
    226       const double wi = (pt / (err_pt * err_z)) * (pt / (err_pt * err_z));
    227       meanpos += z * wi;
    228 
    229       meanerr2 += err_z * err_z * wi;
    230       normpos += wi;
    231       sumpt2 += pt * pt;
    232 
    233       // while we are here store cluster index in tracks
    234       track->ClusterIndex = ivtx;
    235     }
    236 
    237     meantime = meantime / normw;
    238     expv_x2 = expv_x2 / normw;
    239     errtime = TMath::Sqrt((expv_x2 - meantime * meantime) / itr);
    240     meanpos = meanpos / normpos;
    241     meanerr2 = meanerr2 / normpos;
    242     errpos = TMath::Sqrt(meanerr2 / itr);
    243 
    244     candidate->Position.SetXYZT(0.0, 0.0, meanpos * 10.0, meantime * c_light);
    245     candidate->PositionError.SetXYZT(0.0, 0.0, errpos * 10.0, errtime * c_light);
    246     candidate->SumPT2 = sumpt2;
    247     candidate->ClusterNDF = itr;
    248     candidate->ClusterIndex = ivtx;
    249 
    250     fVertexOutputArray->Add(candidate);
    251 
    252     ivtx++;
    253 
    254     if(fVerbose)
    255     {
    256       std::cout << "x,y,z";
    257       std::cout << ",t";
    258       std::cout << "=" << candidate->Position.X() / 10.0 << " " << candidate->Position.Y() / 10.0 << " " << candidate->Position.Z() / 10.0;
    259       std::cout << " " << candidate->Position.T() / c_light;
    260 
    261       std::cout << std::endl;
    262       std::cout << "sumpt2 " << candidate->SumPT2 << endl;
    263 
    264       std::cout << "ex,ey,ez";
    265       std::cout << ",et";
    266       std::cout << "=" << candidate->PositionError.X() / 10.0 << " " << candidate->PositionError.Y() / 10.0 << " " << candidate->PositionError.Z() / 10.0;
    267       std::cout << " " << candidate->PositionError.T() / c_light;
    268       std::cout << std::endl;
    269     }
    270   } // end of cluster loop
    271 
    272   if(fVerbose)
    273   {
    274     std::cout << "PrimaryVertexProducerAlgorithm::vertices candidates =" << ClusterArray->GetEntriesFast() << std::endl;
    275   }
    276 
    277   //TBC maybe this can be done later
    278   // sort vertices by pt**2  vertex (aka signal vertex tagging)
    279   /*if(pvs.size()>1){
    280       sort(pvs.begin(), pvs.end(), VertexHigherPtSquared());
    281     }
    282      */
    283 
    284   delete ClusterArray;
     256    cout << "Tracks added: " << nt << endl;
     257  }
     258  if (nt == 0) return;
     259
     260
     261
     262  vertex_t vtx; // the vertex prototypes
     263  vtx.ZSize = fVertexZSize;
     264  vtx.TSize = fVertexTSize;
     265  // initialize:single vertex at infinite temperature
     266  vtx.addItem(0, 0, 1);
     267
     268  // Fit the vertex at T=inf and return the starting temperature
     269  double beta=beta0(tks, vtx);
     270
     271  if( fVerbose > 1 )
     272  {
     273    cout << "Cluster position at T=inf: z = " << vtx.z[0] << " mm , t = " << vtx.t[0] << " ps" << "  pk = " << vtx.pk[0] << endl;
     274    cout << Form("Beta Start = %2.1e", beta) << endl;
     275  }
     276
     277  if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Ast");
     278
     279  if( fVerbose > 2){cout << "Cool down untill reaching the temperature to finish increasing the number of vertexes" << endl;}
     280
     281  double rho0=0.0;  // start with no outlier rejection
     282
     283  unsigned int last_round = 0;
     284  while(last_round < 2)
     285  {
     286
     287    unsigned int niter=0;
     288    double delta2 = 0;
     289    do  {
     290      delta2 = update(beta, tks, vtx, rho0);
     291
     292      if( fVerbose > 10 ) plot_status(beta, vtx, tks, niter, "Bup");
     293      if (fVerbose > 3)
     294      {
     295        cout << "Update " << niter << " : " << delta2 << endl;
     296      }
     297      niter++;
     298    }
     299    while (delta2 > fD2UpdateLim &&  niter < fMaxIterations);
     300
     301
     302    unsigned int n_it = 0;
     303    while(merge(vtx, fD2Merge) && n_it < fMaxIterations)
     304    {
     305      unsigned int niter=0;
     306      double delta2 = 0;
     307      do  {
     308        delta2 = update(beta, tks, vtx, rho0);
     309        niter++;
     310      }
     311      while (delta2 > fD2UpdateLim &&  niter < fMaxIterations);
     312      n_it++;
     313
     314      if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");
     315    }
     316
     317    beta /= fCoolingFactor;
     318
     319    if( beta < fBetaStop )
     320    {
     321      split(beta, vtx, tks);
     322      if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Asp");
     323    }
     324    else
     325    {
     326      beta = fBetaStop;
     327      last_round++;
     328    }
     329
     330    if(fVerbose > 3)
     331    {
     332      cout << endl << endl << " ----- Beta = " << beta << " --------" << endl;
     333      cout << "Nv: " << vtx.getSize() << endl;
     334    }
     335  }
     336
     337  if( fVerbose > 4)
     338  {
     339    for(unsigned int k = 0; k < vtx.getSize(); k++)
     340    {
     341      cout << Form("Vertex %d next beta_c = %.3f", k, vtx.beta_c[k]) << endl;
     342    }
     343  }
     344
     345  if(fVerbose > 2)  {cout << "Adiabatic switch on of outlayr rejection" << endl;}
     346  rho0 = 1./nt;
     347  const double N_cycles = 10;
     348  for(unsigned int f = 1; f <= N_cycles; f++)
     349  {
     350    unsigned int niter=0;
     351    double delta2 = 0;
     352    do  {
     353      delta2 = update(beta, tks, vtx, rho0 * f/N_cycles);
     354      niter++;
     355    }
     356    while (delta2 > 0.3*fD2UpdateLim &&  niter < fMaxIterations);
     357    if( fVerbose > 10 ) plot_status(beta, vtx, tks, f, "Dadout");
     358  }
     359
     360  do {
     361    beta /= fCoolingFactor;
     362    if(beta > fBetaPurge) beta = fBetaPurge;
     363    unsigned int i_pu = 0;
     364    for(int min_trk = 2; min_trk<=fMinNTrack; min_trk++)
     365    {
     366      while( purge(vtx, tks, rho0, beta, fMinTrackProb, min_trk) )
     367      {
     368        unsigned int niter=0;
     369        double delta2 = 0;
     370        do  {
     371          delta2 = update(beta, tks, vtx, rho0);
     372          niter++;
     373        }
     374        while (delta2 > fD2UpdateLim &&  niter < fMaxIterations);
     375        if( fVerbose > 10 ) plot_status(beta, vtx, tks, i_pu, Form("Eprg%d",min_trk));
     376        i_pu++;
     377      }
     378    }
     379
     380    unsigned int n_it = 0;
     381    while(merge(vtx, fD2Merge) && n_it < fMaxIterations)
     382    {
     383      unsigned int niter=0;
     384      double delta2 = 0;
     385      do  {
     386        delta2 = update(beta, tks, vtx, rho0);
     387        niter++;
     388      }
     389      while (delta2 > fD2UpdateLim &&  niter < fMaxIterations);
     390      n_it++;
     391
     392      if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");
     393    }
     394  } while( beta < fBetaPurge );
     395
     396
     397  if(fVerbose > 2){cout << "Cooldown untill the limit before assigning track to vertices" << endl;}
     398  last_round = 0;
     399  while(last_round < 2)
     400  {
     401    unsigned int niter=0;
     402    double delta2 = 0;
     403    do  {
     404      delta2 = update(beta, tks, vtx, rho0);
     405      niter++;
     406      if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Bup");
     407    }
     408    while (delta2 > 0.3*fD2UpdateLim &&  niter < fMaxIterations);
     409
     410    beta /= fCoolingFactor;
     411    if ( beta >= fBetaMax )
     412    {
     413      beta = fBetaMax;
     414      last_round++;
     415    }
     416  }
     417
     418
     419  // Build the cluster candidates
     420  for(unsigned int k = 0; k < vtx.getSize(); k++)
     421  {
     422    DelphesFactory *factory = GetFactory();
     423    Candidate * candidate = factory->NewCandidate();
     424
     425    candidate->ClusterIndex = k;
     426    candidate->Position.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);
     427    candidate->PositionError.SetXYZT(0.0, 0.0, fVertexZSize , fVertexTSize*1E-9*c_light);
     428    candidate->SumPT2 = 0;
     429    candidate->SumPt = 0;
     430    candidate->ClusterNDF = 0;
     431
     432    clusters.Add(candidate);
     433  }
     434
     435
     436  // Assign each track to the most probable vertex
     437  double Z_init = rho0 * exp(-beta * fMuOutlayer * fMuOutlayer); // Add fDtCutOff here toghether  with this
     438  vector<double> pk_exp_mBetaE = Compute_pk_exp_mBetaE(beta, vtx, tks, Z_init);
     439  for(unsigned int i = 0; i< tks.getSize(); i++)
     440  {
     441    if(tks.w[i] <= 0) continue;
     442
     443    double p_max = 0;
     444    unsigned int k_max = 0;
     445
     446    for(unsigned int k = 0; k < vtx.getSize(); k++)
     447    {
     448      unsigned int idx = k*nt + i;
     449      if(pk_exp_mBetaE[idx] == 0 || tks.Z[i] == 0 || vtx.pk[k] == 0)
     450      {
     451        continue;
     452      }
     453
     454      double pv_max = vtx.pk[k] / (vtx.pk[k] + rho0 * exp(-beta * fMuOutlayer* fMuOutlayer));
     455      double p = pk_exp_mBetaE[idx] / tks.Z[i];
     456
     457      p /= pv_max;
     458
     459      if(p > p_max)
     460      {
     461        p_max = p;
     462        k_max = k;
     463      }
     464    }
     465
     466    if(p_max > fMinTrackProb)
     467    {
     468      tks.tt[i]->ClusterIndex = k_max;
     469      tks.tt[i]->InitialPosition.SetT(1E-9*vtx.t[k_max]*c_light);
     470      tks.tt[i]->InitialPosition.SetZ(vtx.z[k_max]);
     471
     472      ((Candidate *) clusters.At(k_max))->AddCandidate(tks.tt[i]);
     473      ((Candidate *) clusters.At(k_max))->SumPT2 += tks.tt[i]->Momentum.Pt()*tks.tt[i]->Momentum.Pt();
     474      ((Candidate *) clusters.At(k_max))->SumPt += tks.tt[i]->Momentum.Pt();
     475      ((Candidate *) clusters.At(k_max))->ClusterNDF += 1;
     476    }
     477    else
     478    {
     479      tks.tt[i]->ClusterIndex = -1;
     480      tks.tt[i]->InitialPosition.SetT(1E3*1000000*c_light);
     481      tks.tt[i]->InitialPosition.SetZ(1E8);
     482    }
     483    fTrackOutputArray->Add(tks.tt[i]);
     484  }
     485
     486  if(fVerbose > 10) plot_status_end(vtx, tks);
     487
    285488}
    286489
    287490//------------------------------------------------------------------------------
    288 
    289 void VertexFinderDA4D::clusterize(const TObjArray &tracks, TObjArray &clusters)
    290 {
    291   if(fVerbose)
    292   {
    293     cout << "###################################################" << endl;
    294     cout << "# VertexFinderDA4D::clusterize   nt=" << tracks.GetEntriesFast() << endl;
    295     cout << "###################################################" << endl;
    296   }
    297 
    298   vector<Candidate *> pv = vertices();
    299 
    300   if(fVerbose)
    301   {
    302     cout << "# VertexFinderDA4D::clusterize   pv.size=" << pv.size() << endl;
    303   }
    304   if(pv.size() == 0)
    305   {
    306     return;
    307   }
    308 
    309   // convert into vector of candidates
    310   //TObjArray *ClusterArray = pv.begin()->GetCandidates();
    311   //Candidate *aCluster = static_cast<Candidate*>(&(pv.at(0)));
    312   Candidate *aCluster = pv.at(0);
    313 
    314   // fill into clusters and merge
    315 
    316   if(fVerbose)
    317   {
    318     std::cout << '\t' << 0;
    319     std::cout << ' ' << (*pv.begin())->Position.Z() / 10.0 << ' ' << (*pv.begin())->Position.T() / c_light << std::endl;
    320   }
    321 
    322   for(vector<Candidate *>::iterator k = pv.begin() + 1; k != pv.end(); k++)
    323   {
    324     if(fVerbose)
    325     {
    326       std::cout << '\t' << std::distance(pv.begin(), k);
    327       std::cout << ' ' << (*k)->Position.Z() << ' ' << (*k)->Position.T() << std::endl;
    328     }
    329 
    330     // TBC - check units here
    331     if(std::abs((*k)->Position.Z() - (*(k - 1))->Position.Z()) / 10.0 > (2 * fVertexSpaceSize) || std::abs((*k)->Position.T() - (*(k - 1))->Position.Z()) / c_light > 2 * 0.010)
    332     {
    333       // close a cluster
    334       clusters.Add(aCluster);
    335       //aCluster.clear();
    336     }
    337     //for(unsigned int i=0; i<k->GetCandidates().GetEntriesFast(); i++){
    338     aCluster = *k;
    339     //}
    340   }
    341   clusters.Add(aCluster);
    342 
    343   if(fVerbose)
    344   {
    345     std::cout << "# VertexFinderDA4D::clusterize clusters.size=" << clusters.GetEntriesFast() << std::endl;
    346   }
     491// Definition of the distance metrci between track and vertex
     492double VertexFinderDA4D::Energy(double t_z, double v_z, double dz2_o, double t_t, double v_t, double dt2_o)
     493{
     494  return (t_z - v_z)*(t_z - v_z)* dz2_o + (t_t - v_t)*(t_t - v_t)*dt2_o;
    347495}
    348496
    349497//------------------------------------------------------------------------------
    350 
    351 vector<Candidate *> VertexFinderDA4D::vertices()
    352 {
     498// Fill tks with the input candidates array
     499void VertexFinderDA4D::fill(tracks_t &tks)
     500{
     501  tks.sum_w_o_dt2 = 0;
     502  tks.sum_w_o_dz2 = 0;
     503  tks.sum_w = 0;
     504
    353505  Candidate *candidate;
    354   UInt_t clusterIndex = 0;
    355   vector<Candidate *> clusters;
    356 
    357   vector<track_t> tks;
    358   track_t tr;
    359   Double_t z, dz, t, l, dt, d0, d0error;
    360 
    361   // loop over input tracks
     506
    362507  fItInputArray->Reset();
    363   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    364   {
    365     //TBC everything in cm
    366     z = candidate->DZ / 10;
    367     tr.z = z;
    368     dz = candidate->ErrorDZ / 10;
    369     tr.dz2 = dz * dz // track error
    370       //TBC: beamspot size induced error, take 0 for now.
    371       // + (std::pow(beamspot.BeamWidthX()*cos(phi),2.)+std::pow(beamspot.BeamWidthY()*sin(phi),2.))/std::pow(tantheta,2.) // beam-width induced
    372       + fVertexSpaceSize * fVertexSpaceSize; // intrinsic vertex size, safer for outliers and short lived decays
    373 
    374     // TBC: the time is in ns for now TBC
    375     //t = candidate->Position.T()/c_light;
    376     t = candidate->InitialPosition.T() / c_light;
    377     l = candidate->L / c_light;
     508  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
     509  {
     510    unsigned int discard = 0;
     511
    378512    double pt = candidate->Momentum.Pt();
    379     double eta = candidate->Momentum.Eta();
    380     double phi = candidate->Momentum.Phi();
    381 
    382     tr.pt = pt;
    383     tr.eta = eta;
    384     tr.phi = phi;
    385     tr.t = t; //
    386     tr.dtz = 0.;
    387     dt = candidate->ErrorT / c_light;
    388     tr.dt2 = dt * dt + fVertexTimeSize * fVertexTimeSize; // the ~injected~ timing error plus a small minimum vertex size in time
    389     if(fD0CutOff > 0)
    390     {
    391 
    392       d0 = TMath::Abs(candidate->D0) / 10.0;
    393       d0error = candidate->ErrorD0 / 10.0;
    394 
    395       tr.pi = 1. / (1. + exp((d0 * d0) / (d0error * d0error) - fD0CutOff * fD0CutOff)); // reduce weight for high ip tracks
     513    if(pt<fPtMin || pt>fPtMax) discard = 1;
     514
     515    // ------------- Compute cloasest approach Z ----------------
     516    double z = candidate->DZ; // [mm]
     517
     518    candidate->Zd = candidate->DZ; //Set the cloasest approach z
     519    if(fabs(z) > 3*fDzCutOff) discard = 1;
     520
     521    // ------------- Compute cloasest approach T ----------------
     522    //Asumme pion mass which is the most common particle
     523    double M = 0.139570;
     524    candidate->Mass = M;
     525    double p = pt * sqrt(1 + candidate->CtgTheta*candidate->CtgTheta);
     526    double e = sqrt(p*p + M*M);
     527
     528    double t = candidate->Position.T()*1.E9/c_light; // from [mm] to [ps]
     529    if(t <= -9999) discard = 1;                    // Means that the time information has not been added
     530
     531    // DEBUG Here backpropagete for the whole length and not noly for z. Could improve resolution
     532    // double bz = pt * candidate->CtgTheta/e;
     533    // t += (z - candidate->Position.Z())*1E9/(c_light*bz);
     534
     535    // Use full path Length
     536    t -= candidate->L*1E9/(c_light*p/e);
     537
     538    candidate->Td = t*1E-9*c_light;
     539    if(fabs(t) > 3*fDtCutOff) discard = 1;
     540
     541    // auto genp = (Candidate*) candidate->GetCandidates()->At(0);
     542    // cout << "Eta: " << candidate->Position.Eta() << endl;
     543    // cout << genp->Momentum.Pt() << " -- " << candidate->Momentum.Pt() << endl;
     544    // cout << genp->Momentum.Pz() << " -- " << candidate->Momentum.Pz() << endl;
     545    // cout << genp->Momentum.P() << " -- " << p << endl;
     546    // cout << genp->Momentum.E() << " -- " << e << endl;
     547    // cout << Form("bz_true: %.4f -- bz_gen: %.4f", genp->Momentum.Pz()/genp->Momentum.E(), bz) << endl;
     548
     549    double dz2_o = candidate->ErrorDZ*candidate->ErrorDZ;
     550    dz2_o += fVertexZSize*fVertexZSize;
     551    // when needed add beam spot width (x-y)?? mha?
     552    dz2_o = 1/dz2_o; //Multipling is faster than dividing all the times
     553
     554    double dt2_o = candidate->ErrorT*1.E9/c_light; // [ps]
     555    dt2_o *= dt2_o;
     556    dt2_o += fVertexTSize*fVertexTSize; // [ps^2]
     557    // Ideally we should also add the induced uncertantiy from dz, z_out, pt, ctgthetaand all the other thing used above (total around 5ps). For the moment we compensae using a high value for vertex time.
     558    dt2_o = 1/dt2_o;
     559
     560    double w;
     561    if(fD0CutOff > 0 && candidate->ErrorD0 > 0)
     562    {
     563      double d0_sig = fabs(candidate->D0/candidate->ErrorD0);
     564      w = exp(d0_sig*d0_sig - fD0CutOff*fD0CutOff);
     565      w = 1./(1. + w);
     566      if (w < 1E-4) discard = 1;
    396567    }
    397568    else
    398569    {
    399       tr.pi = 1.;
    400     }
    401     tr.tt = &(*candidate);
    402     tr.Z = 1.;
    403 
    404     // TBC now putting track selection here (> fPTMin)
    405     if(tr.pi > 1e-3 && tr.pt > fMinPT)
    406     {
    407       tks.push_back(tr);
    408     }
    409   }
    410 
    411   //print out input tracks
    412 
    413   if(fVerbose)
    414   {
    415     std::cout << " start processing vertices ..." << std::endl;
    416     std::cout << " Found " << tks.size() << " input tracks" << std::endl;
    417     //loop over input tracks
    418 
    419     for(std::vector<track_t>::const_iterator it = tks.begin(); it != tks.end(); it++)
    420     {
    421       double z = it->z;
    422       double pt = it->pt;
    423       double eta = it->eta;
    424       double phi = it->phi;
    425       double t = it->t;
    426 
    427       std::cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", z: " << z << ", t: " << t << std::endl;
    428     }
    429   }
    430 
    431   unsigned int nt = tks.size();
    432   double rho0 = 0.0; // start with no outlier rejection
    433 
    434   if(tks.empty()) return clusters;
    435 
    436   vector<vertex_t> y; // the vertex prototypes
    437 
    438   // initialize:single vertex at infinite temperature
    439   vertex_t vstart;
    440   vstart.z = 0.;
    441   vstart.t = 0.;
    442   vstart.pk = 1.;
    443   y.push_back(vstart);
    444   int niter = 0; // number of iterations
    445 
    446   // estimate first critical temperature
    447   double beta = beta0(fBetaMax, tks, y, fCoolingFactor);
    448   niter = 0;
    449   while((update1(beta, tks, y) > 1.e-6) && (niter++ < fMaxIterations))
    450   {
    451   }
    452 
    453   // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
    454   while(beta < fBetaMax)
    455   {
    456 
    457     if(fUseTc)
    458     {
    459       update1(beta, tks, y);
    460       while(merge(y, beta))
     570      w = 1;
     571    }
     572    candidate->VertexingWeight = w;
     573
     574
     575    if(discard)
     576    {
     577      candidate->ClusterIndex = -1;
     578      candidate->InitialPosition.SetT(1E3*1000000*c_light);
     579      candidate->InitialPosition.SetZ(1E8);
     580      fTrackOutputArray->Add(candidate);
     581    }
     582    else
     583    {
     584      tks.sum_w_o_dt2 += w * dt2_o;
     585      tks.sum_w_o_dz2 += w * dz2_o;
     586      tks.sum_w += w;
     587      tks.addItem(z, t, dz2_o, dt2_o, &(*candidate), w, candidate->PID); //PROVA: rimuovi &(*---)
     588    }
     589
     590  }
     591
     592  if(fVerbose > 1)
     593  {
     594    cout << "----->Filled tracks" << endl;
     595    cout << "M        z           dz        t            dt        w" << endl;
     596    for(unsigned int i = 0; i < tks.getSize(); i++)
     597    {
     598      cout << Form("%d\t%1.1e\t%1.1e\t%1.1e\t%1.1e\t%1.1e", tks.PID[i], tks.z[i], 1/sqrt(tks.dz2_o[i]), tks.t[i], 1/sqrt(tks.dt2_o[i]), tks.w[i]) << endl;
     599    }
     600  }
     601
     602  return;
     603}
     604
     605//------------------------------------------------------------------------------
     606// Compute higher phase transition temperature
     607double VertexFinderDA4D::beta0(tracks_t & tks, vertex_t &vtx)
     608{
     609  if(vtx.getSize() != 1)
     610  {
     611    throw std::invalid_argument( "Unexpected number of vertices" );
     612  }
     613
     614  unsigned int nt = tks.getSize();
     615
     616  //Set vertex position at T=inf as the weighted average of the tracks
     617  double sum_wz = 0, sum_wt = 0;
     618  for(unsigned int i = 0; i < nt; i++)
     619  {
     620    sum_wz += tks.w[i] * tks.z[i] * tks.dz2_o[i];
     621    sum_wt += tks.w[i] * tks.t[i] * tks.dt2_o[i];
     622  }
     623  vtx.t[0] = sum_wt / tks.sum_w_o_dt2;
     624  vtx.z[0] = sum_wz / tks.sum_w_o_dz2;
     625
     626  // Compute the posterior distribution covariance matrix elements
     627  double s_zz = 0, s_tt = 0, s_tz = 0;
     628  for(unsigned int i = 0; i < nt; i++)
     629  {
     630    double dz = (tks.z[i] - vtx.z[0]) * tks.dz_o[i];
     631    double dt = (tks.t[i] - vtx.t[0]) * tks.dt_o[i];
     632
     633    s_zz += tks.w[i] * dz * dz;
     634    s_tt += tks.w[i] * dt * dt;
     635    s_tz += tks.w[i] * dt * dz;
     636  }
     637  s_tt /= tks.sum_w;
     638  s_zz /= tks.sum_w;
     639  s_tz /= tks.sum_w;
     640
     641  // Copute the max eighenvalue
     642  double beta_c = (s_tt - s_zz)*(s_tt - s_zz) + 4*s_tz*s_tz;
     643  beta_c = 1. / (s_tt + s_zz + sqrt(beta_c));
     644
     645  double out;
     646  if (beta_c < fBetaMax)
     647  {
     648    // Cool down up to a step before the phase transition
     649    out = beta_c * sqrt(fCoolingFactor);
     650  }
     651  else
     652  {
     653    out = fBetaMax * fCoolingFactor;
     654  }
     655
     656  return out;
     657}
     658
     659//------------------------------------------------------------------------------
     660// Compute the new vertexes position and mass (probability) -- mass constrained annealing without noise
     661// Compute and store the posterior covariance matrix elements
     662// Returns the squared sum of changes of vertexex position normalized by the vertex size declared in the init
     663double VertexFinderDA4D::update(double beta, tracks_t &tks, vertex_t &vtx, double rho0)
     664{
     665  unsigned int nt = tks.getSize();
     666  unsigned int nv = vtx.getSize();
     667
     668  //initialize sums
     669  double Z_init = rho0 * exp(-beta * fMuOutlayer * fMuOutlayer);
     670
     671  // Compute all the energies (aka distances) and normalization partition function
     672  vector<double> pk_exp_mBetaE = Compute_pk_exp_mBetaE(beta, vtx, tks, Z_init);
     673
     674  double sum_pk = 0;
     675  double delta2_max = 0;
     676  for (unsigned int k = 0; k < nv; k++)
     677  {
     678    // Compute the new vertex positions and masses
     679    double pk_new = 0;
     680    double sw_z = 0, sw_t = 0;
     681    // Compute the posterior covariance matrix Elements
     682    double szz = 0, stt = 0, stz = 0;
     683    double sum_wt = 0, sum_wz = 0;
     684    double sum_ptt = 0, sum_pzz = 0, sum_ptz = 0;
     685
     686
     687    for (unsigned int i = 0; i < nt; i++)
     688    {
     689      unsigned int idx = k*nt + i;
     690
     691      if(pk_exp_mBetaE[idx] == 0 || tks.Z[i] == 0)
    461692      {
    462         update1(beta, tks, y);
    463       }
    464       split(beta, tks, y);
    465       beta = beta / fCoolingFactor;
     693        continue;
     694      }
     695
     696      double p_ygx = pk_exp_mBetaE[idx] / tks.Z[i];      //p(y|x), Gibbs distribution
     697      if(std::isnan(p_ygx) || std::isinf(p_ygx) || p_ygx > 1)
     698      {
     699        cout << Form("%1.6e    %1.6e", pk_exp_mBetaE[idx], tks.Z[i]);
     700        throw std::invalid_argument(Form("p_ygx is %.8f", p_ygx));
     701      }
     702      pk_new += tks.w[i] * p_ygx;
     703
     704      double wt = tks.w[i] * p_ygx * tks.dt2_o[i];
     705      sw_t += wt * tks.t[i];
     706      sum_wt += wt;
     707
     708      double wz = tks.w[i] * p_ygx * tks.dz2_o[i];
     709      sw_z += wz * tks.z[i];
     710      sum_wz += wz;
     711
     712      // Add the track contribution to the covariance matrix
     713      double p_xgy = p_ygx * tks.w[i] / vtx.pk[k];
     714      double dt = (tks.t[i] - vtx.t[k]) * tks.dt_o[i];
     715      double dz = (tks.z[i] - vtx.z[k]) * tks.dz_o[i];
     716
     717      double wtt = p_xgy * tks.dt2_o[i];
     718      double wzz = p_xgy * tks.dz2_o[i];
     719      double wtz = p_xgy * tks.dt_o[i] * tks.dz_o[i];
     720
     721      stt += wtt * dt * dt;
     722      szz += wzz * dz * dz;
     723      stz += wtz * dt * dz;
     724
     725      sum_ptt += wtt;
     726      sum_pzz += wzz;
     727      sum_ptz += wtz;
     728    }
     729    if(pk_new == 0)
     730    {
     731      vtx.removeItem(k);
     732      k--;
     733      // throw std::invalid_argument(Form("pk_new is %.8f", pk_new));
    466734    }
    467735    else
    468736    {
    469       beta = beta / fCoolingFactor;
    470       splitAll(y);
    471     }
    472 
    473     // make sure we are not too far from equilibrium before cooling further
    474     niter = 0;
    475     while((update1(beta, tks, y) > 1.e-6) && (niter++ < fMaxIterations))
    476     {
    477     }
    478   }
    479 
    480   if(fUseTc)
    481   {
    482     // last round of splitting, make sure no critical clusters are left
    483     update1(beta, tks, y);
    484     while(merge(y, beta))
    485     {
    486       update1(beta, tks, y);
    487     }
    488     unsigned int ntry = 0;
    489     while(split(beta, tks, y) && (ntry++ < 10))
    490     {
    491       niter = 0;
    492       while((update1(beta, tks, y) > 1.e-6) && (niter++ < fMaxIterations))
     737      pk_new /= tks.sum_w;
     738      sum_pk += pk_new;
     739
     740      stt /= sum_ptt;
     741      szz /= sum_pzz;
     742      stz /= sum_ptz;
     743
     744      double new_t = sw_t/sum_wt;
     745      double new_z = sw_z/sum_wz;
     746      if(std::isnan(new_z) || std::isnan(new_t))
    493747      {
    494       }
    495       merge(y, beta);
    496       update1(beta, tks, y);
    497     }
     748        cout << endl << endl;
     749        cout << Form("t: %.3e   /   %.3e", sw_t, sum_wt) << endl;
     750        cout << Form("z: %.3e   /   %.3e", sw_z, sum_wz) << endl;
     751        cout << "pk " << k << "  " << vtx.pk[k] << endl;
     752        throw std::invalid_argument("new_z is nan");
     753      }
     754
     755      double z_displ = (new_z - vtx.z[k])/fVertexZSize;
     756      double t_displ = (new_t - vtx.t[k])/fVertexTSize;
     757      double delta2 = z_displ*z_displ + t_displ*t_displ;
     758
     759      if (delta2 > delta2_max) delta2_max =  delta2;
     760
     761      vtx.z[k] = new_z;
     762      vtx.t[k] = new_t;
     763      vtx.pk[k] = pk_new;
     764      vtx.szz[k] = szz;
     765      vtx.stt[k] = stt;
     766      vtx.stz[k] = stz;
     767    }
     768  }
     769
     770  if(fabs((sum_pk - 1.) > 1E-4))
     771  {
     772    cout << "sum_pk " << sum_pk << endl;
     773    for (unsigned int k = 0; k < nv; k++)
     774    {
     775      cout << Form("%d: %1.4e", k, vtx.pk[k]) << endl;
     776    }
     777    throw std::invalid_argument("Sum of masses not unitary");
     778  }
     779  // if(fVerbose > 3)
     780  // {
     781  //   cout << "===Update over" << endl;
     782  //   for (unsigned int k = 0; k < nv; k++)
     783  //   {
     784  //     cout << k << endl;
     785  //     cout << "z: " << vtx.z[k] << " , t: " << vtx.t[k] << " , p: " << vtx.pk[k] << endl;
     786  //     cout << " | " << vtx.szz[k] << "   " << vtx.stz[k] << "|" << endl;
     787  //     cout << " | " << vtx.stz[k] << "   " << vtx.stt[k] << "|" << endl << endl;
     788  //   }
     789  //   cout << "=======" << endl;
     790  // }
     791
     792  return delta2_max;
     793}
     794
     795//------------------------------------------------------------------------------
     796// Split critical vertices (beta_c < beta)
     797// Returns true if at least one cluster was split
     798bool VertexFinderDA4D::split(double &beta, vertex_t &vtx, tracks_t & tks)
     799{
     800  bool split = false;
     801
     802  auto pair_bc_k = vtx.ComputeAllBeta_c(fVerbose);
     803
     804  // If minimum beta_c is higher than beta, no split is necessaire
     805  if( pair_bc_k.first > beta )
     806  {
     807    split = false;
    498808  }
    499809  else
    500810  {
    501     // merge collapsed clusters
    502     while(merge(y, beta))
    503     {
    504       update1(beta, tks, y);
    505     }
    506     if(fVerbose)
    507     {
    508       cout << "dump after 1st merging " << endl;
    509       dump(beta, y, tks);
    510     }
    511   }
    512 
    513   // switch on outlier rejection
    514   rho0 = 1. / nt;
    515   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    516   {
    517     k->pk = 1.;
    518   } // democratic
    519   niter = 0;
    520   while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-8) && (niter++ < fMaxIterations))
    521   {
    522   }
    523   if(fVerbose)
    524   {
    525     cout << "rho0=" << rho0 << " niter=" << niter << endl;
    526     dump(beta, y, tks);
    527   }
    528 
    529   // merge again  (some cluster split by outliers collapse here)
    530   while(merge(y))
    531   {
    532   }
    533   if(fVerbose)
    534   {
    535     cout << "dump after 2nd merging " << endl;
    536     dump(beta, y, tks);
    537   }
    538 
    539   // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices
    540   while(beta <= fBetaStop)
    541   {
    542     while(purge(y, tks, rho0, beta, fDzCutOff))
    543     {
    544       niter = 0;
    545       while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-6) && (niter++ < fMaxIterations))
     811    const unsigned int nv = vtx.getSize();
     812    for(unsigned int k = 0; k < nv; k++)
     813    {
     814      if( fVerbose > 3 )
    546815      {
    547       }
    548     }
    549     beta /= fCoolingFactor;
    550     niter = 0;
    551     while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-6) && (niter++ < fMaxIterations))
    552     {
    553     }
    554   }
    555 
    556   //   // new, one last round of cleaning at T=Tstop
    557   //   while(purge(y,tks,rho0, beta)){
    558   //     niter=0; while((update2(beta, tks,y,rho0, fDzCutOff) > 1.e-6)  && (niter++ < fMaxIterations)){  }
    559   //   }
    560 
    561   if(fVerbose)
    562   {
    563     cout << "Final result, rho0=" << rho0 << endl;
    564     dump(beta, y, tks);
    565   }
    566 
    567   // select significant tracks and use a TransientVertex as a container
    568   //GlobalError dummyError;
    569 
    570   // ensure correct normalization of probabilities, should make double assginment reasonably impossible
    571   for(unsigned int i = 0; i < nt; i++)
    572   {
    573     tks[i].Z = rho0 * exp(-beta * (fDzCutOff * fDzCutOff));
    574     for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    575     {
    576       tks[i].Z += k->pk * exp(-beta * Eik(tks[i], *k));
    577     }
    578   }
    579 
    580   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    581   {
    582 
    583     DelphesFactory *factory = GetFactory();
    584     candidate = factory->NewCandidate();
    585 
    586     //cout<<"new vertex"<<endl;
    587     //GlobalPoint pos(0, 0, k->z);
    588     double time = k->t;
    589     double z = k->z;
    590     //vector< reco::TransientTrack > vertexTracks;
    591     //double max_track_time_err2 = 0;
    592     double mean = 0.;
    593     double expv_x2 = 0.;
    594     double normw = 0.;
    595     for(unsigned int i = 0; i < nt; i++)
    596     {
    597       const double invdt = 1.0 / std::sqrt(tks[i].dt2);
    598       if(tks[i].Z > 0)
     816        cout << "vtx " << k << "  beta_c = " << vtx.beta_c[k] << endl;
     817      }
     818      if(vtx.beta_c[k] <= beta)
    599819      {
    600         double p = k->pk * exp(-beta * Eik(tks[i], *k)) / tks[i].Z;
    601         if((tks[i].pi > 0) && (p > 0.5))
     820        double z_old = vtx.z[k];
     821        double t_old = vtx.t[k];
     822        double pk_old = vtx.pk[k];
     823
     824        // Compute splitting direction: given by the max eighenvalue eighenvector
     825        double zn = (vtx.szz[k] - vtx.stt[k])*(vtx.szz[k] - vtx.stt[k]) + 4*vtx.stz[k]*vtx.stz[k];
     826        zn = vtx.szz[k] - vtx.stt[k] + sqrt(zn);
     827        double tn = 2*vtx.stz[k];
     828        double norm = hypot(zn, tn);
     829        tn /= norm;
     830        zn /= norm;
     831
     832        // Estimate subcluster positions and weight
     833        double p1=0, z1=0, t1=0, wz1=0, wt1=0;
     834        double p2=0, z2=0, t2=0, wz2=0, wt2=0;
     835        const unsigned int nt = tks.getSize();
     836        for(unsigned int i=0; i<nt; ++i)
    602837        {
    603           //std::cout << "pushing back " << i << ' ' << tks[i].tt << std::endl;
    604           //vertexTracks.push_back(*(tks[i].tt)); tks[i].Z=0;
    605 
    606           candidate->AddCandidate(tks[i].tt);
    607           tks[i].Z = 0;
    608 
    609           mean += tks[i].t * invdt * p;
    610           expv_x2 += tks[i].t * tks[i].t * invdt * p;
    611           normw += invdt * p;
    612         } // setting Z=0 excludes double assignment
    613       }
    614     }
    615 
    616     mean = mean / normw;
    617     expv_x2 = expv_x2 / normw;
    618     const double time_var = expv_x2 - mean * mean;
    619     const double crappy_error_guess = std::sqrt(time_var);
    620     /*GlobalError dummyErrorWithTime(0,
    621                                    0,0,
    622                                    0,0,0,
    623                                    0,0,0,crappy_error_guess);*/
    624     //TransientVertex v(pos, time, dummyErrorWithTime, vertexTracks, 5);
    625 
    626     candidate->ClusterIndex = clusterIndex++;
    627     ;
    628     candidate->Position.SetXYZT(0.0, 0.0, z * 10.0, time * c_light);
    629 
    630     // TBC - fill error later ...
    631     candidate->PositionError.SetXYZT(0.0, 0.0, 0.0, crappy_error_guess * c_light);
    632 
    633     clusterIndex++;
    634     clusters.push_back(candidate);
    635   }
    636 
    637   return clusters;
    638 }
    639 
    640 //------------------------------------------------------------------------------
    641 
    642 static double Eik(const track_t &t, const vertex_t &k)
    643 {
    644   return std::pow(t.z - k.z, 2.) / t.dz2 + std::pow(t.t - k.t, 2.) / t.dt2;
    645 }
    646 
    647 //------------------------------------------------------------------------------
    648 
    649 static void dump(const double beta, const vector<vertex_t> &y, const vector<track_t> &tks0)
    650 {
    651   // copy and sort for nicer printout
    652   vector<track_t> tks;
    653   for(vector<track_t>::const_iterator t = tks0.begin(); t != tks0.end(); t++)
    654   {
    655     tks.push_back(*t);
    656   }
    657   std::stable_sort(tks.begin(), tks.end(), recTrackLessZ1);
    658 
    659   cout << "-----DAClusterizerInZT::dump ----" << endl;
    660   cout << " beta=" << beta << endl;
    661   cout << "                                                               z= ";
    662   cout.precision(4);
    663   for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
    664   {
    665     //cout  <<  setw(8) << fixed << k->z;
    666   }
    667   cout << endl
    668        << "                                                               t= ";
    669   for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
    670   {
    671     //cout  <<  setw(8) << fixed << k->t;
    672   }
    673   //cout << endl << "T=" << setw(15) << 1./beta <<"                                             Tc= ";
    674   for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
    675   {
    676     //cout  <<  setw(8) << fixed << k->Tc ;
    677   }
    678 
    679   cout << endl
    680        << "                                                              pk=";
    681   double sumpk = 0;
    682   for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
    683   {
    684     //cout <<  setw(8) <<  setprecision(3) <<  fixed << k->pk;
    685     sumpk += k->pk;
    686   }
    687   cout << endl;
    688 
    689   double E = 0, F = 0;
    690   cout << endl;
    691   cout << "----       z +/- dz        t +/- dt        ip +/-dip       pt    phi  eta    weights  ----" << endl;
    692   cout.precision(4);
    693   for(unsigned int i = 0; i < tks.size(); i++)
    694   {
    695     if(tks[i].Z > 0)
    696     {
    697       F -= log(tks[i].Z) / beta;
    698     }
    699     double tz = tks[i].z;
    700     double tt = tks[i].t;
    701     //cout <<  setw (3)<< i << ")" <<  setw (8) << fixed << setprecision(4)<<  tz << " +/-" <<  setw (6)<< sqrt(tks[i].dz2)
    702     //     << setw(8) << fixed << setprecision(4) << tt << " +/-" << setw(6) << std::sqrt(tks[i].dt2)  ;
    703 
    704     double sump = 0.;
    705     for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
    706     {
    707       if((tks[i].pi > 0) && (tks[i].Z > 0))
    708       {
    709         //double p=pik(beta,tks[i],*k);
    710         double p = k->pk * std::exp(-beta * Eik(tks[i], *k)) / tks[i].Z;
    711         if(p > 0.0001)
     838          if (tks.Z[i] > 0)
     839          {
     840            double lr = (tks.t[i] - vtx.t[k]) * tn + (tks.z[i]-vtx.z[k]) * zn;
     841            // winner-takes-all, usually overestimates splitting
     842            double tl = lr < 0 ? 1.: 0.;
     843            double tr = 1. - tl;
     844
     845            // soften it, especially at low T
     846            // double arg = lr * sqrt(beta * ( zn*zn*tks.dz2_o[i] + tn*tn*tks.dt2_o[i] ) );
     847            // if(abs(arg) < 20)
     848            // {
     849            //   double t = exp(-arg);
     850            //   tl = t/(t+1.);
     851            //   tr = 1/(t+1.);
     852            // }
     853
     854            double p = vtx.pk[k] * tks.w[i];
     855            p *= exp(-beta * Energy(tks.z[i], vtx.z[k], tks.dz2_o[i], tks.t[i], vtx.t[k], tks.dt2_o[i])) / tks.Z[i];
     856            double wt = p*tks.dt2_o[i];
     857            double wz = p*tks.dz2_o[i];
     858            p1 += p*tl;  z1 += wz*tl*tks.z[i]; t1 += wt*tl*tks.t[i]; wz1 += wz*tl; wt1 += wt*tl;
     859            p2 += p*tr;  z2 += wz*tr*tks.z[i]; t2 += wt*tr*tks.t[i]; wz2 += wz*tr; wt2 += wt*tr;
     860          }
     861        }
     862
     863        if(wz1 > 0  && wt1 > 0 && wz2 > 0 && wt2 > 0)
    712864        {
    713           //cout <<  setw (8) <<  setprecision(3) << p;
     865          t1 /= wt1;
     866          z1 /= wz1;
     867          t2 /= wt2;
     868          z2 /= wz2;
     869
     870          if( fVerbose > 3 )
     871          {
     872            double aux = (z1-z2)*(z1-z2)/(fVertexZSize*fVertexZSize) + (t1-t2)*(t1-t2)/(fVertexTSize*fVertexTSize);
     873            cout << "weighted split:  delta = " << sqrt(aux) << endl;
     874          }
    714875        }
    715876        else
    716877        {
    717           cout << "    .   ";
     878          continue;
     879          // plot_split_crush(zn, tn, vtx, tks, k);
     880          // throw std::invalid_argument( "0 division" );
    718881        }
    719         E += p * Eik(tks[i], *k);
    720         sump += p;
    721       }
    722       else
     882
     883        while(vtx.NearestCluster(t1, z1) != k || vtx.NearestCluster(t2, z2) != k)
     884        {
     885          t1 = 0.5 * (t1 + t_old);
     886          z1 = 0.5 * (z1 + z_old);
     887          t2 = 0.5 * (t2 + t_old);
     888          z2 = 0.5 * (z2 + z_old);
     889        }
     890
     891        // Compute final distance and split if the distance is enough
     892        double delta2 = (z1-z2)*(z1-z2)/(fVertexZSize*fVertexZSize) + (t1-t2)*(t1-t2)/(fVertexTSize*fVertexTSize);
     893        if(delta2 > fD2Merge)
     894        {
     895          split = true;
     896          vtx.t[k] = t1;
     897          vtx.z[k] = z1;
     898          vtx.pk[k] = p1 * pk_old/(p1+p2);
     899
     900          double new_t = t2;
     901          double new_z = z2;
     902          double new_pk = p2 * pk_old/(p1+p2);
     903
     904          vtx.addItem(new_z, new_t, new_pk);
     905
     906          if( fVerbose > 3 )
     907          {
     908            cout << "===Split happened on vtx " << k << endl;
     909            cout << "OLD     z: " << z_old << " , t: " << t_old << " , pk: " << pk_old << endl;
     910            cout << "NEW+    z: " << vtx.z[k] << " , t: " << vtx.t[k] << " , pk: " << vtx.pk[k] << endl;
     911            cout << "NEW-    z: " << new_z << " , t: " << new_t << " , pk: " << new_pk <<  endl;
     912          }
     913        }
     914      }
     915    }
     916  }
     917  return split;
     918}
     919
     920
     921//------------------------------------------------------------------------------
     922// Merge vertexes closer than declared dimensions
     923bool VertexFinderDA4D::merge(vertex_t & vtx, double d2_merge = 2)
     924{
     925  bool merged = false;
     926
     927  if(vtx.getSize() < 2) return merged;
     928
     929  bool last_merge = false;
     930  do {
     931    double min_d2 = d2_merge;
     932    unsigned int k1_min, k2_min;
     933    for(unsigned int k1 = 0; k1 < vtx.getSize(); k1++)
     934    {
     935      for(unsigned int k2 = k1+1; k2 < vtx.getSize();k2++)
    723936      {
    724         cout << "        ";
    725       }
    726     }
    727     cout << endl;
    728   }
    729   cout << endl
    730        << "T=" << 1 / beta << " E=" << E << " n=" << y.size() << "  F= " << F << endl
    731        << "----------" << endl;
     937        double d2_tmp = vtx.DistanceSquare(k1, k2);
     938        if(d2_tmp < min_d2)
     939        {
     940          min_d2 = d2_tmp;
     941          k1_min = k1;
     942          k2_min = k2;
     943        }
     944      }
     945    }
     946
     947    if(min_d2 < d2_merge)
     948    {
     949      vtx.mergeItems(k1_min, k2_min);
     950      last_merge = true;
     951      merged = true;
     952    }
     953    else last_merge = false;
     954  } while(last_merge);
     955
     956  return merged;
     957}
     958
     959// -----------------------------------------------------------------------------
     960// Compute all the energies and set the partition function normalization for each track
     961vector<double> VertexFinderDA4D::Compute_pk_exp_mBetaE(double beta, vertex_t &vtx, tracks_t &tks, double Z_init)
     962{
     963  unsigned int nt = tks.getSize();
     964  unsigned int nv = vtx.getSize();
     965
     966  vector<double> pk_exp_mBetaE(nt * nv);
     967  for (unsigned int k = 0; k < nv; k++)
     968  {
     969    for (unsigned int i = 0; i < nt; i++)
     970    {
     971      if(k == 0) tks.Z[i] = Z_init;
     972
     973      double aux = Energy(tks.z[i], vtx.z[k], tks.dz2_o[i], tks.t[i], vtx.t[k], tks.dt2_o[i]);
     974      aux = vtx.pk[k] * exp(-beta * aux);
     975      // if(aux < 1E-10) continue;
     976      tks.Z[i] += aux;
     977
     978      unsigned int idx = k*nt + i;
     979      pk_exp_mBetaE[idx] = aux;
     980    }
     981  }
     982  return pk_exp_mBetaE;
    732983}
    733984
    734985//------------------------------------------------------------------------------
    735 
    736 static double update1(double beta, vector<track_t> &tks, vector<vertex_t> &y)
    737 {
    738   //update weights and vertex positions
    739   // mass constrained annealing without noise
    740   // returns the squared sum of changes of vertex positions
    741 
    742   unsigned int nt = tks.size();
    743 
    744   //initialize sums
    745   double sumpi = 0;
    746   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); ++k)
    747   {
    748     k->sw = 0.;
    749     k->swz = 0.;
    750     k->swt = 0.;
    751     k->se = 0.;
    752     k->swE = 0.;
    753     k->Tc = 0.;
    754   }
    755 
    756   // loop over tracks
    757   for(unsigned int i = 0; i < nt; i++)
    758   {
    759 
    760     // update pik and Zi
    761     double Zi = 0.;
    762     for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); ++k)
    763     {
    764       k->ei = std::exp(-beta * Eik(tks[i], *k)); // cache exponential for one track at a time
    765       Zi += k->pk * k->ei;
    766     }
    767     tks[i].Z = Zi;
    768 
    769     // normalization for pk
    770     if(tks[i].Z > 0)
    771     {
    772       sumpi += tks[i].pi;
    773       // accumulate weighted z and weights for vertex update
    774       for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); ++k)
     986// Eliminate clusters with only one significant/unique track
     987bool VertexFinderDA4D::purge(vertex_t & vtx, tracks_t & tks, double & rho0, const double beta, double min_prob, double min_trk)
     988{
     989  const unsigned int nv = vtx.getSize();
     990  const unsigned int nt = tks.getSize();
     991
     992  if (nv < 2)
     993    return false;
     994
     995  double sumpmin = nt;
     996  unsigned int k0 = nv;
     997
     998  int nUnique = 0;
     999  double sump = 0;
     1000
     1001  double Z_init = rho0 * exp(-beta * fMuOutlayer * fMuOutlayer); // Add fDtCutOff here toghether  with this
     1002  vector<double> pk_exp_mBetaE = Compute_pk_exp_mBetaE(beta, vtx, tks, Z_init);
     1003
     1004  for (unsigned int k = 0; k < nv; ++k) {
     1005
     1006    nUnique = 0;
     1007    sump = 0;
     1008
     1009    double pmax = vtx.pk[k] / (vtx.pk[k] + rho0 * exp(-beta * fMuOutlayer* fMuOutlayer));
     1010    double pcut = min_prob * pmax;
     1011
     1012    for (unsigned int i = 0; i < nt; ++i) {
     1013      unsigned int idx = k*nt + i;
     1014
     1015      if(pk_exp_mBetaE[idx] == 0 || tks.Z[i] == 0)
    7751016      {
    776         k->se += tks[i].pi * k->ei / Zi;
    777         const double w = k->pk * tks[i].pi * k->ei / (Zi * (tks[i].dz2 * tks[i].dt2));
    778         k->sw += w;
    779         k->swz += w * tks[i].z;
    780         k->swt += w * tks[i].t;
    781         k->swE += w * Eik(tks[i], *k);
    782       }
    783     }
    784     else
    785     {
    786       sumpi += tks[i].pi;
    787     }
    788 
    789   } // end of track loop
    790 
    791   // now update z and pk
    792   double delta = 0;
    793   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    794   {
    795     if(k->sw > 0)
    796     {
    797       const double znew = k->swz / k->sw;
    798       const double tnew = k->swt / k->sw;
    799       delta += std::pow(k->z - znew, 2.) + std::pow(k->t - tnew, 2.);
    800       k->z = znew;
    801       k->t = tnew;
    802       k->Tc = 2. * k->swE / k->sw;
    803     }
    804     else
    805     {
    806       // cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl
    807       k->Tc = -1;
    808     }
    809 
    810     k->pk = k->pk * k->se / sumpi;
    811   }
    812 
    813   // return how much the prototypes moved
    814   return delta;
    815 }
    816 
    817 //------------------------------------------------------------------------------
    818 
    819 static double update2(double beta, vector<track_t> &tks, vector<vertex_t> &y, double &rho0, double dzCutOff)
    820 {
    821   // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise
    822   // returns the squared sum of changes of vertex positions
    823 
    824   unsigned int nt = tks.size();
    825 
    826   //initialize sums
    827   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    828   {
    829     k->sw = 0.;
    830     k->swz = 0.;
    831     k->swt = 0.;
    832     k->se = 0.;
    833     k->swE = 0.;
    834     k->Tc = 0.;
    835   }
    836 
    837   // loop over tracks
    838   for(unsigned int i = 0; i < nt; i++)
    839   {
    840 
    841     // update pik and Zi and Ti
    842     double Zi = rho0 * std::exp(-beta * (dzCutOff * dzCutOff)); // cut-off (eventually add finite size in time)
    843     //double Ti = 0.; // dt0*std::exp(-beta*fDtCutOff);
    844     for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    845     {
    846       k->ei = std::exp(-beta * Eik(tks[i], *k)); // cache exponential for one track at a time
    847       Zi += k->pk * k->ei;
    848     }
    849     tks[i].Z = Zi;
    850 
    851     // normalization
    852     if(tks[i].Z > 0)
    853     {
    854       // accumulate weighted z and weights for vertex update
    855       for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    856       {
    857         k->se += tks[i].pi * k->ei / Zi;
    858         double w = k->pk * tks[i].pi * k->ei / (Zi * (tks[i].dz2 * tks[i].dt2));
    859         k->sw += w;
    860         k->swz += w * tks[i].z;
    861         k->swt += w * tks[i].t;
    862         k->swE += w * Eik(tks[i], *k);
    863       }
    864     }
    865 
    866   } // end of track loop
    867 
    868   // now update z
    869   double delta = 0;
    870   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    871   {
    872     if(k->sw > 0)
    873     {
    874       const double znew = k->swz / k->sw;
    875       const double tnew = k->swt / k->sw;
    876       delta += std::pow(k->z - znew, 2.) + std::pow(k->t - tnew, 2.);
    877       k->z = znew;
    878       k->t = tnew;
    879       k->Tc = 2 * k->swE / k->sw;
    880     }
    881     else
    882     {
    883       // cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;
    884       k->Tc = 0;
    885     }
    886   }
    887 
    888   // return how much the prototypes moved
    889   return delta;
    890 }
    891 
    892 //------------------------------------------------------------------------------
    893 
    894 static bool merge(vector<vertex_t> &y)
    895 {
    896   // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise
    897 
    898   if(y.size() < 2) return false;
    899 
    900   for(vector<vertex_t>::iterator k = y.begin(); (k + 1) != y.end(); k++)
    901   {
    902     if(std::abs((k + 1)->z - k->z) < 1.e-3 && std::abs((k + 1)->t - k->t) < 1.e-3)
    903     { // with fabs if only called after freeze-out (splitAll() at highter T)
    904       double rho = k->pk + (k + 1)->pk;
    905       if(rho > 0)
    906       {
    907         k->z = (k->pk * k->z + (k + 1)->z * (k + 1)->pk) / rho;
    908         k->t = (k->pk * k->t + (k + 1)->t * (k + 1)->pk) / rho;
    909       }
    910       else
    911       {
    912         k->z = 0.5 * (k->z + (k + 1)->z);
    913         k->t = 0.5 * (k->t + (k + 1)->t);
    914       }
    915       k->pk = rho;
    916 
    917       y.erase(k + 1);
    918       return true;
    919     }
    920   }
    921 
    922   return false;
    923 }
    924 
    925 //------------------------------------------------------------------------------
    926 
    927 static bool merge(vector<vertex_t> &y, double &beta)
    928 {
    929   // merge clusters that collapsed or never separated,
    930   // only merge if the estimated critical temperature of the merged vertex is below the current temperature
    931   // return true if vertices were merged, false otherwise
    932   if(y.size() < 2) return false;
    933 
    934   for(vector<vertex_t>::iterator k = y.begin(); (k + 1) != y.end(); k++)
    935   {
    936     if(std::abs((k + 1)->z - k->z) < 2.e-3 && std::abs((k + 1)->t - k->t) < 2.e-3)
    937     {
    938       double rho = k->pk + (k + 1)->pk;
    939       double swE = k->swE + (k + 1)->swE - k->pk * (k + 1)->pk / rho * (std::pow((k + 1)->z - k->z, 2.) + std::pow((k + 1)->t - k->t, 2.));
    940       double Tc = 2 * swE / (k->sw + (k + 1)->sw);
    941 
    942       if(Tc * beta < 1)
    943       {
    944         if(rho > 0)
    945         {
    946           k->z = (k->pk * k->z + (k + 1)->z * (k + 1)->pk) / rho;
    947           k->t = (k->pk * k->t + (k + 1)->t * (k + 1)->pk) / rho;
    948         }
    949         else
    950         {
    951           k->z = 0.5 * (k->z + (k + 1)->z);
    952           k->t = 0.5 * (k->t + (k + 1)->t);
    953         }
    954         k->pk = rho;
    955         k->sw += (k + 1)->sw;
    956         k->swE = swE;
    957         k->Tc = Tc;
    958         y.erase(k + 1);
    959         return true;
    960       }
    961     }
    962   }
    963 
    964   return false;
    965 }
    966 
    967 //------------------------------------------------------------------------------
    968 
    969 static bool purge(vector<vertex_t> &y, vector<track_t> &tks, double &rho0, const double beta, const double dzCutOff)
    970 {
    971   // eliminate clusters with only one significant/unique track
    972   if(y.size() < 2) return false;
    973 
    974   unsigned int nt = tks.size();
    975   double sumpmin = nt;
    976   vector<vertex_t>::iterator k0 = y.end();
    977   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    978   {
    979     int nUnique = 0;
    980     double sump = 0;
    981     double pmax = k->pk / (k->pk + rho0 * exp(-beta * dzCutOff * dzCutOff));
    982     for(unsigned int i = 0; i < nt; i++)
    983     {
    984       if(tks[i].Z > 0)
    985       {
    986         double p = k->pk * std::exp(-beta * Eik(tks[i], *k)) / tks[i].Z;
    987         sump += p;
    988         if((p > 0.9 * pmax) && (tks[i].pi > 0))
    989         {
    990           nUnique++;
    991         }
    992       }
    993     }
    994 
    995     if((nUnique < 2) && (sump < sumpmin))
    996     {
     1017        continue;
     1018      }
     1019
     1020      double p = pk_exp_mBetaE[idx] / tks.Z[i];
     1021      sump += p;
     1022      if( ( p > pcut ) & ( tks.w[i] > 0 ) ) nUnique++;
     1023    }
     1024
     1025    if ((nUnique < min_trk) && (sump < sumpmin)) {
    9971026      sumpmin = sump;
    9981027      k0 = k;
    9991028    }
    1000   }
    1001 
    1002   if(k0 != y.end())
    1003   {
    1004     //cout << "eliminating prototype at " << k0->z << "," << k0->t << " with sump=" << sumpmin << endl;
    1005     //rho0+=k0->pk;
    1006     y.erase(k0);
     1029
     1030  }
     1031
     1032  if (k0 != nv) {
     1033    if (fVerbose > 5) {
     1034      std::cout  << Form("eliminating prototype at z = %.3f mm, t = %.0f ps", vtx.z[k0], vtx.t[k0]) << " with sump=" << sumpmin
     1035                 << "  rho*nt =" << vtx.pk[k0]*nt
     1036                 << endl;
     1037    }
     1038    vtx.removeItem(k0);
    10071039    return true;
    1008   }
    1009   else
    1010   {
     1040  } else {
    10111041    return false;
    10121042  }
    10131043}
    10141044
    1015 //------------------------------------------------------------------------------
    1016 
    1017 static double beta0(double betamax, vector<track_t> &tks, vector<vertex_t> &y, const double coolingFactor)
    1018 {
    1019 
    1020   double T0 = 0; // max Tc for beta=0
    1021   // estimate critical temperature from beta=0 (T=inf)
    1022   unsigned int nt = tks.size();
    1023 
    1024   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    1025   {
    1026 
    1027     // vertex fit at T=inf
    1028     double sumwz = 0.;
    1029     double sumwt = 0.;
    1030     double sumw = 0.;
    1031     for(unsigned int i = 0; i < nt; i++)
    1032     {
    1033       double w = tks[i].pi / (tks[i].dz2 * tks[i].dt2);
    1034       sumwz += w * tks[i].z;
    1035       sumwt += w * tks[i].t;
    1036       sumw += w;
    1037     }
    1038     k->z = sumwz / sumw;
    1039     k->t = sumwt / sumw;
    1040 
    1041     // estimate Tcrit, eventually do this in the same loop
    1042     double a = 0, b = 0;
    1043     for(unsigned int i = 0; i < nt; i++)
    1044     {
    1045       double dx = tks[i].z - (k->z);
    1046       double dt = tks[i].t - (k->t);
    1047       double w = tks[i].pi / (tks[i].dz2 * tks[i].dt2);
    1048       a += w * (std::pow(dx, 2.) / tks[i].dz2 + std::pow(dt, 2.) / tks[i].dt2);
    1049       b += w;
    1050     }
    1051     double Tc = 2. * a / b; // the critical temperature of this vertex
    1052     if(Tc > T0) T0 = Tc;
    1053   } // vertex loop (normally there should be only one vertex at beta=0)
    1054 
    1055   if(T0 > 1. / betamax)
    1056   {
    1057     return betamax / pow(coolingFactor, int(std::log(T0 * betamax) / std::log(coolingFactor)) - 1);
    1058   }
    1059   else
    1060   {
    1061     // ensure at least one annealing step
    1062     return betamax / coolingFactor;
    1063   }
    1064 }
    1065 
    1066 //------------------------------------------------------------------------------
    1067 
    1068 static bool split(double beta, vector<track_t> &tks, vector<vertex_t> &y)
    1069 {
    1070   // split only critical vertices (Tc >~ T=1/beta   <==>   beta*Tc>~1)
    1071   // an update must have been made just before doing this (same beta, no merging)
    1072   // returns true if at least one cluster was split
    1073 
    1074   const double epsilon = 1e-3; // split all single vertices by 10 um
    1075   bool split = false;
    1076 
    1077   // avoid left-right biases by splitting highest Tc first
    1078 
    1079   std::vector<std::pair<double, unsigned int> > critical;
    1080   for(unsigned int ik = 0; ik < y.size(); ik++)
    1081   {
    1082     if(beta * y[ik].Tc > 1.)
    1083     {
    1084       critical.push_back(make_pair(y[ik].Tc, ik));
    1085     }
    1086   }
    1087   std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int> >());
    1088 
    1089   for(unsigned int ic = 0; ic < critical.size(); ic++)
    1090   {
    1091     unsigned int ik = critical[ic].second;
    1092     // estimate subcluster positions and weight
    1093     double p1 = 0, z1 = 0, t1 = 0, w1 = 0;
    1094     double p2 = 0, z2 = 0, t2 = 0, w2 = 0;
    1095     //double sumpi=0;
    1096     for(unsigned int i = 0; i < tks.size(); i++)
    1097     {
    1098       if(tks[i].Z > 0)
    1099       {
    1100         //sumpi+=tks[i].pi;
    1101         double p = y[ik].pk * exp(-beta * Eik(tks[i], y[ik])) / tks[i].Z * tks[i].pi;
    1102         double w = p / (tks[i].dz2 * tks[i].dt2);
    1103         if(tks[i].z < y[ik].z)
    1104         {
    1105           p1 += p;
    1106           z1 += w * tks[i].z;
    1107           t1 += w * tks[i].t;
    1108           w1 += w;
    1109         }
    1110         else
    1111         {
    1112           p2 += p;
    1113           z2 += w * tks[i].z;
    1114           t2 += w * tks[i].t;
    1115           w2 += w;
    1116         }
    1117       }
    1118     }
    1119     if(w1 > 0)
    1120     {
    1121       z1 = z1 / w1;
    1122       t1 = t1 / w1;
     1045
     1046// -----------------------------------------------------------------------------
     1047// Plot status
     1048void VertexFinderDA4D::plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it, const char* flag)
     1049{
     1050  vector<int> vtx_color = {2,4,8,1,5,6,9,14,46,3};
     1051  while(vtx.getSize() > vtx_color.size()) vtx_color.push_back(40);
     1052
     1053  vector<double> t_PV, dt_PV, z_PV, dz_PV;
     1054  vector<double> t_PU, dt_PU, z_PU, dz_PU;
     1055
     1056  double ETot = 0;
     1057  vector<double> pk_exp_mBetaE = Compute_pk_exp_mBetaE(beta, vtx, tks, 0);
     1058
     1059  for(unsigned int i = 0; i < tks.getSize(); i++)
     1060  {
     1061    for(unsigned int k = 0; k < vtx.getSize(); k++)
     1062    {
     1063      unsigned int idx = k*tks.getSize() + i;
     1064      if(pk_exp_mBetaE[idx] == 0) continue;
     1065
     1066      double p_ygx = pk_exp_mBetaE[idx] / tks.Z[i];
     1067
     1068      ETot += tks.w[i] * p_ygx * Energy(tks.z[i], vtx.z[k], tks.dz2_o[i], tks.t[i], vtx.t[k], tks.dt2_o[i]);
     1069    }
     1070
     1071    if(tks.tt[i]->IsPU)
     1072    {
     1073      t_PU.push_back(tks.t[i]);
     1074      dt_PU.push_back(1./tks.dt_o[i]);
     1075      z_PU.push_back(tks.z[i]);
     1076      dz_PU.push_back(1./tks.dz_o[i]);
    11231077    }
    11241078    else
    11251079    {
    1126       z1 = y[ik].z - epsilon;
    1127       t1 = y[ik].t - epsilon;
    1128     }
    1129     if(w2 > 0)
    1130     {
    1131       z2 = z2 / w2;
    1132       t2 = t2 / w2;
    1133     }
    1134     else
    1135     {
    1136       z2 = y[ik].z + epsilon;
    1137       t2 = y[ik].t + epsilon;
    1138     }
    1139 
    1140     // reduce split size if there is not enough room
    1141     if((ik > 0) && (y[ik - 1].z >= z1))
    1142     {
    1143       z1 = 0.5 * (y[ik].z + y[ik - 1].z);
    1144       t1 = 0.5 * (y[ik].t + y[ik - 1].t);
    1145     }
    1146     if((ik + 1 < y.size()) && (y[ik + 1].z <= z2))
    1147     {
    1148       z2 = 0.5 * (y[ik].z + y[ik + 1].z);
    1149       t2 = 0.5 * (y[ik].t + y[ik + 1].t);
    1150     }
    1151 
    1152     // split if the new subclusters are significantly separated
    1153     if((z2 - z1) > epsilon || std::abs(t2 - t1) > epsilon)
    1154     {
    1155       split = true;
    1156       vertex_t vnew;
    1157       vnew.pk = p1 * y[ik].pk / (p1 + p2);
    1158       y[ik].pk = p2 * y[ik].pk / (p1 + p2);
    1159       vnew.z = z1;
    1160       vnew.t = t1;
    1161       y[ik].z = z2;
    1162       y[ik].t = t2;
    1163       y.insert(y.begin() + ik, vnew);
    1164 
    1165       // adjust remaining pointers
    1166       for(unsigned int jc = ic; jc < critical.size(); jc++)
    1167       {
    1168         if(critical[jc].second > ik)
    1169         {
    1170           critical[jc].second++;
    1171         }
    1172       }
    1173     }
    1174   }
    1175 
    1176   //  stable_sort(y.begin(), y.end(), clusterLessZ);
    1177   return split;
    1178 }
    1179 
    1180 //------------------------------------------------------------------------------
    1181 
    1182 void splitAll(vector<vertex_t> &y)
    1183 {
    1184 
    1185   const double epsilon = 1e-3; // split all single vertices by 10 um
    1186   const double zsep = 2 * epsilon; // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
    1187   const double tsep = 2 * epsilon; // check t as well
    1188 
    1189   vector<vertex_t> y1;
    1190 
    1191   for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
    1192   {
    1193     if(((k == y.begin()) || (k - 1)->z < k->z - zsep) && (((k + 1) == y.end()) || (k + 1)->z > k->z + zsep))
    1194     {
    1195       // isolated prototype, split
    1196       vertex_t vnew;
    1197       vnew.z = k->z - epsilon;
    1198       vnew.t = k->t - epsilon;
    1199       (*k).z = k->z + epsilon;
    1200       (*k).t = k->t + epsilon;
    1201       vnew.pk = 0.5 * (*k).pk;
    1202       (*k).pk = 0.5 * (*k).pk;
    1203       y1.push_back(vnew);
    1204       y1.push_back(*k);
    1205     }
    1206     else if(y1.empty() || (y1.back().z < k->z - zsep) || (y1.back().t < k->t - tsep))
    1207     {
    1208       y1.push_back(*k);
    1209     }
    1210     else
    1211     {
    1212       y1.back().z -= epsilon;
    1213       y1.back().t -= epsilon;
    1214       k->z += epsilon;
    1215       k->t += epsilon;
    1216       y1.push_back(*k);
    1217     }
    1218   } // vertex loop
    1219 
    1220   y = y1;
    1221 }
     1080      t_PV.push_back(tks.t[i]);
     1081      dt_PV.push_back(1./tks.dt_o[i]);
     1082      z_PV.push_back(tks.z[i]);
     1083      dz_PV.push_back(1./tks.dz_o[i]);
     1084    }
     1085  }
     1086
     1087
     1088  ETot /= tks.sum_w;
     1089  fEnergy_rec.push_back(ETot);
     1090  fBeta_rec.push_back(beta);
     1091  fNvtx_rec.push_back(vtx.getSize());
     1092
     1093  double t_min = TMath::Min(  TMath::MinElement(t_PV.size(), &t_PV[0]), TMath::MinElement(t_PU.size(), &t_PU[0])  );
     1094  t_min = TMath::Min(t_min, TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - fVertexTSize;
     1095  double t_max = TMath::Max(  TMath::MaxElement(t_PV.size(), &t_PV[0]), TMath::MaxElement(t_PU.size(), &t_PU[0])  );
     1096  t_max = TMath::Max(t_max, TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + fVertexTSize;
     1097
     1098  double z_min = TMath::Min(  TMath::MinElement(z_PV.size(), &z_PV[0]), TMath::MinElement(z_PU.size(), &z_PU[0])  );
     1099  z_min = TMath::Min(z_min, TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
     1100  double z_max = TMath::Max(  TMath::MaxElement(z_PV.size(), &z_PV[0]), TMath::MaxElement(z_PU.size(), &z_PU[0])  );
     1101  z_max = TMath::Max(z_max, TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
     1102
     1103  auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
     1104
     1105  TGraphErrors* gr_PVtks = new TGraphErrors(t_PV.size(), &t_PV[0], &z_PV[0], &dt_PV[0], &dz_PV[0]);
     1106  gr_PVtks->SetTitle(Form("Clustering space - #beta = %.6f", beta));
     1107  gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
     1108  gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
     1109  gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
     1110  gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
     1111  gr_PVtks->SetMarkerStyle(4);
     1112  gr_PVtks->SetMarkerColor(8);
     1113  gr_PVtks->SetLineColor(8);
     1114  gr_PVtks->Draw("APE1");
     1115
     1116  TGraphErrors* gr_PUtks = new TGraphErrors(t_PU.size(), &t_PU[0], &z_PU[0], &dt_PU[0], &dz_PU[0]);
     1117  gr_PUtks->SetMarkerStyle(3);
     1118  gr_PUtks->Draw("PE1");
     1119
     1120  TGraph* gr_vtx = new TGraph(vtx.getSize(), &(vtx.t[0]), &(vtx.z[0]));
     1121  gr_vtx->SetMarkerStyle(28);
     1122  gr_vtx->SetMarkerColor(2);
     1123  gr_vtx->SetMarkerSize(2.);
     1124  gr_vtx->Draw("PE1");
     1125
     1126  fItInputGenVtx->Reset();
     1127  TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
     1128  Candidate *candidate;
     1129  unsigned int k = 0;
     1130  while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
     1131  {
     1132    gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
     1133    k++;
     1134  }
     1135  gr_genvtx->SetMarkerStyle(33);
     1136  gr_genvtx->SetMarkerColor(6);
     1137  gr_genvtx->SetMarkerSize(2.);
     1138  gr_genvtx->Draw("PE1");
     1139
     1140  // auto leg = new TLegend(0.1, 0.1);
     1141  // leg->AddEntry(gr_PVtks, "PV tks", "ep");
     1142  // leg->AddEntry(gr_PUtks, "PU tks", "ep");
     1143  // leg->AddEntry(gr_vtx, "Cluster center", "p");
     1144  // leg->Draw();
     1145
     1146  c_2Dspace->SetGrid();
     1147  c_2Dspace->SaveAs(fFigFolderPath + Form("/c_2Dspace_beta%010.0f-%s%d.png", 1E7*beta, flag, n_it));
     1148
     1149  delete c_2Dspace;
     1150}
     1151
     1152// -----------------------------------------------------------------------------
     1153// Plot status at the end
     1154void VertexFinderDA4D::plot_status_end(vertex_t &vtx, tracks_t &tks)
     1155{
     1156  unsigned int nv = vtx.getSize();
     1157
     1158  // Define colors in a meaningfull way
     1159  vector<int> MyPalette(nv);
     1160
     1161  const int Number = 3;
     1162  double Red[Number]    = { 1.00, 0.00, 0.00};
     1163  double Green[Number]  = { 0.00, 1.00, 0.00};
     1164  double Blue[Number]   = { 1.00, 0.00, 1.00};
     1165  double Length[Number] = { 0.00, 0.50, 1.00 };
     1166  int FI = TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nv);
     1167  for (unsigned int i=0;i<nv;i++) MyPalette[i] = FI+i;
     1168
     1169  TCanvas * c_out = new TCanvas("c_out", "c_out", 800, 600);
     1170  double t_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 2*fVertexTSize;
     1171  double t_max = TMath::Max(TMath::MaxElement(tks.getSize(), &tks.t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 2*fVertexTSize;
     1172
     1173  double z_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 15;
     1174  double z_max = TMath::Max( TMath::MaxElement(tks.getSize(), &tks.z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 15;
     1175
     1176  // Draw tracks
     1177  for(unsigned int i = 0; i < tks.getSize(); i++)
     1178  {
     1179    double dt[] = {1./tks.dt_o[i]};
     1180    double dz[] = {1./tks.dz_o[i]};
     1181    TGraphErrors* gr = new TGraphErrors(1, &(tks.t[i]), &(tks.z[i]), dt, dz);
     1182
     1183    gr->SetNameTitle(Form("gr%d",i), Form("gr%d",i));
     1184
     1185    int marker = tks.tt[i]->IsPU? 1 : 4;
     1186    gr->SetMarkerStyle(marker);
     1187
     1188    int idx = tks.tt[i]->ClusterIndex;
     1189    int color = idx>=0 ? MyPalette[idx] : 13;
     1190    gr->SetMarkerColor(color);
     1191    gr->SetLineColor(color);
     1192
     1193    int line_style = idx>=0 ? 1 : 3;
     1194    gr->SetLineStyle(line_style);
     1195
     1196    if(i==0)
     1197    {
     1198      gr->SetTitle(Form("Clustering space - Tot Vertexes = %d", nv));
     1199      gr->GetXaxis()->SetTitle("t CA [ps]");
     1200      gr->GetXaxis()->SetLimits(t_min, t_max);
     1201      gr->GetYaxis()->SetTitle("z CA [mm]");
     1202      gr->GetYaxis()->SetRangeUser(z_min, z_max);
     1203      gr->Draw("APE1");
     1204    }
     1205    else gr->Draw("PE1");
     1206  }
     1207
     1208  // Draw vertices
     1209  for(unsigned int k = 0; k < vtx.getSize(); k++)
     1210  {
     1211    TGraph* gr = new TGraph(1, &(vtx.t[k]), &(vtx.z[k]));
     1212
     1213    gr->SetNameTitle(Form("grv%d",k), Form("grv%d",k));
     1214
     1215    gr->SetMarkerStyle(41);
     1216    gr->SetMarkerSize(2.);
     1217    gr->SetMarkerColor(MyPalette[k]);
     1218
     1219    gr->Draw("P");
     1220  }
     1221
     1222  fItInputGenVtx->Reset();
     1223  TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
     1224  TGraph* gr_genPV = new TGraph(1);
     1225  Candidate *candidate;
     1226  unsigned int k = 0;
     1227  while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
     1228  {
     1229    if(k == 0 ) {
     1230      gr_genPV->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
     1231    }
     1232    else gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
     1233
     1234    k++;
     1235  }
     1236  gr_genvtx->SetMarkerStyle(20);
     1237  gr_genvtx->SetMarkerColorAlpha(kBlack, 0.8);
     1238  gr_genvtx->SetMarkerSize(.8);
     1239  gr_genvtx->Draw("PE1");
     1240  gr_genPV->SetMarkerStyle(33);
     1241  gr_genPV->SetMarkerColorAlpha(kBlack, 1);
     1242  gr_genPV->SetMarkerSize(2.5);
     1243  gr_genPV->Draw("PE1");
     1244
     1245  // auto note =  new TLatex();
     1246  // note->DrawLatexNDC(0.5, 0.8, Form("#splitline{Vertexes Reco = %d }{Vertexes gen = %d}", vtx.getSize(), k) );
     1247
     1248  c_out->SetGrid();
     1249  c_out->SaveAs(fFigFolderPath + Form("/c_final.root"));
     1250  delete c_out;
     1251}
     1252
     1253// -----------------------------------------------------------------------------
     1254// Plot splitting
     1255void VertexFinderDA4D::plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx)
     1256{
     1257  vector<double> t, dt, z, dz;
     1258
     1259  for(unsigned int i = 0; i < tks.getSize(); i++)
     1260  {
     1261      t.push_back(tks.t[i]);
     1262      dt.push_back(1./tks.dt_o[i]);
     1263      z.push_back(tks.z[i]);
     1264      dz.push_back(1./tks.dz_o[i]);
     1265  }
     1266
     1267
     1268  double t_min = TMath::Min(TMath::MinElement(t.size(), &t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 50;
     1269  double t_max = TMath::Max(TMath::MaxElement(t.size(), &t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 50;
     1270
     1271  double z_min = TMath::Min(TMath::MinElement(z.size(), &z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
     1272  double z_max = TMath::Max(TMath::MaxElement(z.size(), &z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
     1273
     1274  auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
     1275
     1276  TGraphErrors* gr_PVtks = new TGraphErrors(t.size(), &t[0], &z[0], &dt[0], &dz[0]);
     1277  gr_PVtks->SetTitle(Form("Clustering space"));
     1278  gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
     1279  gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
     1280  gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
     1281  gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
     1282  gr_PVtks->SetMarkerStyle(4);
     1283  gr_PVtks->SetMarkerColor(1);
     1284  gr_PVtks->SetLineColor(1);
     1285  gr_PVtks->Draw("APE1");
     1286
     1287  TGraph* gr_vtx = new TGraph(1, &(vtx.t[i_vtx]), &(vtx.z[i_vtx]));
     1288  gr_vtx->SetMarkerStyle(28);
     1289  gr_vtx->SetMarkerColor(2);
     1290  gr_vtx->SetMarkerSize(2.);
     1291  gr_vtx->Draw("PE1");
     1292
     1293  double t_pos[] = {vtx.t[i_vtx], vtx.t[i_vtx]+100};
     1294  double t_neg[] = {vtx.t[i_vtx], vtx.t[i_vtx]-100};
     1295  double z_pos[] = {vtx.z[i_vtx], vtx.z[i_vtx]+(zn/tn)*100};
     1296  double z_neg[] = {vtx.z[i_vtx], vtx.z[i_vtx]-(zn/tn)*100};
     1297
     1298  TGraph* gr_pos = new TGraph(2, &t_pos[0], &z_pos[0]);
     1299  gr_pos->SetLineColor(8);
     1300  gr_pos->SetMarkerColor(8);
     1301  gr_pos->Draw("PL");
     1302  TGraph* gr_neg = new TGraph(2, &t_neg[0], &z_neg[0]);
     1303  gr_neg->SetLineColor(4);
     1304  gr_neg->SetMarkerColor(4);
     1305  gr_neg->Draw("PL");
     1306
     1307
     1308  c_2Dspace->SetGrid();
     1309  c_2Dspace->SaveAs(fFigFolderPath + Form("/crush_splitting.png"));
     1310
     1311  delete c_2Dspace;
     1312}

  • modules/VertexFinderDA4D.h

    rc614dd7 r84a097e  
    66 *  Cluster vertices from tracks using deterministic annealing and timing information
    77 *
    8  *  \authors M. Selvaggi, L. Gray
     8 *  Author O. Cerri
    99 *
    1010 */
    … …  
    1313
    1414#include <vector>
     15#include <iostream>
    1516
    1617class TObjArray;
    … …  
    2021class VertexFinderDA4D: public DelphesModule
    2122{
    22 public:
    23   VertexFinderDA4D();
    24   ~VertexFinderDA4D();
    25 
    26   void Init();
    27   void Process();
    28   void Finish();
    29 
    30   void clusterize(const TObjArray &tracks, TObjArray &clusters);
    31   std::vector<Candidate *> vertices();
    32 
    33 private:
    34   Bool_t fVerbose;
    35   Double_t fMinPT;
    36 
    37   Float_t fVertexSpaceSize;
    38   Float_t fVertexTimeSize;
    39   Bool_t fUseTc;
    40   Float_t fBetaMax;
    41   Float_t fBetaStop;
    42   Double_t fCoolingFactor;
    43   Int_t fMaxIterations;
    44   Double_t fDzCutOff;
    45   Double_t fD0CutOff;
    46   Double_t fDtCutOff; // for when the beamspot has time
    47 
    48   TObjArray *fInputArray;
    49   TIterator *fItInputArray;
    50 
    51   TObjArray *fOutputArray;
    52   TObjArray *fVertexOutputArray;
    53 
    54   ClassDef(VertexFinderDA4D, 1)
     23  public:
     24
     25    class tracks_t
     26    {
     27      public:
     28        std::vector<double> z;      // z-coordinate at point of closest approach to the beamline [mm]
     29        std::vector<double> r;      // z-coordinate at point of closest approach to the beamline [mm]
     30
     31        std::vector<double> t;      // t-coordinate at point of closest approach to the beamline  [ps]
     32
     33        std::vector<double> dz2_o;    //1 over the squared error of z(pca)
     34        std::vector<double> dz_o;
     35        std::vector<double> dt2_o;    //1 over the squared error of t(pca)
     36        std::vector<double> dt_o;
     37
     38        std::vector<Candidate*> tt; // a pointer to the Candidate Track
     39        std::vector<double> Z;      // Normalization of P(y|x) = p_k * exp(-beta*E_ik)/Z_i
     40
     41        std::vector<double> w;     // track weight
     42
     43        // std::vector<double> pt;
     44        // std::vector<double> pz;
     45
     46        std::vector<int> PID;
     47
     48        double sum_w_o_dz2 = 0;
     49        double sum_w_o_dt2 = 0;
     50        double sum_w = 0;
     51
     52        tracks_t(){}
     53        ~tracks_t(){}
     54
     55        void addItem( double new_z, double new_t, double new_dz2_o, double new_dt2_o, Candidate * new_tt, double new_w, int new_PID)
     56        {
     57          z.push_back( new_z );
     58          t.push_back( new_t );
     59          dz2_o.push_back( new_dz2_o );
     60          dz_o.push_back( sqrt(new_dz2_o) );
     61          dt2_o.push_back( new_dt2_o );
     62          dt_o.push_back( sqrt(new_dt2_o) );
     63          tt.push_back( new_tt );
     64
     65          w.push_back( new_w );
     66          Z.push_back( 1.0 );
     67
     68          PID.push_back(new_PID);
     69        }
     70
     71        unsigned int getSize()
     72        {
     73          return z.size();
     74        }
     75
     76    };
     77
     78    class vertex_t
     79    {
     80      public:
     81        std::vector<double> z;  //           z coordinate
     82        std::vector<double> t;  //           t coordinate
     83        std::vector<double> pk; //          vertex weight for "constrained" clustering
     84
     85        // Elements of the covariance matrix of the posterior distribution
     86        std::vector<double> szz, stt, stz;
     87        std::vector<double> beta_c;
     88
     89        double ZSize;
     90        double TSize;
     91
     92        // Used in the tracks-cluster assignment
     93        double sum_pt;
     94        std::vector<unsigned int> i_tracks;
     95
     96
     97        vertex_t(){}
     98        ~vertex_t(){}
     99
     100        void addItem( double new_z, double new_t, double new_pk)
     101        {
     102          z.push_back( new_z);
     103          t.push_back( new_t);
     104          pk.push_back( new_pk);
     105
     106          szz.push_back(0);
     107          stt.push_back(0);
     108          stz.push_back(0);
     109          beta_c.push_back(0);
     110        }
     111
     112        void removeItem( unsigned int i )
     113        {
     114          z.erase( z.begin() + i );
     115          t.erase( t.begin() + i );
     116          pk.erase( pk.begin() + i );
     117
     118          szz.erase(szz.begin() + i);
     119          stt.erase(stt.begin() + i);
     120          stz.erase(stz.begin() + i);
     121          beta_c.erase(beta_c.begin() + i);
     122        }
     123
     124        void mergeItems(unsigned int k1, unsigned int k2)
     125        {
     126          double new_pk = pk[k1] + pk[k2];
     127          double new_z = (z[k1]*pk[k1] + z[k2]*pk[k2])/new_pk;
     128          double new_t = (t[k1]*pk[k1] + t[k2]*pk[k2])/new_pk;
     129
     130          removeItem(k2);
     131          z[k1] = new_z;
     132          t[k1] = new_t;
     133          pk[k1] = new_pk;
     134        }
     135
     136        unsigned int getSize() const
     137        {
     138          return z.size();
     139        }
     140
     141        double DistanceSquare(unsigned int k1, unsigned int k2)
     142        {
     143          double dz = (z[k1] - z[k2])/ZSize;
     144          double dt = (t[k1] - t[k2])/TSize;
     145
     146          return dz*dz + dt*dt;
     147        }
     148
     149        unsigned int NearestCluster(double t_, double z_)
     150        {
     151          unsigned int k_min = 0;
     152          double d2_min = 0;
     153
     154          for(unsigned int k = 0; k < getSize(); k++)
     155          {
     156            double dt = (t[k] - t_)/TSize;
     157            double dz = (z[k] - z_)/ZSize;
     158            double d2 = dt*dt + dz*dz;
     159            if(k == 0 || d2 < d2_min)
     160            {
     161              k_min = k;
     162              d2_min = d2;
     163            }
     164          }
     165
     166          return k_min;
     167        }
     168
     169        std::pair<double, unsigned int> ComputeAllBeta_c(unsigned int verbose = 0)
     170        {
     171          unsigned int nv = getSize();
     172          unsigned int k_min = 0;
     173          double beta_c_min = 0;
     174          for(unsigned int k = 0; k < nv; k++)
     175          {
     176            if(szz[k] == 0 || stt[k] == 0)
     177            {
     178              beta_c[k] = 1000;
     179              if(verbose>5)
     180              {
     181                std::cout << std::endl << Form("Varianca too small! beta_c set to 1000\nVertex %d:  t=%1.2e, z=%1.2e, pk=%1.2e, szz=%1.2e, stt=%1.2e, stz=%1.2e", k, t[k], z[k], pk[k], szz[k], stt[k], stz[k]);
     182                // throw std::invalid_argument( "Attempting to compute beta c for uninitialized vertex" );
     183              }
     184            }
     185
     186            double aux = (stt[k] - szz[k])*(stt[k] - szz[k]) + 4*stz[k]*stz[k];
     187            aux = 1. / (stt[k] + szz[k] + sqrt(aux));
     188
     189            if(k == 0 || aux < beta_c_min)
     190            {
     191              beta_c_min = aux;
     192              k_min = k;
     193            }
     194
     195            beta_c[k] = aux;
     196          }
     197
     198          std::pair<double, unsigned int> out(beta_c_min, k_min);
     199          return out;
     200        }
     201    };
     202
     203    VertexFinderDA4D();
     204    ~VertexFinderDA4D();
     205
     206    void Init();
     207    void Process();
     208    void Finish();
     209
     210    void clusterize(TObjArray &clusters);
     211
     212    // Define the distance metric between tracks and vertices
     213    double Energy(double t_z, double v_z, double dz2_o, double t_t, double v_t, double dt2_o);
     214
     215    // Fill the tracks structure from the input array
     216    void fill(tracks_t &tks);
     217
     218    // Compute higher phase transition temperature
     219    double beta0(tracks_t & tks, vertex_t &vtx);
     220
     221    // Compute the new vertexes position and mass
     222    double update(double beta, tracks_t &tks, vertex_t &vtx, double rho0);
     223
     224    // If a vertex has beta_c lower than beta, split it
     225    bool split(double &beta,  vertex_t & vtx, tracks_t & tks);
     226
     227    // Merge vertexes closer than declared dimensions
     228    bool merge(vertex_t & vtx, double d2_merge);
     229
     230    // Eliminate clusters with only one significant/unique track
     231    bool purge(vertex_t & vtx, tracks_t & tks, double & rho0, const double beta, double min_prob, double min_trk);
     232
     233    // Compute all the energies and set the partition function normalization for each track
     234    std::vector<double> Compute_pk_exp_mBetaE(double beta, vertex_t &vtx, tracks_t &tks, double Z_init);
     235
     236    // Plot status of tracks and Vertices
     237    void plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it = 0, const char* flag ="");
     238
     239    // Plot status at the end of the fitting
     240    void plot_status_end(vertex_t &vtx, tracks_t &tks);
     241
     242    // Plot Crush
     243    void plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx);
     244
     245
     246  private:
     247
     248    unsigned int fVerbose;
     249
     250    UInt_t fMaxIterations;
     251    UInt_t fMaxVertexNumber;
     252
     253    Float_t fBetaMax;
     254    Float_t fBetaStop;
     255    Float_t fBetaPurge;
     256
     257
     258    Float_t fVertexZSize;
     259    Float_t fVertexTSize;
     260
     261    Double_t fCoolingFactor;
     262
     263    Double_t fDzCutOff;
     264    Double_t fD0CutOff;
     265    Double_t fDtCutOff;
     266    Double_t fPtMin;
     267    Double_t fPtMax;
     268
     269    Double_t fD2UpdateLim;
     270    Double_t fD2Merge;
     271    Double_t fMuOutlayer;
     272    Double_t fMinTrackProb;
     273    Int_t fMinNTrack;
     274
     275    TObjArray *fInputArray;
     276    TIterator *fItInputArray;
     277
     278    TObjArray *fTrackOutputArray;
     279    TObjArray *fVertexOutputArray;
     280
     281    ClassDef(VertexFinderDA4D, 1)
     282
     283    //Used to keep track of number of verices in verbose mode
     284    std::vector<double> fEnergy_rec;
     285    std::vector<double> fBeta_rec;
     286    std::vector<double> fNvtx_rec;
     287
     288    // Debug purpose only
     289    TObjArray *fInputGenVtx;
     290    TIterator *fItInputGenVtx;
     291
     292    TString fFigFolderPath = "";
    55293};
    56294

  • readers/DelphesPythia8.cpp

    rc614dd7 r84a097e  
    122122
    123123    pdgParticle = pdg->GetParticle(pid);
    124     candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge() / 3.0) : -999;
     124
     125    if(pdgParticle)
     126    {
     127      candidate->Charge = int(pdgParticle->Charge()/3.0);
     128    }
     129    else if(abs(pid) == 1000612) candidate->Charge = TMath::Sign(1, pid);
     130    else if(abs(pid) == 1000622) candidate->Charge = 0;
     131    else if(abs(pid) == 1000632) candidate->Charge = TMath::Sign(1, pid);
     132    else if(abs(pid) == 1000642) candidate->Charge = 0;
     133    else if(abs(pid) == 1000652) candidate->Charge = TMath::Sign(1, pid);
     134    else if(abs(pid) == 1006113) candidate->Charge = 0;
     135    else if(abs(pid) == 1006211) candidate->Charge = TMath::Sign(1, pid);
     136    else if(abs(pid) == 1006213) candidate->Charge = TMath::Sign(1, pid);
     137    else if(abs(pid) == 1006223) candidate->Charge = TMath::Sign(2, pid);
     138    else if(abs(pid) == 1006311) candidate->Charge = 0;
     139    else if(abs(pid) == 1006313) candidate->Charge = 0;
     140    else if(abs(pid) == 1006321) candidate->Charge = TMath::Sign(1, pid);
     141    else if(abs(pid) == 1006323) candidate->Charge = TMath::Sign(1, pid);
     142    else if(abs(pid) == 1006333) candidate->Charge = 0;
     143    else {
     144      candidate->Charge = -999;
     145      cout << "[WARNING] Unrecognised particle. PID: " << pid << endl;
     146    }
     147
    125148    candidate->Mass = mass;
    126149
Note: See TracChangeset for help on using the changeset viewer.