Changeset a5bbe8a in git for modules/IdentificationMap.cc

Timestamp:

Dec 21, 2014, 7:48:30 PM (10 years ago)

Author:

Pavel Demin <pavel.demin@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

d4b9697

Parents:

e5767b57

Message:

fix probability calculation in IdentificationMap

File:

• modules/IdentificationMap.cc (modified) (6 diffs)

modules/IdentificationMap.cc

@@ -69 +69 @@
 void IdentificationMap::Init()
 {
-  // read efficiency formula
-
-
   TMisIDMap::iterator itEfficiencyMap;
   ExRootConfParam param;
@@ -87 +84 @@
     formula->Compile(param[i*3 + 2].GetString());
     pdg = param[i*3].GetInt();
-    fEfficiencyMap.insert(make_pair(pdg,make_pair(param[i*3 + 1].GetInt(),formula)));
-
-   // cout<<param[i*3].GetInt()<<","<<param[i*3+1].GetInt()<<","<<param[i*3 + 2].GetString()<<endl;
-
+    fEfficiencyMap.insert(make_pair(pdg, make_pair(param[i*3 + 1].GetInt(), formula)));
   }
 
@@ -100 +94 @@
     formula->Compile("1.0");
 
-    fEfficiencyMap.insert(make_pair(0,make_pair(0,formula)));
+    fEfficiencyMap.insert(make_pair(0, make_pair(0, formula)));
   }
 
@@ -126 +120 @@
     if(formula) delete formula;
   }
-
 }
 
@@ -138 +131 @@
   pair <TMisIDMap::iterator, TMisIDMap::iterator> range;
   DelphesFormula *formula;
-  Int_t pdgIn, pdgOut, charge;
+  Int_t pdgCodeIn, pdgCodeOut, charge;
 
-  Double_t P, Pi;
-
- // cout<<"------------ New Event ------------"<<endl;
+  Double_t p, r, total;
 
   fItInputArray->Reset();
@@ -152 +143 @@
     phi = candidatePosition.Phi();
     pt = candidateMomentum.Pt();
-    pdgIn = candidate->PID;
+    pdgCodeIn = candidate->PID;
     charge = candidate->Charge;
 
-   // cout<<"------------ New Candidate ------------"<<endl;
-   // cout<<candidate->PID<<"   "<<pt<<","<<eta<<","<<phi<<endl;
+    // first check that PID of this particle is specified in the map
+    // otherwise, look for PID = 0
 
-    P = 1.0;
+    itEfficiencyMap = fEfficiencyMap.find(pdgCodeIn);
 
-    //first check that PID of this particle is specified in cfg, if not set look for PID=0
-
-    itEfficiencyMap = fEfficiencyMap.find(pdgIn);
-
-    range = fEfficiencyMap.equal_range(pdgIn);
-    if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgIn);
+    range = fEfficiencyMap.equal_range(pdgCodeIn);
+    if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgCodeIn);
     if(range.first == range.second) range = fEfficiencyMap.equal_range(0);
 
-    //loop over submap for this pid
-    for (TMisIDMap::iterator it=range.first; it!=range.second; ++it)
+    r = gRandom->Uniform();
+    total = 0.0;
+
+    // loop over sub-map for this PID
+    for(TMisIDMap::iterator it = range.first; it != range.second; ++it)
     {
+      formula = (it->second).second;
+      pdgCodeOut = (it->second).first;
 
-      formula = (it->second).second;
-      pdgOut = (it->second).first;
+      p = formula->Eval(pt, eta);
 
-      Pi = formula->Eval(pt, eta);
-
-      // check that sum of probabilities does not exceed 1.
-      P = (P - Pi)/P;
-
-      if( P < 0.0 ) continue;
-      else
+      if(total <= r && r < total + p)
       {
-
-       //randomly assign a PID to particle according to map
-       Double_t rndm = gRandom->Uniform();
-
-       if(rndm > P)
-       {
-         //change PID of particle
-         if(pdgOut != 0) candidate->PID = charge*pdgOut;
-         fOutputArray->Add(candidate);
-         break;
-       }
+        // change PID of particle
+        if(pdgCodeOut != 0) candidate->PID = charge*pdgCodeOut;
+        fOutputArray->Add(candidate);
+        break;
       }
 
+      total += p;
     }
-
-   }
-
+  }
 }