#ifdef __CINT__
-#pragma link off all globals;
+#pragma link off all glols;
#pragma link off all classes;
#pragma link off all functions;
#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
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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)
+
--- /dev/null
+#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
--- /dev/null
+#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))
+}
+/******************************************************************/
--- /dev/null
+#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
--- /dev/null
+#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);
+
+}
+/**************************************************************************/
+
--- /dev/null
+#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
--- /dev/null
+#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);
+// }
+//}
--- /dev/null
+#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
//
// 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
//
///////////////////////////////////////////////////////////
ClassImp(AliAnalysis)
-Int_t AliAnalysis::fgkDebug = 0;
-
AliAnalysis::AliAnalysis()
{
//ctor
#include <TTask.h>
-class AliESD;
+class AliAOD;
class AliStack;
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)
};
#include "TObject.h"
-class AliESD;
+class AliAOD;
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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
/*********************************************************/
-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;
}
#include "TObject.h"
-class AliESD;
+class AliAOD;
class TObjArray;
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;
#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;
}
}
/*********************************************************/
+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;
}
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;
}
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;
}
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;
}
v2 = TMath::Hypot(ssin,scos);
}
+/*********************************************************/
#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)
};
--- /dev/null
+#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);
+// }
+//}
--- /dev/null
+#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
--- /dev/null
+#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");
+ }
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
#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
}
AliRunAnalysis::~AliRunAnalysis()
{
//dtor
- delete fDirs;
- delete fAnalysies;
+ delete fEventCut;
+ delete fReader;
delete fEventCut;
}
/*********************************************************/
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;
}
//
///////////////////////////////////////////////////////////
-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
--- /dev/null
+//_________________________________
+////////////////////////////////////////////////////////////
+// //
+// 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;
+}
--- /dev/null
+#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
--- /dev/null
+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;
+}
-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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff