AliAOD and AODParticle (T.Kuhr) - Readers, AODStdParticle and Cuts (P.Skowronski...
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Apr 2004 06:12:30 +0000 (06:12 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Apr 2004 06:12:30 +0000 (06:12 +0000)
32 files changed:
ANALYSIS/ANALYSISLinkDef.h
ANALYSIS/AliAOD.cxx [new file with mode: 0644]
ANALYSIS/AliAOD.h [new file with mode: 0644]
ANALYSIS/AliAODParticle.cxx [new file with mode: 0644]
ANALYSIS/AliAODParticle.h [new file with mode: 0644]
ANALYSIS/AliAODParticleCut.cxx [new file with mode: 0644]
ANALYSIS/AliAODParticleCut.h [new file with mode: 0644]
ANALYSIS/AliAODRun.cxx [new file with mode: 0644]
ANALYSIS/AliAODRun.h [new file with mode: 0644]
ANALYSIS/AliAODStdParticle.cxx [new file with mode: 0644]
ANALYSIS/AliAODStdParticle.h [new file with mode: 0644]
ANALYSIS/AliAnalysis.cxx
ANALYSIS/AliAnalysis.h
ANALYSIS/AliBaseEventCut.h
ANALYSIS/AliClusterMap.cxx [new file with mode: 0644]
ANALYSIS/AliClusterMap.h [new file with mode: 0644]
ANALYSIS/AliEventCut.cxx
ANALYSIS/AliEventCut.h
ANALYSIS/AliFlowAnalysis.cxx
ANALYSIS/AliFlowAnalysis.h
ANALYSIS/AliParticle.cxx [new file with mode: 0644]
ANALYSIS/AliParticle.h [new file with mode: 0644]
ANALYSIS/AliReader.cxx [new file with mode: 0644]
ANALYSIS/AliReader.h [new file with mode: 0644]
ANALYSIS/AliReaderESD.cxx [new file with mode: 0644]
ANALYSIS/AliReaderESD.h [new file with mode: 0644]
ANALYSIS/AliRunAnalysis.cxx
ANALYSIS/AliRunAnalysis.h
ANALYSIS/AliTrackPoints.cxx [new file with mode: 0644]
ANALYSIS/AliTrackPoints.h [new file with mode: 0644]
ANALYSIS/analysis.C [new file with mode: 0644]
ANALYSIS/libANALYSIS.pkg

index 0f37990..8d59d57 100644 (file)
@@ -1,6 +1,6 @@
 #ifdef __CINT__
 
-#pragma link off all globals;
+#pragma link off all glols;
 #pragma link off all classes;
 #pragma link off all functions;
 
@@ -9,9 +9,43 @@
 #pragma link C++ class AliRunAnalysis+;
 #pragma link C++ class AliAnalysis+;
 
+#pragma link C++ class AliAOD+;
+#pragma link C++ class AliAODParticle+;
+#pragma link C++ class AliAODStdParticle+;
+
+#pragma link C++ class AliAODRun+;
+
+#pragma link C++ class AliTrackPoints+;
+#pragma link C++ class AliClusterMap+;
+
+#pragma link C++ class AliReader+;
+#pragma link C++ class AliReaderESD+;
+
+
 #pragma link C++ class AliFlowAnalysis+;
 
 #pragma link C++ class AliEventCut+;
 
+#pragma link C++ class AliAODParticleCut-;
+#pragma link C++ class AliAODEmptyParticleCut-;
+#pragma link C++ class AliAODBaseCut+;
+
+#pragma link C++ class AliAODMomentumCut+;
+#pragma link C++ class AliAODPtCut+;
+#pragma link C++ class AliAODEnergyCut+;
+#pragma link C++ class AliAODRapidityCut+;
+#pragma link C++ class AliAODPseudoRapidityCut+;
+#pragma link C++ class AliAODPxCut+;
+#pragma link C++ class AliAODPyCut+;
+#pragma link C++ class AliAODPzCut+;
+#pragma link C++ class AliAODPhiCut+;
+#pragma link C++ class AliAODThetaCut+;
+#pragma link C++ class AliAODVxCut+;
+#pragma link C++ class AliAODVyCut+;
+#pragma link C++ class AliAODVzCut+;
+#pragma link C++ class AliAODPIDCut+;
+#pragma link C++ class AliAODLogicalOperCut-;
+#pragma link C++ class AliAODAndCut+;
+#pragma link C++ class AliAODOrCut+;
 
 #endif
diff --git a/ANALYSIS/AliAOD.cxx b/ANALYSIS/AliAOD.cxx
new file mode 100644 (file)
index 0000000..7e3de10
--- /dev/null
@@ -0,0 +1,76 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+/////////////////////////////////////////////////////////////
+//
+// base class for AOD containers
+//
+/////////////////////////////////////////////////////////////
+
+#include <TParticle.h>
+#include "AliAOD.h"
+
+#include "AliAODStdParticle.h"
+
+ClassImp(AliAOD)
+
+/**************************************************************************/
+
+void  AliAOD::AddParticle(TParticle* part, Int_t idx)
+{
+  //Adds TParticle to event
+  if (part == 0x0) 
+   {
+     Error("AddParticle(TParticle*,Int_t)","pointer to particle is NULL");
+     return;
+   }
+  AddParticle( new AliAODStdParticle(*part,idx) );
+}
+/**************************************************************************/
+
+void  AliAOD::AddParticle(Int_t pdg, Int_t idx,
+                          Double_t px, Double_t py, Double_t pz, Double_t etot,
+                          Double_t vx, Double_t vy, Double_t vz, Double_t time)
+{
+  //adds particle to event
+  AddParticle(new  AliAODStdParticle(pdg,idx,px,py,pz,etot,vx,vy,vz,time));
+}
+/**************************************************************************/
+
+void AliAOD::SwapParticles(Int_t i, Int_t j)
+{
+//swaps particles positions; used by AliHBTEvent::Blend
+  if ( (i<0) || (i>=GetNumberOfParticles()) ) return;
+  if ( (j<0) || (j>=GetNumberOfParticles()) ) return;
+
+  AliAODParticle* tmp = (AliAODParticle*)fParticles.At(i);
+  fParticles.AddAt(fParticles.At(j),i);
+  fParticles.AddAt(tmp,j);
+}
+/**************************************************************************/
+
+void  AliAOD::Reset()
+{
+  //deletes all particles from the event
+   for(Int_t i =0; i<GetNumberOfParticles(); i++)
+    {
+      for (Int_t j = i+1; j<GetNumberOfParticles(); j++)
+        if ( fParticles.At(j) == fParticles.At(i) ) fParticles.RemoveAt(j);
+      delete fParticles.RemoveAt(i);
+    }
+//   fRandomized = kFALSE;
+}
diff --git a/ANALYSIS/AliAOD.h b/ANALYSIS/AliAOD.h
new file mode 100644 (file)
index 0000000..0ab8c2f
--- /dev/null
@@ -0,0 +1,41 @@
+#ifndef ALIAOD_H
+#define ALIAOD_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+/////////////////////////////////////////////////////////////
+//
+// base class for AOD containers
+//
+/////////////////////////////////////////////////////////////
+
+#include <TObject.h>
+#include <TObjArray.h>
+#include "AliAODParticle.h"
+
+class TParticle;
+
+class AliAOD: public TObject {
+public:
+  AliAOD(){}
+  virtual ~AliAOD() { Reset(); }
+
+  virtual TObjArray*       GetParticles() {return &fParticles;};
+  virtual Int_t            GetNumberOfParticles() const  {return fParticles.GetEntriesFast();}
+  virtual AliAODParticle*  GetParticle(Int_t index) const {return (AliAODParticle*) fParticles[index];}
+  virtual void             AddParticle(AliAODParticle* particle)  {fParticles.Add(particle);};
+  virtual void             AddParticle(TParticle* part, Int_t idx); //adds particle to the event
+  virtual void             AddParticle(Int_t pdg, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot,
+                                       Double_t vx, Double_t vy, Double_t vz, Double_t time);
+  
+  virtual void             Reset();
+  void                     SwapParticles(Int_t i, Int_t j);//swaps particles positions; used by AliReader::Blend
+private:
+  TObjArray                fParticles;   // array of AOD particles, AliAOD is owner of particles
+
+  ClassDef(AliAOD,1)  // base class for AOD containers
+};
+
+#endif
diff --git a/ANALYSIS/AliAODParticle.cxx b/ANALYSIS/AliAODParticle.cxx
new file mode 100644 (file)
index 0000000..36cf642
--- /dev/null
@@ -0,0 +1,39 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+/////////////////////////////////////////////////////////////
+//
+// base class for AOD particles
+//
+// Most of the methods defined in this base class (or even all of them)
+// should be redefined in derived classes. Some methods like Pt() or
+// Vx() have default implementation, but they are not optimised for
+// performance. It is likely that they can be (and should be) implemented
+// in the derived classes in a way which gives faster access.
+//
+// Algorithms and analysis classes should as much as possible use the
+// interface of the base class instead of special methods of derived
+// classes. This allows to run the code on all types of particles.
+//
+/////////////////////////////////////////////////////////////
+
+#include "AliAODParticle.h"
+
+Int_t AliAODParticle::fgDebug = 0;
+
+ClassImp(AliAODParticle)
+
diff --git a/ANALYSIS/AliAODParticle.h b/ANALYSIS/AliAODParticle.h
new file mode 100644 (file)
index 0000000..cc1d032
--- /dev/null
@@ -0,0 +1,81 @@
+#ifndef ALIAODPARTICLE_H
+#define ALIAODPARTICLE_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+/////////////////////////////////////////////////////////////
+//
+// base class for AOD particles
+//
+/////////////////////////////////////////////////////////////
+
+#include <TObject.h>
+#include <TLorentzVector.h>
+#include <TVector3.h>
+
+#include "AliAnalysis.h"
+
+class AliTrackPoints;
+class AliClusterMap;
+
+class AliAODParticle : public TObject {
+public:
+   AliAODParticle(){}
+   virtual ~AliAODParticle(){}
+  // kinematics
+  virtual TLorentzVector   FourMomentum() const = 0;
+  virtual TVector3         Momentum() const {return FourMomentum().Vect();};
+  virtual Double_t         Mass() const {return FourMomentum().M();};
+  virtual Double_t         E() const {return FourMomentum().E();};
+  virtual Double_t         P() const {return FourMomentum().P();};
+  virtual Double_t         Pt() const {return FourMomentum().Pt();};
+  virtual Double_t         Px() const {return FourMomentum().Px();};
+  virtual Double_t         Py() const {return FourMomentum().Py();};
+  virtual Double_t         Pz() const {return FourMomentum().Pz();};
+  virtual Double_t         Phi() const {return FourMomentum().Phi();};
+  virtual Double_t         Theta() const {return FourMomentum().Theta();};
+  virtual Double_t         Eta() const {return FourMomentum().Eta();};
+  virtual Double_t         Y() const {return FourMomentum().Rapidity();};
+
+  // PID
+  virtual Double_t         Charge() const = 0;
+  virtual Double_t         GetProbability(Int_t pdg) const = 0;
+  virtual Int_t            GetMostProbable() const = 0;
+  
+  virtual Int_t            GetPdgCode() const = 0;//We need to assume some PID (f.e. energy calculation) 
+                                                  //sotimes one track can apear in analysis twise (e.g. ones as pion ones as kaon)
+  
+  // vertices
+  virtual TVector3         ProductionVertex() const = 0;
+  virtual Double_t         Vx() const {return ProductionVertex().X();};
+  virtual Double_t         Vy() const {return ProductionVertex().Y();};
+  virtual Double_t         Vz() const {return ProductionVertex().Z();};
+  virtual AliAODParticle*  Mother() const {return NULL;};
+  virtual Bool_t           HasDecayVertex() const {return kFALSE;};
+  virtual TVector3         DecayVertex() const {return TVector3();};
+  virtual Int_t            NumberOfDaughters() const {return 0;};
+  virtual AliAODParticle*  Daughter(Int_t /*index*/) const {return NULL;};
+
+
+  // type information
+  virtual Bool_t           IsSimulated() {return kFALSE;};
+  virtual Bool_t           IsTrack() {return kFALSE;};
+  virtual Bool_t           IsCluster() {return kFALSE;};
+
+  //HBT specific 
+  virtual AliTrackPoints*  GetTrackPoints() const {return 0x0;}
+  virtual AliClusterMap*   GetClusterMap() const {return 0x0;}
+  virtual void             Print() const = 0;
+
+  static void    SetDebug(Int_t dbg=1){fgDebug=dbg;}
+  static Int_t   GetDebug(){return fgDebug;}
+
+private:
+  static Int_t fgDebug;//! debug level for all the analysis package
+
+  ClassDef(AliAODParticle,1)  // base class for AOD particles
+};
+
+#endif
diff --git a/ANALYSIS/AliAODParticleCut.cxx b/ANALYSIS/AliAODParticleCut.cxx
new file mode 100644 (file)
index 0000000..aa4a40e
--- /dev/null
@@ -0,0 +1,487 @@
+#include "AliAODParticleCut.h"
+//__________________________________________________________________________
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// class AliAODParticleCut                                                //
+//                                                                        //
+// Classes for single particle cuts                                       //
+// User should use only AliAODParticleCut, 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 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 And Function by a proper method     //
+//                                                                        //
+//                                                                        //
+// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html   //
+// resonsible: Piotr Skowronski@cern.ch                                   //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#include <Riostream.h>
+
+
+ClassImp(AliAODParticleCut)
+const Int_t AliAODParticleCut::fgkMaxCuts = 50;
+/******************************************************************/
+
+AliAODParticleCut::AliAODParticleCut():
+ fCuts(new AliAODBaseCut* [fgkMaxCuts]),//last property in the property enum => defines number of properties
+ fNCuts(0),
+ fPID(0)
+{
+  //default ctor
+}
+/******************************************************************/
+
+AliAODParticleCut::AliAODParticleCut(const AliAODParticleCut& in):
+ TObject(in)
+{
+  //cpy ctor
+  fCuts = new AliAODBaseCut* [fgkMaxCuts];//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] = (AliAODBaseCut*)in.fCuts[i]->Clone();//create new object (clone) and rember pointer to it
+   }
+}
+/******************************************************************/
+AliAODParticleCut& AliAODParticleCut::operator=(const AliAODParticleCut& in)
+{
+  //assigment operator
+  Info("operator=","operator=operator=operator=operator=\noperator=operator=operator=operator=");
+  for (Int_t i = 0;i<fNCuts;i++)
+   {
+     delete fCuts[i];
+   }
+
+  fNCuts = in.fNCuts;
+  fPID  = in.fPID;
+  for (Int_t i = 0;i<fNCuts;i++)
+   {
+     fCuts[i] = (AliAODBaseCut*)in.fCuts[i]->Clone();//create new object (clone) and rember pointer to it
+   }
+  return *this;
+}    
+
+/******************************************************************/
+AliAODParticleCut::~AliAODParticleCut()
+{
+  //dtor
+  for (Int_t i = 0;i<fNCuts;i++)
+   {
+     delete fCuts[i];
+   }
+  delete []fCuts;
+} 
+/******************************************************************/
+
+Bool_t AliAODParticleCut::Pass(AliAODParticle* p) const
+{
+//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)) )
+     {
+//       fCuts[i]->Print();
+       return kTRUE; //if one of the cuts rejects, then reject
+     }
+   }
+  return kFALSE;
+}
+/******************************************************************/
+
+void AliAODParticleCut::AddBasePartCut(AliAODBaseCut* basecut)
+{
+  //adds the base pair cut (cut on one value)
+   if (!basecut) return;
+   if( fNCuts == (fgkMaxCuts-1) )
+    {
+      Warning("AddBasePartCut","Not enough place for another cut");
+      return;
+    }
+   fCuts[fNCuts++]=basecut;
+}
+
+/******************************************************************/
+AliAODBaseCut* AliAODParticleCut::FindCut(AliAODCutProperty property)
+{
+ //returns pointer to the cut checking the given 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 AliAODParticleCut::SetMomentumRange(Double_t min, Double_t max)
+{
+  //Sets momentum range
+  AliAODMomentumCut* cut= (AliAODMomentumCut*)FindCut(kAODP);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODMomentumCut(min,max);
+}
+/******************************************************************/
+
+
+void AliAODParticleCut::SetPtRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODPtCut* cut= (AliAODPtCut*)FindCut(kAODPt);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODPtCut(min,max);
+
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetEnergyRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODEnergyCut* cut= (AliAODEnergyCut*)FindCut(kAODE);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODEnergyCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetRapidityRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODRapidity);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODRapidityCut(min,max);
+
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetPseudoRapidityRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODPseudoRapidity);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODPseudoRapidityCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetPxRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODPx);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODPxCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetPyRange(Double_t min, Double_t max)
+{  
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODPy);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODPyCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetPzRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODPz);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODPzCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetPhiRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODPhi);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODPhiCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetThetaRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODTheta);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODThetaCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetVxRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODVx);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODVxCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetVyRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODVy);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODVyCut(min,max);
+}
+/******************************************************************/
+
+void AliAODParticleCut::SetVzRange(Double_t min, Double_t max)
+{
+  //name self descriptive
+  AliAODBaseCut* cut = FindCut(kAODVz);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliAODVzCut(min,max);
+}
+
+/******************************************************************/
+void AliAODParticleCut::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()) 
+    {
+      Int_t i;
+      for (i = 0;i<fNCuts;i++) delete fCuts[i];
+      b.ReadVersion(&R__s, &R__c);
+      TObject::Streamer(b);
+      b >> fPID;
+      b >> fNCuts;
+      for (i = 0;i<fNCuts;i++) 
+       {
+         b >> fCuts[i];
+       }
+       b.CheckByteCount(R__s, R__c,AliAODParticleCut::IsA());
+    } 
+   else 
+    {
+     R__c = b.WriteVersion(AliAODParticleCut::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 AliAODParticleCut::Print(void) const
+{
+  //prints all information about the cut to stdout
+  cout<<"Printing AliAODParticleCut, 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(AliAODEmptyParticleCut)
+void AliAODEmptyParticleCut::Streamer(TBuffer &b)
+ {
+  //stramer
+  AliAODParticleCut::Streamer(b);
+ }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp(AliAODBaseCut)
+void AliAODBaseCut::Print(void) const
+{
+  // prints the information anout the base cut to stdout
+  cout<<"fMin="<<fMin <<", fMax=" <<fMax<<"    ";
+  PrintProperty();
+}
+/******************************************************************/
+
+void AliAODBaseCut::PrintProperty(void) const
+{
+ //prints the property name 
+ switch (fProperty)
+  {
+   case  kAODP: 
+     cout<<"kAODP"; break;
+   case  kAODPt: 
+     cout<<"kAODPt"; break;
+   case  kAODE: 
+     cout<<"kAODE"; break;
+   case  kAODRapidity: 
+     cout<<"kAODRapidity"; break;
+   case  kAODPseudoRapidity: 
+     cout<<"kAODPseudoRapidity"; break;
+   case  kAODPx: 
+     cout<<"kAODPx"; break;
+   case  kAODPy: 
+     cout<<"kAODPy"; break;
+   case  kAODPz: 
+     cout<<"kAODPz"; break;   
+   case  kAODPhi: 
+     cout<<"kAODPhi"; break;
+   case  kAODTheta: 
+     cout<<"kAODTheta"; break;
+   case  kAODVx: 
+     cout<<"kAODVx"; break;
+   case  kAODVy: 
+     cout<<"kAODVy"; break;
+   case  kAODVz: 
+     cout<<"kAODVz"; break;
+   case  kAODPid: 
+     cout<<"kAODPid"; break;
+   case  kAODNone: 
+     cout<<"kAODNone"; break;
+   default: 
+     cout<<"Property Not Found";
+  }
+ cout<<endl;
+}
+ClassImp( AliAODMomentumCut )
+
+ClassImp( AliAODPtCut )
+ClassImp( AliAODEnergyCut )
+ClassImp( AliAODRapidityCut )
+ClassImp( AliAODPseudoRapidityCut )
+ClassImp( AliAODPxCut )
+ClassImp( AliAODPyCut )
+ClassImp( AliAODPzCut )
+ClassImp( AliAODPhiCut )
+ClassImp( AliAODThetaCut )
+ClassImp( AliAODVxCut )
+ClassImp( AliAODVyCut )
+ClassImp( AliAODVzCut )
+
+ClassImp( AliAODPIDCut )
+
+void AliAODPIDCut::Print(void) const
+{
+  cout<<"PID "<<fPID<<" ";
+  AliAODBaseCut::Print();
+}
+
+ClassImp( AliAODLogicalOperCut )
+
+AliAODLogicalOperCut::AliAODLogicalOperCut():
+ AliAODBaseCut(-10e10,10e10,kAODNone),
+ fFirst(new AliAODDummyBaseCut),
+ fSecond(new AliAODDummyBaseCut)
+{
+ //ctor
+}
+/******************************************************************/
+
+AliAODLogicalOperCut::AliAODLogicalOperCut(AliAODBaseCut* first, AliAODBaseCut* second):
+ AliAODBaseCut(-10e10,10e10,kAODNone),
+ fFirst((first)?(AliAODBaseCut*)first->Clone():0x0),
+ fSecond((second)?(AliAODBaseCut*)second->Clone():0x0)
+{
+  //ctor
+  if ( (fFirst && fSecond) == kFALSE) 
+   {
+     Fatal("AliAODLogicalOperCut","One of parameters is NULL!");
+   }
+}
+/******************************************************************/
+
+AliAODLogicalOperCut::~AliAODLogicalOperCut()
+{
+  //destructor
+  delete fFirst;
+  delete fSecond;
+}
+/******************************************************************/
+
+Bool_t AliAODLogicalOperCut::AliAODDummyBaseCut::Pass(AliAODParticle* /*part*/)  const
+{
+  //checks if particles passes properties defined by this cut
+  Warning("Pass","You are using dummy base cut! Probobly some logical cut is not set up properly");
+  return kFALSE;//accept
+}
+/******************************************************************/
+
+void AliAODLogicalOperCut::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()) 
+   {
+     delete fFirst;
+     delete fSecond;
+     fFirst  = 0x0;
+     fSecond = 0x0;
+
+     b.ReadVersion(&R__s, &R__c);
+     TObject::Streamer(b);
+     b >> fFirst;
+     b >> fSecond;
+     b.CheckByteCount(R__s, R__c,AliAODLogicalOperCut::IsA());
+   } 
+  else 
+   {
+     R__c = b.WriteVersion(AliAODLogicalOperCut::IsA(), kTRUE);
+     TObject::Streamer(b);
+     b << fFirst;
+     b << fSecond;
+     b.SetByteCount(R__c, kTRUE);
+  }
+}
+
+/******************************************************************/
+ClassImp(AliAODOrCut)
+
+Bool_t AliAODOrCut::Pass(AliAODParticle * p) const
+{
+  //returns true when rejected 
+  //AND operation is a little bit misleading but is correct
+  //User wants to build logical cuts with natural (positive) logic
+  //while AODAN use inernally reverse (returns true when rejected)
+  if (fFirst->Pass(p) && fSecond->Pass(p)) return kTRUE;//rejected (both rejected, returned kTRUE)
+  return kFALSE;//accepted, at least one accepted (returned kFALSE)
+}
+/******************************************************************/
+
+ClassImp(AliAODAndCut)
+
+Bool_t AliAODAndCut::Pass(AliAODParticle * p)  const
+{
+  //returns true when rejected 
+  //OR operation is a little bit misleading but is correct
+  //User wants to build logical cuts with natural (positive) logic
+  //while AODAN use inernally reverse (returns true when rejected)
+  if (fFirst->Pass(p) || fSecond->Pass(p)) return kTRUE;//rejected (any of two rejected(returned kTRUE) )
+  return kFALSE;//accepted (both accepted (returned kFALSE))
+}
+/******************************************************************/
diff --git a/ANALYSIS/AliAODParticleCut.h b/ANALYSIS/AliAODParticleCut.h
new file mode 100644 (file)
index 0000000..ae7aac2
--- /dev/null
@@ -0,0 +1,390 @@
+#ifndef ALIAODPARTICLECUT_H
+#define ALIAODPARTICLECUT_H
+//__________________________________________________________________________
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// class AliAODParticleCut                                                //
+//                                                                        //
+// Classes for single particle cuts                                       //
+// User should use only AliAODParticleCut, 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 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 And Function by a proper method     //
+//                                                                        //
+//                                                                        //
+// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html   //
+// resonsible: Piotr Skowronski@cern.ch                                   //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+
+#include <TObject.h>
+#include "AliAODParticle.h"
+
+
+class AliAODEmptyParticleCut;
+class AliAODParticleCut;
+class AliAODBaseCut;
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+enum AliAODCutProperty
+ {
+//codes particle property
+  kAODP,  //Momentum
+  kAODPt, //Transverse momentum
+  kAODE,  //Energy
+  kAODRapidity, //
+  kAODPseudoRapidity,
+  kAODPx, //X coAnddinate of the momentum
+  kAODPy, //Y coAnddinate of the momentum
+  kAODPz, //Z coAnddinate of the momentum
+  kAODPhi,//angle
+  kAODTheta,//angle
+  kAODVx,  // vertex X coAnddinate
+  kAODVy,  // vertex Y coAnddinate
+  kAODVz,  // vertex Z coAnddinate
+  kAODPid, // vertex Z coAnddinate
+//_____________________________
+  kAODNone
+ };
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliAODParticleCut: public TObject
+{
+//Class describing cut on particle
+  public:
+    AliAODParticleCut();
+    AliAODParticleCut(const AliAODParticleCut& in);
+    virtual ~AliAODParticleCut();
+    AliAODParticleCut& operator = (const AliAODParticleCut& in);
+    
+    virtual Bool_t Pass(AliAODParticle* p) const;
+    Bool_t IsEmpty() const {return kFALSE;}
+    
+    void AddBasePartCut(AliAODBaseCut* basecut);
+    
+    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) const;
+  protected:
+     
+    AliAODBaseCut* FindCut(AliAODCutProperty property);
+
+    AliAODBaseCut ** fCuts;//! Array with cuts
+    Int_t fNCuts; //number of base cuts stored in fCuts
+
+    Int_t fPID; //particle PID  - if=0 (rootino) all pids are accepted
+          
+  private:
+    static const Int_t fgkMaxCuts; //Size of the fCuts array
+
+    ClassDef(AliAODParticleCut,1)
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliAODEmptyParticleCut:  public AliAODParticleCut
+{
+//Empty - it passes possitively all particles - it means returns always False
+//Class describing cut on particles
+  public:
+    AliAODEmptyParticleCut(){};
+    virtual ~AliAODEmptyParticleCut(){};
+    
+    Bool_t Pass(AliAODParticle*) const {return kFALSE;} //accept everything
+    Bool_t IsEmpty() const {return kTRUE;}
+
+    ClassDef(AliAODEmptyParticleCut,1)
+};
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliAODBaseCut: public TObject
+ {
+   //This class defines the range of some property - pure virtual
+   //Property is coded by AliAODCutTypes type
+   
+   public:
+     
+     AliAODBaseCut(Double_t min = 0.0, Double_t max = 0.0,AliAODCutProperty prop = kAODNone):
+                   fProperty(prop),fMin(min),fMax(max){}
+
+     virtual           ~AliAODBaseCut(){}
+     
+     virtual Bool_t    Pass(AliAODParticle *p) const;
+     
+     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;}
+     
+     AliAODCutProperty GetProperty() const {return fProperty;}
+     virtual void Print(void) const;
+     
+   protected:
+     virtual Double_t  GetValue(AliAODParticle *) const = 0;
+
+     AliAODCutProperty fProperty; //property that this cut describes
+     Double_t fMin;//minimum value
+     Double_t fMax;//maximum value
+     
+   private:
+     void PrintProperty(void) const;
+     ClassDef(AliAODBaseCut,1)
+   
+ };
+
+inline Bool_t
+AliAODBaseCut::Pass(AliAODParticle *p) const
+{
+  //cjecks if particle property fits in range
+  if ( (GetValue(p) < fMin) || (GetValue(p) > fMax ) ) return kTRUE; //rejected
+  else return kFALSE; //accepted
+}
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliAODMomentumCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODMomentumCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODP){}
+    virtual ~AliAODMomentumCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->P();}
+    ClassDef(AliAODMomentumCut,1)
+ };
+
+class AliAODPtCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODPtCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPt){}
+    virtual ~AliAODPtCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Pt();}
+    ClassDef(AliAODPtCut,1)
+ };
+
+
+class AliAODEnergyCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODEnergyCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODE){}
+    virtual ~AliAODEnergyCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const {return p->E();}
+    ClassDef(AliAODEnergyCut,1)
+ };
+
+class AliAODRapidityCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODRapidityCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODRapidity){}
+    virtual ~AliAODRapidityCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Y();}
+    ClassDef(AliAODRapidityCut,1)
+ };
+
+class AliAODPseudoRapidityCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODPseudoRapidityCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPseudoRapidity){}
+    virtual ~AliAODPseudoRapidityCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Eta();}
+    ClassDef(AliAODPseudoRapidityCut,1)
+ };
+
+class AliAODPxCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODPxCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPx){}
+    virtual ~AliAODPxCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Px();}
+    ClassDef(AliAODPxCut,1)
+ };
+
+class AliAODPyCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODPyCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPy){}
+    virtual ~AliAODPyCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Py();}
+    ClassDef(AliAODPyCut,1)
+ };
+
+
+class AliAODPzCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODPzCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPz){}
+    virtual ~AliAODPzCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Pz();}
+    ClassDef(AliAODPzCut,1)
+ };
+
+class AliAODPhiCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODPhiCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPhi){}
+    virtual ~AliAODPhiCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Phi();}
+    ClassDef(AliAODPhiCut,1)
+  
+ };
+
+class AliAODThetaCut: public AliAODBaseCut
+ {
+  public: 
+    AliAODThetaCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODTheta){}
+    virtual ~AliAODThetaCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Theta();}
+    ClassDef(AliAODThetaCut,1)
+  
+ };
+
+class AliAODVxCut: public AliAODBaseCut
+ {
+ //Cut of the X coAnddinate of the vertex position
+  public: 
+    AliAODVxCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODVx){}
+    virtual ~AliAODVxCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Vx();} //retruns value of the vertex
+    ClassDef(AliAODVxCut,1)
+  
+ };
+
+
+class AliAODVyCut: public AliAODBaseCut
+ {
+ //Cut of the X coAnddinate of the vertex position
+  public: 
+    AliAODVyCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODVy){}
+    virtual ~AliAODVyCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Vy();} //retruns value of the vertex
+    ClassDef(AliAODVyCut,1)
+  
+ };
+
+class AliAODVzCut: public AliAODBaseCut
+ {
+ //Cut of the X coAnddinate of the vertex position
+  public: 
+    AliAODVzCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODVz){}
+    virtual ~AliAODVzCut(){}
+  protected:
+    Double_t  GetValue(AliAODParticle * p)const{return p->Vz();} //retruns value of the vertex
+    
+    ClassDef(AliAODVzCut,1)
+  
+ };
+
+class AliAODPIDCut:  public AliAODBaseCut
+ {
+   public:
+     AliAODPIDCut():AliAODBaseCut(0.0,0.0,kAODPid),fPID(0){}
+     AliAODPIDCut(Int_t pid, Double_t min = 0.0, Double_t max = 1.0):AliAODBaseCut(min,max,kAODPid),fPID(pid){}
+     virtual ~AliAODPIDCut(){}
+     
+     void SetPID(Int_t pid){fPID = pid;}
+     void Print(void) const;
+   protected:
+     Double_t  GetValue(AliAODParticle * p)const{return p->GetProbability(fPID);}
+     Int_t     fPID; //pid of particle that the pid is set 
+     ClassDef(AliAODPIDCut,1)
+ };
+//___________________________________________________
+/////////////////////////////////////////////////////
+//                                                 //
+// class AliAODLogicalOperCut                      //
+//                                                 //
+// This cut is base class fAnd class that perfAndms  //
+// logical operations on cuts                      //
+//                                                 //
+/////////////////////////////////////////////////////
+class AliAODLogicalOperCut:  public AliAODBaseCut
+ {
+   public:
+     AliAODLogicalOperCut();
+     AliAODLogicalOperCut(AliAODBaseCut* first, AliAODBaseCut* second);
+     virtual   ~AliAODLogicalOperCut();
+   protected:
+     Double_t  GetValue(AliAODParticle * /*part*/) const {MayNotUse("GetValue");return 0.0;}
+     
+     AliAODBaseCut* fFirst;   //second cut
+     AliAODBaseCut* fSecond;  //first cut
+   private:  
+    class  AliAODDummyBaseCut: public AliAODBaseCut 
+     {
+       Double_t  GetValue(AliAODParticle * /*part*/) const {return 0.0;}
+       Bool_t    Pass(AliAODParticle* /*part*/) const;
+     };
+     
+    ClassDef(AliAODLogicalOperCut,1)
+ };
+
+class AliAODOrCut: public AliAODLogicalOperCut
+{
+   public:
+     AliAODOrCut(){}
+     AliAODOrCut(AliAODBaseCut* first, AliAODBaseCut* second):AliAODLogicalOperCut(first,second){}
+     virtual   ~AliAODOrCut(){}
+     Bool_t    Pass(AliAODParticle *p) const;
+     ClassDef(AliAODOrCut,1)
+};
+
+class AliAODAndCut: public AliAODLogicalOperCut
+{
+   public:
+     AliAODAndCut(){}
+     AliAODAndCut(AliAODBaseCut* first, AliAODBaseCut* second):AliAODLogicalOperCut(first,second){}
+     virtual   ~AliAODAndCut(){}
+     Bool_t    Pass(AliAODParticle *p) const;
+     ClassDef(AliAODAndCut,1)
+};
+
+#endif
diff --git a/ANALYSIS/AliAODRun.cxx b/ANALYSIS/AliAODRun.cxx
new file mode 100644 (file)
index 0000000..1f96274
--- /dev/null
@@ -0,0 +1,86 @@
+#include "AliAODRun.h"
+//____________________
+///////////////////////////////////////////////////////
+//                                                   //
+// AliAODRun                                         //
+//                                                   //
+// Class storing and managing events                 //
+//                                                   //
+// Piotr.Skowronski@cern.ch                          //
+// http://alisoft.cern.ch/people/skowron/analyzer    //
+//                                                   //
+///////////////////////////////////////////////////////
+
+#include <TObjArray.h>
+
+ClassImp(AliAODRun)
+/**************************************************************************/ 
+
+AliAODRun::AliAODRun()
+{ 
+ //contructor
+  fEvents = new TObjArray();//create array for AliAODs
+  if(!fEvents) Fatal("AliAODRun::AliAODRun","Can not allocate memory");
+  fEvents->SetOwner(); //array is an owner: when is deleted or cleared it deletes objects that it contains
+}
+/**************************************************************************/
+
+AliAODRun::~AliAODRun()
+{
+  //destructor
+  delete fEvents;//delete array with events
+}
+/**************************************************************************/
+
+void AliAODRun::Reset()
+ { 
+   fEvents->Clear();//clear an array with events. 
+                    //All events are deleted because array is an owner (set in constructor)
+ }
+/**************************************************************************/
+
+void AliAODRun::AddParticle(Int_t event, AliAODParticle* 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 AliAOD, event);
+ GetEvent(event)->AddParticle(part);
+}
+/**************************************************************************/
+
+void AliAODRun::AddParticle(Int_t event, TParticle* part, Int_t idx)
+{
+ //if there is no event of this number, crate it and add to the collection
+ if(!GetEvent(event))  fEvents->AddAtAndExpand(new AliAOD, event);
+ GetEvent(event)->AddParticle(part,idx);
+}
+/**************************************************************************/
+
+void AliAODRun::AddParticle(Int_t event, Int_t pdg, Int_t idx,
+                            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 AliAOD, event);
+ GetEvent(event)->AddParticle(pdg,idx,px,py,pz,etot,vx,vy,vz,time);
+}
+/**************************************************************************/ 
+
+void AliAODRun::SetEvent(Int_t number, AliAOD* event)
+{
+  //adds an event to the run
+  if (event == 0x0)
+   {
+     delete fEvents->RemoveAt(number);
+     return;
+   }
+  AliAOD* ev = GetEvent(number);
+  if (ev == event) return;
+  
+  delete fEvents->RemoveAt(number);
+  fEvents->AddAtAndExpand(event, number);
+  
+}
+/**************************************************************************/ 
+
diff --git a/ANALYSIS/AliAODRun.h b/ANALYSIS/AliAODRun.h
new file mode 100644 (file)
index 0000000..bbd7d81
--- /dev/null
@@ -0,0 +1,94 @@
+#ifndef ALIAODRUN_H
+#define ALIAODRUN_H
+//____________________
+///////////////////////////////////////////////////////
+//                                                   //
+// AliAODRun                                         //
+//                                                   //
+// Class storing and managing set events             //
+// designed for fast acces                           //
+//                                                   //
+// Piotr.Skowronski@cern.ch                          //
+// http://alisoft.cern.ch/people/skowron/analyzer    //
+//                                                   //
+///////////////////////////////////////////////////////
+
+
+#include "AliAOD.h"
+#include <TObjArray.h>
+
+class AliAODParticle;
+class TParticle;
+
+class AliAODRun: public TObject
+ {
+  public:
+    AliAODRun();
+    virtual ~AliAODRun();
+
+    void            AddParticle(Int_t event, AliAODParticle* part); //inerface to AliAOD::AddParticle(AliAODParticle*) 
+    void            AddParticle(Int_t event, TParticle* part, Int_t idx);//inerface to AliAOD::AddParticle(TParticle*) 
+    
+    //inerface to AliAOD::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, Int_t idx,
+                                Double_t px, Double_t py, Double_t pz, Double_t etot,
+                                Double_t vx, Double_t vy, Double_t vz, Double_t time); 
+    
+    void            SetEvent(Int_t number, AliAOD* event);
+    AliAODParticle* GetParticle(Int_t event, Int_t n); //returns nth particle from event
+    AliAOD*    GetEvent(Int_t event) const; //returns AliAOD 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 AliAODs
+  private:
+    
+  public:
+    ClassDef(AliAODRun,1)
+ };
+/**************************************************************************/
+
+inline
+AliAOD* AliAODRun::GetEvent(Int_t event) const
+{
+//returns pointer to AliAOD number "event"
+  //check if array is enough big - protect from error massage from array "Out Of Bounds"
+  if (event>=fEvents->GetSize()) return 0x0;//WARNING, that line relies 
+                                            //on index of first object in TObjArray is 0
+                   //== LowerBound = 0
+  return (AliAOD*)fEvents->At(event);
+}
+/**************************************************************************/
+inline
+AliAODParticle* AliAODRun::GetParticle(Int_t event, Int_t n) 
+{
+ //returns nth particle from event number event
+  AliAOD* e = GetEvent(event);
+  return (e)?e->GetParticle(n):0x0;
+}
+
+/**************************************************************************/
+
+inline
+Int_t AliAODRun::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 AliAODRun::GetNumberOfParticlesInEvent(Int_t event) const
+{
+//returns number of Particles in event 
+  AliAOD* e = GetEvent(event);
+  return (e)?e->GetNumberOfParticles():0x0;
+}
+#endif
diff --git a/ANALYSIS/AliAODStdParticle.cxx b/ANALYSIS/AliAODStdParticle.cxx
new file mode 100644 (file)
index 0000000..628a7cf
--- /dev/null
@@ -0,0 +1,410 @@
+#include "AliAODStdParticle.h"
+//___________________________________________________________
+/////////////////////////////////////////////////////////////
+//
+// class AliAODStdParticle
+//
+// 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
+// Additionaly adjusted for HBT Analysies purposes
+// + pointer to Track Points
+// + pointer to Cluster Map(s)
+//
+// Piotr.Skowronski@cern.ch
+//
+/////////////////////////////////////////////////////////////
+#include <TParticle.h>
+#include "AliTrackPoints.h"
+#include "AliClusterMap.h"
+
+ClassImp(AliAODStdParticle)
+
+//______________________________________________________________________________
+AliAODStdParticle::AliAODStdParticle():  
+ fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
+ fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0),
+ fTrackPoints(0x0),fClusterMap(0x0)
+{//empty particle
+}
+//______________________________________________________________________________
+
+AliAODStdParticle::AliAODStdParticle(Int_t pdg, Int_t idx,
+               Double_t px, Double_t py, Double_t pz, Double_t etot,
+               Double_t vx, Double_t vy, Double_t vz, Double_t time):  
+  fPdgIdx(0), fIdxInEvent(idx),fNPids(0),fPids(0x0),fPidProb(0x0),
+  fCalcMass(0), 
+  fPx(px), fPy(py),fPz(pz),fE(etot), 
+  fVx(vx), fVy(vy),fVz(vz),fVt(time),
+  fTrackPoints(0x0),fClusterMap(0x0)
+{
+//mormal constructor
+  SetPdgCode(pdg);
+  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);
+  }
+}
+//______________________________________________________________________________
+
+AliAODStdParticle::AliAODStdParticle(Int_t pdg, Float_t prob, Int_t idx, 
+                               Double_t px, Double_t py, Double_t pz, Double_t etot,
+                               Double_t vx, Double_t vy, Double_t vz, Double_t time):
+  fPdgIdx(0), fIdxInEvent(idx),fNPids(0),fPids(0x0),fPidProb(0x0),
+  fCalcMass(0), 
+  fPx(px), fPy(py),fPz(pz),fE(etot), 
+  fVx(vx), fVy(vy),fVz(vz),fVt(time),
+  fTrackPoints(0x0),fClusterMap(0x0)
+{
+//mormal constructor
+  SetPdgCode(pdg,prob);
+  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);
+  }
+}
+//______________________________________________________________________________
+AliAODStdParticle::AliAODStdParticle(const AliAODStdParticle& in):
+   AliAODParticle(in),
+   fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
+   fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
+   fCalcMass(in.GetCalcMass()),
+   fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.Energy()), 
+   fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T()),
+   fTrackPoints(0x0), fClusterMap(0x0)
+{
+ //Copy constructor
+ for(Int_t i = 0; i<fNPids; i++)
+  {
+    fPids[i] =  in.fPids[i];
+    fPidProb[i] = in.fPidProb[i];
+  }
+ if (in.fTrackPoints)
+   fTrackPoints = (AliTrackPoints*)in.fTrackPoints->Clone();
+ if (in.fClusterMap)
+   fClusterMap = (AliClusterMap*)in.fClusterMap->Clone();
+}
+
+//______________________________________________________________________________
+AliAODStdParticle::AliAODStdParticle(const TParticle &p,Int_t idx):
+   fPdgIdx(0), fIdxInEvent(idx),
+   fNPids(0),fPids(0x0),fPidProb(0x0),
+   fCalcMass(p.GetCalcMass()),
+   fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()), 
+   fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T()),
+   fTrackPoints(0x0),fClusterMap(0x0)
+{
+ //all copied in the initialization
+ SetPdgCode(p.GetPdgCode());
+}
+//______________________________________________________________________________
+
+AliAODStdParticle::~AliAODStdParticle()
+{
+//dtor  
+  delete [] fPids;
+  delete [] fPidProb;
+  delete fTrackPoints;
+  delete fClusterMap;
+}
+//______________________________________________________________________________
+
+AliAODStdParticle& AliAODStdParticle::operator=(const AliAODStdParticle& in)
+{
+//assigment operator
+  
+  fNPids = in.fNPids;
+  delete [] fPids;
+  delete [] fPidProb;
+  Int_t* fPids = new Int_t[fNPids];
+  Float_t* fPidProb = new Float_t[fNPids];
+  for (Int_t i = 0; i < fNPids;i++)
+   {
+     fPids[i]    = in.fPids[i];
+     fPidProb[i] = in.fPidProb[i];
+   }
+   
+  fPdgIdx = in.fPdgIdx; 
+  fIdxInEvent = in.fIdxInEvent;  
+  fCalcMass = in.GetCalcMass();
+  fPx = in.Px();
+  fPy = in.Py();
+  fPz = in.Pz();
+  fE = in.Energy(); 
+  fVx = in.Vx();
+  fVy = in.Vy();
+  fVz = in.Vz();
+  fVt = in.T();
+  
+  delete fTrackPoints;
+  fTrackPoints = (in.fTrackPoints)?(AliTrackPoints*)fTrackPoints->Clone():0x0;
+  
+  delete fClusterMap;
+  fClusterMap =  (in.fClusterMap)?(AliClusterMap*)in.fClusterMap->Clone():0x0;
+  
+  return *this;
+}
+//______________________________________________________________________________
+
+void AliAODStdParticle::SetPdgCode(Int_t pdg,Float_t prob)
+{
+  SetPIDprobability(pdg,prob);
+  fPdgIdx = GetPidSlot(pdg);
+}
+
+//______________________________________________________________________________
+void AliAODStdParticle::SetPIDprobability(Int_t pdg, Float_t prob)
+{
+//Sets another pdg code and corresponding probabilty
+//Ids are set in decreasing order
+//Check if total prbaility is not ivercoming unity is performed
+//in case, warning is printed
+  if (GetDebug() > 9) Info("SetPIDprobability","Setting PID %d prob %f",pdg,prob);
+
+  Float_t totprob = 0.0;//sums up probabilities
+  Int_t idx = GetPidSlot(pdg);
+  Int_t i;
+  if (idx > -1) 
+   {
+     fPidProb[idx] = prob;
+     for (i = 0; i < fNPids;i++) totprob+=fPidProb[i];
+     if (totprob > (1.0+0.000001))
+       {
+         Warning("SetPIDprobability","Total probability greater than unity (%f)",totprob);
+       }
+     if (GetDebug() > 9) 
+      {
+        Info("SetPIDprobability","Current Total probability: %f",totprob);
+      }
+     return;
+   }
+    
+  Int_t currentpid = GetPdgCode();
+  fNPids++;
+  Float_t* aPidProbNew = new Float_t[fNPids];
+  Int_t* aPidsNew = new Int_t[fNPids];
+  
+  for (i = 0; i < fNPids-1;i++)//find a slot
+   {
+     if ( fPidProb[i] > prob)
+      {
+        if (GetDebug()>9) Info("SetPID","Copying entry %d",i);
+        aPidProbNew[i] = fPidProb[i];
+        aPidsNew[i] = fPids[i];
+        totprob+=fPidProb[i];
+      }
+     else break;
+   }
+
+  if (GetDebug() > 9) Info("SetPID","Setting new PID on entry %d",i);
+  aPidProbNew[i] = prob;
+  aPidsNew[i] = pdg;
+  totprob+=prob;
+  
+
+  for (Int_t j = fNPids-1; j > i ;j--)//copy rest of old araays 
+   {
+     if (GetDebug() > 9) Info("SetPID","Copying from old entry %d to new entry %d",j-1,j);
+     aPidProbNew[j] = fPidProb[j-1];
+     aPidsNew[j] = fPids[j-1];
+     totprob+=fPidProb[j-1];
+   }
+
+  delete [] fPidProb;
+  delete [] fPids;
+  
+  fPidProb = aPidProbNew;
+  fPids = aPidsNew;
+  
+  fPdgIdx = GetPidSlot(currentpid);
+  if (fPdgIdx == -1) fPdgIdx = 0;
+  
+  if (totprob > (1.0+0.000001))//place for numerical error
+   {
+     Warning("SetId","Total probability is greater than unity (%f)!!!",totprob);
+     Print();
+   }
+}
+//______________________________________________________________________________
+
+Float_t AliAODStdParticle::GetPIDprobability(Int_t pdg) const
+{
+//Returns probability that this particle is the type of pdg
+  Int_t idx = GetPidSlot(pdg);
+  if (idx < 0) return 0.0;//such pid was not specified for this particle
+  return fPidProb[idx];
+}
+//______________________________________________________________________________
+
+const Char_t* AliAODStdParticle::GetName() const 
+{
+  //returns name of this particle 
+   static char def[4] = "XXX";
+   const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
+   if (ap) return ap->GetName();
+   else    return def;
+}
+//______________________________________________________________________________
+
+Int_t AliAODStdParticle::GetPidSlot(Int_t pdg) const
+{
+ //returns position of the given PID in fPids (and fPidProb) array.
+ if (fPids == 0x0) return -1;
+ for (Int_t i = 0; i< fNPids; i++)
+  {
+   if (fPids[i] == pdg) return i;
+  }
+ return -1;
+}
+//______________________________________________________________________________
+
+Int_t AliAODStdParticle::GetNthPid(Int_t idx) const
+{
+  //returns PID sitting on slot idx in fPids
+  if ( (idx < 0) || (idx >= fNPids) )
+   {
+     Error("GetNthPid","Out Of Bounds");
+     return 0;
+   }
+  return fPids[idx];
+}
+//______________________________________________________________________________
+
+Float_t AliAODStdParticle::GetNthPidProb(Int_t idx) const
+{
+  //returns PID sitting on slot idx in fPidProb
+  if ( (idx < 0) || (idx >= fNPids) )
+   {
+     Error("GetNthPid","Out Of Bounds");
+     return 0;
+   }
+  return fPidProb[idx];
+}
+//______________________________________________________________________________
+
+void AliAODStdParticle::Print() const
+{
+//prints information about particle
+  printf("____________________________________________________\n");
+  printf("Idx: %d  PID: %d  Name: ",fIdxInEvent,GetPdgCode());
+  TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
+  if (pdgp)
+   {
+     printf("%s Mass: %f\n",pdgp->GetName(),pdgp->Mass());
+   }
+  else
+   {
+     printf("Not known\n");
+   }
+  
+  printf("Px: %+f Py: %+f Pz: %+f E: %+f Calculated Mass: %f\nVx: %+f Vy: %+f Vz: %+f T: %+f\n",
+          Px(),Py(),Pz(),Energy(),GetCalcMass(),Vx(),Vy(),Vz(),T());
+
+  for (Int_t i = 0; i < fNPids; i++)
+   {
+     printf("# %d  PID: %d  Probability %f name ",i,fPids[i],fPidProb[i]);
+     const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPids[i]);
+     if (ap)
+      {
+        printf("%s Mass %f\n",ap->GetName(),ap->Mass());
+      }
+     else
+      {
+        printf("Not known\n");
+      }
+   }
+}
+
+//______________________________________________________________________________
+
+//void AliAODStdParticle::Streamer(TBuffer &b)
+//{
+//     // Stream all objects in the array to or from the I/O buffer.
+//   UInt_t R__s, R__c;
+//   Int_t i;
+//   if (b.IsReading()) 
+//    {
+//      delete [] fPids;
+//      delete [] fPidProb;
+//      
+//      Version_t v = b.ReadVersion(&R__s, &R__c);
+//      if (v == 1)
+//       {
+//         AliAODStdParticle::Class()->ReadBuffer(b, this);
+//      }
+//      else
+//       {
+//        TObject::Streamer(b);
+//       b >> fPdgIdx;
+//       b >> fIdxInEvent;
+//       
+//       b >> fNPids;
+//       Int_t* fPids = new Int_t[fNPids];
+//        Float_t* fPidProb = new Float_t[fNPids];
+//        for (i = 0;i<fNPids;i++) 
+//         {
+//           b >> fPids[i];
+//         }
+//        for (i = 0;i<fNPids;i++) 
+//        {
+//          b >> fPidProb[i];
+//         }
+//        b >> fCalcMass;
+//
+//        b >> fPx;
+//       b >> fPy;
+//       b >> fPz;
+//        b >> fE;
+//
+//       b >> fVx;
+//        b >> fVy;
+//        b >> fVz;
+//       b >> fVt;
+//       Info("Streamer","Read data");
+//        Print();
+//       }
+//
+//      b.CheckByteCount(R__s, R__c,AliAODStdParticle::IsA());
+//    } 
+//  else 
+//   {
+//     R__c = b.WriteVersion(AliAODStdParticle::IsA(), kTRUE);
+//     TObject::Streamer(b);
+//     Info("Streamer","Read data");
+//     Print();
+//
+//     b << fPdgIdx;
+//     b << fIdxInEvent;
+//     b << fNPids;
+//     for (i = 0;i<fNPids;i++) 
+//      {
+//        b << fPids[i];
+//      }
+//      {
+//      {
+//     for (i = 0;i<fNPids;i++) 
+//     {
+//        b << fPidProb[i];
+//      }
+//     b << fCalcMass;
+//
+//     b << fPx;
+//     b << fPy;
+//     b << fPz;
+//     b << fE;
+//
+//     b << fVx;
+//     b << fVy;
+//     b << fVz;
+//     b << fVt;
+//
+//    b.SetByteCount(R__c, kTRUE);
+//   }
+//}
diff --git a/ANALYSIS/AliAODStdParticle.h b/ANALYSIS/AliAODStdParticle.h
new file mode 100644 (file)
index 0000000..1462436
--- /dev/null
@@ -0,0 +1,161 @@
+#ifndef ALIAODSTDPARTICLE_H
+#define ALIAODSTDPARTICLE_H
+//___________________________________________________________
+/////////////////////////////////////////////////////////////
+//
+// class AliAODStdParticle
+//
+// Ali HBT Particle: simplified class TParticle
+// Simplified in order to minimize the size of the object
+//  - we want to keep a lot of such a objects in memory
+// Additionaly adjusted for HBT Analysies purposes
+// + pointer to Track Points
+// + pointer to Cluster Map(s)
+//
+// Piotr.Skowronski@cern.ch
+//
+/////////////////////////////////////////////////////////////
+
+#include "AliAODParticle.h"
+//#include <TLorentzVector.h>
+//#include <TMath.h>
+#include <TDatabasePDG.h>
+
+
+class TParticle;
+class AliTrackPoints;
+class AliClusterMap;
+
+class AliAODStdParticle: public AliAODParticle
+{
+public:
+                                // ****** constructors and destructor
+  AliAODStdParticle();
+  AliAODStdParticle(const AliAODStdParticle& in); 
+  AliAODStdParticle(Int_t pdg, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot,
+                 Double_t vx, Double_t vy, Double_t vz, Double_t time);
+
+  AliAODStdParticle(Int_t pdg, Float_t prob, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot,
+                 Double_t vx, Double_t vy, Double_t vz, Double_t time);
+
+  AliAODStdParticle(const TParticle& p, Int_t idx);
+
+  virtual ~AliAODStdParticle();
+  
+  AliAODStdParticle& operator=(const AliAODStdParticle& in); 
+  
+  TLorentzVector FourMomentum() const {TLorentzVector v(fPx,fPy,fPz,fE);return v;}
+  TVector3       ProductionVertex() const {TVector3 v(fVx,fVy,fVz);return v;}
+  
+  void           SetPIDprobability(Int_t pdg, Float_t prob = 1.0);
+  Float_t        GetPIDprobability(Int_t pdg) const;
+  Double_t       GetProbability(Int_t pdg) const {return GetPIDprobability(pdg);}
+  Int_t          GetMostProbable() const { return (fPids)?fPids[0]:0;}
+  
+  Int_t          GetPdgCode      () const { return (fPids)?fPids[fPdgIdx]:0;}
+
+  Float_t        GetPidProb      () const { return (fPidProb)?fPidProb[fPdgIdx]:0;}
+  
+  Int_t          GetUID          () const { return fIdxInEvent;}
+  Int_t          GetNumberOfPids () const { return fNPids;}
+  Int_t          GetNthPid         (Int_t idx) const;
+  Float_t        GetNthPidProb     (Int_t idx) const;
+      
+  void           SetPdgCode(Int_t pdg, Float_t prob = 1.0);
+  Double_t       GetCalcMass     () const { return fCalcMass; }
+  Double_t       Mass            () const { return (GetPDG())?GetPDG()->Mass():-1.;}
+  
+  
+  TParticlePDG*  GetPDG          () const {return TDatabasePDG::Instance()->GetParticle(GetPdgCode());}
+  Double_t       Charge          () const { return (GetPDG())?GetPDG()->Charge():1.e8;}
+  
+  Int_t          Beauty          ()  { return GetPDG()->Beauty(); }
+  Int_t          Charm           ()  { return GetPDG()->Charm(); }
+  Int_t          Strangeness     ()  { return GetPDG()->Strangeness();}
+  void ProductionVertex(TLorentzVector &v) const { 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) const { 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 TMath::Pi()+TMath::ATan2(-fPy,-fPx); }
+
+  Double_t         Theta () const { return (fPz==0)?TMath::PiOver2():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());}
+  void           SetCalcMass(Double_t mass) {fCalcMass = mass;}
+  
+  void           SetUID(Int_t id){fIdxInEvent = id;}
+  
+  const Char_t*  GetName() const; 
+  void           Print() const;
+  
+  void           SetTrackPoints(AliTrackPoints* tpts){fTrackPoints = tpts;}
+  AliTrackPoints* GetTrackPoints() const {return fTrackPoints;}
+  void           SetClusterMap(AliClusterMap* cm){fClusterMap = cm;}
+  AliClusterMap* GetClusterMap() const {return fClusterMap;}
+  
+  
+protected:
+  Int_t          GetPidSlot(Int_t pdg) const;//returns position of the given PID in fPids (and fPidProb) array.
+
+private:
+  Char_t         fPdgIdx;               // index of PDG code of the particle in fPids
+  Int_t          fIdxInEvent;           // index of a particle: the same particle can appear in the event
+                                        //  many times with different pid's. Idx allows to check that they are the same particles
+  Int_t          fNPids;                // number of non-zero proboble Pids
+  Int_t         *fPids;                 // [fNPids] Array with PIDs
+  Float_t       *fPidProb;              // [fNPids] PIDs probabilities
+  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
+
+  AliTrackPoints* fTrackPoints;      // track positions along trajectory - used by anti-merging cut
+  AliClusterMap*  fClusterMap;       // bit map of cluters occupation; 1 if has cluter on given layer/padrow/...
+    
+  ClassDef(AliAODStdParticle,3)  // TParticle vertex particle information
+};
+
+#endif
index 5dc2f68..1a2fd96 100644 (file)
@@ -4,9 +4,11 @@
 //
 // class AliAnalysis
 //
