Resonance analysis (A.Pulvirenti)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Jun 2006 20:36:43 +0000 (20:36 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Jun 2006 20:36:43 +0000 (20:36 +0000)
14 files changed:
Makefile
PWG2/AliRsnAnalysis.cxx [new file with mode: 0644]
PWG2/AliRsnAnalysis.h [new file with mode: 0644]
PWG2/AliRsnDaughter.cxx [new file with mode: 0644]
PWG2/AliRsnDaughter.h [new file with mode: 0644]
PWG2/AliRsnDaughterCut.cxx [new file with mode: 0644]
PWG2/AliRsnDaughterCut.h [new file with mode: 0644]
PWG2/AliRsnEvent.cxx [new file with mode: 0644]
PWG2/AliRsnEvent.h [new file with mode: 0644]
PWG2/AliRsnReader.cxx [new file with mode: 0644]
PWG2/AliRsnReader.h [new file with mode: 0644]
PWG2/PWG2LinkDef.h [new file with mode: 0644]
PWG2/libPWG2.pkg [new file with mode: 0644]
build/module.dep

index d369900..0844026 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -103,6 +103,10 @@ ifeq ($(findstring PWG0,$(MAKECMDGOALS)),PWG0)
 ALIROOTMODULES += PWG0
 endif
 
+ifeq ($(findstring PWG2,$(MAKECMDGOALS)),PWG2)
+ALIROOTMODULES += PWG2
+endif
+
 ifeq ($(findstring PWG3,$(MAKECMDGOALS)),PWG3)
 ALIROOTMODULES += PWG3
 endif
diff --git a/PWG2/AliRsnAnalysis.cxx b/PWG2/AliRsnAnalysis.cxx
new file mode 100644 (file)
index 0000000..4eedede
--- /dev/null
@@ -0,0 +1,494 @@
+/**************************************************************************
+ * 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 ("???");
+       }
+}
+//--------------------------------------------------------------------------------------------------------
diff --git a/PWG2/AliRsnAnalysis.h b/PWG2/AliRsnAnalysis.h
new file mode 100644 (file)
index 0000000..15576bd
--- /dev/null
@@ -0,0 +1,113 @@
+/**************************************************************************
+ * 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
diff --git a/PWG2/AliRsnDaughter.cxx b/PWG2/AliRsnDaughter.cxx
new file mode 100644 (file)
index 0000000..cc85cf6
--- /dev/null
@@ -0,0 +1,222 @@
+/**************************************************************************
+ * 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 &copy) : 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;
+}
diff --git a/PWG2/AliRsnDaughter.h b/PWG2/AliRsnDaughter.h
new file mode 100644 (file)
index 0000000..124a412
--- /dev/null
@@ -0,0 +1,100 @@
+/**************************************************************************
+ * 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 &copy);
+                               
+       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
diff --git a/PWG2/AliRsnDaughterCut.cxx b/PWG2/AliRsnDaughterCut.cxx
new file mode 100644 (file)
index 0000000..f79b3ab
--- /dev/null
@@ -0,0 +1,67 @@
+/**************************************************************************
+ * 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;
+}
diff --git a/PWG2/AliRsnDaughterCut.h b/PWG2/AliRsnDaughterCut.h
new file mode 100644 (file)
index 0000000..a72bf52
--- /dev/null
@@ -0,0 +1,75 @@
+/**************************************************************************
+ * 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
diff --git a/PWG2/AliRsnEvent.cxx b/PWG2/AliRsnEvent.cxx
new file mode 100644 (file)
index 0000000..d8751fb
--- /dev/null
@@ -0,0 +1,257 @@
+/**************************************************************************
+ * 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;
+               }
+       }
+}
+//--------------------------------------------------------------------------------------------------------
diff --git a/PWG2/AliRsnEvent.h b/PWG2/AliRsnEvent.h
new file mode 100644 (file)
index 0000000..11390ff
--- /dev/null
@@ -0,0 +1,62 @@
+/**************************************************************************
+ * 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
diff --git a/PWG2/AliRsnReader.cxx b/PWG2/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/AliRsnReader.h b/PWG2/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
+
diff --git a/PWG2/PWG2LinkDef.h b/PWG2/PWG2LinkDef.h
new file mode 100644 (file)
index 0000000..a02fcd6
--- /dev/null
@@ -0,0 +1,15 @@
+#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
diff --git a/PWG2/libPWG2.pkg b/PWG2/libPWG2.pkg
new file mode 100644 (file)
index 0000000..e4fca0d
--- /dev/null
@@ -0,0 +1,10 @@
+SRCS= AliRsnDaughter.cxx AliRsnDaughterCut.cxx AliRsnEvent.cxx \
+       AliRsnReader.cxx AliRsnAnalysis.cxx
+
+HDRS= $(SRCS:.cxx=.h) 
+
+DHDR:=PWG2LinkDef.h
+
+EXPORT:=
+
+EINCLUDE:= PYTHIA6
index e3fd6ef..990d47f 100644 (file)
@@ -46,3 +46,6 @@ TRD/module.mk:                 TRD/libTRDbase.pkg TRD/libTRDsim.pkg TRD/libTRDrec.pkg TRD/libT
 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