]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliStack.cxx
Changes needed on Mac (F.Carminati)
[u/mrichter/AliRoot.git] / STEER / AliStack.cxx
index 8b5128b5d64dfe961b25bbbe509905817bf4f85c..fbf97c7e11875caff8881c2e844b0a6926602daf 100644 (file)
@@ -31,6 +31,7 @@
 #include <TParticlePDG.h>
 #include <TTree.h>
 
+#include "AliLog.h"
 #include "AliHit.h"
 #include "AliModule.h"
 #include "AliRun.h"
@@ -104,7 +105,7 @@ AliStack::AliStack(const AliStack& st):
 //_______________________________________________________________________
 void AliStack::Copy(TObject&) const
 {
-  Fatal("Copy","Not implemented!\n");
+  AliFatal("Not implemented!");
 }
 
 //_______________________________________________________________________
@@ -169,8 +170,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 !");
     }
 }
 
@@ -205,7 +206,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++]) 
@@ -221,6 +223,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) {
@@ -230,7 +234,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 { 
@@ -296,11 +300,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<fNtrack; i++) {
     // Preset map, to be removed later
     if(i<=fHgwmk) map[i]=i ; 
@@ -365,11 +369,12 @@ void AliStack::PurifyKine()
   TList* hitLists = gAlice->GetMCApp()->GetHitLists();
   TIter next(hitLists);
   TCollection *hitList;
+  
   while((hitList = dynamic_cast<TCollection*>(next()))) {
     TIter nexthit(hitList);
     AliHit *hit;
     while((hit = dynamic_cast<AliHit*>(nexthit()))) {
-      hit->SetTrack(map[hit->GetTrack()]);
+       hit->SetTrack(map[hit->GetTrack()]);
     }
   }
 
@@ -402,11 +407,143 @@ void AliStack::PurifyKine()
 
    Int_t toshrink = fNtrack-fHgwmk-1;
    fLoadPoint-=toshrink;
+
+   
    for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(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<TParticle*>(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<TCollection*>(next()))) {
+         TIter nexthit(hitList);
+         AliHit *hit;
+         while((hit = dynamic_cast<AliHit*>(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<AliModule*>(nextmod()))) {
+         detector->RemapTrackHitIDs(map.GetArray());
+         detector->RemapTrackReferencesIDs(map.GetArray());
+      }
+      gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray());
+  } // new particles poduced
 }
 
 Bool_t AliStack::KeepPhysics(TParticle* part)
@@ -472,7 +609,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 
      {
@@ -559,10 +696,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;      
@@ -584,9 +721,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);
@@ -636,7 +773,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;
   }  
 }
@@ -780,6 +917,8 @@ TParticle* AliStack::GetNextParticle()
   
   // nothing to be tracked
   fCurrent = -1;
+  
   return particle;  
 }
 //__________________________________________________________________________________________
@@ -796,7 +935,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();
@@ -807,10 +946,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
@@ -822,12 +958,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) 
@@ -837,35 +973,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::GetKineBranchName());
   if(branch == 0x0)
    {
     branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &fParticleBuffer, 4000);
-    if (AliLoader::GetDebug()) Info("ConnectTree","Creating Branch in Tree");
+    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");
 }
 //__________________________________________________________________________________________
 
@@ -893,7 +1026,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;
      }