Changeset 7dac4ea in git for examples/ExamplePVtxFit.C

Timestamp:

May 17, 2021, 6:05:45 PM (3 years ago)

Author:

GitHub <noreply@…>

Branches:

master

Children:

1c1c9c2, 33a6b3a, f8e61b2

Parents:

4afb18d (diff), 4acf2fd (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

git-author:

Michele Selvaggi <michele.selvaggi@…> (05/17/21 18:05:45)

git-committer:

GitHub <noreply@…> (05/17/21 18:05:45)

Message:

Merge pull request #95 from fbedesch/master

Fix for backward tracks and vertex fitting improvements

File:

• examples/ExamplePVtxFit.C (modified) (6 diffs)

examples/ExamplePVtxFit.C

@@ -14 +14 @@
 
 #endif
+
 
 
@@ -45 +46 @@
                 // Load selected branches with data from specified event
                 treeReader->ReadEntry(entry);
+                // 
                 Int_t Ntr = 0;  // # of tracks from primary vertex
                 Int_t NtrG = branchTrack->GetEntries();
                 TVectorD** pr = new TVectorD * [NtrG];
                 TMatrixDSym** cv = new TMatrixDSym * [NtrG];
+                //
+                // True vertex
+                Double_t xpv, ypv, zpv;
                 // If event contains at least 1 track
                 //
@@ -60 +65 @@
                                 // Get associated generated particle
                                 GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
+                                TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);
+                                TVector3 pg(gp->Px,gp->Py,gp->Pz);
+                                Double_t Q = (Double_t)gp->Charge;
+                                Double_t Bz = 2.0;
+                                TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz);
+                                TVectorD genPar = TrkUtil::ParToMm(genParM);
                                 //
                                 // Position of origin in mm
@@ -65 +76 @@
                                 Double_t y = gp->Y;
                                 Double_t z = gp->Z;
+                                Bool_t prim = kTRUE;    // Is primary?
+                                Int_t mp = gp->M1;      // Mother
+                                while(mp>0){
+                                        GenParticle* gm = 
+                                        (GenParticle*)branchGenPart->At(mp);
+                                        Double_t xm = gm->X;
+                                        Double_t ym = gm->Y;
+                                        Double_t zm = gm->Z;
+                                        if(x!=xm || y!=ym || z!=zm){
+                                                prim = kFALSE;
+                                                break;
+                                        }else mp = gm->M1;
+                                }
                                 //
                                 // group tracks originating from the primary vertex
-                                if (x == 0.0 && y == 0.0)
+                                if (prim)
                                 {
+                                        //std::cout<<"Event: "<<entry<<", track: "<<it;
+                                        //std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl; 
+                                        xpv = x;
+                                        ypv = y;
+                                        zpv = z; 
                                         //
                                         // Reconstructed track parameters
@@ -75 +104 @@
                                         Double_t obsC = trk->C;
                                         Double_t obsZ0 = trk->DZ;
+                                        //std::cout<<"Z0 track = "<< obsZ0
+                                        //<<", gen Z0 = "<<genPar(3)
+                                        //<<", gen cot = "<<genPar(4)<<std::endl;
                                         Double_t obsCtg = trk->CtgTheta;
                                         Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
@@ -96 +128 @@
                         Double_t Chi2 = Vtx->GetVtxChi2();
                         Double_t Ndof = 2 * (Double_t)Ntr - 3;
-                        Double_t PullX = xvtx(0) / TMath::Sqrt(covX(0, 0));
-                        Double_t PullY = xvtx(1) / TMath::Sqrt(covX(1, 1));
-                        Double_t PullZ = xvtx(2) / TMath::Sqrt(covX(2, 2));
+                        Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
+                        Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
+                        Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
+                        //std::cout<<"**** True  vertex (x, y, z) "<<xpv<<", "
+                        //<<ypv<<", "<<zpv<<std::endl;
+                        //std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
+                        //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
                         //
                         // Fill histograms