Fork me on GitHub

Changeset 53f4746 in git


Ignore:
Timestamp:
Apr 30, 2021, 4:30:36 PM (4 years ago)
Author:
Franco BEDESCHI <bed@…>
Branches:
master
Children:
a47edcc
Parents:
952bbbc
Message:

Cluster counting example and VertexFit update

Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

  • examples/Pythia8/ee_zh.cmnd

    r952bbbc r53f4746  
    11
    2 Main:numberOfEvents = 10000         ! number of events to generate
     2Main:numberOfEvents = 10         ! number of events to generate
    33
    44Beams:idA = 11                   ! first beam, e- = -11
    … …  
    77
    88! Vertex smearing :
    9 Beams:allowVertexSpread = on
    10 Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
    11 Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
    12 Beams:sigmaVertexZ = 0.64      !  0.64 mm
    13 ï¿ŒBeams:sigmaVertexT = 0.64      !  0.64 mm
     9Beams:allowVertexSpread = off
     10!Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
     11!Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
     12!Beams:sigmaVertexZ = 20.64      !  0.64 mm
     13!ï¿ŒBeams:sigmaVertexT = 0.64      !  0.64 mm
    1414ￌ
    1515
    … …  
    2020! 5) Force the Z decays to muons
    212123:onMode = off
    22 23:onIfAny = 11 -11
     2223:onIfAny = 13 -13

  • external/TrackCovariance/ObsTrk.cc

    r952bbbc r53f4746  
    1414ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)
    1515{
    16         SetBfield(B);
     16        SetB(B);
    1717        fGC = GC;
    1818        fGenX = x;
    … …  
    5656ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC)
    5757{
    58         SetBfield(B);
     58        SetB(B);
    5959        fGC = GC;
    6060        fGenX.SetXYZ(x[0],x[1],x[2]);

  • external/TrackCovariance/TrkUtil.cc

    r952bbbc r53f4746  
    22#include <iostream>
    33#include <algorithm>
     4#include <TSpline.h>
    45
    56// Constructor
    … …  
    9697        if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl;
    9798        //
    98         return ParToP(Par,fBz);
     99        return ParToP(Par, fBz);
    99100}
    100101//
    … …  
    233234}
    234235//
    235 
    236 //
    237 void TrkUtil::SetBfield(Double_t Bz)
    238 {
    239         fBz = Bz;
    240 }
    241 
    242236// Setup chamber volume
    243237void TrkUtil::SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax)
    … …  
    274268                //      << ", C= " << C << ", z0= " << z0 << ", ct= " << ct << std::endl;
    275269                //
    276                 // Track length per unit phase change
    277                 Double_t Scale = sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C));
     270                // Track length per unit phase change 
     271                Double_t Scale = sqrt(1.0 + ct * ct) / (2.0 * TMath::Abs(C));
    278272                //
    279273                // Find intersections with chamber boundaries
    280274                //
    281                 Double_t phRin = 0.0;                   // phase of inner cylinder
    282                 Double_t phRin2= 0.0;                   // phase of inner cylinder intersection (2nd branch)
     275                Double_t phRin = 0.0;                   // phase of inner cylinder 
     276                Double_t phRin2 = 0.0;                  // phase of inner cylinder intersection (2nd branch)
    283277                Double_t phRhi = 0.0;                   // phase of outer cylinder intersection
    284278                Double_t phZmn = 0.0;                   // phase of left wall intersection
    285279                Double_t phZmx = 0.0;                   // phase of right wall intersection
    286280                //  ... with inner cylinder
    287                 Double_t Rtop = TMath::Abs((1.0 + C*D) / C);
     281                Double_t Rtop = TMath::Abs((1.0 + C * D) / C);
    288282
    289283                if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment ***
    290284                {
    291                         Double_t ph = 2 * asin(C*sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D)));
    292                         Double_t z = z0 + ct*ph / (2.0*C);
     285                        Double_t ph = 2 * asin(C * sqrt((fRmin * fRmin - D * D) / (1.0 + 2.0 * C * D)));
     286                        Double_t z = z0 + ct * ph / (2.0 * C);
    293287
    294288                        //std::cout << "Rin intersection: ph = " << ph<<", z= "<<z << std::endl;
    295289
    296                         if (z < fZmax && z > fZmin)     phRin = TMath::Abs(ph); // Intersection inside chamber volume
     290                        if (z < fZmax && z > fZmin)     phRin = TMath::Abs(ph); // Intersection inside chamber volume   
    297291                        //
    298292                        // Include second branch of loopers
    299293                        Double_t Pi = 3.14159265358979323846;
    300                         Double_t ph2 = 2*Pi - TMath::Abs(ph);
     294                        Double_t ph2 = 2 * Pi - TMath::Abs(ph);
    301295                        if (ph < 0)ph2 = -ph2;
    302296                        z = z0 + ct * ph2 / (2.0 * C);
    … …  
    306300                if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment ***
    307301                {
    308                         Double_t ph = 2 * asin(C*sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D)));
    309                         Double_t z = z0 + ct*ph / (2.0*C);
    310                         if (z < fZmax && z > fZmin)     phRhi = TMath::Abs(ph); // Intersection inside chamber volume
     302                        Double_t ph = 2 * asin(C * sqrt((fRmax * fRmax - D * D) / (1.0 + 2.0 * C * D)));
     303                        Double_t z = z0 + ct * ph / (2.0 * C);
     304                        if (z < fZmax && z > fZmin)     phRhi = TMath::Abs(ph); // Intersection inside chamber volume   
    311305                }
    312306                //  ... with left wall
    … …  
    314308                if (Zdir > 0.0)
    315309                {
    316                         Double_t ph = 2.0*C*Zdir;
    317                         Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C));
    318                         if (Rint < fRmax && Rint > fRmin)       phZmn = TMath::Abs(ph); // Intersection inside chamber volume
     310                        Double_t ph = 2.0 * C * Zdir;
     311                        Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C));
     312                        if (Rint < fRmax && Rint > fRmin)       phZmn = TMath::Abs(ph); // Intersection inside chamber volume   
    319313                }
    320314                //  ... with right wall
    … …  
    322316                if (Zdir > 0.0)
    323317                {
    324                         Double_t ph = 2.0*C*Zdir;
    325                         Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C));
    326                         if (Rint < fRmax && Rint > fRmin)       phZmx = TMath::Abs(ph); // Intersection inside chamber volume
     318                        Double_t ph = 2.0 * C * Zdir;
     319                        Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C));
     320                        if (Rint < fRmax && Rint > fRmin)       phZmx = TMath::Abs(ph); // Intersection inside chamber volume   
    327321                }
    328322                //
    … …  
    342336                {
    343337                        dPhase = ph_arr[iPos + 2] - ph_arr[iPos + 1];
    344                         tLength = dPhase*Scale;
     338                        tLength = dPhase * Scale;
    345339                }
    346340        }
    … …  
    349343//
    350344// Return number of ionization clusters
    351 Bool_t TrkUtil::IonClusters(Double_t &Ncl, Double_t mass, TVectorD Par)
     345Bool_t TrkUtil::IonClusters(Double_t& Ncl, Double_t mass, TVectorD Par)
    352346{
    353347        //
    … …  
    386380                        TVector3 p = ParToP(Par);
    387381                        bg = p.Mag() / mass;
    388                         muClu = Nclusters(bg)*tLen;                             // Avg. number of clusters
     382                        muClu = Nclusters(bg) * tLen;                           // Avg. number of clusters
    389383
    390384                        Ncl = gRandom->PoissonD(muClu);                 // Actual number of clusters
    … …  
    392386
    393387        }
    394 //
     388        //
    395389        return Signal;
    396390}
    … …  
    413407        //
    414408        //
    415         /*
    416         std::vector<double> bg{ 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,
    417         12., 15., 20., 50., 100., 200., 500., 1000. };
    418         // He 90 - Isobutane 10
    419         std::vector<double> ncl_He_Iso{ 42.94, 23.6,18.97,12.98,12.2,12.13,
    420         12.24,12.73,13.03,13.29,13.63,14.08,15.56,16.43,16.8,16.95,16.98 };
    421         //
    422         // pure He
    423         std::vector<double> ncl_He{ 11.79,6.5,5.23,3.59,3.38,3.37,3.4,3.54,3.63,
    424         3.7,3.8,3.92,4.33,4.61,4.78,4.87,4.89 };
    425         //
    426         // Argon 50 - Ethane 50
    427         std::vector<double> ncl_Ar_Eth{ 130.04,71.55,57.56,39.44,37.08,36.9,
    428         37.25,38.76,39.68,40.49,41.53,42.91,46.8,48.09,48.59,48.85,48.93 };
    429         //
    430         // pure Argon
    431         std::vector<double> ncl_Ar{ 88.69,48.93,39.41,27.09,25.51,25.43,25.69,
    432         26.78,27.44,28.02,28.77,29.78,32.67,33.75,34.24,34.57,34.68 };
    433         //
    434         Int_t nPoints = (Int_t)bg.size();
    435         bg.push_back(10000.);
    436         std::vector<double> ncl;
    437         switch (Opt)
    438         {
    439         case 0: ncl = ncl_He_Iso;                       // He-Isobutane
    440                 break;
    441         case 1: ncl = ncl_He;                           // pure He
    442                 break;
    443         case 2: ncl = ncl_Ar_Eth;                       // Argon - Ethane
    444                 break;
    445         case 3: ncl = ncl_Ar;                           // pure Argon
    446                 break;
    447         }
    448         ncl.push_back(ncl[nPoints - 1]);
    449         */
    450409        const Int_t Npt = 18;
    451410        Double_t bg[Npt] = { 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,
    … …  
    469428        //
    470429        Double_t ncl[Npt];
    471         switch (Opt)
    472         {
    473                 case 0: std::copy(ncl_He_Iso, ncl_He_Iso + Npt, ncl);   // He-Isobutane
     430        switch (Opt)
     431        {
     432        case 0: std::copy(ncl_He_Iso, ncl_He_Iso + Npt, ncl);   // He-Isobutane
    474433                break;
    475                 case 1: std::copy(ncl_He, ncl_He + Npt, ncl);           // pure He
     434        case 1: std::copy(ncl_He, ncl_He + Npt, ncl);           // pure He
    476435                break;
    477                 case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl);   // Argon - Ethane
     436        case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl);   // Argon - Ethane
    478437                break;
    479                 case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl);           // pure Argon
     438        case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl);           // pure Argon
    480439                break;
    481         }
    482         //
    483         Int_t ilow = 0;
    484         while (begam > bg[ilow])ilow++;
    485         ilow--;
    486         //std::cout << "ilow= " << ilow << ", low = " << bg[ilow] << ", val = " << begam
    487         //      << ", high = " << bg[ilow + 1] << std::endl;
    488         //
    489         Int_t ind[3] = { ilow, ilow + 1, ilow + 2 };
    490         TVectorD y(3);
    491         for (Int_t i = 0; i < 3; i++)y(i) = ncl[ind[i]];
    492         TVectorD x(3);
    493         for (Int_t i = 0; i < 3; i++)x(i) = bg[ind[i]];
    494         TMatrixD Xval(3, 3);
    495         for (Int_t i = 0; i < 3; i++)Xval(i, 0) = 1.0;
    496         for (Int_t i = 0; i < 3; i++)Xval(i, 1) = x(i);
    497         for (Int_t i = 0; i < 3; i++)Xval(i, 2) = x(i) * x(i);
    498         //std::cout << "Xval:" << std::endl; Xval.Print();
    499         Xval.Invert();
    500         TVectorD coeff = Xval * y;
    501         Double_t interp = coeff[0] + coeff[1] * begam + coeff[2] * begam * begam;
    502         //std::cout << "val1= (" <<x(0)<<", "<< y(0) << "), val2= ("
    503         //      <<x(1)<<", "<< y(1) << "), val3= ("
    504         //      <<x(2)<<", "<< y(2)
    505         //      << "), result= (" <<begam<<", "<< interp<<")" << std::endl;
    506         //
    507         //if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << std::endl;
    508         if (begam < bg[0]) interp = 0.0;
    509         //std::cout << "bg= " << begam << ", Opt= " << Opt <<", interp = "<<interp<< std::endl;
    510         return 100*interp;
    511 }
    512 //
    513 Double_t TrkUtil::funcNcl(Double_t *xp, Double_t *par){
     440        }
     441        //
     442        Double_t interp = 0.0;
     443        TSpline3* sp3 = new TSpline3("sp3", bg, ncl, Npt);
     444        if (begam > bg[0] && begam < bg[Npt - 1]) interp = sp3->Eval(begam);
     445        if(begam < bg[0]) interp = bg[0];
     446        if(begam > bg[Npt-1]) interp = bg[Npt-1];
     447        return 100 * interp;
     448}
     449//
     450Double_t TrkUtil::funcNcl(Double_t* xp, Double_t* par) {
    514451        Double_t bg = xp[0];
    515452        return Nclusters(bg);

  • external/TrackCovariance/TrkUtil.h

    r952bbbc r53f4746  
    2525        // Service routines
    2626        //
     27        void SetB(Double_t Bz){ fBz = Bz; }
    2728        TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q);
    2829        TVector3 ParToP(TVectorD Par);
    … …  
    7071        // Cluster counting in gas
    7172        //
    72         void SetBfield(Double_t Bz);
     73        void SetBfield(Double_t Bz){ fBz = Bz; }
    7374        // Define gas volume (units = meters)
    7475        void SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax);

  • external/TrackCovariance/VertexFit.cc

    r952bbbc r53f4746  
    1818        fxCst.ResizeTo(3);
    1919        fCovCst.ResizeTo(3, 3);
     20        fCovCstInv.ResizeTo(3, 3);
    2021        fXv.ResizeTo(3);
    2122        fcovXv.ResizeTo(3, 3);
    … …  
    3132        fxCst.ResizeTo(3);
    3233        fCovCst.ResizeTo(3, 3);
     34        fCovCstInv.ResizeTo(3, 3);
    3335        fXv.ResizeTo(3);
    3436        fcovXv.ResizeTo(3, 3);
    … …  
    3840                TVectorD pr = *trkPar[i];
    3941                fPar.push_back(new TVectorD(pr));
     42                fParNew.push_back(new TVectorD(pr));
    4043                TMatrixDSym cv = *trkCov[i];
    4144                fCov.push_back(new TMatrixDSym(cv));
     45                fCovNew.push_back(new TMatrixDSym(cv));
    4246        }
    4347        fChi2List.ResizeTo(fNtr);
    … …  
    5458        fxCst.ResizeTo(3);
    5559        fCovCst.ResizeTo(3, 3);
     60        fCovCstInv.ResizeTo(3, 3);
    5661        fXv.ResizeTo(3);
    5762        fcovXv.ResizeTo(3, 3);
    … …  
    6166        {
    6267                fPar.push_back(new TVectorD(track[i]->GetObsPar()));
     68                fParNew.push_back(new TVectorD(track[i]->GetObsPar()));
    6369                fCov.push_back(new TMatrixDSym(track[i]->GetCov()));
     70                fCovNew.push_back(new TMatrixDSym(track[i]->GetCov()));
    6471        }
    6572}
    … …  
    7380        {
    7481                if (fx0i[i])  { fx0i[i]->Clear();               delete fx0i[i]; }
    75                 if (fai[i])   { fai[i]->Clear();                        delete fai[i]; }
    76                 if (fDi[i])   { fDi[i]->Clear();                        delete fDi[i]; }
    77                 if (fWi[i])   { fWi[i]->Clear();                        delete fWi[i]; }
     82                if (fai[i])   { fai[i]->Clear();                delete fai[i]; }
     83                if (fdi[i])   { fdi[i]->Clear();                delete fdi[i]; }
     84                if (fAti[i])  { fAti[i]->Clear();               delete fAti[i]; }
     85                if (fDi[i])   { fDi[i]->Clear();                delete fDi[i]; }
     86                if (fWi[i])   { fWi[i]->Clear();                delete fWi[i]; }
    7887                if (fWinvi[i]){ fWinvi[i]->Clear();     delete fWinvi[i]; }
    7988        }
    … …  
    8190        fx0i.clear();
    8291        fai.clear();
     92        fdi.clear();
     93        fAti.clear();
    8394        fDi.clear();
    8495        fWi.clear();
    … …  
    8899{
    89100        fxCst.Clear();         
    90         fCovCst.Clear();               
     101        fCovCst.Clear();
     102        fCovCstInv.Clear();
    91103        fXv.Clear();           
    92104        fcovXv.Clear();         
    … …  
    95107        for (Int_t i = 0; i < fNtr; i++)
    96108        {
    97                 fPar[i]->Clear();       delete fPar[i];
    98                 fCov[i]->Clear();       delete fCov[i];
     109                fPar[i]->Clear();               delete fPar[i];
     110                fParNew[i]->Clear();    delete fParNew[i];
     111                fCov[i]->Clear();               delete fCov[i];
     112                fCovNew[i]->Clear();    delete fCovNew[i];
    99113        }
    100114        fPar.clear();
    101         fCov.clear();   
     115        fParNew.clear();
     116        fCov.clear();
     117        fCovNew.clear();               
    102118        //
    103119        ResetWrkArrays();
    … …  
    424440        //
    425441        // Get track parameters and their covariance
    426         TVectorD par = *fPar[i];
     442        TVectorD par = *fParNew[i];
    427443        TMatrixDSym Cov = *fCov[i];
    428444        //
    … …  
    433449        //
    434450        TMatrixD A = Fill_A(par, fs);                           // A = dx/da = derivatives wrt track parameters
     451        TMatrixD At(TMatrixD::kTransposed, A);          // A transposed
     452        fAti.push_back(new TMatrixD(At));                       // Store A'
     453        TVectorD di = A * (par - *fPar[i]);             // x-shift
     454        fdi.push_back(new TVectorD(di));                        // Store x-shift
    435455        TMatrixDSym Winv = Cov.Similarity(A);           // W^-1 = A*C*A'
    436         fWinvi.push_back(new TMatrixDSym(Winv));                // Store W^-1 matrix
     456        fWinvi.push_back(new TMatrixDSym(Winv));        // Store W^-1 matrix
    437457        //
    438458        TMatrixDSym W = RegInv(Winv);                           // W = (A*C*A')^-1
    … …  
    456476{
    457477        //std::cout << "VertexFitter: just in" << std::endl;
    458         if (fNtr < 2)
    459         {
     478        if (fNtr < 2 && !fVtxCst){
    460479                std::cout << "VertexFit::VertexFitter - Method called with less than 2 tracks - Aborting " << std::endl;
    461480                std::exit(1);
    … …  
    469488        Double_t Chi2 = 0;
    470489        TVectorD x0 = fXv;      // If previous fit done
    471         if (fRold < 0.0)for (Int_t i = 0; i < 3; i++)x0(i) = 1000.;     // Set to arbitrary large value if not
     490        if (fRold < 0.0) {
     491                // External constraint
     492                if (fVtxCst) fRold = TMath::Sqrt(fxCst(0) * fxCst(0) + fxCst(1) * fxCst(1));
     493                // No constraint
     494                else for (Int_t i = 0; i < 3; i++)x0(i) = 1000.;        // Set to arbitrary large value
     495        }
    472496        //
    473497        // Starting vertex radius approximation
    … …  
    475499        Double_t R = fRold;                                             // Use previous fit if available
    476500        if (R < 0.0) R = StartRadius();                 // Rough vertex estimate
     501        //std::cout << "Start radius: " << R << std::endl;
    477502        //
    478503        // Iteration properties
    … …  
    529554                        H += Ds;
    530555                        // update constant term
    531                         TVectorD xs = *fx0i[i];
     556                        TVectorD xs = *fx0i[i] - *fdi[i];
     557                        TVectorD xx0 = *fx0i[i];
     558                        /*
     559                        std::cout << "Iter. " << Ntry << ", trk " << i << ", x= "
     560                                << xx0(0) << ", " << xx0(1) << ", " << xx0(2)<<
     561                                ", ph0= "<<par(1)<< std::endl;
     562                        */
    532563                        cterm += Ds * xs;
    533564                }                               // End loop on tracks
     565                // Some additions in case of external constraints
     566                if (fVtxCst) {
     567                        H += fCovCstInv;
     568                        cterm += fCovCstInv * fxCst;
     569                        DW1D += fCovCstInv;
     570                }
    534571                //
    535572                // update vertex position
    … …  
    544581                for (Int_t i = 0; i < fNtr; i++)
    545582                {
    546                         TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x);
     583                        TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x - *fdi[i]);
    547584                        TMatrixDSym Wm1 = *fWinvi[i];
    548585                        fChi2List(i) = Wm1.Similarity(lambda);
    … …  
    551588                        TVectorD b = (*fWi[i]) * (x - (*fx0i[i]));
    552589                        for (Int_t j = 0; j < 3; j++)ffi[i] += a(j) * b(j) / fa2i[i];
     590                        TVectorD newPar = *fPar[i] - ((*fCov[i]) * (*fAti[i])) * lambda;
     591                        fParNew[i] = new TVectorD(newPar);
    553592                }
     593                // Add external constraint to Chi2
     594                if (fVtxCst) Chi2 += fCovCstInv.Similarity(x - fxCst);
    554595                //
    555596                TVectorD dx = x - x0;
    … …  
    562603                // Store result
    563604                //
    564                 fXv = x;                                // Vertex position
     605                fXv = x;                        // Vertex position
    565606                fcovXv = covX;          // Vertex covariance
    566607                fChi2 = Chi2;           // Vertex fit Chi2
     608                //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2) << std::endl;
    567609        }               // end of iteration loop
    568610        //
    … …  
    603645void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov)  // Add gaussian vertex constraint
    604646{
    605         std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
     647        //std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
     648        fVtxCst = kTRUE;                                // Vertex constraint flag
     649        fxCst = xv;                                             // Constraint value
     650        fCovCst = cov;                                  // Constraint covariance
     651        fCovCstInv = cov;
     652        fCovCstInv.Invert();                            // Its inverse
     653        //
     654        // Set starting vertex as external constraint
     655        fXv = fxCst;
     656        fcovXv = fCovCst;
    606657}
    607658//
    … …  
    613664        fPar.push_back(par);                    // add new track
    614665        fCov.push_back(Cov);
     666        fParNew.push_back(par);                 // add new track
     667        fCovNew.push_back(Cov);
    615668        //
    616669        // Reset previous vertex temp arrays
    … …  
    628681        fPar.erase(fPar.begin() + iTrk);                // Remove track
    629682        fCov.erase(fCov.begin() + iTrk);
     683        fParNew.erase(fParNew.begin() + iTrk);          // Remove track
     684        fCovNew.erase(fCovNew.begin() + iTrk);
    630685        //
    631686        // Reset previous vertex temp arrays

  • external/TrackCovariance/VertexFit.h

    r952bbbc r53f4746  
    2323        Int_t fNtr;                                     // Number of tracks
    2424        std::vector<TVectorD*> fPar;            // Input parameter array
    25         std::vector<TMatrixDSym*> fCov;// Input parameter covariances
     25        std::vector<TVectorD*> fParNew;         // Updated parameter array
     26        std::vector<TMatrixDSym*> fCov;         // Input parameter covariances
     27        std::vector<TMatrixDSym*> fCovNew;      // Updated parameter covariances
    2628        // Constraints
    2729        Bool_t fVtxCst;                         // Vertex constraint flag
    2830        TVectorD fxCst;                         // Constraint value
    29         TMatrixDSym fCovCst;                    // Constraint covariance
     31        TMatrixDSym fCovCst;                    // Constraint
     32        TMatrixDSym fCovCstInv;                 // Inverse of constraint covariance
    3033        //
    3134        // Results
    … …  
    3942        // Work arrays
    4043        std::vector<Double_t> ffi;                              // Fit phases
    41         std::vector<TVectorD*> fx0i;                            // Track expansion points
     44        std::vector<TVectorD*> fx0i;                    // Track expansion points
    4245        std::vector<TVectorD*> fai;                             // dx/dphi
     46        std::vector<TVectorD*> fdi;                             // x-shift
    4347        std::vector<Double_t> fa2i;                             // a'Wa
     48        std::vector<TMatrixD*> fAti;                    // A transposed
    4449        std::vector<TMatrixDSym*> fDi;                  // W-WBW
    4550        std::vector<TMatrixDSym*> fWi;                  // (ACA')^-1
Note: See TracChangeset for help on using the changeset viewer.