--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+// Implementation of the ITS PID class
+// Very naive one... Should be made better by the detector experts...
+// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-----------------------------------------------------------------
+
+#include "AliITSpidESD.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+
+ClassImp(AliITSpidESD)
+
+//_________________________________________________________________________
+AliITSpidESD::AliITSpidESD(Double_t *param)
+{
+ //
+ // The main constructor
+ //
+ fMIP=param[0];
+ fRes=param[1];
+ fRange=param[2];
+}
+
+Double_t AliITSpidESD::Bethe(Double_t bg) {
+ //
+ // This is the Bethe-Bloch function normalised to 1 at the minimum
+ //
+ Double_t bg2=bg*bg;
+ Double_t bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2));
+ return bethe/11.091;
+}
+
+//_________________________________________________________________________
+Int_t AliITSpidESD::MakePID(AliESD *event)
+{
+ //
+ // This function calculates the "detector response" PID probabilities
+ //
+ static const Double_t masses[]={
+ 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
+ };
+ Int_t ntrk=event->GetNumberOfTracks();
+ for (Int_t i=0; i<ntrk; i++) {
+ AliESDtrack *t=event->GetTrack(i);
+ if ((t->GetStatus()&AliESDtrack::kITSin )==0)
+ if ((t->GetStatus()&AliESDtrack::kITSout)==0) continue;
+ Int_t ns=AliESDtrack::kSPECIES;
+ Double_t p[10];
+ for (Int_t j=0; j<ns; j++) {
+ Double_t mass=masses[j];
+ Double_t mom=t->GetP();
+ Double_t dedx=t->GetITSsignal()/fMIP;
+ Double_t bethe=Bethe(mom/mass);
+ Double_t sigma=fRes*bethe;
+ if (TMath::Abs(dedx-bethe) > fRange*sigma) {
+ p[j]=0.;
+ continue;
+ }
+ p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma));
+ }
+ t->SetITSpid(p);
+ }
+ return 0;
+}
--- /dev/null
+#ifndef ALIITSpIDESD_H
+#define ALIITSpIDESD_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//-------------------------------------------------------
+// ITS PID class
+// A very naive design... Should be made better by the detector experts...
+// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-------------------------------------------------------
+#include <Rtypes.h>
+
+class AliESD;
+
+class AliITSpidESD {
+public:
+ AliITSpidESD(Double_t *param);
+ Int_t MakePID(AliESD *event);
+ static Double_t Bethe(Double_t bg);
+private:
+ Double_t fMIP; // dEdx for MIP
+ Double_t fRes; // relative dEdx resolution
+ Double_t fRange; // one particle type PID range (in sigmas)
+ ClassDef(AliITSpidESD,1) // ITS PID class
+};
+
+#endif
+
+
#pragma link C++ class AliITSclusterSSD+;
#pragma link C++ class AliITSpackageSSD+;
#pragma link C++ class AliITSPid+;
+#pragma link C++ class AliITSpidESD+;
#pragma link C++ class AliITStrackV2Pid+;
// Classes used for Tracking
#pragma link C++ class AliITSTrackV1+;
AliITSFindClustersV2.cxx \
AliITSRiemannFit.cxx \
AliITSFDigitizer.cxx \
- AliITSDDLRawData.cxx
+ AliITSDDLRawData.cxx AliITSpidESD.cxx
# AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
HDRS:= $(SRCS:.cxx=.h)
//********************************************************************
-// Example (very basic for the moment) of the data analysis
+// Example (very naive for the moment) of the data analysis
// using the ESD classes
-//
+// It demonstrates the idea of the "combined PID".
// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
//********************************************************************
#endif
Int_t AliESDanalysis(Int_t nev=1) {
- TH2F *tpcHist=new TH2F("tpcHist","dE/dX vs momentum",50,0.,2.,50,0.,400.);
- TH2F *itsHist=new TH2F("itsHits","dE/dX vs momentum",50,0.,2.,50,0.,200.);
+ TH2F *tpcHist=
+ new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
+ tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
+ tpcHist->SetMarkerStyle(8);
+ tpcHist->SetMarkerSize(0.3);
+
+ TH2F *elHist=new TH2F("elHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
+ elHist->SetMarkerStyle(8);
+ elHist->SetMarkerSize(0.3);
+
+ TH2F *piHist=new TH2F("piHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
+ piHist->SetMarkerColor(2);
+ piHist->SetMarkerStyle(8);
+ piHist->SetMarkerSize(0.3);
+
+ TH2F *kaHist=new TH2F("kaHist","dE/dX vs momentum",100,0.,4.,100,0.,500.);
+ kaHist->SetMarkerColor(4);
+ kaHist->SetMarkerStyle(8);
+ kaHist->SetMarkerSize(0.3);
+
+ TH2F *prHist=new
+ TH2F("prHist","Classification into e, pi, K and p",100,0.,4.,100,0.,500.);
+ prHist->SetXTitle("p (GeV/c)"); prHist->SetYTitle("dE/dx (Arb. Units)");
+ prHist->SetMarkerColor(6);
+ prHist->SetMarkerStyle(8);
+ prHist->SetMarkerSize(0.3);
TFile *ef=TFile::Open("AliESDs.root");
if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
TKey *key=0;
TIter next(ef->GetListOfKeys());
+ //****** Tentative particle type "concentrations"
+ Double_t c[5]={0.05, 0., 0.85, 0.10, 0.05};
+
//******* The loop over events
while ((key=(TKey*)next())!=0) {
cerr<<"Processing event number : "<<n++<<endl;
//****** The loop over tracks
while (ntrk--) {
AliESDtrack *t=event->GetTrack(ntrk);
+
Double_t p=t->GetP();
+
if (t->GetStatus()&AliESDtrack::kTPCin) {
Double_t dedx=t->GetTPCsignal();
tpcHist->Fill(p,dedx,1);
}
- if (t->GetStatus()&AliESDtrack::kITSin) {
- Double_t dedx=t->GetITSsignal();
- itsHist->Fill(p,dedx,1);
+
+ if (t->GetStatus()&AliESDtrack::kESDpid) {
+ Double_t dedx=t->GetTPCsignal();
+ Double_t r[10]; t->GetESDpid(r);
+ Double_t rc=0.;
+ Int_t i;
+ for (i=0; i<AliESDtrack::kSPECIES; i++) rc+=(c[i]*r[i]);
+ if (rc==0.) continue;
+
+ //Here we apply Bayes' formula
+ Double_t w[10];
+ for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rc;
+
+ if (w[4]>w[3] && w[4]>w[2] && w[4]>w[0]) prHist->Fill(p,dedx,1);
+ if (w[3]>w[4] && w[3]>w[2] && w[3]>w[0]) kaHist->Fill(p,dedx,1);
+ if (w[2]>w[3] && w[2]>w[4] && w[2]>w[0]) piHist->Fill(p,dedx,1);
+ if (w[0]>w[3] && w[0]>w[2] && w[0]>w[4]) elHist->Fill(p,dedx,1);
}
+
}
delete event;
}
c1->cd(1);
tpcHist->Draw();
c1->cd(2);
- itsHist->Draw();
+ prHist->Draw();
+ kaHist->Draw("same");
+ piHist->Draw("same");
+ elHist->Draw("same");
ef->Close();
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+// Implementation of the combined PID class
+//
+// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-----------------------------------------------------------------
+
+#include "AliESDpid.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+
+ClassImp(AliESDpid)
+
+//_________________________________________________________________________
+Int_t AliESDpid::MakePID(AliESD *event)
+{
+ Int_t ntrk=event->GetNumberOfTracks();
+ for (Int_t i=0; i<ntrk; i++) {
+ Int_t ns=AliESDtrack::kSPECIES;
+ Double_t p[10]={1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
+
+ AliESDtrack *t=event->GetTrack(i);
+
+ if ((t->GetStatus()&AliESDtrack::kITSpid )!=0) {
+ Double_t d[10];
+ t->GetITSpid(d);
+ for (Int_t j=0; j<ns; j++) p[j]*=d[j];
+ }
+
+ if ((t->GetStatus()&AliESDtrack::kTPCpid )!=0) {
+ Double_t d[10];
+ t->GetTPCpid(d);
+ for (Int_t j=0; j<ns; j++) p[j]*=d[j];
+ }
+
+ t->SetESDpid(p);
+ }
+ return 0;
+}
--- /dev/null
+#ifndef ALIESDPID_H
+#define ALIESDPID_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//-------------------------------------------------------
+// Combined PID class
+//
+// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-------------------------------------------------------
+#include <Rtypes.h>
+
+class AliESD;
+
+class AliESDpid {
+public:
+ AliESDpid(){}
+ static Int_t MakePID(AliESD *event);
+private:
+ ClassDef(AliESDpid,1) // TPC PID class
+};
+
+#endif
+
+
#include "TStopwatch.h"
#include "AliESD.h"
+ #include "AliESDpid.h"
+ #include "AliTPCpidESD.h"
#include "AliTPCParam.h"
#include "AliTPCtracker.h"
#include "AliITSgeom.h"
#include "AliITStrackerV2.h"
+ #include "AliITSpidESD.h"
#endif
Int_t AliESDtest(Int_t nev=1) {
//An instance of the TPC tracker
AliTPCtracker tpcTracker(par);
+ //An instance of the TPC PID maker
+ Double_t parTPC[]={47.,0.1,3.};
+ AliTPCpidESD tpcPID(parTPC);
//File with the ITS clusters
TFile *itscf=TFile::Open("AliITSclustersV2.root");
//An instance of the ITS tracker
AliITStrackerV2 itsTracker(geom);
-
+
+ //An instance of the ITS PID maker
+ Double_t parITS[]={34.,0.12,3.};
+ AliITSpidESD itsPID(parITS);
+
TFile *ef=TFile::Open("AliESDs.root","new");
if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
rc+=itsTracker.PropagateBack(event);
itsTracker.UnloadClusters();
+ itsPID.MakePID(event);
rc+=tpcTracker.PropagateBack(event);
tpcTracker.UnloadClusters();
+ tpcPID.MakePID(event);
+
+ //Here is the combined PID
+ AliESDpid::MakePID(event);
if (rc==0) {
Char_t ename[100];
for (Int_t i=0; i<fTPCncls; i++) idx[i]=fTPCindex[i];
return fTPCncls;
}
+
+//_______________________________________________________________________
+void AliESDtrack::SetTPCpid(const Double_t *p) {
+ for (Int_t i=0; i<kSPECIES; i++) fTPCr[i]=p[i];
+ SetStatus(AliESDtrack::kTPCpid);
+}
+
+//_______________________________________________________________________
+void AliESDtrack::GetTPCpid(Double_t *p) const {
+ for (Int_t i=0; i<kSPECIES; i++) p[i]=fTPCr[i];
+}
+
+//_______________________________________________________________________
+void AliESDtrack::SetESDpid(const Double_t *p) {
+ for (Int_t i=0; i<kSPECIES; i++) fR[i]=p[i];
+ SetStatus(AliESDtrack::kESDpid);
+}
+
+//_______________________________________________________________________
+void AliESDtrack::GetESDpid(Double_t *p) const {
+ for (Int_t i=0; i<kSPECIES; i++) p[i]=fR[i];
+}
+
Bool_t UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags);
void SetIntegratedLength(Double_t l) {fTrackLength=l;}
void SetIntegratedTimes(const Double_t *times);
+ void SetESDpid(const Double_t *p);
+ void GetESDpid(Double_t *p) const;
ULong_t GetStatus() const {return fFlags;}
Int_t GetLabel() const {return fLabel;}
void GetPxPyPz(Double_t *p) const;
void GetXYZ(Double_t *r) const;
+ void SetTPCpid(const Double_t *p);
+ void GetTPCpid(Double_t *p) const;
Float_t GetTPCsignal() const {return fTPCsignal;}
Int_t GetTPCclusters(UInt_t *idx) const;
+ void SetITSpid(const Double_t *p) {;}
+ void GetITSpid(Double_t *p) const {;}
Float_t GetITSsignal() const {return fITSsignal;}
Int_t GetITSclusters(UInt_t *idx) const;
enum {
- kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSPID=0x0008,
- kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCPID=0x0080,
- kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDPID=0x0800,
- kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFPID=0x8000,
+ kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
+ kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
+ kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDpid=0x0800,
+ kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFpid=0x8000,
+ kESDpid=0x40000000,
kTIME=0x80000000
};
enum {kSPECIES=5}; // Number of particle species recognized by the PID
Int_t fTPCncls; // number of clusters assigned in the TPC
UInt_t fTPCindex[180]; //! indices of the assigned TPC clusters
Float_t fTPCsignal; // detector's PID signal
- Float_t fTPCr[kSPECIES]; //! "detector response probabilities" (for the PID)
+ Float_t fTPCr[kSPECIES]; // "detector response probabilities" (for the PID)
// TRD related track information
// TOF related track information
#pragma link C++ class AliESD+;
#pragma link C++ class AliESDtrack+;
#pragma link C++ class AliESDvertex+;
+#pragma link C++ class AliESDpid+;
#pragma link C++ class AliTrackMap+;
#pragma link C++ class AliTrackMapper+;
#pragma link C++ class AliCollisionGeometry+;
AliGausCorr.cxx AliTrackReference.cxx AliESD.cxx \
AliTrackMap.cxx AliTrackMapper.cxx AliCollisionGeometry.cxx \
AliMemoryWatcher.cxx AliBarrelTrack.cxx \
-AliESDtrack.cxx AliESDvertex.cxx
+AliESDtrack.cxx AliESDvertex.cxx AliESDpid.cxx
HDRS:= $(SRCS:.cxx=.h)
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+// Implementation of the TPC PID class
+// Very naive one... Should be made better by the detector experts...
+// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-----------------------------------------------------------------
+
+#include "AliTPCpidESD.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+
+ClassImp(AliTPCpidESD)
+
+//_________________________________________________________________________
+AliTPCpidESD::AliTPCpidESD(Double_t *param)
+{
+ //
+ // The main constructor
+ //
+ fMIP=param[0];
+ fRes=param[1];
+ fRange=param[2];
+}
+
+Double_t AliTPCpidESD::Bethe(Double_t bg) {
+ //
+ // This is the Bethe-Bloch function normalised to 1 at the minimum
+ //
+ Double_t bg2=bg*bg;
+ Double_t bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2));
+ return bethe/11.091;
+}
+
+//_________________________________________________________________________
+Int_t AliTPCpidESD::MakePID(AliESD *event)
+{
+ //
+ // This function calculates the "detector response" PID probabilities
+ //
+ static const Double_t masses[]={
+ 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
+ };
+ Int_t ntrk=event->GetNumberOfTracks();
+ for (Int_t i=0; i<ntrk; i++) {
+ AliESDtrack *t=event->GetTrack(i);
+ if ((t->GetStatus()&AliESDtrack::kTPCin )==0)
+ if ((t->GetStatus()&AliESDtrack::kTPCout)==0) continue;
+ Int_t ns=AliESDtrack::kSPECIES;
+ Double_t p[10];
+ for (Int_t j=0; j<ns; j++) {
+ Double_t mass=masses[j];
+ Double_t mom=t->GetP();
+ Double_t dedx=t->GetTPCsignal()/fMIP;
+ Double_t bethe=Bethe(mom/mass);
+ Double_t sigma=fRes*bethe;
+ if (TMath::Abs(dedx-bethe) > fRange*sigma) {
+ p[j]=0.;
+ continue;
+ }
+ p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma));
+ }
+ t->SetTPCpid(p);
+ }
+ return 0;
+}
--- /dev/null
+#ifndef ALITPCpIDESD_H
+#define ALITPCpIDESD_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//-------------------------------------------------------
+// TPC PID class
+// A very naive design... Should be made better by the detector experts...
+// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-------------------------------------------------------
+#include <Rtypes.h>
+
+class AliESD;
+
+class AliTPCpidESD {
+public:
+ AliTPCpidESD(Double_t *param);
+ Int_t MakePID(AliESD *event);
+ static Double_t Bethe(Double_t bg);
+private:
+ Double_t fMIP; // dEdx for MIP
+ Double_t fRes; // relative dEdx resolution
+ Double_t fRange; // one particle type PID range (in sigmas)
+ ClassDef(AliTPCpidESD,1) // TPC PID class
+};
+
+#endif
+
+
#pragma link C++ class AliTPCKalmanSegment+;
#pragma link C++ class AliTPCPid+;
+#pragma link C++ class AliTPCpidESD+;
#pragma link C++ class AliTPCtrackPid+;
AliTPCpolyTrack.cxx \
AliTPCBuffer.cxx AliTPCBuffer160.cxx \
AliTPCCompression.cxx AliTPCDDLRawData.cxx \
- AliTPCHuffman.cxx AliTPCPid.cxx AliTPCtrackPid.cxx
+ AliTPCHuffman.cxx AliTPCPid.cxx AliTPCtrackPid.cxx AliTPCpidESD.cxx
HDRS:= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff