ATLAS Offline Software
fbtNtupleTester.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include <AsgMessaging/MessageCheck.h>
6 #include <AsgTools/StandaloneToolHandle.h>
7 
8 #include <iostream>
9 #include <sstream>
10 #include <memory>
11 #include <vector>
12 #include <string>
13 #include <sstream>
14 #include <cmath>
15 #include <numeric>
16 
17 #include <TFile.h>
18 #include <TChain.h>
19 #include <TString.h>
20 #include <TSystem.h>
21 
22 #ifdef XAOD_STANDALONE
23 #include "xAODRootAccess/Init.h"
24 #include "xAODRootAccess/TEvent.h"
25 #include "xAODRootAccess/TStore.h"
26 #endif
27 
28 #include "AsgTools/AnaToolHandle.h"
29 #include "xAODEventInfo/EventInfo.h"
30 #include "xAODEventInfo/EventAuxInfo.h"
31 #include "xAODBase/IParticle.h"
32 #include "xAODEgamma/Electron.h"
33 #include "xAODMuon/Muon.h"
34 #include "FakeBkgTools/ApplyFakeFactor.h"
35 
36 #ifndef READ_TREE_ADDRESSES
37 #define READ_TREE_ADDRESSES(typeName, name) \
38  typeName name = 0; chain->SetBranchAddress(#name, &name)
39 #endif
40 
41 #ifndef READ_TREE_ADDRESSES_VEC
42 #define READ_TREE_ADDRESSES_VEC(typeName, name) \
43  typeName* name = 0; chain->SetBranchAddress(#name, &name)
44 #endif
45 
46 
47 
48 // messaging
49 ANA_MSG_HEADER(Test)
50 ANA_MSG_SOURCE(Test, "fbtNtupleTester")
51 using namespace Test;
52 
53 // an example input file is /afs/cern.ch/user/j/jreicher/public/fbtTestData.root
54 int main(int argc, char* argv[])
55 {
56  ANA_CHECK_SET_TYPE (int); // makes ANA_CHECK return ints if exiting function
57 
58  if(argc < 2){
59  std::cerr << "No input file specified! Exiting." << std::endl;
60  return EXIT_FAILURE;
61  }
62 
63  TString fileNameFullPath = argv[1];
64  TString fileName = fileNameFullPath;
65 
66 #ifdef XAOD_STANDALONE
67 
68 
69  TFile *inFile = TFile::Open(fileName, "READ");
70  if(!inFile){
71  std::cerr << "Invalid root file! Exiting." << std::endl;
72  return EXIT_FAILURE;
73  }
74 
75  // Name of the tree in the file
76  TString sample = "tree_NoSys";
77 
78  std::cout << "Running sample: " << sample << std::endl;
79 
80  // Create instance of the tools and set the FF file names + the histogram names
81  CP::ApplyFakeFactor* ffTool = new CP::ApplyFakeFactor("ApplyFF");
82 
83  // Path to the file of FF histograms
84  std::vector<std::string> ff_file = { "FakeBkgTools/testValuesFF.root" };
85 
86  ANA_CHECK( ffTool->setProperty("InputFiles", ff_file) );
87  ANA_CHECK( ffTool->setProperty("EnergyUnit", "GeV") );
88  ANA_CHECK( ffTool->setProperty("SkipUncertainties", true) );
89  ANA_CHECK( ffTool->setProperty("OutputLevel", MSG::FATAL) );
90 
91  // ntuple values are in GeV, but xAOD::IParticles that
92  // we create will be in MeV
93  float convertToMeV = 1000.;
94 
95  // Initialize tools
96  ANA_CHECK( ffTool->initialize() );
97  TChain* chain = new TChain("chain");
98  std::vector<TString> datasets = {};
99 
100  datasets = {sample};
101 
102  for(TString dataset : datasets){
103  chain->Add(fileNameFullPath+"/"+dataset);
104  }
105 
106  READ_TREE_ADDRESSES(unsigned long long, EventNumber);
107  READ_TREE_ADDRESSES(int, DatasetNumber);
108  READ_TREE_ADDRESSES(int, nLep_base);
109  READ_TREE_ADDRESSES(int, nLep_signal);
110  READ_TREE_ADDRESSES(int, nLep_antiID);
111 
112  READ_TREE_ADDRESSES(float, lep1Pt);
113  READ_TREE_ADDRESSES(float, lep1Eta);
114  READ_TREE_ADDRESSES(float, lep1Phi);
115  READ_TREE_ADDRESSES(float, lep1D0Sig);
116  READ_TREE_ADDRESSES(float, lep1Z0SinTheta);
117  READ_TREE_ADDRESSES(int, lep1Flavor);
118  READ_TREE_ADDRESSES(int, lep1Charge);
119  READ_TREE_ADDRESSES(int, lep1PassOR);
120  READ_TREE_ADDRESSES(int, lep1PassBL);
121  READ_TREE_ADDRESSES(bool, lep1Loose);
122  READ_TREE_ADDRESSES(bool, lep1Medium);
123  READ_TREE_ADDRESSES(bool, lep1TruthMatched);
124  READ_TREE_ADDRESSES(bool, lep1Signal);
125  READ_TREE_ADDRESSES(bool, lep1AntiID);
126  READ_TREE_ADDRESSES(int, lep1Type);
127 
128  READ_TREE_ADDRESSES(float, lep2Pt);
129  READ_TREE_ADDRESSES(float, lep2Eta);
130  READ_TREE_ADDRESSES(float, lep2Phi);
131  READ_TREE_ADDRESSES(float, lep2D0Sig);
132  READ_TREE_ADDRESSES(float, lep2Z0SinTheta);
133  READ_TREE_ADDRESSES(int, lep2Flavor);
134  READ_TREE_ADDRESSES(int, lep2Charge);
135  READ_TREE_ADDRESSES(bool, lep2PassOR);
136  READ_TREE_ADDRESSES(bool, lep2PassBL);
137  READ_TREE_ADDRESSES(bool, lep2Loose);
138  READ_TREE_ADDRESSES(bool, lep2Medium);
139  READ_TREE_ADDRESSES(bool, lep2TruthMatched);
140  READ_TREE_ADDRESSES(bool, lep2Signal);
141  READ_TREE_ADDRESSES(bool, lep2AntiID);
142  READ_TREE_ADDRESSES(int, lep2Type);
143 
144  READ_TREE_ADDRESSES(float, lep3Pt);
145  READ_TREE_ADDRESSES(float, lep3Eta);
146  READ_TREE_ADDRESSES(float, lep3Phi);
147  READ_TREE_ADDRESSES(float, lep3D0Sig);
148  READ_TREE_ADDRESSES(float, lep3Z0SinTheta);
149  READ_TREE_ADDRESSES(int, lep3Flavor);
150  READ_TREE_ADDRESSES(int, lep3Charge);
151  READ_TREE_ADDRESSES(bool, lep3PassOR);
152  READ_TREE_ADDRESSES(bool, lep3PassBL);
153  READ_TREE_ADDRESSES(bool, lep3Loose);
154  READ_TREE_ADDRESSES(bool, lep3Medium);
155  READ_TREE_ADDRESSES(bool, lep3TruthMatched);
156  READ_TREE_ADDRESSES(bool, lep3Signal);
157  READ_TREE_ADDRESSES(bool, lep3AntiID);
158  READ_TREE_ADDRESSES(int, lep3Type);
159 
160  READ_TREE_ADDRESSES(double, FFWeight);
161 
162  // Turn on all branches. Our cloned tree will share
163  // the references to the variables which we wish to modify!
164  chain->SetBranchStatus("*",1);
165 
166  TFile* outFile_nominal = new TFile("antiIDTimesFF_"+sample+"_nominal.root","RECREATE");
167  TTree* outTree_nominal = chain->CloneTree(0);
168  outTree_nominal->SetTitle("fakes");
169  outTree_nominal->SetName("fakes");
170 
171  // We need to make a TEvent and TStore object so that the FBT can look
172  // for the evtStore()
173  xAOD::TEvent* tEvent = new xAOD::TEvent();
174  xAOD::TStore* tStore = new xAOD::TStore();
175 
176  // For creating an EventInfo object that the tool will look for
177  std::unique_ptr<xAOD::EventInfo> eventInfo = std::make_unique<xAOD::EventInfo>();
178  std::unique_ptr<xAOD::EventAuxInfo> eventAuxInfo = std::make_unique<xAOD::EventAuxInfo>();
179 
180  eventInfo->setStore(eventAuxInfo.get());
181  ANA_CHECK( ffTool->evtStore()->record(eventInfo.get(), "EventInfo") );
182  ANA_CHECK( ffTool->evtStore()->record(eventAuxInfo.get(), "EventInfoAux.") );;
183 
184  // we can clear this vector in each loop, rather than allocate
185  // the memory for each event
186  xAOD::IParticleContainer particles;
187 
188  static const SG::Decorator<char> tightDecor("Tight");
189  static const SG::Decorator<unsigned int> runNumberDecor("RandomRunNumber");
190 
191  // Loop over events!
192  int nEventsProcessed(0);
193  int nEventsPassed(0);
194 
195  Long64_t nentries = chain->GetEntries();
196  for(Long64_t i(0); i < nentries; ++i){
197  std::cout << "processed " << nEventsProcessed << " events" << std::endl;
198  ++nEventsProcessed;
199 
200  int getEntry = chain->GetEntry(i);
201  if(getEntry < 1){
202  std::cerr << "GetEntry() error!! Either an I/O issue, or an invalid entry. Aborting." << std::endl;
203  abort();
204  }
205 
206  // just as an example for an EventInfo decor that one could set...
207  runNumberDecor(*eventInfo) = 507;
208 
209  std::vector<float> lepPt = {lep1Pt, lep2Pt, lep3Pt};
210  std::vector<float> lepEta = {lep1Eta, lep2Eta, lep3Eta};
211  std::vector<float> lepPhi = {lep1Phi, lep2Phi, lep3Phi};
212  std::vector<float> lepD0Sig = {lep1D0Sig, lep2D0Sig, lep3D0Sig};
213  std::vector<float> lepZ0SinTheta = {lep1Z0SinTheta, lep2Z0SinTheta, lep3Z0SinTheta};
214  std::vector<int> lepFlavor = {lep1Flavor, lep2Flavor, lep3Flavor};
215  std::vector<int> lepCharge = {lep1Charge, lep2Charge, lep3Charge};
216  std::vector<int> lepPassOR = {lep1PassOR, lep2PassOR, lep3PassOR};
217  std::vector<int> lepPassBL = {lep1PassBL, lep2PassBL, lep3PassBL};
218  std::vector<int> lepLoose = {lep1Loose, lep2Loose, lep3Loose};
219  std::vector<int> lepMedium = {lep1Medium, lep2Medium, lep3Medium};
220  std::vector<int> lepSignal = {lep1Signal, lep2Signal, lep3Signal};
221  std::vector<int> lepTruthMatched = {lep1TruthMatched, lep2TruthMatched, lep3TruthMatched};
222 
224  FFWeight = 1;
225  nLep_antiID = 0;
226  lep1AntiID = 0;
227  lep2AntiID = 0;
228  lep3AntiID = 0;
229 
230  // Right now, we have up to 3 leptons saved in our ntuples
231  int nTotalLeps = 3;
232 
233  // Look at ID and isolation to fill whether a given lepton passes
234  // the signal ID or antiID
235  std::vector<int> lepPassSigID(nTotalLeps, 0);
236  std::vector<int> lepPassAntiID(nTotalLeps, 0);
237 
238  // clear the vector of IParticle objects
239  for(auto part : particles){
240  delete part;
241  }
242  particles.clear();
243 
244  // Loop over the events to assign ID vs. antiID leptons
245  for(int i = 0; i < nTotalLeps; ++i){
246 
247  bool passesSignal = false;
248  bool passesAntiID = false;
249  xAOD::Type::ObjectType lepType = xAOD::Type::Other;
250 
251  // Fill ID and SF vectors below:
252  if( (lepFlavor).at(i) == 1 ){ // electrons
253 
254  //double z0sinTheta = fabs( (lepZ0SinTheta).at(i) );
255 
256  passesSignal = (lepSignal).at(i);
257  passesAntiID = !(lepSignal).at(i);
258  //passesAntiID = (lepLoose.at(i) && lepPassBL.at(i) && (z0sinTheta < 0.5) && (lepPassOR).at(i)) && !(lepSignal.at(i));
259 
260  // Now fill the vectors
261  lepPassSigID.at(i) = passesSignal;
262  lepPassAntiID.at(i) = passesAntiID;
263 
264  lepType = xAOD::Type::Electron;
265  }
266  else if( (lepFlavor).at(i) == 2){ // muons
267 
268  //double z0sinTheta = fabs( (lepZ0SinTheta).at(i) );
269 
270  passesSignal = (lepSignal).at(i);
271  passesAntiID = !(lepSignal).at(i);
272  //passesAntiID = (lepMedium.at(i) && (z0sinTheta < 0.5) && !(lepSignal).at(i));
273 
274  // Now fill the vectors
275  lepPassSigID.at(i) = passesSignal;
276  lepPassAntiID.at(i) = passesAntiID;
277 
278  lepType = xAOD::Type::Muon;
279  }
280 
281  // if this lepton isn't relevant, we can skip
282  if(!passesSignal && !passesAntiID) continue;
283 
284  // FIXME need to be sure all pT/eta/phi values are in MeV, and then just need to
285  // unpack the FFs into MeV
286  if (lepType == xAOD::Type::Electron){
287  xAOD::Electron* particle = new xAOD::Electron();
288  particle->makePrivateStore();
289  particle->setP4(lepPt.at(i)*convertToMeV,lepEta.at(i),lepPhi.at(i),0.511);
290  particle->setCharge(lepCharge.at(i));
291  tightDecor(*particle) = passesSignal;
292  particles.push_back(static_cast<xAOD::IParticle*>(particle));
293  }
294  else if(lepType == xAOD::Type::Muon){
295  xAOD::Muon* particle = new xAOD::Muon();
296  particle->makePrivateStore();
297  particle->setP4(lepPt.at(i)*convertToMeV,lepEta.at(i),lepPhi.at(i));
298  particle->setCharge(lepCharge.at(i));
299  tightDecor(*particle) = passesSignal;
300  particles.push_back(static_cast<xAOD::IParticle*>(particle));
301  }
302  else{
303  std::cerr << "invalid lepton type!" << std::endl;
304  continue;
305  }
306 
307  // Here we'll redefine "signal" lepton for our outputs, because
308  // we're about to promote our antiID leptons to signal leptons.
309  if(passesAntiID){
310 
311  ++nLep_signal;
312  ++nLep_antiID;
313 
314  if(i == 0){
315  lep1Signal = 1;
316  lep1AntiID = 1;
317  }
318  else if(i == 1){
319  lep2Signal = 1;
320  lep2AntiID = 1;
321  }
322  else if(i == 2){
323  lep3Signal = 1;
324  lep3AntiID = 1;
325  }
326  }
327 
328  } // end loop over nTotalLeps
329 
330  // Event is not relevant for fake estimation!
331  if(nLep_antiID == 0) continue;
332 
333  // Now actually compute the weight!
334  float fbtWeight;
335  ANA_CHECK( ffTool->addEvent(particles) );;
336  ANA_CHECK( ffTool->getEventWeight(fbtWeight, "3T", ">=1F[T]") );
337 
338  // FIXME note that the relevant MC events will need their
339  // FFWeight to be multiplied by -1 here, since these are used
340  // for the prompt MC subtraction
341  FFWeight = fbtWeight;
342  std::cout << "nLep_antiID " << nLep_antiID << std::endl;
343  std::cout << "FFWeight " << FFWeight << std::endl;
344 
345  // Fill the tree
346  outTree_nominal->Fill();
347 
348  // just to keep track
349  ++nEventsPassed;
350 
351  } // end loop over events
352 
353  std::cout << "nEventsProcessed is " << nEventsProcessed << " and nEventsPassed is " << nEventsPassed << std::endl;
354 
355  outFile_nominal->Write();
356 
357  delete chain;
358  delete ffTool;
359  delete outFile_nominal;
360  delete tStore;
361  delete tEvent;
362 
363 #endif
364 
365  return EXIT_SUCCESS;
366 }
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:76
xAOD::Electron
Electron_v1 Electron
Definition of the current "egamma version".
Definition: Event/xAOD/xAODEgamma/xAODEgamma/Electron.h:17
runLayerRecalibration.chain
chain
Definition: runLayerRecalibration.py:175
READ_TREE_ADDRESSES
#define READ_TREE_ADDRESSES(typeName, name)
Definition: fbtNtupleTester.cxx:37
IParticle.h
python.Constants.FATAL
int FATAL
Definition: Control/AthenaCommon/python/Constants.py:19
ObjectType
ObjectType
Definition: BaseObject.h:11
Muon.h
makeTOC.inFile
string inFile
Definition: makeTOC.py:5
ANA_CHECK
#define ANA_CHECK(EXP)
check whether the given expression was successful
Definition: Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:324
CP::ApplyFakeFactor
Definition: ApplyFakeFactor.h:22
sct_calib_tf.EventNumber
int EventNumber
Definition: sct_calib_tf.py:29
xAOD::IParticle
Class providing the definition of the 4-vector interface.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:41
CP::BaseLinearFakeBkgTool::getEventWeight
virtual StatusCode getEventWeight(float &weight, const std::string &selection, const std::string &process) override final
returns an event weight addEvent() must have been called before hand.
Definition: BaseLinearFakeBkgTool.cxx:39
xAOD::Muon_v1
Class describing a Muon.
Definition: Muon_v1.h:38
ANA_MSG_HEADER
#define ANA_MSG_HEADER(NAME)
for standalone code this creates a new message category
Definition: Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:113
AthCommonDataStore::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
PlotCalibFromCool.nentries
nentries
Definition: PlotCalibFromCool.py:798
FortranAlgorithmOptions.fileName
fileName
Definition: FortranAlgorithmOptions.py:13
SG::Decorator< char >
FullCPAlgorithmsTest_eljob.sample
sample
Definition: FullCPAlgorithmsTest_eljob.py:116
lumiFormat.i
int i
Definition: lumiFormat.py:85
SG::AuxElement::setStore
void setStore(const SG::IConstAuxStore *store)
Set the store associated with this object.
Definition: AuxElement.cxx:221
LArCellNtuple.argv
argv
Definition: LArCellNtuple.py:152
MessageCheck.h
macros for messaging and checking status codes
TEvent.h
main
int main(int argc, char *argv[])
Definition: fbtNtupleTester.cxx:54
Init.h
dataset
Definition: dataset.h:27
DQHistogramMergeRegExp.argc
argc
Definition: DQHistogramMergeRegExp.py:20
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
AnaToolHandle.h
ANA_MSG_SOURCE
#define ANA_MSG_SOURCE(NAME, TITLE)
the source code part of ANA_MSG_SOURCE
Definition: Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:133
EventAuxInfo.h
xAOD::Muon
Muon_v1 Muon
Reference the current persistent version:
Definition: Event/xAOD/xAODMuon/xAODMuon/Muon.h:13
xAOD::Electron_v1
Definition: Electron_v1.h:34
EventInfo.h
Muon
struct TBPatternUnitContext Muon
xAOD::TStore
A relatively simple transient store for objects created in analysis.
Definition: TStore.h:44
ANA_CHECK_SET_TYPE
#define ANA_CHECK_SET_TYPE(TYPE)
set the type for ANA_CHECK to report failures
Definition: Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:314
LArG4FSStartPointFilter.particles
list particles
Definition: LArG4FSStartPointFilter.py:84
Electron.h
StandaloneToolHandle.h
CP::BaseFakeBkgTool::addEvent
virtual StatusCode addEvent(const xAOD::IParticleContainer &particles, float extraWeight=1.f) override final
supply list of leptons / global variables, internal counters incremented Does not return anything; ev...
Definition: BaseFakeBkgTool.cxx:238
ApplyFakeFactor.h
TStore.h
xAOD::TEvent
Tool for accessing xAOD files outside of Athena.
Definition: Control/xAODRootAccess/xAODRootAccess/TEvent.h:84
CP::ApplyFakeFactor::initialize
virtual StatusCode initialize() override
Dummy implementation of the initialisation function.
Definition: ApplyFakeFactor.cxx:32
Generated on Thu Nov 7 2024 21:15:12 for ATLAS Offline Software by   doxygen 1.8.18