- Timestamp:
- Dec 21, 2014, 7:48:30 PM (10 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- d4b9697
- Parents:
- e5767b57
- Location:
- modules
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/IdentificationMap.cc
re5767b57 ra5bbe8a 69 69 void IdentificationMap::Init() 70 70 { 71 // read efficiency formula72 73 74 71 TMisIDMap::iterator itEfficiencyMap; 75 72 ExRootConfParam param; … … 87 84 formula->Compile(param[i*3 + 2].GetString()); 88 85 pdg = param[i*3].GetInt(); 89 fEfficiencyMap.insert(make_pair(pdg,make_pair(param[i*3 + 1].GetInt(),formula))); 90 91 // cout<<param[i*3].GetInt()<<","<<param[i*3+1].GetInt()<<","<<param[i*3 + 2].GetString()<<endl; 92 86 fEfficiencyMap.insert(make_pair(pdg, make_pair(param[i*3 + 1].GetInt(), formula))); 93 87 } 94 88 … … 100 94 formula->Compile("1.0"); 101 95 102 fEfficiencyMap.insert(make_pair(0, make_pair(0,formula)));96 fEfficiencyMap.insert(make_pair(0, make_pair(0, formula))); 103 97 } 104 98 … … 126 120 if(formula) delete formula; 127 121 } 128 129 122 } 130 123 … … 138 131 pair <TMisIDMap::iterator, TMisIDMap::iterator> range; 139 132 DelphesFormula *formula; 140 Int_t pdg In, pdgOut, charge;133 Int_t pdgCodeIn, pdgCodeOut, charge; 141 134 142 Double_t P, Pi; 143 144 // cout<<"------------ New Event ------------"<<endl; 135 Double_t p, r, total; 145 136 146 137 fItInputArray->Reset(); … … 152 143 phi = candidatePosition.Phi(); 153 144 pt = candidateMomentum.Pt(); 154 pdg In = candidate->PID;145 pdgCodeIn = candidate->PID; 155 146 charge = candidate->Charge; 156 147 157 // cout<<"------------ New Candidate ------------"<<endl;158 // cout<<candidate->PID<<" "<<pt<<","<<eta<<","<<phi<<endl;148 // first check that PID of this particle is specified in the map 149 // otherwise, look for PID = 0 159 150 160 P = 1.0;151 itEfficiencyMap = fEfficiencyMap.find(pdgCodeIn); 161 152 162 //first check that PID of this particle is specified in cfg, if not set look for PID=0 163 164 itEfficiencyMap = fEfficiencyMap.find(pdgIn); 165 166 range = fEfficiencyMap.equal_range(pdgIn); 167 if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgIn); 153 range = fEfficiencyMap.equal_range(pdgCodeIn); 154 if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgCodeIn); 168 155 if(range.first == range.second) range = fEfficiencyMap.equal_range(0); 169 156 170 //loop over submap for this pid 171 for (TMisIDMap::iterator it=range.first; it!=range.second; ++it) 157 r = gRandom->Uniform(); 158 total = 0.0; 159 160 // loop over sub-map for this PID 161 for(TMisIDMap::iterator it = range.first; it != range.second; ++it) 172 162 { 163 formula = (it->second).second; 164 pdgCodeOut = (it->second).first; 173 165 174 formula = (it->second).second; 175 pdgOut = (it->second).first; 166 p = formula->Eval(pt, eta); 176 167 177 Pi = formula->Eval(pt, eta); 178 179 // check that sum of probabilities does not exceed 1. 180 P = (P - Pi)/P; 181 182 if( P < 0.0 ) continue; 183 else 168 if(total <= r && r < total + p) 184 169 { 185 186 //randomly assign a PID to particle according to map 187 Double_t rndm = gRandom->Uniform(); 188 189 if(rndm > P) 190 { 191 //change PID of particle 192 if(pdgOut != 0) candidate->PID = charge*pdgOut; 193 fOutputArray->Add(candidate); 194 break; 195 } 170 // change PID of particle 171 if(pdgCodeOut != 0) candidate->PID = charge*pdgCodeOut; 172 fOutputArray->Add(candidate); 173 break; 196 174 } 197 175 176 total += p; 198 177 } 199 200 } 201 178 } 202 179 } 203 180 -
modules/IdentificationMap.h
re5767b57 ra5bbe8a 49 49 private: 50 50 51 typedef std::multimap< Int_t, std::pair< Int_t, DelphesFormula * > > TMisIDMap; //!51 typedef std::multimap< Int_t, std::pair< Int_t, DelphesFormula * > > TMisIDMap; //! 52 52 53 TMisIDMap fEfficiencyMap; 53 TMisIDMap fEfficiencyMap; //! 54 54 55 55 TIterator *fItInputArray; //!
Note:
See TracChangeset
for help on using the changeset viewer.