267 std:: cout <<
"----------------- Event " <<
evth_->
getEventNumber() <<
" -----------------" << std::endl;
273 if (
debug_) std::cout <<
"Check pbc_configs" << std::endl;
276 int check_run = std::stoi(run.key());
277 if(check_run > run_number)
280 closest_run = check_run;
284 bpc_configs_[std::to_string(closest_run)][
"beamspot_y"],
285 bpc_configs_[std::to_string(closest_run)][
"beamspot_z"]};
291 double apMass = -0.9;
293 double apEnergy = -0.9;
295 double vdMass = -0.9;
297 double vdEnergy = -0.9;
299 if (
debug_) std::cout <<
"plot trigger info" << std::endl;
307 if (
debug_) std::cout <<
"plot vtx N" << std::endl;
311 for(
int i = 0; i <
mcParts_->size(); i++)
313 if(
mcParts_->at(i)->getPDG() == 622)
315 apMass =
mcParts_->at(i)->getMass();
316 apZ =
mcParts_->at(i)->getVertexPosition().at(2);
317 apEnergy =
mcParts_->at(i)->getEnergy();
319 if(
mcParts_->at(i)->getPDG() == 625)
321 vdMass =
mcParts_->at(i)->getMass();
322 vdZ =
mcParts_->at(i)->getVertexPosition().at(2);
323 vdEnergy =
mcParts_->at(i)->getEnergy();
330 std::vector<Vertex*> selected_vtxs;
331 bool passVtxPresel =
false;
334 std::cout<<
"Number of vertices found in event: "<<
vtxs_->size()<<std::endl;
338 for (
int i_vtx = 0; i_vtx <
vtxs_->size(); i_vtx++ ) {
352 bool foundParts =
_ah->GetParticlesFromVtx(vtx,ele,pos);
354 if(
debug_) std::cout<<
"NewVertexAnaProcessor::WARNING::Found vtx without ele/pos. Skip."<<std::endl;
358 if (
debug_) std::cout <<
"got parts" << std::endl;
363 std::cout<<
"Check Ele/Pos Track momenta"<<std::endl;
364 std::cout<<ele_trk.
getP()<<
" "<<pos_trk.
getP()<<std::endl;
393 std::cout<<
"Corrected Ele/Pos Track momenta"<<std::endl;
394 std::cout<<ele_trk.
getP()<<
" "<<pos_trk.
getP()<<std::endl;
397 double invm_smear = 1.;
399 double unsmeared_prod = ele_trk.
getP()*pos_trk.
getP();
402 double smeared_prod = ele_trk.
getP()*pos_trk.
getP();
403 invm_smear = sqrt(smeared_prod/unsmeared_prod);
413 if (
debug_) std::cout <<
"got tracks" << std::endl;
420 TLorentzVector p_ele;
422 TLorentzVector p_pos;
429 if (
debug_) std::cout <<
"start selection" << std::endl;
458 double botClusTime = 0.0;
463 if (!
vtxSelector->passCutLt(
"botCluTime_lt", botClusTime, weight))
466 if (!
vtxSelector->passCutGt(
"botCluTime_gt", botClusTime, weight))
470 if (!
vtxSelector->passCutLt(
"eleposCluTimeDiff_lt",fabs(corr_eleClusterTime - corr_posClusterTime),weight))
517 if (!
vtxSelector->passCutLt(
"eleMom_lt",ele_mom.Mag(),weight))
521 if (!
vtxSelector->passCutGt(
"eleMom_gt",ele_mom.Mag(),weight))
525 if (!
vtxSelector->passCutGt(
"posMom_gt",pos_mom.Mag(),weight))
533 if (!
vtxSelector->passCutGt(
"eleN2Dhits_gt",ele2dHits,weight)) {
542 if (!
vtxSelector->passCutGt(
"posN2Dhits_gt",pos2dHits,weight)) {
557 if (!
vtxSelector->passCutLt(
"chi2unc_lt",vtx->getChi2(),weight))
561 if (!
vtxSelector->passCutLt(
"maxVtxMom_lt",(ele_mom+pos_mom).Mag(),weight))
566 if (!
vtxSelector->passCutGt(
"minVtxMom_gt",(ele_mom+pos_mom).Mag(),weight))
569 if (
debug_) std::cout <<
"fill 1D Vertex" << std::endl;
577 if (
debug_) std::cout <<
"fill track histos" << std::endl;
578 double ele_pos_dt = corr_eleClusterTime - corr_posClusterTime;
579 double psum = ele_mom.Mag()+pos_mom.Mag();
583 _vtx_histos->Fill1DHisto(
"ele_track_n2dhits_h", ele2dHits, weight);
584 _vtx_histos->Fill1DHisto(
"pos_track_n2dhits_h", pos2dHits, weight);
585 _vtx_histos->Fill1DHisto(
"vtx_Psum_h", p_ele.P()+p_pos.P(), weight);
586 _vtx_histos->Fill1DHisto(
"vtx_Esum_h", ele_E + pos_E, weight);
587 _vtx_histos->Fill1DHisto(
"vtx_smear_InvM_h", invm_smear*(vtx->getInvMass()), weight);
588 _vtx_histos->Fill1DHisto(
"ele_pos_clusTimeDiff_h", (corr_eleClusterTime - corr_posClusterTime), weight);
603 _vtx_histos->Fill2DHisto(
"ele_pos_clusTimeDiff_v_pSum_hh",ele_mom.Mag()+pos_mom.Mag(), ele_pos_dt, weight);
629 passVtxPresel =
true;
631 selected_vtxs.push_back(vtx);
637 _vtx_histos->Fill1DHisto(
"n_vertices_h",selected_vtxs.size());
649 if (
debug_) std::cout <<
"start regions" << std::endl;
654 Vertex* goodVtx =
nullptr;
655 std::vector<Vertex*> goodVtxs;
660 for (
auto vtx : selected_vtxs) {
669 _ah->GetParticlesFromVtx(vtx,ele,pos);
712 double invm_smear = 1.;
714 double unsmeared_prod = ele_trk.
getP()*pos_trk.
getP();
717 double smeared_prod = ele_trk.
getP()*pos_trk.
getP();
718 invm_smear = sqrt(smeared_prod/unsmeared_prod);
725 TLorentzVector p_ele;
727 TLorentzVector p_pos;
731 std::vector<int> ele_hit_layers = ele_trk.
getHitLayers();
734 int ele_lastlayer = 0;
735 for(
int i=0; i<ele_hit_layers.size(); i++)
737 int layer = ele_hit_layers.at(i);
738 ele_lastlayer = layer;
739 if (layer == 0) ele_Si0++;
740 if (layer == 1) ele_Si1++;
743 std::vector<int> pos_hit_layers = pos_trk.
getHitLayers();
746 int pos_lastlayer = 0;
747 for(
int i=0; i<pos_hit_layers.size(); i++)
749 int layer = pos_hit_layers.at(i);
750 pos_lastlayer = layer;
751 if (layer == 0) pos_Si0++;
752 if (layer == 1) pos_Si1++;
760 std::cout<<
"Check on ele_Track"<<std::endl;
764 bool foundL1ele =
false;
765 bool foundL2ele =
false;
766 _ah->InnermostLayerCheck(&ele_trk, foundL1ele, foundL2ele);
770 std::cout<<
"Check on pos_Track"<<std::endl;
773 bool foundL1pos =
false;
774 bool foundL2pos =
false;
776 _ah->InnermostLayerCheck(&pos_trk, foundL1pos, foundL2pos);
779 std::cout<<
"Check on pos_Track"<<std::endl;
780 std::cout<<
"Innermost:"<<foundL1pos<<
" Second Innermost:"<<foundL2pos<<std::endl;
815 double botClusTime = 0.0;
827 if (!
_reg_vtx_selectors[region]->passCutLt(
"eleposCluTimeDiff_lt",fabs(corr_eleClusterTime - corr_posClusterTime),weight))
909 if (!
_reg_vtx_selectors[region]->passCutLt(
"maxVtxMom_lt",(ele_mom+pos_mom).Mag(),weight))
913 if (!
_reg_vtx_selectors[region]->passCutGt(
"minVtxMom_gt",(ele_mom+pos_mom).Mag(),weight))
919 if (!
_reg_vtx_selectors[region]->passCutEq(
"L1Requirement_eq",(
int)(foundL1ele&&foundL1pos),weight))
923 if (!
_reg_vtx_selectors[region]->passCutEq(
"L2Requirement_eq",(
int)(foundL2ele&&foundL2pos),weight))
990 std::map<int, int> nHits4part;
991 for(
int i =0; i < ele_trk.
getMcpHits().size(); i++)
993 int partID = ele_trk.
getMcpHits().at(i).second;
994 if ( nHits4part.find(partID) == nHits4part.end() )
997 nHits4part[partID] = 1;
1002 nHits4part[partID]++;
1009 for (std::map<int,int>::iterator it=nHits4part.begin(); it!=nHits4part.end(); ++it)
1011 if(it->second > maxNHits)
1013 maxNHits = it->second;
1019 int isRadEle = -999;
1020 int isRecEle = -999;
1023 trueEleP.SetXYZ(-999,-999,-999);
1024 truePosP.SetXYZ(-999,-999,-999);
1026 float trueEleE = -1;
1027 float truePosE = -1;
1028 for(
int i = 0; i <
mcParts_->size(); i++)
1030 int momPDG =
mcParts_->at(i)->getMomPDG();
1033 std::vector<double> lP =
mcParts_->at(i)->getMomentum();
1034 trueEleP.SetXYZ(lP[0],lP[1],lP[2]);
1035 trueEleE =
mcParts_->at(i)->getEnergy();
1040 std::vector<double> lP =
mcParts_->at(i)->getMomentum();
1041 truePosP.SetXYZ(lP[0],lP[1],lP[2]);
1042 truePosE =
mcParts_->at(i)->getEnergy();
1045 if(trueEleP.X() != -999 && truePosP.X() != -999){
1046 truePsum = trueEleP.Mag() + trueEleP.Mag();
1047 trueEsum = trueEleE + truePosE;
1050 if(
mcParts_->at(i)->getID() != maxID)
continue;
1053 if(momPDG == 623) isRecEle = 1;
1056 double momRatio = recEleP.Mag() / trueEleP.Mag();
1057 double momAngle = trueEleP.Angle(recEleP) * TMath::RadToDeg();
1058 if (!
_reg_vtx_selectors[region]->passCutLt(
"momRatio_lt", momRatio, weight))
continue;
1059 if (!
_reg_vtx_selectors[region]->passCutGt(
"momRatio_gt", momRatio, weight))
continue;
1060 if (!
_reg_vtx_selectors[region]->passCutLt(
"momAngle_lt", momAngle, weight))
continue;
1062 if (!
_reg_vtx_selectors[region]->passCutEq(
"isRadEle_eq", isRadEle, weight))
continue;
1063 if (!
_reg_vtx_selectors[region]->passCutEq(
"isNotRadEle_eq", isRadEle, weight))
continue;
1064 if (!
_reg_vtx_selectors[region]->passCutEq(
"isRecEle_eq", isRecEle, weight))
continue;
1069 goodVtxs.push_back(vtx);
1076 _reg_vtx_histos[region]->Fill1DHisto(
"n_vertices_h", nGoodVtx, weight);
1079 for(std::vector<Vertex*>::iterator it = goodVtxs.begin(); it != goodVtxs.end(); it++){
1086 if (!vtx || !
_ah->GetParticlesFromVtx(vtx,ele,pos))
1127 double invm_smear = 1.;
1129 double unsmeared_prod = ele_trk.
getP()*pos_trk.
getP();
1132 double smeared_prod = ele_trk.
getP()*pos_trk.
getP();
1133 invm_smear = sqrt(smeared_prod/unsmeared_prod);
1137 std::vector<int> ele_hit_layers = ele_trk.
getHitLayers();
1140 int ele_lastlayer = 0;
1141 for(
int i=0; i<ele_hit_layers.size(); i++)
1143 int layer = ele_hit_layers.at(i);
1144 ele_lastlayer = layer;
1145 if (layer == 0) ele_Si0++;
1146 if (layer == 1) ele_Si1++;
1149 std::vector<int> pos_hit_layers = pos_trk.
getHitLayers();
1152 int pos_lastlayer = 0;
1153 for(
int i=0; i<pos_hit_layers.size(); i++)
1155 int layer = pos_hit_layers.at(i);
1156 pos_lastlayer = layer;
1157 if (layer == 0) pos_Si0++;
1158 if (layer == 1) pos_Si1++;
1162 std::vector<float> vtx_cov = vtx->getCovariance();
1163 float cxx = vtx_cov.at(0);
1164 float cyx = vtx_cov.at(1);
1165 float cyy = vtx_cov.at(2);
1166 float czx = vtx_cov.at(3);
1167 float czy = vtx_cov.at(4);
1168 float czz = vtx_cov.at(5);
1172 int L1L2hitCode = 0;
1179 if (!
_reg_vtx_selectors[region]->passCutLt(
"hitCode_lt",((
double)L1L2hitCode)-0.5, weight))
continue;
1180 if (!
_reg_vtx_selectors[region]->passCutGt(
"hitCode_gt",((
double)L1L2hitCode)+0.5, weight))
continue;
1182 _reg_vtx_histos[region]->Fill1DHisto(
"hitCode_h", L1L2hitCode,weight);
1183 _reg_vtx_histos[region]->Fill1DHisto(
"L1hitCode_h", L1hitCode,weight);
1184 _reg_vtx_histos[region]->Fill1DHisto(
"L2hitCode_h", L2hitCode,weight);
1189 bool hasL1ele =
false;
1190 bool hasL2ele =
false;
1191 _ah->InnermostLayerCheck(&ele_trk, hasL1ele, hasL2ele);
1193 bool hasL1pos =
false;
1194 bool hasL2pos =
false;
1195 _ah->InnermostLayerCheck(&pos_trk, hasL1pos, hasL2pos);
1214 double ele_pos_dt = corr_eleClusterTime - corr_posClusterTime;
1215 double psum = ele_mom.Mag()+pos_mom.Mag();
1237 TLorentzVector p_ele;
1239 TLorentzVector p_pos;
1251 _reg_vtx_histos[region]->Fill1DHisto(
"ele_pos_clusTimeDiff_h", (corr_eleClusterTime - corr_posClusterTime), weight);
1252 _reg_vtx_histos[region]->Fill1DHisto(
"ele_track_n2dhits_h", ele2dHits, weight);
1253 _reg_vtx_histos[region]->Fill1DHisto(
"pos_track_n2dhits_h", pos2dHits, weight);
1254 _reg_vtx_histos[region]->Fill1DHisto(
"vtx_Psum_h", p_ele.P()+p_pos.P(), weight);
1256 _reg_vtx_histos[region]->Fill1DHisto(
"vtx_smear_InvM_h", invm_smear*(vtx->getInvMass()), weight);
1271 _reg_vtx_histos[region]->Fill2DHisto(
"vtx_Psum_vs_true_Psum_hh",p_ele.P()+p_pos.P(), truePsum, weight);
1272 _reg_vtx_histos[region]->Fill1DHisto(
"true_vtx_psum_h",truePsum,weight);
1275 double reconz = vtx->getZ();
1276 double ele_trk_z0 = ele_trk.
getZ0();
1277 double ele_trk_z0err = ele_trk.
getZ0Err();
1278 double pos_trk_z0 = pos_trk.
getZ0();
1279 double pos_trk_z0err = pos_trk.
getZ0Err();
1285 double vtx_proj_x = -999.9;
1286 double vtx_proj_y = -999.9;
1287 double vtx_proj_x_sig = -999.9;
1288 double vtx_proj_y_sig = -999.9;
1289 double vtx_proj_sig = -999.9;
1292 vtx_proj_x, vtx_proj_y, vtx_proj_x_sig, vtx_proj_y_sig, vtx->getX(), vtx->getY(),
1293 reconz, vtx->getP().X(), vtx->getP().Y(), vtx->getP().Z());
1295 _reg_vtx_histos[region]->Fill2DHisto(
"unc_vtx_x_v_unc_vtx_y_hh", vtx->getX(), vtx->getY());
1296 _reg_vtx_histos[region]->Fill2DHisto(
"unc_vtx_proj_x_v_unc_vtx_proj_y_hh", vtx_proj_x, vtx_proj_y);
1297 _reg_vtx_histos[region]->Fill2DHisto(
"unc_vtx_proj_x_y_significance_hh", vtx_proj_x_sig, vtx_proj_y_sig);
1298 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_vtx_proj_significance_hh", vtx_proj_sig, reconz);
1299 _reg_vtx_histos[region]->Fill2DHisto(
"vtx_track_reconz_v_Z0err_hh", ele_trk_z0err, reconz);
1300 _reg_vtx_histos[region]->Fill2DHisto(
"vtx_track_reconz_v_Z0err_hh", pos_trk_z0err, reconz);
1301 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_z0_hh", ele_trk_z0, reconz);
1302 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_z0_hh", pos_trk_z0, reconz);
1306 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_cxx_hh", cxx, reconz);
1307 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_cyy_hh", cyy, reconz);
1308 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_czz_hh", czz, reconz);
1309 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_cyx_hh", cyx, reconz);
1310 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_czx_hh", czx, reconz);
1311 _reg_vtx_histos[region]->Fill2DHisto(
"recon_z_v_czy_hh", czy, reconz);
1324 _reg_vtx_histos[region]->Fill2DHisto(
"ele_pos_clusTimeDiff_v_pSum_hh",ele_mom.Mag()+pos_mom.Mag(), ele_pos_dt, weight);
1355 vtxPosSvt.SetX(vtx->getX());
1356 vtxPosSvt.SetY(vtx->getY());
1357 vtxPosSvt.SetZ(vtx->getZ());
1358 vtxPosSvt.RotateY(-0.0305);
1363 _reg_tuples[region]->setVariableValue(
"ap_true_vtx_z", apZ);
1364 _reg_tuples[region]->setVariableValue(
"ap_true_vtx_mass", apMass);
1365 _reg_tuples[region]->setVariableValue(
"ap_true_vtx_energy", apEnergy);
1366 _reg_tuples[region]->setVariableValue(
"vd_true_vtx_z", vdZ);
1367 _reg_tuples[region]->setVariableValue(
"vd_true_vtx_mass", vdMass);
1368 _reg_tuples[region]->setVariableValue(
"vd_true_vtx_energy", vdEnergy);
1369 _reg_tuples[region]->setVariableValue(
"hitCode",
float(L1L2hitCode));
1370 _reg_tuples[region]->setVariableValue(
"L1hitCode",
float(L1hitCode));
1371 _reg_tuples[region]->setVariableValue(
"L2hitCode",
float(L2hitCode));
1374 _reg_tuples[region]->setVariableValue(
"unc_vtx_mass", vtx->getInvMass());
1375 _reg_tuples[region]->setVariableValue(
"unc_vtx_z" , vtxPosSvt.Z());
1376 _reg_tuples[region]->setVariableValue(
"unc_vtx_chi2", vtx->getChi2());
1377 _reg_tuples[region]->setVariableValue(
"unc_vtx_psum", p_ele.P()+p_pos.P());
1378 _reg_tuples[region]->setVariableValue(
"unc_vtx_px", vtx->getP().X());
1379 _reg_tuples[region]->setVariableValue(
"unc_vtx_py", vtx->getP().Y());
1380 _reg_tuples[region]->setVariableValue(
"unc_vtx_pz", vtx->getP().Z());
1381 _reg_tuples[region]->setVariableValue(
"unc_vtx_x", vtx->getX());
1382 _reg_tuples[region]->setVariableValue(
"unc_vtx_y", vtx->getY());
1383 _reg_tuples[region]->setVariableValue(
"unc_vtx_proj_x", vtx_proj_x);
1384 _reg_tuples[region]->setVariableValue(
"unc_vtx_proj_y", vtx_proj_y);
1385 _reg_tuples[region]->setVariableValue(
"unc_vtx_proj_x_sig", vtx_proj_x_sig);
1386 _reg_tuples[region]->setVariableValue(
"unc_vtx_proj_y_sig", vtx_proj_y_sig);
1387 _reg_tuples[region]->setVariableValue(
"unc_vtx_proj_sig", vtx_proj_sig);
1388 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_pos_clust_dt", corr_eleClusterTime - corr_posClusterTime);
1390 _reg_tuples[region]->setVariableValue(
"unc_vtx_cxx", cxx);
1391 _reg_tuples[region]->setVariableValue(
"unc_vtx_cyy", cyy);
1392 _reg_tuples[region]->setVariableValue(
"unc_vtx_czz", czz);
1393 _reg_tuples[region]->setVariableValue(
"unc_vtx_cyx", cyx);
1394 _reg_tuples[region]->setVariableValue(
"unc_vtx_czy", czy);
1395 _reg_tuples[region]->setVariableValue(
"unc_vtx_czx", czx);
1396 _reg_tuples[region]->setVariableValue(
"unc_vtx_deltaZ", deltaZ);
1399 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_p", ele_trk.
getP());
1401 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_d0", ele_trk.
getD0());
1402 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_phi0", ele_trk.
getPhi());
1405 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_z0", ele_trk.
getZ0());
1407 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_clust_dt", ele_trk.
getTrackTime() - corr_eleClusterTime);
1413 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_nhits",ele2dHits);
1414 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_lastlayer",ele_lastlayer);
1415 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_si0",ele_Si0);
1416 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_track_si1",ele_Si1);
1418 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_p", pos_trk.
getP());
1420 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_d0", pos_trk.
getD0());
1421 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_phi0", pos_trk.
getPhi());
1424 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_z0", pos_trk.
getZ0());
1426 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_clust_dt", pos_trk.
getTrackTime() - corr_posClusterTime);
1432 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_nhits",pos2dHits);
1433 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_lastlayer",pos_lastlayer);
1434 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_si0",pos_Si0);
1435 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_track_si1",pos_Si1);
1440 _reg_tuples[region]->setVariableValue(
"unc_vtx_ele_clust_corr_t",corr_eleClusterTime);
1444 _reg_tuples[region]->setVariableValue(
"unc_vtx_pos_clust_corr_t",corr_posClusterTime);