LBNEPerfectFocusProcess.cc

Go to the documentation of this file.

//
// For perfect focusing.  The OO way.  (Will probably end-up 
// increase tracking time for everybody..Hoefully not too much ..) 
// 
#include "LBNEPerfectFocusProcess.hh"
#include "LBNETrackInformation.hh"
#include "G4ios.hh"
#include "LBNEDetectorConstruction.hh"
#include "LBNERunManager.hh"

LBNEPerfectFocusParticleChange::LBNEPerfectFocusParticleChange() {
 ;
}
LBNEPerfectFocusParticleChange::~LBNEPerfectFocusParticleChange() {
;
}
LBNEPerfectFocusProcess::LBNEPerfectFocusProcess(const G4String &aName, G4ProcessType   aType):
G4VDiscreteProcess(aName, aType)
 {
  enablePostStepDoIt = true;
  const G4RunManager *pRunManager = G4RunManager::GetRunManager();
  const LBNEDetectorConstruction *LBNEDetCr = reinterpret_cast<const LBNEDetectorConstruction*>(pRunManager->GetUserDetectorConstruction());
  fZCoordForPerfectFocusing = LBNEDetCr->GetZCoordForPerfectFocusing();
  const double rr = LBNEDetCr->GetRCoordOutOfTarget();
  fRCoordOutOfTargetSq = rr*rr;
  fParticleChange = 0;
//   std::cerr << "  LBNEPerfectFocusProcess::LBNEPerfectFocusProcess, got here!  fZCoordForPerfectFocusing = " 
//             << fZCoordForPerfectFocusing << " And quit " << std::endl; exit(2);
}
LBNEPerfectFocusProcess::~LBNEPerfectFocusProcess() {; }
//
// Define the change of direction as Post step, if not already done for this track. 
// The new particle will be considered as a secondary, preserving the particle ancestry for the dk2nu Ntuple. 
//

LBNEPerfectFocusParticleChange* LBNEPerfectFocusProcess::PostStepDoIt (const G4Track &aTrack, const G4Step &stepData) {


   delete fParticleChange; // start by deleting the previous instance of the particle change.. 
   fParticleChange = 0;
   const G4DynamicParticle *pDef = aTrack.GetDynamicParticle();
   const int pId = pDef->GetDefinition()->GetPDGEncoding();
   if ((std::abs(pId) != 211) && (std::abs(pId) !=321)) return 0; // Should not occur
   const double totalP = pDef->GetTotalMomentum(); 
   const LBNETrackInformation *oldTrInfo = reinterpret_cast<const LBNETrackInformation*>(aTrack.GetUserInformation());
   if (oldTrInfo->GetPFFlag() == 1) return 0;
   
   // Now stack it's replacement. Define it as a secondary. 
   //
   LBNETrackInformation *newTrInfo = new LBNETrackInformation();
   newTrInfo->SetDecayCode(oldTrInfo->GetDecayCode());
   newTrInfo->SetTgen(oldTrInfo->GetTgen());
   newTrInfo->SetNImpWt(oldTrInfo->GetNImpWt());
   G4ThreeVector aNewPos(stepData.GetPostStepPoint()->GetPosition());
   //
   // Displace the track radially and in Z to avoid having it going through the perfect focusing again. 
   //
   aNewPos[2] += 1.*CLHEP::mm;
   const double phi = std::atan2(aNewPos[1], aNewPos[0]);
   aNewPos[0] += std::cos(phi)*1.*CLHEP::mm;
   aNewPos[1] += std::sin(phi)*1.*CLHEP::mm;
   newTrInfo->SetPFFlag(1);
   newTrInfo->SetDecayCode(oldTrInfo->GetDecayCode());
   newTrInfo->SetTgen(oldTrInfo->GetTgen());
   newTrInfo->SetNImpWt(oldTrInfo->GetNImpWt());
   //
   G4ThreeVector aNewP(0., 0., totalP);
   G4DynamicParticle *aNewPart =  new G4DynamicParticle(pDef->GetDefinition(), aNewP);
   G4Track *aNewTrack = new G4Track(aNewPart, aTrack.GetLocalTime(), aNewPos);
//   aNewTrack->SetUserInformation(newTrInfo); // It will get lost, G4 tracking Kernel does not bother to copy the User information. 
// So, we coCLHEP::mment this out... Sad... 
   delete newTrInfo;
   // 
   fParticleChange = new LBNEPerfectFocusParticleChange();
   fParticleChange->AddSecondary(aNewTrack);
   fParticleChange->ProposeTrackStatus( fStopAndKill );   
   // Check the flag.. 
   // 
//   std::cerr << " LBNEPerfectFocusProcess::PostStepDoIt .. Total P " << totalP 
//           << " Position " << aNewPos << " old Tr ID  " << aTrack.GetTrackID() 
//         << " New Id " << aNewTrack->GetTrackID() 
//         << " zAttrigger " << fZAtTrigger << " rAttrigger " << std::sqrt(fRAtTriggerSq) << std::endl;
//   const LBNETrackInformation *newTrInfoCheck = reinterpret_cast<const LBNETrackInformation*>(aNewTrack->GetUserInformation());
//   std::cerr << " Pointer for Tr User Info " <<  newTrInfoCheck << " New track PFFlag " << newTrInfoCheck->GetPFFlag() << std::endl;
   return fParticleChange;

}
 
