musrDecay.cc

Go to the documentation of this file.

#include"musrDecay.hh"

#include "G4DynamicParticle.hh"
#include "G4DecayProducts.hh"
#include "G4DecayTable.hh"
#include "G4PhysicsLogVector.hh"
#include "G4ParticleChangeForDecay.hh"
#include "G4VExtDecayer.hh"
#include "G4ThreeVector.hh"
#include "musrMuonDecayChannel.hh"
#include <iomanip>
#include "G4ios.hh"
#include "globals.hh"

G4VParticleChange* musrDecay::DecayIt(const G4Track& aTrack, const G4Step& aStep)
{




  //  G4cerr<<"METHOD CALLED WHEN SHOULD NOT BE AT ALL!!!!\n";//enable for geant code
  // The DecayIt() method returns by pointer a particle-change object.
  // Units are expressed in GEANT4 internal units.
  // get particle 
  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
  G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();


  // if particle is not a muon
  if(aParticle->GetDefinition()->GetParticleName()!="mu+")
        {
          G4ParticleChangeForDecay* k;
                    k = (G4ParticleChangeForDecay*)G4Decay::DecayIt(aTrack,aStep);
           return k ;
        }
  // m_ExtDecayer = 0;//G4Decay::GetExtDecayer();

  G4double p_time = aTrack.GetDynamicParticle()->GetPreAssignedDecayProperTime();
  G4double tp_time = aTrack.GetProperTime();
  m_RemainderLifeTime =p_time-tp_time;
  
  // check if  the particle is stable
  if (aParticleDef->GetPDGStable()) return &pParticleChangeForDecay ;
 

  //check if thePreAssignedDecayProducts exists
  const G4DecayProducts* o_products = (aParticle->GetPreAssignedDecayProducts());
  G4bool isPreAssigned = (o_products != 0);   
  G4DecayProducts* products = 0;

  // decay table
  G4DecayTable   *decaytable = aParticleDef->GetDecayTable();
 
  // check if external decayer exists
  G4bool isExtDecayer = (decaytable == 0);

  // Error due to NO Decay Table 
  if ( (decaytable == 0) && !isExtDecayer &&!isPreAssigned ){
#ifdef G4VERBOSE
    if (GetVerboseLevel()>0) {
      G4cerr <<  "G4Decay::DoIt  : decay table not defined  for ";
      G4cerr << aParticle->GetDefinition()->GetParticleName()<< G4endl;
    }
#endif
    pParticleChangeForDecay.SetNumberOfSecondaries(0);
    // Kill the parent particle
    pParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
    pParticleChangeForDecay.ProposeLocalEnergyDeposit(0.0); 
    
    ClearNumberOfInteractionLengthLeft();
    return &pParticleChangeForDecay ;
  }

  if (isPreAssigned) {
    // copy decay products 
    products = new G4DecayProducts(*o_products); 
  }  else {
 

    //+++++++++++++++++++++++++++++++++++ACCORDING TO DECAY TABLE++++++++++++++++++++++++++
   // decay acoording to decay table
    // choose a decay channel
    G4VDecayChannel *decaychannel = decaytable->SelectADecayChannel();
    if (decaychannel == 0 ){
      // decay channel not found
      G4Exception("G4Decay::DoIt  : can not determine decay channel ");
    } else {
      G4int temp = decaychannel->GetVerboseLevel();
      // execute DecayIt() 
#ifdef G4VERBOSE
      if (GetVerboseLevel()>1) {
        G4cerr << "G4Decay::DoIt  : selected decay channel  addr:" << decaychannel <<G4endl;
        decaychannel->SetVerboseLevel(GetVerboseLevel());
      }
#endif

      //  CHECK FOR POLARIASATION AND ASSIGN IT TO THE DECAY METHOD
      
      G4ThreeVector parent_polarization = aParticle->GetPolarization();
      
      //-----------------------------INFORMATION---------------------------------
      //         G4cout <<"LEMuSRDecay MESSAGE:: polarization is " << parent_polarization <<" .\n";
      if(aParticle->GetDefinition()->GetParticleName()=="mu+")
        {
          musrMuonDecayChannel *muondecaychannel =new musrMuonDecayChannel("mu+",1.00);
          if(parent_polarization == G4ThreeVector(0,0,0))
            {  
              products = muondecaychannel->DecayIt(aParticle->GetMass());
            }
          if(parent_polarization != G4ThreeVector(0,0,0))
            {  
              products = muondecaychannel->DecayItPolarized(aParticle->GetMass(),parent_polarization);// has to be included in G4VDecayChannel.
            }
        }
      
      else products = decaychannel->DecayIt(aParticle->GetMass()) ;
      
      
      
#ifdef G4VERBOSE
      if (GetVerboseLevel()>1) {
        decaychannel->SetVerboseLevel(temp);
      }
#endif
#ifdef G4VERBOSE
      // for debug
      //if (! products->IsChecked() ) products->DumpInfo();
#endif
    }
  }

  // get parent particle information ...................................
  G4double   ParentEnergy  = aParticle->GetTotalEnergy();
  G4ThreeVector ParentDirection(aParticle->GetMomentumDirection());

  //boost all decay products to laboratory frame
  G4double energyDeposit = 0.0;
  G4double finalGlobalTime = aTrack.GetGlobalTime();
  if (aTrack.GetTrackStatus() == fStopButAlive ){
    // AtRest case
    finalGlobalTime += m_RemainderLifeTime;
    energyDeposit += aParticle->GetKineticEnergy();
    if (isPreAssigned) products->Boost( ParentEnergy, ParentDirection);
  } else {
    // PostStep case
    if (!isExtDecayer) products->Boost( ParentEnergy, ParentDirection);
  }
  //add products in pParticleChangeForDecay
  G4int numberOfSecondaries = products->entries();

  pParticleChangeForDecay.SetNumberOfSecondaries(numberOfSecondaries);

  //      G4cerr << "G4Decay::DoIt  : Decay vertex :";
  //    G4cerr << " Time: " << finalGlobalTime/ns << "[ns]";
  //    G4cerr << " X:" << (aTrack.GetPosition()).x() /cm << "[cm]";
  //    G4cerr << " Y:" << (aTrack.GetPosition()).y() /cm << "[cm]";
  //    G4cerr << " Z:" << (aTrack.GetPosition()).z() /cm << "[cm]";
  //    G4cerr << G4endl;
    //    products->DumpInfo();

  G4int index;
  G4ThreeVector currentPosition;
  const G4TouchableHandle thand = aTrack.GetTouchableHandle();

  for (index=0; index < numberOfSecondaries; index++)
  {
     // get current position of the track
     currentPosition = aTrack.GetPosition();
     // create a new track object
     G4Track* secondary = new G4Track( products->PopProducts(),
                                      finalGlobalTime ,
                                      currentPosition );
     // switch on good for tracking flag
     secondary->SetGoodForTrackingFlag();
     secondary->SetTouchableHandle(thand);
     // add the secondary track in the List
     pParticleChangeForDecay.AddSecondary(secondary);
  }
  delete products;

  // Kill the parent particle
  pParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
  pParticleChangeForDecay.ProposeLocalEnergyDeposit(energyDeposit); 
  pParticleChangeForDecay.ProposeGlobalTime( finalGlobalTime );
  // reset NumberOfInteractionLengthLeft
  ClearNumberOfInteractionLengthLeft();

  return &pParticleChangeForDecay ;     

} 

---