104 "Validating crystal and crate time calibrations using electron clusters in events with lots of tracks and clusters") ;
117 "Validating crystal and crate time calibrations using electron clusters in events with lots of tracks and clusters") ;
130 m_dbg_tree_run =
new TTree(
"tree_run",
"Storing crate time constants") ;
141 B2INFO(
"eclBhabhaTimeCalibrationValidationCollector: Experiment = " <<
m_EventMetaData->getExperiment() <<
147 double binSize = 2000.0 / pow(2, 12);
148 double halfBinSize = binSize / 2.0;
152 double nBinsNeg = floor((
m_timeAbsMax - halfBinSize) / binSize);
153 double min_t = -nBinsNeg * binSize - halfBinSize;
154 int nbins = nBinsNeg * 2 + 1;
155 double max_t = min_t + nbins * binSize;
159 const Int_t N_E_BIN_EDGES = 64;
160 const Int_t N_E_BINS = N_E_BIN_EDGES - 1;
161 Double_t energyBinEdges[N_E_BIN_EDGES] = {0, 0.05, 0.051, 0.052, 0.053, 0.054, 0.055, 0.056, 0.057, 0.058, 0.059, 0.06, 0.062, 0.064, 0.066, 0.068, 0.07, 0.075, 0.08, 0.085, 0.09, 0.095, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.25, 2.5, 2.8, 3.2, 3.6, 4, 4.4, 4.8, 5.2, 5.6, 6, 6.4, 6.8, 7.2, 7.6, 8};
164 auto cutflow =
new TH1F(
"cutflow",
" ;Cut label number ;Number of events passing cut", 10, 0, 10) ;
167 auto clusterTime =
new TH1F(
"clusterTime",
";Electron ECL cluster time [ns]; number of ECL clusters", nbins, min_t, max_t) ;
170 auto clusterTime_cid =
new TH2F(
"clusterTime_cid",
175 auto clusterTime_run =
new TH2F(
"clusterTime_run",
176 ";Run number ;Electron ECL cluster time [ns]", 7000, 0, 7000, nbins, min_t, max_t) ;
180 auto clusterTimeClusterE =
new TH2F(
"clusterTimeClusterE",
181 ";Electron cluster energy [GeV];Electron cluster time [ns]", N_E_BINS, energyBinEdges, nbins, min_t, max_t) ;
184 auto dt99_clusterE =
new TH2F(
"dt99_clusterE",
185 ";Electron cluster energy [GeV];dt99 [ns]", N_E_BINS, energyBinEdges, nbins, 0, max_t) ;
189 auto eventT0 =
new TH1F(
"eventT0",
";event t0 [ns]; number of events", nbins, min_t, max_t) ;
192 auto eventT0Detector =
new TH1F(
"eventT0Detector",
193 "detector used for eventT0 (SVD=2, CDC=4, TOP=8, ECL=32);detector number; number of events", 48, 0, 48) ;
196 auto clusterTimeE0E1diff =
new TH1F(
"clusterTimeE0E1diff",
197 ";ECL cluster time of max E electron - ECL cluster time of 2nd max E electron [ns]; number of electron ECL cluster time differences",
198 nbins, min_t, max_t) ;
214 int cutIndexPassed = 0;
216 B2DEBUG(22,
"Cutflow: no cuts: index = " << cutIndexPassed);
225 B2WARNING(
"SoftwareTriggerResult required to select bhabha event is not found");
231 const std::map<std::string, int>& fresults =
m_TrgResult->getResults();
232 if (fresults.find(
"software_trigger_cut&skim&accept_bhabha") == fresults.end()) {
233 B2WARNING(
"Can't find required bhabha trigger identifier");
237 const bool eBhabha = (
m_TrgResult->getResult(
"software_trigger_cut&skim&accept_bhabha") ==
239 B2DEBUG(22,
"eBhabha (trigger passed) = " << eBhabha);
250 B2DEBUG(22,
"Cutflow: Trigger cut passed: index = " << cutIndexPassed);
262 if (ECLchannelMapHasChanged) {
263 B2INFO(
"eclBhabhaTimeCalibrationValidationCollectorModule::collect() " <<
LogVar(
"ECLchannelMapHasChanged",
264 ECLchannelMapHasChanged));
266 B2FATAL(
"eclBhabhaTimeCalibrationValidationCollectorModule::collect() : Can't initialize eclChannelMapper!");
274 B2DEBUG(29,
"Finished checking if previous crystal time payload has changed");
281 B2DEBUG(25,
"eclBhabhaTimeCalibrationValidationCollector:: loaded ECLCrateTimeOffset from the database"
291 vector<float> Crate_time_ns(52, 0.0);
292 vector<float> Crate_time_unc_ns(52, 0.0);
297 Crate_time_ns[crateID_temp - 1] =
m_CrateTime[crysID] * TICKS_TO_NS;
298 Crate_time_unc_ns[crateID_temp - 1] =
m_CrateTimeUnc[crysID] * TICKS_TO_NS;
306 int tempCrysID = eclCalDigit.getCellId() - 1;
307 m_EperCrys[tempCrysID] = eclCalDigit.getEnergy();
313 double evt_t0 = -1000 ;
314 double evt_t0_unc = -1000 ;
315 int evt_t0_detector = 0;
327 evt_t0_unc =
m_eventT0->getEventT0Uncertainty() ;
328 if (
m_eventT0->isSVDEventT0()) {evt_t0_detector += 2;}
329 if (
m_eventT0->isCDCEventT0()) {evt_t0_detector += 4;}
330 if (
m_eventT0->isTOPEventT0()) {evt_t0_detector += 8;}
331 if (
m_eventT0->isECLEventT0()) {evt_t0_detector += 32;}
334 B2DEBUG(26,
"Found event t0") ;
345 double maxp[2] = {0., 0.};
346 int maxiTrk[2] = { -1, -1};
347 int nTrkAll =
tracks.getEntries() ;
356 for (
int iTrk = 0 ; iTrk < nTrkAll ; iTrk++) {
360 if (not tempTrackFit) {continue ;}
364 double z0 = tempTrackFit->
getZ0() ;
365 double d0 = tempTrackFit->
getD0() ;
404 if (charge > 0) {icharge = 1;}
405 if (p > maxp[icharge]) {
407 maxiTrk[icharge] = iTrk;
416 if (nTrkTight != 2) {
422 B2DEBUG(22,
"Cutflow: Two tight tracks: index = " << cutIndexPassed);
425 if (nTrkLoose != 2) {
431 B2DEBUG(22,
"Cutflow: No additional loose tracks: index = " << cutIndexPassed) ;
436 bool oppositelyChargedTracksPassed = maxiTrk[0] != -1 && maxiTrk[1] != -1;
437 if (!oppositelyChargedTracksPassed) {
443 B2DEBUG(22,
"Cutflow: Oppositely charged tracks: index = " << cutIndexPassed);
451 double trkEClustLab[2] = {0., 0.};
452 double trkEClustCOM[2] = {0., 0.};
455 ROOT::Math::PxPyPzEVector trkp4Lab[2];
456 ROOT::Math::PxPyPzEVector trkp4COM[2];
459 int numClustersPerTrack[2] = { 0, 0 };
462 vector<double> goodClustTimes ;
463 vector<double> goodClust_dt99 ;
464 vector<double> goodClustE ;
465 vector<int> goodClustMaxEcrys_cid ;
467 for (
int icharge = 0; icharge < 2; icharge++) {
468 if (maxiTrk[icharge] > -1) {
469 B2DEBUG(22,
"looping over the 2 max pt tracks");
472 if (not tempTrackFit) {continue ;}
475 trkp4COM[icharge] = boostrotate.
rotateLabToCms() * trkp4Lab[icharge];
476 trkpLab[icharge] = trkp4Lab[icharge].P();
477 trkpCOM[icharge] = trkp4COM[icharge].P();
483 auto eclClusterRelationsFromTracks =
tracks[maxiTrk[icharge]]->getRelationsTo<
ECLCluster>();
484 for (
unsigned int clusterIdx = 0; clusterIdx < eclClusterRelationsFromTracks.size(); clusterIdx++) {
486 B2DEBUG(22,
"Looking at clusters. index = " << clusterIdx);
487 auto cluster = eclClusterRelationsFromTracks[clusterIdx];
490 numClustersPerTrack[icharge]++;
492 double electronTime = cluster->getTime();
493 bool badElectronTime = cluster->hasFailedFitTime();
494 bool badElectronTimeResolution = cluster->hasFailedTimeResolution();
496 (!badElectronTime) &&
497 (!badElectronTimeResolution)) {
498 trkEClustLab[icharge] += eClust ;
499 short cid = cluster->getMaxECellId() ;
500 goodClustMaxEcrys_cid.push_back(cid) ;
501 goodClustTimes.push_back(electronTime) ;
502 goodClust_dt99.push_back(cluster->getDeltaTime99()) ;
503 goodClustE.push_back(eClust);
507 trkEClustCOM[icharge] = trkEClustLab[icharge] * trkpCOM[icharge] / trkpLab[icharge];
512 E_DIV_p[icharge] = trkEClustCOM[icharge] / trkpCOM[icharge];
519 B2DEBUG(26,
"Extracted time information and E/p for the tracks") ;
530 long unsigned int numGoodElectronClusters_cut = 2 ;
531 if (goodClustTimes.size() != numGoodElectronClusters_cut) {
537 B2DEBUG(22,
"Cutflow: Exactly " << numGoodElectronClusters_cut
538 <<
" good clusters connected to tracks: index = " << cutIndexPassed);
544 bool E_DIV_p_instance_passed[2] = {
false,
false};
545 double E_DIV_p_CUT = 0.7;
546 for (
int icharge = 0; icharge < 2; icharge++) {
547 E_DIV_p_instance_passed[icharge] = E_DIV_p[icharge] > E_DIV_p_CUT;
549 if (!E_DIV_p_instance_passed[0] || !E_DIV_p_instance_passed[1]) {
555 B2DEBUG(22,
"Cutflow: E_i/p_i > " << E_DIV_p_CUT <<
": index = " << cutIndexPassed);
560 double invMassTrk = (trkp4Lab[0] + trkp4Lab[1]).M();
561 double invMass_CUT = 0.9;
563 bool invMassCutsPassed = invMassTrk > (invMass_CUT * boostrotate.
getCMSEnergy());
564 if (!invMassCutsPassed) {
570 B2DEBUG(22,
"Cutflow: m(track 1+2) > " << invMass_CUT <<
"*E_COM = " << invMass_CUT <<
" * " << boostrotate.
getCMSEnergy() <<
571 " : index = " << cutIndexPassed);
574 B2DEBUG(22,
"Event passed all cuts");
581 bool isSVDt0 =
m_eventT0->isSVDEventT0();
582 bool isCDCt0 =
m_eventT0->isCDCEventT0();
583 bool isECLt0 =
m_eventT0->isECLEventT0();
584 string t0Detector =
"UNKNOWN... WHY?";
587 }
else if (isCDCt0) {
589 }
else if (isECLt0) {
593 B2DEBUG(26,
"t0 = " << evt_t0 <<
" ns. t0 is from SVD?=" << isSVDt0 <<
", t0 is from CDC?=" << isCDCt0 <<
", t0 is from ECL?=" <<
594 isECLt0 <<
" t0 from " <<
599 for (
long unsigned int i = 0 ; i < goodClustTimes.size() ; i++) {
601 getObjectPtr<TH2F>(
"clusterTime_cid")->Fill(goodClustMaxEcrys_cid[i] + 0.001, goodClustTimes[i], 1) ;
603 getObjectPtr<TH2F>(
"clusterTimeClusterE")->Fill(goodClustE[i], goodClustTimes[i], 1) ;
624 B2DEBUG(26,
"Filled cluster tree") ;
628 if (goodClustE[0] > goodClustE[1]) {
629 tDiff = goodClustTimes[0] - goodClustTimes[1];
631 tDiff = goodClustTimes[1] - goodClustTimes[0];
654 for (
int icrate = 1; icrate <= 52; icrate++) {
666 B2DEBUG(26,
"Filled event tree") ;