From a5556ea5dcb96ed6c8fb42e9637b4e5992497f70 Mon Sep 17 00:00:00 2001 From: skowron Date: Mon, 5 Apr 2004 06:12:30 +0000 Subject: [PATCH] AliAOD and AODParticle (T.Kuhr) - Readers, AODStdParticle and Cuts (P.Skowronski) implemented + test macro. First attemts to analysis package --- ANALYSIS/ANALYSISLinkDef.h | 36 +- ANALYSIS/AliAOD.cxx | 76 ++++ ANALYSIS/AliAOD.h | 41 ++ ANALYSIS/AliAODParticle.cxx | 39 ++ ANALYSIS/AliAODParticle.h | 81 ++++ ANALYSIS/AliAODParticleCut.cxx | 487 ++++++++++++++++++++++ ANALYSIS/AliAODParticleCut.h | 390 ++++++++++++++++++ ANALYSIS/AliAODRun.cxx | 86 ++++ ANALYSIS/AliAODRun.h | 94 +++++ ANALYSIS/AliAODStdParticle.cxx | 410 +++++++++++++++++++ ANALYSIS/AliAODStdParticle.h | 161 ++++++++ ANALYSIS/AliAnalysis.cxx | 10 +- ANALYSIS/AliAnalysis.h | 9 +- ANALYSIS/AliBaseEventCut.h | 6 +- ANALYSIS/AliClusterMap.cxx | 192 +++++++++ ANALYSIS/AliClusterMap.h | 42 ++ ANALYSIS/AliEventCut.cxx | 10 +- ANALYSIS/AliEventCut.h | 4 +- ANALYSIS/AliFlowAnalysis.cxx | 126 +++++- ANALYSIS/AliFlowAnalysis.h | 18 +- ANALYSIS/AliParticle.cxx | 411 +++++++++++++++++++ ANALYSIS/AliParticle.h | 156 ++++++++ ANALYSIS/AliReader.cxx | 434 ++++++++++++++++++++ ANALYSIS/AliReader.h | 101 +++++ ANALYSIS/AliReaderESD.cxx | 711 +++++++++++++++++++++++++++++++++ ANALYSIS/AliReaderESD.h | 130 ++++++ ANALYSIS/AliRunAnalysis.cxx | 256 ++---------- ANALYSIS/AliRunAnalysis.h | 32 +- ANALYSIS/AliTrackPoints.cxx | 589 +++++++++++++++++++++++++++ ANALYSIS/AliTrackPoints.h | 49 +++ ANALYSIS/analysis.C | 56 +++ ANALYSIS/libANALYSIS.pkg | 12 +- 32 files changed, 4986 insertions(+), 269 deletions(-) create mode 100644 ANALYSIS/AliAOD.cxx create mode 100644 ANALYSIS/AliAOD.h create mode 100644 ANALYSIS/AliAODParticle.cxx create mode 100644 ANALYSIS/AliAODParticle.h create mode 100644 ANALYSIS/AliAODParticleCut.cxx create mode 100644 ANALYSIS/AliAODParticleCut.h create mode 100644 ANALYSIS/AliAODRun.cxx create mode 100644 ANALYSIS/AliAODRun.h create mode 100644 ANALYSIS/AliAODStdParticle.cxx create mode 100644 ANALYSIS/AliAODStdParticle.h create mode 100644 ANALYSIS/AliClusterMap.cxx create mode 100644 ANALYSIS/AliClusterMap.h create mode 100644 ANALYSIS/AliParticle.cxx create mode 100644 ANALYSIS/AliParticle.h create mode 100644 ANALYSIS/AliReader.cxx create mode 100644 ANALYSIS/AliReader.h create mode 100644 ANALYSIS/AliReaderESD.cxx create mode 100644 ANALYSIS/AliReaderESD.h create mode 100644 ANALYSIS/AliTrackPoints.cxx create mode 100644 ANALYSIS/AliTrackPoints.h create mode 100644 ANALYSIS/analysis.C diff --git a/ANALYSIS/ANALYSISLinkDef.h b/ANALYSIS/ANALYSISLinkDef.h index 0f37990eb93..8d59d57f353 100644 --- a/ANALYSIS/ANALYSISLinkDef.h +++ b/ANALYSIS/ANALYSISLinkDef.h @@ -1,6 +1,6 @@ #ifdef __CINT__ -#pragma link off all globals; +#pragma link off all glols; #pragma link off all classes; #pragma link off all functions; @@ -9,9 +9,43 @@ #pragma link C++ class AliRunAnalysis+; #pragma link C++ class AliAnalysis+; +#pragma link C++ class AliAOD+; +#pragma link C++ class AliAODParticle+; +#pragma link C++ class AliAODStdParticle+; + +#pragma link C++ class AliAODRun+; + +#pragma link C++ class AliTrackPoints+; +#pragma link C++ class AliClusterMap+; + +#pragma link C++ class AliReader+; +#pragma link C++ class AliReaderESD+; + + #pragma link C++ class AliFlowAnalysis+; #pragma link C++ class AliEventCut+; +#pragma link C++ class AliAODParticleCut-; +#pragma link C++ class AliAODEmptyParticleCut-; +#pragma link C++ class AliAODBaseCut+; + +#pragma link C++ class AliAODMomentumCut+; +#pragma link C++ class AliAODPtCut+; +#pragma link C++ class AliAODEnergyCut+; +#pragma link C++ class AliAODRapidityCut+; +#pragma link C++ class AliAODPseudoRapidityCut+; +#pragma link C++ class AliAODPxCut+; +#pragma link C++ class AliAODPyCut+; +#pragma link C++ class AliAODPzCut+; +#pragma link C++ class AliAODPhiCut+; +#pragma link C++ class AliAODThetaCut+; +#pragma link C++ class AliAODVxCut+; +#pragma link C++ class AliAODVyCut+; +#pragma link C++ class AliAODVzCut+; +#pragma link C++ class AliAODPIDCut+; +#pragma link C++ class AliAODLogicalOperCut-; +#pragma link C++ class AliAODAndCut+; +#pragma link C++ class AliAODOrCut+; #endif diff --git a/ANALYSIS/AliAOD.cxx b/ANALYSIS/AliAOD.cxx new file mode 100644 index 00000000000..7e3de102970 --- /dev/null +++ b/ANALYSIS/AliAOD.cxx @@ -0,0 +1,76 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +///////////////////////////////////////////////////////////// +// +// base class for AOD containers +// +///////////////////////////////////////////////////////////// + +#include +#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 +#include +#include "AliAODParticle.h" + +class TParticle; + +class AliAOD: public TObject { +public: + AliAOD(){} + virtual ~AliAOD() { Reset(); } + + virtual TObjArray* GetParticles() {return &fParticles;}; + virtual Int_t GetNumberOfParticles() const {return fParticles.GetEntriesFast();} + virtual AliAODParticle* GetParticle(Int_t index) const {return (AliAODParticle*) fParticles[index];} + virtual void AddParticle(AliAODParticle* particle) {fParticles.Add(particle);}; + virtual void AddParticle(TParticle* part, Int_t idx); //adds particle to the event + virtual void AddParticle(Int_t pdg, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot, + Double_t vx, Double_t vy, Double_t vz, Double_t time); + + virtual void Reset(); + void SwapParticles(Int_t i, Int_t j);//swaps particles positions; used by AliReader::Blend +private: + TObjArray fParticles; // array of AOD particles, AliAOD is owner of particles + + ClassDef(AliAOD,1) // base class for AOD containers +}; + +#endif diff --git a/ANALYSIS/AliAODParticle.cxx b/ANALYSIS/AliAODParticle.cxx new file mode 100644 index 00000000000..36cf642fdb5 --- /dev/null +++ b/ANALYSIS/AliAODParticle.cxx @@ -0,0 +1,39 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +///////////////////////////////////////////////////////////// +// +// base class for AOD particles +// +// Most of the methods defined in this base class (or even all of them) +// should be redefined in derived classes. Some methods like Pt() or +// Vx() have default implementation, but they are not optimised for +// performance. It is likely that they can be (and should be) implemented +// in the derived classes in a way which gives faster access. +// +// Algorithms and analysis classes should as much as possible use the +// interface of the base class instead of special methods of derived +// classes. This allows to run the code on all types of particles. +// +///////////////////////////////////////////////////////////// + +#include "AliAODParticle.h" + +Int_t AliAODParticle::fgDebug = 0; + +ClassImp(AliAODParticle) + diff --git a/ANALYSIS/AliAODParticle.h b/ANALYSIS/AliAODParticle.h new file mode 100644 index 00000000000..cc1d0325ae6 --- /dev/null +++ b/ANALYSIS/AliAODParticle.h @@ -0,0 +1,81 @@ +#ifndef ALIAODPARTICLE_H +#define ALIAODPARTICLE_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +///////////////////////////////////////////////////////////// +// +// base class for AOD particles +// +///////////////////////////////////////////////////////////// + +#include +#include +#include + +#include "AliAnalysis.h" + +class AliTrackPoints; +class AliClusterMap; + +class AliAODParticle : public TObject { +public: + AliAODParticle(){} + virtual ~AliAODParticle(){} + // kinematics + virtual TLorentzVector FourMomentum() const = 0; + virtual TVector3 Momentum() const {return FourMomentum().Vect();}; + virtual Double_t Mass() const {return FourMomentum().M();}; + virtual Double_t E() const {return FourMomentum().E();}; + virtual Double_t P() const {return FourMomentum().P();}; + virtual Double_t Pt() const {return FourMomentum().Pt();}; + virtual Double_t Px() const {return FourMomentum().Px();}; + virtual Double_t Py() const {return FourMomentum().Py();}; + virtual Double_t Pz() const {return FourMomentum().Pz();}; + virtual Double_t Phi() const {return FourMomentum().Phi();}; + virtual Double_t Theta() const {return FourMomentum().Theta();}; + virtual Double_t Eta() const {return FourMomentum().Eta();}; + virtual Double_t Y() const {return FourMomentum().Rapidity();}; + + // PID + virtual Double_t Charge() const = 0; + virtual Double_t GetProbability(Int_t pdg) const = 0; + virtual Int_t GetMostProbable() const = 0; + + virtual Int_t GetPdgCode() const = 0;//We need to assume some PID (f.e. energy calculation) + //sotimes one track can apear in analysis twise (e.g. ones as pion ones as kaon) + + // vertices + virtual TVector3 ProductionVertex() const = 0; + virtual Double_t Vx() const {return ProductionVertex().X();}; + virtual Double_t Vy() const {return ProductionVertex().Y();}; + virtual Double_t Vz() const {return ProductionVertex().Z();}; + virtual AliAODParticle* Mother() const {return NULL;}; + virtual Bool_t HasDecayVertex() const {return kFALSE;}; + virtual TVector3 DecayVertex() const {return TVector3();}; + virtual Int_t NumberOfDaughters() const {return 0;}; + virtual AliAODParticle* Daughter(Int_t /*index*/) const {return NULL;}; + + + // type information + virtual Bool_t IsSimulated() {return kFALSE;}; + virtual Bool_t IsTrack() {return kFALSE;}; + virtual Bool_t IsCluster() {return kFALSE;}; + + //HBT specific + virtual AliTrackPoints* GetTrackPoints() const {return 0x0;} + virtual AliClusterMap* GetClusterMap() const {return 0x0;} + virtual void Print() const = 0; + + static void SetDebug(Int_t dbg=1){fgDebug=dbg;} + static Int_t GetDebug(){return fgDebug;} + +private: + static Int_t fgDebug;//! debug level for all the analysis package + + ClassDef(AliAODParticle,1) // base class for AOD particles +}; + +#endif diff --git a/ANALYSIS/AliAODParticleCut.cxx b/ANALYSIS/AliAODParticleCut.cxx new file mode 100644 index 00000000000..aa4a40e6e41 --- /dev/null +++ b/ANALYSIS/AliAODParticleCut.cxx @@ -0,0 +1,487 @@ +#include "AliAODParticleCut.h" +//__________________________________________________________________________ +//////////////////////////////////////////////////////////////////////////// +// // +// class AliAODParticleCut // +// // +// Classes for single particle cuts // +// User should use only AliAODParticleCut, eventually // +// EmptyCut which passes all particles // +// There is all interface for setting cuts on all particle properties // +// The main method is Pass - which returns // +// True to reject particle // +// False in case it meets all the criteria of the given cut // +// // +// User should create (and also destroy) cuts himself // +// and then pass them to the Analysis And Function by a proper method // +// // +// // +// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html // +// resonsible: Piotr Skowronski@cern.ch // +// // +//////////////////////////////////////////////////////////////////////////// + +#include + + +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;iClone();//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;iClone();//create new object (clone) and rember pointer to it + } + return *this; +} + +/******************************************************************/ +AliAODParticleCut::~AliAODParticleCut() +{ + //dtor + for (Int_t i = 0;iGetPdgCode() != fPID) && ( fPID != 0)) return kTRUE; + + for (Int_t i = 0;iPass(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;iGetProperty() == 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> fPID; + b >> fNCuts; + for (i = 0;i> 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;iPrint(); + } +} + +/******************************************************************/ +/******************************************************************/ + +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="<Clone():0x0), + fSecond((second)?(AliAODBaseCut*)second->Clone():0x0) +{ + //ctor + if ( (fFirst && fSecond) == kFALSE) + { + Fatal("AliAODLogicalOperCut","One of parameters is NULL!"); + } +} +/******************************************************************/ + +AliAODLogicalOperCut::~AliAODLogicalOperCut() +{ + //destructor + delete fFirst; + delete fSecond; +} +/******************************************************************/ + +Bool_t AliAODLogicalOperCut::AliAODDummyBaseCut::Pass(AliAODParticle* /*part*/) const +{ + //checks if particles passes properties defined by this cut + Warning("Pass","You are using dummy base cut! Probobly some logical cut is not set up properly"); + return kFALSE;//accept +} +/******************************************************************/ + +void AliAODLogicalOperCut::Streamer(TBuffer &b) +{ + // Stream all objects in the array to or from the I/O buffer. + UInt_t R__s, R__c; + if (b.IsReading()) + { + delete fFirst; + delete fSecond; + fFirst = 0x0; + fSecond = 0x0; + + b.ReadVersion(&R__s, &R__c); + TObject::Streamer(b); + b >> fFirst; + b >> fSecond; + b.CheckByteCount(R__s, R__c,AliAODLogicalOperCut::IsA()); + } + else + { + R__c = b.WriteVersion(AliAODLogicalOperCut::IsA(), kTRUE); + TObject::Streamer(b); + b << fFirst; + b << fSecond; + b.SetByteCount(R__c, kTRUE); + } +} + +/******************************************************************/ +ClassImp(AliAODOrCut) + +Bool_t AliAODOrCut::Pass(AliAODParticle * p) const +{ + //returns true when rejected + //AND operation is a little bit misleading but is correct + //User wants to build logical cuts with natural (positive) logic + //while AODAN use inernally reverse (returns true when rejected) + if (fFirst->Pass(p) && fSecond->Pass(p)) return kTRUE;//rejected (both rejected, returned kTRUE) + return kFALSE;//accepted, at least one accepted (returned kFALSE) +} +/******************************************************************/ + +ClassImp(AliAODAndCut) + +Bool_t AliAODAndCut::Pass(AliAODParticle * p) const +{ + //returns true when rejected + //OR operation is a little bit misleading but is correct + //User wants to build logical cuts with natural (positive) logic + //while AODAN use inernally reverse (returns true when rejected) + if (fFirst->Pass(p) || fSecond->Pass(p)) return kTRUE;//rejected (any of two rejected(returned kTRUE) ) + return kFALSE;//accepted (both accepted (returned kFALSE)) +} +/******************************************************************/ diff --git a/ANALYSIS/AliAODParticleCut.h b/ANALYSIS/AliAODParticleCut.h new file mode 100644 index 00000000000..ae7aac2e47d --- /dev/null +++ b/ANALYSIS/AliAODParticleCut.h @@ -0,0 +1,390 @@ +#ifndef ALIAODPARTICLECUT_H +#define ALIAODPARTICLECUT_H +//__________________________________________________________________________ +//////////////////////////////////////////////////////////////////////////// +// // +// class AliAODParticleCut // +// // +// Classes for single particle cuts // +// User should use only AliAODParticleCut, eventually // +// EmptyCut which passes all particles // +// There is all interface for setting cuts on all particle properties // +// The main method is Pass - which returns // +// True to reject particle // +// False in case it meets all the criteria of the given cut // +// // +// User should create (and also destroy) cuts himself // +// and then pass them to the Analysis And Function by a proper method // +// // +// // +// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html // +// resonsible: Piotr Skowronski@cern.ch // +// // +//////////////////////////////////////////////////////////////////////////// + + +#include +#include "AliAODParticle.h" + + +class AliAODEmptyParticleCut; +class AliAODParticleCut; +class AliAODBaseCut; + + +/******************************************************************/ +/******************************************************************/ +/******************************************************************/ + +enum AliAODCutProperty + { +//codes particle property + kAODP, //Momentum + kAODPt, //Transverse momentum + kAODE, //Energy + kAODRapidity, // + kAODPseudoRapidity, + kAODPx, //X coAnddinate of the momentum + kAODPy, //Y coAnddinate of the momentum + kAODPz, //Z coAnddinate of the momentum + kAODPhi,//angle + kAODTheta,//angle + kAODVx, // vertex X coAnddinate + kAODVy, // vertex Y coAnddinate + kAODVz, // vertex Z coAnddinate + kAODPid, // vertex Z coAnddinate +//_____________________________ + kAODNone + }; + +/******************************************************************/ +/******************************************************************/ +/******************************************************************/ + +class AliAODParticleCut: public TObject +{ +//Class describing cut on particle + public: + AliAODParticleCut(); + AliAODParticleCut(const AliAODParticleCut& in); + virtual ~AliAODParticleCut(); + AliAODParticleCut& operator = (const AliAODParticleCut& in); + + virtual Bool_t Pass(AliAODParticle* p) const; + Bool_t IsEmpty() const {return kFALSE;} + + void AddBasePartCut(AliAODBaseCut* basecut); + + Int_t GetPID() const { return fPID;} + void SetPID(Int_t pid){fPID=pid;} + void SetMomentumRange(Double_t min, Double_t max); + void SetPRange(Double_t min, Double_t max){SetMomentumRange(min,max);} + void SetPtRange(Double_t min, Double_t max); + void SetEnergyRange(Double_t min, Double_t max); + void SetRapidityRange(Double_t min, Double_t max); + void SetYRange(Double_t min, Double_t max){SetRapidityRange(min,max);} + void SetPseudoRapidityRange(Double_t min, Double_t max); + void SetPxRange(Double_t min, Double_t max); + void SetPyRange(Double_t min, Double_t max); + void SetPzRange(Double_t min, Double_t max); + void SetPhiRange(Double_t min, Double_t max); + void SetThetaRange(Double_t min, Double_t max); + void SetVxRange(Double_t min, Double_t max); + void SetVyRange(Double_t min, Double_t max); + void SetVzRange(Double_t min, Double_t max); + + void Print(void) const; + protected: + + AliAODBaseCut* FindCut(AliAODCutProperty property); + + AliAODBaseCut ** fCuts;//! Array with cuts + Int_t fNCuts; //number of base cuts stored in fCuts + + Int_t fPID; //particle PID - if=0 (rootino) all pids are accepted + + private: + static const Int_t fgkMaxCuts; //Size of the fCuts array + + ClassDef(AliAODParticleCut,1) +}; +/******************************************************************/ +/******************************************************************/ +/******************************************************************/ + +class AliAODEmptyParticleCut: public AliAODParticleCut +{ +//Empty - it passes possitively all particles - it means returns always False +//Class describing cut on particles + public: + AliAODEmptyParticleCut(){}; + virtual ~AliAODEmptyParticleCut(){}; + + Bool_t Pass(AliAODParticle*) const {return kFALSE;} //accept everything + Bool_t IsEmpty() const {return kTRUE;} + + ClassDef(AliAODEmptyParticleCut,1) + +}; + +/******************************************************************/ +/******************************************************************/ +/******************************************************************/ + +class AliAODBaseCut: public TObject + { + //This class defines the range of some property - pure virtual + //Property is coded by AliAODCutTypes type + + public: + + AliAODBaseCut(Double_t min = 0.0, Double_t max = 0.0,AliAODCutProperty prop = kAODNone): + fProperty(prop),fMin(min),fMax(max){} + + virtual ~AliAODBaseCut(){} + + virtual Bool_t Pass(AliAODParticle *p) const; + + void SetRange(Double_t min, Double_t max){fMin = min; fMax = max;} + + void SetMinimum(Double_t min){fMin = min;} + void SetMaximum(Double_t max){fMax = max;} + + Double_t GetMinimum() const {return fMin;} + Double_t GetMaximum() const {return fMax;} + + AliAODCutProperty GetProperty() const {return fProperty;} + virtual void Print(void) const; + + protected: + virtual Double_t GetValue(AliAODParticle *) const = 0; + + AliAODCutProperty fProperty; //property that this cut describes + Double_t fMin;//minimum value + Double_t fMax;//maximum value + + private: + void PrintProperty(void) const; + ClassDef(AliAODBaseCut,1) + + }; + +inline Bool_t +AliAODBaseCut::Pass(AliAODParticle *p) const +{ + //cjecks if particle property fits in range + if ( (GetValue(p) < fMin) || (GetValue(p) > fMax ) ) return kTRUE; //rejected + else return kFALSE; //accepted +} +/******************************************************************/ +/******************************************************************/ +/******************************************************************/ + + +class AliAODMomentumCut: public AliAODBaseCut + { + public: + AliAODMomentumCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODP){} + virtual ~AliAODMomentumCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->P();} + ClassDef(AliAODMomentumCut,1) + }; + +class AliAODPtCut: public AliAODBaseCut + { + public: + AliAODPtCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPt){} + virtual ~AliAODPtCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Pt();} + ClassDef(AliAODPtCut,1) + }; + + +class AliAODEnergyCut: public AliAODBaseCut + { + public: + AliAODEnergyCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODE){} + virtual ~AliAODEnergyCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const {return p->E();} + ClassDef(AliAODEnergyCut,1) + }; + +class AliAODRapidityCut: public AliAODBaseCut + { + public: + AliAODRapidityCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODRapidity){} + virtual ~AliAODRapidityCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Y();} + ClassDef(AliAODRapidityCut,1) + }; + +class AliAODPseudoRapidityCut: public AliAODBaseCut + { + public: + AliAODPseudoRapidityCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPseudoRapidity){} + virtual ~AliAODPseudoRapidityCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Eta();} + ClassDef(AliAODPseudoRapidityCut,1) + }; + +class AliAODPxCut: public AliAODBaseCut + { + public: + AliAODPxCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPx){} + virtual ~AliAODPxCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Px();} + ClassDef(AliAODPxCut,1) + }; + +class AliAODPyCut: public AliAODBaseCut + { + public: + AliAODPyCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPy){} + virtual ~AliAODPyCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Py();} + ClassDef(AliAODPyCut,1) + }; + + +class AliAODPzCut: public AliAODBaseCut + { + public: + AliAODPzCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPz){} + virtual ~AliAODPzCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Pz();} + ClassDef(AliAODPzCut,1) + }; + +class AliAODPhiCut: public AliAODBaseCut + { + public: + AliAODPhiCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODPhi){} + virtual ~AliAODPhiCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Phi();} + ClassDef(AliAODPhiCut,1) + + }; + +class AliAODThetaCut: public AliAODBaseCut + { + public: + AliAODThetaCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODTheta){} + virtual ~AliAODThetaCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Theta();} + ClassDef(AliAODThetaCut,1) + + }; + +class AliAODVxCut: public AliAODBaseCut + { + //Cut of the X coAnddinate of the vertex position + public: + AliAODVxCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODVx){} + virtual ~AliAODVxCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Vx();} //retruns value of the vertex + ClassDef(AliAODVxCut,1) + + }; + + +class AliAODVyCut: public AliAODBaseCut + { + //Cut of the X coAnddinate of the vertex position + public: + AliAODVyCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODVy){} + virtual ~AliAODVyCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Vy();} //retruns value of the vertex + ClassDef(AliAODVyCut,1) + + }; + +class AliAODVzCut: public AliAODBaseCut + { + //Cut of the X coAnddinate of the vertex position + public: + AliAODVzCut(Double_t min = 0.0, Double_t max = 0.0):AliAODBaseCut(min,max,kAODVz){} + virtual ~AliAODVzCut(){} + protected: + Double_t GetValue(AliAODParticle * p)const{return p->Vz();} //retruns value of the vertex + + ClassDef(AliAODVzCut,1) + + }; + +class AliAODPIDCut: public AliAODBaseCut + { + public: + AliAODPIDCut():AliAODBaseCut(0.0,0.0,kAODPid),fPID(0){} + AliAODPIDCut(Int_t pid, Double_t min = 0.0, Double_t max = 1.0):AliAODBaseCut(min,max,kAODPid),fPID(pid){} + virtual ~AliAODPIDCut(){} + + void SetPID(Int_t pid){fPID = pid;} + void Print(void) const; + protected: + Double_t GetValue(AliAODParticle * p)const{return p->GetProbability(fPID);} + Int_t fPID; //pid of particle that the pid is set + ClassDef(AliAODPIDCut,1) + }; +//___________________________________________________ +///////////////////////////////////////////////////// +// // +// class AliAODLogicalOperCut // +// // +// This cut is base class fAnd class that perfAndms // +// logical operations on cuts // +// // +///////////////////////////////////////////////////// +class AliAODLogicalOperCut: public AliAODBaseCut + { + public: + AliAODLogicalOperCut(); + AliAODLogicalOperCut(AliAODBaseCut* first, AliAODBaseCut* second); + virtual ~AliAODLogicalOperCut(); + protected: + Double_t GetValue(AliAODParticle * /*part*/) const {MayNotUse("GetValue");return 0.0;} + + AliAODBaseCut* fFirst; //second cut + AliAODBaseCut* fSecond; //first cut + private: + class AliAODDummyBaseCut: public AliAODBaseCut + { + Double_t GetValue(AliAODParticle * /*part*/) const {return 0.0;} + Bool_t Pass(AliAODParticle* /*part*/) const; + }; + + ClassDef(AliAODLogicalOperCut,1) + }; + +class AliAODOrCut: public AliAODLogicalOperCut +{ + public: + AliAODOrCut(){} + AliAODOrCut(AliAODBaseCut* first, AliAODBaseCut* second):AliAODLogicalOperCut(first,second){} + virtual ~AliAODOrCut(){} + Bool_t Pass(AliAODParticle *p) const; + ClassDef(AliAODOrCut,1) +}; + +class AliAODAndCut: public AliAODLogicalOperCut +{ + public: + AliAODAndCut(){} + AliAODAndCut(AliAODBaseCut* first, AliAODBaseCut* second):AliAODLogicalOperCut(first,second){} + virtual ~AliAODAndCut(){} + Bool_t Pass(AliAODParticle *p) const; + ClassDef(AliAODAndCut,1) +}; + +#endif diff --git a/ANALYSIS/AliAODRun.cxx b/ANALYSIS/AliAODRun.cxx new file mode 100644 index 00000000000..1f962743910 --- /dev/null +++ b/ANALYSIS/AliAODRun.cxx @@ -0,0 +1,86 @@ +#include "AliAODRun.h" +//____________________ +/////////////////////////////////////////////////////// +// // +// AliAODRun // +// // +// Class storing and managing events // +// // +// Piotr.Skowronski@cern.ch // +// http://alisoft.cern.ch/people/skowron/analyzer // +// // +/////////////////////////////////////////////////////// + +#include + +ClassImp(AliAODRun) +/**************************************************************************/ + +AliAODRun::AliAODRun() +{ + //contructor + fEvents = new TObjArray();//create array for AliAODs + if(!fEvents) Fatal("AliAODRun::AliAODRun","Can not allocate memory"); + fEvents->SetOwner(); //array is an owner: when is deleted or cleared it deletes objects that it contains +} +/**************************************************************************/ + +AliAODRun::~AliAODRun() +{ + //destructor + delete fEvents;//delete array with events +} +/**************************************************************************/ + +void AliAODRun::Reset() + { + fEvents->Clear();//clear an array with events. + //All events are deleted because array is an owner (set in constructor) + } +/**************************************************************************/ + +void AliAODRun::AddParticle(Int_t event, AliAODParticle* part) +{ + //Adds particle to event + //if there is no event of this number, crate it and add to the collection + if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliAOD, event); + + GetEvent(event)->AddParticle(part); +} +/**************************************************************************/ + +void AliAODRun::AddParticle(Int_t event, TParticle* part, Int_t idx) +{ + //if there is no event of this number, crate it and add to the collection + if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliAOD, event); + GetEvent(event)->AddParticle(part,idx); +} +/**************************************************************************/ + +void AliAODRun::AddParticle(Int_t event, Int_t pdg, Int_t idx, + Double_t px, Double_t py, Double_t pz, Double_t etot, + Double_t vx, Double_t vy, Double_t vz, Double_t time) +{ + //if there is no event of this number, crate it and add to the collection + if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliAOD, event); + GetEvent(event)->AddParticle(pdg,idx,px,py,pz,etot,vx,vy,vz,time); +} +/**************************************************************************/ + +void AliAODRun::SetEvent(Int_t number, AliAOD* event) +{ + //adds an event to the run + if (event == 0x0) + { + delete fEvents->RemoveAt(number); + return; + } + AliAOD* ev = GetEvent(number); + if (ev == event) return; + + delete fEvents->RemoveAt(number); + fEvents->AddAtAndExpand(event, number); + +} +/**************************************************************************/ + diff --git a/ANALYSIS/AliAODRun.h b/ANALYSIS/AliAODRun.h new file mode 100644 index 00000000000..bbd7d817db3 --- /dev/null +++ b/ANALYSIS/AliAODRun.h @@ -0,0 +1,94 @@ +#ifndef ALIAODRUN_H +#define ALIAODRUN_H +//____________________ +/////////////////////////////////////////////////////// +// // +// AliAODRun // +// // +// Class storing and managing set events // +// designed for fast acces // +// // +// Piotr.Skowronski@cern.ch // +// http://alisoft.cern.ch/people/skowron/analyzer // +// // +/////////////////////////////////////////////////////// + + +#include "AliAOD.h" +#include + +class AliAODParticle; +class TParticle; + +class AliAODRun: public TObject + { + public: + AliAODRun(); + virtual ~AliAODRun(); + + void AddParticle(Int_t event, AliAODParticle* part); //inerface to AliAOD::AddParticle(AliAODParticle*) + void AddParticle(Int_t event, TParticle* part, Int_t idx);//inerface to AliAOD::AddParticle(TParticle*) + + //inerface to AliAOD::AddParticle(Int_t.Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t) + void AddParticle(Int_t event, Int_t pdg, Int_t idx, + Double_t px, Double_t py, Double_t pz, Double_t etot, + Double_t vx, Double_t vy, Double_t vz, Double_t time); + + void SetEvent(Int_t number, AliAOD* event); + AliAODParticle* GetParticle(Int_t event, Int_t n); //returns nth particle from event + AliAOD* GetEvent(Int_t event) const; //returns AliAOD number "event" + + Int_t GetNumberOfEvents() const; //returns number of events + Int_t GetNumberOfParticlesInEvent(Int_t event) const; //returns number of particles in event number "event" + void Reset();//clears all events in the array (deletes) + protected: + TObjArray* fEvents;//!Array containig AliAODs + private: + + public: + ClassDef(AliAODRun,1) + }; + + +/**************************************************************************/ + +inline +AliAOD* AliAODRun::GetEvent(Int_t event) const +{ +//returns pointer to AliAOD number "event" + //check if array is enough big - protect from error massage from array "Out Of Bounds" + if (event>=fEvents->GetSize()) return 0x0;//WARNING, that line relies + //on index of first object in TObjArray is 0 + //== LowerBound = 0 + return (AliAOD*)fEvents->At(event); +} +/**************************************************************************/ +inline +AliAODParticle* AliAODRun::GetParticle(Int_t event, Int_t n) +{ + //returns nth particle from event number event + AliAOD* e = GetEvent(event); + return (e)?e->GetParticle(n):0x0; +} + +/**************************************************************************/ + +inline +Int_t AliAODRun::GetNumberOfEvents() const + { +//returns number of events in collection + + return fEvents->GetEntriesFast(); //there may be empty slots but we do not care + //Analysis checks it if return is not NULL + } +/**************************************************************************/ + +inline +Int_t AliAODRun::GetNumberOfParticlesInEvent(Int_t event) const +{ +//returns number of Particles in event + AliAOD* e = GetEvent(event); + return (e)?e->GetNumberOfParticles():0x0; +} + +#endif diff --git a/ANALYSIS/AliAODStdParticle.cxx b/ANALYSIS/AliAODStdParticle.cxx new file mode 100644 index 00000000000..628a7cffe08 --- /dev/null +++ b/ANALYSIS/AliAODStdParticle.cxx @@ -0,0 +1,410 @@ +#include "AliAODStdParticle.h" +//___________________________________________________________ +///////////////////////////////////////////////////////////// +// +// class AliAODStdParticle +// +// Ali HBT Particle: simplified class TParticle +// Simplified in order to minimize the size of object +// - we want to keep a lot of such a objects in memory +// Additionaly adjusted for HBT Analysies purposes +// + pointer to Track Points +// + pointer to Cluster Map(s) +// +// Piotr.Skowronski@cern.ch +// +///////////////////////////////////////////////////////////// +#include +#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; iClone(); + 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> fPids[i]; +// } +// for (i = 0;i> 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 +//#include +#include + + +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 diff --git a/ANALYSIS/AliAnalysis.cxx b/ANALYSIS/AliAnalysis.cxx index 5dc2f68c274..1a2fd96f9e3 100644 --- a/ANALYSIS/AliAnalysis.cxx +++ b/ANALYSIS/AliAnalysis.cxx @@ -4,9 +4,11 @@ // // class AliAnalysis // -// Base class for analysis -// -// +// Base class for analysis. +// Each inheriting calss must define 3 methods: +// - Init() : that is called before event processing +// - ProcessEvent(AliESD*,AliStack*) +// - // Piotr.Skowronski@cern.ch // /////////////////////////////////////////////////////////// @@ -14,8 +16,6 @@ ClassImp(AliAnalysis) -Int_t AliAnalysis::fgkDebug = 0; - AliAnalysis::AliAnalysis() { //ctor diff --git a/ANALYSIS/AliAnalysis.h b/ANALYSIS/AliAnalysis.h index e61e8be28da..4b1e572046f 100644 --- a/ANALYSIS/AliAnalysis.h +++ b/ANALYSIS/AliAnalysis.h @@ -14,7 +14,7 @@ #include -class AliESD; +class AliAOD; class AliStack; class AliAnalysis: public TTask @@ -25,17 +25,12 @@ class AliAnalysis: public TTask virtual ~AliAnalysis(); virtual Int_t Init() = 0; - virtual Int_t ProcessEvent(AliESD* esd, AliStack* stack = 0x0) = 0; + virtual Int_t ProcessEvent(AliAOD* aodrec, AliAOD* aodsim = 0x0) = 0; virtual Int_t Finish() = 0; - - static Int_t GetDebug() {return fgkDebug;} - static void SetDebug(Int_t level) {fgkDebug = level;} - protected: private: - static Int_t fgkDebug;//! debug level ClassDef(AliAnalysis,1) }; diff --git a/ANALYSIS/AliBaseEventCut.h b/ANALYSIS/AliBaseEventCut.h index 1561ab28b99..d9923a3f1a4 100644 --- a/ANALYSIS/AliBaseEventCut.h +++ b/ANALYSIS/AliBaseEventCut.h @@ -13,7 +13,7 @@ #include "TObject.h" -class AliESD; +class AliAOD; class AliBaseEventCut: public TObject { @@ -21,9 +21,9 @@ class AliBaseEventCut: public TObject AliBaseEventCut(); virtual ~AliBaseEventCut(){} - virtual Bool_t Pass(AliESD* esd) const;//returns kTRUE if rejected + virtual Bool_t Pass(AliAOD* aod) const;//returns kTRUE if rejected protected: - virtual Double_t GetValue(AliESD* esd) const = 0; + virtual Double_t GetValue(AliAOD* aod) const = 0; Double_t fMin;//Minimum value Double_t fMax;//Maximum value diff --git a/ANALYSIS/AliClusterMap.cxx b/ANALYSIS/AliClusterMap.cxx new file mode 100644 index 00000000000..fafdba536aa --- /dev/null +++ b/ANALYSIS/AliClusterMap.cxx @@ -0,0 +1,192 @@ +#include "AliClusterMap.h" + +//_________________________________________________ +/////////////////////////////////////////////////// +// +// class AliClusterMap +// +// class that describes cluster occupation at TPC +// Each padraw has a corresponding bit in fPadRawMap +// +// +// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html +// Piotr.Skowronski@cern.ch +// +/////////////////////////////////////////////////////////////////////////// + + +#include "AliESDtrack.h" +#include "AliTPCtrack.h" +#include "AliAODParticle.h" +#include +const Int_t AliClusterMap::fNPadRows = 159; + +AliClusterMap::AliClusterMap(): + fPadRawMap(fNPadRows) +{ +//ctor + +} +/***********************************************************************/ +AliClusterMap::AliClusterMap(AliESDtrack* track): + fPadRawMap( (track)?track->GetTPCClusterMap():fNPadRows ) +{ + //ctor + + if (AliAODParticle::GetDebug() > 2) + { + Info("AliClusterMap(AliESDtrack*)",""); + Print(); + } +} +/***********************************************************************/ + +AliClusterMap::AliClusterMap(AliTPCtrack* track): + fPadRawMap(fNPadRows) +{ + //ctor + + //Does not work since indeces in the claster index array + //in the TPC track does not correspond to the padraw segmatation + if (AliAODParticle::GetDebug() > 9) + Info("AliClusterMap", + "#####################################################################"); + if (track == 0x0) + { + Error("AliClusterMap","Pointer to TPC track is NULL"); + return; + } + Int_t prevrow = -1; + Int_t i = 0; + for ( ; i < track->GetNumberOfClusters(); i++) + { + Int_t idx = track->GetClusterIndex(i); + Int_t sect = (idx&0xff000000)>>24; + Int_t row = (idx&0x00ff0000)>>16; + if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors + if (AliAODParticle::GetDebug() > 9) + Info("AliClusterMap","Cl.idx is %d, sect %d, row %d",idx,sect,row); + + fPadRawMap.SetBitNumber(row,kTRUE); + + //Fill the gap between previous row and this row with 0 bits + if (prevrow < 0) + { + prevrow = row;//if previous bit was not assigned yet == this is the first one + } + else + { //we don't know the order (inner to outer or reverse) + //just to be save in case it is going to change + Int_t n = 0, m = 0; + if (prevrow < row) + { + n = prevrow; + m = row; + } + else + { + n = row; + m = prevrow; + } + for (Int_t j = n+1; j < m; j++) + { + fPadRawMap.SetBitNumber(j,kFALSE); + } + prevrow = row; + } + } + + if (AliAODParticle::GetDebug() > 2) + { + Info("AliClusterMap(AliTPCtrack*)",""); + Print(); + } +} +/***********************************************************************/ + +void AliClusterMap::Print() const +{ +//Prints the bit map + TString msg; + for ( Int_t i = 0; i < fNPadRows; i++) + { + if ( fPadRawMap.TestBitNumber(i) ) + { + msg+="1"; + } + else + { + msg+="0"; + } + } + Info("AliClusterMap","BitMap is\n %s",msg.Data()); + +} + +/***********************************************************************/ + +Float_t AliClusterMap::GetOverlapFactor(const AliClusterMap& clmap) const +{ + //Returns quality factor FQ = Sum(An)/Sum(clusters) + // | -1; if both tracks have a cluster on padrow n + //An = < 0; if neither track has a cluster on padrow n + // | 1; if only one trackhas a cluster on padrow n + // Returned value ranges between + // -0.5 (low probability that these tracks are a split track) + // and + // 1.0 (high probability that these tracks are a split track) + TString msg1; + TString msg2; + + Int_t nh = 0; + Int_t an = 0; + for ( Int_t i = 0; i < fNPadRows; i++) + { + Bool_t x = HasClAtPadRow(i); + Bool_t y = clmap.HasClAtPadRow(i); + + if (x) msg1+="1";else msg1+="0"; + if (y) msg2+="1";else msg2+="0"; + + if (x && y)//both have clasters + { + an--; + nh+=2; + } + else + { + + if (x || y)//only one have cluters + { + an++; + nh++; + } + } + } + + + Float_t retval = 0.0; + if (nh > 0) retval = ((Float_t)an)/((Float_t)nh); + else Warning("GetOverlapFactor","Number of counted cluters is 0."); + + if (AliAODParticle::GetDebug() > 2) + { + Info("GetOverlapFactor","Splitting Quality Factor is %f. SumAn = %d, SumClusters %d",retval,an,nh); + if (retval == 1.0) + { + Print(); + Info("AliClusterMap","BitMap is\n %s\n",msg1.Data()); + clmap.Print(); + Info("AliClusterMap","BitMap is\n %s\n\n\n\n",msg2.Data()); + } + if (retval == -.5) + { + Print(); + Info("AliClusterMap","BitMap is\n %s\n",msg1.Data()); + clmap.Print(); + Info("AliClusterMap","BitMap is\n %s\n\n\n\n",msg2.Data()); + } + } + + return retval; +} diff --git a/ANALYSIS/AliClusterMap.h b/ANALYSIS/AliClusterMap.h new file mode 100644 index 00000000000..f7d8774afe0 --- /dev/null +++ b/ANALYSIS/AliClusterMap.h @@ -0,0 +1,42 @@ +#ifndef ALICLUSTERMAP_H +#define ALICLUSTERMAP_H +//_________________________________________________ +/////////////////////////////////////////////////// +// +// class AliClusterMap +// +// class that describes cluster occupation at TPC +// Each padraw has a corresponding bit in fPadRawMap +// +// +// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html +// Piotr.Skowronski@cern.ch +// +/////////////////////////////////////////////////////////////////////////// + + +#include +#include + +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 diff --git a/ANALYSIS/AliEventCut.cxx b/ANALYSIS/AliEventCut.cxx index e7f1eb0216b..32e943f421a 100644 --- a/ANALYSIS/AliEventCut.cxx +++ b/ANALYSIS/AliEventCut.cxx @@ -30,14 +30,20 @@ AliEventCut::~AliEventCut() /*********************************************************/ -Bool_t AliEventCut::Pass(AliESD* esd) const +Bool_t AliEventCut::Pass(AliAOD* aod) const { //returns kTRUE if rejected + if (aod == 0x0) + { + Error("Pass","Pointer to AOD is NULL. Not passed the cut"); + return kFALSE; + } + TIter iter(fBaseCuts); AliBaseEventCut* becut; while (( becut = (AliBaseEventCut*)iter() )) { - if (becut->Pass(esd)) return kTRUE; + if (becut->Pass(aod)) return kTRUE; } return kFALSE; } diff --git a/ANALYSIS/AliEventCut.h b/ANALYSIS/AliEventCut.h index 6ab0369b86d..21e822b6337 100644 --- a/ANALYSIS/AliEventCut.h +++ b/ANALYSIS/AliEventCut.h @@ -12,7 +12,7 @@ #include "TObject.h" -class AliESD; +class AliAOD; class TObjArray; class AliEventCut: public TObject @@ -21,7 +21,7 @@ class AliEventCut: public TObject AliEventCut(); virtual ~AliEventCut(); - virtual Bool_t Pass(AliESD* esd) const;//returns kTRUE if rejected + virtual Bool_t Pass(AliAOD* aod) const;//returns kTRUE if rejected protected: TObjArray* fBaseCuts; diff --git a/ANALYSIS/AliFlowAnalysis.cxx b/ANALYSIS/AliFlowAnalysis.cxx index 36d029e897a..2d5184e26c6 100644 --- a/ANALYSIS/AliFlowAnalysis.cxx +++ b/ANALYSIS/AliFlowAnalysis.cxx @@ -17,26 +17,54 @@ #include #include +#include +#include +#include + #include #include ClassImp(AliFlowAnalysis) +AliFlowAnalysis::AliFlowAnalysis(): + fPartCut(0x0) +{ + //ctor +} +/*********************************************************/ + +AliFlowAnalysis::~AliFlowAnalysis() +{ + //dtor + delete fPartCut; +} +/*********************************************************/ + Int_t AliFlowAnalysis::Init() { //Initilizes anaysis Info("Init",""); return 0; } - /*********************************************************/ -Int_t AliFlowAnalysis::ProcessEvent(AliESD* esd, AliStack* stack) +Int_t AliFlowAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim) { - Info("ProcessEvent","Stack address %#x",stack); - Double_t psi = GetEventPlane(esd); - Info("ProcessEvent","Event plane is %f",psi); + Info("ProcessEvent","Sim AOD address %#x",aodsim); + Double_t psi = 0, v2 = 0; + if (aodrec) + { + GetFlow(aodrec,v2,psi); + Info("ProcessEvent","Reconstructed Event: Event plane is %f, V2 is %f",psi,v2); + } + + if (aodsim) + { + GetFlow(aodsim,v2,psi); + Info("ProcessEvent","Simulated Event: Event plane is %f, V2 is %f",psi,v2); + } + return 0; } @@ -50,12 +78,90 @@ Int_t AliFlowAnalysis::Finish() } /*********************************************************/ +Double_t AliFlowAnalysis::GetEventPlane(AliAOD* aod) +{ + //returns event plane in degrees + if (aod == 0x0) + { + Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL"); + return -1.0; + } + + Double_t psi; + Int_t mult = aod->GetNumberOfParticles(); + + Double_t ssin = 0, scos = 0; + + for (Int_t i=0; iGetParticle(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; iGetParticle(i); + if (aodtrack == 0x0) + { + Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i); + continue; + } + if (fPartCut) + if (fPartCut->Pass(aodtrack)) + continue; + + Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px()); + ssin += TMath::Sin( 2.0 * (phi - psi)); + scos += TMath::Cos( 2.0 * (phi - psi)); + } + + v2 = TMath::Hypot(ssin,scos); +} + + +/*********************************************************/ + Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd) { //returns event plane if (esd == 0x0) { - ::Error("GetFlow","Pointer to ESD is NULL"); + ::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL"); return -1.0; } @@ -69,7 +175,7 @@ Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd) AliESDtrack* esdtrack = esd->GetTrack(i); if (esdtrack == 0x0) { - ::Error("GetEventPlane","Can not get track %d", i); + ::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i); continue; } @@ -87,13 +193,14 @@ Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd) return psi; } +/*********************************************************/ void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi) { //returns flow parameters: v2 and event plane if (esd == 0x0) { - ::Error("GetFlow","Pointer to ESD is NULL"); + ::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL"); return; } @@ -107,7 +214,7 @@ void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi) AliESDtrack* esdtrack = esd->GetTrack(i); if (esdtrack == 0x0) { - ::Error("GetEventPlane","Can not get track %d", i); + ::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i); continue; } @@ -121,3 +228,4 @@ void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi) v2 = TMath::Hypot(ssin,scos); } +/*********************************************************/ diff --git a/ANALYSIS/AliFlowAnalysis.h b/ANALYSIS/AliFlowAnalysis.h index 5cf49400bdf..f20dcb3fe14 100644 --- a/ANALYSIS/AliFlowAnalysis.h +++ b/ANALYSIS/AliFlowAnalysis.h @@ -16,24 +16,30 @@ #include "AliAnalysis.h" class AliESD; +class AliAOD; class AliStack; - +class AliAODParticleCut; + class AliFlowAnalysis: public AliAnalysis { public: - AliFlowAnalysis(){} - ~AliFlowAnalysis(){} + AliFlowAnalysis(); + virtual ~AliFlowAnalysis(); Int_t Init(); - Int_t ProcessEvent(AliESD* esd, AliStack* stack = 0x0); + Int_t ProcessEvent(AliAOD* aodrec, AliAOD* aodsim = 0x0); Int_t Finish(); + void SetParticleCut(AliAODParticleCut* pcut){fPartCut = pcut;} static Double_t GetEventPlane(AliESD* esd); static void GetFlow(AliESD* esd,Double_t& v2,Double_t& psi); - protected: + Double_t GetEventPlane(AliAOD* aod); + void GetFlow(AliAOD* aod,Double_t& v2,Double_t& psi); + protected: + private: - + AliAODParticleCut* fPartCut;//Particle Cut ClassDef(AliFlowAnalysis,1) }; diff --git a/ANALYSIS/AliParticle.cxx b/ANALYSIS/AliParticle.cxx new file mode 100644 index 00000000000..eeced89f53d --- /dev/null +++ b/ANALYSIS/AliParticle.cxx @@ -0,0 +1,411 @@ +#include "AliHBTParticle.h" +//___________________________________________________________ +///////////////////////////////////////////////////////////// +// +// class AliHBTParticle +// +// Ali HBT Particle: simplified class TParticle +// Simplified in order to minimize the size of object +// - we want to keep a lot of such a objects in memory +// Additionaly adjusted for HBT Analysies purposes +// + pointer to Track Points +// + pointer to Cluster Map(s) +// +// Piotr.Skowronski@cern.ch +// +///////////////////////////////////////////////////////////// +#include +#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; iClone(); + 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> fPids[i]; +// } +// for (i = 0;i> 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 +#include +#include +#include + + +class TParticle; +class AliHBTTrackPoints; +class AliHBTClusterMap; + +class AliHBTParticle : public TObject +{ +public: + // ****** constructors and destructor + AliHBTParticle(); + AliHBTParticle(const AliHBTParticle& in); + + AliHBTParticle(Int_t pdg, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot, + Double_t vx, Double_t vy, Double_t vz, Double_t time); + + AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot, + Double_t vx, Double_t vy, Double_t vz, Double_t time); + + AliHBTParticle(const TParticle& p,Int_t idx); + + virtual ~AliHBTParticle(); + + AliHBTParticle& operator=(const AliHBTParticle& in); + + void SetPIDprobability(Int_t pdg, Float_t prob = 1.0); + Float_t GetPIDprobability(Int_t pdg) const; + + Int_t GetPdgCode () const { return (fPids)?fPids[fPdgIdx]:0;} + Int_t GetPid () const { return GetPdgCode();} + Float_t GetPidProb () const { return (fPidProb)?fPidProb[fPdgIdx]:0;} + + Int_t GetUID () const { return fIdxInEvent;} + Int_t GetNumberOfPids () const { return fNPids;} + Int_t GetNthPid (Int_t idx) const; + Float_t GetNthPidProb (Int_t idx) const; + + void SetPdgCode(Int_t pdg, Float_t prob = 1.0); + Double_t GetCalcMass () const { return fCalcMass; } + Double_t GetMass () { return (GetPDG())?GetPDG()->Mass():-1.;} + + TParticlePDG* GetPDG (){return TDatabasePDG::Instance()->GetParticle(GetPdgCode());} + + Int_t Beauty () { return GetPDG()->Beauty(); } + Int_t Charm () { return GetPDG()->Charm(); } + Int_t Strangeness () { return GetPDG()->Strangeness();} + void ProductionVertex(TLorentzVector &v) const { v.SetXYZT(fVx,fVy,fVz,fVt);} + + + Double_t Vx () const { return fVx;} + Double_t Vy () const { return fVy;} + Double_t Vz () const { return fVz;} + Double_t T () const { return fVt;} + + Double_t Px () const { return fPx; } //X coordinate of the momentum + Double_t Py () const { return fPy; } //Y coordinate of the momentum + Double_t Pz () const { return fPz; } //Z coordinate of the momentum + Double_t P () const //momentum + { return TMath::Sqrt(fPx*fPx+fPy*fPy+fPz*fPz); } + + void Momentum(TLorentzVector &v) const { v.SetPxPyPzE(fPx,fPy,fPz,fE);} + + Double_t Pt () const //transverse momentum + { return TMath::Sqrt(fPx*fPx+fPy*fPy); } + Double_t Energy() const { return fE; } + + //Pseudo Rapidity + Double_t Eta () const { if (P() != fPz) return 0.5*TMath::Log((P()+fPz)/(P()-fPz)); + else return 1.e30;} + + //Rapidity + Double_t Y () const { if (fE != fPz) return 0.5*TMath::Log((fE+fPz)/(fE-fPz)); + else return 1.e30;} + + Double_t Phi () const { return TMath::Pi()+TMath::ATan2(-fPy,-fPx); } + + Double_t Theta () const { return (fPz==0)?TMath::PiOver2():TMath::ACos(fPz/P()); } + + // setters + + void SetMomentum(Double_t px, Double_t py, Double_t pz, Double_t e) + {fPx=px; fPy=py; fPz=pz; fE=e;} + void SetMomentum(const TLorentzVector& p) + {SetMomentum(p.Px(),p.Py(),p.Pz(),p.Energy());} + + void SetProductionVertex(Double_t vx, Double_t vy, Double_t vz, Double_t t) + {fVx=vx; fVy=vy; fVz=vz; fVt=t;} + void SetProductionVertex(const TLorentzVector& v) + {SetProductionVertex(v.X(),v.Y(),v.Z(),v.T());} + void SetCalcMass(Double_t mass) {fCalcMass = mass;} + + void SetUID(Int_t id){fIdxInEvent = id;} + + const Char_t* GetName() const; + void Print() const; + + void SetTrackPoints(AliHBTTrackPoints* tpts){fTrackPoints = tpts;} + AliHBTTrackPoints* GetTrackPoints() const {return fTrackPoints;} + void SetClusterMap(AliHBTClusterMap* cm){fClusterMap = cm;} + AliHBTClusterMap* GetClusterMap() const {return fClusterMap;} + + static void SetDebug(Int_t dbg=1){fgDebug=dbg;} + static Int_t GetDebug(){return fgDebug;} + +protected: + Int_t GetPidSlot(Int_t pdg) const;//returns position of the given PID in fPids (and fPidProb) array. + +private: + Char_t fPdgIdx; // index of PDG code of the particle in fPids + Int_t fIdxInEvent; // index of a particle: the same particle can appear in the event + // many times with different pid's. Idx allows to check that they are the same particles + Int_t fNPids; // number of non-zero proboble Pids + Int_t *fPids; // [fNPids] Array with PIDs + Float_t *fPidProb; // [fNPids] PIDs probabilities + Double_t fCalcMass; // Calculated mass + + Double_t fPx; // x component of momentum + Double_t fPy; // y component of momentum + Double_t fPz; // z component of momentum + Double_t fE; // Energy + + Double_t fVx; // x of production vertex + Double_t fVy; // y of production vertex + Double_t fVz; // z of production vertex + Double_t fVt; // t of production vertex + + AliHBTTrackPoints* fTrackPoints; // track positions along trajectory - used by anti-merging cut + AliHBTClusterMap* fClusterMap; // bit map of cluters occupation; 1 if has cluter on given layer/padrow/... + + static Int_t fgDebug; //debug printout level + ClassDef(AliHBTParticle,3) // TParticle vertex particle information +}; + +#endif diff --git a/ANALYSIS/AliReader.cxx b/ANALYSIS/AliReader.cxx new file mode 100644 index 00000000000..1b5743c5576 --- /dev/null +++ b/ANALYSIS/AliReader.cxx @@ -0,0 +1,434 @@ +#include "AliReader.h" +//_________________________________________________________________________ +/////////////////////////////////////////////////////////////////////////// +// +// class AliReader +// +// Reader Base class (reads particles and tracks and +// puts it to the AliAODRun objects +// +// Provides functionality for both buffering and non-buffering reading +// This can be switched on/off via method SetEventBuffering(bool) +// The main method that inheriting classes need to implement is ReadNext() +// that read next event in queue. +// The others are: +// Bool_t ReadsRec() const; specifies if reader is able to read simulated particles +// Bool_t ReadsSim() const; specifies if reader is able to read reconstructed tracks +// void Rewind(); rewind reading to the beginning +// +// reading of next event is triggered via method Next() +// +// This class provides full functionality for reading from many sources +// User can provide TObjArray of TObjString (SetDirs method or via parameter +// in constructor) which desribes paths of directories to search data in. +// If none specified current directory is searched. +// +// Piotr.Skowronski@cern.ch +/////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include + +#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; iGetEntriesFast(); 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; iGetEntriesFast(); i++) + { + AliAODParticleCut &cut = *((AliAODParticleCut*)fCuts->At(i)); + //if some of cuts accepts all particles or some accepts particles of this type, accept + if ( (cut.GetPID() == 0) || (cut.GetPID() == pid) ) return kFALSE; + } + return kTRUE; +} +/*************************************************************************************/ + +TString& AliReader::GetDirName(Int_t entry) +{ +//returns directory name of next one to read + TString* retval;//return value + if (fDirs == 0x0) + { + retval = new TString("."); + return *retval; + } + + if ( (entry>fDirs->GetEntries()) || (entry<0))//if out of bounds return empty string + { //note that entry==0 is accepted even if array is empty (size=0) + Error("GetDirName","Name out of bounds"); + retval = new TString(); + return *retval; + } + + if (fDirs->GetEntries() == 0) + { + retval = new TString("."); + return *retval; + } + + TClass *objclass = fDirs->At(entry)->IsA(); + TClass *stringclass = TObjString::Class(); + + TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry)); + + if(dir == 0x0) + { + Error("GetDirName","Object in TObjArray is not a TObjString or its descendant"); + retval = new TString(); + return *retval; + } + if (gDebug > 0) Info("GetDirName","Returned ok %s",dir->String().Data()); + return dir->String(); +} +/*************************************************************************************/ + +void AliReader::Blend() +{ + //randomly change positions of the particles after reading + //is used to check if some distributions (of many particle properties) + //depend on the order of particles + //(tracking gives particles Pt sorted) + + if (fEventSim == 0x0) return; + + for (Int_t i = 2; i < fEventSim->GetNumberOfParticles(); i++) + { + Int_t with = gRandom->Integer(i); + fEventSim->SwapParticles(i,with); + if (fEventRec) fEventRec->SwapParticles(i,with); + } +} +/*************************************************************************************/ + +void AliReader::WriteTrackCounter() const +{ + //writes the counter histogram + + if (fTrackCounter) fTrackCounter->Write(0,TObject::kOverwrite); + else + { + Warning("WriteTrackCounter","Counter is NULL"); + } +} diff --git a/ANALYSIS/AliReader.h b/ANALYSIS/AliReader.h new file mode 100644 index 00000000000..bf891a48bf3 --- /dev/null +++ b/ANALYSIS/AliReader.h @@ -0,0 +1,101 @@ +#ifndef ALIREADER_H +#define ALIREADER_H +//_________________________________________________________________________ +/////////////////////////////////////////////////////////////////////////// +// +// class AliReader +// +// Reader Base class +// Reads particles and tracks and +// puts it to the AliAOD objects and eventuall buffers in AliAODRuns +// +// Piotr.Skowronski@cern.ch +// +/////////////////////////////////////////////////////////////////////////// + +#include +#include + +class AliAODRun; +class AliAOD; +class AliAODParticleCut; +class AliAODParticle; +class TString; +class TH1I; + +class AliReader: public TNamed +{ + public: + AliReader(); + AliReader(TObjArray*); + AliReader(const AliReader& in); + virtual ~AliReader(); + + AliReader& operator=(const AliReader& in); + + virtual Int_t Next(); + virtual void Rewind() = 0; // + + virtual Bool_t ReadsSim() const = 0; //specifies if reader is able to read simulated particles + virtual Bool_t ReadsRec() const = 0;//specifies if reader is able to read reconstructed tracks + + void AddParticleCut(AliAODParticleCut* cut); + + virtual Int_t Read(AliAODRun* particles, AliAODRun *tracks); + + virtual AliAOD* GetEventRec() const {return fEventRec;}// + virtual AliAOD* GetEventSim() const {return fEventSim;}//can not be const because position randomizer overloads it + + virtual AliAOD* GetEventRec(Int_t n); + virtual AliAOD* GetEventSim(Int_t n); + + virtual Int_t GetNumberOfRecEvents(); + virtual Int_t GetNumberOfSimEvents(); + + void SetDirs(TObjArray* dirs){fDirs = dirs;} //sets array directories names + void SetEventBuffering(Bool_t flag){fBufferEvents = flag;}//switches on/off buffering - read data are kept in local buffer + void SetBlend(Bool_t flag = kTRUE){fBlend=flag;} //set blending - randomizing particle order + virtual Int_t GetNumberOfDirs() const {return (fDirs)?fDirs->GetEntries():0;} + void ReadEventsFromTo(Int_t first,Int_t last){fFirst = first; fLast = last;} + virtual TH1I* GetTrackCounter() const {return fTrackCounter;} + virtual void WriteTrackCounter() const; + + protected: + + TObjArray* fCuts;//array with particle cuts + TObjArray* fDirs;//arry with directories to read data from + + Int_t fCurrentEvent;//! number of current event in current directory + Int_t fCurrentDir;//! number of current directory + + Int_t fNEventsRead;//!total + + AliAOD* fEventRec; //! tracks read from current event + AliAOD* fEventSim; //! particles read from current event + + AliAODRun* fRunSim; //!simulated particles + AliAODRun* fRunRec; //!reconstructed tracks + + Bool_t fIsRead;//!flag indicating if the data are already read + Bool_t fBufferEvents;//flag indicating if the data should be bufferred + + Bool_t fBlend;// flag indicating if randomly change positions of the particles after reading + + Int_t fFirst;//first event to return (all are before are skipped) + Int_t fLast;//last + + TH1I* fTrackCounter; //histogram with number of tracks read + + virtual Int_t ReadNext() = 0; //this methods reads next event and put result in fTracksEvent and/or fParticlesEvent + Bool_t Pass(AliAODParticle* p); + Bool_t Pass(Int_t pid); + void Blend(); + + TString& GetDirName(Int_t entry); + + private: + + ClassDef(AliReader,1)// +}; + +#endif diff --git a/ANALYSIS/AliReaderESD.cxx b/ANALYSIS/AliReaderESD.cxx new file mode 100644 index 00000000000..4d2e75bea8a --- /dev/null +++ b/ANALYSIS/AliReaderESD.cxx @@ -0,0 +1,711 @@ +#include "AliReaderESD.h" +//____________________________________________________________________ +////////////////////////////////////////////////////////////////////// +// // +// class AliReaderESD // +// // +// reader for ALICE Event Summary Data (ESD). // +// // +// Piotr.Skowronski@cern.ch // +// // +////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#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(esdobj); + + TString esdname = "ESD"; + esdname+=fCurrentEvent; + AliESD* esd = dynamic_cast(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;iGetTrack(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 2) + Info("ReadNext","Particle rejected since total bayessian PID probab. is zero."); + continue; + } + + for (Int_t s=0; s 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;s4) + + 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; smaxprob) + { + maxprob = w[s]; + spec = s; + } + } + firstspecie = spec; + lastspecie = spec + 1; + } + + for (Int_t s = firstspecie; s 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; kSetPIDprobability(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; sGetNumberOfParticles(), fEventSim->GetNumberOfParticles(), + fNEventsRead,fCurrentEvent,fCurrentDir); + fTrackCounter->Fill(fEventRec->GetNumberOfParticles()); + return 0; +} + +/**********************************************************/ + +void AliReaderESD::Rewind() +{ + //rewinds reading + delete fKeyIterator; + delete fFile; + fFile = 0x0; + fKeyIterator = 0x0; + delete fRunLoader; + fRunLoader = 0x0; + fCurrentDir = 0; + fNEventsRead = 0; + if (fTrackCounter) fTrackCounter->Reset(); +} +/**********************************************************/ + +TFile* AliReaderESD::OpenFile(Int_t n) +{ +//opens fFile with kine tree + + const TString& dirname = GetDirName(n); + if (dirname == "") + { + Error("OpenFiles","Can not get directory name"); + return 0x0; + } + TString filename = dirname +"/"+ fESDFileName; + TFile *ret = TFile::Open(filename.Data()); + + if ( ret == 0x0) + { + Error("OpenFiles","Can't open fFile %s",filename.Data()); + return 0x0; + } + if (!ret->IsOpen()) + { + Error("OpenFiles","Can't open fFile %s",filename.Data()); + return 0x0; + } + + if (fReadSim ) + { + fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName); + if (fRunLoader == 0x0) + { + Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data()); + delete ret; + return 0x0; + } + + fRunLoader->LoadHeader(); + if (fRunLoader->LoadKinematics()) + { + Error("Next","Error occured while loading kinematics."); + delete fRunLoader; + delete ret; + return 0x0; + } + } + + return ret; +} +/**********************************************************/ + +Int_t AliReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron +{ + //returns pdg code from the PID index + //ask jura about charge + switch (spec) + { + case kESDElectron: + return kPositron; + break; + case kESDMuon: + return kMuonPlus; + break; + case kESDPion: + return kPiPlus; + break; + case kESDKaon: + return kKPlus; + break; + case kESDProton: + return kProton; + break; + default: + ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec); + break; + } + return 0; +} +/********************************************************************/ +Bool_t AliReaderESD::CheckTrack(AliESDtrack* t) const +{ + //Performs check of the track + + if ( (t->GetConstrainedChi2() < fChi2Min) || (t->GetConstrainedChi2() > fChi2Max) ) return kTRUE; + + if ( (t->GetTPCclusters(0x0) < fNTPCClustMin) || (t->GetTPCclusters(0x0) > fNTPCClustMax) ) return kTRUE; + + if (t->GetTPCclusters(0x0) > 0) + { + Float_t chisqpercl = t->GetTPCchi2()/((Double_t)t->GetTPCclusters(0x0)); + if ( (chisqpercl < fTPCChi2PerClustMin) || (chisqpercl > fTPCChi2PerClustMax) ) return kTRUE; + } + + Double_t cc[15]; + t->GetConstrainedExternalCovariance(cc); + + if ( (cc[0] < fC00Min) || (cc[0] > fC00Max) ) return kTRUE; + if ( (cc[2] < fC11Min) || (cc[2] > fC11Max) ) return kTRUE; + if ( (cc[5] < fC22Min) || (cc[5] > fC22Max) ) return kTRUE; + if ( (cc[9] < fC33Min) || (cc[9] > fC33Max) ) return kTRUE; + if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE; + + + t->GetInnerExternalCovariance(cc); + + if ( (cc[0] < fTPCC00Min) || (cc[0] > fTPCC00Max) ) return kTRUE; + if ( (cc[2] < fTPCC11Min) || (cc[2] > fTPCC11Max) ) return kTRUE; + if ( (cc[5] < fTPCC22Min) || (cc[5] > fTPCC22Max) ) return kTRUE; + if ( (cc[9] < fTPCC33Min) || (cc[9] > fTPCC33Max) ) return kTRUE; + if ( (cc[14] < fTPCC44Min) || (cc[14] > fTPCC44Max) ) return kTRUE; + + return kFALSE; + +} +/********************************************************************/ + +void AliReaderESD::SetChi2Range(Float_t min, Float_t max) +{ + //sets range of Chi2 per Cluster + fChi2Min = min; + fChi2Max = max; +} +/********************************************************************/ + +void AliReaderESD::SetTPCNClustersRange(Int_t min,Int_t max) +{ + //sets range of Number Of Clusters that tracks have to have + fNTPCClustMin = min; + fNTPCClustMax = max; +} +/********************************************************************/ + +void AliReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max) +{ + //sets range of Chi2 per Cluster + fTPCChi2PerClustMin = min; + fTPCChi2PerClustMax = max; +} +/********************************************************************/ + +void AliReaderESD::SetC00Range(Float_t min, Float_t max) +{ + //Sets range of C00 parameter of covariance matrix of the track + //it defines uncertainty of the momentum + fC00Min = min; + fC00Max = max; +} +/********************************************************************/ + +void AliReaderESD::SetC11Range(Float_t min, Float_t max) +{ + //Sets range of C11 parameter of covariance matrix of the track + //it defines uncertainty of the momentum + fC11Min = min; + fC11Max = max; +} +/********************************************************************/ + +void AliReaderESD::SetC22Range(Float_t min, Float_t max) +{ + //Sets range of C22 parameter of covariance matrix of the track + //it defines uncertainty of the momentum + fC22Min = min; + fC22Max = max; +} +/********************************************************************/ + +void AliReaderESD::SetC33Range(Float_t min, Float_t max) +{ + //Sets range of C33 parameter of covariance matrix of the track + //it defines uncertainty of the momentum + fC33Min = min; + fC33Max = max; +} +/********************************************************************/ + +void AliReaderESD::SetC44Range(Float_t min, Float_t max) +{ + //Sets range of C44 parameter of covariance matrix of the track + //it defines uncertainty of the momentum + fC44Min = min; + fC44Max = max; +} diff --git a/ANALYSIS/AliReaderESD.h b/ANALYSIS/AliReaderESD.h new file mode 100644 index 00000000000..0b6fe2dc180 --- /dev/null +++ b/ANALYSIS/AliReaderESD.h @@ -0,0 +1,130 @@ +#ifndef AliReaderESD_H +#define AliReaderESD_H +//___________________________________________________________________________ +///////////////////////////////////////////////////////////////////////////// +// // +// Multi file reader for ESD // +// // +// This reader reads tracks from Event Summary Data // +// do not read particles // +// Piotr.Skowronski@cern.ch // +// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html // +// // +///////////////////////////////////////////////////////////////////////////// + +#include "AliReader.h" +#include +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 diff --git a/ANALYSIS/AliRunAnalysis.cxx b/ANALYSIS/AliRunAnalysis.cxx index 321b4e592f4..20d90431320 100644 --- a/ANALYSIS/AliRunAnalysis.cxx +++ b/ANALYSIS/AliRunAnalysis.cxx @@ -26,15 +26,17 @@ #include #include "AliEventCut.h" +#include "AliReader.h" +#include "AliAODParticle.h" ClassImp(AliRunAnalysis) AliRunAnalysis::AliRunAnalysis(): TTask("RunAnalysis","Alice Analysis Manager") , - fDirs(), + fReader(0x0), fEventCut(0x0), - fFileName("AliESDs.root"), - fReadKinematics(kFALSE) + fCutOnSim(kFALSE), + fCutOnRec(kTRUE) { //ctor } @@ -43,8 +45,8 @@ AliRunAnalysis::AliRunAnalysis(): AliRunAnalysis::~AliRunAnalysis() { //dtor - delete fDirs; - delete fAnalysies; + delete fEventCut; + delete fReader; delete fEventCut; } /*********************************************************/ @@ -52,249 +54,75 @@ AliRunAnalysis::~AliRunAnalysis() Int_t AliRunAnalysis::Run() { //makes analysis - - Int_t currentdir = 0; - Int_t ndirs; - if (fDirs) //if array with directories is supplied by user - { - ndirs = fDirs->GetEntries(); //get the number if directories - } - else - { - ndirs = 0; //if the array is not supplied read only from current directory - } + if (fReader == 0x0) + { + Error("Run","Reader is not set"); + return 1; + } /******************************/ /* Init Event */ /******************************/ - if (fAnalysies == 0x0) + for (Int_t an = 0; an < fAnalysies.GetEntries(); an++) { - Info("Run","No analysis present"); - return 0; - } - - for (Int_t an = 0; an < fAnalysies->GetEntries(); an++) - { - AliAnalysis* analysis = (AliAnalysis*)fAnalysies->At(an); + AliAnalysis* analysis = (AliAnalysis*)fAnalysies.At(an); analysis->Init(); } - - do + + while (fReader->Next() == kFALSE) { - TFile* file = OpenFile(currentdir); - if (file == 0x0) - { - Error("Run","Cannot get File for dir no. %d",currentdir); - currentdir++; - continue; - } - AliStack* stack = 0x0; - AliRunLoader* rl = 0x0; - if (fReadKinematics) - { - const TString& dirname = GetDirName(currentdir); - if (dirname == "") - { - Error("Run","Can not get directory name"); - return 0x0; - } - TString filename = dirname +"/galice.root"; - - rl = AliRunLoader::Open(filename); - if ( rl == 0x0 ) - { - Error("Run","Can't get Run Loader from dir %s",filename.Data()); - delete file; - currentdir++; - continue; - } - if( rl->LoadHeader() ) - { - Error("Run","Error while loading Header from dir %s",filename.Data()); - delete file; - delete rl; - currentdir++; - continue; - } - if( rl->LoadKinematics() ) - { - Error("Run","Error while loading Kinematics from dir %s",filename.Data()); - delete file; - delete rl; - currentdir++; - continue; - } - } + AliAOD* eventsim = fReader->GetEventSim(); + AliAOD* eventrec = fReader->GetEventRec(); - file->cd(); - TIter keyiter(file->GetListOfKeys()); - TKey* key; - while (( key = (TKey*)keyiter.Next() )) - { - if (key == 0x0) - { - if (GetDebug() > 2 ) - { - Info("Run","No more keys."); - } - break; - } - - TObject* esdobj = key->ReadObj(); - if (esdobj == 0x0) - { - if (GetDebug() > 2 ) - { - Info("ReadNext","Key read NULL. Key Name is %s",key->GetName()); - key->Dump(); - } - continue; - } - if (GetDebug() > 9 ) esdobj->Dump(); - AliESD* esd = dynamic_cast(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; } diff --git a/ANALYSIS/AliRunAnalysis.h b/ANALYSIS/AliRunAnalysis.h index 127084ce603..4b369a4f24e 100644 --- a/ANALYSIS/AliRunAnalysis.h +++ b/ANALYSIS/AliRunAnalysis.h @@ -12,38 +12,36 @@ // /////////////////////////////////////////////////////////// -class AliEventCut; -class TObjArray; -class TFile; - -#include #include - +#include #include "AliAnalysis.h" +class AliEventCut; +class TFile; +class AliReader; + class AliRunAnalysis: public TTask { public: AliRunAnalysis(); virtual ~AliRunAnalysis(); - Int_t Run(); - void Add(AliAnalysis* a); - void ReadKinematics(Bool_t flag){fReadKinematics = flag;} + Int_t Run(); + void Add(AliAnalysis* a); + void SetReader(AliReader* reader){fReader = reader;} + + const char* GetName(){return "RunAnalysis";} - Int_t GetDebug() {return AliAnalysis::GetDebug();} - void SetDirs(TObjArray* dirs){fDirs = dirs;} //sets array directories names; protected: - TObjArray* fAnalysies;//arry with analysies - TObjArray* fDirs;//arry with directories to read data from + TObjArray fAnalysies;//arry with analysies + AliReader* fReader;//arry with directories to read data from AliEventCut* fEventCut;//event cut - TString fFileName;//name of the file with ESDs - Bool_t fReadKinematics; + Bool_t fCutOnSim;//flag indicating that event cut is performed on simulated particles + Bool_t fCutOnRec;//flag indicating that event cut is performed on reconstructed tracks - TString& GetDirName(Int_t entry); - TFile* OpenFile(Int_t n); + Bool_t Pass(AliAOD* recevent, AliAOD* simevent); private: void SetName(const char *){}//change SetName to be private diff --git a/ANALYSIS/AliTrackPoints.cxx b/ANALYSIS/AliTrackPoints.cxx new file mode 100644 index 00000000000..3878afd3a73 --- /dev/null +++ b/ANALYSIS/AliTrackPoints.cxx @@ -0,0 +1,589 @@ +//_________________________________ +//////////////////////////////////////////////////////////// +// // +// class AliTrackPoints // +// // +// used by Anti-Merging cut // +// contains set of poits the lay on track trajectory // +// according to reconstructed track parameters - // +// NOT CLUSTERS POSITIONS!!! // +// Anti-Merging cut is applied only on tracks coming from // +// different events (that are use to fill deniminators) // +// // +//////////////////////////////////////////////////////////// + +#include +#include +#include + +#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 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; i9) + { +// 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(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;iPositionAt(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;iPositionAt(i,x,y,z); + hxy->Fill(x,y); + hxz->Fill(x,z); + printf("Rdemanded %f\n",r); + printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y)); + + //BUT they are local!!!! + t->PropagateTo(r); +// Double_t phi = t->Phi(); + Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius + + Double_t alpha = t->GetAlpha(); + Double_t salpha = TMath::Sin(alpha); + Double_t calpha = TMath::Cos(alpha); + x = t->GetX()*calpha - t->GetY()*salpha; + y = t->GetX()*salpha + t->GetY()*calpha; + z = t->GetZ(); + + printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl); + + printf("tpz - tz %f\n",z-t->GetZ()); + printf("\n"); + hxyt->Fill(x,y); + hxzt->Fill(x,z); + + } + + rl->LoadTrackRefs(); + TTree* treeTR = rl->TreeTR(); + b = treeTR->GetBranch("TPC"); + + TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100); + AliTrackReference* tref; + b->SetAddress(&trackrefs); + + Int_t tlab = TMath::Abs(t->GetLabel()); + + Int_t netr = (Int_t)treeTR->GetEntries(); + printf("Found %d entries in TR tree\n",netr); + + for (Int_t e = 0; e < netr; e++) + { + treeTR->GetEntry(e); + tref = (AliTrackReference*)trackrefs->At(0); + if (tref == 0x0) continue; + if (tref->GetTrack() != tlab) continue; + + printf("Found %d entries in TR array\n",trackrefs->GetEntries()); + + for (Int_t i = 0; i < trackrefs->GetEntries(); i++) + { + tref = (AliTrackReference*)trackrefs->At(i); + if (tref->GetTrack() != tlab) continue; + x = tref->X(); + y = tref->Y(); + z = tref->Z(); + printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z()); + + hxzTR->Fill(x,z); + hxyTR->Fill(x,y); + for (Int_t j = 1; j < 10; j++) + { + hxyTR->Fill(x, y+j*0.1); + hxyTR->Fill(x, y-j*0.1); + hxyTR->Fill(x+j*0.1,y); + hxyTR->Fill(x-j*0.1,y); + + hxzTR->Fill(x,z-j*0.1); + hxzTR->Fill(x,z+j*0.1); + hxzTR->Fill(x-j*0.1,z); + hxzTR->Fill(x+j*0.1,z); + } + } + break; + } + hxz->Draw(""); +// hxzt->Draw("same"); + hxzTR->Draw("same"); + + delete rl; +} diff --git a/ANALYSIS/AliTrackPoints.h b/ANALYSIS/AliTrackPoints.h new file mode 100644 index 00000000000..74bf907643c --- /dev/null +++ b/ANALYSIS/AliTrackPoints.h @@ -0,0 +1,49 @@ +#ifndef AliTrackPoints_H +#define AliTrackPoints_H +//_________________________________ +//////////////////////////////////////////////////////////// +// // +// class AliTrackPoints // +// // +// used by Anti-Merging cut // +// contains set of poits the lay on track trajectory // +// according to reconstructed track parameters - // +// NOT CLUSTERS POSITIONS!!! // +// Anti-Merging cut is applied only on tracks coming from // +// different events (that are use to fill deniminators) // +// // +//////////////////////////////////////////////////////////// +#include + +class AliTPCtrack; +class AliESDtrack; + +class AliTrackPoints: public TObject +{ + public: + AliTrackPoints(); + AliTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr=30, Float_t r0 = 84.1); //min TPC R = 84.1; max TPC R = 246.6cm, + AliTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr=30,Float_t r0 = 84.1); //min TPC R = 84.1; max TPC R = 246.6cm, + AliTrackPoints(const AliTrackPoints& in); + + virtual ~AliTrackPoints(); + AliTrackPoints& operator=(const AliTrackPoints& in); + + Double_t AvarageDistance(const AliTrackPoints& tr); + void PositionAt(Int_t n, Float_t &x, Float_t &y, Float_t &z); + Int_t GetDebug() const {return fgDebug;} + void SetDebug(Int_t deblevel){fgDebug = deblevel;} + static void TestTPC(Int_t entr); + static void TestESD(Int_t entr,const char* fname = "AliESDs.root"); + protected: + void MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha); + private: + Int_t fN;//number of points + Float_t* fX;//[fN]positions at x + Float_t* fY;//[fN]positions at y + Float_t* fZ;//[fN] positions at z +// Float_t* fR;//! [fN] radii + static Int_t fgDebug;//! debug level + ClassDef(AliTrackPoints,1) +}; +#endif diff --git a/ANALYSIS/analysis.C b/ANALYSIS/analysis.C new file mode 100644 index 00000000000..090b82b0c00 --- /dev/null +++ b/ANALYSIS/analysis.C @@ -0,0 +1,56 @@ +void analysis(Int_t first = -1, Int_t last = -1, const char* directory=".") +{ + + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libANALYSIS"); + + AliAODParticle::SetDebug(0); + AliRunAnalysis* analysis = new AliRunAnalysis(); + + ::Info("analysis.C","Setting dirs"); + TObjArray* dirs=0; + if ( ((first >= 0) && (last >= 0) ) && ( (last-first)>=0 ) ) + {//read from many dirs dirs + char buff[50]; + dirs = new TObjArray(last-first+1); + for (Int_t i = first; i<=last; i++) + { +// sprintf(buff,"%s/%s/%s/%05.5d",basedir,field,serie,i); + printf("%s/%d\n",directory,i); + sprintf(buff,"%s/%d",directory,i); + TObjString *odir= new TObjString(buff); + dirs->Add(odir); + } + } + + AliReaderESD* reader = new AliReaderESD(dirs); + reader->ReadSimulatedData(kTRUE); + reader->SetReadMostProbableOnly(kTRUE); + +/* + //example PID cuts + AliAODParticleCut* partcut = new AliAODParticleCut(); + partcut->SetPID(kPiPlus);//here we define the incarnation + AliAODPIDCut* pidcut = new AliAODPIDCut(kPiPlus,0.5);//accept all particles types that have PID prob > 50% + partcut->AddBasePartCut(pidcut);// + reader->AddParticleCut(partcut);//This guy makes a copy of a cut for himself so we can modify it here + + partcut->SetPID(kPiMinus);//here we define that particle has incarnation PiMinus + pidcut->SetPID(kPiMinus);//here we define to check if PID probability of being kPiMinus is greater thann 0.5 (number defined few lines above) + reader->AddParticleCut(partcut); + + pidcut->SetPID(kKPlus); + pidcut->SetPID(kKPlus); + reader->AddParticleCut(partcut); + + pidcut->SetPID(kKMinus); + pidcut->SetPID(kKMinus); + reader->AddParticleCut(partcut); +*/ + + AliFlowAnalysis* flow = new AliFlowAnalysis(); + analysis->SetReader(reader); + + analysis->Add(flow); + analysis->Run(); + delete analysis; +} diff --git a/ANALYSIS/libANALYSIS.pkg b/ANALYSIS/libANALYSIS.pkg index ac70c18b93a..a336363d355 100644 --- a/ANALYSIS/libANALYSIS.pkg +++ b/ANALYSIS/libANALYSIS.pkg @@ -1,7 +1,13 @@ -SRCS= AliRunAnalysis.cxx AliAnalysis.cxx AliEventCut.cxx \ - AliD0toKpi.cxx AliD0toKpiAnalysis.cxx AliFlowAnalysis.cxx - +SRCS= AliAOD.cxx \ + AliAODParticle.cxx AliAODStdParticle.cxx \ + AliAODParticleCut.cxx AliAODRun.cxx \ + AliRunAnalysis.cxx AliAnalysis.cxx AliEventCut.cxx \ + AliReader.cxx AliReaderESD.cxx\ + AliTrackPoints.cxx AliClusterMap.cxx \ + AliD0toKpi.cxx AliD0toKpiAnalysis.cxx AliFlowAnalysis.cxx \ HDRS= $(SRCS:.cxx=.h) DHDR:=ANALYSISLinkDef.h + +EINCLUDE:= TPC CONTAINERS ITS -- 2.43.0