ALIROOTMODULES += PWG0
endif
+ifeq ($(findstring PWG2,$(MAKECMDGOALS)),PWG2)
+ALIROOTMODULES += PWG2
+endif
+
ifeq ($(findstring PWG3,$(MAKECMDGOALS)),PWG3)
ALIROOTMODULES += PWG3
endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnAnalysis
+// Reconstruction and analysis of a binary resonance
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#include <Riostream.h>
+
+#include <TH1.h>
+#include <TTree.h>
+#include <TRefArray.h>
+#include <TParticle.h>
+#include <TObjectTable.h>
+#include <TDatabasePDG.h>
+
+#include "AliRsnDaughter.h"
+#include "AliRsnDaughterCut.h"
+#include "AliRsnEvent.h"
+#include "AliRsnAnalysis.h"
+
+ClassImp(AliRsnAnalysis)
+
+//--------------------------------------------------------------------------------------------------------
+AliRsnAnalysis::AliRsnAnalysis()
+//
+// Constructor
+// Initializes all pointers and collections to NULL.
+//
+{
+ fMixHistograms = 0;
+ fHistograms = 0;
+ fEventsTree = 0;
+ fMixPairDefs = 0;
+ fPairDefs = 0;
+ fPairCuts = 0;
+
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i] = 0;
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnAnalysis::AddCutPair(AliRsnDaughterCut *cut)
+//
+// Add a cut on pairs.
+// This cut is global for all pairs.
+//
+{
+ if (!cut->IsPairCut()) {
+ Warning("AddCutPair", "This is a single cut, cannot be added");
+ return;
+ }
+
+ if (!fPairCuts) fPairCuts = new TObjArray(0);
+
+ fPairCuts->AddLast(cut);
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnAnalysis::AddCutSingle(AliPID::EParticleType type, AliRsnDaughterCut *cut)
+//
+// Add a cut on single particles.
+// This cut must be specified for each particle type.
+//
+{
+ if (cut->IsPairCut()) {
+ Warning("AddCutSingle", "This is a pair cut, cannot be added");
+ return;
+ }
+
+ if (type >= AliPID::kElectron && type <= AliPID::kProton) {
+ if (!fCuts[type]) fCuts[type] = new TObjArray(0);
+ fCuts[type]->AddLast(cut);
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnAnalysis::AddMixPairDef
+(AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2)
+//
+// Adds a new pair definition to create a new histogram,
+// for the event mixing step.
+// If the pair defs array is NULL, it is initialized here.
+//
+{
+ if (!fMixPairDefs) fMixPairDefs = new TObjArray(0);
+ fMixPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, 0, kFALSE) );
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnAnalysis::AddPairDef
+(AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2, Bool_t onlyTrue)
+//
+// Adds a new pair definition to create a new histogram,
+// for the signal evaluation (same event) step.
+// If the pair defs array is NULL, it is initialized here.
+//
+{
+ if (!fPairDefs) fPairDefs = new TObjArray(0);
+ fPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, fTrueMotherPDG, onlyTrue) );
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnAnalysis::Clear(Option_t* /* option */)
+//
+// Clear heap
+//
+{
+ fHistograms->Clear("C");
+ fPairDefs->Clear("C");
+
+ fMixHistograms->Clear("C");
+ fMixPairDefs->Clear("C");
+
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i]->Clear("C");
+ fPairCuts->Clear("C");
+}
+//--------------------------------------------------------------------------------------------------------
+Stat_t AliRsnAnalysis::EventMix(Int_t nmix, Int_t multDiffMax, Double_t vzDiffMax, Bool_t compareTotMult)
+//
+// Performs event mixing.
+// It takes the array of fMixPairDefs and stores results in fMixHistograms.
+// First parameter defines how many events must be mixed with each one.
+// Events to be mixed are chosen using the other parameters:
+//
+// - multDiffMax defines the maximum allowed difference in particle multiplicity,
+// a) when last argument is kFALSE (default) the multiplicity comparison is made with respect
+// to the particles of 'type 2' in the pair
+// b) when last argument is kTRUE, the multiplicity comparison is made with respect of total
+// particle multiplicity in the events
+//
+// - vzDiffMax defines maximum allowed difference in Z coordinate of prim. vertex.
+//
+// If one wants to exchange the particle types, one has to add both combinations of particles
+// as different pair defs.
+//
+// EXAMPLE:
+// analysis->AddMixPairDef(AliRsnEvent::kPion, '+', AliRsnEvent::kKaon, '-');
+// analysis->AddMixPairDef(AliRsnEvent::kKaon, '-', AliRsnEvent::kPion, '+');
+//
+{
+ // allocate the histograms array
+ Int_t i, npairdefs = (Int_t)fMixPairDefs->GetEntries();
+ fMixHistograms = new TObjArray(npairdefs);
+ TObjArray &histos = *fMixHistograms;
+ AliPairDef *pd = 0;
+ for (i = 0; i < npairdefs; i++) {
+ pd = (AliPairDef*)fMixPairDefs->At(i);
+ histos[i] = new TH1D(Form("hmix_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax);
+ }
+
+ // Link events branch
+ AliRsnEvent *event2 = new AliRsnEvent;
+ fEventsTree->SetBranchAddress("events", &event2);
+
+ // define variables to store data about particles
+ Stat_t nPairs = 0;
+ Int_t iev, ievmix, nEvents = (Int_t)fEventsTree->GetEntries();
+
+ // loop on events
+ Int_t nmixed, mult1, mult2, diffMult;
+ Double_t vz1, vz2, diffVz;
+ for (iev = 0; iev < nEvents; iev++) {
+
+ // get event
+ event2->Clear("DELETE");
+ fEventsTree->GetEntry(iev);
+
+ // copy this event into 'event 1'
+ AliRsnEvent *event1 = new AliRsnEvent(*event2);
+
+ // get Z position of primary vertex
+ vz1 = event1->GetPrimaryVertexZ();
+
+ // if total multiplicities must be used
+ // it is computed here
+ if (compareTotMult) {
+ mult1 = event1->GetMultiplicity(kTRUE);
+ }
+ else {
+ mult1 = 0;
+ }
+
+ // message
+ if (iev % 10 == 0) cout << "\rMixing event " << iev << flush;
+
+ // loop on pair definitions
+ for (i = 0; i < npairdefs; i++) {
+
+ // get pair definition
+ pd = (AliPairDef*)fMixPairDefs->At(i);
+
+ // get histogram reference
+ TH1D *histogram = (TH1D*)histos[i];
+
+ // get multiplicity of particle type '2' in event '1'
+ if (!mult1) {
+ mult1 = (Int_t)event1->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries();
+ }
+
+ // loop on events for mixing
+ nmixed = 0;
+ ievmix = iev;
+ while (nmixed < nmix) {
+
+ // get other event (it starts from iev + 1, and
+ // loops again to 0 if reachs end of tree)
+ ievmix++;
+ if (ievmix >= nEvents) ievmix -= nEvents;
+
+ // if event-2-mix is equal to 'iev', then
+ // a complete loop has been done, and no more
+ // mixing can be done, even if they are too few
+ // then the loop breaks anyway
+ if (ievmix == iev) break;
+
+ // get other event
+ event2->Clear("DELETE");
+ fEventsTree->GetEntry(ievmix);
+
+ // get event stats related to event 2
+ vz2 = event2->GetPrimaryVertexZ();
+ if (compareTotMult) {
+ mult2 = event2->GetMultiplicity(kTRUE);
+ }
+ else {
+ mult2 = event2->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries();
+ }
+
+ // fill histogram if checks are passed
+ diffMult = TMath::Abs(mult1 - mult2);
+ diffVz = TMath::Abs(vz1 - vz2);
+ if (diffVz <= vzDiffMax && diffMult <= multDiffMax) {
+ //cout << ievmix << " " << flush;
+ nPairs += Compute(pd, histogram, event1, event2);
+ nmixed++;
+ }
+ }
+ }
+
+ event1->Clear("DELETE");
+ delete event1;
+ event1 = 0;
+
+ } // end loop on events
+ cout << endl;
+
+ return nPairs;
+}
+//--------------------------------------------------------------------------------------------------------
+Stat_t AliRsnAnalysis::Process()
+//
+// Reads the list 'fPairDefs', and builds an inv-mass histogram for each definition.
+// For each event, particle of 'type 1' are combined with particles of 'type 2' as
+// defined in each pair definition specified in the list, taking the list of both types
+// from the same event.
+// This method is supposed to evaluate the 'signal' of the resonance, containing the peak.
+// It can also be used when one wants to evaluate the 'background' with the 'wrong sign'
+// particle pairs.
+//
+{
+ // allocate the histograms array in order to contain
+ // as many objects as the number of pair definitionss
+ Int_t i, npairdefs = (Int_t)fPairDefs->GetEntries();
+ fHistograms = new TObjArray(npairdefs);
+ TObjArray &histos = *fHistograms;
+
+ // loop over pair defs in order to retrieve the particle species chosen
+ // which are used to set an unique name for the output histogram
+ // There will be a direct correspondence between the inder of pairdef in its array
+ // and the corresponding histogram in the 'fHistograms' list.
+ AliPairDef *pd = 0;
+ for (i = 0; i < npairdefs; i++) {
+ pd = (AliPairDef*)fPairDefs->At(i);
+ histos[i] = new TH1D(Form("h_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax);
+ }
+
+ // Link events branch
+ AliRsnEvent *event = new AliRsnEvent;
+ fEventsTree->SetBranchAddress("events", &event);
+
+ // define variables to store data about particles
+ Stat_t nPairs = 0;
+ Int_t iev, nEvents = (Int_t)fEventsTree->GetEntries();
+
+ // loop on events
+ for (iev = 0; iev < nEvents; iev++) {
+
+ // get event
+ event->Clear("DELETE");
+ fEventsTree->GetEntry(iev);
+ if (iev % 1000 == 0) cout << "\rProcessing event " << iev << flush;
+
+ // loop over the collection of pair defs
+ for (i = 0; i < npairdefs; i++) {
+ pd = (AliPairDef*)fPairDefs->At(i);
+ TH1D *histogram = (TH1D*)histos[i];
+ // invoke AliPairDef method to fill histogram
+ nPairs += Compute(pd, histogram, event, event);
+ }
+
+ } // end loop on events
+ cout << endl;
+
+ return nPairs;
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnAnalysis::WriteHistograms()
+//
+// Writes histograms in current directory
+//
+{
+ TH1D *histogram;
+ TObjArrayIter iter(fHistograms);
+ while ( (histogram = (TH1D*)iter.Next()) ) {
+ histogram->Write();
+ }
+
+ if (fMixHistograms) {
+ TObjArrayIter iterMix(fMixHistograms);
+ while ( (histogram = (TH1D*)iterMix.Next()) ) {
+ histogram->Write();
+ }
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+Stat_t AliRsnAnalysis::Compute
+(AliPairDef *pd, TH1D* &histogram, AliRsnEvent *event1, AliRsnEvent *event2)
+//
+// Adds to the specified histogram the invariant mass spectrum calculated taking
+// particles of type 1 from event 1, and particles of type 2 from event 2.
+// Events can be equal (signal) or different (background with event mixing).
+//
+{
+ // define two 'cursor' objects
+ AliRsnDaughter *track1 = 0, *track2 = 0;
+
+ // define iterators for the two collections
+ if (!event1 || !event1->GetTracks(pd->GetSign1(), pd->GetParticle1())) {
+ Error("Compute", "Null pointer for particle 1 collection");
+ return 0.0;
+ }
+ if (!event2 || !event2->GetTracks(pd->GetSign2(), pd->GetParticle2())) {
+ Error("Compute", "Null pointer for particle 2 collection");
+ return 0.0;
+ }
+ TObjArrayIter iter1(event1->GetTracks(pd->GetSign1(), pd->GetParticle1()));
+ TObjArrayIter iter2(event2->GetTracks(pd->GetSign2(), pd->GetParticle2()));
+
+ // define temporary variables for better code readability
+ Stat_t nPairs = 0;
+
+ // loop on particle of type 1 (in event 1)
+ while ( (track1 = (AliRsnDaughter*)iter1.Next()) ) {
+
+ if (fRejectFakes && (track1->GetLabel() < 0)) continue;
+ if (!SingleCutCheck(pd->GetParticle1(), track1)) continue;
+
+ iter2.Reset();
+
+ // loop on particles of type 2 (in event 2)
+ while ( (track2 = (AliRsnDaughter*)iter2.Next()) ) {
+
+ if (fRejectFakes && (track2->GetLabel() < 0)) continue;
+ if (event1 == event2 && track1->GetIndex() == track2->GetIndex()) continue;
+ if (!SingleCutCheck(pd->GetParticle2(), track2)) continue;
+ if (!PairCutCheck(track1, track2)) continue;
+
+ /*
+ // check
+ Char_t sign1 = (track1->GetSign() > 0) ? '+' : '-';
+ Char_t sign2 = (track2->GetSign() > 0) ? '+' : '-';
+ Info("Compute", "Particle 1: PDG = %d, sign = %c --- Particle 2: PDG = %d, sign = %c", track1->GetTruePDG(), sign1, track2->GetTruePDG(), sign2);
+ */
+
+ // total 4-momentum
+ track1->SetMass(pd->GetMass1());
+ track2->SetMass(pd->GetMass2());
+ AliRsnDaughter sum = (*track1) + (*track2);
+
+ // if the choice to get ONLY true pairs is selected, a check is made on the mothers
+ if (pd->GetOnlyTrue()) {
+ // the 'sum' AliRsnDaughter object is initialized with
+ // the PDG code of the common mother (if it is the case)
+ if (TMath::Abs(sum.GetMotherPDG()) == TMath::Abs(fTrueMotherPDG)) {
+ histogram->Fill(sum.GetMass());
+ nPairs++;
+ }
+ }
+ else {
+ histogram->Fill(sum.GetMass());
+ nPairs++;
+ }
+
+ } // end loop 2
+
+ } // end loop 1
+
+ return nPairs;
+}
+//--------------------------------------------------------------------------------------------------------
+Bool_t AliRsnAnalysis::SingleCutCheck(Int_t itype, AliRsnDaughter *track)
+//
+// Checks a track against single particle cuts (if defined)
+//
+{
+ if (!fCuts[itype]) return kTRUE;
+
+ TObjArrayIter iter(fCuts[itype]);
+ AliRsnDaughterCut *cut = 0;
+ while ( (cut = (AliRsnDaughterCut*)iter.Next()) ) {
+ if (!cut->Pass(track)) return kFALSE;
+ }
+
+ return kTRUE;
+}
+//--------------------------------------------------------------------------------------------------------
+Bool_t AliRsnAnalysis::PairCutCheck(AliRsnDaughter *track1, AliRsnDaughter *track2)
+//
+// Checks a pair against pair cuts (if defined)
+//
+{
+ if (!fPairCuts) return kTRUE;
+
+ TObjArrayIter iter(fPairCuts);
+ AliRsnDaughterCut *cut = 0;
+ while ( (cut = (AliRsnDaughterCut*)iter.Next()) ) {
+ if (!cut->Pass(track1, track2)) return kFALSE;
+ }
+
+ return kTRUE;
+}
+//--------------------------------------------------------------------------------------------------------
+AliRsnAnalysis::AliPairDef::AliPairDef
+(AliPID::EParticleType p1, Char_t sign1,
+ AliPID::EParticleType p2, Char_t sign2, Int_t pdgMother, Bool_t onlyTrue)
+//
+// Constructor for nested class
+//
+{
+ fOnlyTrue = onlyTrue;
+ fTrueMotherPDG = 0;
+ if (fOnlyTrue) fTrueMotherPDG = pdgMother;
+
+ fSign1 = sign1;
+ fParticle1 = p1;
+
+ fSign2 = sign2;
+ fParticle2 = p2;
+
+ // instance a PDG database to recovery true masses of particles
+ TDatabasePDG *db = TDatabasePDG::Instance();
+ Int_t pdg1 = AliPID::ParticleCode((Int_t)p1);
+ Int_t pdg2 = AliPID::ParticleCode((Int_t)p2);
+ fMass1 = db->GetParticle(pdg1)->Mass();
+ fMass2 = db->GetParticle(pdg2)->Mass();
+
+ // set name according to the chosen particles
+ fName.Append(Form("%s(%c)_%s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2));
+ fTitle.Append(Form("Inv. mass for %s(%c) and %s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2));
+ if (onlyTrue) {
+ fName.Append("_true");
+ fTitle.Append(" (true pairs)");
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+Text_t* AliRsnAnalysis::AliPairDef::ParticleName(AliPID::EParticleType part)
+//
+// [PRIVATE]
+// Returns the name of the particle in text format
+//
+{
+ if (part == AliPID::kElectron) return ("E");
+ else if (part == AliPID::kMuon) return ("Mu");
+ else if (part == AliPID::kPion) return ("Pi");
+ else if (part == AliPID::kKaon) return ("K");
+ else if (part == AliPID::kProton) return ("P");
+ else {
+ Warning("ParticleName", "Unrecognized value of EParticle argument");
+ return ("???");
+ }
+}
+//--------------------------------------------------------------------------------------------------------
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnAnalysis
+// Reconstruction and analysis of K* Rsn
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#ifndef AliRsnANALYSIS_H
+#define AliRsnANALYSIS_H
+
+#include <TNamed.h>
+
+#include "AliPID.h"
+
+class TH1D;
+class TObjArray;
+class AliRsnEvent;
+class AliRsnDaughter;
+class AliRsnDaughterCut;
+
+class AliRsnAnalysis : public TObject
+{
+public:
+
+ AliRsnAnalysis();
+ virtual ~AliRsnAnalysis() {Clear();}
+
+ void AddCutPair(AliRsnDaughterCut *cut);
+ void AddCutSingle(AliPID::EParticleType type, AliRsnDaughterCut *cut);
+ void AddMixPairDef(AliPID::EParticleType p1, Char_t s1, AliPID::EParticleType p2, Char_t s2);
+ void AddPairDef(AliPID::EParticleType p1, Char_t s1, AliPID::EParticleType p2, Char_t s2, Bool_t onlyTrue = kFALSE);
+ void Clear(Option_t *option = "");
+ Stat_t EventMix(Int_t nmix = 5, Int_t multDiffMax = 10, Double_t vzDiffMax = 0.01, Bool_t compareTotalMult = kFALSE);
+ Stat_t Process();
+ void SetBins(Int_t nbins, Double_t min, Double_t max) {fNBins=nbins;fHistoMin=min;fHistoMax=max;}
+ void SetEventsTree(TTree *tree) {fEventsTree = tree;}
+ void SetRejectFakes(Bool_t doit=kTRUE) {fRejectFakes = doit;}
+ void SetTrueMotherPDG(Int_t pdg) {fTrueMotherPDG = pdg;}
+ void WriteHistograms();
+
+private:
+
+ class AliPairDef : public TNamed
+ {
+ public:
+
+ AliPairDef(AliPID::EParticleType p1, Char_t s1, AliPID::EParticleType p2, Char_t s2, Int_t pdg, Bool_t onlyTrue = kFALSE);
+
+ virtual ~AliPairDef() { }
+
+ Char_t GetSign1() {return fSign1;}
+ Char_t GetSign2() {return fSign2;}
+ Bool_t GetOnlyTrue() {return fOnlyTrue;}
+ AliPID::EParticleType GetParticle1() {return fParticle1;}
+ AliPID::EParticleType GetParticle2() {return fParticle2;}
+ Double_t GetMass1() {return fMass1;}
+ Double_t GetMass2() {return fMass2;}
+
+ void SetSign1(Char_t value) {fSign1 = value;}
+ void SetSign2(Char_t value) {fSign2 = value;}
+ void SetOnlyTrue(Bool_t value = kTRUE) {fOnlyTrue = value;}
+ void SetTrueMotherPDG(Int_t pdg) {fTrueMotherPDG = pdg;}
+ void SetParticle1(AliPID::EParticleType p) {fParticle1 = p;}
+ void SetParticle2(AliPID::EParticleType p) {fParticle2 = p;}
+
+ Text_t* ParticleName(AliPID::EParticleType part);
+
+ private:
+
+ Bool_t fOnlyTrue; // flag to be used for spectra of true pairs
+ Int_t fTrueMotherPDG; // PDG code of true mother (if requested)
+
+ Double_t fMass1; //
+ Char_t fSign1; // info about particle 1
+ AliPID::EParticleType fParticle1; //
+
+ Double_t fMass2; //
+ Char_t fSign2; // info about particle 2
+ AliPID::EParticleType fParticle2; //
+ };
+
+ Stat_t Compute(AliPairDef *pd, TH1D* &h, AliRsnEvent *ev1, AliRsnEvent *ev2);
+ Bool_t SingleCutCheck(Int_t itype, AliRsnDaughter *track);
+ Bool_t PairCutCheck(AliRsnDaughter *track1, AliRsnDaughter *track2);
+
+ Bool_t fRejectFakes; // reject particles labeled as fake
+
+ Int_t fNBins; // number of histogram bins
+ Double_t fHistoMin; // minimum of the histograms
+ Double_t fHistoMax; // maximum of the histograms
+
+ Int_t fTrueMotherPDG; // PDG code of true mother (used to create 'true' histos)
+
+ TObjArray *fMixPairDefs; // list of pair definitions for histograms (event mixing)
+ TObjArray *fMixHistograms; //! list of invmass histograms created (event mixing)
+
+ TObjArray *fPairDefs; // list of pair definitions for histograms
+ TObjArray *fHistograms; //! list of invmass histograms created
+
+ TObjArray *fCuts[AliPID::kSPECIES]; //! list of single particle cuts for each particle type
+ TObjArray *fPairCuts; //! list of pair cuts
+ TTree *fEventsTree; //! TTree of events (can not be created here, must be passed)
+
+ // Rsn analysis implementation
+ ClassDef(AliRsnAnalysis,1)
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnDaughter
+//
+// A simple object which describes a reconstructed track
+// with some references to its related generated particle
+// and some facilities which could help in composing its
+// 4-momentum, for resonance study.
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#include <Riostream.h>
+
+#include <TParticle.h>
+#include <TParticlePDG.h>
+
+#include "AliESDtrack.h"
+#include "AliRsnDaughter.h"
+
+ClassImp(AliRsnDaughter)
+
+//--------------------------------------------------------------------------------------------------------
+AliRsnDaughter::AliRsnDaughter()
+//
+// Default constructor.
+// Its unique argument defines how many PID weights are allowed.
+// Actually, it should be the same as AliESDtrack::kSPECIES (=5).
+//
+{
+ fSign = (Char_t)0;
+ fPDG = (UShort_t)0;
+ fIndex = (UShort_t)0;
+
+ fP[0] = fP[1] = fP[2] = 0.0;
+ fV[0] = fV[1] = fV[2] = 0.0;
+ fMass = 0.0;
+
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ fPIDwgt[i] = 0.0;
+ }
+
+ fLabel = -1;
+ fTruePDG = 0;
+ fMother = -1;
+ fMotherPDG = 0;
+}
+//--------------------------------------------------------------------------------------------------------
+AliRsnDaughter::AliRsnDaughter(const AliRsnDaughter ©) : TObject(copy)
+//
+// Copy constructor
+//
+{
+ fSign = copy.fSign;
+ fPDG = copy.fPDG;
+ fIndex = copy.fIndex;
+
+ Int_t i;
+ for (i = 0; i < 3; i++) {
+ fP[i] = copy.fP[i];
+ fV[i] = copy.fV[i];
+ }
+ fMass = copy.fMass;
+
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ fPIDwgt[i] = copy.fPIDwgt[i];
+ }
+
+ fLabel = copy.fLabel;
+ fTruePDG = copy.fTruePDG;
+ fMother = copy.fMother;
+ fMotherPDG = copy.fMotherPDG;
+}
+//--------------------------------------------------------------------------------------------------------
+Bool_t AliRsnDaughter::Adopt(const AliESDtrack* esdTrack, Bool_t checkITSRefit)
+//
+// Copies reconstructed data from an AliESDtrack:
+//
+// - charge sign
+// - momentum
+// - point of closest approach to primary vertex
+// - ESD pid weights
+// - track label (AliESDtrack::GetLabel())
+//
+// Makes the following checks:
+//
+// - if argument 'checkITSRefit' is TRUE and the "ITS refit" flag is FALSE
+// in the track, the "fIsOK" flag of (this) is set to "false" (track should be rejected)
+// - if the passed label is negative, track is considered as "fake", and
+// this info is kept to allow fake tracks exclusion when doing analysis
+//
+{
+ // check for refit in the ITS (if requested)
+ if (checkITSRefit) {
+ if ( !(esdTrack->GetStatus() & AliESDtrack::kITSrefit) ) {
+ return kFALSE;
+ }
+ }
+
+ // get sign and number of species allowed for PID
+ fSign = (Char_t)esdTrack->GetSign();
+
+ // get (and check) momentum
+ esdTrack->GetPxPyPz(fP);
+ if (fP[0] == 0.0 || fP[1] == 0.0 || fP[2] == 0.0) {
+ return kFALSE;
+ }
+
+ // get (and check) vertex
+ esdTrack->GetXYZ(fV);
+ if (fV[0] == 0.0 || fV[1] == 0.0 || fV[2] == 0.0) {
+ return kFALSE;
+ }
+
+ // get label
+ // (other kinematics informations are set to default and meaningless values)
+ fLabel = esdTrack->GetLabel();
+
+ // get PID weights
+ esdTrack->GetESDpid(fPIDwgt);
+
+ return kTRUE;
+}
+//--------------------------------------------------------------------------------------------------------
+Bool_t AliRsnDaughter::Adopt(TParticle* particle)
+//
+// Copies data from a generated particle:
+//
+// - PDG code
+// - charge sign
+// - momentum
+// - production vertex
+// - GEANT label of mother track
+//
+// When an AliRsnDaughter is copied from a TParticle, it is
+// considered always good for analysis and never fake.
+//
+{
+ // get particle sign form the sign of PDG code
+ Int_t pdg = particle->GetPdgCode();
+ if (TMath::Abs(pdg) < 20) {
+ if (pdg > 0) fSign = -1; else fSign = 1;
+ }
+ else if (TMath::Abs(pdg) < 3000) {
+ if (pdg > 0) fSign = 1; else fSign = -1;
+ }
+
+ // get momentum
+ fP[0] = particle->Px();
+ fP[1] = particle->Py();
+ fP[2] = particle->Pz();
+
+ // get vertex
+ fV[0] = particle->Vx();
+ fV[1] = particle->Vy();
+ fV[2] = particle->Vz();
+
+ // set simulation data
+ fPDG = fTruePDG = (Short_t)particle->GetPdgCode();
+ fMother = particle->GetFirstMother();
+ fMotherPDG = 0;
+
+ return kTRUE;
+}
+//--------------------------------------------------------------------------------------------------------
+AliRsnDaughter operator+(AliRsnDaughter t1, AliRsnDaughter t2)
+//
+// Sum operator overloading.
+// Builds a new AliRsnDaughter object with the sum of momenta of two particles.
+//
+{
+ // create new AliRsnDaughter with default useless values for datamembers
+ AliRsnDaughter out;
+
+ // if the summed particles are daughters of the same resonance
+ // their common mother label becomes the label of the sum
+ Int_t mum1 = t1.GetMother();
+ Int_t mum2 = t2.GetMother();
+ if (mum1 == mum2) {
+ out.SetLabel(mum1);
+ out.SetMotherPDG(t1.GetMotherPDG());
+ }
+ else {
+ out.SetLabel(-1);
+ out.SetMotherPDG(0);
+ }
+
+ // compute total 4-momentum
+ Double_t Etot = t1.GetEnergy() + t2.GetEnergy();
+ Double_t pxTot = t1.GetPx() + t2.GetPx();
+ Double_t pyTot = t1.GetPy() + t2.GetPy();
+ Double_t pzTot = t1.GetPz() + t2.GetPz();
+ Double_t mass = TMath::Sqrt(Etot*Etot - pxTot*pxTot - pyTot*pyTot - pzTot*pzTot);
+
+ //TLorentzVector v1 = track1.Get4Momentum();
+ //TLorentzVector v2 = track2.Get4Momentum();
+ //TLorentzVector sm = v1 + v2;
+ //Double_t pxTot = sum.X();
+ //Double_t pyTot = sum.Y();
+ //Double_t pzTot = sum.Z();
+ //Double_t mass = sm.M();
+
+ out.SetPxPyPz(pxTot, pyTot, pzTot);
+ out.SetMass(mass);
+
+ return out;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnDaughter
+//
+// A simple object which describes a reconstructed track
+// with some references to its related generated particle
+// and some facilities which could help in composing its
+// 4-momentum, for resonance study.
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#ifndef ALIRSNDAUGHTER_H
+#define ALIRSNDAUGHTER_H
+
+#include "AliPID.h"
+#include <TVector3.h>
+#include <TParticle.h>
+#include <TLorentzVector.h>
+
+class AliESDtrack;
+
+class AliRsnDaughter : public TObject
+{
+public:
+ AliRsnDaughter();
+ AliRsnDaughter(const AliRsnDaughter ©);
+
+ virtual ~AliRsnDaughter() { }
+
+ Bool_t Adopt(TParticle* particle);
+ Bool_t Adopt(const AliESDtrack* track, Bool_t checkRefit = kTRUE);
+ TVector3 Get3Momentum() const {TVector3 v(fP[0],fP[1],fP[2]); return v;}
+ TLorentzVector Get4Momentum() const {TLorentzVector v(fP[0],fP[1],fP[2],GetEnergy()); return v;}
+ Double_t GetEnergy() const {return TMath::Sqrt(fMass*fMass + GetP2());}
+ UShort_t GetIndex() const {return fIndex;}
+ Int_t GetLabel() const {return fLabel;}
+ Double_t GetMass() const {return fMass;}
+ Int_t GetMother() const {return fMother;}
+ Short_t GetMotherPDG() const {return fMotherPDG;}
+ UShort_t GetPDG() const {return fPDG;}
+ Double_t GetPIDweight(Int_t i) const {return ( (i>=0&&i<AliPID::kSPECIES)?fPIDwgt[i]:-1.0 );}
+ Char_t GetSign() const {return fSign;}
+ Double_t GetP2() const {return fP[0]*fP[0] + fP[1]*fP[1] + fP[2]*fP[2];}
+ Double_t GetP() const {return TMath::Sqrt(GetP2());}
+ Double_t GetPx() const {return fP[0];}
+ Double_t GetPy() const {return fP[1];}
+ Double_t GetPz() const {return fP[2];}
+ Double_t GetPt() const {return TMath::Sqrt(fP[0]*fP[0] + fP[1]*fP[1]);}
+ Short_t GetTruePDG() const {return fTruePDG;}
+ TVector3 GetVertex() const {TVector3 v(fV[0],fV[1],fV[2]); return v;}
+ Double_t GetVx() const {return fV[0];}
+ Double_t GetVy() const {return fV[1];}
+ Double_t GetVz() const {return fV[2];}
+ Double_t GetVt() const {return TMath::Sqrt(fV[0]*fV[0] + fV[1]*fV[1]);}
+ void SetIndex(UShort_t value) {fIndex = value;}
+ void SetIndex(Int_t value) {fIndex = (UShort_t)value;}
+ void SetLabel(Int_t l) {fLabel = l;}
+ void SetMass(Double_t m) {fMass = m;}
+ void SetMother(Int_t l) {fMother = l;}
+ void SetMotherPDG(Short_t pdg) {fMotherPDG = pdg;}
+ void SetPDG(UShort_t pdg) {fPDG = TMath::Abs(pdg);}
+ void SetPDG(Int_t pdg) {fPDG = (UShort_t)TMath::Abs(pdg);}
+ void SetPIDweights(const Double_t *pid) {Int_t i;for(i=0;i<AliPID::kSPECIES;i++)fPIDwgt[i]=pid[i];}
+ void SetPxPyPz(Double_t px,Double_t py,Double_t pz) {fP[0]=px;fP[1]=py;fP[2]=pz;}
+ void SetSign(Char_t value) {fSign = value;}
+ void SetSign(Int_t value) {fSign = (Char_t)value;}
+ void SetTruePDG(Short_t pdg) {fTruePDG = pdg;}
+ void SetVxVyVz(Double_t vx,Double_t vy,Double_t vz) {fV[0]=vx;fV[1]=vy;fV[2]=vz;}
+
+ friend AliRsnDaughter operator+(AliRsnDaughter t1, AliRsnDaughter t2);
+
+private:
+
+ Char_t fSign; // charge sign
+ UShort_t fPDG; // assigned PDG code from PID (0=undefined)
+ UShort_t fIndex; // reference index in AliESD container
+
+ Double_t fP[3]; // vector momentum
+ Double_t fV[3]; // production vertex
+ Double_t fMass; // mass
+
+ Double_t fPIDwgt[AliPID::kSPECIES]; // particle PID weights
+
+ // The following data are useful for simulated events only
+ // and can be left blank when using real or 'realistic' data
+
+ Int_t fLabel; // GEANT label of corresponding particle
+ Short_t fTruePDG; // PDG code of corresponding particle
+ Int_t fMother; // GEANT label of mother particle
+ Short_t fMotherPDG; // PDG code of mother particle
+
+ ClassDef(AliRsnDaughter,1)
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnDaughterCut
+//
+// Implementation of track cuts for analysis.
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#include <Riostream.h>
+
+#include "AliRsnDaughter.h"
+#include "AliRsnDaughterCut.h"
+
+ClassImp(AliRsnDaughterCut)
+
+//--------------------------------------------------------------------------------------------------------
+Bool_t AliRsnDaughterCut::Pass(AliRsnDaughter* /*track1*/, AliRsnDaughter* /*track2*/)
+//
+// Virtual method for cut passing.
+// This function must be overridden and return kTRUE when cut is passed.
+//
+{
+ TObject::Error("Pass", "This method must be overridden!");
+ return kFALSE;
+}
+//--------------------------------------------------------------------------------------------------------
+Bool_t AliRsnDaughterCutPtSingle::Pass(AliRsnDaughter *track1, AliRsnDaughter* /*track2*/)
+//
+// Cut on single track momentum.
+//
+{
+ if (!track1) return kFALSE;
+ if (track1->GetPt() < fPtMin) return kFALSE;
+ if (track1->GetPt() > fPtMax) return kFALSE;
+
+ return kTRUE;
+}
+//--------------------------------------------------------------------------------------------------------
+Bool_t AliRsnDaughterCutPtPair::Pass(AliRsnDaughter *track1, AliRsnDaughter *track2)
+//
+// Cut on single track momentum.
+//
+{
+ if (!track1 || !track2) return kFALSE;
+
+ AliRsnDaughter sum = (*track1) + (*track2);
+
+ if (sum.GetPt() < fPtMin) return kFALSE;
+ if (sum.GetPt() > fPtMax) return kFALSE;
+
+ return kTRUE;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnDaughterCut
+//
+// Implementation of various cut which can be applied
+// during resonance analysis.
+// First is the virtual base class.
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#ifndef ALIRSNDAUGHTERCUT_H
+#define ALIRSNDAUGHTERCUT_H
+
+class AliRsnDaughter;
+
+class AliRsnDaughterCut : public TObject
+{
+public:
+ AliRsnDaughterCut() { }
+ virtual ~AliRsnDaughterCut() { }
+
+ Bool_t IsPairCut() {return fPairCut;}
+ virtual Bool_t Pass(AliRsnDaughter *track1, AliRsnDaughter *track2 = 0);
+
+protected:
+
+ Bool_t fPairCut;
+
+ ClassDef(AliRsnDaughterCut,1)
+};
+
+//-------------------------------------------------------------------------
+
+class AliRsnDaughterCutPtSingle : public AliRsnDaughterCut
+{
+public:
+ AliRsnDaughterCutPtSingle(Double_t min, Double_t max) {fPairCut=kFALSE;fPtMin=min;fPtMax=max;}
+ virtual ~AliRsnDaughterCutPtSingle() { }
+
+ virtual Bool_t Pass(AliRsnDaughter *track1, AliRsnDaughter *track2 = 0);
+
+protected:
+
+ Double_t fPtMin; // smallest allowed Pt
+ Double_t fPtMax; // largest allowed Pt
+
+ ClassDef(AliRsnDaughterCutPtSingle,1)
+};
+
+//-------------------------------------------------------------------------
+
+class AliRsnDaughterCutPtPair : public AliRsnDaughterCut
+{
+public:
+ AliRsnDaughterCutPtPair(Double_t min, Double_t max) {fPairCut=kTRUE;fPtMin=min;fPtMax=max;}
+ virtual ~AliRsnDaughterCutPtPair() { }
+
+ virtual Bool_t Pass(AliRsnDaughter *track1, AliRsnDaughter *track2);
+
+protected:
+
+ Double_t fPtMin; // smallest allowed Pt
+ Double_t fPtMax; // largest allowed Pt
+
+ ClassDef(AliRsnDaughterCutPtPair,1)
+};
+
+//-------------------------------------------------------------------------
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnEvent
+// Simple collection of reconstructed tracks, selected from an ESD event
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#include <Riostream.h>
+
+#include <TTree.h>
+#include <TString.h>
+#include <TParticle.h>
+#include <TObjArray.h>
+#include <TRefArray.h>
+#include <TClonesArray.h>
+
+#include "AliESD.h"
+#include "AliStack.h"
+#include "AliESDtrack.h"
+#include "AliRsnDaughter.h"
+#include "AliRsnEvent.h"
+
+ClassImp(AliRsnEvent)
+
+//--------------------------------------------------------------------------------------------------------
+AliRsnEvent::AliRsnEvent() :
+ fIsESD(kTRUE),
+ fPath(""),
+ fPVx(0.0),
+ fPVy(0.0),
+ fPVz(0.0),
+ fMultiplicity(-1)
+//
+// Default constructor
+//
+{
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ fPos[i] = NULL;
+ fNeg[i] = NULL;
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+AliRsnEvent::AliRsnEvent(const AliRsnEvent &event) :
+ TObject((TObject)event),
+ fIsESD(event.fIsESD),
+ fPath(event.fPath),
+ fPVx(event.fPVx),
+ fPVy(event.fPVy),
+ fPVz(event.fPVz),
+ fMultiplicity(event.fMultiplicity)
+//
+// Copy constructor.
+// Creates new instances of all collections to store a copy of all objects.
+//
+{
+ // clone tracks collections
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ fPos[i] = 0;
+ fNeg[i] = 0;
+ if (event.fPos[i]) fPos[i] = (TClonesArray*)event.fPos[i]->Clone();
+ if (event.fNeg[i]) fNeg[i] = (TClonesArray*)event.fNeg[i]->Clone();
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnEvent::AddTrack(AliRsnDaughter track)
+//
+// Stores a track into the correct array
+//
+{
+ Int_t pdg, sign, ifirst, ilast;
+
+ // if sign is zero, track is not stored
+ sign = (Int_t)track.GetSign();
+ if (!sign) return;
+
+ // if PDG code is assigned, track is stored in the corresponding collection
+ // otherwise, a copy of track is is stored in each collection (undefined track)
+ pdg = track.GetPDG();
+ if (pdg != 0) {
+ ifirst = PDG2Enum(pdg);
+ ilast = ifirst;
+ if (ifirst < AliPID::kElectron || ifirst > AliPID::kProton) return;
+ }
+ else {
+ ifirst = AliPID::kElectron;
+ ilast = AliPID::kProton;
+ }
+
+ // track is stored
+ Int_t i, index;
+ for (i = ifirst; i <= ilast; i++) {
+ if (sign > 0) {
+ index = fPos[i]->GetEntries();
+ TClonesArray &array = *fPos[i];
+ new(array[index]) AliRsnDaughter(track);
+ }
+ else {
+ index = fNeg[i]->GetEntries();
+ TClonesArray &array = *fNeg[i];
+ new(array[index]) AliRsnDaughter(track);
+ }
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnEvent::Clear(Option_t *option)
+//
+// Clears list of tracks and references.
+// If the string "DELETE" is specified, the collection classes
+// are also cleared from heap.
+//
+{
+ // evaluate option
+ TString opt(option);
+ Bool_t deleteCollections = opt.Contains("DELETE", TString::kIgnoreCase);
+
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ if (fPos[i]) fPos[i]->Delete();
+ if (fNeg[i]) fNeg[i]->Delete();
+ if (deleteCollections) {
+ delete fPos[i];
+ delete fNeg[i];
+ fPos[i] = 0;
+ fNeg[i] = 0;
+ }
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+Int_t AliRsnEvent::GetMultiplicity(Bool_t recalc)
+//
+// Computes multiplicity.
+// If it is already computed (fMultiplicity > -1), it returns that value,
+// unless one sets the argument to kTRUE.
+//
+{
+ if (fMultiplicity < 0) recalc = kTRUE;
+ if (recalc) {
+ fMultiplicity = 0;
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ fMultiplicity += (Int_t)fPos[i]->GetEntries();
+ fMultiplicity += (Int_t)fNeg[i]->GetEntries();
+ }
+ }
+
+ return fMultiplicity;
+}
+//--------------------------------------------------------------------------------------------------------
+const char * AliRsnEvent::GetOriginFileName()
+//
+// Returns the path where input file was stored
+//
+{
+ TString str(fPath);
+ if (fIsESD) {
+ str.Append("/AliESDs.root");
+ }
+ else {
+ str.Append("/galice.root");
+ }
+
+ return str.Data();
+}
+//--------------------------------------------------------------------------------------------------------
+TClonesArray* AliRsnEvent::GetTracks(Char_t sign, AliPID::EParticleType type)
+//
+// Returns the particle collection specified in argument
+//
+{
+ Int_t itype = (Int_t)type;
+ if (itype >= 0 && itype < AliPID::kSPECIES) {
+ if (sign == '+') return fPos[type]; else return fNeg[type];
+ }
+ else {
+ return NULL;
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnEvent::Init()
+//
+// Action 1: define default values for some data members (including pointers).
+// Action 2: if 'ntracks' > 0 allocates memory to store tracks.
+//
+{
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ fPos[i] = new TClonesArray("AliRsnDaughter", 0);
+ fNeg[i] = new TClonesArray("AliRsnDaughter", 0);
+ fPos[i]->BypassStreamer(kFALSE);
+ fNeg[i]->BypassStreamer(kFALSE);
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+Int_t AliRsnEvent::PDG2Enum(Int_t pdgcode)
+//
+// Converts a PDG code into the correct slot in the EParticleType enumeration in AliPID
+//
+{
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ if (AliPID::ParticleCode((AliPID::EParticleType)i) == TMath::Abs(pdgcode)) {
+ return i;
+ }
+ }
+
+ return -1;
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnEvent::PrintTracks()
+//
+// Print data for particles stored in this event
+//
+{
+ cout << endl;
+
+ AliRsnDaughter *track = 0;
+
+ Int_t type;
+ for (type = 0; type < AliPID::kSPECIES; type++) {
+ TObjArrayIter iterPos(fPos[type]);
+ cout << "Positive " << AliPID::ParticleName(type) << "s" << endl;
+ if (!fPos[type]) cout << "NOT INITIALIZED" << endl;
+ while ( (track = (AliRsnDaughter*)iterPos.Next()) ) {
+ cout << "Index, Sign, PDG = ";
+ cout << track->GetIndex() << " ";
+ cout << (Int_t)track->GetSign() << " ";
+ cout << track->GetPDG() << endl;
+ }
+ TObjArrayIter iterNeg(fNeg[type]);
+ cout << "Negative " << AliPID::ParticleName(type) << "s" << endl;
+ if (!fNeg[type]) cout << "NOT INITIALIZED" << endl;
+ while ( (track = (AliRsnDaughter*)iterNeg.Next()) ) {
+ cout << "Index, Sign, PDG = ";
+ cout << track->GetIndex() << " ";
+ cout << (Int_t)track->GetSign() << " ";
+ cout << track->GetPDG() << endl;
+ }
+ }
+}
+//--------------------------------------------------------------------------------------------------------
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnEvent
+// Simple collection of reconstructed tracks, selected from an ESD event
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#ifndef ALIRSNEVENT_H
+#define ALIRSNEVENT_H
+
+#include <TMath.h>
+#include "AliPID.h"
+
+class AliRsnDaughter;
+
+class AliRsnEvent : public TObject
+{
+public:
+ AliRsnEvent();
+ AliRsnEvent(const AliRsnEvent& copy);
+
+ virtual ~AliRsnEvent() {Clear("DELETE");}
+
+ void AddTrack(AliRsnDaughter track);
+ void Clear(Option_t *option = "");
+ Int_t GetMultiplicity(Bool_t recalc = kFALSE);
+ const char* GetOriginFileName();
+ Double_t GetPrimaryVertexX() {return fPVx;}
+ Double_t GetPrimaryVertexY() {return fPVy;}
+ Double_t GetPrimaryVertexZ() {return fPVz;}
+ void GetPrimaryVertex(Double_t &x, Double_t &y, Double_t &z) {x=fPVx;y=fPVy;z=fPVz;}
+ TClonesArray* GetTracks(Char_t sign, AliPID::EParticleType type);
+ void Init();
+ Int_t PDG2Enum(Int_t pdgcode);
+ void PrintTracks();
+ void SetESD(Bool_t yesno = kTRUE) {fIsESD=yesno;}
+ void SetPath(TString path) {fPath=path;}
+ void SetPrimaryVertex(Double_t x, Double_t y, Double_t z) {fPVx=x;fPVy=y;fPVz=z;}
+
+private:
+
+ Bool_t fIsESD; // if true, it is ESD event, otherwise it comes from Kine
+ TString fPath; // complete path where input event file is stored
+
+ Double_t fPVx; //
+ Double_t fPVy; // primary vertex
+ Double_t fPVz; //
+
+ Int_t fMultiplicity; // global event multiplicity
+
+ TClonesArray *fPos[AliPID::kSPECIES]; // collections of positive particles
+ TClonesArray *fNeg[AliPID::kSPECIES]; // collections of negative particles
+
+ ClassDef(AliRsnEvent,1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnReader
+//
+// Reader for conversion of ESD or Kinematics output into AliRsnEvent
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#include <Riostream.h>
+
+#include <TH1.h>
+#include <TH3.h>
+#include <TFile.h>
+#include <TChain.h>
+#include <TArrayF.h>
+#include <TParticle.h>
+#include <TRandom.h>
+#include <TObjString.h>
+#include <TObjectTable.h>
+#include <TOrdCollection.h>
+
+#include "AliRun.h"
+#include "AliESD.h"
+#include "AliStack.h"
+#include "AliHeader.h"
+#include "AliESDtrack.h"
+#include "AliRunLoader.h"
+#include "AliGenEventHeader.h"
+#include "AliGenPythiaEventHeader.h"
+
+#include "AliRsnDaughter.h"
+#include "AliRsnEvent.h"
+#include "AliRsnReader.h"
+
+ClassImp(AliRsnReader)
+
+//--------------------------------------------------------------------------------------------------------
+AliRsnReader::AliRsnReader()
+//
+// Constructor.
+// Initializes all working parameters to default values:
+// - ESD particle identification
+// - rejection of non-ITS-refitted tracks
+// - maximum distance allowed from primary vertex = 3.0 cm (beam pipe)
+//
+{
+ fPIDMethod = kESDPID;
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = 1.0;
+ fPtLimit4PID = 4.0;
+ fProbThreshold = 0.0;
+ fMaxRadius = 3.0; // the beam pipe
+
+ fEvents = 0;
+ for (i = AliPID::kElectron; i <= AliPID::kProton; i++) {
+ fEffPos[i] = 0;
+ fEffNeg[i] = 0;
+ }
+ fResPt = 0;
+ fResLambda = 0;
+ fResPhi = 0;
+}
+//--------------------------------------------------------------------------------------------------------
+AliRsnReader::AliRsnReader(const AliRsnReader ©) : TObject(copy)
+//
+// Copy constructor.
+// Initializes all working parameters to same values of another simlar object.
+// TTree data member is not created.
+//
+{
+ fPIDMethod = copy.fPIDMethod;
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = copy.fPrior[i];
+ fPtLimit4PID = copy.fPtLimit4PID;
+ fProbThreshold = copy.fProbThreshold;
+ fMaxRadius = copy.fMaxRadius;
+
+ fEvents = 0;
+ for (i = AliPID::kElectron; i <= AliPID::kProton; i++) {
+ fEffPos[i] = copy.fEffPos[i];
+ fEffNeg[i] = copy.fEffNeg[i];
+ }
+
+ fResPt = copy.fResPt;
+ fResLambda = copy.fResLambda;
+ fResPhi = copy.fResPhi;
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnReader::Clear(Option_t *option)
+//
+// Clear collection of filenames.
+// If requested with the option "DELETELIST",
+// the collection object is also deleted.
+//
+{
+ TString opt(option);
+
+ if (!opt.CompareTo("TREE", TString::kIgnoreCase)) {
+ fEvents->Reset();
+ if (!opt.CompareTo("DELTREE", TString::kIgnoreCase)) {
+ delete fEvents;
+ fEvents = 0;
+ }
+ }
+
+ Int_t i;
+ for (i = AliPID::kElectron; i <= AliPID::kProton; i++) {
+ fEffPos[i] = 0;
+ fEffNeg[i] = 0;
+ }
+ fResPt = 0;
+ fResLambda = 0;
+ fResPhi = 0;
+}
+//--------------------------------------------------------------------------------------------------------
+Bool_t AliRsnReader::EffSim(Int_t pdg, Double_t pt, Double_t eta, Double_t z)
+//
+// If efficiency histogram are present, they are used to simulate efficiency.
+// Pt, Eta and Z are used to find reconstruction efficiency value to be used
+// and PDG is used to select efficiency histogram for a given particle.
+// An extraction is done, and it is supposed that particle must be accepted
+// only when this function retunrs kTRUE (= successful extraction).
+//
+{
+ // find particle sign from PDG code
+ Int_t sign;
+ if (TMath::Abs(pdg) >= 211) {
+ if (pdg > 0) sign = 1; else sign = -1;
+ }
+ else {
+ if (pdg > 0) sign = -1; else sign = 1;
+ }
+
+ // convert PDG code into one value in AliPID::kSPECIES
+ // (if returned value is not a charged particle, procedure is skipped)
+ Int_t index = FindType(pdg);
+ TH3D *eff = 0;
+ if (index >= AliPID::kElectron && index <= AliPID::kProton) {
+ if (sign > 0) eff = fEffPos[index]; else eff = fEffNeg[index];
+ }
+
+ // if no efficiency histogram is defined, method returns a 'fail' result
+ if (!eff) return kFALSE;
+
+ // otherwise, a random extraction is made
+ Int_t ibin = eff->FindBin(z, eta, pt);
+ Double_t ref = (Double_t)eff->GetBinContent(ibin);
+ Double_t ran = gRandom->Rndm();
+ return (ran <= ref);
+}
+//--------------------------------------------------------------------------------------------------------
+Double_t* AliRsnReader::GetPIDprobabilities(AliRsnDaughter track)
+//
+// Computes PID probabilites starting from priors and weights
+//
+{
+ Double_t *prob = new Double_t[AliPID::kSPECIES];
+
+ Int_t i;
+
+ // step 1 - compute the normalization factor
+ Double_t sum = 0.0;
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ prob[i] = fPrior[i] * track.GetPIDweight(i);
+ sum += prob[i];
+ }
+ if (sum <= 0.0) return 0;
+
+ // step 2 - normalize PID weights by the relative prior probability
+ for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
+ prob[i] /= sum;
+ }
+
+ return prob;
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnReader::Identify(AliRsnDaughter &track)
+//
+// Identifies a track according to method selected
+//
+{
+ Bool_t doESDPID = (fPIDMethod == kESDPID);
+ Bool_t keepRecSign = (fPIDMethod == kPerfectPIDwithRecSign);
+
+ if (doESDPID) {
+ // when doing ESD PID it is supposed that charge sign
+ // comes from ESD track and it is not modified
+ Double_t pt = track.GetPt();
+ if (pt <= fPtLimit4PID) {
+ Double_t *prob = GetPIDprobabilities(track);
+ if (!prob) track.SetPDG(0);
+ Int_t idx[AliPID::kSPECIES];
+ TMath::Sort(AliPID::kSPECIES, prob, idx);
+ Int_t imax = idx[0];
+ Double_t maxprob = prob[imax];
+ if (maxprob >= fProbThreshold) {
+ track.SetPDG((UShort_t)AliPID::ParticleCode(imax));
+ }
+ delete [] prob;
+ }
+ else {
+ track.SetPDG(0);
+ }
+ }
+ else {
+ Short_t truePDG = track.GetTruePDG();
+ track.SetPDG((UShort_t)TMath::Abs(truePDG));
+ if (!keepRecSign) {
+ if (TMath::Abs(truePDG) <= 13) {
+ if (truePDG > 0) track.SetSign(-1); else track.SetSign(1);
+ }
+ else {
+ if (truePDG > 0) track.SetSign(1); else track.SetSign(-1);
+ }
+ }
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+TTree* AliRsnReader::ReadParticles(const char *path, Option_t *option)
+//
+// Opens the file "galice.root" and kinematics in the path specified,
+// loads Kinematics informations and fills and AliRsnEvent object
+// with data coming from simulation.
+// If required, an efficiency simulation can be done which cuts off some tracks
+// depending on a MonteCarlo extraction weighted on the efficiency histograms
+// passed to the class.
+// Moreover, a smearing can be done if requested, using the resolution histograms
+// passed to the class.
+// Allowed options:
+// - "E" --> do efficiency simulation
+// - "P" --> do momentum smearing
+//
+{
+ fPIDMethod = kPerfectPID;
+
+ TTree *events = new TTree("selection", "AliRsnEvents");
+ TTree::SetBranchStyle(1);
+ AliRsnEvent *event = new AliRsnEvent;
+ TBranch *branch = events->Branch("events", "AliRsnEvent", &event, 32000, 1);
+ branch->SetAutoDelete(kFALSE);
+
+ // get path
+ TString strPath(path);
+ if (strPath.Length() < 1) return 0;
+ strPath += '/';
+
+ // evaluate options
+ TString opt(option);
+ opt.ToUpper();
+ Bool_t doEffSim = (opt.Contains("E"));
+ Bool_t doSmearing = (opt.Contains("P"));
+
+ // initialize run loader
+ AliRunLoader *runLoader = OpenRunLoader(path);
+ if (!runLoader) return 0;
+ Int_t nEvents = runLoader->GetNumberOfEvents();
+ TArrayF vertex(3);
+
+ // loop on events
+ Int_t i, iev, index = 0, imum, nParticles, nStoredParticles = 0;
+ Double_t vtot, px, py, pz, arg, eta;
+ TParticle *particle = 0, *mum = 0;
+ for (iev = 0; iev < nEvents; iev++) {
+
+ // load given event
+ cout << "\rEvent " << iev << " of " << nEvents << flush;
+ runLoader->GetEvent(iev);
+
+ // primary vertex
+ runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
+ //cout << "Process type: " << ((AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader())->ProcessType() << endl;
+
+ // kinematics
+ AliStack *stack = runLoader->Stack();
+ if (!stack) continue;
+
+ // get number of particles to collect
+ nParticles = stack->GetNtrack();
+ if (!nParticles) return 0;
+
+ // create new AliRsnEvent object
+ event->Clear("DELETE");
+ event->SetESD(kFALSE);
+ event->SetPath(strPath);
+ event->Init();
+
+ // loop on tracks
+ index = 0;
+ for (i = 0; i < nParticles; i++) {
+
+ // get particle
+ particle = stack->Particle(i);
+ if (!particle) continue;
+
+ // check against maximum allowed distance from primary vertex
+ vtot = (particle->Vx() - vertex[0])*(particle->Vx() - vertex[0]);
+ vtot += (particle->Vy() - vertex[1])*(particle->Vy() - vertex[1]);
+ vtot += (particle->Vz() - vertex[2])*(particle->Vz() - vertex[2]);
+ vtot = TMath::Sqrt(vtot);
+ if (vtot > fMaxRadius) continue;
+
+ // efficiency selection
+ if (doEffSim) {
+ arg = 0.5 * TMath::ATan2(particle->Pz(), particle->Pt());
+ if (arg > 0.) eta = -TMath::Log(arg); else continue;
+ Bool_t result = EffSim(particle->GetPdgCode(), (Double_t)vertex[2], eta, (Double_t)particle->Pt());
+ if (result == kFALSE) continue;
+ }
+
+ // smear kinematic parameters (if requested)
+ px = particle->Px();
+ py = particle->Py();
+ pz = particle->Pz();
+ if (doSmearing) SmearMomentum(px, py, pz);
+
+ // add to collection of tracks
+ nStoredParticles++;
+ AliRsnDaughter *track = new AliRsnDaughter;
+ if (!track->Adopt(particle)) {
+ delete track;
+ continue;
+ }
+ track->SetIndex(index);
+ track->SetLabel(i);
+ track->SetPxPyPz(px, py, pz);
+ imum = particle->GetFirstMother();
+ if (imum >= 0) {
+ mum = stack->Particle(imum);
+ track->SetMotherPDG( (Short_t)mum->GetPdgCode() );
+ }
+ Identify(*track);
+
+ event->AddTrack(*track);
+ delete track;
+
+ } // end of loop over tracks
+
+ // compute total multiplicity
+ event->GetMultiplicity();
+
+ // link to events tree and fill
+ events->Fill();
+ }
+
+ runLoader->UnloadKinematics();
+ delete runLoader;
+
+ return events;
+}
+//--------------------------------------------------------------------------------------------------------
+TTree* AliRsnReader::ReadTracks(const char *path, Option_t *option)
+//
+// Reads AliESD in a given path, with and "experimental" method.
+// When using this method, no Kinematics information is assumed
+// and particle identification is always done with the bayesian method.
+// No Kinematics informations are stored.
+// Allowed options (actually):
+// - "R" : reject tracks which do not have the flag "kITSRefit" set to TRUE
+// - "F" : reject 'fake' tracks (negative label)
+//
+{
+ fPIDMethod = kESDPID;
+
+ TTree *events = new TTree("selection", "AliRsnEvents");
+ TTree::SetBranchStyle(1);
+ AliRsnEvent *event = new AliRsnEvent;
+ TBranch *branch = events->Branch("events", "AliRsnEvent", &event, 32000, 1);
+ branch->SetAutoDelete(kFALSE);
+
+ // get path
+ TString strPath(path);
+ if (strPath.Length() < 1) return 0;
+ strPath += '/';
+
+ // evaluate options
+ TString opt(option);
+ opt.ToUpper();
+ Bool_t checkITSrefit = (opt.Contains("R"));
+ Bool_t rejectFakes = (opt.Contains("F"));
+
+ // opens ESD file
+ TFile *fileESD = TFile::Open(Form("%s/AliESDs.root", strPath.Data()));
+ if (!fileESD) return 0;
+ if (fileESD->IsOpen() == kFALSE) return 0;
+ TTree* treeESD = (TTree*)fileESD->Get("esdTree");
+ AliESD *esd = 0;
+ treeESD->SetBranchAddress("ESD", &esd);
+ Int_t nev = (Int_t)treeESD->GetEntries();
+
+ // loop on events
+ Int_t i, nSelTracks = 0;
+ Double_t vertex[3];
+ for (i = 0; i < nev; i++) {
+
+ // message
+ cout << "\rEvent " << i << flush;
+ treeESD->GetEntry(i);
+
+ // primary vertex
+ vertex[0] = esd->GetVertex()->GetXv();
+ vertex[1] = esd->GetVertex()->GetYv();
+ vertex[2] = esd->GetVertex()->GetZv();
+
+ // create new AliRsnEvent object
+ event->Clear("DELETE");
+ event->Init();
+ event->SetPath(strPath);
+ event->SetESD();
+
+ // get number of tracks
+ Int_t ntracks = esd->GetNumberOfTracks();
+ if (!ntracks) continue;
+
+ // store tracks from ESD
+ Int_t index, label;
+ Double_t vtot, v[3];
+ for (index = 0; index < ntracks; index++) {
+
+ // get track
+ AliESDtrack *esdTrack = esd->GetTrack(index);
+
+ // check against vertex constraint
+ esdTrack->GetXYZ(v);
+ vtot = (v[0] - vertex[0])*(v[0] - vertex[0]);
+ vtot += (v[1] - vertex[1])*(v[1] - vertex[1]);
+ vtot += (v[2] - vertex[2])*(v[2] - vertex[2]);
+ vtot = TMath::Sqrt(vtot);
+ if (vtot > fMaxRadius) continue;
+
+ // check for ITS refit
+ if (checkITSrefit) {
+ if (!(esdTrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
+ }
+
+ // check for fakes
+ label = esdTrack->GetLabel();
+ if (rejectFakes && (label < 0)) continue;
+
+ // create AliRsnDaughter (and make Bayesian PID)
+ AliRsnDaughter track;
+ if (!track.Adopt(esdTrack, checkITSrefit)) continue;
+ track.SetIndex(index);
+ track.SetLabel(label);
+ Identify(track);
+
+ // store in TClonesArray
+ event->AddTrack(track);
+ nSelTracks++;
+ }
+
+ // compute total multiplicity
+ event->GetMultiplicity();
+
+ // link to events tree and fill
+ events->Fill();
+ }
+
+ fileESD->Close();
+ return events;
+}
+//--------------------------------------------------------------------------------------------------------
+TTree* AliRsnReader::ReadTracksAndParticles(const char *path, Option_t *option)
+//
+// Reads AliESD in a given path, getting also informations from Kinematics.
+// In this case, the PID method used is the one selected with apposite setter.
+// Allowed options (actually):
+// - "R" : reject tracks which do not have the flag "kITSRefit" set to TRUE
+// - "F" : reject 'fake' tracks (negative label)
+// - "M" : use 'true' momentum instead of reconstructed one
+//
+{
+ TTree *events = new TTree("selection", "AliRsnEvents");
+ TTree::SetBranchStyle(1);
+ AliRsnEvent *event = new AliRsnEvent;
+ TBranch *branch = events->Branch("events", "AliRsnEvent", &event, 32000, 1);
+ branch->SetAutoDelete(kFALSE);
+
+ // get path
+ TString strPath(path);
+ if (strPath.Length() < 1) return 0;
+
+ // evaluate options
+ TString opt(option);
+ opt.ToUpper();
+ Bool_t checkITSrefit = (opt.Contains("R"));
+ Bool_t rejectFakes = (opt.Contains("F"));
+ Bool_t copyMomentum = (opt.Contains("M"));
+
+ // opens ESD file
+ TFile *fileESD = TFile::Open(Form("%s/AliESDs.root", strPath.Data()));
+ if (!fileESD) return 0;
+ if (fileESD->IsOpen() == kFALSE) return 0;
+ TTree* treeESD = (TTree*)fileESD->Get("esdTree");
+ AliESD *esd = 0;
+ treeESD->SetBranchAddress("ESD", &esd);
+ Int_t nevRec = (Int_t)treeESD->GetEntries();
+
+ // initialize run loader
+ AliRunLoader *runLoader = OpenRunLoader(path);
+ if (!runLoader) return 0;
+ Int_t nevSim = runLoader->GetNumberOfEvents();
+
+ // check number of reconstructed and simulated events
+ if ( (nevSim != 0 && nevRec != 0) && (nevSim != nevRec) ) {
+ cerr << "Count mismatch: sim = " << nevSim << " -- rec = " << nevRec << endl;
+ return 0;
+ }
+ else if (nevSim == 0 && nevRec == 0) {
+ cerr << "Count error: sim = rec = 0" << endl;
+ return 0;
+ }
+
+ // loop on events
+ Int_t i, procType, ntracks, nSelTracks = 0;
+ Double_t vertex[3];
+ for (i = 0; i < nevRec; i++) {
+
+ // get event
+ cout << "\rEvent " << i << " " << flush;
+ treeESD->GetEntry(i);
+ runLoader->GetEvent(i);
+
+ // reject event if it is diffractive
+ procType = ((AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader())->ProcessType();
+ if (procType == 92 || procType == 93 || procType == 94) {
+ cout << "Skipping diffractive event" << endl;
+ continue;
+ }
+
+ // get particle stack
+ AliStack *stack = runLoader->Stack();
+
+ // primary vertex
+ vertex[0] = esd->GetVertex()->GetXv();
+ vertex[1] = esd->GetVertex()->GetYv();
+ vertex[2] = esd->GetVertex()->GetZv();
+
+ // multiplicity
+ ntracks = esd->GetNumberOfTracks();
+ if (!ntracks) {
+ Warning("ReadTracksAndParticles", "Event %d has no tracks!", i);
+ continue;
+ }
+
+ // create new AliRsnEvent object
+ event->Clear("DELETE");
+ event->Init();
+ event->SetPath(strPath);
+ event->SetESD();
+
+ // store tracks from ESD
+ Int_t index, label;
+ Double_t vtot, v[3];
+ for (index = 0; index < ntracks; index++) {
+
+ // get track
+ AliESDtrack *esdTrack = esd->GetTrack(index);
+
+ // check against vertex constraint
+ esdTrack->GetXYZ(v);
+ vtot = (v[0] - vertex[0])*(v[0] - vertex[0]);
+ vtot += (v[1] - vertex[1])*(v[1] - vertex[1]);
+ vtot += (v[2] - vertex[2])*(v[2] - vertex[2]);
+ vtot = TMath::Sqrt(vtot);
+ if (vtot > fMaxRadius) continue;
+
+ // check for fakes
+ label = esdTrack->GetLabel();
+ if (rejectFakes) {
+ if (label < 0) continue;
+ }
+
+ // create AliRsnDaughter (and make Bayesian PID)
+ AliRsnDaughter track;
+ if (!track.Adopt(esdTrack, checkITSrefit)) continue;
+ track.SetIndex(index);
+
+ // retrieve particle and get Kine info
+ TParticle *part = stack->Particle(TMath::Abs(label));
+ track.SetTruePDG(part->GetPdgCode());
+ Int_t mother = part->GetFirstMother();
+ track.SetMother(mother);
+ if (mother >= 0) {
+ TParticle *mum = stack->Particle(mother);
+ track.SetMotherPDG(mum->GetPdgCode());
+ }
+ if (copyMomentum) track.SetPxPyPz(part->Px(), part->Py(), part->Pz());
+
+ // identification
+ Identify(track);
+
+ // store in TClonesArray
+ event->AddTrack(track);
+ nSelTracks++;
+ }
+
+ // compute total multiplicity
+ event->GetMultiplicity();
+
+ // link to events tree and fill
+ events->Fill();
+ }
+
+ runLoader->UnloadKinematics();
+ delete runLoader;
+ fileESD->Close();
+
+ return events;
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnReader::SetEfficiencyHistogram(AliPID::EParticleType type, TH3D *h, Bool_t pos)
+//
+// Sets efficiency histogram for a given particle species.
+// If third argument is "true", histo is assigned to positive particles,
+// otherwise it is assigned to negative particles.
+//
+{
+ if (type >= AliPID::kElectron && type < AliPID::kPhoton) {
+ if (pos) fEffPos[type] = h; else fEffNeg[type] = h;
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnReader::SetPriorProbabilities(Double_t *prior)
+//
+// Set prior probabilities to be used in case of ESD PID.
+//
+{
+ if (!prior) return;
+
+ Int_t i = 0;
+ for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = prior[i];
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnReader::SetPriorProbability(AliPID::EParticleType type, Double_t value)
+//
+// Sets prior probability referred to a single particle species.
+//
+{
+ if (type >= AliPID::kElectron && type < AliPID::kPhoton) fPrior[type] = value;
+}
+//--------------------------------------------------------------------------------------------------------
+void AliRsnReader::SmearMomentum(Double_t &px, Double_t &py, Double_t &pz)
+//
+// Use resolution histograms to do smearing of momentum
+// (to introduce reconstruction effects)
+//
+{
+ Double_t pt = TMath::Sqrt(px*px + py*py);
+ Double_t lambda = TMath::ATan2(pz, pt);
+ Double_t phi = TMath::ATan2(py, px);
+
+ Int_t ibin;
+ Double_t ref;
+ if (fResPt) {
+ ibin = fResPt->FindBin(pt);
+ ref = (Double_t)fResPt->GetBinContent(ibin);
+ pt *= 1.0 + ref;
+ }
+ if (fResLambda) {
+ ibin = fResLambda->FindBin(pt);
+ ref = (Double_t)fResLambda->GetBinContent(ibin);
+ lambda *= 1.0 + ref;
+ }
+ if (fResPhi) {
+ ibin = fResPhi->FindBin(pt);
+ ref = (Double_t)fResPhi->GetBinContent(ibin);
+ phi *= 1.0 + ref;
+ }
+
+ px = pt * TMath::Cos(phi);
+ py = pt * TMath::Sin(phi);
+ pz = pt * TMath::Tan(lambda);
+}
+//--------------------------------------------------------------------------------------------------------
+AliPID::EParticleType AliRsnReader::FindType(Int_t pdg)
+//
+// Finds the enum value corresponding to a PDG code
+//
+{
+ pdg = TMath::Abs(pdg);
+ switch (pdg) {
+ case 11: return AliPID::kElectron; break;
+ case 13: return AliPID::kMuon; break;
+ case 211: return AliPID::kPion; break;
+ case 321: return AliPID::kKaon; break;
+ case 2212: return AliPID::kProton; break;
+ default : return AliPID::kPhoton;
+ }
+}
+//--------------------------------------------------------------------------------------------------------
+AliRunLoader* AliRsnReader::OpenRunLoader(const char *path)
+//
+// Open the Run loader with events in a given path
+//
+{
+ // clear gALICE
+ if (gAlice) {
+ delete gAlice;
+ gAlice = 0;
+ }
+
+ // initialize run loader
+ TString name(path);
+ name += "/galice.root";
+ AliRunLoader *runLoader = AliRunLoader::Open(name.Data());
+ if (runLoader) {
+ runLoader->LoadgAlice();
+ gAlice = runLoader->GetAliRun();
+ runLoader->LoadKinematics();
+ runLoader->LoadHeader();
+ }
+
+ return runLoader;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// Class AliRsnEventReader
+//
+// Reader for conversion of ESD or Kinematics output into AliRsnEvent
+//
+// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
+//-------------------------------------------------------------------------
+
+#ifndef AliRSNREADER_H
+#define AliRSNREADER_H
+
+#include "AliPID.h"
+
+class TH3D;
+class TH1D;
+class TOrdCollection;
+class TTree;
+class AliRunLoader;
+
+class AliRsnReader : public TObject
+{
+public:
+
+ enum EPIDMethod {
+ kNoPID = 0, // no PID performed
+ kPerfectPID = 1, // use Kinematics to simulate a 100% PID efficiency
+ kESDPID = 2, // use experimental ESD weights
+ kPerfectPIDwithRecSign = 3 // get particle type from Kine and charge sign from reconstruction
+ };
+
+ AliRsnReader();
+ AliRsnReader(const AliRsnReader& copy);
+ virtual ~AliRsnReader() {Clear("DELTREE");}
+
+ void Clear(Option_t *option = "");
+ Bool_t EffSim(Int_t pdg, Double_t pt, Double_t eta, Double_t z);
+ TTree* GetEvents() const {return fEvents;}
+ Double_t* GetPIDprobabilities(AliRsnDaughter track);
+ void Identify(AliRsnDaughter &track);
+ TTree* ReadParticles(const char *path, Option_t *option="");
+ TTree* ReadTracks(const char *path, Option_t *option="R");
+ TTree* ReadTracksAndParticles(const char *path, Option_t *option="R");
+ void SetEfficiencyHistogram(AliPID::EParticleType type, TH3D *h, Bool_t pos=kTRUE);
+ void SetMaxRadius(Double_t value) {fMaxRadius=value;}
+ void SetPIDMethod(AliRsnReader::EPIDMethod pm) {fPIDMethod=pm;}
+ void SetPriorProbabilities(Double_t *prior);
+ void SetPriorProbability(AliPID::EParticleType type, Double_t p);
+ void SetProbabilityThreshold(Double_t p) {fProbThreshold=p;}
+ void SetPtLimit4PID(Double_t value) {fPtLimit4PID=value;}
+ void SetResolutionPt(TH1D *histo) {fResPt=histo;}
+ void SetResolutionLambda(TH1D *histo) {fResLambda=histo;}
+ void SetResolutionPhi(TH1D *histo) {fResPhi=histo;}
+ void SmearMomentum(Double_t &px, Double_t &py, Double_t &pz);
+
+protected:
+
+ AliPID::EParticleType FindType(Int_t pdg);
+ AliRunLoader* OpenRunLoader(const char *path);
+
+ Bool_t fDoEffSim; // set to true to introduce efficiency effects
+ Bool_t fDoSmearing; // when selecting only particles, set this to true
+ // to smear particle momentum
+
+ EPIDMethod fPIDMethod; // PID method
+ Double_t fPrior[AliPID::kSPECIES]; // prior probabilities (in case of REAL pid)
+ Double_t fPtLimit4PID; // maximum transverse momentum to accept realistic PID
+ Double_t fProbThreshold; // minimum required probability to accept realistic PID
+ Double_t fMaxRadius; // maximum allowed distance from primary vertex
+
+ TTree *fEvents; //! tree of read events
+ TH3D *fEffPos[AliPID::kSPECIES]; // tracking efficiency (in pt, eta and Zv) for positive particles
+ TH3D *fEffNeg[AliPID::kSPECIES]; // tracking efficiency (in pt, eta and Zv) for negative particles
+ TH1D *fResPt; // relative Pt resolution
+ TH1D *fResLambda; // relative Lambda resolution
+ TH1D *fResPhi; // relative Phi resolution
+
+ // Rsn event reader implementation
+ ClassDef(AliRsnReader,1)
+};
+
+#endif
+
--- /dev/null
+#ifdef __CINT__
+
+#pragma link off all glols;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliRsnDaughter+;
+#pragma link C++ class AliRsnDaughterCut+;
+#pragma link C++ class AliRsnDaughterCutPtSingle+;
+#pragma link C++ class AliRsnDaughterCutPtPair+;
+#pragma link C++ class AliRsnEvent+;
+#pragma link C++ class AliRsnReader+;
+#pragma link C++ class AliRsnAnalysis+;
+
+#endif
--- /dev/null
+SRCS= AliRsnDaughter.cxx AliRsnDaughterCut.cxx AliRsnEvent.cxx \
+ AliRsnReader.cxx AliRsnAnalysis.cxx
+
+HDRS= $(SRCS:.cxx=.h)
+
+DHDR:=PWG2LinkDef.h
+
+EXPORT:=
+
+EINCLUDE:= PYTHIA6
VZERO/module.mk: VZERO/libVZERObase.pkg VZERO/libVZEROsim.pkg VZERO/libVZEROrec.pkg
ZDC/module.mk: ZDC/libZDCbase.pkg ZDC/libZDCrec.pkg ZDC/libZDCsim.pkg
EVE/module.mk: EVE/libReve.pkg EVE/binreve.pkg EVE/libAlieve.pkg EVE/binalieve.pkg
+PWG0/module.mk: PWG0/libPWG0base.pkg PWG0/libPWG0selectors.pkg
+PWG2/module.mk: PWG2/libPWG2.pkg
+PWG3/module.mk: PWG3/libPWG3base.pkg
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff