File indexing completed on 2024-04-06 12:11:15
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <memory>
0020
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "DataFormats/Math/interface/Point3D.h"
0023 #include "DataFormats/Math/interface/Vector.h"
0024 #include "DataFormats/Math/interface/Vector3D.h"
0025 #include "FWCore/Framework/interface/ConsumesCollector.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/Framework/interface/stream/EDProducer.h"
0029 #include "FWCore/ParameterSet/interface/FileInPath.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/isFinite.h"
0032 #include "FWCore/Utilities/interface/StreamID.h"
0033 #include "FastSimDataFormats/CTPPSFastSim/interface/CTPPSFastRecHit.h"
0034 #include "FastSimDataFormats/CTPPSFastSim/interface/CTPPSFastRecHitContainer.h"
0035 #include "FastSimDataFormats/CTPPSFastSim/interface/CTPPSFastTrack.h"
0036 #include "FastSimDataFormats/CTPPSFastSim/interface/CTPPSFastTrackContainer.h"
0037 #include "FastSimulation/CTPPSFastGeometry/interface/CTPPSToFDetector.h"
0038 #include "FastSimulation/CTPPSFastGeometry/interface/CTPPSTrkDetector.h"
0039 #include "Utilities/PPS/interface/PPSUnitConversion.h"
0040 #include "Utilities/PPS/interface/PPSUtilities.h"
0041
0042 #include "TLorentzVector.h"
0043
0044
0045 #include "H_Parameters.h"
0046 #include "H_BeamLine.h"
0047 #include "H_RecRPObject.h"
0048 #include "H_BeamParticle.h"
0049
0050 class CTPPSFastTrackingProducer : public edm::stream::EDProducer<> {
0051 public:
0052 explicit CTPPSFastTrackingProducer(const edm::ParameterSet&);
0053
0054 private:
0055 void produce(edm::Event&, const edm::EventSetup&) override;
0056
0057
0058
0059 typedef std::vector<CTPPSFastRecHit> CTPPSFastRecHitContainer;
0060 edm::EDGetTokenT<CTPPSFastRecHitContainer> _recHitToken;
0061 void ReadRecHits(edm::Handle<CTPPSFastRecHitContainer>&);
0062 void FastReco(int Direction, H_RecRPObject* station);
0063 void Reconstruction();
0064 void ReconstructArm(
0065 H_RecRPObject* pps_station, double x1, double y1, double x2, double y2, double& tx, double& ty, double& eloss);
0066 void MatchCellId(int cellId, std::vector<int> vrecCellId, std::vector<double> vrecTof, bool& match, double& recTof);
0067 bool SearchTrack(int,
0068 int,
0069 int Direction,
0070 double& xi,
0071 double& t,
0072 double& partP,
0073 double& pt,
0074 double& thx,
0075 double& thy,
0076 double& x0,
0077 double& y0,
0078 double& xt,
0079 double& yt,
0080 double& X1d,
0081 double& Y1d,
0082 double& X2d,
0083 double& Y2d);
0084 void TrackerStationClear();
0085 void TrackerStationStarting();
0086 void ProjectToToF(const double x1, const double y1, const double x2, const double y2, double& xt, double& yt) {
0087 xt = ((fz_timing - fz_tracker2) * (x2 - x1) / (fz_tracker2 - fz_tracker1)) + x2;
0088 yt = ((fz_timing - fz_tracker2) * (y2 - y1) / (fz_tracker2 - fz_tracker1)) + y2;
0089 };
0090
0091 bool SetBeamLine();
0092
0093 std::unique_ptr<H_BeamLine> m_beamlineCTPPS1;
0094 std::unique_ptr<H_BeamLine> m_beamlineCTPPS2;
0095 std::unique_ptr<H_RecRPObject> pps_stationF;
0096 std::unique_ptr<H_RecRPObject> pps_stationB;
0097
0098 std::string beam1filename;
0099 std::string beam2filename;
0100
0101
0102 double lengthctpps;
0103 bool m_verbosity;
0104 double fBeamEnergy;
0105 double fBeamMomentum;
0106 bool fCrossAngleCorr;
0107 double fCrossingAngleBeam1;
0108 double fCrossingAngleBeam2;
0109
0110 std::unique_ptr<CTPPSTrkStation> TrkStation_F;
0111 std::unique_ptr<CTPPSTrkStation> TrkStation_B;
0112 std::unique_ptr<CTPPSTrkDetector> det1F;
0113 std::unique_ptr<CTPPSTrkDetector> det1B;
0114 std::unique_ptr<CTPPSTrkDetector> det2F;
0115 std::unique_ptr<CTPPSTrkDetector> det2B;
0116 std::unique_ptr<CTPPSToFDetector> detToF_F;
0117 std::unique_ptr<CTPPSToFDetector> detToF_B;
0118
0119 std::vector<CTPPSFastTrack> theCTPPSFastTrack;
0120
0121 CTPPSFastTrack track;
0122
0123 std::vector<int> recCellId_F, recCellId_B;
0124 std::vector<double> recTof_F, recTof_B;
0125
0126 double fz_tracker1, fz_tracker2, fz_timing;
0127 double fTrackerWidth, fTrackerHeight, fTrackerInsertion, fBeamXRMS_Trk1, fBeamXRMS_Trk2, fTrk1XOffset, fTrk2XOffset;
0128 std::vector<double> fToFCellWidth;
0129 double fToFCellHeight, fToFPitchX, fToFPitchY;
0130 int fToFNCellX, fToFNCellY;
0131 double fToFInsertion, fBeamXRMS_ToF, fToFXOffset, fTimeSigma, fImpParcut;
0132 };
0133
0134
0135
0136 CTPPSFastTrackingProducer::CTPPSFastTrackingProducer(const edm::ParameterSet& iConfig)
0137 : m_verbosity(false), fBeamMomentum(0.), fCrossAngleCorr(false), fCrossingAngleBeam1(0.), fCrossingAngleBeam2(0.) {
0138
0139 produces<edm::CTPPSFastTrackContainer>("CTPPSFastTrack");
0140 using namespace edm;
0141 _recHitToken = consumes<CTPPSFastRecHitContainer>(iConfig.getParameter<edm::InputTag>("recHitTag"));
0142 m_verbosity = iConfig.getParameter<bool>("Verbosity");
0143
0144
0145
0146 lengthctpps = iConfig.getParameter<double>("BeamLineLengthCTPPS");
0147 beam1filename = iConfig.getParameter<string>("Beam1");
0148 beam2filename = iConfig.getParameter<string>("Beam2");
0149 fBeamEnergy = iConfig.getParameter<double>("BeamEnergy");
0150 fBeamMomentum = sqrt(fBeamEnergy * fBeamEnergy - PPSTools::ProtonMassSQ);
0151 fCrossingAngleBeam1 = iConfig.getParameter<double>("CrossingAngleBeam1");
0152 fCrossingAngleBeam2 = iConfig.getParameter<double>("CrossingAngleBeam2");
0153
0154 if (fCrossingAngleBeam1 != 0 || fCrossingAngleBeam2 != 0)
0155 fCrossAngleCorr = true;
0156
0157
0158 fz_tracker1 = iConfig.getParameter<double>("Z_Tracker1");
0159 fz_tracker2 = iConfig.getParameter<double>("Z_Tracker2");
0160 fz_timing = iConfig.getParameter<double>("Z_Timing");
0161
0162 fTrackerWidth = iConfig.getParameter<double>("TrackerWidth");
0163 fTrackerHeight = iConfig.getParameter<double>("TrackerHeight");
0164 fTrackerInsertion = iConfig.getParameter<double>("TrackerInsertion");
0165 fBeamXRMS_Trk1 = iConfig.getParameter<double>("BeamXRMS_Trk1");
0166 fBeamXRMS_Trk2 = iConfig.getParameter<double>("BeamXRMS_Trk2");
0167 fTrk1XOffset = iConfig.getParameter<double>("Trk1XOffset");
0168 fTrk2XOffset = iConfig.getParameter<double>("Trk2XOffset");
0169 fToFCellWidth = iConfig.getUntrackedParameter<std::vector<double> >("ToFCellWidth");
0170 fToFCellHeight = iConfig.getParameter<double>("ToFCellHeight");
0171 fToFPitchX = iConfig.getParameter<double>("ToFPitchX");
0172 fToFPitchY = iConfig.getParameter<double>("ToFPitchY");
0173 fToFNCellX = iConfig.getParameter<int>("ToFNCellX");
0174 fToFNCellY = iConfig.getParameter<int>("ToFNCellY");
0175 fToFInsertion = iConfig.getParameter<double>("ToFInsertion");
0176 fBeamXRMS_ToF = iConfig.getParameter<double>("BeamXRMS_ToF");
0177 fToFXOffset = iConfig.getParameter<double>("ToFXOffset");
0178 fTimeSigma = iConfig.getParameter<double>("TimeSigma");
0179 fImpParcut = iConfig.getParameter<double>("ImpParcut");
0180
0181 if (!SetBeamLine()) {
0182 if (m_verbosity)
0183 LogDebug("CTPPSFastTrackingProducer") << "CTPPSFastTrackingProducer: WARNING: lengthctpps= " << lengthctpps;
0184 return;
0185 }
0186
0187
0188
0189
0190
0191 det1F = std::make_unique<CTPPSTrkDetector>(
0192 fTrackerWidth, fTrackerHeight, fTrackerInsertion * fBeamXRMS_Trk1 + fTrk1XOffset);
0193 det2F = std::make_unique<CTPPSTrkDetector>(
0194 fTrackerWidth, fTrackerHeight, fTrackerInsertion * fBeamXRMS_Trk2 + fTrk2XOffset);
0195 det1B = std::make_unique<CTPPSTrkDetector>(
0196 fTrackerWidth, fTrackerHeight, fTrackerInsertion * fBeamXRMS_Trk1 + fTrk1XOffset);
0197 det2B = std::make_unique<CTPPSTrkDetector>(
0198 fTrackerWidth, fTrackerHeight, fTrackerInsertion * fBeamXRMS_Trk2 + fTrk2XOffset);
0199
0200
0201 std::vector<double> vToFCellWidth;
0202 vToFCellWidth.reserve(8);
0203 for (int i = 0; i < 8; i++) {
0204 vToFCellWidth.push_back(fToFCellWidth[i]);
0205 }
0206 double pos_tof = fToFInsertion * fBeamXRMS_ToF + fToFXOffset;
0207 detToF_F = std::make_unique<CTPPSToFDetector>(
0208 fToFNCellX, fToFNCellY, vToFCellWidth, fToFCellHeight, fToFPitchX, fToFPitchY, pos_tof, fTimeSigma);
0209 detToF_B = std::make_unique<CTPPSToFDetector>(
0210 fToFNCellX, fToFNCellY, vToFCellWidth, fToFCellHeight, fToFPitchX, fToFPitchY, pos_tof, fTimeSigma);
0211
0212 }
0213
0214
0215 void CTPPSFastTrackingProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0216 using namespace edm;
0217 TrackerStationStarting();
0218 Handle<CTPPSFastRecHitContainer> recHits;
0219 iEvent.getByToken(_recHitToken, recHits);
0220 recCellId_F.clear();
0221 recCellId_B.clear();
0222 recTof_F.clear();
0223 recTof_B.clear();
0224 ReadRecHits(recHits);
0225 Reconstruction();
0226 TrackerStationClear();
0227
0228 std::unique_ptr<CTPPSFastTrackContainer> output_tracks(new CTPPSFastTrackContainer);
0229 for (std::vector<CTPPSFastTrack>::const_iterator i = theCTPPSFastTrack.begin(); i != theCTPPSFastTrack.end(); i++) {
0230 output_tracks->push_back(*i);
0231 }
0232
0233 iEvent.put(std::move(output_tracks), "CTPPSFastTrack");
0234 }
0235
0236
0237 void CTPPSFastTrackingProducer::TrackerStationClear() {
0238 TrkStation_F->first.clear();
0239 TrkStation_F->second.clear();
0240 TrkStation_B->first.clear();
0241 TrkStation_B->second.clear();
0242 }
0243
0244 void CTPPSFastTrackingProducer::TrackerStationStarting() {
0245 det1F->clear();
0246 det1B->clear();
0247 det2F->clear();
0248 det2B->clear();
0249 detToF_F->clear();
0250 detToF_B->clear();
0251 }
0252
0253
0254 void CTPPSFastTrackingProducer::ReadRecHits(edm::Handle<CTPPSFastRecHitContainer>& recHits) {
0255
0256
0257
0258
0259
0260
0261
0262
0263 for (unsigned int irecHits = 0; irecHits < recHits->size(); ++irecHits) {
0264 const CTPPSFastRecHit* recHitDet = &(*recHits)[irecHits];
0265 unsigned int detlayerId = recHitDet->detUnitId();
0266 double x = recHitDet->entryPoint().x();
0267 double y = recHitDet->entryPoint().y();
0268 double z = recHitDet->entryPoint().z();
0269 float tof = recHitDet->tof();
0270 if (detlayerId == 2014314496)
0271 det1F->AddHit(detlayerId, x, y, z);
0272 else if (detlayerId == 2014838784)
0273 det2F->AddHit(detlayerId, x, y, z);
0274 else if (detlayerId == 2031091712)
0275 det1B->AddHit(detlayerId, x, y, z);
0276 else if (detlayerId == 2031616000)
0277 det2B->AddHit(detlayerId, x, y, z);
0278 else if (detlayerId == 2046820352) {
0279 detToF_F->AddHit(x, y, tof);
0280 recCellId_F.push_back(detToF_F->findCellId(x, y));
0281 recTof_F.push_back(tof);
0282 } else if (detlayerId == 2063597568) {
0283 detToF_B->AddHit(x, y, tof);
0284 recCellId_B.push_back(detToF_B->findCellId(x, y));
0285 recTof_B.push_back(tof);
0286 }
0287
0288 }
0289
0290 TrkStation_F = std::make_unique<CTPPSTrkStation>(*det1F, *det2F);
0291 TrkStation_B = std::make_unique<CTPPSTrkStation>(*det1B, *det2B);
0292 }
0293
0294 void CTPPSFastTrackingProducer::Reconstruction() {
0295 theCTPPSFastTrack.clear();
0296 int Direction;
0297 Direction = 1;
0298 FastReco(Direction, &*pps_stationF);
0299 Direction = -1;
0300 FastReco(Direction, &*pps_stationB);
0301 }
0302
0303 bool CTPPSFastTrackingProducer::SearchTrack(int i,
0304 int j,
0305 int Direction,
0306 double& xi,
0307 double& t,
0308 double& partP,
0309 double& pt,
0310 double& thx,
0311 double& thy,
0312 double& x0,
0313 double& y0,
0314 double& xt,
0315 double& yt,
0316 double& X1d,
0317 double& Y1d,
0318 double& X2d,
0319 double& Y2d) {
0320
0321 xi = 0;
0322 t = 0;
0323 partP = 0;
0324 pt = 0;
0325 x0 = 0.;
0326 y0 = 0.;
0327 xt = 0.;
0328 yt = 0.;
0329 X1d = 0.;
0330 Y1d = 0.;
0331 X2d = 0.;
0332 Y2d = 0.;
0333 CTPPSTrkDetector* det1 = nullptr;
0334 CTPPSTrkDetector* det2 = nullptr;
0335 H_RecRPObject* station = nullptr;
0336
0337 if (Direction > 0) {
0338 det1 = &(TrkStation_F->first);
0339 det2 = &(TrkStation_F->second);
0340 station = &*pps_stationF;
0341 } else {
0342 det1 = &(TrkStation_B->first);
0343 det2 = &(TrkStation_B->second);
0344 station = &*pps_stationB;
0345 }
0346 if (det1->ppsNHits_ <= i || det2->ppsNHits_ <= j)
0347 return false;
0348
0349 double x1 = det1->ppsX_.at(i);
0350 double y1 = det1->ppsY_.at(i);
0351 double x2 = det2->ppsX_.at(j);
0352 double y2 = det2->ppsY_.at(j);
0353 double eloss;
0354
0355
0356 ReconstructArm(
0357 station, Direction * x1, y1, Direction * x2, y2, thx, thy, eloss);
0358 thx *= -Direction;
0359
0360
0361 if (edm::isNotFinite(eloss) || edm::isNotFinite(thx) || edm::isNotFinite(thy))
0362 return false;
0363
0364
0365 if (m_verbosity)
0366 LogDebug("CTPPSFastTrackingProducer::SearchTrack:") << "thx " << thx << " thy " << thy << " eloss " << eloss;
0367
0368
0369 x0 = -Direction * station->getX0() * um_to_cm;
0370 y0 = station->getY0() * um_to_cm;
0371 double ImpPar = sqrt(x0 * x0 + y0 * y0);
0372 if (ImpPar > fImpParcut)
0373 return false;
0374 if (eloss < 0. || eloss > fBeamEnergy)
0375 return false;
0376
0377
0378 double theta = sqrt(thx * thx + thy * thy) * urad;
0379 xi = eloss / fBeamEnergy;
0380 double energy = fBeamEnergy * (1. - xi);
0381 partP = sqrt(energy * energy - PPSTools::ProtonMassSQ);
0382 t = -2. * (PPSTools::ProtonMassSQ - fBeamEnergy * energy + fBeamMomentum * partP * cos(theta));
0383 pt = partP * theta;
0384 if (xi < 0. || xi > 1. || t < 0. || t > 10. || pt <= 0.) {
0385 xi = 0.;
0386 t = 0.;
0387 partP = 0.;
0388 pt = 0.;
0389 x0 = 0.;
0390 y0 = 0.;
0391 return false;
0392 }
0393
0394 ProjectToToF(x1, y1, x2, y2, xt, yt);
0395 X1d = x1;
0396 Y1d = y1;
0397 X2d = x2;
0398 Y2d = y2;
0399 return true;
0400 }
0401
0402 void CTPPSFastTrackingProducer::ReconstructArm(
0403 H_RecRPObject* pps_station, double x1, double y1, double x2, double y2, double& tx, double& ty, double& eloss) {
0404 tx = 0.;
0405 ty = 0.;
0406 eloss = 0.;
0407 if (!pps_station)
0408 return;
0409 x1 *= mm_to_um;
0410 x2 *= mm_to_um;
0411 y1 *= mm_to_um;
0412 y2 *= mm_to_um;
0413 pps_station->setPositions(x1, y1, x2, y2);
0414 double energy = pps_station->getE(AM);
0415 if (edm::isNotFinite(energy))
0416 return;
0417 tx = pps_station->getTXIP();
0418 ty = pps_station->getTYIP();
0419 eloss = pps_station->getE();
0420 }
0421
0422 void CTPPSFastTrackingProducer::MatchCellId(
0423 int cellId, std::vector<int> vrecCellId, std::vector<double> vrecTof, bool& match, double& recTof) {
0424 for (unsigned int i = 0; i < vrecCellId.size(); i++) {
0425 if (cellId == vrecCellId.at(i)) {
0426 match = true;
0427 recTof = vrecTof.at(i);
0428 continue;
0429 }
0430 }
0431 }
0432
0433 void CTPPSFastTrackingProducer::FastReco(int Direction, H_RecRPObject* station) {
0434 double theta = 0.;
0435 double xi, t, partP, pt, phi, x0, y0, thx, thy, xt, yt, X1d, Y1d, X2d, Y2d;
0436 CTPPSTrkDetector* Trk1 = nullptr;
0437 CTPPSTrkDetector* Trk2 = nullptr;
0438 double pos_tof = fToFInsertion * fBeamXRMS_ToF + fToFXOffset;
0439 int cellId = 0;
0440 std::vector<double> vToFCellWidth;
0441 vToFCellWidth.reserve(8);
0442 for (int i = 0; i < 8; i++) {
0443 vToFCellWidth.push_back(fToFCellWidth[i]);
0444 }
0445 CTPPSToFDetector* ToF = new CTPPSToFDetector(
0446 fToFNCellX, fToFNCellY, vToFCellWidth, fToFCellHeight, fToFPitchX, fToFPitchY, pos_tof, fTimeSigma);
0447 if (Direction > 0) {
0448 Trk1 = &(TrkStation_F->first);
0449 Trk2 = &(TrkStation_F->second);
0450 } else {
0451 Trk1 = &(TrkStation_B->first);
0452 Trk2 = &(TrkStation_B->second);
0453 }
0454
0455
0456 for (int i = 0; i < (int)Trk1->ppsNHits_; i++) {
0457 for (int j = 0; j < (int)Trk2->ppsNHits_; j++) {
0458 if (SearchTrack(i, j, Direction, xi, t, partP, pt, thx, thy, x0, y0, xt, yt, X1d, Y1d, X2d, Y2d)) {
0459
0460 cellId = ToF->findCellId(xt, yt);
0461 double recTof = 0.;
0462 bool matchCellId = false;
0463 if (Direction > 0) {
0464 theta = sqrt(thx * thx + thy * thy) * urad;
0465 MatchCellId(cellId, recCellId_F, recTof_F, matchCellId, recTof);
0466 } else if (Direction < 0) {
0467 theta = TMath::Pi() - sqrt(thx * thx + thy * thy) * urad;
0468 MatchCellId(cellId, recCellId_B, recTof_B, matchCellId, recTof);
0469 }
0470 phi = atan2(thy, thx);
0471
0472 double px = partP * sin(theta) * cos(phi);
0473 double py = partP * sin(theta) * sin(phi);
0474 double pz = partP * cos(theta);
0475 double e = sqrt(partP * partP + PPSTools::ProtonMassSQ);
0476 TLorentzVector p(px, py, pz, e);
0477
0478 if (fCrossAngleCorr) {
0479 PPSTools::LorentzBoost(p, "MC", {fCrossingAngleBeam1, fCrossingAngleBeam2, fBeamMomentum, fBeamEnergy});
0480 }
0481
0482 PPSTools::Get_t_and_xi(const_cast<TLorentzVector*>(&p), t, xi, {fBeamMomentum, fBeamEnergy});
0483 double pxx = p.Px();
0484 double pyy = p.Py();
0485 double pzz = p.Pz();
0486 math::XYZVector momentum(pxx, pyy, pzz);
0487 math::XYZPoint vertex(x0, y0, 0);
0488
0489 track.setp(momentum);
0490 track.setvertex(vertex);
0491 track.sett(t);
0492 track.setxi(xi);
0493 track.setx1(X1d);
0494 track.sety1(Y1d);
0495 track.setx2(X2d);
0496 track.sety2(Y2d);
0497 if (matchCellId) {
0498 track.setcellid(cellId);
0499 track.settof(recTof);
0500 } else {
0501 track.setcellid(0);
0502 track.settof(0.);
0503 }
0504 theCTPPSFastTrack.push_back(track);
0505 }
0506 }
0507 }
0508 }
0509
0510 bool CTPPSFastTrackingProducer::SetBeamLine() {
0511 edm::FileInPath b1(beam1filename.c_str());
0512 edm::FileInPath b2(beam2filename.c_str());
0513 if (lengthctpps <= 0)
0514 return false;
0515 m_beamlineCTPPS1 = std::make_unique<H_BeamLine>(-1, lengthctpps + 0.1);
0516 m_beamlineCTPPS1->fill(b2.fullPath(), 1, "IP5");
0517 m_beamlineCTPPS2 = std::make_unique<H_BeamLine>(1, lengthctpps + 0.1);
0518 m_beamlineCTPPS2->fill(b1.fullPath(), 1, "IP5");
0519 m_beamlineCTPPS1->offsetElements(120, 0.097);
0520 m_beamlineCTPPS2->offsetElements(120, -0.097);
0521 pps_stationF = std::make_unique<H_RecRPObject>(fz_tracker1, fz_tracker2, *m_beamlineCTPPS1);
0522 pps_stationB = std::make_unique<H_RecRPObject>(fz_tracker1, fz_tracker2, *m_beamlineCTPPS2);
0523 return true;
0524 }
0525
0526 DEFINE_FWK_MODULE(CTPPSFastTrackingProducer);