Changeset b750b0a in git for examples/ExamplePVtxFind.C

Timestamp:

Jan 19, 2022, 10:11:39 AM (3 years ago)

Author:

Franco BEDESCHI <bed@…>

Branches:

master

Children:

78ce8d1, 7bca620

Parents:

bd376e3

Message:

Fix CovarianceMatrix scaling error. Vertexing improvements.

File:

• examples/ExamplePVtxFind.C (modified) (16 diffs)

examples/ExamplePVtxFind.C

@@ -1 +1 @@
 /*
 Example of using vertex fitter class to fit primary vertex
-assumed to be generated in (0,0,0)
+assumed to be generated with Pythia8/ee_zh_smr-shf.cmn
 */
 
@@ -44 +44 @@
         TH1D* hTrFoundG = new TH1D("hTrFoundG", "Found GOOD primary tracks", 41, -0.5, 40.5);
         TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD  primary tracks", 41, -0.5, 40.5);
+        TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5);
         //
         // Loop over all events
@@ -87 +88 @@
                                 pr[it] = new TVectorD(obsPar);
                                 cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
+                                // Check covariance
+                                /*
+                                TMatrixDSymEigen Eign(*cv[it]);
+                                TVectorD lambda = Eign.GetEigenValues();
+                                Bool_t good = kTRUE;
+                                for(Int_t i=0; i<5; i++) if(lambda(i)<0.0) good = kFALSE;
+                                if(!good){
+                                        std::cout<<"Cov["<<it<<"]"<<" eigenvalues: "
+                                        <<lambda(0)<<", "<<lambda(1)<<", "<<lambda(2)<<", "
+                                        <<lambda(3)<<", "<<lambda(4)<<std::endl;
+                                        TMatrixDSym Ctmp = *cv[it];
+                                        std::cout<<"Error D (dir - cov) "<<trk->ErrorD0
+                                        <<" - "<<TMath::Sqrt(Ctmp(0,0))<<std::endl;
+                                }
+                                */
                                 //std::cout << "Track loaded it= " << it << std::endl;
                         //
@@ -120 +136 @@
                         }               // End loop on tracks
                 }
+                std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;
                 //
                 // Find primary vertex
@@ -132 +149 @@
                 covpvc(1, 1) = 2.55e-05 * 2.55e-05;
                 covpvc(2, 2) = 0.64 * 0.64;
+                if(Nprim == 0){
+                        xpv = xpvc(0);
+                        ypv = xpvc(1);
+                        zpv = xpvc(2);
+                }
                 //
                 //
@@ -137 +159 @@
                 Int_t nSkim = 0;
                 Int_t* nSkimmed = new Int_t[NtrG];
-                VertexFit* Vskim = new VertexFit();
-                Vskim->AddVtxConstraint(xpvc, covpvc);
-                Double_t MaxChi2 = 2.0+3*sqrt(2.);
+                TVectorD** PrSk = new TVectorD * [1];
+                TMatrixDSym** CvSk = new TMatrixDSym * [1];
+                Double_t MaxChi2 = 2.0+3*2.;
                 for (Int_t n = 0; n < NtrG; n++) {
-                        Vskim->AddTrk(pr[n], cv[n]);
+                        PrSk[0] = new TVectorD(*pr[n]);
+                        CvSk[0] = new TMatrixDSym(*cv[n]);
+                        VertexFit* Vskim = new VertexFit(1,PrSk, CvSk);
+                        Vskim->AddVtxConstraint(xpvc, covpvc);
                         Double_t Chi2One = Vskim->GetVtxChi2();
                         //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
@@ -149 +174 @@
                                 nSkim++;
                         }
-                        Vskim->RemoveTrk(0);
-                }
+                        // Cleanup
+                        delete Vskim;
+                }
+                delete PrSk[0];
+                delete CvSk[0];
+                delete[] PrSk;
+                delete[] CvSk;
+                //
                 // Load tracks for primary fit
                 Int_t MinTrk = 1;       // Minumum # tracks for vertex fit
