diff --git a/src/libraries/ANALYSIS/DEventWriterROOT.cc b/src/libraries/ANALYSIS/DEventWriterROOT.cc index deaa335bcb..c31e280f65 100644 --- a/src/libraries/ANALYSIS/DEventWriterROOT.cc +++ b/src/libraries/ANALYSIS/DEventWriterROOT.cc @@ -768,12 +768,19 @@ void DEventWriterROOT::Create_Branches_ChargedHypotheses(DTreeBranchRegister& lo locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "Energy_ECAL"), locArraySizeString, dInitNumNeutralArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "NumBlocks_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "E1E9_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "E9E25_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "SumU_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "SumV_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + //SHOWER MATCHING: locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaPhi"), locArraySizeString, dInitNumTrackArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaZ"), locArraySizeString, dInitNumTrackArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackFCAL_DOCA"), locArraySizeString, dInitNumTrackArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackFCAL_DeltaT"), locArraySizeString, dInitNumTrackArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackECAL_DOCA"), locArraySizeString, dInitNumTrackArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackECAL_DeltaT"), locArraySizeString, dInitNumTrackArraySize); //DIRC: if(DIRC_OUTPUT) { @@ -857,13 +864,20 @@ void DEventWriterROOT::Create_Branches_NeutralHypotheses(DTreeBranchRegister& lo locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "Energy_CCAL"), locArraySizeString, dInitNumNeutralArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "Energy_ECAL"), locArraySizeString, dInitNumNeutralArraySize); if(ECAL_VERBOSE_OUTPUT) { + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "E1E9_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "E9E25_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "SumU_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "SumV_ECAL"), locArraySizeString, dInitNumNeutralArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "NumBlocks_ECAL"), locArraySizeString, dInitNumNeutralArraySize); + } //NEARBY TRACKS locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaPhi"), locArraySizeString, dInitNumNeutralArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaZ"), locArraySizeString, dInitNumNeutralArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackFCAL_DOCA"), locArraySizeString, dInitNumNeutralArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackFCAL_DeltaT"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackECAL_DOCA"), locArraySizeString, dInitNumNeutralArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "TrackECAL_DeltaT"), locArraySizeString, dInitNumNeutralArraySize); //PHOTON PID INFO //Computed using DVertex (best estimate of reaction vertex using all "good" tracks) @@ -1884,8 +1898,8 @@ void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned locECALShower = locECALShowerMatchParams->dECALShower; } - //shared_ptr locECALSingleHitMatchParams - // = locChargedTrackHypothesis->Get_ECALSingleHitMatchParams(); + shared_ptr locECALSingleHitMatchParams + = locChargedTrackHypothesis->Get_ECALSingleHitMatchParams(); //IDENTIFIERS locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackID"), locTrackTimeBased->candidateid, locArrayIndex); @@ -2007,6 +2021,16 @@ void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "NumBlocks_FCAL"), locNumBlocksFCAL, locArrayIndex); double locNumBlocksECAL = (locECALShower != NULL) ? locECALShower->nBlocks : 0.0; + + double locE1E9ECAL = (locECALShower != NULL) ? locECALShower->E1E9 : 0.0; + double locE9E25ECAL = (locECALShower != NULL) ? locECALShower->E9E25 : 0.0; + double locSumUECAL = (locECALShower != NULL) ? locECALShower->sumU : 0.0; + double locSumVECAL = (locECALShower != NULL) ? locECALShower->sumV : 0.0; + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "E1E9_ECAL"), locE1E9ECAL, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "E9E25_ECAL"), locE9E25ECAL, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "SumU_ECAL"), locSumUECAL, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "SumV_ECAL"), locSumVECAL, locArrayIndex); + if (locECALSingleHitMatchParams!=nullptr) locNumBlocksECAL=1.; locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "NumBlocks_ECAL"), locNumBlocksECAL, locArrayIndex); //TIMING INFO @@ -2030,6 +2054,31 @@ void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaPhi"), locTrackBCAL_DeltaPhi, locArrayIndex); locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaZ"), locTrackBCAL_DeltaZ, locArrayIndex); + //SHOWER MATCHING: ECAL + double locDOCAToShower_ECAL = 999.0; + double locDeltaT_TrackToShower_ECAL = 999.0; + if(locECALShowerMatchParams!=nullptr) { + locDOCAToShower_ECAL = locECALShowerMatchParams->dDOCAToShower; + locDeltaT_TrackToShower_ECAL = locECALShower->t + -(locChargedTrackHypothesis->t0() + + locECALShowerMatchParams->dFlightTime); + } + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackECAL_DOCA"), locDOCAToShower_ECAL, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackECAL_DeltaT"), locDeltaT_TrackToShower_ECAL, locArrayIndex); + + // Single hit matching: ECAL + if(locECALSingleHitMatchParams!=nullptr){ + double locDOCAToECALHit = locECALSingleHitMatchParams->dDOCAToHit; + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackECAL_DOCA"), locDOCAToECALHit, locArrayIndex); + + locDeltaT_TrackToShower_ECAL = locECALSingleHitMatchParams->dTHit + - ( locChargedTrackHypothesis->t0() + + locECALSingleHitMatchParams->dFlightTime); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackECAL_DeltaT"), locDeltaT_TrackToShower_ECAL, locArrayIndex); + + } + + //SHOWER MATCHING: FCAL double locDOCAToShower_FCAL = 999.0; double locDeltaT_TrackToShower_FCAL = 999.0; @@ -2217,7 +2266,7 @@ void DEventWriterROOT::Fill_NeutralHypo(DTreeFillData* locTreeFillData, unsigned double locE9E25FCAL = (locFCALShower != NULL) ? locFCALShower->getE9E25() : 0.0; double locSumUFCAL = (locFCALShower != NULL) ? locFCALShower->getSumU() : 0.0; double locSumVFCAL = (locFCALShower != NULL) ? locFCALShower->getSumV() : 0.0; - double locNumBlocksFCAL = (locFCALShower != NULL) ? locFCALShower->getNumBlocks() : 0.0; + double locNumBlocksFCAL = (locFCALShower != NULL) ? locFCALShower->getNumBlocks() : 0.0; locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "E1E9_FCAL"), locE1E9FCAL, locArrayIndex); locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "E9E25_FCAL"), locE9E25FCAL, locArrayIndex); locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "SumU_FCAL"), locSumUFCAL, locArrayIndex); @@ -2227,6 +2276,14 @@ void DEventWriterROOT::Fill_NeutralHypo(DTreeFillData* locTreeFillData, unsigned if(ECAL_VERBOSE_OUTPUT) { double locNumBlocksECAL = (locECALShower != NULL) ? locECALShower->nBlocks : 0.0; + double locE1E9ECAL = (locECALShower != NULL) ? locECALShower->E1E9 : 0.0; + double locE9E25ECAL = (locECALShower != NULL) ? locECALShower->E9E25 : 0.0; + double locSumUECAL = (locECALShower != NULL) ? locECALShower->sumU : 0.0; + double locSumVECAL = (locECALShower != NULL) ? locECALShower->sumV : 0.0; + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "E1E9_ECAL"), locE1E9ECAL, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "E9E25_ECAL"), locE9E25ECAL, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "SumU_ECAL"), locSumUECAL, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "SumV_ECAL"), locSumVECAL, locArrayIndex); locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "NumBlocks_ECAL"), locNumBlocksECAL, locArrayIndex); } @@ -2248,6 +2305,20 @@ void DEventWriterROOT::Fill_NeutralHypo(DTreeFillData* locTreeFillData, unsigned locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaPhi"), locNearestTrackBCALDeltaPhi, locArrayIndex); locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaZ"), locNearestTrackBCALDeltaZ, locArrayIndex); + //Track DOCA to Shower - FCAL + double locDistanceToNearestTrack_ECAL = 999.0; + double locDeltaT_TrackToShower_ECAL = 999.0; + if(locECALShower != NULL) + { + if(!locDetectorMatches->Get_DistanceToNearestTrack(locECALShower, locDistanceToNearestTrack_ECAL)) + locDistanceToNearestTrack_ECAL = 999.0; + if(locDistanceToNearestTrack_ECAL > 999.0) + locDistanceToNearestTrack_ECAL = 999.0; + locDeltaT_TrackToShower_ECAL = locECALShower->t - ( locNeutralParticleHypothesis->t0() + locECALShower->timeTrack ); + } + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackECAL_DOCA"), locDistanceToNearestTrack_ECAL, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "TrackECAL_DeltaT"), locDeltaT_TrackToShower_ECAL, locArrayIndex); + //Track DOCA to Shower - FCAL double locDistanceToNearestTrack_FCAL = 999.0; double locDeltaT_TrackToShower_FCAL = 999.0; diff --git a/src/libraries/ECAL/DECALShower.h b/src/libraries/ECAL/DECALShower.h index a834098c3a..43e0fefcb9 100644 --- a/src/libraries/ECAL/DECALShower.h +++ b/src/libraries/ECAL/DECALShower.h @@ -24,6 +24,7 @@ struct DECALShower : public JObject { TMatrixFSym ExyztCovariance; bool isNearBorder; int nBlocks; + double sumU,sumV,docaTrack,timeTrack; float EErr() const { return sqrt(ExyztCovariance(0,0)); } float xErr() const { return sqrt(ExyztCovariance(1,1)); } @@ -79,6 +80,10 @@ struct DECALShower : public JObject { summary.add(ZTcorr(), "ZTcorr", "%5.3f"); summary.add(nBlocks,"Number of blocks","%d"); summary.add(isNearBorder,"Near border?","%d"); + summary.add(sumU,"U moment","%5.3f"); + summary.add(sumV,"V moment","%5.3f"); + summary.add(docaTrack,"Distance to nearest track","%5.3f"); + summary.add(timeTrack,"Time for nearest track","%5.3f"); } }; diff --git a/src/libraries/ECAL/DECALShower_factory.cc b/src/libraries/ECAL/DECALShower_factory.cc index e12366ba64..1d16598566 100644 --- a/src/libraries/ECAL/DECALShower_factory.cc +++ b/src/libraries/ECAL/DECALShower_factory.cc @@ -15,6 +15,7 @@ #include #include #include +#include #include //------------------ @@ -53,6 +54,10 @@ void DECALShower_factory::BeginRun(const std::shared_ptr& event) auto geom = geo_manager->GetDGeometry(runnumber); geom->GetECALZ(mECALz); mECALzBack=mECALz+20.; + if (geom->HaveInsert()){ + event->GetSingle(dECALGeom); + } + geom->GetTargetZ(mVertexZ); // Get calibration constant auto jcalib = app->GetService()->GetJCalibration(runnumber); @@ -76,6 +81,7 @@ void DECALShower_factory::Process(const std::shared_ptr& event) { vectorclusters; event->Get(clusters); + auto wbtracks=event->Get(); for (size_t i=0;i& event) cov(3,0)=cov(0,3)=X0_over_E*cov(0,0); shower->ExyztCovariance.ResizeTo(5,5); shower->ExyztCovariance=cov; - + + // Find position of closest track to the shower + double min_distance=1e6,timeTrack=1e6; + DVector3 proj_pos_for_shower_shape(0,0,shower->pos.z()); + for (const auto& wbtrack:wbtracks){ + DVector3 proj_pos,proj_mom; + double flight_time=0.; + if (!wbtrack->GetProjection(SYS_ECAL,proj_pos,&proj_mom,&flight_time)) continue; + proj_pos+=((shower->pos.z()-proj_pos.z())/proj_mom.z())*proj_mom; + DVector3 diff=shower->pos-proj_pos; + double d=diff.Perp(); + if (dposition().z()-mVertexZ)/SPEED_OF_LIGHT+flight_time; + } + } + // Compute a couple of shower shaper parameters + GetUV(cluster,shower->pos,proj_pos_for_shower_shape,shower->sumU, + shower->sumV); + shower->docaTrack=min_distance; + shower->timeTrack=timeTrack; shower->AddAssociatedObject(cluster); Insert(shower); @@ -148,3 +175,29 @@ double DECALShower_factory::GetCorrectedZ(double E) const { return mECALz+dZmax; } + +// Compute the energy-weighted second moment of the shower along and +// perpendicular to an axis formed by pointing from the shower to nearest +// track. This code mimics similar code in DFCALShower_factory. +void DECALShower_factory::GetUV(const DECALCluster *cluster, + const DVector3 &showerPos, + const DVector3 &trackPos,double &sumU, + double &sumV) const{ + DVector3 u = ( showerPos-trackPos).Unit(); + DVector3 zhat(0,0,1 ); + DVector3 v = u.Cross(zhat); + sumU=0.; + sumV=0.; + double sumE=0.; + + auto hits=cluster->Get(); + for (const auto& hit:hits){ + DVector2 blockPos=dECALGeom->positionOnFace(hit->row,hit->column); + DVector3 diff(blockPos.X()-showerPos.x(),blockPos.Y()-showerPos.y(),0.); + sumU+=hit->E*pow(u.Dot(diff),2); + sumV+=hit->E*pow(v.Dot(diff),2); + sumE+=hit->E; + } + sumU/=sumE; + sumV/=sumE; +} diff --git a/src/libraries/ECAL/DECALShower_factory.h b/src/libraries/ECAL/DECALShower_factory.h index 72d8423f27..8cc36a9cfe 100644 --- a/src/libraries/ECAL/DECALShower_factory.h +++ b/src/libraries/ECAL/DECALShower_factory.h @@ -11,6 +11,8 @@ #include #include "DECALCluster.h" #include "DECALShower.h" +#include "DECALHit.h" +#include "DECALGeometry.h" class DECALShower_factory:public JFactoryT{ public: @@ -26,6 +28,8 @@ class DECALShower_factory:public JFactoryT{ double GetCorrectedEnergy(double E) const; double GetCorrectedZ(double E) const; + void GetUV(const DECALCluster *cluster,const DVector3 &showerPos, + const DVector3 &trackPos,double &sumU,double &sumV) const; double mECALz,mECALzBack; double SHOWER_ENERGY_THRESHOLD,ECAL_C_EFFECTIVE; @@ -33,6 +37,9 @@ class DECALShower_factory:public JFactoryT{ double E_CORRECTION_PAR1,E_CORRECTION_PAR2,E_CORRECTION_PAR3; double E_CORRECTION_PAR4; bool ENABLE_ENERGY_CORRECTION; + + const DECALGeometry *dECALGeom=NULL; + double mVertexZ; // assume center of target, so it does not vary event by event }; #endif // _DECALShower_factory_ diff --git a/src/libraries/HDDM/DEventSourceREST.cc b/src/libraries/HDDM/DEventSourceREST.cc index dce30c4916..bcdc660673 100644 --- a/src/libraries/HDDM/DEventSourceREST.cc +++ b/src/libraries/HDDM/DEventSourceREST.cc @@ -1167,6 +1167,18 @@ bool DEventSourceREST::Extract_DECALShower(hddm_r::HDDM *record, shower->E1E9=locEcalShowerPropertiesIterator->getE1E9(); shower->E9E25=locEcalShowerPropertiesIterator->getE9E25(); } + const hddm_r::EcalShowerMorePropertiesList& locEcalShowerMorePropertiesList = iter->getEcalShowerMorePropertiesList(); + hddm_r::EcalShowerMorePropertiesList::iterator locEcalShowerMorePropertiesIterator = locEcalShowerMorePropertiesList.begin(); + shower->sumU=0.; + shower->sumV=0.; + shower->docaTrack=0.; + shower->timeTrack=0.; + if(locEcalShowerMorePropertiesIterator != locEcalShowerMorePropertiesList.end()) { + shower->sumU=locEcalShowerMorePropertiesIterator->getSumU(); + shower->sumV=locEcalShowerMorePropertiesIterator->getSumV(); + shower->docaTrack=locEcalShowerMorePropertiesIterator->getDocaTrack(); + shower->timeTrack=locEcalShowerMorePropertiesIterator->getTimeTrack(); + } data.push_back(shower); } @@ -2094,6 +2106,15 @@ bool DEventSourceREST::Extract_DDetectorMatches(const std::shared_ptrSet_DistanceToNearestTrack(locBCALShowers[locShowerIndex], locDeltaPhi, locDeltaZ); } + const hddm_r::EcalDOCAtoTrackList &ecaldocaList = iter->getEcalDOCAtoTracks(); + hddm_r::EcalDOCAtoTrackList::iterator ecaldocaIter = ecaldocaList.begin(); + for(; ecaldocaIter != ecaldocaList.end(); ++ecaldocaIter) + { + size_t locShowerIndex = ecaldocaIter->getShower(); + double locDOCA = ecaldocaIter->getDoca(); + locDetectorMatches->Set_DistanceToNearestTrack(locECALShowers[locShowerIndex], locDOCA); + } + const hddm_r::FcalDOCAtoTrackList &fcaldocaList = iter->getFcalDOCAtoTracks(); hddm_r::FcalDOCAtoTrackList::iterator fcaldocaIter = fcaldocaList.begin(); for(; fcaldocaIter != fcaldocaList.end(); ++fcaldocaIter) diff --git a/src/libraries/HDDM/DEventWriterREST.cc b/src/libraries/HDDM/DEventWriterREST.cc index ddc6c27801..9f182aa18c 100644 --- a/src/libraries/HDDM/DEventWriterREST.cc +++ b/src/libraries/HDDM/DEventWriterREST.cc @@ -282,6 +282,11 @@ bool DEventWriterREST::Write_RESTEvent(const std::shared_ptr& locE hddm_r::EcalShowerPropertiesList locEcalShowerPropertiesList = ecal().addEcalShowerPropertiesList(1); locEcalShowerPropertiesList().setE1E9(ecalshowers[i]->E1E9); locEcalShowerPropertiesList().setE9E25(ecalshowers[i]->E9E25); + hddm_r::EcalShowerMorePropertiesList locEcalShowerMorePropertiesList = ecal().addEcalShowerMorePropertiesList(1); + locEcalShowerMorePropertiesList().setSumU(ecalshowers[i]->sumU); + locEcalShowerMorePropertiesList().setSumV(ecalshowers[i]->sumV); + locEcalShowerMorePropertiesList().setDocaTrack(ecalshowers[i]->docaTrack); + locEcalShowerMorePropertiesList().setTimeTrack(ecalshowers[i]->timeTrack); } // push any DFCALShower objects to the output record for (size_t i=0; i < fcalshowers.size(); i++) @@ -994,6 +999,17 @@ bool DEventWriterREST::Write_RESTEvent(const std::shared_ptr& locE fcalDocaList().setShower(loc_j); fcalDocaList().setDoca(locDistance); } + + for(size_t loc_j = 0; loc_j < ecalshowers.size(); ++loc_j) + { + double locDistance = 0.0; + if(!locDetectorMatches[loc_i]->Get_DistanceToNearestTrack(ecalshowers[loc_j], locDistance)) + continue; + + hddm_r::EcalDOCAtoTrackList ecalDocaList = matches().addEcalDOCAtoTracks(1); + ecalDocaList().setShower(loc_j); + ecalDocaList().setDoca(locDistance); + } } // write the resulting record to the output stream diff --git a/src/libraries/HDDM/rest.xml b/src/libraries/HDDM/rest.xml index 7023011c7f..e3b9eef222 100644 --- a/src/libraries/HDDM/rest.xml +++ b/src/libraries/HDDM/rest.xml @@ -36,6 +36,9 @@ lunit="cm" tunit="ns" Eunit="GeV"> + +