-// Base class for analysis
-//
-//
+// Base class for analysis.
+// Each inheriting calss must define 3 methods:
+//   - Init() : that is called before event processing
+//   - ProcessEvent(AliESD*,AliStack*)
+//   - 
 // Piotr.Skowronski@cern.ch
 //
 ///////////////////////////////////////////////////////////
@@ -14,8 +16,6 @@
 
 ClassImp(AliAnalysis)
 
-Int_t AliAnalysis::fgkDebug = 0;
-
 AliAnalysis::AliAnalysis()
 {
  //ctor
index e61e8be..4b1e572 100644 (file)
@@ -14,7 +14,7 @@
 
 #include <TTask.h>
 
-class AliESD;
+class AliAOD;
 class AliStack;
  
 class AliAnalysis: public TTask
@@ -25,17 +25,12 @@ class AliAnalysis: public TTask
     virtual ~AliAnalysis();
     
     virtual Int_t Init() = 0;
-    virtual Int_t ProcessEvent(AliESD* esd, AliStack* stack = 0x0) = 0;
+    virtual Int_t ProcessEvent(AliAOD* aodrec, AliAOD* aodsim = 0x0) = 0;
     virtual Int_t Finish() = 0;
     
-    
-    static Int_t GetDebug() {return fgkDebug;}
-    static void  SetDebug(Int_t level) {fgkDebug = level;}
-    
   protected:
     
   private:
-    static Int_t fgkDebug;//! debug level
     ClassDef(AliAnalysis,1)
 };
 
