DecayConfig.C

Go to the documentation of this file.

/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *              GNU Lesser General Public Licence (LGPL) version 3,             *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
void DecayConfig()
{

    // This script uses the external decayer TPythia6Decayer in place of the
    // concrete Monte Carlo native decay mechanisms only for the
    // specific types of decays defined below.

    // Access the external decayer singleton and initialize it
    TPythia6Decayer* decayer = TPythia6Decayer::Instance();
    // The following just tells pythia6 to not decay particles only to
    // certain channels.

    decayer->SetForceDecay(TPythia6Decayer::kAll);
    // example:  Force the J/PSI decay channel e+e-
    //        Int_t products[2];
    //        Int_t mult[2];
    //        Int_t npart=2;

    // decay products
    //        products[0]=11;
    //        products[1]=-11;
    // multiplicity
    //        mult[0]=1;
    //        mult[1]=1;
    // force the decay channel
    //        decayer->ForceParticleDecay(443,products,mult,npart);

    decayer->Init();

    // Tell the concrete monte carlo to use the external decayer.  The
    // external decayer will be used for:
    // i)particle decays not defined in concrete monte carlo, or
    // ii)particles for which the concrete monte carlo is told
    //   to use the external decayer for its type via:
    //     TVirtualMC::GetMC()->SetUserDecay(pdgId);
    //   If this is invoked, the external decayer will be used for particles
    //   of type pdgId even if the concrete monte carlo has a decay mode
    //   already defined for that particle type.
    TVirtualMC::GetMC()->SetExternalDecayer(decayer);

    TPythia6& pythia6 = *(TPythia6::Instance());

    // The pythia6 decayer is used in place of the concrete Monte Carlo
    // decay for the particles type mu+/-,pi+/-, K+/-, K0L in order to preserve
    // the decay product neutrino flavor, which is otherwise not preserved in
    // Geant3 decays.
    const Int_t npartnf = 9;
    // mu-,mu+,pi+,pi-,K+,K-,K0L, Xi-
    Int_t pdgnf[npartnf] = {13, -13, 211, -211, 321, -321, 130, 3312, 443};
    for (Int_t ipartnf = 0; ipartnf < npartnf; ipartnf++) {
        Int_t ipdg = pdgnf[ipartnf];

        if (TString(TVirtualMC::GetMC()->GetName()) == "TGeant3")
            TVirtualMC::GetMC()->SetUserDecay(ipdg);   // Force the decay to be done w/external decayer

        pythia6.SetMDCY(pythia6.Pycomp(ipdg), 1, 1);   // Activate decay in pythia
    }

    // The following will print the decay modes
    pythia6.Pyupda(1, 6);

    // rho0 (113), rho+ (213), rho- (-213) and
    // D+(411) ,D-(-411),D0(421),D0bar(-421) have decay modes defined in
    // TGeant3::DefineParticles, but for these particles
    // those decay modes are overridden to make use of pythia6.
    const Int_t nparthq = 3;
    // rho0,rho+,rho-,D+,D-,D0,D0bar
    // Int_t pdghq[nparthq] = {113,213,-213,411,-411,421,-421};
    Int_t pdghq[nparthq] = {421, 3122, -3122};
    for (Int_t iparthq = 0; iparthq < nparthq; iparthq++) {
        Int_t ipdg = pdghq[iparthq];
        if (TString(TVirtualMC::GetMC()->GetName()) == "TGeant3")
            TVirtualMC::GetMC()->SetUserDecay(ipdg);   // Force the decay to be done w/external decayer
        pythia6.SetMDCY(pythia6.Pycomp(ipdg), 1, 1);   // Activate decay in pythia
    }
    // Set pi0 to be stable in pythia6 so that Geant3 can handle decay.
    // In general, TGeant3 is set up through TGeant3gu::gudcay to pass
    // all pythia6 decay products back to the G3 transport mechanism if they
    // have a lifetime > 1.E-15 sec for further transport.
    // Since the pi0 lifetime is less than this, if pi0 is produced as a decay
    // product in pythia6, e.g. KL0 -> pi0 pi+ pi-, the pi0 will be immediately
    // decayed by pythia6 to 2 gammas, and the KL0 decay product list passed
    // back to the transport mechanism will be "gamma gamma pi+ pi-", i.e.
    // the pi0 will not be visible in the list of secondaries passed back to
    // the transport mechanism and will not be pushed to the stack for possible
    // storage to the stdhep output array.
    // To avoid this, the pi0 is set to stable in pythia6, and its decay
    // will be handled by Geant3.
    // pythia6.SetMDCY(pythia6.Pycomp(111),1,0);
    //}
}