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...
true by default - command line options sets to false so should override what is in the config file
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
new code - can extract vtx name, pt, any extra options that we want, but also chop off everythiung after :post
track selectors so we can select multiple times with different filters if we want (simpler then looping over vectors each time
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 ...
Determine what sort of matching is required ...
I don't prefer range based for loops ...
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
leave this in for the moment ...
both legs have the same reference ...
different references for each leg ...
tag leg reference ...
do we want to filter on the RoI properties? If so, if the RoI fails the cuts, then skip this roi
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 ...
check configuration is provided ...
if requesting an object match, remove any tracks which correspond to an object below the object PT threshold
do for all vertices now ...
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 ...
532 std::cerr <<
"Error: no config file specified\n" << std::endl;
546 std::cout <<
"$0 :: compiled " << __DATE__ <<
" " << __TIME__ << std::endl;
550 std::string histofilename(
"");
552 std::string datafile =
"";
554 std::vector<std::string> refChains(0);
556 std::vector<TrackFilter*> refFilters(0);
561 std::string refChain =
"";
565 std::vector<std::string> testChains;
567 std::string binningConfigFile =
"";
569 bool useoldrms =
true;
572 std::string vertexSelection =
"";
573 std::string vertexSelection_rec =
"";
575 unsigned Nentries = 0;
577 for (
int i=1 ;
i<
argc ;
i++ ) {
578 if ( std::string(
argv[
i])==
"-h" || std::string(
argv[
i])==
"--help" ) {
581 else if ( std::string(
argv[
i])==
"-e" || std::string(
argv[
i])==
"--entries" ) {
585 else if ( std::string(
argv[
i])==
"-o" || std::string(
argv[
i])==
"-f" || std::string(
argv[
i])==
"--file" ) {
587 histofilename =
argv[
i];
588 if ( histofilename.find(
".root")==std::string::npos ) histofilename +=
".root";
590 else if ( std::string(
argv[
i])==
"-r" || std::string(
argv[
i])==
"--refChain" ) {
595 if (refChain.find(
"+") != string::npos){
596 std::istringstream iss(refChain);
598 while (std::getline(iss, token,
'+')){
599 refChains.push_back(token);
600 refFilters.push_back(0);
604 refChains.push_back(
argv[
i]);
605 refFilters.push_back(0);
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" ) {
612 testChains.push_back(
argv[
i]);
614 else if ( std::string(
argv[
i])==
"-p" || std::string(
argv[
i])==
"--pdgId" ) {
618 else if ( std::string(
argv[
i])==
"--vr" ) {
620 vertexSelection =
argv[
i];
622 else if ( std::string(
argv[
i])==
"--vt" ) {
624 vertexSelection_rec =
argv[
i];
626 else if ( std::string(
argv[
i])==
"-b" || std::string(
argv[
i])==
"--binConfig" ) {
628 binningConfigFile = std::string(
argv[
i]);
630 else if ( std::string(
argv[
i]).
find(
'-')==0 ) {
632 std::cerr <<
"unknown option " <<
argv[
i] << std::endl;
642 if ( datafile==
"" ) {
643 std::cerr <<
"no config file specifed\n" << endl;
647 std::cout <<
time_str() << std::endl;
656 bool oldrms95 =
true;
657 inputdata.declareProperty(
"OldRMS95", oldrms95 );
658 std::cout <<
"setting Resplot old rms95 " << oldrms95 << std::endl;
662 std::cout <<
"setting Resplot old rms95 " << useoldrms << std::endl;
668 std::cout <<
"Not fitting resplots " << std::endl;
674 inputdata.declareProperty(
"NFiles", nfiles );
677 if ( histofilename==
"" ){
678 if (
inputdata.isTagDefined(
"outputFile") ) histofilename =
inputdata.GetString(
"outputFile");
680 std::cerr <<
"Error: no output file defined\n" << std::endl;
686 TFile foutput( histofilename.c_str(),
"recreate" );
687 if (!foutput.IsOpen()) {
688 std::cerr <<
"Error: could not open output file\n" << std::endl;
694 std::cout <<
"writing output to " << histofilename << std::endl;
696 TH1D* hevent =
new TH1D(
"event",
"event", 1000, 10000, 80000 );
697 Resplot* hcorr =
new Resplot(
"correlation", 21, -0.5, 20.5, 75, 0, 600);
715 bool expectBL =
false;
729 double massMax = 130;
734 bool rotate_testtracks =
false;
736 if (
inputdata.isTagDefined(
"RotateTestTracks") ) rotate_testtracks = (
inputdata.GetValue(
"RotateTestTracks") ? true : false );
738 bool truthMatch =
false;
740 if (
inputdata.isTagDefined(
"TruthMatch") ) truthMatch = (
inputdata.GetValue(
"TruthMatch") ? true : false );
745 std::vector<size_t> refPtOrd_indices;
746 bool use_pt_ordered_ref =
false;
748 if (
inputdata.isTagDefined(
"UsePtOrderedRefTracks") ) {
749 std::vector<double> refPtOrd_indices_tmp;
751 use_pt_ordered_ref =
true;
752 refPtOrd_indices_tmp = (
inputdata.GetVector(
"UsePtOrderedRefTracks") );
754 std::cout <<
"using PT ordered reference tracks: " << refPtOrd_indices_tmp << std::endl;
756 for (
size_t i=refPtOrd_indices_tmp.size();
i-- ; ) {
757 refPtOrd_indices.push_back( refPtOrd_indices_tmp.at(
i) );
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 );
792 if ( pdgId==0 &&
inputdata.isTagDefined(
"pdgId") ) pdgId =
inputdata.GetValue(
"pdgId");
794 if (
inputdata.isTagDefined(
"InvMassMax") ) massMax =
inputdata.GetValue(
"InvMassMax");
795 if (
inputdata.isTagDefined(
"InvMassMin") ) massMin =
inputdata.GetValue(
"InvMassMin");
808 std::string useMatcher =
"DeltaR";
809 if (
inputdata.isTagDefined(
"UseMatcher") ) useMatcher =
inputdata.GetString(
"UseMatcher");
813 if ( refChain==
"" ) {
814 if (
inputdata.isTagDefined(
"refChain") ) {
815 refChain =
inputdata.GetString(
"refChain");
816 refChains.push_back(refChain);
817 refFilters.push_back(0);
820 std::cerr <<
"Error: no reference chain defined\n" << std::endl;
830 if (
inputdata.isTagDefined(
"Exclude") ) {
835 if (refChains.size() == 0){
836 std::cerr <<
"Error: refChains is empty\n" <<std::endl;
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") );
849 std::vector<ChainString> chainConfig;
851 for (
size_t ic=0 ;
ic<testChains.size() ;
ic++ ) {
853 testChains[
ic] = chainConfig.back().pre();
864 bool select_roi =
true;
870 bool use_custom_ref_roi =
false;
871 const int custRefRoi_nParams = 3;
873 double custRefRoi_params[custRefRoi_nParams] = {-999., -999., -999.};
874 std::vector<std::string> custRefRoi_chainList;
875 std::set<std::string> customRoi_chains;
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");
881 if (
inputdata.isTagDefined(
"customRefRoi_chainList") ) custRefRoi_chainList =
inputdata.GetStringVector(
"customRefRoi_chainList");
883 for (
unsigned ic=0 ;
ic<custRefRoi_chainList.size() ;
ic++ ) customRoi_chains.insert( custRefRoi_chainList[
ic] );
885 for (
int param_idx=0 ; param_idx<custRefRoi_nParams ; param_idx++ ) {
886 if ( custRefRoi_params[param_idx] != -999 ) {
888 use_custom_ref_roi =
true;
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;
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;
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;
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;
905 if ( !custRefRoi_chainList.empty() ) {
906 std::cout <<
"**** \t****" << std::endl;
907 std::cout <<
"**** Applying custom RoI only to specified chains\t****" << std::endl;
909 std::cout <<
"**** \t****" << std::endl;
914 if (
inputdata.isTagDefined(
"SelectRoi") ) {
915 select_roi = (
inputdata.GetValue(
"SelectRoi")!=0 ? true : false );
919 std::cout <<
"**** ****" << std::endl;
920 std::cout <<
"**** RoI filtering of reference tracks is disabled ****" << std::endl;
921 std::cout <<
"**** ****" << std::endl;
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] );
941 else if (
inputdata.isTagDefined(
"LumiBlocks") ) {
954 if ( vertexSelection ==
"" ) {
955 if (
inputdata.isTagDefined(
"VertexSelection") ) vertexSelection =
inputdata.GetString(
"VertexSelection");
958 bool bestPTVtx =
false;
959 bool bestPT2Vtx =
false;
962 if ( vertexSelection!=
"" ) {
963 if ( vertexSelection==
"BestPT" ) bestPTVtx =
true;
964 else if ( vertexSelection==
"BestPT2" ) bestPT2Vtx =
true;
973 if ( vertexSelection_rec ==
"" ) {
974 if (
inputdata.isTagDefined(
"VertexSelectionRec") ) vertexSelection_rec =
inputdata.GetString(
"VertexSelectionRec");
977 bool bestPTVtx_rec =
false;
978 bool bestPT2Vtx_rec =
false;
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 );
987 std::cout <<
"vertexSelection: " << vertexSelection << std::endl;
988 std::cout <<
"vertexSelection_rec: " << vertexSelection_rec << std::endl;
996 bool useBestVertex =
false;
997 if (
inputdata.isTagDefined(
"useBestVertex") ) useBestVertex = (
inputdata.GetValue(
"useBestVertex") ? 1 : 0 );
999 bool useSumPtVertex =
true;
1000 if (
inputdata.isTagDefined(
"useSumPtVertex") ) useSumPtVertex = (
inputdata.GetValue(
"useSumPtVertex") ? 1 : 0 );
1002 int MinVertices = 1;
1003 if (
inputdata.isTagDefined(
"MinVertices") ) MinVertices =
inputdata.GetValue(
"MinVertices");
1010 bool useVertexTracks =
false;
1011 if (
inputdata.isTagDefined(
"UseVertexTracks") ) useVertexTracks = (
inputdata.GetValue(
"UseVertexTracks") > 0 );
1014 std::string vertex_refname =
"Vertex";
1015 if (
inputdata.isTagDefined(
"VertexReference") ) vertex_refname +=
":" +
inputdata.GetString(
"VertexReference");
1018 int NVtxTrackCut = 2;
1019 if (
inputdata.isTagDefined(
"NVtxTrackCut") ) NVtxTrackCut =
inputdata.GetValue(
"NVtxTrackCut");
1023 bool event_selector_flag =
false;
1029 if ( es.size() ) event_selector_flag =
true;
1032 std::vector<double> beamTest;
1033 std::vector<double> beamRef;
1036 bool correctBeamlineRef =
false;
1037 bool correctBeamlineTest =
false;
1039 if (
inputdata.isTagDefined(
"CorrectBeamlineRef") ) correctBeamlineRef = (
inputdata.GetValue(
"CorrectBeamlineRef") == 0 ? false : true );
1040 if (
inputdata.isTagDefined(
"CorrectBeamlineTest") ) correctBeamlineTest = (
inputdata.GetValue(
"CorrectBeamlineTest") == 0 ? false : true );
1043 if (
inputdata.isTagDefined(
"BeamTest") ) beamTest =
inputdata.GetVector(
"BeamTest");
1045 if (
inputdata.isTagDefined(
"BeamTestx") ) beamTest.push_back(
inputdata.GetValue(
"BeamTestx"));
1046 if (
inputdata.isTagDefined(
"BeamTesty") ) beamTest.push_back(
inputdata.GetValue(
"BeamTesty"));
1052 if (
inputdata.isTagDefined(
"BeamRefx") ) beamRef.push_back(
inputdata.GetValue(
"BeamRefx"));
1053 if (
inputdata.isTagDefined(
"BeamRefy") ) beamRef.push_back(
inputdata.GetValue(
"BeamRefy"));
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;
1064 if ( beamTest.size()>0 ) correctBeamlineTest =
true;
1065 if ( beamRef.size()>0 ) correctBeamlineRef =
true;
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;
1072 if ( beamRef.size()>0 ) std::cout <<
"beamref " << beamRef << std::endl;
1073 if ( beamTest.size()>0 ) std::cout <<
"beamtest " << beamTest << std::endl;
1082 double a0vrec = 1000;
1083 double z0vrec = 2000;
1088 bool initialiseFirstEvent =
false;
1089 if (
inputdata.isTagDefined(
"InitialiseFirstEvent") ) initialiseFirstEvent =
inputdata.GetValue(
"InitialiseFirstEvent");
1096 bool doPurity =
false;
1097 if (
inputdata.isTagDefined(
"doPurity") ) doPurity = (
inputdata.GetValue(
"doPurity")==0 ? false : true );
1102 bool monitorZBeam =
false;
1103 if (
inputdata.isTagDefined(
"MonitorinZBeam") ) monitorZBeam = (
inputdata.GetValue(
"MonitorZBeam")==0 ? false : true );
1105 std::cout <<
"dbg " << __LINE__ << std::endl;
1111 if ( binningConfigFile!=
"" ) binningConfig =
new ReadCards( binningConfigFile );
1123 if ( binningConfig!=&
inputdata )
delete binningConfig;
1128 int selectcharge = 0;
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;
1147 std::cout <<
"a0v: " << a0v << std::endl;
1148 std::cout <<
"z0v: " << z0v << std::endl;
1153 npix, nsct, -1, nbl,
1155 npixholes, nsctholes, nsiholes, expectBL );
1157 if ( selectcharge!=0 ) filter_offline.chargeSelection( selectcharge );
1158 if ( pTMax>
pT ) filter_offline.maxpT( pTMax );
1185 Filter_Track filter_onlinekine( eta_rec, 1000, 0., 2000,
pT, -1, npix, nsct, -1, -2, -2);
1187 Filter_Combined filter_online( &filter_onlinekine, &filter_onlinevertex );
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 );
1199 if (
inputdata.isTagDefined(
"Filter" ) ) {
1201 std::cout <<
"Filter: " <<
inputdata.GetString(
"Filter") <<
" : " <<
filter << std::endl;
1202 if (
filter.head()==
"Offline2017" ) {
1203 std::string filter_type =
filter.tail();
1205 filter_offCP =
new Filter_Combined ( filter_offline2017, &filter_vertex);
1206 refFilter = filter_offCP;
1208 else if (
filter.head()==
"OfflineR22" ) {
1209 std::string filter_type =
filter.tail();
1211 filter_offCP =
new Filter_Combined ( filter_offlineR22, &filter_vertex);
1212 refFilter = filter_offCP;
1215 std::cerr <<
"unimplemented Filter requested: " <<
filter.head() << std::endl;
1220 if ( !(refFilter =
getFilter( refChains[0], pdgId, &filter_off, &filter_muon, &filter_truth ) ) ) {
1221 std::cerr <<
"unknown reference chain defined" << std::endl;
1226 refFilters.push_back(refFilter);
1229 std::map<std::string,TIDA::Reference>
ref;
1231 std::vector<NtupleTrackSelector*> refSelectors;
1238 for (
size_t ic=0 ;
ic<refChains.size() ;
ic++ ) {
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;
1245 refFilter = refFilters[
ic];
1247 else refFilters[
ic] = refFilter;
1255 if (pdgId==0) truthFilter = &filter_off;
1256 else truthFilter = &filter_truth;
1262 std::cout <<
"filter_passthrough" << std::endl;
1264 Filter_Track filter_passthrough( 10, 1000, 0., 2000, pT_rec, npix_rec, nsct_rec, 1, -2, -2, -2);
1269 std::cout <<
"using tracks: " << refChain <<
" for reference sample" << std::endl;
1275 std::vector<std::string>& test_chains = testChains;
1280 std::map<std::string,TrackAnalysis*> analysis;
1286 std::vector<TrackAnalysis*>
analyses;
1287 analyses.reserve(test_chains.size());
1289 std::cout <<
"booking " << test_chains.size() <<
" analyses" << std::endl;
1291 for (
unsigned i=0 ;
i<test_chains.size() ;
i++ ) {
1295 std::vector<std::string> chainnames;
1297 chainnames.push_back(chainname);
1305 std::string probe_extra = probe.
extra();
1307 if ( probe_extra.empty() ) probe_extra = probe.
postvalue(
"extra");
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;
1315 size_t p = probe_extra.find(
"_probe");
1317 if (
p!=std::string::npos ) {
1319 std::string probe_ref = refChains[0];
1321 if ( !probe.
postvalue(
"ref").empty() ) {
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;
1334 std::string probe_key = probe_extra.erase(
p, 6);
1336 for (
unsigned j=0 ; j<test_chains.size(); ++j) {
1338 if (
i==j )
continue;
1342 if (
tag.head() != probe.
head() )
continue;
1346 std::string tag_extra =
tag.extra();
1348 if ( tag_extra.empty() ) tag_extra =
tag.postvalue(
"extra");
1349 if ( tag_extra.find(
"_tag")==std::string::npos )
continue;
1352 std::string tag_key = tag_extra.erase( tag_extra.find(
"_tag"), 4) ;
1355 if ( tag_key != probe_key )
continue;
1357 if (
tag.element() == probe.
element() )
continue;
1359 std::string tag_ref = refChains[0];
1367 if ( !
tag.postvalue(
"ref").empty() ) {
1369 tag_ref =
tag.postvalue(
"ref");
1372 std::cout <<
"tag ref: " << tag_ref << std::endl;
1374 if (
ref.find(tag_ref)==
ref.end() ) {
1385 TnP_tool =
new TagNProbe( tag_ref, probe_ref, massMin, massMax);
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;
1402 if ( TnP_tool )
replace( chainname,
"/",
"_" );
1409 std::string vtxTool = chainConfig[
i].postvalue(
"rvtx");
1411 if ( vtxTool!=
"" ) {
1416 analy_conf->
store().
insert( anal_confvtx,
"rvtx" );
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;
1442 if ( analysis.find( chainname )==analysis.end() ) {
1443 analysis.insert( std::map<std::string,TrackAnalysis*>::value_type( chainname, analy_conf ) );
1447 std::cerr <<
"WARNING: Duplicated chain"
1449 <<
"---------------------------------"
1450 <<
"---------------------------------"
1451 <<
"---------------------------------" << std::endl;
1455 std::cout <<
"analysis: " << chainname <<
"\t" << analy_conf
1457 <<
"---------------------------------"
1458 <<
"---------------------------------"
1459 <<
"---------------------------------" << std::endl;
1463 std::cout <<
"purity " << (chainnames[0]+
"-purity") << std::endl;
1467 analysis[chainnames[0]+
"-purity"] = analp;
1473 std::cout <<
"main() finished looping" << std::endl;
1477 bool fullyContainTracks =
false;
1479 if (
inputdata.isTagDefined(
"FullyContainTracks") ) {
1480 fullyContainTracks = (
inputdata.GetValue(
"FullyContainTracks")==0 ? false : true );
1483 bool containTracksPhi =
true;
1485 if (
inputdata.isTagDefined(
"ContainTracksPhi") ) {
1486 containTracksPhi = (
inputdata.GetValue(
"ContainTracksPhi")==0 ? false : true );
1491 dynamic_cast<Filter_Combined*
>(refFilter)->containtracks(fullyContainTracks);
1492 dynamic_cast<Filter_Combined*
>(refFilter)->containtracksPhi(containTracksPhi);
1509 bool filterRoi =
false;
1512 bool roicomposite =
false;
1515 if (
inputdata.isTagDefined(
"FilterRoi" ) ) {
1519 std::vector<double> filter_values =
inputdata.GetVector(
"FilterRoi" );
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);
1527 RoiFilter roiFilter( roieta, roicomposite, roimult );
1534 else if ( useMatcher ==
"DeltaRZ" || useMatcher ==
"DeltaRZSinTheta" ) {
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");
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");
1550 else if ( useMatcher ==
"Truth" ) {
1563 else if ( useMatcher ==
"DeltaRZ" || useMatcher ==
"DeltaRZSinTheta" ) {
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");
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");
1579 else if ( useMatcher ==
"Truth" ) {
1603 if (
inputdata.isTagDefined(
"DataSets") ) {
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;
1610 std::vector<std::string> filenames_ =
d.datafiles();
1611 std::cout <<
"\tdataset contains " << filenames_.size() <<
" files" << std::endl;
1617 std::cerr <<
"no input data specified" << std::endl;
1627 bool show_release =
true;
1629 std::vector<std::string> release_data;
1631 std::string release_data_save =
"";
1633 if ( show_release ){
1639 TFile* finput = TFile::Open(
filenames[
i].c_str() );
1642 if ( finput==0 || !finput->IsOpen() || finput->IsZombie() ) {
1643 std::cerr <<
"Error: could not open input file: " <<
filenames[
i] << std::endl;
1647 TTree* dataTree = (TTree*)finput->Get(
"dataTree");
1648 TString* releaseData =
new TString(
"");
1651 dataTree->SetBranchAddress(
"ReleaseMetaData", &releaseData);
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;
1660 release_data_save = release_data.back();
1664 if ( finput )
delete finput;
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());
1677 TTree* dataTree =
new TTree(
"dataTree",
"dataTree");
1678 TString* releaseData =
new TString(
"");
1681 dataTree->Branch(
"ReleaseMetaData",
"TString", &releaseData);
1683 for (
unsigned ird=0 ; ird<release_data.size() ; ird++ ) {
1684 *releaseData = release_data[ird];
1688 dataTree->Write(
"", TObject::kOverwrite);
1706 if ( Nentries==0 &&
inputdata.isTagDefined(
"Nentries") ) {
1712 typedef std::pair<int,double> zpair;
1713 std::vector<zpair> refz;
1714 std::vector<zpair> testz;
1716 std::vector<double> beamline_ref;
1717 std::vector<double> beamline_test;
1722 std::cout <<
"opening files" << std::endl;
1726 int pregrl_events = 0;
1727 int grl_counter = 0;
1730 std::cout <<
"starting event loop " <<
time_str() << std::endl;
1734 if ( nfiles!=0 && nfiles<max_files ) max_files = nfiles;
1743 if ( finput==0 || !finput->IsOpen() || finput->IsZombie() ) {
1744 std::cerr <<
"Error: could not open output file " <<
filenames[
ifile] << std::endl;
1748 TTree*
data = (TTree*)finput->Get(
"tree");
1751 std::cerr <<
"Error: cannot open TTree: " <<
filenames[
ifile] << std::endl;
1759 data->SetBranchAddress(
"TIDA::Event",&track_ev);
1765 unsigned cNentries =
data->GetEntries();
1771 for (
unsigned int i=0;
skip &&
run &&
i<cNentries ;
i++ ) {
1793 if ( !ingrl )
continue;
1799 if ( event_selector_flag && !es.in(
event ) )
continue;
1801 if ( mintime>
ts ) mintime =
ts;
1802 if ( maxtime<
ts ) maxtime =
ts;
1806 std::cout <<
"breaking out " <<
run << std::endl;
1816 hevent->Fill(
event );
1819 if ( (cNentries<10) ||
i%(cNentries/10)==0 ||
i%1000==0 ||
debugPrintout ) {
1820 std::cout <<
"run " << track_ev->
run_number()
1823 <<
"\tchains " << track_ev->
chains().size()
1825 std::cout <<
"\t : processed " <<
i <<
" events so far (" <<
int((1000*
i)/cNentries)*0.1 <<
"%)\t" <<
time_str() << std::endl;
1832 if ( nfiles>0 ) pfiles = nfiles;
1835 std::cout <<
"file entries=" <<
data->GetEntries();
1837 if (
data->GetEntries()<100 ) std::cout <<
" ";
1838 if (
data->GetEntries()<1000 ) std::cout <<
" ";
1839 if (
data->GetEntries()<10000 ) std::cout <<
" ";
1844 std::cout <<
"run " << track_ev->
run_number()
1847 <<
"\tchains " << track_ev->
chains().size()
1850 std::cout <<
"\t : processed " <<
ifile <<
" files so far (" <<
int((1
e3*
ifile)/pfiles)*0.1 <<
"%)\t" <<
time_str() << std::endl;
1863 truthTracks.clear();
1864 refPurityTracks.clear();
1868 mit->second.selector()->clear();
1874 const std::vector<TIDA::Chain>&
chains = track_ev->
chains();
1881 truthTracks.selectTracks(
chains[
ic][0].tracks() );
1888 for (
const std::string& rc : refChains){
1891 offTracks.selectTracks(
chains[
ic][0].tracks() );
1893 beamline_ref =
chains[
ic][0].user();
1902 std::vector<TIDA::Vertex> vertices;
1906 if ( vtxchain && vtxchain->
size()>0 ) {
1908 const std::vector<TIDA::Vertex>& mv = vtxchain->
at(0).
vertices();
1913 if (
debugPrintout ) std::cout <<
"vertices:\n" << mv << std::endl;
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;
1920 vtx_temp.selectTracks( offTracks.tracks() );
1921 for (
unsigned itr=0; itr<vtx_temp.tracks().
size(); itr++) {
1923 if ( bestPTVtx ) selection_ += std::fabs(tr->
pT());
1924 else if ( bestPT2Vtx ) selection_ += std::fabs(tr->
pT())*std::fabs(tr->
pT());
1931 if ( selectvtx!=-1 ) {
1932 vertices.push_back( mv[selectvtx] );
1935 else if ( vtxind>=0 ) {
1936 if (
size_t(vtxind)<mv.size() ) {
1937 vertices.push_back( mv[vtxind] );
1941 for (
size_t iv=0 ; iv<mv.size() ; iv++ ) {
1942 vertices.push_back( mv[iv] );
1947 filter_vertex.setVertex( vertices );
1953 for (
unsigned iv=0 ; iv<mv.size() ; iv++ ) {
1954 int Ntracks = mv[iv].Ntracks();
1955 if ( Ntracks>NVtxTrackCut ) {
1972 bool foundReference =
false;
1981 for (
const std::string& rc : refChains ) {
1990 std::string refname = mitr->first;
1999 *mitr->second.tom() = rtom;
2013 if ( !foundReference )
continue;
2016 std::cout <<
"reference chain:\n" << *refchain << std::endl;
2030 if ( analitr==analysis.end() )
continue;
2033 std::cout <<
"test chain:\n" <<
chain << std::endl;
2038 std::vector<TIDA::Roi*>
rois;
2058 if ( TnP_tool->
type0() == TnP_tool->
type1() ) {
2060 rois = TnP_tool->
GetRois( track_ev->
chains(), selector0, filter0, invmass, invmass_obj, rtom );
2074 for (
size_t ia=0 ; ia<selector1->
tracks().size() ; ia++ ) {
2076 std::cout << ia <<
"\t" << *
track << std::endl;
2077 std::cout <<
"\t" << rtom1->
object(
track->id()) << std::endl;
2080 std::cout << *TnP_tool << std::endl;
2083 rois = TnP_tool->
GetRois( track_ev->
chains(), selector0, filter0, selector1, filter1, invmass, invmass_obj, rtom, rtom1 );
2093 if (
false ) std::cout <<
"++++++++++++ rois size = " <<
rois.size() <<
" +++++++++++" << std::endl;
2096 for (
unsigned int ir=0 ;
ir<
rois.size() ;
ir++ ) {
2106 testTracks.selectTracks( troi.
tracks() );
2109 std::vector<TIDA::Track*> testp = testTracks.tracks();
2114 if ( filterRoi && !roiFilter.filter( roi ) )
continue;
2117 const std::vector<TIDA::Vertex>& mvt = troi.
vertices();
2119 std::vector<TIDA::Vertex> vertices_test;
2124 if ( bestPTVtx_rec || bestPT2Vtx_rec ) {
2128 for (
unsigned iv=0 ; iv<mvt.size() ; iv++ ) {
2129 double selection_ = 0.0;
2131 vtx_temp.selectTracks( testp );
2132 for (
unsigned itr=0; itr<vtx_temp.tracks().
size(); itr++) {
2134 if ( bestPTVtx ) selection_ += std::fabs(tr->
pT());
2135 else if ( bestPT2Vtx ) selection_ += std::fabs(tr->
pT())*std::fabs(tr->
pT());
2142 if ( selectvtx!=-1 ) {
2144 if ( useVertexTracks ) selected.selectTracks( testp );
2145 vertices_test.push_back(selected);
2148 else if ( vtxind_rec!=-1 ) {
2149 if (
unsigned(vtxind_rec)<mvt.size() ) {
2151 if ( useVertexTracks ) selected.selectTracks( testp );
2152 vertices_test.push_back( selected );
2156 for (
unsigned iv=0 ; iv<mvt.size() ; iv++ ) {
2158 if ( useVertexTracks ) selected.selectTracks( testp );
2159 vertices_test.push_back( selected );
2164 beamline_test =
rois[
ir]->user();
2171 if ( correctBeamlineTest ) {
2172 if ( beamTest.size()==2 ) analy_track->
setBeamTest( beamTest[0], beamTest[1] );
2175 if ( !
inputdata.isTagDefined(
"BeamTest") ) {
2176 if ( beamline_test.size()==2 ) analy_track->
setBeamTest( beamline_test[0], beamline_test[1] );
2182 if ( correctBeamlineRef ) {
2183 if ( beamRef.size()==2 ) analy_track->
setBeamRef( beamRef[0], beamRef[1] );
2186 if ( !
inputdata.isTagDefined(
"BeamRef") ) {
2187 if ( beamline_ref.size()==2 ) analy_track->
setBeamRef( beamline_ref[0], beamline_ref[1] );
2206 bool customRefRoi_thisChain =
false;
2208 if ( use_custom_ref_roi ) {
2209 if ( customRoi_chains.size() ) {
2210 if ( customRoi_chains.find(
chain.name() )!=customRoi_chains.end() ) customRefRoi_thisChain =
true;
2212 else customRefRoi_thisChain =
true;
2215 if ( use_custom_ref_roi && customRefRoi_thisChain ) {
2216 refRoi =
makeCustomRefRoi( roi, custRefRoi_params[0], custRefRoi_params[1], custRefRoi_params[2] );
2232 std::vector<TIDA::Track*> refp_vec = refTracks.tracks( refFilter );
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;
2242 truth_matcher->
match( refp_tmp, truth );
2244 std::vector<TIDA::Track*> refp_matched;
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] );
2252 refp_vec = refp_matched;
2256 if ( use_pt_ordered_ref ) {
2259 size_t nRefTracks = refp_vec.size();
2261 std::vector<TIDA::Track*> refp_chosenPtOrd;
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;
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));
2285 refp_vec = refp_chosenPtOrd;
2296 if ( ptconfig!=
"" ) {
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 );
2310 if ( d0config!=
"" ) {
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 );
2332 if ( objectET !=
"" ) {
2334 double ETconfig =
std::atof( objectET.c_str() );
2340 while ( itr!=refp_vec.end() ) {
2343 if ( tobj==0 || tobj->
pt()<ETconfig )
2344 itr=refp_vec.erase( itr );
2352 const std::vector<TIDA::Track*>& refp = refp_vec;
2369 std::vector<TIDA::Vertex> vertices_roi;
2377 vertices_roi.clear();
2379 const std::vector<TIDA::Vertex>& mv = vertices;
2384 for (
unsigned iv=0 ; iv<mv.size() ; iv++ ) {
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;
2397 if ( roi.zedMinus()<=vx.
z() && roi.zedPlus()>=vx.
z() ) accept_vertex =
true;
2400 if ( !accept_vertex )
continue;
2406 if ( useVertexTracks ) {
2409 vertex_roi.selectTracks( refp );
2410 trackcount = vertex_roi.Ntracks();
2411 if ( trackcount>=ntracks && trackcount>0 ) {
2412 vertices_roi.push_back( vertex_roi );
2417 for (
unsigned itr=0; itr<refp.size(); itr++){
2420 double dzsintheta = std::fabs( (vx.
z()-tr->
z0()) *
std::sin(theta_) );
2421 if( dzsintheta < 1.5 ) trackcount++;
2426 if ( trackcount>=ntracks && trackcount>0 ) {
2428 vertices_roi.push_back( vertex_roi );
2436 if ( rotate_testtracks )
for (
size_t i=testp.size() ;
i-- ; ) testp[
i]->rotate();
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]) );
2449 matcher->
match( refp, testp);
2451 if ( tom.
status() ) analitr->second->execute( refp, testp, matcher, &tom );
2452 else analitr->second->execute( refp, testp, matcher );
2455 analitr->second->store().find( vtxanal,
"rvtx" );
2459 if ( vtxanal && ( refp.size() >=
size_t(ntracks) ) ) {
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;
2482 TrackAssociator::map_type::const_iterator titr = matcher->TrackAssociator::matched().begin();
2483 TrackAssociator::map_type::const_iterator
tend = matcher->TrackAssociator::matched().end();
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;
2492 std::cout <<
"completed : " <<
chain.name() <<
"\n";
2493 std::cout <<
"-----------------------------------" << std::endl;
2498 if ( _matcher->size()<refp.size() ) {
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;
2512 const std::vector<TIDA::Track*>& refpp = refPurityTracks.tracks( refFilter );
2514 testPurityTracks.clear();
2516 testPurityTracks.selectTracks( troi.
tracks() );
2517 std::vector<TIDA::Track*> testpp = testPurityTracks.tracks();
2519 matcher->
match(refpp, testpp);
2524 if ( analitrp == analysis.end() )
continue;
2527 analitrp->second->execute( refpp, testpp, matcher );
2548 <<
"\ttimes " << mintime <<
" " << maxtime
2549 <<
"\t( grl: " << grl_counter <<
" / " << pregrl_events <<
" )" << std::endl;
2551 if ( monitorZBeam )
zbeam zb( refz, testz );
2565 analyses[
i]->store().find( vtxanal,
"rvtx" );
2566 if ( vtxanal ) vtxanal->
finalise();
2577 mit->second.Clean();
2581 if ( ex_matcher )
delete ex_matcher;
2583 std::cout <<
"done" << std::endl;