index 1561ab2..d9923a3 100644 (file)
@@ -13,7 +13,7 @@
 
 #include "TObject.h"
 
-class AliESD;
+class AliAOD;
 
 class AliBaseEventCut: public TObject
 {
@@ -21,9 +21,9 @@ class AliBaseEventCut: public TObject
     AliBaseEventCut();
     virtual ~AliBaseEventCut(){}
     
-    virtual Bool_t Pass(AliESD* esd) const;//returns kTRUE if rejected
+    virtual Bool_t Pass(AliAOD* aod) const;//returns kTRUE if rejected
   protected:
-    virtual Double_t GetValue(AliESD* esd) const = 0;
+    virtual Double_t GetValue(AliAOD* aod) const = 0;
     
     Double_t fMin;//Minimum value
     Double_t fMax;//Maximum value
diff --git a/ANALYSIS/AliClusterMap.cxx b/ANALYSIS/AliClusterMap.cxx
new file mode 100644 (file)
index 0000000..fafdba5
--- /dev/null
@@ -0,0 +1,192 @@
+#include "AliClusterMap.h"
+
+//_________________________________________________
+///////////////////////////////////////////////////
+//
+// class AliClusterMap
+//
+// class that describes cluster occupation at TPC
+// Each padraw has a corresponding bit in fPadRawMap
+// 
+//
+// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////////////////////
+
+
+#include "AliESDtrack.h"
+#include "AliTPCtrack.h"
+#include "AliAODParticle.h"
+#include <TString.h>
+const Int_t AliClusterMap::fNPadRows = 159;
+
+AliClusterMap::AliClusterMap():
+ fPadRawMap(fNPadRows)
+{
+//ctor
+
+}
+/***********************************************************************/
+AliClusterMap::AliClusterMap(AliESDtrack* track):
+ fPadRawMap( (track)?track->GetTPCClusterMap():fNPadRows )
+{
+ //ctor
+ if (AliAODParticle::GetDebug() > 2)
+  { 
+    Info("AliClusterMap(AliESDtrack*)","");
+    Print();
+  } 
+} 
+/***********************************************************************/
+
+AliClusterMap::AliClusterMap(AliTPCtrack* track):
+ fPadRawMap(fNPadRows)
+{
+ //ctor
+ //Does not work since indeces in the claster index array 
+ //in the TPC track does not correspond to the padraw segmatation
+ if (AliAODParticle::GetDebug() > 9) 
+   Info("AliClusterMap",
+      "#####################################################################"); 
+ if (track == 0x0)
+  {
+    Error("AliClusterMap","Pointer to TPC track is NULL");
+    return;
+  }
+ Int_t prevrow = -1;
+ Int_t i = 0;
+ for ( ; i < track->GetNumberOfClusters(); i++)
+  {
+    Int_t idx = track->GetClusterIndex(i);
+    Int_t sect = (idx&0xff000000)>>24;
+    Int_t row = (idx&0x00ff0000)>>16;
+    if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
+    if (AliAODParticle::GetDebug() > 9)  
+      Info("AliClusterMap","Cl.idx is %d, sect %d, row %d",idx,sect,row);
+      
+    fPadRawMap.SetBitNumber(row,kTRUE);
+    
+    //Fill the gap between previous row and this row with 0 bits
+    if (prevrow < 0) 
+     {
+       prevrow = row;//if previous bit was not assigned yet == this is the first one
+     }
+    else
+     { //we don't know the order (inner to outer or reverse)
+       //just to be save in case it is going to change
+       Int_t n = 0, m = 0;
+       if (prevrow < row)
+        {
+          n = prevrow;
+          m = row;
+        }
+       else
+        {
+          n = row;
+          m = prevrow;
+        }
+       for (Int_t j = n+1; j < m; j++)
+        {
+          fPadRawMap.SetBitNumber(j,kFALSE);
+        }
+       prevrow = row; 
+     }
+  }
+  
+ if (AliAODParticle::GetDebug() > 2)
+  { 
+    Info("AliClusterMap(AliTPCtrack*)","");
+    Print();
+  } 
+}
+/***********************************************************************/
+
+void AliClusterMap::Print() const
+{
+//Prints the bit map 
+  TString msg;
+  for ( Int_t i = 0; i < fNPadRows; i++)
+   {
+     if ( fPadRawMap.TestBitNumber(i) )
+      {
+        msg+="1";
+      }
+     else
+      {
+        msg+="0";
+      }
+   }
+  Info("AliClusterMap","BitMap is\n  %s",msg.Data());
+  
+}
+
+/***********************************************************************/
+
+Float_t AliClusterMap::GetOverlapFactor(const AliClusterMap& clmap) const
+{
+  //Returns quality factor FQ = Sum(An)/Sum(clusters)
+  //      | -1; if both tracks have a cluster on padrow n
+  //An = <  0; if neither track has a cluster on padrow n
+  //     |  1; if only one trackhas a cluster on padrow n
+  // Returned value ranges between 
+  //  -0.5 (low probability that these tracks are a split track)
+  //  and
+  //   1.0 (high probability that these tracks are a split track)
+  TString msg1;
+  TString msg2;
+  
+  Int_t nh = 0;
+  Int_t an = 0;
+  for ( Int_t i = 0; i < fNPadRows; i++)
+   {
+     Bool_t x = HasClAtPadRow(i);
+     Bool_t y = clmap.HasClAtPadRow(i);
+     
+     if (x) msg1+="1";else msg1+="0";
+     if (y) msg2+="1";else msg2+="0";
+     
+     if (x && y)//both have clasters
+      {
+       an--;
+       nh+=2;
+      }
+     else 
+      {
+       
+       if (x || y)//only one have cluters
+        {
+          an++;
+          nh++;
+        }
+      }
+   }
+  
+  
+  Float_t retval = 0.0;
+  if (nh > 0) retval = ((Float_t)an)/((Float_t)nh);
+  else Warning("GetOverlapFactor","Number of counted cluters is 0.");
+  
+  if (AliAODParticle::GetDebug() > 2)
+   {
+     Info("GetOverlapFactor","Splitting Quality Factor is %f. SumAn = %d, SumClusters %d",retval,an,nh); 
+     if (retval == 1.0) 
+      { 
+        Print();
+        Info("AliClusterMap","BitMap is\n  %s\n",msg1.Data());
+        clmap.Print(); 
+        Info("AliClusterMap","BitMap is\n  %s\n\n\n\n",msg2.Data());
+      }
+     if (retval == -.5) 
+      { 
+        Print();
+        Info("AliClusterMap","BitMap is\n  %s\n",msg1.Data());
+        clmap.Print(); 
+        Info("AliClusterMap","BitMap is\n  %s\n\n\n\n",msg2.Data());
+      }
+   } 
+  return retval;
+}
diff --git a/ANALYSIS/AliClusterMap.h b/ANALYSIS/AliClusterMap.h
new file mode 100644 (file)
index 0000000..f7d8774
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALICLUSTERMAP_H
+#define ALICLUSTERMAP_H
+//_________________________________________________
+///////////////////////////////////////////////////
+//
+// class AliClusterMap
+//
+// class that describes cluster occupation at TPC
+// Each padraw has a corresponding bit in fPadRawMap
+// 
+//
+// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////////////////////
+
+
+#include <TObject.h>
+#include <TBits.h>
+
+class AliTPCtrack;
+class AliESDtrack;
+class TBits;
+
+class AliClusterMap: public TObject
+{
+  public:
+   AliClusterMap();
+   AliClusterMap(AliTPCtrack* track);
+   AliClusterMap(AliESDtrack* track);
+   virtual ~AliClusterMap(){}
+   Float_t GetOverlapFactor(const AliClusterMap& clmap) const;
+   Bool_t  HasClAtPadRow(Int_t i) const { return fPadRawMap.TestBitNumber(i);}
+   void    Print() const;
+  protected:
+  private:
+   TBits    fPadRawMap;//bit vector of length 150 correspondind to total number of padraws in TPC
+   static const Int_t fNPadRows;
+   ClassDef(AliClusterMap,1)
+};
+
+#endif 
index e7f1eb0..32e943f 100644 (file)
@@ -30,14 +30,20 @@ AliEventCut::~AliEventCut()
 
 /*********************************************************/
 
-Bool_t AliEventCut::Pass(AliESD* esd) const
+Bool_t AliEventCut::Pass(AliAOD* aod) const
 {
   //returns kTRUE if rejected
+  if (aod == 0x0)
+   {
+     Error("Pass","Pointer to AOD is NULL. Not passed the cut");
+     return kFALSE;
+   }
+   
   TIter iter(fBaseCuts);
   AliBaseEventCut* becut;
   while (( becut = (AliBaseEventCut*)iter() ))
    {
-     if (becut->Pass(esd)) return kTRUE;
+     if (becut->Pass(aod)) return kTRUE;
    }
   return kFALSE;
 }
index 6ab0369..21e822b 100644 (file)
@@ -12,7 +12,7 @@
 
 #include "TObject.h"
 
-class AliESD;
+class AliAOD;
 class TObjArray;
 
 class AliEventCut: public TObject
@@ -21,7 +21,7 @@ class AliEventCut: public TObject
     AliEventCut();
     virtual ~AliEventCut();
     
-    virtual Bool_t Pass(AliESD* esd) const;//returns kTRUE if rejected
+    virtual Bool_t Pass(AliAOD* aod) const;//returns kTRUE if rejected
     
   protected:
     TObjArray* fBaseCuts;
index 36d029e..2d5184e 100644 (file)
 #include <TParticle.h>
 
 #include <AliStack.h>
+#include <AliAOD.h>
+#include <AliAODParticle.h>
+#include <AliAODParticleCut.h>
+
 #include <AliESDtrack.h>
 #include <AliESD.h>
 
 ClassImp(AliFlowAnalysis)
 
+AliFlowAnalysis::AliFlowAnalysis():
+ fPartCut(0x0)
+{
+ //ctor
+}
+/*********************************************************/
+
+AliFlowAnalysis::~AliFlowAnalysis()
+{
+ //dtor
+  delete fPartCut;
+}
+/*********************************************************/
+
 Int_t AliFlowAnalysis::Init()
 {
   //Initilizes anaysis
   Info("Init","");
   return 0; 
 }
-
 /*********************************************************/
 
-Int_t AliFlowAnalysis::ProcessEvent(AliESD* esd, AliStack* stack)
+Int_t AliFlowAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
 {
  
-  Info("ProcessEvent","Stack address %#x",stack);
-  Double_t psi = GetEventPlane(esd);
-  Info("ProcessEvent","Event plane is %f",psi);
+  Info("ProcessEvent","Sim AOD address %#x",aodsim);
+  Double_t psi = 0, v2 = 0;
+  if (aodrec)
+   {
+     GetFlow(aodrec,v2,psi);
+     Info("ProcessEvent","Reconstructed Event: Event plane is %f, V2 is %f",psi,v2);
+   }  
+
+  if (aodsim)
+   {
+     GetFlow(aodsim,v2,psi);
+     Info("ProcessEvent","Simulated Event: Event plane is %f, V2 is %f",psi,v2);
+   }  
+  
   return 0;
   
 }
@@ -50,12 +78,90 @@ Int_t AliFlowAnalysis::Finish()
 }
 /*********************************************************/
 
