+//
+// Decay the pi0
+// and put decay products on the stack
+//
+ static TClonesArray *particles;
+ if(!particles) particles = new TClonesArray("TParticle",1000);
+//
+ const Float_t kMass = TDatabasePDG::Instance()->GetParticle(kPi0)->Mass();
+ Float_t e = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]+ kMass * kMass);
+//
+// Decay the pi0
+ TLorentzVector pmom(p[0], p[1], p[2], e);
+ fDecayer->Decay(kPi0, &pmom);
+
+//
+// Put decay particles on the stack
+//
+ Float_t polar[3] = {0., 0., 0.};
+ Int_t np = fDecayer->ImportParticles(particles);
+ fNpartProd += (np-1);
+ Int_t nt;
+ for (Int_t i = 1; i < np; i++)
+ {
+ TParticle* iParticle = (TParticle *) particles->At(i);
+ p[0] = iParticle->Px();
+ p[1] = iParticle->Py();
+ p[2] = iParticle->Pz();
+ Int_t part = iParticle->GetPdgCode();
+
+ PushTrack(fTrackIt, fNt, part, p, orig, polar, 0, kPDecay, nt, fParentWeight);
+ KeepTrack(nt);
+ }
+ fNt = nt;