X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStack.cxx;h=724e5310fd15c0457c04bc970d87553d5a05221f;hb=bc77d57e77b56b0d14f0643090c1e792eb530be2;hp=03bfdf84e4d5497446de648222bae9c2f082007f;hpb=e191bb57c96a66bf4f60381b9a3f8c511bc6188e;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx index 03bfdf84e4d..724e5310fd1 100644 --- a/STEER/AliStack.cxx +++ b/STEER/AliStack.cxx @@ -17,18 +17,21 @@ /////////////////////////////////////////////////////////////////////////////// // // -// Particles stack class +// Particles stack class // +// Implements the TMCVirtualStack of the Virtual Monte Carlo // +// Holds the particles transported during simulation // +// Is used to compare results of reconstruction with simulation // +// Author A.Morsch // // // /////////////////////////////////////////////////////////////////////////////// -#include -#include #include #include #include #include +#include "AliLog.h" #include "AliHit.h" #include "AliModule.h" #include "AliRun.h" @@ -44,6 +47,7 @@ AliStack::AliStack(): fParticleMap(0), fParticleFileMap(0), fParticleBuffer(0), + fCurrentTrack(0), fTreeK(0), fNtrack(0), fNprimary(0), @@ -64,6 +68,7 @@ AliStack::AliStack(Int_t size, const char* evfoldname): fParticleMap(new TObjArray(size)), fParticleFileMap(0), fParticleBuffer(0), + fCurrentTrack(0), fTreeK(0), fNtrack(0), fNprimary(0), @@ -85,13 +90,15 @@ AliStack::AliStack(const AliStack& st): fParticleMap(0), fParticleFileMap(0), fParticleBuffer(0), + fCurrentTrack(0), fTreeK(0), fNtrack(0), fNprimary(0), fCurrent(-1), fCurrentPrimary(-1), fHgwmk(0), - fLoadPoint(0) + fLoadPoint(0), + fEventFolderName() { // // Copy constructor @@ -102,7 +109,7 @@ AliStack::AliStack(const AliStack& st): //_______________________________________________________________________ void AliStack::Copy(TObject&) const { - Fatal("Copy","Not implemented!\n"); + AliFatal("Not implemented!"); } //_______________________________________________________________________ @@ -167,8 +174,8 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2], mech, ntr, weight, is); } else { - Warning("PushTrack", "Particle type %d not defined in PDG Database !\n", pdg); - Warning("PushTrack", "Particle skipped !\n"); + AliWarning(Form("Particle type %d not defined in PDG Database !", pdg)); + AliWarning("Particle skipped !"); } } @@ -203,7 +210,8 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Int_t kFirstDaughter=-1; const Int_t kLastDaughter=-1; - + + TClonesArray &particles = *fParticles; TParticle* particle = new(particles[fLoadPoint++]) @@ -219,6 +227,8 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, // Declare that the daughter information is valid particle->SetBit(kDaughtersBit); // Add the particle to the stack + + fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!! if(parent>=0) { @@ -228,7 +238,7 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); } else { - printf("Error in AliStack::PushTrack: Parent %d does not exist\n",parent); + AliError(Form("Parent %d does not exist",parent)); } } else { @@ -294,11 +304,11 @@ void AliStack::PurifyKine() TArrayI map(particles.GetLast()+1); // Save in Header total number of tracks before compression - // If no tracks generated return now if(fHgwmk+1 == fNtrack) return; // First pass, invalid Daughter information + for(i=0; iGetMCApp()->GetHitLists(); TIter next(hitLists); TCollection *hitList; + while((hitList = dynamic_cast(next()))) { TIter nexthit(hitList); AliHit *hit; while((hit = dynamic_cast(nexthit()))) { - hit->SetTrack(map[hit->GetTrack()]); + hit->SetTrack(map[hit->GetTrack()]); } } @@ -400,11 +412,140 @@ void AliStack::PurifyKine() Int_t toshrink = fNtrack-fHgwmk-1; fLoadPoint-=toshrink; + + for(i=fLoadPoint; iRemoveAt(i); fNtrack=nkeep; fHgwmk=nkeep-1; - // delete [] map; +} + +void AliStack::ReorderKine() +{ +// +// In some transport code children might not come in a continuous sequence. +// In this case the stack has to be reordered in order to establish the +// mother daughter relation using index ranges. +// + if(fHgwmk+1 == fNtrack) return; + + // + // Howmany secondaries have been produced ? + Int_t nNew = fNtrack - fHgwmk - 1; + + if (nNew > 0) { + Int_t i, j; + TObjArray &particles = *fParticleMap; + TArrayI map1(nNew); + // + // Copy pointers to temporary array + TParticle** tmp = new TParticle*[nNew]; + + for (i = 0; i < nNew; i++) { + if (particles.At(fHgwmk + 1 + i)) { + tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i)); + } else { + tmp[i] = 0x0; + } + map1[i] = -99; + } + + + // + // Reset LoadPoint + // + fLoadPoint = fHgwmk + 1; + // + // Re-Push particles into stack + // The outer loop is over parents, the inner over children. + // -1 refers to the primary particle + // + for (i = -1; i < nNew-1; i++) { + Int_t ipa; + TParticle* parP; + if (i == -1) { + ipa = tmp[0]->GetFirstMother(); + parP =dynamic_cast(particles.At(ipa)); + } else { + ipa = (fHgwmk + 1 + i); + // Skip deleted particles + if (!tmp[i]) continue; + // Skip particles without children + if (tmp[i]->GetFirstDaughter() == -1) continue; + parP = tmp[i]; + } + // Reset daughter information + + Int_t idaumin = parP->GetFirstDaughter() - fHgwmk - 1; + Int_t idaumax = parP->GetLastDaughter() - fHgwmk - 1; + parP->SetFirstDaughter(-1); + parP->SetLastDaughter(-1); + for (j = idaumin; j <= idaumax; j++) { + // Skip deleted particles + if (!tmp[j]) continue; + // Skip particles already handled + if (map1[j] != -99) continue; + Int_t jpa = tmp[j]->GetFirstMother(); + // Check if daughter of current parent + if (jpa == ipa) { + particles[fLoadPoint] = tmp[j]; + // Re-establish daughter information + parP->SetLastDaughter(fLoadPoint); + if (parP->GetFirstDaughter() == -1) parP->SetFirstDaughter(fLoadPoint); + // Set Mother information + if (i != -1) { + tmp[j]->SetFirstMother(map1[i]); + } + // Build the map + map1[j] = fLoadPoint; + // Increase load point + fLoadPoint++; + } + } // children + } // parents + + delete[] tmp; + + // + // Build map for remapping of hits + // + TArrayI map(fNtrack); + for (i = 0; i < fNtrack; i ++) { + if (i <= fHgwmk) { + map[i] = i; + } else{ + map[i] = map1[i - fHgwmk -1]; + } + } + + // Now loop on all registered hit lists + + TList* hitLists = gAlice->GetMCApp()->GetHitLists(); + TIter next(hitLists); + TCollection *hitList; + + while((hitList = dynamic_cast(next()))) { + TIter nexthit(hitList); + AliHit *hit; + while((hit = dynamic_cast(nexthit()))) { + hit->SetTrack(map[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 *detector; + while((detector = dynamic_cast(nextmod()))) { + detector->RemapTrackHitIDs(map.GetArray()); + detector->RemapTrackReferencesIDs(map.GetArray()); + } + gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray()); + } // new particles poduced } Bool_t AliStack::KeepPhysics(TParticle* part) @@ -470,7 +611,7 @@ void AliStack::FinishEvent() // When all primaries were filled no particle!=0 // should be left => to be removed later. - if (allFilled) printf("Why != 0 part # %d?\n",i); + if (allFilled) AliWarning(Form("Why != 0 part # %d?\n",i)); } else { @@ -557,10 +698,10 @@ void AliStack::SetHighWaterMark(Int_t) // // Set high water mark for last track in event // - + + fHgwmk = fNtrack-1; fCurrentPrimary=fHgwmk; - // Set also number of primary tracks fNprimary = fHgwmk+1; fNtrack = fHgwmk+1; @@ -582,9 +723,9 @@ TParticle* AliStack::Particle(Int_t i) // and the fParticleFileMap[i] give the same; // give the fatal error if not if (entry != fParticleFileMap[i]) { - Fatal("Particle", + AliFatal(Form( "!! The algorithmic way and map are different: !!\n entry: %d map: %d", - entry, fParticleFileMap[i]); + entry, fParticleFileMap[i])); } TreeK()->GetEntry(entry); @@ -634,7 +775,7 @@ Int_t AliStack::GetCurrentParentTrackNumber() const if (current) return current->GetFirstMother(); else { - Warning("GetCurrentParentTrackNumber", "Current track not found in the stack"); + AliWarning("Current track not found in the stack"); return -1; } } @@ -778,6 +919,8 @@ TParticle* AliStack::GetNextParticle() // nothing to be tracked fCurrent = -1; + + return particle; } //__________________________________________________________________________________________ @@ -794,7 +937,7 @@ TTree* AliStack::TreeK() AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName); if (rl == 0x0) { - Fatal("TreeK","Can not get RunLoader from event folder named %s",fEventFolderName.Data()); + AliFatal(Form("Can not get RunLoader from event folder named %s",fEventFolderName.Data())); return 0x0;//pro forma } fTreeK = rl->TreeK(); @@ -805,10 +948,7 @@ TTree* AliStack::TreeK() else { //don't panic - could be Lego - if (AliLoader::GetDebug()) - { - Warning("TreeK","Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data()); - } + AliWarning(Form("Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data())); } } return fTreeK;//never reached @@ -820,12 +960,12 @@ void AliStack::ConnectTree() // // Creates branch for writing particles // - if (AliLoader::GetDebug()) Info("ConnectTree","Connecting TreeK"); + AliDebug(1, "Connecting TreeK"); if (fTreeK == 0x0) { if (TreeK() == 0x0) { - Fatal("ConnectTree","Parameter is NULL");//we don't like such a jokes + AliFatal("Parameter is NULL");//we don't like such a jokes return; } return;//in this case TreeK() calls back this method (ConnectTree) @@ -835,35 +975,32 @@ void AliStack::ConnectTree() // Create a branch for particles - if (AliLoader::GetDebug()) - Info("ConnectTree","Tree name is %s",fTreeK->GetName()); + AliDebug(2, Form("Tree name is %s",fTreeK->GetName())); if (fTreeK->GetDirectory()) { - if (AliLoader::GetDebug()) - Info("ConnectTree","and dir is %s",fTreeK->GetDirectory()->GetName()); + AliDebug(2, Form("and dir is %s",fTreeK->GetDirectory()->GetName())); } else - Warning("ConnectTree","DIR IS NOT SET !!!"); + AliWarning("DIR IS NOT SET !!!"); - TBranch *branch=fTreeK->GetBranch(AliRunLoader::fgkKineBranchName); + TBranch *branch=fTreeK->GetBranch(AliRunLoader::GetKineBranchName()); if(branch == 0x0) { - branch = fTreeK->Branch(AliRunLoader::fgkKineBranchName, "TParticle", &fParticleBuffer, 4000); - if (AliLoader::GetDebug()) Info("ConnectTree","Creating Branch in Tree"); + branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &fParticleBuffer, 4000); + AliDebug(2, "Creating Branch in Tree"); } else { - if (AliLoader::GetDebug()) Info("ConnectTree","Branch Found in Tree"); + AliDebug(2, "Branch Found in Tree"); branch->SetAddress(&fParticleBuffer); } if (branch->GetDirectory()) { - if (AliLoader::GetDebug()) - Info("ConnectTree","Branch Dir Name is %s",branch->GetDirectory()->GetName()); + AliDebug(1, Form("Branch Dir Name is %s",branch->GetDirectory()->GetName())); } else - Warning("ConnectTree","Branch Dir is NOT SET"); + AliWarning("Branch Dir is NOT SET"); } //__________________________________________________________________________________________ @@ -891,7 +1028,7 @@ Bool_t AliStack::GetEvent() if (TreeK() == 0x0) //forces connecting { - Error("GetEvent","cannot find Kine Tree for current event\n"); + AliError("cannot find Kine Tree for current event"); return kFALSE; }