Fork me on GitHub

source: git/modules/PhotonID.cc@ 09ffaee

Last change on this file since 09ffaee was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

  • Property mode set to 100644
File size: 6.4 KB
Line 
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 or 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/** \class PhotonID
20 *
21 * Applies complex photon Id. Reconstructed photon candidtes are first separated into matched and non-matched to gen particles.
22 * Non-matched pass the "fake" efficiency. Matched photons get further splitted into isolated and non-isolated (user can choose criterion for isolation)
23 * Isolated photons pass the "prompt" efficiency while the non-isolated pass the "non-prompt" efficiency
24 *
25 * \author M. Selvaggi CERN
26 *
27 */
28
29#include "modules/PhotonID.h"
30
31#include "classes/DelphesClasses.h"
32#include "classes/DelphesFactory.h"
33#include "classes/DelphesFormula.h"
34
35#include "ExRootAnalysis/ExRootClassifier.h"
36#include "ExRootAnalysis/ExRootFilter.h"
37#include "ExRootAnalysis/ExRootResult.h"
38
39#include "TDatabasePDG.h"
40#include "TFormula.h"
41#include "TLorentzVector.h"
42#include "TMath.h"
43#include "TObjArray.h"
44#include "TRandom3.h"
45#include "TString.h"
46
47#include <algorithm>
48#include <iostream>
49#include <sstream>
50#include <stdexcept>
51
52using namespace std;
53
54//------------------------------------------------------------------------------
55
56PhotonID::PhotonID() :
57 fPromptFormula(0), fNonPromptFormula(0), fFakeFormula(0), fItInputPhotonArray(0), fItInputGenArray(0)
58{
59 fPromptFormula = new DelphesFormula;
60 fNonPromptFormula = new DelphesFormula;
61 fFakeFormula = new DelphesFormula;
62}
63
64//------------------------------------------------------------------------------
65
66PhotonID::~PhotonID()
67{
68 if(fPromptFormula) delete fPromptFormula;
69 if(fNonPromptFormula) delete fNonPromptFormula;
70 if(fFakeFormula) delete fFakeFormula;
71}
72
73//------------------------------------------------------------------------------
74
75void PhotonID::Init()
76{
77
78 // read PhotonID formulae
79 fPromptFormula->Compile(GetString("PromptFormula", "1.0"));
80 fNonPromptFormula->Compile(GetString("NonPromptFormula", "1.0"));
81 fFakeFormula->Compile(GetString("FakeFormula", "1.0"));
82
83 // import input arrays
84 fInputPhotonArray = ImportArray(GetString("InputPhotonArray", "PhotonIsolation/photons"));
85 fItInputPhotonArray = fInputPhotonArray->MakeIterator();
86
87 // use filtered collection for speed
88 fInputGenArray = ImportArray(GetString("InputGenArray", "GenParticleFilter/filteredParticles"));
89 fItInputGenArray = fInputGenArray->MakeIterator();
90
91 // min pt to be considered, make sure this threshold is higher than threshold in particle filter
92 fPTMin = GetDouble("PTMin", 10.0);
93
94 // to be tuned, since FS and delphes have different isolation profiles
95 fRelIsoMax = GetDouble("fRelIsoMax", 0.3);
96
97 // create output array
98 fOutputArray = ExportArray(GetString("OutputArray", "photons"));
99}
100
101//------------------------------------------------------------------------------
102
103void PhotonID::Finish()
104{
105 if(fItInputPhotonArray) delete fItInputPhotonArray;
106 if(fItInputGenArray) delete fItInputGenArray;
107}
108
109//------------------------------------------------------------------------------
110
111void PhotonID::Process()
112{
113 Candidate *candidate, *mother;
114 Double_t pt, eta, phi, e;
115 Double_t relIso;
116 Bool_t isolated;
117
118 //cout<< "---- new event ---------"<<endl;
119
120 fItInputPhotonArray->Reset();
121 while((candidate = static_cast<Candidate *>(fItInputPhotonArray->Next())))
122 {
123
124 mother = candidate;
125 candidate = static_cast<Candidate *>(candidate->Clone());
126 candidate->AddCandidate(mother);
127
128 const TLorentzVector &candidatePosition = candidate->Position;
129 const TLorentzVector &candidateMomentum = candidate->Momentum;
130 eta = candidatePosition.Eta();
131 phi = candidatePosition.Phi();
132 pt = candidateMomentum.Pt();
133 e = candidateMomentum.E();
134
135 if(pt < fPTMin) continue;
136
137 //cout<< " ---- photon -----: "<<pt<<","<<eta<<","<<phi<<endl;
138
139 // find out if photon matches does not match photon in gen collection and apply fae efficiency
140 if(isFake(candidate))
141 {
142 //cout<<" Fake!"<<endl;
143
144 if(gRandom->Uniform() > fFakeFormula->Eval(pt, eta, phi, e)) continue;
145 //cout<<" passed"<<endl;
146 candidate->Status = 3;
147 fOutputArray->Add(candidate);
148 }
149
150 // if matches photon in gen collection
151 else
152 {
153 relIso = candidate->IsolationVar;
154 isolated = (relIso < 0.3);
155 //cout<<" Prompt!: "<<relIso<<endl;
156
157 // if isolated apply prompt formula
158 if(isolated)
159 {
160 //cout<<" isolated!: "<<relIso<<endl;
161 if(gRandom->Uniform() > fPromptFormula->Eval(pt, eta, phi, e)) continue;
162 //cout<<" passed"<<endl;
163 candidate->Status = 1;
164 fOutputArray->Add(candidate);
165 }
166
167 // if non-isolated apply non-prompt formula
168 else
169 {
170 //cout<<" non-isolated!: "<<relIso<<endl;
171 if(gRandom->Uniform() > fNonPromptFormula->Eval(pt, eta, phi, e)) continue;
172 //cout<<" passed"<<endl;
173 candidate->Status = 2;
174 fOutputArray->Add(candidate);
175 }
176 }
177 }
178}
179
180//------------------------------------------------------------------------------
181
182Bool_t PhotonID::isFake(const Candidate *obj)
183{
184
185 const TLorentzVector &mom_rec = obj->Momentum;
186
187 Bool_t matches = false;
188 fItInputGenArray->Reset();
189 Candidate *gen;
190
191 while((gen = static_cast<Candidate *>(fItInputGenArray->Next())))
192 {
193 const TLorentzVector &mom_gen = gen->Momentum;
194 Int_t status = gen->Status;
195 Int_t pdgCode = TMath::Abs(gen->PID);
196 Float_t dPtOverPt = TMath::Abs((mom_gen.Pt() - mom_rec.Pt()) / mom_rec.Pt());
197 Float_t deltaR = mom_gen.DeltaR(mom_rec);
198
199 if(status != 1) continue;
200 if(pdgCode != 22) continue;
201 if(dPtOverPt > 0.5) continue;
202 if(deltaR > 0.1) continue;
203
204 matches = true;
205 break;
206 }
207
208 return !matches;
209}
Note: See TracBrowser for help on using the repository browser.