ATLAS Offline Software
Loading...
Searching...
No Matches
Classes | Functions | Variables
rmain.cxx File Reference
#include <cstdio>
#include <cmath>
#include <fstream>
#include <iostream>
#include <vector>
#include <map>
#include <set>
#include <sstream>
#include <execinfo.h>
#include <signal.h>
#include <unistd.h>
#include "TChain.h"
#include "TFile.h"
#include "TTree.h"
#include "TDirectory.h"
#include "TrigInDetAnalysis/Track.h"
#include "TrigInDetAnalysis/TIDAEvent.h"
#include "TrigInDetAnalysis/TrackSelector.h"
#include "TrigInDetAnalysis/TIDAVertex.h"
#include "TrigInDetAnalysisUtils/Associator_BestMatch.h"
#include "TrigInDetAnalysisUtils/Filters.h"
#include "TrigInDetAnalysisUtils/Filter_Offline2017.h"
#include "TrigInDetAnalysisUtils/Filter_OfflineR22.h"
#include "TrigInDetAnalysisExample/NtupleTrackSelector.h"
#include "TrigInDetAnalysisExample/ChainString.h"
#include "TrigInDetAnalysisUtils/Associator_TruthMatch.h"
#include "TrigInDetAnalysis/Efficiency1D.h"
#include "TrigInDetAnalysis/TIDARoiDescriptor.h"
#include "TrigInDetAnalysis/TrigObjectMatcher.h"
#include "ReadCards.h"
#include "utils.h"
#include "RoiFilter.h"
#include "ConfAnalysis.h"
#include "PurityAnalysis.h"
#include "ConfVtxAnalysis.h"
#include "TIDAReference.h"
#include "lumiList.h"
#include "lumiParser.h"
#include "dataset.h"
#include "event_selector.h"
#include "BinConfig.h"
#include "zbeam.h"
#include "computils.h"
#include "globals.h"
#include "TrigInDetAnalysisUtils/TagNProbe.h"

Go to the source code of this file.

Classes

struct  event_list

Functions

void copyReleaseInfo (TTree *tree, TFile *foutdir)
void handler (int sig)
 signal handler
std::string time_str ()
int atoi_check (const std::string &s)
 convert string to integer with check if successful
int GetChainAuthor (std::string chainName)
 Return selected author, expects format L2_e20_medium:TrigL2SiTrackFinder:3 where 3 is the track author NB: this isn't actually needed at all.
TIDARoiDescriptor makeCustomRefRoi (const TIDARoiDescriptor &roi, double etaHalfWidth=-999, double phiHalfWidth=-999, double zedHalfWidth=-999)
bool trackPtGrtr (TIDA::Track *trackA, TIDA::Track *trackB)
template<typename T>
std::ostream & operator<< (std::ostream &s, const std::vector< T * > &v)
template<typename T>
std::ostream & operator<< (std::ostream &s, const std::vector< T > &v)
const std::vector< TIDA::Track > ibl_filter (const std::vector< TIDA::Track > &tv)
const std::vector< TIDA::Track > replaceauthor (const std::vector< TIDA::Track > &tv, int a0=5, int a1=4)
int usage (const std::string &name, int status)
template<typename T>
std::vector< T * > pointers (std::vector< T > &v)
bool SelectObjectET (const TrackTrigObject &t)
bool SelectObjectETovPT (const TrackTrigObject &tobj, TIDA::Track *t=0)
TrackFilter * getFilter (const std::string &refname, int pdgId, TrackFilter *foff, TrackFilter *fmu, TrackFilter *ftruth)
 This is awful code, passing in lots of filter pointers just so that they can be assigned neatly ?
bool GetRefTracks (const std::string &rc, const std::string &exclude, const std::vector< TIDA::Chain > &chains, NtupleTrackSelector &refTracks, TrackAssociator *ex_matcher, TrigObjectMatcher &tom)
int main (int argc, char **argv)

Variables

bool PRINT_BRESIDUALS
 stack trace headers
BinConfig g_binConfig
BinConfig electronBinConfig
BinConfig muonBinConfig
BinConfig tauBinConfig
BinConfig bjetBinConfig
BinConfig cosmicBinConfig
bool debugPrintout = false
double ETmin = 0
double ETovPTmin = 0
 this is a swiss knife function - by default if ET/PT > 0 such that fabs(ET/PT) > 0 is always true and we can always call this function - PT can never be 0 for a particle in the detector so can use this for both ET/PT selection and raw ET selection

Detailed Description

Author
mark sutton
Date
Fri 11 Jan 2019 07:41:26 CET

Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration

Definition in file rmain.cxx.

Function Documentation

◆ atoi_check()

int atoi_check ( const std::string & s)

convert string to integer with check if successful

Definition at line 129 of file rmain.cxx.

129 {
130 int i = std::atoi( s.c_str() );
131 char check[128];
132 std::sprintf( check, "%d", i );
133 if ( std::string(check)!=s ) return -1;
134 else return i;
135}

◆ copyReleaseInfo()

void copyReleaseInfo ( TTree * tree,
TFile * foutdir )

◆ GetChainAuthor()

int GetChainAuthor ( std::string chainName)

Return selected author, expects format L2_e20_medium:TrigL2SiTrackFinder:3 where 3 is the track author NB: this isn't actually needed at all.

can parse the entire name properly with this technique - NB: should move to some meaningful name for strategies "A", "B", etc

chops tokens off the front, up to the ":"

Definition at line 144 of file rmain.cxx.

144 {
145
148 std::string s1 = chop(chainName,":");
149 std::string s2 = chop(chainName,":");
150 std::string s3 = chop(chainName,":");
151
152 // if only s3=="" then no third token, if s2=="", no second token etc
153 // return atoi(s3.c_str());
154
155 // std::cout << ">" << s1 << "< >" << s2 << "< >" << s3 << "<" << std::endl;
156 if ( s3=="0" ) return 0;
157 if ( s3=="1" ) return 1;
158 if ( s3=="2" ) return 2;
159 // }
160 return -1;
161
162}
chop
std::string chop(std::string &s1, const std::string &s2)
Definition hcg.cxx:161

◆ getFilter()

TrackFilter * getFilter ( const std::string & refname,
int pdgId,
TrackFilter * foff,
TrackFilter * fmu,
TrackFilter * ftruth )

This is awful code, passing in lots of filter pointers just so that they can be assigned neatly ?

This section of the code needs a bit of a rethink getFilter( refname, &filter_off, &filter_muon, &filter_truth ); Fixme: maybe use a switch statement or something in the future

Definition at line 404 of file rmain.cxx.

407 {
408
409 std::cout << "getFilter(): refname " << refname << std::endl;
410
411 if ( refname=="Offline" ) return foff;
412 else if ( refname=="InDetLargeD0TrackParticles" ) return foff;
413 else if ( contains( refname, "Electrons") ) return foff;
414 else if ( contains( refname, "LRTElectrons") ) return foff;
415 else if ( contains( refname, "Muons" ) ) return fmu;
416 else if ( contains( refname, "MuonsLRT" ) ) return fmu;
417 else if ( contains( refname, "Taus" ) ) return foff; // tau ref chains
418 else if ( contains( refname, "1Prong" ) ) return foff; // tau ref chains
419 else if ( contains( refname, "3Prong" ) ) return foff; // tau ref chains
420 else if ( refname=="Truth" && pdgId!=0 ) return ftruth;
421 else if ( refname=="Truth" && pdgId==0 ) return foff;
422 else {
423 std::cerr << "unknown reference chain defined" << std::endl;
424 return 0;
425 }
426}
contains
bool contains(const std::string &s, const std::string &regx)
does a string contain the substring
Definition hcg.cxx:114

◆ GetRefTracks()

bool GetRefTracks ( const std::string & rc,
const std::string & exclude,
const std::vector< TIDA::Chain > & chains,
NtupleTrackSelector & refTracks,
TrackAssociator * ex_matcher,
TrigObjectMatcher & tom )

get any objects to exclude matched from the reference selection

match between the reference tracks and the objects to be excluded, and rebuild the reference tracks without the matches

if no match then add back to the standard reefrence

get objects if requested

Definition at line 429 of file rmain.cxx.

433 {
434
435 bool foundReference = false;
436
437 for ( size_t ic=chains.size() ; ic-- ; ) {
438
439 if ( chains[ic].name()==rc ) {
440
441 foundReference = true;
442
443 //Get tracks from within reference roi
444
446
447 if ( exclude.empty() ) {
448 refTracks.selectTracks( chains[ic].rois()[0].tracks() );
449 }
450 else {
451
452 std::vector<TIDA::Track> tmptracks = chains[ic][0].tracks();
453
454 for ( size_t ix=chains.size() ; ix-- ; ) {
455
456 if ( chains[ix].name()==exclude ) {
457
458 std::vector<TIDA::Track> extracks = chains[ix][0].tracks();
459
460 std::vector<TIDA::Track*> refp;
461 std::vector<TIDA::Track*> refx;
462
463 for ( size_t it=tmptracks.size() ; it-- ; ) refp.push_back( &tmptracks[it] );
464 for ( size_t it=extracks.size() ; it-- ; ) refx.push_back( &extracks[it] );
465
468
469 ex_matcher->match( refp, refx );
470
471 std::vector<TIDA::Track*> refp_ex;
472
473 for ( size_t it=refp.size() ; it-- ; ) {
475 if ( ex_matcher->matched(refp[it])==0 ) refp_ex.push_back(refp[it]);
476 }
477
478 refTracks.clear();
479 refTracks.selectTracks( refp_ex );
480
481 if ( debugPrintout ) {
482
483 std::cout << "\nexclude: " << refp.size() << "\t" << refp_ex.size() << std::endl;
484 std::cout << "exclude:\n" << extracks << std::endl;
485 std::cout << "reference tracks: " << std::endl;
486
487 size_t it0 = refp_ex.size();
488
489 for ( size_t it=0 ; it<refp.size() && it0>0 ; it++ ) {
490 if ( refp[it]==refp_ex[it0-1] ) std::cout << it << "\t" << *refp[it] << "\t" << *refp_ex[--it0] << std::endl;
491 else std::cout << "\n" << it << "\t" << *refp[it] << "\t----\n" << std::endl;
492 }
493
494 }
495
496 break;
497 }
498
499 }
500 }
501
503
504 if ( chains[ic].rois()[0].objects().size()>0 ) {
505 tom = TrigObjectMatcher( &refTracks, chains[ic].rois()[0].objects(), SelectObjectETovPT );
506 }
507
508 break;
509 }
510 }
511
512 return foundReference;
513
514}
rc
static Double_t rc
Definition LArPhysWaveHECTool.cxx:37
NtupleTrackSelector::clear
virtual void clear() override
the ntple selector manages the tracks itself, so we have to do an explicit delete for each one to pre...
Definition NtupleTrackSelector.h:71
NtupleTrackSelector::selectTracks
void selectTracks(const std::vector< TIDA::Track > &tracks)
extract all the tracks from a vector of Tracks
Definition NtupleTrackSelector.h:48
TIDA::Associator::matched
virtual const S * matched(T *t)
Definition TIDAAssociator.h:45
TIDA::Associator::match
virtual void match(const std::vector< T * > &s1, const std::vector< S * > &s2)=0
exclude
std::set< std::string > exclude
list of directories to be excluded
Definition hcg.cxx:98
grepfile.ic
int ic
Definition grepfile.py:33
SelectObjectETovPT
bool SelectObjectETovPT(const TrackTrigObject &tobj, TIDA::Track *t=0)
Definition rmain.cxx:388
debugPrintout
bool debugPrintout
Definition rmain.cxx:95

◆ handler()

void handler ( int sig)

signal handler

Definition at line 99 of file rmain.cxx.

99 {
100 void *array[10];
101
102 // get void*'s for all entries on the stack
103 size_t size = backtrace(array, 10);
104
105 // print out all the frames to stderr
106 std::cout << "Error: signal %d:\n" << sig << std::endl;
107 backtrace_symbols_fd(array, size, STDERR_FILENO);
108 std::exit(1);
109}
check_log.backtrace
list backtrace
Definition check_log.py:58

◆ ibl_filter()

const std::vector< TIDA::Track > ibl_filter ( const std::vector< TIDA::Track > & tv)

Definition at line 223 of file rmain.cxx.

223 {
224
225 static TH2D* h = 0;
226 static TH2D* h2 = 0;
227
228 static int ic = 0;
229
230
231 if ( h==0 ) {
232 h = new TH2D( "hzvphi", "hzvphi", 150, -300, 300, 150, -M_PI, M_PI );
233 h2 = new TH2D( "hzvphi2", "hzvphi", 150, -300, 300, 150, -M_PI, M_PI );
234 }
235
236 for ( size_t i=tv.size() ; i-- ; ) {
237
238 if ( tv[i].author()!=5 ) break;
239
240 double eta = tv[i].eta();
241
242 double theta = 2*std::atan(-std::exp(eta));
243
244 double ribl = 34.2; // actually 32.26 - 36.21 mm
245
246 double z = tv[i].z0() + ribl/std::tan(theta);
247
248 if ( !tv[i].expectBL() ) {
249 std::cout << "missing IBL: phi: " << tv[i].phi() << "\tz: " << z << " (" << eta << " " << theta*180/M_PI << ")" << std::endl;
250 if ( h ) h->Fill( z, tv[i].phi() );
251 }
252 else {
253 if ( h2 ) h2->Fill( z, tv[i].phi() );
254 }
255 }
256
257 if ( ic>=500 ) {
258 ic = 0;
259 if ( h ) {
260 h->DrawCopy();
261 gPad->Print("zphimap.pdf");
262 h2->DrawCopy();
263 gPad->Print("zphimap2.pdf");
264 }
265 }
266
267 ic++;
268
269 return tv;
270}
M_PI
#define M_PI
Definition ActiveFraction.h:11
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
theta
Scalar theta() const
theta method
Definition AmgMatrixBasePlugin.h:75
z
#define z
h
Header file for AthHistogramAlgorithm.
CaloCondBlobAlgs_fillNoiseFromASCII.author
str author
Definition CaloCondBlobAlgs_fillNoiseFromASCII.py:25

◆ main()

int main ( int argc,
char ** argv )

set some configuration parameters that can be set from the command line - these over-ride those from the configuration file Fixme: which is best? set them from the command line first, and then ignore them in the config file, or set from the config file, and then override them with command line arguments later? Theres only one way to find out...

get a filename from the command line if present

unknown option

true by default - command line options sets to false so should override what is in the config file

open output file

set up the filters etc

essentially no limit

essentially no limit

here we set a pTMax value less than pT, then only set the max pT in the filter if we read a pTMax value greater than pT

only if not set from the command line

only if not set from the command line

get the test chains

NB: don't override command line args

new code - can extract vtx name, pt, any extra options that we want, but also chop off everythiung after :post

now any additional config parameters for the chain are available

print soime debugging output

select only tracks within the roi?

read the (xml?) GRL

else get the list from the dat file directly

reference vertex selection

test vertex selection

option to use updated vertex matching with tracks

is this option needed any longer ???

read the binning config from a separate file if required

set the tags in front of the histogram stuff

clean up

track filters

