Corrections for gcc 4.0
[u/mrichter/AliRoot.git] / HBTAN / AliHBTAnalysis.cxx
index 56e4bc7..b330210 100644 (file)
-
 #include "AliHBTAnalysis.h"
+//_________________________________________________________
+////////////////////////////////////////////////////////////////////////////
+//
+// class AliHBTAnalysis
+//
+// Central Object Of HBTAnalyser: 
+// This class performs main looping within HBT Analysis
+// User must plug a reader of Type AliReader
+// User plugs in coorelation and monitor functions
+// as well as monitor functions
+//
+// HBT Analysis Tool, which is integral part of AliRoot,
+// ALICE Off-Line framework:
+//
+// Piotr.Skowronski@cern.ch
+// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
+//
+////////////////////////////////////////////////////////////////////////////
+//_________________________________________________________
 
-#include <iostream.h>
+#include <Riostream.h>
+#include <AliESD.h>
 
-#include "AliHBTRun.h"
-#include "AliHBTReader.h"
-#include "AliHBTParticle.h"
-#include "AliHBTParticleCut.h"
-#include "AliHBTPair.h"
-#include "AliHBTPairCut.h"
-#include "AliHBTFunction.h"
+#include <TSystem.h>
+#include <TFile.h>
 
-#include <TList.h>
+#include "AliAOD.h"
+#include "AliAODParticle.h"
+#include "AliAODPairCut.h"
+#include "AliEventCut.h"
 
+#include "AliEventBuffer.h"
 
+#include "AliReader.h"
+#include "AliHBTPair.h"
+#include "AliHBTFunction.h"
+#include "AliHBTMonitorFunction.h"
+#include "AliLog.h"
 
 ClassImp(AliHBTAnalysis)
 
 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
-const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
+const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
+const Int_t  AliHBTAnalysis::fgkDefaultBufferSize = 5;
 
-AliHBTAnalysis::AliHBTAnalysis()
+AliHBTAnalysis::AliHBTAnalysis():
+  fProcEvent(0x0),
+  fReader(0x0),
+  fNTrackFunctions(0),
+  fNParticleFunctions(0),
+  fNParticleAndTrackFunctions(0),
+  fNTrackMonitorFunctions(0),
+  fNParticleMonitorFunctions(0), 
+  fNParticleAndTrackMonitorFunctions(0),
+  fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
+  fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
+  fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
+  fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),    
+  fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),    
+  fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),    
+  fBkgEventCut(0x0),
+  fPartBuffer(0x0),
+  fTrackBuffer(0x0),
+  fBufferSize(2),
+  fDisplayMixingInfo(fgkDefaultMixingInfo),
+  fIsOwner(kFALSE),
+  fProcessOption(kSimulatedAndReconstructed),
+  fNoCorrfctns(kFALSE),
+  fOutputFileName(0x0),
+  fVertexX(0.0),
+  fVertexY(0.0),
+  fVertexZ(0.0)
  {
-   fReader = 0x0;
-   
-   fTrackFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
-   fParticleFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
-   fParticleAndTrackFunctions = new AliHBTTwoPairFctn* [fgkFctnArraySize];
-   
-   fNTrackFunctions = 0;
-   fNParticleFunctions = 0;
-   fNParticleAndTrackFunctions = 0;
-   
-   fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
+   //default constructor
    
  }
 /*************************************************************************************/ 
 
+AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
+  AliAnalysis(in),
+  fProcEvent(0x0),
+  fReader(0x0),
+  fNTrackFunctions(0),
+  fNParticleFunctions(0),
+  fNParticleAndTrackFunctions(0),
+  fNTrackMonitorFunctions(0),
+  fNParticleMonitorFunctions(0), 
+  fNParticleAndTrackMonitorFunctions(0),
+  fTrackFunctions(0x0),
+  fParticleFunctions(0x0),
+  fParticleAndTrackFunctions(0x0),
+  fParticleMonitorFunctions(0x0),
+  fTrackMonitorFunctions(0x0),
+  fParticleAndTrackMonitorFunctions(0x0),
+  fBkgEventCut(0x0),
+  fPartBuffer(0x0),
+  fTrackBuffer(0x0),
+  fBufferSize(fgkDefaultBufferSize),
+  fDisplayMixingInfo(fgkDefaultMixingInfo),
+  fIsOwner(kFALSE),
+  fProcessOption(kSimulatedAndReconstructed),
+  fNoCorrfctns(kFALSE),
+  fOutputFileName(0x0),
+  fVertexX(0.0),
+  fVertexY(0.0),
+  fVertexZ(0.0)
+ {
+//copy constructor
+   Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
+ }
+/*************************************************************************************/ 
+AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
+ {
+//operator =
+   Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
+   return *this;
+ }
+/*************************************************************************************/ 
 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
+    
+   if (fIsOwner)
+    {
+      AliDebug(5,"Is Owner: Attempting to delete functions");
+      DeleteFunctions();
+      if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
+    }
    delete [] fTrackFunctions;
    delete [] fParticleFunctions;
    delete [] fParticleAndTrackFunctions;
    
-   delete fPairCut; // always have an copy of an object - we create we dstroy
+   delete [] fParticleMonitorFunctions; 
+   delete [] fTrackMonitorFunctions; 
+   delete [] fParticleAndTrackMonitorFunctions;
+
+   delete fBkgEventCut;
+   delete fOutputFileName;
  }
 
 /*************************************************************************************/ 
 
-void AliHBTAnalysis::Process(Option_t* option)
+Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
 {
- //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("Process","The reader is not set");
-   return;
-  }
+  //Processes one event
+  if (fProcEvent == 0x0)
+   {
+     Error("ProcessEvent","Analysis <<%s>> option not specified.",GetName());
+     return 1;
+   }
+  if ( Rejected(aodrec,aodsim) ) 
+   {
+//     Double_t x = 0, y = 0, z = 0;
+//     aodrec->GetPrimaryVertex(x,y,z);
+//     Info("ProcessEvent","Event has vertex at %f %f %f",x,y,z);
+     Info("ProcessEvent","Nch is %d",aodsim->GetNumberOfCharged());
+     Info("ProcessEvent","%s: Event cut rejected this event",GetName());
+     return 0;
+   } 
+  
+  //Move event to the apparent vertex -> must be after the event cut
+  //It is important for any cut that use any spacial coordiantes, 
+  //f.g. anti-merging cut in order to preserve the same bahavior of variables (f.g. distance between tracks)
+  Double_t dvx = 0, dvy = 0, dvz = 0;
+  if (aodrec)
+   {
+     Double_t pvx,pvy,pvz;
+     aodrec->GetPrimaryVertex(pvx,pvy,pvz);
+     
+     dvx = fVertexX - pvx;
+     dvy = fVertexY - pvy;
+     dvz = fVertexZ - pvz;
+     aodrec->Move(dvx,dvy,dvz);
+   }  
+  
+  Int_t result = (this->*fProcEvent)(aodrec,aodsim);
+
+  if (aodrec) aodrec->Move(-dvx,-dvy,-dvz);//move event back to the same spacial coordinates
+  
+  return  result;
+}
+/*************************************************************************************/ 
+
+Int_t AliHBTAnalysis::Finish()
+{
+//Finishes analysis
+  WriteFunctions();
+  return 0;
+}
+/*************************************************************************************/ 
+
+void AliHBTAnalysis::DeleteFunctions()
+{
+ //Deletes all functions added to analysis
  
+ UInt_t ii;
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+  { 
+    if (AliVAODParticle::GetDebug()>5)
+     {
+       Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]);
+       Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name());
+     }
+    delete fParticleFunctions[ii];
+  } 
+ fNParticleFunctions = 0;
+                
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+  { 
+    if (AliVAODParticle::GetDebug()>5)
+     {
+       Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]);
+       Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name());
+     }
+    delete fTrackFunctions[ii];
+  }  
+ fNTrackFunctions = 0;
  
