Non-buffering readers implemented, proper changes in analysis. Compiler warnings...
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Aug 2003 15:54:32 +0000 (15:54 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Aug 2003 15:54:32 +0000 (15:54 +0000)
31 files changed:
HBTAN/AliHBTAnalysis.cxx
HBTAN/AliHBTEvent.cxx
HBTAN/AliHBTEvent.h
HBTAN/AliHBTEventBuffer.cxx [new file with mode: 0644]
HBTAN/AliHBTEventBuffer.h [new file with mode: 0644]
HBTAN/AliHBTFunction.cxx
HBTAN/AliHBTPairCut.cxx
HBTAN/AliHBTParticle.cxx
HBTAN/AliHBTParticle.h
HBTAN/AliHBTParticleCut.cxx
HBTAN/AliHBTParticleCut.h
HBTAN/AliHBTPositionRandomizer.cxx
HBTAN/AliHBTPositionRandomizer.h
HBTAN/AliHBTReader.cxx
HBTAN/AliHBTReader.h
HBTAN/AliHBTReaderESD.cxx
HBTAN/AliHBTReaderESD.h
HBTAN/AliHBTReaderITSv2.cxx
HBTAN/AliHBTReaderITSv2.h
HBTAN/AliHBTReaderInternal.cxx
HBTAN/AliHBTReaderInternal.h
HBTAN/AliHBTReaderKineTree.cxx
HBTAN/AliHBTReaderKineTree.h
HBTAN/AliHBTReaderTPC.cxx
HBTAN/AliHBTReaderTPC.h
HBTAN/AliHBTRun.cxx
HBTAN/AliHBTRun.h
HBTAN/AliHBTTwoTrackEffFctn.h
HBTAN/AliHBTWriteInternFormat.C
HBTAN/HBTAnalysisLinkDef.h
HBTAN/libHBTAN.pkg

index 22449d4a99131e63b8e6b470d37f9db331e5557e..943a40d67a5e81fc486070708a81a0cbbd18280d 100644 (file)
@@ -8,6 +8,7 @@
 #include "AliHBTPairCut.h"
 #include "AliHBTFunction.h"
 #include "AliHBTMonitorFunction.h"
+#include "AliHBTEventBuffer.h"
  
 #include <TBenchmark.h>
 #include <TList.h>
@@ -221,26 +222,10 @@ void AliHBTAnalysis::Process(Option_t* option)
  
  Bool_t nonid = IsNonIdentAnalysis();
  
+ Init();
  if(oT && oP)
   { 
-    if (fReader->GetNumberOfPartEvents() <1)
-     {
-       Error("Process","There is no Particles. Maybe change the option?");
-       return;
-     }
-    if (fReader->GetNumberOfTrackEvents() <1)
-     {
-       Error("Process","There is no Tracks. Maybe change the option?");
-       return;
-     }
-    
-    if ( RunCoherencyCheck() )
-      {
-        Error("Process",
-              "Coherency check not passed. Maybe change the option?\n");
-        return;
-      }
-    Init();
     if (nonid) ProcessTracksAndParticlesNonIdentAnal();
     else ProcessTracksAndParticles();
     return;
@@ -248,12 +233,6 @@ void AliHBTAnalysis::Process(Option_t* option)
 
  if(oT)
   {
-    if (fReader->GetNumberOfTrackEvents() <1)
-     {
-       Error("Process","There is no data to analyze.");
-       return;
-     }
-    Init();
     if (nonid) ProcessTracksNonIdentAnal();
     else ProcessTracks();
     return;
@@ -261,12 +240,6 @@ void AliHBTAnalysis::Process(Option_t* option)
  
  if(oP)
   {
-    if (fReader->GetNumberOfPartEvents() <1)
-     {
-       Error("Process","There is no data to analyze.");
-       return;
-     }
-    Init();
     if (nonid) ProcessParticlesNonIdentAnal();
     else ProcessParticles();
     return;
@@ -285,45 +258,96 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
   AliHBTParticle * track1, * track2;
   
   AliHBTEvent * trackEvent, *partEvent;
+  AliHBTEvent * trackEvent1 = 0x0,*partEvent1 = 0x0;
   AliHBTEvent * trackEvent2,*partEvent2;
   
 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
   
-  Int_t nev = fReader->GetNumberOfTrackEvents();
-  
-  /***************************************/
-  /******   Looping same events   ********/
-  /******   filling numerators    ********/
-  /***************************************/
+//  Int_t nev = fReader->GetNumberOfTrackEvents();
   AliHBTPair * trackpair = new AliHBTPair();
   AliHBTPair * partpair = new AliHBTPair();
 
   AliHBTPair * tmptrackpair;//temprary pointers to pairs
   AliHBTPair * tmppartpair;
 
+  AliHBTEventBuffer partbuffer(fBufferSize);
+  AliHBTEventBuffer trackbuffer(fBufferSize);
+  
+  
   register UInt_t ii;
   
-  for (Int_t i = 0;i<nev;i++)
+  Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
+  
+  fReader->Rewind();
+  Int_t i = -1;
+  while (fReader->Next() == kFALSE)
     {
-      partEvent= fReader->GetParticleEvent(i);
-      trackEvent = fReader->GetTrackEvent(i);
+      i++;
+      partEvent= fReader->GetParticleEvent();
+      trackEvent = fReader->GetTrackEvent();
       
-      if (!partEvent) continue;
-      
-      //N = 0;
+      if ( !partEvent || !trackEvent ) 
+       {
+         Error("ProcessTracksAndParticles","Can not get event");
+         return;
+       }
+       
+      if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
+       {
+         Fatal("ProcessTracksAndParticles",
+               "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
+               i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
+         return;
+       }
       
-      for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+      if(partEvent1 == 0x0) 
        {
+         partEvent1 = new AliHBTEvent();
+         partEvent1->SetOwner(kTRUE);
 
-         if ( (j%fDisplayMixingInfo) == 0) 
+         trackEvent1 = new AliHBTEvent();
+         trackEvent1->SetOwner(kTRUE);
+       }
+      else
+       {
+         partEvent1->Reset();
+         trackEvent1->Reset();
+       }
+       
+      for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+       {
+         /***************************************/
+         /******   Looping same events   ********/
+         /******   filling numerators    ********/
+         /***************************************/
+         if ( (j%fDisplayMixingInfo) == 0)
             Info("ProcessTracksAndParticles",
                  "Mixing particle %d from event %d with particles from event %d",j,i,i);
 
          part1= partEvent->GetParticle(j);
-         track1= trackEvent->GetParticle(j);       
-
-         if (fPairCut->GetFirstPartCut()->Pass(part1)) continue;
+         track1= trackEvent->GetParticle(j);
+         
+         if( part1->GetPdgCode() != track1->GetPdgCode() )
+          {
+            Fatal("ProcessTracksAndParticles",
+                  "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
+                  i,j, part1->GetPdgCode(),track1->GetPdgCode() );
+            return;
+          }
+         
+         Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
+         
+         if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
+          {
+            //accepted by any cut
+            // we have to copy because reader keeps only one event
+            
+            partEvent1->AddParticle(new AliHBTParticle(*part1));
+            trackEvent1->AddParticle(new AliHBTParticle(*track1));
+          }
 
+         if (firstcut) continue;
+         
          for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
            fParticleMonitorFunctions[ii]->Process(part1);
          for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
@@ -331,8 +355,7 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
          for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
            fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
 
-         if ( (fNParticleFunctions == 0) && (fNTrackFunctions ==0) && (fNParticleAndTrackFunctions == 0))
-           continue; 
+         if (nocorrfctns) continue;
         
          for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
           {
@@ -343,12 +366,12 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
             track2= trackEvent->GetParticle(k);       
             trackpair->SetParticles(track1,track2);
 
-            if(fPairCut->Pass(partpair) ) //check pair cut
-              { //do not meets crietria of the pair cut, try with swaped pairs
+            if(fPairCut->Pass(partpair) ) //check pair cut 
+              { //do not meets crietria of the pair cut, try with swapped pairs
                 if( fPairCut->Pass(partpair->GetSwapedPair()) )
                   continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
                 else 
-               { //swaped pair meets all the criteria
+                 { //swaped pair meets all the criteria
                    tmppartpair = partpair->GetSwapedPair();
                    tmptrackpair = trackpair->GetSwapedPair();
                  }
@@ -366,58 +389,27 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
                  
             for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
                    fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
-           }
-       }
-    }
-
-
-  /***** Filling denominators    *********/
-  /***************************************/
-  if (fBufferSize != 0)
-   for (Int_t i = 0;i<nev-1;i++)   //In each event (but last) ....
-    {
-  
-      if ((fNParticleFunctions == 0) && (fNTrackFunctions ==0) && (fNParticleAndTrackFunctions == 0))
-        continue;  
-
-      partEvent= fReader->GetParticleEvent(i);
-      if (!partEvent) continue;
-      
-      trackEvent = fReader->GetTrackEvent(i); 
-      
-      for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
-       {
-           part1= partEvent->GetParticle(j);
-           track1= trackEvent->GetParticle(j);
-
-           if (fPairCut->GetFirstPartCut()->Pass(part1)) continue;
-           Int_t maxeventnumber;
-           
-           if ( ((i+fBufferSize) >= nev) ||( fBufferSize < 0) ) //if buffer size is negative 
-                                                                //or current event+buffersize is greater
-                                                                //than max nuber of events
-             {
-               maxeventnumber = nev; //set the max event number 
-             }
-           else 
-             {
-               maxeventnumber = i+fBufferSize; //set the current event number + fBufferSize
-             }
-           for (Int_t k = i+1; k<maxeventnumber;k++)  // ... Loop over next event
-            {
-             
-             partEvent2= fReader->GetParticleEvent(k);
-             if (!partEvent2) continue;
-              
-             trackEvent2 = fReader->GetTrackEvent(k);
-             
-             if ( (j%fDisplayMixingInfo) == 0) 
-                 Info("ProcessTracksAndParticles",
-                      "Mixing particle %d from event %d with particles from event %d",j,i,k);
+           //end of 2nd loop over particles from the same event  
+          }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+          
+         /***************************************/
+         /***** Filling denominators    *********/
+         /***************************************/
+         if (fBufferSize == 0) continue;
+         
+         partbuffer.ResetIter();
+         trackbuffer.ResetIter();
+         Int_t m = 0;
+         while (( partEvent2 = partbuffer.Next() ))
+          {
+            trackEvent2 = trackbuffer.Next();
             
-             for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
+            m++;
+            if ( (j%fDisplayMixingInfo) == 0)
+               Info("ProcessTracksAndParticles",
+                    "Mixing particle %d from event %d with particles from event %d",j,i,i-m);
+         
+            for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
               {
                 part2= partEvent2->GetParticle(l);
                 partpair->SetParticles(part1,part2);
@@ -449,10 +441,14 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
                 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
                   fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
-            }//for (Int_t k = i+1; k<maxeventnumber;k++)  // ... Loop over next event
-       }
-    } 
-  /***************************************/
+
+          }
+        //end of loop over particles from first event
+       }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+      partEvent1  = partbuffer.Push(partEvent1);
+      trackEvent1 = trackbuffer.Push(trackEvent1);
+     //end of loop over events  
+    }//while (fReader->Next() == kFALSE)
 } 
 /*************************************************************************************/
 
@@ -462,37 +458,60 @@ void AliHBTAnalysis::ProcessTracks()
 //the loops are splited
   AliHBTParticle * track1, * track2;
   AliHBTEvent * trackEvent;
+  AliHBTEvent * trackEvent1 = 0x0;
   AliHBTEvent * trackEvent2;
 
-  UInt_t ii;
-  Int_t nev = fReader->GetNumberOfTrackEvents();
+  register UInt_t ii;
 
   AliHBTPair * trackpair = new AliHBTPair();
   AliHBTPair * tmptrackpair; //temporary pointer 
   
-  /***************************************/
-  /******   Looping same events   ********/
-  /******   filling numerators    ********/
-  /***************************************/
-  for (Int_t i = 0;i<nev;i++)
+  AliHBTEventBuffer trackbuffer(fBufferSize);
+  
+  fReader->Rewind();
+  Int_t i = -1;
+  while (fReader->Next() == kFALSE)
     {
-      trackEvent = fReader->GetTrackEvent(i);
+      i++;
+      trackEvent = fReader->GetTrackEvent();
       if (!trackEvent) continue;
       
+      if(trackEvent1 == 0x0) 
+       {
+         trackEvent1 = new AliHBTEvent();
+         trackEvent1->SetOwner(kTRUE);
+       }
+      else
+       {
+         trackEvent1->Reset();
+       }
+      
       for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
        {
+         /***************************************/
+         /******   Looping same events   ********/
+         /******   filling numerators    ********/
+         /***************************************/
          if ( (j%fDisplayMixingInfo) == 0) 
                Info("ProcessTracks",
                     "Mixing particle %d from event %d with particles from event %d",j,i,i);
          
-         track1= trackEvent->GetParticle(j);       
-         if (fPairCut->GetFirstPartCut()->Pass(track1)) continue;
+         track1= trackEvent->GetParticle(j);
+         Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(track1);
+         
+         if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(track1) == kFALSE) )
+          {
+            //accepted by any cut
+            // we have to copy because reader keeps only one event
+            trackEvent1->AddParticle(new AliHBTParticle(*track1));
+          }
+
+         if (firstcut) continue;
 
          for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
            fTrackMonitorFunctions[ii]->Process(track1);
 
-         if ( fNTrackFunctions ==0 )
-           continue; 
+         if ( fNTrackFunctions ==0 ) continue; 
         
          for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
           {
@@ -512,70 +531,45 @@ void AliHBTAnalysis::ProcessTracks()
            for(ii = 0;ii<fNTrackFunctions;ii++)
                fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
           }
-       }
-    }
-
-  /***************************************/
-  /***** Filling diff histogram *********/
-  /***************************************/
-  if (fBufferSize != 0)
-   for (Int_t i = 0;i<nev-1;i++)   //In each event (but last) ....
-    {
-      if ( fNTrackFunctions ==0 )
-        continue; 
-
-      trackEvent = fReader->GetTrackEvent(i);
-      if (!trackEvent) continue;
-
-      for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
-       {
-         track1= trackEvent->GetParticle(j);
-         if (fPairCut->GetFirstPartCut()->Pass(track1)) continue;
-
-         Int_t maxeventnumber;
-           
-         if ( ((i+fBufferSize) >= nev) ||( fBufferSize < 0) ) //if buffer size is negative 
-                                                              //or current event+buffersize is greater
-                                                              //than max nuber of events
-          {
-            maxeventnumber = nev; //set the max event number 
-          }
-         else 
-          {
-            maxeventnumber = i+fBufferSize; //set the current event number + fBufferSize
-          }
-         for (Int_t k = i+1; k<maxeventnumber;k++)  // ... Loop over next event
+         /***************************************/
+         /***** Filling denominators    *********/
+         /***************************************/
+          
+         if (fBufferSize == 0) continue;
+         
+         trackbuffer.ResetIter();
+         
+         while (( trackEvent2 = trackbuffer.Next() ))
           {
-            trackEvent2 = fReader->GetTrackEvent(k);
-            if (!trackEvent2) continue;
-             
-            if ( (j%fDisplayMixingInfo) == 0) 
-                 Info("ProcessTracks",
-                      "Mixing particle %d from event %d with particles from event %d",j,i,k);
-            
             for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++)   //  ... on all particles
-             {
-               track2= trackEvent2->GetParticle(l);
-               trackpair->SetParticles(track1,track2);
-                
-               if(fPairCut->Pass(trackpair)) //check pair cut
-                { //do not meets crietria of the 
-                  if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
-                  else tmptrackpair = trackpair->GetSwapedPair();
-                }
-               else
-                {
+              {
+
+                track2= trackEvent2->GetParticle(l);
+                trackpair->SetParticles(track1,track2);
+
+                if( fPairCut->Pass(trackpair) ) //check pair cut
+                  { //do not meets crietria of the 
+                    if( fPairCut->Pass(trackpair->GetSwapedPair()) )
+                      continue;
+                    else 
+                     {
+                       tmptrackpair = trackpair->GetSwapedPair();
+                     }
+                  }
+                else
+                 {//meets criteria of the pair cut
                   tmptrackpair = trackpair;
-                }
-               for(ii = 0;ii<fNTrackFunctions;ii++)
-                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
-               
+                 }
+                 
+                for(ii = 0;ii<fNTrackFunctions;ii++)
+                  fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
+                 
              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
-          }//for (Int_t k = i+1; k<maxeventnumber;k++)  // ... Loop over next event
+            
+          }
        }
-    } 
-  /***************************************/
+      trackEvent1 = trackbuffer.Push(trackEvent1); 
+    }//while (fReader->Next() == kFALSE)
 }
 
 /*************************************************************************************/
@@ -585,131 +579,118 @@ void AliHBTAnalysis::ProcessParticles()
 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
 //the loops are splited
   AliHBTParticle * part1, * part2;
-  
   AliHBTEvent * partEvent;
+  AliHBTEvent * partEvent1 = 0x0;
   AliHBTEvent * partEvent2;
 
+  register UInt_t ii;
+
   AliHBTPair * partpair = new AliHBTPair();
-  AliHBTPair * tmppartpair; //temporary pointer to the pair
+  AliHBTPair * tmppartpair; //temporary pointer 
   
-  Int_t nev = fReader->GetNumberOfPartEvents();
+  AliHBTEventBuffer partbuffer(fBufferSize);
   
-  /***************************************/
-  /******   Looping same events   ********/
-  /******   filling numerators    ********/
-  /***************************************/
-  for (Int_t i = 0;i<nev;i++)
+  fReader->Rewind();
+  Int_t i = -1;
+  while (fReader->Next() == kFALSE)
     {
-      partEvent= fReader->GetParticleEvent(i);
+      i++;
+      partEvent = fReader->GetParticleEvent();
       if (!partEvent) continue;
-      //N = 0;
+      
+      if(partEvent1 == 0x0) 
+       {
+         partEvent1 = new AliHBTEvent();
+         partEvent1->SetOwner(kTRUE);
+       }
+      else
+       {
+         partEvent1->Reset();
+       }
       
       for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
        {
+         /***************************************/
+         /******   Looping same events   ********/
+         /******   filling numerators    ********/
+         /***************************************/
          if ( (j%fDisplayMixingInfo) == 0) 
-                 Info("ProcessParticles",
-                      "Mixing particle %d from event %d with particles from event %d",j,i,i);
-
+               Info("ProcessParts",
+                    "Mixing particle %d from event %d with particles from event %d",j,i,i);
+         
          part1= partEvent->GetParticle(j);
-         if (fPairCut->GetFirstPartCut()->Pass(part1)) continue;
-        
-         UInt_t zz;
-         for(zz = 0; zz<fNParticleMonitorFunctions; zz++)
-           fParticleMonitorFunctions[zz]->Process(part1);
-
-         if ( fNParticleFunctions ==0 )
-           continue; 
-
-         for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+         Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
+         
+         if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
           {
-            part2= partEvent->GetParticle(k);
-            if (part1->GetUID() == part2->GetUID()) continue; //this is the same particle but different incarnation (PID)
-            partpair->SetParticles(part1,part2);
-            
-            if( fPairCut->Pass(partpair) ) //check pair cut
-              { //do not meets crietria of the pair cut, try with swaped pairs
-                if(  fPairCut->Pass(partpair->GetSwapedPair() )  ) 
-                  continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
-                else 
-                 { //swaped pair meets all the criteria
-                   tmppartpair = partpair->GetSwapedPair();
-                 }
-              }
-            else
-              {
-                tmppartpair = partpair;
-              }
-            
-            UInt_t ii;
-            for(ii = 0;ii<fNParticleFunctions;ii++)
-                   fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
-           }
-       }
-    }
+            //accepted by any cut
+            // we have to copy because reader keeps only one event
+            partEvent1->AddParticle(new AliHBTParticle(*part1));
+          }
 
-  /***************************************/
-  /***** Filling diff histogram *********/
-  /***************************************/
-  if (fBufferSize != 0) //less then 0 mix everything, == 0 do not mix denominator
-   for (Int_t i = 0;i<nev-1;i++)   //In each event (but last)....
-    {
-      if ( fNParticleFunctions ==0 )
-        continue; 
+         if (firstcut) continue;
 
-      partEvent= fReader->GetParticleEvent(i);
-      if (!partEvent) continue;
+         for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+           fParticleMonitorFunctions[ii]->Process(part1);
 
-      for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
-       {
-           part1= partEvent->GetParticle(j);
-           if (fPairCut->GetFirstPartCut()->Pass(part1)) continue;
-           Int_t maxeventnumber;
+         if ( fNParticleFunctions == 0 ) continue; 
+        
+         for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+          {
+           part2= partEvent->GetParticle(k);
+           if (part1->GetUID() == part2->GetUID()) continue;
 
-           if ( ((i+fBufferSize) >= nev) ||( fBufferSize < 0) ) //if buffer size is negative 
-                                                                //or current event+buffersize is greater
-                                                                //than max nuber of events
-            {
-             maxeventnumber = nev; //set the max event number 
+           partpair->SetParticles(part1,part2);
+           if(fPairCut->Pass(partpair)) //check pair cut
+            { //do not meets crietria of the 
+              if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
+              else tmppartpair = partpair->GetSwapedPair();
             }
-           else 
+           else
             {
-             maxeventnumber = i+fBufferSize; //set the current event number + fBufferSize
+              tmppartpair = partpair;
             }
-           
-           for (Int_t k = i+1; k<maxeventnumber;k++)  // ... Loop over next event
-            {
-             
-             partEvent2= fReader->GetParticleEvent(k);
-             if (!partEvent2) continue;
-             
-             if ( (j%fDisplayMixingInfo) == 0) 
-                Info("ProcessParticles",
-                     "Mixing particle %d from event %d with particles from event %d",j,i,k);
-            
-             for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
+           for(ii = 0;ii<fNParticleFunctions;ii++)
+               fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
+          }
+         /***************************************/
+         /***** Filling denominators    *********/
+         /***************************************/
+          
+         if (fBufferSize == 0) continue;
+         
+         partbuffer.ResetIter();
+         
+         while (( partEvent2 = partbuffer.Next() ))
+          {
+            for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
               {
+
                 part2= partEvent2->GetParticle(l);
                 partpair->SetParticles(part1,part2);
-                
-                if(fPairCut->Pass(partpair)) //check pair cut
+
+                if( fPairCut->Pass(partpair) ) //check pair cut
                   { //do not meets crietria of the 
-                    if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
-                    else tmppartpair = partpair->GetSwapedPair();
+                    if( fPairCut->Pass(partpair->GetSwapedPair()) )
+                      continue;
+                    else 
+                     {
+                       tmppartpair = partpair->GetSwapedPair();
+                     }
                   }
                 else
-                  {
-                    tmppartpair = partpair;
-                  }
-                UInt_t ii;
+                 {//meets criteria of the pair cut
+                  tmppartpair = partpair;
+                 }
+                 
                 for(ii = 0;ii<fNParticleFunctions;ii++)
                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
-               
-              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
-            }//for (Int_t k = i+1; k<maxeventnumber;k++)  // ... Loop over next event
+                 
+             }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
+          }
        }
-    } 
-  /***************************************/
-
+      partEvent1 = partbuffer.Push(partEvent1); 
+    }//while (fReader->Next() == kFALSE)
 }
 /*************************************************************************************/
 
@@ -911,50 +892,56 @@ void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
   AliHBTEvent * trackEvent2=0x0,*partEvent2=0x0;
   AliHBTEvent * trackEvent3=0x0,*partEvent3=0x0;
 
-  AliHBTEvent * rawtrackEvent, *rawpartEvent;
+  AliHBTEvent * rawtrackEvent, *rawpartEvent;//this we get from Reader
   
-  Int_t nev = fReader->GetNumberOfTrackEvents();
-
   AliHBTPair * trackpair = new AliHBTPair();
   AliHBTPair * partpair = new AliHBTPair();
 
-  TList tbuffer;
-  TList pbuffer;
-  Int_t ninbuffer = 0;
-  UInt_t ii;
+  AliHBTEventBuffer partbuffer(fBufferSize);
+  AliHBTEventBuffer trackbuffer(fBufferSize);
+  
+  register UInt_t ii;
 
   trackEvent1 = new AliHBTEvent();
   partEvent1 = new AliHBTEvent();
   trackEvent1->SetOwner(kFALSE);
   partEvent1->SetOwner(kFALSE);;
   
+  Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
+  
+  fReader->Rewind();
+  
   Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
   Info("ProcessTracksAndParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
   Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
   
-  for (Int_t i = 0;i<nev;i++)
+  for (Int_t i = 0;;i++)//infinite loop
     {
-      rawpartEvent  = fReader->GetParticleEvent(i);
-      rawtrackEvent = fReader->GetTrackEvent(i);
+      if (fReader->Next()) break; //end when no more events available
+      
+      rawpartEvent  = fReader->GetParticleEvent();
+      rawtrackEvent = fReader->GetTrackEvent();
       if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error
-
+      
+      if ( rawpartEvent->GetNumberOfParticles() != rawtrackEvent->GetNumberOfParticles() )
+       {
+         Fatal("ProcessTracksAndParticlesNonIdentAnal",
+               "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
+               i,rawpartEvent->GetNumberOfParticles() , rawtrackEvent->GetNumberOfParticles());
+         return;
+       }
+      
       /********************************/
       /*      Filtering out           */
       /********************************/
-      if ( (fBufferSize != 0) || ( (partEvent2==0x0)||(trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
-       if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
-        {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
-         partEvent2  = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
-         trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
-         ninbuffer--;
-        }
-       else
-        {
+      if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
+       {
          partEvent2  = new AliHBTEvent();
          trackEvent2 = new AliHBTEvent();
-         partEvent2->SetOwner(kFALSE);
-         trackEvent2->SetOwner(kFALSE);
-        }
+         partEvent2->SetOwner(kTRUE);
+         trackEvent2->SetOwner(kTRUE);
+       }
+       
       FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
       
       for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
@@ -965,7 +952,15 @@ void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
 
          part1= partEvent1->GetParticle(j);
          track1= trackEvent1->GetParticle(j);
-
+         
+         if( part1->GetPdgCode() != track1->GetPdgCode() )
+          {
+            Fatal("ProcessTracksAndParticlesNonIdentAnal",
+                  "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
+                  i,j, part1->GetPdgCode(),track1->GetPdgCode() );
+            return;
+          }
+         
          for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
            fParticleMonitorFunctions[ii]->Process(part1);
          for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
@@ -973,11 +968,14 @@ void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
          for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
            fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
 
+         if (nocorrfctns) continue;
+
          /***************************************/
          /******   filling numerators    ********/
          /****** (particles from event2) ********/
          /***************************************/
-         for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
+         
+         for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups 
           {
             part2= partEvent2->GetParticle(k);
             if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
@@ -1007,13 +1005,15 @@ void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
      /***************************************/
      /***** Filling denominators    *********/
      /***************************************/
-     TIter piter(&pbuffer);
-     TIter titer(&tbuffer);
+     partbuffer.ResetIter();
+     trackbuffer.ResetIter();
+     
      Int_t nmonitor = 0;
      
-     while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
+     while ( (partEvent3 = partbuffer.Next() ) != 0x0)
       {
-        trackEvent3 = (AliHBTEvent*)titer.Next();
+        trackEvent3 = trackbuffer.Next();
+        
         if ( (j%fDisplayMixingInfo) == 0) 
           Info("ProcessTracksAndParticlesNonIdentAnal",
                "Mixing particle %d from event %d with particles from event %d",j,i,i-(++nmonitor));
@@ -1046,25 +1046,20 @@ void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
          }//while event2
        }//for over particles in event1
      
-     if ( fBufferSize == 0) continue;//do not mix diff histograms-> do not add to buffer list
-     pbuffer.AddFirst(partEvent2);
-     tbuffer.AddFirst(trackEvent2);
-     partEvent2 = 0x0;
-     trackEvent2 = 0x0;
-     ninbuffer++;
+     partEvent2 = partbuffer.Push(partEvent2);
+     trackEvent2 = trackbuffer.Push(trackEvent2);
 
   }//end of loop over events (1)
- pbuffer.SetOwner();  //to clean stored events by the way of buffer destruction
- tbuffer.SetOwner();  
+  
+ partbuffer.SetOwner(kTRUE);
+ trackbuffer.SetOwner(kTRUE);
  delete partEvent1;
  delete trackEvent1;
+ delete partEvent2;
+ delete trackEvent2;
  delete partpair;
  delete trackpair;
- if ( fBufferSize == 0)
-  {//in that case we did not add these events to list
-    delete partEvent2;
-    delete trackEvent2;
-  }
 
 }
 /*************************************************************************************/  
@@ -1078,43 +1073,39 @@ void AliHBTAnalysis::ProcessTracksNonIdentAnal()
   AliHBTEvent * trackEvent2=0x0;
   AliHBTEvent * trackEvent3=0x0;
 
-
   AliHBTEvent * rawtrackEvent;
   
-  Int_t nev = fReader->GetNumberOfTrackEvents();
-
   AliHBTPair * trackpair = new AliHBTPair();
+  AliHBTEventBuffer trackbuffer(fBufferSize);
 
-  TList tbuffer;
-  Int_t ninbuffer = 0;
   register UInt_t ii;
 
   trackEvent1 = new AliHBTEvent();
   trackEvent1->SetOwner(kFALSE);
   
+  fReader->Rewind();
+  
   Info("ProcessTracksNonIdentAnal","**************************************");
   Info("ProcessTracksNonIdentAnal","*****  NON IDENT MODE ****************");
   Info("ProcessTracksNonIdentAnal","**************************************");
   
-  for (Int_t i = 0;i<nev;i++)
+  
+  for (Int_t i = 0;;i++)//infinite loop
     {
-      rawtrackEvent = fReader->GetTrackEvent(i);
+      if (fReader->Next()) break; //end when no more events available
+      rawtrackEvent = fReader->GetTrackEvent();
+      
       if (rawtrackEvent == 0x0)  continue;//in case of any error
 
       /********************************/
       /*      Filtering out           */
       /********************************/
-      if ( (fBufferSize != 0) || ( trackEvent2==0x0) )//in case fBufferSize == 0 and pointers are created do not eneter
-       if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
-        {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
-         trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
-         ninbuffer--;
-        }
-       else
-        {
+      if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
+       {
          trackEvent2 = new AliHBTEvent();
-         trackEvent2->SetOwner(kFALSE);
-        }
+         trackEvent2->SetOwner(kTRUE);
+       }
+       
       FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
       
       for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
@@ -1127,7 +1118,9 @@ void AliHBTAnalysis::ProcessTracksNonIdentAnal()
 
          for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
            fTrackMonitorFunctions[ii]->Process(track1);
-
+         
+         if (fNTrackFunctions == 0x0) continue;
+         
          /***************************************/
          /******   filling numerators    ********/
          /****** (particles from event2) ********/
@@ -1154,10 +1147,10 @@ void AliHBTAnalysis::ProcessTracksNonIdentAnal()
      /***************************************/
      /***** Filling denominators    *********/
      /***************************************/
-     TIter titer(&tbuffer);
+     trackbuffer.ResetIter();
      Int_t nmonitor = 0;
      
-     while ( (trackEvent3 = (AliHBTEvent*)titer.Next()) != 0x0)
+     while ( (trackEvent3 = trackbuffer.Next() ) != 0x0)
       {
         
         if ( (j%fDisplayMixingInfo) == 0) 
@@ -1184,17 +1177,14 @@ void AliHBTAnalysis::ProcessTracksNonIdentAnal()
          }//while event2
        }//for over particles in event1
      
-     if ( fBufferSize == 0) continue;//do not mix diff histograms     
-     tbuffer.AddFirst(trackEvent2);
-     trackEvent2 = 0x0;
-     ninbuffer++;
+     trackEvent2 = trackbuffer.Push(trackEvent2);
 
   }//end of loop over events (1)
 
+ trackbuffer.SetOwner(kTRUE);
  delete trackpair;
  delete trackEvent1;
- if (fBufferSize == 0) delete trackEvent2;
- tbuffer.SetOwner();  
+ delete trackEvent2;
 }
 /*************************************************************************************/  
 
@@ -1209,39 +1199,36 @@ void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
 
   AliHBTEvent * rawpartEvent = 0x0;
 
-  Int_t nev = fReader->GetNumberOfPartEvents();
-
   AliHBTPair * partpair = new AliHBTPair();
+  AliHBTEventBuffer partbuffer(fBufferSize);
 
-  TList pbuffer;
-  Int_t ninbuffer = 0;
+  register UInt_t ii;
 
   partEvent1 = new AliHBTEvent();
-  partEvent1->SetOwner(kFALSE);;
+  partEvent1->SetOwner(kFALSE);
+  
+  fReader->Rewind();
   
   Info("ProcessParticlesNonIdentAnal","**************************************");
   Info("ProcessParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
   Info("ProcessParticlesNonIdentAnal","**************************************");
 
-  for (Int_t i = 0;i<nev;i++)
+  for (Int_t i = 0;;i++)//infinite loop
     {
-      rawpartEvent  = fReader->GetParticleEvent(i);
+      if (fReader->Next()) break; //end when no more events available
+      
+      rawpartEvent  = fReader->GetParticleEvent();
       if ( rawpartEvent == 0x0  ) continue;//in case of any error
 
       /********************************/
       /*      Filtering out           */
       /********************************/
-      if ( (fBufferSize != 0) || ( partEvent2==0x0) )//in case fBufferSize == 0 and pointers are created do not eneter
-       if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
-        {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
-         partEvent2  = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
-         ninbuffer--; 
-        }
-       else
-        {
+      if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
+       {
          partEvent2  = new AliHBTEvent();
-         partEvent2->SetOwner(kFALSE);
-        }
+         partEvent2->SetOwner(kTRUE);
+       }
+       
       FilterOut(partEvent1, partEvent2, rawpartEvent);
       
       for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
@@ -1252,10 +1239,11 @@ void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
 
          part1= partEvent1->GetParticle(j);
 
-         UInt_t zz;
-         for(zz = 0; zz<fNParticleMonitorFunctions; zz++)
-           fParticleMonitorFunctions[zz]->Process(part1);
-
+         for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+           fParticleMonitorFunctions[ii]->Process(part1);
+         
+         if (fNParticleFunctions == 0) continue;
+         
          /***************************************/
          /******   filling numerators    ********/
          /****** (particles from event2) ********/
@@ -1272,21 +1260,18 @@ void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
               }
             else
              {//meets criteria of the pair cut
-              UInt_t ii;
               for(ii = 0;ii<fNParticleFunctions;ii++)
-                {
                   fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
-                }
              }
            }
      if ( fBufferSize == 0) continue;//do not mix diff histograms
      /***************************************/
      /***** Filling denominators    *********/
      /***************************************/
-     TIter piter(&pbuffer);
+     partbuffer.ResetIter();
      Int_t nmonitor = 0;
 
-     while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
+     while ( (partEvent3 = partbuffer.Next() ) != 0x0)
       {
         if ( (j%fDisplayMixingInfo) == 0) 
             Info("ProcessParticlesNonIdentAnal",
@@ -1303,7 +1288,6 @@ void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
               }
             else
              {//meets criteria of the pair cut
-              UInt_t ii;
               for(ii = 0;ii<fNParticleFunctions;ii++)
                {
                  fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
@@ -1312,16 +1296,13 @@ void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
            }// for particles event2
          }//while event2
        }//for over particles in event1
-     if ( fBufferSize == 0) continue;//do not mix diff histograms
-     pbuffer.AddFirst(partEvent2);
-     partEvent2 = 0x0; 
-     ninbuffer++;
-
+    partEvent2 = partbuffer.Push(partEvent2);
   }//end of loop over events (1)
+  
+ partbuffer.SetOwner(kTRUE);
  delete partpair;
  delete partEvent1;
- if ( fBufferSize == 0) delete partEvent2;
- pbuffer.SetOwner();//to delete stered events.
+ delete partEvent2;
 }
 
 /*************************************************************************************/  
@@ -1369,8 +1350,8 @@ void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, Ali
      
      if (in2)
       {
-        outpart2->AddParticle(part);
-        outtrack2->AddParticle(track);
+        outpart2->AddParticle(new AliHBTParticle(*part));
+        outtrack2->AddParticle(new AliHBTParticle(*track));
         continue;
       }
    }
index de7cdf29fd8cde7e64ba99aedb0b29732576b96d..fd2e0937575dab478bee650cf57092e31950ec4c 100644 (file)
@@ -36,7 +36,42 @@ AliHBTEvent::AliHBTEvent():
    }
  }
 /**************************************************************************/ 
+AliHBTEvent::AliHBTEvent(const AliHBTEvent& source):
+ TObject(source),
+ fSize(source.fSize),
+ fParticles(new AliHBTParticle* [fSize]),
+ fNParticles(source.fNParticles),
+ fOwner(source.fNParticles),
+ fRandomized(source.fRandomized)
+{
+//copy constructor
+  for(Int_t i =0; i<fNParticles; i++)
+   {
+      fParticles[i] = new AliHBTParticle( *(source.fParticles[i]) );
+   }
+  
+}
+/**************************************************************************/ 
 
+AliHBTEvent& AliHBTEvent::operator=(const AliHBTEvent source)
+{
+  // assigment operator
+  Reset();
+  if (fParticles) delete [] fParticles;
+  fSize = source.fSize;
+  fParticles = new AliHBTParticle* [fSize];
+  fNParticles = source.fNParticles;
+  fOwner = source.fNParticles;
+  fRandomized = source.fRandomized;
+      
+  for(Int_t i =0; i<fNParticles; i++)
+   {
+      fParticles[i] = new AliHBTParticle( *(source.fParticles[i]) );
+   }
+  return *this;
+}
+/**************************************************************************/ 
 AliHBTEvent::~AliHBTEvent()
  {
 //destructor   
@@ -54,7 +89,11 @@ void  AliHBTEvent::Reset()
   if(fParticles && fOwner)
     {
       for(Int_t i =0; i<fNParticles; i++)
-        delete fParticles[i];
+       {
+         for (Int_t j = i+1; j<fNParticles; j++)
+           if (fParticles[j] == fParticles[i]) fParticles[j] = 0x0;
+         delete fParticles[i];
+       }
     }
    fNParticles = 0;
 } 
