diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx index 26708fd2b58f26226ce1384732ed2461534ddd51..aa3e0beaa217e5e1c723028679d1b66519ce30e5 100644 --- a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx +++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx @@ -282,6 +282,191 @@ std::unique_ptr<ZDCDataAnalyzer> ZdcAnalysisTool::initializeLHCf2022() } +std::unique_ptr<ZDCDataAnalyzer> ZdcAnalysisTool::initializepp2023() +{ + + m_deltaTSample = 3.125; + m_numSample = 24; + + ZDCDataAnalyzer::ZDCModuleIntArray peak2ndDerivMinSamples = {{{0, 9, 9, 9}, {0, 9, 10, 8}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinThresholdsHG, peak2ndDerivMinThresholdsLG; + ZDCDataAnalyzer::ZDCModuleFloatArray deltaT0CutLow, deltaT0CutHigh, chisqDivAmpCut; + ZDCDataAnalyzer::ZDCModuleBoolArray fixTau1Arr, fixTau2Arr; + + ZDCDataAnalyzer::ZDCModuleFloatArray tau1 = {{{0, 1.1, 1.1, 1.1}, + {0, 1.1, 1.1, 1.1}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray tau2 = {{{6, 5, 5, 5}, {5.5, 5.5, 5.5, 5.5}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray t0HG = {{{0, 26.3, 26.5, 26.8}, {32, 32, 32, 32}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray t0LG = {{{0, 31.1, 28.1, 27.0}, {0, 26.6, 26.3, 25.3}}}; + + for (size_t side : {0, 1}) { + for (size_t module : {0, 1, 2, 3}) { + fixTau1Arr[side][module] = true; + fixTau2Arr[side][module] = false; + + peak2ndDerivMinThresholdsHG[side][module] = -35; + peak2ndDerivMinThresholdsLG[side][module] = -16; + + deltaT0CutLow[side][module] = -10; + deltaT0CutHigh[side][module] = 10; + chisqDivAmpCut[side][module] = 20; + } + } + + ATH_MSG_DEBUG( "LHCF2022: delta t cut, value low = " << deltaT0CutLow[0][0] << ", high = " << deltaT0CutHigh[0][0] ); + + ZDCDataAnalyzer::ZDCModuleFloatArray HGOverFlowADC = {{{{4000, 4000, 4000, 4000}}, {{4000, 4000, 4000, 4000}}}}; + ZDCDataAnalyzer::ZDCModuleFloatArray HGUnderFlowADC = {{{{1, 1, 1, 1}}, {{1, 1, 1, 1}}}}; + ZDCDataAnalyzer::ZDCModuleFloatArray LGOverFlowADC = {{{{4000, 4000, 4000, 4000}}, {{4000, 4000, 4000, 4000}}}}; + + // For the LHCf run, use low gain samples + // + m_lowGainOnly = true; + + // Construct the data analyzer + // + std::unique_ptr<ZDCDataAnalyzer> zdcDataAnalyzer (new ZDCDataAnalyzer(MakeMessageFunction(), + m_numSample, m_deltaTSample, + m_presample, "FermiExpLHCf", + peak2ndDerivMinSamples, + peak2ndDerivMinThresholdsHG, + peak2ndDerivMinThresholdsLG, + m_lowGainOnly)); + + zdcDataAnalyzer->SetPeak2ndDerivMinTolerances(2); + zdcDataAnalyzer->SetADCOverUnderflowValues(HGOverFlowADC, HGUnderFlowADC, LGOverFlowADC); + zdcDataAnalyzer->SetTauT0Values(fixTau1Arr, fixTau2Arr, tau1, tau2, t0HG, t0LG); + zdcDataAnalyzer->SetCutValues(chisqDivAmpCut, chisqDivAmpCut, deltaT0CutLow, deltaT0CutHigh, deltaT0CutLow, deltaT0CutHigh); + + zdcDataAnalyzer->SetGainFactorsHGLG(0.1, 1); // a gain adjustment of unity applied to LG ADC, 0.1 to HG ADC values + + ZDCDataAnalyzer::ZDCModuleFloatArray noiseSigmasLG = {{{2, 2, 2, 2}, {2, 2, 2, 2}}}; + ZDCDataAnalyzer::ZDCModuleFloatArray noiseSigmasHG = {{{20, 20, 20, 20}, {20, 20, 20, 20}}}; + + zdcDataAnalyzer->SetNoiseSigmas(noiseSigmasHG, noiseSigmasLG); + + // Enable two-pass analysis + // + ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinRepassHG = {{{-12, -12, -12, -12}, + {-12, -12, -12, -12}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinRepassLG = {{{-8, -8, -8, -8}, + {-8, -8, -8, -8}}}; + + zdcDataAnalyzer->enableRepass(peak2ndDerivMinRepassHG, peak2ndDerivMinRepassLG); + + // Set the amplitude fit range limits + // + zdcDataAnalyzer->SetFitMinMaxAmpValues(5, 2, 5000, 5000); + + + m_rpdDataAnalyzer.push_back(std::make_unique<RPDDataAnalyzer>(MakeMessageFunction(), "rpdA", 4, 4, m_numSample)); + m_rpdDataAnalyzer.push_back(std::make_unique<RPDDataAnalyzer>(MakeMessageFunction(), "rpdC", 4, 4, m_numSample)); + + return zdcDataAnalyzer; + +} + +std::unique_ptr<ZDCDataAnalyzer> ZdcAnalysisTool::initializePbPb2023() +{ + // Key configuration parameters needed for the data analyzer construction + // + m_deltaTSample = 3.125; + m_numSample = 24; + + ZDCDataAnalyzer::ZDCModuleIntArray peak2ndDerivMinSamples; + ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinThresholdsHG, peak2ndDerivMinThresholdsLG; + + for (size_t side : {0, 1}) { + for (size_t module : {0, 1, 2, 3}) { + peak2ndDerivMinSamples[side][module] = 9; + peak2ndDerivMinThresholdsHG[side][module] = -35; + peak2ndDerivMinThresholdsLG[side][module] = -16; + } + } + + m_lowGainOnly = false; + + // Construct the data analyzer + // + std::unique_ptr<ZDCDataAnalyzer> zdcDataAnalyzer (new ZDCDataAnalyzer(MakeMessageFunction(), + m_numSample, m_deltaTSample, + m_presample, "FermiExpRun3", + peak2ndDerivMinSamples, + peak2ndDerivMinThresholdsHG, + peak2ndDerivMinThresholdsLG, + m_lowGainOnly)); + zdcDataAnalyzer->set2ndDerivStep(2); + zdcDataAnalyzer->SetPeak2ndDerivMinTolerances(2); + zdcDataAnalyzer->SetGainFactorsHGLG(10, 1); // a gain adjustment of 10 applied to LG ADC, 1 to HG ADC values + + // Now set cuts and default fit parameters + // + ZDCDataAnalyzer::ZDCModuleFloatArray deltaT0CutLow, deltaT0CutHigh, chisqDivAmpCut; + ZDCDataAnalyzer::ZDCModuleBoolArray fixTau1Arr, fixTau2Arr; + + ZDCDataAnalyzer::ZDCModuleFloatArray tau1 = {{{1.1, 1.1, 1.1, 1.1}, + {1.1, 1.1, 1.1, 1.1}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray tau2 = {{{5.5, 5.5, 5.5, 5.5}, {5.5, 5.5, 5.5, 5.5}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray t0HG = {{{27, 27, 27, 27}, {27, 27, 27, 27}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray t0LG = {{{29, 29, 29, 29}, {29, 29, 29, 29}}}; + + for (size_t side : {0, 1}) { + for (size_t module : {0, 1, 2, 3}) { + fixTau1Arr[side][module] = true; + fixTau2Arr[side][module] = true; + + peak2ndDerivMinThresholdsHG[side][module] = -30; + peak2ndDerivMinThresholdsLG[side][module] = -12; + + deltaT0CutLow[side][module] = -m_deltaTCut; + deltaT0CutHigh[side][module] = m_deltaTCut; + chisqDivAmpCut[side][module] = m_ChisqRatioCut; + } + } + + ATH_MSG_DEBUG( "PbPb2023: delta t cut, value low = " << deltaT0CutLow[0][0] << ", high = " << deltaT0CutHigh[0][0] ); + + ZDCDataAnalyzer::ZDCModuleFloatArray HGOverFlowADC = {{{{3500, 3500, 3500, 3500}}, {{3500, 3500, 3500, 3500}}}}; + ZDCDataAnalyzer::ZDCModuleFloatArray HGUnderFlowADC = {{{{1, 1, 1, 1}}, {{1, 1, 1, 1}}}}; + ZDCDataAnalyzer::ZDCModuleFloatArray LGOverFlowADC = {{{{4000, 4000, 4000, 4000}}, {{4000, 4000, 4000, 4000}}}}; + + zdcDataAnalyzer->SetADCOverUnderflowValues(HGOverFlowADC, HGUnderFlowADC, LGOverFlowADC); + zdcDataAnalyzer->SetTauT0Values(fixTau1Arr, fixTau2Arr, tau1, tau2, t0HG, t0LG); + zdcDataAnalyzer->SetCutValues(chisqDivAmpCut, chisqDivAmpCut, deltaT0CutLow, deltaT0CutHigh, deltaT0CutLow, deltaT0CutHigh); + + ZDCDataAnalyzer::ZDCModuleFloatArray noiseSigmasLG = {{{2, 2, 2, 2}, {2, 2, 2, 2}}}; + ZDCDataAnalyzer::ZDCModuleFloatArray noiseSigmasHG = {{{1, 1, 1, 1}, {1, 1, 1, 1}}}; + + zdcDataAnalyzer->SetNoiseSigmas(noiseSigmasHG, noiseSigmasLG); + + // Enable two-pass analysis + // + ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinRepassHG = {{{-12, -12, -12, -12}, + {-12, -12, -12, -12}}}; + + ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinRepassLG = {{{-8, -8, -8, -8}, + {-8, -8, -8, -8}}}; + + zdcDataAnalyzer->enableRepass(peak2ndDerivMinRepassHG, peak2ndDerivMinRepassLG); + + // Set the amplitude fit range limits + // + zdcDataAnalyzer->SetFitMinMaxAmpValues(5, 2, 5000, 5000); + + m_rpdDataAnalyzer.push_back(std::make_unique<RPDDataAnalyzer>(MakeMessageFunction(), "rpdA", 4, 4, m_numSample)); + m_rpdDataAnalyzer.push_back(std::make_unique<RPDDataAnalyzer>(MakeMessageFunction(), "rpdC", 4, 4, m_numSample)); + + return zdcDataAnalyzer; +} + std::unique_ptr<ZDCDataAnalyzer> ZdcAnalysisTool::initializeDefault() { // We rely completely on the default parameters specified in the job properties to control: @@ -905,6 +1090,12 @@ StatusCode ZdcAnalysisTool::initialize() else if (m_configuration == "LHCf2022") { m_zdcDataAnalyzer = initializeLHCf2022(); } + else if (m_configuration == "pp2023") { + m_zdcDataAnalyzer = initializepp2023(); + } + else if (m_configuration == "PbPb2023") { + m_zdcDataAnalyzer = initializePbPb2023(); + } else { ATH_MSG_ERROR("Unknown configuration: " << m_configuration); return StatusCode::FAILURE; @@ -1024,7 +1215,11 @@ StatusCode ZdcAnalysisTool::recoZdcModules(const xAOD::ZdcModuleContainer& modul ATH_MSG_DEBUG("Processing modules"); for (const auto zdcModule : moduleContainer) { + + ATH_MSG_DEBUG(ZdcModuleToString(*zdcModule)); + if (zdcModule->zdcType() == 1) { + ATH_MSG_DEBUG("RPD side " << zdcModule->zdcSide() << " chan " << zdcModule->zdcChannel() ); if (m_LHCRun < 3) continue; // type == 1 -> pixel data in runs 2 and 3, skip // This is RPD data in Run 3 @@ -1128,15 +1323,15 @@ StatusCode ZdcAnalysisTool::recoZdcModules(const xAOD::ZdcModuleContainer& modul int side = (zdcModule->zdcSide() == -1) ? 0 : 1 ; int mod = zdcModule->zdcModule(); - if (zdcModule->zdcType() == 1) { + if (zdcModule->zdcType() == 1 && m_LHCRun==3) { // this is the RPD if (m_writeAux) { int rpdChannel = zdcModule->zdcChannel() - 1; zdcModule->auxdecor<float>("RPDChannelAmplitude" + m_auxSuffix) = m_rpdDataAnalyzer.at(side)->getChSumADC(rpdChannel); zdcModule->auxdecor<unsigned int>("RPDChannelMaxSample" + m_auxSuffix) = m_rpdDataAnalyzer.at(side)->getChMaxADCSample(rpdChannel); } - } else { - // this is the ZDC + } else if (zdcModule->zdcType() == 0) { + // this is the main ZDC if (m_writeAux) { if (m_doCalib) { float calibEnergy = m_zdcDataAnalyzer->GetModuleCalibAmplitude(side, mod); diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h index 1dbabcbe38f5a504b83573749ff3ca59fbc460ad..6be9acf1859f8f440919c6c9e2c18cd686a11f54 100644 --- a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h +++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h @@ -9,6 +9,7 @@ #include "AsgDataHandles/ReadHandleKey.h" #include "xAODForward/ZdcModuleContainer.h" +#include "xAODForward/ZdcModuleToString.h" #include "xAODTrigL1Calo/TriggerTowerContainer.h" @@ -113,6 +114,8 @@ private: std::unique_ptr<ZDCDataAnalyzer> initializepPb2016(); std::unique_ptr<ZDCDataAnalyzer> initializePbPb2018(); std::unique_ptr<ZDCDataAnalyzer> initializeLHCf2022(); + std::unique_ptr<ZDCDataAnalyzer> initializepp2023(); + std::unique_ptr<ZDCDataAnalyzer> initializePbPb2023(); StatusCode configureNewRun(unsigned int runNumber); diff --git a/ForwardDetectors/ZDC/ZdcNtuple/Root/ZdcNtuple.cxx b/ForwardDetectors/ZDC/ZdcNtuple/Root/ZdcNtuple.cxx index 409e4137318d47ef25196200d0ed0570d58c6648..b27569d4a8b743dd503f383238aff1352c85d3d9 100644 --- a/ForwardDetectors/ZDC/ZdcNtuple/Root/ZdcNtuple.cxx +++ b/ForwardDetectors/ZDC/ZdcNtuple/Root/ZdcNtuple.cxx @@ -49,6 +49,8 @@ ZdcNtuple :: ZdcNtuple (const std::string& name, ISvcLocator *pSvcLocator) declareProperty("auxSuffix", auxSuffix = "", "comment"); declareProperty("nsamplesZdc", nsamplesZdc = 24, "comment"); declareProperty("lhcf2022", lhcf2022 = true,"comment"); + declareProperty("lhcf2022afp", lhcf2022afp = true,"comment"); + declareProperty("lhcf2022zdc", lhcf2022zdc = true,"comment"); declareProperty("zdcConfig", zdcConfig = "PbPb2018", "argument to configure ZdcAnalysisTool"); declareProperty("doZdcCalib", doZdcCalib = false, "perform ZDC energy calibration"); @@ -239,7 +241,56 @@ StatusCode ZdcNtuple :: initialize () m_outputTree->Branch("edgeGapA", &t_edgeGapA, "edgeGapA/F"); m_outputTree->Branch("edgeGapC", &t_edgeGapC, "edgeGapC/F"); + m_outputTree->Branch("ntrk", &t_ntrk, "ntrk/i"); + if (enableTracks) + { + m_outputTree->Branch("trk_pt", "vector<float>", &t_trk_pt); + m_outputTree->Branch("trk_eta", "vector<float>", &t_trk_eta); + m_outputTree->Branch("trk_theta", "vector<float>", &t_trk_theta); + m_outputTree->Branch("trk_phi", "vector<float>", &t_trk_phi); + m_outputTree->Branch("trk_e", "vector<float>", &t_trk_e); + m_outputTree->Branch("trk_index", "vector<int>", &t_trk_index); + m_outputTree->Branch("trk_d0", "vector<float>", &t_trk_d0); + m_outputTree->Branch("trk_z0", "vector<float>", &t_trk_z0); + m_outputTree->Branch("trk_vz", "vector<float>", &t_trk_vz); + m_outputTree->Branch("trk_vtxz", "vector<float>", &t_trk_vtxz); + m_outputTree->Branch("trk_pixeldEdx", "vector<float>", &t_trk_pixeldEdx); + m_outputTree->Branch("trk_charge", "vector<int8_t>", &t_trk_charge); + m_outputTree->Branch("trk_quality", "vector<int16_t>", &t_trk_quality); + m_outputTree->Branch("trk_nPixHits", "vector<uint8_t>", &t_trk_nPixHits); + m_outputTree->Branch("trk_nSctHits", "vector<uint8_t>", &t_trk_nSctHits); + m_outputTree->Branch("trk_nPixDead", "vector<uint8_t>", &t_trk_nPixDead); + m_outputTree->Branch("trk_nSctDead", "vector<uint8_t>", &t_trk_nSctDead); + m_outputTree->Branch("trk_nPixHoles", "vector<uint8_t>", &t_trk_nPixHoles); + m_outputTree->Branch("trk_nSctHoles", "vector<uint8_t>", &t_trk_nSctHoles); + m_outputTree->Branch("trk_nTrtHits", "vector<uint8_t>", &t_trk_nTrtHits); + m_outputTree->Branch("trk_nTrtOutliers", "vector<uint8_t>", &t_trk_nTrtOutliers); + m_outputTree->Branch("trk_inPixHits", "vector<uint8_t>", &t_trk_inPixHits); + m_outputTree->Branch("trk_exPixHits", "vector<uint8_t>", &t_trk_exPixHits); + m_outputTree->Branch("trk_ninPixHits", "vector<uint8_t>", &t_trk_ninPixHits); + m_outputTree->Branch("trk_nexPixHits", "vector<uint8_t>", &t_trk_nexPixHits); + } + + } + + if( lhcf2022 || lhcf2022zdc || lhcf2022afp){ + m_outputTree->Branch("nProtons", &nProtons, "nProtons/i"); + m_outputTree->Branch("proton_pt", "vector<double>", &proton_pt); + m_outputTree->Branch("proton_eta", "vector<double>", &proton_eta); + m_outputTree->Branch("proton_phi", "vector<double>", &proton_phi); + m_outputTree->Branch("proton_e", "vector<double>", &proton_e); + m_outputTree->Branch("proton_side", "vector<int>", &proton_side); + m_outputTree->Branch("proton_eLoss", "vector<double>", &proton_eLoss); + m_outputTree->Branch("proton_t", "vector<double>", &proton_t); + m_outputTree->Branch("proton_track_stationID", "vector<vector<int>>", &proton_track_stationID); + m_outputTree->Branch("proton_track_nClusters", "vector<vector<int>>", &proton_track_nClusters); + m_outputTree->Branch("proton_track_xLocal", "vector<vector<float>>", &proton_track_xLocal); + m_outputTree->Branch("proton_track_yLocal", "vector<vector<float>>", &proton_track_yLocal); + m_outputTree->Branch("proton_track_zLocal", "vector<vector<float>>", &proton_track_zLocal); + m_outputTree->Branch("proton_track_xSlope", "vector<vector<float>>", &proton_track_xSlope); + m_outputTree->Branch("proton_track_ySlope", "vector<vector<float>>", &proton_track_ySlope); } + } ANA_MSG_DEBUG("Anti-howdy from Initialize!"); @@ -435,6 +486,10 @@ StatusCode ZdcNtuple :: execute () ANA_CHECK(evtStore()->retrieve( m_truthParticleContainer, "TruthParticles")); } + if( lhcf2022||lhcf2022zdc||lhcf2022afp ){ + ANA_CHECK(evtStore()->retrieve( m_afpProtons, "AFPProtonContainer")); + processProtons(); + } // if trigger enabled, only write out events which pass one of them, unless using MC if (enableTrigger && !passTrigger && !m_isMC && writeOnlyTriggers) return StatusCode::SUCCESS; @@ -478,6 +533,22 @@ void ZdcNtuple::processZdcNtupleFromModules() t_ZdcModuleFitT0[iside][imod] = 0; t_ZdcModuleBkgdMaxFraction[iside][imod] = 0; t_ZdcModuleAmpError[iside][imod] = 0; t_ZdcModuleMinDeriv2nd[iside][imod] = 0; t_ZdcModulePresample[iside][imod] = 0; t_ZdcModulePreSampleAmp[iside][imod] = 0; t_ZdcLucrodTriggerAmp[iside][imod] = 0;t_ZdcModuleMaxADC[iside][imod] = 0; + + if (enableOutputSamples) + { + for (int ig=0;ig<2;ig++) + { + for (int id=0;id<2;id++) + { + for (int isamp=0;isamp<((int)nsamplesZdc);isamp++) + { + if (nsamplesZdc==7) t_raw7[iside][imod][ig][id][isamp]=0; + if (nsamplesZdc==15) t_raw15[iside][imod][ig][id][isamp]=0; + if (nsamplesZdc==24) t_raw24[iside][imod][ig][id][isamp]=0; + } + } + } + } } } @@ -499,7 +570,8 @@ void ZdcNtuple::processZdcNtupleFromModules() t_ZdcAmp[iside] = zdcSum->auxdataConst<float>("UncalibSum"+auxSuffix); t_ZdcAmpErr[iside] = zdcSum->auxdataConst<float>("UncalibSumErr"+auxSuffix); - t_ZdcLucrodTriggerSideAmp[iside] = zdcSum->auxdataConst<uint16_t>("LucrodTriggerSideAmp"); + if (zdcSum->isAvailable<uint16_t>("LucrodTriggerSideAmp")) + t_ZdcLucrodTriggerSideAmp[iside] = zdcSum->auxdataConst<uint16_t>("LucrodTriggerSideAmp"); ANA_MSG_VERBOSE("processZdcNtupleFromModules: ZdcSum energy = " << t_ZdcEnergy[iside]); t_ZdcTime[iside] = zdcSum->auxdataConst<float>("AverageTime"+auxSuffix); @@ -537,8 +609,10 @@ void ZdcNtuple::processZdcNtupleFromModules() t_ZdcModuleMinDeriv2nd[iside][imod] = zdcMod->auxdataConst<float>("MinDeriv2nd" + auxSuffix); t_ZdcModulePresample[iside][imod] = zdcMod->auxdataConst<float>("Presample" + auxSuffix); t_ZdcModulePreSampleAmp[iside][imod] = zdcMod->auxdataConst<float>("PreSampleAmp" + auxSuffix); - t_ZdcLucrodTriggerAmp[iside][imod] = zdcMod->auxdataConst<uint16_t>("LucrodTriggerAmp"); - t_ZdcModuleMaxADC[iside][imod] = zdcMod->auxdataConst<float>("MaxADC"); + if (zdcMod->isAvailable<uint16_t>("LucrodTriggerAmp")) + t_ZdcLucrodTriggerAmp[iside][imod] = zdcMod->auxdataConst<uint16_t>("LucrodTriggerAmp"); + if (zdcMod->isAvailable<float>("MaxADC")) + t_ZdcModuleMaxADC[iside][imod] = zdcMod->auxdataConst<float>("MaxADC"); if (enableOutputSamples) { @@ -1328,6 +1402,66 @@ void ZdcNtuple::processClusters() return; } +void ZdcNtuple::processProtons(){ + + proton_pt.clear(); + proton_eta.clear(); + proton_phi.clear(); + proton_e.clear(); + proton_side.clear(); + proton_eLoss.clear(); + proton_t.clear(); + + proton_track_stationID.clear(); + proton_track_nClusters.clear(); + proton_track_xLocal.clear(); + proton_track_yLocal.clear(); + proton_track_zLocal.clear(); + proton_track_xSlope.clear(); + proton_track_ySlope.clear(); + + proton_track_stationID.resize(m_afpProtons->size(), std::vector<int>(2)); + proton_track_nClusters.resize(m_afpProtons->size(), std::vector<int>(2)); + proton_track_xLocal.resize(m_afpProtons->size(), std::vector<float>(2)); + proton_track_yLocal.resize(m_afpProtons->size(), std::vector<float>(2)); + proton_track_zLocal.resize(m_afpProtons->size(), std::vector<float>(2)); + proton_track_xSlope.resize(m_afpProtons->size(), std::vector<float>(2)); + proton_track_ySlope.resize(m_afpProtons->size(), std::vector<float>(2)); + + nProtons = 0; + + for(const auto * proton: *m_afpProtons){ + + proton_pt.push_back(proton->pt()); + proton_eta.push_back(proton->eta()); + proton_phi.push_back(proton->phi()); + proton_e.push_back(proton->e()); + proton_side.push_back(proton->side()); + + proton_eLoss.push_back((6800.-proton->e())/6800.); + p_scat.SetPtEtaPhiE(proton->pt(), proton->eta(), proton->phi(), proton->e()); + (signbit(proton->eta())) ? p_beam.SetPxPyPzE(0.0, 0.0, -6800.0, 6800.0) : p_beam.SetPxPyPzE(0.0, 0.0, 6800.0,\ + 6800.0); + + proton_t.push_back( (p_beam - p_scat)*(p_beam - p_scat)); + + for(int i=0; i< int(proton->nTracks()); i++){ + + proton_track_stationID.at(nProtons).at(i) = proton->track(i)->stationID(); + proton_track_nClusters.at(nProtons).at(i) = proton->track(i)->nClusters(); + proton_track_xLocal.at(nProtons).at(i) = proton->track(i)->xLocal(); + proton_track_yLocal.at(nProtons).at(i) = proton->track(i)->yLocal(); + proton_track_zLocal.at(nProtons).at(i) = proton->track(i)->zLocal(); + proton_track_xSlope.at(nProtons).at(i) = proton->track(i)->xSlope(); + proton_track_ySlope.at(nProtons).at(i) = proton->track(i)->ySlope(); + + } + + nProtons++; + } + + return; +} uint32_t ZdcNtuple::acceptEvent() { @@ -1393,22 +1527,30 @@ void ZdcNtuple::setupTriggerHistos() else { if (lhcf2022) - { - triggers.push_back("HLT_noalg_L1LHCF"); - triggers.push_back("HLT_noalg_ZDCPEB_L1ZDC_OR"); - triggers.push_back("HLT_noalg_ZDCPEB_L1LHCF"); - triggers.push_back("HLT_noalg_L1ZDC_OR"); - triggers.push_back("HLT_noalg_L1ZDC_XOR_E2"); - triggers.push_back("HLT_noalg_L1ZDC_XOR_E1_E3"); - triggers.push_back("HLT_noalg_L1ZDC_A_AND_C"); - triggers.push_back("HLT_mb_sptrk_L1ZDC_OR"); - triggers.push_back("HLT_mb_sptrk_L1ZDC_XOR_E2"); - triggers.push_back("HLT_mb_sptrk_L1ZDC_XOR_E1_E3"); - triggers.push_back("HLT_mb_sptrk_L1ZDC_A_AND_C"); - triggers.push_back("HLT_mb_sp100_trk30_hmt_L1ZDC_XOR_E2"); - triggers.push_back("HLT_mb_sp100_trk30_hmt_L1ZDC_XOR_E1_E3"); - triggers.push_back("HLT_mb_sp100_trk30_hmt_L1ZDC_A_AND_C"); - } + { + triggers.push_back("HLT_noalg_L1LHCF"); + } + if (lhcf2022afp) + { + triggers.push_back("HLT_noalg_AFPPEB_L1AFP_A"); + triggers.push_back("HLT_noalg_AFPPEB_L1AFP_C"); + } + if (lhcf2022zdc) + { + triggers.push_back("HLT_noalg_ZDCPEB_L1ZDC_OR"); + triggers.push_back("HLT_noalg_ZDCPEB_L1LHCF"); + triggers.push_back("HLT_noalg_L1ZDC_OR"); + triggers.push_back("HLT_noalg_L1ZDC_XOR_E2"); + triggers.push_back("HLT_noalg_L1ZDC_XOR_E1_E3"); + triggers.push_back("HLT_noalg_L1ZDC_A_AND_C"); + triggers.push_back("HLT_mb_sptrk_L1ZDC_OR"); + triggers.push_back("HLT_mb_sptrk_L1ZDC_XOR_E2"); + triggers.push_back("HLT_mb_sptrk_L1ZDC_XOR_E1_E3"); + triggers.push_back("HLT_mb_sptrk_L1ZDC_A_AND_C"); + triggers.push_back("HLT_mb_sp100_trk30_hmt_L1ZDC_XOR_E2"); + triggers.push_back("HLT_mb_sp100_trk30_hmt_L1ZDC_XOR_E1_E3"); + triggers.push_back("HLT_mb_sp100_trk30_hmt_L1ZDC_A_AND_C"); + } } } diff --git a/ForwardDetectors/ZDC/ZdcNtuple/ZdcNtuple/ZdcNtuple.h b/ForwardDetectors/ZDC/ZdcNtuple/ZdcNtuple/ZdcNtuple.h index cedc7af8e379a7f17b064c5ad78801eb9115ebbd..9774a4e24e799957cb4ca5b80c2697289ded8daf 100644 --- a/ForwardDetectors/ZDC/ZdcNtuple/ZdcNtuple/ZdcNtuple.h +++ b/ForwardDetectors/ZDC/ZdcNtuple/ZdcNtuple/ZdcNtuple.h @@ -37,6 +37,10 @@ #include "ZdcAnalysis/IZdcAnalysisTool.h" +#include "xAODForward/AFPProtonContainer.h" +#include "xAODForward/AFPTrackContainer.h" +#include <TLorentzVector.h> + class ZdcNtuple : public EL::AnaAlgorithm { public: @@ -67,6 +71,8 @@ public: std::string auxSuffix; // what to add to name the new data, when reprocessing size_t nsamplesZdc; // nsamples expected by ZDC tool bool lhcf2022; + bool lhcf2022zdc; + bool lhcf2022afp; bool doZdcCalib; std::string zdcConfig; @@ -98,6 +104,7 @@ public: const xAOD::EnergySumRoI* m_lvl1EnergySumRoI; const xAOD::TruthParticleContainer* m_truthParticleContainer; const xAOD::TriggerTowerContainer* m_TTcontainer; + const xAOD::AFPProtonContainer* m_afpProtons; int m_nTriggers; @@ -296,6 +303,27 @@ public: float t_clusetaMax; float t_clusphiMax; + // AFP protons + + int nProtons; + std::vector<double> proton_pt; + std::vector<double> proton_eta; + std::vector<double> proton_phi; + std::vector<double> proton_e; + std::vector<int> proton_side; + std::vector<double> proton_eLoss; + std::vector<double> proton_t; + std::vector<std::vector<int>> proton_track_stationID; + std::vector<std::vector<float>> proton_track_xLocal; + std::vector<std::vector<float>> proton_track_yLocal; + std::vector<std::vector<float>> proton_track_zLocal; + std::vector<std::vector<float>> proton_track_xSlope; + std::vector<std::vector<float>> proton_track_ySlope; + std::vector<std::vector<int>> proton_track_nClusters; + + TLorentzVector p_beam; + TLorentzVector p_scat; + // end of Histograms // this is a standard constructor @@ -314,6 +342,7 @@ public: void reprocessZdcModule(const xAOD::ZdcModule* zdcMod, bool flipdelay=0); void processTriggerTowers(); void processGaps(); + void processProtons(); double dR(const double eta1, const double phi1, const double eta2, const double phi2);