CMS 3D CMS Logo

CTPPSFastTrackingProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: FastSimulation/CTPPSFastTrackingProducer
4 // Class: CTPPSFastTrackingProducer
5 //
13 //
14 // Original Author: Sandro Fonseca De Souza
15 // Created: Thu, 29 Sep 2016 16:13:41 GMT
16 //
17 //
18 
19 
20 // system include files
21 #include <memory>
22 
23 // user include files
26 #include "TLorentzVector.h"
27 // hector includes
28 #include "H_Parameters.h"
29 #include "H_BeamLine.h"
30 #include "H_RecRPObject.h"
31 #include "H_BeamParticle.h"
33 // constructors and destructor
34 //
36  m_verbosity(false), fBeamMomentum(0.), fCrossAngleCorr(false), fCrossingAngleBeam1(0.),fCrossingAngleBeam2(0.)
37 {
38  //register your products
39  produces<edm::CTPPSFastTrackContainer>("CTPPSFastTrack");
40  using namespace edm;
41  _recHitToken = consumes<CTPPSFastRecHitContainer>(iConfig.getParameter<edm::InputTag>("recHitTag"));
42  m_verbosity = iConfig.getParameter<bool>("Verbosity");
43  // User definitons
44 
45  // Read beam parameters needed for Hector reconstruction
46  lengthctpps = iConfig.getParameter<double>("BeamLineLengthCTPPS" );
47  beam1filename = iConfig.getParameter<string>("Beam1");
48  beam2filename = iConfig.getParameter<string>("Beam2");
49  fBeamEnergy = iConfig.getParameter<double>("BeamEnergy"); // beam energy in GeV
50  fBeamMomentum = sqrt(fBeamEnergy*fBeamEnergy - PPSTools::ProtonMassSQ);
51  fCrossingAngleBeam1 = iConfig.getParameter<double>("CrossingAngleBeam1");
52  fCrossingAngleBeam2 = iConfig.getParameter<double>("CrossingAngleBeam2");
53 
54  if (fCrossingAngleBeam1!=0||fCrossingAngleBeam2!=0) fCrossAngleCorr = true;
55  //Read detectors positions and parameters
56 
57  fz_tracker1 = iConfig.getParameter<double>("Z_Tracker1");
58  fz_tracker2 = iConfig.getParameter<double>("Z_Tracker2");
59  fz_timing = iConfig.getParameter<double>("Z_Timing");
60  //
61  fTrackerWidth = iConfig.getParameter<double>("TrackerWidth");
62  fTrackerHeight = iConfig.getParameter<double>("TrackerHeight");
63  fTrackerInsertion = iConfig.getParameter<double>("TrackerInsertion");
64  fBeamXRMS_Trk1 = iConfig.getParameter<double>("BeamXRMS_Trk1");
65  fBeamXRMS_Trk2 = iConfig.getParameter<double>("BeamXRMS_Trk2");
66  fTrk1XOffset = iConfig.getParameter<double>("Trk1XOffset");
67  fTrk2XOffset = iConfig.getParameter<double>("Trk2XOffset");
68  fToFCellWidth = iConfig.getUntrackedParameter<std::vector<double> >("ToFCellWidth");
69  fToFCellHeight = iConfig.getParameter<double>("ToFCellHeight");
70  fToFPitchX = iConfig.getParameter<double>("ToFPitchX");
71  fToFPitchY = iConfig.getParameter<double>("ToFPitchY");
72  fToFNCellX = iConfig.getParameter<int>("ToFNCellX");
73  fToFNCellY = iConfig.getParameter<int>("ToFNCellY");
74  fToFInsertion = iConfig.getParameter<double>("ToFInsertion");
75  fBeamXRMS_ToF = iConfig.getParameter<double>("BeamXRMS_ToF");
76  fToFXOffset = iConfig.getParameter<double>("ToFXOffset");
77  fTimeSigma = iConfig.getParameter<double>("TimeSigma");
78  fImpParcut = iConfig.getParameter<double>("ImpParcut");
79 
80 
81  if(!SetBeamLine() ) {
82  if ( m_verbosity ) LogDebug("CTPPSFastTrackingProducer") << "CTPPSFastTrackingProducer: WARNING: lengthctpps= " << lengthctpps;
83  return;
84  }
85 
86  // Create a particle to get the beam energy from the beam file
87  // Take care: the z inside the station is in meters
88  //
89  //Tracker Detector Description
90  det1F =std::unique_ptr<CTPPSTrkDetector>(new CTPPSTrkDetector(fTrackerWidth,fTrackerHeight,fTrackerInsertion*fBeamXRMS_Trk1+fTrk1XOffset));
91  det2F =std::unique_ptr<CTPPSTrkDetector>(new CTPPSTrkDetector(fTrackerWidth,fTrackerHeight,fTrackerInsertion*fBeamXRMS_Trk2+fTrk2XOffset));
92  det1B =std::unique_ptr<CTPPSTrkDetector>(new CTPPSTrkDetector(fTrackerWidth,fTrackerHeight,fTrackerInsertion*fBeamXRMS_Trk1+fTrk1XOffset));
93  det2B =std::unique_ptr<CTPPSTrkDetector>(new CTPPSTrkDetector(fTrackerWidth,fTrackerHeight,fTrackerInsertion*fBeamXRMS_Trk2+fTrk2XOffset));
94 
95  //Timing Detector Description
96  std::vector<double> vToFCellWidth;
97  for (int i = 0 ; i < 8 ; i++){
98  vToFCellWidth.push_back(fToFCellWidth[i]);
99  }
100  double pos_tof = fToFInsertion*fBeamXRMS_ToF+fToFXOffset;
101  detToF_F =std::unique_ptr<CTPPSToFDetector>(new CTPPSToFDetector(fToFNCellX,fToFNCellY,vToFCellWidth,fToFCellHeight,fToFPitchX,fToFPitchY,pos_tof,fTimeSigma));
102  detToF_B =std::unique_ptr<CTPPSToFDetector>(new CTPPSToFDetector(fToFNCellX,fToFNCellY,vToFCellWidth,fToFCellHeight,fToFPitchX,fToFPitchY,pos_tof,fTimeSigma));
103 //
104 }
106 {
107  for (std::map<unsigned int,H_BeamParticle*>::iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it ) {
108  delete (*it).second;
109  }
110 }
111 // ------------ method called to produce the data ------------
112  void
114 {
115  using namespace edm;
118  iEvent.getByToken( _recHitToken,recHits);
119  recCellId_F.clear(); recCellId_B.clear();
120  recTof_F.clear(); recTof_B.clear();
121  ReadRecHits(recHits);
122  Reconstruction();
124 
125  std::unique_ptr<CTPPSFastTrackContainer> output_tracks(new CTPPSFastTrackContainer);
126  int n = 0;
127  for ( std::vector<CTPPSFastTrack>::const_iterator i = theCTPPSFastTrack.begin();
128  i != theCTPPSFastTrack.end(); i++ ) {
129  output_tracks->push_back(*i);
130  n += 1;
131  }
132 
133  iEvent.put(std::move(output_tracks),"CTPPSFastTrack");
134 }//end
136 {
138 
139 }
142 {
144 }
145 
148 
149  TrkStation_F->first.clear(); TrkStation_F->second.clear();
150  TrkStation_B->first.clear(); TrkStation_B->second.clear();
151 }
154  det1F->clear();
155  det1B->clear();
156  det2F->clear();
157  det2B->clear();
158  detToF_F->clear();
159  detToF_B->clear();
160 
161 }
162 
163 
165  void
167 {
168  // DetId codification for PSimHit taken from CTPPSPixel- It will be replaced by CTPPSDetId
169  // 2014314496 -> Tracker1 zPositive
170  // 2014838784 -> Tracker2 zPositive
171  // 2046820352 -> Timing zPositive
172  // 2031091712 -> Tracker1 zNegative
173  // 2031616000 -> Tracker2 zNegative
174  // 2063597568 -> Timing zNegative
175 
176  for (unsigned int irecHits = 0; irecHits < recHits->size(); ++irecHits)
177  {
178  const CTPPSFastRecHit* recHitDet = &(*recHits)[irecHits];
179  unsigned int detlayerId = recHitDet->detUnitId();
180  double x = recHitDet->entryPoint().x();
181  double y = recHitDet->entryPoint().y();
182  double z = recHitDet->entryPoint().z();
183  float tof = recHitDet->tof();
184  if(detlayerId == 2014314496) det1F->AddHit(detlayerId,x,y,z);
185  else if(detlayerId == 2014838784) det2F->AddHit(detlayerId,x,y,z);
186  else if(detlayerId == 2031091712) det1B->AddHit(detlayerId,x,y,z);
187  else if(detlayerId == 2031616000) det2B->AddHit(detlayerId,x,y,z);
188  else if(detlayerId == 2046820352) { detToF_F->AddHit(x,y,tof); recCellId_F.push_back(detToF_F->findCellId(x,y)); recTof_F.push_back(tof); }
189  else if(detlayerId == 2063597568) { detToF_B->AddHit(x,y,tof); recCellId_B.push_back(detToF_B->findCellId(x,y)); recTof_B.push_back(tof); }
190 
191  }//LOOP TRK
192  //creating Stations
193  TrkStation_F = std::unique_ptr<CTPPSTrkStation>(new std::pair<CTPPSTrkDetector,CTPPSTrkDetector>(*det1F,*det2F));
194  TrkStation_B = std::unique_ptr<CTPPSTrkStation>(new std::pair<CTPPSTrkDetector,CTPPSTrkDetector>(*det1B,*det2B));
195 } // end function
196 
198 {
199  theCTPPSFastTrack.clear();
200  int Direction;
201  Direction=1; //cms positive Z / forward
202  FastReco(Direction,&*pps_stationF);
203  Direction=-1; //cms negative Z / backward
204  FastReco(Direction,&*pps_stationB);
205 }//end Reconstruction
206 
207 bool CTPPSFastTrackingProducer::SearchTrack(int i,int j,int Direction,double& xi,double& t,double& partP,double& pt,double& thx,double& thy,double& x0, double& y0, double& xt, double& yt, double& X1d, double& Y1d, double& X2d, double& Y2d)
208 {
209  // Given 1 hit in Tracker1 and 1 hit in Tracker2 try to make a track with Hector
210  double theta=0.;
211  xi = 0; t=0; partP=0; pt=0; x0=0.;y0=0.;xt =0.;yt =0.;X1d=0.;Y1d=0.;X2d=0.;Y2d=0.;
212  CTPPSTrkDetector* det1 = nullptr;
213  CTPPSTrkDetector* det2 = nullptr;
214  H_RecRPObject* station = nullptr;
215  // Separate in forward and backward stations according to direction
216  if (Direction>0) {
217  det1=&(TrkStation_F->first);det2=&(TrkStation_F->second);
218  station = &*pps_stationF;
219  } else {
220  det1=&(TrkStation_B->first);det2=&(TrkStation_B->second);
221  station = &*pps_stationB;
222  }
223  if (det1->ppsNHits_<=i||det2->ppsNHits_<=j) return false;
224 //
225  double x1 = det1->ppsX_.at(i); double y1 = det1->ppsY_.at(i);
226  double x2 = det2->ppsX_.at(j); double y2 = det2->ppsY_.at(j);
227  double eloss;
228 
229  //thx and thy are returned in microrad
230  ReconstructArm(station,Direction*x1,y1,Direction*x2,y2,thx,thy,eloss); // Pass the hits in the LHC ref. frame
231  thx*=-Direction;// invert to the CMS ref frame
232 
233  // Protect for unphysical results
234  if (std::isnan(eloss)||std::isinf(eloss)||
235  std::isnan(thx) || std::isinf(thx) ||
236  std::isnan(thy) || std::isinf(thy)) return false;
237  //
238 
239  if ( m_verbosity ) LogDebug("CTPPSFastTrackingProducer::SearchTrack:") << "thx " << thx << " thy " << thy << " eloss " << eloss;
240 
241  // Get the start point of the reconstructed track near the origin made by Hector in the CMS ref. frame
242  x0 = -Direction*station->getX0()*um_to_cm;
243  y0 = station->getY0()*um_to_cm;
244  double ImpPar=sqrt(x0*x0+y0*y0);
245  if (ImpPar>fImpParcut) return false;
246  if (eloss<0.||eloss>fBeamEnergy) return false;
247 //
248  // Calculate the reconstructed track parameters
249  theta = sqrt(thx*thx+thy*thy)*urad;
250  xi = eloss/fBeamEnergy;
251  double energy= fBeamEnergy*(1.-xi);
252  partP = sqrt(energy*energy-PPSTools::ProtonMassSQ);
253  t = -2.*(PPSTools::ProtonMassSQ- fBeamEnergy*energy + fBeamMomentum*partP*cos(theta));
254  pt = sqrt(pow(partP*thx*urad,2)+pow(partP*thy*urad,2));
255  if (xi<0.||xi>1.||t<0.||t>10.||pt<=0.) {
256  xi = 0.; t=0.; partP=0.; pt=0.; theta=0.; x0=0.;y0=0.;
257  return false; // unphysical values
258  }
259  //Try to include the timing detector in the track
260  ProjectToToF(x1, y1, x2, y2, xt, yt); // the projections is done in the CMS ref frame
261  X1d = x1;
262  Y1d = y1;
263  X2d = x2;
264  Y2d = y2;
265  return true;
266 }//end SearchTrack
267 
268 void CTPPSFastTrackingProducer::ReconstructArm(H_RecRPObject* pps_station, double x1, double y1, double x2, double y2, double& tx, double& ty, double& eloss)
269 {
270  tx=0.;
271  ty=0.;
272  eloss=0.;
273  if (!pps_station) return;
274  x1*=mm_to_um;
275  x2*=mm_to_um;
276  y1*=mm_to_um;
277  y2*=mm_to_um;
278  pps_station->setPositions(x1,y1,x2,y2);
279  double energy = pps_station->getE(AM); // dummy call needed to calculate some Hector internal parameter
280  if (std::isnan(energy)||std::isinf(energy)) return;
281  tx = pps_station->getTXIP(); // change orientation to CMS
282  ty = pps_station->getTYIP();
283  eloss = pps_station->getE();
284 }
285 
286 void CTPPSFastTrackingProducer::MatchCellId(int cellId, std::vector<int> vrecCellId, std::vector<double> vrecTof, bool& match, double& recTof){
287  for (unsigned int i = 0 ; i < vrecCellId.size(); i++){
288  if(cellId == vrecCellId.at(i)) {
289  match = true;
290  recTof = vrecTof.at(i);
291  continue;
292  }
293  }
294 }
295 
296 void CTPPSFastTrackingProducer::FastReco(int Direction,H_RecRPObject* station)
297 {
298  double theta = 0.;
299  double xi,t,partP,pt,phi,x0,y0,thx,thy,xt,yt,X1d,Y1d,X2d,Y2d;
300  CTPPSTrkDetector* Trk1 = nullptr;
301  CTPPSTrkDetector* Trk2 = nullptr;
302  double pos_tof = fToFInsertion*fBeamXRMS_ToF+fToFXOffset;
303  int cellId = 0;
304  std::vector<double> vToFCellWidth;
305  for (int i = 0 ; i < 8 ; i++){
306  vToFCellWidth.push_back(fToFCellWidth[i]);
307  }
309  if (Direction>0) {
310  Trk1=&(TrkStation_F->first);Trk2=&(TrkStation_F->second);
311  } else {
312  Trk1=&(TrkStation_B->first);Trk2=&(TrkStation_B->second);
313  }
314  // Make a track from EVERY pair of hits combining Tracker1 and Tracker2.
315  // The tracks may not be independent as 1 hit may belong to more than 1 track.
316  for(int i=0;i<(int)Trk1->ppsNHits_;i++) {
317  for(int j=0;j<(int)Trk2->ppsNHits_;j++){
318  if (SearchTrack(i,j,Direction,xi,t,partP,pt,thx,thy,x0,y0,xt,yt,X1d,Y1d,X2d,Y2d)) {
319  // Check if the hitted timing cell matches the reconstructed track
320  cellId = ToF->findCellId(xt,yt);
321  double recTof = 0.;
322  bool matchCellId = false;
323  if (Direction > 0 ) {
324  theta = sqrt(thx*thx+thy*thy)*urad;
325  MatchCellId(cellId, recCellId_F, recTof_F, matchCellId, recTof);
326  }
327  else if (Direction<0) {
328  theta = TMath::Pi() - sqrt(thx*thx+thy*thy)*urad;
329  MatchCellId(cellId, recCellId_B, recTof_B, matchCellId, recTof);
330  }
331  phi = atan2(thy,thx); // at this point, thx is already in the cms ref. frame
332 
333  double px = partP*sin(theta)*cos(phi);
334  double py = partP*sin(theta)*sin(phi);
335  double pz = partP*cos(theta);
336  double e = sqrt(partP*partP+PPSTools::ProtonMassSQ);
337  TLorentzVector p(px,py,pz,e);
338  // Invert the Lorentz boost made to take into account the crossing angle during simulation
339  if (fCrossAngleCorr) {
342  }
343  //Getting the Xi and t (squared four momentum transferred) of the reconstructed track
344  PPSTools::Get_t_and_xi(const_cast<TLorentzVector*>(&p),t,xi, {fBeamMomentum,fBeamEnergy});
345  double pxx = p.Px(); double pyy = p.Py(); double pzz = p.Pz();
346  math::XYZVector momentum (pxx,pyy,pzz);
347  math::XYZPoint vertex (x0,y0,0);
348 
349  track.setp(momentum);
350  track.setvertex(vertex);
351  track.sett(t);
352  track.setxi(xi);
353  track.setx1(X1d);
354  track.sety1(Y1d);
355  track.setx2(X2d);
356  track.sety2(Y2d);
357  if (matchCellId) {
358  track.setcellid(cellId);
359  track.settof(recTof);
360  }
361  else {
362  track.setcellid(0);
363  track.settof(0.);
364  }
365  theCTPPSFastTrack.push_back(track);
366  }
367  }
368  }
369 }//end FastReco
370 
371 // ------------ method called once each stream before processing any runs, lumis or events ------------
372 
374 {
375 }
376 
377 // ------------ method called once each stream after processing all runs, lumis and events ------------
378 
380 {
381 }
383 {
384  edm::FileInPath b1(beam1filename.c_str());
385  edm::FileInPath b2(beam2filename.c_str());
386  if (lengthctpps<=0) return false;
387  m_beamlineCTPPS1 = std::unique_ptr<H_BeamLine>(new H_BeamLine( -1, lengthctpps + 0.1)); // (direction, length)
388  m_beamlineCTPPS1->fill( b2.fullPath(), 1, "IP5");
389  m_beamlineCTPPS2 = std::unique_ptr<H_BeamLine>(new H_BeamLine( 1, lengthctpps + 0.1 )); //
390  m_beamlineCTPPS2->fill( b1.fullPath(), 1, "IP5");
391  m_beamlineCTPPS1->offsetElements( 120, 0.097 );
392  m_beamlineCTPPS2->offsetElements( 120,-0.097 );
393  pps_stationF = std::unique_ptr<H_RecRPObject>(new H_RecRPObject(fz_tracker1,fz_tracker2,*m_beamlineCTPPS1));
394  pps_stationB = std::unique_ptr<H_RecRPObject>(new H_RecRPObject(fz_tracker1,fz_tracker2,*m_beamlineCTPPS2));
395  return true;
396 }
397 //define this as a plug-in
#define LogDebug(id)
void LorentzBoost(H_BeamParticle &h_p, int dir, const std::string &frame, FullBeamInfo const &bi)
void sett(float t)
const double Pi
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
void beginStream(edm::StreamID) override
std::vector< double > ppsX_
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:125
std::unique_ptr< CTPPSTrkDetector > det2B
int findCellId(double x, double y)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
std::unique_ptr< H_BeamLine > m_beamlineCTPPS1
virtual void endEvent(edm::Event &event, const edm::EventSetup &eventSetup)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void FastReco(int Direction, H_RecRPObject *station)
Geom::Theta< T > theta() const
T y() const
Definition: PV3DBase.h:63
void setvertex(const Point &vertex)
static const double urad
void sety1(float y1)
void setx2(float x2)
void ReconstructArm(H_RecRPObject *pps_station, double x1, double y1, double x2, double y2, double &tx, double &ty, double &eloss)
std::unique_ptr< H_RecRPObject > pps_stationB
std::unique_ptr< CTPPSToFDetector > detToF_F
float tof() const
deprecated name for timeOfFlight()
std::unique_ptr< CTPPSTrkDetector > det2F
std::unique_ptr< H_RecRPObject > pps_stationF
void setp(const Vector &momentum)
int iEvent
Definition: GenABIO.cc:224
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void setx1(float x1)
CTPPSFastTrackingProducer(const edm::ParameterSet &)
std::unique_ptr< CTPPSTrkDetector > det1B
bool isnan(float x)
Definition: math.h:13
virtual void beginEvent(edm::Event &event, const edm::EventSetup &eventSetup)
T sqrt(T t)
Definition: SSEVec.h:18
unsigned int detUnitId() const
T z() const
Definition: PV3DBase.h:64
std::unique_ptr< CTPPSTrkStation > TrkStation_B
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
void setxi(float xi)
const double ProtonMassSQ
Definition: PPSUtilities.h:33
std::unique_ptr< CTPPSToFDetector > detToF_B
std::map< unsigned int, H_BeamParticle * > m_beamPart
std::unique_ptr< H_BeamLine > m_beamlineCTPPS2
void ProjectToToF(const double x1, const double y1, const double x2, const double y2, double &xt, double &yt)
std::vector< double > ppsY_
void MatchCellId(int cellId, std::vector< int > vrecCellId, std::vector< double > vrecTof, bool &match, double &recTof)
std::unique_ptr< CTPPSTrkStation > TrkStation_F
std::vector< CTPPSFastTrack > theCTPPSFastTrack
const GeomDet * recHitDet(const TrackingRecHit &hit, const TrackingGeometry *geom)
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
static const double um_to_cm
edm::EDGetTokenT< CTPPSFastRecHitContainer > _recHitToken
Local3DPoint entryPoint() const
Entry point in the local Det frame.
std::unique_ptr< CTPPSTrkDetector > det1F
void settof(float tof)
void sety2(float y2)
static const double mm_to_um
HLT enums.
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
void Get_t_and_xi(const TLorentzVector *proton, double &t, double &xi, LimitedBeamInfo const &bi)
Definition: PPSUtilities.cc:65
T x() const
Definition: PV3DBase.h:62
void produce(edm::Event &, const edm::EventSetup &) override
std::vector< CTPPSFastTrack > CTPPSFastTrackContainer
void setcellid(unsigned int cellid)
void ReadRecHits(edm::Handle< CTPPSFastRecHitContainer > &)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
def move(src, dest)
Definition: eostools.py:511
Definition: event.py:1
bool SearchTrack(int, int, int Direction, double &xi, double &t, double &partP, double &pt, double &thx, double &thy, double &x0, double &y0, double &xt, double &yt, double &X1d, double &Y1d, double &X2d, double &Y2d)