- const char *oT = strstr(option,"Tracks");
- const char *oP = strstr(option,"Particles");
+ for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+  { 
+    if (AliVAODParticle::GetDebug()>5)
+     {
+       Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
+       Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
+     } 
+    delete fParticleAndTrackFunctions[ii];
+  }  
+ fNParticleAndTrackFunctions = 0;
  
- if(oT && oP)
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
   { 
-    if (fReader->GetNumberOfPartEvents() <1)
+    if (AliVAODParticle::GetDebug()>5)
      {
-       Error("Process","There is no Particles. Maybe change the option?");
-       return;
+       Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
+       Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
      }
-    if (fReader->GetNumberOfTrackEvents() <1)
+    delete fParticleMonitorFunctions[ii];
+  } 
+ fNParticleMonitorFunctions = 0;
+   
+ for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+  { 
+    if (AliVAODParticle::GetDebug()>5)
      {
-       Error("Process","There is no Tracks. Maybe change the option?");
-       return;
+       Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
+       Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
      }
-    
-    if ( RunCoherencyCheck() )
-      {
-        Error("Process",
-              "Coherency check not passed. Maybe change the option?\n");
-        return;
-      }
-    ProcessTracksAndParticles();
-    return;
-  }
-
- if(oT)
-  {
-    if (fReader->GetNumberOfTrackEvents() <1)
+    delete fTrackMonitorFunctions[ii];
+  } 
+ fNTrackMonitorFunctions = 0;
+   
+ for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
+  { 
+    if (AliVAODParticle::GetDebug()>5)
      {
-       Error("Process","There is no data to analyze.");
-       return;
+      Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
+      Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
      }
-    ProcessTracks();
-    return;
-  }
+    delete fParticleAndTrackMonitorFunctions[ii];
+  } 
+ fNParticleAndTrackMonitorFunctions = 0;
  
- if(oP)
+}
+/*************************************************************************************/ 
+
+Int_t AliHBTAnalysis::Init()
+{
+//Initializeation method
+//calls Init for all functions
+
+ //Depending on option and pair cut assigns proper analysis method
+ Bool_t nonid = IsNonIdentAnalysis();
+ switch(fProcessOption)
   {
-    if (fReader->GetNumberOfPartEvents() <1)
-     {
-       Error("Process","There is no data to analyze.");
-       return;
-     }
-    ProcessParticles();
-    return;
+    case kReconstructed:
+      if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecNonId;
+      else fProcEvent = &AliHBTAnalysis::ProcessRec;
+      SetCutsOnRec();
+      break;
+
+    case kSimulated:
+      if (nonid) fProcEvent = &AliHBTAnalysis::ProcessSimNonId;
+      else fProcEvent = &AliHBTAnalysis::ProcessSim;
+      SetCutsOnSim();
+      break;
+
+    case kSimulatedAndReconstructed:
+      if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecAndSimNonId;
+      else fProcEvent = &AliHBTAnalysis::ProcessRecAndSim;
+      break;
   }
+  
+ if (fPartBuffer == 0x0)  fPartBuffer = new AliEventBuffer (fBufferSize);
+ else fPartBuffer->Reset();
+
+ if (fTrackBuffer == 0x0)  fTrackBuffer = new AliEventBuffer (fBufferSize);
+ else fTrackBuffer->Reset();
+
+
+ fNoCorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
  
-}
+ UInt_t ii;
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+   fParticleFunctions[ii]->Init();
+                
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+   fTrackFunctions[ii]->Init();
 
+ for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+   fParticleAndTrackFunctions[ii]->Init();
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+   fParticleMonitorFunctions[ii]->Init();
+   
+ for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+   fTrackMonitorFunctions[ii]->Init();
+   
+ for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
+   fParticleAndTrackMonitorFunctions[ii]->Init();
+   
+ return 0;
+}
 /*************************************************************************************/ 
 
