+
+ //Update the energy deposit tables
+ Int_t i;
+ for(i=0;i<fEventEnergy.GetSize();i++)
+ {
+ fSummEnergy[i]+=fEventEnergy[i];
+ fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
+ }
+
+ AliRunLoader *runloader=AliRunLoader::Instance();
+
+ AliHeader* header = runloader->GetHeader();
+ AliStack* stack = runloader->Stack();
+ if ( (header == 0x0) || (stack == 0x0) )
+ {//check if we got header and stack. If not cry and exit aliroot
+ AliFatal("Can not get the stack or header from LOADER");
+ return;//never reached
+ }
+ // Update Header information
+ header->SetNprimary(stack->GetNprimary());
+ header->SetNtrack(stack->GetNtrack());
+
+ // Write out the kinematics
+ if (!AliSimulation::Instance()->Lego()) stack->FinishEvent();
+
+ // Synchronize the TreeTR with TreeK
+ if (fTmpTreeTR) ReorderAndExpandTreeTR();
+
+ // Write out the event Header information
+ TTree* treeE = runloader->TreeE();
+ if (treeE)
+ {
+ header->SetStack(stack);
+ treeE->Fill();
+ }
+ else
+ {
+ AliError("Can not get TreeE from RL");
+ }
+
+ if(AliSimulation::Instance()->Lego() == 0x0)
+ {
+ runloader->WriteKinematics("OVERWRITE");
+ runloader->WriteTrackRefs("OVERWRITE");
+ runloader->WriteHits("OVERWRITE");
+ }
+
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, " FINISHING EVENT ");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+}
+
+//_______________________________________________________________________
+void AliMC::Init()
+{
+ // MC initialization
+
+ //=================Create Materials and geometry
+ gMC->Init();
+ // Set alignable volumes for the whole geometry (with old root)
+#if ROOT_VERSION_CODE < 331527
+ SetAllAlignableVolumes();
+#endif
+ //Read the cuts for all materials
+ ReadTransPar();
+ //Build the special IMEDIA table
+ MediaTable();
+
+ //Compute cross-sections
+ gMC->BuildPhysics();
+
+ //Initialise geometry deposition table
+ fEventEnergy.Set(gMC->NofVolumes()+1);
+ fSummEnergy.Set(gMC->NofVolumes()+1);
+ fSum2Energy.Set(gMC->NofVolumes()+1);
+
+ // Register MC in configuration
+ AliConfig::Instance()->Add(gMC);
+
+}
+
+//_______________________________________________________________________
+void AliMC::MediaTable()
+{
+ //
+ // Built media table to get from the media number to
+ // the detector id
+ //
+
+ Int_t kz, nz, idt, lz, i, k, ind;
+ // Int_t ibeg;
+ TObjArray &dets = *gAlice->Detectors();
+ AliModule *det;
+ Int_t ndets=gAlice->GetNdets();
+ //
+ // For all detectors
+ for (kz=0;kz<ndets;kz++) {
+ // If detector is defined
+ if((det=dynamic_cast<AliModule*>(dets[kz]))) {
+ TArrayI &idtmed = *(det->GetIdtmed());
+ for(nz=0;nz<100;nz++) {
+
+ // Find max and min material number
+ if((idt=idtmed[nz])) {
+ det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
+ det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
+ }
+ }
+ if(det->LoMedium() > det->HiMedium()) {
+ det->LoMedium() = 0;
+ det->HiMedium() = 0;
+ } else {
+ if(det->HiMedium() > fImedia->GetSize()) {
+ AliError(Form("Increase fImedia from %d to %d",
+ fImedia->GetSize(),det->HiMedium()));
+ return;
+ }
+ // Tag all materials in rage as belonging to detector kz
+ for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
+ (*fImedia)[lz]=kz;
+ }
+ }
+ }
+ }
+ //
+ // Print summary table
+ AliInfo("Tracking media ranges:");
+ ToAliInfo(
+ for(i=0;i<(ndets-1)/6+1;i++) {
+ for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
+ ind=i*6+k;
+ det=dynamic_cast<AliModule*>(dets[ind]);
+ if(det)
+ printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
+ det->HiMedium());
+ else
+ printf(" %6s: %3d -> %3d;","NULL",0,0);
+ }
+ printf("\n");
+ }
+ )
+}
+
+//_______________________________________________________________________
+void AliMC::ReadTransPar()
+{
+ //
+ // Read filename to set the transport parameters
+ //
+
+
+ const Int_t kncuts=10;
+ const Int_t knflags=12;
+ const Int_t knpars=kncuts+knflags;
+ const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
+ "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
+ "BREM","COMP","DCAY","DRAY","HADR","LOSS",
+ "MULS","PAIR","PHOT","RAYL","STRA"};
+ char line[256];
+ char detName[7];
+ char* filtmp;
+ Float_t cut[kncuts];
+ Int_t flag[knflags];
+ Int_t i, itmed, iret, ktmed, kz;
+ FILE *lun;
+ //
+ // See whether the file is there
+ filtmp=gSystem->ExpandPathName(fTransParName.Data());
+ lun=fopen(filtmp,"r");
+ delete [] filtmp;
+ if(!lun) {
+ AliWarning(Form("File %s does not exist!",fTransParName.Data()));
+ return;
+ }
+ //
+ while(1) {
+ // Initialise cuts and flags
+ for(i=0;i<kncuts;i++) cut[i]=-99;
+ for(i=0;i<knflags;i++) flag[i]=-99;
+ itmed=0;
+ memset(line,0,256);
+ // Read up to the end of line excluded
+ iret=fscanf(lun,"%255[^\n]",line);
+ if(iret<0) {
+ //End of file
+ fclose(lun);
+ return;
+ }
+ // Read the end of line
+ fscanf(lun,"%*c");
+ if(!iret) continue;
+ if(line[0]=='*') continue;
+ // Read the numbers
+ iret=sscanf(line,"%s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d %d",
+ detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
+ &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
+ &flag[8],&flag[9],&flag[10],&flag[11]);
+ if(!iret) continue;
+ if(iret<0) {
+ //reading error
+ AliWarning(Form("Error reading file %s",fTransParName.Data()));
+ continue;
+ }
+ // Check that the module exist
+ AliModule *mod = gAlice->GetModule(detName);
+ if(mod) {
+ // Get the array of media numbers
+ TArrayI &idtmed = *mod->GetIdtmed();
+ // Check that the tracking medium code is valid
+ if(0<=itmed && itmed < 100) {
+ ktmed=idtmed[itmed];
+ if(!ktmed) {
+ AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName()));
+ continue;
+ }
+ // Set energy thresholds
+ for(kz=0;kz<kncuts;kz++) {
+ if(cut[kz]>=0) {
+ AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
+ kpars[kz],cut[kz],itmed,mod->GetName()));
+ gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
+ }
+ }
+ // Set transport mechanisms
+ for(kz=0;kz<knflags;kz++) {
+ if(flag[kz]>=0) {
+ AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
+ kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
+ gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
+ }
+ }
+ } else {
+ AliWarning(Form("Invalid medium code %d",itmed));
+ continue;
+ }
+ } else {
+ AliDebug(1, Form("%s not present",detName));
+ continue;
+ }
+ }
+}
+
+//_______________________________________________________________________
+void AliMC::SetTransPar(const char *filename)
+{
+ //
+ // Sets the file name for transport parameters
+ //
+ fTransParName = filename;
+}
+
+//_______________________________________________________________________
+void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
+{
+ //
+ // Add a hit to detector id
+ //
+ TObjArray &dets = *gAlice->Modules();
+ if(dets[id]) static_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
+}
+
+//_______________________________________________________________________
+void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
+{
+ //
+ // Add digit to detector id
+ //
+ TObjArray &dets = *gAlice->Modules();
+ if(dets[id]) static_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
+}
+
+//_______________________________________________________________________
+Int_t AliMC::GetCurrentTrackNumber() const {
+ //
+ // Returns current track
+ //
+ return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber();
+}
+
+//_______________________________________________________________________
+void AliMC::DumpPart (Int_t i) const
+{
+ //
+ // Dumps particle i in the stack
+ //
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader->Stack())
+ runloader->Stack()->DumpPart(i);
+}
+
+//_______________________________________________________________________
+void AliMC::DumpPStack () const
+{
+ //
+ // Dumps the particle stack
+ //
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader->Stack())
+ runloader->Stack()->DumpPStack();
+}
+
+//_______________________________________________________________________
+Int_t AliMC::GetNtrack() const {
+ //
+ // Returns number of tracks in stack
+ //
+ Int_t ntracks = -1;
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader->Stack())
+ ntracks = runloader->Stack()->GetNtrack();
+ return ntracks;
+}
+
+//_______________________________________________________________________
+Int_t AliMC::GetPrimary(Int_t track) const
+{
+ //
+ // return number of primary that has generated track
+ //
+ Int_t nprimary = -999;
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader->Stack())
+ nprimary = runloader->Stack()->GetPrimary(track);
+ return nprimary;
+}
+
+//_______________________________________________________________________
+TParticle* AliMC::Particle(Int_t i) const
+{
+ // Returns the i-th particle from the stack taking into account
+ // the remaping done by PurifyKine
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader)
+ if (runloader->Stack())
+ return runloader->Stack()->Particle(i);
+ return 0x0;
+}
+
+//_______________________________________________________________________
+const TObjArray* AliMC::Particles() const {
+ //
+ // Returns pointer to Particles array
+ //
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader)
+ if (runloader->Stack())
+ return runloader->Stack()->Particles();
+ return 0x0;
+}
+
+//_______________________________________________________________________
+void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom,
+ const Float_t *vpos, const Float_t *polar, Float_t tof,
+ TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
+{
+// Delegate to stack
+//
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader)
+ if (runloader->Stack())
+ runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
+ mech, ntr, weight, is);
+}
+
+//_______________________________________________________________________
+void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
+ Double_t px, Double_t py, Double_t pz, Double_t e,
+ Double_t vx, Double_t vy, Double_t vz, Double_t tof,
+ Double_t polx, Double_t poly, Double_t polz,
+ TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
+{
+ // Delegate to stack
+ //
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader)
+ if (runloader->Stack())
+ runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
+ polx, poly, polz, mech, ntr, weight, is);
+}
+
+//_______________________________________________________________________
+void AliMC::SetHighWaterMark(Int_t nt) const
+{
+ //
+ // Set high water mark for last track in event
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader)
+ if (runloader->Stack())
+ runloader->Stack()->SetHighWaterMark(nt);
+}
+
+//_______________________________________________________________________
+void AliMC::KeepTrack(Int_t track) const
+{
+ //
+ // Delegate to stack
+ //
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader)
+ if (runloader->Stack())
+ runloader->Stack()->KeepTrack(track);
+}
+
+//_______________________________________________________________________
+void AliMC::FlagTrack(Int_t track) const
+{
+ // Delegate to stack
+ //
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader)
+ if (runloader->Stack())
+ runloader->Stack()->FlagTrack(track);
+}
+
+//_______________________________________________________________________
+void AliMC::SetCurrentTrack(Int_t track) const
+{
+ //
+ // Set current track number
+ //
+ AliRunLoader * runloader = AliRunLoader::Instance();
+ if (runloader)
+ if (runloader->Stack())
+ runloader->Stack()->SetCurrentTrack(track);
+}
+
+//_______________________________________________________________________
+AliTrackReference* AliMC::AddTrackReference(Int_t label, Int_t id)
+{
+ //
+ // add a trackrefernce to the list
+ Int_t primary = GetPrimary(label);
+ Particle(primary)->SetBit(kKeepBit);
+
+ Int_t nref = fTmpTrackReferences.GetEntriesFast();
+ return new(fTmpTrackReferences[nref]) AliTrackReference(label, id);
+}
+
+
+
+//_______________________________________________________________________
+void AliMC::ResetTrackReferences()
+{
+ //
+ // Reset all references
+ //
+ fTmpTrackReferences.Clear();
+}
+
+//_______________________________________________________________________
+void AliMC::RemapTrackReferencesIDs(const Int_t *map)
+{
+ //
+ // Remapping track reference
+ // Called at finish primary
+ //
+
+ Int_t nEntries = fTmpTrackReferences.GetEntries();
+ for (Int_t i=0; i < nEntries; i++){
+ AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTmpTrackReferences.UncheckedAt(i));
+ if (ref) {
+ Int_t newID = map[ref->GetTrack()];
+ if (newID>=0) ref->SetTrack(newID);
+ else {
+ ref->SetBit(kNotDeleted,kFALSE);
+ fTmpTrackReferences.RemoveAt(i);
+ }
+ } // if ref
+ }
+ 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", &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();
+ // make branch for central track references
+ TClonesArray* pRef = &fTrackReferences;
+ treeTR->Branch("TrackReferences", &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");