Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/VertexFinderDA4D.h

    r341014c 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
Note: See TracChangeset for help on using the changeset viewer.