This commit was generated by cvs2svn to compensate for changes in r4472,
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 12 Nov 2001 08:12:54 +0000 (08:12 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 12 Nov 2001 08:12:54 +0000 (08:12 +0000)
which included commits to RCS files with non-trunk default branches.

34 files changed:
HBTAN/.rootrc [new file with mode: 0644]
HBTAN/AliHBTAnalysis.cxx [new file with mode: 0644]
HBTAN/AliHBTAnalysis.h [new file with mode: 0644]
HBTAN/AliHBTCorrelFctn.cxx [new file with mode: 0644]
HBTAN/AliHBTCorrelFctn.h [new file with mode: 0644]
HBTAN/AliHBTEvent.cxx [new file with mode: 0644]
HBTAN/AliHBTEvent.h [new file with mode: 0644]
HBTAN/AliHBTFunction.cxx [new file with mode: 0644]
HBTAN/AliHBTFunction.h [new file with mode: 0644]
HBTAN/AliHBTPair.cxx [new file with mode: 0644]
HBTAN/AliHBTPair.h [new file with mode: 0644]
HBTAN/AliHBTPairCut.cxx [new file with mode: 0644]
HBTAN/AliHBTPairCut.h [new file with mode: 0644]
HBTAN/AliHBTParticle.cxx [new file with mode: 0644]
HBTAN/AliHBTParticle.h [new file with mode: 0644]
HBTAN/AliHBTParticleCut.cxx [new file with mode: 0644]
HBTAN/AliHBTParticleCut.h [new file with mode: 0644]
HBTAN/AliHBTQResolutionFctns.cxx [new file with mode: 0644]
HBTAN/AliHBTQResolutionFctns.h [new file with mode: 0644]
HBTAN/AliHBTReader.cxx [new file with mode: 0644]
HBTAN/AliHBTReader.h [new file with mode: 0644]
HBTAN/AliHBTReaderITSv1.cxx [new file with mode: 0644]
HBTAN/AliHBTReaderITSv1.h [new file with mode: 0644]
HBTAN/AliHBTReaderITSv2.cxx [new file with mode: 0644]
HBTAN/AliHBTReaderITSv2.h [new file with mode: 0644]
HBTAN/AliHBTReaderPPprod.cxx [new file with mode: 0644]
HBTAN/AliHBTReaderPPprod.h [new file with mode: 0644]
HBTAN/AliHBTReaderTPC.cxx [new file with mode: 0644]
HBTAN/AliHBTReaderTPC.h [new file with mode: 0644]
HBTAN/AliHBTRun.cxx [new file with mode: 0644]
HBTAN/AliHBTRun.h [new file with mode: 0644]
HBTAN/HBTAnalysisLinkDef.h [new file with mode: 0644]
HBTAN/Makefile [new file with mode: 0644]
HBTAN/libHBTAnalysis.pkg [new file with mode: 0644]

diff --git a/HBTAN/.rootrc b/HBTAN/.rootrc
new file mode 100644 (file)
index 0000000..413ae4c
--- /dev/null
@@ -0,0 +1,7 @@
+
+Unix.*.Root.MacroPath:      .:$(ALICE_ROOT)/macros:$(ALICE_ROOT)
+
+Root.Html.OutputDir: $(ALICE_ROOT)/html/roothtml/HBTANALYSIS
+Root.Html.SourceDir: ./
+Root.Html.Author:  Piotr Krzysztof Skowronski
+Root.Html.Root:  http://root.cern.ch/root/html
diff --git a/HBTAN/AliHBTAnalysis.cxx b/HBTAN/AliHBTAnalysis.cxx
new file mode 100644 (file)
index 0000000..fa30c16
--- /dev/null
@@ -0,0 +1,634 @@
+
+#include "AliHBTAnalysis.h"
+
+#include <iostream.h>
+
+#include "AliHBTRun.h"
+#include "AliHBTReader.h"
+#include "AliHBTParticle.h"
+#include "AliHBTParticleCut.h"
+#include "AliHBTPair.h"
+#include "AliHBTPairCut.h"
+#include "AliHBTFunction.h"
+
+#include <TList.h>
+
+
+
+ClassImp(AliHBTAnalysis)
+
+const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
+const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
+
+AliHBTAnalysis::AliHBTAnalysis()
+ {
+   fReader = 0x0;
+   
+   fTrackFunctions = new AliHBTTwoPartFctn* [fgkFctnArraySize];
+   fParticleFunctions = new AliHBTTwoPartFctn* [fgkFctnArraySize];
+   fParticleAndTrackFunctions = new AliHBTFourPartFctn* [fgkFctnArraySize];
+   
+   fNTrackFunctions = 0;
+   fNParticleFunctions = 0;
+   fNParticleAndTrackFunctions = 0;
+   
+   fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
+   
+ }
+
+AliHBTAnalysis::~AliHBTAnalysis()
+ {
+ //destructor
+ //note that we do not delete functions itself
+ // they should be deleted by whom where created
+ //we only store pointers, and we use "new" only for pointers array
+   delete [] fTrackFunctions;
+   delete [] fParticleFunctions;
+   delete [] fParticleAndTrackFunctions;
+   
+   delete fPairCut; // always have an copy of an object - we create we dstroy
+ }
+
+/*************************************************************************************/ 
+void AliHBTAnalysis::Process(Option_t* option)
+{
+ //default option  = "TracksAndParticles"
+ //Main method of the HBT Analysis Package
+ //It triggers reading with the global cut (default is an empty cut)
+ //Than it checks options and data which are read
+ //if everything is OK, then it calls one of the looping methods
+ //depending on tfReaderhe option
+ //These methods differs on what thay are looping on
+ //
+ //        METHOD                           OPTION
+ //--------------------------------------------------------------------
+ //ProcessTracksAndParticles    -     "TracksAndParticles"
+ //     DEFAULT
+ //     it loops over both, tracks(reconstructed) and particles(simulated)
+ //     all function gethered in all 3 lists are called for each (double)pair
+ //                             
+ //ProcessTracks                -          "Tracks" 
+ //     it loops only on tracks(reconstructed),
+ //     functions ONLY from fTrackFunctions list are called
+ //
+ //ProcessParticles             -         "Particles"
+ //     it loops only on particles(simulated),
+ //     functions ONLY from fParticleAndTrackFunctions list are called
+ //
+ //
+ if (!fReader) 
+  {
+   Error("AliHBTAnalysis::Process","The reader is not set");
+   return;
+  }
+ const char *oT = strstr(option,"Tracks");
+ const char *oP = strstr(option,"Particles");
+ if(oT && oP)
+  { 
+    if (fReader->GetNumberOfPartEvents() <1)
+     {
+       Error("AliHBTAnalysis::Process","There is no Particles. Maybe change the option?");
+       return;
+     }
+    if (fReader->GetNumberOfTrackEvents() <1)
+     {
+       Error("AliHBTAnalysis::Process","There is no Tracks. Maybe change the option?");
+       return;
+     }
+    
+    if ( RunCoherencyCheck() )
+      {
+        Error("AliHBTAnalysis::Process",
+              "Coherency check not passed. Maybe change the option?\n");
+        return;
+      }
+    ProcessTracksAndParticles();
+    return;
+  }
+
+ if(oT)
+  {
+    ProcessTracks();
+    return;
+  }
+ if(oP)
+  {
+    ProcessParticles();
+    return;
+  }
+}
+
+/*************************************************************************************/ 
+void AliHBTAnalysis::ProcessTracksAndParticles()
+{
+
+//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+  
+  
+  AliHBTParticle * part1, * part2;
+  AliHBTParticle * track1, * track2;
+  
+  AliHBTEvent * trackEvent, *partEvent;
+  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    ********/
+  /***************************************/
+  AliHBTPair * trackpair = new AliHBTPair();
+  AliHBTPair * partpair = new AliHBTPair();
+  
+  for (Int_t i = 0;i<Nev;i++)
+    {
+      partEvent= fReader->GetParticleEvent(i);
+      trackEvent = fReader->GetTrackEvent(i);
+      
+      if (!partEvent) continue;
+      
+      //N = 0;
+      
+      for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+       {
+         if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+
+         part1= partEvent->GetParticle(j);
+         track1= trackEvent->GetParticle(j);   
+        
+         for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+          {
+            part2= partEvent->GetParticle(k);
+            partpair->SetParticles(part1,part2);
+           
+            track2= trackEvent->GetParticle(k);        
+            trackpair->SetParticles(track1,track2);
+
+            if(fPairCut->Pass(partpair) || (fPairCut->Pass(trackpair))) //check pair cut
+              { //do not meets crietria of the pair cut, try with swaped pairs
+                if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) ) 
+                  continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
+                else 
+                 { //swaped pair meets all the criteria
+                   partpair = partpair->GetSwapedPair();
+                   trackpair = trackpair->GetSwapedPair();
+                 }
+              }
+            UInt_t ii;
+            for(ii = 0;ii<fNParticleFunctions;ii++)
+                   fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
+                
+            for(ii = 0;ii<fNTrackFunctions;ii++)
+                   fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
+                 
+            for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+                   fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
+           }
+       }
+    }
+
+  /***************************************/
+  /***** Filling denominators    *********/
+  /***************************************/
+  for (Int_t i = 0;i<Nev;i++)   //In each event ....
+    {
+      
+      partEvent= fReader->GetParticleEvent(i);
+      if (!partEvent) continue;
+      
+      trackEvent = fReader->GetTrackEvent(i); 
+      
+//      N=0;
+      
+      for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
+       {
+//         if (N>MAXCOMB) break;
+           
+           part1= partEvent->GetParticle(j);
+
+           track1= trackEvent->GetParticle(j);
+//         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
+           Int_t NNN;
+  
+           if ( (i+2) < Nev) NNN = i+2;
+           else NNN = Nev;
+           for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
+            {
+             
+              partEvent2= fReader->GetParticleEvent(k);
+              if (!partEvent2) continue;
+              
+              trackEvent2 = fReader->GetTrackEvent(k);
+             
+             if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
+            
+             for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
+              {
+               
+                // if (N>MAXCOMB) break;
+                
+                part2= partEvent2->GetParticle(l);
+                partpair->SetParticles(part1,part2);
+
+                track2= trackEvent2->GetParticle(l);
+                trackpair->SetParticles(track1,track2);
+
+                if( (fPairCut->Pass(partpair)) || (fPairCut->Pass(trackpair)) ) //check pair cut
+                  { //do not meets crietria of the 
+                    if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
+          continue;
+       else 
+        {
+          partpair = partpair->GetSwapedPair();
+          trackpair = trackpair->GetSwapedPair();
+        }
+                  }
+                UInt_t ii;
+                for(ii = 0;ii<fNParticleFunctions;ii++)
+                  fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
+                 
+                for(ii = 0;ii<fNTrackFunctions;ii++)
+                  fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
+                 
+                for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+                  fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
+           
+               
+              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
+            }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
+       }
+    
+    } 
+
+  /***************************************/
+  
+   
+} 
+/*************************************************************************************/
+
+void AliHBTAnalysis::ProcessTracks()
+{
+  //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+  AliHBTParticle * track1, * track2;
+  
+  AliHBTEvent * trackEvent;
+  AliHBTEvent * trackEvent2;
+  
+//  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    ********/
+  /***************************************/
+  AliHBTPair * trackpair = new AliHBTPair();
+  
+  for (Int_t i = 0;i<Nev;i++)
+    {
+      trackEvent = fReader->GetTrackEvent(i);
+      if (!trackEvent) continue;
+      //N = 0;
+      
+      for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
+       {
+         if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+
+         track1= trackEvent->GetParticle(j);   
+        
+         for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
+          {
+            track2= trackEvent->GetParticle(k);        
+            trackpair->SetParticles(track1,track2);
+            if(fPairCut->Pass(trackpair)) //check pair cut
+              { //do not meets crietria of the 
+                if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
+                else trackpair = trackpair->GetSwapedPair();
+              }
+            
+            UInt_t ii;
+                
+            for(ii = 0;ii<fNTrackFunctions;ii++)
+                   fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
+                 
+           
+           }
+       }
+    }
+
+  /***************************************/
+  /***** Filling diff histogram *********/
+  /***************************************/
+  for (Int_t i = 0;i<Nev;i++)   //In each event ....
+    {
+      trackEvent = fReader->GetTrackEvent(i);
+      if (!trackEvent) continue;
+//      N=0;
+      
+      for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
+       {
+//         if (N>MAXCOMB) break;
+           
+           track1= trackEvent->GetParticle(j);
+//         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
+           Int_t NNN;
+  
+           if ( (i+2) < Nev) NNN = i+2;
+           else NNN = Nev;
+           for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
+            {
+             
+             trackEvent2 = fReader->GetTrackEvent(k);
+             if (!trackEvent2) continue;
+             
+             if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
+            
+             for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++)   //  ... on all particles
+              {
+               
+                // if (N>MAXCOMB) break;
+                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 trackpair = trackpair->GetSwapedPair();
+                  }
+                UInt_t ii;
+                for(ii = 0;ii<fNTrackFunctions;ii++)
+                  fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
+               
+              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
+            }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
+       }
+    
+    } 
+
+  /***************************************/
+  
+
+}
+
+/*************************************************************************************/
+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 * partEvent2;
+  
+//  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
+  
+  Int_t Nev = fReader->GetNumberOfPartEvents();
+  
+  /***************************************/
+  /******   Looping same events   ********/
+  /******   filling numerators    ********/
+  /***************************************/
+  AliHBTPair * partpair = new AliHBTPair();
+  
+  for (Int_t i = 0;i<Nev;i++)
+    {
+      partEvent= fReader->GetParticleEvent(i);
+      if (!partEvent) continue;
+      //N = 0;
+      
+      for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+       {
+         if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+
+         part1= partEvent->GetParticle(j);
+        
+         for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+          {
+            part2= partEvent->GetParticle(k);
+            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
+                   partpair = partpair->GetSwapedPair();
+                 }
+              }
+
+            UInt_t ii;
+                
+            for(ii = 0;ii<fNParticleFunctions;ii++)
+                   fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
+           }
+       }
+    }
+
+  /***************************************/
+  /***** Filling diff histogram *********/
+  /***************************************/
+  for (Int_t i = 0;i<Nev;i++)   //In each event ....
+    {
+      partEvent= fReader->GetParticleEvent(i);
+      if (!partEvent) continue;
+
+//      N=0;
+      
+      for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
+       {
+//         if (N>MAXCOMB) break;
+           
+           part1= partEvent->GetParticle(j);
+//         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
+           Int_t NNN;
+  
+           if ( (i+2) < Nev) NNN = i+2; //loop over next event
+           else NNN = Nev;
+           for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
+            {
+             
+             partEvent2= fReader->GetParticleEvent(k);
+             if (!partEvent2) continue;
+             
+             if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
+            
+             for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
+              {
+               
+                // if (N>MAXCOMB) break;
+                part2= partEvent2->GetParticle(l);
+                partpair->SetParticles(part1,part2);
+                
+                if(fPairCut->Pass(partpair)) //check pair cut
+                  { //do not meets crietria of the 
+                    if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
+       else partpair = partpair->GetSwapedPair();
+                  }
+                UInt_t ii;
+                for(ii = 0;ii<fNParticleFunctions;ii++)
+                  fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
+               
+              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
+            }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
+       }
+    
+    } 
+
+  /***************************************/
+  
+
+}
+
+/*************************************************************************************/
+
+
+void AliHBTAnalysis::Write()
+{
+  UInt_t ii;
+  for(ii = 0;ii<fNParticleFunctions;ii++)
+    fParticleFunctions[ii]->Write();
+                 
+  for(ii = 0;ii<fNTrackFunctions;ii++)
+    fTrackFunctions[ii]->Write();
+                 
+  for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+   fParticleAndTrackFunctions[ii]->Write();
+}
+/*************************************************************************************/
+void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
+{
+  if (cut == 0x0)
+   {
+     Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
+   }
+  delete fPairCut;
+  fPairCut = (AliHBTPairCut*)cut->Clone();
+}
+
+/*************************************************************************************/
+void AliHBTAnalysis::AddTrackFunction(AliHBTTwoPartFctn* f)
+{
+  if (f == 0x0) return;
+  if (fNTrackFunctions == fgkFctnArraySize)
+   {
+    Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
+   }
+  fTrackFunctions[fNTrackFunctions] = f;
+  fNTrackFunctions++;
+}
+void AliHBTAnalysis::AddParticleFunction(AliHBTTwoPartFctn* f)
+{
+  if (f == 0x0) return;
+  
+  if (fNParticleFunctions == fgkFctnArraySize)
+   {
+    Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
+   }
+  fParticleFunctions[fNParticleFunctions] = f;
+  fNParticleFunctions++;
+  
+  
+}
+void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTFourPartFctn* f)
+{
+  if (f == 0x0) return;
+  if (fNParticleAndTrackFunctions == fgkFctnArraySize)
+   {
+    Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
+   }  
+  fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
+  fNParticleAndTrackFunctions++;
+}
+
+
+/*************************************************************************************/ 
+
+
+/*************************************************************************************/  
+Bool_t AliHBTAnalysis::RunCoherencyCheck()
+{
+ //Checks if both HBTRuns are similar
+ //return true if error found
+ //if they seem to be OK return false
+ Int_t i;  
+ cout<<"Checking HBT Runs Coherency"<<endl;
+ cout<<"Number of events ..."; fflush(stdout);
+ if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same  number of events
+  {
+    cout<<"OK. "<<fReader->GetNumberOfTrackEvents()<<"  found"<<endl;
+  }
+ else
+  { //if not the same - ERROR
+   Error("AliHBTAnalysis::RunCoherencyCheck()",
+         "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
+         fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
+   return kTRUE;
+  }
+ cout<<"Checking number of Particles AND Particles Types in each event ...";fflush(stdout);
+ AliHBTEvent *partEvent;
+ AliHBTEvent *trackEvent;
+ for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
+  {
+    partEvent= fReader->GetParticleEvent(i); //gets the "ith" event 
+    trackEvent = fReader->GetTrackEvent(i);
+    
+    if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
+    if ( (partEvent == 0x0) || (partEvent == 0x0) )
+     {
+       Error("AliHBTAnalysis::RunCoherencyCheck()",
+             "One event is NULL and the other one not. Event Number %d",i);
+       return kTRUE;    
+     }
+    
+    if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
+     {
+       Error("AliHBTAnalysis::RunCoherencyCheck()",
+             "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
+             i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
+       return kTRUE;
+     }
+    else
+     for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
+      {
+        if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
+         {
+           Error("AliHBTAnalysis::RunCoherencyCheck()",
+                 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
+                 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
+           return kTRUE;
+           
+         }
+      }
+  }
+ cout<<"  Done"<<endl;
+ cout<<"  Everything looks OK"<<endl;
+ return kFALSE;
+}
+
+
+/*************************************************************************************/  
+
+/*************************************************************************************/  
+
+
diff --git a/HBTAN/AliHBTAnalysis.h b/HBTAN/AliHBTAnalysis.h
new file mode 100644 (file)
index 0000000..67dbf4a
--- /dev/null
@@ -0,0 +1,83 @@
+#ifndef ALIHBTANALYSIS_H
+#define ALIHBTANALYSIS_H
+
+#include <TObject.h>
+
+
+
+class AliHBTParticleCut;
+class AliHBTCut;
+class AliHBTPairCut;
+class AliHBTPair;
+
+class AliHBTRun;
+class AliHBTReader;
+class AliHBTTwoPartFctn;      
+class AliHBTFourPartFctn;
+
+
+class TList;
+
+class AliHBTAnalysis: public TObject
+ {
+   public:
+     AliHBTAnalysis();
+
+     virtual ~AliHBTAnalysis();
+
+     virtual void Process(Option_t* option = "TracksAndParticles");
+     
+
+     void SetGlobalPairCut(AliHBTPairCut* cut);
+     
+     void AddTrackFunction(AliHBTTwoPartFctn*);
+     void AddParticleFunction(AliHBTTwoPartFctn*);
+     void AddParticleAndTrackFunction(AliHBTFourPartFctn*);
+     
+     void AddResolutionFunction(AliHBTFourPartFctn* f){AddParticleAndTrackFunction(f);}
+     
+     void SetReader(AliHBTReader* r){fReader = r;}
+     
+     void Write();
+   protected:
+     
+     Bool_t RunCoherencyCheck();
+     
+     
+     AliHBTReader* fReader;
+     
+     virtual void ProcessTracks();
+     virtual void ProcessParticles();
+     virtual void ProcessTracksAndParticles();
+     
+     
+     AliHBTTwoPartFctn**  fTrackFunctions; //array of pointers to functions that analyze rekonstructed tracks
+     AliHBTTwoPartFctn**  fParticleFunctions; //array of pointers to functions that analyze generated particles
+     AliHBTFourPartFctn** fParticleAndTrackFunctions; //array of pointers to functions that analyze both 
+                                        //reconstructed tracks and generated particles
+               //i.e. - resolution analyzers
+     UInt_t fNTrackFunctions;
+     UInt_t fNParticleFunctions;
+     UInt_t fNParticleAndTrackFunctions;
+               
+     /**********************************************/
+     /* Control parameters  */
+
+      AliHBTPairCut *fPairCut;
+      
+   // AliHBTCut *fParticleCut; 
+     /**********************************************/
+     
+     
+   private:
+     static const Int_t fgkHbtAnalyzeAll;
+     static const UInt_t fgkFctnArraySize;
+/*********************************************/   
+   public:
+     ClassDef(AliHBTAnalysis,0)
+ };
+
+
+
+
+#endif
diff --git a/HBTAN/AliHBTCorrelFctn.cxx b/HBTAN/AliHBTCorrelFctn.cxx
new file mode 100644 (file)
index 0000000..aa76170
--- /dev/null
@@ -0,0 +1,30 @@
+#include "AliHBTCorrelFctn.h"
+
+
+
+ClassImp(AliHBTQInvCorrelFctn)
+
+//Corroleation function is created from dividing two histograms of QInvariant:
+//  of particles from the same evnt
+//by 
+//  of particles from different events
+
+TH1* AliHBTQInvCorrelFctn::GetResult()
+{
+ return GetRatio(GetDenominator()->GetMaximum()/GetNumerator()->GetMaximum());
+}
+
+
+ClassImp(AliHBTInvMassCorrelFctn)
+
+AliHBTInvMassCorrelFctn::
+AliHBTInvMassCorrelFctn(Int_t nbins, Double_t maxXval, Double_t minXval):
+                        AliHBTTwoPartFctn1D(nbins,maxXval,minXval)
+{
+  Rename("InvMass CF","Invariant Mass Correlation Function");
+}
+
+TH1* AliHBTInvMassCorrelFctn::GetResult()
+{
+ return GetRatio(GetDenominator()->GetMaximum()/GetNumerator()->GetMaximum());
+}
diff --git a/HBTAN/AliHBTCorrelFctn.h b/HBTAN/AliHBTCorrelFctn.h
new file mode 100644 (file)
index 0000000..b34919c
--- /dev/null
@@ -0,0 +1,47 @@
+#ifndef ALIHBTCORRELFUNCTION_H
+#define ALIHBTCORRELFUNCTION_H
+
+#include "AliHBTFunction.h"
+#include "AliHBTParticle.h"
+//Set of functions:
+//   Q Invaraint Correlation Function
+//   Invariant Mass Function
+//
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+//Piotr.Skowronski@cern.ch
+
+
+class AliHBTQInvCorrelFctn: public AliHBTTwoPartFctn1D
+{
+//Q Invaraint Correlation Function
+//1D two particle function 
+ public:
+   AliHBTQInvCorrelFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0):
+                        AliHBTTwoPartFctn1D(nbins,maxXval,minXval){}
+   virtual ~AliHBTQInvCorrelFctn(){};
+   TH1* GetResult();
+ protected:
+   Double_t GetValue(AliHBTPair * pair){return pair->GetQInv();}
+  public:
+    ClassDef(AliHBTQInvCorrelFctn,1)
+};
+
+
+class AliHBTInvMassCorrelFctn: public AliHBTTwoPartFctn1D
+{
+//   Invariant Mass Function 
+ public:
+   AliHBTInvMassCorrelFctn(Int_t nbins = 2000, Double_t maxXval = 2., Double_t minXval = 0.0);
+   virtual ~AliHBTInvMassCorrelFctn(){};
+   TH1* GetResult();
+ protected:
+   Double_t GetValue(AliHBTPair * pair) { return pair->GetInvMass();}
+  public:
+    ClassDef(AliHBTInvMassCorrelFctn,1)
+};
+
+
+
+#endif
diff --git a/HBTAN/AliHBTEvent.cxx b/HBTAN/AliHBTEvent.cxx
new file mode 100644 (file)
index 0000000..215b24f
--- /dev/null
@@ -0,0 +1,99 @@
+#include "AliHBTEvent.h"
+#include "AliHBTParticle.h"
+
+ClassImp(AliHBTEvent)
+
+const UInt_t AliHBTEvent::fgkInitEventSize = 10000;
+
+
+
+/**************************************************************************/ 
+AliHBTEvent::AliHBTEvent()
+ {
+    if(fgkInitEventSize<1) 
+     {
+      Fatal("AliHBTEvent::AliHBTEvent()",
+            "fgkInitEventSize has a stiupid value (%d). Change it to positive number and recompile",
+             fgkInitEventSize);
+      
+     }
+    fSize=fgkInitEventSize;
+    fParticles = new AliHBTParticle* [fSize];
+    fNParticles = 0;
+ }
+/**************************************************************************/ 
+
+AliHBTEvent::~AliHBTEvent()
+ {
+   this->Reset();//delete all particles
+   if(fParticles)
+    { 
+     delete [] fParticles; //and delete array itself
+    }
+   fParticles = 0x0;
+ }
+/**************************************************************************/ 
+void  AliHBTEvent::Reset()
+{
+  //deletes all particles from the event
+  if(fParticles)
+    {
+      for(Int_t i =0; i<fNParticles; i++)
+        delete fParticles[i];
+    }
+   fNParticles = 0;
+} 
+
+AliHBTParticle* AliHBTEvent::GetParticleSafely(Int_t n)
+ {
+   if( (n<0) || (fNParticles<=n) ) return 0x0;
+   else return fParticles[n];
+   
+ }
+/**************************************************************************/ 
+
+void  AliHBTEvent:: AddParticle(AliHBTParticle* hbtpart)
+ {
+   //Adds new perticle to the event
+   if ( fNParticles+1 >= fSize) Expand(); //if there is no space in array, expand it
+   fParticles[fNParticles++] = hbtpart; //add a pointer
+ }
+
+/**************************************************************************/ 
+void  AliHBTEvent::AddParticle(TParticle* part)
+ {
+   AddParticle( new AliHBTParticle(*part) );
+ }
+/**************************************************************************/ 
+void  AliHBTEvent::
+AddParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+            Double_t vx, Double_t vy, Double_t vz, Double_t time)
+ {
+   AddParticle(new  AliHBTParticle(pdg,px,py,pz,etot,vx,vy,vz,time) );
+ }
+/**************************************************************************/ 
+
+void AliHBTEvent::Expand()
+ {
+ //expands the array with pointers to particles
+ //about the size defined in fgkInitEventSize
+  fSize+=fgkInitEventSize;
+  AliHBTParticle** tmpParticles = new AliHBTParticle* [fSize]; //create new array of pointers
+  //check if we got memory / if not Abort
+  if (!tmpParticles) Fatal("AliHBTEvent::Expand()","No more space in memory");
+  
+  
+  for(Int_t i = 0; i<fNParticles ;i++)
+   //copy pointers to the new array
+   {
+     tmpParticles[i] = fParticles[i];
+   }
+  delete [] fParticles; //delete old array
+   fParticles = tmpParticles; //copy new pointer to the array of pointers to particles
+  
+ }
diff --git a/HBTAN/AliHBTEvent.h b/HBTAN/AliHBTEvent.h
new file mode 100644 (file)
index 0000000..d83e588
--- /dev/null
@@ -0,0 +1,56 @@
+#ifndef ALIHBTEvent_H
+#define ALIHBTEvent_H
+//This class sters HBT perticles for one event
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+#include <TObject.h>
+
+class AliHBTParticle;
+class TParticle;
+class AliHBTEvent: public TObject
+ {
+  public:
+    AliHBTEvent();
+    virtual ~AliHBTEvent();
+    const static UInt_t fgkInitEventSize; //initial number of the array
+                                         //if expanded, this size is used also
+    AliHBTParticle* GetParticle(Int_t n);  //gets particle 
+    AliHBTParticle* GetParticleSafely(Int_t n); //gets particle with index check
+    
+    void    AddParticle(AliHBTParticle*); //adds particle to the event
+    void    AddParticle(TParticle*); //adds particle to the event
+    void    AddParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+                        Double_t vx, Double_t vy, Double_t vz, Double_t time);
+    
+    Int_t   GetNumberOfParticles() const;
+    void   Reset(); //deletes all entries
+  protected:
+    AliHBTParticle ** fParticles; //!array of pointers to the particles
+    Int_t fNParticles; //!number of particles in Event
+    Int_t fSize;       //!current size of the array
+    
+    void  Expand();    //expands the array if necessary
+  private:
+    
+  public:
+    ClassDef(AliHBTEvent,1)
+ };
+/**************************************************************************/ 
+inline 
+AliHBTParticle* AliHBTEvent::GetParticle(Int_t n)
+ {
+ //Gets particle without boundary check
+   return fParticles[n];
+ }
+
+/**************************************************************************/ 
+inline 
+Int_t  AliHBTEvent::GetNumberOfParticles() const
+ {
+ //reurns number of particles in this event
+   return fNParticles;
+ }
+#endif
diff --git a/HBTAN/AliHBTFunction.cxx b/HBTAN/AliHBTFunction.cxx
new file mode 100644 (file)
index 0000000..a80843f
--- /dev/null
@@ -0,0 +1,335 @@
+#include "AliHBTFunction.h"
+/******************************************************************/
+/*
+Piotr Krzysztof Skowronski
+Piotr.Skowronski@cern.ch
+Base classes for HBT functions
+
+             function
+              /    \
+             /      \
+            /        \
+           /          \
+          /            \
+         /              \
+        /                \
+    two part          four part 
+    /   |   \         /   |   \
+   /    |    \       /    |    \
+  1D   2D    3D     1D   2D    3D
+
+ four particle functions are intendent to be resolution functions:
+ it is mecessary to have simulated particle pair corresponding to given
+ recontructed track pair in order to calculate function simualted value 
+ and recontructed value to be further histogrammed
+*/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTFunction )
+
+AliHBTFunction::AliHBTFunction()
+{
+ fPairCut = new AliHBTEmptyPairCut(); //dummy cut
+}
+
+void AliHBTFunction::
+Write()
+ {
+   if (GetNumerator()) GetNumerator()->Write();
+   if (GetDenominator()) GetDenominator()->Write();
+   TH1* res = GetResult();
+   if (res) GetResult()->Write();
+ }
+/******************************************************************/
+
+TH1* AliHBTFunction::
+GetRatio(Double_t normfactor)
+ {
+   TString str = fName + " ratio";
+   TH1 *result = (TH1*)GetNumerator()->Clone(str.Data());
+   
+   result->SetTitle(str.Data());
+   result->Sumw2();
+   
+   result->Divide(GetNumerator(),GetDenominator(),normfactor);
+   
+   return result;
+   
+ }
+/******************************************************************/
+void AliHBTFunction::SetPairCut(AliHBTPairCut* cut)
+{
+//Sets new Pair Cut. Old one is deleted
+//Note that it is created new object instead of simple pointer set
+//I do not want to have pointer 
+//to object created somewhere else
+//because in that case I could not believe that 
+//it would always exist (sb could delete it)
+//so we have always own copy
+
+ if(!cut) 
+   {
+     Error("AliHBTFunction::SetPairCut","argument is NULL");
+     return;
+   }
+ delete fPairCut;
+ fPairCut = (AliHBTPairCut*)cut->Clone();
+}
+
+/******************************************************************/
+
+void AliHBTFunction::
+Rename(const Char_t * name)
+ {
+ //renames the function and histograms
+  SetName(name);
+  SetTitle(name);
+  
+  TString numstr = fName + " Numerator";  //title and name of the 
+                                           //numerator histogram
+  TString denstr = fName + " Denominator";//title and name of the 
+                                           //denominator histogram
+  GetNumerator()->SetName(numstr.Data());
+  GetNumerator()->SetTitle(numstr.Data());
+  
+  GetDenominator()->SetName(denstr.Data());
+  GetDenominator()->SetTitle(denstr.Data());
+  
+ }
+
+void AliHBTFunction::
+Rename(const Char_t * name, const Char_t * title)
+ {
+ //renames and retitle the function and histograms
+  SetName(name);
+  SetTitle(title);
+  
+  TString numstrn = fName + " Numerator";  //name of the 
+                                           //numerator histogram
+
+  TString numstrt = fTitle + " Numerator";  //title of the 
+                                           //numerator histogram
+                  
+  TString denstrn = fName + " Denominator";//name of the 
+                                           //denominator histogram
+
+  TString denstrt = fTitle + " Denominator";//title of the 
+                                           //denominator histogram
+                  
+  GetNumerator()->SetName(numstrn.Data());
+  GetNumerator()->SetTitle(numstrt.Data());
+  
+  GetDenominator()->SetName(denstrn.Data());
+  GetDenominator()->SetTitle(denstrt.Data());
+
+
+ }
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTTwoPartFctn )
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTFourPartFctn)
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTTwoPartFctn1D )
+
+AliHBTTwoPartFctn1D::
+AliHBTTwoPartFctn1D(Int_t nbins, Double_t maxXval, Double_t minXval)
+ {
+ //Constructor of Two Part One Dimentional Function 
+ // nbins: number of bins in histograms - default 100
+ // maxXval and minXval: range of histgram(s) default 0 - 0.15 (GeV)
+   TString numstr = fName + " Numerator";  //title and name of the 
+                                           //numerator histogram
+   TString denstr = fName + " Denominator";//title and name of the 
+                                           //denominator histogram
+         
+   fNumerator   = new TH1D(numstr.Data(),numstr.Data(),nbins,minXval,maxXval);
+   fDenominator = new TH1D(denstr.Data(),denstr.Data(),nbins,minXval,maxXval);
+   
+   fNumerator->Sumw2();
+   fDenominator->Sumw2();
+   
+ }
+/******************************************************************/
+AliHBTTwoPartFctn1D::~AliHBTTwoPartFctn1D()
+{
+  delete fNumerator;
+  delete fDenominator;
+}
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTTwoPartFctn2D )
+
+AliHBTTwoPartFctn2D::
+AliHBTTwoPartFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval , 
+                    Int_t nYbins, Double_t maxYval, Double_t minYval)
+
+{
+   TString numstr = fName + " Numerator";  //title and name of the 
+                                           //numerator histogram
+   TString denstr = fName + " Denominator";//title and name of the 
+                                           //denominator histogram
+         
+   fNumerator   = new TH2D(numstr.Data(),numstr.Data(),
+                           nXbins,minXval,maxXval,
+              nYbins,minYval,maxYval);
+              
+   fDenominator = new TH2D(denstr.Data(),denstr.Data(),
+                           nXbins,minXval,maxXval,
+              nYbins,minYval,maxYval);
+   
+   fNumerator->Sumw2();
+   fDenominator->Sumw2();
+
+}        
+AliHBTTwoPartFctn2D::~AliHBTTwoPartFctn2D()
+{
+  delete fNumerator;
+  delete fDenominator;
+}
+void AliHBTTwoPartFctn2D::ProcessSameEventParticles(AliHBTPair* pair)
+{
+  pair = CheckPair(pair);
+  if(pair) 
+   { 
+     Double_t x,y;
+     GetValues(pair,x,y);
+     fNumerator->Fill(y,x);
+   }
+}
+
+void AliHBTTwoPartFctn2D::ProcessDiffEventParticles(AliHBTPair* pair)
+{
+  pair = CheckPair(pair);
+  if(pair) 
+   { 
+     Double_t x,y;
+     GetValues(pair,x,y);
+     fDenominator->Fill(y,x);
+   }
+
+}
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTTwoPartFctn3D)
+
+AliHBTTwoPartFctn3D::
+AliHBTTwoPartFctn3D(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                    Int_t nYbins, Double_t maxYval, Double_t minYval, 
+                    Int_t nZbins, Double_t maxZval, Double_t minZval)
+
+{
+   TString numstr = fName + " Numerator";  //title and name of the 
+                                           //numerator histogram
+   TString denstr = fName + " Denominator";//title and name of the 
+                                           //denominator histogram
+         
+   fNumerator   = new TH3D(numstr.Data(),numstr.Data(),
+                           nXbins,minXval,maxXval,
+              nYbins,minYval,maxYval,
+              nZbins,minZval,maxZval);
+              
+   fDenominator = new TH3D(denstr.Data(),denstr.Data(),
+                           nXbins,minXval,maxXval,
+              nYbins,minYval,maxYval,
+              nZbins,minZval,maxZval);
+   
+   fNumerator->Sumw2();
+   fDenominator->Sumw2();
+
+}        
+
+
+AliHBTTwoPartFctn3D::~AliHBTTwoPartFctn3D()
+{
+  delete fNumerator;
+  delete fDenominator;
+}
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTFourPartFctn2D)
+
+
+AliHBTFourPartFctn2D::
+AliHBTFourPartFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval , 
+                    Int_t nYbins, Double_t maxYval, Double_t minYval)
+
+{
+   TString numstr = fName + " Numerator";  //title and name of the 
+                                           //numerator histogram
+   TString denstr = fName + " Denominator";//title and name of the 
+                                           //denominator histogram
+         
+   fNumerator   = new TH2D(numstr.Data(),numstr.Data(),
+                           nXbins,minXval,maxXval,
+              nYbins,minYval,maxYval);
+              
+   fDenominator = new TH2D(denstr.Data(),denstr.Data(),
+                           nXbins,minXval,maxXval,
+              nYbins,minYval,maxYval);
+   
+   fNumerator->Sumw2();
+   fDenominator->Sumw2();
+
+}        
+AliHBTFourPartFctn2D::~AliHBTFourPartFctn2D()
+{
+  delete fNumerator;
+  delete fDenominator;
+}
+void AliHBTFourPartFctn2D::
+ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
+{
+  partpair  = CheckPair(partpair);
+  trackpair = CheckPair(trackpair);
+  if( partpair && trackpair) 
+   { 
+     Double_t x,y;
+     GetValues(trackpair,partpair,x,y);
+     fNumerator->Fill(y,x);
+   }
+}
+
+void AliHBTFourPartFctn2D::
+ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
+{
+  partpair  = CheckPair(partpair);
+  trackpair = CheckPair(trackpair);
+  if( partpair && trackpair)
+   { 
+     Double_t x,y;
+     GetValues(trackpair,partpair,x,y);
+     fDenominator->Fill(y,x);
+   }
+
+}
+
diff --git a/HBTAN/AliHBTFunction.h b/HBTAN/AliHBTFunction.h
new file mode 100644 (file)
index 0000000..d45571b
--- /dev/null
@@ -0,0 +1,234 @@
+//Piotr Skowronski@cern.ch
+
+#ifndef ALIHBTFUNCTION_H
+#define ALIHBTFUNCTION_H
+
+#include "AliHBTParticleCut.h"
+#include "AliHBTPairCut.h"
+#include "AliHBTPair.h"
+
+#include <TH2.h>
+#include <TH3.h>
+
+
+class AliHBTAnalysis;
+
+class AliHBTFunction: public TNamed
+//Abstract base class for HBT functions
+{
+  public:
+    AliHBTFunction();
+    virtual ~AliHBTFunction(){}
+    
+    virtual TH1* GetNumerator() =0;
+    virtual TH1* GetDenominator() =0;
+    virtual TH1* GetResult() = 0;
+    
+    virtual void Write();
+    
+    TH1* GetRatio(Double_t normfactor = 1.0);
+    void Rename(const Char_t * name); //renames the function and histograms ==title is the same that name
+    void Rename(const Char_t * name, const Char_t * title); //renames and retitle the function and histograms
+    
+    void SetPairCut(AliHBTPairCut*);
+    
+    virtual AliHBTPair* CheckPair(AliHBTPair* pair);
+    
+  protected:
+    
+    AliHBTPairCut*      fPairCut;
+    
+  public:  
+   ClassDef(AliHBTFunction,1)
+};
+inline AliHBTPair* AliHBTFunction::CheckPair(AliHBTPair* pair)
+{
+  //check if pair and both particles meets the cut criteria
+  if(fPairCut->Pass(pair)) //if the pair is BAD
+   {//it is BAD 
+    pair = pair->GetSwapedPair();
+    if(fPairCut->Pass(pair)) //so try reverse combination
+     { 
+       return 0x0;//it is BAD as well - so return
+     }
+   }
+  return pair; 
+}
+
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+class AliHBTTwoPartFctn: public AliHBTFunction
+{
+  public:
+    AliHBTTwoPartFctn(){}
+    virtual ~AliHBTTwoPartFctn(){}
+    
+    virtual void ProcessSameEventParticles(AliHBTPair* pair) = 0;
+    virtual void ProcessDiffEventParticles(AliHBTPair* pair) = 0;
+    
+    
+  protected:
+  public:  
+   ClassDef(AliHBTTwoPartFctn,1)
+  
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+class AliHBTFourPartFctn: public AliHBTFunction
+{
+  public:
+    AliHBTFourPartFctn(){};
+    virtual ~AliHBTFourPartFctn(){};
+    
+    virtual void 
+    ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) = 0;
+    virtual void 
+    ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) = 0;
+            
+  protected:
+  public:  
+   ClassDef(AliHBTFourPartFctn,1)
+  
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+class AliHBTTwoPartFctn1D: public AliHBTTwoPartFctn
+{
+ public:
+  AliHBTTwoPartFctn1D(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0);
+  virtual ~AliHBTTwoPartFctn1D();
+  
+  
+  TH1* GetNumerator(){return fNumerator;}
+  TH1* GetDenominator(){return fDenominator;}
+
+  void ProcessSameEventParticles(AliHBTPair* pair);
+  void ProcessDiffEventParticles(AliHBTPair* pair);
+  
+ protected:
+  //retruns velue to be histogrammed
+  virtual Double_t GetValue(AliHBTPair* pair) = 0; 
+  
+  TH1D* fNumerator;
+  TH1D* fDenominator;
+  
+ public:
+  ClassDef(AliHBTTwoPartFctn1D,1)
+};
+
+inline void 
+AliHBTTwoPartFctn1D::ProcessSameEventParticles(AliHBTPair* pair)
+{
+ //Fills the numerator
+   pair = CheckPair(pair);
+   if(pair) fNumerator->Fill(GetValue(pair));
+}
+  
+  
+inline void
+AliHBTTwoPartFctn1D::ProcessDiffEventParticles(AliHBTPair* pair)
+ {
+  //fills denumerator
+   pair = CheckPair(pair);
+   if(pair) fDenominator->Fill(GetValue(pair));
+
+  }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTTwoPartFctn2D: public AliHBTTwoPartFctn
+{
+ public:
+  AliHBTTwoPartFctn2D(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                      Int_t nYbins = 200, Double_t maxYval = .15, Double_t minYval =-0.15);
+  ~AliHBTTwoPartFctn2D();
+  
+  TH1* GetNumerator(){return fNumerator;}
+  TH1* GetDenominator(){return fDenominator;}
+  
+  void ProcessSameEventParticles(AliHBTPair* pair);
+  void ProcessDiffEventParticles(AliHBTPair* pair);
+  virtual void GetValues(AliHBTPair* pair, Double_t&, Double_t&) = 0;
+  
+ protected:
+  TH2D* fNumerator;
+  TH2D* fDenominator;
+  
+ public:
+  ClassDef(AliHBTTwoPartFctn2D,1)
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTTwoPartFctn3D: public AliHBTTwoPartFctn
+{
+ public:
+  AliHBTTwoPartFctn3D(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                      Int_t nYbins = 200, Double_t maxYval = .15, Double_t minYval =-0.15, 
+                      Int_t nZbins = 200, Double_t maxZval = .15, Double_t minZval =-0.15);
+           
+  virtual ~AliHBTTwoPartFctn3D();
+
+  TH1* GetNumerator(){return fNumerator;}
+  TH1* GetDenominator(){return fDenominator;}
+
+ protected:
+  TH3D* fNumerator;
+  TH3D* fDenominator;
+ public:
+  ClassDef(AliHBTTwoPartFctn3D,1)
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+class AliHBTFourPartFctn2D: public AliHBTFourPartFctn
+{
+ public:
+  AliHBTFourPartFctn2D(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                       Int_t nYbins = 200, Double_t maxYval = .15, Double_t minYval =-0.15);
+  ~AliHBTFourPartFctn2D();
+  
+  TH1* GetNumerator(){return fNumerator;}
+  TH1* GetDenominator(){return fDenominator;}
+  
+  void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair);
+  void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair);
+  virtual void GetValues(AliHBTPair*,AliHBTPair*, Double_t&, Double_t&) = 0;
+  
+ protected:
+  TH2D* fNumerator;
+  TH2D* fDenominator;
+  
+ public:
+  ClassDef(AliHBTFourPartFctn2D,1)
+};
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+
+#endif
diff --git a/HBTAN/AliHBTPair.cxx b/HBTAN/AliHBTPair.cxx
new file mode 100644 (file)
index 0000000..7eacc80
--- /dev/null
@@ -0,0 +1,128 @@
+#include "AliHBTPair.h"
+#include "AliHBTParticle.h"
+
+ClassImp(AliHBTPair)
+
+//value of rev 
+/************************************************************************/
+AliHBTPair::AliHBTPair(Bool_t rev)
+ {
+//value of rev defines if it is Swaped
+//if you pass kTRUE swpaped pair will NOT be created
+//though you wont be able to get the swaped pair from this pair
+
+  if(!rev) fSwapedPair = new AliHBTPair(kTRUE); //if false create swaped pair
+  else fSwapedPair = 0x0; //if true set the pointer to NULL
+  
+ }
+/************************************************************************/
+Double_t AliHBTPair::GetInvMass()
+{
+  if(fInvMassNotCalc)
+   {
+     CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method 
+     
+     if(fInvMassSqr<0)  fInvMass = TMath::Sqrt(-fInvMassSqr);
+     else fInvMass = TMath::Sqrt(fInvMassSqr); 
+     
+     fInvMassNotCalc = kFALSE;
+   }
+  return fInvMass;
+}
+/************************************************************************/
+Double_t AliHBTPair::GetQSideCMSLC()
+{
+  //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
+  if (fQSideCMSLCNotCalc)
+   {
+    fQSideCMSLC = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
+    fQSideCMSLCNotCalc = kFALSE;
+   }
+  return fQSideCMSLC;
+}
+/************************************************************************/
+Double_t AliHBTPair::GetQOutCMSLC()
+{
+ if(fQOutCMSLCNotCalc)
+  {
+   CalculateSums();
+   CalculateDiffs();
+   Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
+   fQOutCMSLC = 0.5*k2/GetKt();
+   fQOutCMSLCNotCalc = kFALSE;
+  }
+ return fQOutCMSLC;
+}
+/************************************************************************/
+Double_t AliHBTPair::GetQLongCMSLC()
+{
+ if (fQLongCMSLCNotCalc)
+  {
+    CalculateSums();
+    CalculateDiffs();
+    Double_t beta = fPzSum/fESum;
+    Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
+    fQLongCMSLC = gamma*( fPzDiff - beta*fEDiff );
+    fQLongCMSLCNotCalc = kFALSE;
+  }
+ return fQLongCMSLC; 
+}
+/************************************************************************/
+Double_t AliHBTPair::GetKt()
+{
+  if(fKtNotCalc)
+   { 
+     CalculateSums();
+     fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
+     fKtNotCalc = kFALSE;
+   }
+  return fKt;
+}
+/************************************************************************/
+
+Double_t AliHBTPair::GetKStar()
+{
+  if (fKStarNotCalc)
+   { 
+    
+    CalculateSums();
+
+    Double_t Ptrans = fPxSum*fPxSum + fPySum*fPySum;
+    Double_t Mtrans = fESum*fESum - fPzSum*fPzSum;
+    Double_t Pinv =   TMath::Sqrt(Mtrans - Ptrans);
+
+    Double_t Q = ( fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/Pinv;
+    
+    CalculateQInvL();
+    
+    Q = TMath::Sqrt( Q*Q - fQInvL);
+    fKStar = Q/2;
+    fKStarNotCalc = kFALSE;
+   }
+  return fKStar;
+}
+/************************************************************************/
+
+Double_t AliHBTPair::GetQInv()
+{
+  if(fQInvNotCalc)
+   {
+    CalculateQInvL();
+
+    if (fQInvL < 0.) fQInv =  TMath::Sqrt(-fQInvL);
+    else             fQInv = -TMath::Sqrt( fQInvL);
+
+    fQInvNotCalc = kFALSE;
+   }
+   
+  return fQInv;
+}
+
+
+
+
+
+
+
+
diff --git a/HBTAN/AliHBTPair.h b/HBTAN/AliHBTPair.h
new file mode 100644 (file)
index 0000000..1d87ec9
--- /dev/null
@@ -0,0 +1,199 @@
+#ifndef ALIHBTPAIR_H
+#define ALIHBTPAIR_H
+
+#include <TObject.h>
+
+
+#include "AliHBTParticle.h"
+//class implements pair of particles and taking care of caluclation (almost)
+//all of pair properties (Qinv, InvMass,...)
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+class AliHBTPair: public TObject
+{
+ public:
+   AliHBTPair(Bool_t rev = kFALSE); //contructor
+   ~AliHBTPair(){}
+   void SetParticles(AliHBTParticle*,AliHBTParticle*); //sets particles in the pair
+   AliHBTPair* GetSwapedPair() {return fSwapedPair;} //returns pair with swapped particles
+   
+   AliHBTParticle* Particle1() const {return fPart1;} //returns pointer to first particle
+   AliHBTParticle* Particle2() const {return fPart2;} //returns pointer to decond particle
+   
+   //Center Mass System - Longitudinally Comoving
+   
+   Double_t GetInvMass(); //returns invariant mass of the pair
+   
+   Double_t GetQSideCMSLC(); //returns Q Side CMS longitudionally co-moving
+   Double_t GetQOutCMSLC(); //returns Q out CMS longitudionally co-moving
+   Double_t GetQLongCMSLC(); //returns Q Long CMS longitudionally co-moving
+   
+   
+   Double_t GetKt();  //returns K transverse
+   Double_t GetKStar();
+   
+   Double_t GetQInv(); //returns Q invariant
+   Double_t GetQSide(); //returns Q side
+   Double_t GetQLong(); //returns Q long 
+   Double_t GetQOut(); //returns Q out
+   
+   
+ protected:
+   AliHBTParticle* fPart1;  //pointer to first particle
+   AliHBTParticle* fPart2;  //pointer to second particle
+  
+   AliHBTPair* fSwapedPair; //pointer to swapped pair
+   
+/************************************************************/
+/************CMS (LC) Q's   *********************************/
+/************************************************************/
+   //Center Mass System - Longitudinally Comoving
+   
+   Double_t fQSideCMSLC;  //value of Q side CMS longitudially co-moving
+   Bool_t   fQSideCMSLCNotCalc; //flag indicating if fQSideCMSLC is already calculated for this pair
+   
+   Double_t fQOutCMSLC; //value of Q out CMS longitudially co-moving
+   Bool_t   fQOutCMSLCNotCalc;//flag indicating if fQOutCMSLC is already calculated for this pair
+   
+   Double_t fQLongCMSLC; //value of Q long CMS longitudially co-moving
+   Bool_t   fQLongCMSLCNotCalc;//flag indicating if fQLongCMSLC is already calculated for this pair
+/************************************************************/
+/************************************************************/
+   Double_t fQInv;  //half of differnece of 4-momenta
+   Bool_t   fQInvNotCalc;//flag indicating if fQInv is already calculated for this pair
+   
+   Double_t fInvMass;  //invariant mass
+   Bool_t   fInvMassNotCalc;//flag indicating if fInvMass is already calculated for this pair
+   
+   Double_t fKt; //K == sum vector of particle's momenta. Kt transverse component
+   Bool_t   fKtNotCalc;//flag indicating if fKt is already calculated for this pair
+   
+   Double_t fKStar;
+   Bool_t   fKStarNotCalc;
+   
+   Double_t fPInv;  //invariant momentum
+   Double_t fQSide; //Q Side
+   Double_t fOut;//Q Out
+   Double_t fQLong;//Q Long
+
+   
+   Double_t fInvMassSqr;//squre of invariant mass
+   Bool_t   fMassSqrNotCalc; //
+   void     CalculateInvMassSqr();
+   
+   Double_t fQInvL;
+   Bool_t   fQInvLNotCalc;
+   void     CalculateQInvL();
+   
+   Double_t fPxSum;
+   Double_t fPySum;
+   Double_t fPzSum;
+   Double_t fESum;
+   Bool_t   fSumsNotCalc;
+   void     CalculateSums();
+   
+   Double_t fPxDiff;
+   Double_t fPyDiff;
+   Double_t fPzDiff;
+   Double_t fEDiff;
+   Bool_t   fDiffsNotCalc;
+   void     CalculateDiffs();
+   
+   
+   /***************************************************/
+   void CalculateBase();
+   Bool_t fChanged;
+   
+   
+ private:
+ public:
+  ClassDef(AliHBTPair,1)
+};
+/****************************************************************/
+inline
+void AliHBTPair::SetParticles(AliHBTParticle* p1,AliHBTParticle* p2)
+{
+ //sets the particle to the pair
+ fPart1 = p1;
+ fPart2 = p2;
+ if (fSwapedPair) //if we have Swaped (so we are not)
+   fSwapedPair->SetParticles(p2,p1); //set particles for him too
+ // Resel all calculations (flags)
+
+ fChanged           = kTRUE;
+ fSumsNotCalc       = kTRUE;
+ fDiffsNotCalc      = kTRUE;
+ fMassSqrNotCalc    = kTRUE;
+ fInvMassNotCalc    = kTRUE;
+ fQInvNotCalc       = kTRUE;
+ fQSideCMSLCNotCalc = kTRUE;
+ fQOutCMSLCNotCalc  = kTRUE;
+ fQLongCMSLCNotCalc = kTRUE;
+ fKtNotCalc         = kTRUE;
+ fKStarNotCalc      = kTRUE;
+ fQInvLNotCalc      = kTRUE;
+ //and do nothing until will be asked for
+}
+/****************************************************************/
+inline
+Double_t AliHBTPair::GetQSide()
+{
+ return 0.0;
+}
+/****************************************************************/
+inline 
+void AliHBTPair::CalculateInvMassSqr()
+ {
+  if (fMassSqrNotCalc)
+   {
+     CalculateSums();
+     Double_t fPart12s= (fPxSum*fPxSum) + (fPySum*fPySum) + (fPzSum*fPzSum);
+     fInvMassSqr=fESum*fESum-fPart12s;
+
+     fMassSqrNotCalc = kFALSE;
+   }
+ }
+/****************************************************************/
+inline 
+void AliHBTPair::CalculateQInvL()
+ {
+  if (fQInvLNotCalc)
+  {
+   CalculateDiffs();
+   fQInvL = fEDiff*fEDiff - ( fPxDiff*fPxDiff + fPyDiff*fPyDiff + fPzDiff*fPzDiff );
+   fQInvLNotCalc = kFALSE;
+  }
+ }
+/****************************************************************/ 
+inline 
+void AliHBTPair::CalculateSums()
+ {
+   if(fSumsNotCalc)
+    {
+     fPxSum = fPart1->Px()+fPart2->Px();
+     fPySum = fPart1->Py()+fPart2->Py();
+     fPzSum = fPart1->Pz()+fPart2->Pz();
+     fESum  = fPart1->Energy() + fPart2->Energy();
+     fSumsNotCalc = kFALSE;
+    }
+ }
+/****************************************************************/
+inline 
+void AliHBTPair::CalculateDiffs()
+ {
+   if(fDiffsNotCalc)
+    {
+     fPxDiff = fPart1->Px()-fPart2->Px();
+     fPyDiff = fPart1->Py()-fPart2->Py();
+     fPzDiff = fPart1->Pz()-fPart2->Pz();
+     fEDiff  = fPart1->Energy() - fPart2->Energy();
+     fDiffsNotCalc = kFALSE;
+    }
+ }
+
+#endif
diff --git a/HBTAN/AliHBTPairCut.cxx b/HBTAN/AliHBTPairCut.cxx
new file mode 100644 (file)
index 0000000..75723c4
--- /dev/null
@@ -0,0 +1,197 @@
+#include "AliHBTPairCut.h"
+#include "AliHBTPair.h"
+#include <iostream.h>
+
+
+ClassImp(AliHBTPairCut)
+const Int_t AliHBTPairCut::fkgMaxCuts = 50;
+/**********************************************************/
+
+AliHBTPairCut::AliHBTPairCut()
+{
+//constructor
+  fFirstPartCut = new AliHBTEmptyParticleCut(); //empty cuts
+  fSecondPartCut= new AliHBTEmptyParticleCut(); //empty cuts
+    
+  fCuts = new AliHbtBasePairCut*[fkgMaxCuts];
+  fNCuts = 0;
+}
+/**********************************************************/
+
+AliHBTPairCut::AliHBTPairCut(const AliHBTPairCut& in)
+{
+ //copy constructor
+ fCuts = new AliHbtBasePairCut*[fkgMaxCuts];
+ fNCuts = in.fNCuts;
+
+ fFirstPartCut = (AliHBTParticleCut*)in.fFirstPartCut->Clone();
+ fSecondPartCut = (AliHBTParticleCut*)in.fSecondPartCut->Clone();
+ for (Int_t i = 0;i<fNCuts;i++)
+   {
+     fCuts[i] = (AliHbtBasePairCut*)in.fCuts[i]->Clone();//create new object (clone) and rember pointer to it
+   }
+}
+/**********************************************************/
+
+AliHBTPairCut::~AliHBTPairCut()
+{
+//destructor
+  for (Int_t i = 0;i<fNCuts;i++)
+   {
+     delete fCuts[i];
+   }
+  delete []fCuts;
+} 
+/**********************************************************/
+
+/**********************************************************/
+
+void AliHBTPairCut::AddBasePairCut(AliHbtBasePairCut* basecut)
+ {
+ //adds the base pair cut (cut on one value)
+   if (!basecut) return;
+   if( fNCuts == (fkgMaxCuts-1) )
+    {
+      Warning("AddBasePairCut","Not enough place for another cut");
+      return;
+    }
+   fCuts[fNCuts++]=basecut;
+ }
+/**********************************************************/
+
+Bool_t AliHBTPairCut::Pass(AliHBTPair* pair)
+{
+//methods which checks if given pair meets all criteria of the cut
+//if it meets returns FALSE
+//if NOT   returns    TRUE
+ if(!pair) 
+  {
+    Warning("Pass()","No Pasaran! We never accept NULL pointers");
+    return kTRUE;
+  }
+ //check particle's cuts
+ if( (   fFirstPartCut->Pass( pair->Particle1() )   ) || 
+     (   fSecondPartCut->Pass(pair->Particle2() )   )   )
+   {  
+     return kTRUE;
+   }
+ //examine all base pair cuts
+ for (Int_t i = 0;i<fNCuts;i++)
+   {
+    if ( (fCuts[i]->Pass(pair)) ) return kTRUE; //if one of the cuts reject, then reject
+   }
+   
+// cout<<"passed "<<pair->Particle1()->GetPdgCode()<<"  "<<pair->Particle2()->GetPdgCode()<<endl;
+ return kFALSE;
+}
+/**********************************************************/
+
+void AliHBTPairCut::SetFirstPartCut(AliHBTParticleCut* cut)
+{
+ if(!cut) 
+   {
+     Error("SetFirstPartCut","argument is NULL");
+     return;
+   }
+ delete fFirstPartCut;
+ fFirstPartCut = (AliHBTParticleCut*)cut->Clone();
+
+}
+/**********************************************************/
+
+void AliHBTPairCut::SetSecondPartCut(AliHBTParticleCut* cut)
+{
+ if(!cut) 
+   {
+     Error("SetSecondPartCut","argument is NULL");
+     return;
+   }
+ delete fSecondPartCut;
+ fSecondPartCut = (AliHBTParticleCut*)cut->Clone();
+}
+/**********************************************************/
+
+void AliHBTPairCut::SetPartCut(AliHBTParticleCut* cut)
+{
+//sets the the same cut on both particles
+ if(!cut) 
+   {
+     Error("SetFirstPartCut","argument is NULL");
+     return;
+   }
+ delete fFirstPartCut;
+ fFirstPartCut = (AliHBTParticleCut*)cut->Clone();
+ delete fSecondPartCut; //even if null should not be harmful
+ fSecondPartCut = fFirstPartCut;
+}
+/**********************************************************/
+
+void AliHBTPairCut::SetQInvRange(Double_t min, Double_t max)
+{
+  AliHBTQInvCut* cut= (AliHBTQInvCut*)FindCut(kHbtPairCutPropQInv);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTQInvCut(min,max);
+
+}
+
+AliHbtBasePairCut* AliHBTPairCut::FindCut(AliHBTPairCutProperty property)
+{
+ for (Int_t i = 0;i<fNCuts;i++)
+  {
+    if (fCuts[i]->GetProperty() == property) 
+       return fCuts[i]; //we found the cut we were searching for
+  }
+ return 0x0; //we did not found this cut
+  
+}
+/**********************************************************/
+
+void AliHBTPairCut::Streamer(TBuffer &b)
+{
+     // Stream all objects in the array to or from the I/O buffer.
+
+   UInt_t R__s, R__c;
+   if (b.IsReading()) 
+    {
+      Version_t v = b.ReadVersion(&R__s, &R__c);
+      TObject::Streamer(b);
+      b >> fFirstPartCut;
+      b >> fSecondPartCut;
+      b >> fNCuts;
+      for (Int_t i = 0;i<fNCuts;i++)
+       {
+        b >> fCuts[i];
+       }
+      b.CheckByteCount(R__s, R__c,AliHBTPairCut::IsA());
+    } 
+   else 
+    {
+     R__c = b.WriteVersion(AliHBTPairCut::IsA(), kTRUE);
+     TObject::Streamer(b);
+     b << fFirstPartCut;
+     b << fSecondPartCut;
+     b << fNCuts;
+     for (Int_t i = 0;i<fNCuts;i++)
+      {
+       b << fCuts[i];
+      }
+     b.SetByteCount(R__c, kTRUE);
+   }
+}
+
+ClassImp(AliHBTEmptyPairCut)
+
+void AliHBTEmptyPairCut::Streamer(TBuffer &b)
+ {
+   AliHBTPairCut::Streamer(b);
+ }
+
+ClassImp(AliHbtBasePairCut)
+
+ClassImp(AliHBTQInvCut)
diff --git a/HBTAN/AliHBTPairCut.h b/HBTAN/AliHBTPairCut.h
new file mode 100644 (file)
index 0000000..2a8e42a
--- /dev/null
@@ -0,0 +1,139 @@
+//Piotr Skowronski@cern.ch
+//Class implemnts cut on the pair of particles
+//
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+#ifndef ALIHBTPAIRCUT_H
+#define ALIHBTPAIRCUT_H
+
+
+#include "AliHBTParticleCut.h"
+#include "AliHBTPair.h"
+
+class AliHbtBasePairCut;
+
+enum AliHBTPairCutProperty
+ {
+  kHbtPairCutPropQInv, //Q invariant
+  kHbtPairCutPropNone
+ };
+
+class AliHBTPairCut: public TObject
+{
+  public:
+    AliHBTPairCut();
+    AliHBTPairCut(const AliHBTPairCut&);
+    
+    virtual ~AliHBTPairCut();
+    virtual Bool_t Pass(AliHBTPair*);
+    
+    void SetFirstPartCut(AliHBTParticleCut*);  //sets the cut on the first particle
+    void SetSecondPartCut(AliHBTParticleCut*); //sets the cut on the first particle
+    
+    void SetPartCut(AliHBTParticleCut*);//sets the the same cut on both particles
+     
+    void AddBasePairCut(AliHbtBasePairCut*);
+    
+    void SetQInvRange(Double_t min, Double_t max);
+    
+  protected:
+    AliHBTParticleCut*      fFirstPartCut;
+    AliHBTParticleCut*      fSecondPartCut;
+
+    AliHbtBasePairCut** fCuts; //!
+    Int_t fNCuts;
+       
+       
+    AliHbtBasePairCut* FindCut(AliHBTPairCutProperty);
+  private:
+    static const Int_t fkgMaxCuts;
+  public:
+    ClassDef(AliHBTPairCut,1)
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTEmptyPairCut:  public AliHBTPairCut
+{
+//Empty - it passes possitively all particles - it means returns always False
+//Class describing cut on pairs of particles
+  public:
+    AliHBTEmptyPairCut(){};
+    AliHBTEmptyPairCut(const AliHBTEmptyPairCut&){}; 
+    virtual ~AliHBTEmptyPairCut(){};
+    
+    Bool_t Pass(AliHBTPair*) {return kFALSE;} //accpept everything
+
+    ClassDef(AliHBTEmptyPairCut,1)
+};
+
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHbtBasePairCut: public TObject
+ {
+   //This class defines the range of some property - pure virtual
+   //Property is coded by AliHBTCutTypes type
+   
+   public:
+     
+     AliHbtBasePairCut(Double_t min = 0.0, Double_t max = 0.0, AliHBTPairCutProperty prop= kHbtPairCutPropNone):
+                      fMin(min),fMax(max),fProperty(prop){}
+
+     virtual   ~AliHbtBasePairCut(){}
+     
+     Bool_t    Pass(AliHBTPair*);
+     
+     void      SetRange(Double_t min, Double_t max){fMin = min; fMax = max;}
+
+     void      SetMinimum(Double_t min){fMin = min;}
+     void      SetMaximum(Double_t max){fMax = max;}
+     
+     Double_t  GetMinimum() const {return fMin;}
+     Double_t  GetMaximum() const {return fMax;}
+     
+     AliHBTPairCutProperty GetProperty() const {return fProperty;}
+     
+   protected:
+     virtual Double_t  GetValue(AliHBTPair*) = 0;
+
+     Double_t fMin;
+     Double_t fMax;
+     
+     AliHBTPairCutProperty fProperty;
+     
+   private:
+   public:
+     ClassDef(AliHbtBasePairCut,1)
+   
+ };
+
+inline Bool_t
+AliHbtBasePairCut::Pass(AliHBTPair* pair) 
+ {
+   Double_t value = GetValue(pair);
+   if ( (value > fMin) && (value <fMax ) ) return kFALSE; //accepted
+   else return kTRUE; //rejected
+ }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+class AliHBTQInvCut: public AliHbtBasePairCut
+ {
+   public:
+    AliHBTQInvCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBasePairCut(min,max,kHbtPairCutPropQInv){}
+    virtual ~AliHBTQInvCut(){}
+   protected:
+    virtual Double_t  GetValue(AliHBTPair* pair){return pair->GetQInv();}
+   private:
+   public:
+     ClassDef(AliHBTQInvCut,1)
+ };
+
+#endif
diff --git a/HBTAN/AliHBTParticle.cxx b/HBTAN/AliHBTParticle.cxx
new file mode 100644 (file)
index 0000000..a1d8d84
--- /dev/null
@@ -0,0 +1,55 @@
+//Simplified TParticle class
+
+
+#include "AliHBTParticle.h"
+
+#include <TParticle.h>
+
+ClassImp(AliHBTParticle)
+
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle():  
+                fPdgCode(0), fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0)
+{//empty particle
+}
+
+
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+                     Double_t vx, Double_t vy, Double_t vz, Double_t time):
+  fPdgCode(pdg), fPx(px), fPy(py),fPz(pz),fE(etot), 
+                 fVx(vx), fVy(vy),fVz(vz),fVt(time)
+{
+//mormal constructor
+  
+  if (GetPDG()) {
+     fCalcMass    = GetPDG()->Mass();
+  } else {
+     Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
+     if (a2 >= 0) fCalcMass =  TMath::Sqrt(a2);
+     else         fCalcMass = -TMath::Sqrt(-a2);
+  }
+}
+
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle(const TParticle &p):
+   fPdgCode(p.GetPdgCode()), 
+   fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()), 
+   fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T())
+{
+ //all copied in the initialization
+}
+
+//______________________________________________________________________________
+const Char_t* AliHBTParticle::GetName() const 
+{
+   static char def[4] = "XXX";
+   const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPdgCode);
+   if (ap) return ap->GetName();
+   else    return def;
+}
+
+
+//______________________________________________________________________________
+
+
diff --git a/HBTAN/AliHBTParticle.h b/HBTAN/AliHBTParticle.h
new file mode 100644 (file)
index 0000000..c1dc4e7
--- /dev/null
@@ -0,0 +1,104 @@
+#ifndef ALIHBTPARTICLE
+#define ALIHBTPARTICLE
+
+//Ali HBT Particle: simplified class TParticle
+//Simplified in order to minimize the size of object
+//  - we want to keep a lot of such a objects in memory
+
+#include <TObject.h>
+#include <TLorentzVector.h>
+#include <TMath.h>
+#include <TDatabasePDG.h>
+
+class TParticle;
+
+class AliHBTParticle : public TObject
+{
+
+public:
+                                // ****** constructors and destructor
+  AliHBTParticle();
+
+  AliHBTParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+                 Double_t vx, Double_t vy, Double_t vz, Double_t time);
+
+  AliHBTParticle(const TParticle& p);
+
+  virtual ~AliHBTParticle(){};
+
+
+  Int_t          GetPdgCode      () const { return fPdgCode; }
+  Double_t       GetCalcMass     () const { return fCalcMass; }
+  Double_t       GetMass         ()       { return GetPDG()->Mass();}
+
+  TParticlePDG*  GetPDG          (){return TDatabasePDG::Instance()->GetParticle(fPdgCode);}
+
+  Int_t          Beauty          ()  { return GetPDG()->Beauty(); }
+  Int_t          Charm           ()  { return GetPDG()->Charm(); }
+  Int_t          Strangeness     ()  { return GetPDG()->Strangeness();}
+  void ProductionVertex(TLorentzVector &v) { v.SetXYZT(fVx,fVy,fVz,fVt);}
+
+
+  Double_t         Vx    () const { return fVx;}
+  Double_t         Vy    () const { return fVy;}
+  Double_t         Vz    () const { return fVz;}
+  Double_t         T     () const { return fVt;}
+
+  Double_t         Px    () const { return fPx; } //X coordinate of the momentum
+  Double_t         Py    () const { return fPy; } //Y coordinate of the momentum
+  Double_t         Pz    () const { return fPz; } //Z coordinate of the momentum
+  Double_t         P     () const                 //momentum
+    { return TMath::Sqrt(fPx*fPx+fPy*fPy+fPz*fPz); }
+  
+  void Momentum(TLorentzVector &v) { v.SetPxPyPzE(fPx,fPy,fPz,fE);}
+    
+  Double_t         Pt    () const  //transverse momentum
+    { return TMath::Sqrt(fPx*fPx+fPy*fPy); }
+  Double_t         Energy() const { return fE; }
+  
+                                   //Pseudo Rapidity
+  Double_t         Eta   () const { if (P() != fPz) return 0.5*TMath::Log((P()+fPz)/(P()-fPz)); 
+                                    else return 1.e30;}
+
+                                   //Rapidity
+  Double_t         Y     () const { if (fE  != fPz) return 0.5*TMath::Log((fE+fPz)/(fE-fPz));
+                                    else return 1.e30;}
+
+  Double_t         Phi   () const { return kPI+TMath::ATan2(-fPy,-fPx); }
+
+  Double_t         Theta () const { return (fPz==0)?kPI/2:TMath::ACos(fPz/P()); }
+
+                                // setters
+
+
+  void           SetMomentum(Double_t px, Double_t py, Double_t pz, Double_t e)
+                             {fPx=px; fPy=py; fPz=pz; fE=e;}
+  void           SetMomentum(const TLorentzVector& p)
+                             {SetMomentum(p.Px(),p.Py(),p.Pz(),p.Energy());}
+
+  void           SetProductionVertex(Double_t vx, Double_t vy, Double_t vz, Double_t t)
+                             {fVx=vx; fVy=vy; fVz=vz; fVt=t;}
+  void           SetProductionVertex(const TLorentzVector& v)
+                             {SetProductionVertex(v.X(),v.Y(),v.Z(),v.T());}
+  const Char_t*  GetName() const;
+  
+protected:
+
+  Int_t          fPdgCode;              // PDG code of the particle
+  Double_t       fCalcMass;             // Calculated mass
+
+  Double_t       fPx;                   // x component of momentum
+  Double_t       fPy;                   // y component of momentum
+  Double_t       fPz;                   // z component of momentum
+  Double_t       fE;                    // Energy
+
+  Double_t       fVx;                   // x of production vertex
+  Double_t       fVy;                   // y of production vertex
+  Double_t       fVz;                   // z of production vertex
+  Double_t       fVt;                   // t of production vertex
+
+public:
+  ClassDef(AliHBTParticle,1)  // TParticle vertex particle information
+};
+
+#endif
diff --git a/HBTAN/AliHBTParticleCut.cxx b/HBTAN/AliHBTParticleCut.cxx
new file mode 100644 (file)
index 0000000..45a4fc9
--- /dev/null
@@ -0,0 +1,315 @@
+#include "AliHBTParticleCut.h"
+
+#include <iostream.h>
+
+ClassImp(AliHBTParticleCut)
+const Int_t AliHBTParticleCut::fkgMaxCuts = 50;
+/******************************************************************/
+
+AliHBTParticleCut::AliHBTParticleCut()
+ {
+   fCuts = new AliHbtBaseCut* [fkgMaxCuts];//last property in the property
+                                         //property enum => defines number of properties
+   fNCuts = 0;
+   fPID  = 0;
+ }
+/******************************************************************/
+
+AliHBTParticleCut::AliHBTParticleCut(const AliHBTParticleCut& in)
+{
+  fCuts = new AliHbtBaseCut* [fkgMaxCuts];//last property in the property
+                                         //property enum => defines number of properties
+  fNCuts = in.fNCuts;
+  fPID  = in.fPID;
+  for (Int_t i = 0;i<fNCuts;i++)
+   {
+     fCuts[i] = (AliHbtBaseCut*)in.fCuts[i]->Clone();//create new object (clone) and rember pointer to it
+   }
+}
+/******************************************************************/
+
+AliHBTParticleCut::~AliHBTParticleCut()
+{
+  for (Int_t i = 0;i<fNCuts;i++)
+   {
+     delete fCuts[i];
+   }
+  delete []fCuts;
+} 
+/******************************************************************/
+
+Bool_t AliHBTParticleCut::Pass(AliHBTParticle* p)
+{
+//method checks all the cuts that are set (in the list)
+//If any of the baseCuts rejects particle False(rejection) is returned
+ if(!p) 
+  {
+    Warning("Pass()","No Pasaran! We never accept NULL pointers");
+    return kTRUE;
+  }
+ if( (p->GetPdgCode() != fPID) && ( fPID != 0)) return kTRUE;
+ for (Int_t i = 0;i<fNCuts;i++)
+   {
+    if ( (fCuts[i]->Pass(p)) ) return kTRUE; //if one of the cuts rejects, then reject
+   }
+  return kFALSE;
+}
+/******************************************************************/
+
+void AliHBTParticleCut::AddBasePartCut(AliHbtBaseCut* basecut)
+{
+  //adds the base pair cut (cut on one value)
+   if (!basecut) return;
+   if( fNCuts == (fkgMaxCuts-1) )
+    {
+      Warning("AddBasePartCut","Not enough place for another cut");
+      return;
+    }
+   fCuts[fNCuts++]=basecut;
+}
+
+/******************************************************************/
+AliHbtBaseCut* AliHBTParticleCut::FindCut(AliHBTCutProperty property)
+{
+ for (Int_t i = 0;i<fNCuts;i++)
+  {
+    if (fCuts[i]->GetProperty() == property) 
+       return fCuts[i]; //we found the cut we were searching for
+  }
+ return 0x0; //we did not found this cut
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetMomentumRange(Double_t min, Double_t max)
+{
+  AliHBTMomentumCut* cut= (AliHBTMomentumCut*)FindCut(kHbtP);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTMomentumCut(min,max);
+}
+/******************************************************************/
+
+
+void AliHBTParticleCut::SetPtRange(Double_t min, Double_t max)
+{
+  AliHBTPtCut* cut= (AliHBTPtCut*)FindCut(kHbtPt);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTPtCut(min,max);
+
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetEnergyRange(Double_t min, Double_t max)
+{
+  AliHBTEnergyCut* cut= (AliHBTEnergyCut*)FindCut(kHbtE);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTEnergyCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetRapidityRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtRapidity);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTRapidityCut(min,max);
+
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPseudoRapidityRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtPseudoRapidity);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTPseudoRapidityCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPxRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtPx);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTPxCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPyRange(Double_t min, Double_t max)
+{  
+  AliHbtBaseCut* cut = FindCut(kHbtPy);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTPyCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPzRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtPz);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTPzCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPhiRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtPhi);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTPhiCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetThetaRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtTheta);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTThetaCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetVxRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtVx);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTVxCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetVyRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtVy);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTVyCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetVzRange(Double_t min, Double_t max)
+{
+  AliHbtBaseCut* cut = FindCut(kHbtVz);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTVzCut(min,max);
+}
+
+/******************************************************************/
+void AliHBTParticleCut::Streamer(TBuffer &b)
+{
+     // Stream all objects in the array to or from the I/O buffer.
+
+   UInt_t R__s, R__c;
+   if (b.IsReading()) 
+    {
+      Version_t v = b.ReadVersion(&R__s, &R__c);
+      TObject::Streamer(b);
+      b >> fPID;
+      b >> fNCuts;
+      for (Int_t i = 0;i<fNCuts;i++)
+       {
+        b >> fCuts[i];
+       }
+       b.CheckByteCount(R__s, R__c,AliHBTParticleCut::IsA());
+    } 
+   else 
+    {
+     R__c = b.WriteVersion(AliHBTParticleCut::IsA(), kTRUE);
+     TObject::Streamer(b);
+     b << fPID;
+     b << fNCuts;
+     for (Int_t i = 0;i<fNCuts;i++)
+      {
+       b << fCuts[i];
+      }
+     b.SetByteCount(R__c, kTRUE);
+   }
+}
+
+void AliHBTParticleCut::Print(void)
+{
+  cout<<"Printing AliHBTParticleCut, this = "<<this<<endl;
+  cout<<"fPID  "<<fPID<<endl;
+  cout<<"fNCuts  "<<fNCuts <<endl;
+  for (Int_t i = 0;i<fNCuts;i++)
+      {
+       cout<<"  fCuts["<<i<<"]  "<<fCuts[i]<<endl<<"   ";
+       fCuts[i]->Print();
+      }
+}
+
+/******************************************************************/
+/******************************************************************/
+
+ClassImp(AliHBTEmptyParticleCut)
+void AliHBTEmptyParticleCut::Streamer(TBuffer &b)
+ {
+  AliHBTParticleCut::Streamer(b);
+ }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp(AliHbtBaseCut)
+void AliHbtBaseCut::Print(void)
+{
+  cout<<"fMin="<<fMin <<", fMax=" <<fMax<<"    ";
+  PrintProperty();
+}
+void AliHbtBaseCut::PrintProperty(void)
+{
+ switch (fProperty)
+  {
+   case  kHbtP: 
+     cout<<"kHbtP"; break;
+   case  kHbtPt: 
+     cout<<"kHbtPt"; break;
+   case  kHbtE: 
+     cout<<"kHbtE"; break;
+   case  kHbtRapidity: 
+     cout<<"kHbtRapidity"; break;
+   case  kHbtPseudoRapidity: 
+     cout<<"kHbtPseudoRapidity"; break;
+   case  kHbtPx: 
+     cout<<"kHbtPx"; break;
+   case  kHbtPy: 
+     cout<<"kHbtPy"; break;
+   case  kHbtPz: 
+     cout<<"kHbtPz"; break;   
+   case  kHbtPhi: 
+     cout<<"kHbtPhi"; break;
+   case  kHbtTheta: 
+     cout<<"kHbtTheta"; break;
+   case  kHbtVx: 
+     cout<<"kHbtVx"; break;
+   case  kHbtVy: 
+     cout<<"kHbtVy"; break;
+   case  kHbtVz: 
+     cout<<"kHbtVz"; break;
+   case  kHbtNone: 
+     cout<<"kHbtNone"; break;
+   default: 
+     cout<<"Property Not Found";
+  }
+ cout<<endl;
+}
+ClassImp( AliHBTMomentumCut )
+
+ClassImp( AliHBTPtCut )
+ClassImp( AliHBTEnergyCut )
+ClassImp( AliHBTRapidityCut )
+ClassImp( AliHBTPseudoRapidityCut )
+ClassImp( AliHBTPxCut )
+ClassImp( AliHBTPyCut )
+ClassImp( AliHBTPzCut )
+ClassImp( AliHBTPhiCut )
+ClassImp( AliHBTThetaCut )
+ClassImp( AliHBTVxCut )
+ClassImp( AliHBTVyCut )
+ClassImp( AliHBTVzCut )
+
+
diff --git a/HBTAN/AliHBTParticleCut.h b/HBTAN/AliHBTParticleCut.h
new file mode 100644 (file)
index 0000000..89020d4
--- /dev/null
@@ -0,0 +1,340 @@
+//Piotr Skowronski@cern.ch
+//Classes for single particle cuts
+//User should use only AliHBTParticleCut, eventually EmptyCut which passes all particles
+//There is all interface for setting cuts on all particle properties
+//The main method is Pass - which returns 
+//        True in order to reject particle
+//        False in case it meets all the criteria of the given cut
+//
+//User should create (and also destroy) cuts himself 
+//and then pass them to the Analysis or Function by a proper method
+//
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+
+#ifndef ALIHBTPARTICLECUT_H
+#define ALIHBTPARTICLECUT_H
+
+#include <TObject.h>
+#include "AliHBTParticle.h"
+
+
+class AliHBTEmptyParticleCut;
+class AliHBTParticleCut;
+class AliHBTPairCut;
+class AliHBTPair;
+class AliHbtBaseCut;
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+enum AliHBTCutProperty
+ {
+//codes particle property
+  kHbtP,  //Momentum
+  kHbtPt, //Transverse momentum
+  kHbtE,  //Energy
+  kHbtRapidity, //
+  kHbtPseudoRapidity,
+  kHbtPx, //X coordinate of the momentum
+  kHbtPy, //Y coordinate of the momentum
+  kHbtPz, //Z coordinate of the momentum
+  kHbtPhi,//angle
+  kHbtTheta,//angle
+  kHbtVx,  // vertex X coordinate
+  kHbtVy,  // vertex Y coordinate
+  kHbtVz,  // vertex Z coordinate
+//_____________________________
+  kHbtNone 
+ };
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTParticleCut: public TObject
+{
+//Class describing cut on pairs of particles
+  public:
+    AliHBTParticleCut();
+    AliHBTParticleCut(const AliHBTParticleCut&);
+    virtual ~AliHBTParticleCut();
+
+    virtual Bool_t Pass(AliHBTParticle*);
+    
+    void AddBasePartCut(AliHbtBaseCut*);
+    
+    Int_t GetPID() const { return fPID;}
+    void SetPID(Int_t pid){fPID=pid;}
+    void SetMomentumRange(Double_t min, Double_t max);
+    void SetPRange(Double_t min, Double_t max){SetMomentumRange(min,max);}
+    void SetPtRange(Double_t min, Double_t max);
+    void SetEnergyRange(Double_t min, Double_t max);
+    void SetRapidityRange(Double_t min, Double_t max);
+    void SetYRange(Double_t min, Double_t max){SetRapidityRange(min,max);}
+    void SetPseudoRapidityRange(Double_t min, Double_t max);
+    void SetPxRange(Double_t min, Double_t max);
+    void SetPyRange(Double_t min, Double_t max);
+    void SetPzRange(Double_t min, Double_t max);
+    void SetPhiRange(Double_t min, Double_t max);
+    void SetThetaRange(Double_t min, Double_t max);
+    void SetVxRange(Double_t min, Double_t max);
+    void SetVyRange(Double_t min, Double_t max);
+    void SetVzRange(Double_t min, Double_t max);
+    
+    void Print(void);
+  protected:
+     
+     AliHbtBaseCut* FindCut(AliHBTCutProperty);
+     
+     AliHbtBaseCut ** fCuts;//! Array with cuts
+     Int_t fNCuts;
+
+     Int_t fPID; //particle PID  - if=0 (rootino) all pids are accepted
+          
+  private:
+     static const Int_t fkgMaxCuts;
+  public:
+    ClassDef(AliHBTParticleCut,1)
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTEmptyParticleCut:  public AliHBTParticleCut
+{
+//Empty - it passes possitively all particles - it means returns always False
+//Class describing cut on pairs of particles
+  public:
+    AliHBTEmptyParticleCut(){};
+    virtual ~AliHBTEmptyParticleCut(){};
+    
+    Bool_t Pass(AliHBTParticle*) {return kFALSE;} //accpept everything
+
+    ClassDef(AliHBTEmptyParticleCut,1)
+};
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHbtBaseCut: public TObject
+ {
+   //This class defines the range of some property - pure virtual
+   //Property is coded by AliHBTCutTypes type
+   
+   public:
+     
+     AliHbtBaseCut(Double_t min = 0.0, Double_t max = 0.0,AliHBTCutProperty prop = kHbtNone):
+                   fProperty(prop),fMin(min),fMax(max){}
+
+     virtual   ~AliHbtBaseCut(){}
+     
+     Bool_t    Pass(AliHBTParticle *);
+     
+     void      SetRange(Double_t min, Double_t max){fMin = min; fMax = max;}
+
+     void      SetMinimum(Double_t min){fMin = min;}
+     void      SetMaximum(Double_t max){fMax = max;}
+     
+     Double_t  GetMinimum() const {return fMin;}
+     Double_t  GetMaximum() const {return fMax;}
+     
+     AliHBTCutProperty GetProperty() const {return fProperty;}
+     void Print(void);     
+   protected:
+     virtual Double_t  GetValue(AliHBTParticle *) = 0;
+
+     AliHBTCutProperty fProperty;
+     Double_t fMin;
+     Double_t fMax;
+   private:
+     void PrintProperty(void);
+   public:
+     ClassDef(AliHbtBaseCut,1)
+   
+ };
+
+inline Bool_t
+AliHbtBaseCut::Pass(AliHBTParticle *p)
+ {
+   if ( (GetValue(p) > fMin) && (GetValue(p) <fMax ) ) return kFALSE; //accepted
+   else return kTRUE; //rejected
+ }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTMomentumCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTMomentumCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtP){}
+    virtual ~AliHBTMomentumCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->P();}
+  public:  
+    ClassDef(AliHBTMomentumCut,1)
+  
+ };
+
+
+class AliHBTPtCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTPtCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPt){}
+    virtual ~AliHBTPtCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Pt();}
+  public: 
+    ClassDef(AliHBTPtCut,1)
+  
+ };
+
+
+class AliHBTEnergyCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTEnergyCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtE){}
+    virtual ~AliHBTEnergyCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Energy();}
+  public: 
+    ClassDef(AliHBTEnergyCut,1)
+  
+ };
+
+class AliHBTRapidityCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTRapidityCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtRapidity){}
+    virtual ~AliHBTRapidityCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Y();}
+  public:   
+    ClassDef(AliHBTRapidityCut,1)
+  
+ };
+
+class AliHBTPseudoRapidityCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTPseudoRapidityCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPseudoRapidity){}
+    virtual ~AliHBTPseudoRapidityCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Eta();}
+  public:   
+    ClassDef(AliHBTPseudoRapidityCut,1)
+  
+ };
+
+class AliHBTPxCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTPxCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPx){}
+    virtual ~AliHBTPxCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Px();}
+  public:   
+    ClassDef(AliHBTPxCut,1)
+  
+ };
+
+class AliHBTPyCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTPyCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPy){}
+    virtual ~AliHBTPyCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Py();}
+  public:   
+    ClassDef(AliHBTPyCut,1)
+  
+ };
+
+
+class AliHBTPzCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTPzCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPz){}
+    virtual ~AliHBTPzCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Pz();}
+  public:   
+    ClassDef(AliHBTPzCut,1)
+  
+ };
+
+class AliHBTPhiCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTPhiCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPhi){}
+    virtual ~AliHBTPhiCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Phi();}
+  public:    
+    ClassDef(AliHBTPhiCut,1)
+  
+ };
+
+class AliHBTThetaCut: public AliHbtBaseCut
+ {
+  public: 
+    AliHBTThetaCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtTheta){}
+    virtual ~AliHBTThetaCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Theta();}
+  public:  
+    ClassDef(AliHBTThetaCut,1)
+  
+ };
+
+class AliHBTVxCut: public AliHbtBaseCut
+ {
+ //Cut of the X coordinate of the vertex position
+  public: 
+    AliHBTVxCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtVx){}
+    virtual ~AliHBTVxCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Vx();} //retruns value of the vertex
+  public:   
+    ClassDef(AliHBTVxCut,1)
+  
+ };
+
+
+class AliHBTVyCut: public AliHbtBaseCut
+ {
+ //Cut of the X coordinate of the vertex position
+  public: 
+    AliHBTVyCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtVy){}
+    virtual ~AliHBTVyCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Vy();} //retruns value of the vertex
+  public:   
+    ClassDef(AliHBTVyCut,1)
+  
+ };
+
+class AliHBTVzCut: public AliHbtBaseCut
+ {
+ //Cut of the X coordinate of the vertex position
+  public: 
+    AliHBTVzCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtVz){}
+    virtual ~AliHBTVzCut(){}
+  protected:
+    Double_t  GetValue(AliHBTParticle * p){return p->Vz();} //retruns value of the vertex
+  public:   
+    ClassDef(AliHBTVzCut,1)
+  
+ };
+
+
+
+
+
+#endif
diff --git a/HBTAN/AliHBTQResolutionFctns.cxx b/HBTAN/AliHBTQResolutionFctns.cxx
new file mode 100644 (file)
index 0000000..1e04955
--- /dev/null
@@ -0,0 +1,144 @@
+#include "AliHBTQResolutionFctns.h"
+
+AliHBTQOutResolVSQInvFctn* xxx = new AliHBTQOutResolVSQInvFctn();
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQOutResolVSQInvFctn )
+AliHBTQOutResolVSQInvFctn::
+AliHBTQOutResolVSQInvFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                          Int_t nYbins, Double_t maxYval, Double_t minYval):
+                AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QOutResolVSQInv","Q_{Out} Resolution vs. Q_{Inv}");
+}
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQSideResolVSQInvFctn )
+
+AliHBTQSideResolVSQInvFctn::
+AliHBTQSideResolVSQInvFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                          Int_t nYbins, Double_t maxYval, Double_t minYval):
+                AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QSideResolVSQInv","Q_{Side} Resolution vs. Q_{Inv}");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQLongResolVSQInvFctn )
+
+AliHBTQLongResolVSQInvFctn::
+AliHBTQLongResolVSQInvFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                           Int_t nYbins, Double_t maxYval, Double_t minYval):
+                           AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QLongResolVSQInv","Q_{Long} Resolution vs. Q_{Inv}");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQInvResolVSKtFctn )
+
+AliHBTQInvResolVSKtFctn::
+AliHBTQInvResolVSKtFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                        Int_t nYbins, Double_t maxYval, Double_t minYval):
+                           AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QInvResolVSKt","Q_{Inv} Resolution vs. K_{t}");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQOutResolVSKtFctn )
+
+AliHBTQOutResolVSKtFctn::
+AliHBTQOutResolVSKtFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                           Int_t nYbins, Double_t maxYval, Double_t minYval):
+                           AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QOutResolVSKt","Q_{Out} Resolution vs. K_{t} ");
+}
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQSideResolVSKtFctn )
+
+AliHBTQSideResolVSKtFctn::
+AliHBTQSideResolVSKtFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                           Int_t nYbins, Double_t maxYval, Double_t minYval):
+                           AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QSideResolVSKt","Q_{Side} Resolution vs. K_{t} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQLongResolVSKtFctn )
+
+AliHBTQLongResolVSKtFctn::
+AliHBTQLongResolVSKtFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                           Int_t nYbins, Double_t maxYval, Double_t minYval):
+                           AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QLongResolVSKt","Q_{Long} Resolution vs. K_{t} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQOutResolVSQOutFctn)
+
+AliHBTQOutResolVSQOutFctn::
+AliHBTQOutResolVSQOutFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                           Int_t nYbins, Double_t maxYval, Double_t minYval):
+                           AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QOutResolVSQOut","Q_{Out} Resolution vs. Q_{Out} ");
+}
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQSideResolVSQSideFctn )
+
+AliHBTQSideResolVSQSideFctn::
+AliHBTQSideResolVSQSideFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                           Int_t nYbins, Double_t maxYval, Double_t minYval):
+                           AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QSideResolVSQSide","Q_{Side} Resolution vs. Q_{Side} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQLongResolVSQLongFctn )
+
+AliHBTQLongResolVSQLongFctn::
+AliHBTQLongResolVSQLongFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, 
+                           Int_t nYbins, Double_t maxYval, Double_t minYval):
+                           AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QLongResolVSQLong","Q_{Long} Resolution vs. Q_{Long} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+
+
diff --git a/HBTAN/AliHBTQResolutionFctns.h b/HBTAN/AliHBTQResolutionFctns.h
new file mode 100644 (file)
index 0000000..2238a77
--- /dev/null
@@ -0,0 +1,240 @@
+#ifndef ALIHBTQOUTVSQINVRESOLFCTN_H
+#define ALIHBTQOUTVSQINVRESOLFCTN_H
+//General Remark:
+//CMSLC means
+//Center of Mass System Longitudially Co-moving
+
+
+//this class creates resolution function of Qout 
+//(difference of simulated pair Qout and recontructed pair)
+//in function of QInv
+//it inherits from AliHBTFourPartFctn2D
+//  it needs two pairs to compare
+//  and is two dimentional: numerator and denominator are TH2D
+
+class AliHBTQOutResolVSQInvFctn;  //QOutCMSLC  Res   VS   QInvCMSLC 
+class AliHBTQSideResolVSQInvFctn; //QSideCMSLC Res   VS   QInvCMSLC 
+class AliHBTQLongResolVSQInvFctn; //QLongCMSLC Res   VS   QInvCMSLC 
+
+class AliHBTQInvResolVSKtFctn;    //QInvCMSLC  Res   VS   Kt
+class AliHBTQOutResolVSKtFctn;    //QOutCMSLC  Res   VS   Kt
+class AliHBTQSideResolVSKtFctn;   //QSideCMSLC Res   VS   Kt
+class AliHBTQLongResolVSKtFctn;   //QLongCMSLC Res   VS   Kt
+
+class AliHBTQOutResolVSQOutFctn;  //QOutCMSLC  Res   VS   QOut
+class AliHBTQSideResolVSQSideFctn;//QSideCMSLC Res   VS   QSide
+class AliHBTQLongResolVSQLongFctn;//QLongCMSLC Res   VS   QLong
+
+
+#include "AliHBTFunction.h"
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQOutResolVSQInvFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQOutResolVSQInvFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .15, Double_t minYval =-0.15);
+   
+   virtual ~AliHBTQOutResolVSQInvFctn(){}
+   
+   TH1* GetResult(){return fNumerator;}  
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     x = partpair->GetQOutCMSLC() - trackpair->GetQOutCMSLC();
+     y = partpair->GetQInv();
+    }
+  protected:
+  private: 
+  public:
+    ClassDef(AliHBTQOutResolVSQInvFctn,1)
+ };
+
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQSideResolVSQInvFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQSideResolVSQInvFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+   virtual ~AliHBTQSideResolVSQInvFctn(){}
+
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair,  Double_t& x, Double_t& y)
+    {
+     x = partpair->GetQSideCMSLC() - trackpair->GetQSideCMSLC();
+     y = partpair->GetQInv();
+    }
+   TH1* GetResult(){return fNumerator;} 
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQSideResolVSQInvFctn,1)
+ };
+
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQLongResolVSQInvFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQLongResolVSQInvFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+   virtual ~AliHBTQLongResolVSQInvFctn(){}
+
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     x = partpair->GetQLongCMSLC() - trackpair->GetQLongCMSLC();
+     y = partpair->GetQInv();
+    }
+   TH1* GetResult(){return fNumerator;} 
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQLongResolVSQInvFctn,1)
+ };
+
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQInvResolVSKtFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQInvResolVSKtFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+   virtual ~AliHBTQInvResolVSKtFctn(){};
+
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     x = partpair->GetQInv() - trackpair->GetQInv();
+     y = partpair->GetKt();
+    }
+   TH1* GetResult(){return fNumerator;} 
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQInvResolVSKtFctn,1)
+ };
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQOutResolVSKtFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQOutResolVSKtFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .15, Double_t minYval =-0.15);
+   virtual ~AliHBTQOutResolVSKtFctn(){}
+   TH1* GetResult(){return GetNumerator();}
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     x = partpair->GetQOutCMSLC() - trackpair->GetQOutCMSLC();
+     y = partpair->GetKt();
+    }
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQOutResolVSKtFctn,1)
+ };
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQSideResolVSKtFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQSideResolVSKtFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                            Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+   virtual ~AliHBTQSideResolVSKtFctn(){}
+   TH1* GetResult(){return GetNumerator();}
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     x = partpair->GetQSideCMSLC() - trackpair->GetQSideCMSLC();
+     y = partpair->GetKt();
+    }
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQSideResolVSKtFctn,1)
+ };
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQLongResolVSKtFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQLongResolVSKtFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+   virtual ~AliHBTQLongResolVSKtFctn(){}
+
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     x = partpair->GetQLongCMSLC() - trackpair->GetQLongCMSLC();
+     y = partpair->GetKt();
+    }
+   TH1* GetResult(){return fNumerator;}  
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQLongResolVSKtFctn,1)
+ };
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQOutResolVSQOutFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQOutResolVSQOutFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .15, Double_t minYval =-0.15);
+   virtual ~AliHBTQOutResolVSQOutFctn(){}
+
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     y = partpair->GetQOutCMSLC();
+     x = y - trackpair->GetQOutCMSLC();
+    }
+   TH1* GetResult(){return fNumerator;}  
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQOutResolVSQOutFctn,1)
+ };
+
+/***********************************************************************/
+/***********************************************************************/
+
+class AliHBTQSideResolVSQSideFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQSideResolVSQSideFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .15, Double_t minYval =-0.15);
+   virtual ~AliHBTQSideResolVSQSideFctn(){}
+
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     y = partpair->GetQSideCMSLC();
+     x = y - trackpair->GetQSideCMSLC();
+    }
+   TH1* GetResult(){return fNumerator;}  
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQSideResolVSQSideFctn,1)
+ };
+
+
+/***********************************************************************/
+/***********************************************************************/
+
+class AliHBTQLongResolVSQLongFctn: public AliHBTFourPartFctn2D
+ {
+  public: 
+   AliHBTQLongResolVSQLongFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0, 
+                             Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+   virtual ~AliHBTQLongResolVSQLongFctn(){}
+
+   void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+    {
+     y = partpair->GetQLongCMSLC();
+     x = y - trackpair->GetQLongCMSLC();
+    }
+   TH1* GetResult(){return fNumerator;}  
+  protected:
+  private:
+  public:
+    ClassDef(AliHBTQLongResolVSQLongFctn,1)
+ };
+
+
+
+#endif
diff --git a/HBTAN/AliHBTReader.cxx b/HBTAN/AliHBTReader.cxx
new file mode 100644 (file)
index 0000000..979bfa0
--- /dev/null
@@ -0,0 +1,82 @@
+#include "AliHBTReader.h"
+
+#include "AliHBTParticleCut.h"
+
+
+ClassImp(AliHBTReader)
+//pure virtual
+
+/*************************************************************************************/
+
+AliHBTReader::AliHBTReader()
+{
+//constructor
+ fCuts = new TObjArray();
+}
+
+/*************************************************************************************/
+
+AliHBTReader::~AliHBTReader()
+{
+//destructor
+  fCuts->SetOwner();
+  delete fCuts;
+}
+
+/*************************************************************************************/
+
+void AliHBTReader::AddParticleCut(AliHBTParticleCut* cut)
+{
+ //sets the new cut 
+  if (!cut) //if cut is NULL return with error
+   {
+    Error("AddParticleType","NULL pointers are not accepted any more.\nIf You want to accept all particles of this type, set an empty cut ");
+    return;
+   }
+  AliHBTParticleCut *c = (AliHBTParticleCut*)cut->Clone();
+  fCuts->Add(c);
+}
+
+/*************************************************************************************/
+
+Bool_t AliHBTReader::Pass(AliHBTParticle* p)
+ {
+ //Method examines whether particle meets all cut and particle type criteria
+  
+   if(p==0x0)//of corse we not pass NULL pointers
+    {
+     Warning("Pass()","No Pasaran! We never accept NULL pointers");
+     return kTRUE;
+    }
+   //if no particle is specified, we pass all particles
+   //excluding NULL pointers, of course
+   
+  for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)   
+   {
+     AliHBTParticleCut &cut = *((AliHBTParticleCut*)fCuts->At(i));
+     if(!cut.Pass(p)) return kFALSE;  //accepted
+   }
+   
+  return kTRUE;//not accepted
+
+ }
+/*************************************************************************************/
+
+Bool_t  AliHBTReader::Pass(Int_t pid)
+{
+//this method checks if any of existing cuts accepts this pid particles
+//or any cuts accepts all particles
+
+ if(pid == 0)
+  return kTRUE;
+  
+ for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)   
+   {
+     AliHBTParticleCut &cut = *((AliHBTParticleCut*)fCuts->At(i));
+     //if some of cuts accepts all particles or some accepts particles of this type, accept
+     if ( (cut.GetPID() == 0) || (cut.GetPID() == pid) ) return kFALSE; 
+   }
+ return kTRUE;
+}
+/*************************************************************************************/
diff --git a/HBTAN/AliHBTReader.h b/HBTAN/AliHBTReader.h
new file mode 100644 (file)
index 0000000..3622244
--- /dev/null
@@ -0,0 +1,46 @@
+#ifndef ALIHBTREADER_H
+#define ALIHBTREADER_H
+
+#include <TObject.h>
+
+//Reader Base class (reads particles and tracks and
+//puts it to the AliHBTRun object
+//Piotr.Skowronski@cern.ch
+
+class AliHBTRun;
+class AliHBTEvent;
+class AliHBTParticleCut;
+class TObjArray;
+class AliHBTParticle;
+
+class AliHBTReader: public TObject
+
+{
+  public:
+    AliHBTReader();
+    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;
+    
+    
+    void AddParticleCut(AliHBTParticleCut* cut);
+    
+  protected:
+    
+    TObjArray *fCuts;//
+
+    Bool_t Pass(AliHBTParticle*);
+    Bool_t Pass(Int_t pid);
+    
+  private:
+  
+  public:
+    ClassDef(AliHBTReader,1)
+};
+
+#endif
diff --git a/HBTAN/AliHBTReaderITSv1.cxx b/HBTAN/AliHBTReaderITSv1.cxx
new file mode 100644 (file)
index 0000000..f6db7b3
--- /dev/null
@@ -0,0 +1,188 @@
+
+#include "AliHBTReader.h"
+
+ClassImp(AliHBTReaderITSv1)
+
+AliHBTReaderITSv1::AliHBTReaderITSv1(const Char_t* goodtracksfilename):
+                 fGoodITSTracksFileName(goodtracksfilename)
+
+ {
+     fParticles = new AliHBTRun();
+     fTracks    = new AliHBTRun();
+     fIsRead = kFALSE;
+ }
+AliHBTReaderITSv1::AliHBTReaderITSv1()
+{
+   delete fParticles;
+   delete fTracks;
+}
+AliHBTEvent* AliHBTReaderITSv1::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+AliHBTEvent* AliHBTReaderITSv1::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fTracks->GetEvent(n);
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderITSv1::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fParticles->GetNumberOfEvents();
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderITSv1::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+  if (!fIsRead) Read(fParticles,fTracks);
+  return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderITSv1::Read(AliHBTRun* particles, AliHBTRun *tracks)
+ {
+   AliGoodTracksITSv1 *goodITStracks = new AliGoodTracksITSv1(fGoodITSTracksFileName);
+   
+   if (!goodITStracks)
+    {
+     Error("AliHBTReaderITSv1::Read","Exiting due to problems with opening files. Errorcode %d",i);
+     return 1;
+    }
+   for (Int_t currentEvent = 0; currentEvent < goodITStracks->fNevents; currentEvent++)
+    {
+      for(Int_t i =0; i<goodITStracks->fGoodInEvent[currentEvent]; i++)
+       {
+         const struct GoodTrackITSv1 & gt = goodTPCTracks->GetTrack(currentEvent,i);
+         
+         if(Pass(gt.code)) continue;
+         
+         Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
+         Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
+         
+         AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
+         if(Pass(part)) { delete part; continue;}
+         
+         Double_t tEtot = TMath::Sqrt( gt.pxg*gt.pxg + gt.pyg*gt.pyg + gt.pzg*gt.pzg + mass*mass);
+         
+         AliHBTParticle* track = new AliHBTParticle(gt.code, gt.pxg, gt.pyg , gt.pzg, tEtot, 0., 0., 0., 0.);
+         if(Pass(track)) { delete  track;continue;}
+            
+         particles->AddParticle(currentEvent,part);
+         tracks->AddParticle(currentEvent,track);
+       }
+    }
+  delete goodITStracks;
+  fIsRead = kTRUE;
+  return 0;
+ }
+/********************************************************************/
+
+/********************************************************************/
+/********************************************************************/
+/********************************************************************/
+
+
+AliGoodTracksITSv1::~AliGoodTracksITSv1()
+{
+ delete [] fGoodInEvent;
+ for (Int_t i = 0;i<fNevents;i++)
+   delete [] fData[i];
+ delete [] fData;
+}
+/********************************************************************/
+AliGoodTracksITSv1::AliGoodTracksITSv1(const TString& infilename)
+{
+
+  cout<<"AliGoodTracksITSv1::AliGoodTracksITSv1()  ....\n";
+
+
+  fNevents = 0;
+  Int_t maxevents = 500;
+  
+  ifstream in(infilename.Data());
+
+  if(!in)
+    {
+      cerr<<"Can not open file with Good ITSv1 Tracks named:"<<infilename.Data()<<endl;
+      delete this;
+      return;
+    }
+
+  
+  fGoodInEvent = new Int_t[maxevents];
+  fData = new struct GoodTrack* [maxevents];
+
+  Int_t i;
+  Int_t lastevent = 0;
+  fGoodInEvent[0] =0;
+  fData[0] = new struct GoodTrack[50000];
+
+  Int_t evno;
+  while(in>>evno)
+   {
+    if(lastevent>evno)
+     {
+       for(i = lastevent;i<=evno;i++)
+        {
+            fGoodInEvent[i] =0;
+            fData[i] = new struct GoodTrack[50000];
+        }
+      lastevent = evno;
+     }
+    if(fGoodInEvent[evno]>=50000)
+     {
+      cerr<<"AliGoodTracksITSv1::AliGoodTracksITSv1() : Not enough place in the array\n";
+      continue;
+     }
+    in>>fData[evno][fGoodInEvent[evno]].lab;
+    in>>fData[evno][fGoodInEvent[evno]].code;
+    in>>fData[evno][fGoodInEvent[evno]].px;
+    in>>fData[evno][fGoodInEvent[evno]].py;
+    in>>fData[evno][fGoodInEvent[evno]].pz;
+    in>>fData[evno][fGoodInEvent[evno]].x;
+    in>>fData[evno][fGoodInEvent[evno]].y;
+    in>>fData[evno][fGoodInEvent[evno]].z;
+    
+ /* cout<<evno<<" ";
+  cout<<fData[evno][fGoodInEvent[evno]].lab;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
+  cout<<"\n";
+ */ 
+  fGoodInEvent[evno]++;
+ }
+ fNevents = evno+1;
+ in.close();
+ cout<<"AliGoodTracksITSv1::AliGoodTracksITSv1()  ....  Done\n";
+}
+
+
+
+const GoodTrack& AliGoodTracksITSv1::GetTrack(Int_t event, Int_t n) const
+ {
+  
+  if( (event>fNevents) || (event<0))
+   {
+     gROOT->Fatal("AliGoodTracksITSv1::GetTrack","No such Event %d",event);
+   }
+  if( (n>fGoodInEvent[event]) || (n<0))
+   {
+     gROOT->Fatal("AliGoodTracksITSv1::GetTrack","No such Good TPC Track %d",n);
+   }
+  return fData[event][n];
+
+ }
diff --git a/HBTAN/AliHBTReaderITSv1.h b/HBTAN/AliHBTReaderITSv1.h
new file mode 100644 (file)
index 0000000..2e7fbf5
--- /dev/null
@@ -0,0 +1,62 @@
+#ifndef ALIHBTREADERITSV1_H
+#define ALIHBTREADERITSV1_H
+
+#include "AliHBTReader.h"
+
+#include <TString.h>
+
+
+class AliHBTReaderITSv1: public AliHBTReader
+{
+  public:    
+    AliHBTReaderITS(const Char_t* goodtracksfilename = "itsgood_tracks");
+    virtual ~AliHBTReaderITS();
+    
+    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);//returns pointer to event with particles
+    Int_t GetNumberOfPartEvents();//returns number of particle events
+    Int_t GetNumberOfTrackEvents();//returns number of track events
+    
+  protected:
+    TString fGoodITSTracksFileName;
+    
+    AliHBTRun* fParticles; //!simulated particles
+    AliHBTRun* fTracks; //!reconstructed tracks (particles)
+    
+    Bool_t fIsRead;//flag indicating if the data are already read
+    
+  private:
+  public:
+    ClassDef(AliHBTReaderITS,1)
+};
+
+struct GoodTrackITSv1 //good tracks produced by ITSComparison V1
+{
+  Int_t fEventN; //event number
+  Int_t lab;
+  Int_t code;
+  Float_t px,py,pz;
+  Float_t x,y,z;
+  Float_t pxg,pyg,pzg,ptg;
+  Bool_t flag;
+};
+
+class AliGoodTracksITSv1
+ { 
+   //container for good tracks ITS tracking V1
+   //this class is for internal use only
+   friend class AliHBTReaderITSv1;
+   
+   private:
+     AliGoodTracksITSv1(const TString& infilename = TString("itsgood_tracks"));
+     ~AliGoodTracksITSv1();
+   
+     const GoodTrackITSv1& GetTrack(Int_t event, Int_t n) const;
+
+     Int_t  fNevents;  //Number of events
+     Int_t* fGoodInEvent; //Numbers of good track in event
+     struct GoodTrack **fData;
+ };
+
diff --git a/HBTAN/AliHBTReaderITSv2.cxx b/HBTAN/AliHBTReaderITSv2.cxx
new file mode 100644 (file)
index 0000000..fd13375
--- /dev/null
@@ -0,0 +1,5 @@
+
+#include "AliHBTReaderITSv2.h"
+
+ClassImp(AliHBTReaderITSv2)
+
diff --git a/HBTAN/AliHBTReaderITSv2.h b/HBTAN/AliHBTReaderITSv2.h
new file mode 100644 (file)
index 0000000..f24b820
--- /dev/null
@@ -0,0 +1,25 @@
+#ifndef ALIHBTREADERITSV2_H
+#define ALIHBTREADERITSV2_H
+
+#include "AliHBTReader.h"
+
+#include <TString.h>
+
+
+class AliHBTReaderITSv2: public AliHBTReaderTPC
+{
+  public:    
+    AliHBTReaderITSv2(const Char_t* trackfilename = "AliITStracksV2.root",
+                    const Char_t* clusterfilename = "AliITSclustersV2.root",
+       const Char_t* goodtracksfilename = "good_tracks_its",
+       const Char_t* galicefilename = "");
+       
+    virtual ~AliHBTReaderITSv2();
+    
+    Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+    
+  protected:    
+  private:
+  public:
+    ClassDef(AliHBTReaderITSv2,1)
+};
diff --git a/HBTAN/AliHBTReaderPPprod.cxx b/HBTAN/AliHBTReaderPPprod.cxx
new file mode 100644 (file)
index 0000000..6a69cf9
--- /dev/null
@@ -0,0 +1,392 @@
+#include "AliHBTReaderPPprod.h"
+
+#include <iostream.h>
+#include <fstream.h>
+#include <TString.h>
+#include <TTree.h>
+#include <TFile.h>
+
+#include <AliTPCtrack.h>
+#include <AliTPCParam.h>
+#include <AliTPCtracker.h>
+
+#include "AliRun.h"
+#include "AliHBTRun.h"
+#include "AliHBTEvent.h"
+#include "AliHBTParticle.h"
+#include "AliHBTParticleCut.h"
+
+
+ClassImp(AliHBTReaderPPprod)
+
+AliHBTReaderPPprod::
+ AliHBTReaderPPprod(const Char_t* trackfilename,const Char_t* clusterfilename,
+                 const Char_t* goodtracksfilename,const Char_t* galicefilename):
+                 fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
+                 fGAliceFileName(galicefilename),
+                 fGoodTPCTracksFileName(goodtracksfilename)
+{
+  //constructor, only file names are set
+  //Defaults:
+  //  trackfilename = "AliTPCtracks.root"
+  //  clusterfilename = "AliTPCclusters.root"
+  //  goodtracksfilename = "good_tracks_tpc"
+  //  galicefilename = ""  - this means: Do not open gAlice file - 
+  //                         just leave the global pointer untached
+  
+  fParticles = new AliHBTRun();
+  fTracks    = new AliHBTRun();
+
+  fTracksFile   = 0x0;  //files are opened during reading only
+  fClustersFile = 0x0;
+  
+  fIsRead = kFALSE;
+}
+/********************************************************************/
+
+AliHBTReaderPPprod::~AliHBTReaderPPprod()
+ {
+   delete fParticles;
+   delete fTracks;
+ }
+/********************************************************************/
+
+AliHBTEvent* AliHBTReaderPPprod::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+AliHBTEvent* AliHBTReaderPPprod::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fTracks->GetEvent(n);
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderPPprod::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fParticles->GetNumberOfEvents();
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderPPprod::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+  if (!fIsRead) Read(fParticles,fTracks);
+  return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+
+Int_t AliHBTReaderPPprod::Read(AliHBTRun* particles, AliHBTRun *tracks)
+ {
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+  Int_t i; //iterator and some temprary values
+  Int_t Nevents;
+  if (!particles) //check if an object is instatiated
+   {
+     Error("AliHBTReaderPPprod::Read"," particles object must instatiated before passing it to the reader");
+   }
+  if (!tracks)  //check if an object is instatiated
+   {
+     Error("AliHBTReaderPPprod::Read"," tracks object must instatiated before passing it to the reader");
+   }
+  particles->Reset();//clear runs == delete all old events
+  tracks->Reset();
+    
+  if( (i=OpenFiles()) )
+   {
+     Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
+     return i;
+   }
+  
+  AliGoodTracksPP *goodTPCTracks = new AliGoodTracksPP(fGoodTPCTracksFileName);
+  if (!goodTPCTracks)
+   {
+     Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
+     return 1;
+   }
+
+  Nevents = 100;
+
+  fClustersFile->cd();
+  AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
+  if (!TPCParam) 
+    { 
+     Error("AliHBTReaderPPprod::Read","TPC parameters have not been found !\n");
+     return 1;
+    }
+
+  TObjArray *tarray = new TObjArray(5000);
+  tarray->SetOwner();//  this causes memory leak, but in some cases deleting is infinite loop
+  
+  for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
+   {
+     cout<<"Reading Event "<<currentEvent<<endl;
+     /**************************************/
+      /**************************************/
+       /**************************************/ 
+         fTracksFile->cd();
+         Char_t  treename[100];
+         sprintf(treename,"TreeT_TPC_%d",currentEvent);
+   
+         TTree *tracktree=0;
+         
+         tracktree=(TTree*)fTracksFile->Get(treename);
+         if (!tracktree) 
+          {
+            Error("AliHBTReaderPPprod::Read","Can't get a tree with TPC tracks !\n"); 
+            return 1;
+          }
+   
+         TBranch *trackbranch=tracktree->GetBranch("tracks");
+         if (!trackbranch) 
+          {
+            Error("AliHBTReaderPPprod::Read","Can't get a branch with TPC tracks !\n"); 
+            return 2;
+          }
+         Int_t NTPCtracks=(Int_t)tracktree->GetEntries();
+         cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
+         //Copy tracks to array
+         
+         AliTPCtrack *iotrack=0;
+         
+         fClustersFile->cd();
+         AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);
+         if (!tracker) 
+          {
+            Error("AliHBTReaderPPprod::Read","Can't get a tracker !\n"); 
+            return 3;
+          }
+         tracker->LoadInnerSectors();
+         tracker->LoadOuterSectors();
+   
+         for (i=0; i<NTPCtracks; i++)
+          {
+            iotrack=new AliTPCtrack;
+            trackbranch->SetAddress(&iotrack);
+            tracktree->GetEvent(i);
+            tracker->CookLabel(iotrack,0.1);
+            tarray->AddLast(iotrack);
+          }
+        
+         
+         fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
+         fTracksFile->Delete("tracks");
+         
+         delete tracker;
+         
+         tracker = 0x0;
+         trackbranch = 0x0;
+         tracktree = 0x0;
+
+         Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent];
+   
+         Double_t xk;
+         Double_t par[5];
+         Float_t phi, lam, pt;
+         Int_t label;
+         Bool_t found;
+   
+         for (i=0; i<ngood; i++)
+          { 
+            const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i);
+            
+            if(Pass(gt.code)) continue;
+            
+            label = gt.lab;
+            found = kFALSE; //guard in case we don't find track with such a label
+            for (Int_t j=0;j<NTPCtracks;j++)
+              {
+                iotrack = (AliTPCtrack*)tarray->At(j);
+                if (iotrack->GetLabel() == label) 
+                  {
+                    found = kTRUE;
+                    break;
+                  }
+              }  
+            if(!found) 
+              {
+                Warning("Read",
+                "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
+                continue; //put comunicate on the screen and continue loop
+              }
+        
+            Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
+            Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
+            
+            AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
+            if(Pass(part)) continue;
+            
+         
+            iotrack->PropagateTo(gt.x);
+            iotrack->GetExternalParameters(xk,par);
+            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 x coordinate of momentum
+            Double_t tpz = pt * lam;
+            
+            Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);
+            
+            AliHBTParticle* track = new AliHBTParticle(gt.code, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+            if(Pass(track)) continue;
+            
+            particles->AddParticle(currentEvent,part);
+            tracks->AddParticle(currentEvent,track);
+
+          }
+         tarray->Clear();
+       /**************************************/
+      /**************************************/
+     /**************************************/  
+   }
+  
+  CloseFiles();
+  delete tarray;
+  delete goodTPCTracks;
+  fIsRead = kTRUE;
+  return 0;
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderPPprod::OpenFiles()
+{
+   fTracksFile = 0;
+   fTracksFile=TFile::Open("AliTPCtracks.root");
+   if (!fTracksFile->IsOpen()) 
+     {
+       Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCtracks.root");
+       return 1;
+     }
+   
+   fClustersFile = 0;
+   
+   fClustersFile=TFile::Open("AliTPCclusters.root");
+   if (!fClustersFile->IsOpen()) 
+    {
+      Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCclusters.root");
+      return 2;
+    }
+
+    
+
+ return 0; 
+}
+  
+
+
+/********************************************************************/
+  
+void AliHBTReaderPPprod::CloseFiles()
+{
+  fTracksFile->Close();
+  fClustersFile->Close();
+}
+
+/********************************************************************/
+
+/********************************************************************/
+/********************************************************************/
+/********************************************************************/
+
+
+AliGoodTracksPP::~AliGoodTracksPP()
+{
+ delete [] fGoodInEvent;
+ for (Int_t i = 0;i<fNevents;i++)
+   delete [] fData[i];
+ delete [] fData;
+}
+/********************************************************************/
+AliGoodTracksPP::AliGoodTracksPP(const TString& infilename)
+{
+
+  fNevents = 100;
+  ifstream in(infilename.Data());
+
+  if(!in)
+    {
+      cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
+      delete this;
+      return;
+    }
+
+  
+  fGoodInEvent = new Int_t[fNevents];
+  fData = new struct GoodTrack* [fNevents];
+
+  Int_t i;
+  for( i = 0;i<fNevents;i++)
+   {
+    fGoodInEvent[i] =0;
+    fData[i] = new struct GoodTrack[500];
+   }
+  Float_t tmp;
+  Int_t evno;
+  while(in>>evno)
+   {
+    if(fGoodInEvent[evno]>=500)
+     {
+      cerr<<"AliGoodTracksPP::AliGoodTracksPP() : Not enough place in the array\n";
+      continue;
+     }
+    in>>fData[evno][fGoodInEvent[evno]].lab;
+    in>>fData[evno][fGoodInEvent[evno]].code;
+    in>>fData[evno][fGoodInEvent[evno]].px;
+    in>>fData[evno][fGoodInEvent[evno]].py;
+    in>>fData[evno][fGoodInEvent[evno]].pz;
+    in>>fData[evno][fGoodInEvent[evno]].x;
+    in>>fData[evno][fGoodInEvent[evno]].y;
+    in>>fData[evno][fGoodInEvent[evno]].z;
+    
+    in>>tmp;
+    in>>tmp;
+    in>>tmp;
+    
+ /* cout<<evno<<" ";
+  cout<<fData[evno][fGoodInEvent[evno]].lab;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
+  cout<<"\n";
+ */ 
+  fGoodInEvent[evno]++;
+ }
+ in.close();
+ cout<<"AliGoodTracksPP::AliGoodTracksPP()  ....  Done\n";
+}
+
+
+
+const GoodTrack& AliGoodTracksPP::GetTrack(Int_t event, Int_t n) const
+ {
+  
+  if( (event>fNevents) || (event<0))
+   {
+     gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Event %d",event);
+   }
+  if( (n>fGoodInEvent[event]) || (n<0))
+   {
+     gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Good TPC Track %d",n);
+   }
+  return fData[event][n];
+
+ }
diff --git a/HBTAN/AliHBTReaderPPprod.h b/HBTAN/AliHBTReaderPPprod.h
new file mode 100644 (file)
index 0000000..cf981e6
--- /dev/null
@@ -0,0 +1,78 @@
+#ifndef ALIHBTREADERPPPROD_H
+#define ALIHBTREADERPPPROD_H
+
+#include "AliHBTReader.h"
+#include "AliHBTReaderTPC.h"
+
+//This reader reads tracks AliTPCtracks.root
+//                  particles form tpc_good_tracks 
+//I am aware that this file is temporary however we do not have any other PID yet
+//Piotr.Skowronski@cern.ch
+
+#include <TString.h>
+class TFile;
+
+class AliHBTReaderPPprod: public AliHBTReader
+{
+  public:
+    AliHBTReaderPPprod(const Char_t* trackfilename = "AliTPCtracks.root",
+                       const Char_t* clusterfilename = "AliTPCclusters.root", 
+                       const Char_t* goodtracksfilename = "good_tracks_tpc",
+                       const Char_t* galicefilename = "");
+
+    virtual ~AliHBTReaderPPprod();
+    
+    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);//returns pointer to event with particles 
+    Int_t GetNumberOfPartEvents(); //returns number of particle events
+    Int_t GetNumberOfTrackEvents();//returns number of track events
+    
+  protected:
+    //in the future this class is will read global tracking
+
+    
+    Int_t OpenFiles(); //opens files to be read
+    void CloseFiles(); //close files
+    
+    AliHBTRun* fParticles; //!simulated particles
+    AliHBTRun* fTracks; //!reconstructed tracks (particles)
+
+    TString fTrackFileName; //name of the file with tracks
+    TString fClusterFileName;//name of the file with clusters
+    TString fGAliceFileName;//name of the file with galice.root
+    TString fGoodTPCTracksFileName; //name of text file with good tracks
+    
+    TFile *fTracksFile; //file with tracks
+    TFile *fClustersFile;//file with clusters
+    
+    Bool_t fIsRead; //flag indicating if the data are already read
+    
+    
+  private:
+  public:
+    ClassDef(AliHBTReaderPPprod,1)
+};
+
+
+class AliGoodTracksPP
+ { 
+   //container for good tracks
+   //this class is for internal use only
+   
+   friend class AliHBTReaderPPprod;
+   
+   private:
+     AliGoodTracksPP(const TString& infilename = TString("good_tracks_tpc")); //constructor
+     ~AliGoodTracksPP(); //dctor
+   
+     const GoodTrack& GetTrack(Int_t event, Int_t n) const; //returns reference to the nth good track in event "event"
+
+     Int_t  fNevents;  //Number of events
+     Int_t* fGoodInEvent; //Numbers of good track in event
+     struct GoodTrack **fData;
+ };
+
+
+#endif
diff --git a/HBTAN/AliHBTReaderTPC.cxx b/HBTAN/AliHBTReaderTPC.cxx
new file mode 100644 (file)
index 0000000..ccfa508
--- /dev/null
@@ -0,0 +1,427 @@
+#include "AliHBTReaderTPC.h"
+
+#include <iostream.h>
+#include <fstream.h>
+#include <TString.h>
+#include <TTree.h>
+#include <TFile.h>
+
+#include <AliTPCtrack.h>
+#include <AliTPCParam.h>
+#include <AliTPCtracker.h>
+
+#include "AliRun.h"
+#include "AliHBTRun.h"
+#include "AliHBTEvent.h"
+#include "AliHBTParticle.h"
+#include "AliHBTParticleCut.h"
+
+
+ClassImp(AliHBTReaderTPC)
+//reader for TPC tracking
+//needs galice.root, AliTPCtracks.root, AliTPCclusters.root, good_tracks_tpc 
+//
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+//Piotr.Skowronski@cern.ch
+
+AliHBTReaderTPC::
+ AliHBTReaderTPC(const Char_t* trackfilename,const Char_t* clusterfilename,
+                 const Char_t* goodtracksfilename,const Char_t* galicefilename):
+                 fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
+                 fGAliceFileName(galicefilename),
+                 fGoodTPCTracksFileName(goodtracksfilename)
+{
+  //constructor, only file names are set
+  //Defaults:
+  //  trackfilename = "AliTPCtracks.root"
+  //  clusterfilename = "AliTPCclusters.root"
+  //  goodtracksfilename = "good_tracks_tpc"
+  //  galicefilename = ""  - this means: Do not open gAlice file - 
+  //                         just leave the global pointer untached
+  
+  fParticles = new AliHBTRun();
+  fTracks    = new AliHBTRun();
+
+  fTracksFile   = 0x0;  //files are opened during reading only
+  fClustersFile = 0x0;
+  
+  fIsRead = kFALSE;
+}
+/********************************************************************/
+
+AliHBTReaderTPC::~AliHBTReaderTPC()
+ {
+ //desctructor
+   delete fParticles;
+   delete fTracks;
+ }
+/********************************************************************/
+
+AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fTracks->GetEvent(n);
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+   if (!fIsRead) Read(fParticles,fTracks);
+   return fParticles->GetNumberOfEvents();
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+  if (!fIsRead) Read(fParticles,fTracks);
+  return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+
+Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
+ {
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+  Int_t i; //iterator and some temprary values
+  Int_t Nevents;
+  if (!particles) //check if an object is instatiated
+   {
+     Error("AliHBTReaderTPC::Read"," particles object must instatiated before passing it to the reader");
+   }
+  if (!tracks)  //check if an object is instatiated
+   {
+     Error("AliHBTReaderTPC::Read"," tracks object must instatiated before passing it to the reader");
+   }
+  particles->Reset();//clear runs == delete all old events
+  tracks->Reset();
+    
+  if( (i=OpenFiles()) )
+   {
+     Error("AliHBTReaderTPC::Read","Exiting due to problems with opening files. Errorcode %d",i);
+     return i;
+   }
+  
+  AliGoodTracks *goodTPCTracks = new AliGoodTracks(fGoodTPCTracksFileName);
+  if (!goodTPCTracks)
+   {
+     Error("AliHBTReaderTPC::Read","Exiting due to problems with opening files. Errorcode %d",i);
+     return 1;
+   }
+  
+    
+  if (gAlice->TreeE())//check if tree E exists
+   {
+    Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
+    cout<<"Found "<<Nevents<<endl;
+   }
+  else
+   {//if not return an error
+     Error("AliHBTReaderPPprod::Read","Can not find Header tree (TreeE) in gAlice");
+     return 1;
+   }
+  
+  fClustersFile->cd();//set cluster file active 
+  AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
+  if (!TPCParam) 
+    { 
+     Error("AliHBTReaderTPC::Read","TPC parameters have not been found !\n");
+     return 1;
+    }
+
+  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
+  
+  for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
+   {
+     cout<<"Reading Event "<<currentEvent<<endl;
+     /**************************************/
+      /**************************************/
+       /**************************************/ 
+         fTracksFile->cd();//set track file active
+         
+         Char_t  treename[100];
+         sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
+   
+         TTree *tracktree=0;
+         
+         tracktree=(TTree*)fTracksFile->Get(treename);//get the tree 
+         if (!tracktree) //check if we got the tree
+          {//if not return with error
+            Error("AliHBTReaderTPC::Read","Can't get a tree with TPC tracks !\n"); 
+            return 1;
+          }
+   
+         TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
+         if (!trackbranch) ////check if we got the branch
+          {//if not return with error
+            Error("AliHBTReaderTPC::Read","Can't get a branch with TPC tracks !\n"); 
+            return 2;
+          }
+         Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks 
+         cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
+         //Copy tracks to array
+         
+         AliTPCtrack *iotrack=0;
+         
+         fClustersFile->cd();//set cluster file active 
+         AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);//create the tacker for this event
+         if (!tracker) //check if it has created succeffuly
+          {//if not return with error
+            Error("AliHBTReaderTPC::Read","Can't get a tracker !\n"); 
+            return 3;
+          }
+         tracker->LoadInnerSectors();
+         tracker->LoadOuterSectors();
+   
+         for (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
+          }
+         
+         fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
+         fTracksFile->Delete("tracks");//delete branch from memmory
+         delete tracker; //delete tracker
+         
+         tracker = 0x0;
+         trackbranch = 0x0;
+         tracktree = 0x0;
+
+         Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent]; //number of good tracks in the current event
+   
+         Double_t xk;
+         Double_t par[5];
+         Float_t phi, lam, pt;//angles and transverse momentum
+         Int_t label; //label of the current track
+         Bool_t found; //flag indicated wether we managed to match good_tpc_track with track
+   
+         for (i=0; i<ngood; i++) //loop over all good tracks
+          { 
+            const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i); //get ith goog track
+            
+            if(Pass(gt.code)) continue; //check if we are intersted with particles of this type 
+                                        //if not take next partilce
+            
+            label = gt.lab;
+            found = kFALSE; //guard in case we don't find track with such a label
+            for (Int_t j=0;j<NTPCtracks;j++)//lopp over all tpc tracks
+              {
+                iotrack = (AliTPCtrack*)tarray->At(j);
+                if (iotrack->GetLabel() == label) //if the label is the same 
+                  {
+                    found = kTRUE; //we found the track
+                    break;
+                  }
+              }  
+            if(!found) //check if we found the track
+              {
+                Warning("Read",
+                "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
+                continue; //put comunicate on the screen and continue loop
+              }
+        
+            Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();//CMS mass of this particle 
+            Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass); //particle total energy
+            
+            AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
+            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(gt.x);
+            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 tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
+            
+            AliHBTParticle* track = new AliHBTParticle(gt.code, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+            if(Pass(track)) { delete  track;continue;}//check if meets all criteria of any of our cuts
+                                                      //if it does not delete it and take next good track
+            
+            particles->AddParticle(currentEvent,part);//put track and particle on the run
+            tracks->AddParticle(currentEvent,track);
+
+          }
+         tarray->Clear(); //clear the array
+    
+       /**************************************/
+      /**************************************/
+     /**************************************/  
+   }
+  
+  //save environment (resouces) --
+  //clean your place after the work
+  CloseFiles(); 
+  delete tarray;
+  delete goodTPCTracks;
+  fIsRead = kTRUE;
+  return 0;
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderTPC::OpenFiles()
+{
+ //opens all the files
+   fTracksFile = 0;
+   fTracksFile=TFile::Open(fTrackFileName.Data());
+   if (!fTracksFile->IsOpen()) 
+     {
+       Error("AliHBTReaderTPC::OpenFiles","Can't open file with tacks named ",fTrackFileName.Data());
+       return 1;
+     }
+   
+   fClustersFile = 0;
+   
+   fClustersFile=TFile::Open(fClusterFileName.Data());
+   if (!fClustersFile->IsOpen()) 
+    {
+      Error("AliHBTReaderTPC::OpenFiles","Can't open file with TPC clusters named ",fClusterFileName.Data());
+      return 2;
+    }
+
+ return 0; 
+}
+  
+
+
+/********************************************************************/
+  
+void AliHBTReaderTPC::CloseFiles()
+{
+  //closes the files
+  fTracksFile->Close();
+  fClustersFile->Close();
+}
+
+/********************************************************************/
+
+/********************************************************************/
+/********************************************************************/
+/********************************************************************/
+
+
+AliGoodTracks::~AliGoodTracks()
+{
+//destructor
+ delete [] fGoodInEvent;
+ for (Int_t i = 0;i<fNevents;i++)
+   delete [] fData[i];
+ delete [] fData;
+}
+/********************************************************************/
+AliGoodTracks::AliGoodTracks(const TString& infilename)
+{
+
+  cout<<"AliGoodTracks::AliGoodTracks()  ....\n";
+  if(!gAlice) 
+    {
+      cerr<<"There is no gAlice"<<endl;
+      delete this;
+      return;
+    }
+  
+  if (!gAlice->TreeE())
+   {
+     cerr<<"Can not find Header tree (TreeE) in gAlice"<<endl;
+     delete this;
+     return;
+   }
+   
+  fNevents = (Int_t)gAlice->TreeE()->GetEntries();
+  //fNevents = 100;
+  cout<<fNevents<<" FOUND"<<endl;
+  ifstream in(infilename.Data());
+
+  if(!in)
+    {
+      cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
+      delete this;
+      return;
+    }
+
+  
+  fGoodInEvent = new Int_t[fNevents];
+  fData = new struct GoodTrack* [fNevents];
+
+  Int_t i;
+  for( i = 0;i<fNevents;i++)
+   {
+    fGoodInEvent[i] =0;
+    fData[i] = new struct GoodTrack[50000];
+   }
+
+  Int_t evno;
+  while(in>>evno)
+   {
+    if(fGoodInEvent[evno]>=50000)
+     {
+      cerr<<"AliGoodTracks::AliGoodTracks() : Not enough place in the array\n";
+      continue;
+     }
+    in>>fData[evno][fGoodInEvent[evno]].lab;
+    in>>fData[evno][fGoodInEvent[evno]].code;
+    in>>fData[evno][fGoodInEvent[evno]].px;
+    in>>fData[evno][fGoodInEvent[evno]].py;
+    in>>fData[evno][fGoodInEvent[evno]].pz;
+    in>>fData[evno][fGoodInEvent[evno]].x;
+    in>>fData[evno][fGoodInEvent[evno]].y;
+    in>>fData[evno][fGoodInEvent[evno]].z;
+    
+ /* cout<<evno<<" ";
+  cout<<fData[evno][fGoodInEvent[evno]].lab;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
+  cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
+  cout<<"\n";
+ */ 
+  fGoodInEvent[evno]++;
+ }
+ in.close();
+ cout<<"AliGoodTracks::AliGoodTracks()  ....  Done\n";
+}
+
+
+
+const GoodTrack& AliGoodTracks::GetTrack(Int_t event, Int_t n) const
+ {
+  
+  if( (event>fNevents) || (event<0))
+   {
+     gAlice->Fatal("AliGoodTracks::GetTrack","No such Event %d",event);
+   }
+  if( (n>fGoodInEvent[event]) || (n<0))
+   {
+     gAlice->Fatal("AliGoodTracks::GetTrack","No such Good TPC Track %d",n);
+   }
+  return fData[event][n];
+
+ }
diff --git a/HBTAN/AliHBTReaderTPC.h b/HBTAN/AliHBTReaderTPC.h
new file mode 100644 (file)
index 0000000..f98d7e3
--- /dev/null
@@ -0,0 +1,85 @@
+#ifndef ALIHBTREADERTPC_H
+#define ALIHBTREADERTPC_H
+
+#include "AliHBTReader.h"
+
+
+//This reader reads tracks AliTPCtracks.root
+//                  particles form tpc_good_tracks 
+//I am aware that this file is temporary however we do not have any other PID
+//Piotr.Skowronski@cern.ch
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+#include <TString.h>
+class TFile;
+
+class AliHBTReaderTPC: public AliHBTReader
+{
+  public:
+    AliHBTReaderTPC(const Char_t* trackfilename = "AliTPCtracks.root",
+                    const Char_t* clusterfilename = "AliTPCclusters.root",
+       const Char_t* goodtracksfilename = "good_tracks_tpc",
+       const Char_t* galicefilename = "");
+       
+    virtual ~AliHBTReaderTPC();
+    
+    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);//returns pointer to event with particles 
+    Int_t GetNumberOfPartEvents();//returns number of particle events
+    Int_t GetNumberOfTrackEvents();//returns number of track events
+    
+  protected:
+    //in the future this class is will read global tracking
+
+    
+    Int_t OpenFiles();//opens files to be read
+    void CloseFiles();//close files
+    
+    AliHBTRun* fParticles; //!simulated particles
+    AliHBTRun* fTracks; //!reconstructed tracks (particles)
+
+    TString fTrackFileName;//name of the file with tracks
+    TString fClusterFileName;//name of the file with clusters
+    TString fGAliceFileName;//name of the file with galice.root
+    TString fGoodTPCTracksFileName;//name of text file with good tracks
+    
+    TFile *fTracksFile;//file with tracks
+    TFile *fClustersFile;//file with clusters
+    
+    Bool_t fIsRead;//flag indicating if the data are already read
+    
+    
+  private:
+  public:
+    ClassDef(AliHBTReaderTPC,1)
+};
+
+
+struct GoodTrack //data of good tracks produced by AliTPCComparison.C  
+ {
+  Int_t lab;
+  Int_t code;
+  Float_t px,py,pz;
+  Float_t x,y,z;
+ };
+
+class AliGoodTracks
+ { 
+   //this class is for internal use only
+   friend class AliHBTReaderTPC;
+   
+   private:
+     AliGoodTracks(const TString& infilename = TString("good_tracks_tpc"));
+     ~AliGoodTracks();
+   
+     const GoodTrack& GetTrack(Int_t event, Int_t n) const;
+
+     Int_t  fNevents;  //Number of events
+     Int_t* fGoodInEvent; //Numbers of good track in event
+     struct GoodTrack **fData;
+ };
+
+
+#endif
diff --git a/HBTAN/AliHBTRun.cxx b/HBTAN/AliHBTRun.cxx
new file mode 100644 (file)
index 0000000..4c8687e
--- /dev/null
@@ -0,0 +1,65 @@
+#include "AliHBTRun.h"
+
+#include <TObjArray.h>
+
+
+ClassImp(AliHBTRun)
+/**************************************************************************/ 
+
+AliHBTRun::AliHBTRun()
+  { 
+  //contructor
+  
+    fEvents = new TObjArray();//create array for AliHBTEvents
+    if(!fEvents) Fatal("AliHBTRun::AliHBTRun","Can not allocate memory");
+    fEvents->SetOwner(); //array is an owner: when is deleted or cleared it deletes objects that it contains
+  }
+/**************************************************************************/
+AliHBTRun::~AliHBTRun()
+  {
+    //destructor
+    
+    delete fEvents;//delete array with events
+  }
+
+
+/**************************************************************************/
+void AliHBTRun::Reset()
+ { 
+   fEvents->Clear();//clear an array with events. 
+                    //All events are deleted because array is an owner (set in constructor)
+ }
+/**************************************************************************/
+
+
+void AliHBTRun::AddParticle(Int_t event, AliHBTParticle* part)
+{
+ //Adds particle to event
+ //if there is no event of this number, crate it and add to the collection
+ if(!GetEvent(event))  fEvents->AddAtAndExpand(new AliHBTEvent, event);
+ GetEvent(event)->AddParticle(part);
+}
+/**************************************************************************/
+
+
+void AliHBTRun::AddParticle(Int_t event, TParticle* part)
+{
+ //if there is no event of this number, crate it and add to the collection
+ if(!GetEvent(event))  fEvents->AddAtAndExpand(new AliHBTEvent, event);
+ GetEvent(event)->AddParticle(part);
+}
+/**************************************************************************/
+
+
+void AliHBTRun::AddParticle(Int_t event, Int_t pdg, 
+                            Double_t px, Double_t py, Double_t pz, Double_t etot,
+                            Double_t vx, Double_t vy, Double_t vz, Double_t time)
+{
+ //if there is no event of this number, crate it and add to the collection
+ if(!GetEvent(event))  fEvents->AddAtAndExpand(new AliHBTEvent, event);
+ GetEvent(event)->AddParticle(pdg,px,py,pz,etot,vx,vy,vz,time);
+}
+
+/**************************************************************************/ 
diff --git a/HBTAN/AliHBTRun.h b/HBTAN/AliHBTRun.h
new file mode 100644 (file)
index 0000000..c3408f1
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef ALIHBTRUN_H
+#define ALIHBTRUN_H
+
+#include "AliHBTEvent.h"
+#include <TObjArray.h>
+
+//class describing set of events (the run)
+//designed for fast acces
+//Piotr.Skowronski@cern.ch
+
+class AliHBTParticle;
+
+class AliHBTEvent;
+
+class AliHBTRun: public TObject
+ {
+  public:
+    AliHBTRun();
+    virtual ~AliHBTRun();
+
+    void            AddParticle(Int_t event, AliHBTParticle* part); //inerface to AliHBTEvent::AddParticle(AliHBTParticle*) 
+    void            AddParticle(Int_t event, TParticle* part);//inerface to AliHBTEvent::AddParticle(TParticle*) 
+    
+    //inerface to AliHBTEvent::AddParticle(Int_t.Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t)
+    void            AddParticle(Int_t event, Int_t pdg, 
+                                Double_t px, Double_t py, Double_t pz, Double_t etot,
+                                Double_t vx, Double_t vy, Double_t vz, Double_t time); 
+    
+    AliHBTParticle* GetParticle(Int_t event, Int_t n); //returns nth particle from event
+    AliHBTEvent*    GetEvent(Int_t event) const; //returns AliHBTEvent number "event"
+    
+    Int_t           GetNumberOfEvents() const; //returns number of events
+    Int_t              GetNumberOfParticlesInEvent(Int_t event) const; //returns number of particles in event number "event"
+    void            Reset();//clears all events in the array (deletes)
+  protected:
+    TObjArray* fEvents;//!Array containig AliHBTEvents
+  private:
+    
+  public:
+    ClassDef(AliHBTRun,1)
+ };
+/**************************************************************************/
+
+inline
+AliHBTEvent* AliHBTRun::GetEvent(Int_t event) const
+{
+//returns pointer to AliHBTEvent number "event"
+  
+  return (AliHBTEvent*)fEvents->At(event);
+}
+/**************************************************************************/
+inline
+AliHBTParticle* AliHBTRun::GetParticle(Int_t event, Int_t n) 
+{
+ //returns nth particle from event number event
+  return GetEvent(event)->GetParticle(n);
+}
+
+/**************************************************************************/
+
+inline
+Int_t AliHBTRun::GetNumberOfEvents() const
+ {
+//returns number of events in collection
+
+   return fEvents->GetEntriesFast(); //there may be empty slots but we do not care
+                                     //Analysis checks it if return is not NULL
+ }
+/**************************************************************************/
+
+inline
+Int_t AliHBTRun::GetNumberOfParticlesInEvent(Int_t event) const
+{
+//returns number of Particles in event 
+  return GetEvent(event)->GetNumberOfParticles();
+}
+#endif
diff --git a/HBTAN/HBTAnalysisLinkDef.h b/HBTAN/HBTAnalysisLinkDef.h
new file mode 100644 (file)
index 0000000..0e2ca8c
--- /dev/null
@@ -0,0 +1,65 @@
+#ifdef __CINT__
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+#pragma link C++ class AliHBTAnalysis+;
+#pragma link C++ class AliHBTParticle+;
+#pragma link C++ class AliHBTEvent+;
+#pragma link C++ class AliHBTRun+;
+#pragma link C++ class AliHBTFunction+;
+#pragma link C++ class AliHBTTwoPartFctn+;
+#pragma link C++ class AliHBTFourPartFctn+;
+#pragma link C++ class AliHBTTwoPartFctn1D+;
+#pragma link C++ class AliHBTTwoPartFctn2D+;
+#pragma link C++ class AliHBTTwoPartFctn3D+;
+
+#pragma link C++ class AliHBTFourPartFctn2D+;
+
+#pragma link C++ class AliHBTPair+;
+#pragma link C++ class AliHBTParticleCut-;
+#pragma link C++ class AliHBTEmptyParticleCut-;
+#pragma link C++ class AliHBTPairCut-;
+#pragma link C++ class AliHBTEmptyPairCut-;
+
+#pragma link C++ class AliHbtBaseCut+;
+#pragma link C++ class AliHbtBasePairCut+;
+#pragma link C++ class AliHBTQInvCut+;
+
+#pragma link C++ class  AliHBTMomentumCut+;
+#pragma link C++ class  AliHBTPtCut+;
+#pragma link C++ class  AliHBTEnergyCut+;
+#pragma link C++ class  AliHBTRapidityCut+;
+#pragma link C++ class  AliHBTPseudoRapidityCut+;
+#pragma link C++ class  AliHBTPxCut+;
+#pragma link C++ class  AliHBTPyCut+;
+#pragma link C++ class  AliHBTPzCut+;
+#pragma link C++ class  AliHBTPhiCut+;
+#pragma link C++ class  AliHBTThetaCut+;
+#pragma link C++ class  AliHBTVxCut+;
+#pragma link C++ class  AliHBTVyCut+;
+#pragma link C++ class  AliHBTVzCut+;
+
+
+#pragma link C++ class AliHBTReader+;
+#pragma link C++ class AliHBTReaderTPC+;
+#pragma link C++ class AliHBTReaderPPprod+;
+
+#pragma link C++ class AliHBTQInvCorrelFctn+;
+#pragma link C++ class AliHBTInvMassCorrelFctn+;
+
+#pragma link C++ class AliHBTQOutResolVSQInvFctn+;
+#pragma link C++ class AliHBTQSideResolVSQInvFctn+;
+#pragma link C++ class AliHBTQLongResolVSQInvFctn+;
+#pragma link C++ class AliHBTQInvResolVSKtFctn+;
+#pragma link C++ class AliHBTQOutResolVSKtFctn+;
+#pragma link C++ class AliHBTQSideResolVSKtFctn+;
+#pragma link C++ class AliHBTQLongResolVSKtFctn+;
+#pragma link C++ class AliHBTQOutResolVSQOutFctn+;
+#pragma link C++ class AliHBTQSideResolVSQSideFctn+;
+#pragma link C++ class AliHBTQLongResolVSQLongFctn+;
+
+
+
+#endif
diff --git a/HBTAN/Makefile b/HBTAN/Makefile
new file mode 100644 (file)
index 0000000..9c4692d
--- /dev/null
@@ -0,0 +1,88 @@
+############################### HBTAnalysis Makefile ##################################
+
+# Include machine specific definitions
+
+include $(ALICE_ROOT)/conf/GeneralDef
+include $(ALICE_ROOT)/conf/MachineDef.$(ALICE_TARGET)
+
+
+OPT = -g
+PACKAGE = HBTAnalysis
+
+# C++ sources
+
+SRCS          = AliHBTAnalysis.cxx AliHBTFunction.cxx \
+                AliHBTEvent.cxx AliHBTRun.cxx \
+               AliHBTParticle.cxx AliHBTParticleCut.cxx \
+               AliHBTPair.cxx AliHBTPairCut.cxx\
+                AliHBTCorrelFctn.cxx \
+                AliHBTReader.cxx AliHBTReaderTPC.cxx\
+                AliHBTQResolutionFctns.cxx AliHBTReaderPPprod.cxx
+
+# Fortran sources
+
+FSRCS         = 
+
+# C++ Headers
+
+HDRS          = $(SRCS:.cxx=.h) HBTAnalysisLinkDef.h
+
+# Library dictionary
+
+DICT          = HBTAnalysisCint.cxx
+DICTH         = $(DICT:.cxx=.h)
+DICTO         = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(DICT))
+
+# FORTRAN Objectrs
+
+FOBJS         = $(patsubst %.F,tgt_$(ALICE_TARGET)/%.o,$(FSRCS))
+
+# C Objects
+
+COBJS         = $(patsubst %.c,tgt_$(ALICE_TARGET)/%.o,$(CSRCS))
+
+# C++ Objects
+
+OBJS          = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(SRCS)) $(DICTO)
+
+# C++ compilation flags
+
+INCLUDES      = -I$(ALICE_ROOT)/TPC/ -I$(ALICE_ROOT)/CONTAINERS/
+
+CXXFLAGS      = $(CXXOPTS) -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include/ $(INCLUDES)
+# FORTRAN compilation flags
+
+FFLAGS      = $(FOPT) -I$(ALICE_ROOT)/GEANT321
+
+##### TARGETS #####
+# Target
+
+SLIBRARY       = $(LIBDIR)/libHBTAnalysis.$(SL)
+ALIBRARY       = $(LIBDIR)/libHBTAnalysis.a
+
+default:       $(SLIBRARY)
+
+
+$(LIBDIR)/libHBTAnalysis.$(SL):                $(OBJS)
+
+$(DICT):                       $(HDRS)
+
+DEPINC += -I$(ALICE_ROOT)/GEANT321
+
+depend:                                $(SRCS) $(FSRCS)
+
+TOCLEAN                = $(OBJS) $(FOBJS) *Cint.cxx *Cint.h
+
+CHECKS        = $(patsubst %.cxx,check/%.viol,$(SRCS))
+
+############################### General Macros ################################
+
+include $(ALICE_ROOT)/conf/GeneralMacros
+
+############################ Dependencies #####################################
+
+-include tgt_$(ALICE_TARGET)/Make-depend 
+
+
diff --git a/HBTAN/libHBTAnalysis.pkg b/HBTAN/libHBTAnalysis.pkg
new file mode 100644 (file)
index 0000000..3382ba5
--- /dev/null
@@ -0,0 +1,14 @@
+SRCS          =   \
+AliHBTAnalysis.cxx    AliHBTPair.cxx         AliHBTQResolutionFctns.cxx  AliHBTReaderPPprod.cxx\
+AliHBTCorrelFctn.cxx  AliHBTPairCut.cxx      AliHBTReader.cxx            AliHBTReaderTPC.cxx\
+AliHBTEvent.cxx       AliHBTParticle.cxx     AliHBTRun.cxx\
+AliHBTFunction.cxx    AliHBTParticleCut.cxx  
+
+
+HDRS:= $(SRCS:.cxx=.h) 
+
+DHDR= HBTAnalysisLinkDef.h
+
+EINCLUDE:= TPC CONTAINERS 
+
+