index be7ddb799f3ac5f767fbf3228a5db0b160545473..1523be2eda194131ad5024ee083e12dc4628f869 100644 (file)
@@ -21,7 +21,11 @@ class AliHBTEvent: public TObject
  {
   public:
     AliHBTEvent();
+    AliHBTEvent(const AliHBTEvent& source);
     virtual ~AliHBTEvent();
+    
+    AliHBTEvent& operator=(const AliHBTEvent source);
+        
     static const UInt_t fgkInitEventSize; //initial number of the array
                                           //if expanded, this size is used also
     AliHBTParticle* GetParticle(Int_t n);  //gets particle 
diff --git a/HBTAN/AliHBTEventBuffer.cxx b/HBTAN/AliHBTEventBuffer.cxx
new file mode 100644 (file)
index 0000000..a97478e
--- /dev/null
@@ -0,0 +1,36 @@
+#include "AliHBTEventBuffer.h"
+
+ClassImp(AliHBTEventBuffer)
+
+//______________________________________________________
+////////////////////////////////////////////////////////
+//
+// class AliHBTEventBuffer
+//
+// FIFO type event buffer
+
+AliHBTEventBuffer::AliHBTEventBuffer():
+ fSize(-1),fEvents(),fIter(&fEvents)
+{
+  //ctor
+}
+/***********************************************************/
+AliHBTEventBuffer::AliHBTEventBuffer(Int_t size):
+ fSize(size),fEvents(),fIter(&fEvents)
+{
+  //ctor
+}
+
+AliHBTEvent* AliHBTEventBuffer::Push(AliHBTEvent* event)
+{
+  //adds a new event, and returns old of do not fit in size
+  if (fSize == 0) return event;
+  
+  AliHBTEvent* ret = 0x0;
+  
+  if (fSize == fEvents.GetSize()) 
+    ret = dynamic_cast<AliHBTEvent*>(fEvents.Remove(fEvents.Last()));
+  if (event) fEvents.AddFirst(event);
+  return ret;
+}
+
diff --git a/HBTAN/AliHBTEventBuffer.h b/HBTAN/AliHBTEventBuffer.h
new file mode 100644 (file)
index 0000000..0d25535
--- /dev/null
@@ -0,0 +1,31 @@
+#ifndef ALIHBTEVENTBUFFER_H
+#define ALIHBTEVENTBUFFER_H
+
+#include <TObject.h>
+#include <TList.h>
+#include "AliHBTEvent.h"
+
+class AliHBTEventBuffer: public TObject
+{
+  public:
+    AliHBTEventBuffer();
+    AliHBTEventBuffer(Int_t size);
+    virtual ~AliHBTEventBuffer(){}
+    
+    AliHBTEvent* Push(AliHBTEvent* event);//adds a new event, and returns old of do not fit in size
+    AliHBTEvent* RemoveLast(){return dynamic_cast<AliHBTEvent*>(fEvents.Remove(fEvents.Last()));}
+    void         ResetIter(){fIter.Reset();}
+    AliHBTEvent* Next(){return dynamic_cast<AliHBTEvent*>( fIter.Next() );}
+    void         SetSize(Int_t size){fSize = size;}
+    Int_t        GetSize() const {return fSize;}
+    void         SetOwner(Bool_t flag) {fEvents.SetOwner(flag);}
+  protected:
+  private:
+    Int_t  fSize;//size of buffer; if 0 infinite size
+    TList  fEvents;//list with arrays
+    TIter  fIter;//iterator
+    ClassDef(AliHBTEventBuffer,1)
+};
+
+
+#endif
index fd1403e74846769b2566d537abbb2a785a4939aa..275bc9a9bd2acabf294e6767a1f617a4110baa82 100644 (file)
@@ -687,10 +687,10 @@ void AliHBTFunction3D::BuildHistos(Int_t nxbins, Float_t xmax, Float_t xmin,
                                            //denominator histogram
          
    fNumerator   = new TH3D(numstr.Data(),numstr.Data(),
-                           nxbins,xmin,xmin,nybins,ymin,ymin,nzbins,zmin,zmin);
+                           nxbins,xmin,xmax,nybins,ymin,ymax,nzbins,zmin,zmax);
               
    fDenominator = new TH3D(denstr.Data(),denstr.Data(),
-                           nxbins,xmin,xmin,nybins,ymin,ymin,nzbins,zmin,zmin);
+                           nxbins,xmin,xmax,nybins,ymin,ymax,nzbins,zmin,zmax);
    
    fNumerator->Sumw2();
    fDenominator->Sumw2();
index 319358f4b764fb8bb146a94e9d1df8ce90f880f9..95901b28762ca38838c7189ca5e0538463ac73de 100644 (file)
@@ -26,7 +26,8 @@ AliHBTPairCut::AliHBTPairCut():
 }
 /**********************************************************/
 
-AliHBTPairCut::AliHBTPairCut(const AliHBTPairCut& in)
+AliHBTPairCut::AliHBTPairCut(const AliHBTPairCut& in):
+ TNamed(in)
 {
   //copy constructor
   fCuts = new AliHbtBasePairCut*[fgkMaxCuts];
index 199c2c7ca34c0cdb26e6d4435b6ca04dcc8f4c22..babbee6a1e4cfae948172ac395683ed8510acc40 100644 (file)
@@ -54,6 +54,7 @@ AliHBTParticle::AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx,
 }
 //______________________________________________________________________________
 AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
+   TObject(in),
    fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
    fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
    fCalcMass(in.GetCalcMass()),
@@ -81,6 +82,45 @@ AliHBTParticle::AliHBTParticle(const TParticle &p,Int_t idx):
 }
 //______________________________________________________________________________
 
+AliHBTParticle::~AliHBTParticle()
+{
+//dtor  
+  delete [] fPids;
+  delete [] fPidProb;
+}
+//______________________________________________________________________________
+
+AliHBTParticle& AliHBTParticle::operator=(const AliHBTParticle& in)
+{
+//assigment operator
+  
+  fNPids = in.fNPids;
+  delete [] fPids;
+  delete [] fPidProb;
+  Int_t* fPids = new Int_t[fNPids];
+  Float_t* fPidProb = new Float_t[fNPids];
+  for (Int_t i = 0; i < fNPids;i++)
+   {
+     fPids[i]    = in.fPids[i];
+     fPidProb[i] = in.fPidProb[i];
+   }
+   
+  fPdgIdx = in.fPdgIdx; 
+  fIdxInEvent = in.fIdxInEvent;  
+  fCalcMass = in.GetCalcMass();
+  fPx = in.Px();
+  fPy = in.Py();
+  fPz = in.Pz();
+  fE = in.Energy(); 
+  fVx = in.Vx();
+  fVy = in.Vy();
+  fVz = in.Vz();
+  fVt = in.T();
+  
+  return *this;
+}
+//______________________________________________________________________________
+
 void AliHBTParticle::SetPdgCode(Int_t pdg,Float_t prob)
 {
   SetPIDprobability(pdg,prob);
index 434ca025a80c1eff36cd5e5ed9d81c85727e5a46..3e7d1cc22f699afa7255cf162e9f4abbe7380942 100644 (file)
@@ -28,8 +28,10 @@ public:
 
   AliHBTParticle(const TParticle& p,Int_t idx);
 
-  virtual ~AliHBTParticle(){};
-
+  virtual ~AliHBTParticle();
+  
+  AliHBTParticle& operator=(const AliHBTParticle& in); 
+  
   void           SetPIDprobability(Int_t pdg, Float_t prob = 1.0);
   Float_t        GetPIDprobability(Int_t pdg);
   
index 3c2232910505bfad50d8d7831abaf804305dfd3d..838781f211cba7654d67832c6838f6e8c4698299 100644 (file)
@@ -15,7 +15,8 @@ AliHBTParticleCut::AliHBTParticleCut()
  }
 /******************************************************************/
 
-AliHBTParticleCut::AliHBTParticleCut(const AliHBTParticleCut& in)
+AliHBTParticleCut::AliHBTParticleCut(const AliHBTParticleCut& in):
+ TObject()
 {
   fCuts = new AliHbtBaseCut* [fkgMaxCuts];//last property in the property
                                          //property enum => defines number of properties
@@ -347,7 +348,7 @@ AliHBTLogicalOperCut::~AliHBTLogicalOperCut()
 }
 /******************************************************************/
 
-Bool_t AliHBTLogicalOperCut::AliHBTDummyBaseCut::Pass(AliHBTParticle*p)
+Bool_t AliHBTLogicalOperCut::AliHBTDummyBaseCut::Pass(AliHBTParticle* /*part*/)
 {
   Warning("Pass","You are using dummy base cut! Probobly some logical cut is not set up properly");
   return kFALSE;//accept
index c15612c4dcf1ac2d15104214819c8a1e2e24f2ec..7a33c3ec75f9442912b4483d1015d6eec4900c4b 100644 (file)
@@ -342,15 +342,15 @@ class AliHBTLogicalOperCut:  public AliHbtBaseCut
      AliHBTLogicalOperCut(AliHbtBaseCut* first, AliHbtBaseCut* second);
      virtual   ~AliHBTLogicalOperCut();
    protected:
-     Double_t  GetValue(AliHBTParticle * p){MayNotUse("GetValue");return 0.0;}
+     Double_t  GetValue(AliHBTParticle * /*part*/){MayNotUse("GetValue");return 0.0;}
      
      AliHbtBaseCut* fFirst;   //second cut
      AliHbtBaseCut* fSecond;  //first cut
    private:  
     class  AliHBTDummyBaseCut: public AliHbtBaseCut 
      {
-       Double_t  GetValue(AliHBTParticle * p){return 0.0;}
-       Bool_t    Pass(AliHBTParticle*p);
+       Double_t  GetValue(AliHBTParticle * /*part*/){return 0.0;}
+       Bool_t    Pass(AliHBTParticle* /*part*/);
      };
      
     ClassDef(AliHBTLogicalOperCut,1)
index b3efd2678f366fff9c44418904aeeae6fd259b0c..a76923ff0dbd0b5074919a640851e6543444399a 100644 (file)
@@ -150,7 +150,7 @@ AliHBTRndmGaussBall::AliHBTRndmGaussBall(Float_t rx, Float_t ry, Float_t rz):
 }
 /*********************************************************************/
 
-void AliHBTRndmGaussBall::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*p)
+void AliHBTRndmGaussBall::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*/*particle*/)
 {
 //randomizez gauss for each coordinate separately
   x = gRandom->Gaus(0.0,fRx);
@@ -165,13 +165,13 @@ void AliHBTRndmGaussBall::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTPa
 //                                                                   //
 ///////////////////////////////////////////////////////////////////////
 
-void AliHBTRndmCyllSurf::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*p)
+void AliHBTRndmCyllSurf::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle* particle)
 {
    Double_t r = fR + gRandom->Gaus(0.0, 1.0);
-   Double_t sf = r/p->Pt();//scaling factor for position transformation ->
+   Double_t sf = r/particle->Pt();//scaling factor for position transformation ->
                              //we move direction of string momentum but legth defined by r
-   x = sf*p->Px();
-   y = sf*p->Py();
+   x = sf*particle->Px();
+   y = sf*particle->Py();
    z = gRandom->Uniform(-fL,fL);
   
 }
index accb4df6fa7007b089f5cede0328f068a899130b..16cbb9037e957c6369193148d9dfa8f58b92c1f7 100644 (file)
@@ -64,7 +64,7 @@ class AliHBTRndmGaussBall: public AliHBTRndm
    AliHBTRndmGaussBall(Float_t r);
    AliHBTRndmGaussBall(Float_t rx, Float_t ry, Float_t rz);
    ~AliHBTRndmGaussBall(){}
-   void Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*p);
+   void Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*/*particle*/);
   private:
    Float_t fRx;
    Float_t fRy;
@@ -79,7 +79,7 @@ class AliHBTRndmCyllSurf: public AliHBTRndm
    AliHBTRndmCyllSurf(Float_t r, Float_t l):fR(r),fL(l){}
    ~AliHBTRndmCyllSurf(){}
    
-   void Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*p);
+   void Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle* particle);
   private:
    Float_t fR;
    Float_t fL;
index dc0435c29751775a7efadbb8d05bccb94fe8b188..b98e1f9813c0ca19a8131be3d377ead592537a1b 100644 (file)
@@ -4,29 +4,49 @@
 #include <TObjString.h>
 #include <TObjArray.h>
 #include <TClass.h>
-#include <Riostream.h>
 
 #include "AliHBTParticleCut.h"
-
+#include "AliHBTEvent.h"
+#include "AliHBTRun.h"
  
 ClassImp(AliHBTReader)
 //pure virtual
-
+    
 /*************************************************************************************/
 
-AliHBTReader::AliHBTReader()
+AliHBTReader::AliHBTReader():
+ fCuts(new TObjArray()),
+ fDirs(0x0),
+ fCurrentEvent(0),
+ fCurrentDir(0),
+ fNEventsRead(0),
+ fTracksEvent(0x0),
+ fParticlesEvent(0x0),
+ fParticles(0x0),
+ fTracks(0x0),
+ fIsRead(kFALSE),
+ fBufferEvents(kFALSE)
 {
 //constructor
- fCuts = new TObjArray();
- fDirs = 0x0;
 }
+/*************************************************************************************/
 
+AliHBTReader::AliHBTReader(TObjArray* dirs):
+ fCuts(new TObjArray()),
+ fDirs(dirs),
+ fCurrentEvent(0),
+ fCurrentDir(0),
+ fNEventsRead(0),
+ fTracksEvent(0x0),
+ fParticlesEvent(0x0),
+ fParticles(0x0),
+ fTracks(0x0),
+ fIsRead(kFALSE),
+ fBufferEvents(kFALSE)
+{
+//ctor with array of directories to read as parameter
+}
 /*************************************************************************************/
-AliHBTReader::AliHBTReader(TObjArray* dirs)
- {
-  fCuts = new TObjArray();
-  fDirs = dirs;
- }
 
 AliHBTReader::~AliHBTReader()
 {
@@ -36,8 +56,26 @@ AliHBTReader::~AliHBTReader()
    fCuts->SetOwner();
    delete fCuts;
   }
+ delete fParticlesEvent;
+ delete fTracksEvent;
 }
+/*************************************************************************************/
 
+Int_t AliHBTReader::Next()
+{
+//moves to next event
+  if ( ReadNext() == kTRUE)
+     return kTRUE;
+  
+  if (fBufferEvents)
+   {
+     if ( ReadsTracks() && fTracksEvent) 
+       fTracks->SetEvent(fNEventsRead-1,fTracksEvent);
+     if ( ReadsParticles() && fParticlesEvent)
+       fParticles->SetEvent(fNEventsRead-1,fParticlesEvent);
+   }
+  return kFALSE;
+}
 /*************************************************************************************/
 
 void AliHBTReader::AddParticleCut(AliHBTParticleCut* cut)
@@ -52,7 +90,137 @@ void AliHBTReader::AddParticleCut(AliHBTParticleCut* cut)
   AliHBTParticleCut *c = (AliHBTParticleCut*)cut->Clone();
   fCuts->Add(c);
 }
+/********************************************************************/
+
+AliHBTEvent* AliHBTReader::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+  if (ReadsParticles() == kFALSE)
+   {
+     Error("GetParticleEvent","This reader is not able to provide simulated particles.");
+     return 0;
+   } 
+   
+  if (!fIsRead) 
+   { 
+    if (ReadsParticles() && (fParticles == 0x0)) fParticles = new AliHBTRun();
+    if (ReadsTracks() && (fTracks == 0x0)) fTracks = new AliHBTRun();
+    
+    if (Read(fParticles,fTracks))
+     {
+       Error("GetParticleEvent","Error in reading");
+       return 0x0;
+     }
+    else fIsRead = kTRUE;
+   }
+  return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+
+AliHBTEvent* AliHBTReader::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+  if (ReadsTracks() == kFALSE)
+   {
+     Error("GetTrackEvent","This reader is not able to provide recosntructed tracks.");
+     return 0;
+   } 
+  if (!fIsRead) 
+   {
+    if (ReadsParticles() && (fParticles == 0x0)) fParticles = new AliHBTRun();
+    if (ReadsTracks() && (fTracks == 0x0)) fTracks = new AliHBTRun();
+    
+    if(Read(fParticles,fTracks))
+     {
+       Error("GetTrackEvent","Error in reading");
+       return 0x0;
+     }
+    else fIsRead = kTRUE;
+   }
+  return fTracks->GetEvent(n);
+ }
+/********************************************************************/
 
+Int_t AliHBTReader::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+  if (ReadsParticles() == kFALSE)
+   {
+     Error("GetNumberOfPartEvents","This reader is not able to provide simulated particles.");
+     return 0;
+   } 
+   
+  if (!fIsRead) 
+   {
+    if (ReadsParticles() && (fParticles == 0x0)) fParticles = new AliHBTRun();
+    if (ReadsTracks() && (fTracks == 0x0)) fTracks = new AliHBTRun();
+    
+    if (Read(fParticles,fTracks))
+     {
+       Error("GetNumberOfPartEvents","Error in reading");
+       return 0;
+     }
+    else fIsRead = kTRUE;
+   }
+   return fParticles->GetNumberOfEvents();
+ }
+/********************************************************************/
+Int_t AliHBTReader::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+  if (ReadsTracks() == kFALSE)
+   {
+     Error("GetNumberOfTrackEvents","This reader is not able to provide recosntructed tracks.");
+     return 0;
+   } 
+  if (!fIsRead)
+   {
+     if (ReadsParticles() && (fParticles == 0x0)) fParticles = new AliHBTRun();
+     if (ReadsTracks() && (fTracks == 0x0)) fTracks = new AliHBTRun();
+     
+     if(Read(fParticles,fTracks))
+      {
+        Error("GetNumberOfTrackEvents","Error in reading");
+        return 0;
+      }
+     else fIsRead = kTRUE;
+   }
+  return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+Int_t AliHBTReader::Read(AliHBTRun* particles, AliHBTRun *tracks)
+{
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+  Info("Read","");
+  
+  if (!particles) //check if an object is instatiated
+   {
+     Error("Read"," particles object must be instatiated before passing it to the reader");
+     return 1;
+   }
+  if (!tracks)  //check if an object is instatiated
+   {
+     Error("Read"," tracks object must be instatiated before passing it to the reader");
+     return 1;
+   }
+  particles->Reset();//clear runs == delete all old events
+  tracks->Reset();
+  
+  Rewind();
+  
+  Int_t i = 0;
+  while(Next() == kFALSE)
+   {
+     if (ReadsTracks()) tracks->SetEvent(i,fTracksEvent);
+     if (ReadsParticles()) particles->SetEvent(i,fParticlesEvent);
+     i++;
+   }
+  return 0;
+}      
 /*************************************************************************************/
 
 Bool_t AliHBTReader::Pass(AliHBTParticle* p)
@@ -98,39 +266,40 @@ Bool_t  AliHBTReader::Pass(Int_t pid)
 /*************************************************************************************/
 
 TString& AliHBTReader::GetDirName(Int_t entry)
- {
-   TString* retval;//return value
-   if (fDirs ==  0x0)
-    {
-      retval = new TString(".");
-      return *retval;
-    }
-   
-   if ( (entry>fDirs->GetEntries()) || (entry<0))//if out of bounds return empty string
-    {                                            //note that entry==0 is accepted even if array is empty (size=0)
-      Error("GetDirName","Name out of bounds");
-      retval = new TString();
-      return *retval;
-    }
-   
-   if (fDirs->GetEntries() == 0)
-    { 
-      retval = new TString(".");
-      return *retval;
-    }
-   
-   TClass *objclass = fDirs->At(entry)->IsA();
-   TClass *stringclass = TObjString::Class();
-   
-   TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry));
-   
-   if(dir == 0x0)
-    {
-      Error("GetDirName","Object in TObjArray is not a TObjString or its descendant");
-      retval = new TString();
-      return *retval;
-    }
-   if (gDebug > 0) cout<<"Returned ok "<<dir->String().Data()<<endl;
-   return dir->String();
- }
+{
+//returns directory name of next one to read
+  TString* retval;//return value
+  if (fDirs ==  0x0)
+   {
+     retval = new TString(".");
+     return *retval;
+   }
+
+  if ( (entry>fDirs->GetEntries()) || (entry<0))//if out of bounds return empty string
+   {                                            //note that entry==0 is accepted even if array is empty (size=0)
+     Error("GetDirName","Name out of bounds");
+     retval = new TString();
+     return *retval;
+   }
+
+  if (fDirs->GetEntries() == 0)
+   { 
+     retval = new TString(".");
+     return *retval;
+   }
+
+  TClass *objclass = fDirs->At(entry)->IsA();
+  TClass *stringclass = TObjString::Class();
+
+  TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry));
+
+  if(dir == 0x0)
+   {
+     Error("GetDirName","Object in TObjArray is not a TObjString or its descendant");
+     retval = new TString();
+     return *retval;
+   }
+  if (gDebug > 0) Info("GetDirName","Returned ok %s",dir->String().Data());
+  return dir->String();
+}
 
index f2d5d7d39b2fb8c73d114f73e3d9f2a30084e7ec..a121fe2ba3357a142b9d749dd0885226a7eb7ef0 100644 (file)
@@ -2,6 +2,7 @@
 #define ALIHBTREADER_H
 
 #include <TNamed.h>
+#include <TObjArray.h>
 
 //Reader Base class (reads particles and tracks and
 //puts it to the AliHBTRun objects
@@ -10,7 +11,6 @@
 class AliHBTRun;
 class AliHBTEvent;
 class AliHBTParticleCut;
-class TObjArray;
 class AliHBTParticle;
 class TString;
  
@@ -20,26 +20,50 @@ class AliHBTReader: public TNamed
     AliHBTReader();
     AliHBTReader(TObjArray*);
     virtual ~AliHBTReader();
-    //in the future this class is will read global tracking
-    virtual Int_t Read(AliHBTRun* particles, AliHBTRun *tracks) = 0;
     
-    virtual AliHBTEvent* GetParticleEvent(Int_t) = 0;
-    virtual AliHBTEvent* GetTrackEvent(Int_t) = 0;
-    virtual Int_t GetNumberOfPartEvents() = 0;
-    virtual Int_t GetNumberOfTrackEvents() = 0;
+    virtual Int_t        Next();
+    virtual void         Rewind() = 0;
     
-    void AddParticleCut(AliHBTParticleCut* cut);
+    virtual Bool_t       ReadsTracks() const = 0;
+    virtual Bool_t       ReadsParticles() const = 0;
     
-    void SetDirs(TObjArray* dirs){fDirs = dirs;} //sets array directories names
-
-    virtual AliHBTEvent* GetNextParticleEvent(){return 0x0;}
-    virtual AliHBTEvent* GetNextTrackEvent(){return 0x0;}
-
+    void                 AddParticleCut(AliHBTParticleCut* cut);
+    
+    virtual Int_t        Read(AliHBTRun* particles, AliHBTRun *tracks);
+    
+    virtual AliHBTEvent* GetParticleEvent() const {return fParticlesEvent;}
+    virtual AliHBTEvent* GetTrackEvent() const {return fTracksEvent;}
+    
+    virtual AliHBTEvent* GetParticleEvent(Int_t);
+    virtual AliHBTEvent* GetTrackEvent(Int_t);
+    
+    virtual Int_t        GetNumberOfPartEvents();
+    virtual Int_t        GetNumberOfTrackEvents();
+    
+    void                 SetDirs(TObjArray* dirs){fDirs = dirs;} //sets array directories names
+    void                 SetEventBuffering(Bool_t flag){fBufferEvents = flag;}
+    
+    virtual Int_t GetNumberOfDirs() const {return (fDirs)?fDirs->GetEntries():0;}
   protected:
     
-    TObjArray *fCuts;//array with particle cuts
-    TObjArray *fDirs;//arry with directories to read data from
-
+    TObjArray*    fCuts;//array with particle cuts
+    TObjArray*    fDirs;//arry with directories to read data from
+    
+    Int_t         fCurrentEvent;//!  number of current event in current directory
+    Int_t         fCurrentDir;//! number of current directory
+    
+    Int_t         fNEventsRead;//!total 
+        
+    AliHBTEvent*  fTracksEvent;    //! tracks read from current event
+    AliHBTEvent*  fParticlesEvent; //! particles read from current event
+    
+    AliHBTRun*    fParticles; //!simulated particles
+    AliHBTRun*    fTracks; //!reconstructed tracks (particles)
+    
+    Bool_t        fIsRead;//!flag indicating if the data are already read
+    Bool_t        fBufferEvents;//flag indicating if the data should be bufferred
+    
+    virtual Int_t ReadNext() = 0; //this methods reads next event and put result in fTracksEvent and/or fParticlesEvent
     Bool_t Pass(AliHBTParticle*);
     Bool_t Pass(Int_t pid);
 
index 0dbffd2c6b1a2df3c5a797cd173ec179a4c3a999..e0fdfe330bd1805e9c4786eecad9b14090f7fd2f 100644 (file)
@@ -7,6 +7,8 @@
 #include <TFile.h>
 #include <TParticle.h>
 
+#include <AliRunLoader.h>
+#include <AliStack.h>
 #include <AliESDtrack.h>
 #include <AliESD.h>
 
 #include "AliHBTParticle.h"
 #include "AliHBTParticleCut.h"
 
+AliHBTReaderESD r;
+
 ClassImp(AliHBTReaderESD)
 
-AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename):
- fParticles(new AliHBTRun()),
- fTracks(new AliHBTRun()),
+AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
  fESDFileName(esdfilename),
- fIsRead(kFALSE)
+ fGAlFileName(galfilename),
+ fFile(0x0),
+ fRunLoader(0x0),
+ fReadParticles(kFALSE)
 {
   //cosntructor
   if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
@@ -29,12 +34,13 @@ AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename):
 }
 /********************************************************************/
   
-AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename):
+AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
  AliHBTReader(dirs), 
- fParticles(new AliHBTRun()),
- fTracks(new AliHBTRun()),
  fESDFileName(esdfilename),
- fIsRead(kFALSE)
+ fGAlFileName(galfilename),
+ fFile(0x0),
+ fRunLoader(0x0),
+ fReadParticles(kFALSE)
 {
   //cosntructor
   if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
@@ -45,239 +51,215 @@ AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename):
 AliHBTReaderESD::~AliHBTReaderESD()
 {
  //desctructor
-  delete fParticles;
-  delete fTracks;
+  delete fFile;
+  delete fRunLoader;
 }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderESD::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
-   if (!fIsRead) 
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetParticleEvent","Error in reading");
-       return 0x0;
-     }
-   return fParticles->GetEvent(n);
- }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderESD::GetTrackEvent(Int_t n)
- {
- //returns Nth event with reconstructed tracks
-   if (!fIsRead) 
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetTrackEvent","Error in reading");
-       return 0x0;
-     }
-   return fTracks->GetEvent(n);
- }
-/********************************************************************/
-
-Int_t AliHBTReaderESD::GetNumberOfPartEvents()
- {
- //returns number of events of particles
-   if (!fIsRead) 
-    if ( Read(fParticles,fTracks))
-     {
-       Error("GetNumberOfPartEvents","Error in reading");
-       return 0;
-     }
-   return fParticles->GetNumberOfEvents();
- }
-
-/********************************************************************/
-Int_t AliHBTReaderESD::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
-  if (!fIsRead)
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetNumberOfTrackEvents","Error in reading");
-       return 0;
-     }
-  return fTracks->GetNumberOfEvents();
- }
-/********************************************************************/
-
-
-Int_t AliHBTReaderESD::Read(AliHBTRun* particles, AliHBTRun *tracks)
+/**********************************************************/
+Int_t AliHBTReaderESD::ReadNext()
 {
- //reads data and puts put to the particles and tracks objects
- //reurns 0 if everything is OK
- //
-  Info("Read","");
-  Int_t totalNevents = 0;
-  Int_t currentdir = 0; //number of current directory name is fDirs array
+//reads next event from fFile
+  
+  //****** Tentative particle type "concentrations"
+  static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
   
   Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
   Double_t w[kNSpecies];
   Double_t mom[3];//momentum
   Double_t pos[3];//position
-   //****** Tentative particle type "concentrations"
-  const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
   
+  if (AliHBTParticle::fgDebug)
+    Info("ReadNext","Entered");
+    
   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
   if (pdgdb == 0x0)
    {
-     Error("Read","Can not get PDG Database Instance.");
+     Error("Next","Can not get PDG Database Instance.");
      return 1;
    }
-  if (!particles) //check if an object is instatiated
-   {
-     Error("Read"," particles object must instatiated before passing it to the reader");
-     return 1;
-   }
-  if (!tracks)  //check if an object is instatiated
-   {
-     Error("Read"," tracks object must instatiated before passing it to the reader");
-     return 1;
-   }
-  particles->Reset();//clear runs == delete all old events
-  tracks->Reset();
-
-  Int_t Ndirs;
-  if (fDirs) //if array with directories is supplied by user
-   {
-     Ndirs = fDirs->GetEntries(); //get the number if directories
-   }
-  else
-   {
-     Ndirs = 0; //if the array is not supplied read only from current directory
-   }
-
+  
+  if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
+  if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
+  
+  fParticlesEvent->Reset();
+  fTracksEvent->Reset();
+        
   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
    {
-     TFile* file = OpenFile(currentdir);
-     if (file == 0x0)
-       {
-         Error("Read","Cannot get File for dir no. %d",currentdir);
-         currentdir++;
-         continue;
-       } 
-       
-     for(Int_t currentEvent = 0; ;currentEvent++)//loop over all events
+     if (fFile == 0x0)
       {
-       TString esdname;
-       esdname+=currentEvent;
-       AliESD* esd = dynamic_cast<AliESD*>(file->Get(esdname));
-       if (esd == 0x0)
+       fFile = OpenFile(fCurrentDir);//rl is opened here
+       if (fFile == 0x0)
          {
-           if (AliHBTParticle::fgDebug > 2 )
-            {
-              Info("Read","Can not find ESD object named %s.",esdname.Data());
-            }
-           currentdir++;
-           break;//we have to assume there is no more ESD objects in the file
-         } 
-       
-       Info("Read","Reading Event %d",currentEvent);
+           Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
+           fCurrentDir++;
+           continue;
+         }
+       fCurrentEvent = 0;
+      } 
+      
+     //try to read
+     TString esdname;
+     esdname+=fCurrentEvent;
+     AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
+     if (esd == 0x0)
+      {
+        if (AliHBTParticle::fgDebug > 2 )
+          {
+            Info("ReadNext","Can not find ESD object named %s.",esdname.Data());
+          }
+        fCurrentDir++;
+        delete fFile;//we have to assume there is no more ESD objects in the fFile
+        delete fRunLoader;
+        fFile = 0x0;
+        continue;
+      }
+     AliStack* stack = 0x0;
+     if (fReadParticles && fRunLoader)
+      {
+        fRunLoader->GetEvent(fCurrentEvent);
+        stack = fRunLoader->Stack();
+      }
+     Info("ReadNext","Reading Event %d",fCurrentEvent);
   
-       Int_t ntr = esd->GetNumberOfTracks();
-       Info("Read","Found %d tracks.",ntr);
-       for (Int_t i = 0;i<ntr; i++)
-        {
-          AliESDtrack *esdtrack = esd->GetTrack(i);
-          if (esdtrack == 0x0)
-           {
-             Error("Read","Can not get track %d", i);
-             continue;
-           }
-          if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) 
-           {
-             if (AliHBTParticle::fgDebug > 2) 
-               Info("Read","Particle skipped: PID BIT is not set.");
-             continue;
-           }
-
-          esdtrack->GetESDpid(pidtable);
-          esdtrack->GetPxPyPz(mom);
-          esdtrack->GetXYZ(pos);
-          Double_t rc=0.;
+     Int_t ntr = esd->GetNumberOfTracks();
+     Info("ReadNext","Found %d tracks.",ntr);
+     for (Int_t i = 0;i<ntr; i++)
+      {
+        AliESDtrack *esdtrack = esd->GetTrack(i);
+        if (esdtrack == 0x0)
+         {
+           Error("Next","Can not get track %d", i);
+           continue;
+         }
+        if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) 
+         {
+           if (AliHBTParticle::fgDebug > 2) 
+             Info("ReadNext","Particle skipped: PID BIT is not set.");
+           continue;
+         }
 
-           //Here we apply Bayes' formula
-          for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
-          if (rc==0.0) 
-           {
-             if (AliHBTParticle::fgDebug > 2) 
-               Info("Read","Particle rejected since total bayessian PID probab. is zero.");
-             continue;
-           }
+        esdtrack->GetESDpid(pidtable);
+        esdtrack->GetPxPyPz(mom);
+        esdtrack->GetXYZ(pos);
+        
+        //Particle from kinematics
+        AliHBTParticle* particle = 0;
+        Bool_t keeppart = kFALSE;
+        if ( fReadParticles && stack  )
+         {
+           if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track 
+           TParticle *p = stack->Particle(esdtrack->GetLabel());
+           if (p==0x0) 
+            {
+              Error("ReadNext","Can not find track with such label.");
+              continue;
+            }
+           particle = new AliHBTParticle(*p,i);
+         }
+         
+        //Here we apply Bayes' formula
+        Double_t rc=0.;
+        for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
+        if (rc==0.0) 
+         {
+           if (AliHBTParticle::fgDebug > 2) 
+             Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
+           continue;
+         }
 
-          for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
+        for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
 
-          if (AliHBTParticle::fgDebug > 4)
-           { 
-             Info("Read","###########################################################################");
-             Info("Read","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
-             Info("Read","Position: %f %f %f",pos[0],pos[1],pos[2]);
-             Info("Read","Radius: %f ,rc: %f",TMath::Hypot(pos[0],pos[1]),rc);
-             TString msg("Pid list got from track:");
-             for (Int_t s = 0;s<kNSpecies;s++)
-              {
-                msg+="\n    ";
-                msg+=s;
-                msg+="(";
-                msg+=GetSpeciesPdgCode((ESpecies)s);
-                msg+="): ";
-                msg+=w[s];
-                msg+=" (";
-                msg+=pidtable[s];
-                msg+=")";
-              }
-             Info("Read","%s",msg.Data());
-           }//if (AliHBTParticle::fgDebug>4)
-           
-          for (Int_t s = 0; s<kNSpecies; s++)
-           {
-             if (w[s] == 0.0) continue;
+        if (AliHBTParticle::fgDebug > 4)
+         { 
+           Info("ReadNext","###########################################################################");
+           Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
+           Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
+           TString msg("Pid list got from track:");
+           for (Int_t s = 0;s<kNSpecies;s++)
+            {
+              msg+="\n    ";
+              msg+=s;
+              msg+="(";
+              msg+=GetSpeciesPdgCode((ESpecies)s);
+              msg+="): ";
+              msg+=w[s];
+              msg+=" (";
+              msg+=pidtable[s];
+              msg+=")";
+            }
+           Info("ReadNext","%s",msg.Data());
+         }//if (AliHBTParticle::fgDebug>4)
 
-             Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
-             if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type 
+        for (Int_t s = 0; s<kNSpecies; s++)
+         {
+           if (w[s] == 0.0) continue;
 
-             Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
-             Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
+           Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
+           if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type 
 
-             AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, 
-                                                        mom[0], mom[1], mom[2], tEtot,
-                                                        pos[0], pos[1], pos[2], 0.);
-             //copy probabilitis of other species (if not zero)
-             for (Int_t k = 0; k<kNSpecies; k++)
-              {
-                if (k == s) continue;
-                if (w[k] == 0.0) continue;
-                track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
-              }
-             if(Pass(track))//check if meets all criteria of any of our cuts
-                            //if it does not delete it and take next good track
-               { 
-                 delete track;
-                 continue;
-               }
-             tracks->AddParticle(totalNevents,track);
-             if (AliHBTParticle::fgDebug > 4 )
-              {
-                Info("Read","\n\nAdding Particle with incarnation %d",pdgcode);
-                track->Print();
-              }
-           }//for (Int_t s = 0; s<kNSpecies; s++)
+           Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
+           Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
 
-        }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
-       totalNevents++;
-      }//for(Int_t currentEvent = 0; ;currentEvent++) -- for over ESD in single file
+           AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, 
+                                                      mom[0], mom[1], mom[2], tEtot,
+                                                      pos[0], pos[1], pos[2], 0.);
+           //copy probabilitis of other species (if not zero)
+           for (Int_t k = 0; k<kNSpecies; k++)
+            {
+              if (k == s) continue;
+              if (w[k] == 0.0) continue;
+              track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
+            }
+           if(Pass(track))//check if meets all criteria of any of our cuts
+                          //if it does not delete it and take next good track
+            { 
+              delete track;
+              continue;
+            }
+           
+           fTracksEvent->AddParticle(track);
+           if (particle) fParticlesEvent->AddParticle(particle);
+           keeppart = kTRUE;
+               
+           if (AliHBTParticle::fgDebug > 4 )
+            {
+              Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
+              track->Print();
+            }
+         }//for (Int_t s = 0; s<kNSpecies; s++)
+         
+        if (keeppart == kFALSE) delete particle;//particle was not stored in event
+        
+      }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
+     
+     Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
+            fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
+            fNEventsRead,fCurrentEvent,fCurrentDir);
       