-void AliHBTAnalysis::ProcessTracksAndParticles()
+void AliHBTAnalysis::ResetFunctions()
 {
-
-//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
-//the loops are splited
+//In case fOwner is true, deletes all functions
+//in other case, just set number of analysis to 0
+ if (fIsOwner) DeleteFunctions();
+ else
+  {
+    fNParticleFunctions = 0;
+    fNTrackFunctions = 0;
+    fNParticleAndTrackFunctions = 0;
+    fNParticleMonitorFunctions = 0;
+    fNTrackMonitorFunctions = 0;
+    fNParticleAndTrackMonitorFunctions = 0;
+  }
+}
+/*************************************************************************************/ 
+Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim)
+{
+//Does analysis for both tracks and particles
+//mainly for resolution study and analysies with weighting algirithms
   
+// cut on particles only -- why?
+// - PID: when we make resolution analysis we want to take only tracks with correct PID
+// We need cut on tracks because there are data characteristic
   
-  AliHBTParticle * part1, * part2;
-  AliHBTParticle * track1, * track2;
+  AliVAODParticle * part1, * part2;
+  AliVAODParticle * track1, * track2;
   
-  AliHBTEvent * trackEvent, *partEvent;
-  AliHBTEvent * trackEvent2,*partEvent2;
+  AliAOD * trackEvent = aodrec, *partEvent = aodsim;
+  AliAOD* trackEvent1 = new AliAOD();
+  AliAOD* partEvent1 = new AliAOD();
+
+  AliAOD * 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();
+//  Int_t nev = fReader->GetNumberOfTrackEvents();
+  static AliHBTPair tpair;
+  static AliHBTPair ppair;
 
+  AliHBTPair* trackpair = &tpair;
+  AliHBTPair* partpair = &ppair;
   AliHBTPair * tmptrackpair;//temprary pointers to pairs
   AliHBTPair * tmppartpair;
   
-  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%100) == 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++)
+  register UInt_t ii;
+  
+  
+
+  if ( !partEvent || !trackEvent ) 
+   {
+     Error("ProcessRecAndSim","<<%s>> Can not get event",GetName());
+     return 1;
+   }
+
+  if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
+   {
+     Error("ProcessRecAndSim",
+           "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d). Skipping Event.",
+            partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
+     return 2;
+   }
+
+
+  for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+   {
+     /***************************************/
+     /******   Looping same events   ********/
+     /******   filling numerators    ********/
+     /***************************************/
+     if ( (j%fDisplayMixingInfo) == 0)
+        Info("ProcessRecAndSim",
+             "Mixing particle %d with particles from the same event",j);
+
+     part1= partEvent->GetParticle(j);
+     track1= trackEvent->GetParticle(j);
+     
+     Bool_t firstcut = (this->*fkPass1)(part1,track1);
+     if (fBufferSize != 0) 
+       if ( (firstcut == kFALSE) || ( (this->*fkPass2)(part1,track1) == kFALSE ) )
+        {
+          //accepted by any cut
+          // we have to copy because reader keeps only one event
+
+          partEvent1->AddParticle(part1);
+          trackEvent1->AddParticle(track1);
+        }
+
+     if (firstcut) continue;
+
+     for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+       fParticleMonitorFunctions[ii]->Process(part1);
+     for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+       fTrackMonitorFunctions[ii]->Process(track1);
+     for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
+       fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
+
+     if (fNoCorrfctns) continue;
+
+     for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+      {
+        part2= partEvent->GetParticle(k);
+        if (part1->GetUID() == part2->GetUID()) continue;
+        partpair->SetParticles(part1,part2);
+
+        track2= trackEvent->GetParticle(k);
+        trackpair->SetParticles(track1,track2);
+
+        if( (this->*fkPass)(partpair,trackpair) ) //check pair cut 
+          { //do not meets crietria of the pair cut, try with swapped pairs
+            if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
+              continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
+            else 
+             { //swaped pair meets all the criteria
+               tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
+               tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
+             }
+          }
+        else
+         {//meets criteria of the pair cut
+           tmptrackpair = trackpair;
+           tmppartpair = partpair;
+         }
+
+        for(ii = 0;ii<fNParticleFunctions;ii++)
+               fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
+
+        for(ii = 0;ii<fNTrackFunctions;ii++)
+               fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
+
+        for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+               fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
+       //end of 2nd loop over particles from the same event  
+      }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+
+     /***************************************/
+     /***** Filling denominators    *********/
+     /***************************************/
+     if (fBufferSize == 0) continue;
+
+     fPartBuffer->ResetIter();
+     fTrackBuffer->ResetIter();
+     Int_t m = 0;
+     while (( partEvent2 = fPartBuffer->Next() ))
+      {
+        trackEvent2 = fTrackBuffer->Next();
+
+        m++;
+        if ( (j%fDisplayMixingInfo) == 0)
+           Info("ProcessRecAndSim",
+                "Mixing particle %d from current event with particles from event %d",j,-m);
+
+        for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
           {
-            part2= partEvent->GetParticle(k);
+            part2= partEvent2->GetParticle(l);
             partpair->SetParticles(part1,part2);
-           
-            track2= trackEvent->GetParticle(k);        
+
+            track2= trackEvent2->GetParticle(l);
             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
+            if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
+              { //do not meets crietria of the 
+                if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
+                  continue;
                 else 
-                 { //swaped pair meets all the criteria
-                   tmppartpair = partpair->GetSwapedPair();
-                   tmptrackpair = trackpair->GetSwapedPair();
-                   
+                 {
+                   tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
+                   tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
                  }
               }
             else
              {//meets criteria of the pair cut
-               tmptrackpair = trackpair;
-               tmppartpair = partpair;
+              tmptrackpair = trackpair;
+              tmppartpair = partpair;
              }
-            UInt_t ii;
+
             for(ii = 0;ii<fNParticleFunctions;ii++)
-                   fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
-                
+              fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
+
             for(ii = 0;ii<fNTrackFunctions;ii++)
-                   fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
-                 
+              fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
+
             for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
-                   fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
-           }
-       }
-    }
+              fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
+          }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
 
-  /***************************************/
-  /***** 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);
+      }
+    //end of loop over particles from first event
+   }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+  delete fPartBuffer->Push(partEvent1);
+  delete fTrackBuffer->Push(trackEvent1);
+ //end of loop over events  
+  return 0;
+}
+/*************************************************************************************/ 
 
-           track1= trackEvent->GetParticle(j);
-//         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
-           Int_t NNN;
+Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim)
+{
+  //Does analysis of simulated data
+  AliVAODParticle * part1, * part2;
   
-           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%100) == 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 
-        {
-          tmppartpair = partpair->GetSwapedPair();
-          tmptrackpair = trackpair->GetSwapedPair();
-        }
-                  }
-                else
-                 {//meets criteria of the pair cut
-                  tmptrackpair = trackpair;
-                  tmppartpair = partpair;
-                 }
-                UInt_t ii;
-                for(ii = 0;ii<fNParticleFunctions;ii++)
-                  fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
-                 
-                for(ii = 0;ii<fNTrackFunctions;ii++)
-                  fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
-                 
-                for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
-                  fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
-           
-               
-              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
-            }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
-       }
-    
-    } 
+  AliAOD* partEvent = aodsim;
+  AliAOD* partEvent1 = new AliAOD();
 
-  /***************************************/
+  AliAOD* partEvent2;
   
-   
-} 
-/*************************************************************************************/
+  AliHBTPair ppair;
 
-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;
+  AliHBTPair* partpair = &ppair;
+  AliHBTPair * tmppartpair;
   
-  AliHBTEvent * trackEvent;
-  AliHBTEvent * trackEvent2;
+  register UInt_t ii;
   
