Changeset 9666d3b in git for modules/StatusPidFilter.cc

Timestamp:

Apr 11, 2018, 4:19:18 PM (7 years ago)

Author:

Ulrike Schnoor <ulrike.schnoor@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

6ba7847

Parents:

6f96f62 (diff), a670493 (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.

Message:

merged the upstream/master

File:

• modules/StatusPidFilter.cc (modified) (2 diffs)

modules/StatusPidFilter.cc

@@ -51 +51 @@
 using namespace std;
 
+namespace {
+  // integer power (faster than TMath::Pow() + cast to integer)
+  int ipow(int base, int exp)
+  {
+      int result = 1;
+      while (exp)
+      {
+          if (exp & 1)
+              result *= base;
+          exp >>= 1;
+          base *= base;
+      }
+
+      return result;
+  }
+
+  // standalone function to extract the i-th digit from a number (counting from 0 = rightmost, etc..)
+  int digit(int val, int i)
+  {
+      int y    = ipow(10,i);
+      int z    = val/y;
+      int val2 = val / (y * 10);
+      return (z - val2*10);
+  }
+
+  //  return the first two digits if this is a "fundamental" particle
+  //  ID = 100 is a special case (internal generator ID's are 81-100)
+  //  also, 101 and 102 are now used (by HepPID) for geantinos
+  int fundamentalID(int pdgCode)
+  {
+      pdgCode = abs(pdgCode);
+      if( ( digit(pdgCode, 9) == 1 ) && ( digit(pdgCode, 8) == 0 ) ) { return 0; }
+      if( digit(pdgCode, 2) == 0 && digit(pdgCode, 3) == 0) {
+          return pdgCode%10000;
+      } else if( pdgCode <= 102 ) {
+          return pdgCode;
+      } else {
+          return 0;
+      }
+  }
+
+  bool hasBottom(int pdgCode)
+  {
+      if ( (pdgCode/10000000) > 0)
+        return false;
+      if (pdgCode <= 100)
+        return false;
+      if (fundamentalID(pdgCode) <= 100 && fundamentalID(pdgCode) > 0)
+        return false;
+      if( digit(pdgCode, 3) == 5 || digit(pdgCode, 2) == 5 || digit(pdgCode, 1) == 5 )
+        return true;      
+      return false;
+  }
+
+  bool isTauDaughter(int pdgCode, int M1, const TObjArray *fInputArray){
+    //not needed, just to speed up the code - can be further refined but gives only negligible improvement:
+    if ( pdgCode==15 || pdgCode<11 || (pdgCode > 22 && pdgCode < 100) || pdgCode>1000 )
+      return false;
+
+    if ( M1 < 0 ) 
+      return false;
+
+    Candidate *mother;
+    mother = static_cast<Candidate*>(fInputArray->At( M1 ));
+    if ( TMath::Abs(mother->PID) == 15 )
+      return true;
+
+    return false;
+  }
+
+}
+
+
 //------------------------------------------------------------------------------
 
@@ -117 +190 @@
     if(pdgCode > 22 && pdgCode < 43) pass = kTRUE;
 
-    if(!pass || candidate->Momentum.Pt() < fPTMin) continue;
+    //Stable photons 
+    if(pdgCode == 22 && status==1) pass = kTRUE;
+
+    // logic ported from HepPDF: http://lcgapp.cern.ch/project/simu/HepPDT/HepPDT.2.05.02/html/ParticleID_8cc-source.html#l00081
+    bool is_b_hadron = hasBottom(pdgCode);
+    bool is_b_quark  = (pdgCode == 5);
+
+    bool is_tau_daughter = isTauDaughter(pdgCode, candidate->M1, fInputArray);
+
+    if (is_b_hadron)
+      pass = kTRUE;
+
+    if (is_tau_daughter)
+      pass = kTRUE;
+
+    // fPTMin not applied to b_hadrons / b_quarks to allow for b-enriched sample stitching
+    // fPTMin not applied to tau decay products to allow visible-tau four momentum determination
+    if(!pass || (candidate->Momentum.Pt() < fPTMin && !(is_b_hadron || is_b_quark || is_tau_daughter)) ) continue;
 
     fOutputArray->Add(candidate);