-   }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
+     fCurrentEvent++;
+     fNEventsRead++;
+     delete esd;
+     return 0;//success -> read one event
+   }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
+   
+  return 1; //no more directories to read
+}
+/**********************************************************/
+void AliHBTReaderESD::Rewind()
+{
+  delete fFile;
+  fFile = 0;
+  delete fRunLoader;
+  fCurrentDir = 0;
+  fNEventsRead = 0;
+  fCurrentEvent++;
   
-  fIsRead = kTRUE;
-  return 0;
-}      
+}
 /**********************************************************/
-   
+
 TFile* AliHBTReaderESD::OpenFile(Int_t n)
 {
-//opens file with kine tree
+//opens fFile with kine tree
 
  const TString& dirname = GetDirName(n);
  if (dirname == "")
@@ -290,15 +272,33 @@ TFile* AliHBTReaderESD::OpenFile(Int_t n)
 
  if ( ret == 0x0)
   {
-    Error("OpenFiles","Can't open file %s",filename.Data());
+    Error("OpenFiles","Can't open fFile %s",filename.Data());
     return 0x0;
   }
  if (!ret->IsOpen())
   {
-    Error("OpenFiles","Can't open file  %s",filename.Data());
+    Error("OpenFiles","Can't open fFile  %s",filename.Data());
     return 0x0;
   }
  
+ if (fReadParticles)
+  {
+   fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
+   if (fRunLoader == 0x0)
+    {
+      Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
+      delete ret;
+      return 0x0;
+    }
+   fRunLoader->LoadHeader();
+   if (fRunLoader->LoadKinematics())
+    {
+      Error("Next","Error occured while loading kinematics.");
+      delete fRunLoader;
+      delete ret;
+      return 0x0;
+    }
+  }
  return ret;
 }
 /**********************************************************/
index ca1c86aa1c6cdbd9ca9bc641f5c0b47166b0a946..b5c888f8734ca127644ec152454eb8b7c93f75ad 100644 (file)
 
 #include <TString.h>
 class TFile;
+class AliRunLoader;
 
 class AliHBTReaderESD: public AliHBTReader
 {
   public:
-    AliHBTReaderESD(const Char_t* esdfilename = "AliESDs.root");
+    AliHBTReaderESD(const Char_t* esdfilename = "AliESDs.root", const Char_t* galfilename = "galice.root");
 
-    AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename = "AliESDs.root");
+    AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename = "AliESDs.root", const Char_t* galfilename = "galice.root");
 
     virtual ~AliHBTReaderESD();
     
-    Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+    void          Rewind();
     
-    AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
-    AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles 
-    Int_t        GetNumberOfPartEvents();//returns number of particle events
-    Int_t        GetNumberOfTrackEvents();//returns number of track events
+    void          ReadParticles(Bool_t flag){fReadParticles = flag;}
+    Bool_t        ReadsTracks() const {return kTRUE;}
+    Bool_t        ReadsParticles() const {return fReadParticles;}
     
     enum ESpecies {kESDElectron = 0, kESDMuon, kESDPion, kESDKaon, kESDProton, kNSpecies};
-    static Int_t GetSpeciesPdgCode(ESpecies spec);//skowron
-  protected:
-      
-    TFile*       OpenFile(Int_t evno);//opens files to be read for given event
-    void         CloseFiles(TFile*);//close files
+    static Int_t  GetSpeciesPdgCode(ESpecies spec);//skowron
     
-    AliHBTRun*   fParticles; //!simulated particles
-    AliHBTRun*   fTracks; //!reconstructed tracks (particles)
-
-    TString      fESDFileName;//name of the file with tracks
-    Bool_t       fIsRead;//!flag indicating if the data are already read
+  protected:
+    Int_t         ReadNext();
+    TFile*        OpenFile(Int_t evno);//opens files to be read for given event
+    void          CloseFiles(TFile*);//close files
+
+    TString       fESDFileName;//name of the file with tracks
+    TString       fGAlFileName;//name of the file with tracks
+    TFile*        fFile;//! pointer to current ESD file
+    AliRunLoader* fRunLoader;//!Run Loader
+    Bool_t        fReadParticles;//flag indicating wether to read particles from kinematics
     
   private:
-    ClassDef(AliHBTReaderESD,1)
+    ClassDef(AliHBTReaderESD,2)
 };
 
 
index dd56bda19c255b13845a40411766e33c512d8099..2018438327760f51ceed98b7c25de9ef7494e487 100644 (file)
@@ -1,22 +1,18 @@
-
 #include "AliHBTReaderITSv2.h"
 
-#include <Riostream.h>
-#include <Riostream.h>
+#include <Varargs.h>
+
 #include <TString.h>
 #include <TObjString.h>
 #include <TTree.h>
-#include <TFile.h>
 #include <TParticle.h>
 
 #include <AliRun.h>
 #include <AliRunLoader.h>
 #include <AliLoader.h>
+#include <AliStack.h>
 #include <AliMagF.h>
-#include <AliITS.h>
 #include <AliITStrackV2.h>
-#include <AliITStrackerV2.h>
-#include <AliITSgeom.h>
 
 #include "AliHBTRun.h"
 #include "AliHBTEvent.h"
 
 ClassImp(AliHBTReaderITSv2)
 
-AliHBTReaderITSv2::AliHBTReaderITSv2():fFileName("galice.root")
+AliHBTReaderITSv2::AliHBTReaderITSv2():
+ fFileName("galice.root"),
+ fRunLoader(0x0),
+ fITSLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
 {
   //constructor, 
   //Defaults:
   //  galicefilename = "galice.root"
-  fParticles = 0x0;
-  fTracks    = 0x0;
-  fIsRead = kFALSE;
 }
 /********************************************************************/
 
-AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):fFileName(galicefilename)
+AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fITSLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
 {
   //constructor, 
   //Defaults:
   //  galicefilename = "galice.root"
-  fParticles = new AliHBTRun();
-  fTracks    = new AliHBTRun();
-  fIsRead = kFALSE;
 }
 /********************************************************************/
 
 AliHBTReaderITSv2::AliHBTReaderITSv2(TObjArray* dirs, const Char_t* galicefilename): 
-       AliHBTReader(dirs), fFileName(galicefilename)
+ AliHBTReader(dirs),
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fITSLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
 {
   //constructor, 
   //Defaults:
   //  galicefilename = "galice.root"
   
-  fParticles = new AliHBTRun();
-  fTracks    = new AliHBTRun();
-  fIsRead = kFALSE;
 }
 /********************************************************************/
-
-AliHBTReaderITSv2::~AliHBTReaderITSv2()
- {
-   if (fParticles) delete fParticles;
-   if (fTracks) delete fTracks;
- }
-/********************************************************************/
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
- if (!fIsRead) 
-  { 
-    if (fParticles == 0x0) fParticles = new AliHBTRun();
-    if (fTracks == 0x0) fTracks    = new AliHBTRun();
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetParticleEvent","Error in reading");
-       return 0x0;
-     }
-  }
- return (fParticles)?fParticles->GetEvent(n):0x0;
+void AliHBTReaderITSv2::Rewind()
+{
+  //rewinds reading
+  delete fRunLoader;
+  fRunLoader = 0x0;
+  fCurrentDir = 0;
+  fNEventsRead= 0;
 }
 /********************************************************************/
 
-AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n)
- {
- //returns Nth event with reconstructed tracks
- if (!fIsRead) 
-  { 
-    if (fParticles == 0x0) fParticles = new AliHBTRun();
-    if (fTracks == 0x0) fTracks    = new AliHBTRun();
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetTrackEvent","Error in reading");
-       return 0x0;
-     }
-   }
-  return (fTracks)?fTracks->GetEvent(n):0x0;
- }
-/********************************************************************/
-
-Int_t AliHBTReaderITSv2::GetNumberOfPartEvents()
- {
- //returns number of events of particles
-  if (!fIsRead)
-   {
-    if (fParticles == 0x0) fParticles = new AliHBTRun();
-    if (fTracks == 0x0) fTracks    = new AliHBTRun();
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetNumberOfPartEvents","Error in reading");
-       return 0;
-     }
-   }
-  return (fParticles)?fParticles->GetNumberOfEvents():0;
- }
-
-/********************************************************************/
-Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
-  if (!fIsRead) 
-   {
-    if (fParticles == 0x0) fParticles = new AliHBTRun();
-    if (fTracks == 0x0) fTracks    = new AliHBTRun();
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetNumberOfTrackEvents","Error in reading");
-       return 0;
-     }
-   }
-  return (fTracks)?fTracks->GetNumberOfEvents():0;
- }
-/********************************************************************/
+AliHBTReaderITSv2::~AliHBTReaderITSv2()
+{
+  //dtor
+  delete fRunLoader;
+}
 /********************************************************************/
-
  
-Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
+Int_t AliHBTReaderITSv2::ReadNext()
 {
-//reads data
- Int_t Nevents = 0; //number of events found in given directory
- Int_t Ndirs; //number of the directories to be read
- Int_t Ntracks; //number of tracks in current event
- Int_t currentdir = 0; //number of events in the current directory 
- Int_t totalNevents = 0; //total number of events read from all directories up to now
+//reads data from next event
+  
  register Int_t i = 0; //iterator
  
 // AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
- TTree *tracktree; // tree for tracks
+ TTree *tracktree = 0x0; // tree for tracks
+ AliITStrackV2 *iotrack = 0x0;
  
  Double_t xk;
  Double_t par[5]; //Kalman track parameters
  Float_t phi, lam, pt;//angles and transverse momentum
  Int_t label; //label of the current track
-
- AliITStrackV2 *iotrack = 0x0; //buffer track for reading data from tree
-
- if (!particles) //check if an object is instatiated
-  {
-    Error("Read"," particles object must instatiated before passing it to the reader");
-  }
- if (!tracks)  //check if an object is instatiated
-  {
-    Error("Read"," tracks object must instatiated before passing it to the reader");
-  }
- particles->Reset();//clear runs == delete all old events
- tracks->Reset();
-
- if (fDirs) //if array with directories is supplied by user
-  {
-    Ndirs = fDirs->GetEntries(); //get the number if directories
-  }
- else
-  {
-    Ndirs = 0; //if the array is not supplied read only from current directory
-  }
  
-// cout<<"Found "<<Ndirs<<" directory entries"<<endl;
  
+ if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
+ if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
+
+ fParticlesEvent->Reset();
+ fTracksEvent->Reset();
  do //do while is good even if Ndirs==0 (than read from current directory)
    {
-    TString filename = GetDirName(currentdir);
-    if (filename.IsNull())
+    if (fRunLoader == 0x0) 
+      if (OpenNextFile()) continue;//directory counter is increased inside in case of error
+
+    if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
      {
-       Error("Read","Can not get directory name");
-       currentdir++;
-       continue;
+       //read next directory
+       delete fRunLoader;//close current session
+       fRunLoader = 0x0;//assure pointer is null
+       fCurrentDir++;//go to next dir
+       continue;//directory counter is increased inside in case of error
      }
-    filename = filename +"/"+ fFileName;
-    AliRunLoader* rl = AliRunLoader::Open(filename);
-    if( rl == 0x0)
+     
+    Info("ReadNext","Reading Event %d",fCurrentEvent);
+     
+    fRunLoader->GetEvent(fCurrentEvent);
+    
+    tracktree=fITSLoader->TreeT();
+    if (!tracktree) 
      {
-       Error("Read","Exiting due to problems with opening files.");
-       currentdir++;
+       Error("ReadNext","Can't get a tree with ITS tracks"); 
+       fCurrentEvent++;
        continue;
      }
-    
-    rl->LoadHeader();
-    rl->LoadKinematics();
-    rl->LoadgAlice();
-    gAlice = rl->GetAliRun();
-    AliITS* its = (AliITS*)gAlice->GetModule("ITS");
-    
-    AliLoader* itsl = rl->GetLoader("ITSLoader");
-    
-    if ((its == 0x0) || ( itsl== 0x0))
+      
+    TBranch *tbranch=tracktree->GetBranch("tracks");
+    if (!tbranch) 
      {
-       Error("Read","Can not found ITS in this run");
-       delete rl;
-       rl = 0x0;
-       currentdir++;
+       Error("ReadNext","Can't get a branch with ITS tracks"); 
+       fCurrentEvent++;
        continue;
      }
-    Nevents = rl->GetNumberOfEvents();
-    if (Nevents > 0)//check if tree E exists
+
+    AliStack* stack = fRunLoader->Stack();
+    if (stack == 0x0)
      {
-      Info("Read","________________________________________________________");
-      Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
-      Float_t mf;
-      if (fUseMagFFromRun)
-       {
-         mf = gAlice->Field()->SolenoidField();
-         Info("Read","Setting Magnetic Field from run: B=%fT",mf/10.);
-       }
-      else
-       {
-         Info("Read","Setting Own Magnetic Field: B=%fT",fMagneticField);
-         mf = fMagneticField*10.;
-       }
-      AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
-      if (iotrack == 0x0) iotrack = new AliITStrackV2();
-     }
-    else
-     {//if not return an error
-       Error("Read","No events in this run");
-       delete rl;
-       rl = 0x0;
-       currentdir++;
+       Error("ReadNext","Can not get stack for current event",fCurrentEvent);
+       fCurrentEvent++;
        continue;
      }
+     
+    //must be here because on the beginning conv. const. is not set yet 
+    if (iotrack == 0x0) iotrack = new AliITStrackV2(); //buffer track for reading data from tree
     
-    AliITSgeom *geom= its->GetITSgeom();
-    if (!geom) 
+    Int_t ntr = (Int_t)tracktree->GetEntries();
+    
+    for (i=0; i < ntr; i++) //loop over all tpc tracks
      { 
-       Error("Read","Can't get the ITS geometry!"); 
-       delete rl;
-       rl = 0x0;
-       currentdir++;
-       continue;
-     }
+       tbranch->SetAddress(&iotrack);
+       tracktree->GetEvent(i);
 
-    itsl->LoadTracks();
+       label=iotrack->GetLabel();
+       if (label < 0) 
+        {
+          continue;
+        }
 
-    for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
-     {
-       cout<<"Reading Event "<<currentEvent<<endl;
-       rl->GetEvent(currentEvent);
-       tracktree=itsl->TreeT();
-       
-       if (!tracktree) 
-         {
-           Error("Read","Can't get a tree with ITS tracks"); 
-           continue;
-         }
-       TBranch *tbranch=tracktree->GetBranch("tracks");
-       Ntracks=(Int_t)tracktree->GetEntries();
+       TParticle *p = stack->Particle(label);
+       if(p == 0x0) continue; //if returned pointer is NULL
+       if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
 
-       Int_t accepted = 0;
-       Int_t tpcfault = 0;
-       Int_t itsfault = 0;
-       for (i=0; i<Ntracks; i++) //loop over all tpc tracks
-        { 
-          if(i%100 == 0)cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"\r";
-          
-          tbranch->SetAddress(&iotrack);
-          tracktree->GetEvent(i);
+       if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
+                                           //if not take next partilce
 
-          label=iotrack->GetLabel();
-          if (label < 0) 
-           {
-             tpcfault++;
-             continue;
-           }
+       AliHBTParticle* part = new AliHBTParticle(*p,i);
+       if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
+                                               //if it does not delete it and take next good track
 
-          TParticle *p = (TParticle*)gAlice->Particle(label);
-          if(p == 0x0) continue; //if returned pointer is NULL
-          if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
+       iotrack->PropagateTo(3.,0.0028,65.19);
+       iotrack->PropagateToVertex();
 
-          if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
-                                              //if not take next partilce
-            
-          AliHBTParticle* part = new AliHBTParticle(*p,i);
-          if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
-                                                  //if it does not delete it and take next good track
+       iotrack->GetExternalParameters(xk,par);     //get properties of the track
+       phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
+       if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+       if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+       lam=par[3]; 
+       pt=1.0/TMath::Abs(par[4]);
 
-          iotrack->PropagateTo(3.,0.0028,65.19);
-          iotrack->PropagateToVertex();
-          iotrack->GetExternalParameters(xk,par);     //get properties of the track
-          phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
-          if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
-          if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
-          lam=par[3]; 
-          pt=1.0/TMath::Abs(par[4]);
-            
-          Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
-          Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
-          Double_t tpz = pt * lam; //track z coordinate of momentum
-           
-          Double_t mass = p->GetMass();
-          Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
-            
-          AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
-          if(Pass(track))//check if meets all criteria of any of our cuts
-                         //if it does not delete it and take next good track
-           { 
-            delete track;
-            delete part;
-            continue;
-           }
-          particles->AddParticle(totalNevents,part);//put track and particle on the run
-          tracks->AddParticle(totalNevents,track);
-          accepted++;
-        }//end of loop over tracks in the event
+       Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
+       Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
+       Double_t tpz = pt * lam; //track z coordinate of momentum
+
+       Double_t mass = p->GetMass();
+       Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
+
+       AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+       if(Pass(track))//check if meets all criteria of any of our cuts
+                      //if it does not delete it and take next good track
+        { 
+         delete track;
+         delete part;
+         continue;
+        }
+        
+       fParticlesEvent->AddParticle(part);
+       fTracksEvent->AddParticle(track);
+     }//end of loop over tracks in the event
        
-       totalNevents++;
-       cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"   its faults: "<<itsfault<<endl;
+    Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
+            fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
+            fNEventsRead,fCurrentEvent,fCurrentDir);
      
-     }//end of loop over events in current directory
-    delete rl;
-    currentdir++;
-   }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
+    fCurrentEvent++;
+    fNEventsRead++;
+    delete iotrack;
+    return 0;
+   }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
 
  delete iotrack;
- fIsRead = kTRUE;
- return 0;
+ return 1;
 }
 
 /********************************************************************/