-                TVectorD** PrFit = new TVectorD * [nSkim];
-                TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
                 if (nSkim >= MinTrk) {
+                        TVectorD** PrFit = new TVectorD * [nSkim];
+                        TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
                         for (Int_t n = 0; n < nSkim; n++) {
                                 PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
                                 CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
-                                TVectorD test = *PrFit[n];
-                                //std::cout << "Test *PrFit[" << n << "](0) = " << test(0) << std::endl;
-                        }
-                }
-                delete[] nSkimmed;
-                delete Vskim;
-                Int_t Nfound = nSkim;
-                if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
-                if (nSkim >= MinTrk) {
+                        }
+                        delete[] nSkimmed;
+                        Int_t Nfound = nSkim;
+                        if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
+                        const Int_t MaxFound = 100; Double_t Chi2LL[MaxFound]; Double_t *Chi2L = Chi2LL;
                         VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
                         //std::cout << "Vertex fit created " << std::endl;
@@ -173 +200 @@
                         //
                         // Remove tracks with large chi2
-                        Double_t MaxChi2Fit = 9.;
+                        Double_t MaxChi2Fit = 4.5;
                         Bool_t Done = kFALSE;
                         while (!Done) {
@@ -180 +207 @@
                                 TVectorD Chi2List = Vtx->GetVtxChi2List();      // Get contributions to Chi2
                                 //std::cout << "After Chi2List.  " << std::endl; Chi2List.Print();
-                                Double_t* Chi2L = new Double_t[Nfound];
+                                //Double_t* Chi2L = new Double_t[Nfound];
                                 Chi2L = Chi2List.GetMatrixArray();
                                 Int_t iMax = TMath::LocMax(Nfound, Chi2L);
@@ -186 +213 @@
                                 Double_t Chi2Mx = Chi2L[iMax];
                                 //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
-                                if (Chi2Mx > MaxChi2Fit && Nfound > 2) {
+                                if (Chi2Mx > MaxChi2Fit && Nfound > 1) {
                                         //std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
                                         Vtx->RemoveTrk(iMax);
@@ -195 +222 @@
                                         Done = kTRUE;
                                 }
-                                //std::cout << "Done =  " << Done << std::endl;
-                                //delete[] Chi2L;
-                                //std::cout << "Array Chi2L removed " << std::endl;
                         }
                         //
@@ -203 +227 @@
                         //
                         // Require minimum number of tracks in vertex
-                        Int_t Nmin = 2;
+                        Int_t Nmin = 1;
                         if (Nfound >= Nmin) {
                                 TVectorD xvtx = Vtx->GetVtx();
-                                //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
+                                std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
                                 TMatrixDSym covX = Vtx->GetVtxCov();
                                 Double_t Chi2 = Vtx->GetVtxChi2();
@@ -222 +246 @@
                                 hTrPrim->Fill(Nprim);
                                 hTrFound->Fill((Double_t)Nfound);
+                                hTrDiff->Fill((Double_t)Nfound-Nprim);
                                 //std::cout << "Histograms filled " << std::endl;
                         }
                         //
+                        // Clean
                         delete Vtx;
+                        for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
+                        for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
+                        delete[] PrFit;
+                        delete[] CvFit;
                 }
 
@@ -233 +263 @@
                 for (Int_t i = 0; i < NtrG; i++) delete pr[i];
                 for (Int_t i = 0; i < NtrG; i++) delete cv[i];
-                for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
-                for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
                 delete[] pr;
                 delete[] cv;
-                delete[] PrFit;
-                delete[] CvFit;
         }
         //
@@ -254 +280 @@
         //
         TCanvas* CnvN = new TCanvas("CnvN", "Primary tracks found", 100, 100, 900, 500);
-        CnvN->Divide(2, 1);
+        CnvN->Divide(2, 2);
         CnvN->cd(1);
         hTrPrim->Draw();
@@ -260 +286 @@
         hTrFound->SetLineColor(kRed);
         hTrFound->Draw();
+        CnvN->cd(3);
+        hTrDiff->Draw();
 }