+Double_t AliFlowAnalysis::GetEventPlane(AliAOD* aod)
+{
+  //returns event plane in degrees
+  if (aod == 0x0)
+   {
+     Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL");
+     return -1.0;
+   }
+
+  Double_t psi;
+  Int_t mult = aod->GetNumberOfParticles();
+  
+  Double_t ssin = 0, scos = 0;
+
+  for (Int_t i=0; i<mult; i++) 
+   {
+     AliAODParticle* aodtrack = aod->GetParticle(i);
+     if (aodtrack == 0x0)
+      {
+        Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
+        continue;
+      }
+     
+     if (fPartCut)
+      if (fPartCut->Pass(aodtrack))
+        continue;
+
+     Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px()); 
+     
+     ssin += TMath::Sin( 2.0 * phi );
+     scos += TMath::Cos( 2.0 * phi );
+   }
+
+  psi = atan2 (ssin, scos) / 2.0;
+  psi = psi * 180. / TMath::Pi(); 
+  
+  return psi;
+
+}
+/*********************************************************/
+
+void AliFlowAnalysis::GetFlow(AliAOD* aod,Double_t& v2,Double_t& psi)
+{
+//returns flow parameters: v2 and event plane
+  if (aod == 0x0)
+   {
+     Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL");
+     return;
+   }
+   
+  psi = GetEventPlane(aod);
+  Int_t mult = aod->GetNumberOfParticles();
+  
+  Double_t ssin = 0, scos = 0;
+
+  for (Int_t i=0; i<mult; i++) 
+   {
+     AliAODParticle* aodtrack = aod->GetParticle(i);
+     if (aodtrack == 0x0)
+      {
+        Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
+        continue;
+      }
+     if (fPartCut)
+      if (fPartCut->Pass(aodtrack))
+        continue;
+      
+     Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px()); 
+     ssin += TMath::Sin( 2.0 * (phi - psi));
+     scos += TMath::Cos( 2.0 * (phi - psi));
+   }
+   
+  v2 = TMath::Hypot(ssin,scos);
+}
+
+
+/*********************************************************/
+
 Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd)
 {
   //returns event plane
   if (esd == 0x0)
    {
-     ::Error("GetFlow","Pointer to ESD is NULL");
+     ::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL");
      return -1.0;
    }
 
@@ -69,7 +175,7 @@ Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd)
      AliESDtrack* esdtrack = esd->GetTrack(i);
      if (esdtrack == 0x0)
       {
-        ::Error("GetEventPlane","Can not get track %d", i);
+        ::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
         continue;
       }
       
@@ -87,13 +193,14 @@ Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd)
   return psi;
 
 }
+/*********************************************************/
 
 void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi)
 {
 //returns flow parameters: v2 and event plane
   if (esd == 0x0)
    {
-     ::Error("GetFlow","Pointer to ESD is NULL");
+     ::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL");
      return;
    }
    
@@ -107,7 +214,7 @@ void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi)
      AliESDtrack* esdtrack = esd->GetTrack(i);
      if (esdtrack == 0x0)
       {
-        ::Error("GetEventPlane","Can not get track %d", i);
+        ::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
         continue;
       }
       
@@ -121,3 +228,4 @@ void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi)
    
   v2 = TMath::Hypot(ssin,scos);
 }
+/*********************************************************/
index 5cf4940..f20dcb3 100644 (file)
 #include "AliAnalysis.h"
 
 class AliESD;
+class AliAOD;
 class AliStack;
-
+class AliAODParticleCut;
+    
 class AliFlowAnalysis: public AliAnalysis
 { 
   public: 
-     AliFlowAnalysis(){}
-     ~AliFlowAnalysis(){}
+     AliFlowAnalysis();
+     virtual ~AliFlowAnalysis();
 
     Int_t Init();
-    Int_t ProcessEvent(AliESD* esd, AliStack* stack = 0x0);
+    Int_t ProcessEvent(AliAOD* aodrec, AliAOD* aodsim = 0x0);
     Int_t Finish();
    
+    void SetParticleCut(AliAODParticleCut* pcut){fPartCut = pcut;}
     static Double_t GetEventPlane(AliESD* esd);
     static void     GetFlow(AliESD* esd,Double_t& v2,Double_t& psi);
-  protected:
 
+    Double_t        GetEventPlane(AliAOD* aod);
+    void            GetFlow(AliAOD* aod,Double_t& v2,Double_t& psi);
+  protected:
+    
   private:
+    AliAODParticleCut* fPartCut;//Particle Cut
     ClassDef(AliFlowAnalysis,1)
 };
 
