Changeset a47edcc in git

Timestamp:

May 1, 2021, 2:15:02 PM (4 years ago)

Author:

Franco BEDESCHI <bed@…>

Branches:

master

Children:

46d3442

Parents:

53f4746

Message:

Fixed issue with backward generated tracks

Files:

• examples/ExamplePVtxFit.C (modified) (6 diffs)
• external/TrackCovariance/TrkUtil.cc (modified) (3 diffs)
• external/TrackCovariance/TrkUtil.h (modified) (3 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

external/TrackCovariance/TrkUtil.cc

@@ -46 +46 @@
         Double_t C = a / (2 * pt);                      // Half curvature
         //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl;
-        Double_t r2 = x.Perp2();
+        Double_t r2 = x(0) * x(0) + x(1) * x(1);
         Double_t cross = x(0) * p(1) - x(1) * p(0);
         Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2);
@@ -61 +61 @@
         Double_t st = asin(B) / C;
         Double_t ct = p(2) / pt;
-        Double_t z0 = x(2) - ct * st;
+        Double_t z0;
+        Double_t dot = x(0) * p(0) + x(1) * p(1);
+        if (dot > 0.0) z0 = x(2) - ct * st;
+        else z0 = x(2) + ct * st;
         //
         Par(3) = z0;            // Store z0
@@ -443 +446 @@
         TSpline3* sp3 = new TSpline3("sp3", bg, ncl, Npt);
         if (begam > bg[0] && begam < bg[Npt - 1]) interp = sp3->Eval(begam);
-        if(begam < bg[0]) interp = bg[0];
-        if(begam > bg[Npt-1]) interp = bg[Npt-1];
         return 100 * interp;
 }

external/TrackCovariance/TrkUtil.h

@@ -25 +25 @@
         // Service routines
         //
-        void SetB(Double_t Bz){ fBz = Bz; }
+        void SetB(Double_t Bz) { fBz = Bz; }
         TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q);
         TVector3 ParToP(TVectorD Par);
@@ -54 +54 @@
                 Double_t c = 2.99792458e8;      // speed of light m/sec
                 //return TMath::C()*1.0e-9;     // Incompatible with root5
-                return c*1.0e-9;                // Reduced speed of light
+                return c*1.0e-9;                // Reduced speed of light       
         }
         //
@@ -71 +71 @@
         // Cluster counting in gas
         //
-        void SetBfield(Double_t Bz){ fBz = Bz; }
-        // Define gas volume (units = meters)
+        void SetBfield(Double_t Bz) { fBz = Bz; }
+        // Define gas volume (units = meters) 
         void SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax);
         // Gas mixture selection