+Int_t AliHBTReaderITSv2::OpenNextFile()
+{
+  //opens next file
+  TString filename = GetDirName(fCurrentDir);
+  if (filename.IsNull())
+   {
+     DoOpenError("Can not get directory name");
+     return 1;
+   }
+  filename = filename +"/"+ fFileName;
+  fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
+  if( fRunLoader == 0x0)
+   {
+     DoOpenError("Can not open session.");
+     return 1;
+   }
+
+  if (fRunLoader->GetNumberOfEvents() <= 0)
+   {
+     DoOpenError("There is no events in this directory.");
+     return 1;
+   }
+
+  if (fRunLoader->LoadKinematics())
+   {
+     DoOpenError("Error occured while loading kinematics.");
+     return 1;
+   }
+  fITSLoader = fRunLoader->GetLoader("ITSLoader");
+  if ( fITSLoader == 0x0)
+   {
+     DoOpenError("Exiting due to problems with opening files.");
+     return 1;
+   }
+   
+  Info("OpenNextSession","________________________________________________________");
+  Info("OpenNextSession","Found %d event(s) in directory %s",
+        fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
+  Float_t mf;
+  if (fUseMagFFromRun)
+   {
+     if (fRunLoader->LoadgAlice())
+      {
+        DoOpenError("Error occured while loading AliRun.");
+        return 1;
+      }
+     mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
+     Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
+     fRunLoader->UnloadgAlice();
+   }
+  else
+   {
+     Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
+     if (fMagneticField == 0x0)
+      {
+        Fatal("OpenNextSession","Magnetic field can not be 0.");
+        return 1;//pro forma
+      }
+     mf = fMagneticField*10.;
+   }
+  AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
+
+  if (fITSLoader->LoadTracks())
+   {
+     DoOpenError("Error occured while loading TPC tracks.");
+     return 1;
+   }
+  
+  fCurrentEvent = 0;
+  return 0;
+}
 /********************************************************************/
 
+void AliHBTReaderITSv2::DoOpenError( const char *va_(fmt), ...)
+{
+  // Does error display and clean-up in case error caught on Open Next Session
 
+   va_list ap;
+   va_start(ap,va_(fmt));
+   Error("OpenNextFile", va_(fmt), ap);
+   va_end(ap);
+   
+   delete fRunLoader;
+   fRunLoader = 0x0;
+   fITSLoader = 0x0;
+   fCurrentDir++;
+}
+
+/********************************************************************/
index 75d6d5ea9f6e19211ab2c05516a6aa18ceef8fdd..e52a0e59ff08d3e13bdf43be21a1002919fafbbd 100644 (file)
@@ -3,9 +3,10 @@
 
 #include "AliHBTReader.h"
 
-#include <TString.h>
 
-class TFile;
+class AliLoader;
+class AliRunLoader;
+class TString;
 
 class AliHBTReaderITSv2: public AliHBTReader
 {
@@ -17,29 +18,29 @@ class AliHBTReaderITSv2: public AliHBTReader
 
     virtual ~AliHBTReaderITSv2();
     
-    Int_t Read(AliHBTRun*, AliHBTRun*);//reads tracks and particles and puts them in runs
+    void          Rewind();
     
-    AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
-    AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles 
-    Int_t GetNumberOfPartEvents();//returns number of particle events
-    Int_t GetNumberOfTrackEvents();//returns number of track events
-
-    void SetMagneticField(Float_t mf){fMagneticField=mf;}
-    void UseMagneticFieldFromRun(Bool_t flag = kTRUE){fUseMagFFromRun=flag;}
+    Bool_t        ReadsTracks() const {return kTRUE;}
+    Bool_t        ReadsParticles() const {return kTRUE;}
     
-  protected:
+    void          SetMagneticField(Float_t mf){fMagneticField=mf;}
+    void          UseMagneticFieldFromRun(Bool_t flag = kTRUE){fUseMagFFromRun=flag;}
     
-    AliHBTRun* fParticles; //!simulated particles
-    AliHBTRun* fTracks; //!reconstructed tracks (particles)
+  protected:
     
-    TString fFileName;//name of the file with galice.root
-  
-    Bool_t fIsRead;//!flag indicating if the data are already read
-
-    Float_t    fMagneticField;//magnetic field value that was enforced while reading
-    Bool_t     fUseMagFFromRun;//flag indicating if using field specified in gAlice (kTRUE)
+    Int_t         ReadNext();//reads tracks and particles and puts them in runs
+    Int_t         OpenNextFile();
+    void          DoOpenError( const char *va_(fmt), ...);
+        
+    TString       fFileName;//name of the file with galice.root
+    AliRunLoader* fRunLoader;//!Run Loader
+    AliLoader*    fITSLoader;//! ITS Loader
+        
+    Float_t       fMagneticField;//magnetic field value that was enforced while reading
+    Bool_t        fUseMagFFromRun;//flag indicating if using field specified in gAlice (kTRUE)
                                // or enforece other defined by fMagneticField
-    ClassDef(AliHBTReaderITSv2,1)
+    
+    ClassDef(AliHBTReaderITSv2,2)
 };
 
 #endif
index c874fa514b83432cb3444890f8fb20faecd90f93..109869a682bf1334a14672730b56cb9f99dc1c27 100644 (file)
 #include "AliHBTParticle.h"
 #include "AliHBTParticleCut.h"
 
+AliHBTReaderInternal test;
 
 ClassImp(AliHBTReaderInternal)
 /********************************************************************/
 
-AliHBTReaderInternal::AliHBTReaderInternal()
+AliHBTReaderInternal::AliHBTReaderInternal():
+ fFileName(),
+ fPartBranch(0x0),
+ fTrackBranch(0x0),
+ fTree(0x0),
+ fFile(0x0)
 {
 //Defalut constructor
-  fParticles = 0x0; 
-  fTracks = 0x0;
-  fIsRead = kFALSE;
 }
 /********************************************************************/
 
-AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):fFileName(filename)
+AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):
+ fFileName(filename),
+ fPartBranch(0x0),
+ fTrackBranch(0x0),
+ fTree(0x0),
+ fFile(0x0)
 { 
 //constructor 
 //filename - name of file to open
- fParticles = new AliHBTRun();
- fTracks    = new AliHBTRun();
- fIsRead = kFALSE;
 }
 /********************************************************************/
 
 AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
-  AliHBTReader(dirs),fFileName(filename)
+ AliHBTReader(dirs),
+ fFileName(filename),
+ fPartBranch(0x0),
+ fTrackBranch(0x0),
+ fTree(0x0),
+ fFile(0x0)
 { 
 //ctor
 //dirs contains strings with directories to look data in
 //filename - name of file to open
- fParticles = new AliHBTRun();
- fTracks    = new AliHBTRun();
- fIsRead = kFALSE;
 }
 /********************************************************************/
 
 AliHBTReaderInternal::~AliHBTReaderInternal()
  {
  //desctructor
-   delete fParticles;
-   delete fTracks;
+   delete fTree;
+   delete fFile;
  }
 /********************************************************************/
-
-AliHBTEvent* AliHBTReaderInternal::GetParticleEvent(Int_t n)
- {
- //returns nth event with simulated particles
-   if (!fIsRead) 
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetParticleEvent","Error in reading");
-       return 0x0;
-     }
-   return fParticles->GetEvent(n);
- }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderInternal::GetTrackEvent(Int_t n)
- {
- //returns nth event with reconstructed tracks
-   if (!fIsRead) 
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetTrackEvent","Error in reading");
-       return 0x0;
-     }
-   return fTracks->GetEvent(n);
- }
-/********************************************************************/
-
-Int_t AliHBTReaderInternal::GetNumberOfPartEvents()
- {
- //returns number of events of particles
-   if (!fIsRead) 
-    if ( Read(fParticles,fTracks))
-     {
-       Error("GetNumberOfPartEvents","Error in reading");
-       return 0;
-     }
-   return fParticles->GetNumberOfEvents();
- }
-
-/********************************************************************/
-Int_t AliHBTReaderInternal::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
-  if (!fIsRead)
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetNumberOfTrackEvents","Error in reading");
-       return 0;
-     }
-  return fTracks->GetNumberOfEvents();
- }
+void AliHBTReaderInternal::Rewind()
+{
+  delete fTree;
+  fTree = 0x0;
+  delete fFile;
+  fFile = 0x0;
+  fCurrentDir = 0;
+  fNEventsRead= 0;
+}
 /********************************************************************/
 
-Int_t AliHBTReaderInternal::Read(AliHBTRun* particles, AliHBTRun *tracks)
+Int_t AliHBTReaderInternal::ReadNext()
  {
  //reads data and puts put to the particles and tracks objects
  //reurns 0 if everything is OK
  //
-  Info("Read","");
+  Info("ReadNext","");
   Int_t i; //iterator and some temprary values
-  Int_t totalnevents = 0; //total number of read events
-  Int_t nevents = 0;
-  Int_t currentdir = 0;
-  Int_t Ndirs;
   Int_t counter;
-  TFile *aFile;//file with tracks
   AliHBTParticle* tpart = 0x0, *ttrack = 0x0;
+  
   TDatabasePDG* pdgdb = TDatabasePDG::Instance();  
   if (pdgdb == 0x0)
    {
-     Error("Read","Can not get PDG Particles Data Base");
-     return 1;
-   }
-  if (!particles) //check if an object is instatiated
-   {
-     Error("Read"," particles object must instatiated before passing it to the reader");
-     return 1;
-   }
-  if (!tracks)  //check if an object is instatiated
-   {
-     Error("Read"," tracks object must instatiated before passing it to the reader");
+     Error("ReadNext","Can not get PDG Particles Data Base");
      return 1;
    }
-  particles->Reset();//clear runs == delete all old events
-  tracks->Reset();
-  if (fDirs) //if array with directories is supplied by user
-   {
-     Ndirs = fDirs->GetEntries(); //get the number if directories
-   }
-  else
-   {
-     Ndirs = 0; //if the array is not supplied read only from current directory
-   }
-
-  TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
-  TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
-//  pbuffer->BypassStreamer(kFALSE);
-//  tbuffer->BypassStreamer(kFALSE);
+   
+  if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
+  if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
+  
+  fParticlesEvent->Reset();
+  fTracksEvent->Reset();
   
   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
    {
-    
-    if( (i=OpenFile(aFile, currentdir)) )
+    if (fTree == 0x0)
+     if( (i=OpenNextFile()) )
+      {
+        Error("ReadNext","Skippimg directory due to problems with opening files. Errorcode %d",i);
+        fCurrentDir++;
+        continue;
+      }
+    if (fCurrentEvent == (Int_t)fTree->GetEntries())
      {
-       Error("Read","Skippimg directory due to problems with opening files. Errorcode %d",i);
-       currentdir++;
+       delete fTree;
+       fTree = 0x0;
+       delete fFile;
+       fFile = 0x0;
+       fPartBranch = 0x0;
+       fTrackBranch= 0x0;
+       fCurrentDir++;
        continue;
      }
    /***************************/
    /***************************/
    /***************************/
     
-     TTree* tree = (TTree*)aFile->Get("data");
-     if (tree == 0x0)
-      {
-       Error("Read","Can not get the tree");
-       currentdir++;
-       continue;
-      }
-     
-     TBranch *trackbranch=tree->GetBranch("tracks");//get the branch with tracks
-     if (trackbranch == 0x0) ////check if we got the branch
-       {//if not return with error
-         Warning("Read","Can't find a branch with tracks !\n"); 
-       }
-     else
-      {
-        trackbranch->SetAddress(&tbuffer);
-      }
+        
+    Info("ReadNext","Event %d",fCurrentEvent);
+    fTree->GetEvent(fCurrentEvent);
 
-     TBranch *partbranch=tree->GetBranch("particles");//get the branch with particles
-     if (partbranch == 0x0) ////check if we got the branch
-       {//if not return with error
-         Warning("Read","Can't find a branch with particles !\n"); 
-       }
-     else
-      {
-        partbranch->SetAddress(&pbuffer);
-      }
-     
-     nevents = (Int_t)tree->GetEntries();
-     Info("Read","________________________________________________________");
-     Info("Read","Found %d event(s) in directory %s",nevents,GetDirName(currentdir).Data());
-     
-     for (Int_t currentEvent =0; currentEvent<nevents;currentEvent++)
-      {
-       Info("Read","Event %d",currentEvent);
-       tree->GetEvent(currentEvent);
-       
-       counter = 0;  
-       if (partbranch && trackbranch)
-        {
-           for(i = 0; i < pbuffer->GetEntries(); i++)
+    counter = 0;  
+    if (fPartBranch && fTrackBranch)
+    {
+       for(i = 0; i < fPartBuffer->GetEntries(); i++)
+         {
+           tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
+           ttrack =  dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
+
+           if( tpart == 0x0 ) continue; //if returned pointer is NULL
+
+           if (ttrack->GetUID() != tpart->GetUID())
              {
-               tpart = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
-               ttrack =  dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
-
-               if( tpart == 0x0 ) continue; //if returned pointer is NULL
-               
-               if (ttrack->GetUID() != tpart->GetUID())
-                 {
-                   Error("Read","Sth. is wrong: Track and Particle has different UID.");
-                   Error("Read","They probobly do not correspond to each other.");
-                 }
-               
-               for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+               Error("ReadNext","Sth. is wrong: Track and Particle has different UID.");
+               Error("ReadNext","They probobly do not correspond to each other.");
+             }
+
+           for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+            {
+              if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
+              if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
+                                          //if not take next partilce
+              AliHBTParticle* part = new AliHBTParticle(*tpart);
+              part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
+              if( Pass(part) )
                 {
-                  if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
-                  if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
-                                              //if not take next partilce
-                  AliHBTParticle* part = new AliHBTParticle(*tpart);
-                  part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
-                  if( Pass(part) )
-                    {
-         delete part;
-         continue; 
-       }
-                  AliHBTParticle* track = new AliHBTParticle(*ttrack);
-                  
-                  particles->AddParticle(totalnevents,part);//put track and particle on the run
-                  tracks->AddParticle(totalnevents,track);
-                  counter++;
+                  delete part;
+                  continue; 
                 }
-             }
-            Info("Read","   Read: %d particles and tracks.",counter);
-        }
-       else
-        {  
-         if (partbranch)
-          {
-            Info("Read","Found %d particles in total.",pbuffer->GetEntries()); 
-            for(i = 0; i < pbuffer->GetEntries(); i++)
-             { 
-               tpart = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
-               if(tpart == 0x0) continue; //if returned pointer is NULL
-               
-               for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+              AliHBTParticle* track = new AliHBTParticle(*ttrack);
+
+              fParticlesEvent->AddParticle(part);
+              fTracksEvent->AddParticle(track);
+
+              counter++;
+            }
+         }
+        Info("ReadNext","   Read: %d particles and tracks.",counter);
+    }
+   else
+    {  
+     if (fPartBranch)
+      {
+        Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());     
+        for(i = 0; i < fPartBuffer->GetEntries(); i++)
+         { 
+           tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
+           if(tpart == 0x0) continue; //if returned pointer is NULL
+
+           for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+            {
+              if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
+              if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
+
+              AliHBTParticle* part = new AliHBTParticle(*tpart);
+              part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
+              if( Pass(part) )
                 {
-                  if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
-                  if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
-               
-                  AliHBTParticle* part = new AliHBTParticle(*tpart);
-                  part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
-                  if( Pass(part) )
-                    {
-         delete part;
-         continue; 
-       }
-                  particles->AddParticle(totalnevents,part);//put track and particle on the run
-                  counter++;
+                  delete part;
+                  continue; 
                 }
-             }
-            Info("Read","   Read: %d particles.",counter);
-          }
-         else Info("Read","   Read: 0 particles.");
-         
-         if (trackbranch)
-          {
-            for(i = 0; i < tbuffer->GetEntries(); i++)
-             {
-               tpart = dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
-               if(tpart == 0x0) continue; //if returned pointer is NULL
-               for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+              fParticlesEvent->AddParticle(part);
+              counter++;
+            }
+         }
+        Info("ReadNext","   Read: %d particles.",counter);
+      }
+     else Info("ReadNext","   Read: 0 particles.");
+
+     if (fTrackBranch)
+      {
+        for(i = 0; i < fTrackBuffer->GetEntries(); i++)
+         {
+           tpart = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
+           if(tpart == 0x0) continue; //if returned pointer is NULL
+           for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+            {
+              if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
+              if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
+
+              AliHBTParticle* part = new AliHBTParticle(*tpart);
+              part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
+              if( Pass(part) )
                 {
-                  if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
-                  if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
-               
-                  AliHBTParticle* part = new AliHBTParticle(*tpart);
-                  part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
-                  if( Pass(part) )
-                    {
-         delete part;
-         continue; 
-       }
-                  tracks->AddParticle(totalnevents,part);//put track and particle on the run
-                  counter++;
+                  delete part;
+                  continue; 
                 }
-             }
-            Info("Read","   Read: %d tracks",counter);
-          }
-         else Info("Read","   Read: 0 tracks.");
-        }
-        totalnevents++;
-       }
-
-    /***************************/
-    /***************************/
-    /***************************/
-    currentdir++;
-    delete tree;
-    aFile->Close();
-    delete aFile;
-    aFile = 0x0;
-
-   }while(currentdir < Ndirs);
+              fTracksEvent->AddParticle(part);
+              counter++;
+            }
+         }
+        Info("ReadNext","   Read: %d tracks",counter);
+      }
+     else Info("ReadNext","   Read: 0 tracks.");
+    }
+    fPartBuffer->Delete();
+    fTrackBuffer->Delete();
+    
+    fCurrentEvent++;
+    fNEventsRead++;
+    
+    return 0;
 
-  delete pbuffer;
-  delete tbuffer;
-  fIsRead = kTRUE;
-  return 0;
+   }while(fCurrentEvent < GetNumberOfDirs());
+
+  return 1;//no more directories to read
 }
 /********************************************************************/
 
-Int_t AliHBTReaderInternal::OpenFile(TFile*& aFile,Int_t event)
+Int_t AliHBTReaderInternal::OpenNextFile()
 {
-
-   const TString& dirname = GetDirName(event); 
+  //open file in current directory
+   const TString& dirname = GetDirName(fCurrentEvent); 
    if (dirname == "")
     {
-      Error("OpenFile","Can not get directory name");
+      Error("OpenNextFile","Can not get directory name");
       return 4;
     }
    
    TString filename = dirname +"/"+ fFileName;
-   aFile = TFile::Open(filename.Data());
-   if ( aFile  == 0x0 ) 
+   fFile = TFile::Open(filename.Data());
+   if ( fFile  == 0x0 ) 
      {
-       Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
+       Error("OpenNextFile","Can't open file named %s",filename.Data());
        return 1;
      }
-   if (!aFile->IsOpen())
+   if (fFile->IsOpen() == kFALSE)
      {
-       Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
+       Error("OpenNextFile","Can't open filenamed %s",filename.Data());
        return 1;
      }
+   
+   fTree = (TTree*)fFile->Get("data");
+   if (fTree == 0x0)
+    {
+     Error("OpenNextFile","Can not get the tree.");
+     return 1;
+    }
+
+    
+   fTrackBranch = fTree->GetBranch("tracks");//get the branch with tracks
+   if (fTrackBranch == 0x0) ////check if we got the branch
+     {//if not return with error
+       Info("OpenNextFile","Can't find a branch with tracks !\n"); 
+     }
+   else
+    {
+      if (fTrackBuffer == 0x0)
+        fTrackBuffer = new TClonesArray("AliHBTParticle",15000);
+      fTrackBranch->SetAddress(&fTrackBuffer);
+    }
+
+   fPartBranch = fTree->GetBranch("particles");//get the branch with particles
+   if (fPartBranch == 0x0) ////check if we got the branch
+     {//if not return with error
+       Info("OpenNextFile","Can't find a branch with particles !\n"); 
+     }
+   else
+    {
+      if (fPartBuffer == 0x0)
+        fPartBuffer = new TClonesArray("AliHBTParticle",15000);
+      fPartBranch->SetAddress(&fPartBuffer);
+    }
+
+   Info("OpenNextFile","________________________________________________________");
+   Info("OpenNextFile","Found %d event(s) in directory %s",
+         (Int_t)fTree->GetEntries(),GetDirName(fCurrentEvent).Data());
+   
+   fCurrentEvent = 0;
+
    return 0; 
 }
 /********************************************************************/
@@ -370,8 +323,6 @@ Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile, Bool
 
   TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
   TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
