#include "AliTPCtracker.h"
#include "AliITSgeom.h"
#include "AliITStrackerV2.h"
+ #include "AliTRDtracker.h"
+ #include "AliTRDPartID.h"
#include "AliITSpidESD.h"
#endif
-Int_t AliESDtest(Int_t nev=1) {
+Int_t AliESDtest(Int_t nev=1,
+ const char* fileNameITSClusters = "its.clusters.root",
+ const char* fileNameTPCClusters = "tpc.clusters.root",
+ const char* fileNameTRDClusters = "trd.clusters.root") {
+
//File with the TPC clusters
- TFile *tpccf=TFile::Open("AliTPCclusters.root");
+ TFile *tpccf=TFile::Open(fileNameTPCClusters);
if (!tpccf->IsOpen()) {
- cerr<<"Can't open AliTPCclusters.root !\n";
+ cerr<<"Can't open "<<fileNameTPCClusters<<" !\n";
return 2;
}
AliTPCParam *par=(AliTPCParam*)tpccf->Get("75x40_100x60_150x60");
AliTPCpidESD tpcPID(parTPC);
//File with the ITS clusters
- TFile *itscf=TFile::Open("AliITSclustersV2.root");
+ TFile *itscf=TFile::Open(fileNameITSClusters);
if (!itscf->IsOpen()) {
- cerr<<"Can't open AliITSclustersV2.root !\n";
+ cerr<<"Can't open "<<fileNameITSClusters<<".root !\n";
return 4;
}
AliITSgeom *geom=(AliITSgeom*)itscf->Get("AliITSgeom");
//An instance of the ITS PID maker
Double_t parITS[]={34.,0.12,3.};
AliITSpidESD itsPID(parITS);
-
- TFile *ef=TFile::Open("AliESDs.root","new");
+
+ //File with the TRD clusters
+ TFile *trdcf=TFile::Open(fileNameTRDClusters);
+ if (!trdcf->IsOpen()) {
+ cerr<<"Can't open "<<fileNameTRDClusters<<".root !\n";
+ return 6;
+ }
+
+ //An instance of the TRD tracker
+ AliTRDtracker trdTracker(trdcf);
+
+ //An instance of the TRD PID maker
+ TFile* pidFile = TFile::Open("pid.root");
+ if (!pidFile->IsOpen()) {
+ cerr << "Can't get pid.root !\n";
+ return 7;
+ }
+ AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID");
+ if (!trdPID) {
+ cerr << "Can't get PID object !\n";
+ return 8;
+ }
+
+ TFile *ef=TFile::Open("AliESDs.root","RECREATE");
if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
TStopwatch timer;
rc+=itsTracker.PropagateBack(event);
itsTracker.UnloadClusters();
+
itsPID.MakePID(event);
rc+=tpcTracker.PropagateBack(event);
tpcTracker.UnloadClusters();
+
tpcPID.MakePID(event);
+ trdTracker.SetEventNumber(i);
+ trdcf->cd();
+ trdTracker.LoadClusters();
+
+ rc+=trdTracker.PropagateBack(event);
+ trdTracker.UnloadClusters();
+
+ for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
+ AliESDtrack* track = event->GetTrack(iTrack);
+ trdPID->MakePID(track);
+ }
+
//Here is the combined PID
AliESDpid::MakePID(event);
if (rc==0) {
Char_t ename[100];
sprintf(ename,"%d",i);
+ ef->cd();
if (!event->Write(ename)) rc++;
}
if (rc) {
}
timer.Stop(); timer.Print();
+ pidFile->Close();
+ trdcf->Close();
delete geom;
itscf->Close();
delete par;
else if (mass<0.4) fR[2]=1.; // the ITS reconstruction
else fR[3]=1.;} //
break;
+ case kTRDin: case kTRDout: case kTRDrefit:
+ fTRDncls=t->GetNumberOfClusters();
+ fTRDchi2=t->GetChi2();
+ fTRDsignal=t->GetPIDsignal();
+ break;
default:
Error("UpdateTrackParams()","Wrong flag !\n");
return kFALSE;
for (Int_t i=0; i<kSPECIES; i++) p[i]=fTPCr[i];
}
+//_______________________________________________________________________
+void AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
+{
+ fTRDr[iSpecies] = p;
+}
+
+Float_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
+{
+ return fTRDr[iSpecies];
+}
+
//_______________________________________________________________________
void AliESDtrack::SetESDpid(const Double_t *p) {
for (Int_t i=0; i<kSPECIES; i++) fR[i]=p[i];
Float_t GetITSsignal() const {return fITSsignal;}
Int_t GetITSclusters(UInt_t *idx) const;
+ Float_t GetTRDsignal() const {return fTRDsignal;}
+ void SetTRDpid(Int_t iSpecies, Float_t p);
+ Float_t GetTRDpid(Int_t iSpecies) const;
+
enum {
kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
kTIME=0x80000000
};
enum {kSPECIES=5}; // Number of particle species recognized by the PID
-
+
protected:
ULong_t fFlags; // Reconstruction status flags
Int_t fLabel; // Track label
Float_t fTPCr[kSPECIES]; // "detector response probabilities" (for the PID)
// TRD related track information
+ Float_t fTRDchi2; // chi2 in the TRD
+ Int_t fTRDncls; // number of clusters assigned in the TRD
+ Float_t fTRDsignal; // detector's PID signal
+ Float_t fTRDr[kSPECIES]; //! "detector response probabilities" (for the PID)
+
// TOF related track information
// HMPID related track information
Double_t par[5];
Double_t localX = GetX();
GetExternalParameters(localX, par);
+ if (par[2] > 1.) par[2] = 1.;
+ if (par[2] < -1.) par[2] = -1.;
Double_t phi = TMath::ASin(par[2]) + GetAlpha();
while (phi < 0) phi += TMath::TwoPi();
while (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
Double_t par[5];
Double_t localX = GetX();
GetExternalParameters(localX, par);
- return 1. / TMath::Abs(par[4] * TMath::Sin(TMath::ATan(par[3])));
+ return 1. / TMath::Abs(par[4] * TMath::Cos(TMath::ATan(par[3])));
}
//_______________________________________________________________________
void AliKalmanTrack::StartTimeIntegral()
--- /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. *
+ **************************************************************************/
+
+#include "AliTRDPartID.h"
+#include "AliESDtrack.h"
+#include "TProfile.h"
+#include "TF1.h"
+
+#include <TMath.h>
+
+ClassImp(AliTRDPartID)
+
+
+AliTRDPartID::AliTRDPartID()
+{
+// create a TRD particle identifier
+
+ fBetheBloch = NULL;
+ fRes = 0.2;
+ fRange = 3.;
+}
+
+AliTRDPartID::AliTRDPartID(TF1* betheBloch, Double_t res, Double_t range)
+{
+// create a TRD particle identifier with custom settings
+
+ fBetheBloch = betheBloch;
+ fRes = res;
+ fRange = range;
+}
+
+AliTRDPartID::~AliTRDPartID()
+{
+ if (fBetheBloch) delete fBetheBloch;
+}
+
+
+Bool_t AliTRDPartID::MakePID(AliESDtrack* track)
+{
+// This function calculates the "detector response" PID probabilities
+
+ static const Double_t masses[]={
+ 0.000511, 0.105658, 0.139570, 0.493677, 0.938272
+ };
+
+ if (((track->GetStatus()&AliESDtrack::kTRDin) == 0) &&
+ ((track->GetStatus()&AliESDtrack::kTRDout) == 0)) return kFALSE;
+ Double_t momentum = track->GetP();
+ if (momentum < 0.001) return kFALSE;
+
+ // get the probability densities
+ Double_t pSum = 0;
+ for (Int_t iSpecies = 0; iSpecies < AliESDtrack::kSPECIES; iSpecies++) {
+ Double_t expectedSignal = fBetheBloch->Eval(momentum/masses[iSpecies]);
+ Double_t expectedError = fRes * expectedSignal;
+ Double_t measuredSignal = track->GetTRDsignal();
+ if (TMath::Abs(measuredSignal - expectedSignal) > fRange * expectedError) {
+ track->SetTRDpid(iSpecies, 0.);
+ } else {
+ Double_t p = TMath::Gaus(measuredSignal, expectedSignal, expectedError);
+ pSum += p;
+ track->SetTRDpid(iSpecies, p);
+ }
+ }
+
+ // "normalize" the probability densities
+ if (pSum <= 0) return kFALSE;
+ for (Int_t iSpecies = 0; iSpecies < AliESDtrack::kSPECIES; iSpecies++) {
+ track->SetTRDpid(iSpecies, track->GetTRDpid(iSpecies) / pSum);
+ }
+
+ return kTRUE;
+}
+
+
+void AliTRDPartID::FitBetheBloch(TProfile* dEdxVsBetaGamma)
+{
+// determine the parameters of the bethe bloch function
+
+ if (fBetheBloch) delete fBetheBloch;
+ fBetheBloch = new TF1("fBetheBlochTRD", fcnBetheBloch,
+ 0.001, 100000., 5);
+ fBetheBloch->SetParameters(1, 10, 0.00002, 2, 2);
+ fBetheBloch->SetParLimits(2, 0., 0.01);
+ fBetheBloch->SetParLimits(3, 0., 10.);
+ fBetheBloch->SetParLimits(4, 0., 10.);
+ fBetheBloch->SetFillStyle(0);
+ fBetheBloch->SetLineColor(kRed);
+ dEdxVsBetaGamma->Fit(fBetheBloch, "N");
+}
+
+TF1* AliTRDPartID::CreateBetheBloch(Double_t mass)
+{
+// create a function for expected dE/dx vs p
+
+ TF1* result = new TF1("betheBlochMass", fcnBetheBlochMass,
+ 0.001, 100000., 6);
+ result->SetParameter(0, mass);
+ if (fBetheBloch) {
+ for (Int_t iPar = 0; iPar < 5; iPar++) {
+ result->SetParameter(iPar+1, fBetheBloch->GetParameter(iPar));
+ result->SetParError(iPar+1, fBetheBloch->GetParError(iPar));
+ }
+ }
+ result->SetFillStyle(0);
+ return result;
+}
+
+Double_t AliTRDPartID::fcnBetheBloch(Double_t* xx, Double_t* par)
+{
+// parametrized bethe bloch function (xx[0] = beta*gamma = p/m):
+//
+// p0/beta^p3 * [ p1 - log(p2 + 1/(beta*gamma)^p4) - beta^p3 ]
+//
+
+ Double_t betaGamma2 = xx[0] * xx[0];
+ Double_t beta2 = betaGamma2 / (1. + betaGamma2);
+ Double_t betaPar3 = TMath::Power(beta2, par[3]/2.);
+ return par[0]/betaPar3 * (par[1] - TMath::Log(TMath::Abs(par[2] + TMath::Power(betaGamma2, -par[4]/2.))) - betaPar3);
+}
+
+Double_t AliTRDPartID::fcnBetheBlochMass(Double_t* xx, Double_t* par)
+{
+// parametrized bethe bloch function: xx[0] = p, p0 = mass
+
+ Double_t betaGamma = xx[0] / par[0];
+ return fcnBetheBloch(&betaGamma, &par[1]);
+}
--- /dev/null
+#ifndef ALITRDPARTID_H
+#define ALITRDPARTID_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+#include <TObject.h>
+
+class AliESDtrack;
+class TProfile;
+class TF1;
+
+
+class AliTRDPartID: public TObject {
+ public:
+ AliTRDPartID();
+ AliTRDPartID(TF1* betheBloch, Double_t res, Double_t range);
+ virtual ~AliTRDPartID();
+
+ Bool_t MakePID(AliESDtrack* track);
+
+ void FitBetheBloch(TProfile* dEdxVsBetaGamma);
+ inline TF1* GetBetheBloch() {return fBetheBloch;};
+ TF1* CreateBetheBloch(Double_t mass);
+
+ private:
+ static Double_t fcnBetheBloch(Double_t* xx, Double_t* par);
+ static Double_t fcnBetheBlochMass(Double_t* xx, Double_t* par);
+
+ TF1* fBetheBloch; // parametrized bethe bloch function
+ Double_t fRes; // relative dE/dx resolution
+ Double_t fRange; // cut off in standard deviations
+
+ ClassDef(AliTRDPartID,1) // TRD PID class
+};
+
+#endif
+
+
/*
$Log$
+Revision 1.19 2003/05/22 10:46:46 hristov
+Using access methods instead of data members
+
Revision 1.18 2003/04/10 10:36:54 hristov
Code for unified TPC/TRD tracking (S.Radomski)
#include "AliTRDcluster.h"
#include "AliTRDtrack.h"
#include "../TPC/AliTPCtrack.h"
+#include "AliESDtrack.h"
ClassImp(AliTRDtrack)
}
}
//_____________________________________________________________________________
+AliTRDtrack::AliTRDtrack(const AliESDtrack& t)
+ :AliKalmanTrack() {
+ //
+ // Constructor from AliESDtrack
+ //
+
+ SetLabel(t.GetLabel());
+ SetChi2(0.);
+ SetMass(t.GetMass());
+ SetNumberOfClusters(0);
+ // WARNING: cluster indices are NOT copied !!!
+
+ fdEdx=0;
+
+ fLhElectron = 0.0;
+ fNWrong = 0;
+ fNRotate = 0;
+
+ fAlpha = t.GetAlpha();
+ if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
+ else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
+
+ Double_t x, p[5]; t.GetExternalParameters(x,p);
+
+ fX=x;
+
+ x = GetConvConst();
+
+ fY=p[0];
+ fZ=p[1];
+ fT=p[3];
+ fC=p[4]/x;
+ fE=fC*fX - p[2];
+
+ //Conversion of the covariance matrix
+ Double_t c[15]; t.GetExternalCovariance(c);
+
+ c[10]/=x; c[11]/=x; c[12]/=x; c[13]/=x; c[14]/=x*x;
+
+ Double_t c22=fX*fX*c[14] - 2*fX*c[12] + c[5];
+ Double_t c32=fX*c[13] - c[8];
+ Double_t c20=fX*c[10] - c[3], c21=fX*c[11] - c[4], c42=fX*c[14] - c[12];
+
+ fCyy=c[0 ];
+ fCzy=c[1 ]; fCzz=c[2 ];
+ fCey=c20; fCez=c21; fCee=c22;
+ fCty=c[6 ]; fCtz=c[7 ]; fCte=c32; fCtt=c[9 ];
+ fCcy=c[10]; fCcz=c[11]; fCce=c42; fCct=c[13]; fCcc=c[14];
+
+ // Initialization [SR, GSI, 18.02.2003]
+ for(Int_t i=0; i<kMAX_CLUSTERS_PER_TRACK; i++) {
+ fdQdl[i] = 0;
+ fIndex[i] = 0;
+ }
+}
+//_____________________________________________________________________________
void AliTRDtrack::GetBarrelTrack(AliBarrelTrack *track) {
//
class AliTRDcluster;
class AliTPCtrack;
+class AliESDtrack;
const unsigned kMAX_CLUSTERS_PER_TRACK=210;
const Double_t cc[15], Double_t xr, Double_t alpha);
AliTRDtrack(const AliTRDtrack& t);
AliTRDtrack(const AliKalmanTrack& t, Double_t alpha);
+ AliTRDtrack(const AliESDtrack& t);
Int_t Compare(const TObject *o) const;
void CookdEdx(Double_t low=0.05, Double_t up=0.70);
void GetCovariance(Double_t cov[15]) const;
Double_t GetdEdx() const {return fdEdx;}
+ Double_t GetPIDsignal() const {return GetdEdx();}
Double_t GetEta() const {return fE;}
void GetExternalCovariance(Double_t cov[15]) const ;
/*
$Log$
+Revision 1.26 2003/04/10 10:36:54 hristov
+Code for unified TPC/TRD tracking (S.Radomski)
+
Revision 1.25 2003/03/19 17:14:11 hristov
Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
#include "AliTRDgeometryDetail.h"
#include "AliTRDcluster.h"
#include "AliTRDtrack.h"
+#include "AliTRDPartID.h"
#include "../TPC/AliTPCtrack.h"
#include "AliTRDtracker.h"
}
else {
in->cd();
- in->ls();
+// in->ls();
fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
fPar = (AliTRDparameter*) in->Get("TRDparameter");
- fGeom->Dump();
+// fGeom->Dump();
}
if(fGeom) {
}
+//_____________________________________________________________________________
+Int_t AliTRDtracker::PropagateBack(AliESD* event) {
+ //
+ // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
+ // backpropagated by the TPC tracker. Each seed is first propagated
+ // to the TRD, and then its prolongation is searched in the TRD.
+ // If sufficiently long continuation of the track is found in the TRD
+ // the track is updated, otherwise it's stored as originaly defined
+ // by the TPC tracker.
+ //
+
+ Int_t found=0;
+ Float_t foundMin = 40;
+
+ Int_t n = event->GetNumberOfTracks();
+ for (Int_t i=0; i<n; i++) {
+ AliESDtrack* seed=event->GetTrack(i);
+ ULong_t status=seed->GetStatus();
+ if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
+ if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
+
+ Int_t lbl = seed->GetLabel();
+ AliTRDtrack *track = new AliTRDtrack(*seed);
+ track->SetSeedLabel(lbl);
+ fNseeds++;
+
+ Int_t expectedClr = FollowBackProlongation(*track);
+
+ Int_t foundClr = track->GetNumberOfClusters();
+ if (foundClr >= foundMin) {
+ if(foundClr >= 2) {
+ track->CookdEdx();
+// CookLabel(track, 1-fLabelFraction);
+ UseClusters(track);
+ }
+
+ // Propagate to outer reference plane [SR, GSI, 18.02.2003]
+// track->PropagateTo(364.8); why?
+
+ seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
+ found++;
+ }
+
+ }
+
+ cerr<<"Number of seeds: "<<fNseeds<<endl;
+ cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
+
+ return 0;
+
+}
+
//---------------------------------------------------------------------------
Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
if (fTrSec[s]->GetLayer(nr)->IsSensitive() !=
fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) {
- if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack);
+// if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack);
}
if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() &&
#include "AliTracker.h"
#include "TObjArray.h"
#include "AliBarrelTrack.h"
+#include "AliESD.h"
class TFile;
+class TTree;
class TParticle;
class TParticlePDG;
Int_t Clusters2Tracks(const TFile *in, TFile *out);
Int_t PropagateBack(const TFile *in, TFile *out);
+ Int_t PropagateBack(AliESD* event);
//Int_t Refit(const TFile *in, TFile *out);
Int_t LoadClusters() {LoadEvent(); return 0;};
void UnloadClusters() {UnloadEvent();};
- AliCluster *GetCluster(Int_t index) const { return (AliCluster*) fClusters->UncheckedAt(index); };
+ AliCluster *GetCluster(Int_t index) const { if (index >= fNclusters) return NULL; return (AliCluster*) fClusters->UncheckedAt(index); };
virtual void CookLabel(AliKalmanTrack *t,Float_t wrong) const;
virtual void UseClusters(const AliKalmanTrack *t, Int_t from=0) const;
#pragma link C++ class AliTRDarrayI+;
#pragma link C++ class AliTRDarrayF+;
#pragma link C++ class AliTRDparameter+;
+#pragma link C++ class AliTRDPartID+;
#endif
AliTRDsimple.cxx \
AliTRDsimpleMC.cxx \
AliTRDsimpleGen.cxx \
- AliTRDparameter.cxx
+ AliTRDparameter.cxx \
+ AliTRDPartID.cxx
HDRS= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff