Changes in modules/IdentificationMap.cc [95aa610:cab38f6] in git

File:

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

modules/IdentificationMap.cc

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