Changes in / [22138b0:bae4fe8] in git

Files:

• cards/trackResolutionCMS.tcl (added)
• cards/trkCountingBTaggingCMS.tcl (added)
• modules/TrackCountingBTagging.cc (modified) (3 diffs)
• modules/TrackCountingBTagging.h (modified) (1 diff)
• modules/TrackSmearing.cc (modified) (3 diffs)
• modules/TrackSmearing.h (modified) (1 diff)

modules/TrackCountingBTagging.cc

@@ -72 +72 @@
   fNtracks = GetInt("Ntracks", 3);
 
+  fUse3D = GetBool("Use3D", false);
+
   // import input array(s)
 
@@ -95 +97 @@
   Candidate *jet, *track;
 
-  Double_t jpx, jpy;
+  Double_t jpx, jpy, jpz;
   Double_t dr, tpt;
-  Double_t xd, yd, d0, dd0, ip, sip;
+  Double_t xd, yd, zd, d0, dd0, dz, ddz, sip;
 
   Int_t sign;
@@ -110 +112 @@
     jpx = jetMomentum.Px();
     jpy = jetMomentum.Py();
+    jpz = jetMomentum.Pz();
 
     // loop over all input tracks
     fItTrackInputArray->Reset();
     count = 0;
-    while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
+    // stop once we have enough tracks
+    while((track = static_cast<Candidate*>(fItTrackInputArray->Next())) and count < fNtracks)
     {
       const TLorentzVector &trkMomentum = track->Momentum;
+      tpt = trkMomentum.Pt();
+      if(tpt < fPtMin) continue;
       
+      d0 = TMath::Abs(track->D0);
+      if(d0 > fIPmax) continue;
+
       dr = jetMomentum.DeltaR(trkMomentum);
-      tpt = trkMomentum.Pt();
+      if(dr > fDeltaR) continue;
+
       xd = track->Xd;
       yd = track->Yd;
-      d0 = TMath::Hypot(xd, yd);
-      dd0 = track->ErrorD0;
+      zd = track->Zd;
+      dd0 = TMath::Abs(track->ErrorD0);
+      dz = TMath::Abs(track->DZ);
+      ddz = TMath::Abs(track->ErrorDZ);
 
-      if(tpt < fPtMin) continue;
-      if(dr > fDeltaR) continue;
-      if(d0 > fIPmax) continue;
-
-      sign = (jpx*xd + jpy*yd > 0.0) ? 1 : -1;
-
-      ip = sign*d0;
-      sip = ip / TMath::Abs(dd0);
+      if(fUse3D){
+        sign = (jpx*xd + jpy*yd + jpz*zd > 0.0) ? 1 : -1;
+        //add transverse and longitudinal significances in quadrature
+        sip = sign * TMath::Sqrt( TMath::Power(d0 / dd0, 2) + TMath::Power(dz / ddz, 2) );
+      }
+      else {
+        sign = (jpx*xd + jpy*yd > 0.0) ? 1 : -1;
+        sip = sign * d0 / TMath::Abs(dd0);
+      }
 
       if(sip > fSigMin) count++;

modules/TrackCountingBTagging.h

@@ -54 +54 @@
   Double_t fSigMin;
   Int_t    fNtracks;
+  Bool_t   fUse3D;
 
   TIterator *fItTrackInputArray; //!

modules/TrackSmearing.cc

@@ -61 +61 @@
 void TrackSmearing::Init()
 {
+  fBz = GetDouble("Bz", 0.0);
+
   // read resolution formula
 
@@ -158 +160 @@
   Double_t pt, eta, d0, d0Error, trueD0, dz, dzError, trueDZ, p, pError, trueP, ctgTheta, ctgThetaError, trueCtgTheta, phi, phiError, truePhi;
   Double_t x, y, z, t, px, py, pz, theta;
+  Double_t q, r;
+  Double_t x_c, y_c, r_c, phi_0;
+  Double_t rcu, rc2, xd, yd, zd;
+  const Double_t c_light = 2.99792458E8;
   TProfile2D *d0ErrorHist = NULL,
              *dzErrorHist = NULL,
@@ -349 +355 @@
     y = candidate->InitialPosition.Y ();
     candidate->InitialPosition.SetZ (z + ((pz * (px * (x - beamSpotPosition.X ()) + py * (y - beamSpotPosition.Y ())) + pt * pt * (dz - z)) / (pt * pt)));
-  
+    z = candidate->InitialPosition.Z ();
+
     candidate->InitialPosition.SetT(t);
-    
+
+    // update closest approach
+    x *= 1.0E-3;
+    y *= 1.0E-3;
+    z *= 1.0E-3;
+
+    q = candidate->Charge;
+
+    r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
+    phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
+
+    // 2. helix axis coordinates
+    x_c = x + r*TMath::Sin(phi_0);
+    y_c = y - r*TMath::Cos(phi_0);
+    r_c = TMath::Hypot(x_c, y_c);
+
+    rcu = TMath::Abs(r);
+    rc2 = r_c*r_c;
+
+    // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
+    xd = x_c*x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;
+    xd = (rc2 > 0.0) ? xd / rc2 : -999;
+    yd = y_c*(-rcu*r_c + rc2);
+    yd = (rc2 > 0.0) ? yd / rc2 : -999;
+    zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
+
+    candidate->Xd = xd*1.0E3;
+    candidate->Yd = yd*1.0E3;
+    candidate->Zd = zd*1.0E3;
+
     if (fApplyToPileUp || !candidate->IsPU)
     {

modules/TrackSmearing.h

@@ -32 +32 @@
 
   Double_t ptError (const Double_t, const Double_t, const Double_t, const Double_t);
+
+  Double_t fBz;
 
   DelphesFormula *fD0Formula; //!