diff --git a/ANALYSIS/AliParticle.cxx b/ANALYSIS/AliParticle.cxx
new file mode 100644 (file)
index 0000000..eeced89
--- /dev/null
@@ -0,0 +1,411 @@
+#include "AliHBTParticle.h"
+//___________________________________________________________
+/////////////////////////////////////////////////////////////
+//
+// class 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
+// Additionaly adjusted for HBT Analysies purposes
+// + pointer to Track Points
+// + pointer to Cluster Map(s)
+//
+// Piotr.Skowronski@cern.ch
+//
+/////////////////////////////////////////////////////////////
+#include <TParticle.h>
+#include "AliHBTTrackPoints.h"
+#include "AliHBTClusterMap.h"
+
+ClassImp(AliHBTParticle)
+
+Int_t AliHBTParticle::fgDebug = 0;
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle():  
+ fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
+ fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0),
+ fTrackPoints(0x0),fClusterMap(0x0)
+{//empty particle
+}
+//______________________________________________________________________________
+
+AliHBTParticle::AliHBTParticle(Int_t pdg, Int_t idx,
+               Double_t px, Double_t py, Double_t pz, Double_t etot,
+               Double_t vx, Double_t vy, Double_t vz, Double_t time):  
+  fPdgIdx(0), fIdxInEvent(idx),fNPids(0),fPids(0x0),fPidProb(0x0),
+  fCalcMass(0), 
+  fPx(px), fPy(py),fPz(pz),fE(etot), 
+  fVx(vx), fVy(vy),fVz(vz),fVt(time),
+  fTrackPoints(0x0),fClusterMap(0x0)
+{
+//mormal constructor
+  SetPdgCode(pdg);
+  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(Int_t pdg, Float_t prob, Int_t idx, 
+                               Double_t px, Double_t py, Double_t pz, Double_t etot,
+                               Double_t vx, Double_t vy, Double_t vz, Double_t time):
+  fPdgIdx(0), fIdxInEvent(idx),fNPids(0),fPids(0x0),fPidProb(0x0),
+  fCalcMass(0), 
+  fPx(px), fPy(py),fPz(pz),fE(etot), 
+  fVx(vx), fVy(vy),fVz(vz),fVt(time),
+  fTrackPoints(0x0),fClusterMap(0x0)
+{
+//mormal constructor
+  SetPdgCode(pdg,prob);
+  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 AliHBTParticle& in):
+   TObject(in),
+   fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
+   fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
+   fCalcMass(in.GetCalcMass()),
+   fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.Energy()), 
+   fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T()),
+   fTrackPoints(0x0), fClusterMap(0x0)
+{
+ //Copy constructor
+ for(Int_t i = 0; i<fNPids; i++)
+  {
+    fPids[i] =  in.fPids[i];
+    fPidProb[i] = in.fPidProb[i];
+  }
+ if (in.fTrackPoints)
+   fTrackPoints = (AliHBTTrackPoints*)in.fTrackPoints->Clone();
+ if (in.fClusterMap)
+   fClusterMap = (AliHBTClusterMap*)in.fClusterMap->Clone();
+}
+
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle(const TParticle &p,Int_t idx):
+   fPdgIdx(0), fIdxInEvent(idx),
+   fNPids(0),fPids(0x0),fPidProb(0x0),
+   fCalcMass(p.GetCalcMass()),
+   fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()), 
+   fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T()),
+   fTrackPoints(0x0),fClusterMap(0x0)
+{
+ //all copied in the initialization
+ SetPdgCode(p.GetPdgCode());
+}
+//______________________________________________________________________________
+
+AliHBTParticle::~AliHBTParticle()
+{
+//dtor  
+  delete [] fPids;
+  delete [] fPidProb;
+  delete fTrackPoints;
+  delete fClusterMap;
+}
+//______________________________________________________________________________
+
+AliHBTParticle& AliHBTParticle::operator=(const AliHBTParticle& in)
+{
+//assigment operator
+  
+  fNPids = in.fNPids;
+  delete [] fPids;
+  delete [] fPidProb;
+  Int_t* fPids = new Int_t[fNPids];
+  Float_t* fPidProb = new Float_t[fNPids];
+  for (Int_t i = 0; i < fNPids;i++)
+   {
+     fPids[i]    = in.fPids[i];
+     fPidProb[i] = in.fPidProb[i];
+   }
+   
+  fPdgIdx = in.fPdgIdx; 
+  fIdxInEvent = in.fIdxInEvent;  
+  fCalcMass = in.GetCalcMass();
+  fPx = in.Px();
+  fPy = in.Py();
+  fPz = in.Pz();
+  fE = in.Energy(); 
+  fVx = in.Vx();
+  fVy = in.Vy();
+  fVz = in.Vz();
+  fVt = in.T();
+  
+  delete fTrackPoints;
+  fTrackPoints = (in.fTrackPoints)?(AliHBTTrackPoints*)fTrackPoints->Clone():0x0;
+  
+  delete fClusterMap;
+  fClusterMap =  (in.fClusterMap)?(AliHBTClusterMap*)in.fClusterMap->Clone():0x0;
+  
+  return *this;
+}
+//______________________________________________________________________________
+
+void AliHBTParticle::SetPdgCode(Int_t pdg,Float_t prob)
+{
+  SetPIDprobability(pdg,prob);
+  fPdgIdx = GetPidSlot(pdg);
+}
+
+//______________________________________________________________________________
+void AliHBTParticle::SetPIDprobability(Int_t pdg, Float_t prob)
+{
+//Sets another pdg code and corresponding probabilty
+//Ids are set in decreasing order
+//Check if total prbaility is not ivercoming unity is performed
+//in case, warning is printed
+  if (fgDebug > 9) Info("SetPIDprobability","Setting PID %d prob %f",pdg,prob);
+
+  Float_t totprob = 0.0;//sums up probabilities
+  Int_t idx = GetPidSlot(pdg);
+  Int_t i;
+  if (idx > -1) 
+   {
+     fPidProb[idx] = prob;
+     for (i = 0; i < fNPids;i++) totprob+=fPidProb[i];
+     if (totprob > (1.0+0.000001))
+       {
+         Warning("SetPIDprobability","Total probability greater than unity (%f)",totprob);
+       }
+     if (fgDebug > 9) 
+      {
+        Info("SetPIDprobability","Current Total probability: %f",totprob);
+      }
+     return;
+   }
+    
+  Int_t currentpid = GetPdgCode();
+  fNPids++;
+  Float_t* aPidProbNew = new Float_t[fNPids];
+  Int_t* aPidsNew = new Int_t[fNPids];
+  
+  for (i = 0; i < fNPids-1;i++)//find a slot
+   {
+     if ( fPidProb[i] > prob)
+      {
+        if (fgDebug>9) Info("SetPID","Copying entry %d",i);
+        aPidProbNew[i] = fPidProb[i];
+        aPidsNew[i] = fPids[i];
+        totprob+=fPidProb[i];
+      }
+     else break;
+   }
+
+  if (fgDebug > 9) Info("SetPID","Setting new PID on entry %d",i);
+  aPidProbNew[i] = prob;
+  aPidsNew[i] = pdg;
+  totprob+=prob;
+  
+
+  for (Int_t j = fNPids-1; j > i ;j--)//copy rest of old araays 
+   {
+     if (fgDebug > 9) Info("SetPID","Copying from old entry %d to new entry %d",j-1,j);
+     aPidProbNew[j] = fPidProb[j-1];
+     aPidsNew[j] = fPids[j-1];
+     totprob+=fPidProb[j-1];
+   }
+
+  delete [] fPidProb;
+  delete [] fPids;
+  
+  fPidProb = aPidProbNew;
+  fPids = aPidsNew;
+  
+  fPdgIdx = GetPidSlot(currentpid);
+  if (fPdgIdx == -1) fPdgIdx = 0;
+  
+  if (totprob > (1.0+0.000001))//place for numerical error
+   {
+     Warning("SetId","Total probability is greater than unity (%f)!!!",totprob);
+     Print();
+   }
+}
+//______________________________________________________________________________
+
+Float_t AliHBTParticle::GetPIDprobability(Int_t pdg) const
+{
+//Returns probability that this particle is the type of pdg
+  Int_t idx = GetPidSlot(pdg);
+  if (idx < 0) return 0.0;//such pid was not specified for this particle
+  return fPidProb[idx];
+}
+//______________________________________________________________________________
+
+const Char_t* AliHBTParticle::GetName() const 
+{
+  //returns name of this particle 
+   static char def[4] = "XXX";
+   const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
+   if (ap) return ap->GetName();
+   else    return def;
+}
+//______________________________________________________________________________
+
+Int_t AliHBTParticle::GetPidSlot(Int_t pdg) const
+{
+ //returns position of the given PID in fPids (and fPidProb) array.
+ if (fPids == 0x0) return -1;
+ for (Int_t i = 0; i< fNPids; i++)
+  {
+   if (fPids[i] == pdg) return i;
+  }
+ return -1;
+}
+//______________________________________________________________________________
+
+Int_t AliHBTParticle::GetNthPid(Int_t idx) const
+{
+  //returns PID sitting on slot idx in fPids
+  if ( (idx < 0) || (idx >= fNPids) )
+   {
+     Error("GetNthPid","Out Of Bounds");
+     return 0;
+   }
+  return fPids[idx];
+}
+//______________________________________________________________________________
+
+Float_t AliHBTParticle::GetNthPidProb(Int_t idx) const
+{
+  //returns PID sitting on slot idx in fPidProb
+  if ( (idx < 0) || (idx >= fNPids) )
+   {
+     Error("GetNthPid","Out Of Bounds");
+     return 0;
+   }
+  return fPidProb[idx];
+}
+//______________________________________________________________________________
+
+void AliHBTParticle::Print() const
+{
+//prints information about particle
+  printf("____________________________________________________\n");
+  printf("Idx: %d  PID: %d  Name: ",fIdxInEvent,GetPdgCode());
+  TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
+  if (pdgp)
+   {
+     printf("%s Mass: %f\n",pdgp->GetName(),pdgp->Mass());
+   }
+  else
+   {
+     printf("Not known\n");
+   }
+  
+  printf("Px: %+f Py: %+f Pz: %+f E: %+f Calculated Mass: %f\nVx: %+f Vy: %+f Vz: %+f T: %+f\n",
+          Px(),Py(),Pz(),Energy(),GetCalcMass(),Vx(),Vy(),Vz(),T());
+
+  for (Int_t i = 0; i < fNPids; i++)
+   {
+     printf("# %d  PID: %d  Probability %f name ",i,fPids[i],fPidProb[i]);
+     const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPids[i]);
+     if (ap)
+      {
+        printf("%s Mass %f\n",ap->GetName(),ap->Mass());
+      }
+     else
+      {
+        printf("Not known\n");
+      }
+   }
+}
+
+//______________________________________________________________________________
+
+//void AliHBTParticle::Streamer(TBuffer &b)
+//{
+//     // Stream all objects in the array to or from the I/O buffer.
+//   UInt_t R__s, R__c;
+//   Int_t i;
+//   if (b.IsReading()) 
+//    {
+//      delete [] fPids;
+//      delete [] fPidProb;
+//      
+//      Version_t v = b.ReadVersion(&R__s, &R__c);
+//      if (v == 1)
+//       {
+//         AliHBTParticle::Class()->ReadBuffer(b, this);
+//      }
+//      else
+//       {
+//        TObject::Streamer(b);
+//       b >> fPdgIdx;
+//       b >> fIdxInEvent;
+//       
+//       b >> fNPids;
+//       Int_t* fPids = new Int_t[fNPids];
+//        Float_t* fPidProb = new Float_t[fNPids];
+//        for (i = 0;i<fNPids;i++) 
+//         {
+//           b >> fPids[i];
+//         }
+//        for (i = 0;i<fNPids;i++) 
+//        {
+//          b >> fPidProb[i];
+//         }
+//        b >> fCalcMass;
+//
+//        b >> fPx;
+//       b >> fPy;
+//       b >> fPz;
+//        b >> fE;
+//
+//       b >> fVx;
+//        b >> fVy;
+//        b >> fVz;
+//       b >> fVt;
+//       Info("Streamer","Read data");
+//        Print();
+//       }
+//
+//      b.CheckByteCount(R__s, R__c,AliHBTParticle::IsA());
+//    } 
+//  else 
+//   {
+//     R__c = b.WriteVersion(AliHBTParticle::IsA(), kTRUE);
+//     TObject::Streamer(b);
+//     Info("Streamer","Read data");
+//     Print();
+//
+//     b << fPdgIdx;
+//     b << fIdxInEvent;
+//     b << fNPids;
+//     for (i = 0;i<fNPids;i++) 
+//      {
+//        b << fPids[i];
+//      }
+//      {
+//      {
+//     for (i = 0;i<fNPids;i++) 
+//     {
+//        b << fPidProb[i];
+//      }
+//     b << fCalcMass;
+//
+//     b << fPx;
+//     b << fPy;
+//     b << fPz;
+//     b << fE;
+//
+//     b << fVx;
+//     b << fVy;
+//     b << fVz;
+//     b << fVt;
+//
+//    b.SetByteCount(R__c, kTRUE);
+//   }
+//}
diff --git a/ANALYSIS/AliParticle.h b/ANALYSIS/AliParticle.h
new file mode 100644 (file)
index 0000000..cd126ea
--- /dev/null
@@ -0,0 +1,156 @@
+#ifndef ALIHBTPARTICLE_H
+#define ALIHBTPARTICLE_H
+//___________________________________________________________
+/////////////////////////////////////////////////////////////
+//
+// class 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
+// Additionaly adjusted for HBT Analysies purposes
+// + pointer to Track Points
+// + pointer to Cluster Map(s)
+//
+// Piotr.Skowronski@cern.ch
+//
+/////////////////////////////////////////////////////////////
+
+#include <TObject.h>
+#include <TLorentzVector.h>
+#include <TMath.h>
+#include <TDatabasePDG.h>
+
+
+class TParticle;
+class AliHBTTrackPoints;
+class AliHBTClusterMap;
+
+class AliHBTParticle : public TObject
+{
+public:
+                                // ****** constructors and destructor
+  AliHBTParticle();
+  AliHBTParticle(const AliHBTParticle& in); 
+  AliHBTParticle(Int_t pdg, Int_t idx, 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(Int_t pdg, Float_t prob, Int_t idx, 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,Int_t idx);
+
+  virtual ~AliHBTParticle();
+  
+  AliHBTParticle& operator=(const AliHBTParticle& in); 
+  
+  void           SetPIDprobability(Int_t pdg, Float_t prob = 1.0);
+  Float_t        GetPIDprobability(Int_t pdg) const;
+  
+  Int_t          GetPdgCode      () const { return (fPids)?fPids[fPdgIdx]:0;}
+  Int_t          GetPid          () const { return GetPdgCode();}
+  Float_t        GetPidProb      () const { return (fPidProb)?fPidProb[fPdgIdx]:0;}
+  
+  Int_t          GetUID          () const { return fIdxInEvent;}
+  Int_t          GetNumberOfPids () const { return fNPids;}
+  Int_t          GetNthPid         (Int_t idx) const;
+  Float_t        GetNthPidProb     (Int_t idx) const;
+      
+  void           SetPdgCode(Int_t pdg, Float_t prob = 1.0);
+  Double_t       GetCalcMass     () const { return fCalcMass; }
+  Double_t       GetMass         ()       { return (GetPDG())?GetPDG()->Mass():-1.;}
+
+  TParticlePDG*  GetPDG          (){return TDatabasePDG::Instance()->GetParticle(GetPdgCode());}
+
+  Int_t          Beauty          ()  { return GetPDG()->Beauty(); }
+  Int_t          Charm           ()  { return GetPDG()->Charm(); }
+  Int_t          Strangeness     ()  { return GetPDG()->Strangeness();}
+  void ProductionVertex(TLorentzVector &v) const { 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) const { 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 TMath::Pi()+TMath::ATan2(-fPy,-fPx); }
+
+  Double_t         Theta () const { return (fPz==0)?TMath::PiOver2():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());}
+  void           SetCalcMass(Double_t mass) {fCalcMass = mass;}
+  
+  void           SetUID(Int_t id){fIdxInEvent = id;}
+  
+  const Char_t*  GetName() const; 
+  void           Print() const;
+  
+  void           SetTrackPoints(AliHBTTrackPoints* tpts){fTrackPoints = tpts;}
+  AliHBTTrackPoints* GetTrackPoints() const {return fTrackPoints;}
+  void           SetClusterMap(AliHBTClusterMap* cm){fClusterMap = cm;}
+  AliHBTClusterMap* GetClusterMap() const {return fClusterMap;}
+  
+  static void    SetDebug(Int_t dbg=1){fgDebug=dbg;}
+  static Int_t   GetDebug(){return fgDebug;}
+  
+protected:
+  Int_t          GetPidSlot(Int_t pdg) const;//returns position of the given PID in fPids (and fPidProb) array.
+
+private:
+  Char_t         fPdgIdx;               // index of PDG code of the particle in fPids
+  Int_t          fIdxInEvent;           // index of a particle: the same particle can appear in the event
+                                        //  many times with different pid's. Idx allows to check that they are the same particles
+  Int_t          fNPids;                // number of non-zero proboble Pids
+  Int_t         *fPids;                 // [fNPids] Array with PIDs
+  Float_t       *fPidProb;              // [fNPids] PIDs probabilities
+  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
+
+  AliHBTTrackPoints* fTrackPoints;      // track positions along trajectory - used by anti-merging cut
+  AliHBTClusterMap*  fClusterMap;       // bit map of cluters occupation; 1 if has cluter on given layer/padrow/...
+    
+  static Int_t   fgDebug; //debug printout level
+  ClassDef(AliHBTParticle,3)  // TParticle vertex particle information
+};
+
+#endif
diff --git a/ANALYSIS/AliReader.cxx b/ANALYSIS/AliReader.cxx
new file mode 100644 (file)
index 0000000..1b5743c
--- /dev/null
@@ -0,0 +1,434 @@
+#include "AliReader.h"
+//_________________________________________________________________________
+///////////////////////////////////////////////////////////////////////////
+//
+// class AliReader
+//
+// Reader Base class (reads particles and tracks and
+// puts it to the AliAODRun objects
+//
+// Provides functionality for both buffering and non-buffering reading
+// This can be switched on/off via method SetEventBuffering(bool)
+// The main method that inheriting classes need to implement is ReadNext()
+// that read next event in queue.
+// The others are:
+// Bool_t  ReadsRec() const; specifies if reader is able to read simulated particles
+// Bool_t  ReadsSim() const; specifies if reader is able to read reconstructed tracks
+// void    Rewind(); rewind reading to the beginning
+//
+// reading of next event is triggered via method Next()
+//
+// This class provides full functionality for reading from many sources
+// User can provide TObjArray of TObjString (SetDirs method or via parameter 
+// in constructor) which desribes paths of directories to search data in.
+// If none specified current directory is searched.
+//
+// Piotr.Skowronski@cern.ch
+///////////////////////////////////////////////////////////////////////////
+
+#include <TString.h>
+#include <TObjString.h>
+#include <TObjArray.h>
+#include <TClass.h>
+#include <TRandom.h>
+#include <TH1.h>
+
+#include "AliAODParticleCut.h"
+#include "AliAOD.h"
+#include "AliAODRun.h"
+ClassImp(AliReader)
+//pure virtual
+    
+/*************************************************************************************/
+
+AliReader::AliReader():
+ fCuts(new TObjArray()),
+ fDirs(0x0),
+ fCurrentEvent(0),
+ fCurrentDir(0),
+ fNEventsRead(0),
+ fEventRec(0x0),
+ fEventSim(0x0),
+ fRunSim(0x0),
+ fRunRec(0x0),
+ fIsRead(kFALSE),
+ fBufferEvents(kFALSE),
+ fBlend(kFALSE),
+ fFirst(0),
+ fLast(0),
+ fTrackCounter(0x0)
+{
+//constructor
+}
+/*************************************************************************************/
+
+AliReader::AliReader(TObjArray* dirs):
+ fCuts(new TObjArray()),
+ fDirs(dirs),
+ fCurrentEvent(0),
+ fCurrentDir(0),
+ fNEventsRead(0),
+ fEventRec(0x0),
+ fEventSim(0x0),
+ fRunSim(0x0),
+ fRunRec(0x0),
+ fIsRead(kFALSE),
+ fBufferEvents(kFALSE),
+ fBlend(kFALSE),
+ fFirst(0),
+ fLast(0),
+ fTrackCounter(0x0)
+{
+//ctor with array of directories to read as parameter
+}
+/*************************************************************************************/
+AliReader::AliReader(const AliReader& in):
+ TNamed(in),
+ fCuts((in.fCuts)?(TObjArray*)in.fCuts->Clone():0x0),
+ fDirs((in.fDirs)?(TObjArray*)in.fDirs->Clone():0x0),
+ fCurrentEvent(0),
+ fCurrentDir(0),
+ fNEventsRead(0),
+ fEventRec(0x0),
+ fEventSim(0x0),
+ fRunSim(0x0),
+ fRunRec(0x0),
+ fIsRead(kFALSE),
+ fBufferEvents(in.fBufferEvents),
+ fBlend(in.fBlend),
+ fFirst(in.fFirst),
+ fLast(in.fLast),
+ fTrackCounter(0x0)
+{
+ //cpy constructor
+}
+
+AliReader::~AliReader()
+{
+//destructor
+ if(fCuts)
+  {
+   fCuts->SetOwner();
+   delete fCuts;
+  }
+ delete fEventSim;
+ delete fEventRec;
+ delete fTrackCounter;
+}
+/*************************************************************************************/
+
+AliReader& AliReader::operator=(const AliReader& in)
+{
+  //Assigment operator
+ if (this == &in) return *this;  
+ TNamed::operator=( (const TNamed&)in );
+  
+ fCuts = (in.fCuts)?(TObjArray*)in.fCuts->Clone():0x0;
+ fDirs = (in.fDirs)?(TObjArray*)in.fDirs->Clone():0x0;
+ fCurrentEvent = 0;
+ fCurrentDir = 0;
+ fNEventsRead = 0;
+ fEventRec = 0x0;
+ fEventSim = 0x0;
+ fRunSim = 0x0;
+ fRunRec = 0x0;
+ fIsRead = kFALSE;
+ fBufferEvents = in.fBufferEvents;
+ fBlend = in.fBlend;
+ fFirst = in.fFirst;
+ fLast = in.fLast;
+ fTrackCounter = 0x0;
+ return *this;  
+}
+/*************************************************************************************/
+
+Int_t AliReader::Next()
+{
+//moves to next event
+
+  //if asked to read up to event nb. fLast, and it is overcome, report no more events
+  if ((fNEventsRead > fLast) && (fLast > 0) ) return kTRUE;
+  
+  if (fTrackCounter == 0x0)//create Track Counter
+   {
+     fTrackCounter = new TH1I("trackcounter","Track Counter",20000,0,20000);
+     fTrackCounter->SetDirectory(0x0);
+   }
+  
+  do //if asked to read from event fFirst, rewind to it
+   {
+    if ( ReadNext() == kTRUE) //if no more evets, return it
+      return kTRUE;
+   }while (fNEventsRead < fFirst);
+   
+  //here we have event
+  
+  if (fBlend) Blend();//Mix particles order 
+  
+  if (fBufferEvents)//store events if buffering is on
+   {
+     if ( ReadsRec() && fEventRec) 
+       fRunRec->SetEvent(fNEventsRead-1-fFirst,fEventRec);
+     if ( ReadsSim() && fEventSim)
+       fRunSim->SetEvent(fNEventsRead-1-fFirst,fEventSim);
+   }
+  return kFALSE;
+}
+/*************************************************************************************/
+
+void AliReader::AddParticleCut(AliAODParticleCut* cut)
+{
+ //sets the new cut. MAKES A COPY OF THE 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;
+   }
+  AliAODParticleCut *c = (AliAODParticleCut*)cut->Clone();
+  fCuts->Add(c);
+}
+/********************************************************************/
+
+AliAOD* AliReader::GetEventSim(Int_t n)
+ {
+ //returns Nth event with simulated particles
+  if (ReadsSim() == kFALSE)
+   {
+     Error("GetParticleEvent","This reader is not able to provide simulated particles.");
+     return 0;
+   } 
+   
+  if (!fIsRead) 
+   { 
+    if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
+    if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
+    
+    if (Read(fRunSim,fRunRec))
+     {
+       Error("GetParticleEvent","Error in reading");
+       return 0x0;
+     }
+    else fIsRead = kTRUE;
+   }
+  return fRunSim->GetEvent(n);
+ }
+/********************************************************************/
+
+AliAOD* AliReader::GetEventRec(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+  if (ReadsRec() == kFALSE)
+   {
+     Error("GetTrackEvent","This reader is not able to provide recosntructed tracks.");
+     return 0;
+   } 
+  if (!fIsRead) 
+   {
+    if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
+    if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
+    
+    if(Read(fRunSim,fRunRec))
+     {
+       Error("GetTrackEvent","Error in reading");
+       return 0x0;
+     }
+    else fIsRead = kTRUE;
+   }
+  return fRunRec->GetEvent(n);
+ }
+/********************************************************************/
+
+Int_t AliReader::GetNumberOfSimEvents()
+ {
+ //returns number of events of particles
+  if (ReadsSim() == kFALSE)
+   {
+     Error("GetNumberOfPartEvents","This reader is not able to provide simulated particles.");
+     return 0;
+   } 
+   
+  if (!fIsRead) 
+   {
+    if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
+    if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
+    
+    if (Read(fRunSim,fRunRec))
+     {
+       Error("GetNumberOfPartEvents","Error in reading");
+       return 0;
+     }
+    else fIsRead = kTRUE;
+   }
+   return fRunSim->GetNumberOfEvents();
+ }
+/********************************************************************/
+Int_t AliReader::GetNumberOfRecEvents()
+ {
+ //returns number of events of tracks
+  if (ReadsRec() == kFALSE)
+   {
+     Error("GetNumberOfTrackEvents","This reader is not able to provide recosntructed tracks.");
+     return 0;
+   } 
+  if (!fIsRead)
+   {
+     if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
+     if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
+     
+     if(Read(fRunSim,fRunRec))
+      {
+        Error("GetNumberOfTrackEvents","Error in reading");
+        return 0;
+      }
+     else fIsRead = kTRUE;
+   }
+  return fRunRec->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+Int_t AliReader::Read(AliAODRun* particles, AliAODRun *tracks)
+{
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+  Info("Read","");
+  
+  if ( ReadsSim() && (particles == 0x0) ) //check if an object is instatiated
+   {
+     Error("Read"," particles object must be instatiated before passing it to the reader");
+     return 1;
+   }
+  if ( ReadsRec() && (tracks == 0x0) )  //check if an object is instatiated
+   {
+     Error("Read"," tracks object must be instatiated before passing it to the reader");
+     return 1;
+   }
+   
+  if (ReadsSim()) particles->Reset();//clear runs == delete all old events
+  if (ReadsRec()) tracks->Reset();
+  
+  Rewind();
+  
+  Int_t i = 0;
+  while(Next() == kFALSE)
+   {
+     if (ReadsRec()) tracks->SetEvent(i,fEventRec);
+     if (ReadsSim()) particles->SetEvent(i,fEventSim);
+     i++;
+   }
+  return 0;
+}      
+/*************************************************************************************/
+
+Bool_t AliReader::Pass(AliAODParticle* 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
+  if ( fCuts->GetEntriesFast() == 0 ) return kFALSE; //if no cut specified accept all particles
+  for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)   
+   {
+     AliAODParticleCut &cut = *((AliAODParticleCut*)fCuts->At(i));
+     if(!cut.Pass(p)) return kFALSE;  //accepted
+   }
+   
+  return kTRUE;//not accepted
+}
+/*************************************************************************************/
+
+Bool_t  AliReader::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;
+
+ if ( fCuts->GetEntriesFast() == 0 ) return kFALSE; //if no cut specified accept all particles
+  
+ for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)   
+   {
+     AliAODParticleCut &cut = *((AliAODParticleCut*)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;
+}
+/*************************************************************************************/
+
+TString& AliReader::GetDirName(Int_t entry)
+{
+//returns directory name of next one to read
+  TString* retval;//return value
+  if (fDirs ==  0x0)
+   {
+     retval = new TString(".");
+     return *retval;
+   }
+
+  if ( (entry>fDirs->GetEntries()) || (entry<0))//if out of bounds return empty string
+   {                                            //note that entry==0 is accepted even if array is empty (size=0)
+     Error("GetDirName","Name out of bounds");
+     retval = new TString();
+     return *retval;
+   }
+
+  if (fDirs->GetEntries() == 0)
+   { 
+     retval = new TString(".");
+     return *retval;
+   }
+
+  TClass *objclass = fDirs->At(entry)->IsA();
+  TClass *stringclass = TObjString::Class();
+
+  TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry));
+
+  if(dir == 0x0)
+   {
+     Error("GetDirName","Object in TObjArray is not a TObjString or its descendant");
+     retval = new TString();
+     return *retval;
+   }
+  if (gDebug > 0) Info("GetDirName","Returned ok %s",dir->String().Data());
+  return dir->String();
+}
+/*************************************************************************************/
+
+void AliReader::Blend()
+{
+  //randomly change positions of the particles after reading
+  //is used to check if some distributions (of many particle properties) 
+  //depend on the order of particles
+  //(tracking gives particles Pt sorted)
+  
+  if (fEventSim == 0x0) return;
+  
+  for (Int_t i = 2; i < fEventSim->GetNumberOfParticles(); i++)
+   {
+     Int_t with = gRandom->Integer(i);
+     fEventSim->SwapParticles(i,with);
+     if (fEventRec) fEventRec->SwapParticles(i,with);
+   }
+}
+/*************************************************************************************/
+
+void AliReader::WriteTrackCounter() const
+{
+ //writes the counter histogram
+ if (fTrackCounter) fTrackCounter->Write(0,TObject::kOverwrite);
+ else 
+  {
+    Warning("WriteTrackCounter","Counter is NULL");
+  }
+}
diff --git a/ANALYSIS/AliReader.h b/ANALYSIS/AliReader.h
new file mode 100644 (file)
index 0000000..bf891a4
--- /dev/null
@@ -0,0 +1,101 @@
+#ifndef ALIREADER_H
+#define ALIREADER_H
+//_________________________________________________________________________
+///////////////////////////////////////////////////////////////////////////
+//
+// class AliReader
+//
+// Reader Base class 
+// Reads particles and tracks and
+// puts it to the AliAOD objects and eventuall buffers in AliAODRuns
+//
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////////////////////
+
+#include <TNamed.h>
+#include <TObjArray.h>
+
+class AliAODRun;
+class AliAOD;
+class AliAODParticleCut;
+class AliAODParticle;
+class TString;
+class TH1I;
+class AliReader: public TNamed
+{
+  public:
+    AliReader();
+    AliReader(TObjArray*);
+    AliReader(const AliReader& in);
+    virtual ~AliReader();
+    
+    AliReader& operator=(const AliReader& in);
+        
+    virtual Int_t        Next();
+    virtual void         Rewind() = 0; //
+    
+    virtual Bool_t       ReadsSim() const = 0; //specifies if reader is able to read simulated particles
+    virtual Bool_t       ReadsRec() const = 0;//specifies if reader is able to read reconstructed tracks
+    
+    void                 AddParticleCut(AliAODParticleCut* cut);
+    
+    virtual Int_t        Read(AliAODRun* particles, AliAODRun *tracks);
+    
+    virtual AliAOD*      GetEventRec() const {return fEventRec;}//
+    virtual AliAOD*      GetEventSim() const {return fEventSim;}//can not be const because position randomizer overloads it
+    
+    virtual AliAOD*      GetEventRec(Int_t n);
+    virtual AliAOD*      GetEventSim(Int_t n);
+    
+    virtual Int_t        GetNumberOfRecEvents();
+    virtual Int_t        GetNumberOfSimEvents();
+    
+    void                 SetDirs(TObjArray* dirs){fDirs = dirs;} //sets array directories names
+    void                 SetEventBuffering(Bool_t flag){fBufferEvents = flag;}//switches on/off buffering - read data are kept in local buffer
+    void                 SetBlend(Bool_t flag = kTRUE){fBlend=flag;} //set blending - randomizing particle order
+    virtual Int_t        GetNumberOfDirs() const {return (fDirs)?fDirs->GetEntries():0;}
+    void                 ReadEventsFromTo(Int_t first,Int_t last){fFirst = first; fLast = last;}
+    virtual TH1I*        GetTrackCounter() const {return fTrackCounter;}
+    virtual void         WriteTrackCounter() const;
+    
+  protected:
+    
+    TObjArray*           fCuts;//array with particle cuts
+    TObjArray*           fDirs;//arry with directories to read data from
+    
+    Int_t                fCurrentEvent;//!  number of current event in current directory
+    Int_t                fCurrentDir;//! number of current directory
+    
+    Int_t                fNEventsRead;//!total 
+        
+    AliAOD*              fEventRec;    //! tracks read from current event
+    AliAOD*              fEventSim;    //! particles read from current event
+    
+    AliAODRun*           fRunSim; //!simulated particles
+    AliAODRun*           fRunRec; //!reconstructed tracks
+    
+    Bool_t               fIsRead;//!flag indicating if the data are already read
+    Bool_t               fBufferEvents;//flag indicating if the data should be bufferred
+    
+    Bool_t               fBlend;// flag indicating if randomly change positions of the particles after reading
+    
+    Int_t                fFirst;//first event to return (all are before are skipped)
+    Int_t                fLast;//last
+
+    TH1I*                fTrackCounter; //histogram with number of tracks read
+    
+    virtual Int_t        ReadNext() = 0; //this methods reads next event and put result in fTracksEvent and/or fParticlesEvent
+    Bool_t               Pass(AliAODParticle* p);
+    Bool_t               Pass(Int_t pid);
+    void                 Blend();
+    
+    TString&             GetDirName(Int_t entry);
+    
+  private:
+  
+    ClassDef(AliReader,1)//
+};
+
+#endif
diff --git a/ANALYSIS/AliReaderESD.cxx b/ANALYSIS/AliReaderESD.cxx
new file mode 100644 (file)
index 0000000..4d2e75b
--- /dev/null
@@ -0,0 +1,711 @@
+#include "AliReaderESD.h"
+//____________________________________________________________________
+//////////////////////////////////////////////////////////////////////
+//                                                                  //
+// class AliReaderESD                                               //
+//                                                                  //
+// reader for ALICE Event Summary Data (ESD).                       //
+//                                                                  //
+// Piotr.Skowronski@cern.ch                                         //
+//                                                                  //
+//////////////////////////////////////////////////////////////////////
+
+#include <TPDGCode.h>
+#include <TString.h>
+#include <TObjString.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TKey.h>
+#include <TParticle.h>
+#include <TH1.h>
+
+#include <AliRun.h>
+#include <AliRunLoader.h>
+#include <AliStack.h>
+#include <AliESDtrack.h>
+#include <AliESD.h>
+
+#include "AliAnalysis.h"
+#include "AliAODRun.h"
+#include "AliAOD.h"
+#include "AliAODStdParticle.h"
+#include "AliAODParticleCut.h"
+#include "AliTrackPoints.h"
+#include "AliClusterMap.h"
+
+ClassImp(AliReaderESD)
+
+AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
+ fESDFileName(esdfilename),
+ fGAlFileName(galfilename),
+ fFile(0x0),
+ fRunLoader(0x0),
+ fKeyIterator(0x0),
+ fReadSim(kFALSE),
+ fCheckParticlePID(kFALSE),
+ fReadMostProbableOnly(kFALSE),
+ fNTrackPoints(0),
+ fdR(0.0),
+ fClusterMap(kFALSE),
+ fNTPCClustMin(0),
+ fNTPCClustMax(150),
+ fTPCChi2PerClustMin(0.0),
+ fTPCChi2PerClustMax(10e5),
+ fChi2Min(0.0),
+ fChi2Max(10e5),
+ fC00Min(0.0),
+ fC00Max(10e5),
+ fC11Min(0.0),
+ fC11Max(10e5),
+ fC22Min(0.0),
+ fC22Max(10e5),
+ fC33Min(0.0),
+ fC33Max(10e5),
+ fC44Min(0.0),
+ fC44Max(10e5),
+ fTPCC00Min(0.0),
+ fTPCC00Max(10e5),
+ fTPCC11Min(0.0),
+ fTPCC11Max(10e5),
+ fTPCC22Min(0.0),
+ fTPCC22Max(10e5),
+ fTPCC33Min(0.0),
+ fTPCC33Max(10e5),
+ fTPCC44Min(0.0),
+ fTPCC44Max(10e5)
+
+{
+  //cosntructor
+  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
+    Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
+}
+/********************************************************************/
+  
+AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
+ AliReader(dirs), 
+ fESDFileName(esdfilename),
+ fGAlFileName(galfilename),
+ fFile(0x0),
+ fRunLoader(0x0),
+ fKeyIterator(0x0),
+ fReadSim(kFALSE),
+ fCheckParticlePID(kFALSE),
+ fReadMostProbableOnly(kFALSE),
+ fNTrackPoints(0),
+ fdR(0.0),
+ fClusterMap(kFALSE),
+ fNTPCClustMin(0),
+ fNTPCClustMax(150),
+ fTPCChi2PerClustMin(0.0),
+ fTPCChi2PerClustMax(10e5),
+ fChi2Min(0.0),
+ fChi2Max(10e5),
+ fC00Min(0.0),
+ fC00Max(10e5),
+ fC11Min(0.0),
+ fC11Max(10e5),
+ fC22Min(0.0),
+ fC22Max(10e5),
+ fC33Min(0.0),
+ fC33Max(10e5),
+ fC44Min(0.0),
+ fC44Max(10e5),
+ fTPCC00Min(0.0),
+ fTPCC00Max(10e5),
+ fTPCC11Min(0.0),
+ fTPCC11Max(10e5),
+ fTPCC22Min(0.0),
+ fTPCC22Max(10e5),
+ fTPCC33Min(0.0),
+ fTPCC33Max(10e5),
+ fTPCC44Min(0.0),
+ fTPCC44Max(10e5)
+{
+  //cosntructor
+  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
+    Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
+}
+/********************************************************************/
+
+AliReaderESD::~AliReaderESD()
+{
+ //desctructor
+  delete fRunLoader;
+  delete fKeyIterator;
+  delete fFile;
+}
+/**********************************************************/
+Int_t AliReaderESD::ReadNext()
+{
+//reads next event from fFile
+  //fRunLoader is for reading Kine
+  
+  if (AliAODParticle::GetDebug())
+    Info("ReadNext","Entered");
+    
+  if (fEventSim == 0x0)  fEventSim = new AliAOD();
+  if (fEventRec == 0x0)  fEventRec = new AliAOD();
+  
+  fEventSim->Reset();
+  fEventRec->Reset();
+        
+  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
+   {
+     if (fFile == 0x0)
+      {
+       fFile = OpenFile(fCurrentDir);//rl is opened here
+       if (fFile == 0x0)
+         {
+           Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
+           fCurrentDir++;
+           continue;
+         }
+       fCurrentEvent = 0;
+       fKeyIterator = new TIter(fFile->GetListOfKeys());  
+//       fFile->Dump();
+//       fFile->GetListOfKeys()->Print();
+      } 
+     TKey* key = (TKey*)fKeyIterator->Next();
+     if (key == 0x0)
+      {
+        if (AliAODParticle::GetDebug() > 2 )
+          {
+            Info("ReadNext","No more keys.");
+          }
+        fCurrentDir++;
+        delete fKeyIterator;
+        fKeyIterator = 0x0;
+        delete fFile;//we have to assume there is no more ESD objects in the fFile
+        fFile = 0x0;
+        delete fRunLoader;
+        fRunLoader = 0x0;
+        continue;
+      }
+     //try to read
+     
+     
+//     TObject* esdobj = key->ReadObj();
+//     if (esdobj == 0x0)
+//      {
+//        if (AliAODParticle::GetDebug() > 2 )
+//          {
+//            Info("ReadNext","Key read NULL. Key Name is %s",key->GetName());
+//            key->Dump();
+//          }
+//        continue;
+//      }
+//     esdobj->Dump();
+//     AliESD* esd = dynamic_cast<AliESD*>(esdobj);
+     
+     TString esdname = "ESD";
+     esdname+=fCurrentEvent;
+     AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
+     if (esd == 0x0)
+      {
+//        if (AliAODParticle::GetDebug() > 2 )
+//          {
+//            Info("ReadNext","This key is not an AliESD object %s",key->GetName());
+//          }
+        if (AliAODParticle::GetDebug() > 2 )
+          {
+            Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
+          }
+        fCurrentDir++;
+        delete fKeyIterator;
+        fKeyIterator = 0x0;
+        delete fFile;//we have to assume there is no more ESD objects in the fFile
+        fFile = 0x0;
+        delete fRunLoader;
+        fRunLoader = 0x0;
+        continue;
+      }
+     
+     ReadESD(esd);
+      
+     fCurrentEvent++;
+     fNEventsRead++;
+     delete esd;
+     return 0;//success -> read one event
+   }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
+   
+  return 1; //no more directories to read
+}
+/**********************************************************/
+
+Int_t AliReaderESD::ReadESD(AliESD* esd)
+{
+  //****** Tentative particle type "concentrations"
+  static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
+  
+  Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
+  Double_t w[kNSpecies];
+  Double_t mom[3];//momentum
+  Double_t pos[3];//position
+  Double_t vertexpos[3];//vertex position
+  //Reads one ESD
+  if (esd == 0x0)
+   {
+     Error("ReadESD","ESD is NULL");
+     return 1;
+   }
+
+  TDatabasePDG* pdgdb = TDatabasePDG::Instance();
+  if (pdgdb == 0x0)
+   {
+     Error("ReadESD","Can not get PDG Database Instance.");
+     return 1;
+   }
+   
+  Float_t mf = esd->GetMagneticField(); 
+
+  if ( (mf == 0.0) && (fNTrackPoints > 0) )
+   {
+      Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
+      return 1;
+   }
+
+  AliStack* stack = 0x0;
+  if (fReadSim && fRunLoader)
+   {
+     fRunLoader->GetEvent(fCurrentEvent);
+     stack = fRunLoader->Stack();
+   }
+
+  const AliESDVertex* vertex = esd->GetVertex();
+  if (vertex == 0x0)
+   {
+     Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
+     vertexpos[0] = 0.0;
+     vertexpos[1] = 0.0;
+     vertexpos[2] = 0.0;
+   }
+  else
+   {
+     vertex->GetXYZ(vertexpos);
+   }
+   
+  if (AliAODParticle::GetDebug() > 0)
+   {
+     Info("ReadESD","Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]);
+   }
+   
+  Info("ReadESD","Reading Event %d",fCurrentEvent);
+
+  Int_t ntr = esd->GetNumberOfTracks();
+  Info("ReadESD","Found %d tracks.",ntr);
+  for (Int_t i = 0;i<ntr; i++)
+   {
+     AliESDtrack *esdtrack = esd->GetTrack(i);
+     if (esdtrack == 0x0)
+      {
+        Error("Next","Can not get track %d", i);
+        continue;
+      }
+
+     //if (esdtrack->HasVertexParameters() == kFALSE) 
+     if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
+      {
+        if (AliAODParticle::GetDebug() > 2) 
+          Info("ReadNext","Particle skipped: Data at vertex not available.");
+        continue;
+      }
+
+     if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE) 
+      {
+        if (AliAODParticle::GetDebug() > 2) 
+          Info("ReadNext","Particle skipped: PID BIT is not set.");
+        continue;
+      }
+
+
+     Double_t extx;
+     Double_t extp[5];
+     esdtrack->GetConstrainedExternalParameters(extx,extp);
+     if (extp[4] == 0.0)
+      {
+        if (AliAODParticle::GetDebug() > 2) 
+          Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
+        continue;
+      } 
+     esdtrack->GetESDpid(pidtable);
+     esdtrack->GetConstrainedPxPyPz(mom);
+     esdtrack->GetConstrainedXYZ(pos);
+     pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point
+     pos[1] -= vertexpos[1];
+     pos[2] -= vertexpos[2];
+
+     Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive
+
+     //Particle from kinematics
+     AliAODStdParticle* particle = 0;
+     Bool_t keeppart = kFALSE;
+     if ( fReadSim && stack  )
+      {
+        if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track 
+        TParticle *p = stack->Particle(esdtrack->GetLabel());
+        if (p==0x0) 
+         {
+           Error("ReadNext","Can not find track with such label.");
+           continue;
+         }
+        if(Pass(p->GetPdgCode())) 
+         {
+           if ( AliAODParticle::GetDebug() > 5 )
+             Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode());
+           continue; //check if we are intersted with particles of this type 
+         }
+//           if(p->GetPdgCode()<0) charge = -1;
+        particle = new AliAODStdParticle(*p,i);
+
+      }
+      
+     if(CheckTrack(esdtrack)) continue;
+      
+     //Here we apply Bayes' formula
+     Double_t rc=0.;
+     for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
+     if (rc==0.0) 
+      {
+        if (AliAODParticle::GetDebug() > 2) 
+          Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
+        continue;
+      }
+
+     for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
+
+     if (AliAODParticle::GetDebug() > 4)
+      { 
+        Info("ReadNext","###########################################################################");
+        Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
+        Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
+        TString msg("Pid list got from track:");
+        for (Int_t s = 0;s<kNSpecies;s++)
+         {
+           msg+="\n    ";
+           msg+=s;
+           msg+="(";
+           msg+=charge*GetSpeciesPdgCode((ESpecies)s);
+           msg+="): ";
+           msg+=w[s];
+           msg+=" (";
+           msg+=pidtable[s];
+           msg+=")";
+         }
+        Info("ReadNext","%s",msg.Data());
+      }//if (AliAODParticle::GetDebug()>4)
+
+      AliTrackPoints* tpts = 0x0;
+      if (fNTrackPoints > 0) 
+       {
+         tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf,fdR);
+       }
+
+      AliClusterMap* cmap = 0x0; 
+      if (  fClusterMap ) 
+       {
+         cmap = new AliClusterMap(esdtrack);
+       }
+     
+     //If this flag fReadMostProbableOnly is false the 
+     //loop over species (see "LOOPPIDS") is over all possible PIDs
+     //in other case the most probablle PID is searched
+     //and the loop is limited to that PID only
+
+     Int_t firstspecie = 0;
+     Int_t lastspecie = kNSpecies;
+     
+     if (fReadMostProbableOnly)
+      { 
+        //find the most probable PID
+        Int_t spec = 0;
+        Float_t maxprob = w[0];
+        for (Int_t s=1; s<AliESDtrack::kSPECIES; s++) 
+         {
+           if (w[s]>maxprob)
+            {
+              maxprob = w[s];
+              spec = s;
+            }
+         }
+        firstspecie = spec;
+        lastspecie = spec + 1;
+      }
+     
+     for (Int_t s = firstspecie; s<lastspecie; s++)//LOOPPIDS
+      {
+        Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
+        Float_t pp = w[s];
+        if (pp == 0.0) 
+         {
+           if ( AliAODParticle::GetDebug() > 5 )
+             Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
+           continue;
+         }
+
+        if(Pass(pdgcode)) 
+         {
+           if ( AliAODParticle::GetDebug() > 5 )
+             Info("ReadNext","PID (%d) did not pass the cut.",pdgcode);
+           continue; //check if we are intersted with particles of this type 
+         }
+
+        Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
+        Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
+
+        AliAODStdParticle* track = new AliAODStdParticle(pdgcode, w[s],i, 
+                                                   mom[0], mom[1], mom[2], tEtot,
+                                                   pos[0], pos[1], pos[2], 0.);
+        //copy probabilitis of other species (if not zero)
+        for (Int_t k = 0; k<kNSpecies; k++)
+         {
+           if (k == s) continue;
+           if (w[k] == 0.0) continue;
+           track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
+         }
+
+        if(Pass(track))//check if meets all criteria of any of our cuts
+                       //if it does not delete it and take next good track
+         { 
+           if ( AliAODParticle::GetDebug() > 4 )
+             Info("ReadNext","Track did not pass the cut");
+           delete track;
+           continue;
+         }
+
+        //Single Particle cuts on cluster map and track points rather do not have sense
+        if (tpts)
+         {
+           track->SetTrackPoints(tpts); 
+         }
+
+        if (cmap) 
+         { 
+           track->SetClusterMap(cmap);
+         }
+
+        fEventRec->AddParticle(track);
+        if (particle) fEventSim->AddParticle(particle);
+        keeppart = kTRUE;
+
+        if (AliAODParticle::GetDebug() > 4 )
+         {
+           Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
+           track->Print();
+           if (particle) particle->Print();
+           Info("ReadNext","\n----------------------------------------------\n");
+         }
+      }//for (Int_t s = 0; s<kNSpecies; s++)
+
+     if (keeppart == kFALSE) 
+      {
+        delete particle;//particle was not stored in event
+        delete tpts;
+        delete cmap;
+      }
+
+   }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
+
+  Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
+         fEventRec->GetNumberOfParticles(), fEventSim->GetNumberOfParticles(),
+         fNEventsRead,fCurrentEvent,fCurrentDir);
+  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
+  return 0;
+}
+
+/**********************************************************/
+
+void AliReaderESD::Rewind()
+{
+  //rewinds reading 
+  delete fKeyIterator;
+  delete fFile;
+  fFile = 0x0;
+  fKeyIterator = 0x0;
+  delete fRunLoader;
+  fRunLoader = 0x0;
+  fCurrentDir = 0;
+  fNEventsRead = 0;
+  if (fTrackCounter) fTrackCounter->Reset();
+}
+/**********************************************************/
+
+TFile* AliReaderESD::OpenFile(Int_t n)
+{
+//opens fFile with kine tree
+
+ const TString& dirname = GetDirName(n);
+ if (dirname == "")
+  {
+   Error("OpenFiles","Can not get directory name");
+   return 0x0;
+  }
+ TString filename = dirname +"/"+ fESDFileName;
+ TFile *ret = TFile::Open(filename.Data()); 
+
+ if ( ret == 0x0)
+  {
+    Error("OpenFiles","Can't open fFile %s",filename.Data());
+    return 0x0;
+  }
+ if (!ret->IsOpen())
+  {
+    Error("OpenFiles","Can't open fFile  %s",filename.Data());
+    return 0x0;
+  }
+ if (fReadSim )
+  {
+   fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
+   if (fRunLoader == 0x0)
+    {
+      Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
+      delete ret;
+      return 0x0;
+    }
+    
+   fRunLoader->LoadHeader();
+   if (fRunLoader->LoadKinematics())
+    {
+      Error("Next","Error occured while loading kinematics.");
+      delete fRunLoader;
+      delete ret;
+      return 0x0;
+    }
+  }
+   
+ return ret;
+}
+/**********************************************************/
+
+Int_t AliReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
+{
+  //returns pdg code from the PID index
+  //ask jura about charge
+  switch (spec)
+   {
+     case kESDElectron:
+       return kPositron;
+       break;
+     case kESDMuon:
+       return kMuonPlus;
+       break;
+     case kESDPion:
+       return kPiPlus;
+       break;
+     case kESDKaon:
+       return kKPlus;
+       break;
+     case kESDProton:
+       return kProton;
+       break;
+     default:
+       ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
+       break;
+   }
+  return 0;
+}
+/********************************************************************/
+Bool_t AliReaderESD::CheckTrack(AliESDtrack* t) const
+{
+  //Performs check of the track
+  
+  if ( (t->GetConstrainedChi2() < fChi2Min) || (t->GetConstrainedChi2() > fChi2Max) ) return kTRUE;
+  
+  if ( (t->GetTPCclusters(0x0) < fNTPCClustMin) || (t->GetTPCclusters(0x0) > fNTPCClustMax) ) return kTRUE;
+
+  if (t->GetTPCclusters(0x0) > 0)
+   {
+     Float_t chisqpercl = t->GetTPCchi2()/((Double_t)t->GetTPCclusters(0x0));
+     if ( (chisqpercl < fTPCChi2PerClustMin) || (chisqpercl > fTPCChi2PerClustMax) ) return kTRUE;
+   }
+
+  Double_t cc[15];
+  t->GetConstrainedExternalCovariance(cc);
+
+  if ( (cc[0]  < fC00Min) || (cc[0]  > fC00Max) ) return kTRUE;
+  if ( (cc[2]  < fC11Min) || (cc[2]  > fC11Max) ) return kTRUE;
+  if ( (cc[5]  < fC22Min) || (cc[5]  > fC22Max) ) return kTRUE;
+  if ( (cc[9]  < fC33Min) || (cc[9]  > fC33Max) ) return kTRUE;
+  if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE;
+
+
+  t->GetInnerExternalCovariance(cc);
+
+  if ( (cc[0]  < fTPCC00Min) || (cc[0]  > fTPCC00Max) ) return kTRUE;
+  if ( (cc[2]  < fTPCC11Min) || (cc[2]  > fTPCC11Max) ) return kTRUE;
+  if ( (cc[5]  < fTPCC22Min) || (cc[5]  > fTPCC22Max) ) return kTRUE;
+  if ( (cc[9]  < fTPCC33Min) || (cc[9]  > fTPCC33Max) ) return kTRUE;
+  if ( (cc[14] < fTPCC44Min) || (cc[14] > fTPCC44Max) ) return kTRUE;
+
+  return kFALSE;
+
+}
+/********************************************************************/
+
+void AliReaderESD::SetChi2Range(Float_t min, Float_t max)
+{
+  //sets range of Chi2 per Cluster
+  fChi2Min = min;
+  fChi2Max = max;
+}
+/********************************************************************/
+
+void AliReaderESD::SetTPCNClustersRange(Int_t min,Int_t max)
+{
+ //sets range of Number Of Clusters that tracks have to have
+ fNTPCClustMin = min;
+ fNTPCClustMax = max;
+}
+/********************************************************************/
+
+void AliReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max)
+{
+  //sets range of Chi2 per Cluster
+  fTPCChi2PerClustMin = min;
+  fTPCChi2PerClustMax = max;
+}
+/********************************************************************/
+
+void AliReaderESD::SetC00Range(Float_t min, Float_t max)
+{
+ //Sets range of C00 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC00Min = min;
+  fC00Max = max;
+}
+/********************************************************************/
+
+void AliReaderESD::SetC11Range(Float_t min, Float_t max)
+{
+ //Sets range of C11 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC11Min = min;
+  fC11Max = max;
+}
+/********************************************************************/
+
+void AliReaderESD::SetC22Range(Float_t min, Float_t max)
+{
+ //Sets range of C22 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC22Min = min;
+  fC22Max = max;
+}
+/********************************************************************/
+
+void AliReaderESD::SetC33Range(Float_t min, Float_t max)
+{
+ //Sets range of C33 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC33Min = min;
+  fC33Max = max;
+}
+/********************************************************************/
+
+void AliReaderESD::SetC44Range(Float_t min, Float_t max)
+{
+ //Sets range of C44 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC44Min = min;
+  fC44Max = max;
+}
diff --git a/ANALYSIS/AliReaderESD.h b/ANALYSIS/AliReaderESD.h
new file mode 100644 (file)
index 0000000..0b6fe2d
--- /dev/null
@@ -0,0 +1,130 @@
+#ifndef AliReaderESD_H
+#define AliReaderESD_H
+//___________________________________________________________________________
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// Multi file reader for ESD                                               //
+//                                                                         //
+// This reader reads tracks from Event Summary Data                        //
+// do not read particles                                                   //
+// Piotr.Skowronski@cern.ch                                                //
+// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html    //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+#include "AliReader.h"
+#include <TString.h>
+class TFile;
+class AliRunLoader;
+class AliESD;
+class AliESDtrack;
+
+class AliReaderESD: public AliReader
+{
+  public:
+    AliReaderESD(const Char_t* esdfilename = "AliESDs.root", const Char_t* galfilename = "galice.root");
+
+    AliReaderESD(TObjArray* dirs,const Char_t* esdfilename = "AliESDs.root", const Char_t* galfilename = "galice.root");
+
+    virtual ~AliReaderESD();
+    
+    void          Rewind();
+    
+    void          ReadSimulatedData(Bool_t flag){fReadSim = flag;}//switches reading MC data
+    Bool_t        ReadsRec() const {return kTRUE;}
+    Bool_t        ReadsSim() const {return fReadSim;}
+    void          SetCheckParticlePID(Bool_t flag){fCheckParticlePID = flag;}
+    void          SetReadMostProbableOnly(Bool_t flag){fReadMostProbableOnly = flag;}
+            
+    void          ReadDataTPC(){}
+    void          ReadDataITS(){}
+
+    void          SetTPCNClustersRange(Int_t min,Int_t max);
+    void          SetTPCChi2PerCluserRange(Float_t min, Float_t max);
+    
+    void          SetChi2Range(Float_t min, Float_t max);
+    void          SetC00Range(Float_t min, Float_t max);
+    void          SetC11Range(Float_t min, Float_t max);
+    void          SetC22Range(Float_t min, Float_t max);
+    void          SetC33Range(Float_t min, Float_t max);
+    void          SetC44Range(Float_t min, Float_t max);
+    void          SetNumberOfTrackPoints(Int_t n = 5,Float_t dr = 30.0) {fNTrackPoints = n; fdR = dr;}
+    Int_t         GetNumberOfTrackPoints() const {return fNTrackPoints;}
+    void          SetClusterMap(Bool_t flag = kTRUE){fClusterMap = flag;}
+
+    
+    enum ESpecies {kESDElectron = 0, kESDMuon, kESDPion, kESDKaon, kESDProton, kNSpecies};
+    static Int_t  GetSpeciesPdgCode(ESpecies spec);//skowron
+    
+    Int_t         ReadESD(AliESD* esd);
+    
+  protected:
+    Int_t         ReadNext();
+    TFile*        OpenFile(Int_t evno);//opens files to be read for given event
+    Bool_t        CheckTrack(AliESDtrack* t) const;
+    
+    TString       fESDFileName;//name of the file with tracks
+    TString       fGAlFileName;//name of the file with tracks
+    TFile*        fFile;//! pointer to current ESD file
+    AliRunLoader* fRunLoader;//!Run Loader
+    TIter*        fKeyIterator;//!iterator over keys in ESD file
+    Bool_t        fReadSim;//flag indicating wether to read particles from kinematics
+    Bool_t        fCheckParticlePID;//flag indicating to perform the check on PID of simulated particle - usefull in resoluion analysis
+    Bool_t        fReadMostProbableOnly;//flag indicating to read ony one incarnation with the highest probability
+    Int_t         fNTrackPoints;//number of track points; if==0 track points are not created
+    Float_t       fdR;//spacing between points (along radius) in cm
+                      //Track Points are needed for Anti-Merging Cut
+    
+    Bool_t        fClusterMap;//Flag indicating if Claster Map should be created for each track
+                              //Claster map is needed for Anti-Splitting Cut
+
+    //Cut Parameters specific to TPC tracks
+        
+    Int_t         fNTPCClustMin;//Number of clusters min value
+    Int_t         fNTPCClustMax;//Number of clusters max value
+    
+    Float_t       fTPCChi2PerClustMin;//Chi^2 per number of clusters min value
+    Float_t       fTPCChi2PerClustMax;//Chi^2 per number of clusters max value
+
+
+    // Required parameters at vertex
+    Float_t       fChi2Min;//Chi^2 min value
+    Float_t       fChi2Max;//Chi^2 max value
+
+    Float_t       fC00Min;//C00 (0th diagonal element of covariance matrix) min value
+    Float_t       fC00Max;//C00 (0th diagonal element of covariance matrix) max value
+            
+    Float_t       fC11Min;//C11 (1th diagonal element of covariance matrix) min value
+    Float_t       fC11Max;//C11 (1th diagonal element of covariance matrix) max value
+    
+    Float_t       fC22Min;//C22 (2th diagonal element of covariance matrix) min value
+    Float_t       fC22Max;//C22 (2th diagonal element of covariance matrix) max value
+    
+    Float_t       fC33Min;//C33 (3th diagonal element of covariance matrix) min value
+    Float_t       fC33Max;//C33 (3th diagonal element of covariance matrix) max value
+    
+    Float_t       fC44Min;//C44 (4th diagonal element of covariance matrix) min value
+    Float_t       fC44Max;//C44 (4th diagonal element of covariance matrix) max value
+
+    // Required parameters at TPC Inner Layer
+    Float_t       fTPCC00Min;//C00 (0th diagonal element of covariance matrix) min value
+    Float_t       fTPCC00Max;//C00 (0th diagonal element of covariance matrix) max value
+            
+    Float_t       fTPCC11Min;//C11 (1th diagonal element of covariance matrix) min value
+    Float_t       fTPCC11Max;//C11 (1th diagonal element of covariance matrix) max value
+    
+    Float_t       fTPCC22Min;//C22 (2th diagonal element of covariance matrix) min value
+    Float_t       fTPCC22Max;//C22 (2th diagonal element of covariance matrix) max value
+    
+    Float_t       fTPCC33Min;//C33 (3th diagonal element of covariance matrix) min value
+    Float_t       fTPCC33Max;//C33 (3th diagonal element of covariance matrix) max value
+    
+    Float_t       fTPCC44Min;//C44 (4th diagonal element of covariance matrix) min value
+    Float_t       fTPCC44Max;//C44 (4th diagonal element of covariance matrix) max value
+    
+  private:
+    ClassDef(AliReaderESD,1)
+};
+
+
+#endif
index 321b4e5..20d9043 100644 (file)
 #include <AliESD.h>
 
 #include "AliEventCut.h"
