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
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;
1333 std::string probe_key = probe_extra.erase(
p, 6);
1335 for (
unsigned j=0 ; j<test_chains.size(); ++j) {
1337 if (
i==j )
continue;
1341 if (
tag.head() != probe.
head() )
continue;
1345 std::string tag_extra =
tag.extra();
1347 if ( tag_extra.empty() ) tag_extra =
tag.postvalue(
"extra");
1348 if ( tag_extra.find(
"_tag")==std::string::npos )
continue;
1351 std::string tag_key = tag_extra.erase( tag_extra.find(
"_tag"), 4) ;
1354 if ( tag_key != probe_key )
continue;
1356 if (
tag.element() == probe.
element() )
continue;
1358 std::string tag_ref = refChains[0];
1366 if ( !
tag.postvalue(
"ref").empty() ) {
1368 tag_ref =
tag.postvalue(
"ref");
1371 std::cout <<
"tag ref: " << tag_ref << std::endl;
1373 if (
ref.find(tag_ref)==
ref.end() ) {
1384 TnP_tool =
new TagNProbe( tag_ref, probe_ref, massMin, massMax);
1386 TnP_tool->
probe(probe);
1387 std::cout <<
"Tag and probe pair found! \n\t Tag : " <<
tag <<
"\n\t Probe : " << probe
1388 <<
"\n\t tag ref: " << tag_ref
1389 <<
"\n\tprobe ref: " << probe_ref
1390 <<
"\n-------------------" << std::endl;
1401 if ( TnP_tool )
replace( chainname,
"/",
"_" );
1408 std::string vtxTool = chainConfig[
i].postvalue(
"rvtx");
1410 if ( vtxTool!=
"" ) {
1415 analy_conf->
store().
insert( anal_confvtx,
"rvtx" );
1432 std::cout <<
"chain:: " << chainname <<
" : size " << chainConfig[
i].values().size() << std::endl;
1433 for (
unsigned ik=chainConfig[
i].
values().
size() ; ik-- ; ) {
1434 std::cout <<
"\tchainconfig: " << ik <<
"\tkey " << chainConfig[
i].keys()[ik] <<
" " << chainConfig[
i].values()[ik] << std::endl;
1441 if ( analysis.find( chainname )==analysis.end() ) {
1442 analysis.insert( std::map<std::string,TrackAnalysis*>::value_type( chainname, analy_conf ) );
1446 std::cerr <<
"WARNING: Duplicated chain"
1448 <<
"---------------------------------"
1449 <<
"---------------------------------"
1450 <<
"---------------------------------" << std::endl;
1454 std::cout <<
"analysis: " << chainname <<
"\t" << analy_conf
1456 <<
"---------------------------------"
1457 <<
"---------------------------------"
1458 <<
"---------------------------------" << std::endl;
1462 std::cout <<
"purity " << (chainnames[0]+
"-purity") << std::endl;
1466 analysis[chainnames[0]+
"-purity"] = analp;
1472 std::cout <<
"main() finished looping" << std::endl;
1476 bool fullyContainTracks =
false;
1478 if (
inputdata.isTagDefined(
"FullyContainTracks") ) {
1479 fullyContainTracks = (
inputdata.GetValue(
"FullyContainTracks")==0 ? false : true );
1482 bool containTracksPhi =
true;
1484 if (
inputdata.isTagDefined(
"ContainTracksPhi") ) {
1485 containTracksPhi = (
inputdata.GetValue(
"ContainTracksPhi")==0 ? false : true );
1490 dynamic_cast<Filter_Combined*
>(refFilter)->containtracks(fullyContainTracks);
1491 dynamic_cast<Filter_Combined*
>(refFilter)->containtracksPhi(containTracksPhi);
1508 bool filterRoi =
false;
1511 bool roicomposite =
false;
1514 if (
inputdata.isTagDefined(
"FilterRoi" ) ) {
1518 std::vector<double> filter_values =
inputdata.GetVector(
"FilterRoi" );
1520 if ( filter_values.size()>0 ) roieta = filter_values[0];
1521 if ( filter_values.size()>1 ) roicomposite = ( filter_values[1]==0 ? false : true );
1522 if ( filter_values.size()>2 ) roimult =
int(filter_values[2]+0.5);
1526 RoiFilter roiFilter( roieta, roicomposite, roimult );
1533 else if ( useMatcher ==
"DeltaRZ" || useMatcher ==
"DeltaRZSinTheta" ) {
1537 if (
inputdata.isTagDefined(
"Matcher_deta" ) ) deta =
inputdata.GetValue(
"Matcher_deta");
1538 if (
inputdata.isTagDefined(
"Matcher_dphi" ) ) dphi =
inputdata.GetValue(
"Matcher_dphi");
1539 if (
inputdata.isTagDefined(
"Matcher_dzed" ) ) dzed =
inputdata.GetValue(
"Matcher_dzed");
1544 else if ( useMatcher ==
"pT_2" ) {
1545 double pTmatchLim_2 = 1.0;
1546 if (
inputdata.isTagDefined(
"Matcher_pTLim_2") ) pTmatchLim_2 =
inputdata.GetValue(
"Matcher_pTLim_2");
1549 else if ( useMatcher ==
"Truth" ) {
1562 else if ( useMatcher ==
"DeltaRZ" || useMatcher ==
"DeltaRZSinTheta" ) {
1566 if (
inputdata.isTagDefined(
"Matcher_deta" ) ) deta =
inputdata.GetValue(
"Matcher_deta");
1567 if (
inputdata.isTagDefined(
"Matcher_dphi" ) ) dphi =
inputdata.GetValue(
"Matcher_dphi");
1568 if (
inputdata.isTagDefined(
"Matcher_dzed" ) ) dzed =
inputdata.GetValue(
"Matcher_dzed");
1573 else if ( useMatcher ==
"pT_2" ) {
1574 double pTmatchLim_2 = 1.0;
1575 if (
inputdata.isTagDefined(
"Matcher_pTLim_2") ) pTmatchLim_2 =
inputdata.GetValue(
"Matcher_pTLim_2");
1578 else if ( useMatcher ==
"Truth" ) {
1602 if (
inputdata.isTagDefined(
"DataSets") ) {
1604 std::cout <<
"fetching dataset details" << std::endl;
1605 std::vector<std::string> datasets =
inputdata.GetStringVector(
"DataSets");
1606 for (
unsigned int ids=0 ;
ids<datasets.size() ;
ids++ ) {
1607 std::cout <<
"\tdataset " << datasets[
ids] << std::endl;
1609 std::vector<std::string> filenames_ =
d.datafiles();
1610 std::cout <<
"\tdataset contains " << filenames_.size() <<
" files" << std::endl;
1616 std::cerr <<
"no input data specified" << std::endl;
1626 bool show_release =
true;
1628 std::vector<std::string> release_data;
1630 std::string release_data_save =
"";
1632 if ( show_release ){
1638 TFile* finput = TFile::Open(
filenames[
i].c_str() );
1641 if ( finput==0 || !finput->IsOpen() || finput->IsZombie() ) {
1642 std::cerr <<
"Error: could not open input file: " <<
filenames[
i] << std::endl;
1646 TTree* dataTree = (TTree*)finput->Get(
"dataTree");
1647 TString* releaseData =
new TString(
"");
1650 dataTree->SetBranchAddress(
"ReleaseMetaData", &releaseData);
1652 for (
unsigned int i=0;
i<dataTree->GetEntries() ;
i++ ) {
1653 dataTree->GetEntry(
i);
1654 release_data.push_back( releaseData->Data() );
1655 if ( release_data_save != release_data.back() ) {
1656 std::cout <<
"main() release data: " << release_data.back() <<
" : " << *releaseData << std::endl;
1659 release_data_save = release_data.back();
1663 if ( finput )
delete finput;
1668 if ( !release_data.empty() ) {
1669 std::sort(release_data.begin(), release_data.end());
1670 release_data.erase(std::unique(release_data.begin(), release_data.end()), release_data.end());
1676 TTree* dataTree =
new TTree(
"dataTree",
"dataTree");
1677 TString* releaseData =
new TString(
"");
1680 dataTree->Branch(
"ReleaseMetaData",
"TString", &releaseData);
1682 for (
unsigned ird=0 ; ird<release_data.size() ; ird++ ) {
1683 *releaseData = release_data[ird];
1687 dataTree->Write(
"", TObject::kOverwrite);
1705 if ( Nentries==0 &&
inputdata.isTagDefined(
"Nentries") ) {
1711 typedef std::pair<int,double> zpair;
1712 std::vector<zpair> refz;
1713 std::vector<zpair> testz;
1715 std::vector<double> beamline_ref;
1716 std::vector<double> beamline_test;
1721 std::cout <<
"opening files" << std::endl;
1725 int pregrl_events = 0;
1726 int grl_counter = 0;
1729 std::cout <<
"starting event loop " <<
time_str() << std::endl;
1733 if ( nfiles!=0 && nfiles<max_files ) max_files = nfiles;
1742 if ( finput==0 || !finput->IsOpen() || finput->IsZombie() ) {
1743 std::cerr <<
"Error: could not open output file " <<
filenames[
ifile] << std::endl;
1747 TTree*
data = (TTree*)finput->Get(
"tree");
1750 std::cerr <<
"Error: cannot open TTree: " <<
filenames[
ifile] << std::endl;
1758 data->SetBranchAddress(
"TIDA::Event",&track_ev);
1764 unsigned cNentries =
data->GetEntries();
1770 for (
unsigned int i=0;
skip &&
run &&
i<cNentries ;
i++ ) {
1792 if ( !ingrl )
continue;
1798 if ( event_selector_flag && !es.in(
event ) )
continue;
1800 if ( mintime>
ts ) mintime =
ts;
1801 if ( maxtime<
ts ) maxtime =
ts;
1805 std::cout <<
"breaking out " <<
run << std::endl;
1815 hevent->Fill(
event );
1818 if ( (cNentries<10) ||
i%(cNentries/10)==0 ||
i%1000==0 ||
debugPrintout ) {
1819 std::cout <<
"run " << track_ev->
run_number()
1822 <<
"\tchains " << track_ev->
chains().size()
1824 std::cout <<
"\t : processed " <<
i <<
" events so far (" <<
int((1000*
i)/cNentries)*0.1 <<
"%)\t" <<
time_str() << std::endl;
1831 if ( nfiles>0 ) pfiles = nfiles;
1834 std::cout <<
"file entries=" <<
data->GetEntries();
1836 if (
data->GetEntries()<100 ) std::cout <<
" ";
1837 if (
data->GetEntries()<1000 ) std::cout <<
" ";
1838 if (
data->GetEntries()<10000 ) std::cout <<
" ";
1843 std::cout <<
"run " << track_ev->
run_number()
1846 <<
"\tchains " << track_ev->
chains().size()
1849 std::cout <<
"\t : processed " <<
ifile <<
" files so far (" <<
int((1
e3*
ifile)/pfiles)*0.1 <<
"%)\t" <<
time_str() << std::endl;
1862 truthTracks.clear();
1863 refPurityTracks.clear();
1867 mit->second.selector()->clear();
1873 const std::vector<TIDA::Chain>&
chains = track_ev->
chains();
1880 truthTracks.selectTracks(
chains[
ic][0].tracks() );
1887 for (
const std::string& rc : refChains){
1890 offTracks.selectTracks(
chains[
ic][0].tracks() );
1892 beamline_ref =
chains[
ic][0].user();
1901 std::vector<TIDA::Vertex> vertices;
1905 if ( vtxchain && vtxchain->
size()>0 ) {
1907 const std::vector<TIDA::Vertex>& mv = vtxchain->
at(0).
vertices();
1912 if (
debugPrintout ) std::cout <<
"vertices:\n" << mv << std::endl;
1914 if ( bestPTVtx || bestPT2Vtx ) {
1915 for (
size_t iv=0 ; iv<mv.size() ; iv++ ) {
1916 if ( mv[iv].Ntracks()==0 )
continue;
1917 double selection_ = 0.0;
1919 vtx_temp.selectTracks( offTracks.tracks() );
1920 for (
unsigned itr=0; itr<vtx_temp.tracks().
size(); itr++) {
1922 if ( bestPTVtx ) selection_ += std::fabs(tr->
pT());
1923 else if ( bestPT2Vtx ) selection_ += std::fabs(tr->
pT())*std::fabs(tr->
pT());
1930 if ( selectvtx!=-1 ) {
1931 vertices.push_back( mv[selectvtx] );
1934 else if ( vtxind>=0 ) {
1935 if (
size_t(vtxind)<mv.size() ) {
1936 vertices.push_back( mv[vtxind] );
1940 for (
size_t iv=0 ; iv<mv.size() ; iv++ ) {
1941 vertices.push_back( mv[iv] );
1946 filter_vertex.setVertex( vertices );
1952 for (
unsigned iv=0 ; iv<mv.size() ; iv++ ) {
1953 int Ntracks = mv[iv].Ntracks();
1954 if ( Ntracks>NVtxTrackCut ) {
1971 bool foundReference =
false;
1980 for (
const std::string& rc : refChains ) {
1989 std::string refname = mitr->first;
1998 *mitr->second.tom() = rtom;
2012 if ( !foundReference )
continue;
2015 std::cout <<
"reference chain:\n" << *refchain << std::endl;
2029 if ( analitr==analysis.end() )
continue;
2032 std::cout <<
"test chain:\n" <<
chain << std::endl;
2037 std::vector<TIDA::Roi*>
rois;
2057 if ( TnP_tool->
type0() == TnP_tool->
type1() ) {
2059 rois = TnP_tool->
GetRois( track_ev->
chains(), selector0, filter0, invmass, invmass_obj, rtom );
2073 for (
size_t ia=0 ; ia<selector1->
tracks().size() ; ia++ ) {
2075 std::cout << ia <<
"\t" << *
track << std::endl;
2076 std::cout <<
"\t" << rtom1->
object(
track->id()) << std::endl;
2079 std::cout << *TnP_tool << std::endl;
2082 rois = TnP_tool->
GetRois( track_ev->
chains(), selector0, filter0, selector1, filter1, invmass, invmass_obj, rtom, rtom1 );
2092 if (
false ) std::cout <<
"++++++++++++ rois size = " <<
rois.size() <<
" +++++++++++" << std::endl;
2095 for (
unsigned int ir=0 ;
ir<
rois.size() ;
ir++ ) {
2105 testTracks.selectTracks( troi.
tracks() );
2108 std::vector<TIDA::Track*> testp = testTracks.tracks();
2113 if ( filterRoi && !roiFilter.filter( roi ) )
continue;
2116 const std::vector<TIDA::Vertex>& mvt = troi.
vertices();
2118 std::vector<TIDA::Vertex> vertices_test;
2123 if ( bestPTVtx_rec || bestPT2Vtx_rec ) {
2127 for (
unsigned iv=0 ; iv<mvt.size() ; iv++ ) {
2128 double selection_ = 0.0;
2130 vtx_temp.selectTracks( testp );
2131 for (
unsigned itr=0; itr<vtx_temp.tracks().
size(); itr++) {
2133 if ( bestPTVtx ) selection_ += std::fabs(tr->
pT());
2134 else if ( bestPT2Vtx ) selection_ += std::fabs(tr->
pT())*std::fabs(tr->
pT());
2141 if ( selectvtx!=-1 ) {
2143 if ( useVertexTracks ) selected.selectTracks( testp );
2144 vertices_test.push_back(selected);
2147 else if ( vtxind_rec!=-1 ) {
2148 if (
unsigned(vtxind_rec)<mvt.size() ) {
2150 if ( useVertexTracks ) selected.selectTracks( testp );
2151 vertices_test.push_back( selected );
2155 for (
unsigned iv=0 ; iv<mvt.size() ; iv++ ) {
2157 if ( useVertexTracks ) selected.selectTracks( testp );
2158 vertices_test.push_back( selected );
2163 beamline_test =
rois[
ir]->user();
2170 if ( correctBeamlineTest ) {
2171 if ( beamTest.size()==2 ) analy_track->
setBeamTest( beamTest[0], beamTest[1] );
2174 if ( !
inputdata.isTagDefined(
"BeamTest") ) {
2175 if ( beamline_test.size()==2 ) analy_track->
setBeamTest( beamline_test[0], beamline_test[1] );
2181 if ( correctBeamlineRef ) {
2182 if ( beamRef.size()==2 ) analy_track->
setBeamRef( beamRef[0], beamRef[1] );
2185 if ( !
inputdata.isTagDefined(
"BeamRef") ) {
2186 if ( beamline_ref.size()==2 ) analy_track->
setBeamRef( beamline_ref[0], beamline_ref[1] );
2205 bool customRefRoi_thisChain =
false;
2207 if ( use_custom_ref_roi ) {
2208 if ( customRoi_chains.size() ) {
2209 if ( customRoi_chains.find(
chain.name() )!=customRoi_chains.end() ) customRefRoi_thisChain =
true;
2211 else customRefRoi_thisChain =
true;
2214 if ( use_custom_ref_roi && customRefRoi_thisChain ) {
2215 refRoi =
makeCustomRefRoi( roi, custRefRoi_params[0], custRefRoi_params[1], custRefRoi_params[2] );
2231 std::vector<TIDA::Track*> refp_vec = refTracks.tracks( refFilter );
2236 if ( select_roi )
dynamic_cast<Filter_Combined*
>(truthFilter)->setRoi(&roi);
2237 const std::vector<TIDA::Track*>& truth = truthTracks.tracks(truthFilter);
2238 const std::vector<TIDA::Track*>& refp_tmp = refp_vec;
2241 truth_matcher->
match( refp_tmp, truth );
2243 std::vector<TIDA::Track*> refp_matched;
2246 for (
unsigned i=0 ;
i<refp_vec.size() ;
i++ ) {
2247 if ( truth_matcher->
matched( refp_vec[
i] ) ) refp_matched.push_back( refp_vec[
i] );
2251 refp_vec = refp_matched;
2255 if ( use_pt_ordered_ref ) {
2258 size_t nRefTracks = refp_vec.size();
2260 std::vector<TIDA::Track*> refp_chosenPtOrd;
2264 for (
size_t sidx=refPtOrd_indices.size() ; sidx-- ; ) {
2265 if ( refPtOrd_indices.at(sidx)>nRefTracks ) {
2266 std::cout <<
"WARNING: for event " <<
event
2267 <<
", pT ordered reference track at vector position " << refPtOrd_indices.at(sidx)
2268 <<
" requested but not found" << std::endl;
2274 for (
size_t sidx=refPtOrd_indices.size() ; sidx-- ; ) {
2275 for (
size_t idx=0 ;
idx<refp_vec.size() ;
idx++ ) {
2276 if (
idx == refPtOrd_indices.at(sidx) ) {
2277 refp_chosenPtOrd.push_back(refp_vec.at(
idx));
2284 refp_vec = refp_chosenPtOrd;
2295 if ( ptconfig!=
"" ) {
2298 std::vector<TIDA::Track*> reft; reft.reserve(refp_vec.size());
2299 for ( std::vector<TIDA::Track*>::const_iterator itr=refp_vec.begin() ; itr!=refp_vec.end() ; ++itr ) {
2300 if ( std::fabs((*itr)->pT())>=
pt ) reft.push_back( *itr );
2309 if ( d0config!=
"" ) {
2312 std::vector<TIDA::Track*> reft; reft.reserve(refp_vec.size());
2313 for ( std::vector<TIDA::Track*>::const_iterator itr=refp_vec.begin() ; itr!=refp_vec.end() ; ++itr ) {
2314 if ( std::fabs((*itr)->a0())<=
d0 ) reft.push_back( *itr );
2331 if ( objectET !=
"" ) {
2333 double ETconfig =
std::atof( objectET.c_str() );
2339 while ( itr!=refp_vec.end() ) {
2342 if ( tobj==0 || tobj->
pt()<ETconfig )
2343 itr=refp_vec.erase( itr );
2351 const std::vector<TIDA::Track*>& refp = refp_vec;
2368 std::vector<TIDA::Vertex> vertices_roi;
2376 vertices_roi.clear();
2378 const std::vector<TIDA::Vertex>& mv = vertices;
2383 for (
unsigned iv=0 ; iv<mv.size() ; iv++ ) {
2389 bool accept_vertex =
false;
2390 if ( roi.composite() ) {
2391 for (
size_t ir=0 ;
ir<roi.size() ;
ir++ ) {
2392 if ( roi[
ir]->zedMinus()<=vx.
z() && roi[
ir]->zedPlus()>=vx.
z() ) accept_vertex =
true;
2396 if ( roi.zedMinus()<=vx.
z() && roi.zedPlus()>=vx.
z() ) accept_vertex =
true;
2399 if ( !accept_vertex )
continue;
2405 if ( useVertexTracks ) {
2408 vertex_roi.selectTracks( refp );
2409 trackcount = vertex_roi.Ntracks();
2410 if ( trackcount>=ntracks && trackcount>0 ) {
2411 vertices_roi.push_back( vertex_roi );
2416 for (
unsigned itr=0; itr<refp.size(); itr++){
2419 double dzsintheta = std::fabs( (vx.
z()-tr->
z0()) *
std::sin(theta_) );
2420 if( dzsintheta < 1.5 ) trackcount++;
2425 if ( trackcount>=ntracks && trackcount>0 ) {
2427 vertices_roi.push_back( vertex_roi );
2435 if ( rotate_testtracks )
for (
size_t i=testp.size() ;
i-- ; ) testp[
i]->rotate();
2441 if ( monitorZBeam ) {
2442 if ( beamline_ref.size()>2 && beamline_test.size()>2 ) {
2443 refz.push_back( zpair(
lb, beamline_ref[2]) );
2444 testz.push_back( zpair(
lb, beamline_test[2]) );
2448 matcher->
match( refp, testp);
2450 if ( tom.
status() ) analitr->second->execute( refp, testp, matcher, &tom );
2451 else analitr->second->execute( refp, testp, matcher );
2454 analitr->second->store().find( vtxanal,
"rvtx" );
2458 if ( vtxanal && ( refp.size() >=
size_t(ntracks) ) ) {
2477 std::cout <<
"\nselected tracks:" <<
chain.name() << std::endl;
2478 std::cout <<
"ref tracks refp.size() " << refp.size() <<
"\n" << refp << std::endl;
2479 std::cout <<
"test tracks testp.size() " << testp.size() <<
"\n" << testp << std::endl;
2481 TrackAssociator::map_type::const_iterator titr = matcher->TrackAssociator::matched().begin();
2482 TrackAssociator::map_type::const_iterator
tend = matcher->TrackAssociator::matched().end();
2484 std::cout <<
"track matches:\n";
2485 while (titr!=
tend) {
2486 std::cout <<
"\t" <<
im++ <<
"\t" << *titr->first <<
" ->\n\t\t" << *titr->second << std::endl;
2491 std::cout <<
"completed : " <<
chain.name() <<
"\n";
2492 std::cout <<
"-----------------------------------" << std::endl;
2497 if ( _matcher->size()<refp.size() ) {
2499 if ( refp.size()==1 && testp.size()==0 ) {
2500 std::cout << track_ev->
chains()[
ic] <<std::endl;
2501 std::cout <<
"roi\n" << track_ev->
chains()[
ic].rois()[
ir] << endl;
2511 const std::vector<TIDA::Track*>& refpp = refPurityTracks.tracks( refFilter );
2513 testPurityTracks.clear();
2515 testPurityTracks.selectTracks( troi.
tracks() );
2516 std::vector<TIDA::Track*> testpp = testPurityTracks.tracks();
2518 matcher->
match(refpp, testpp);
2523 if ( analitrp == analysis.end() )
continue;
2526 analitrp->second->execute( refpp, testpp, matcher );
2547 <<
"\ttimes " << mintime <<
" " << maxtime
2548 <<
"\t( grl: " << grl_counter <<
" / " << pregrl_events <<
" )" << std::endl;
2550 if ( monitorZBeam )
zbeam zb( refz, testz );
2564 analyses[
i]->store().find( vtxanal,
"rvtx" );
2565 if ( vtxanal ) vtxanal->
finalise();
2576 mit->second.Clean();
2580 if ( ex_matcher )
delete ex_matcher;
2582 std::cout <<
"done" << std::endl;