27 #include <ecl/modules/eclBhabhaTimeCalibrationValidationCollector/eclBhabhaTimeCalibrationValidationCollectorModule.h>
28 #include <framework/dataobjects/EventMetaData.h>
29 #include <framework/gearbox/Const.h>
30 #include <ecl/dbobjects/ECLCrystalCalib.h>
31 #include <ecl/dataobjects/ECLDigit.h>
32 #include <ecl/dataobjects/ECLCalDigit.h>
33 #include <ecl/dataobjects/ECLTrig.h>
34 #include <mdst/dataobjects/ECLCluster.h>
35 #include <mdst/dataobjects/Track.h>
36 #include <mdst/dataobjects/HitPatternCDC.h>
37 #include <tracking/dataobjects/RecoTrack.h>
38 #include <ecl/digitization/EclConfiguration.h>
40 #include <analysis/utility/PCmsLabTransform.h>
41 #include <analysis/ClusterUtility/ClusterUtils.h>
42 #include <boost/optional.hpp>
55 REG_MODULE(eclBhabhaTimeCalibrationValidationCollector)
63 m_dbg_tree_electronClusters(0),
66 m_CrateTimeDB("ECLCrateTimeOffset")
68 setDescription(
"This module validates the ECL cluster times");
70 addParam(
"timeAbsMax", m_timeAbsMax,
71 "Events with fabs(getTimeFit) > m_timeAbsMax "
72 "are excluded", (
short)80);
74 addParam(
"saveTree", m_saveTree,
75 "If true, TTree 'tree' with more detailed event info is saved in "
76 "the output file specified by HistoManager",
79 addParam(
"looseTrkZ0", m_looseTrkZ0,
"max Z0 for loose tracks (cm)", 10.);
80 addParam(
"tightTrkZ0", m_tightTrkZ0,
"max Z0 for tight tracks (cm)", 2.);
81 addParam(
"looseTrkD0", m_looseTrkD0,
"max D0 for loose tracks (cm)", 2.);
82 addParam(
"tightTrkD0", m_tightTrkD0,
"max D0 for tight tracks (cm)", 0.5);
86 setPropertyFlags(c_ParallelProcessingCertified);
103 m_dbg_tree_electronClusters =
new TTree(
"tree_electronClusters",
104 "Validating crystal and crate time calibrations using electron clusters in events with lots of tracks and clusters") ;
105 m_dbg_tree_electronClusters->Branch(
"EventNum" , &m_tree_evt_num) ->SetTitle(
"Event number") ;
106 m_dbg_tree_electronClusters->Branch(
"cluster_time" , &m_tree_time) ->SetTitle(
"Cluster time t (calibrated), ns") ;
107 m_dbg_tree_electronClusters->Branch(
"clust_E" , &m_E_electron_clust) ->SetTitle(
"Electron type cluster energy, GeV") ;
108 m_dbg_tree_electronClusters->Branch(
"t0" , &m_tree_t0) ->SetTitle(
"T0, ns") ;
109 m_dbg_tree_electronClusters->Branch(
"t0_unc" , &m_tree_t0_unc) ->SetTitle(
"T0 uncertainty, ns") ;
110 m_dbg_tree_electronClusters->Branch(
"runNum" , &m_tree_run) ->SetTitle(
"Run number") ;
111 m_dbg_tree_electronClusters->Branch(
"CrystalCellID" , &m_tree_cid) ->SetTitle(
"Cell ID, 1..8736") ;
112 m_dbg_tree_electronClusters->Branch(
"dt99" , &m_tree_dt99) ->SetTitle(
"Cluster dt99, ns") ;
113 m_dbg_tree_electronClusters->SetAutoSave(10) ;
116 m_dbg_tree_event =
new TTree(
"tree_event",
117 "Validating crystal and crate time calibrations using electron clusters in events with lots of tracks and clusters") ;
118 m_dbg_tree_event->Branch(
"EventNum" , &m_tree_evt_num) ->SetTitle(
"Event number") ;
119 m_dbg_tree_event->Branch(
"t0" , &m_tree_t0) ->SetTitle(
"T0, ns") ;
120 m_dbg_tree_event->Branch(
"t0_unc" , &m_tree_t0_unc) ->SetTitle(
"T0 uncertainty, ns") ;
121 m_dbg_tree_event->Branch(
"runNum" , &m_tree_run) ->SetTitle(
"Run number") ;
122 m_dbg_tree_event->Branch(
"E0" , &m_tree_E0) ->SetTitle(
"Highest E cluster E") ;
123 m_dbg_tree_event->Branch(
"E1" , &m_tree_E1) ->SetTitle(
"2nd highest E cluster E") ;
124 m_dbg_tree_event->Branch(
"time_E0" , &m_tree_time_fromE0) ->SetTitle(
"Cluster time of highest E cluster") ;
125 m_dbg_tree_event->Branch(
"time_E1" , &m_tree_time_fromE1) ->SetTitle(
"Cluster time of 2nd highest E cluster") ;
126 m_dbg_tree_event->SetAutoSave(10) ;
130 m_dbg_tree_run =
new TTree(
"tree_run",
"Storing crate time constants") ;
131 m_dbg_tree_run->Branch(
"runNum" , &m_tree_run) ->SetTitle(
"Run number") ;
132 m_dbg_tree_run->Branch(
"crateid" , &m_tree_crateid) ->SetTitle(
"Crate ID") ;
133 m_dbg_tree_run->Branch(
"tcrate" , &m_tree_tcrate) ->SetTitle(
"Crate time") ;
134 m_dbg_tree_run->Branch(
"tcrate_unc" , &m_tree_tcrate_unc) ->SetTitle(
"Crate time uncertainty") ;
135 m_dbg_tree_run->SetAutoSave(10) ;
142 B2INFO(
"eclBhabhaTimeCalibrationValidationCollector: Experiment = " << evtMetaData->getExperiment() <<
143 " run = " << evtMetaData->getRun()) ;
148 double binSize = 2000.0 / pow(2, 12);
149 double halfBinSize = binSize / 2.0;
153 double nBinsNeg = floor((m_timeAbsMax - halfBinSize) / binSize);
154 double min_t = -nBinsNeg * binSize - halfBinSize;
155 int nbins = nBinsNeg * 2 + 1;
156 double max_t = min_t + nbins * binSize;
160 const Int_t N_E_BIN_EDGES = 64;
161 const Int_t N_E_BINS = N_E_BIN_EDGES - 1;
162 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};
165 auto cutflow =
new TH1F(
"cutflow",
" ;Cut label number ;Number of events passing cut", 10, 0, 10) ;
166 registerObject<TH1F>(
"cutflow", cutflow) ;
168 auto clusterTime =
new TH1F(
"clusterTime",
";Electron ECL cluster time [ns]; number of ECL clusters", nbins, min_t, max_t) ;
169 registerObject<TH1F>(
"clusterTime", clusterTime) ;
171 auto clusterTime_cid =
new TH2F(
"clusterTime_cid",
172 ";crystal Cell ID ;Electron ECL cluster time [ns]", 8736, 1, 8736 + 1, nbins, min_t, max_t) ;
173 registerObject<TH2F>(
"clusterTime_cid", 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) ;
177 registerObject<TH2F>(
"clusterTime_run", clusterTime_run) ;
180 auto clusterTimeClusterE =
new TH2F(
"clusterTimeClusterE",
181 ";Electron cluster energy [GeV];Electron cluster time [ns]", N_E_BINS, energyBinEdges, nbins, min_t, max_t) ;
182 registerObject<TH2F>(
"clusterTimeClusterE", clusterTimeClusterE) ;
184 auto dt99_clusterE =
new TH2F(
"dt99_clusterE",
185 ";Electron cluster energy [GeV];dt99 [ns]", N_E_BINS, energyBinEdges, nbins, 0, max_t) ;
186 registerObject<TH2F>(
"dt99_clusterE", dt99_clusterE) ;
189 auto eventT0 =
new TH1F(
"eventT0",
";event t0 [ns]; number of events", nbins, min_t, max_t) ;
190 registerObject<TH1F>(
"eventT0", eventT0) ;
192 auto clusterTimeE0E1diff =
new TH1F(
"clusterTimeE0E1diff",
193 ";ECL cluster time of max E electron - ECL cluster time of 2nd max E electron [ns]; number of electron ECL cluster time differences",
194 nbins, min_t, max_t) ;
195 registerObject<TH1F>(
"clusterTimeE0E1diff", clusterTimeE0E1diff) ;
199 tracks.isRequired() ;
200 m_eclClusterArray.isRequired() ;
201 m_eclCalDigitArray.isRequired() ;
207 int cutIndexPassed = 0;
208 getObjectPtr<TH1F>(
"cutflow")->Fill(cutIndexPassed);
209 B2DEBUG(22,
"Cutflow: no cuts: index = " << cutIndexPassed);
219 B2DEBUG(29,
"Finished checking if previous crystal time payload has changed");
221 if (m_CrateTimeDB.hasChanged()) {
222 m_CrateTime = m_CrateTimeDB->getCalibVector();
223 m_CrateTimeUnc = m_CrateTimeDB->getCalibUncVector();
226 B2DEBUG(25,
"eclBhabhaTimeCalibrationValidationCollector:: loaded ECLCrateTimeOffset from the database"
227 <<
LogVar(
"IoV", m_CrateTimeDB.getIoV())
228 <<
LogVar(
"Checksum", m_CrateTimeDB.getChecksum()));
236 vector<float> Crate_time_ns(52, 0.0);
237 vector<float> Crate_time_unc_ns(52, 0.0);
240 for (
int crysID = 1; crysID <= 8736; crysID++) {
241 int crateID_temp = crystalMapper->
getCrateID(crysID);
242 Crate_time_ns[crateID_temp - 1] = m_CrateTime[crysID] * TICKS_TO_NS;
243 Crate_time_unc_ns[crateID_temp - 1] = m_CrateTimeUnc[crysID] * TICKS_TO_NS;
249 m_EperCrys.resize(8736);
250 for (
auto& eclCalDigit : m_eclCalDigitArray) {
251 int tempCrysID = eclCalDigit.getCellId() - 1;
252 m_EperCrys[tempCrysID] = eclCalDigit.getEnergy();
258 double evt_t0 = -1000 ;
259 double evt_t0_unc = -1000 ;
262 if (m_eventT0.isOptional()) {
263 if (!m_eventT0.isValid()) {
266 if (!m_eventT0->hasEventT0()) {
270 evt_t0 = m_eventT0->getEventT0() ;
271 evt_t0_unc = m_eventT0->getEventT0Uncertainty() ;
273 B2DEBUG(26,
"Found event t0") ;
286 double maxp[2] = {0., 0.};
287 int maxiTrk[2] = { -1, -1};
288 int nTrkAll = tracks.getEntries() ;
297 for (
int iTrk = 0 ; iTrk < nTrkAll ; iTrk++) {
301 if (not tempTrackFit) {continue ;}
304 short charge = tempTrackFit->getChargeSign() ;
305 double z0 = tempTrackFit->getZ0() ;
306 double d0 = tempTrackFit->getD0() ;
307 int nCDChits = tempTrackFit->getHitPatternCDC().getNHits() ;
308 double p = tempTrackFit->getMomentum().Mag() ;
313 if (fabs(d0) > m_looseTrkD0) {
316 if (fabs(z0) > m_looseTrkZ0) {
334 if (fabs(d0) > m_tightTrkD0) {
337 if (fabs(z0) > m_tightTrkZ0) {
345 if (charge > 0) {icharge = 1;}
346 if (p > maxp[icharge]) {
348 maxiTrk[icharge] = iTrk;
357 if (nTrkTight != 2) {
362 getObjectPtr<TH1F>(
"cutflow")->Fill(cutIndexPassed);
363 B2DEBUG(22,
"Cutflow: Two tight tracks: index = " << cutIndexPassed);
366 if (nTrkLoose != 2) {
371 getObjectPtr<TH1F>(
"cutflow")->Fill(cutIndexPassed) ;
372 B2DEBUG(22,
"Cutflow: No additional loose tracks: index = " << cutIndexPassed) ;
377 bool oppositelyChargedTracksPassed = maxiTrk[0] != -1 && maxiTrk[1] != -1;
378 if (!oppositelyChargedTracksPassed) {
383 getObjectPtr<TH1F>(
"cutflow")->Fill(cutIndexPassed);
384 B2DEBUG(22,
"Cutflow: Oppositely charged tracks: index = " << cutIndexPassed);
392 double trkEClustLab[2] = {0., 0.};
393 double trkEClustCOM[2] = {0., 0.};
396 TLorentzVector trkp4Lab[2];
397 TLorentzVector trkp4COM[2];
400 int numClustersPerTrack[2] = { 0, 0 };
403 vector<double> goodClustTimes ;
404 vector<double> goodClust_dt99 ;
405 vector<double> goodClustE ;
406 vector<int> goodClustMaxEcrys_cid ;
408 for (
int icharge = 0; icharge < 2; icharge++) {
409 if (maxiTrk[icharge] > -1) {
410 B2DEBUG(22,
"looping over the 2 max pt tracks");
413 if (not tempTrackFit) {continue ;}
415 trkp4Lab[icharge] = tempTrackFit->get4Momentum();
416 trkp4COM[icharge] = boostrotate.
rotateLabToCms() * trkp4Lab[icharge];
417 trkpLab[icharge] = trkp4Lab[icharge].Rho();
418 trkpCOM[icharge] = trkp4COM[icharge].Rho();
424 auto eclClusterRelationsFromTracks = tracks[maxiTrk[icharge]]->getRelationsTo<
ECLCluster>();
425 for (
unsigned int clusterIdx = 0; clusterIdx < eclClusterRelationsFromTracks.size(); clusterIdx++) {
427 B2DEBUG(22,
"Looking at clusters. index = " << clusterIdx);
428 auto cluster = eclClusterRelationsFromTracks[clusterIdx];
431 numClustersPerTrack[icharge]++;
433 double electronTime = cluster->getTime();
434 bool badElectronTime = cluster->hasFailedFitTime();
435 bool badElectronTimeResolution = cluster->hasFailedTimeResolution();
436 if ((fabs(electronTime) < m_timeAbsMax) &&
437 (!badElectronTime) &&
438 (!badElectronTimeResolution)) {
439 trkEClustLab[icharge] += eClust ;
440 short cid = cluster->getMaxECellId() ;
441 goodClustMaxEcrys_cid.push_back(cid) ;
442 goodClustTimes.push_back(electronTime) ;
443 goodClust_dt99.push_back(cluster->getDeltaTime99()) ;
444 goodClustE.push_back(eClust);
448 trkEClustCOM[icharge] = trkEClustLab[icharge] * trkpCOM[icharge] / trkpLab[icharge];
453 E_DIV_p[icharge] = trkEClustCOM[icharge] / trkpCOM[icharge];
460 B2DEBUG(26,
"Extracted time information and E/p for the tracks") ;
471 long unsigned int numGoodElectronClusters_cut = 2 ;
472 if (goodClustTimes.size() != numGoodElectronClusters_cut) {
477 getObjectPtr<TH1F>(
"cutflow")->Fill(cutIndexPassed) ;
478 B2DEBUG(22,
"Cutflow: Exactly " << numGoodElectronClusters_cut
479 <<
" good clusters connected to tracks: index = " << cutIndexPassed);
485 bool E_DIV_p_instance_passed[2] = {
false,
false};
486 double E_DIV_p_CUT = 0.7;
487 for (
int icharge = 0; icharge < 2; icharge++) {
488 E_DIV_p_instance_passed[icharge] = E_DIV_p[icharge] > E_DIV_p_CUT;
490 if (!E_DIV_p_instance_passed[0] || !E_DIV_p_instance_passed[1]) {
495 getObjectPtr<TH1F>(
"cutflow")->Fill(cutIndexPassed);
496 B2DEBUG(22,
"Cutflow: E_i/p_i > " << E_DIV_p_CUT <<
": index = " << cutIndexPassed);
501 double invMassTrk = (trkp4Lab[0] + trkp4Lab[1]).M();
502 double invMass_CUT = 0.9;
504 bool invMassCutsPassed = invMassTrk > (invMass_CUT * boostrotate.
getCMSEnergy());
505 if (!invMassCutsPassed) {
510 getObjectPtr<TH1F>(
"cutflow")->Fill(cutIndexPassed);
511 B2DEBUG(22,
"Cutflow: m(track 1+2) > " << invMass_CUT <<
"*E_COM = " << invMass_CUT <<
" * " << boostrotate.
getCMSEnergy() <<
512 " : index = " << cutIndexPassed);
515 B2DEBUG(22,
"Event passed all cuts");
519 getObjectPtr<TH1F>(
"eventT0")->Fill(evt_t0) ;
521 bool isCDCt0 = (
static_cast<EventT0::EventT0Component>(*m_eventT0->getEventT0Component())).detectorSet.contains(Const::CDC);
522 bool isECLt0 = (
static_cast<EventT0::EventT0Component>(*m_eventT0->getEventT0Component())).detectorSet.contains(Const::ECL);
523 string t0Detector =
"UNKNOWN... WHY?";
526 }
else if (isECLt0) {
530 B2DEBUG(26,
"t0 = " << evt_t0 <<
" ns. t0 is from CDC?=" << isCDCt0 <<
", t0 is from ECL?=" << isECLt0 <<
" t0 from " <<
537 for (
long unsigned int i = 0 ; i < goodClustTimes.size() ; i++) {
538 getObjectPtr<TH1F>(
"clusterTime")->Fill(goodClustTimes[i]) ;
539 getObjectPtr<TH2F>(
"clusterTime_cid")->Fill(goodClustMaxEcrys_cid[i] + 0.001, goodClustTimes[i] , 1) ;
540 getObjectPtr<TH2F>(
"clusterTime_run")->Fill(evtMetaData->getRun() + 0.001, goodClustTimes[i] , 1) ;
541 getObjectPtr<TH2F>(
"clusterTimeClusterE")->Fill(goodClustE[i], goodClustTimes[i], 1) ;
542 getObjectPtr<TH2F>(
"dt99_clusterE")->Fill(goodClustE[i], goodClust_dt99[i], 1) ;
548 m_tree_time = goodClustTimes[i] ;
550 m_tree_t0_unc = evt_t0_unc ;
551 m_E_electron_clust = goodClustE[i] ;
552 m_NtightTracks = nTrkTight ;
553 m_tree_evt_num = evtMetaData->getEvent() ;
554 m_tree_run = evtMetaData->getRun() ;
555 m_tree_cid = goodClustMaxEcrys_cid[i] ;
556 m_tree_dt99 = goodClust_dt99[i] ;
558 m_dbg_tree_electronClusters->Fill() ;
562 B2DEBUG(26,
"Filled cluster tree") ;
566 if (goodClustE[0] > goodClustE[1]) {
567 tDiff = goodClustTimes[0] - goodClustTimes[1];
569 tDiff = goodClustTimes[1] - goodClustTimes[0];
572 getObjectPtr<TH1F>(
"clusterTimeE0E1diff")->Fill(tDiff) ;
578 m_tree_t0_unc = evt_t0_unc ;
579 m_tree_evt_num = evtMetaData->getEvent() ;
580 m_tree_run = evtMetaData->getRun() ;
581 m_NtightTracks = nTrkTight ;
582 m_tree_E0 = goodClustE[0] ;
583 m_tree_E1 = goodClustE[1] ;
584 m_tree_time_fromE0 = goodClustTimes[0] ;
585 m_tree_time_fromE1 = goodClustTimes[1] ;
587 m_dbg_tree_event->Fill() ;
590 int runNum = evtMetaData->getRun();
591 if (m_tree_PreviousRun != runNum) {
592 for (
int icrate = 1; icrate <= 52; icrate++) {
593 m_tree_run = runNum ;
594 m_tree_crateid = icrate ;
595 m_tree_tcrate = Crate_time_ns[icrate] ;
596 m_tree_tcrate_unc = Crate_time_unc_ns[icrate] ;
598 m_dbg_tree_run->Fill() ;
600 m_tree_PreviousRun = m_tree_run ;
604 B2DEBUG(26,
"Filled event tree") ;