-//  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
+  if ( !partEvent )
+   {
+     Error("ProcessSim","Can not get event");
+     return 1;
+   }
+
+
+  for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+   {
+   /***************************************/
+   /******   Looping same events   ********/
+   /******   filling numerators    ********/
+   /***************************************/
+   if ( (j%fDisplayMixingInfo) == 0)
+      Info("ProcessSim",
+           "Mixing particle %d with particles from the same event",j);
+
+   part1= partEvent->GetParticle(j);
+
+   Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(part1);
+
+   if (fBufferSize != 0) 
+     if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(part1) == kFALSE ) )
+      {
+        //accepted by any cut
+        // we have to copy because reader keeps only one event
+
+        partEvent1->AddParticle(part1);
+      }
+
+   if (firstcut) continue;
+
+   for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+     fParticleMonitorFunctions[ii]->Process(part1);
+
+   if ( fNParticleFunctions == 0 ) continue;
+
+   for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+    {
+      part2= partEvent->GetParticle(k);
+      if (part1->GetUID() == part2->GetUID()) continue;
+      partpair->SetParticles(part1,part2);
+
+         if(fPairCut->Rejected(partpair)) //check pair cut
+          { //do not meets crietria of the 
+            if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) ) continue;
+            else tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
+          }
+         else
+          {
+            tmppartpair = partpair;
+          }
+
+      for(ii = 0;ii<fNParticleFunctions;ii++)
+             fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
+
+     //end of 2nd loop over particles from the same event  
+    }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+
+   /***************************************/
+   /***** Filling denominators    *********/
+   /***************************************/
+   if (fBufferSize == 0) continue;
+
+   fPartBuffer->ResetIter();
+       Int_t m = 0;
+       while (( partEvent2 = fPartBuffer->Next() ))
+        {
+          m++;
+          if ( (j%fDisplayMixingInfo) == 0)
+             Info("ProcessSim",
+                  "Mixing particle %d from current event with particles from event %d",j,-m);
+          for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
+            {
+
+              part2= partEvent2->GetParticle(l);
+              partpair->SetParticles(part1,part2);
+
+              if( fPairCut->Rejected(partpair) ) //check pair cut
+                { //do not meets crietria of the 
+                  if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) )
+                    continue;
+                  else 
+                   {
+                     tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
+                   }
+                }
+              else
+               {//meets criteria of the pair cut
+                tmppartpair = partpair;
+               }
+
+              for(ii = 0;ii<fNParticleFunctions;ii++)
+                fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
+
+           }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
+        }
+     }
+
+       
+  delete fPartBuffer->Push(partEvent1);
+ //end of loop over events  
+  return 0;
+}
+/*************************************************************************************/ 
+Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/)
+{
+  //Does analysis of reconstructed data
+  AliVAODParticle * track1, * track2;
   
-  Int_t Nev = fReader->GetNumberOfTrackEvents();
+  AliAOD* trackEvent = aodrec;
+  AliAOD* trackEvent1 = new AliAOD();
+
+  AliAOD* trackEvent2;
   
-  /***************************************/
-  /******   Looping same events   ********/
-  /******   filling numerators    ********/
-  /***************************************/
-  AliHBTPair * trackpair = new AliHBTPair();
-  AliHBTPair * tmptrackpair; //temporary pointer 
+  AliHBTPair tpair;
+
+  AliHBTPair* trackpair = &tpair;
+  AliHBTPair * tmptrackpair;
   
-  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%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+  register UInt_t ii;
 
-         track1= trackEvent->GetParticle(j);   
-        
-         for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
+  if ( !trackEvent )
+   {
+     Error("ProcessRecAndSim","Can not get event");
+     return 1;
+   }
+
+
+  for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
+   {
+     /***************************************/
+     /******   Looping same events   ********/
+     /******   filling numerators    ********/
+     /***************************************/
+     if ( (j%fDisplayMixingInfo) == 0)
+        Info("ProcessRec",
+             "Mixing Particle %d with Particles from the same event",j);
+
+     track1= trackEvent->GetParticle(j);
+
+     Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(track1);
+
+     if (fBufferSize != 0) 
+       if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(track1) == kFALSE ) )
+        {
+          //accepted by any cut
+          // we have to copy because reader keeps only one event
+
+          trackEvent1->AddParticle(track1);
+        }
+
+     if (firstcut) continue;
+
+     for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+       fParticleMonitorFunctions[ii]->Process(track1);
+
+     if ( fNTrackFunctions == 0 ) continue;
+
+     for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
+      {
+        track2= trackEvent->GetParticle(k);
+        if (track1->GetUID() == track2->GetUID()) continue;
+        trackpair->SetParticles(track1,track2);
+
+           if(fPairCut->Rejected(trackpair)) //check pair cut
+            { //do not meets crietria of the 
+              if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) ) continue;
+              else tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
+            }
+           else
+            {
+              tmptrackpair = trackpair;
+            }
+
+        for(ii = 0;ii<fNTrackFunctions;ii++)
+               fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
+
+       //end of 2nd loop over Particles from the same event  
+      }//for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
+
+     /***************************************/
+     /***** Filling denominators    *********/
+     /***************************************/
+     if (fBufferSize == 0) continue;
+
+     fTrackBuffer->ResetIter();
+     Int_t m = 0;
+     while (( trackEvent2 = fTrackBuffer->Next() ))
+      {
+        m++;
+        if ( (j%fDisplayMixingInfo) == 0)
+           Info("ProcessRec",
+                "Mixing Particle %d from current event with Particles from event %d",j,-m);
+        for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++)   //  ... on all Particles
           {
-            track2= trackEvent->GetParticle(k);        
+
+            track2= trackEvent2->GetParticle(l);
             trackpair->SetParticles(track1,track2);
-            if(fPairCut->Pass(trackpair)) //check pair cut
+
+            if( fPairCut->Rejected(trackpair) ) //check pair cut
               { //do not meets crietria of the 
-                if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
-                else tmptrackpair = trackpair->GetSwapedPair();
+                if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) )
+                  continue;
+                else 
+                 {
+                   tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
+                 }
               }
             else
-              {
-                tmptrackpair = trackpair;
-              }
-            UInt_t ii;
-                
+             {//meets criteria of the pair cut
+              tmptrackpair = trackpair;
+             }
+
             for(ii = 0;ii<fNTrackFunctions;ii++)
