18 #ifdef XAOD_STANDALONE
23 #include "GaudiKernel/SmartIF.h"
27 #include <boost/program_options.hpp>
28 namespace po = boost::program_options;
55 #define CHECK_RETRIEVE( container , name ) { \
56 if (!event.retrieve( container , name ).isSuccess()){ \
57 Error( APP_NAME , "Could not load event %s from the file!" , name ); \
58 throw std::runtime_error("Container retrieval failed"); \
66 for (
size_t n=0;
n<
p->nChildren();++
n){
76 for (
size_t n=0;
n<
p->nParents();++
n){
89 using namespace asg::msgUserCode;
93 #ifdef XAOD_STANDALONE
95 throw std::runtime_error(
"Cannot initialise xAOD access !");
104 po::options_description
desc(
"Running a simple truth analysis");
107 (
"help,h",
"print usage and exit")
108 (
"output,o", po::value<std::string>(&
outputName),
"Output file name")
109 (
"input,i", po::value<std::string>(&inputName),
"Input file name")
110 (
"nevents", po::value<int>()->default_value(-1),
"number of events to run on (set to -1 to ignore it")
113 po::variables_map vm;
117 if (vm.count(
"help")) {
118 std::cout <<
desc << std::endl;
123 if (vm.count(
"nevents")) {
124 nevents = vm[
"nevents"].as<
int>();
127 Info(
"Reading from %s and writing to %s" , inputName.c_str() ,
outputName.c_str() );
130 std::unique_ptr< TFile >
ifile( TFile::Open( inputName.c_str() ,
"READ" ) );
132 #ifdef XAOD_STANDALONE
143 const size_t nParticleContainers = 9;
144 std::string particleKeyList[nParticleContainers] = {
"TruthElectrons",
"TruthMuons",
"TruthPhotons",
"TruthTaus",
"TruthNeutrinos",
"TruthBSM",
"TruthBottom",
"TruthTop",
"TruthBoson" };
152 for (
size_t n=0;
n<nParticleContainers;++
n) truthParticles[
n] =
nullptr;
159 TH1D* h_metNonInt =
new TH1D(
"MET_NonInt",
"",50,0,500.);
160 TH1D* h_metInt =
new TH1D(
"MET_Int",
"",50,0,200.);
161 TH1D* h_metIntOut =
new TH1D(
"MET_IntOut",
"",50,0,100.);
162 TH1D* h_metIntMuons =
new TH1D(
"MET_IntMuons",
"",50,0,300.);
164 TH1D* h_partPt[nParticleContainers];
165 for (
size_t n=0;
n<nParticleContainers;++
n) h_partPt[
n] =
new TH1D((particleKeyList[
n]+
"_pT").c_str(),
"",50,0.,500.);
166 TH1D* h_partConn[nParticleContainers];
167 for (
size_t n=0;
n<nParticleContainers;++
n) h_partConn[
n] =
new TH1D((particleKeyList[
n]+
"_Connections").c_str(),
"",35,-10,25);
169 TH1D* h_elPtCone =
new TH1D(
"ElPtCone",
"",50,0.,20.);
170 TH1D* h_elEtCone =
new TH1D(
"ElEtCone",
"",50,0.,20.);
171 TH1D* h_muPtCone =
new TH1D(
"MuPtCone",
"",50,0.,20.);
172 TH1D* h_muEtCone =
new TH1D(
"MuEtCone",
"",50,0.,20.);
173 TH1D* h_phPtCone =
new TH1D(
"PhPtCone",
"",50,0.,20.);
174 TH1D* h_phEtCone =
new TH1D(
"PhEtCone",
"",50,0.,20.);
176 TH1D* h_elType =
new TH1D(
"ElType",
"",50,0,50);
177 TH1D* h_elOrig =
new TH1D(
"ElOrigin",
"",50,0,50);
178 TH1D* h_muType =
new TH1D(
"MuType",
"",50,0,50);
179 TH1D* h_muOrig =
new TH1D(
"MuOrigin",
"",50,0,50);
180 TH1D* h_phType =
new TH1D(
"PhType",
"",50,0,50);
181 TH1D* h_phOrig =
new TH1D(
"PhOrigin",
"",50,0,50);
182 TH1D* h_taType =
new TH1D(
"TaType",
"",50,0,50);
183 TH1D* h_taOrig =
new TH1D(
"TaOrigin",
"",50,0,50);
184 TH1D* h_nuType =
new TH1D(
"NuType",
"",50,0,50);
185 TH1D* h_nuOrig =
new TH1D(
"NuOrigin",
"",50,0,50);
187 TH1D* h_elDressedPt =
new TH1D(
"ElDressedPt",
"",50,0.,200.);
188 TH1D* h_muDressedPt =
new TH1D(
"MuDressedPt",
"",50,0.,200.);
189 TH1D* h_taDressedPt =
new TH1D(
"TaDressedPt",
"",50,0.,200.);
191 TH1D* h_jetPt =
new TH1D(
"JetPt",
"",50,0,500.);
192 TH1D* h_jetFTAG =
new TH1D(
"JetFTAG",
"",25,0,25);
193 TH1D* h_jetJeMe =
new TH1D(
"JetJeMe",
"",25,0,25);
194 TH1D* h_jetFull =
new TH1D(
"JetFull",
"",50,0,50);
196 TH1D* h_jetLRPt =
new TH1D(
"JetLRPt",
"",50,0,500.);
197 TH1D* h_jetLRT2 =
new TH1D(
"JetLRT2",
"",50,0,500.);
198 TH1D* h_jetLRD2 =
new TH1D(
"JetLRD2",
"",50,0,500.);
200 TH1D* h_x1 =
new TH1D(
"X1",
"",50,0,1.);
201 TH1D* h_x2 =
new TH1D(
"X2",
"",50,0,1.);
202 TH1D* h_xsec =
new TH1D(
"XSec",
"",50,0,1.);
203 TH1D* h_HTFilt =
new TH1D(
"HTFilt",
"",50,0,1000.);
209 std::vector<TH1D*> h_weights;
212 std::vector<std::string> weightNames;
226 CHECK_RETRIEVE( largeRJets ,
"AntiKt10TruthTrimmedPtFrac5SmallR20Jets" )
228 for (
size_t n=0;
n<nParticleContainers;++
n){
233 Error(
APP_NAME ,
"Could not load event TruthMetaData from the file!" );
234 throw std::runtime_error(
"Could not load event TruthMetaData from the file!");
238 if (truthMeta->
size()>1){
239 Error(
APP_NAME ,
"Truth metadata size: %lu . No one will look past item 0!" , truthMeta->
size() );
240 throw std::runtime_error(
"Truth metadata size >1");
242 if (channelNumber==0){
243 channelNumber = (*truthMeta)[0]->mcChannelNumber();
244 weightNames = std::vector<std::string>((*truthMeta)[0]->weightNames());
245 Info(
APP_NAME ,
"Got channel number %u and %lu weight names" , channelNumber , weightNames.size() );
249 throw std::runtime_error(
"Channel number changed mid-file!");
251 if (weightNames != (*truthMeta)[0]->weightNames() ){
253 for (
size_t n=0;
n<
std::max( weightNames.size() , (*truthMeta)[0]->weightNames().size() );++
n){
254 std::cerr <<
" " <<
n <<
" ";
255 if (
n<weightNames.size()) std::cerr << weightNames[
n] <<
" ";
256 else std::cerr <<
"- ";
257 if (
n<(*truthMeta)[0]->weightNames().size()) std::cerr << (*truthMeta)[0]->weightNames()[
n] <<
" ";
258 else std::cerr <<
"- ";
259 std::cerr << std::endl;
261 throw std::runtime_error(
"Weights have changed!");
264 if (h_weights.size()==0){
266 for (
auto &
name : weightNames ){
267 h_weights.push_back(
new TH1D((
"h_W_"+
name).c_str(),
"",100,-10.,10.) );
269 h_weights.push_back(
new TH1D(
"h_W_nominalTest",
"",100,-10.,10.) );
271 for (
size_t n=0;
n<weightNames.size();++
n) h_weights[
n]->Fill( weightTool->
getWeight(xEventInfo,weightNames[
n]) );
273 h_weights[weightNames.size()]->Fill( weightTool->
getWeight(xEventInfo,
" nominal ") );
280 h_xsec->Fill( (*xTruthEventContainer)[0]->
crossSection() );
281 h_HTFilt->Fill( xEventInfo->
auxdata<
float>(
"GenFiltHT") );
283 h_metNonInt->Fill( (*truthMET)[
"NonInt"]->
met()*0.001 );
284 h_metNonInt->Fill( (*truthMET)[
"Int"]->
met()*0.001 );
285 h_metNonInt->Fill( (*truthMET)[
"IntOut"]->
met()*0.001 );
286 h_metNonInt->Fill( (*truthMET)[
"IntMuons"]->
met()*0.001 );
288 for (
size_t n=0;
n<nParticleContainers;++
n){
289 for (
const auto *
p : *truthParticles[
n]){
290 h_partPt[
n]->Fill(
p->pt() );
295 for (
const auto *
p : *truthParticles[0]){
296 h_elPtCone->Fill(
p->auxdata<
float>(
"ptcone30")*0.001 );
297 h_elEtCone->Fill(
p->auxdata<
float>(
"etcone20")*0.001 );
298 h_elDressedPt->Fill(
p->auxdata<
float>(
"pt_dressed")*0.001 );
299 h_elType->Fill(
p->auxdata<
unsigned int>(
"classifierParticleType") );
300 h_elOrig->Fill(
p->auxdata<
unsigned int>(
"classifierParticleOrigin") );
302 for (
const auto *
p : *truthParticles[1]){
303 h_muPtCone->Fill(
p->auxdata<
float>(
"ptcone30")*0.001 );
304 h_muEtCone->Fill(
p->auxdata<
float>(
"etcone20")*0.001 );
305 h_muDressedPt->Fill(
p->auxdata<
float>(
"pt_dressed")*0.001 );
306 h_muType->Fill(
p->auxdata<
unsigned int>(
"classifierParticleType") );
307 h_muOrig->Fill(
p->auxdata<
unsigned int>(
"classifierParticleOrigin") );
309 for (
const auto *
p : *truthParticles[2]){
310 h_phPtCone->Fill(
p->auxdata<
float>(
"ptcone20")*0.001 );
311 h_phEtCone->Fill(
p->auxdata<
float>(
"etcone20")*0.001 );
312 h_phType->Fill(
p->auxdata<
unsigned int>(
"classifierParticleType") );
313 h_phOrig->Fill(
p->auxdata<
unsigned int>(
"classifierParticleOrigin") );
315 for (
const auto *
p : *truthParticles[3]){
316 h_taDressedPt->Fill(
p->auxdata<
float>(
"pt_vis_dressed")*0.001 );
317 h_taType->Fill(
p->auxdata<
unsigned int>(
"classifierParticleType") );
318 h_taOrig->Fill(
p->auxdata<
unsigned int>(
"classifierParticleOrigin") );
320 for (
const auto *
p : *truthParticles[4]){
321 h_nuType->Fill(
p->auxdata<
unsigned int>(
"classifierParticleType") );
322 h_nuOrig->Fill(
p->auxdata<
unsigned int>(
"classifierParticleOrigin") );
324 for (
const auto * j : *smallRJets){
325 h_jetPt->Fill( j->pt()*0.001 );
326 h_jetFTAG->Fill( j->auxdata<
int>(
"HadronConeExclTruthLabelID") );
327 h_jetJeMe->Fill( j->auxdata<
int>(
"PartonTruthLabelID") );
328 h_jetFull->Fill( j->auxdata<
int>(
"TrueFlavor") );
330 for (
const auto * j : *largeRJets){
331 h_jetLRPt->Fill( j->pt()*0.001 );
332 h_jetLRT2->Fill( j->auxdata<
float>(
"Tau2_wta") );
333 h_jetLRD2->Fill( j->auxdata<
float>(
"D2") );
337 Info(
APP_NAME ,
"Done with event loop" );
340 TFile * oRoot =
new TFile(
outputName.c_str(),
"RECREATE");
345 h_metNonInt->Write();
347 h_metIntOut->Write();
348 h_metIntMuons->Write();
350 for (
size_t n=0;
n<nParticleContainers;++
n) h_partPt[
n]->Write();
351 for (
size_t n=0;
n<nParticleContainers;++
n) h_partConn[
n]->Write();
358 h_elDressedPt->Write();
359 h_muDressedPt->Write();
360 h_taDressedPt->Write();
384 for (
auto *
h : h_weights)
h->Write();
390 #ifndef XAOD_STANDALONE
391 if (app->finalize().isFailure()){
392 Warning(
APP_NAME ,
"Finalization failed?" );
396 Info(
APP_NAME ,
"All done -- goodbye." );