X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStack.cxx;h=2fc8d5a25a7cd0048cac4ae60326339a8765bca1;hb=594d89909ea15cd9abc17048a25a153436b409a2;hp=4d1040f7df396564f92723db7003bf19a4f5f2cd;hpb=fe046ade3136e16900c79722861d96d4594459e0;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx index 4d1040f7df3..2fc8d5a25a7 100644 --- a/STEER/AliStack.cxx +++ b/STEER/AliStack.cxx @@ -17,24 +17,26 @@ /////////////////////////////////////////////////////////////////////////////// // // -// 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 -#include -#include "AliStack.h" -#include "AliRun.h" -#include "AliModule.h" #include "AliHit.h" +#include "AliModule.h" +#include "AliRun.h" +#include "AliMC.h" +#include "AliRunLoader.h" +#include "AliStack.h" ClassImp(AliStack) @@ -50,7 +52,8 @@ AliStack::AliStack(): fCurrent(-1), fCurrentPrimary(-1), fHgwmk(0), - fLoadPoint(0) + fLoadPoint(0), + fEventFolderName(AliConfig::GetDefaultEventFolderName()) { // // Default constructor @@ -58,7 +61,7 @@ AliStack::AliStack(): } //_______________________________________________________________________ -AliStack::AliStack(Int_t size): +AliStack::AliStack(Int_t size, const char* evfoldname): fParticles(new TClonesArray("TParticle",1000)), fParticleMap(new TObjArray(size)), fParticleFileMap(0), @@ -69,7 +72,8 @@ AliStack::AliStack(Int_t size): fCurrent(-1), fCurrentPrimary(-1), fHgwmk(0), - fLoadPoint(0) + fLoadPoint(0), + fEventFolderName(evfoldname) { // // Constructor @@ -98,7 +102,7 @@ AliStack::AliStack(const AliStack& st): } //_______________________________________________________________________ -void AliStack::Copy(AliStack&) const +void AliStack::Copy(TObject&) const { Fatal("Copy","Not implemented!\n"); } @@ -115,7 +119,7 @@ AliStack::~AliStack() delete fParticles; } delete fParticleMap; - delete fTreeK; + //PH??? delete fTreeK; } // @@ -123,9 +127,9 @@ AliStack::~AliStack() // //_____________________________________________________________________________ -void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, - Float_t *vpos, Float_t *polar, Float_t tof, - TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) +void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, + Float_t *vpos, Float_t *polar, Float_t tof, + TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) { // // Load a track on the stack @@ -161,17 +165,17 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, // mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS); - SetTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e, + PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e, vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2], mech, ntr, weight, is); } else { - Warning("SetTrack", "Particle type %d not defined in PDG Database !\n", pdg); - Warning("SetTrack", "Particle skipped !\n"); + Warning("PushTrack", "Particle type %d not defined in PDG Database !\n", pdg); + Warning("PushTrack", "Particle skipped !\n"); } } //_____________________________________________________________________________ -void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, +void AliStack::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, @@ -201,7 +205,8 @@ void AliStack::SetTrack(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++]) @@ -217,6 +222,8 @@ void AliStack::SetTrack(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) { @@ -226,7 +233,7 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); } else { - printf("Error in AliStack::SetTrack: Parent %d does not exist\n",parent); + printf("Error in AliStack::PushTrack: Parent %d does not exist\n",parent); } } else { @@ -242,7 +249,7 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, } //_____________________________________________________________________________ -TParticle* AliStack::GetNextTrack(Int_t& itrack) +TParticle* AliStack::PopNextTrack(Int_t& itrack) { // // Returns next track from stack of particles @@ -262,58 +269,8 @@ TParticle* AliStack::GetNextTrack(Int_t& itrack) return track; } -/* -//_____________________________________________________________________________ -void AliStack::GetNextTrack(Int_t& itrack, 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) -{ - // - // Return next track from stack of particles - // - - - TParticle* track = GetNextParticle(); -// cout << "GetNextTrack():" << fCurrent << fNprimary << endl; - - if (track) { - itrack = fCurrent; - pdg = track->GetPdgCode(); - px = track->Px(); - py = track->Py(); - pz = track->Pz(); - e = track->Energy(); - vx = track->Vx(); - vy = track->Vy(); - vz = track->Vz(); - tof = track->T(); - TVector3 pol; - track->GetPolarisation(pol); - polx = pol.X(); - poly = pol.Y(); - polz = pol.Z(); - track->SetBit(kDoneBit); -// cout << "Filled params" << endl; - } - else - itrack = -1; - - // - // stop and start timer when we start a primary track - Int_t nprimaries = fNprimary; - if (fCurrent >= nprimaries) return; - if (fCurrent < nprimaries-1) { - fTimer.Stop(); - track=(TParticle*) fParticleMap->At(fCurrent+1); - // track->SetProcessTime(fTimer.CpuTime()); - } - fTimer.Start(); -} - -*/ //_____________________________________________________________________________ -TParticle* AliStack::GetPrimaryForTracking(Int_t i) +TParticle* AliStack::PopPrimaryForTracking(Int_t i) { // // Returns i-th primary particle if it is flagged to be tracked, @@ -328,7 +285,6 @@ TParticle* AliStack::GetPrimaryForTracking(Int_t i) return 0; } - //_____________________________________________________________________________ void AliStack::PurifyKine() { @@ -343,11 +299,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; i(particles.At(i)))) { // // Check of this track should be kept for physics reasons - if (KeepPhysics(part)) KeepTrack(i); + if (KeepPhysics(part)) KeepTrack(i); // part->ResetBit(kDaughtersBit); part->SetFirstDaughter(-1); @@ -409,14 +365,15 @@ void AliStack::PurifyKine() } // Now loop on all registered hit lists - TList* hitLists = gAlice->GetHitLists(); + 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()]); + hit->SetTrack(map[hit->GetTrack()]); } } @@ -424,7 +381,7 @@ void AliStack::PurifyKine() // 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; @@ -432,29 +389,162 @@ void AliStack::PurifyKine() detector->RemapTrackHitIDs(map.GetArray()); detector->RemapTrackReferencesIDs(map.GetArray()); } - + // + gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray()); + // Now the output bit, from fHgwmk to nkeep we write everything and we erase if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5)); for (i=fHgwmk+1; i(particles.At(i)); - fParticleFileMap[i]=static_cast(fTreeK->GetEntries()); - fTreeK->Fill(); + fParticleFileMap[i]=static_cast(TreeK()->GetEntries()); + TreeK()->Fill(); particles[i]=fParticleBuffer=0; - } + } for (i=nkeep; 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))->Clone(); + tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i)); + +// if (((TParticle*) (particles.At(fHgwmk + 1 + i)))->TestBit(kKeepBit)) +// tmp[i]->SetBit(kKeepBit); +// if (((TParticle*) (particles.At(fHgwmk + 1 + i)))->TestBit(kDoneBit)) +// tmp[i]->SetBit(kDoneBit); + } 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 + parP->SetFirstDaughter(-1); + parP->SetLastDaughter(-1); + for (j = i + 1; j < nNew; 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) { // @@ -486,22 +576,21 @@ Bool_t AliStack::KeepPhysics(TParticle* part) //_____________________________________________________________________________ void AliStack::FinishEvent() { - // - // Write out the kinematics that was not yet filled - // +// +// Write out the kinematics that was not yet filled +// - // Update event header +// Update event header - if (!fTreeK) { + if (!TreeK()) { // Fatal("FinishEvent", "No kinematics tree is defined."); // Don't panic this is a probably a lego run return; - } CleanParents(); - if(fTreeK->GetEntries() ==0) { + if(TreeK()->GetEntries() ==0) { // set the fParticleFileMap size for the first time fParticleFileMap.Set(fHgwmk+1); } @@ -511,8 +600,8 @@ void AliStack::FinishEvent() for(Int_t i=0; iAt(i))) { fParticleBuffer = dynamic_cast(part); - fParticleFileMap[i]= static_cast(fTreeK->GetEntries()); - fTreeK->Fill(); + fParticleFileMap[i]= static_cast(TreeK()->GetEntries()); + TreeK()->Fill(); //PH (*fParticleMap)[i]=fParticleBuffer=0; fParticleBuffer=0; fParticleMap->AddAt(0,i); @@ -521,19 +610,18 @@ void AliStack::FinishEvent() // should be left => to be removed later. if (allFilled) printf("Why != 0 part # %d?\n",i); } - else { + else + { // // printf("Why = 0 part # %d?\n",i); => We know. // break; - // we don't break now in order to be sure there is no - // particle !=0 left. - // To be removed later and replaced with break. - if(!allFilled) allFilled = kTRUE; - } -// cout << "Nof particles: " << fNtrack << endl; - //Reset(); + // we don't break now in order to be sure there is no + // particle !=0 left. + // To be removed later and replaced with break. + if(!allFilled) allFilled = kTRUE; + } } - //_____________________________________________________________________________ + void AliStack::FlagTrack(Int_t track) { // @@ -573,12 +661,13 @@ void AliStack::Reset(Int_t size) // // Resets stack // - + fNtrack=0; fNprimary=0; fHgwmk=0; fLoadPoint=0; fCurrent = -1; + fTreeK = 0x0; ResetArrays(size); } @@ -606,10 +695,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; @@ -636,7 +725,7 @@ TParticle* AliStack::Particle(Int_t i) entry, fParticleFileMap[i]); } - fTreeK->GetEntry(entry); + TreeK()->GetEntry(entry); new ((*fParticles)[nentries]) TParticle(*fParticleBuffer); fParticleMap->AddAt((*fParticles)[nentries],i); } @@ -672,7 +761,7 @@ Int_t AliStack::TreeKEntry(Int_t id) const } //_____________________________________________________________________________ -Int_t AliStack::CurrentTrackParent() const +Int_t AliStack::GetCurrentParentTrackNumber() const { // // Return number of the parent of the current track @@ -683,7 +772,7 @@ Int_t AliStack::CurrentTrackParent() const if (current) return current->GetFirstMother(); else { - Warning("CurrentTrackParent", "Current track not found in the stack"); + Warning("GetCurrentParentTrackNumber", "Current track not found in the stack"); return -1; } } @@ -725,8 +814,7 @@ void AliStack::DumpPStack () Int_t i; - printf( - "\n\n=======================================================================\n"); + printf("\n\n=======================================================================\n"); for (i=0;i(fParticleMap->At(i)); if ((particle) && (!particle->TestBit(kDoneBit))) { fCurrent=i; - //cout << "GetNextParticle() - secondary " - // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; return particle; } } @@ -825,88 +910,139 @@ TParticle* AliStack::GetNextParticle() fCurrent = fCurrentPrimary; particle = dynamic_cast(fParticleMap->At(fCurrentPrimary--)); if ((particle) && (!particle->TestBit(kDoneBit))) { - //cout << "GetNextParticle() - primary " - // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; return particle; } } // nothing to be tracked fCurrent = -1; - //cout << "GetNextParticle() - none " - // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; + + return particle; } +//__________________________________________________________________________________________ +TTree* AliStack::TreeK() +{ +//returns TreeK + if (fTreeK) + { + return fTreeK; + } + else + { + AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName); + if (rl == 0x0) + { + Fatal("TreeK","Can not get RunLoader from event folder named %s",fEventFolderName.Data()); + return 0x0;//pro forma + } + fTreeK = rl->TreeK(); + if ( fTreeK ) + { + ConnectTree(); + } + else + { + //don't panic - could be Lego + if (AliLoader::GetDebug()) + { + Warning("TreeK","Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data()); + } + } + } + return fTreeK;//never reached +} //__________________________________________________________________________________________ -void AliStack::MakeTree(Int_t event, const char * /*file*/) + +void AliStack::ConnectTree() { // -// Make Kine tree and creates branch for writing particles -// - TBranch *branch=0; - // Make Kinematics Tree - char hname[30]; - if (!fTreeK) { - sprintf(hname,"TreeK%d",event); - fTreeK = new TTree(hname,"Kinematics"); - // Create a branch for particles - branch = fTreeK->Branch("Particles", "TParticle", &fParticleBuffer, 4000); - fTreeK->Write(0,TObject::kOverwrite); - } +// Creates branch for writing particles +// + if (AliLoader::GetDebug()) Info("ConnectTree","Connecting TreeK"); + if (fTreeK == 0x0) + { + if (TreeK() == 0x0) + { + Fatal("ConnectTree","Parameter is NULL");//we don't like such a jokes + return; + } + return;//in this case TreeK() calls back this method (ConnectTree) + //tree after setting fTreeK, the rest was already executed + //it is safe to return now + } + + // Create a branch for particles + + if (AliLoader::GetDebug()) + Info("ConnectTree","Tree name is %s",fTreeK->GetName()); + + if (fTreeK->GetDirectory()) + { + if (AliLoader::GetDebug()) + Info("ConnectTree","and dir is %s",fTreeK->GetDirectory()->GetName()); + } + else + Warning("ConnectTree","DIR IS NOT SET !!!"); + + TBranch *branch=fTreeK->GetBranch(AliRunLoader::GetKineBranchName()); + if(branch == 0x0) + { + branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &fParticleBuffer, 4000); + if (AliLoader::GetDebug()) Info("ConnectTree","Creating Branch in Tree"); + } + else + { + if (AliLoader::GetDebug()) Info("ConnectTree","Branch Found in Tree"); + branch->SetAddress(&fParticleBuffer); + } + if (branch->GetDirectory()) + { + if (AliLoader::GetDebug()) + Info("ConnectTree","Branch Dir Name is %s",branch->GetDirectory()->GetName()); + } + else + Warning("ConnectTree","Branch Dir is NOT SET"); } +//__________________________________________________________________________________________ -//_____________________________________________________________________________ -void AliStack::BeginEvent(Int_t event) + +void AliStack::BeginEvent() { // start a new event -// -// - fNprimary = 0; - fNtrack = 0; - - char hname[30]; - if(fTreeK) { - fTreeK->Reset(); - sprintf(hname,"TreeK%d",event); - fTreeK->SetName(hname); - } + Reset(); } //_____________________________________________________________________________ void AliStack::FinishRun() { // Clean TreeK information - if (fTreeK) { - delete fTreeK; fTreeK = 0; - } } - //_____________________________________________________________________________ -Bool_t AliStack::GetEvent(Int_t event) + +Bool_t AliStack::GetEvent() { // // Get new event from TreeK // Reset/Create the particle stack - if (fTreeK) delete fTreeK; - - // Get Kine Tree from file - char treeName[20]; - sprintf(treeName,"TreeK%d",event); - fTreeK = dynamic_cast(gDirectory->Get(treeName)); + fTreeK = 0x0; - if (fTreeK) - fTreeK->SetBranchAddress("Particles", &fParticleBuffer); - else { - // Error("GetEvent","cannot find Kine Tree for event:%d\n",event); - Warning("GetEvent","cannot find Kine Tree for event:%d\n",event); + if (TreeK() == 0x0) //forces connecting + { + Error("GetEvent","cannot find Kine Tree for current event\n"); return kFALSE; - } -// printf("\n primaries %d", fNprimary); -// printf("\n tracks %d", fNtrack); + } - Int_t size = static_cast(fTreeK->GetEntries()); + Int_t size = (Int_t)TreeK()->GetEntries(); ResetArrays(size); return kTRUE; } +//_____________________________________________________________________________ + +void AliStack::SetEventFolderName(const char* foldname) +{ + //Sets event folder name + fEventFolderName = foldname; +}