-//  pbuffer->BypassStreamer(kFALSE);
-//  tbuffer->BypassStreamer(kFALSE);
 
   TClonesArray &particles = *pbuffer;
   TClonesArray &tracks = *tbuffer;
index 4419f2335ac03d5e437fa73ccbd073b25e5bade7..7f4cf6b9d15e05f4cefd91bead28011ac8910d13 100644 (file)
@@ -13,6 +13,8 @@
 
 #include <TString.h>
 class TFile;
+class TTree;
+class TBranch;
 class TClonesArray;
 
 class AliHBTReaderInternal: public AliHBTReader
@@ -23,23 +25,27 @@ class AliHBTReaderInternal: public AliHBTReader
     AliHBTReaderInternal(TObjArray* dirs, const char *filename);
     virtual ~AliHBTReaderInternal();
     
-    Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
-    static Int_t Write(AliHBTReader* reader,const char* outfile = "data.root", Bool_t multcheck = kFALSE);//reads tracks from runs and writes them to file
+    void          Rewind();
+    
+    Bool_t        ReadsTracks() const {return kTRUE;}
+    Bool_t        ReadsParticles() const {return kTRUE;}
     
-    AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
-    AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles 
-    Int_t GetNumberOfPartEvents();//returns number of particle events
-    Int_t GetNumberOfTrackEvents();//returns number of track events
+    static Int_t Write(AliHBTReader* reader,const char* outfile = "data.root", Bool_t multcheck = kFALSE);//reads tracks from runs and writes them to file
     
   protected:
-    AliHBTRun* fParticles; //!simulated particles
-    AliHBTRun* fTracks; //!reconstructed tracks (particles)
-    Bool_t     fIsRead;//!flag indicating if the data are already read    
-    TString    fFileName;//name of the file with tracks
-
-    Int_t      OpenFile(TFile*& aFile,Int_t event);//opens file to be read for given event 
+    
+    TString       fFileName;//name of the file with tracks
+    TBranch*      fPartBranch;//!branch with particles
+    TBranch*      fTrackBranch;//!branch with tracks
+    TTree*        fTree;//!tree
+    TFile*        fFile;//!file
+    TClonesArray* fPartBuffer;//!buffer array that tree is read to
+    TClonesArray* fTrackBuffer;//!
+        
+    Int_t         ReadNext();//reads next event
+    Int_t         OpenNextFile();//opens file to be read for given event 
     static Bool_t FindIndex(TClonesArray* arr,Int_t idx);
     
-    ClassDef(AliHBTReaderInternal,1)
+    ClassDef(AliHBTReaderInternal,2)
 };
 #endif 
index 8160be6b1622d5c3ad95bbe67e00a1daa8add2f1..d5b32de096b0d81835b818702cd7b123a3beb256 100644 (file)
@@ -1,6 +1,5 @@
 #include "AliHBTReaderKineTree.h"
 
-#include <Riostream.h>
 #include <TString.h>
 #include <TObjString.h>
 #include <TTree.h>
@@ -21,114 +20,86 @@ ClassImp(AliHBTReaderKineTree)
 /**********************************************************/
 const TString AliHBTReaderKineTree::fgkEventFolderName("HBTReaderKineTree");
 
-AliHBTReaderKineTree::AliHBTReaderKineTree():fFileName("galice.root")
+AliHBTReaderKineTree::AliHBTReaderKineTree():
+ fFileName("galice.root"),
+ fRunLoader(0x0)
 {
-  fParticles = new AliHBTRun();
-  fIsRead = kFALSE;
+  //ctor
 }
+/**********************************************************/
 
-AliHBTReaderKineTree::
-AliHBTReaderKineTree(TString& fname):fFileName(fname)
+AliHBTReaderKineTree::AliHBTReaderKineTree(TString& fname):
+ fFileName(fname),
+ fRunLoader(0x0)
 {
-  fParticles = new AliHBTRun();
-  fIsRead = kFALSE;
+  //ctor
 }
-
-
 /**********************************************************/
-AliHBTReaderKineTree::
-AliHBTReaderKineTree(TObjArray* dirs,const Char_t *filename):AliHBTReader(dirs),fFileName(filename)
+
+AliHBTReaderKineTree::AliHBTReaderKineTree(TObjArray* dirs,const Char_t *filename):
+ AliHBTReader(dirs),
+ fFileName(filename),
+ fRunLoader(0x0)
 {
-  fParticles = new AliHBTRun();
-  fIsRead = kFALSE;
+  //ctor
 }
 /**********************************************************/
 
-AliHBTEvent* AliHBTReaderKineTree::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
-   if (!fIsRead) 
-    if(Read(fParticles,0x0))
-     {
-       Error("GetParticleEvent","Error in reading");
-       return 0x0;
-     }
-
-   return fParticles->GetEvent(n);
- }
-
-/********************************************************************/
-
-Int_t AliHBTReaderKineTree::GetNumberOfPartEvents()
- {
- //returns number of events of particles
-   if (!fIsRead)
-    if(Read(fParticles,0x0))
-     {
-       Error("GetNumberOfPartEvents","Error in reading");
-       return 0;
-     }
-   return fParticles->GetNumberOfEvents();
- }
-
+AliHBTReaderKineTree::~AliHBTReaderKineTree()
+{
+  //dtor
+  delete fRunLoader;
+}
+/**********************************************************/
 
+void AliHBTReaderKineTree::Rewind()
+{
+  delete fRunLoader;
+  fRunLoader = 0x0;
+  fCurrentDir = 0;
+  fNEventsRead= 0;  
+}
 /**********************************************************/
-Int_t AliHBTReaderKineTree::Read(AliHBTRun* particles, AliHBTRun* /*tracks*/)
+
+Int_t AliHBTReaderKineTree::ReadNext()
 {
  //Reads Kinematics Tree
   
  Info("Read","");
- if (!particles) //check if an object is instatiated
-   {
-     Error("Read"," particles object must instatiated before passing it to the reader");
-     return 1;
-   }  
- particles->Reset();//clear runs == delete all old events
- Int_t Ndirs;//total number of directories specified in fDirs array
- Int_t Nevents; //number of events read in current directory
- Int_t totalNevents = 0; //total number of read events 
- Int_t currentdir = 0; //number of current directory name is fDirs array
-
- if (fDirs) //if array with directories is supplied by user
-  {
-    Ndirs = fDirs->GetEntries(); //get the number if directories
-  }
- else
-  {
-    Ndirs = 0; //if the array is not supplied read only from current directory
-  }
+ if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
+ fParticlesEvent->Reset();
 
  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
   { 
-    Info("Read","________________________________________________________");
-    AliRunLoader * rl = OpenFile(currentdir);
-
-    if (rl == 0x0)
+    if (fRunLoader == 0x0) 
+      if (OpenNextFile()) 
+        { 
+          fCurrentDir++;
+          continue;
+        }
+    
+    if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
      {
-      Error("Read","Cannot open session");
-      return 2;
+       //read next directory
+       delete fRunLoader;//close current session
+       fRunLoader = 0x0;//assure pointer is null
+       fCurrentDir++;//go to next dir
+       continue;//directory counter is increased inside in case of error
      }
-    rl->LoadHeader();
-    rl->LoadKinematics("READ");
-    Nevents = rl->GetNumberOfEvents();
+     
+    Info("ReadNext","Reading Event %d",fCurrentEvent);
+    
+    fRunLoader->GetEvent(fCurrentEvent);
     
-    for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
+    AliStack* stack = fRunLoader->Stack();
+    if (!stack)
      {
-//      cout<<"processing event "<<currentEvent<<" in current dir."<<endl;
-      rl->GetEvent(currentEvent);
-      AliStack* stack = rl->Stack();
-      if (!stack)
-       {
-         Error("Read","Can not get stack for event %d",currentEvent);
-         continue;
-       }
-      Int_t npart = stack->GetNtrack();
-      Int_t nnn = 0;
-      for (Int_t i = 0;i<npart; i++)
-       {
-         
+       Error("ReadNext","Can not get stack for event %d",fCurrentEvent);
+       continue;
+     }
+    Int_t npart = stack->GetNtrack();
+    for (Int_t i = 0;i<npart; i++)
+      {
          TParticle * p = stack->Particle(i);
 //         if (p->GetFirstMother() >= 0) continue; do not apply with pythia etc
          
@@ -138,44 +109,56 @@ Int_t AliHBTReaderKineTree::Read(AliHBTRun* particles, AliHBTRun* /*tracks*/)
          AliHBTParticle* part = new AliHBTParticle(*p,i);
          if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
                                                   //if it does not delete it and take next good track
-         particles->AddParticle(totalNevents,part);//put track and particle on the run
-
-         if ( (nnn%100) == 0) 
-          {
-            cout<<nnn<<"\r";
-            fflush(0);
-          }
-         nnn++;
-       }
-      Info("Read","Total read %d",nnn);
-      totalNevents++;
-     }
-     delete rl;
-     currentdir++;
-  }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
- fIsRead = kTRUE;
- return 0;
+         fParticlesEvent->AddParticle(part);//put particle in event
+      }
+    Info("ReadNext","Read %d particles from event %d (event %d in dir %d).",
+                     fParticlesEvent->GetNumberOfParticles(),
+                     fNEventsRead,fCurrentEvent,fCurrentDir);
+      
+    fCurrentEvent++;
+    fNEventsRead++;
+    return 0;
+  }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
+  
+ return 1;
 }
+/**********************************************************/
 
-AliRunLoader* AliHBTReaderKineTree::OpenFile(Int_t n)
+Int_t AliHBTReaderKineTree::OpenNextFile()
 {
 //opens file with kine tree
-
- const TString& dirname = GetDirName(n);
+ Info("OpenNextFile","________________________________________________________");
+ const TString& dirname = GetDirName(fCurrentEvent);
  if (dirname == "")
   {
-   Error("OpenFiles","Can not get directory name");
-   return 0x0;
+   Error("OpenNextFile","Can not get directory name");
+   return 1;
   }
  TString filename = dirname +"/"+ fFileName;
 
AliRunLoader *ret = AliRunLoader::Open(filename.Data(),fgkEventFolderName,"READ"); 
fRunLoader = AliRunLoader::Open(filename.Data(),fgkEventFolderName,"READ"); 
 
- if ( ret == 0x0)
+ if ( fRunLoader == 0x0)
   {
-    Error("OpenFiles","Can't open session from file %s",filename.Data());
+    Error("OpenNextFile","Can't open session from file %s",filename.Data());
     return 0x0;
   }
- return ret;
+  
+ if (fRunLoader->GetNumberOfEvents() <= 0)
+  {
+    Error("OpenNextFile","There is no events in this directory.");
+    delete fRunLoader;
+    fRunLoader = 0x0;
+    return 1;
+  }
+  
+ if (fRunLoader->LoadKinematics())
+  {
+    Error("OpenNextFile","Error occured while loading kinematics.");
+    return 1;
+  }
+  
+ fCurrentEvent = 0;
+ return 0;
 }
index 9c08ba769b30caae7d9559a023757b121a029498..0c47fe1561e1a11b620b4fedcdf91ffefa08ec31 100644 (file)
@@ -16,30 +16,26 @@ class AliHBTReaderKineTree: public AliHBTReader
     AliHBTReaderKineTree(TString&);
     AliHBTReaderKineTree(TObjArray*,const Char_t *filename="galice.root");
 
-    virtual ~AliHBTReaderKineTree(){}
+    virtual ~AliHBTReaderKineTree();
     
-    Int_t        Read(AliHBTRun* particles, AliHBTRun* /*tracks*/);//reads tracks and particles and puts them in runs
-
-    AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
-    AliHBTEvent* GetTrackEvent(Int_t){return 0x0;}//returns pointer to event with particles
-    Int_t        GetNumberOfPartEvents();//returns number of particle events
-    Int_t        GetNumberOfTrackEvents(){return 0;}//returns number of track events
-
+    void          Rewind();
+    
+    Bool_t        ReadsTracks() const {return kFALSE;}
+    Bool_t        ReadsParticles() const {return kTRUE;}
     
    protected:
-    TString       fFileName;
-    AliHBTRun*    fParticles; //!simulated particles
-
-    AliRunLoader* OpenFile(Int_t);
-
-    Bool_t fIsRead;//!flag indicating if the data are already read
+    Int_t         ReadNext();//reads tracks and particles and puts them in runs
+    Int_t         OpenNextFile();
+   
+    TString       fFileName;//file name 
+    AliRunLoader* fRunLoader;
     
     static const TString fgkEventFolderName;
     
    private:
    
    public:
-     ClassDef(AliHBTReaderKineTree,1)
+     ClassDef(AliHBTReaderKineTree,2)
  };
 
 #endif
index a07dd50f5f1936fdb0161b0fd58c1af0c640fd07..c2aa96289b71d0e588d68b76b1c99096c581ff92 100644 (file)
@@ -4,6 +4,8 @@
 #include <TFile.h>
 #include <TParticle.h>
 
+#include <Varargs.h>
+
 #include <AliRun.h>
 #include <AliLoader.h>
 #include <AliStack.h>
 #include <AliTPCtrack.h>
 #include <AliTPCParam.h>
 #include <AliTPCtracker.h>
+#include <AliTPCLoader.h>
 
 #include "AliHBTRun.h"
 #include "AliHBTEvent.h"
 #include "AliHBTParticle.h"
 #include "AliHBTParticleCut.h"
 
-
 ClassImp(AliHBTReaderTPC)
 //______________________________________________
 //
@@ -28,348 +30,303 @@ ClassImp(AliHBTReaderTPC)
 //
 //more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
 //Piotr.Skowronski@cern.ch
-AliHBTReaderTPC::AliHBTReaderTPC():fFileName("galice.root")
+AliHBTReaderTPC::AliHBTReaderTPC():
+ fFileName("galice.root"),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
 {
   //constructor, 
   //Defaults:
   //  galicefilename = ""  - this means: Do not open gAlice file - 
   //                         just leave the global pointer untouched
   
-  fParticles = 0x0;
-  fTracks    = 0x0;
-  fIsRead = kFALSE;
 }
 
-AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):fFileName(galicefilename)
+AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
 {
   //constructor, 
   //Defaults:
   //  galicefilename = ""  - this means: Do not open gAlice file - 
   //                         just leave the global pointer untouched
   
-  fParticles = new AliHBTRun();
-  fTracks    = new AliHBTRun();
-  fIsRead = kFALSE;
 }
 /********************************************************************/
 AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
-                  AliHBTReader(dirs), fFileName(galicefilename)
-
+ AliHBTReader(dirs), 
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
 {
   //constructor, 
   //Defaults:
   //  galicefilename = ""  - this means: Do not open gAlice file - 
   //                         just leave the global pointer untached
   
-  fParticles = new AliHBTRun();
-  fTracks    = new AliHBTRun();
-  fIsRead = kFALSE;
 }
 /********************************************************************/
 
 AliHBTReaderTPC::~AliHBTReaderTPC()
- {
+{
  //desctructor
-   delete fParticles;
-   delete fTracks;
- }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
-   if (!fIsRead)
-    {
-     if (fParticles == 0x0) fParticles = new AliHBTRun();
-     if (fTracks == 0x0) fTracks    = new AliHBTRun();
-      
-     if(Read(fParticles,fTracks))
-      {
-        Error("GetParticleEvent","Error in reading");
-        return 0x0;
-      }
-    }
-   return (fParticles)?fParticles->GetEvent(n):0x0;
- }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
- {
- //returns Nth event with reconstructed tracks
-   if (!fIsRead) 
-    {
-      if (fParticles == 0x0) fParticles = new AliHBTRun();
-      if (fTracks == 0x0) fTracks    = new AliHBTRun();
-      if(Read(fParticles,fTracks))
-       {
-         Error("GetTrackEvent","Error in reading");
-         return 0x0;
-       }
-    }
-   return (fTracks)?fTracks->GetEvent(n):0x0;
- }
-/********************************************************************/
-
-Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
- {
- //returns number of events of particles
-   if (!fIsRead) 
-    {
-      if (fParticles == 0x0) fParticles = new AliHBTRun();
-      if (fTracks == 0x0) fTracks    = new AliHBTRun();
-      if ( Read(fParticles,fTracks))
-       {
-         Error("GetNumberOfPartEvents","Error in reading");
-         return 0;
-       }
-    }
-   return (fParticles)?fParticles->GetNumberOfEvents():0;
- }
-
+   delete fRunLoader;
+}
 /********************************************************************/
-Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
-  if (!fIsRead)
-   {
-    if (fParticles == 0x0) fParticles = new AliHBTRun();
-    if (fTracks == 0x0) fTracks    = new AliHBTRun();
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetNumberOfTrackEvents","Error in reading");
-       return 0;
-     }
-    }
-  return (fTracks)?fTracks->GetNumberOfEvents():0;
- }
+void AliHBTReaderTPC::Rewind()
+{
+  delete fRunLoader;
+  fRunLoader = 0x0;
+  fCurrentDir = 0;
+  fNEventsRead= 0;
+}
 /********************************************************************/
 