truth articles

reminder: this is the new constructor - it now has too many parameters

filters for true selection for efficiency

include chi2 probability cut

track selectors so we can select multiple times with different filters if we want (simpler then looping over vectors each time

strings for which chains to analyse

create a map of the name to analysis string for selecting the correct analysis from the chain name

create the analyses and initialise them - this will create a root directory and book all the histograms

try not exiting if ref chain is different from default, temporarily ...

now we know that there is a corresponding tag, we can add the probe ...

if matching tag found then initialise tag and probe object and store tag and probe chains in there this will be passed into the ConfAnalysis, which will delete it when necessary
could perhaps be done with a unique_ptr

useful debug - leave this here ...

don't use the vertex name as the directory, since then we cannot easily plot on the same plot

hooray !! Set any additional aspects to the analysis up here so should be able to set individual, chain related variables eg set individual pt limits, different selection of reference objects, whether to run a vertex analysis ...

When setting different pt cuts, then an additional track selector will need to be created, this can then be added to the TIDA::FeatireStore of the analysis, and retrieved each time it is required

still needed for the moment ...

track selectors for efficiencies

do we want to filter the RoIs ?

Determine what sort of matching is required ...

default to deltaR matcher track matcher for best fit deltaR matcher

extra matcher for additionally matching reference to truth

default to deltaR matcher track matcher for best fit deltaR matcher

track selectors for purities

copy the release data to the output file

so we can specify the number of entries

we like, rather than run on all of them

check whether in good lumi block range

check whether it's in the event selector list

get the reference tracks (or perhaps even the true tracks !!!

I don't prefer range based for loops ...

select the reference offline vertices

always push back the vector - if required there will be only one vertex on it

calculate number of "vertex tracks"

do we really want this cut ???

get the tracks from the reference chain - shouldn't be one of the test chains, since there will be a 1-to-1 mapping but might be a useful check

get reference tracks and the revference TrigObjectMatcher if one is found in the reference chain

leave this in for the moment ...

reset the tom for this event

if expecting an Electron or a Tau and we don't have the TrigObjectMatcher then we don't have any actual offline electrons or taus so skip this event - code needs to be expanded a bit before being used ...

find the analysis for this chain - is there a matching analysis?

if no matching analysis then continue

these are the rois to process

both legs have the same reference ...

different references for each leg ...

tag leg reference ...

debug printout

get the rois and filter on them if required

trigger tracks already restricted by roi - so no roi filtering required

do we want to filter on the RoI properties? If so, if the RoI fails the cuts, then skip this roi

select the test sample (trigger) vertices

here we set the roi for the filter so we can request only those tracks inside the roi

what is all this logic doing ? It looks needlessly convoluted, and should at the very least have some proper explanation of what it is supposed to be doing

get the truth particles ...

truth_matcher match against current reference selection

which truth tracks have a matching reference track ?

check configuration is provided ...

remove any tracks below the pt threshold if one is specifed for the analysis

remove any tracks outside the d0 region if one is specifed for the analysis

if requesting an object match, remove any tracks which correspond to an object below the object PT threshold

only bother is objects actual exists

now found all the tracks within the RoI - now if required find the the count how many of these reference tracks are on each of the offline vertices

do for all vertices now ...

select the reference offline vertices

don't add vertices with no matching tracks - remember the tracks are filtered by Roi already so some vertices may have no tracks in the Roi - ntracks set to 0 by default

NB: because ntracks = 0 by default, this second clause should always be true unless ntracks has been set to some value

AAAAAARGH!!! because you cannot cast vector<T> to const vector<const T> we first need to copy the actual elements from the const vector, to a normal vector. This is because if we take the address of elements of a const vector<T> they will by of type const T*, so we can only add them to a vector<const T*> and all our functions are using const vector<T>&, so we would need to duplicate all the functions to allow over riding with vector<T*> and vector<const T*> to get this nonsense to work

so we now use a handy wrapper function to do the conversion for us ...

run the analysis for this chain

???

Definition at line 521 of file rmain.cxx.

522{
523 signal( SIGABRT, handler );
524 signal( SIGFPE, handler );
525 signal( SIGILL, handler );
526 signal( SIGSEGV, handler );
527 signal( SIGTERM, handler );
528
529 // ROOT::Cintex::Cintex::Enable();
530
531 if ( argc<2 ) {
532 std::cerr << "Error: no config file specified\n" << std::endl;
533 return usage(argv[0], -1);
534 }
535
536 // event_list elist("events.txt");
537
545
546 std::cout << "$0 :: compiled " << __DATE__ << " " << __TIME__ << std::endl;
547
549
550 std::string histofilename("");
551
552 std::string datafile = "";
553
554 std::vector<std::string> refChains(0);
555
556 std::vector<TrackFilter*> refFilters(0);
557
558 TrackFilter* refFilter = 0;
559 TrackFilter* truthFilter = 0;
560
561 std::string refChain = "";
562
563 int pdgId = 0;
564
565 std::vector<std::string> testChains;
566
567 std::string binningConfigFile = "";
568
569 bool useoldrms = true;
570 bool nofit = false;
571
572 std::string vertexSelection = "";
573 std::string vertexSelection_rec = "";
574
575 unsigned Nentries = 0;
576
577 for ( int i=1 ; i<argc ; i++ ) {
578 if ( std::string(argv[i])=="-h" || std::string(argv[i])=="--help" ) {
579 return usage(argv[0], 0);
580 }
581 else if ( std::string(argv[i])=="-e" || std::string(argv[i])=="--entries" ) {
582 if ( ++i>=argc ) return usage(argv[0], -1);
583 Nentries = std::atoi(argv[i]);
584 }
585 else if ( std::string(argv[i])=="-o" || std::string(argv[i])=="-f" || std::string(argv[i])=="--file" ) {
586 if ( ++i>=argc ) return usage(argv[0], -1);
587 histofilename = argv[i];
588 if ( histofilename.find(".root")==std::string::npos ) histofilename += ".root";
589 }
590 else if ( std::string(argv[i])=="-r" || std::string(argv[i])=="--refChain" ) {
591 if ( ++i>=argc ) return usage(argv[0], -1);
592 refChain = argv[i];
593
594 // Merge multiple references
595 if (refChain.find("+") != string::npos){
596 std::istringstream iss(refChain);
597 std::string token;
598 while (std::getline(iss, token, '+')){ // tokenize string based on '+' delimeter
599 refChains.push_back(token);
600 refFilters.push_back(0);
601 }
602 }
603 else {
604 refChains.push_back(argv[i]); // standard single reference
605 refFilters.push_back(0);
606 }
607 }
608 else if ( std::string(argv[i])=="--rms" ) useoldrms = false;
609 else if ( std::string(argv[i])=="-n" || std::string(argv[i])=="--nofit" ) nofit = true;
610 else if ( std::string(argv[i])=="-t" || std::string(argv[i])=="--testChain" ) {
611 if ( ++i>=argc ) return usage(argv[0], -1);
612 testChains.push_back(argv[i]);
613 }
614 else if ( std::string(argv[i])=="-p" || std::string(argv[i])=="--pdgId" ) {
615 if ( ++i>=argc ) return usage(argv[0], -1);
616 pdgId = atoi(argv[i]);
617 }
618 else if ( std::string(argv[i])=="--vr" ) {
619 if ( ++i>=argc ) return usage(argv[0], -1);
620 vertexSelection = argv[i];
621 }
622 else if ( std::string(argv[i])=="--vt" ) {
623 if ( ++i>=argc ) return usage(argv[0], -1);
624 vertexSelection_rec = argv[i];
625 }
626 else if ( std::string(argv[i])=="-b" || std::string(argv[i])=="--binConfig" ) {
627 if ( ++i>=argc ) return usage(argv[0], -1);
628 binningConfigFile = std::string(argv[i]);
629 }
630 else if ( std::string(argv[i]).find('-')==0 ) {
632 std::cerr << "unknown option " << argv[i] << std::endl;
633 return usage(argv[0], -1);
634 }
635 else {
636 datafile = argv[i];
637 }
638 }
639
640
641
642 if ( datafile=="" ) {
643 std::cerr << "no config file specifed\n" << endl;
644 return usage(argv[0], -1);
645 }
646
647 std::cout << time_str() << std::endl;
648
649 ReadCards inputdata(datafile);
650
651 inputdata.print();
652
655 if ( useoldrms ) {
656 bool oldrms95 = true;
657 inputdata.declareProperty( "OldRMS95", oldrms95 );
658 std::cout << "setting Resplot old rms95 " << oldrms95 << std::endl;
659 Resplot::setoldrms95( oldrms95 );
660 }
661 else {
662 std::cout << "setting Resplot old rms95 " << useoldrms << std::endl;
663 Resplot::setoldrms95( useoldrms );
664 }
665
666
667 if ( nofit ) {
668 std::cout << "Not fitting resplots " << std::endl;
669 Resplot::setnofit( nofit );
670 }
671
672
673 unsigned nfiles = 0;
674 inputdata.declareProperty( "NFiles", nfiles );
675
676
677 if ( histofilename=="" ){
678 if ( inputdata.isTagDefined("outputFile") ) histofilename = inputdata.GetString("outputFile");
679 else {
680 std::cerr << "Error: no output file defined\n" << std::endl;
681 return usage(argv[0], -1);
682 }
683 }
684
686 TFile foutput( histofilename.c_str(), "recreate" );
687 if (!foutput.IsOpen()) {
688 std::cerr << "Error: could not open output file\n" << std::endl;
689 return -1;
690 }
691
692 TDirectory* foutdir = gDirectory;
693
694 std::cout << "writing output to " << histofilename << std::endl;
695
696 TH1D* hevent = new TH1D( "event", "event", 1000, 10000, 80000 );
697 Resplot* hcorr = new Resplot( "correlation", 21, -0.5, 20.5, 75, 0, 600);
698
700
701 double pT = 1000;
702 double eta = 2.5;
703 double zed = 2000;
704
705 double a0min=0.;
706
707 int npix = 1;
708 int nsct = 6;
709 int nbl = -1;
710
711 int nsiholes = 2;
712 int npixholes = 20;
713 int nsctholes = 20;
714
715 bool expectBL = false;
716
717 double chi2prob = 0;
718
719 int npix_rec = -2;
720 int nsct_rec = -2;
721
722 double pT_rec = 0;
723 double eta_rec = 5;
724
725 double Rmatch = 0.1;
726
727 int ntracks = 0;
728
729 double massMax = 130; // default invariant mass window for Tag and Probe Analyses only
730 double massMin = 50;
731
732 //bool printflag = false; // JK removed (unused)
733
734 bool rotate_testtracks = false;
735
736 if ( inputdata.isTagDefined("RotateTestTracks") ) rotate_testtracks = ( inputdata.GetValue("RotateTestTracks") ? true : false );
737
738 bool truthMatch = false;
739
740 if ( inputdata.isTagDefined("TruthMatch") ) truthMatch = ( inputdata.GetValue("TruthMatch") ? true : false );
741
742 // CK: Option to specify which ref tracks to use from ref track vector, ordered by pT
743 // User parameter is a vector of ref track vector indices
744 // e.g. specifyPtOrderedRefTracks = {0, 1, 5} will use the first leading, second leading, and sixth leading track
745 std::vector<size_t> refPtOrd_indices;
746 bool use_pt_ordered_ref = false;
747
748 if ( inputdata.isTagDefined("UsePtOrderedRefTracks") ) {
749 std::vector<double> refPtOrd_indices_tmp;
750
751 use_pt_ordered_ref = true;
752 refPtOrd_indices_tmp = ( inputdata.GetVector("UsePtOrderedRefTracks") );
753
754 std::cout << "using PT ordered reference tracks: " << refPtOrd_indices_tmp << std::endl;
755
756 for ( size_t i=refPtOrd_indices_tmp.size(); i-- ; ) {
757 refPtOrd_indices.push_back( refPtOrd_indices_tmp.at(i) );
758 }
759 }
760
761 if ( inputdata.isTagDefined("pT") ) pT = inputdata.GetValue("pT");
762 if ( inputdata.isTagDefined("ET") ) ETmin = inputdata.GetValue("ET");
763 if ( inputdata.isTagDefined("ETovPT") ) ETovPTmin = inputdata.GetValue("ETovPT");
764
767 double pTMax = pT-1;
768 if ( inputdata.isTagDefined("pTMax") ) pTMax = inputdata.GetValue("pTMax");
769
770
771 if ( inputdata.isTagDefined("NMod") ) NMod = inputdata.GetValue("NMod");
772
773
774 if ( inputdata.isTagDefined("eta") ) eta = inputdata.GetValue("eta");
775 if ( inputdata.isTagDefined("zed") ) zed = inputdata.GetValue("zed");
776 else if ( inputdata.isTagDefined("z0") ) zed = inputdata.GetValue("z0");
777 if ( inputdata.isTagDefined("npix") ) npix = inputdata.GetValue("npix");
778 if ( inputdata.isTagDefined("nsiholes") ) nsiholes = inputdata.GetValue("nsiholes");
779 if ( inputdata.isTagDefined("npixholes") ) npixholes = inputdata.GetValue("npixholes");
780 if ( inputdata.isTagDefined("nsctholes") ) nsctholes = inputdata.GetValue("nsctholes");
781 if ( inputdata.isTagDefined("expectBL") ) expectBL = ( inputdata.GetValue("expectBL") > 0.5 ? true : false );
782 if ( inputdata.isTagDefined("nsct") ) nsct = inputdata.GetValue("nsct");
783 if ( inputdata.isTagDefined("nbl") ) nbl = inputdata.GetValue("nbl");
784 if ( inputdata.isTagDefined("chi2prob") ) chi2prob = inputdata.GetValue("chi2prob");
785
786 if ( inputdata.isTagDefined("ntracks") ) ntracks = inputdata.GetValue("ntracks")+0.5; // rounding necessary ?
787
788
789 if ( inputdata.isTagDefined("chi2prob") ) chi2prob = inputdata.GetValue("chi2prob");
790
792 if ( pdgId==0 && inputdata.isTagDefined("pdgId") ) pdgId = inputdata.GetValue("pdgId");
793
794 if ( inputdata.isTagDefined("InvMassMax") ) massMax = inputdata.GetValue("InvMassMax");
795 if ( inputdata.isTagDefined("InvMassMin") ) massMin = inputdata.GetValue("InvMassMin");
796
797 if ( inputdata.isTagDefined("npix_rec") ) npix_rec = inputdata.GetValue("npix_rec");
798 if ( inputdata.isTagDefined("nsct_rec") ) nsct_rec = inputdata.GetValue("nsct_rec");
799
800 if ( inputdata.isTagDefined("pT_rec") ) pT_rec = inputdata.GetValue("pT_rec");
801 if ( inputdata.isTagDefined("eta_rec") ) eta_rec = inputdata.GetValue("eta_rec");
802
803 if ( inputdata.isTagDefined("a0") ) a0 = inputdata.GetValue("a0");
804 if ( inputdata.isTagDefined("a0min") ) a0min = inputdata.GetValue("a0min");
805
806 if ( inputdata.isTagDefined("Rmatch") ) Rmatch = inputdata.GetValue("Rmatch");
807
808 std::string useMatcher = "DeltaR";
809 if ( inputdata.isTagDefined("UseMatcher") ) useMatcher = inputdata.GetString("UseMatcher");
810
811
813 if ( refChain=="" ) {
814 if ( inputdata.isTagDefined("refChain") ) {
815 refChain = inputdata.GetString("refChain");
816 refChains.push_back(refChain);
817 refFilters.push_back(0);
818 }
819 else {
820 std::cerr << "Error: no reference chain defined\n" << std::endl;
821 // return usage(argv[0], -1);
822 return -1;
823 }
824 }
825
826 std::string exclude = "";
827
828 TrackAssociator* ex_matcher = 0;
829
830 if ( inputdata.isTagDefined("Exclude") ) {
831 exclude = inputdata.GetString("Exclude");
832 ex_matcher = new Associator_BestDeltaRZSinThetaMatcher( "deltaRZ", 0.01, 0.01, 2 );
833 }
834
835 if (refChains.size() == 0){
836 std::cerr << "Error: refChains is empty\n" <<std::endl;
837 return -1;
838 }
839
841 if ( testChains.size()==0 ) {
842 if ( inputdata.isTagDefined("testChains") ) testChains = inputdata.GetStringVector("testChains");
843 else if ( inputdata.isTagDefined("testChain") ) testChains.push_back( inputdata.GetString("testChain") );
844 }
845
848
849 std::vector<ChainString> chainConfig;
850
851 for ( size_t ic=0 ; ic<testChains.size() ; ic++ ) {
852 chainConfig.push_back( ChainString( testChains[ic] ) );
853 testChains[ic] = chainConfig.back().pre();
854 }
855
856
858
859
861 //if ( inputdata.isTagDefined("printflag") ) printflag = ( inputdata.GetValue("printflag") ? 1 : 0 ); // JK removed (unused)
862
864 bool select_roi = true;
865
866
867 // CK: use a custom filter RoI for ref tracks using refFilter?
868 // Note that if the reference chain uses a non-full-scan RoI, the custom filter RoI could be wider than the
869 // ref chain RoI
870 bool use_custom_ref_roi = false;
871 const int custRefRoi_nParams = 3;
872
873 double custRefRoi_params[custRefRoi_nParams] = {-999., -999., -999.};
874 std::vector<std::string> custRefRoi_chainList; // Vector of chain names to apply the custom roi to
875 std::set<std::string> customRoi_chains;
876
877 if ( inputdata.isTagDefined("customRefRoi_etaHalfWidth") ) custRefRoi_params[0] = inputdata.GetValue("customRefRoi_etaHalfWidth");
878 if ( inputdata.isTagDefined("customRefRoi_phiHalfWidth") ) custRefRoi_params[1] = inputdata.GetValue("customRefRoi_phiHalfWidth");
879 if ( inputdata.isTagDefined("customRefRoi_zedHalfWidth") ) custRefRoi_params[2] = inputdata.GetValue("customRefRoi_zedHalfWidth");
880
881 if ( inputdata.isTagDefined("customRefRoi_chainList") ) custRefRoi_chainList = inputdata.GetStringVector("customRefRoi_chainList");
882
883 for ( unsigned ic=0 ; ic<custRefRoi_chainList.size() ; ic++ ) customRoi_chains.insert( custRefRoi_chainList[ic] );
884
885 for ( int param_idx=0 ; param_idx<custRefRoi_nParams ; param_idx++ ) {
886 if ( custRefRoi_params[param_idx] != -999 ) {
887 select_roi = true; // In case select_roi is ever set to default to false
888 use_custom_ref_roi = true;
889 }
890 }
891
892 if ( use_custom_ref_roi ) {
893 std::cout << "**** \t****" << std::endl;
894 std::cout << "**** Custom RoI will be used to filter ref. tracks\t****" << std::endl;
895
896 if ( custRefRoi_params[0] != -999. ) std::cout << "**** etaHalfWidth = " << custRefRoi_params[0] << "\t\t\t\t****" << std::endl;
897 else std::cout << "**** etaHalfWidth = value used in trigger RoI\t****" << std::endl;
898
899 if ( custRefRoi_params[1] != -999. ) std::cout << "**** phiHalfWidth = " << custRefRoi_params[1] << "\t\t\t\t****" << std::endl;
900 else std::cout << "**** phiHalfWidth = value used in trigger RoI\t****" << std::endl;
901
902 if ( custRefRoi_params[2] != -999. ) std::cout << "**** zedHalfWidth = " << custRefRoi_params[2] << "\t\t\t\t****" << std::endl;
903 else std::cout << "**** zedHalfWidth = value used in trigger RoI\t****" << std::endl;
904
905 if ( !custRefRoi_chainList.empty() ) {
906 std::cout << "**** \t****" << std::endl;
907 std::cout << "**** Applying custom RoI only to specified chains\t****" << std::endl;
908 }
909 std::cout << "**** \t****" << std::endl;
910 }
911
912 // Checking for SelectRoi after any other options that will set select_roi, to ensure that the value set
913 // here will be used
914 if ( inputdata.isTagDefined("SelectRoi") ) {
915 select_roi = ( inputdata.GetValue("SelectRoi")!=0 ? true : false );
916 }
917
918 if ( !select_roi ) {
919 std::cout << "**** ****" << std::endl;
920 std::cout << "**** RoI filtering of reference tracks is disabled ****" << std::endl;
921 std::cout << "**** ****" << std::endl;
922 }
923
924 // bool selectfake_roi = false; // JK removed (unused)
925 // if ( inputdata.isTagDefined("SelectFakeRoi") ) {
926 // selectfake_roi = ( inputdata.GetValue("SelectFakeRoi")!=0 ? true : false ); // JK removed (unused)
927 // }
928
929
930 std::vector<double> lumiblocks;
931 lumiParser goodrunslist;
932
933
934 if ( inputdata.isTagDefined("GRL") ) {
936 std::vector<std::string> grlvector = inputdata.GetStringVector("GRL");
937 std::cout << "Reading GRL from: " << grlvector << std::endl;
938 for ( size_t igrl=0 ; igrl<grlvector.size() ; igrl++ ) goodrunslist.read( grlvector[igrl] );
939 // std::cout << goodrunslist << std::endl;
940 }
941 else if ( inputdata.isTagDefined("LumiBlocks") ) {
943 lumiblocks = inputdata.GetVector("LumiBlocks");
944
945 for (unsigned int i=0 ; i<lumiblocks.size()-2 ; i+=3 ){
946 goodrunslist.addRange( lumiblocks[i], lumiblocks[i+1], lumiblocks[i+2] );
947 }
948 }
949
950
952
953
954 if ( vertexSelection == "" ) {
955 if ( inputdata.isTagDefined("VertexSelection") ) vertexSelection = inputdata.GetString("VertexSelection");
956 }
957
958 bool bestPTVtx = false;
959 bool bestPT2Vtx = false;
960 int vtxind = -1;
961
962 if ( vertexSelection!="" ) {
963 if ( vertexSelection=="BestPT" ) bestPTVtx = true;
964 else if ( vertexSelection=="BestPT2" ) bestPT2Vtx = true;
965 else vtxind = atoi_check( vertexSelection );
966 }
967
968
969
970
972
973 if ( vertexSelection_rec == "" ) {
974 if ( inputdata.isTagDefined("VertexSelectionRec") ) vertexSelection_rec = inputdata.GetString("VertexSelectionRec");
975 }
976
977 bool bestPTVtx_rec = false;
978 bool bestPT2Vtx_rec = false;
979 int vtxind_rec = -1;
980
981 if ( vertexSelection_rec!="" ) {
982 if ( vertexSelection_rec=="BestPT" ) bestPTVtx_rec = true;
983 else if ( vertexSelection_rec=="BestPT2" ) bestPT2Vtx_rec = true;
984 else vtxind_rec = atoi_check( vertexSelection_rec );
985 }
986
987 std::cout << "vertexSelection: " << vertexSelection << std::endl;
988 std::cout << "vertexSelection_rec: " << vertexSelection_rec << std::endl;
989
990
991#if 0
992
995
996 bool useBestVertex = false;
997 if ( inputdata.isTagDefined("useBestVertex") ) useBestVertex = ( inputdata.GetValue("useBestVertex") ? 1 : 0 );
998
999 bool useSumPtVertex = true;
1000 if ( inputdata.isTagDefined("useSumPtVertex") ) useSumPtVertex = ( inputdata.GetValue("useSumPtVertex") ? 1 : 0 );
1001
1002 int MinVertices = 1;
1003 if ( inputdata.isTagDefined("MinVertices") ) MinVertices = inputdata.GetValue("MinVertices");
1004
1006
1007#endif
1008
1010 bool useVertexTracks = false;
1011 if ( inputdata.isTagDefined("UseVertexTracks") ) useVertexTracks = ( inputdata.GetValue("UseVertexTracks") > 0 );
1012
1013 // option to vary default "PrimaryVertices" offline reference collection
1014 std::string vertex_refname = "Vertex";
1015 if ( inputdata.isTagDefined("VertexReference") ) vertex_refname += ":" + inputdata.GetString("VertexReference");
1016
1018 int NVtxTrackCut = 2;
1019 if ( inputdata.isTagDefined("NVtxTrackCut") ) NVtxTrackCut = inputdata.GetValue("NVtxTrackCut");
1020
1021
1022 std::vector<double> event_list;
1023 bool event_selector_flag = false;
1024
1025 if ( inputdata.isTagDefined("EventSelector") ) event_list = inputdata.GetVector("EventSelector");
1026
1027 event_selector es(event_list);
1028
1029 if ( es.size() ) event_selector_flag = true;
1030
1031
1032 std::vector<double> beamTest;
1033 std::vector<double> beamRef;
1034
1035
1036 bool correctBeamlineRef = false;
1037 bool correctBeamlineTest = false;
1038
1039 if ( inputdata.isTagDefined("CorrectBeamlineRef") ) correctBeamlineRef = ( inputdata.GetValue("CorrectBeamlineRef") == 0 ? false : true );
1040 if ( inputdata.isTagDefined("CorrectBeamlineTest") ) correctBeamlineTest = ( inputdata.GetValue("CorrectBeamlineTest") == 0 ? false : true );
1041
1042
1043 if ( inputdata.isTagDefined("BeamTest") ) beamTest = inputdata.GetVector("BeamTest");
1044 else {
1045 if ( inputdata.isTagDefined("BeamTestx") ) beamTest.push_back(inputdata.GetValue("BeamTestx"));
1046 if ( inputdata.isTagDefined("BeamTesty") ) beamTest.push_back(inputdata.GetValue("BeamTesty"));
1047 }
1048
1049
1050 if ( inputdata.isTagDefined("BeamRef") ) beamRef = inputdata.GetVector("BeamRef");
1051 else {
1052 if ( inputdata.isTagDefined("BeamRefx") ) beamRef.push_back(inputdata.GetValue("BeamRefx"));
1053 if ( inputdata.isTagDefined("BeamRefy") ) beamRef.push_back(inputdata.GetValue("BeamRefy"));
1054 }
1055
1056
1057
1058 if ( ( beamTest.size()!=0 && beamTest.size()!=2 && beamTest.size()!=3 ) ||
1059 ( beamRef.size()!=0 && beamRef.size()!=2 && beamRef.size()!=3 ) ) {
1060 std::cerr << "incorrectly specified beamline position" << std::endl;
1061 return (-1);
1062 }
1063
1064 if ( beamTest.size()>0 ) correctBeamlineTest = true;
1065 if ( beamRef.size()>0 ) correctBeamlineRef = true;
1066
1067 if ( correctBeamlineRef ) std::cout << "main() correcting beamline for reference tracks" << std::endl;
1068 if ( correctBeamlineTest ) std::cout << "main() correcting beamline for test tracks" << std::endl;
1069
1070
1071
1072 if ( beamRef.size()>0 ) std::cout << "beamref " << beamRef << std::endl;
1073 if ( beamTest.size()>0 ) std::cout << "beamtest " << beamTest << std::endl;
1074
1075 double a0v = 1000;
1076 double z0v = 2000;
1077
1078 if ( inputdata.isTagDefined("a0v") ) a0v = inputdata.GetValue("a0v");
1079 if ( inputdata.isTagDefined("z0v") ) z0v = inputdata.GetValue("z0v");
1080
1081
1082 double a0vrec = 1000;
1083 double z0vrec = 2000;
1084
1085 if ( inputdata.isTagDefined("a0vrec") ) a0vrec = inputdata.GetValue("a0vrec");
1086 if ( inputdata.isTagDefined("z0vrec") ) z0vrec = inputdata.GetValue("z0vrec");
1087
1088 bool initialiseFirstEvent = false;
1089 if ( inputdata.isTagDefined("InitialiseFirstEvent") ) initialiseFirstEvent = inputdata.GetValue("InitialiseFirstEvent");
1090
1091
1092 // set the flag to prinout the missing track list in ConfAnalysis
1093 if ( inputdata.isTagDefined("dumpflag") ) dumpflag = ( inputdata.GetValue("dumpflag")==0 ? false : true );
1094
1095
1096 bool doPurity = false;
1097 if ( inputdata.isTagDefined("doPurity") ) doPurity = ( inputdata.GetValue("doPurity")==0 ? false : true );
1098
1099 if ( inputdata.isTagDefined("DebugPrintout") ) debugPrintout = ( inputdata.GetValue("DebugPrintout")==0 ? false : true );
1100
1101
1102 bool monitorZBeam = false;
1103 if ( inputdata.isTagDefined("MonitorinZBeam") ) monitorZBeam = ( inputdata.GetValue("MonitorZBeam")==0 ? false : true );
1104
1105 std::cout << "dbg " << __LINE__ << std::endl;
1106
1107 ReadCards* binningConfig = &inputdata;
1108
1110
1111 if ( binningConfigFile!="" ) binningConfig = new ReadCards( binningConfigFile );
1112
1114 g_binConfig.set( *binningConfig, "" );
1115 electronBinConfig.set( *binningConfig, "e_" );
1116 muonBinConfig.set( *binningConfig, "mu_" );
1117 tauBinConfig.set( *binningConfig, "tau_" );
1118 bjetBinConfig.set( *binningConfig, "bjet_" );
1119 cosmicBinConfig.set( *binningConfig, "cosmic_" );
1120
1121
1123 if ( binningConfig!=&inputdata ) delete binningConfig; // cppcheck-suppress autovarInvalidDeallocation; false positive
1124
1125
1126 if ( inputdata.isTagDefined("PRINT_BRESIDUALS") ) PRINT_BRESIDUALS = ( inputdata.GetValue("PRINT_BRESIDUALS")==0 ? false : true );
1127
1128 int selectcharge = 0;
1129 if ( inputdata.isTagDefined("Charge") ) selectcharge = inputdata.GetValue("Charge");
1130
1131
1132 std::cout << "using reference " << refChain << std::endl;
1133 if ( refChain.find("Truth") != string::npos ) std::cout << "using pdgId " << pdgId << std::endl;
1134 if ( refChains.size() > 1 ) std::cout<<"Multiple reference chains split to: " << refChains <<std::endl;
1135
1137
1139
1141 // Filter_Track( double etaMax, double d0Max, double z0Max, double pTMin,
1142 // int minPixelHits, int minSctHits, int minSiHits, int minBlayerHits,
1143 // int minStrawHits, int minTrHits, double prob=0 ) :
1144
1146
1147 std::cout << "a0v: " << a0v << std::endl;
1148 std::cout << "z0v: " << z0v << std::endl;
1149
1150 Filter_Vertex filter_vertex( a0v, z0v );
1151
1152 Filter_Track filter_offline( eta, 1000, a0min, zed, pT,
1153 npix, nsct, -1, nbl,
1154 -2, -2, chi2prob,
1155 npixholes, nsctholes, nsiholes, expectBL );
1156
1157 if ( selectcharge!=0 ) filter_offline.chargeSelection( selectcharge );
1158 if ( pTMax>pT ) filter_offline.maxpT( pTMax );
1159
1160 Filter_etaPT filter_etaPT(eta,pT);
1161
1162
1163 // std::cout << "pdgId " << pdgId << std::endl;
1164 //Select only true particles matching pdgId
1165 Filter_pdgId filter_pdgtruth(pdgId); //BP
1166
1167 // select inside-out offline tracking
1168
1169 // Filter_TruthParticle filter_passthrough(&filter_offline);
1170 // Filter_Track filter_track( eta, 1000, 2000, pT, npix, nsct, -1, -1, -2, -2);
1171
1172 Filter_Author filter_inout(0);
1173
1174 int author = 0;
1175
1176 Filter_Author filter_auth(author);
1177
1178 Filter_TrackQuality filter_q(0.01);
1179 Filter_Combined filter_off(&filter_offline, &filter_vertex);
1180
1181 Filter_Combined filter_truth( &filter_pdgtruth, &filter_etaPT);
1182
1183 Filter_Combined filter_muon( &filter_offline, &filter_vertex);
1184
1185 Filter_Track filter_onlinekine( eta_rec, 1000, 0., 2000, pT, -1, npix, nsct, -1, -2, -2);
1186 Filter_Vertex filter_onlinevertex( a0vrec, z0vrec);
1187 Filter_Combined filter_online( &filter_onlinekine, &filter_onlinevertex );
1188
1189 Filter_Track filter_offkinetight( 5, 1000, 0., 2000, pT, -1, 0, 0, -1, -2, -2);
1190 Filter_Combined filter_offtight( &filter_offkinetight, &filter_inout );
1191
1192 Filter_Offline2017* filter_offline2017 = 0;
1193 Filter_OfflineR22* filter_offlineR22 = 0;
1194 Filter_Combined* filter_offCP = 0;
1197
1198
1199 if ( inputdata.isTagDefined("Filter" ) ) {
1200 ChainString filter = inputdata.GetString("Filter");
1201 std::cout << "Filter: " << inputdata.GetString("Filter") << " : " << filter << std::endl;
1202 if ( filter.head()=="Offline2017" ) {
1203 std::string filter_type = filter.tail();
1204 filter_offline2017 = new Filter_Offline2017( pT, filter_type, zed, a0 );
1205 filter_offCP = new Filter_Combined ( filter_offline2017, &filter_vertex);
1206 refFilter = filter_offCP;
1207 }
1208 else if ( filter.head()=="OfflineR22" ) {
1209 std::string filter_type = filter.tail();
1210 filter_offlineR22 = new Filter_OfflineR22( pT, filter_type, zed, a0 );
1211 filter_offCP = new Filter_Combined ( filter_offlineR22, &filter_vertex);
1212 refFilter = filter_offCP;
1213 }
1214 else {
1215 std::cerr << "unimplemented Filter requested: " << filter.head() << std::endl;
1216 return -1;
1217 }
1218 }
1219 else {
1220 if ( !(refFilter = getFilter( refChains[0], pdgId, &filter_off, &filter_muon, &filter_truth ) ) ) {
1221 std::cerr << "unknown reference chain defined" << std::endl;
1222 return (-1);
1223 }
1224 }
1225
1226 refFilters.push_back(refFilter);
1227
1228
1229 std::map<std::string,TIDA::Reference> ref;
1230
1231 std::vector<NtupleTrackSelector*> refSelectors;
1232
1233#if 0
1237
1238 for ( size_t ic=0 ; ic<refChains.size() ; ic++ ) {
1239
1240 if ( refFilter==0 ) {
1241 if ( !(refFilters[ic] = getFilter( refChains[ic], pdgId, &filter_off, &filter_muon, &filter_truth ) ) ) {
1242 std::cerr << "unknown reference chain defined" << std::endl;
1243 return (-1);
1244 }
1245 refFilter = refFilters[ic];
1246 }
1247 else refFilters[ic] = refFilter;
1248
1249 ref.insert( TIDA::ReferenceMap::value_type( refChains[ic], TIDA::Reference( refChains[ic], new NtupleTrackSelector(refFilter), refFilter, new TrigObjectMatcher ) ) );
1250
1251 }
1252
1253# endif
1254
1255 if (pdgId==0) truthFilter = &filter_off;
1256 else truthFilter = &filter_truth;
1257
1258
1259 // use an actual filter requiring at least 1 silicon hit
1260 // to get rid of the EF trt only tracks
1261
1262 std::cout << "filter_passthrough" << std::endl;
1263
1264 Filter_Track filter_passthrough( 10, 1000, 0., 2000, pT_rec, npix_rec, nsct_rec, 1, -2, -2, -2);
1265
1266 TrackFilter* testFilter = &filter_passthrough;
1267
1268
1269 std::cout << "using tracks: " << refChain << " for reference sample" << std::endl;
1270
1271
1273 // std::vector<std::string> chains;
1274
1275 std::vector<std::string>& test_chains = testChains;
1276
1279 //smh: TrackAnalysis is purely abstract, analysis will contain ConfAnalysis
1280 std::map<std::string,TrackAnalysis*> analysis;
1281
1285
1286 std::vector<TrackAnalysis*> analyses;
1287 analyses.reserve(test_chains.size());
1288
1289 std::cout << "booking " << test_chains.size() << " analyses" << std::endl;
1290
1291 for ( unsigned i=0 ; i<test_chains.size() ; i++ ) {
1292
1293 std::string chainname = ChainString(test_chains[i]);
1294
1295 std::vector<std::string> chainnames;
1296
1297 chainnames.push_back(chainname);
1298
1299
1300 // tag and probe object creation and configuration
1301
1302 TagNProbe* TnP_tool = 0;
1303 ChainString probe = chainConfig[i];
1304
1305 std::string probe_extra = probe.extra();
1306
1307 if ( probe_extra.empty() ) probe_extra = probe.postvalue("extra");
1308
1309 if ( probe_extra.find("_tag")!=std::string::npos || probe.extra().find("_tag")!=std::string::npos ) {
1310 std::cout << "rejecting tag chain " << probe << std::endl;
1311 continue;
1312 }
1313
1314 // probe can be the .head() so convert m_chainNames to a ChainString and search the .extra() specifically
1315 size_t p = probe_extra.find("_probe");
1316
1317 if ( p!=std::string::npos ) {
1318
1319 std::string probe_ref = refChains[0];
1320
1321 if ( !probe.postvalue("ref").empty() ) {
1322
1323 probe_ref = probe.postvalue("ref");
1324
1325 if ( refChains[0] != probe_ref ) {
1326 std::cerr << "default and probe chain references do not match: probe ref: " << probe_ref << " ref: " << refChains[0] << std::endl;
1328 // return -1;
1329 }
1330
1331 }
1332
1333
1334 std::string probe_key = probe_extra.erase(p, 6);
1335
1336 for ( unsigned j=0 ; j<test_chains.size(); ++j) {
1337
1338 if ( i==j ) continue;
1339
1340 ChainString tag = chainConfig[j];
1341
1342 if ( tag.head() != probe.head() ) continue;
1343
1344 // if ( tag.tail() != probe.tail() ) continue; // should not enforce this for mu + tau tnp chains
1345
1346 std::string tag_extra = tag.extra();
1347
1348 if ( tag_extra.empty() ) tag_extra = tag.postvalue("extra");
1349 if ( tag_extra.find("_tag")==std::string::npos ) continue;
1350
1351 // need to compare just the 'el1' part of of .extra() so create variables without '_probe/_tag' part
1352 std::string tag_key = tag_extra.erase( tag_extra.find("_tag"), 4) ;
1353
1354 // tag chain must be the same as probe chain but with te=0 and extra=*_tag
1355 if ( tag_key != probe_key ) continue;
1356
1357 if ( tag.element() == probe.element() ) continue;
1358
1359 std::string tag_ref = refChains[0];
1360
1362
1363 TrackFilter* filter = getFilter( probe_ref, pdgId, &filter_off, &filter_muon, &filter_truth );
1364
1365 ref.insert( TIDA::ReferenceMap::value_type( probe_ref, TIDA::Reference( probe_ref, new NtupleTrackSelector(filter), filter, new TrigObjectMatcher ) ) );
1366
1367 if ( !tag.postvalue("ref").empty() ) {
1368
1369 tag_ref = tag.postvalue("ref");
1370 // refChains.push_back( tag_ref);
1371
1372 std::cout << "tag ref: " << tag_ref << std::endl;
1373
1374 if ( ref.find(tag_ref)==ref.end() ) {
1375 TrackFilter* filter = getFilter( tag_ref, pdgId, &filter_off, &filter_muon, &filter_truth );
1376 // ref.insert( std::map<std::string,TIDA::Reference>::value_type( tag_ref, TIDA::Reference( tag_ref, new NtupleTrackSelector(filter) ) ) );
1377 ref.insert( TIDA::ReferenceMap::value_type( tag_ref, TIDA::Reference( tag_ref, new NtupleTrackSelector(filter), filter, new TrigObjectMatcher ) ) );
1378 }
1379 }
1380
1384 // TnP_tool = new TagNProbe( refChains[0], refChains[0], massMin, massMax);
1385 TnP_tool = new TagNProbe( tag_ref, probe_ref, massMin, massMax);
1386 TnP_tool->tag(tag);
1387 TnP_tool->probe(probe);
1388 std::cout << "Tag and probe pair found! \n\t Tag : " << tag << "\n\t Probe : " << probe
1389 << "\n\t tag ref: " << tag_ref
1390 << "\n\tprobe ref: " << probe_ref
1391 << "\n-------------------" << std::endl;
1392
1394 // std::cout << *TnP_tool << std::endl;
1395
1396 break ;
1397 }
1398
1399 }
1400
1401 // Replace "/" with "_" in chain names, for Tag&Probe analysis
1402 if ( TnP_tool ) replace( chainname, "/", "_" );
1403 ConfAnalysis* analy_conf = new ConfAnalysis( chainname, chainConfig[i], TnP_tool );
1404
1405 analy_conf->initialiseFirstEvent(initialiseFirstEvent);
1406 analy_conf->initialise();
1407 analy_conf->setprint(false);
1408
1409 std::string vtxTool = chainConfig[i].postvalue("rvtx");
1410
1411 if ( vtxTool!="" ) {
1414 ConfVtxAnalysis* anal_confvtx = new ConfVtxAnalysis( vtxTool, (vertex_refname!="Vertex") );
1415 // ConfVtxAnalysis* anal_confvtx = new ConfVtxAnalysis( "vertex" );
1416 analy_conf->store().insert( anal_confvtx, "rvtx" );
1417 }
1418
1419
1420 // analy_conf->setprint(true);
1421
1426
1431
1432 if ( chainConfig[i].values().size()>0 ) {
1433 std::cout << "chain:: " << chainname << " : size " << chainConfig[i].values().size() << std::endl;
1434 for ( unsigned ik=chainConfig[i].values().size() ; ik-- ; ) {
1435 std::cout << "\tchainconfig: " << ik << "\tkey " << chainConfig[i].keys()[ik] << " " << chainConfig[i].values()[ik] << std::endl;
1436 }
1437 }
1438
1440 // for (unsigned int ic=0 ; ic<chainnames.size() ; ic++ ) analysis[chainnames[ic]] = analy_conf;
1441
1442 if ( analysis.find( chainname )==analysis.end() ) {
1443 analysis.insert( std::map<std::string,TrackAnalysis*>::value_type( chainname, analy_conf ) );
1444 analyses.push_back(analy_conf);
1445 }
1446 else {
1447 std::cerr << "WARNING: Duplicated chain"
1448 << "\n"
1449 << "---------------------------------"
1450 << "---------------------------------"
1451 << "---------------------------------" << std::endl;
1452 continue;
1453 }
1454
1455 std::cout << "analysis: " << chainname << "\t" << analy_conf
1456 << "\n"
1457 << "---------------------------------"
1458 << "---------------------------------"
1459 << "---------------------------------" << std::endl;
1460
1461 if ( doPurity ) {
1462 PurityAnalysis* analp = new PurityAnalysis(chainnames[0]+"-purity");
1463 std::cout << "purity " << (chainnames[0]+"-purity") << std::endl;
1464 analp->initialise();
1465 analp->setprint(false);
1466 // analp->setprint(true);
1467 analysis[chainnames[0]+"-purity"] = analp;
1468 analyses.push_back(analp);
1469 }
1470
1471 }
1472
1473 std::cout << "main() finished looping" << std::endl;
1474
1476
1477 bool fullyContainTracks = false;
1478
1479 if ( inputdata.isTagDefined("FullyContainTracks") ) {
1480 fullyContainTracks = ( inputdata.GetValue("FullyContainTracks")==0 ? false : true );
1481 }
1482
1483 bool containTracksPhi = true;
1484
1485 if ( inputdata.isTagDefined("ContainTracksPhi") ) {
1486 containTracksPhi = ( inputdata.GetValue("ContainTracksPhi")==0 ? false : true );
1487 }
1488
1489 if ( dynamic_cast<Filter_Combined*>(refFilter) ) {
1490 dynamic_cast<Filter_Combined*>(refFilter)->setDebug(debugPrintout);
1491 dynamic_cast<Filter_Combined*>(refFilter)->containtracks(fullyContainTracks);
1492 dynamic_cast<Filter_Combined*>(refFilter)->containtracksPhi(containTracksPhi);
1493 }
1494
1495 NtupleTrackSelector refTracks( refFilter );
1496 NtupleTrackSelector offTracks( testFilter );
1497 NtupleTrackSelector testTracks( testFilter);
1498
1499 NtupleTrackSelector truthTracks( truthFilter );
1500
1501 // NtupleTrackSelector refTracks( &filter_passthrough );
1502 // NtupleTrackSelector testTracks( refFilter );
1503 // NtupleTrackSelector refTracks( &filter_pdgtruth);
1504 // NtupleTrackSelector testTracks( &filter_off );
1505 // NtupleTrackSelector testTracks(&filter_roi);
1506
1508
1509 bool filterRoi = false;
1510
1511 double roieta = 0;
1512 bool roicomposite = false;
1513 size_t roimult = 1;
1514
1515 if ( inputdata.isTagDefined( "FilterRoi" ) ) {
1516
1517 filterRoi = true;
1518
1519 std::vector<double> filter_values = inputdata.GetVector( "FilterRoi" );
1520
1521 if ( filter_values.size()>0 ) roieta = filter_values[0];
1522 if ( filter_values.size()>1 ) roicomposite = ( filter_values[1]==0 ? false : true );
1523 if ( filter_values.size()>2 ) roimult = int(filter_values[2]+0.5);
1524
1525 }
1526
1527 RoiFilter roiFilter( roieta, roicomposite, roimult );
1528
1530
1531 TrackAssociator* matcher = 0;
1532
1533 if ( useMatcher == "Sigma" ) matcher = new Associator_BestSigmaMatcher("sigma", Rmatch);
1534 else if ( useMatcher == "DeltaRZ" || useMatcher == "DeltaRZSinTheta" ) {
1535 double deta = 0.05;
1536 double dphi = 0.05;
1537 double dzed = 25;
1538 if ( inputdata.isTagDefined("Matcher_deta" ) ) deta = inputdata.GetValue("Matcher_deta");
1539 if ( inputdata.isTagDefined("Matcher_dphi" ) ) dphi = inputdata.GetValue("Matcher_dphi");
1540 if ( inputdata.isTagDefined("Matcher_dzed" ) ) dzed = inputdata.GetValue("Matcher_dzed");
1541
1542 if ( useMatcher == "DeltaRZ" ) matcher = new Associator_BestDeltaRZMatcher( "deltaRZ", deta, dphi, dzed );
1543 else matcher = new Associator_BestDeltaRZSinThetaMatcher( "deltaRZ", deta, dphi, dzed );
1544 }
1545 else if ( useMatcher == "pT_2" ) {
1546 double pTmatchLim_2 = 1.0;
1547 if ( inputdata.isTagDefined("Matcher_pTLim_2") ) pTmatchLim_2 = inputdata.GetValue("Matcher_pTLim_2");
1548 matcher = new Associator_SecondBestpTMatcher("SecpT", pTmatchLim_2);
1549 }
1550 else if ( useMatcher == "Truth" ) {
1551 matcher = new Associator_TruthMatcher();
1552 }
1553 else {
1556 matcher = new Associator_BestDeltaRMatcher("deltaR", Rmatch);
1557 }
1558
1560 TrackAssociator* truth_matcher = 0;
1561 if ( truthMatch ) {
1562 if ( useMatcher == "Sigma" ) truth_matcher = new Associator_BestSigmaMatcher("sigma_truth", Rmatch);
1563 else if ( useMatcher == "DeltaRZ" || useMatcher == "DeltaRZSinTheta" ) {
1564 double deta = 0.05;
1565 double dphi = 0.05;
1566 double dzed = 25;
1567 if ( inputdata.isTagDefined("Matcher_deta" ) ) deta = inputdata.GetValue("Matcher_deta");
1568 if ( inputdata.isTagDefined("Matcher_dphi" ) ) dphi = inputdata.GetValue("Matcher_dphi");
1569 if ( inputdata.isTagDefined("Matcher_dzed" ) ) dzed = inputdata.GetValue("Matcher_dzed");
1570
1571 if ( useMatcher == "DeltaRZ" ) truth_matcher = new Associator_BestDeltaRZMatcher( "deltaRZ_truth", deta, dphi, dzed );
1572 else truth_matcher = new Associator_BestDeltaRZSinThetaMatcher( "deltaRZ_truth", deta, dphi, dzed );
1573 }
1574 else if ( useMatcher == "pT_2" ) {
1575 double pTmatchLim_2 = 1.0;
1576 if ( inputdata.isTagDefined("Matcher_pTLim_2") ) pTmatchLim_2 = inputdata.GetValue("Matcher_pTLim_2");
1577 truth_matcher = new Associator_SecondBestpTMatcher("SecpT_truth", pTmatchLim_2);
1578 }
1579 else if ( useMatcher == "Truth" ) {
1580 truth_matcher = new Associator_TruthMatcher();
1581 }
1582 else {
1585 truth_matcher = new Associator_BestDeltaRMatcher("deltaR_truth", Rmatch);
1586 }
1587 }
1588
1589
1590 // NtupleTrackSelector roiTracks( refFilter );
1591
1592
1594
1595 NtupleTrackSelector refPurityTracks( &filter_inout );
1596 NtupleTrackSelector testPurityTracks( &filter_online );
1597
1598 // get the list of input files
1599
1600 std::vector<std::string> filenames;
1601
1602
1603 if ( inputdata.isTagDefined("DataSets") ) {
1604
1605 std::cout << "fetching dataset details" << std::endl;
1606 std::vector<std::string> datasets = inputdata.GetStringVector("DataSets");
1607 for (unsigned int ids=0 ; ids<datasets.size() ; ids++ ) {
1608 std::cout << "\tdataset " << datasets[ids] << std::endl;
1609 dataset d( datasets[ids] );
1610 std::vector<std::string> filenames_ = d.datafiles();
1611 std::cout << "\tdataset contains " << filenames_.size() << " files" << std::endl;
1612 filenames.insert(filenames.end(), filenames_.begin(),filenames_.end());
1613 }
1614 }
1615 else if ( inputdata.isTagDefined("DataFiles") ) filenames = inputdata.GetStringVector("DataFiles");
1616 else {
1617 std::cerr << "no input data specified" << std::endl;
1618 return (-1);
1619 }
1620
1621
1623
1624 // TString* releaseMetaData = 0;
1625 // data->SetBranchAddress("ReleaseMetaData",&releaseMetaData);
1626
1627 bool show_release = true;
1628
1629 std::vector<std::string> release_data;
1630
1631 std::string release_data_save = "";
1632
1633 if ( show_release ){
1634
1635 bool first = true;
1636
1637 for ( unsigned i=0 ; first && i<filenames.size() ; i++ ) {
1638
1639 TFile* finput = TFile::Open( filenames[i].c_str() );
1640
1641
1642 if ( finput==0 || !finput->IsOpen() || finput->IsZombie() ) {
1643 std::cerr << "Error: could not open input file: " << filenames[i] << std::endl;
1644 exit(-1);
1645 }
1646
1647 TTree* dataTree = (TTree*)finput->Get("dataTree");
1648 TString* releaseData = new TString("");
1649
1650 if ( dataTree ) {
1651 dataTree->SetBranchAddress( "ReleaseMetaData", &releaseData);
1652
1653 for (unsigned int i=0; i<dataTree->GetEntries() ; i++ ) {
1654 dataTree->GetEntry(i);
1655 release_data.push_back( releaseData->Data() );
1656 if ( release_data_save != release_data.back() ) {
1657 std::cout << "main() release data: " << release_data.back() << " : " << *releaseData << std::endl;
1658 }
1659 first = false;
1660 release_data_save = release_data.back();
1661 }
1662 }
1663
1664 if ( finput ) delete finput;
1665 }
1666
1667 // for ( unsigned ird=0 ; ird<release_data.size() ; ird++ ) std::cout << "presort " << ird << " " << release_data[ird] << std::endl;
1668
1669 if ( !release_data.empty() ) {
1670 std::sort(release_data.begin(), release_data.end());
1671 release_data.erase(std::unique(release_data.begin(), release_data.end()), release_data.end());
1672
1673 // for ( unsigned ird=0 ; ird<release_data.size() ; ird++ ) std::cout << "postsort " << ird << " " << release_data[ird] << std::endl;
1674
1675 foutdir->cd();
1676
1677 TTree* dataTree = new TTree("dataTree", "dataTree");
1678 TString* releaseData = new TString("");
1679
1680 if ( dataTree ) {
1681 dataTree->Branch( "ReleaseMetaData", "TString", &releaseData);
1682
1683 for ( unsigned ird=0 ; ird<release_data.size() ; ird++ ) {
1684 *releaseData = release_data[ird];
1685 dataTree->Fill();
1686 }
1687
1688 dataTree->Write("", TObject::kOverwrite);
1689 delete dataTree;
1690 }
1691 }
1692
1693 }
1694
1695 // foutput.Write();
1696 // foutput.Close();
1697
1698 // exit(0);
1699
1700
1701
1702
1703
1704
1705
1706 if ( Nentries==0 && inputdata.isTagDefined("Nentries") ) {
1707 Nentries = unsigned(inputdata.GetValue("Nentries"));
1708 }
1709
1710 unsigned event_counter = 0;
1711
1712 typedef std::pair<int,double> zpair;
1713 std::vector<zpair> refz;
1714 std::vector<zpair> testz;
1715
1716 std::vector<double> beamline_ref;
1717 std::vector<double> beamline_test;
1718
1719 int maxtime = 0;
1720 int mintime = 0;
1721
1722 std::cout << "opening files" << std::endl;
1723
1724 bool run = true;
1725
1726 int pregrl_events = 0;
1727 int grl_counter = 0;
1728
1729
1730 std::cout << "starting event loop " << time_str() << std::endl;
1731
1732
1733 size_t max_files = filenames.size();
1734 if ( nfiles!=0 && nfiles<max_files ) max_files = nfiles;
1735
1736 for ( size_t ifile=0 ; run && ifile<max_files; ifile++ ) {
1737
1738 bool newfile = true;
1739
1740
1741 TFile* finput = TFile::Open( filenames[ifile].c_str() );
1742
1743 if ( finput==0 || !finput->IsOpen() || finput->IsZombie() ) {
1744 std::cerr << "Error: could not open output file " << filenames[ifile] << std::endl;
1745 continue;
1746 }
1747
1748 TTree* data = (TTree*)finput->Get("tree");
1749
1750 if ( !data ) {
1751 std::cerr << "Error: cannot open TTree: " << filenames[ifile] << std::endl;
1752 continue;
1753 }
1754
1755 TIDA::Event* track_ev = new TIDA::Event();
1756
1757 gevent = track_ev;
1758
1759 data->SetBranchAddress("TIDA::Event",&track_ev);
1760
1761
1762 maxtime = track_ev->time_stamp();
1763 mintime = track_ev->time_stamp();
1764
1765 unsigned cNentries = data->GetEntries();
1766
1767 bool skip = true;
1768
1771 for (unsigned int i=0; skip && run && i<cNentries ; i++ ) {
1772
1773 data->GetEntry(i);
1774
1775 r = track_ev->run_number();
1776 ev = track_ev->event_number();
1777 lb = track_ev->lumi_block();
1778 ts = track_ev->time_stamp();
1779
1780 int event = track_ev->event_number();
1781 //int bc = track_ev->bunch_crossing_id();
1782
1783
1784 hipt = false;
1785
1786
1787
1788 bool ingrl = goodrunslist.inRange( r, lb );
1789
1790 pregrl_events++;
1791
1793 if ( !ingrl ) continue;
1794
1795 grl_counter++;
1796
1797
1799 if ( event_selector_flag && !es.in( event ) ) continue;
1800
1801 if ( mintime>ts ) mintime = ts;
1802 if ( maxtime<ts ) maxtime = ts;
1803
1804 if ( Nentries>0 && event_counter>Nentries ) {
1805 run = false;
1806 std::cout << "breaking out " << run << std::endl;
1807 break;
1808 }
1809
1810 event_counter++;
1811
1812 // if ( !elist.find(event) ) continue;
1813
1814 // std::cout << "run " << r << "\tevent " << event << "\tlb " << lb << std::endl;
1815
1816 hevent->Fill( event );
1817
1818 if ( filenames.size()<2 ) {
1819 if ( (cNentries<10) || i%(cNentries/10)==0 || i%1000==0 || debugPrintout ) {
1820 std::cout << "run " << track_ev->run_number()
1821 << "\tevent " << track_ev->event_number()
1822 << "\tlb " << track_ev->lumi_block()
1823 << "\tchains " << track_ev->chains().size()
1824 << "\ttime " << track_ev->time_stamp();
1825 std::cout << "\t : processed " << i << " events so far (" << int((1000*i)/cNentries)*0.1 << "%)\t" << time_str() << std::endl;
1826 // std::cerr << "\tprocessed " << i << " events so far \t" << time_str() << std::endl;
1827 }
1828 }
1829 else if ( newfile ) {
1830
1831 int pfiles = filenames.size();
1832 if ( nfiles>0 ) pfiles = nfiles;
1833
1834
1835 std::cout << "file entries=" << data->GetEntries();
1836
1837 if ( data->GetEntries()<100 ) std::cout << " ";
1838 if ( data->GetEntries()<1000 ) std::cout << " ";
1839 if ( data->GetEntries()<10000 ) std::cout << " ";
1840
1841 std::cout << "\t";
1842
1843
1844 std::cout << "run " << track_ev->run_number()
1845 << "\tevent " << track_ev->event_number()
1846 << "\tlb " << track_ev->lumi_block()
1847 << "\tchains " << track_ev->chains().size()
1848 << "\ttime " << track_ev->time_stamp();
1849
1850 std::cout << "\t : processed " << ifile << " files so far (" << int((1e3*ifile)/pfiles)*0.1 << "%)\t" << time_str() << std::endl;
1851
1852 newfile = false;
1853 }
1854
1855 // if ( printflag ) std::cout << *track_ev << std::endl;
1856
1857 r = track_ev->run_number();
1858
1860
1861 offTracks.clear();
1862 refTracks.clear();
1863 truthTracks.clear();
1864 refPurityTracks.clear();
1865
1867 for ( TIDA::ReferenceMap::iterator mit=ref.begin() ; mit!=ref.end() ; ++mit ) {
1868 mit->second.selector()->clear();
1869 dynamic_cast<Filter_Combined*>(mit->second.filter())->setRoi(0);
1870 }
1871
1872 Nvtxtracks = 0;
1873
1874 const std::vector<TIDA::Chain>& chains = track_ev->chains();
1875
1876 dynamic_cast<Filter_Combined*>(truthFilter)->setRoi(0);
1878 if ( truthMatch ) {
1879 for (unsigned int ic=0 ; ic<chains.size() ; ic++ ) {
1880 if ( chains[ic].name()=="Truth" ) {
1881 truthTracks.selectTracks( chains[ic][0].tracks() );
1882 break;
1883 }
1884 }
1885 }
1886
1888 for (const std::string& rc : refChains){
1889 for (unsigned int ic=0 ; ic<chains.size() ; ic++ ) {
1890 if ( chains[ic].name()==rc ) {
1891 offTracks.selectTracks( chains[ic][0].tracks() );
1892 //extract beamline position values from rois
1893 beamline_ref = chains[ic][0].user();
1894 // std::cout << "beamline: " << chains[ic].name() << " " << beamline_ref << std::endl;
1895 break;
1896 }
1897 }
1898 }
1899
1901
1902 std::vector<TIDA::Vertex> vertices; // keep for now as will be needed later ...
1903
1904 const TIDA::Chain* vtxchain = track_ev->chain(vertex_refname);
1905
1906 if ( vtxchain && vtxchain->size()>0 ) {
1907
1908 const std::vector<TIDA::Vertex>& mv = vtxchain->at(0).vertices();
1909
1910 int selectvtx = -1;
1911 double selection = 0;
1912
1913 if ( debugPrintout ) std::cout << "vertices:\n" << mv << std::endl;
1914
1915 if ( bestPTVtx || bestPT2Vtx ) {
1916 for ( size_t iv=0 ; iv<mv.size() ; iv++ ) {
1917 if ( mv[iv].Ntracks()==0 ) continue;
1918 double selection_ = 0.0;
1919 TIDA::Vertex vtx_temp( mv[iv] );
1920 vtx_temp.selectTracks( offTracks.tracks() );
1921 for (unsigned itr=0; itr<vtx_temp.tracks().size(); itr++) {
1922 TIDA::Track* tr = vtx_temp.tracks().at(itr);
1923 if ( bestPTVtx ) selection_ += std::fabs(tr->pT());
1924 else if ( bestPT2Vtx ) selection_ += std::fabs(tr->pT())*std::fabs(tr->pT());
1925 }
1926 if( selection_>selection ) {
1927 selection = selection_;
1928 selectvtx = iv;
1929 }
1930 }
1931 if ( selectvtx!=-1 ) {
1932 vertices.push_back( mv[selectvtx] );
1933 }
1934 }
1935 else if ( vtxind>=0 ) {
1936 if ( size_t(vtxind)<mv.size() ) {
1937 vertices.push_back( mv[vtxind] );
1938 }
1939 }
1940 else {
1941 for ( size_t iv=0 ; iv<mv.size() ; iv++ ) {
1942 vertices.push_back( mv[iv] );
1943 }
1944 }
1945
1947 filter_vertex.setVertex( vertices );
1948
1950
1951 NvtxCount = 0;
1952
1953 for ( unsigned iv=0 ; iv<mv.size() ; iv++ ) {
1954 int Ntracks = mv[iv].Ntracks();
1955 if ( Ntracks>NVtxTrackCut ) {
1956 Nvtxtracks += Ntracks;
1957 // vertices.push_back( mv[iv] );
1958 NvtxCount++;
1959 }
1960 }
1961 }
1962
1963 // filter_vertex.setVertex( vvtx ) ;
1964
1965 hcorr->Fill( vertices.size(), Nvtxtracks );
1966
1967 dynamic_cast<Filter_Combined*>(refFilter)->setRoi(0);
1968
1972 bool foundReference = false;
1973
1974 const TIDA::Chain* refchain = 0;
1975
1978
1979 TrigObjectMatcher tom;
1980
1981 for ( const std::string& rc : refChains ) {
1982 foundReference |= GetRefTracks( rc, exclude, chains, refTracks, ex_matcher, tom );
1983 }
1984
1986 // bool skip_tnp = false;
1987
1988 for ( TIDA::ReferenceMap::iterator mitr=ref.begin() ; mitr!=ref.end() ; ++mitr ) {
1989
1990 std::string refname = mitr->first;
1991
1992 NtupleTrackSelector& selector = *dynamic_cast<NtupleTrackSelector*>( mitr->second.selector() );
1993
1995 TrigObjectMatcher rtom;
1996
1997 foundReference |= GetRefTracks( refname, exclude, chains, selector, ex_matcher, rtom );
1998
1999 *mitr->second.tom() = rtom;
2000
2004 // if ( refname.find("Tau")!=std::string::npos || refname.find("Electron")!=std::string::npos ) {
2005 // if ( rtom.status()==0 ) skip_tnp = true; /// maybe don't use this at the moment
2006 // }
2007
2008 }
2009
2010
2011
2012
2013 if ( !foundReference ) continue;
2014
2015 if ( debugPrintout ) {
2016 std::cout << "reference chain:\n" << *refchain << std::endl;
2017 }
2018
2019 for ( unsigned ic=0 ; ic<track_ev->chains().size() ; ic++ ) {
2020
2021 TIDA::Chain& chain = track_ev->chains()[ic];
2022
2023 // std::cout << ic << " chain " << chain.name() << " size " << chain.size() << std::endl;
2024
2026 std::map<std::string,TrackAnalysis*>::iterator analitr = analysis.find(chain.name());
2027
2029
2030 if ( analitr==analysis.end() ) continue;
2031
2032 if ( debugPrintout ) {
2033 std::cout << "test chain:\n" << chain << std::endl;
2034 }
2035
2036 ConfAnalysis* cf = dynamic_cast<ConfAnalysis*>(analitr->second);
2037
2038 std::vector<TIDA::Roi*> rois;
2039
2040 // tag and probe object retreival and filling of the roi vector
2041 TagNProbe* TnP_tool = cf->getTnPtool();
2042
2043 if ( TnP_tool ) {
2044
2045 foutdir->cd();
2046 cf->initialiseInternal();
2047 // changes to output directory and books the invariant mass histograms
2048 TH1F* invmass = cf->getHist_invmass();
2049 TH1F* invmass_obj = cf->getHist_invmassObj();
2050
2051 std::map<std::string,TIDA::Reference>::iterator mit0=ref.find(TnP_tool->type0());
2052
2053 TrackSelector* selector0 = mit0->second.selector();
2054 TrackFilter* filter0 = mit0->second.filter();
2055
2056 TrigObjectMatcher* rtom = mit0->second.tom();
2057
2058 if ( TnP_tool->type0() == TnP_tool->type1() ) {
2060 rois = TnP_tool->GetRois( track_ev->chains(), selector0, filter0, invmass, invmass_obj, rtom );
2061 }
2062 else {
2064
2066 std::map<std::string,TIDA::Reference>::iterator mit1=ref.find(TnP_tool->type1());
2067
2068 TrackSelector* selector1 = mit1->second.selector();
2069 TrackFilter* filter1 = mit1->second.filter();
2070
2071 TrigObjectMatcher* rtom1 = mit1->second.tom();
2072
2073#if 0
2074 for ( size_t ia=0 ; ia<selector1->tracks().size() ; ia++ ) {
2075 const TIDA::Track* track = selector1->tracks()[ia];
2076 std::cout << ia << "\t" << *track << std::endl;
2077 std::cout << "\t" << rtom1->object(track->id()) << std::endl;
2078 }
2079
2080 std::cout << *TnP_tool << std::endl;
2081#endif
2082
2083 rois = TnP_tool->GetRois( track_ev->chains(), selector0, filter0, selector1, filter1, invmass, invmass_obj, rtom, rtom1 );
2084 }
2085 }
2086 else {
2087 // if not a tnp analysis then fill rois in the normal way
2088 rois.reserve( chain.size() );
2089 for ( size_t ir=0 ; ir<chain.size() ; ir++ ) rois.push_back( &(chain.rois()[ir]) );
2090 }
2091
2093 if ( false ) std::cout << "++++++++++++ rois size = " << rois.size() << " +++++++++++" << std::endl;
2094
2095 //for (unsigned int ir=0 ; ir<chain.size() ; ir++ ) {
2096 for (unsigned int ir=0 ; ir<rois.size() ; ir++ ) { // changed for tagNprobe
2097
2099
2100 //TIDA::Roi& troi = chain.rois()[ir];
2101 TIDA::Roi& troi = *rois[ir]; // changed for tagNprobe
2102 TIDARoiDescriptor roi( troi.roi() );
2103
2104 testTracks.clear();
2105
2106 testTracks.selectTracks( troi.tracks() );
2107
2109 std::vector<TIDA::Track*> testp = testTracks.tracks();
2110
2113
2114 if ( filterRoi && !roiFilter.filter( roi ) ) continue;
2115
2117 const std::vector<TIDA::Vertex>& mvt = troi.vertices();
2118
2119 std::vector<TIDA::Vertex> vertices_test;
2120
2121 int selectvtx = -1;
2122 double selection = 0;
2123
2124 if ( bestPTVtx_rec || bestPT2Vtx_rec ) {
2125
2126 // const std::vector<TIDA::Track>& recTracks = troi.tracks();
2127
2128 for ( unsigned iv=0 ; iv<mvt.size() ; iv++ ) {
2129 double selection_ = 0.0;
2130 TIDA::Vertex vtx_temp( mvt[iv] );
2131 vtx_temp.selectTracks( testp );
2132 for (unsigned itr=0; itr<vtx_temp.tracks().size(); itr++) {
2133 TIDA::Track* tr = vtx_temp.tracks().at(itr);
2134 if ( bestPTVtx ) selection_ += std::fabs(tr->pT());
2135 else if ( bestPT2Vtx ) selection_ += std::fabs(tr->pT())*std::fabs(tr->pT());
2136 }
2137 if( selection_>selection){
2138 selection = selection_;
2139 selectvtx = iv;
2140 }
2141 }
2142 if ( selectvtx!=-1 ) {
2143 TIDA::Vertex selected( mvt[selectvtx] );
2144 if ( useVertexTracks ) selected.selectTracks( testp );
2145 vertices_test.push_back(selected);
2146 }
2147 }
2148 else if ( vtxind_rec!=-1 ) {
2149 if ( unsigned(vtxind_rec)<mvt.size() ) {
2150 TIDA::Vertex selected( mvt[vtxind] );
2151 if ( useVertexTracks ) selected.selectTracks( testp );
2152 vertices_test.push_back( selected );
2153 }
2154 }
2155 else {
2156 for ( unsigned iv=0 ; iv<mvt.size() ; iv++ ) {
2157 TIDA::Vertex selected( mvt[iv] );
2158 if ( useVertexTracks ) selected.selectTracks( testp );
2159 vertices_test.push_back( selected );
2160 }
2161 }
2162
2163 //extract beamline position values from rois
2164 beamline_test = rois[ir]->user(); // changed for tagNprobe
2165
2166 //set values of track analysis to these so can access elsewhere
2167 for ( size_t i=analyses.size() ; i-- ; ) {
2168
2169 TrackAnalysis* analy_track = analyses[i];
2170
2171 if ( correctBeamlineTest ) {
2172 if ( beamTest.size()==2 ) analy_track->setBeamTest( beamTest[0], beamTest[1] );
2173 // else if ( beamTest.size()==3 ) analy_track->setBeamTest( beamTest[0], beamTest[1], beamTest[2] );
2174 else {
2175 if ( !inputdata.isTagDefined("BeamTest") ) {
2176 if ( beamline_test.size()==2 ) analy_track->setBeamTest( beamline_test[0], beamline_test[1] );
2177 // else if ( beamline_test.size()==3 ) analy_track->setBeamTest( beamline_test[0], beamline_test[1], beamline_test[2] );
2178 }
2179 }
2180 }
2181
2182 if ( correctBeamlineRef ) {
2183 if ( beamRef.size()==2 ) analy_track->setBeamRef( beamRef[0], beamRef[1] );
2184 // else if ( beamRef.size()==3 ) analy_track->setBeamRef( beamRef[0], beamRef[1], beamRef[2] );
2185 else {
2186 if ( !inputdata.isTagDefined("BeamRef") ) {
2187 if ( beamline_ref.size()==2 ) analy_track->setBeamRef( beamline_ref[0], beamline_ref[1] );
2188 // else if ( beamline_ref.size()==3 ) analy_track->setBeamRef( beamline_ref[0], beamline_ref[1], beamline_ref[2] );
2189 }
2190 }
2191 }
2192
2193 }
2194
2197
2198 TIDARoiDescriptor refRoi;
2199
2200 if ( select_roi ) {
2201
2205
2206 bool customRefRoi_thisChain = false;
2207
2208 if ( use_custom_ref_roi ) { // Ideally just want to say ( use_custom_ref_roi && (chain.name() in custRefRoi_chain]sist) )
2209 if ( customRoi_chains.size() ) {
2210 if ( customRoi_chains.find( chain.name() )!=customRoi_chains.end() ) customRefRoi_thisChain = true;
2211 }
2212 else customRefRoi_thisChain = true; // Apply custom ref roi to all chains
2213 }
2214
2215 if ( use_custom_ref_roi && customRefRoi_thisChain ) {
2216 refRoi = makeCustomRefRoi( roi, custRefRoi_params[0], custRefRoi_params[1], custRefRoi_params[2] );
2217 }
2218 else refRoi = roi;
2219
2220 dynamic_cast<Filter_Combined*>(refFilter)->setRoi(&refRoi);
2221 }
2222
2223 // if ( debugPrintout ) {
2224 // std::cout << "refTracks.size() " << refTracks.size() << " before roi filtering" << std::endl;
2225 // std::cout << "filter with roi " << roi << std::endl;
2226 // }
2227
2228
2229 // Now filterng ref tracks by refFilter, and performing any further filtering and selecting,
2230 // before finally creating the const reference object refp
2231
2232 std::vector<TIDA::Track*> refp_vec = refTracks.tracks( refFilter );
2233
2234 // Selecting only truth matched reference tracks
2235 if ( truthMatch ) {
2237 if ( select_roi ) dynamic_cast<Filter_Combined*>(truthFilter)->setRoi(&roi);
2238 const std::vector<TIDA::Track*>& truth = truthTracks.tracks(truthFilter);
2239 const std::vector<TIDA::Track*>& refp_tmp = refp_vec;
2240
2242 truth_matcher->match( refp_tmp, truth );
2243
2244 std::vector<TIDA::Track*> refp_matched;
2245
2247 for ( unsigned i=0 ; i<refp_vec.size() ; i++ ) {
2248 if ( truth_matcher->matched( refp_vec[i] ) ) refp_matched.push_back( refp_vec[i] );
2249 }
2250
2251 refp_vec.clear();
2252 refp_vec = refp_matched;
2253 }
2254
2255 // Choose the pT ordered refp tracks that have been asked for by the user
2256 if ( use_pt_ordered_ref ) {
2257 std::sort( refp_vec.begin(), refp_vec.end(), trackPtGrtr ); // Should this sorting be done to a temporary copied object, instead of the object itself?
2258
2259 size_t nRefTracks = refp_vec.size();
2260
2261 std::vector<TIDA::Track*> refp_chosenPtOrd;
2262
2263 if ( debugPrintout ) {
2264 // Checking if any user specifed track indices are out of bounds for this event
2265 for ( size_t sidx=refPtOrd_indices.size() ; sidx-- ; ) {
2266 if ( refPtOrd_indices.at(sidx)>nRefTracks ) {
2267 std::cout << "WARNING: for event " << event
2268 << ", pT ordered reference track at vector position " << refPtOrd_indices.at(sidx)
2269 << " requested but not found" << std::endl;
2270 break;
2271 }
2272 }
2273 }
2274
2275 for ( size_t sidx=refPtOrd_indices.size() ; sidx-- ; ) {
2276 for ( size_t idx=0 ; idx<refp_vec.size() ; idx++ ) {
2277 if ( idx == refPtOrd_indices.at(sidx) ) {
2278 refp_chosenPtOrd.push_back(refp_vec.at(idx));
2279 break;
2280 }
2281 }
2282 }
2283
2284 refp_vec.clear();
2285 refp_vec = refp_chosenPtOrd; // Note: not necessarily pT ordered.
2286 // Ordered by order of indices the user has passed
2287 // (which is ideally pT ordered e.g. 0, 1, 3)
2288 }
2289
2291 if ( cf ) {
2292
2294
2295 std::string ptconfig = cf->config().postvalue("pt");
2296 if ( ptconfig!="" ) {
2297 double pt = std::atof( ptconfig.c_str() );
2298 if ( pt>0 ) {
2299 std::vector<TIDA::Track*> reft; reft.reserve(refp_vec.size());
2300 for ( std::vector<TIDA::Track*>::const_iterator itr=refp_vec.begin() ; itr!=refp_vec.end() ; ++itr ) {
2301 if ( std::fabs((*itr)->pT())>=pt ) reft.push_back( *itr );
2302 }
2303 refp_vec = reft;
2304 }
2305 }
2306
2308
2309 std::string d0config = cf->config().postvalue("d0");
2310 if ( d0config!="" ) {
2311 double d0 = std::atof( d0config.c_str() );
2312 if ( d0>0 ) {
2313 std::vector<TIDA::Track*> reft; reft.reserve(refp_vec.size());
2314 for ( std::vector<TIDA::Track*>::const_iterator itr=refp_vec.begin() ; itr!=refp_vec.end() ; ++itr ) {
2315 if ( std::fabs((*itr)->a0())<=d0 ) reft.push_back( *itr );
2316 }
2317 refp_vec = reft;
2318 }
2319 }
2320
2321 }
2322
2325
2327
2328 if ( tom.status() ) {
2329
2330 std::string objectET = cf->config().postvalue("ET");
2331
2332 if ( objectET != "" ) {
2333
2334 double ETconfig = std::atof( objectET.c_str() );
2335
2336 if ( ETconfig>0 ) {
2337
2338 std::vector<TIDA::Track*>::iterator itr=refp_vec.begin() ;
2339
2340 while ( itr!=refp_vec.end() ) {
2341 const TrackTrigObject* tobj = tom.object( (*itr)->id() );
2342
2343 if ( tobj==0 || tobj->pt()<ETconfig )
2344 itr=refp_vec.erase( itr );
2345 else
2346 ++itr;
2347 }
2348 }
2349 }
2350 }
2351
2352 const std::vector<TIDA::Track*>& refp = refp_vec;
2353
2354 // if ( debugPrintout ) {
2355 // std::cout << "refp.size() " << refp.size() << " after roi filtering" << std::endl;
2356 // }
2357
2358
2359 // std::cout << "ref tracks refp.size() " << refp.size() << "\n" << refp << std::endl;
2360 // std::cout << "test tracks testp.size() " << testp.size() << "\n" << testp << std::endl;
2361
2362 groi = &roi;
2363
2367
2368 // new vertex class containing tracks, offline
2369 std::vector<TIDA::Vertex> vertices_roi;
2370
2372 // if ( chain.name().find("SuperRoi") ) {
2373 {
2374
2376
2377 vertices_roi.clear();
2378
2379 const std::vector<TIDA::Vertex>& mv = vertices;
2380
2381 // std::cout << "vertex filtering " << mv.size() << std::endl;
2382
2383
2384 for ( unsigned iv=0 ; iv<mv.size() ; iv++ ) {
2385
2386 const TIDA::Vertex& vx = mv[iv];
2387
2388 // reject all vertices that are not in the roi
2389
2390 bool accept_vertex = false;
2391 if ( roi.composite() ) {
2392 for ( size_t ir=0 ; ir<roi.size() ; ir++ ) {
2393 if ( roi[ir]->zedMinus()<=vx.z() && roi[ir]->zedPlus()>=vx.z() ) accept_vertex = true;
2394 }
2395 }
2396 else {
2397 if ( roi.zedMinus()<=vx.z() && roi.zedPlus()>=vx.z() ) accept_vertex = true;
2398 }
2399
2400 if ( !accept_vertex ) continue;
2401
2402 // std::cout << "\t" << iv << "\t" << vx << std::endl;
2403
2404 int trackcount = 0;
2405
2406 if ( useVertexTracks ) {
2407 // refp contains roi filtered tracks, vx contains ids of tracks belonging to vertex
2408 TIDA::Vertex vertex_roi( vx );
2409 vertex_roi.selectTracks( refp );
2410 trackcount = vertex_roi.Ntracks();
2411 if ( trackcount>=ntracks && trackcount>0 ) {
2412 vertices_roi.push_back( vertex_roi );
2413 }
2414 }
2415 else {
2416 // old track count method still in use?
2417 for (unsigned itr=0; itr<refp.size(); itr++){
2418 TIDA::Track* tr = refp[itr];
2419 double theta_ = 2*std::atan(std::exp(-tr->eta()));
2420 double dzsintheta = std::fabs( (vx.z()-tr->z0()) * std::sin(theta_) );
2421 if( dzsintheta < 1.5 ) trackcount++;
2422 }
2426 if ( trackcount>=ntracks && trackcount>0 ) {
2427 const TIDA::Vertex& vertex_roi( vx );
2428 vertices_roi.push_back( vertex_roi );
2429 }
2430 }
2431
2432 }
2433 }
2434 // else vertices_roi = vertices;
2435
2436 if ( rotate_testtracks ) for ( size_t i=testp.size() ; i-- ; ) testp[i]->rotate();
2437
2438 foutdir->cd();
2439
2440 // do analysing
2441
2442 if ( monitorZBeam ) {
2443 if ( beamline_ref.size()>2 && beamline_test.size()>2 ) {
2444 refz.push_back( zpair( lb, beamline_ref[2]) );
2445 testz.push_back( zpair( lb, beamline_test[2]) );
2446 }
2447 }
2448
2449 matcher->match( refp, testp);
2450
2451 if ( tom.status() ) analitr->second->execute( refp, testp, matcher, &tom );
2452 else analitr->second->execute( refp, testp, matcher );
2453
2454 ConfVtxAnalysis* vtxanal = 0;
2455 analitr->second->store().find( vtxanal, "rvtx" );
2456
2459 if ( vtxanal && ( refp.size() >= size_t(ntracks) ) ) {
2460
2468
2470
2471 if ( vertices_roi.size()>0 ) vtxanal->execute( pointers(vertices_roi), pointers(vertices_test), track_ev );
2472
2473 }
2474
2475
2476 if ( debugPrintout ) {
2477 // std::cout << "-----------------------------------\n\nselected tracks:" << chain.name() << std::endl;
2478 std::cout << "\nselected tracks:" << chain.name() << std::endl;
2479 std::cout << "ref tracks refp.size() " << refp.size() << "\n" << refp << std::endl;
2480 std::cout << "test tracks testp.size() " << testp.size() << "\n" << testp << std::endl;
2481
2482 TrackAssociator::map_type::const_iterator titr = matcher->TrackAssociator::matched().begin();
2483 TrackAssociator::map_type::const_iterator tend = matcher->TrackAssociator::matched().end();
2484 int im=0;
2485 std::cout << "track matches:\n";
2486 while (titr!=tend) {
2487 std::cout << "\t" << im++ << "\t" << *titr->first << " ->\n\t\t" << *titr->second << std::endl;
2488 ++titr;
2489 }
2490
2491
2492 std::cout << "completed : " << chain.name() << "\n";
2493 std::cout << "-----------------------------------" << std::endl;
2494
2495 }
2496
2497#if 0
2498 if ( _matcher->size()<refp.size() ) {
2499
2500 if ( refp.size()==1 && testp.size()==0 ) {
2501 std::cout << track_ev->chains()[ic] <<std::endl;
2502 std::cout << "roi\n" << track_ev->chains()[ic].rois()[ir] << endl;
2503 }
2504
2505 }
2506#endif
2507
2509
2510 if ( doPurity ) {
2511
2512 const std::vector<TIDA::Track*>& refpp = refPurityTracks.tracks( refFilter );
2513
2514 testPurityTracks.clear();
2515
2516 testPurityTracks.selectTracks( troi.tracks() );
2517 std::vector<TIDA::Track*> testpp = testPurityTracks.tracks();
2518
2519 matcher->match(refpp, testpp);
2520
2521
2522 std::map<std::string,TrackAnalysis*>::iterator analitrp = analysis.find(chain.name()+"-purity");
2523
2524 if ( analitrp == analysis.end() ) continue;
2525
2526
2527 analitrp->second->execute( refpp, testpp, matcher );
2528
2529
2530
2531 }
2532
2533 } // loop through rois
2534
2535 } // loop through chanines - this block is indented incorrectly
2536
2537 } // loop through nentries
2538
2539 delete track_ev;
2540 delete data;
2541 delete finput;
2542
2543 // std::cout << "run: " << run << std::endl;
2544
2545 }// loop through files
2546
2547 std::cout << "done " << time_str() << "\tprocessed " << event_counter << " events"
2548 << "\ttimes " << mintime << " " << maxtime
2549 << "\t( grl: " << grl_counter << " / " << pregrl_events << " )" << std::endl;
2550
2551 if ( monitorZBeam ) zbeam zb( refz, testz );
2552
2553 foutdir->cd();
2554
2555 // hcorr->Finalise(Resplot::FitPoisson);
2556
2557 hcorr->Finalise(Resplot::FitNull95);
2558 hcorr->Write();
2559
2560 for ( int i=analyses.size() ; i-- ; ) {
2561 // std::cout << "finalise analysis chain " << analyses[i]->name() << std::endl;
2562 analyses[i]->finalise();
2563
2564 ConfVtxAnalysis* vtxanal = 0;
2565 analyses[i]->store().find( vtxanal, "rvtx" );
2566 if ( vtxanal ) vtxanal->finalise();
2567
2568 delete analyses[i];
2569 }
2570
2571 foutput.Write();
2572 foutput.Close();
2573
2574
2575 // for ( std::map<std::string,TIDA::Reference>::iterator mit=ref.begin() ; mit!=ref.end() ; ++mit++ ) {
2576 for ( TIDA::ReferenceMap::iterator mit=ref.begin() ; mit!=ref.end() ; ++mit++ ) {
2577 mit->second.Clean();
2578 }
2579
2580
2581 if ( ex_matcher ) delete ex_matcher;
2582
2583 std::cout << "done" << std::endl;
2584
2585 return 0;
2586}
ref
const boost::regex ref(r_ef)
Nvtxtracks
int Nvtxtracks
Definition globals.cxx:17
cosmicBinConfig
BinConfig cosmicBinConfig("cosmic")
electronBinConfig
BinConfig electronBinConfig("electron")
NvtxCount
int NvtxCount
Definition globals.cxx:18
muonBinConfig
BinConfig muonBinConfig("muon")
PRINT_BRESIDUALS
bool PRINT_BRESIDUALS
stack trace headers
Definition ConfAnalysis.cxx:29
g_binConfig
BinConfig g_binConfig("standard")
bjetBinConfig
BinConfig bjetBinConfig("bjet")
tauBinConfig
BinConfig tauBinConfig("tau")
gevent
TIDA::Event * gevent
Definition globals.cxx:13
data
char data[hepevt_bytes_allocation_ATLAS]
Definition HepEvt.cxx:11
TrackAssociator
TIDA::Associator< TIDA::Track > TrackAssociator
Definition TrackAssociator.h:26
Associator_BestDeltaRMatcher
Definition Associator_BestMatch.h:80
Associator_BestDeltaRZMatcher
Definition Associator_BestMatch.h:109
Associator_BestDeltaRZSinThetaMatcher
Definition Associator_BestMatch.h:152
Associator_BestSigmaMatcher
Definition Associator_BestMatch.h:197
Associator_SecondBestpTMatcher
Definition Associator_BestMatch.h:58
Associator_TruthMatcher
Definition Associator_TruthMatch.h:44
ChainString::extra
const std::string & extra() const
Definition ChainString.h:38
ChainString::head
const std::string & head() const
Definition ChainString.h:33
ChainString::postvalue
std::string postvalue(const std::string &key) const
same here regarding returning a reference
Definition ChainString.h:54
ChainString::element
const std::string & element() const
Definition ChainString.h:37
ConfAnalysis::getHist_invmass
virtual TH1F * getHist_invmass()
Definition ConfAnalysis.h:126
ConfAnalysis::config
const ChainString & config() const
Definition ConfAnalysis.h:118
ConfAnalysis::getTnPtool
virtual TagNProbe * getTnPtool()
Definition ConfAnalysis.h:124
ConfAnalysis::setprint
void setprint(bool p)
Definition ConfAnalysis.h:112
ConfAnalysis::initialise
virtual void initialise()
book all the histograms
Definition ConfAnalysis.cxx:59
ConfAnalysis::getHist_invmassObj
virtual TH1F * getHist_invmassObj()
Definition ConfAnalysis.h:128
ConfAnalysis::initialiseInternal
virtual void initialiseInternal()
Definition ConfAnalysis.cxx:66
ConfAnalysis::initialiseFirstEvent
void initialiseFirstEvent(bool b=true)
Definition ConfAnalysis.h:116
ConfVtxAnalysis::execute
void execute(const std::vector< TIDA::Vertex * > &vtx0, const std::vector< TIDA::Vertex * > &vtx1, const TIDA::Event *tevt=0)
Definition ConfVtxAnalysis.cxx:186
Filter_Author
Definition Filters.h:32
Filter_Combined
Definition Filters.h:225
Filter_Combined::setRoi
void setRoi(TIDARoiDescriptor *r)
Definition Filters.h:236
Filter_Offline2017
Definition Filter_Offline2017.h:28
Filter_OfflineR22
Definition Filter_OfflineR22.h:28
Filter_TrackQuality
Definition Filter_TrackQuality.h:27
Filter_Track
sadly need to include a root dependency, but no matter - the TIDA::Track class itself inherets from T...
Definition Filter_Track.h:41
Filter_Vertex
Definition Filters.h:155
Filter_etaPT
Definition Filter_etaPT.h:27
Filter_pdgId
Definition Filters.h:53
PurityAnalysis::setprint
void setprint(bool p)
Definition PurityAnalysis.h:73
PurityAnalysis::initialise
virtual void initialise()
book all the histograms
Definition PurityAnalysis.cxx:34
ReadCards
Get tag-value pairs from a file.
Definition ReadCards.h:50
Resplot::Finalise
int Finalise(double a=-999, double b=-999, TF1 *(*func)(TH1D *s, double a, double b)=Resplot::FitGaussian)
Definition Resplot.cxx:387
Resplot::Write
Int_t Write(const char *=0, Int_t=0, Int_t=0) const
Hooray, this stupidity is to overwride both the const and non-const TObject Write methods Fixme: shou...
Definition Resplot.h:454
Resplot::FitNull95
static TF1 * FitNull95(TH1D *s, double a=0, double b=0)
Definition Resplot.cxx:1672
Resplot::setoldrms95
static bool setoldrms95(bool b)
Definition Resplot.h:555
Resplot::Fill
int Fill(double x, double y, double w=1)
Definition Resplot.h:281
Resplot::setnofit
static bool setnofit(bool b)
Definition Resplot.h:557
TIDARoiDescriptor
Describes the Region of Ineterest geometry It has basically 8 parameters.
Definition TIDARoiDescriptor.h:40
TIDA::Chain::size
unsigned size() const
number of rois
Definition TIDAChain.h:54
TIDA::Chain::at
TIDA::Roi & at(int i)
Definition TIDAChain.h:63
TIDA::Event::chain
const TIDA::Chain * chain(const std::string &s) const
Definition TIDAEvent.cxx:57
TIDA::Event::lumi_block
void lumi_block(unsigned lb)
Definition TIDAEvent.h:44
TIDA::Event::event_number
void event_number(unsigned long long e)
Definition TIDAEvent.h:43
TIDA::Event::time_stamp
void time_stamp(unsigned t)
Definition TIDAEvent.h:45
TIDA::Event::run_number
void run_number(unsigned r)
accessors
Definition TIDAEvent.h:42
TIDA::Event::chains
const std::vector< TIDA::Chain > & chains() const
Definition TIDAEvent.h:76
TIDA::FeatureStore::insert
void insert(T *t, const std::string &key)
Definition TIDAFeatureStore.h:44
TIDA::ReferenceMap::iterator
std::map< std::string, Reference >::iterator iterator
Definition TIDAReference.h:72
TIDA::ReferenceMap::value_type
std::map< std::string, Reference >::value_type value_type
Definition TIDAReference.h:74
TIDA::Roi::roi
const TIDARoiDescriptor & roi() const
access the roi information
Definition TIDARoi.h:42
TIDA::Roi::tracks
const std::vector< TIDA::Track > & tracks() const
Definition TIDARoi.h:52
TIDA::Roi::vertices
const std::vector< TIDA::Vertex > & vertices() const
access the vertices
Definition TIDARoi.h:56
TIDA::Vertex::z
double z() const
Definition TIDAVertex.h:53
TagNProbe::probe
void probe(const std::string &chainName)
Definition TagNProbe.h:46
TagNProbe::GetRois
std::vector< TIDA::Roi * > GetRois(std::vector< TIDA::Chain > &chains, const TrackSelector *selector, TrackFilter *filter, T *hmass, T *hmass_obj, TrigObjectMatcher *tom=0) const
Definition TagNProbe.h:60
TagNProbe::tag
void tag(const std::string &chainName)
getters and setters
Definition TagNProbe.h:45
TagNProbe::type0
const std::string & type0() const
Definition TagNProbe.h:51
TagNProbe::type1
const std::string & type1() const
Definition TagNProbe.h:52
TrackAnalysis::store
TIDA::FeatureStore & store()
Definition TrackAnalysis.h:133
TrackAnalysis::setBeamTest
void setBeamTest(double x, double y, double z=0)
Definition TrackAnalysis.h:110
TrackAnalysis::setBeamRef
void setBeamRef(double x, double y, double z=0)
set the beamline positions
Definition TrackAnalysis.h:109
TrackTrigObject::pt
double pt() const
Definition TrackTrigObject.h:41
TrigObjectMatcher::object
const TrackTrigObject * object(unsigned long track_id)
Definition TrigObjectMatcher.h:59
event_selector
Definition event_selector.h:20
zbeam
Definition zbeam.h:24
ir
int ir
counter of the current depth
Definition fastadd.cxx:49
selection
const std::string selection
Definition fbtTestBasics.cxx:75
usage
StatusCode usage()
Definition fbtTestToyMC.cxx:1168
dumpflag
bool dumpflag
Definition globals.cxx:30
groi
TIDARoiDescriptor * groi
all these externals initialised in globals.cxx
Definition globals.cxx:14
NMod
int NMod
Definition globals.cxx:20
lb
int lb
Definition globals.cxx:23
ts
int ts
Definition globals.cxx:24
hipt
bool hipt
Definition globals.cxx:29
r
int r
Definition globals.cxx:22
a0
double a0
Definition globals.cxx:27
ev
int ev
Definition globals.cxx:25
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition hcg.cxx:310
find
std::string find(const std::string &s)
return a remapped string
Definition hcg.cxx:138
CaloCondBlobAlgs_fillNoiseFromASCII.tag
str tag
Definition CaloCondBlobAlgs_fillNoiseFromASCII.py:23
CxxUtils::atoi
int atoi(std::string_view str)
Helper functions to unpack numbers decoded in string into integers and doubles The strings are requir...
Definition Control/CxxUtils/Root/StringUtils.cxx:85
DeMoScan.first
bool first
Definition DeMoScan.py:534
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition LArNewCalib_DelayDump_OFC_Cali.py:69
PowhegPythia8EvtGen_H2a4X_ctauY.newfile
str newfile
Definition PowhegPythia8EvtGen_H2a4X_ctauY.py:13
RivetAnalysis_JO.analyses
str analyses
Definition RivetAnalysis_JO.py:30
example.goodrunslist
goodrunslist
Definition example.py:26
grepfile.filenames
list filenames
Definition grepfile.py:34
hist_file_dump.d
d
Definition hist_file_dump.py:142
python.TrigEgammaMonitorHelper.TH1F
TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition TrigEgammaMonitorHelper.py:24
python.atlas_oh.im
im
Definition atlas_oh.py:167
python.utility.LHE.event_counter
event_counter(input_file_pattern)
Count total number of events in input files.
Definition LHE.py:89
run
Definition run.py:1
std::unique
DataModel_detail::iterator< DVL > unique(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of unique for DataVector/List.
Definition DVL_algorithms.h:135
std::sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554
test_pyathena.pt
pt
Definition test_pyathena.py:11
trackPtGrtr
bool trackPtGrtr(TIDA::Track *trackA, TIDA::Track *trackB)
Definition rmain.cxx:205
getFilter
TrackFilter * getFilter(const std::string &refname, int pdgId, TrackFilter *foff, TrackFilter *fmu, TrackFilter *ftruth)
This is awful code, passing in lots of filter pointers just so that they can be assigned neatly ?
Definition rmain.cxx:404
pointers
std::vector< T * > pointers(std::vector< T > &v)
Definition rmain.cxx:367
handler
void handler(int sig)
signal handler
Definition rmain.cxx:99
ETmin
double ETmin
Definition rmain.cxx:375
makeCustomRefRoi
TIDARoiDescriptor makeCustomRefRoi(const TIDARoiDescriptor &roi, double etaHalfWidth=-999, double phiHalfWidth=-999, double zedHalfWidth=-999)
Definition rmain.cxx:165
GetRefTracks
bool GetRefTracks(const std::string &rc, const std::string &exclude, const std::vector< TIDA::Chain > &chains, NtupleTrackSelector &refTracks, TrackAssociator *ex_matcher, TrigObjectMatcher &tom)
Definition rmain.cxx:429
atoi_check
int atoi_check(const std::string &s)
convert string to integer with check if successful
Definition rmain.cxx:129
ETovPTmin
double ETovPTmin
this is a swiss knife function - by default if ET/PT > 0 such that fabs(ET/PT) > 0 is always true and...
Definition rmain.cxx:386
time_str
std::string time_str()
Definition rmain.cxx:116
event_list
Definition rmain.cxx:299

◆ makeCustomRefRoi()

TIDARoiDescriptor makeCustomRefRoi ( const TIDARoiDescriptor & roi,
double etaHalfWidth = -999,
double phiHalfWidth = -999,
double zedHalfWidth = -999 )

Definition at line 165 of file rmain.cxx.

168 {
169
170 double roi_eta = roi.eta();
171 double roi_phi = roi.phi();
172 double roi_zed = roi.zed();
173
174 double etaMinus = roi.etaMinus();
175 double etaPlus = roi.etaPlus();
176
177 double zedMinus = roi.zedMinus();
178 double zedPlus = roi.zedPlus();
179
180 double phiMinus = roi.phiMinus();
181 double phiPlus = roi.phiPlus();
182
183 if ( etaHalfWidth != -999 ) {
184 etaMinus = roi_eta - etaHalfWidth;
185 etaPlus = roi_eta + etaHalfWidth;
186 }
187
188 if ( phiHalfWidth != -999 ) {
189 phiMinus = roi_phi - phiHalfWidth; // !!! careful! will this always wrap correctly?
190 phiPlus = roi_phi + phiHalfWidth; // !!! careful! will this always wrap correctly?
191 }
192
193 if ( zedHalfWidth != -999 ) {
194 zedMinus = roi_zed - zedHalfWidth;
195 zedPlus = roi_zed + zedHalfWidth;
196 }
197
198 return TIDARoiDescriptor( roi_eta, etaMinus, etaPlus,
199 roi_phi, phiMinus, phiPlus, // do wrap phi here (done within TIDARoiDescriptor already)
200 roi_zed, zedMinus, zedPlus );
201}
TIDARoiDescriptor::etaMinus
double etaMinus() const
Definition TIDARoiDescriptor.h:135
TIDARoiDescriptor::zedPlus
double zedPlus() const
Definition TIDARoiDescriptor.h:142
TIDARoiDescriptor::etaPlus
double etaPlus() const
Definition TIDARoiDescriptor.h:136
TIDARoiDescriptor::zedMinus
double zedMinus() const
Definition TIDARoiDescriptor.h:141
TIDARoiDescriptor::phiPlus
double phiPlus() const
Definition TIDARoiDescriptor.h:139
TIDARoiDescriptor::phiMinus
double phiMinus() const
Definition TIDARoiDescriptor.h:138
IDTrig_MC23a_preInclude.zedHalfWidth
zedHalfWidth
Definition IDTrig_MC23a_preInclude.py:10

◆ operator<<() [1/2]

template<typename T>
std::ostream & operator<< ( std::ostream & s,
const std::vector< T * > & v )

Definition at line 209 of file rmain.cxx.

209 {
210 for ( size_t i=0 ; i<v.size() ; i++ ) s << "\t" << *v[i] << "\n";
211 return s;
212}

◆ operator<<() [2/2]

template<typename T>
std::ostream & operator<< ( std::ostream & s,
const std::vector< T > & v )

Definition at line 216 of file rmain.cxx.

216 {
217 if ( v.size()<5 ) for ( unsigned i=0 ; i<v.size() ; i++ ) s << "\t" << v[i];
218 else for ( unsigned i=0 ; i<v.size() ; i++ ) s << "\n\t" << v[i];
219 return s;
220}

◆ pointers()

template<typename T>
std::vector< T * > pointers ( std::vector< T > & v)

this is slow - all this copying

Definition at line 367 of file rmain.cxx.

367 {
369 std::vector<T*> vp(v.size(),0);
370 for ( unsigned i=v.size() ; i-- ; ) vp[i] = &v[i];
371 return vp;
372}

◆ replaceauthor()

const std::vector< TIDA::Track > replaceauthor ( const std::vector< TIDA::Track > & tv,
int a0 = 5,
int a1 = 4 )

Definition at line 273 of file rmain.cxx.

273 {
274
275 if ( a0==a1 ) return tv;
276
277 std::vector<TIDA::Track> tr;
278 tr.reserve(tv.size());
279
280 for ( size_t i=0 ; i<tv.size() ; i++ ) {
281 const TIDA::Track& t = tv[i];
282 int a = t.author();
283 if ( a==a0 ) a=a1;
284 tr.push_back( TIDA::Track( t.eta(), t.phi(), t.z0(), t.a0(), t.pT(), t.chi2(), t.dof(),
285 t.deta(), t.dphi(), t.dz0(), t.da0(), t.dpT(),
286 t.bLayerHits(), t.pixelHits(), t.sctHits(), t.siHits(),
287 t.strawHits(), t.trHits(),
288 t.hitPattern(), t.multiPattern(), a, t.hasTruth(),
289 t.barcode(), t.match_barcode(), t.expectBL(), t.id() ) );
290
291 }
292
293 return tr;
294
295}
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
read_hist_ntuple.t
t
Definition read_hist_ntuple.py:5

◆ SelectObjectET()

bool SelectObjectET ( const TrackTrigObject & t)

Definition at line 377 of file rmain.cxx.

377{ return std::fabs(t.pt())>ETmin; }

◆ SelectObjectETovPT()

bool SelectObjectETovPT ( const TrackTrigObject & tobj,
TIDA::Track * t = 0 )

Definition at line 388 of file rmain.cxx.

388 {
389 bool ETselection = std::fabs(tobj.pt())>=ETmin;
390 bool ETovPTselection = true;
391 if ( t ) ETovPTselection = std::fabs(tobj.pt()/t->pT())>=ETovPTmin;
392 return ETselection&ETovPTselection;
393}

◆ time_str()

std::string time_str ( )

Definition at line 116 of file rmain.cxx.

116 {
117 time_t t;
118 time(&t);
119 std::string s(ctime(&t));
120 // std::string::size_type pos = s.find("\n");
121 // if ( pos != std::string::npos )
122 return s.substr(0,s.find('\n'));
123 // return s;
124}
CaloSwCorrections.time
time(flags, cells_name, *args, **kw)
Definition CaloSwCorrections.py:242

◆ trackPtGrtr()

bool trackPtGrtr ( TIDA::Track * trackA,
TIDA::Track * trackB )

Definition at line 205 of file rmain.cxx.

205{ return ( std::fabs(trackA->pT()) > std::fabs(trackB->pT()) ); }

◆ usage()

int usage ( const std::string & name,
int status )

Definition at line 332 of file rmain.cxx.

332 {
333 std::ostream& s = std::cout;
334
335 s << "Usage: " << name << " <config filename> [OPTIONS]" << std::endl;
336 s << "\nOptions: \n";
337 s << " -o, -f, --file value\toutput filename, \n";
338 s << " -b, --binConfig value\tconfig file for histogram configuration, \n";
339 s << " -r, --refChain value\treference chain. + separated keys will be merged, \n";
340 s << " -t, --testChain value\ttest chain, \n";
341 s << " -p, --pdgId value\tpdg ID of truth particle if requiring truth particle processing,\n";
342 s << " --vt value\tuse value as the test vertex selector - overrides value in the config file,\n";
343 s << " --vr value\tuse value as the reference vertex selector - overrides value in the config file,\n";
344 s << " -n, --nofit \ttest do not fit resplots, \n";
345 s << " --rms \ttest force new rms95 errors, \n";
346 s << " -e, --entries value\ttest only run over value entries, \n";
347 // s << " -a, --all \tadd all grids (default)\n";
348 s << " -h, --help \tthis help\n";
349 // s << "\nSee " << PACKAGE_URL << " for more details\n";
350 s << "\nReport bugs to sutt@cern.ch";
351 s << std::endl;
352
353 return status;
354}
merge.status
status
Definition merge.py:16

Variable Documentation

◆ bjetBinConfig

BinConfig bjetBinConfig
extern

◆ cosmicBinConfig

BinConfig cosmicBinConfig
extern

◆ debugPrintout

bool debugPrintout = false

Definition at line 95 of file rmain.cxx.

◆ electronBinConfig

BinConfig electronBinConfig
extern

◆ ETmin

double ETmin = 0

Definition at line 375 of file rmain.cxx.

◆ ETovPTmin

double ETovPTmin = 0

this is a swiss knife function - by default if ET/PT > 0 such that fabs(ET/PT) > 0 is always true and we can always call this function - PT can never be 0 for a particle in the detector so can use this for both ET/PT selection and raw ET selection

Definition at line 386 of file rmain.cxx.

◆ g_binConfig

BinConfig g_binConfig
extern

◆ muonBinConfig

BinConfig muonBinConfig
extern

◆ PRINT_BRESIDUALS

bool PRINT_BRESIDUALS
extern

stack trace headers

globals for communicating with *Analyses TagNProbe class

Definition at line 29 of file ConfAnalysis.cxx.

◆ tauBinConfig

BinConfig tauBinConfig
extern
Generated on for ATLAS Offline Software by doxygen 1.14.0