+ fTmpTrackReferences.Compress();
+}
+
+//_______________________________________________________________________
+void AliMC::FixParticleDecaytime()
+{
+ //
+ // Fix the particle decay time according to rmin and rmax for decays
+ //
+
+ TLorentzVector p;
+ gMC->TrackMomentum(p);
+ Double_t tmin, tmax;
+ Double_t b;
+
+ // Transverse velocity
+ Double_t vt = p.Pt() / p.E();
+
+ if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) { // [kG]
+
+ // Radius of helix
+
+ Double_t rho = p.Pt() / 0.0003 / b; // [cm]
+
+ // Revolution frequency
+
+ Double_t omega = vt / rho;
+
+ // Maximum and minimum decay time
+ //
+ // Check for curlers first
+ const Double_t kOvRhoSqr2 = 1./(rho*TMath::Sqrt(2.));
+ if (fRDecayMax * kOvRhoSqr2 > 1.) return;
+
+ //
+
+ tmax = TMath::ACos((1.-fRDecayMax*kOvRhoSqr2)*(1.+fRDecayMax*kOvRhoSqr2)) / omega; // [ct]
+ tmin = TMath::ACos((1.-fRDecayMin*kOvRhoSqr2)*(1.+fRDecayMin*kOvRhoSqr2)) / omega; // [ct]
+ } else {
+ tmax = fRDecayMax / vt; // [ct]
+ tmin = fRDecayMin / vt; // [ct]
+ }
+
+ //
+ // Dial t using the two limits
+ Double_t t = tmin + (tmax - tmin) * gRandom->Rndm(); // [ct]
+ //
+ //
+ // Force decay time in transport code
+ //
+ gMC->ForceDecayTime(t / 2.99792458e10);
+}
+
+void AliMC::MakeTmpTrackRefsTree()
+{
+ // Make the temporary track reference tree
+ fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
+ fTmpTreeTR = new TTree("TreeTR", "Track References");
+ TClonesArray* pRef = &fTmpTrackReferences;
+ fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &pRef, 4000);
+}
+
+//_______________________________________________________________________
+void AliMC::ReorderAndExpandTreeTR()
+{
+//
+// Reorder and expand the temporary track reference tree in order to match the kinematics tree
+//
+
+ AliRunLoader *rl = AliRunLoader::Instance();
+//
+// TreeTR
+ AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
+ rl->MakeTrackRefsContainer();
+ TTree * treeTR = rl->TreeTR();
+ if (treeTR){
+ // make branch for central track references
+ TBranch *branch;
+ TClonesArray* pRef = &fTrackReferences;
+ branch = treeTR->Branch("TrackReferences", &pRef);
+ branch->SetAddress(&pRef);
+ }
+
+ AliStack* stack = rl->Stack();
+ Int_t np = stack->GetNprimary();
+ Int_t nt = fTmpTreeTR->GetEntries();
+ //
+ // Loop over tracks and find the secondaries with the help of the kine tree
+ Int_t ifills = 0;
+ Int_t it = 0;
+ for (Int_t ip = np - 1; ip > -1; ip--) {
+ TParticle *part = stack->Particle(ip);
+ //printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter());
+
+ // Skip primaries that have not been transported
+ Int_t dau1 = part->GetFirstDaughter();
+ Int_t dau2 = -1;
+ if (!part->TestBit(kTransportBit)) continue;
+ //
+ fTmpTreeTR->GetEntry(it++);
+ Int_t nh = fTmpTrackReferences.GetEntries();
+ // Determine range of secondaries produced by this primary
+ if (dau1 > -1) {
+ Int_t inext = ip - 1;
+ while (dau2 < 0) {
+ if (inext >= 0) {
+ part = stack->Particle(inext);
+ dau2 = part->GetFirstDaughter();
+ if (!(part->TestBit(kTransportBit)) || dau2 == -1 || dau2 < np) {
+// if (dau2 == -1 || dau2 < np) {
+ dau2 = -1;
+ } else {
+ dau2--;
+ }
+ } else {
+ dau2 = stack->GetNtrack() - 1;
+ }
+ inext--;
+ } // find upper bound
+ } // dau2 < 0
+// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
+ //
+ // Loop over reference hits and find secondary label
+ for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
+ for (Int_t ih = 0; ih < nh; ih++) {
+ AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih);
+ Int_t label = tr->Label();
+ // Skip primaries
+ if (label == ip) continue;
+ if (label > dau2 || label < dau1)
+ AliWarning(Form("Track Reference Label out of range !: %5d %5d %5d \n", label, dau1, dau2));
+ if (label == id) {
+ // secondary found
+ Int_t nref = fTrackReferences.GetEntriesFast();
+ new(fTrackReferences[nref]) AliTrackReference(*tr);
+ }
+ } // hits
+ treeTR->Fill();
+ fTrackReferences.Clear();
+ ifills++;
+ } // daughters
+ } // tracks
+ //
+ // Now loop again and write the primaries
+ it = nt - 1;
+ for (Int_t ip = 0; ip < np; ip++) {
+ TParticle* part = stack->Particle(ip);
+// if ((part->GetFirstDaughter() == -1 && part->GetStatusCode() <= 1) || part->GetFirstDaughter() >= np)
+ if (part->TestBit(kTransportBit))
+ {
+ // Skip particles that have not been transported
+ fTmpTreeTR->GetEntry(it--);
+ Int_t nh = fTmpTrackReferences.GetEntries();
+ //
+ // Loop over reference hits and find primary labels
+ for (Int_t ih = 0; ih < nh; ih++) {
+ AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih);
+ Int_t label = tr->Label();
+ if (label == ip) {
+ Int_t nref = fTrackReferences.GetEntriesFast();
+ new(fTrackReferences[nref]) AliTrackReference(*tr);
+ }
+ }
+ }
+ treeTR->Fill();
+ fTrackReferences.Clear();
+ ifills++;
+ } // tracks
+ // Check
+ if (ifills != stack->GetNtrack())
+ AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack()));
+//
+// Clean-up
+ delete fTmpTreeTR;
+ fTmpFileTR->Close();
+ delete fTmpFileTR;
+ fTmpTrackReferences.Clear();
+ gSystem->Exec("rm -rf TrackRefsTmp.root");