//
// Not actually called 
//
G4double LBNEPerfectFocusProcess::PostStepGetPhysicalInteractionLength (const G4Track &aTrack, G4double previousStepSize, 
                                                   G4ForceCondition *condition)
{
//  std::cerr << "  LBNEPerfectFocusProcess::GetPhysicalInteractionLength, at " << aTrack.GetPosition() << std::endl;
  *condition = NotForced;
  return this->GetFreePath(aTrack, previousStepSize);
}

G4double LBNEPerfectFocusProcess::GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition* aCond) {
//  std::cerr << "  LBNEPerfectFocusProcess::GetMeanFreePath, at " << aTrack.GetPosition() << std::endl;
  *aCond = NotForced;
  return this->GetFreePath(aTrack, previousStepSize);
}



G4double LBNEPerfectFocusProcess::GetFreePath(const G4Track &aTrack, G4double previousStepSize) {
  
   // In this context, we have no use for this 2nd parameter in the argument call. 
   // However, one can not change the interface, G4VDiscreteProcess OO design.. We use to flag a crazy 
   // value to avoid a compiler warning. (and burn some CPU cyvle, but it's good for the Grid economy. 
   // Paul Lebrun, Feb. 2015 
  
   if (previousStepSize < -999999) std::cerr << " LBNEPerfectFocusProcess::GetFreePath, crazy prev step size.. " << std::endl;
//   const LBNETrackInformation *oldTrInfo = reinterpret_cast<const LBNETrackInformation*>(aTrack.GetUserInformation());
//   if (previousStepSize < 1.0e-4) return previousStepSize;
//   std::cerr << "  LBNEPerfectFocusProcess::GetMeanFreePath, got here! flag = " << oldTrInfo->GetPFFlag() << std::endl;
//   if (oldTrInfo->GetPFFlag() == 1) return 1.0e103; // we are done,  the change of direction already occured. DOes not currently work!!!
   if (fZCoordForPerfectFocusing < 0. ) {
        const G4ThreeVector posT = aTrack.GetPosition(); 
        const double radSq= posT[0]*posT[0] + posT[1]*posT[1];
        if ((std::abs(radSq - fRCoordOutOfTargetSq) < 0.5*CLHEP::mm) && (posT[2] < 5000.*CLHEP::mm)) {
         // it is unlikely that we will ever have a target longer than 5 meters 
          fRAtTriggerSq = radSq;
          fZAtTrigger = posT[2];
//          std::cerr << "  LBNEPerfectFocusProcess::GetMeanFreePath, at crit. radius " << std::sqrt(rAtTriggerSq) << std::endl;
          return 1.0e-5*CLHEP::mm; 
        } else {
          return 1.0e103;
        }                                               
   } else {
     if(aTrack.GetPosition()[2] > fZCoordForPerfectFocusing && (fZStart < fZCoordForPerfectFocusing)) {
      // If the track has been created downstream of fZCoordForPerfectFocusing, via an other process, mean free path inf. 
        fZAtTrigger = aTrack.GetPosition()[2];
//        std::cerr << "  LBNEPerfectFocusProcess::GetMeanFreePath, At Z " << zAtTrigger << std::endl;
        return 0.00001; 
     } else {
       return 1.0e103;
     }
   }
}

G4bool LBNEPerfectFocusProcess::IsApplicable(const G4ParticleDefinition  &pDef) {
  int pid = pDef.GetPDGEncoding();
  if (std::abs(pid) == 211) return true;
  if (std::abs(pid) == 321) return true;
  return false;
}

void LBNEPerfectFocusProcess::StartTracking(G4Track *aTrack) {
//   const LBNETrackInformation *oldTrInfo = reinterpret_cast<const LBNETrackInformation*>(aTrack->GetUserInformation());
//   std::cerr << "  LBNEPerfectFocusProcess::StartTracking, id " << aTrack->GetTrackID() 
//             << " at  " << aTrack->GetPosition() << " fRCoordOutOfTargetSq " << fRCoordOutOfTargetSq << std::endl;
   fZAtTrigger = -1.0e103;
   const G4ThreeVector posT = aTrack->GetPosition(); 
   fZStart = posT[2];
   fRAtTriggerSq = posT[0]*posT[0] + posT[1]*posT[1];
             
//   exit(2);
}