+#include "AliReader.h"
+#include "AliAODParticle.h"
 
 
 ClassImp(AliRunAnalysis)
 AliRunAnalysis::AliRunAnalysis():
  TTask("RunAnalysis","Alice Analysis Manager") ,
- fDirs(),
+ fReader(0x0),
  fEventCut(0x0),
- fFileName("AliESDs.root"),
- fReadKinematics(kFALSE)
+ fCutOnSim(kFALSE),
+ fCutOnRec(kTRUE)
 {
   //ctor
 }
@@ -43,8 +45,8 @@ AliRunAnalysis::AliRunAnalysis():
 AliRunAnalysis::~AliRunAnalysis()
 {
   //dtor
-  delete fDirs;
-  delete fAnalysies;
+  delete fEventCut;
+  delete fReader;
   delete fEventCut;
 }
 /*********************************************************/
@@ -52,249 +54,75 @@ AliRunAnalysis::~AliRunAnalysis()
 Int_t AliRunAnalysis::Run()
 {
  //makes analysis
- Int_t currentdir = 0;
- Int_t ndirs;
- if (fDirs) //if array with directories is supplied by user
-  {
-    ndirs = fDirs->GetEntries(); //get the number if directories
-  }
- else
-  {
-    ndirs = 0; //if the array is not supplied read only from current directory
-  }
 
+  if (fReader == 0x0)
+   {
+     Error("Run","Reader is not set");
+     return 1;
+   }
  /******************************/ 
  /*  Init Event                */ 
  /******************************/ 
if (fAnalysies == 0x0)
for (Int_t an = 0; an < fAnalysies.GetEntries(); an++)
   {
-    Info("Run","No analysis present");
-    return 0;
-  }
-  
- for (Int_t an = 0; an < fAnalysies->GetEntries(); an++)
-  {
-      AliAnalysis* analysis = (AliAnalysis*)fAnalysies->At(an);
+      AliAnalysis* analysis = (AliAnalysis*)fAnalysies.At(an);
       analysis->Init();
   }
-
- do
+ while (fReader->Next() == kFALSE)
   {
-    TFile* file = OpenFile(currentdir);
-    if (file == 0x0)
-      {
-        Error("Run","Cannot get File for dir no. %d",currentdir);
-        currentdir++;
-        continue;
-      } 
-    AliStack* stack = 0x0;
-    AliRunLoader* rl = 0x0;
-    if (fReadKinematics)
-     {
-       const TString& dirname = GetDirName(currentdir);
-       if (dirname == "")
-        {
-         Error("Run","Can not get directory name");
-         return 0x0;
-        }
-       TString filename = dirname +"/galice.root";
-
-       rl = AliRunLoader::Open(filename);
-       if ( rl == 0x0 )
-        {
-          Error("Run","Can't get Run Loader from dir %s",filename.Data());
-          delete file;
-          currentdir++;
-          continue;
-        }
-       if( rl->LoadHeader() )
-        {
-          Error("Run","Error while loading Header from dir %s",filename.Data());
-          delete file;
-          delete rl;
-          currentdir++;
-          continue;
-        }
-       if( rl->LoadKinematics() )
-        {
-          Error("Run","Error while loading Kinematics from dir %s",filename.Data());
-          delete file;
-          delete rl;
-          currentdir++;
-          continue;
-        }
-     }
+   AliAOD* eventsim = fReader->GetEventSim();
+   AliAOD* eventrec = fReader->GetEventRec();
    
-    file->cd();
-    TIter keyiter(file->GetListOfKeys());
-    TKey* key;
-    while (( key = (TKey*)keyiter.Next() ))
-     {
-      if (key == 0x0)
-       {
-        if (GetDebug() > 2 )
-          {
-            Info("Run","No more keys.");
-          }
-        break;
-       }
-
-      TObject* esdobj = key->ReadObj();
-      if (esdobj == 0x0)
-       {
-         if (GetDebug() > 2 )
-           {
-             Info("ReadNext","Key read NULL. Key Name is %s",key->GetName());
-             key->Dump();
-           }
-         continue;
-       }
-      if (GetDebug() > 9 ) esdobj->Dump();
-      AliESD* esd = dynamic_cast<AliESD*>(esdobj);
-
-      if (esd == 0x0)
-       {
-         if (GetDebug() > 7 )
-           {
-             Info("ReadNext","It is not an ESD object");
-           }
-         delete esdobj;
-         continue;
-       }
-
-      if (fReadKinematics)
-       {
-        TString esdname(esd->GetName());
-        esdname.ReplaceAll("ESD","");
-        Int_t nev = atoi(esdname);
-        Info("Run","ESD name is %s, Event Number is %d",esd->GetName(),nev);
-        if (rl->GetEvent(nev))
-         {
-           Error("Run","Error occured while RunLoader GetEvent %d",nev);
-           delete esd;
-           continue;
-         }
-        stack = rl->Stack();
-        if (stack == 0x0) 
-         {
-           Error("Run","Can not get stack");
-           delete esd;
-           continue;
-         }
-       }
       /******************************/ 
       /*  Event Cut                 */ 
       /******************************/ 
-      if (fEventCut)
+      if ( Pass(eventrec,eventsim) )
        {
-         if (fEventCut->Pass(esd))
-          {
-            if (GetDebug()) Info("Run","Event rejected by Event Cut");
-            delete esd;
-            continue; //Did not pass the 
-          }
+         if (AliAODParticle::GetDebug()) Info("Run","Event rejected by Event Cut");
+         continue; //Did not pass the 
        }
       /******************************/ 
       /*  Process Event             */ 
       /******************************/ 
-      for (Int_t an = 0; an < fAnalysies->GetEntries(); an++)
+      for (Int_t an = 0; an < fAnalysies.GetEntries(); an++)
        {
-           AliAnalysis* analysis = (AliAnalysis*)fAnalysies->At(an);
-           analysis->ProcessEvent(esd,stack);
+           AliAnalysis* analysis = (AliAnalysis*)fAnalysies.At(an);
+           analysis->ProcessEvent(eventsim,eventrec);
        }
-      delete esd;
-     }//end of loop over keys in file
-     
-   delete file;
-   delete rl;
-   currentdir++;
     
-  }while (currentdir < ndirs);//end of loop over directories
+  }//end of loop over events
 
  /******************************/ 
  /*  Finish Event              */ 
  /******************************/ 
