X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliMC.cxx;h=12dd60c4b4e1c939ae41b14bede5b6a3a6df2de6;hb=098cca5293b87b12dcc1564932a3a76117fdfca8;hp=8b14fc7faf0ef773ddf65ec2b08db3b06c704779;hpb=392808e852c26a813bd17ef6aec482eeccd9027d;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliMC.cxx b/STEER/AliMC.cxx index 8b14fc7faf0..12dd60c4b4e 100644 --- a/STEER/AliMC.cxx +++ b/STEER/AliMC.cxx @@ -22,27 +22,33 @@ // Federico.Carminati@cern.ch #include -#include +#include #include +#include +#include #include +#include +#include #include #include #include -#include "AliLog.h" +#include "AliCDBEntry.h" +#include "AliCDBManager.h" +#include "AliCDBStorage.h" #include "AliDetector.h" #include "AliGenerator.h" +#include "AliGeomManager.h" #include "AliHeader.h" +#include "AliHit.h" #include "AliLego.h" +#include "AliLog.h" #include "AliMC.h" -#include "AliMCQA.h" +#include "AliMagF.h" #include "AliRun.h" +#include "AliSimulation.h" #include "AliStack.h" -#include "AliMagF.h" #include "AliTrackReference.h" -#include "AliSimulation.h" -#include "AliGeomManager.h" - ClassImp(AliMC) @@ -59,9 +65,11 @@ AliMC::AliMC() : fDecayPdg(0), fImedia(0), fTransParName("\0"), - fMCQA(0), fHitLists(0), - fTrackReferences(0) + fTmpTreeTR(0), + fTmpFileTR(0), + fTrackReferences(), + fTmpTrackReferences() { //default constructor @@ -82,9 +90,11 @@ AliMC::AliMC(const char *name, const char *title) : fDecayPdg(0), fImedia(new TArrayI(1000)), fTransParName("\0"), - fMCQA(0), fHitLists(new TList()), - fTrackReferences(new TClonesArray("AliTrackReference", 100)) + fTmpTreeTR(0), + fTmpFileTR(0), + fTrackReferences("AliTrackReference", 100), + fTmpTrackReferences("AliTrackReference", 100) { //constructor // Set transport parameters @@ -94,52 +104,14 @@ AliMC::AliMC(const char *name, const char *title) : for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99; } -//_______________________________________________________________________ -AliMC::AliMC(const AliMC &mc) : - TVirtualMCApplication(mc), - fGenerator(0), - fEventEnergy(0), - fSummEnergy(0), - fSum2Energy(0), - fTrRmax(1.e10), - fTrZmax(1.e10), - fRDecayMax(1.e10), - fRDecayMin(-1.), - fDecayPdg(0), - fImedia(0), - fTransParName("\0"), - fMCQA(0), - fHitLists(0), - fTrackReferences(0) -{ - // - // Copy constructor for AliMC - // - mc.Copy(*this); -} - //_______________________________________________________________________ AliMC::~AliMC() { //destructor delete fGenerator; delete fImedia; - delete fMCQA; delete fHitLists; // Delete track references - if (fTrackReferences) { - fTrackReferences->Delete(); - delete fTrackReferences; - fTrackReferences = 0; - } - -} - -//_______________________________________________________________________ -void AliMC::Copy(TObject &) const -{ - //dummy Copy function - AliFatal("Not implemented!"); } //_______________________________________________________________________ @@ -152,18 +124,29 @@ void AliMC::ConstructGeometry() // at InitGeometry(). // - if(gAlice->IsRootGeometry()){ - // Load geometry - const char *geomfilename = gAlice->GetGeometryFileName(); - if(gSystem->ExpandPathName(geomfilename)){ - AliInfo(Form("Loading geometry from file:\n %40s\n\n",geomfilename)); - TGeoManager::Import(geomfilename); + if(gAlice->IsRootGeometry()){ //load geometry either from CDB or from file + if(gAlice->IsGeomFromCDB()){ + AliInfo("Loading geometry from CDB default storage"); + AliCDBPath path("GRP","Geometry","Data"); + AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath()); + if(!entry) AliFatal("Unable to load geometry from CDB!"); + entry->SetOwner(0); + gGeoManager = (TGeoManager*) entry->GetObject(); + if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!"); }else{ - AliInfo(Form("Geometry file %40s not found!\n",geomfilename)); - return; + // Load geometry + const char *geomfilename = gAlice->GetGeometryFileName(); + if(gSystem->ExpandPathName(geomfilename)){ + AliInfo(Form("Loading geometry from file:\n %40s",geomfilename)); + TGeoManager::Import(geomfilename); + }else{ + AliInfo(Form("Geometry file %40s not found!\n",geomfilename)); + return; + } } }else{ // Create modules, materials, geometry + if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry"); TStopwatch stw; TIter next(gAlice->Modules()); AliModule *detector; @@ -185,12 +168,17 @@ Bool_t AliMC::MisalignGeometry() { // Call misalignment code if AliSimulation object was defined. - //Set alignable volumes for the whole geometry - SetAllAlignableVolumes(); + if(!gAlice->IsRootGeometry()){ + //Set alignable volumes for the whole geometry + SetAllAlignableVolumes(); + } // Misalign geometry via AliSimulation instance - if (!AliSimulation::GetInstance()) return kFALSE; + if (!AliSimulation::Instance()) return kFALSE; AliGeomManager::SetGeometry(gGeoManager); - return AliSimulation::GetInstance()->MisalignGeometry(gAlice->GetRunLoader()); + if(!AliGeomManager::CheckSymNamesLUT("ALL")) + AliFatal("Current loaded geometry differs in the definition of symbolic names!"); + + return AliSimulation::Instance()->MisalignGeometry(AliRunLoader::Instance()); } //_______________________________________________________________________ @@ -306,7 +294,6 @@ void AliMC::BeginPrimary() // Reset Hits info ResetHits(); ResetTrackReferences(); - } //_______________________________________________________________________ @@ -320,7 +307,6 @@ void AliMC::PreTrack() if((module = dynamic_cast(dets[i]))) module->PreTrack(); - fMCQA->PreTrack(); } //_______________________________________________________________________ @@ -329,7 +315,6 @@ void AliMC::Stepping() // // Called at every step during transport // - Int_t id = DetFromMate(gMC->CurrentMedium()); if (id < 0) return; @@ -347,22 +332,25 @@ void AliMC::Stepping() // // --- If lego option, do it and leave - if (gAlice->Lego()) - gAlice->Lego()->StepManager(); + if (AliSimulation::Instance()->Lego()) + AliSimulation::Instance()->Lego()->StepManager(); else { Int_t copy; //Update energy deposition tables AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep()); // // write tracke reference for track which is dissapearing - MI - if (gMC->IsTrackDisappeared()) { - if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber()); + + if (gMC->IsTrackDisappeared() && !(gMC->IsTrackAlive())) { + if (gMC->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(), + AliTrackReference::kDisappeared); + + } - + //Call the appropriate stepping routine; AliModule *det = dynamic_cast(gAlice->Modules()->At(id)); if(det && det->StepManagerIsEnabled()) { - if(AliLog::GetGlobalDebugLevel()>0) fMCQA->StepManager(id); det->StepManager(); } } @@ -371,7 +359,7 @@ void AliMC::Stepping() //_______________________________________________________________________ void AliMC::EnergySummary() { - // + //e // Print summary of deposited energy // @@ -380,7 +368,7 @@ void AliMC::EnergySummary() Float_t ed, ed2; Int_t kn, i, left, j, id; const Float_t kzero=0; - Int_t ievent=gAlice->GetRunLoader()->GetHeader()->GetEvent()+1; + Int_t ievent=AliRunLoader::Instance()->GetHeader()->GetEvent()+1; // // Energy loss information if(ievent) { @@ -448,7 +436,7 @@ void AliMC::BeginEvent() AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - AliRunLoader *runloader=gAlice->GetRunLoader(); + AliRunLoader *runloader=AliRunLoader::Instance(); /*******************************/ /* Clean after eventual */ @@ -465,18 +453,19 @@ void AliMC::BeginEvent() // Clean detector information if (runloader->Stack()) - runloader->Stack()->Reset();//clean stack -> tree is unloaded + runloader->Stack()->Reset();//clean stack -> tree is unloaded else - runloader->MakeStack();//or make a new one + runloader->MakeStack();//or make a new one + + + if(AliSimulation::Instance()->Lego() == 0x0) + { + AliDebug(1, "fRunLoader->MakeTree(K)"); + runloader->MakeTree("K"); + } - if(gAlice->Lego() == 0x0) - { - AliDebug(1, "fRunLoader->MakeTree(K)"); - runloader->MakeTree("K"); - } - AliDebug(1, "gMC->SetStack(fRunLoader->Stack())"); - gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here + gMC->SetStack(runloader->Stack());//Was in InitMC - but was moved here //because we don't have guarantee that //stack pointer is not going to change from event to event //since it bellobgs to header and is obtained via RunLoader @@ -484,15 +473,15 @@ void AliMC::BeginEvent() // Reset all Detectors & kinematics & make/reset trees // - runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(), - gAlice->GetEventNrInRun()); + runloader->GetHeader()->Reset(AliCDBManager::Instance()->GetRun(),gAlice->GetEvNumber(), + gAlice->GetEventNrInRun()); // fRunLoader->WriteKinematics("OVERWRITE"); is there any reason to rewrite here since MakeTree does so - if(gAlice->Lego()) - { - gAlice->Lego()->BeginEvent(); - return; - } + if(AliSimulation::Instance()->Lego()) + { + AliSimulation::Instance()->Lego()->BeginEvent(); + return; + } AliDebug(1, "ResetHits()"); @@ -501,32 +490,17 @@ void AliMC::BeginEvent() AliDebug(1, "fRunLoader->MakeTree(H)"); runloader->MakeTree("H"); - AliDebug(1, "fRunLoader->MakeTrackRefsContainer()"); - runloader->MakeTrackRefsContainer();//for insurance + MakeTmpTrackRefsTree(); //create new branches and SetAdresses TIter next(gAlice->Modules()); AliModule *detector; while((detector = (AliModule*)next())) { - AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName())); - detector->MakeBranch("H"); - AliDebug(2, Form("%s->MakeBranchTR()",detector->GetName())); - detector->MakeBranchTR(); - AliDebug(2, Form("%s->SetTreeAddress()",detector->GetName())); - detector->SetTreeAddress(); + AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName())); + detector->MakeBranch("H"); } - // make branch for AliRun track References - TTree * treeTR = runloader->TreeTR(); - if (treeTR){ - // make branch for central track references - if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0); - TBranch *branch; - branch = treeTR->Branch("AliRun",&fTrackReferences); - branch->SetAddress(&fTrackReferences); - } - // } //_______________________________________________________________________ @@ -542,6 +516,32 @@ void AliMC::ResetHits() } } +//_______________________________________________________________________ +void AliMC::ResetDigits() +{ + // + // Reset all Detectors digits + // + TIter next(gAlice->Modules()); + AliModule *detector; + while((detector = dynamic_cast(next()))) { + detector->ResetDigits(); + } +} + +//_______________________________________________________________________ +void AliMC::ResetSDigits() +{ + // + // Reset all Detectors digits + // + TIter next(gAlice->Modules()); + AliModule *detector; + while((detector = dynamic_cast(next()))) { + detector->ResetSDigits(); + } +} + //_______________________________________________________________________ void AliMC::PostTrack() { @@ -560,15 +560,16 @@ void AliMC::FinishPrimary() // // Called at the end of each primary track // - AliRunLoader *runloader=gAlice->GetRunLoader(); + AliRunLoader *runloader=AliRunLoader::Instance(); // static Int_t count=0; // const Int_t times=10; // This primary is finished, purify stack #if ROOT_VERSION_CODE > 262152 - if (!(gMC->SecondariesAreOrdered())) - runloader->Stack()->ReorderKine(); + if (!(gMC->SecondariesAreOrdered())) { + if (runloader->Stack()->ReorderKine()) RemapHits(); + } #endif - runloader->Stack()->PurifyKine(); + if (runloader->Stack()->PurifyKine()) RemapHits(); TIter next(gAlice->Modules()); AliModule *detector; @@ -583,7 +584,42 @@ void AliMC::FinishPrimary() } // Write out track references if any - if (runloader->TreeTR()) runloader->TreeTR()->Fill(); + if (fTmpTreeTR) fTmpTreeTR->Fill(); +} + +void AliMC::RemapHits() +{ +// +// Remaps the track labels of the hits + AliRunLoader *runloader=AliRunLoader::Instance(); + AliStack* stack = runloader->Stack(); + TList* hitLists = GetHitLists(); + TIter next(hitLists); + TCollection *hitList; + + while((hitList = dynamic_cast(next()))) { + TIter nexthit(hitList); + AliHit *hit; + while((hit = dynamic_cast(nexthit()))) { + hit->SetTrack(stack->TrackLabel(hit->GetTrack())); + } + } + + // + // This for detectors which have a special mapping mechanism + // for hits, such as TPC and TRD + // + + + TObjArray* modules = gAlice->Modules(); + TIter nextmod(modules); + AliModule *module; + while((module = (AliModule*) nextmod())) { + AliDetector* det = dynamic_cast (module); + if (det) det->RemapTrackHitIDs(stack->TrackLabelMap()); + } + // + RemapTrackReferencesIDs(stack->TrackLabelMap()); } //_______________________________________________________________________ @@ -595,7 +631,7 @@ void AliMC::FinishEvent() // - if(gAlice->Lego()) gAlice->Lego()->FinishEvent(); + if(AliSimulation::Instance()->Lego()) AliSimulation::Instance()->Lego()->FinishEvent(); TIter next(gAlice->Modules()); AliModule *detector; @@ -611,7 +647,7 @@ void AliMC::FinishEvent() fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i]; } - AliRunLoader *runloader=gAlice->GetRunLoader(); + AliRunLoader *runloader=AliRunLoader::Instance(); AliHeader* header = runloader->GetHeader(); AliStack* stack = runloader->Stack(); @@ -624,9 +660,11 @@ void AliMC::FinishEvent() header->SetNprimary(stack->GetNprimary()); header->SetNtrack(stack->GetNtrack()); - // Write out the kinematics - if (!gAlice->Lego()) stack->FinishEvent(); + if (!AliSimulation::Instance()->Lego()) stack->FinishEvent(); + + // Synchronize the TreeTR with TreeK + if (fTmpTreeTR) ReorderAndExpandTreeTR(); // Write out the event Header information TTree* treeE = runloader->TreeE(); @@ -640,7 +678,7 @@ void AliMC::FinishEvent() AliError("Can not get TreeE from RL"); } - if(gAlice->Lego() == 0x0) + if(AliSimulation::Instance()->Lego() == 0x0) { runloader->WriteKinematics("OVERWRITE"); runloader->WriteTrackRefs("OVERWRITE"); @@ -656,13 +694,6 @@ void AliMC::FinishEvent() AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); } -//_______________________________________________________________________ -void AliMC::Field(const Double_t* x, Double_t* b) const -{ - // Calculates field "b" at point "x" - gAlice->Field(x,b); -} - //_______________________________________________________________________ void AliMC::Init() { @@ -687,9 +718,6 @@ void AliMC::Init() fSummEnergy.Set(gMC->NofVolumes()+1); fSum2Energy.Set(gMC->NofVolumes()+1); - // - fMCQA = new AliMCQA(gAlice->GetNdets()); - // Register MC in configuration AliConfig::Instance()->Add(gMC); @@ -865,18 +893,6 @@ void AliMC::SetTransPar(const char *filename) fTransParName = filename; } -//_______________________________________________________________________ -void AliMC::Browse(TBrowser *b) -{ - // - // Called when the item "Run" is clicked on the left pane - // of the Root browser. - // It displays the Root Trees and all detectors. - // - //detectors are in folders anyway - b->Add(fMCQA,"AliMCQA"); -} - //_______________________________________________________________________ void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const { @@ -902,7 +918,7 @@ Int_t AliMC::GetCurrentTrackNumber() const { // // Returns current track // - return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber(); + return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber(); } //_______________________________________________________________________ @@ -911,7 +927,7 @@ void AliMC::DumpPart (Int_t i) const // // Dumps particle i in the stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) runloader->Stack()->DumpPart(i); } @@ -922,7 +938,7 @@ void AliMC::DumpPStack () const // // Dumps the particle stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) runloader->Stack()->DumpPStack(); } @@ -933,7 +949,7 @@ Int_t AliMC::GetNtrack() const { // Returns number of tracks in stack // Int_t ntracks = -1; - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) ntracks = runloader->Stack()->GetNtrack(); return ntracks; @@ -946,7 +962,7 @@ Int_t AliMC::GetPrimary(Int_t track) const // return number of primary that has generated track // Int_t nprimary = -999; - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) nprimary = runloader->Stack()->GetPrimary(track); return nprimary; @@ -957,7 +973,7 @@ 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 = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) return runloader->Stack()->Particle(i); @@ -965,11 +981,11 @@ TParticle* AliMC::Particle(Int_t i) const } //_______________________________________________________________________ -TObjArray* AliMC::Particles() const { +const TObjArray* AliMC::Particles() const { // // Returns pointer to Particles array // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) return runloader->Stack()->Particles(); @@ -977,13 +993,13 @@ TObjArray* AliMC::Particles() const { } //_______________________________________________________________________ -void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, - Float_t *vpos, Float_t *polar, Float_t tof, +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 = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof, @@ -999,7 +1015,7 @@ void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, { // Delegate to stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof, @@ -1011,7 +1027,7 @@ void AliMC::SetHighWaterMark(Int_t nt) const { // // Set high water mark for last track in event - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->SetHighWaterMark(nt); @@ -1023,7 +1039,7 @@ void AliMC::KeepTrack(Int_t track) const // // Delegate to stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->KeepTrack(track); @@ -1034,7 +1050,7 @@ void AliMC::FlagTrack(Int_t track) const { // Delegate to stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->FlagTrack(track); @@ -1046,23 +1062,22 @@ void AliMC::SetCurrentTrack(Int_t track) const // // Set current track number // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->SetCurrentTrack(track); } //_______________________________________________________________________ -void AliMC::AddTrackReference(Int_t label){ +AliTrackReference* AliMC::AddTrackReference(Int_t label, Int_t id) +{ // // add a trackrefernce to the list - if (!fTrackReferences) { - AliError("Container trackrefernce not active"); - return; - } - Int_t nref = fTrackReferences->GetEntriesFast(); - TClonesArray &lref = *fTrackReferences; - new(lref[nref]) AliTrackReference(label); + Int_t primary = GetPrimary(label); + Particle(primary)->SetBit(kKeepBit); + + Int_t nref = fTmpTrackReferences.GetEntriesFast(); + return new(fTmpTrackReferences[nref]) AliTrackReference(label, id); } @@ -1073,13 +1088,7 @@ void AliMC::ResetTrackReferences() // // Reset all references // - if (fTrackReferences) fTrackReferences->Clear(); - - TIter next(gAlice->Modules()); - AliModule *detector; - while((detector = dynamic_cast(next()))) { - detector->ResetTrackReferences(); - } + fTmpTrackReferences.Clear(); } void AliMC::RemapTrackReferencesIDs(Int_t *map) @@ -1088,20 +1097,20 @@ void AliMC::RemapTrackReferencesIDs(Int_t *map) // Remapping track reference // Called at finish primary // - if (!fTrackReferences) return; - Int_t nEntries = fTrackReferences->GetEntries(); + + Int_t nEntries = fTmpTrackReferences.GetEntries(); for (Int_t i=0; i < nEntries; i++){ - AliTrackReference * ref = dynamic_cast(fTrackReferences->UncheckedAt(i)); + AliTrackReference * ref = dynamic_cast(fTmpTrackReferences.UncheckedAt(i)); if (ref) { Int_t newID = map[ref->GetTrack()]; if (newID>=0) ref->SetTrack(newID); else { ref->SetBit(kNotDeleted,kFALSE); - fTrackReferences->RemoveAt(i); + fTmpTrackReferences.RemoveAt(i); } } // if ref } - fTrackReferences->Compress(); + fTmpTrackReferences.Compress(); } void AliMC::FixParticleDecaytime() @@ -1118,7 +1127,7 @@ void AliMC::FixParticleDecaytime() // Transverse velocity Double_t vt = p.Pt() / p.E(); - if ((b = gAlice->Field()->SolenoidField()) > 0.) { // [kG] + if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) { // [kG] // Radius of helix @@ -1151,3 +1160,130 @@ void AliMC::FixParticleDecaytime() // 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"); +} +