Adding AliRsnReader in the RESONANCES directory
authorpanos <panos@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 6 Jun 2006 16:31:02 +0000 (16:31 +0000)
committerpanos <panos@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 6 Jun 2006 16:31:02 +0000 (16:31 +0000)
PWG2/RESONANCES/AliRsnReader.cxx [new file with mode: 0644]
PWG2/RESONANCES/AliRsnReader.h [new file with mode: 0644]

diff --git a/PWG2/RESONANCES/AliRsnReader.cxx b/PWG2/RESONANCES/AliRsnReader.cxx
new file mode 100644 (file)
index 0000000..5453951
--- /dev/null
@@ -0,0 +1,729 @@
+/**************************************************************************
+ * 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 &copy) : 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;
+}
diff --git a/PWG2/RESONANCES/AliRsnReader.h b/PWG2/RESONANCES/AliRsnReader.h
new file mode 100644 (file)
index 0000000..c70b91f
--- /dev/null
@@ -0,0 +1,87 @@
+/**************************************************************************
+ * 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
+