- for (Int_t an = 0; an < fAnalysies->GetEntries(); an++)
+ for (Int_t an = 0; an < fAnalysies.GetEntries(); an++)
   {
-      AliAnalysis* analysis = (AliAnalysis*)fAnalysies->At(an);
-      analysis->Init();
+      AliAnalysis* analysis = (AliAnalysis*)fAnalysies.At(an);
+      analysis->Finish();
   }
 
  return 0;   
 }
 /*********************************************************/
 
-TFile* AliRunAnalysis::OpenFile(Int_t n)
-{
-//opens file with kine tree
-
- const TString& dirname = GetDirName(n);
- if (dirname == "")
-  {
-   Error("OpenFiles","Can not get directory name");
-   return 0x0;
-  }
- TString filename = dirname +"/"+ fFileName;
- TFile *ret = TFile::Open(filename.Data()); 
-
- if ( ret == 0x0)
-  {
-    Error("OpenFiles","Can't open file %s",filename.Data());
-    return 0x0;
-  }
- if (!ret->IsOpen())
-  {
-    Error("OpenFiles","Can't open file  %s",filename.Data());
-    return 0x0;
-  }
- return ret;
-}
-/*********************************************************/
-
-TString& AliRunAnalysis::GetDirName(Int_t entry)
+void  AliRunAnalysis::Add(AliAnalysis* a)
 {
-//returns directory name of next one to read
-  TString* retval;//return value
-  if (fDirs ==  0x0)
-   {
-     retval = new TString(".");
-     return *retval;
-   }
-
-  if ( (entry>fDirs->GetEntries()) || (entry<0))//if out of bounds return empty string
-   {                                            //note that entry==0 is accepted even if array is empty (size=0)
-     Error("GetDirName","Name out of bounds");
-     retval = new TString();
-     return *retval;
-   }
-
-  if (fDirs->GetEntries() == 0)
-   { 
-     retval = new TString(".");
-     return *retval;
-   }
-
-  TClass *objclass = fDirs->At(entry)->IsA();
-  TClass *stringclass = TObjString::Class();
-
-  TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry));
-
-  if(dir == 0x0)
-   {
-     Error("GetDirName","Object in TObjArray is not a TObjString or its descendant");
-     retval = new TString();
-     return *retval;
-   }
-  if (gDebug > 0) Info("GetDirName","Returned ok %s",dir->String().Data());
-  return dir->String();
+  //adds a to the list of analysis
+  fAnalysies.Add(a);
 }
 /*********************************************************/
 
-void  AliRunAnalysis::Add(AliAnalysis* a)
+Bool_t AliRunAnalysis::Pass(AliAOD* recevent, AliAOD* simevent)
 {
-  //adds a to the list of analysis
-  if (fAnalysies == 0x0) fAnalysies = new TObjArray();
-  fAnalysies->Add(a);
+  //checks the event cut
+  if (fEventCut == 0x0) return kFALSE;
+  
+  if (fCutOnRec)
+    if (fEventCut->Pass(recevent)) return kTRUE;
+    
+  if (fCutOnSim)
+    if (fEventCut->Pass(simevent)) return kTRUE;
+  
+  return kFALSE;
 }
index 127084c..4b369a4 100644 (file)
 //
 ///////////////////////////////////////////////////////////
 
-class AliEventCut;
-class TObjArray;
-class TFile;
-
-#include <TString.h>
 #include <TTask.h>
-
+#include <TObjArray.h>
 #include "AliAnalysis.h"
 
