Fork me on GitHub

Changes in / [7dac4ea:4afb18d] in git


Ignore:
Files:
2 deleted
6 edited

Legend:

Unmodified
Added
Removed
  • examples/ExamplePVtxFind.C

    r7dac4ea r4afb18d  
    5151                // Load selected branches with data from specified event
    5252                treeReader->ReadEntry(entry);
     53                Int_t Ntr = 0;  // # of starting tracks used in primary vertex
    5354                Int_t NtrG = branchTrack->GetEntries();
    54                 if(entry%500 ==0)std::cout << "Event "<<entry<<" opened containing " << NtrG << " tracks" << std::endl;
     55                //std::cout << "Event opened containing " << NtrG << " tracks" << std::endl;
    5556                TVectorD** pr = new TVectorD * [NtrG];
    5657                TMatrixDSym** cv = new TMatrixDSym * [NtrG];
     
    6263                //
    6364                Double_t Nprim = 0.0;
    64                 Double_t xpv = 0.0;             // Init true primary vertex position
    65                 Double_t ypv = 0.0;
    66                 Double_t zpv = 0.0;
    6765                if (branchTrack->GetEntries() > 0)
    6866                {
     
    8280                                //std::cout << "Got track parameters for track " << it << std::endl;
    8381                                //
    84                                 // Load all tracks for vertex fit
    85                                 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
    86                                 TVectorD obsPar(5, oPar);       // Fill observed parameters
    87                                 pr[it] = new TVectorD(obsPar);
    88                                 cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
    89                                 //std::cout << "Track loaded it= " << it << std::endl;
    90                         //
    91                         // Find true primary vertex
     82                                // Load tracks for vertex fit if impact parameters is not ridiculous
     83                                Double_t Dmax = 1.0;    // max is 1 mm
     84                                if (TMath::Abs(obsD0) < Dmax) {
     85                                        Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
     86                                        TVectorD obsPar(5, oPar);       // Fill observed parameters
     87                                        pr[Ntr] = new TVectorD(obsPar);
     88                                        cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
     89                                        Ntr++;
     90                                        //std::cout << "Track loaded Ntr= " << Ntr << std::endl;
     91                                }
     92                                //
     93                                // Get associated generated particle
    9294                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
    9395                                //std::cout << "GenParticle pointer "<<gp << std::endl;
     
    9799                                Double_t y = gp->Y;
    98100                                Double_t z = gp->Z;
    99                                 Bool_t prim = kTRUE;    // Is primary?
    100                                 Int_t mp = gp->M1;      // Mother
    101                                 while (mp > 0) {
    102                                         GenParticle* gm =
    103                                                 (GenParticle*)branchGenPart->At(mp);
    104                                         Double_t xm = gm->X;
    105                                         Double_t ym = gm->Y;
    106                                         Double_t zm = gm->Z;
    107                                         if (x != xm || y != ym || z != zm) {
    108                                                 prim = kFALSE;
    109                                                 break;
    110                                         }
    111                                         else mp = gm->M1;
    112                                 }
    113                                 if (prim) {             // It's a primary track
    114                                         Nprim++;
    115                                         xpv = x;        // Store true primary
    116                                         ypv = y;
    117                                         zpv = z;
    118                                 }
     101                                //std::cout << "Got position of origin " << std::endl;
     102                                //
     103                                // Count tracks originating from the primary vertex
     104                                if (x == 0.0 && y == 0.0) Nprim++;
    119105
    120106                        }               // End loop on tracks
    121107                }
    122108                //
    123                 // Find primary vertex
    124                 //
    125                 //Beam constraint
    126                 TVectorD xpvc(3);
    127                 xpvc(0) = 1.0;
    128                 xpvc(1) = -2.0;
    129                 xpvc(2) = 10.0;
    130                 TMatrixDSym covpvc(3); covpvc.Zero();
    131                 covpvc(0, 0) = 0.0097 * 0.0097;
    132                 covpvc(1, 1) = 2.55e-05 * 2.55e-05;
    133                 covpvc(2, 2) = 0.64 * 0.64;
    134                 //
    135                 //
    136                 // Skim tracks
    137                 Int_t nSkim = 0;
    138                 Int_t* nSkimmed = new Int_t[NtrG];
    139                 VertexFit* Vskim = new VertexFit();
    140                 Vskim->AddVtxConstraint(xpvc, covpvc);
    141                 Double_t MaxChi2 = 2.0+3*sqrt(2.);
    142                 for (Int_t n = 0; n < NtrG; n++) {
    143                         Vskim->AddTrk(pr[n], cv[n]);
    144                         Double_t Chi2One = Vskim->GetVtxChi2();
    145                         //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
    146                         if (Chi2One < MaxChi2) {
    147                                 nSkimmed[nSkim] = n;
    148                                 //std::cout << "nSkimmed[" << nSkim << "] = " << n << std::endl;
    149                                 nSkim++;
    150                         }
    151                         Vskim->RemoveTrk(0);
    152                 }
    153                 // Load tracks for primary fit
    154                 Int_t MinTrk = 1;       // Minumum # tracks for vertex fit
    155                 TVectorD** PrFit = new TVectorD * [nSkim];
    156                 TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
    157                 if (nSkim >= MinTrk) {
    158                         for (Int_t n = 0; n < nSkim; n++) {
    159                                 PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
    160                                 CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
    161                                 TVectorD test = *PrFit[n];
    162                                 //std::cout << "Test *PrFit[" << n << "](0) = " << test(0) << std::endl;
    163                         }
    164                 }
    165                 delete[] nSkimmed;
    166                 delete Vskim;
    167                 Int_t Nfound = nSkim;
    168                 if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
    169                 if (nSkim >= MinTrk) {
    170                         VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
     109                // Fit primary vertex
     110                //
     111                Int_t Nfound = Ntr;
     112                //std::cout << "Found tracks "<<Nfound << std::endl;
     113                Int_t MinTrk = 2;       // Minumum # tracks for vertex fit
     114                Double_t MaxChi2 = 6.;
     115                if (Ntr >= MinTrk) {
     116                        VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
    171117                        //std::cout << "Vertex fit created " << std::endl;
    172                         Vtx->AddVtxConstraint(xpvc, covpvc);
    173118                        //
    174119                        // Remove tracks with large chi2
    175                         Double_t MaxChi2Fit = 9.;
    176120                        Bool_t Done = kFALSE;
    177                         while (!Done) {
     121                        while (!Done){
    178122                                //std::cout << "After while " << std::endl;
    179123                                // Find largest Chi2 contribution
     
    186130                                Double_t Chi2Mx = Chi2L[iMax];
    187131                                //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
    188                                 if (Chi2Mx > MaxChi2Fit && Nfound > 2) {
    189                                         //std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
    190                                         Vtx->RemoveTrk(iMax);
    191                                         //std::cout << "After remove." << std::endl;
     132                                if (Chi2Mx > MaxChi2 && Nfound > 2) {
     133                                        Vtx->RemoveTrk(iMax);
    192134                                        Nfound--;
    193135                                }
     
    203145                        //
    204146                        // Require minimum number of tracks in vertex
    205                         Int_t Nmin = 2;
     147                        Int_t Nmin = 4;
    206148                        if (Nfound >= Nmin) {
    207149                                TVectorD xvtx = Vtx->GetVtx();
     
    209151                                TMatrixDSym covX = Vtx->GetVtxCov();
    210152                                Double_t Chi2 = Vtx->GetVtxChi2();
    211                                 Double_t Ndof = 2 * (Double_t)Nfound;
    212                                 Double_t PullX = (xvtx(0) - xpv) / TMath::Sqrt(covX(0, 0));
    213                                 Double_t PullY = (xvtx(1) - ypv) / TMath::Sqrt(covX(1, 1));
    214                                 Double_t PullZ = (xvtx(2) - zpv) / TMath::Sqrt(covX(2, 2));
     153                                Double_t Ndof = 2 * (Double_t)Nfound - 3;
     154                                Double_t PullX = xvtx(0) / TMath::Sqrt(covX(0, 0));
     155                                Double_t PullY = xvtx(1) / TMath::Sqrt(covX(1, 1));
     156                                Double_t PullZ = xvtx(2) / TMath::Sqrt(covX(2, 2));
    215157                                //
    216158                                // Fill histograms
     
    231173                //
    232174                // Cleanup
    233                 for (Int_t i = 0; i < NtrG; i++) delete pr[i];
    234                 for (Int_t i = 0; i < NtrG; i++) delete cv[i];
    235                 for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
    236                 for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
     175                for (Int_t i = 0; i < Ntr; i++) delete pr[i];
     176                for (Int_t i = 0; i < Ntr; i++) delete cv[i];
    237177                delete[] pr;
    238178                delete[] cv;
    239                 delete[] PrFit;
    240                 delete[] CvFit;
    241179        }
    242180        //
     
    257195        CnvN->cd(1);
    258196        hTrPrim->Draw();
    259         CnvN->cd(2);
     197        CnvN-> cd(2);
    260198        hTrFound->SetLineColor(kRed);
    261199        hTrFound->Draw();
  • examples/ExamplePVtxFit.C

    r7dac4ea r4afb18d  
    1414
    1515#endif
    16 
    1716
    1817
     
    4645                // Load selected branches with data from specified event
    4746                treeReader->ReadEntry(entry);
    48                 //
    4947                Int_t Ntr = 0;  // # of tracks from primary vertex
    5048                Int_t NtrG = branchTrack->GetEntries();
    5149                TVectorD** pr = new TVectorD * [NtrG];
    5250                TMatrixDSym** cv = new TMatrixDSym * [NtrG];
    53                 //
    54                 // True vertex
    55                 Double_t xpv, ypv, zpv;
    5651                // If event contains at least 1 track
    5752                //
     
    6560                                // Get associated generated particle
    6661                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
    67                                 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);
    68                                 TVector3 pg(gp->Px,gp->Py,gp->Pz);
    69                                 Double_t Q = (Double_t)gp->Charge;
    70                                 Double_t Bz = 2.0;
    71                                 TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz);
    72                                 TVectorD genPar = TrkUtil::ParToMm(genParM);
    7362                                //
    7463                                // Position of origin in mm
     
    7665                                Double_t y = gp->Y;
    7766                                Double_t z = gp->Z;
    78                                 Bool_t prim = kTRUE;    // Is primary?
    79                                 Int_t mp = gp->M1;      // Mother
    80                                 while(mp>0){
    81                                         GenParticle* gm =
    82                                         (GenParticle*)branchGenPart->At(mp);
    83                                         Double_t xm = gm->X;
    84                                         Double_t ym = gm->Y;
    85                                         Double_t zm = gm->Z;
    86                                         if(x!=xm || y!=ym || z!=zm){
    87                                                 prim = kFALSE;
    88                                                 break;
    89                                         }else mp = gm->M1;
    90                                 }
    9167                                //
    9268                                // group tracks originating from the primary vertex
    93                                 if (prim)
     69                                if (x == 0.0 && y == 0.0)
    9470                                {
    95                                         //std::cout<<"Event: "<<entry<<", track: "<<it;
    96                                         //std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl;
    97                                         xpv = x;
    98                                         ypv = y;
    99                                         zpv = z;
    10071                                        //
    10172                                        // Reconstructed track parameters
     
    10475                                        Double_t obsC = trk->C;
    10576                                        Double_t obsZ0 = trk->DZ;
    106                                         //std::cout<<"Z0 track = "<< obsZ0
    107                                         //<<", gen Z0 = "<<genPar(3)
    108                                         //<<", gen cot = "<<genPar(4)<<std::endl;
    10977                                        Double_t obsCtg = trk->CtgTheta;
    11078                                        Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
     
    12896                        Double_t Chi2 = Vtx->GetVtxChi2();
    12997                        Double_t Ndof = 2 * (Double_t)Ntr - 3;
    130                         Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
    131                         Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
    132                         Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
    133                         //std::cout<<"**** True  vertex (x, y, z) "<<xpv<<", "
    134                         //<<ypv<<", "<<zpv<<std::endl;
    135                         //std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
    136                         //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
     98                        Double_t PullX = xvtx(0) / TMath::Sqrt(covX(0, 0));
     99                        Double_t PullY = xvtx(1) / TMath::Sqrt(covX(1, 1));
     100                        Double_t PullZ = xvtx(2) / TMath::Sqrt(covX(2, 2));
    137101                        //
    138102                        // Fill histograms
  • examples/Pythia8/ee_zh.cmnd

    r7dac4ea r4afb18d  
    11
    2 ! number of events to generate
    32Main:numberOfEvents = 1000         ! number of events to generate
    43
     
    87
    98! 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:sigmaTime    = 0.64      !  0.64 mm
     14ï¿Œ
    1015
    11 Beams:allowVertexSpread = on
    12 Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
    13 Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
    14 Beams:sigmaVertexZ = 0.64      !  0.64 mm
    15 Beams:sigmaTime    = 0.64      !  0.64 mm
    1616
    1717! Higgsstrahlung process
     
    2020! 5) Force the Z decays to muons
    212123:onMode = off
    22 23:onIfAny = 13 -13
     2223:onIfAny = 11 -11
  • external/TrackCovariance/ObsTrk.cc

    r7dac4ea r4afb18d  
    1414ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)
    1515{
    16         SetB(B);
     16        SetBfield(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         SetB(B);
     58        SetBfield(B);
    5959        fGC = GC;
    6060        fGenX.SetXYZ(x[0],x[1],x[2]);
  • external/TrackCovariance/VertexFit.cc

    r7dac4ea r4afb18d  
    1818        fxCst.ResizeTo(3);
    1919        fCovCst.ResizeTo(3, 3);
    20         fCovCstInv.ResizeTo(3, 3);
    2120        fXv.ResizeTo(3);
    2221        fcovXv.ResizeTo(3, 3);
     
    3231        fxCst.ResizeTo(3);
    3332        fCovCst.ResizeTo(3, 3);
    34         fCovCstInv.ResizeTo(3, 3);
    3533        fXv.ResizeTo(3);
    3634        fcovXv.ResizeTo(3, 3);
     
    4038                TVectorD pr = *trkPar[i];
    4139                fPar.push_back(new TVectorD(pr));
    42                 fParNew.push_back(new TVectorD(pr));
    4340                TMatrixDSym cv = *trkCov[i];
    4441                fCov.push_back(new TMatrixDSym(cv));
    45                 fCovNew.push_back(new TMatrixDSym(cv));
    4642        }
    4743        fChi2List.ResizeTo(fNtr);
     
    5854        fxCst.ResizeTo(3);
    5955        fCovCst.ResizeTo(3, 3);
    60         fCovCstInv.ResizeTo(3, 3);
    6156        fXv.ResizeTo(3);
    6257        fcovXv.ResizeTo(3, 3);
     
    6661        {
    6762                fPar.push_back(new TVectorD(track[i]->GetObsPar()));
    68                 fParNew.push_back(new TVectorD(track[i]->GetObsPar()));
    6963                fCov.push_back(new TMatrixDSym(track[i]->GetCov()));
    70                 fCovNew.push_back(new TMatrixDSym(track[i]->GetCov()));
    7164        }
    7265}
     
    8073        {
    8174                if (fx0i[i])  { fx0i[i]->Clear();               delete fx0i[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]; }
     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]; }
    8778                if (fWinvi[i]){ fWinvi[i]->Clear();     delete fWinvi[i]; }
    8879        }
     
    9081        fx0i.clear();
    9182        fai.clear();
    92         fdi.clear();
    93         fAti.clear();
    9483        fDi.clear();
    9584        fWi.clear();
     
    9988{
    10089        fxCst.Clear();         
    101         fCovCst.Clear();
    102         fCovCstInv.Clear();
     90        fCovCst.Clear();               
    10391        fXv.Clear();           
    10492        fcovXv.Clear();         
     
    10795        for (Int_t i = 0; i < fNtr; i++)
    10896        {
    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];
     97                fPar[i]->Clear();       delete fPar[i];
     98                fCov[i]->Clear();       delete fCov[i];
    11399        }
    114100        fPar.clear();
    115         fParNew.clear();
    116         fCov.clear();
    117         fCovNew.clear();               
     101        fCov.clear();   
    118102        //
    119103        ResetWrkArrays();
     
    440424        //
    441425        // Get track parameters and their covariance
    442         TVectorD par = *fParNew[i];
     426        TVectorD par = *fPar[i];
    443427        TMatrixDSym Cov = *fCov[i];
    444428        //
     
    449433        //
    450434        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
    455435        TMatrixDSym Winv = Cov.Similarity(A);           // W^-1 = A*C*A'
    456         fWinvi.push_back(new TMatrixDSym(Winv));        // Store W^-1 matrix
     436        fWinvi.push_back(new TMatrixDSym(Winv));                // Store W^-1 matrix
    457437        //
    458438        TMatrixDSym W = RegInv(Winv);                           // W = (A*C*A')^-1
     
    476456{
    477457        //std::cout << "VertexFitter: just in" << std::endl;
    478         if (fNtr < 2 && !fVtxCst){
     458        if (fNtr < 2)
     459        {
    479460                std::cout << "VertexFit::VertexFitter - Method called with less than 2 tracks - Aborting " << std::endl;
    480461                std::exit(1);
     
    488469        Double_t Chi2 = 0;
    489470        TVectorD x0 = fXv;      // If previous fit done
    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         }
     471        if (fRold < 0.0)for (Int_t i = 0; i < 3; i++)x0(i) = 1000.;     // Set to arbitrary large value if not
    496472        //
    497473        // Starting vertex radius approximation
     
    499475        Double_t R = fRold;                                             // Use previous fit if available
    500476        if (R < 0.0) R = StartRadius();                 // Rough vertex estimate
    501         //std::cout << "Start radius: " << R << std::endl;
    502477        //
    503478        // Iteration properties
     
    554529                        H += Ds;
    555530                        // update constant term
    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                         */
     531                        TVectorD xs = *fx0i[i];
    563532                        cterm += Ds * xs;
    564533                }                               // 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                 }
    571534                //
    572535                // update vertex position
     
    581544                for (Int_t i = 0; i < fNtr; i++)
    582545                {
    583                         TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x - *fdi[i]);
     546                        TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x);
    584547                        TMatrixDSym Wm1 = *fWinvi[i];
    585548                        fChi2List(i) = Wm1.Similarity(lambda);
     
    588551                        TVectorD b = (*fWi[i]) * (x - (*fx0i[i]));
    589552                        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);
    592553                }
    593                 // Add external constraint to Chi2
    594                 if (fVtxCst) Chi2 += fCovCstInv.Similarity(x - fxCst);
    595554                //
    596555                TVectorD dx = x - x0;
     
    603562                // Store result
    604563                //
    605                 fXv = x;                        // Vertex position
     564                fXv = x;                                // Vertex position
    606565                fcovXv = covX;          // Vertex covariance
    607566                fChi2 = Chi2;           // Vertex fit Chi2
    608                 //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2) << std::endl;
    609567        }               // end of iteration loop
    610568        //
     
    645603void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov)  // Add gaussian vertex constraint
    646604{
    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;
     605        std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
    657606}
    658607//
     
    664613        fPar.push_back(par);                    // add new track
    665614        fCov.push_back(Cov);
    666         fParNew.push_back(par);                 // add new track
    667         fCovNew.push_back(Cov);
    668615        //
    669616        // Reset previous vertex temp arrays
     
    681628        fPar.erase(fPar.begin() + iTrk);                // Remove track
    682629        fCov.erase(fCov.begin() + iTrk);
    683         fParNew.erase(fParNew.begin() + iTrk);          // Remove track
    684         fCovNew.erase(fCovNew.begin() + iTrk);
    685630        //
    686631        // Reset previous vertex temp arrays
  • external/TrackCovariance/VertexFit.h

    r7dac4ea r4afb18d  
    2323        Int_t fNtr;                                     // Number of tracks
    2424        std::vector<TVectorD*> fPar;            // Input parameter array
    25         std::vector<TVectorD*> fParNew;         // Updated parameter array
    26         std::vector<TMatrixDSym*> fCov;         // Input parameter covariances
    27         std::vector<TMatrixDSym*> fCovNew;      // Updated parameter covariances
     25        std::vector<TMatrixDSym*> fCov;// Input parameter covariances
    2826        // Constraints
    2927        Bool_t fVtxCst;                         // Vertex constraint flag
    3028        TVectorD fxCst;                         // Constraint value
    31         TMatrixDSym fCovCst;                    // Constraint
    32         TMatrixDSym fCovCstInv;                 // Inverse of constraint covariance
     29        TMatrixDSym fCovCst;                    // Constraint covariance
    3330        //
    3431        // Results
     
    4239        // Work arrays
    4340        std::vector<Double_t> ffi;                              // Fit phases
    44         std::vector<TVectorD*> fx0i;                    // Track expansion points
     41        std::vector<TVectorD*> fx0i;                            // Track expansion points
    4542        std::vector<TVectorD*> fai;                             // dx/dphi
    46         std::vector<TVectorD*> fdi;                             // x-shift
    4743        std::vector<Double_t> fa2i;                             // a'Wa
    48         std::vector<TMatrixD*> fAti;                    // A transposed
    4944        std::vector<TMatrixDSym*> fDi;                  // W-WBW
    5045        std::vector<TMatrixDSym*> fWi;                  // (ACA')^-1
Note: See TracChangeset for help on using the changeset viewer.