220 B2INFO(
"ECLBhabhaTCollector: Experiment = " <<
m_EventMetaData->getExperiment() <<
230 auto TimevsCrysPrevCrateCalibPrevCrystCalib =
new TH2F(
"TimevsCrysPrevCrateCalibPrevCrystCalib",
231 "Time t psi - ts - tcrate (previous calibs) vs crystal cell ID;crystal cell ID;Time t_psi with previous calib (ns)",
233 registerObject<TH2F>(
"TimevsCrysPrevCrateCalibPrevCrystCalib", TimevsCrysPrevCrateCalibPrevCrystCalib);
235 auto TimevsCratePrevCrateCalibPrevCrystCalib =
new TH2F(
"TimevsCratePrevCrateCalibPrevCrystCalib",
236 "Time t psi - ts - tcrate (previous calibs) vs crate ID;crate ID;Time t_psi previous calib (ns)",
237 52, 1, 52 + 1, nbins, min_t, max_t);
238 registerObject<TH2F>(
"TimevsCratePrevCrateCalibPrevCrystCalib", TimevsCratePrevCrateCalibPrevCrystCalib);
240 auto TimevsCrysNoCalibrations =
new TH2F(
"TimevsCrysNoCalibrations",
245 auto TimevsCrateNoCalibrations =
new TH2F(
"TimevsCrateNoCalibrations",
246 "Time tpsi vs crate ID;crate ID;Time t_psi (ns)", 52, 1, 52 + 1, nbins, min_t, max_t);
249 auto TimevsCrysPrevCrateCalibNoCrystCalib =
new TH2F(
"TimevsCrysPrevCrateCalibNoCrystCalib",
250 "Time tpsi - tcrate (previous calib) vs crystal cell ID;crystal cell ID;Time t_psi including previous crate calib (ns)",
252 registerObject<TH2F>(
"TimevsCrysPrevCrateCalibNoCrystCalib", TimevsCrysPrevCrateCalibNoCrystCalib);
254 auto TimevsCrateNoCrateCalibPrevCrystCalib =
new TH2F(
"TimevsCrateNoCrateCalibPrevCrystCalib",
255 "Time tpsi - ts (previous calib) vs crate ID;crate ID;Time t_psi including previous crystal calib (ns)",
256 52, 1, 52 + 1, nbins, min_t, max_t);
257 registerObject<TH2F>(
"TimevsCrateNoCrateCalibPrevCrystCalib", TimevsCrateNoCrateCalibPrevCrystCalib);
277 auto tcrateDatabase_ns =
new TH1F(
"tcrateDatabase_ns",
";crate id;tcrate derived from database", 52, 1, 52 + 1);
281 auto databaseCounter =
new TH1I(
"databaseCounter",
282 ";A database was read in;Number of times database was saved to histogram", 1, 1, 2);
286 auto numCrystalEntriesPerEvent =
new TH1F(
"numCrystalEntriesPerEvent",
287 ";Number crystal entries;Number of events", 15, 0, 15);
290 auto cutflow =
new TH1F(
"cutflow",
";Cut label number;Number of events passing cut", 20, 0, 20);
293 auto svdEventT0 =
new TH1D(
"svdEventT0",
";EventT0 from SVD", 200, -50, 50);
296 auto maxEcrsytalEnergyFraction =
new TH1F(
"maxEcrsytalEnergyFraction",
297 ";Maximum energy crystal energy / (sum) cluster energy;Number", 22, 0, 1.1);
300 auto refCrysIDzeroingCrate =
new TH1F(
"refCrysIDzeroingCrate",
";cell id;Boolean - is reference crystal",
304 auto SVDEventT0Correction =
new TH1F(
"SVDEventT0Correction",
";;SVD event t0 offset correction [ns]", 1, 1, 2);
316 " ns correction to SVD event t0 will be applied");
324 int cutIndexPassed = 0;
326 B2DEBUG(22,
"Cutflow: no cuts: index = " << cutIndexPassed);
336 B2WARNING(
"SoftwareTriggerResult required to select bhabha event is not found");
342 const std::map<std::string, int>& fresults =
m_TrgResult->getResults();
343 if (fresults.find(
"software_trigger_cut&skim&accept_bhabha") == fresults.end()) {
344 B2WARNING(
"Can't find required bhabha trigger identifier");
348 const bool eBhabha = (
m_TrgResult->getResult(
"software_trigger_cut&skim&accept_bhabha") ==
350 B2DEBUG(22,
"eBhabha (trigger passed) = " << eBhabha);
361 B2DEBUG(22,
"Cutflow: Trigger cut passed: index = " << cutIndexPassed);
371 if (ECLchannelMapHasChanged) {
372 B2INFO(
"ECLBhabhaTCollectorModule::collect() " <<
LogVar(
"ECLchannelMapHasChanged", ECLchannelMapHasChanged));
374 B2FATAL(
"ECLBhabhaTCollectorModule::collect() : Can't initialize eclChannelMapper!");
398 B2DEBUG(29,
"Finished checking if previous crystal time payload has changed");
403 B2DEBUG(29,
"Finished checking if previous crate time payload has changed");
404 B2DEBUG(29,
"m_CrateTime size = " <<
m_CrateTime.size());
408 B2DEBUG(29,
"Finished checking if previous crate time payload has changed");
412 B2DEBUG(29,
"Finished checking if reference crystal cell ids payload has changed");
415 B2DEBUG(25,
"ECLBhabhaTCollector:: loaded ECLCrystalTimeOffset from the database"
418 B2DEBUG(25,
"ECLBhabhaTCollector:: loaded ECLCrateTimeOffset from the database"
421 B2DEBUG(25,
"ECLBhabhaTCollector:: loaded ECLCrystalElectronics from the database"
424 B2DEBUG(25,
"ECLBhabhaTCollector:: loaded ECLCrystalElectronicsTime from the database"
427 B2DEBUG(25,
"ECLBhabhaTCollector:: loaded ECLCrystalFlightTime from the database"
430 B2DEBUG(25,
"ECLBhabhaTCollector:: loaded ECLReferenceCrystalPerCrateCalib from the database"
442 vector<float> Crate_time_ns(52, 0.0);
447 Crate_time_ns[crateID_temp - 1] =
m_CrateTime[crysID] * TICKS_TO_NS;
457 for (
int crateID_temp = 1; crateID_temp <= 52; crateID_temp++) {
470 B2INFO(
"ECLBhabhaTCollector:: ECLCrystalTimeOffset from the database information:"
473 B2INFO(
"First event so print out previous ts values");
483 for (
int crateID_temp = 1; crateID_temp <= 52; crateID_temp++) {
484 getObjectPtr<TH1F>(
"tcrateDatabase_ns")->SetBinContent(crateID_temp + 0.001, Crate_time_ns[crateID_temp - 1]);
501 double evt_t0_unc = -1;
502 double evt_t0_ECL_closestSVD = -1;
503 double evt_t0_ECL_minChi2 = -1;
509 }
else if (!
m_eventT0->hasTemporaryEventT0(Const::EDetector::SVD)) {
516 B2DEBUG(22,
"Cutflow: Event t0 exists: index = " << cutIndexPassed);
519 const auto bestSVDEventT0Candidate =
m_eventT0->getBestSVDTemporaryEventT0();
520 evt_t0 = bestSVDEventT0Candidate->eventT0;
521 evt_t0_unc = bestSVDEventT0Candidate->eventT0Uncertainty;
529 if (
m_eventT0->hasTemporaryEventT0(Const::EDetector::ECL)) {
530 vector<EventT0::EventT0Component> evt_t0_list_ECL =
m_eventT0->getTemporaryEventT0s(Const::EDetector::ECL);
533 double smallest_SVD_ECL_t0_diff = fabs(evt_t0_list_ECL[0].eventT0 - evt_t0);
534 int smallest_SVD_ECL_t0_diff_idx = 0;
535 for (
long unsigned int ECLi = 0; ECLi < evt_t0_list_ECL.size(); ECLi++) {
536 double tempt_ECL_t0 = evt_t0_list_ECL[ECLi].eventT0;
537 if (fabs(tempt_ECL_t0 - evt_t0) < smallest_SVD_ECL_t0_diff) {
538 smallest_SVD_ECL_t0_diff = fabs(tempt_ECL_t0 - evt_t0);
539 smallest_SVD_ECL_t0_diff_idx = ECLi;
543 evt_t0_ECL_closestSVD = evt_t0_list_ECL[smallest_SVD_ECL_t0_diff_idx].eventT0;
544 B2DEBUG(26,
"evt_t0_ECL_closestSVD = " << evt_t0_ECL_closestSVD);
548 double smallest_ECL_t0_minChi2 = evt_t0_list_ECL[0].quality;
549 int smallest_ECL_t0_minChi2_idx = 0;
551 B2DEBUG(26,
"evt_t0_list_ECL[0].quality = " << evt_t0_list_ECL[0].quality
552 <<
", with ECL event t0 = " << evt_t0_list_ECL[0].eventT0);
554 for (
long unsigned int ECLi = 0; ECLi < evt_t0_list_ECL.size(); ECLi++) {
555 B2DEBUG(26,
"evt_t0_list_ECL[" << ECLi <<
"].quality = " << evt_t0_list_ECL[ECLi].quality
556 <<
", with ECL event t0 = " <<
557 evt_t0_list_ECL[ECLi].eventT0);
558 if (evt_t0_list_ECL[ECLi].quality < smallest_ECL_t0_minChi2) {
559 smallest_ECL_t0_minChi2 = evt_t0_list_ECL[ECLi].quality;
560 smallest_ECL_t0_minChi2_idx = ECLi;
564 evt_t0_ECL_minChi2 = evt_t0_list_ECL[smallest_ECL_t0_minChi2_idx].eventT0;
566 B2DEBUG(26,
"evt_t0_ECL_minChi2 = " << evt_t0_ECL_minChi2);
567 B2DEBUG(26,
"smallest_ECL_t0_minChi2_idx = " << smallest_ECL_t0_minChi2_idx);
585 int tempCrysID = eclCalDigit.getCellId() - 1;
586 m_EperCrys[tempCrysID] = eclCalDigit.getEnergy();
593 int tempCrysID = eclDigit.getCellId() - 1;
608 double maxp[2] = {0., 0.};
609 int maxiTrk[2] = { -1, -1};
610 int nTrkAll =
tracks.getEntries();
618 for (
int iTrk = 0; iTrk < nTrkAll; iTrk++) {
622 if (not tempTrackFit) {
continue;}
626 double z0 = tempTrackFit->
getZ0();
627 double d0 = tempTrackFit->
getD0();
682 if (charge > 0) {icharge = 1;}
683 if (p > maxp[icharge]) {
685 maxiTrk[icharge] = iTrk;
694 if (nTrkTight != 2) {
700 B2DEBUG(22,
"Cutflow: Two tight tracks: index = " << cutIndexPassed);
703 if (nTrkLoose != 2) {
709 B2DEBUG(22,
"Cutflow: No additional loose tracks: index = " << cutIndexPassed);
715 bool oppositelyChargedTracksPassed = maxiTrk[0] != -1 && maxiTrk[1] != -1;
716 if (!oppositelyChargedTracksPassed) {
722 B2DEBUG(22,
"Cutflow: Oppositely charged tracks: index = " << cutIndexPassed);
731 double trkEClustLab[2] = {0., 0.};
732 double trkEClustCOM[2] = {0., 0.};
735 ROOT::Math::PxPyPzEVector trkp4Lab[2];
736 ROOT::Math::PxPyPzEVector trkp4COM[2];
739 int crysIDMax[2] = { -1, -1 };
740 double crysEMax[2] = { -1, -1 };
741 double crysE2Max[2] = { -1, -1 };
742 int numClustersPerTrack[2] = { 0, 0 };
744 double clusterTime[2] = {0, 0};
748 vector<double> time_ECLCaldigits_bothClusters;
749 vector<int> cid_ECLCaldigits_bothClusters;
750 vector<double> E_ECLCaldigits_bothClusters;
751 vector<double> amp_ECLDigits_bothClusters;
752 vector<int> chargeID_ECLCaldigits_bothClusters;
754 for (
int icharge = 0; icharge < 2; icharge++) {
755 if (maxiTrk[icharge] > -1) {
756 B2DEBUG(22,
"looping over the 2 max pt tracks");
759 if (not tempTrackFit) {
continue;}
761 trkp4COM[icharge] = boostrotate.
rotateLabToCms() * trkp4Lab[icharge];
762 trkpLab[icharge] = trkp4Lab[icharge].P();
763 trkpCOM[icharge] = trkp4COM[icharge].P();
769 auto eclClusterRelationsFromTracks =
tracks[maxiTrk[icharge]]->getRelationsTo<
ECLCluster>();
770 for (
unsigned int clusterIdx = 0; clusterIdx < eclClusterRelationsFromTracks.size(); clusterIdx++) {
772 B2DEBUG(22,
"Looking at clusters. index = " << clusterIdx);
773 auto cluster = eclClusterRelationsFromTracks[clusterIdx];
774 bool goodClusterType =
false;
778 goodClusterType =
true;
779 numClustersPerTrack[icharge]++;
782 if (goodClusterType) {
784 clusterTime[icharge] = cluster->getTime();
786 auto eclClusterRelations = cluster->getRelationsTo<
ECLCalDigit>(
"ECLCalDigits");
789 for (
unsigned int ir = 0; ir < eclClusterRelations.size(); ir++) {
790 const auto calDigit = eclClusterRelations.object(ir);
791 int tempCrysID = calDigit->getCellId() - 1;
798 if (tempE > crysEMax[icharge]) {
800 crysE2Max[icharge] = crysEMax[icharge];
802 crysEMax[icharge] = tempE;
803 crysIDMax[icharge] = tempCrysID;
806 if (tempE > crysE2Max[icharge] && tempCrysID != crysIDMax[icharge]) {
807 crysE2Max[icharge] = tempE;
814 B2DEBUG(26,
"calDigit(ir" << ir <<
") time = " << calDigit->getTime() <<
"ns , with E = " << tempE <<
" GeV");
815 time_ECLCaldigits_bothClusters.push_back(calDigit->getTime());
816 cid_ECLCaldigits_bothClusters.push_back(tempCrysID);
817 E_ECLCaldigits_bothClusters.push_back(tempE);
818 amp_ECLDigits_bothClusters.push_back(ecl_dig->
getAmp());
819 chargeID_ECLCaldigits_bothClusters.push_back(icharge);
824 trkEClustCOM[icharge] = trkEClustLab[icharge] * trkpCOM[icharge] / trkpLab[icharge];
829 E_DIV_p[icharge] = trkEClustCOM[icharge] / trkpCOM[icharge];
841 int numCrystalsPassingCuts = 0;
843 int crystalIDs[2] = { -1, -1};
844 int crateIDs[2] = { -1, -1};
845 double ts_prevCalib[2] = { -1, -1};
846 double tcrate_prevCalib[2] = { -1, -1};
847 double times_noPreviousCalibrations[2] = { -1, -1};
848 bool crystalCutsPassed[2] = {
false,
false};
849 double crystalEnergies[2] = { -1, -1};
850 double crystalEnergies2[2] = { -1, -1};
852 for (
int iCharge = 0; iCharge < 2; iCharge++) {
853 int crystal_idx = crysIDMax[iCharge];
863 auto amplitude = ecl_dig->
getAmp();
864 crystalEnergies[iCharge] = en;
867 double time = ecl_dig->
getTimeFit() * TICKS_TO_NS - evt_t0;
876 double energyTimeShift =
m_ECLTimeUtil->energyDependentTimeOffsetElectronic(amplitude *
m_Electronics[cid - 1]) * TICKS_TO_NS;
878 B2DEBUG(29,
"cellid = " << cid <<
", amplitude = " << amplitude <<
", time before t(E) shift = " << time <<
879 ", t(E) shift = " << energyTimeShift <<
" ns");
880 time -= energyTimeShift;
893 crystalIDs[iCharge] = cid;
898 tcrate_prevCalib[iCharge] =
m_CrateTime[cid - 1] * TICKS_TO_NS;
899 times_noPreviousCalibrations[iCharge] = time;
902 B2DEBUG(26,
"iCharge = " << iCharge);
903 B2DEBUG(26,
"crateIDs[iCharge] = " << crateIDs[iCharge]);
904 B2DEBUG(26,
"times_noPreviousCalibrations[iCharge] = " << times_noPreviousCalibrations[iCharge]);
905 B2DEBUG(26,
"tcrate_prevCalib[iCharge] = " << tcrate_prevCalib[iCharge]);
906 B2DEBUG(26,
"ts_prevCalib[iCharge] = " << ts_prevCalib[iCharge]);
909 crystalCutsPassed[iCharge] =
true;
913 crystalEnergies2[iCharge] = crysE2Max[iCharge];
951 bool E_DIV_p_instance_passed[2] = {
false,
false};
952 double E_DIV_p_CUT = 0.7;
953 for (
int icharge = 0; icharge < 2; icharge++) {
954 E_DIV_p_instance_passed[icharge] = E_DIV_p[icharge] > E_DIV_p_CUT;
956 if (!E_DIV_p_instance_passed[0] || !E_DIV_p_instance_passed[1]) {
962 B2DEBUG(22,
"Cutflow: E_i/p_i > " << E_DIV_p_CUT <<
": index = " << cutIndexPassed);
968 double invMassTrk = (trkp4Lab[0] + trkp4Lab[1]).M();
969 double invMass_CUT = 0.9;
978 bool invMassCutsPassed = invMassTrk > (invMass_CUT * boostrotate.
getCMSEnergy());
979 if (!invMassCutsPassed) {
985 B2DEBUG(22,
"Cutflow: m(track 1+2) > " << invMass_CUT <<
"*E_COM = " << invMass_CUT <<
" * " << boostrotate.
getCMSEnergy() <<
986 " : index = " << cutIndexPassed);
991 for (
int iCharge = 0; iCharge < 2; iCharge++) {
992 if (crystalCutsPassed[iCharge]) {
993 getObjectPtr<TH2F>(
"TimevsCrysPrevCrateCalibPrevCrystCalib")->Fill((crystalIDs[iCharge]) + 0.001,
994 times_noPreviousCalibrations[iCharge] - ts_prevCalib[iCharge] - tcrate_prevCalib[iCharge], 1);
995 getObjectPtr<TH2F>(
"TimevsCratePrevCrateCalibPrevCrystCalib")->Fill((crateIDs[iCharge]) + 0.001,
996 times_noPreviousCalibrations[iCharge] - ts_prevCalib[iCharge] - tcrate_prevCalib[iCharge], 1);
997 getObjectPtr<TH2F>(
"TimevsCrysNoCalibrations")->Fill((crystalIDs[iCharge]) + 0.001, times_noPreviousCalibrations[iCharge], 1);
998 getObjectPtr<TH2F>(
"TimevsCrateNoCalibrations")->Fill((crateIDs[iCharge]) + 0.001, times_noPreviousCalibrations[iCharge], 1);
999 getObjectPtr<TH2F>(
"TimevsCrysPrevCrateCalibNoCrystCalib")->Fill((crystalIDs[iCharge]) + 0.001,
1000 times_noPreviousCalibrations[iCharge] - tcrate_prevCalib[iCharge], 1);
1001 getObjectPtr<TH2F>(
"TimevsCrateNoCrateCalibPrevCrystCalib")->Fill((crateIDs[iCharge]) + 0.001,
1002 times_noPreviousCalibrations[iCharge] - ts_prevCalib[iCharge], 1);
1005 numCrystalsPassingCuts++;
1013 if (crystalCutsPassed[0] || crystalCutsPassed[1]) {
1017 B2DEBUG(22,
"Cutflow: At least one crystal time and quality cuts passed: index = " << cutIndexPassed);
1024 for (
int iCharge = 0; iCharge < 2; iCharge++) {
1025 if (crystalCutsPassed[iCharge]) {
1030 m_tree_time = times_noPreviousCalibrations[iCharge];
1034 m_tree_E1Etot = crystalEnergies[iCharge] / trkEClustLab[iCharge];
1035 m_tree_E1E2 = crystalEnergies[iCharge] / crystalEnergies2[iCharge];
1036 m_tree_E1p = crystalEnergies[iCharge] / trkpLab[iCharge];
1055 if (crystalCutsPassed[0] && crystalCutsPassed[1] &&
1056 numClustersPerTrack[0] == 1 && numClustersPerTrack[1] == 1) {
1077 B2DEBUG(26,
"SVD evt_t0 = " << evt_t0);
1082 for (
long unsigned int digit_i = 0; digit_i < time_ECLCaldigits_bothClusters.size(); digit_i++) {
1092 ts_prevCalib[chargeID_ECLCaldigits_bothClusters[digit_i]] -
1093 tcrate_prevCalib[chargeID_ECLCaldigits_bothClusters[digit_i]];
1094 m_tree_cid = cid_ECLCaldigits_bothClusters[digit_i];
1105 B2DEBUG(26,
"This was for event number = " <<
m_tree_evtNum);