Changeset 6cdc544 in git for modules/PileUpJetID.cc
- Timestamp:
- Dec 21, 2014, 12:07:11 PM (10 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 38bf1ae
- Parents:
- 1d1f6a4
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/PileUpJetID.cc
r1d1f6a4 r6cdc544 16 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. 17 17 */ 18 18 19 19 20 /** \class PileUpJetID … … 53 54 54 55 PileUpJetID::PileUpJetID() : 55 fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0) 56 fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0) 56 57 { 57 58 … … 85 86 fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers")); 86 87 fItNeutralInputArray = fNeutralInputArray->MakeIterator(); 87 88 88 89 fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices")); 89 90 fItVertexInputArray = fVertexInputArray->MakeIterator(); 90 91 91 92 fZVertexResolution = GetDouble("ZVertexResolution", 0.005)*1.0E3; 92 // create output array(s) 93 94 // create output array(s) 93 95 94 96 fOutputArray = ExportArray(GetString("OutputArray", "jets")); 95 96 97 } 97 98 … … 100 101 void PileUpJetID::Finish() 101 102 { 102 103 103 if(fItJetInputArray) delete fItJetInputArray; 104 104 if(fItTrackInputArray) delete fItTrackInputArray; 105 105 if(fItNeutralInputArray) delete fItNeutralInputArray; 106 106 if(fItVertexInputArray) delete fItVertexInputArray; 107 108 107 } 109 108 … … 114 113 Candidate *candidate, *constituent; 115 114 TLorentzVector momentum, area; 116 Double_t zvtx=0; 117 118 Candidate *trk; 119 120 // find z position of primary vertex 121 115 Int_t i, nc, nn; 116 Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq; 117 Double_t dr, pt, pt_ann[5]; 118 Double_t zvtx = 0.0; 119 120 Candidate *track; 121 122 // find z position of primary vertex 123 122 124 fItVertexInputArray->Reset(); 123 125 while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next()))) … … 125 127 if(!candidate->IsPU) 126 128 { 127 zvtx = candidate->Position.Z();128 break;129 zvtx = candidate->Position.Z(); 130 break; 129 131 } 130 132 } … … 137 139 area = candidate->Area; 138 140 139 float sumpt = 0.; 140 float sumptch = 0.; 141 float sumptchpv = 0.; 142 float sumptchpu = 0.; 143 float sumdrsqptsq = 0.; 144 float sumptsq = 0.; 145 int nc = 0; 146 int nn = 0; 147 float pt_ann[5]; 148 149 for (int i = 0 ; i < 5 ; i++) { 150 pt_ann[i] = 0.; 151 } 152 153 if (fUseConstituents) { 141 sumpt = 0.0; 142 sumptch = 0.0; 143 sumptchpv = 0.0; 144 sumptchpu = 0.0; 145 sumdrsqptsq = 0.0; 146 sumptsq = 0.0; 147 nc = 0; 148 nn = 0; 149 150 for(i = 0; i < 5; ++i) 151 { 152 pt_ann[i] = 0.0; 153 } 154 155 if(fUseConstituents) 156 { 154 157 TIter itConstituents(candidate->GetCandidates()); 155 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) { 156 float pt = constituent->Momentum.Pt(); 157 float dr = candidate->Momentum.DeltaR(constituent->Momentum); 158 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) 159 { 160 pt = constituent->Momentum.Pt(); 161 dr = candidate->Momentum.DeltaR(constituent->Momentum); 158 162 sumpt += pt; 159 163 sumdrsqptsq += dr*dr*pt*pt; 160 164 sumptsq += pt*pt; 161 if (constituent->Charge == 0) { 162 // neutrals 163 nn++; 164 } else { 165 // charged 166 if (constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution) { 167 sumptchpu += pt; 168 } else { 169 sumptchpv += pt; 170 } 171 sumptch += pt; 172 nc++; 173 } 174 for (int i = 0 ; i < 5 ; i++) { 175 if (dr > 0.1*i && dr < 0.1*(i+1)) { 176 pt_ann[i] += pt; 177 } 178 } 179 } 180 } else { 165 if(constituent->Charge == 0) 166 { 167 // neutrals 168 ++nn; 169 } 170 else 171 { 172 // charged 173 if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution) 174 { 175 sumptchpu += pt; 176 } 177 else 178 { 179 sumptchpv += pt; 180 } 181 sumptch += pt; 182 ++nc; 183 } 184 for(i = 0; i < 5; ++i) 185 { 186 if(dr > 0.1*i && dr < 0.1*(i + 1)) 187 { 188 pt_ann[i] += pt; 189 } 190 } 191 } 192 } 193 else 194 { 181 195 // Not using constituents, using dr 182 196 fItTrackInputArray->Reset(); 183 while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) { 184 if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) { 185 float pt = trk->Momentum.Pt(); 186 sumpt += pt; 187 sumptch += pt; 188 if (trk->IsPU && TMath::Abs(trk->Position.Z()-zvtx) > fZVertexResolution) { 189 sumptchpu += pt; 190 } else { 191 sumptchpv += pt; 192 } 193 float dr = candidate->Momentum.DeltaR(trk->Momentum); 194 sumdrsqptsq += dr*dr*pt*pt; 195 sumptsq += pt*pt; 196 nc++; 197 for (int i = 0 ; i < 5 ; i++) { 198 if (dr > 0.1*i && dr < 0.1*(i+1)) { 197 while((track = static_cast<Candidate*>(fItTrackInputArray->Next()))) 198 { 199 if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR) 200 { 201 pt = track->Momentum.Pt(); 202 sumpt += pt; 203 sumptch += pt; 204 if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution) 205 { 206 sumptchpu += pt; 207 } 208 else 209 { 210 sumptchpv += pt; 211 } 212 dr = candidate->Momentum.DeltaR(track->Momentum); 213 sumdrsqptsq += dr*dr*pt*pt; 214 sumptsq += pt*pt; 215 nc++; 216 for(i = 0; i < 5; ++i) 217 { 218 if(dr > 0.1*i && dr < 0.1*(i + 1)) 219 { 199 220 pt_ann[i] += pt; 200 } 201 } 202 } 203 } 221 } 222 } 223 } 224 } 225 204 226 fItNeutralInputArray->Reset(); 205 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) { 206 if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) { 207 float pt = constituent->Momentum.Pt(); 208 sumpt += pt; 209 float dr = candidate->Momentum.DeltaR(constituent->Momentum); 210 sumdrsqptsq += dr*dr*pt*pt; 211 sumptsq += pt*pt; 212 nn++; 213 for (int i = 0 ; i < 5 ; i++) { 214 if (dr > 0.1*i && dr < 0.1*(i+1)) { 215 pt_ann[i] += pt; 216 } 217 } 218 } 219 } 220 } 221 222 if (sumptch > 0.) { 227 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) 228 { 229 if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) 230 { 231 pt = constituent->Momentum.Pt(); 232 sumpt += pt; 233 dr = candidate->Momentum.DeltaR(constituent->Momentum); 234 sumdrsqptsq += dr*dr*pt*pt; 235 sumptsq += pt*pt; 236 nn++; 237 for(i = 0; i < 5; ++i) 238 { 239 if(dr > 0.1*i && dr < 0.1*(i + 1)) 240 { 241 pt_ann[i] += pt; 242 } 243 } 244 } 245 } 246 } 247 248 if(sumptch > 0.0) 249 { 223 250 candidate->Beta = sumptchpu/sumptch; 224 251 candidate->BetaStar = sumptchpv/sumptch; 225 } else { 226 candidate->Beta = -999.; 227 candidate->BetaStar = -999.; 228 } 229 if (sumptsq > 0.) { 252 } 253 else 254 { 255 candidate->Beta = -999.0; 256 candidate->BetaStar = -999.0; 257 } 258 if(sumptsq > 0.0) 259 { 230 260 candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq; 231 } else { 232 candidate->MeanSqDeltaR = -999.; 261 } 262 else 263 { 264 candidate->MeanSqDeltaR = -999.0; 233 265 } 234 266 candidate->NCharged = nc; 235 267 candidate->NNeutrals = nn; 236 if (sumpt > 0.) { 268 if(sumpt > 0.0) 269 { 237 270 candidate->PTD = TMath::Sqrt(sumptsq) / sumpt; 238 for (int i = 0 ; i < 5 ; i++) { 271 for(i = 0; i < 5; ++i) 272 { 239 273 candidate->FracPt[i] = pt_ann[i]/sumpt; 240 274 } 241 } else { 242 candidate->PTD = -999.; 243 for (int i = 0 ; i < 5 ; i++) { 244 candidate->FracPt[i] = -999.; 275 } 276 else 277 { 278 candidate->PTD = -999.0; 279 for(i = 0; i < 5; ++i) 280 { 281 candidate->FracPt[i] = -999.0; 245 282 } 246 283 }
Note:
See TracChangeset
for help on using the changeset viewer.