+class AliEventCut;
+class TFile;
+class AliReader;
+
 class AliRunAnalysis: public TTask
 {
   public: 
     AliRunAnalysis();
     virtual ~AliRunAnalysis();
     
-    Int_t Run();
-    void  Add(AliAnalysis* a);
-    void  ReadKinematics(Bool_t flag){fReadKinematics = flag;}
+    Int_t         Run();
+    void          Add(AliAnalysis* a);
+    void          SetReader(AliReader* reader){fReader = reader;}
+    
+    const char*   GetName(){return "RunAnalysis";}
     
-    Int_t GetDebug() {return AliAnalysis::GetDebug();}
-    void SetDirs(TObjArray* dirs){fDirs = dirs;} //sets array directories names;
   protected:
-    TObjArray*    fAnalysies;//arry with analysies
-    TObjArray*    fDirs;//arry with directories to read data from
+    TObjArray     fAnalysies;//arry with analysies
+    AliReader*    fReader;//arry with directories to read data from
     
     AliEventCut*  fEventCut;//event cut    
     
-    TString       fFileName;//name of the file with ESDs
-    Bool_t        fReadKinematics;
+    Bool_t        fCutOnSim;//flag indicating that event cut is performed on simulated particles 
+    Bool_t        fCutOnRec;//flag indicating that event cut is performed on reconstructed tracks
     
-    TString& GetDirName(Int_t entry);
-    TFile* OpenFile(Int_t n);
+    Bool_t        Pass(AliAOD* recevent, AliAOD* simevent);
     
   private:
     void SetName(const char *){}//change SetName to be private
diff --git a/ANALYSIS/AliTrackPoints.cxx b/ANALYSIS/AliTrackPoints.cxx
new file mode 100644 (file)
index 0000000..3878afd
--- /dev/null
@@ -0,0 +1,589 @@
+//_________________________________
+////////////////////////////////////////////////////////////
+//                                                        //
+// class AliTrackPoints                                //
+//                                                        //
+// used by Anti-Merging cut                               //
+// contains set of poits the lay on track trajectory      //
+// according to reconstructed track parameters -          //
+// NOT CLUSTERS POSITIONS!!!                              //
+// Anti-Merging cut is applied only on tracks coming from //
+// different events (that are use to fill deniminators)   //
+//                                                        //
+////////////////////////////////////////////////////////////
+
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TMath.h>
+
+#include "AliESDtrack.h"
+#include "AliTrackPoints.h"
+#include "AliTPCtrack.h"
+#include "AliTrackReference.h"
+
+ClassImp(AliTrackPoints)
+
+Int_t AliTrackPoints::fgDebug = 0;
+AliTrackPoints::AliTrackPoints():
+ fN(0),
+ fX(0x0),
+ fY(0x0),
+ fZ(0x0)
+{
+  //constructor
+}
+/***************************************************************/
+
+AliTrackPoints::AliTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr, Float_t r0):
+ fN(n),
+ fX(new Float_t[fN]),
+ fY(new Float_t[fN]),
+ fZ(new Float_t[fN])
+{
+  //constructor
+  //mf - magnetic field in kG - needed to calculated curvature out of Pt
+  //r0 - starting radius
+  //dr - calculate points every dr cm, default every 30cm
+  if (track == 0x0)
+   {
+     Error("AliTrackPoints","ESD track is null");
+     fN = 0;
+     delete [] fX;
+     delete [] fY;
+     delete [] fZ;
+     fX = fY = fZ = 0x0;
+     return;
+   }
+
+  if ( ((track->GetStatus() & AliESDtrack::kTPCrefit) == kFALSE)&& 
+       ((track->GetStatus() & AliESDtrack::kTPCin) == kFALSE)  )
+   {
+     //could happend: its stand alone tracking
+     if (GetDebug() > 3) 
+       Warning("AliTrackPoints","This ESD track does not contain TPC information");
+       
+     fN = 0;
+     delete [] fX;
+     delete [] fY;
+     delete [] fZ;
+     fX = fY = fZ = 0x0;
+     
+     return;
+   }
+  
+  Double_t x;
+  Double_t par[5];
+  track->GetInnerExternalParameters(x,par);     //get properties of the track
+  if (par[4] == 0)
+   {
+     Error("AliTrackPoints","This ESD track seem not to contain TPC information (curv is 0)");
+     return;
+   }
+
+  if (mf == 0.0)
+   {
+     Error("AliTrackPoints","Zero Magnetic field passed as parameter.");
+     return;
+   }
+   
+  Double_t alpha = track->GetInnerAlpha();
+  Double_t cc = 1000./0.299792458/mf;//conversion constant
+  Double_t c=par[4]/cc;
+  
+  MakePoints(dr,r0,x,par,c,alpha);
+  
+}
+/***************************************************************/
+
+AliTrackPoints::AliTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
+ fN(n),
+ fX(new Float_t[fN]),
+ fY(new Float_t[fN]),
+ fZ(new Float_t[fN])
+{
+  //constructor
+  //r0 starting radius
+  //dr - calculate points every dr cm, default every 30cm
+  if (track == 0x0)
+   {
+     Error("AliTrackPoints","TPC track is null");
+     fN = 0;
+     delete [] fX;
+     delete [] fY;
+     delete [] fZ;
+     fX = fY = fZ = 0x0;
+     return;
+   }
+  track->PropagateTo(r0);
+
+ //* This formation is now fixed in the following way:                         *
+ //*      external param0:   local Y-coordinate of a track (cm)                *
+ //*      external param1:   local Z-coordinate of a track (cm)                *
+ //*      external param2:   local sine of the track momentum azimuth angle    *
+ //*      external param3:   tangent of the track momentum dip angle           *
+ //*      external param4:   1/pt (1/(GeV/c))                                  *
+    
+  Double_t x = 0;
+  Double_t par[5];
+  track->GetExternalParameters(x,par);     //get properties of the track
+  
+  Double_t alpha = track->GetAlpha();
+  Double_t c=track->GetC();
+  MakePoints(dr,r0,x,par,c,alpha);
+}  
+
+void AliTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha)
+{
+  //Calculates points starting at radius r0
+  //spacing every dr (in radial direction)
+  // according to track parameters
+  // x - position in sector local reference frame. x i parallel to R and sector is symmetric with respect to x
+  // par - track external parameters; array with 5 elements;  look at AliTPCtrack.h or AliESDtrack.h for their meaning
+  // c - track curvature
+  // alpha - sector's rotation angle (phi) == angle needed for local to global transformation
+  
+  Double_t y = par[0];
+  Double_t z0 = par[1];
+      
+  Double_t phi0local = TMath::ATan2(y,x);
+  Double_t phi0global = phi0local + alpha;
+  
+  if (phi0local<0) phi0local+=2*TMath::Pi();
+  if (phi0local>=2.*TMath::Pi()) phi0local-=2*TMath::Pi();
+  
+  if (phi0global<0) phi0global+=2*TMath::Pi();
+  if (phi0global>=2.*TMath::Pi()) phi0global-=2*TMath::Pi();
+  
+  Double_t r = TMath::Hypot(x,y);
+  
+
+  if (GetDebug() > 9) 
+    Info("AliTrackPoints","Radius0 %f, Real Radius %f",r0,r); 
+  
+  if (GetDebug() > 5) 
+    Info("AliTrackPoints","Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local);
+    
+  Double_t eta = x*c - par[2] ;//par[2] = fX*C - eta; eta==fP2 ; C==fP4
+  
+  //this calculattions are assuming new (current) model 
+  Double_t tmp = par[2];
+  tmp = 1. - tmp*tmp;
+  tmp = c*y + TMath::Sqrt(tmp);
+  Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);
+
+  //Here is old model Cold=Cnew/2.
+  Double_t dcasq = dca*dca;
+  Double_t c2 = c/2.;
+  Double_t cst1 = (1.+c2*dca)*dca;//first constant
+  Double_t cst2 = 1. + 2.*c2*dca;//second constant
+      
+  Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
+  Double_t factorZ0 = TMath::ASin(c2*TMath::Sqrt((r*r-dcasq)/cst2))*par[3]/c2;
+  for(Int_t i = 0; i<fN; i++)
+   {
+     Double_t rc = r0 + i*dr;
+     Double_t ftmp = (c2*rc + cst1/rc)/cst2;
+     if (ftmp > 1.0) 
+      {
+        if (GetDebug() > 1) 
+          Warning("AliTrackPoints","ASin argument > 1 %f:",ftmp);
+        ftmp=1.0;
+      }
+     else if (ftmp < -1.0) 
+      {
+        if (GetDebug() > 1) 
+          Warning("AliTrackPoints","ASin argument < -1 %f:",ftmp);
+        ftmp=-1.0;
+      }
+      
+     Double_t factorPhi = TMath::ASin( ftmp );//factor phi od rc
+     Double_t phi = phi0global + factorPhi - factorPhi0;
+     
+     ftmp = (rc*rc-dcasq)/cst2;
+     if (ftmp < 0.0) 
+      {
+        if (GetDebug() > 1) 
+          Warning("AliTrackPoints","Sqrt argument < 0: %f",ftmp);
+        ftmp=0.0;
+      }
+     
+     ftmp = c2*TMath::Sqrt(ftmp);
+     if (ftmp > 1.0) 
+      {
+        if (GetDebug() > 1) 
+          Warning("AliTrackPoints","ASin argument > 1: %f",ftmp);
+        ftmp=1.0;
+      }
+     else if (ftmp < -1.0) 
+      {
+        if (GetDebug() > 1) 
+          Warning("AliTrackPoints","ASin argument < -1: %f",ftmp);
+        ftmp=-1.0;
+      }
+     Double_t factorZ = TMath::ASin(ftmp)*par[3]/c2;
+     fZ[i] = z0 + factorZ - factorZ0;
+     fX[i] = rc*TMath::Cos(phi);
+     fY[i] = rc*TMath::Sin(phi);
+     
+     if ( GetDebug() > 2 )
+      {
+        Info("AliTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
+             fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i]));
+      }
+   }
+}
+/***************************************************************/
+
+AliTrackPoints::~AliTrackPoints()
+{
+  //destructor
+  delete [] fX;
+  delete [] fY;
+  delete [] fZ;
+}
+/***************************************************************/
+
+void AliTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
+{
+  //returns position at point n
+  if ((n<0) || (n>fN))
+   {
+     Error("PositionAt","Point %d out of range",n);
+     return;
+   }
+  
+  x = fX[n];
+  y = fY[n];
+  z = fZ[n];
+  if ( GetDebug() > 1 )
+    {
+      Info("AliTrackPoints","n %d; X %f; Y %f; Z %f",n,x,y,z);
+    }
+}
+/***************************************************************/
+
+Double_t AliTrackPoints::AvarageDistance(const AliTrackPoints& tr)
+{
+  //returns the aritmethic avarage distance between two tracks
+//  Info("AvarageDistance","Entered");
+  if ( (fN <= 0) || (tr.fN <=0) )
+   {
+     if (GetDebug()) Warning("AvarageDistance","One of tracks is empty");
+     return -1;
+   }
+
+  if (fN != tr.fN)
+   {
+     Warning("AvarageDistance","Number of points is not equal");
+     return -1;
+   }
+   
+  Double_t sum = 0;
+  for (Int_t i = 0; i<fN; i++)
+   {
+     if (GetDebug()>9)
+      {
+//       Float_t r1sq = fX[i]*fX[i]+fY[i]*fY[i];
+//       Float_t r2sq = tr.fX[i]*tr.fX[i]+tr.fY[i]*tr.fY[i];
+       Float_t r1sq = TMath::Hypot(fX[i],fY[i]);
+       Float_t r2sq = TMath::Hypot(tr.fX[i],tr.fY[i]);
+       Info("AvarageDistance","radii: %f %f",r1sq,r2sq);
+      } 
+
+      
+     Double_t dx = fX[i]-tr.fX[i];
+     Double_t dy = fY[i]-tr.fY[i];
+     Double_t dz = fZ[i]-tr.fZ[i];
+     sum+=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
+     
+     if (GetDebug()>1)
+      {
+       Info("AvarageDistance","Diff: x ,y z: %f , %f, %f",dx,dy,dz);
+       Info("AvarageDistance","xxyyzz %f %f %f %f %f %f",
+            fX[i],tr.fX[i],fY[i],tr.fY[i],fZ[i],tr.fZ[i]);
+      } 
+   }
+   
+  Double_t retval = sum/((Double_t)fN);
+  if ( GetDebug() )
+    {
+      Info("AvarageDistance","Avarage distance is %f.",retval);
+    }
+  return retval;
+}
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+
+#include "AliRun.h"
+#include "AliESD.h"
+#include "AliRunLoader.h"
+#include "AliTPCtrack.h"
+#include "TTree.h"
+#include "TBranch.h"
+#include "TH2D.h"
+#include "TCanvas.h"
+#include "AliMagF.h"
+
+
+
+
+void AliTrackPoints::TestESD(Int_t entr,const char* fname )
+{
+//Tests creation of track points for ESD tracks
+  delete gAlice;
+  gAlice = 0x0;
+  AliRunLoader* rl = AliRunLoader::Open();
+  rl->LoadgAlice();
+  
+  Float_t mf = rl->GetAliRun()->Field()->SolenoidField();
+  
+  TFile* fFile = TFile::Open(fname);
+  
+  if (fFile == 0x0)
+   {
+     printf("testesd: There is no suche a ESD file\n");
+     return;
+   }
+  AliESD* esd = dynamic_cast<AliESD*>(fFile->Get("0"));
+  AliESDtrack *t = esd->GetTrack(entr);
+  if (t == 0x0)
+   {
+     ::Error("testesd","Can not get track %d",entr);
+     return;
+   }
+
+  
+  Int_t N = 170;
+  AliTrackPoints* tp = new AliTrackPoints(N,t,mf,1.);
+  
+  Float_t xmin = -250;
+  Float_t xmax =  250;
+  
+  Float_t ymin = -250;
+  Float_t ymax =  250;
+
+  Float_t zmin = -250;
+  Float_t zmax =  250;
+  
+  TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
+  TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
+  TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
+
+  TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
+  TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
+  TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
+
+  hxyt->SetDirectory(0x0);   
+  hxy->SetDirectory(0x0);   
+  hxyTR->SetDirectory(0x0);   
+  
+  hxzt->SetDirectory(0x0);   
+  hxz->SetDirectory(0x0);   
+  hxzTR->SetDirectory(0x0);   
+
+  Float_t x,y,z;
+   
+  for (Int_t i = 0;i<N;i++)
+   {  
+      Double_t r = 84.1+i;
+      tp->PositionAt(i,x,y,z);
+      hxy->Fill(x,y);
+      hxz->Fill(x,z);
+      printf("Rdemanded %f\n",r);
+      printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
+      
+   }
+
+  rl->LoadTrackRefs();
+  TTree* treeTR = rl->TreeTR();
+  TBranch* b = treeTR->GetBranch("TPC");
+  
+  TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
+  AliTrackReference* tref;
+  b->SetAddress(&trackrefs);
+  
+  Int_t tlab = TMath::Abs(t->GetLabel());
+
+  Int_t netr = (Int_t)treeTR->GetEntries();
+  printf("Found %d entries in TR tree\n",netr);
+  
+  for (Int_t e = 0; e < netr; e++)
+   { 
+     treeTR->GetEntry(e);
+     tref = (AliTrackReference*)trackrefs->At(0);
+     if (tref == 0x0) continue;
+     if (tref->GetTrack() != tlab) continue;
+    
+     printf("Found %d entries in TR array\n",trackrefs->GetEntries());
+     
+     for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
+      {
+        tref = (AliTrackReference*)trackrefs->At(i);
+        if (tref->GetTrack() != tlab) continue;
+        x = tref->X();
+        y = tref->Y();
+        z = tref->Z();
+        printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
+
+        hxzTR->Fill(x,z);
+        hxyTR->Fill(x,y);
+        for (Int_t j = 1; j < 10; j++)
+         {
+            hxyTR->Fill(x, y+j*0.1);
+            hxyTR->Fill(x, y-j*0.1);
+            hxyTR->Fill(x+j*0.1,y);
+            hxyTR->Fill(x-j*0.1,y);
+
+            hxzTR->Fill(x,z-j*0.1);
+            hxzTR->Fill(x,z+j*0.1);
+            hxzTR->Fill(x-j*0.1,z);
+            hxzTR->Fill(x+j*0.1,z);
+         }
+      }
+      break; 
+    }  
+  hxy->Draw("");
+//  hxzt->Draw("same");
+  hxyTR->Draw("same");
+  
+  delete rl;
+}
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+
+void AliTrackPoints::TestTPC(Int_t entr)
+{
+//Tests creation of track points for TPC tracks
+  delete gAlice;
+  gAlice = 0x0;
+  AliRunLoader* rl = AliRunLoader::Open();
+  AliLoader* l = rl->GetLoader("TPCLoader");
+  rl->LoadgAlice();
+  AliKalmanTrack::SetConvConst(100/0.299792458/0.2/rl->GetAliRun()->Field()->Factor());
+  l->LoadTracks();
+  AliTPCtrack* t = new AliTPCtrack();
+  TBranch* b=l->TreeT()->GetBranch("tracks");
+  b->SetAddress(&t);
+  l->TreeT()->GetEntry(entr);  
+  Int_t N = 160;
+  AliTrackPoints* tp = new AliTrackPoints(N,t,1.);
+  
+  Float_t xmin = -250;
+  Float_t xmax =  250;
+  
+  Float_t ymin = -250;
+  Float_t ymax =  250;
+
+  Float_t zmin = -250;
+  Float_t zmax =  250;
+  
+  TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
+  TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
+  TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
+
+  TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
+  TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
+  TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
+
+  hxyt->SetDirectory(0x0);   
+  hxy->SetDirectory(0x0);   
+  hxyTR->SetDirectory(0x0);   
+  
+  hxzt->SetDirectory(0x0);   
+  hxz->SetDirectory(0x0);   
+  hxzTR->SetDirectory(0x0);   
+
+  Float_t x,y,z;
+   
+  for (Int_t i = 0;i<N;i++)
+   {  
+      Double_t r = 84.1+i;
+      tp->PositionAt(i,x,y,z);
+      hxy->Fill(x,y);
+      hxz->Fill(x,z);
+      printf("Rdemanded %f\n",r);
+      printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
+      
+      //BUT they are local!!!!
+      t->PropagateTo(r);
+//      Double_t phi = t->Phi();
+      Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius 
+      
+      Double_t alpha = t->GetAlpha();
+      Double_t salpha = TMath::Sin(alpha);
+      Double_t calpha = TMath::Cos(alpha);
+      x = t->GetX()*calpha - t->GetY()*salpha;
+      y = t->GetX()*salpha + t->GetY()*calpha;
+      z = t->GetZ();
+      
+      printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl);
+      
+      printf("tpz - tz %f\n",z-t->GetZ());
+      printf("\n");
+      hxyt->Fill(x,y);
+      hxzt->Fill(x,z);
+      
+   }
+
+  rl->LoadTrackRefs();
+  TTree* treeTR = rl->TreeTR();
+  b = treeTR->GetBranch("TPC");
+  
+  TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
+  AliTrackReference* tref;
+  b->SetAddress(&trackrefs);
+  
+  Int_t tlab = TMath::Abs(t->GetLabel());
+
+  Int_t netr = (Int_t)treeTR->GetEntries();
+  printf("Found %d entries in TR tree\n",netr);
+  
+  for (Int_t e = 0; e < netr; e++)
+   { 
+     treeTR->GetEntry(e);
+     tref = (AliTrackReference*)trackrefs->At(0);
+     if (tref == 0x0) continue;
+     if (tref->GetTrack() != tlab) continue;
+    
+     printf("Found %d entries in TR array\n",trackrefs->GetEntries());
+     
+     for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
+      {
+        tref = (AliTrackReference*)trackrefs->At(i);
+        if (tref->GetTrack() != tlab) continue;
+        x = tref->X();
+        y = tref->Y();
+        z = tref->Z();
+        printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
+
+        hxzTR->Fill(x,z);
+        hxyTR->Fill(x,y);
+        for (Int_t j = 1; j < 10; j++)
+         {
+            hxyTR->Fill(x, y+j*0.1);
+            hxyTR->Fill(x, y-j*0.1);
+            hxyTR->Fill(x+j*0.1,y);
+            hxyTR->Fill(x-j*0.1,y);
+
+            hxzTR->Fill(x,z-j*0.1);
+            hxzTR->Fill(x,z+j*0.1);
+            hxzTR->Fill(x-j*0.1,z);
+            hxzTR->Fill(x+j*0.1,z);
+         }
+      }
+      break; 
+    }  
+  hxz->Draw("");
+//  hxzt->Draw("same");
+  hxzTR->Draw("same");
+  
+  delete rl;
+}
diff --git a/ANALYSIS/AliTrackPoints.h b/ANALYSIS/AliTrackPoints.h
new file mode 100644 (file)
index 0000000..74bf907
--- /dev/null
@@ -0,0 +1,49 @@
+#ifndef AliTrackPoints_H
+#define AliTrackPoints_H
+//_________________________________
+////////////////////////////////////////////////////////////
+//                                                        //
+// class AliTrackPoints                                //
+//                                                        //
+// used by Anti-Merging cut                               //
+// contains set of poits the lay on track trajectory      //
+// according to reconstructed track parameters -          //
+// NOT CLUSTERS POSITIONS!!!                              //
+// Anti-Merging cut is applied only on tracks coming from //
+// different events (that are use to fill deniminators)   //
+//                                                        //
+////////////////////////////////////////////////////////////
+#include <TObject.h>
+
+class AliTPCtrack;
+class AliESDtrack;
+
+class AliTrackPoints: public TObject
+{
+  public:
+    AliTrackPoints();
+    AliTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr=30, Float_t r0 = 84.1); //min TPC R  = 84.1; max TPC R =  246.6cm, 
+    AliTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr=30,Float_t r0 = 84.1); //min TPC R  = 84.1; max TPC R =  246.6cm, 
+    AliTrackPoints(const AliTrackPoints& in);
+    
+    virtual ~AliTrackPoints();
+    AliTrackPoints& operator=(const AliTrackPoints& in);
+    
+    Double_t AvarageDistance(const AliTrackPoints& tr);
+    void PositionAt(Int_t n, Float_t &x, Float_t &y, Float_t &z);
+    Int_t GetDebug() const {return fgDebug;}
+    void  SetDebug(Int_t deblevel){fgDebug = deblevel;} 
+    static void TestTPC(Int_t entr);
+    static void TestESD(Int_t entr,const char* fname = "AliESDs.root");
+  protected:
+    void MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha);
+  private:
+    Int_t    fN;//number of points
+    Float_t* fX;//[fN]positions at x
+    Float_t* fY;//[fN]positions at y
+    Float_t* fZ;//[fN] positions at z
+//    Float_t* fR;//! [fN] radii
+    static Int_t fgDebug;//! debug level
+    ClassDef(AliTrackPoints,1)
+};
+#endif
diff --git a/ANALYSIS/analysis.C b/ANALYSIS/analysis.C
new file mode 100644 (file)
index 0000000..090b82b
--- /dev/null
@@ -0,0 +1,56 @@
+void analysis(Int_t first = -1, Int_t last = -1, const char* directory=".")
+{
+
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libANALYSIS");
+
+  AliAODParticle::SetDebug(0);
+  AliRunAnalysis* analysis = new AliRunAnalysis();
+  
+  ::Info("analysis.C","Setting dirs");
+  TObjArray* dirs=0;
+  if ( ((first >= 0) && (last >= 0) ) && ( (last-first)>=0 ) )
+  {//read from many dirs dirs
+    char buff[50];
+    dirs = new TObjArray(last-first+1);
+    for (Int_t i = first; i<=last; i++)
+     {
+//     sprintf(buff,"%s/%s/%s/%05.5d",basedir,field,serie,i);
+       printf("%s/%d\n",directory,i);
+       sprintf(buff,"%s/%d",directory,i);
+       TObjString *odir= new TObjString(buff);
+       dirs->Add(odir);
+     }
+   }
+
+  AliReaderESD* reader = new AliReaderESD(dirs);
+  reader->ReadSimulatedData(kTRUE);
+  reader->SetReadMostProbableOnly(kTRUE);
+
+/*  
+   //example PID cuts
+  AliAODParticleCut* partcut = new  AliAODParticleCut();
+  partcut->SetPID(kPiPlus);//here we define the incarnation 
+  AliAODPIDCut* pidcut = new AliAODPIDCut(kPiPlus,0.5);//accept all particles types that have PID prob > 50%
+  partcut->AddBasePartCut(pidcut);//
+  reader->AddParticleCut(partcut);//This guy makes a copy of a cut for himself so we can modify it here
+  
+  partcut->SetPID(kPiMinus);//here we define that particle has incarnation PiMinus
+  pidcut->SetPID(kPiMinus);//here we define to check if PID probability of being kPiMinus is greater thann 0.5 (number defined few lines above)
+  reader->AddParticleCut(partcut);
+  
+  pidcut->SetPID(kKPlus);
+  pidcut->SetPID(kKPlus);
+  reader->AddParticleCut(partcut);
+  
+  pidcut->SetPID(kKMinus);
+  pidcut->SetPID(kKMinus);
+  reader->AddParticleCut(partcut);
+*/  
+  
+  AliFlowAnalysis* flow = new AliFlowAnalysis();
+  analysis->SetReader(reader);
+  
+  analysis->Add(flow);
+  analysis->Run();
+  delete analysis;
+}
index ac70c18..a336363 100644 (file)
@@ -1,7 +1,13 @@
-SRCS= AliRunAnalysis.cxx AliAnalysis.cxx AliEventCut.cxx \
-      AliD0toKpi.cxx  AliD0toKpiAnalysis.cxx AliFlowAnalysis.cxx
-
+SRCS= AliAOD.cxx \
+      AliAODParticle.cxx AliAODStdParticle.cxx \
+      AliAODParticleCut.cxx AliAODRun.cxx \
+      AliRunAnalysis.cxx AliAnalysis.cxx AliEventCut.cxx \
+      AliReader.cxx AliReaderESD.cxx\
+      AliTrackPoints.cxx AliClusterMap.cxx \
+      AliD0toKpi.cxx  AliD0toKpiAnalysis.cxx AliFlowAnalysis.cxx \
 
 HDRS= $(SRCS:.cxx=.h)
 
 DHDR:=ANALYSISLinkDef.h
+
+EINCLUDE:= TPC CONTAINERS ITS