-
-Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
+Int_t AliHBTReaderTPC::ReadNext()
  {
  //reads data and puts put to the particles and tracks objects
  //reurns 0 if everything is OK
  //
   Info("Read","");
-  Int_t Nevents = 0;
-  Int_t totalNevents = 0;
 
-  if (!particles) //check if an object is instatiated
-   {
-     Error("Read"," particles object must instatiated before passing it to the reader");
-   }
-  if (!tracks)  //check if an object is instatiated
-   {
-     Error("Read"," tracks object must instatiated before passing it to the reader");
-   }
-  particles->Reset();//clear runs == delete all old events
-  tracks->Reset();
  
   TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
   tarray->SetOwner(); //set the ownership of the objects it contains
                       //when array is is deleted or cleared all objects 
                       //that it contains are deleted
-  Int_t currentdir = 0;
 
-  Int_t Ndirs;
-  if (fDirs) //if array with directories is supplied by user
-   {
-     Ndirs = fDirs->GetEntries(); //get the number if directories
-   }
-  else
-   {
-     Ndirs = 0; //if the array is not supplied read only from current directory
-   }
+  if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
+  if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
+
+  fParticlesEvent->Reset();
+  fTracksEvent->Reset();
 
   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
    {
-    TString filename = GetDirName(currentdir);
-    if (filename.IsNull())
-     {
-       Error("Read","Can not get directory name");
-       return 4;
-     }
-    filename = filename +"/"+ fFileName;
-    AliRunLoader* rl = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
-    if( rl == 0x0)
+      
+    if (fRunLoader == 0x0) 
+      if (OpenNextSession()) continue;//directory counter is increased inside in case of error
+
+    if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
      {
-       Error("Read","Can not open session.");
-       currentdir++;
-       continue;
+       //read next directory
+       delete fRunLoader;//close current session
+       fRunLoader = 0x0;//assure pointer is null
+       fCurrentDir++;//go to next dir
+       continue;//directory counter is increased inside in case of error
      }
+     
+    Info("ReadNext","Reading Event %d",fCurrentEvent);
     
-    rl->LoadHeader();
-    rl->LoadKinematics();
-    AliLoader* tpcl = rl->GetLoader("TPCLoader");
-    
-    if ( tpcl== 0x0)
+    fRunLoader->GetEvent(fCurrentEvent);
+    TTree *tracktree = fTPCLoader->TreeT();//get the tree 
+    if (!tracktree) //check if we got the tree
+      {//if not return with error
+         Error("ReadNext","Can't get a tree with TPC tracks !\n"); 
+         continue;
+      }
+   
+    TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
+    if (!trackbranch) ////check if we got the branch
+      {//if not return with error
+        Error("ReadNext","Can't get a branch with TPC tracks !\n"); 
+        continue;
+      }
+    Int_t ntpctracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks 
+    Info("ReadNext","Found %d TPC tracks.",ntpctracks);
+    //Copy tracks to array
+
+    AliTPCtrack *iotrack=0;
+
+    for (Int_t i=0; i<ntpctracks; i++) //loop over all tpc tracks
      {
-       Error("Read","Exiting due to problems with opening files.");
-       currentdir++;
-       continue;
+       iotrack=new AliTPCtrack;   //create new tracks
+       trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
+       tracktree->GetEvent(i); //stream track i to the iotrack
+       tarray->AddLast(iotrack); //put the track in the array
      }
-    Nevents = rl->GetNumberOfEvents();
-    if (Nevents > 0)//check if tree E exists
+
+    Double_t xk;
+    Double_t par[5];
+    Float_t phi, lam, pt;//angles and transverse momentum
+    Int_t label; //label of the current track
+
+    AliStack* stack = fRunLoader->Stack();
+    if (stack == 0x0)
      {
-      Info("Read","________________________________________________________");
-      Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
-      rl->LoadgAlice();
-      Info("Read","Setting Magnetic Field: B=%fT",rl->GetAliRun()->Field()->SolenoidField());
-      AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
-      rl->UnloadgAlice();
-     }
-    else
-     {//if not return an error
-       Error("Read","Can not find Header tree (TreeE) in gAlice");
-       currentdir++;
+       Error("ReadNext","Can not get stack for current event",fCurrentEvent);
+       fCurrentEvent++;
        continue;
      }
+    stack->Particles();
+
+    for (Int_t i=0; i<ntpctracks; i++) //loop over all good tracks
+     { 
+       iotrack = (AliTPCtrack*)tarray->At(i);
+       label = iotrack->GetLabel();
+
+       if (label < 0) continue;
+
+       TParticle *p = (TParticle*)stack->Particle(label);
+
+       if(p == 0x0) continue; //if returned pointer is NULL
+       if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
+
+       if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
+                                   //if not take next partilce
+
+       AliHBTParticle* part = new AliHBTParticle(*p,i);
+       if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
+                                               //if it does not delete it and take next good track
+
+//       iotrack->PropagateTo(3.,0.0028,65.19);
+//       iotrack->PropagateToVertex(36.66,1.2e-3);//orig
+       iotrack->PropagateToVertex(50.,0.0353);
+       
+       iotrack->GetExternalParameters(xk,par);     //get properties of the track
+       phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
+       if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+       if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+       lam=par[3]; 
+       pt=1.0/TMath::Abs(par[4]);
+
+       Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
+       Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
+       Double_t tpz = pt * lam; //track z coordinate of momentum
+
+       Double_t mass = p->GetMass();
+       Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
+
+       AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+       if(Pass(track))//check if meets all criteria of any of our cuts
+                    //if it does not delete it and take next good track
+        { 
+          delete track;
+          delete part;
+          continue;
+        }
+       fParticlesEvent->AddParticle(part);
+       fTracksEvent->AddParticle(track);
+      }
+      
+    Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
+            fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
+            fNEventsRead,fCurrentEvent,fCurrentDir);
     
-   rl->CdGAFile();
-   AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
-   
-   if (!TPCParam) 
-    {
-     TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
-     if (!TPCParam) 
-      { 
-        Error("Read","TPC parameters have not been found !\n");
-        delete rl;
-        rl = 0x0;
-        currentdir++;
-        continue;
+    fCurrentEvent++;
+    fNEventsRead++;
+    delete tarray;
+    return 0;
+   }while(fCurrentDir < GetNumberOfDirs());
+
+  delete tarray;
+  return 1;
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderTPC::OpenNextSession()
+{
+  TString filename = GetDirName(fCurrentDir);
+  if (filename.IsNull())
+   {
+     DoOpenError("Can not get directory name");
+     return 1;
+   }
+  filename = filename +"/"+ fFileName;
+  fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
+  if( fRunLoader == 0x0)
+   {
+     DoOpenError("Can not open session.");
+     return 1;
+   }
+
+  fRunLoader->LoadHeader();
+  if (fRunLoader->GetNumberOfEvents() <= 0)
+   {
+     DoOpenError("There is no events in this directory.");
+     return 1;
+   }
+
+  if (fRunLoader->LoadKinematics())
+   {
+     DoOpenError("Error occured while loading kinematics.");
+     return 1;
+   }
+  fTPCLoader = (AliTPCLoader*)fRunLoader->GetLoader("TPCLoader");
+  if ( fTPCLoader == 0x0)
+   {
+     DoOpenError("Exiting due to problems with opening files.");
+     return 1;
+   }
+  Info("OpenNextSession","________________________________________________________");
+  Info("OpenNextSession","Found %d event(s) in directory %s",
+        fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
+  Float_t mf;
+  if (fUseMagFFromRun)
+   {
+     fRunLoader->LoadgAlice();
+     mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
+     Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
+     fRunLoader->UnloadgAlice();
+   }
+  else
+   {
+     Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
+     if (fMagneticField == 0x0)
+      {
+        Fatal("OpenNextSession","Magnetic field can not be 0.");
+        return 1;//pro forma
       }
-    }
+     mf = fMagneticField*10.;
+   }
+  AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
 
-    tpcl->LoadTracks();
-  
-    for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
-     {
-       Info("Read","Reading Event %d",currentEvent);
-       /**************************************/
-        /**************************************/
-         /**************************************/ 
-         rl->GetEvent(currentEvent);
-         TTree *tracktree = tpcl->TreeT();//get the tree 
-         if (!tracktree) //check if we got the tree
-          {//if not return with error
-            Error("Read","Can't get a tree with TPC tracks !\n"); 
-            continue;
-          }
-   
-         TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
-         if (!trackbranch) ////check if we got the branch
-          {//if not return with error
-            Error("Read","Can't get a branch with TPC tracks !\n"); 
-            continue;
-          }
-         Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks 
-         Info("Read","Found %d TPC tracks.",NTPCtracks);
-         //Copy tracks to array
-         
-         AliTPCtrack *iotrack=0;
-         
-printf("This method is not converted to the NewIO !\n"); //I.B.
-         //AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent,AliConfig::fgkDefaultEventFolderName);//create the tacker for this event
-         AliTPCtracker *tracker = new AliTPCtracker(TPCParam); //I.B.
-         if (!tracker) //check if it has created succeffuly
-          {//if not return with error
-            Error("Read","Can't get a tracker !\n"); 
-            continue;
-          }
-         tracker->LoadClusters(0);//I.Belikov, "0" must be a pointer to a tree
-   
-         for (Int_t i=0; i<NTPCtracks; i++) //loop over all tpc tracks
-          {
-            iotrack=new AliTPCtrack;   //create new tracks
-            trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
-            tracktree->GetEvent(i); //stream track i to the iotrack
-            tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
-                                             //which is the label of corresponding simulated particle 
-            tarray->AddLast(iotrack); //put the track in the array
-          }
-         
-         delete tracker; //delete tracker
-         
-         tracker = 0x0;
-         trackbranch = 0x0;
-         tracktree = 0x0;
-   
-         Double_t xk;
-         Double_t par[5];
-         Float_t phi, lam, pt;//angles and transverse momentum
-         Int_t label; //label of the current track
-
-         rl->Stack()->Particles();
-         
-         for (Int_t i=0; i<NTPCtracks; i++) //loop over all good tracks
-          { 
-            iotrack = (AliTPCtrack*)tarray->At(i);
-            label = iotrack->GetLabel();
-
-            if (label < 0) continue;
-            
-            TParticle *p = (TParticle*)rl->Stack()->Particle(label);
-             
-            if(p == 0x0) continue; //if returned pointer is NULL
-            if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
-           
-            if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
-                                        //if not take next partilce
-            
-            AliHBTParticle* part = new AliHBTParticle(*p,i);
-            if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
-                                                    //if it does not delete it and take next good track
-         
-            iotrack->PropagateToVertex();
-
-            iotrack->GetExternalParameters(xk,par);     //get properties of the track
-            phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
-            if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
-            if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
-            lam=par[3]; 
-            pt=1.0/TMath::Abs(par[4]);
-            
-            Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
-            Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
-            Double_t tpz = pt * lam; //track z coordinate of momentum
-            
-            Double_t mass = p->GetMass();
-            Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
-            
-            AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
-            if(Pass(track))//check if meets all criteria of any of our cuts
-                         //if it does not delete it and take next good track
-             { 
-               delete track;
-               delete part;
-               continue;
-             }
-            particles->AddParticle(totalNevents,part);//put track and particle on the run
-            tracks->AddParticle(totalNevents,track);
-
-          }
-         tarray->Clear(); //clear the array
-         
-        /**************************************/
-       /**************************************/
-      /**************************************/  
-     totalNevents++;
-    }
-  
-    //save environment (resouces) --
-    //clean your place after the work
-    delete rl;
-    currentdir++;
-   }while(currentdir < Ndirs);
 
-  delete tarray;
-  fIsRead = kTRUE;
+  fRunLoader->CdGAFile();
+  AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
+  if (!TPCParam) 
+   {
+    TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
+    if (!TPCParam) 
+     { 
+       DoOpenError("TPC parameters have not been found !\n");
+       return 1;
+     }
+   }
+
+  if (fTPCLoader->LoadTracks())
+   {
+     DoOpenError("Error occured while loading TPC tracks.");
+     return 1;
+   }
   
+  fCurrentEvent = 0;
   return 0;
- }
+}
+/********************************************************************/
 
+void AliHBTReaderTPC::DoOpenError( const char *va_(fmt), ...)
+{
+  // Does error display and clean-up in case error caught on Open Next Session
+
+   va_list ap;
+   va_start(ap,va_(fmt));
+   Error("OpenNextSession", va_(fmt), ap);
+   va_end(ap);
+   
+   delete fRunLoader;
+   fRunLoader = 0x0;
+   fTPCLoader = 0x0;
+   fCurrentDir++;
+}
 
 /********************************************************************/
 /********************************************************************/
index 43b64b531aba81c139dd475e6af85bab4166e0e9..8c6e17eece0fa956618dca0b01f758444cb8681f 100644 (file)
@@ -12,6 +12,8 @@
 
 #include <TString.h>
 class TFile;
+class AliRunLoader;
+class AliTPCLoader;
 
 class AliHBTReaderTPC: public AliHBTReader
 {
@@ -22,25 +24,30 @@ class AliHBTReaderTPC: public AliHBTReader
 
     virtual ~AliHBTReaderTPC();
     
-    Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+    void          Rewind();
     
-    AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
-    AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles 
-    Int_t GetNumberOfPartEvents();//returns number of particle events
-    Int_t GetNumberOfTrackEvents();//returns number of track events
+    Bool_t        ReadsTracks() const {return kTRUE;}
+    Bool_t        ReadsParticles() const {return kTRUE;}
+    
+    void          SetMagneticField(Float_t mf){fMagneticField=mf;}
+    void          UseMagneticFieldFromRun(Bool_t flag = kTRUE){fUseMagFFromRun=flag;}
     
   protected:
     //in the future this class is will read global tracking
+    Int_t         ReadNext();
+    Int_t         OpenNextSession();
+    void          DoOpenError(const char* msgfmt, ...);
     
-    AliHBTRun* fParticles; //!simulated particles
-    AliHBTRun* fTracks; //!reconstructed tracks (particles)
-
-    TString fFileName;//name of the file with galice.root
+    TString       fFileName;//name of the file with galice.root
+    AliRunLoader* fRunLoader;//!RL
+    AliTPCLoader* fTPCLoader;//!TPCLoader
+    Float_t       fMagneticField;//magnetic field value that was enforced while reading
+    Bool_t        fUseMagFFromRun;//flag indicating if using field specified in gAlice (kTRUE)
+                               // or enforece other defined by fMagneticField
 
-    Bool_t fIsRead;//!flag indicating if the data are already read
   private:
   public:
-    ClassDef(AliHBTReaderTPC,2)
+    ClassDef(AliHBTReaderTPC,3)
 };
 
 
index d9e156e9c70f44444070a32692c047a7e3ac456e..04fcec2cf1001d96042b5cba0a93be4e9cc54a06 100644 (file)
@@ -66,5 +66,20 @@ void AliHBTRun::AddParticle(Int_t event, Int_t pdg, Int_t idx,
  if(!GetEvent(event))  fEvents->AddAtAndExpand(new AliHBTEvent, event);
  GetEvent(event)->AddParticle(pdg,idx,px,py,pz,etot,vx,vy,vz,time);
 }
+/**************************************************************************/ 
 
+void AliHBTRun::SetEvent(Int_t number, AliHBTEvent* event)
+{
+  //adds an event to the run
+  // makes an own copy of the event!
+  if (event == 0x0)
+   {
+     delete fEvents->RemoveAt(number);
+     return;
+   }
+  AliHBTEvent* ev = GetEvent(number);
+  if (ev) *ev = *event;
+  else fEvents->AddAtAndExpand(new AliHBTEvent(*event), number);
+}
 /**************************************************************************/ 
+
index 67fe48fbcc76243d4f5c11777a80b2c1603b5efb..672633b3e78c4ed4f46c39bd3c4489a907d2c5eb 100644 (file)
@@ -35,6 +35,7 @@ class AliHBTRun: public TObject
                                 Double_t px, Double_t py, Double_t pz, Double_t etot,
                                 Double_t vx, Double_t vy, Double_t vz, Double_t time); 
     
+    void            SetEvent(Int_t number, AliHBTEvent* event);
     AliHBTParticle* GetParticle(Int_t event, Int_t n); //returns nth particle from event
     AliHBTEvent*    GetEvent(Int_t event) const; //returns AliHBTEvent number "event"
     
index 35470ca1df1d91b566651206f723622da172b131..464de7325f773f75f8cdd814e08b1a1f61d8491c 100644 (file)
@@ -28,8 +28,8 @@ class AliHBTTwoTrackEffFctn3D: public AliHBTOnePairFctn3D
     AliHBTTwoTrackEffFctn3D();
     virtual ~AliHBTTwoTrackEffFctn3D(){}
 
-    void ProcessSameEventParticles(AliHBTPair* pair){}
-    void ProcessDiffEventParticles(AliHBTPair* pair){}
+    void ProcessSameEventParticles(AliHBTPair* /*pair*/){}
+    void ProcessDiffEventParticles(AliHBTPair* /*pair*/){}
 
   protected:
     void GetValues(AliHBTPair*,Double_t&, Double_t&,Double_t&);
index a2bb032f3a30d219ad0c93ef3d6309c6e001dce0..292100b66216a4ca02283a758c74e6702bea214a 100644 (file)
@@ -60,6 +60,7 @@ void AliHBTWriteInternFormat(Option_t* datatype, Option_t* processopt="TracksAnd
    
   AliHBTReader* reader;
   Int_t kine = strcmp(datatype,"Kine");
+  Int_t ESD = strcmp(datatype,"ESD");
   Int_t TPC = strcmp(datatype,"TPC");
   Int_t ITSv1 = strcmp(datatype,"ITSv1");
   Int_t ITSv2 = strcmp(datatype,"ITSv2");
@@ -70,10 +71,16 @@ void AliHBTWriteInternFormat(Option_t* datatype, Option_t* processopt="TracksAnd
     reader = new AliHBTReaderKineTree();
     processopt="Particles"; //this reader by definition reads only simulated particles
    }
+  else if(!ESD)
+   {
+    AliHBTReaderESD* esdreader = new AliHBTReaderESD();
+    esdreader->ReadParticles(kTRUE);
+    reader = esdreader;
+   }
   else if(!TPC)
    {
     cout<<"AliHBTWriteInternFormat.C: Creating Reader TPC .....\n";
-    reader = new AliHBTReaderTPC("tpc.tracks.root","tpc.clusters.root");
+    reader = new AliHBTReaderTPC();
     cout<<"AliHBTWriteInternFormat.C: ..... Created\n";
    }
   else if(!ITSv1)
@@ -83,12 +90,12 @@ void AliHBTWriteInternFormat(Option_t* datatype, Option_t* processopt="TracksAnd
   else if(!ITSv2)
    {
     cout<<"AliHBTWriteInternFormat.C: Creating Reader ITSv2 .....\n";
-    reader = new AliHBTReaderITSv2("its.tracksv2.root","its.clustersv2.root");
+    reader = new AliHBTReaderITSv2();
     cout<<"AliHBTWriteInternFormat.C: ..... Created\n";
    }
   else if(!intern)
    {
-    reader = new AliHBTReaderInternal("TPC.root");
+    reader = new AliHBTReaderInternal("data.root");
    }
   else
    {
@@ -115,7 +122,7 @@ void AliHBTWriteInternFormat(Option_t* datatype, Option_t* processopt="TracksAnd
 
    cout<<"AliHBTWriteInternFormat.C:   P R O C S E S S I N G .....\n\n";
    AliHBTReaderInternal::Write(reader,outfile);
-   cout<<"\n\nAliHBTWriteInternFormat.C:   F I N I S H E D";
+   cout<<"\n\nAliHBTWriteInternFormat.C:   F I N I S H E D\n";
    
    if (dirs) delete dirs;
    delete reader;
index c4b51afd2f3d00c5089b3e357f22bc40fd07e077..d4d25fafc4bb8bac9c79d7262479248426e037eb 100644 (file)
@@ -8,6 +8,7 @@
 #pragma link C++ class AliHBTParticle+;
 #pragma link C++ class AliHBTEvent+;
 #pragma link C++ class AliHBTRun+;
+#pragma link C++ class AliHBTEventBuffer+;
 #pragma link C++ class AliHBTFunction+;
 #pragma link C++ class AliHBTMonitorFunction+;
 
index f64184ed6be1ef97f1a64ca2b3e2e1c33e068f40..c38145c12c96ce005d54536619999734c48caed7 100644 (file)
@@ -12,7 +12,7 @@ AliHBTQResolutionFctns.cxx      AliHBTQDistributionFctns.cxx \
 AliHBTMonDistributionFctns.cxx  AliHBTMonResolutionFctns.cxx \
 AliHBTLLWeights.cxx             AliHBTLLWeightFctn.cxx \
 AliHBTLLWeightsPID.cxx          AliHBTLLWeightTheorFctn.cxx \
-AliHBTPositionRandomizer.cxx
+AliHBTPositionRandomizer.cxx AliHBTEventBuffer.cxx
 
 FSRCS   = fsiini.F  fsiw.F  led_bldata.F  ltran12.F