Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:19 

0001 // -*- C++ -*-
0002 //
0003 // Package:     Forward
0004 // Class  :     TotemSD
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author:
0010 //         Created:  Tue May 16 10:14:34 CEST 2006
0011 //
0012 
0013 // system include files
0014 
0015 // user include files
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 
0019 #include "SimG4Core/Notification/interface/TrackInformation.h"
0020 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0021 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0022 
0023 #include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
0024 #include "SimDataFormats/SimHitMaker/interface/TrackingSlaveSD.h"
0025 
0026 #include "SimG4CMS/Forward/interface/TotemSD.h"
0027 #include "SimG4CMS/Forward/interface/TotemNumberMerger.h"
0028 #include "SimG4CMS/Forward/interface/TotemT1NumberingScheme.h"
0029 #include "SimG4CMS/Forward/interface/TotemT2NumberingSchemeGem.h"
0030 #include "SimG4CMS/Forward/interface/TotemRPNumberingScheme.h"
0031 
0032 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0033 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0034 
0035 #include "G4SDManager.hh"
0036 #include "G4Step.hh"
0037 #include "G4Track.hh"
0038 #include "G4VProcess.hh"
0039 
0040 #include "G4PhysicalConstants.hh"
0041 #include <CLHEP/Units/SystemOfUnits.h>
0042 
0043 //#define EDM_ML_DEBUG
0044 //
0045 // constructors and destructor
0046 //
0047 TotemSD::TotemSD(const std::string& name,
0048                  const SensitiveDetectorCatalog& clg,
0049                  edm::ParameterSet const& p,
0050                  const SimTrackManager* manager)
0051     : SensitiveTkDetector(name, clg),
0052       numberingScheme(nullptr),
0053       hcID(-1),
0054       theHC(nullptr),
0055       theManager(manager),
0056       currentHit(nullptr),
0057       theTrack(nullptr),
0058       currentPV(nullptr),
0059       unitID(0),
0060       previousUnitID(0),
0061       preStepPoint(nullptr),
0062       postStepPoint(nullptr) {
0063   //Parameters
0064   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("TotemSD");
0065   int verbn = m_p.getUntrackedParameter<int>("Verbosity");
0066 
0067   SetVerboseLevel(verbn);
0068 
0069   slave = new TrackingSlaveSD(name);
0070 
0071   if (name == "TotemHitsT1") {
0072     numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT1NumberingScheme(1));
0073   } else if (name == "TotemHitsT2Si") {
0074     numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT2NumberingSchemeGem(3));
0075   } else if (name == "TotemHitsT2Gem") {
0076     numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT2NumberingSchemeGem(4));
0077   } else if (name == "TotemHitsRP") {
0078     numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemRPNumberingScheme(3));
0079   } else {
0080     edm::LogWarning("ForwardSim") << "TotemSD: ReadoutName not supported\n";
0081   }
0082 }
0083 
0084 TotemSD::~TotemSD() {
0085   delete slave;
0086   delete numberingScheme;
0087 }
0088 
0089 bool TotemSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0090   getStepInfo(aStep);
0091   if (!hitExists() && edeposit > 0.) {
0092     createNewHit();
0093     return true;
0094   }
0095   if (!hitExists() && (((unitID == 1111 || unitID == 2222) && ParentId == 0 && ParticleType == 2212))) {
0096     createNewHitEvo();
0097   }
0098   return true;
0099 }
0100 
0101 uint32_t TotemSD::setDetUnitId(const G4Step* aStep) {
0102   return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
0103 }
0104 
0105 void TotemSD::Initialize(G4HCofThisEvent* HCE) {
0106 #ifdef EDM_ML_DEBUG
0107   edm::LogVerbatim("ForwardSim") << "TotemSD : Initialize called for " << GetName();
0108 #endif
0109 
0110   theHC = new TotemG4HitCollection(GetName(), collectionName[0]);
0111   if (hcID < 0)
0112     hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0113   HCE->AddHitsCollection(hcID, theHC);
0114 
0115   tsID = -2;
0116   primID = -2;
0117 }
0118 
0119 void TotemSD::EndOfEvent(G4HCofThisEvent*) {
0120   // here we loop over transient hits and make them persistent
0121   int thehc_entries = theHC->entries();
0122   for (int j = 0; j < thehc_entries && j < 15000; j++) {
0123     TotemG4Hit* aHit = (*theHC)[j];
0124 #ifdef EDM_ML_DEBUG
0125     edm::LogVerbatim("ForwardSim") << "HIT NUMERO " << j << "unit ID = " << aHit->getUnitID()
0126                                    << "\n               enrty z " << aHit->getEntry().z() << "\n               theta   "
0127                                    << aHit->getThetaAtEntry() << "\n";
0128 #endif
0129     Local3DPoint theExitPoint(0, 0, 0);
0130     Local3DPoint Entrata(aHit->getEntry().x(), aHit->getEntry().y(), aHit->getEntry().z());
0131     slave->processHits(PSimHit(Entrata,
0132                                theExitPoint,
0133                                aHit->getPabs(),
0134                                aHit->getTof(),
0135                                aHit->getEnergyLoss(),
0136                                aHit->getParticleType(),
0137                                aHit->getUnitID(),
0138                                aHit->getTrackID(),
0139                                aHit->getThetaAtEntry(),
0140                                aHit->getPhiAtEntry()));
0141   }
0142 }
0143 
0144 void TotemSD::PrintAll() {
0145 #ifdef EDM_ML_DEBUG
0146   edm::LogVerbatim("ForwardSim") << "TotemSD: Collection " << theHC->GetName();
0147 #endif
0148   theHC->PrintAllHits();
0149 }
0150 
0151 void TotemSD::fillHits(edm::PSimHitContainer& cc, const std::string& hname) {
0152   if (slave->name() == hname) {
0153     cc = slave->hits();
0154   }
0155 }
0156 
0157 void TotemSD::update(const BeginOfEvent* i) {
0158 #ifdef EDM_ML_DEBUG
0159   edm::LogVerbatim("ForwardSim") << " Dispatched BeginOfEvent for " << GetName() << " !";
0160 #endif
0161   clearHits();
0162 }
0163 
0164 void TotemSD::clearHits() { slave->Initialize(); }
0165 
0166 G4ThreeVector TotemSD::setToLocal(const G4ThreeVector& global) {
0167   G4ThreeVector localPoint;
0168   const G4VTouchable* touch = preStepPoint->GetTouchable();
0169   localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(global);
0170   return localPoint;
0171 }
0172 
0173 void TotemSD::getStepInfo(const G4Step* aStep) {
0174   preStepPoint = aStep->GetPreStepPoint();
0175   postStepPoint = aStep->GetPostStepPoint();
0176   theTrack = aStep->GetTrack();
0177   hitPoint = preStepPoint->GetPosition();
0178   currentPV = preStepPoint->GetPhysicalVolume();
0179 
0180   // double weight = 1;
0181   G4String name = currentPV->GetName();
0182   name.assign(name, 0, 4);
0183   G4String particleType = theTrack->GetDefinition()->GetParticleName();
0184   edeposit = aStep->GetTotalEnergyDeposit();
0185 
0186   tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
0187   tSliceID = (int)tSlice;
0188   unitID = setDetUnitId(aStep);
0189 #ifdef EDM_ML_DEBUG
0190   edm::LogVerbatim("ForwardSim") << "UNITa " << unitID;
0191 #endif
0192   primaryID = theTrack->GetTrackID();
0193 
0194   Posizio = hitPoint;
0195   Pabs = aStep->GetPreStepPoint()->GetMomentum().mag() / CLHEP::GeV;
0196   Tof = aStep->GetPostStepPoint()->GetGlobalTime() / CLHEP::nanosecond;
0197 
0198   Eloss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0199   ParticleType = theTrack->GetDefinition()->GetPDGEncoding();
0200 
0201   ThetaAtEntry = aStep->GetPreStepPoint()->GetPosition().theta() / CLHEP::deg;
0202   PhiAtEntry = aStep->GetPreStepPoint()->GetPosition().phi() / CLHEP::deg;
0203 
0204   ParentId = theTrack->GetParentID();
0205   Vx = theTrack->GetVertexPosition().x();
0206   Vy = theTrack->GetVertexPosition().y();
0207   Vz = theTrack->GetVertexPosition().z();
0208 }
0209 
0210 bool TotemSD::hitExists() {
0211   if (primaryID < 1) {
0212     edm::LogWarning("ForwardSim") << "***** TotemSD error: primaryID = " << primaryID << " maybe detector name changed";
0213   }
0214 
0215   // Update if in the same detector, time-slice and for same track
0216   //  if (primaryID == primID && tSliceID == tsID && unitID==previousUnitID) {
0217   if (tSliceID == tsID && unitID == previousUnitID) {
0218     updateHit();
0219     return true;
0220   }
0221 
0222   // Reset entry point for new primary
0223   if (primaryID != primID)
0224     resetForNewPrimary();
0225 
0226   //look in the HitContainer whether a hit with the same primID, unitID,
0227   //tSliceID already exists:
0228 
0229   bool found = false;
0230   int thehc_entries = theHC->entries();
0231   for (int j = 0; j < thehc_entries && !found; j++) {
0232     TotemG4Hit* aPreviousHit = (*theHC)[j];
0233     if (aPreviousHit->getTrackID() == primaryID && aPreviousHit->getTimeSliceID() == tSliceID &&
0234         aPreviousHit->getUnitID() == unitID) {
0235       currentHit = aPreviousHit;
0236       found = true;
0237       break;
0238     }
0239   }
0240 
0241   if (found) {
0242     updateHit();
0243     return true;
0244   } else {
0245     return false;
0246   }
0247 }
0248 
0249 void TotemSD::createNewHit() {
0250 #ifdef EDM_ML_DEBUG
0251   edm::LogVerbatim("ForwardSim") << "TotemSD CreateNewHit for PV " << currentPV->GetName()
0252                                  << " PVid = " << currentPV->GetCopyNo() << " Unit " << unitID << "\n primary "
0253                                  << primaryID << " time slice " << tSliceID << " For Track  " << theTrack->GetTrackID()
0254                                  << " which is a " << theTrack->GetDefinition()->GetParticleName();
0255 
0256   if (theTrack->GetTrackID() == 1) {
0257     edm::LogVerbatim("ForwardSim") << " of energy " << theTrack->GetTotalEnergy();
0258   } else {
0259     edm::LogVerbatim("ForwardSim") << " daughter of part. " << theTrack->GetParentID();
0260   }
0261 
0262   if (theTrack->GetCreatorProcess() != nullptr)
0263     edm::LogVerbatim("ForwardSim") << theTrack->GetCreatorProcess()->GetProcessName();
0264   else
0265     edm::LogVerbatim("ForwardSim") << "NO process";
0266 #endif
0267 
0268   currentHit = new TotemG4Hit;
0269   currentHit->setTrackID(primaryID);
0270   currentHit->setTimeSlice(tSlice);
0271   currentHit->setUnitID(unitID);
0272   currentHit->setIncidentEnergy(incidentEnergy);
0273 
0274   currentHit->setPabs(Pabs);
0275   currentHit->setTof(Tof);
0276   currentHit->setEnergyLoss(Eloss);
0277   currentHit->setParticleType(ParticleType);
0278   currentHit->setThetaAtEntry(ThetaAtEntry);
0279   currentHit->setPhiAtEntry(PhiAtEntry);
0280 
0281   currentHit->setEntry(Posizio.x(), Posizio.y(), Posizio.z());
0282 
0283   currentHit->setParentId(ParentId);
0284   currentHit->setVx(Vx);
0285   currentHit->setVy(Vy);
0286   currentHit->setVz(Vz);
0287 
0288   updateHit();
0289 
0290   storeHit(currentHit);
0291 }
0292 
0293 void TotemSD::createNewHitEvo() {
0294   // edm::LogVerbatim("ForwardSim") << "INSIDE CREATE NEW HIT EVO ";
0295 
0296   currentHit = new TotemG4Hit;
0297   currentHit->setTrackID(primaryID);
0298   currentHit->setTimeSlice(tSlice);
0299   currentHit->setUnitID(unitID);
0300   currentHit->setIncidentEnergy(incidentEnergy);
0301 
0302   currentHit->setPabs(Pabs);
0303   currentHit->setTof(Tof);
0304   currentHit->setEnergyLoss(Eloss);
0305   currentHit->setParticleType(ParticleType);
0306   currentHit->setThetaAtEntry(ThetaAtEntry);
0307   currentHit->setPhiAtEntry(PhiAtEntry);
0308 
0309   //  edm::LogVerbatim("ForwardSim") << Posizio.x() << " " << Posizio.y() << " " << Posizio.z();
0310 
0311   currentHit->setParentId(ParentId);
0312   currentHit->setVx(Vx);
0313   currentHit->setVy(Vy);
0314   currentHit->setVz(Vz);
0315 
0316   G4ThreeVector _PosizioEvo;
0317   int flagAcc = 0;
0318   _PosizioEvo = posizioEvo(Posizio, Vx, Vy, Vz, Pabs, flagAcc);
0319 
0320   if (flagAcc == 1) {
0321     currentHit->setEntry(_PosizioEvo.x(), _PosizioEvo.y(), _PosizioEvo.z());
0322 
0323     updateHit();
0324 
0325     storeHit(currentHit);
0326   }
0327   // edm::LogVerbatim("ForwardSim") << "STORED HIT IN: " << unitID;
0328 }
0329 
0330 G4ThreeVector TotemSD::posizioEvo(
0331     const G4ThreeVector& Pos, double vx, double vy, double vz, double pabs, int& accettanza) {
0332   accettanza = 0;
0333   //Pos.xyz() in mm
0334   G4ThreeVector PosEvo;
0335   double ThetaX = std::atan((Pos.x() - vx) / (Pos.z() - vz));
0336   double ThetaY = std::atan((Pos.y() - vy) / (Pos.z() - vz));
0337   double X_at_0 = (vx - ((Pos.x() - vx) / (Pos.z() - vz)) * vz) / 1000.;
0338   double Y_at_0 = (vy - ((Pos.y() - vy) / (Pos.z() - vz)) * vz) / 1000.;
0339 
0340   //  double temp_evo_X;
0341   //  double temp_evo_Y;
0342   //  double temp_evo_Z;
0343   //  temp_evo_Z = std::abs(Pos.z())/Pos.z()*220000.;
0344 
0345   //csi=-dp/d
0346   double csi = std::abs((7000. - pabs) / 7000.);
0347 
0348   // all in m
0349   const int no_rp = 4;
0350   double x_par[no_rp + 1];
0351   double y_par[no_rp + 1];
0352   //rp z position
0353   double rp[no_rp] = {141., 149., 198., 220.};
0354   //{lx0,mlx} for each rp; Lx=lx0+mlx*csi
0355   double leffx[][2] = {{122.5429, -46.9312}, {125.4194, -49.1849}, {152.6, -81.157}, {98.8914, -131.8390}};
0356   //{ly0,mly} for each rp; Ly=ly0+mly*csi
0357   double leffy[][2] = {{124.2314, -55.4852}, {127.7825, -57.4503}, {179.455, -76.274}, {273.0931, -40.4626}};
0358   //{vx0,mvx0} for each rp; vx=vx0+mvx*csi
0359   double avx[][2] = {{0.515483, -1.0123}, {0.494122, -1.0534}, {0.2217, -1.483}, {0.004633, -1.0719}};
0360   //{vy0,mvy0} for each rp; vy=vy0+mvy*csi
0361   double avy[][2] = {{0.371418, -1.6327}, {0.349035, -1.6955}, {0.0815, -2.59}, {0.007592, -4.0841}};
0362   //{D0,md,a,b} for each rp; D=D0+(md+a*thetax)*csi+b*thetax
0363   double ddx[][4] = {{-0.082336, -0.092513, 112.3436, -82.5029},
0364                      {-0.086927, -0.097670, 114.9513, -82.9835},
0365                      {-0.092117, -0.0915, 180.6236, -82.443},
0366                      {-0.050470, 0.058837, 208.1106, 20.8198}};
0367   // {10sigma_x+0.5mm,10sigma_y+0.5mm}
0368   double detlim[][2] = {{0, 0}, {0.0028, 0.0021}, {0, 0}, {0.0008, 0.0013}};
0369   //{rmax,dmax}
0370   double pipelim[][2] = {{0.026, 0.026}, {0.04, 0.04}, {0.0226, 0.0177}, {0.04, 0.04}};
0371 
0372   for (int j = 0; j < no_rp; j++) {
0373     //y=Ly*thetay+vy*y0
0374     //x=Lx*thetax+vx*x0-csi*D
0375     y_par[j] = ThetaY * (leffy[j][0] + leffy[j][1] * csi) + (avy[j][0] + avy[j][1] * csi) * Y_at_0;
0376     x_par[j] = ThetaX * (leffx[j][0] + leffx[j][1] * csi) + (avx[j][0] + avx[j][1] * csi) * X_at_0 -
0377                csi * (ddx[j][0] + (ddx[j][1] + ddx[j][2] * ThetaX) * csi + ddx[j][3] * ThetaX);
0378   }
0379 
0380   //pass TAN@141
0381   if (std::abs(y_par[0]) < pipelim[0][1] && sqrt((y_par[0] * y_par[0]) + (x_par[0] * x_par[0])) < pipelim[0][0]) {
0382     //pass 149
0383     if ((sqrt((y_par[1] * y_par[1]) + (x_par[1] * x_par[1])) < pipelim[1][0]) &&
0384         (std::abs(y_par[1]) > detlim[1][1] || x_par[1] > detlim[1][0])) {
0385       accettanza = 1;
0386     }
0387   }
0388 
0389   //pass TAN@141
0390   if (std::abs(y_par[0]) < pipelim[0][1] && sqrt((y_par[0]) * (y_par[0]) + (x_par[0]) * (x_par[0])) < pipelim[0][0]) {
0391     //pass Q5@198
0392     if (std::abs(y_par[2]) < pipelim[2][1] && sqrt((y_par[2] * y_par[2]) + (x_par[2] * x_par[2])) < pipelim[2][0]) {
0393       //pass 220
0394       if ((sqrt((y_par[3] * y_par[3]) + (x_par[3] * x_par[3])) < pipelim[3][0]) &&
0395           (std::abs(y_par[3]) > detlim[3][1] || x_par[3] > detlim[3][0])) {
0396         accettanza = 1;
0397 
0398         PosEvo.setX(1000 * x_par[3]);
0399         PosEvo.setY(1000 * y_par[3]);
0400         PosEvo.setZ(1000 * rp[3]);
0401         if (Pos.z() < vz)
0402           PosEvo.setZ(-1000 * rp[3]);
0403       }
0404     }
0405   }
0406   /*
0407   edm::LogVerbatim("ForwardSim") << "\n"
0408              << "ACCETTANZA: "<<accettanza << "\n" 
0409              << "CSI: "<< csi << "\n"
0410              << "Theta_X: " << ThetaX << "\n"
0411              << "Theta_Y: " << ThetaY << "\n"
0412              << "X_at_0: "<< X_at_0 << "\n"
0413              << "Y_at_0: "<< Y_at_0 << "\n" 
0414              << "x_par[3]: "<< x_par[3] << "\n"
0415              << "y_par[3]: "<< y_par[3] << "\n"
0416              << "pos " << Pos.x() << " " << Pos.y() << " " 
0417              << Pos.z() << "\n" << "V "<< vx << " " << vy << " "
0418              << vz << "\n"
0419 */
0420   // --------------
0421   return PosEvo;
0422 }
0423 
0424 void TotemSD::updateHit() {
0425   //
0426   if (Eloss > 0.) {
0427 #ifdef EDM_ML_DEBUG
0428     edm::LogVerbatim("ForwardSim") << "G4TotemT1SD updateHit: add eloss " << Eloss << "\nCurrentHit=" << currentHit
0429                                    << ", PostStepPoint=" << postStepPoint->GetPosition();
0430 #endif
0431 
0432     currentHit->setEnergyLoss(Eloss);
0433   }
0434 
0435   // buffer for next steps:
0436   tsID = tSliceID;
0437   primID = primaryID;
0438   previousUnitID = unitID;
0439 }
0440 
0441 void TotemSD::storeHit(TotemG4Hit* hit) {
0442   if (primID < 0)
0443     return;
0444   if (hit == nullptr) {
0445     edm::LogWarning("ForwardSim") << "TotemSD: hit to be stored is NULL !!";
0446     return;
0447   }
0448 
0449   theHC->insert(hit);
0450 }
0451 
0452 void TotemSD::resetForNewPrimary() {
0453   entrancePoint = setToLocal(hitPoint);
0454   incidentEnergy = preStepPoint->GetKineticEnergy();
0455 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!