FP420/src/FP420SD.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: FP420SD.cc
0003 // Date: 02.2006
0004 // Description: Sensitive Detector class for FP420
0005 // Modifications:
0006 ///////////////////////////////////////////////////////////////////////////////
0007
0008 #include "SimG4Core/SensitiveDetector/interface/FrameRotation.h"
0009 #include "SimG4Core/Notification/interface/TrackInformation.h"
0010 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0011 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0012
0013 #include "SimDataFormats/SimHitMaker/interface/TrackingSlaveSD.h"
0014 #include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
0015 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0016
0017 #include "SimG4CMS/FP420/interface/FP420SD.h"
0018 #include "SimG4CMS/FP420/interface/FP420G4Hit.h"
0019 #include "SimG4CMS/FP420/interface/FP420G4HitCollection.h"
0020 #include "SimG4CMS/FP420/interface/FP420NumberingScheme.h"
0021
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025
0026 #include "SimG4Core/Notification/interface/TrackInformation.h"
0027 #include "SimG4Core/Notification/interface/SimTrackManager.h"
0028
0029 #include "G4Track.hh"
0030 #include "G4SDManager.hh"
0031 #include "G4VProcess.hh"
0032 #include "G4EventManager.hh"
0033 #include "G4Step.hh"
0034 #include "G4ParticleTable.hh"
0035
0036 #include <CLHEP/Units/SystemOfUnits.h>
0037
0038 #include <string>
0039 #include <vector>
0040 #include <iostream>
0041
0042 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0043
0044 //#define debug
0045 //-------------------------------------------------------------------
0046 FP420SD::FP420SD(const std::string& name,
0047                  const SensitiveDetectorCatalog& clg,
0048                  edm::ParameterSet const& p,
0049                  const SimTrackManager* manager)
0050     : SensitiveTkDetector(name, clg),
0051       numberingScheme(nullptr),
0052       hcID(-1),
0053       theHC(nullptr),
0054       theManager(manager),
0055       currentHit(nullptr),
0056       theTrack(nullptr),
0057       currentPV(nullptr),
0058       unitID(0),
0059       previousUnitID(0),
0060       preStepPoint(nullptr),
0061       postStepPoint(nullptr),
0062       eventno(0) {
0063   //Parameters
0064   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("FP420SD");
0065   int verbn = m_p.getUntrackedParameter<int>("Verbosity");
0066   //int verbn = 1;
0067
0068   SetVerboseLevel(verbn);
0069
0070   slave = new TrackingSlaveSD(name);
0071
0072   if (name == "FP420SI") {
0073     if (verbn > 0) {
0074       edm::LogInfo("FP420Sim") << "name = FP420SI and  new FP420NumberingSchem";
0075     }
0076     numberingScheme = new FP420NumberingScheme();
0077   } else {
0078     edm::LogWarning("FP420Sim") << "FP420SD: ReadoutName not supported\n";
0079   }
0080 }
0081
0082 FP420SD::~FP420SD() {
0083   delete slave;
0084   delete numberingScheme;
0085 }
0086
0087 double FP420SD::getEnergyDeposit(G4Step* aStep) { return aStep->GetTotalEnergyDeposit(); }
0088
0089 void FP420SD::Initialize(G4HCofThisEvent* HCE) {
0090 #ifdef debug
0091   LogDebug("FP420Sim") << "FP420SD : Initialize called for " << name << std::endl;
0092 #endif
0093
0094   theHC = new FP420G4HitCollection(GetName(), collectionName[0]);
0095   if (hcID < 0)
0096     hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0097   HCE->AddHitsCollection(hcID, theHC);
0098
0099   tsID = -2;
0100   //  primID = -2;
0101   primID = 0;
0102
0103   ////    slave->Initialize();
0104 }
0105
0106 bool FP420SD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0107   if (aStep == nullptr) {
0108     return true;
0109   } else {
0110     GetStepInfo(aStep);
0111     //   LogDebug("FP420Sim") << edeposit <<std::endl;
0112
0113     //AZ
0114 #ifdef debug
0115     LogDebug("FP420Sim") << "FP420SD :  number of hits = " << theHC->entries() << std::endl;
0116 #endif
0117
0118     if (HitExists() == false && edeposit > 0. && theHC->entries() < 200) {
0119       CreateNewHit();
0120       return true;
0121     }
0122   }
0123   return true;
0124 }
0125
0126 void FP420SD::GetStepInfo(G4Step* aStep) {
0127   preStepPoint = aStep->GetPreStepPoint();
0128   postStepPoint = aStep->GetPostStepPoint();
0129   theTrack = aStep->GetTrack();
0130   hitPoint = preStepPoint->GetPosition();
0131   currentPV = preStepPoint->GetPhysicalVolume();
0132   hitPointExit = postStepPoint->GetPosition();
0133
0134   hitPointLocal = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0135   hitPointLocalExit = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPointExit);
0136
0137   G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
0138   if (particleCode == emPDG || particleCode == epPDG || particleCode == gammaPDG) {
0139     edepositEM = getEnergyDeposit(aStep);
0140     edepositHAD = 0.;
0141   } else {
0142     edepositEM = 0.;
0143     edepositHAD = getEnergyDeposit(aStep);
0144   }
0145   edeposit = aStep->GetTotalEnergyDeposit();
0146   tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
0147   tSliceID = (int)tSlice;
0148   unitID = setDetUnitId(aStep);
0149 #ifdef debug
0150   LogDebug("FP420Sim") << "unitID=" << unitID << std::endl;
0151 #endif
0152   primaryID = theTrack->GetTrackID();
0153   //  Position     = hitPoint;
0154   Pabs = aStep->GetPreStepPoint()->GetMomentum().mag() / CLHEP::GeV;
0155   //Tof          = 1400. + aStep->GetPostStepPoint()->GetGlobalTime()/nanosecond;
0156   Tof = aStep->GetPostStepPoint()->GetGlobalTime() / CLHEP::nanosecond;
0157   Eloss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0158   ParticleType = theTrack->GetDefinition()->GetPDGEncoding();
0159   ThetaAtEntry = aStep->GetPreStepPoint()->GetPosition().theta() / CLHEP::deg;
0160   PhiAtEntry = aStep->GetPreStepPoint()->GetPosition().phi() / CLHEP::deg;
0161
0162   ParentId = theTrack->GetParentID();
0163   Vx = theTrack->GetVertexPosition().x();
0164   Vy = theTrack->GetVertexPosition().y();
0165   Vz = theTrack->GetVertexPosition().z();
0166   X = hitPoint.x();
0167   Y = hitPoint.y();
0168   Z = hitPoint.z();
0169 }
0170
0171 uint32_t FP420SD::setDetUnitId(const G4Step* aStep) {
0172   return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
0173 }
0174
0175 G4bool FP420SD::HitExists() {
0176   if (primaryID < 1) {
0177     edm::LogWarning("FP420Sim") << "***** FP420SD error: primaryID = " << primaryID << " maybe detector name changed";
0178   }
0179
0180   // Update if in the same detector, time-slice and for same track
0181   //  if (primaryID == primID && tSliceID == tsID && unitID==previousUnitID) {
0182   if (tSliceID == tsID && unitID == previousUnitID) {
0183     //AZ:
0184     UpdateHit();
0185     return true;
0186   }
0187   // Reset entry point for new primary
0188   if (primaryID != primID)
0189     ResetForNewPrimary();
0190
0191   //look in the HitContainer whether a hit with the same primID, unitID,
0192   //tSliceID already exists:
0193
0194   G4bool found = false;
0195
0196   //    LogDebug("FP420Sim") << "FP420SD: HCollection=  " << theHC->entries()    <<std::endl;
0197   int nhits = theHC->entries();
0198   for (int j = 0; j < nhits && !found; j++) {
0199     FP420G4Hit* aPreviousHit = (*theHC)[j];
0200     if (aPreviousHit->getTrackID() == primaryID && aPreviousHit->getTimeSliceID() == tSliceID &&
0201         aPreviousHit->getUnitID() == unitID) {
0202       //AZ:
0203       currentHit = aPreviousHit;
0204       found = true;
0205     }
0206   }
0207
0208   if (found) {
0209     //AZ:
0210     UpdateHit();
0211     return true;
0212   } else {
0213     return false;
0214   }
0215 }
0216
0217 void FP420SD::ResetForNewPrimary() {
0218   entrancePoint = SetToLocal(hitPoint);
0219   exitPoint = SetToLocalExit(hitPointExit);
0220   incidentEnergy = preStepPoint->GetKineticEnergy();
0221 }
0222
0223 void FP420SD::StoreHit(FP420G4Hit* hit) {
0224   //  if (primID<0) return;
0225   if (hit == nullptr) {
0226     edm::LogWarning("FP420Sim") << "FP420SD: hit to be stored is NULL !!";
0227     return;
0228   }
0229
0230   theHC->insert(hit);
0231 }
0232
0233 void FP420SD::CreateNewHit() {
0234 #ifdef debug
0235   //       << " MVid = " << currentPV->GetMother()->GetCopyNo()
0236   LogDebug("FP420Sim") << "FP420SD CreateNewHit for"
0237                        << " PV " << currentPV->GetName() << " PVid = " << currentPV->GetCopyNo() << " Unit " << unitID
0238                        << std::endl;
0239   LogDebug("FP420Sim") << " primary " << primaryID << " time slice " << tSliceID << " For Track  "
0240                        << theTrack->GetTrackID() << " which is a " << theTrack->GetDefinition()->GetParticleName();
0241
0242   if (theTrack->GetTrackID() == 1) {
0243     LogDebug("FP420Sim") << " of energy " << theTrack->GetTotalEnergy();
0244   } else {
0245     LogDebug("FP420Sim") << " daughter of part. " << theTrack->GetParentID();
0246   }
0247
0248   LogDebug("FP420Sim") << " and created by ";
0249   if (theTrack->GetCreatorProcess() != NULL)
0250     LogDebug("FP420Sim") << theTrack->GetCreatorProcess()->GetProcessName();
0251   else
0252     LogDebug("FP420Sim") << "NO process";
0253   LogDebug("FP420Sim") << std::endl;
0254 #endif
0255
0256   currentHit = new FP420G4Hit;
0257   currentHit->setTrackID(primaryID);
0258   currentHit->setTimeSlice(tSlice);
0259   currentHit->setUnitID(unitID);
0260   currentHit->setIncidentEnergy(incidentEnergy);
0261
0262   currentHit->setPabs(Pabs);
0263   currentHit->setTof(Tof);
0264   currentHit->setEnergyLoss(Eloss);
0265   currentHit->setParticleType(ParticleType);
0266   currentHit->setThetaAtEntry(ThetaAtEntry);
0267   currentHit->setPhiAtEntry(PhiAtEntry);
0268
0269   // currentHit->setEntry(entrancePoint);
0270   currentHit->setEntry(hitPoint);
0271
0272   currentHit->setEntryLocalP(hitPointLocal);
0273   currentHit->setExitLocalP(hitPointLocalExit);
0274
0275   currentHit->setParentId(ParentId);
0276   currentHit->setVx(Vx);
0277   currentHit->setVy(Vy);
0278   currentHit->setVz(Vz);
0279
0280   currentHit->setX(X);
0281   currentHit->setY(Y);
0282   currentHit->setZ(Z);
0283   //AZ:12.10.2007
0284   //  UpdateHit();
0285   // buffer for next steps:
0286   tsID = tSliceID;
0287   primID = primaryID;
0288   previousUnitID = unitID;
0289
0290   StoreHit(currentHit);
0291 }
0292
0293 void FP420SD::UpdateHit() {
0294   //
0295   if (Eloss > 0.) {
0296     currentHit->addEnergyDeposit(edepositEM, edepositHAD);
0297
0298 #ifdef debug
0299     LogDebug("FP420Sim") << "updateHit: add eloss " << Eloss << std::endl;
0300     LogDebug("FP420Sim") << "CurrentHit=" << currentHit << ", PostStepPoint=" << postStepPoint->GetPosition()
0301                          << std::endl;
0302 #endif
0303     //AZ
0304     //   currentHit->setEnergyLoss(Eloss);
0305     currentHit->addEnergyLoss(Eloss);
0306   }
0307
0308   // buffer for next steps:
0309   tsID = tSliceID;
0310   primID = primaryID;
0311   previousUnitID = unitID;
0312 }
0313
0314 G4ThreeVector FP420SD::SetToLocal(const G4ThreeVector& global) {
0315   const G4VTouchable* touch = preStepPoint->GetTouchable();
0316   theEntryPoint = touch->GetHistory()->GetTopTransform().TransformPoint(global);
0317   return theEntryPoint;
0318 }
0319
0320 G4ThreeVector FP420SD::SetToLocalExit(const G4ThreeVector& globalPoint) {
0321   const G4VTouchable* touch = postStepPoint->GetTouchable();
0322   theExitPoint = touch->GetHistory()->GetTopTransform().TransformPoint(globalPoint);
0323   return theExitPoint;
0324 }
0325
0326 void FP420SD::EndOfEvent(G4HCofThisEvent*) {
0327   // here we loop over transient hits and make them persistent
0328
0329   //  if(theHC->entries() > 100){
0330   //    LogDebug("FP420Sim") << "FP420SD: warning!!! Number of hits exceed 100 and =" << theHC->entries() << "\n";
0331   //  }
0332   //  for (int j=0; j<theHC->entries() && j<100; j++) {
0333   int nhitsHPS240 = 0;
0334   int nhitsFP420 = 0;
0335   int nhits = theHC->entries();
0336   for (int j = 0; j < nhits; j++) {
0337     FP420G4Hit* aHit = (*theHC)[j];
0338     if ((fabs(aHit->getTof()) > 780. && fabs(aHit->getTof()) < 840.))
0339       ++nhitsHPS240;
0340     if ((fabs(aHit->getTof()) > 1380. && fabs(aHit->getTof()) < 1450.))
0341       ++nhitsFP420;
0342     //    if(fabs(aHit->getTof()) < 1700.) {
0343     if ((fabs(aHit->getTof()) > 780. && fabs(aHit->getTof()) < 840. && nhitsHPS240 < 200.) ||
0344         (fabs(aHit->getTof()) > 1380. && fabs(aHit->getTof()) < 1450. && nhitsFP420 < 200.)) {
0345 #ifdef ddebug
0346       //    LogDebug("FP420SD") << " FP420Hit " << j << " " << *aHit << std::endl;
0347       LogDebug("FP420Sim") << "hit number" << j << "unit ID = " << aHit->getUnitID() << "\n";
0348       LogDebug("FP420Sim") << "entry z " << aHit->getEntry().z() << "\n";
0349       LogDebug("FP420Sim") << "entr theta " << aHit->getThetaAtEntry() << "\n";
0350 #endif
0351
0352       //    Local3DPoint locExitPoint(0,0,0);
0353       //  Local3DPoint locEntryPoint(aHit->getEntry().x(),
0354       //     aHit->getEntry().y(),
0355       //     aHit->getEntry().z());
0356       Local3DPoint locExitPoint(aHit->getExitLocalP().x(), aHit->getExitLocalP().y(), aHit->getExitLocalP().z());
0357       Local3DPoint locEntryPoint(aHit->getEntryLocalP().x(), aHit->getEntryLocalP().y(), aHit->getEntryLocalP().z());
0358       // implicit conversion (slicing) to PSimHit!!!
0359       slave->processHits(PSimHit(locEntryPoint,
0360                                  locExitPoint,             //entryPoint(), exitPoint()  Local3DPoint
0361                                  aHit->getPabs(),          // pabs()  float
0362                                  aHit->getTof(),           // tof() float
0363                                  aHit->getEnergyLoss(),    // energyLoss() float
0364                                  aHit->getParticleType(),  // particleType()   int
0365                                  aHit->getUnitID(),        // detUnitId() unsigned int
0366                                  aHit->getTrackID(),       // trackId() unsigned int
0367                                  aHit->getThetaAtEntry(),  //  thetaAtEntry()   float
0368                                  aHit->getPhiAtEntry()));  //  phiAtEntry()   float
0369
0370       //PSimHit( const Local3DPoint& entry, const Local3DPoint& exit,
0371       //       float pabs, float tof, float eloss, int particleType,
0372       //       unsigned int detId, unsigned int trackId,
0373       //       float theta, float phi, unsigned short processType=0) :
0374
0375       //  LocalVector direction = hit.exitPoint() - hit.entryPoint();
0376       //hit.energyLoss
0377
0378       /*
0379           aHit->getEM(),               -
0380           aHit->getHadr(),             -
0381           aHit->getIncidentEnergy(),   -
0382           aHit->getTimeSlice(),       -
0383           aHit->getEntry(),     -
0384           aHit->getParentId(),
0385           aHit->getEntryLocalP(),  -
0386           aHit->getExitLocalP(),   -
0387           aHit->getX(),    -
0388           aHit->getY(),   -
0389           aHit->getZ(),   -
0390           aHit->getVx(),  -
0391           aHit->getVy(),  -
0392           aHit->getVz()));  -
0393       */
0394     }  //if Tof<1600. if nhits<100
0395   }  // for loop on hits
0396
0397   Summarize();
0398 }
0399
0400 void FP420SD::Summarize() {}
0401
0402 void FP420SD::clear() {}
0403
0404 void FP420SD::DrawAll() {}
0405
0406 void FP420SD::PrintAll() {
0407   LogDebug("FP420Sim") << "FP420SD: Collection " << theHC->GetName() << "\n";
0408   theHC->PrintAllHits();
0409 }
0410
0411 //void FP420SD::SetNumberingScheme(FP420NumberingScheme* scheme){
0412 //
0413 //  if (numberingScheme)
0414 //    delete numberingScheme;
0415 //  numberingScheme = scheme;
0416 //
0417 //}
0418
0419 void FP420SD::fillHits(edm::PSimHitContainer& cc, const std::string& hname) {
0420   if (slave->name() == hname) {
0421     cc = slave->hits();
0422   }
0423 }
0424
0425 void FP420SD::update(const BeginOfEvent* i) {
0426   LogDebug("ForwardSim") << " Dispatched BeginOfEvent for " << GetName() << " !";
0427   clearHits();
0428   eventno = (*i)()->GetEventID();
0429 }
0430
0431 void FP420SD::update(const BeginOfRun*) {
0432   G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
0433   G4String particleName;
0434   emPDG = theParticleTable->FindParticle(particleName = "e-")->GetPDGEncoding();
0435   epPDG = theParticleTable->FindParticle(particleName = "e+")->GetPDGEncoding();
0436   gammaPDG = theParticleTable->FindParticle(particleName = "gamma")->GetPDGEncoding();
0437 }
0438
0439 void FP420SD::update(const ::EndOfEvent*) {}
0440
0441 void FP420SD::clearHits() {
0442   //AZ:
0443   slave->Initialize();
0444 }