| 1 | /*
|
|---|
| 2 | * Delphes: a framework for fast simulation of a generic collider experiment
|
|---|
| 3 | * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
|
|---|
| 4 | *
|
|---|
| 5 | * This program is free software: you can redistribute it and/or modify
|
|---|
| 6 | * it under the terms of the GNU General Public License as published by
|
|---|
| 7 | * the Free Software Foundation, either version 3 of the License, or
|
|---|
| 8 | * (at your option) any later version.
|
|---|
| 9 | *
|
|---|
| 10 | * This program is distributed in the hope that it will be useful,
|
|---|
| 11 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 12 | * MERCHANTABILITY || FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 13 | * GNU General Public License for more details.
|
|---|
| 14 | *
|
|---|
| 15 | * You should have received a copy of the GNU General Public License
|
|---|
| 16 | * along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 17 | */
|
|---|
| 18 |
|
|---|
| 19 |
|
|---|
| 20 | /** \class JetFlavorAssociation
|
|---|
| 21 | *
|
|---|
| 22 | * Find origin of jet && evaluate jet flavor
|
|---|
| 23 | *
|
|---|
| 24 | * \author P. Demin - UCL, Louvain-la-Neuve
|
|---|
| 25 | *
|
|---|
| 26 | */
|
|---|
| 27 |
|
|---|
| 28 | #include "modules/JetFlavorAssociation.h"
|
|---|
| 29 |
|
|---|
| 30 | #include "classes/DelphesClasses.h"
|
|---|
| 31 | #include "classes/DelphesFactory.h"
|
|---|
| 32 | #include "classes/DelphesFormula.h"
|
|---|
| 33 |
|
|---|
| 34 | #include "ExRootAnalysis/ExRootResult.h"
|
|---|
| 35 | #include "ExRootAnalysis/ExRootFilter.h"
|
|---|
| 36 | #include "ExRootAnalysis/ExRootClassifier.h"
|
|---|
| 37 |
|
|---|
| 38 | #include "TMath.h"
|
|---|
| 39 | #include "TString.h"
|
|---|
| 40 | #include "TFormula.h"
|
|---|
| 41 | #include "TRandom3.h"
|
|---|
| 42 | #include "TObjArray.h"
|
|---|
| 43 | #include "TDatabasePDG.h"
|
|---|
| 44 | #include "TLorentzVector.h"
|
|---|
| 45 |
|
|---|
| 46 | #include <algorithm>
|
|---|
| 47 | #include <stdexcept>
|
|---|
| 48 | #include <iostream>
|
|---|
| 49 | #include <sstream>
|
|---|
| 50 |
|
|---|
| 51 | using namespace std;
|
|---|
| 52 |
|
|---|
| 53 | //------------------------------------------------------------------------------
|
|---|
| 54 |
|
|---|
| 55 | class PartonClassifier: public ExRootClassifier
|
|---|
| 56 | {
|
|---|
| 57 | public:
|
|---|
| 58 |
|
|---|
| 59 | PartonClassifier() {}
|
|---|
| 60 | Int_t GetCategory(TObject *object);
|
|---|
| 61 | Double_t fEtaMax, fPTMin;
|
|---|
| 62 | };
|
|---|
| 63 |
|
|---|
| 64 | //------------------------------------------------------------------------------
|
|---|
| 65 | // https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/PartonSelector.cc
|
|---|
| 66 |
|
|---|
| 67 | Int_t PartonClassifier::GetCategory(TObject *object)
|
|---|
| 68 | {
|
|---|
| 69 | // select parton in the parton list
|
|---|
| 70 |
|
|---|
| 71 | Candidate *parton = static_cast<Candidate *>(object);
|
|---|
| 72 | const TLorentzVector &momentum = parton->Momentum;
|
|---|
| 73 | Int_t pdgCode;
|
|---|
| 74 |
|
|---|
| 75 | // inside the eta && momentum range (be a little bit larger that the tracking coverage
|
|---|
| 76 | if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
|
|---|
| 77 |
|
|---|
| 78 | pdgCode = TMath::Abs(parton->PID);
|
|---|
| 79 |
|
|---|
| 80 | if(parton->Status == -1) return -1;
|
|---|
| 81 | if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
|
|---|
| 82 | if(parton->Status == 3 || parton->Status == 2) return 0; // if status 3 return
|
|---|
| 83 |
|
|---|
| 84 | return 0;
|
|---|
| 85 | }
|
|---|
| 86 |
|
|---|
| 87 | //------------------------------------------------------------------------------
|
|---|
| 88 |
|
|---|
| 89 | class ParticleLHEFClassifier: public ExRootClassifier
|
|---|
| 90 | {
|
|---|
| 91 | public:
|
|---|
| 92 |
|
|---|
| 93 | ParticleLHEFClassifier() {}
|
|---|
| 94 | Int_t GetCategory(TObject *object);
|
|---|
| 95 | Double_t fEtaMax, fPTMin;
|
|---|
| 96 | };
|
|---|
| 97 |
|
|---|
| 98 | Int_t ParticleLHEFClassifier::GetCategory(TObject *object)
|
|---|
| 99 | {
|
|---|
| 100 | // select parton in the parton list
|
|---|
| 101 |
|
|---|
| 102 | Candidate *particleLHEF = static_cast<Candidate *>(object);
|
|---|
| 103 | const TLorentzVector &momentum = particleLHEF->Momentum;
|
|---|
| 104 | Int_t pdgCode;
|
|---|
| 105 |
|
|---|
| 106 | // inside the eta && momentum range (be a little bit larger that the tracking coverage
|
|---|
| 107 | if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
|
|---|
| 108 |
|
|---|
| 109 | pdgCode = TMath::Abs(particleLHEF->PID);
|
|---|
| 110 | if(particleLHEF->Status == -1) return -1;
|
|---|
| 111 | if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
|
|---|
| 112 | if(particleLHEF->Status != 1) return -1; // if status 3 return
|
|---|
| 113 |
|
|---|
| 114 | return 0;
|
|---|
| 115 | }
|
|---|
| 116 |
|
|---|
| 117 | //------------------------------------------------------------------------------
|
|---|
| 118 |
|
|---|
| 119 | JetFlavorAssociation::JetFlavorAssociation() :
|
|---|
| 120 | fPartonClassifier(0), fPartonFilter(0), fParticleLHEFFilter(0),
|
|---|
| 121 | fItPartonInputArray(0), fItParticleInputArray(0),
|
|---|
| 122 | fItParticleLHEFInputArray(0), fItJetInputArray(0)
|
|---|
| 123 | {
|
|---|
| 124 | fPartonClassifier = new PartonClassifier;
|
|---|
| 125 | fParticleLHEFClassifier = new ParticleLHEFClassifier;
|
|---|
| 126 | }
|
|---|
| 127 |
|
|---|
| 128 | //------------------------------------------------------------------------------
|
|---|
| 129 |
|
|---|
| 130 | JetFlavorAssociation::~JetFlavorAssociation()
|
|---|
| 131 | {
|
|---|
| 132 | if(fPartonClassifier) delete fPartonClassifier;
|
|---|
| 133 | if(fParticleLHEFClassifier) delete fParticleLHEFClassifier;
|
|---|
| 134 | }
|
|---|
| 135 |
|
|---|
| 136 | //------------------------------------------------------------------------------
|
|---|
| 137 |
|
|---|
| 138 | void JetFlavorAssociation::Init()
|
|---|
| 139 | {
|
|---|
| 140 | ExRootConfParam param;
|
|---|
| 141 |
|
|---|
| 142 | fDeltaR = GetDouble("DeltaR", 0.5);
|
|---|
| 143 |
|
|---|
| 144 | fPartonClassifier->fPTMin = GetDouble("PartonPTMin", 0.0);
|
|---|
| 145 | fPartonClassifier->fEtaMax = GetDouble("PartonEtaMax", 2.5);
|
|---|
| 146 |
|
|---|
| 147 | fParticleLHEFClassifier->fPTMin = GetDouble("PartonPTMin", 0.0);
|
|---|
| 148 | fParticleLHEFClassifier->fEtaMax = GetDouble("PartonEtaMax", 2.5);
|
|---|
| 149 |
|
|---|
| 150 | // import input array(s)
|
|---|
| 151 | fPartonInputArray = ImportArray(GetString("PartonInputArray", "Delphes/partons"));
|
|---|
| 152 | fItPartonInputArray = fPartonInputArray->MakeIterator();
|
|---|
| 153 |
|
|---|
| 154 | fParticleInputArray = ImportArray(GetString("ParticleInputArray", "Delphes/allParticles"));
|
|---|
| 155 | fItParticleInputArray = fParticleInputArray->MakeIterator();
|
|---|
| 156 |
|
|---|
| 157 | fParticleLHEFInputArray = ImportArray(GetString("ParticleLHEFInputArray", "Delphes/allParticlesLHEF"));
|
|---|
| 158 | fItParticleLHEFInputArray = fParticleLHEFInputArray->MakeIterator();
|
|---|
| 159 |
|
|---|
| 160 | fPartonFilter = new ExRootFilter(fPartonInputArray);
|
|---|
| 161 | fParticleLHEFFilter = new ExRootFilter(fParticleLHEFInputArray);
|
|---|
| 162 |
|
|---|
| 163 | fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
|
|---|
| 164 | fItJetInputArray = fJetInputArray->MakeIterator();
|
|---|
| 165 | }
|
|---|
| 166 |
|
|---|
| 167 | //------------------------------------------------------------------------------
|
|---|
| 168 |
|
|---|
| 169 | void JetFlavorAssociation::Finish()
|
|---|
| 170 | {
|
|---|
| 171 | if(fPartonFilter) delete fPartonFilter;
|
|---|
| 172 | if(fParticleLHEFFilter) delete fParticleLHEFFilter;
|
|---|
| 173 |
|
|---|
| 174 | if(fItJetInputArray) delete fItJetInputArray;
|
|---|
| 175 | if(fItParticleLHEFInputArray) delete fItParticleLHEFInputArray;
|
|---|
| 176 | if(fItParticleInputArray) delete fItParticleInputArray;
|
|---|
| 177 | if(fItPartonInputArray) delete fItPartonInputArray;
|
|---|
| 178 | }
|
|---|
| 179 |
|
|---|
| 180 | //------------------------------------------------------------------------------
|
|---|
| 181 |
|
|---|
| 182 | void JetFlavorAssociation::Process(){
|
|---|
| 183 |
|
|---|
| 184 | Candidate *jet;
|
|---|
| 185 | TObjArray *partonArray;
|
|---|
| 186 | TObjArray *partonLHEFArray;
|
|---|
| 187 |
|
|---|
| 188 | // select quark && gluons
|
|---|
| 189 | fPartonFilter->Reset();
|
|---|
| 190 | partonArray = fPartonFilter->GetSubArray(fPartonClassifier, 0); // get the filtered parton array
|
|---|
| 191 |
|
|---|
| 192 | if(partonArray == 0) return;
|
|---|
| 193 | TIter itPartonArray(partonArray);
|
|---|
| 194 |
|
|---|
| 195 | fParticleLHEFFilter->Reset();
|
|---|
| 196 | partonLHEFArray = fParticleLHEFFilter->GetSubArray(fParticleLHEFClassifier, 0); // get the filtered parton array
|
|---|
| 197 |
|
|---|
| 198 | if(partonLHEFArray == 0) return;
|
|---|
| 199 | TIter itPartonLHEFArray(partonLHEFArray);
|
|---|
| 200 |
|
|---|
| 201 | // loop over all input jets
|
|---|
| 202 | fItJetInputArray->Reset();
|
|---|
| 203 | while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
|
|---|
| 204 | {
|
|---|
| 205 | // get standard flavor
|
|---|
| 206 | GetAlgoFlavor(jet, itPartonArray, itPartonLHEFArray);
|
|---|
| 207 | GetPhysicsFlavor(jet, itPartonArray, itPartonLHEFArray);
|
|---|
| 208 | }
|
|---|
| 209 | }
|
|---|
| 210 |
|
|---|
| 211 | //------------------------------------------------------------------------------
|
|---|
| 212 | // Standard definition of jet flavor in
|
|---|
| 213 | // https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/JetPartonMatcher.cc?v=CMSSW_7_3_0_pre1
|
|---|
| 214 |
|
|---|
| 215 | void JetFlavorAssociation::GetAlgoFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
|
|---|
| 216 | {
|
|---|
| 217 | float maxPt = 0;
|
|---|
| 218 | bool isGoodParton = true;
|
|---|
| 219 | int daughterCounter = 0;
|
|---|
| 220 | Candidate *parton, *partonLHEF;
|
|---|
| 221 | Candidate *tempParton = 0, *tempPartonHighestPt = 0;
|
|---|
| 222 | int pdgCode, pdgCodeMax = -1;
|
|---|
| 223 |
|
|---|
| 224 | itPartonArray.Reset();
|
|---|
| 225 | while((parton = static_cast<Candidate *>(itPartonArray.Next())))
|
|---|
| 226 | {
|
|---|
| 227 | // default delphes method
|
|---|
| 228 | pdgCode = TMath::Abs(parton->PID);
|
|---|
| 229 | if(TMath::Abs(parton->PID) == 21) pdgCode = 0;
|
|---|
| 230 | if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
|
|---|
| 231 | {
|
|---|
| 232 | if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
|
|---|
| 233 | }
|
|---|
| 234 |
|
|---|
| 235 | isGoodParton = true;
|
|---|
| 236 | itPartonLHEFArray.Reset();
|
|---|
| 237 | while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
|
|---|
| 238 | {
|
|---|
| 239 | if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.001 &&
|
|---|
| 240 | parton->PID == partonLHEF->PID &&
|
|---|
| 241 | partonLHEF->Charge == parton->Charge)
|
|---|
| 242 | {
|
|---|
| 243 | isGoodParton = false;
|
|---|
| 244 | break;
|
|---|
| 245 | }
|
|---|
| 246 |
|
|---|
| 247 | if(!isGoodParton) continue;
|
|---|
| 248 |
|
|---|
| 249 | // check the daugheter
|
|---|
| 250 | daughterCounter = 0;
|
|---|
| 251 | if(parton->D1 != -1 || parton->D2 != -1)
|
|---|
| 252 | {
|
|---|
| 253 | // partons are only quarks || gluons
|
|---|
| 254 | int daughterFlavor1 = -1;
|
|---|
| 255 | int daughterFlavor2 = -1;
|
|---|
| 256 | if(parton->D1 != -1) daughterFlavor1 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D1))->PID);
|
|---|
| 257 | if(parton->D2 != -1) daughterFlavor2 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D2))->PID);
|
|---|
| 258 | if((daughterFlavor1 == 1 || daughterFlavor1 == 2 || daughterFlavor1 == 3 || daughterFlavor1 == 4 || daughterFlavor1 == 5 || daughterFlavor1 == 21)) daughterCounter++;
|
|---|
| 259 | if((daughterFlavor2 == 1 || daughterFlavor2 == 2 || daughterFlavor2 == 3 || daughterFlavor2 == 4 || daughterFlavor1 == 5 || daughterFlavor2 == 21)) daughterCounter++;
|
|---|
| 260 | }
|
|---|
| 261 | if(daughterCounter > 0) continue;
|
|---|
| 262 | if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
|
|---|
| 263 | {
|
|---|
| 264 | // if not yet found && pdgId is a c, take as c
|
|---|
| 265 | if(TMath::Abs(parton->PID) == 4) tempParton = parton;
|
|---|
| 266 | if(TMath::Abs(parton->PID) == 5) tempParton = parton;
|
|---|
| 267 | if(parton->Momentum.Pt() > maxPt)
|
|---|
| 268 | {
|
|---|
| 269 | maxPt = parton->Momentum.Pt();
|
|---|
| 270 | tempPartonHighestPt = parton;
|
|---|
| 271 | }
|
|---|
| 272 | }
|
|---|
| 273 | }
|
|---|
| 274 | }
|
|---|
| 275 |
|
|---|
| 276 | if(!tempParton) tempParton = tempPartonHighestPt;
|
|---|
| 277 | jet->FlavorAlgo = tempParton ? TMath::Abs(tempParton->PID) : 0;
|
|---|
| 278 |
|
|---|
| 279 | if(pdgCodeMax == 0) pdgCodeMax = 21;
|
|---|
| 280 | if(pdgCodeMax == -1) pdgCodeMax = 0;
|
|---|
| 281 |
|
|---|
| 282 | jet->Flavor = pdgCodeMax;
|
|---|
| 283 | }
|
|---|
| 284 |
|
|---|
| 285 | //------------------------------------------------------------------------------
|
|---|
| 286 |
|
|---|
| 287 | void JetFlavorAssociation::GetPhysicsFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
|
|---|
| 288 | {
|
|---|
| 289 | int partonCounter = 0;
|
|---|
| 290 | float biggerConeSize = 0.7;
|
|---|
| 291 | float dist;
|
|---|
| 292 | bool isGoodCandidate;
|
|---|
| 293 | int contaminatingFlavor = 0;
|
|---|
| 294 | int motherCounter = 0;
|
|---|
| 295 | Candidate *parton, *partonLHEF, *mother1, *mother2;
|
|---|
| 296 | Candidate *tempParton = 0;
|
|---|
| 297 | vector<Candidate *> contaminations;
|
|---|
| 298 | vector<Candidate *>::iterator itContaminations;
|
|---|
| 299 |
|
|---|
| 300 | contaminations.clear();
|
|---|
| 301 |
|
|---|
| 302 | itPartonLHEFArray.Reset();
|
|---|
| 303 | while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
|
|---|
| 304 | {
|
|---|
| 305 | dist = jet->Momentum.DeltaR(partonLHEF->Momentum); // take the DR
|
|---|
| 306 |
|
|---|
| 307 | if(partonLHEF->Status == 1 && dist <= fDeltaR)
|
|---|
| 308 | {
|
|---|
| 309 | tempParton = partonLHEF;
|
|---|
| 310 | partonCounter++;
|
|---|
| 311 | }
|
|---|
| 312 | }
|
|---|
| 313 |
|
|---|
| 314 | itPartonArray.Reset();
|
|---|
| 315 | itPartonLHEFArray.Reset();
|
|---|
| 316 | while((parton = static_cast<Candidate *>(itPartonArray.Next())))
|
|---|
| 317 | {
|
|---|
| 318 | dist = jet->Momentum.DeltaR(parton->Momentum); // take the DR
|
|---|
| 319 | isGoodCandidate = true;
|
|---|
| 320 | while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
|
|---|
| 321 | {
|
|---|
| 322 | if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.01 &&
|
|---|
| 323 | parton->PID == partonLHEF->PID &&
|
|---|
| 324 | partonLHEF->Charge == parton->Charge)
|
|---|
| 325 | {
|
|---|
| 326 | isGoodCandidate = false;
|
|---|
| 327 | break;
|
|---|
| 328 | }
|
|---|
| 329 | }
|
|---|
| 330 |
|
|---|
| 331 | if(!isGoodCandidate) continue;
|
|---|
| 332 |
|
|---|
| 333 | if(parton->D1 != -1 || parton->D2 != -1)
|
|---|
| 334 | {
|
|---|
| 335 | if((TMath::Abs(parton->PID) < 4 || TMath::Abs(parton->PID) == 21)) continue;
|
|---|
| 336 | if(dist < biggerConeSize) contaminations.push_back(parton);
|
|---|
| 337 | }
|
|---|
| 338 | }
|
|---|
| 339 |
|
|---|
| 340 | if(partonCounter != 1)
|
|---|
| 341 | {
|
|---|
| 342 | jet->FlavorPhys = 0;
|
|---|
| 343 | }
|
|---|
| 344 | else if(contaminations.size() == 0)
|
|---|
| 345 | {
|
|---|
| 346 | jet->FlavorPhys = TMath::Abs(tempParton->PID);
|
|---|
| 347 | }
|
|---|
| 348 | else if(contaminations.size() > 0)
|
|---|
| 349 | {
|
|---|
| 350 | jet->FlavorPhys = TMath::Abs(tempParton->PID);
|
|---|
| 351 |
|
|---|
| 352 | for(itContaminations = contaminations.begin(); itContaminations != contaminations.end(); ++itContaminations)
|
|---|
| 353 | {
|
|---|
| 354 | parton = *itContaminations;
|
|---|
| 355 | contaminatingFlavor = TMath::Abs(parton->PID);
|
|---|
| 356 | motherCounter = 0;
|
|---|
| 357 | if(parton->M1 != -1) motherCounter++;
|
|---|
| 358 | if(parton->M2 != -1) motherCounter++;
|
|---|
| 359 |
|
|---|
| 360 | if(parton->M1 != -1)
|
|---|
| 361 | {
|
|---|
| 362 | mother1 = static_cast<Candidate *>(fParticleInputArray->At(parton->M1));
|
|---|
| 363 | if(mother1 && motherCounter > 0 && mother1->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
|
|---|
| 364 | }
|
|---|
| 365 | if(parton->M2 != -1)
|
|---|
| 366 | {
|
|---|
| 367 | mother2 = static_cast<Candidate *>(fParticleInputArray->At(parton->M2));
|
|---|
| 368 | if(mother2 && motherCounter > 0 && mother2->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
|
|---|
| 369 | }
|
|---|
| 370 | // mother is the initialParton --> OK
|
|---|
| 371 | if(TMath::Abs(tempParton->PID) == 4)
|
|---|
| 372 | {
|
|---|
| 373 | // keep association --> the initialParton is a c --> the contaminated parton is a c
|
|---|
| 374 | if(contaminatingFlavor == 4) continue;
|
|---|
| 375 | jet->FlavorPhys = 0; // all the other cases reject!
|
|---|
| 376 | break;
|
|---|
| 377 | }
|
|---|
| 378 | }
|
|---|
| 379 | }
|
|---|
| 380 | }
|
|---|