-                   fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
-                 
-           
-           }
+              fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
+
+         }//for(Int_t l = 0; l<N2;l++)   //  ... on all Particles
+      }
+   }
+  delete fTrackBuffer->Push(trackEvent1);
+ //end of loop over events  
+  return 0;
+}
+/*************************************************************************************/ 
+     
+Int_t AliHBTAnalysis::ProcessRecAndSimNonId(AliAOD* aodrec, AliAOD* aodsim)
+{
+//Analyzes both reconstructed and simulated data
+  if (aodrec == 0x0) 
+   {
+     Error("ProcessTracksAndParticlesNonIdentAnal","Reconstructed event is NULL");
+     return 1;
+   }
+
+  if (aodsim == 0x0) 
+   {
+     Error("ProcessTracksAndParticlesNonIdentAnal","Simulated event is NULL");
+     return 1;
+   }
+
+  if ( aodrec->GetNumberOfParticles() != aodsim->GetNumberOfParticles() )
+   {
+     Error("ProcessTracksAndParticlesNonIdentAnal",
+           "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
+           aodsim->GetNumberOfParticles() , aodrec->GetNumberOfParticles());
+     return 2;
+   }
+
+  AliVAODParticle * part1, * part2;
+  AliVAODParticle * track1, * track2;
+
+  static AliAOD aodrec1;
+  static AliAOD aodsim1;
+
+  AliAOD * trackEvent1=&aodrec1,*partEvent1=&aodsim1;//Particle that passes first particle cut, this event
+  trackEvent1->Reset();
+  partEvent1->Reset();
+  AliAOD * trackEvent2=0x0,*partEvent2=0x0;//Particle that passes second particle cut, this event
+  AliAOD * trackEvent3=0x0,*partEvent3=0x0;//Particle that passes second particle cut, events from buffer
+
+  AliAOD* rawtrackEvent = aodrec;//this we get from Reader
+  AliAOD* rawpartEvent = aodsim;//this we get from Reader
+
+  static AliHBTPair tpair;
+  static AliHBTPair ppair;
+  
+  AliHBTPair* trackpair = &tpair;
+  AliHBTPair* partpair = &ppair;
+
+  register UInt_t ii;
+  
+  /********************************/
+  /*      Filtering out           */
+  /********************************/
+  if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
+   {
+     partEvent2  = new AliAOD();
+     trackEvent2 = new AliAOD();
+   }
+
+  FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
+
+  for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
+   {
+     if ( (j%fDisplayMixingInfo) == 0) 
+        Info("ProcessTracksAndParticlesNonIdentAnal",
+             "Mixing particle %d from current event with particles from current event",j);
+
+     part1= partEvent1->GetParticle(j);
+     track1= trackEvent1->GetParticle(j);
+
+
+     for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+       fParticleMonitorFunctions[ii]->Process(part1);
+     for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+       fTrackMonitorFunctions[ii]->Process(track1);
+     for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
+       fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
+
+     if (fNoCorrfctns) continue;
+
+     /***************************************/
+     /******   filling numerators    ********/
+     /****** (particles from event2) ********/
+     /***************************************/
+
+     for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups 
+      {
+        part2= partEvent2->GetParticle(k);
+        if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
+        partpair->SetParticles(part1,part2);
+
+        track2= trackEvent2->GetParticle(k);
+        trackpair->SetParticles(track1,track2);
+
+        if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
+         { //do not meets crietria of the pair cut
+          continue; 
+         }
+        else
+         {//meets criteria of the pair cut
+          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 diff histogram *********/
-  /***************************************/
-  for (Int_t i = 0;i<Nev;i++)   //In each event ....
+   if ( fBufferSize == 0) continue;//do not mix diff histograms
+   /***************************************/
+   /***** Filling denominators    *********/
+   /***************************************/
+   fPartBuffer->ResetIter();
+   fTrackBuffer->ResetIter();
+
+   Int_t nmonitor = 0;
+
+   while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
     {
-      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);
+      trackEvent3 = fTrackBuffer->Next();
+
+      if ( (j%fDisplayMixingInfo) == 0) 
+        Info("ProcessTracksAndParticlesNonIdentAnal",
+             "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
+
+      for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
+        {
+          part2= partEvent3->GetParticle(k);
+          partpair->SetParticles(part1,part2);
+
+          track2= trackEvent3->GetParticle(k);
+          trackpair->SetParticles(track1,track2);
+
+          if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
+           { //do not meets crietria of the pair cut
+            continue; 
+           }
+          else
+           {//meets criteria of the pair cut
+            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 particles event2
+       }//while event2
+   }//for over particles in event1
+
+ delete fPartBuffer->Push(partEvent2);
+ delete fTrackBuffer->Push(trackEvent2);
  
-//         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;
+ return 0;
+}
+/*************************************************************************************/ 
+Int_t AliHBTAnalysis::ProcessSimNonId(AliAOD* /*aodrec*/, AliAOD* aodsim)
+{
+//does analysis of simulated (MC) data in non-identical mode
+//i.e. when particles selected by first part. cut are a disjunctive set than particles
+//passed by the second part. cut
+ if (aodsim == 0x0) 
+  {
+    return 1;
+  }
  
-           for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
-            {
-             
-             trackEvent2 = fReader->GetTrackEvent(k);
-             if (!trackEvent2) continue;
-             
-             if ( (j%100) == 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 tmptrackpair = trackpair->GetSwapedPair();
-                  }
-                else
-                  {
-                    tmptrackpair = trackpair;
-                  }
-                UInt_t ii;
-                for(ii = 0;ii<fNTrackFunctions;ii++)
-                  fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
-               
-              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
-            }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
-       }
-    
-    } 
+  AliVAODParticle * part1, * part2;
+
+  static AliAOD aodsim1;
+
+  AliAOD* partEvent1=&aodsim1;//Particle that passes first particle cut, this event
+  partEvent1->Reset();
+  AliAOD* partEvent2=0x0;//Particle that passes second particle cut, this event
+  AliAOD* partEvent3=0x0;//Particle that passes second particle cut, events from buffer
+
+  AliAOD* rawpartEvent = aodsim;//this we get from Reader
 
-  /***************************************/
+  static AliHBTPair ppair;
   
+  AliHBTPair* partpair = &ppair;
 
-}
+  register UInt_t ii;
+  
+  /********************************/
+  /*      Filtering out           */
+  /********************************/
+  if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
+   {
+     partEvent2  = new AliAOD();
+   }
 
-/*************************************************************************************/
+  FilterOut(partEvent1, partEvent2, rawpartEvent);
 
-void AliHBTAnalysis::ProcessParticles()
+  for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
+   {
+     if ( (j%fDisplayMixingInfo) == 0) 
+        Info("ProcessParticlesNonIdentAnal",
+             "Mixing particle %d from current event with particles from current event",j);
+
+     part1= partEvent1->GetParticle(j);
+
+
+     for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+       fParticleMonitorFunctions[ii]->Process(part1);
+
+     if (fNParticleFunctions == 0) continue;
+
+     /***************************************/
+     /******   filling numerators    ********/
+     /****** (particles from event2) ********/
+     /***************************************/
+
+     for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups 
+      {
+        part2= partEvent2->GetParticle(k);
+        if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
+        partpair->SetParticles(part1,part2);
+
+
+        if(fPairCut->PassPairProp(partpair) ) //check pair cut
+          { //do not meets crietria of the pair cut
+              continue; 
+          }
+        else
+         {//meets criteria of the pair cut
+          for(ii = 0;ii<fNParticleFunctions;ii++)
+              fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
+         }
+       }
+
+ if ( fBufferSize == 0) continue;//do not mix diff histograms
+ /***************************************/
+ /***** Filling denominators    *********/
+ /***************************************/
+ fPartBuffer->ResetIter();
+
+ Int_t nmonitor = 0;
+
+ while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
+  {
+
+    if ( (j%fDisplayMixingInfo) == 0) 
+      Info("ProcessParticlesNonIdentAnal",
+           "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
+
+    for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
+      {
+        part2= partEvent3->GetParticle(k);
+        partpair->SetParticles(part1,part2);
+
+
+        if(fPairCut->PassPairProp(partpair) ) //check pair cut
+          { //do not meets crietria of the pair cut
+              continue; 
+          }
+        else
+         {//meets criteria of the pair cut
+          for(ii = 0;ii<fNParticleFunctions;ii++)
+           {
+             fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
+           }
+         }
+       }// for particles event2
+     }//while event2
+   }//for over particles in event1
+
+ delete fPartBuffer->Push(partEvent2);
+ return 0;
+}
+/*************************************************************************************/ 
+Int_t AliHBTAnalysis::ProcessRecNonId(AliAOD* aodrec, AliAOD* /*aodsim*/)
 {
-  //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;
+//Analyzes both reconstructed and simulated data
+ if (aodrec == 0x0) 
+  {
+    return 1;
+  }
+  AliVAODParticle * track1, * track2;
+
+  static AliAOD aodrec1;
+  AliAOD * trackEvent1=&aodrec1;//Particle that passes first particle cut, this event
+  trackEvent1->Reset();
+  AliAOD * trackEvent2=0x0;//Particle that passes second particle cut, this event
+  AliAOD * trackEvent3=0x0;//Particle that passes second particle cut, events from buffer
+  AliAOD* rawtrackEvent = aodrec;//this we get from Reader
 
-  AliHBTPair * partpair = new AliHBTPair();
-  AliHBTPair * tmppartpair; //temporary pointer to the pair
+  static AliHBTPair tpair;
   
-//  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
+  AliHBTPair* trackpair = &tpair;
+
+  register UInt_t ii;
   
-  Int_t Nev = fReader->GetNumberOfPartEvents();
   
-  /***************************************/
-  /******   Looping same events   ********/
-  /******   filling numerators    ********/
-  /***************************************/
-  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%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+  /********************************/
+  /*      Filtering out           */
+  /********************************/
+  if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
+   {
+     trackEvent2 = new AliAOD();
+   }
 
-         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
-                   tmppartpair = partpair->GetSwapedPair();
-                 }
-              }
-            else
-              {
-                tmppartpair = partpair;
-              }
+  FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
 
-            UInt_t ii;
-                
-            for(ii = 0;ii<fNParticleFunctions;ii++)
-                   fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
-           }
+  for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
+   {
+     if ( (j%fDisplayMixingInfo) == 0) 
+        Info("ProcessTracksNonIdentAnal",
+             "Mixing particle %d from current event with particles from current event",j);
+
+     track1= trackEvent1->GetParticle(j);
+
+
+     for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+       fTrackMonitorFunctions[ii]->Process(track1);
+
+     if (fNTrackFunctions == 0x0) continue;
+
+     /***************************************/
+     /******   filling numerators    ********/
+     /****** (particles from event2) ********/
+     /***************************************/
+
+     for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups 
+      {
+        track2= trackEvent2->GetParticle(k);
+        if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
+        trackpair->SetParticles(track1,track2);
+
+
+        if( fPairCut->PassPairProp(trackpair)) //check pair cut
+         { //do not meets crietria of the pair cut
+           continue; 
+         }
+        else
+         {//meets criteria of the pair cut
+           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 ....
-    {
-      partEvent= fReader->GetParticleEvent(i);
-      if (!partEvent) continue;
+ if ( fBufferSize == 0) continue;//do not mix diff histograms
+ /***************************************/
+ /***** Filling denominators    *********/
+ /***************************************/
+ fTrackBuffer->ResetIter();
 
-//      N=0;
-      
-      for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
-       {
-//         if (N>MAXCOMB) break;
-           
-           part1= partEvent->GetParticle(j);
+ Int_t nmonitor = 0;
+
+ while ( (trackEvent3 = fTrackBuffer->Next() ) != 0x0)
+  {
+    if ( (j%fDisplayMixingInfo) == 0) 
+      Info("ProcessTracksNonIdentAnal",
+           "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
+
+    for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
+      {
+        track2= trackEvent3->GetParticle(k);
+        trackpair->SetParticles(track1,track2);
+
+        if( fPairCut->PassPairProp(trackpair)) //check pair cut
+         { //do not meets crietria of the pair cut
+           continue; 
+         }
+        else
+         {//meets criteria of the pair cut
+           for(ii = 0;ii<fNTrackFunctions;ii++)
+               fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
+         }
+       }// for particles event2
+     }//while event2
+   }//for over particles in event1
+
+ delete fTrackBuffer->Push(trackEvent2);
  
-//         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;
+ return 0;
+}
+/*************************************************************************************/ 
+
+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("Process","The reader is not set");
+   return;
+  }
  
-           for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
-            {
-             
-             partEvent2= fReader->GetParticleEvent(k);
-             if (!partEvent2) continue;
-             
-             if ( (j%100) == 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 tmppartpair = partpair->GetSwapedPair();
-                  }
-                else
-                  {
-                    tmppartpair = partpair;
-                  }
-                UInt_t ii;
-                for(ii = 0;ii<fNParticleFunctions;ii++)
-                  fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
-               
-              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
-            }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
-       }
-    
-    } 
+ const char *oT = strstr(option,"Tracks");
+ const char *oP = strstr(option,"Particles");
+ Bool_t nonid = IsNonIdentAnalysis();
+ Init();
+ if(oT && oP)
+  { 
+    if (nonid) ProcessTracksAndParticlesNonIdentAnal();
+    else ProcessTracksAndParticles();
+    return;
+  }
+
+ if(oT)
+  {
+    if (nonid) ProcessTracksNonIdentAnal();
+    else ProcessTracks();
+    return;
+  }
+ if(oP)
+  {
+    if (nonid) ProcessParticlesNonIdentAnal();
+    else ProcessParticles();
+    return;
+  }
+}
+/*************************************************************************************/ 
 
-  /***************************************/
+void AliHBTAnalysis::ProcessTracksAndParticles()
+{
+//Makes analysis for both tracks and particles
+//mainly for resolution study and analysies with weighting algirithms
+//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+  
+// cut on particles only -- why?
+// - PID: when we make resolution analysis we want to take only tracks with correct PID
+// We need cut on tracks because there are data characteristic to 
   
+  AliAOD * trackEvent, *partEvent;
+
+  fReader->Rewind();
+  while (fReader->Next() == kFALSE)
+    {
+      partEvent = fReader->GetEventSim();
+      trackEvent = fReader->GetEventRec();
+      ProcessRecAndSim(trackEvent,partEvent);
+
+    }//while (fReader->Next() == kFALSE)
+    
+} 
+/*************************************************************************************/
 
+void AliHBTAnalysis::ProcessTracks()
+{
+//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+  AliAOD * trackEvent;
+  fReader->Rewind();
+  while (fReader->Next() == kFALSE)
+    {
+      trackEvent = fReader->GetEventRec();
+      ProcessRec(trackEvent,0x0);
+    }//while (fReader->Next() == kFALSE)
 }
 
 /*************************************************************************************/
 
+void AliHBTAnalysis::ProcessParticles()
+{
+//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+  AliAOD * partEvent;
+  fReader->Rewind();
+  while (fReader->Next() == kFALSE)
+    {
+      partEvent = fReader->GetEventSim();
+      ProcessSim(0x0,partEvent);
+    }//while (fReader->Next() == kFALSE)
+}
+/*************************************************************************************/
 
 void AliHBTAnalysis::WriteFunctions()
 {
+//Calls Write for all defined functions in analysis
+//== writes all results
+  TFile* oututfile = 0x0;
+  if (fOutputFileName)
+   {
+     oututfile = TFile::Open(*fOutputFileName,"update");
+   }
   UInt_t ii;
   for(ii = 0;ii<fNParticleFunctions;ii++)
-    fParticleFunctions[ii]->Write();
+   {
+     if (AliVAODParticle::GetDebug()>5)
+      {
+        Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
+        Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
+      }
+     fParticleFunctions[ii]->Write();
+   }
                  
   for(ii = 0;ii<fNTrackFunctions;ii++)
-    fTrackFunctions[ii]->Write();
+   {
+     if (AliVAODParticle::GetDebug()>5)
+      {
+        Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
+        Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
+      }
+     fTrackFunctions[ii]->Write();
+   }
                  
   for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
-   fParticleAndTrackFunctions[ii]->Write();
+   {
+     if (AliVAODParticle::GetDebug()>5)
+      {
+        Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
+        Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
+      } 
+     fParticleAndTrackFunctions[ii]->Write();
+   }
+
+  for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
+   {
+     if (AliVAODParticle::GetDebug()>5)
+      {
+        Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
+        Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
+      }
+     fParticleMonitorFunctions[ii]->Write();
+   }
+
+  for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
+   {
+     if (AliVAODParticle::GetDebug()>5)
+      {
+        Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
+        Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
+      }
+     fTrackMonitorFunctions[ii]->Write();
+   }
+
+  for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
+   {
+     if (AliVAODParticle::GetDebug()>5)
+      {
+       Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
+       Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
+      }
+     fParticleAndTrackMonitorFunctions[ii]->Write();
+   }
+  delete oututfile; 
 }
 /*************************************************************************************/
 
-void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
+void AliHBTAnalysis::SetOutputFileName(const char* fname)
 {
+  //Sets fiele name where to dump results, 
+  //if not specified reults are written to gDirectory
+  if (fname == 0x0)
+   {
+     delete fOutputFileName;
+     fOutputFileName = 0x0;
+     return;
+   }
+  if ( strcmp(fname,"") == 0 ) 
+   {
+     delete fOutputFileName;
+     fOutputFileName = 0x0;
+     return;
+   }
+  if (fOutputFileName == 0x0) fOutputFileName = new TString(fname);
+  else *fOutputFileName = fname;
+}
+/*************************************************************************************/
+
+void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
+{
+//Sets the global cut
   if (cut == 0x0)
    {
      Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
    }
   delete fPairCut;
-  fPairCut = (AliHBTPairCut*)cut->Clone();
+  fPairCut = (AliAODPairCut*)cut->Clone();
 }
 
 /*************************************************************************************/
 
 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
 {
+//Adds track function
   if (f == 0x0) return;
   if (fNTrackFunctions == fgkFctnArraySize)
    {
@@ -583,6 +1393,7 @@ void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
 
 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
 {
+//adds particle function
   if (f == 0x0) return;
   
   if (fNParticleFunctions == fgkFctnArraySize)
@@ -591,11 +1402,12 @@ void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
    }
   fParticleFunctions[fNParticleFunctions] = f;
   fNParticleFunctions++;
-  
-  
 }
+/*************************************************************************************/ 
+
 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
 {
+//add resolution function
   if (f == 0x0) return;
   if (fNParticleAndTrackFunctions == fgkFctnArraySize)
    {
@@ -604,11 +1416,50 @@ void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
   fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
   fNParticleAndTrackFunctions++;
 }
+/*************************************************************************************/ 
+
+void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
+{
+//add particle monitoring function
+  if (f == 0x0) return;
+
+  if (fNParticleMonitorFunctions == fgkFctnArraySize)
+    {
+      Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
+   }
+  fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
+  fNParticleMonitorFunctions++;
+}
+/*************************************************************************************/ 
 
+void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
+{
+//add track monitoring function
+  if (f == 0x0) return;
 
+  if (fNTrackMonitorFunctions == fgkFctnArraySize)
+   {
+    Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
+   }
+  fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
+  fNTrackMonitorFunctions++;
+}
 /*************************************************************************************/ 
 
+void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
+{
+//add resolution monitoring function
+  if (f == 0x0) return;
+  if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
+    {
+      Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
+    }
+  fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
+  fNParticleAndTrackMonitorFunctions++;
+}
+
 
+/*************************************************************************************/ 
 /*************************************************************************************/  
 
 Bool_t AliHBTAnalysis::RunCoherencyCheck()
@@ -616,43 +1467,46 @@ 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;
-  }
+ Info("RunCoherencyCheck","Checking HBT Runs Coherency");
+
+//When we use non-buffering reader this is a big waste of time -> We need to read all data to check it
+//and reader is implemented safe in this case anyway
+// Info("RunCoherencyCheck","Number of events ...");
+// if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same  number of events
+//  {
+//    Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
+//  }
+// else
+//  { //if not the same - ERROR
+//   Error("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);
+ Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
  
- AliHBTEvent *partEvent;
- AliHBTEvent *trackEvent;
+ AliAOD *partEvent;
+ AliAOD *trackEvent;
  for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
   {
-    partEvent= fReader->GetParticleEvent(i); //gets the "ith" event 
-    trackEvent = fReader->GetTrackEvent(i);
+    partEvent= fReader->GetEventSim(i); //gets the "ith" event 
+    trackEvent = fReader->GetEventRec(i);
     
     if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
     if ( (partEvent == 0x0) || (partEvent == 0x0) )
      {
-       Error("AliHBTAnalysis::RunCoherencyCheck()",
+       Error("RunCoherencyCheck",
              "One event is NULL and the other one not. Event Number %d",i);
        return kTRUE;    
      }
     
     if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
      {
-       Error("AliHBTAnalysis::RunCoherencyCheck()",
+       Error("RunCoherencyCheck",
              "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
              i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
        return kTRUE;
@@ -662,7 +1516,7 @@ Bool_t AliHBTAnalysis::RunCoherencyCheck()
       {
         if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
          {
-           Error("AliHBTAnalysis::RunCoherencyCheck()",
+           Error("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;
@@ -670,15 +1524,220 @@ Bool_t AliHBTAnalysis::RunCoherencyCheck()
          }
       }
   }
- cout<<"  Done"<<endl;
- cout<<"  Everything looks OK"<<endl;
+ Info("RunCoherencyCheck","  Done");
+ Info("RunCoherencyCheck","  Everything looks OK");
+*/ 
  return kFALSE;
 }
 
+/*************************************************************************************/  
+void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
+{
+//Performs analysis for both, tracks and particles
+  AliAOD* rawtrackEvent, * rawpartEvent;
+  fReader->Rewind();
+  Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
+  Info("ProcessTracksAndParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
+  Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
 
+  for (Int_t i = 0;;i++)//infinite loop
+    {
+      if (fReader->Next()) break; //end when no more events available
+      
+      rawpartEvent  = fReader->GetEventSim();
+      rawtrackEvent = fReader->GetEventRec();
+      
+      ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
+  }//end of loop over events (1)
+}
 /*************************************************************************************/  
  
+void AliHBTAnalysis::ProcessTracksNonIdentAnal()
+{
+//Process Tracks only with non identical mode
+  AliAOD * rawtrackEvent;
+  fReader->Rewind();
+  
+  Info("ProcessTracksNonIdentAnal","**************************************");
+  Info("ProcessTracksNonIdentAnal","*****  NON IDENT MODE ****************");
+  Info("ProcessTracksNonIdentAnal","**************************************");
+  
+  for (Int_t i = 0;;i++)//infinite loop
+    {
+      if (fReader->Next()) break; //end when no more events available
+      rawtrackEvent = fReader->GetEventRec();
+      ProcessRecNonId(rawtrackEvent,0x0);
+  }//end of loop over events (1)
+}
+/*************************************************************************************/  
+
+void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
+{
+//process paricles only with non identical mode
+  AliAOD * rawpartEvent = 0x0;
+  fReader->Rewind();
+  
+  Info("ProcessParticlesNonIdentAnal","**************************************");
+  Info("ProcessParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
+  Info("ProcessParticlesNonIdentAnal","**************************************");
+
+  for (Int_t i = 0;;i++)//infinite loop
+    {
+      if (fReader->Next()) break; //end when no more events available
+      
+      rawpartEvent  = fReader->GetEventSim();
+      ProcessSimNonId(0x0,rawpartEvent);
+  }//end of loop over events (1)
+}
 
 /*************************************************************************************/  
+void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
+                     AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
+{
+ //Puts particles accepted as a first particle by global cut in out1
+ //and as a second particle in out2
+
+  AliVAODParticle* part, *track;
+
+  outpart1->Reset();
+  outpart2->Reset();
+  outtrack1->Reset();
+  outtrack2->Reset();
+  
+  Bool_t in1, in2;
+  
+  for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
+   {
+     in1 = in2 = kTRUE;
+     part = inpart->GetParticle(i);
+     track = intrack->GetParticle(i);
+     
+     if ( ((this->*fkPass1)(part,track))  ) in1 = kFALSE; //if part  is rejected by cut1, in1 is false
+     if ( ((this->*fkPass2)(part,track))  ) in2 = kFALSE; //if part  is rejected by cut2, in2 is false
+     
+     if (gDebug)//to be removed in real analysis     
+     if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
+      {
+      //Particle accpted by both cuts
+       Error("FilterOut","Particle accepted by both cuts");
+       continue;
+      }
 
+     if (in1)
+      {
+        outpart1->AddParticle(part);
+        outtrack1->AddParticle(track);
+        continue;
+      }
+     
+     if (in2)
+      {
+        outpart2->AddParticle(part);
+        outtrack2->AddParticle(track);
+        continue;
+      }
+   }
+}
+/*************************************************************************************/  
+
+void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
+{
+ //Puts particles accepted as a first particle by global cut in out1
+ //and as a second particle in out2
+  AliVAODParticle* part;
+  
+  out1->Reset();
+  out2->Reset();
+  
+  AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
+  AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
+  
+  Bool_t in1, in2;
+  
+  for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
+   {
+     in1 = in2 = kTRUE;
+     part = in->GetParticle(i);
+     
+     if ( cut1->Rejected(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
+     if ( cut2->Rejected(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
+     
+     if (gDebug)//to be removed in real analysis     
+     if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
+      {
+      //Particle accpted by both cuts
+       Error("FilterOut","Particle accepted by both cuts");
+       continue;
+      }
+
+     if (in1)
+      { 
+        out1->AddParticle(part);
+        continue;
+      }
+     
+     if (in2)
+      {
+        out2->AddParticle(part);
+        continue;
+      }
+   }
+}
+/*************************************************************************************/ 
+
+Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
+{
+ //checks if it is possible to use special analysis for non identical particles
+ //it means - in global pair cut first particle id is different than second one 
+ //and both are different from 0
+ //in the future is possible to perform more sophisticated check 
+ //if cuts have excluding requirements
+ if (fPairCut->IsEmpty()) 
+   return kFALSE;
+ if (fPairCut->GetFirstPartCut()->IsEmpty()) 
+   return kFALSE;
+
+ if (fPairCut->GetSecondPartCut()->IsEmpty()) 
+   return kFALSE;
+ Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
+ Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
+
+ if ( (id1==0) || (id2==0) || (id1==id2) ) 
+   return kFALSE;
+ return kTRUE;
+}
+/*************************************************************************************/ 
+
+void AliHBTAnalysis::SetApparentVertex(Double_t x, Double_t y, Double_t z)
+{ 
+  //Sets apparent vertex
+  // All events have to be moved to the same vertex position in order to
+  // to be able to comare any space positions (f.g. anti-merging)
+  // This method defines this position
+  
+  fVertexX = x;
+  fVertexY = y;
+  fVertexZ = z;
+}
+/*************************************************************************************/ 
+
+void AliHBTAnalysis::PressAnyKey()
+{ 
+ //small utility function that helps to make comfortable macros
+  char c;
+  int nread = -1;
+  fcntl(0,  F_SETFL, O_NONBLOCK);
+  ::Info("","Press Any Key to continue ...");
+  while (nread<1) 
+   {
+     nread = read(0, &c, 1);
+     gSystem->ProcessEvents();
+   }
+}