Update of TRD code
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 7 May 2001 08:08:06 +0000 (08:08 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 7 May 2001 08:08:06 +0000 (08:08 +0000)
18 files changed:
TRD/AliTRDcluster.cxx
TRD/AliTRDcluster.h
TRD/AliTRDdataArray.h
TRD/AliTRDdataArrayF.cxx
TRD/AliTRDdataArrayF.h
TRD/AliTRDdataArrayI.cxx
TRD/AliTRDdataArrayI.h
TRD/AliTRDgeometry.cxx
TRD/AliTRDgeometry.h
TRD/AliTRDgeometryFull.h
TRD/AliTRDhit.h
TRD/AliTRDpid.cxx [new file with mode: 0644]
TRD/AliTRDpid.h [new file with mode: 0644]
TRD/AliTRDsim.h
TRD/AliTRDtrack.cxx
TRD/AliTRDtrack.h
TRD/AliTRDtracker.cxx
TRD/AliTRDtracker.h

index 25f0008c0310c0a4c9482f56bde43e2f721f95f6..b70e5b6b8301f845bc1d27bbf91790695a776ea0 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.3  2000/12/08 16:07:02  cblume
+Update of the tracking by Sergei
+
 Revision 1.2  2000/10/06 16:49:46  cblume
 Made Getters const
 
@@ -27,38 +30,48 @@ Add the tracking code
 #include "AliTRDgeometry.h"
 #include "AliTRDrecPoint.h"
 
-
 ClassImp(AliTRDcluster)
  
 //_____________________________________________________________________________
-AliTRDcluster::AliTRDcluster() {
-  //default constructor
+AliTRDcluster::AliTRDcluster() 
+{
+  //
+  // Default constructor
+  //
+
+  fDetector  = 0;
+  fTimeBin   = 0;
+  fTracks[0] = 0;
+  fTracks[1] = 0;
+  fTracks[2] = 0; 
+  fY         = 0;
+  fZ         = 0;
+  fQ         = 0;
+  fSigmaY2   = 0;
+  fSigmaZ2   = 0;
 
-  fDetector = fTimeBin = 0;
-  fTracks[0]=fTracks[1]=fTracks[2]=0; 
-  fY=fZ=fQ=fSigmaY2=fSigmaZ2=0.;
 }
 
 //_____________________________________________________________________________
-AliTRDcluster::AliTRDcluster(AliTRDrecPoint *rp)
+AliTRDcluster::AliTRDcluster(const AliTRDrecPoint &p)
 {
   //
-  // constructor from AliTRDrecPoint
+  // Constructor from AliTRDrecPoint
   //
 
-  fDetector   = rp->GetDetector();
-  fTimeBin    = rp->GetLocalTimeBin();
+  fDetector   = p.GetDetector();
+  fTimeBin    = p.GetLocalTimeBin();
 
-  fTracks[0]  = rp->GetTrackIndex(0);
-  fTracks[1]  = rp->GetTrackIndex(1);
-  fTracks[2]  = rp->GetTrackIndex(2);
+  fTracks[0]  = p.GetTrackIndex(0);
+  fTracks[1]  = p.GetTrackIndex(1);
+  fTracks[2]  = p.GetTrackIndex(2);
 
-  fQ          = rp->GetEnergy();
+  fQ          = p.GetEnergy();
 
-  fY          = rp->GetY();
-  fZ          = rp->GetZ();
-  fSigmaY2    = rp->GetSigmaY2();
-  fSigmaZ2    = rp->GetSigmaZ2();  
+  fY          = p.GetY();
+  fZ          = p.GetZ();
+  fSigmaY2    = p.GetSigmaY2();
+  fSigmaZ2    = p.GetSigmaZ2();  
 
   fSigmaY2    = 0.2;
   fSigmaZ2    = 5.;  
@@ -66,25 +79,126 @@ AliTRDcluster::AliTRDcluster(AliTRDrecPoint *rp)
 }
 
 //_____________________________________________________________________________
-AliTRDcluster::AliTRDcluster(AliTRDcluster *cl)
+AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
 {
   //
   // Copy constructor 
   //
 
-  fDetector   = cl->GetDetector();
-  fTimeBin    = cl->GetLocalTimeBin();
+  ((AliTRDcluster &) c).Copy(*this);
+
+}
+
+//_____________________________________________________________________________
+AliTRDcluster::~AliTRDcluster()
+{
+  //
+  // AliTRDcluster destructor
+  //
+
+}
+
+//_____________________________________________________________________________
+AliTRDcluster &AliTRDcluster::operator=(const AliTRDcluster &c)
+{
+  //
+  // Assignment operator
+  //
+
+  if (this != &c) ((AliTRDcluster &) c).Copy(*this);
+  return *this;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDcluster::Copy(TObject &c)
+{
+  //
+  // Copy function
+  //
+
+  ((AliTRDcluster &) c).fDetector   = fDetector;
+  ((AliTRDcluster &) c).fTimeBin    = fTimeBin;
 
-  fTracks[0]  = cl->GetTrackIndex(0);
-  fTracks[1]  = cl->GetTrackIndex(1);
-  fTracks[2]  = cl->GetTrackIndex(2);
+  ((AliTRDcluster &) c).fTracks[0]  = fTracks[0];
+  ((AliTRDcluster &) c).fTracks[1]  = fTracks[1];
+  ((AliTRDcluster &) c).fTracks[2]  = fTracks[2];
 
-  fQ          = cl->GetQ();
+  ((AliTRDcluster &) c).fQ          = fQ;
+
+  ((AliTRDcluster &) c).fY          = fY;
+  ((AliTRDcluster &) c).fZ          = fZ;
+  ((AliTRDcluster &) c).fSigmaY2    = fSigmaY2;
+  ((AliTRDcluster &) c).fSigmaZ2    = fSigmaZ2;  
+
+}
+
+
+//_____________________________________________________________________________
+void AliTRDcluster::AddTrackIndex(Int_t *track)
+{
+  //
+  // Adds track index. Currently assumed that track is an array of
+  // size 9, and up to 3 track indexes are stored in fTracks[3].
+  // Indexes are sorted according to:
+  //  1) index of max number of appearances is stored first
+  //  2) if two or more indexes appear equal number of times, the lowest
+  //     ones are stored first;
+  //
 
-  fY          = cl->GetY();
-  fZ          = cl->GetZ();
-  fSigmaY2    = cl->GetSigmaY2();
-  fSigmaZ2    = cl->GetSigmaZ2();  
+  const Int_t size = 9;
+
+  Int_t entries[size][2], i, j, index;
+
+  Bool_t index_added;
+
+  for (i=0; i<size; i++) {
+    entries[i][0]=-1;
+    entries[i][1]=0;
+  }
+
+  for (Int_t k=0; k<size; k++) {
+    index=track[k];
+    index_added=kFALSE; j=0;
+    if (index >= 0) {
+      while ( (!index_added) && ( j < size ) ) {
+        if ((entries[j][0]==index) || (entries[j][1]==0)) {
+          entries[j][0]=index;
+          entries[j][1]=entries[j][1]+1;
+          index_added=kTRUE;
+        }
+        j++;
+      }
+    }
+  }
+
+  // sort by number of appearances and index value
+  Int_t swap=1, tmp0, tmp1;
+  while ( swap > 0) {
+    swap=0;
+    for(i=0; i<(size-1); i++) {
+      if ((entries[i][0] >= 0) && (entries[i+1][0] >= 0)) {
+        if ((entries[i][1] < entries[i+1][1]) ||
+            ((entries[i][1] == entries[i+1][1]) &&
+             (entries[i][0] > entries[i+1][0]))) {
+               tmp0=entries[i][0];
+               tmp1=entries[i][1];
+               entries[i][0]=entries[i+1][0];
+               entries[i][1]=entries[i+1][1];
+               entries[i+1][0]=tmp0;
+               entries[i+1][1]=tmp1;
+               swap++;
+        }
+      }
+    }
+  }
+
+  // set track indexes
+  for(i=0; i<3; i++) {
+    fTracks[i] = entries[i][0];
+  }
+
+  return;
 
 }
 
index 5951d52e0d7c090a4c482e76603bbc849ca1f482..032f073b5d40022f7ac6c5aab480a97504ef1f38 100644 (file)
@@ -14,40 +14,57 @@ class AliTRDcluster : public TObject {
  public:
 
   AliTRDcluster();
-  AliTRDcluster(AliTRDrecPoint *rp);
-  AliTRDcluster(AliTRDcluster *cl);
-  
-  Int_t   GetDetector() const           { return fDetector; };
-  Int_t   GetLocalTimeBin() const       { return fTimeBin; }
-
-  Float_t GetSigmaY2() const            { return fSigmaY2; }
-  Float_t GetSigmaZ2() const            { return fSigmaZ2; }
-  Float_t GetY() const                  { return fY; }
-  Float_t GetZ() const                  { return fZ; }
-  Float_t GetQ() const                  { return fQ; }
-
-  Int_t   IsUsed() const                { return (fQ<0) ? 1 : 0; }
-  void    Use()                         { fQ=-fQ; }
-  Int_t   GetTrackIndex(Int_t i) const  { return fTracks[i]; }
-
-  void    SetSigmaY2(Float_t s)         { fSigmaY2 = s; }
-  void    SetSigmaZ2(Float_t s)         { fSigmaZ2 = s; }
-
+  AliTRDcluster(const AliTRDcluster &c);
+  AliTRDcluster(const AliTRDrecPoint &p);
+  virtual ~AliTRDcluster();
+  AliTRDcluster &operator=(const AliTRDcluster &c);
+
+  virtual void    Copy(TObject &c);
+  virtual void    AddTrackIndex(Int_t *i);
+   
+          Int_t   GetDetector() const             { return fDetector; }
+          Int_t   GetLocalTimeBin() const         { return fTimeBin;  }
+
+          Float_t GetSigmaY2() const              { return fSigmaY2; }
+          Float_t GetSigmaZ2() const              { return fSigmaZ2; }
+          Float_t GetY() const                    { return fY; }
+          Float_t GetZ() const                    { return fZ; }
+          Float_t GetQ() const                    { return fQ; }
+
+          Int_t   IsUsed() const                  { return (fQ < 0) ? 1 : 0; }
+          void    Use()                           { fQ = -fQ; }
+          Int_t   GetTrackIndex(Int_t i) const    { return fTracks[i]; }
+
+          Bool_t  FromUnfolding() const           { return TestBit(kUnfold); }
+
+          void    SetDetector(Int_t d)            { fDetector  = d; }
+          void    SetLocalTimeBin(Int_t t)        { fTimeBin   = t; }
+          void    SetQ(Float_t q)                 { fQ         = q; }
+          void    SetY(Float_t y)                 { fY         = y; }
+          void    SetZ(Float_t z)                 { fZ         = z; }
+          void    SetTrackIndex(Int_t i, Int_t t) { fTracks[i] = t; } 
+          void    SetSigmaY2(Float_t s)           { fSigmaY2   = s; }
+          void    SetSigmaZ2(Float_t s)           { fSigmaZ2   = s; }
+
+          void    SetUnfolding()                  { SetBit(kUnfold); }
 
  protected:
 
-  Int_t    fDetector;    // TRD detector number
-  Int_t    fTimeBin;     // Time bin number within the detector
+  enum {
+    kUnfold = 0x00000001    // Cluster results from unfolding procedure
+  };
 
-  Int_t    fTracks[3];   // labels of overlapped tracks
-  Float_t  fQ;           // amplitude 
-  Float_t  fY;           // local Rphi coordinate (cm) within tracking sector
-  Float_t  fZ;           // local Z coordinate (cm) within tracking sector
-  Float_t  fSigmaY2;     // Y variance (cm)
-  Float_t  fSigmaZ2;     // Z variance (cm)  
+  Int_t    fDetector;       // TRD detector number
+  Int_t    fTimeBin;        // Time bin number within the detector
 
+  Int_t    fTracks[3];      // labels of overlapped tracks
+  Float_t  fQ;              // amplitude 
+  Float_t  fY;              // local Rphi coordinate (cm) within tracking sector
+  Float_t  fZ;              // local Z coordinate (cm) within tracking sector
+  Float_t  fSigmaY2;        // Y variance (cm)
+  Float_t  fSigmaZ2;        // Z variance (cm)  
 
-  ClassDef(AliTRDcluster,1) // Reconstructed point for the TRD
+  ClassDef(AliTRDcluster,1) // Cluster for the TRD
  
 };
 
index 28c1dfb9fba5d5f44cf24dfdba0d5b150355f99e..fbee032ae5838d70ed51784b2ee3bae618a072d7 100644 (file)
@@ -33,12 +33,16 @@ class AliTRDdataArray : public AliTRDsegmentID {
   virtual Int_t  GetNCol() const               { return fNcol;       };
   virtual Int_t  GetNtime() const              { return fNtime;      };
           Int_t  GetIndex(Int_t row, Int_t col, Int_t time) const;
+          Int_t  GetIndexUnchecked(Int_t row, Int_t col, Int_t time) const
+           { return time * fNrow*fNcol + GetIdx1Unchecked(row,col); };
           Int_t  GetBufType() const            { return fBufType;    };
   virtual Int_t  GetNelems() const             { return fNelems;     };
 
  protected:
 
           Int_t  GetIdx1(Int_t row, Int_t col) const;
+  inline  Int_t  GetIdx1Unchecked(Int_t row, Int_t col) const
+                                               { return row + col * fNrow; };
   inline  Bool_t CheckBounds(const char *where, Int_t idx1, Int_t idx2);
   inline  Bool_t OutOfBoundsError(const char *where, Int_t idx1, Int_t idx2);
  
index 43c9e75afdd30c1b1756b8183e68d548bf10b129..45638635668278ef63ddf3989ff01c9583637950 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.8  2000/11/23 14:34:08  cblume
+Fixed bug in expansion routine of arrays (initialize buffers properly)
+
 Revision 1.7  2000/11/20 08:56:07  cblume
 Cleanup of data arrays
 
@@ -233,6 +236,18 @@ Float_t AliTRDdataArrayF::GetData(Int_t row, Int_t col, Int_t time) const
 
 }
 
+//_____________________________________________________________________________
+Float_t AliTRDdataArrayF::GetDataFast(Int_t idx1, Int_t idx2) const
+{
+  //
+  // Returns the data value at a given position of the array
+  // No boundary checking
+  //
+
+  return fElements->At(fIndex->At(idx2)+idx1);
+
+}
+
 //_____________________________________________________________________________
 void AliTRDdataArrayF::Compress(Int_t bufferType, Float_t threshold)
 {
@@ -639,17 +654,6 @@ Float_t AliTRDdataArrayF::GetData1(Int_t idx1, Int_t idx2) const
 
 }
 
-//____________________________________________________________________________
-Float_t AliTRDdataArrayF::GetDataFast(Int_t idx1, Int_t idx2) const
-{
-  //
-  // Returns the value at a given position in the array
-  //
-
-  return fElements->At(fIndex->At(idx2) + idx1); 
-
-}
-
 //_____________________________________________________________________________
 void AliTRDdataArrayF::SetData(Int_t row, Int_t col, Int_t time, Float_t value)
 {
@@ -675,20 +679,14 @@ void AliTRDdataArrayF::SetData(Int_t row, Int_t col, Int_t time, Float_t value)
 }
 
 //_____________________________________________________________________________
-void  AliTRDdataArrayF::SetDataFast(Int_t idx1, Int_t idx2, Float_t value)
+void AliTRDdataArrayF::SetDataFast(Int_t idx1, Int_t idx2, Float_t value)
 {
   //
-  // Set the value at a given position in the array
+  // Sets the data value at a given position of the array
+  // No boundary checking
   //
 
-  if ((idx1 < 0) || (idx1 >= fNdim1) || 
-      (idx2 < 0) || (idx2 >= fNdim2)) { 
-    TObject::Error("SetDataFast"
-                  ,"idx1 %d  idx2 %d out of bounds (size: %d x %d, this: 0x%08x)"
-                  ,idx1,idx2,fNdim1,fNdim2,this);
-  }
-
-  (*fElements)[fIndex->fArray[idx2] + idx1] = value; 
+  (*fElements)[fIndex->fArray[idx2]+idx1] = value;
 
 }
 
index d8fcd799e150d73663f5bb7601289c49a58783b8..221083b3b48570c075d7c1692e9e40967c8545ee 100644 (file)
@@ -40,9 +40,15 @@ class AliTRDdataArrayF : public AliTRDdataArray {
   virtual void    Reset();
 
           void    SetData(Int_t row, Int_t col, Int_t time, Float_t value);
+          void    SetDataUnchecked(Int_t row, Int_t col, Int_t time, Float_t value)
+                                  { SetDataFast(GetIdx1Unchecked(row,col),time,value); };
+
   virtual void    SetThreshold(Float_t threshold) { fThreshold = threshold; };
 
   virtual Float_t GetData(Int_t row, Int_t col, Int_t time) const;
+          Float_t GetDataUnchecked(Int_t row, Int_t col, Int_t time) const
+                                  { return GetDataFast(GetIdx1Unchecked(row,col),time); };
+                                                  
   virtual Float_t GetThreshold() const            { return fThreshold;  };
 
   virtual Int_t   GetSize();
@@ -51,8 +57,8 @@ class AliTRDdataArrayF : public AliTRDdataArray {
 
  protected:
 
-  inline  void    SetDataFast(Int_t idx1, Int_t idx2, Float_t value); 
-  inline  Float_t GetDataFast(Int_t idx1, Int_t idx2) const; 
+  inline  void    SetDataFast(Int_t idx1, Int_t idx2, Float_t v);  
+  inline  Float_t GetDataFast(Int_t idx1, Int_t idx2) const;
 
   Float_t         GetData1(Int_t idx1, Int_t idx2) const; 
   void            Expand1(); 
index 1ed83e5bacbaae7d06008ca5667424ddf9284a21..1f97dc48d77eebdd3e0e0887f7ec13994f1b4f56 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.8  2000/11/23 14:34:08  cblume
+Fixed bug in expansion routine of arrays (initialize buffers properly)
+
 Revision 1.7  2000/11/20 08:56:07  cblume
 Cleanup of data arrays
 
@@ -232,6 +235,18 @@ Int_t AliTRDdataArrayI::GetData(Int_t row, Int_t col, Int_t time) const
 
 }
 
+//_____________________________________________________________________________
+Int_t AliTRDdataArrayI::GetDataFast(Int_t idx1, Int_t idx2) const
+{
+  //
+  // Returns the data value at a given position of the array
+  // No boundary checking
+  //
+
+  return fElements->At(fIndex->At(idx2)+idx1);
+
+}
+
 //_____________________________________________________________________________
 void AliTRDdataArrayI::Compress(Int_t bufferType, Int_t threshold)
 {
@@ -634,17 +649,6 @@ Int_t AliTRDdataArrayI::GetData1(Int_t idx1, Int_t idx2) const
 
 }
 
-//_____________________________________________________________________________
-Int_t AliTRDdataArrayI::GetDataFast(Int_t idx1, Int_t idx2) const
-{
-  //
-  // Returns the value at a given position in the array
-  //
-  
-  return fElements->At(fIndex->At(idx2) + idx1); 
-
-}
-
 //_____________________________________________________________________________
 void AliTRDdataArrayI::SetData(Int_t row, Int_t col, Int_t time, Int_t value)
 {
@@ -670,20 +674,14 @@ void AliTRDdataArrayI::SetData(Int_t row, Int_t col, Int_t time, Int_t value)
 }
 
 //_____________________________________________________________________________
-void  AliTRDdataArrayI::SetDataFast(Int_t idx1, Int_t idx2, Int_t value)
+void AliTRDdataArrayI::SetDataFast(Int_t idx1, Int_t idx2, Int_t value)
 {
   //
-  // Set the value at a given position in the array
+  // Sets the data value at a given position of the array
+  // No boundary checking
   //
 
-  if ((idx1 < 0) || (idx1 >= fNdim1) || 
-      (idx2 < 0) || (idx2 >= fNdim2)) { 
-    TObject::Error("SetDataFast"
-                  ,"idx1 %d  idx2 %d out of bounds (size: %d x %d, this: 0x%08x)"
-                  ,idx1,idx2,fNdim1,fNdim2,this);
-  }
-
-  (*fElements)[fIndex->fArray[idx2] + idx1] = value; 
+  (*fElements)[fIndex->fArray[idx2]+idx1] = value;
 
 }
 
index 7f1008ba37906fb01fbe38576830447896ac1f1b..789f645e9462af1a71976573b75d85dea1430009 100644 (file)
@@ -36,9 +36,15 @@ class AliTRDdataArrayI : public AliTRDdataArray {
   virtual void   Reset();
 
           void   SetData(Int_t row, Int_t col, Int_t time, Int_t value);
+          void   SetDataUnchecked(Int_t row, Int_t col, Int_t time, Int_t value)
+                                 { SetDataFast(GetIdx1Unchecked(row,col),time,value); };
+
   virtual void   SetThreshold(Int_t threshold) { fThreshold = threshold; };
 
   virtual Int_t  GetData(Int_t row, Int_t col, Int_t time) const;
+          Int_t  GetDataUnchecked(Int_t row, Int_t col, Int_t time) const
+                                 { return GetDataFast(GetIdx1Unchecked(row,col),time); };
+
   virtual Int_t  GetThreshold() const          { return fThreshold;  };
 
   virtual Int_t  GetSize();
@@ -47,8 +53,8 @@ class AliTRDdataArrayI : public AliTRDdataArray {
 
  protected:
 
-  inline  void   SetDataFast(Int_t idx1, Int_t idx2, Int_t value); 
-  inline  Int_t  GetDataFast(Int_t idx1, Int_t idx2) const; 
+  inline  void   SetDataFast(Int_t idx1, Int_t idx2, Int_t value);
+  inline  Int_t  GetDataFast(Int_t idx1, Int_t idx2) const;
 
   Int_t          GetData1(Int_t idx1, Int_t idx2) const; 
   void           Expand1(); 
index 737a78c337d35277d629f29e9e86d582f0145948..574d242763b0f2bf98a70a8060e46191224ac61b 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.9  2001/03/27 12:48:33  cblume
+Correct for volume overlaps
+
 Revision 1.8  2001/03/13 09:30:35  cblume
 Update of digitization. Moved digit branch definition to AliTRD
 
@@ -210,15 +213,20 @@ void AliTRDgeometry::Init()
   //     +----------------------------+      +------>
   //                                             z
   //                                             
-  // IMPORTANT: time bin 0 is now the one closest to the readout !!!
+  // IMPORTANT: time bin 0 is now the first one in the drift region 
+  // closest to the readout !!!
   //
 
   // The pad column (rphi-direction)  
   SetNColPad(96);
 
-  // The time bucket. Default is 100 ns timbin size
+  // The number of time bins. Default is 100 ns timbin size
   SetNTimeBin(15);
 
+  // Additional time bins before and after the drift region.
+  // Default is to only sample the drift region
+  SetExpandTimeBin(0,0);
+
   // The rotation matrix elements
   Float_t phi = 0;
   for (isect = 0; isect < fgkNsect; isect++) {
@@ -237,7 +245,7 @@ void AliTRDgeometry::Init()
 }
 
 //_____________________________________________________________________________
-void AliTRDgeometry::SetNColPad(Int_t npad)
+void AliTRDgeometry::SetNColPad(const Int_t npad)
 {
   //
   // Redefines the number of pads in column direction
@@ -252,10 +260,12 @@ void AliTRDgeometry::SetNColPad(Int_t npad)
 }
 
 //_____________________________________________________________________________
-void AliTRDgeometry::SetNTimeBin(Int_t nbin)
+void AliTRDgeometry::SetNTimeBin(const Int_t nbin)
 {
   //
-  // Redefines the number of time bins
+  // Redefines the number of time bins in the drift region.
+  // The time bin width is defined by the length of the
+  // drift region divided by <nbin>.
   //
 
   fTimeMax     = nbin;
@@ -476,9 +486,9 @@ Bool_t AliTRDgeometry::Local2Global(Int_t iplan, Int_t icham, Int_t isect
   Float_t  rot[3];
 
   // calculate (x,y,z) position in rotated chamber
-  rot[0] = time0 - timeSlice * fTimeBinSize;
-  rot[1] = col0  + padCol    * fColPadSize[iplan];
-  rot[2] = row0  + padRow    * fRowPadSize[iplan][icham][isect];
+  rot[0] = time0 - (timeSlice - fTimeBefore) * fTimeBinSize;
+  rot[1] = col0  + padCol                    * fColPadSize[iplan];
+  rot[2] = row0  + padRow                    * fRowPadSize[iplan][icham][isect];
 
   // Rotate back to original position
   return RotateBack(idet,rot,global);
index ffa844b266844226a8dfde8f0d20288edcefacfa..9fb9f95024a1194ba205388ede575f72cd5294dc 100644 (file)
@@ -3,7 +3,11 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-/* $Id$ */
+/* $Id: AliTRDgeometry.h,v 1.8 2001/02/14 18:22:26 cbl
+ public:
+
+  enum { kNplan = 6, kNcham = 5, kNsect = 18, kNdet = 540 };
+ume Exp $ */
 
 #include "AliGeometry.h"
 
@@ -38,43 +42,58 @@ class AliTRDgeometry : public AliGeometry {
                              / fgkSheight * (fgkCheight + fgkCspace); };
   static  Float_t  Cheight() { return fgkCheight; };
   static  Float_t  Cspace()  { return fgkCspace;  };
+  static  Float_t  Ccframe() { return fgkCcframe; };
+  static  Float_t  SeThick() { return fgkSeThick; };
   static  Float_t  MyThick() { return fgkMyThick; };
   static  Float_t  DrThick() { return fgkDrThick; };
+  static  Float_t  AmThick() { return fgkAmThick; };
   static  Float_t  RaThick() { return fgkRaThick; };
+  static  Float_t  DrZpos()  { return fgkDrZpos;  };
 
   virtual void     SetPHOShole() = 0;
   virtual void     SetRICHhole() = 0;
 
-  virtual void     SetNRowPad(Int_t p, Int_t c, Int_t npad) {};
-  virtual void     SetNColPad(Int_t npad);
-  virtual void     SetNTimeBin(Int_t nbin);
+  virtual void     SetNRowPad(const Int_t p, const Int_t c, const Int_t npad) {};
+  virtual void     SetNColPad(const Int_t npad);
+  virtual void     SetNTimeBin(const Int_t nbin);
+  virtual void     SetExpandTimeBin(const Int_t nbefore, const Int_t nafter)
+                                                                  { fTimeBefore = nbefore;
+                                                                    fTimeAfter  = nafter; };
 
   virtual Bool_t   GetPHOShole() const = 0;
   virtual Bool_t   GetRICHhole() const = 0;
 
-  virtual Int_t    GetDetector(Int_t p, Int_t c, Int_t s) const;
-  virtual Int_t    GetPlane(Int_t d)   const;
-  virtual Int_t    GetChamber(Int_t d) const;
-  virtual Int_t    GetSector(Int_t d)  const;
+  virtual Int_t    GetDetector(const Int_t p, const Int_t c, const Int_t s) const;
+  virtual Int_t    GetPlane(const Int_t d)   const;
+  virtual Int_t    GetChamber(const Int_t d) const;
+  virtual Int_t    GetSector(const Int_t d)  const;
 
-  virtual Float_t  GetChamberWidth(Int_t p)                 const { return fCwidth[p]; };
+          Float_t  GetChamberWidth(const Int_t p)           const { return fCwidth[p]; };
    
-  virtual Int_t    GetRowMax(Int_t p, Int_t c, Int_t s)     const { return fRowMax[p][c][s]; };
-  virtual Int_t    GetColMax(Int_t p)                       const { return fColMax[p];       };
-  virtual Int_t    GetTimeMax()                             const { return fTimeMax;         };
-  virtual Float_t  GetRow0(Int_t p, Int_t c, Int_t s)       const { return fRow0[p][c][s]; };
-  virtual Float_t  GetCol0(Int_t p)                         const { return fCol0[p];       };
-  virtual Float_t  GetTime0(Int_t p)                        const { return fTime0[p];      };
-
-  virtual Float_t  GetRowPadSize(Int_t p, Int_t c, Int_t s) const { return fRowPadSize[p][c][s]; };
-  virtual Float_t  GetColPadSize(Int_t p)                   const { return fColPadSize[p];       };
-  virtual Float_t  GetTimeBinSize()                         const { return fTimeBinSize;         };
+          Int_t    GetRowMax(const Int_t p, const Int_t c, const Int_t s)     
+                                                            const { return fRowMax[p][c][s]; };
+          Int_t    GetColMax(const Int_t p)                 const { return fColMax[p];       };
+          Int_t    GetTimeMax()                             const { return fTimeMax;         };
+          Int_t    GetTimeBefore()                          const { return fTimeBefore;      }; 
+          Int_t    GetTimeAfter()                           const { return fTimeAfter;       }; 
+          Int_t    GetTimeTotal()                           const { return fTimeMax 
+                                                                         + fTimeBefore 
+                                                                         + fTimeAfter; };
+
+          Float_t  GetRow0(const Int_t p, const Int_t c, const Int_t s)       
+                                                            const { return fRow0[p][c][s]; };
+          Float_t  GetCol0(const Int_t p)                   const { return fCol0[p];       };
+          Float_t  GetTime0(const Int_t p)                  const { return fTime0[p];      };
+
+          Float_t  GetRowPadSize(const Int_t p, const Int_t c, const Int_t s) 
+                                                            const { return fRowPadSize[p][c][s]; };
+          Float_t  GetColPadSize(const Int_t p)             const { return fColPadSize[p];       };
+          Float_t  GetTimeBinSize()                         const { return fTimeBinSize;         };
 
   virtual void     GetGlobal(const AliRecPoint *p, TVector3 &pos, TMatrix &mat) const; 
   virtual void     GetGlobal(const AliRecPoint *p, TVector3 &pos) const;   
 
-  static  Double_t GetAlpha() { return 2 * 3.14159265358979323846 / fgkNsect; }; 
+  static  Double_t GetAlpha()  { return 2 * 3.14159265358979323846 / fgkNsect; }; 
 
  protected:
 
@@ -130,7 +149,9 @@ class AliTRDgeometry : public AliGeometry {
 
   Int_t                fRowMax[kNplan][kNcham][kNsect];     // Number of pad-rows
   Int_t                fColMax[kNplan];                     // Number of pad-columns
-  Int_t                fTimeMax;                            // Number of time buckets
+  Int_t                fTimeMax;                            // Number of timebins in the drift region
+  Int_t                fTimeBefore;                         // Number of timebins before the drift region
+  Int_t                fTimeAfter;                          // Number of timebins after the drift region
 
   Float_t              fCwidth[kNplan];                     // Width of the chambers
 
index da9e6164c4cb7a12dbbc20d78ea584e06e44ba65..cd03aa127e3301444057e379662e84118cc1aa47 100644 (file)
@@ -27,6 +27,10 @@ class AliTRDgeometryFull : public AliTRDgeometry {
           Bool_t  GetPHOShole() const { return fPHOShole;  };
           Bool_t  GetRICHhole() const { return fRICHhole;  };
 
+  virtual Float_t GetChamberLengthI(Int_t p) { return fClengthI[p];  }; 
+  virtual Float_t GetChamberLengthM(Int_t p) { return fClengthM1[p]; }; 
+  virtual Float_t GetChamberLengthO(Int_t p) { return fClengthO1[p]; }; 
+
  protected:
 
   Bool_t          fPHOShole;                       // Switch for the hole in front of the PHOS
index 5b7bd4cda42f84b29cf665a8bef6b9aff671db1e..46749f066a3751c8c77d93c907fa23cd9900aa32 100644 (file)
@@ -20,15 +20,32 @@ class AliTRDhit : public AliHit {
   AliTRDhit(Int_t shunt, Int_t track, Int_t det, Float_t *hits, Int_t q);
   virtual ~AliTRDhit();
 
-          Int_t GetDetector() const { return fDetector; };
-          Int_t GetCharge() const   { return fQ;        };
+          Int_t  GetDetector() const         { return fDetector; };
+          Int_t  GetCharge() const           { return fQ;        };
+
+          Bool_t FromDrift() const           { return TestBit(kDrift);         };
+          Bool_t FromAmplification() const   { return TestBit(kAmplification); };
+          Bool_t FromTRphoton() const        { return TestBit(kTRphoton);      };
+          Bool_t FromTest() const            { return TestBit(kTest);          };
+
+          void   SetDrift()                  { SetBit(kDrift);         };
+          void   SetAmplification()          { SetBit(kAmplification); };
+          void   SetTRphoton()               { SetBit(kTRphoton);      };
+          void   SetTest()                   { SetBit(kTest);          }; 
 
  protected:
 
-  UShort_t     fDetector;   // TRD detector number
-  Short_t      fQ;          // Charge created by a hit. TR signals are negative.
+  enum {
+    kDrift         = 0x00000001,    // Hit is from the drift region
+    kAmplification = 0x00000002,    // Hit is from the amplification region
+    kTRphoton      = 0x00000004,    // Hit is from a TR photon
+    kTest          = 0x00000008     // Hit is a special test hit
+  };
+
+  UShort_t     fDetector;           // TRD detector number
+  Short_t      fQ;                  // Charge created by a hit. TR signals are negative.
 
-  ClassDef(AliTRDhit,3)     // Hit for the Transition Radiation Detector
+  ClassDef(AliTRDhit,3)             // Hit for the Transition Radiation Detector
 
 };
 
diff --git a/TRD/AliTRDpid.cxx b/TRD/AliTRDpid.cxx
new file mode 100644 (file)
index 0000000..1629e8a
--- /dev/null
@@ -0,0 +1,706 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//   The TRD particle identification class                                   //
+//                                                                           //
+//   Its main purposes are:                                                  //
+//      - Creation and bookkeeping of the propability distributions          //
+//      - Assignment of a e/pi propability to a given track                  //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <stdlib.h>
+#include <math.h>
+
+#include <TROOT.h>
+#include <TH1.h>
+#include <TObjArray.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TParticle.h>
+
+#include "AliRun.h"
+#include "AliTRD.h"
+#include "AliTRDpid.h"
+#include "AliTRDcluster.h"
+#include "AliTRDtrack.h"
+#include "AliTRDtracker.h"
+#include "AliTRDgeometry.h"
+
+ClassImp(AliTRDpid)
+
+//_____________________________________________________________________________
+AliTRDpid::AliTRDpid():TNamed()
+{
+  //
+  // AliTRDpid default constructor
+  // 
+
+  fNMom   = 0;
+  fMinMom = 0;
+  fMaxMom = 0;
+  fWidMom = 0;
+
+  fQHist        = NULL;
+  fLQHist       = NULL;
+  fTrackArray   = NULL;
+  fClusterArray = NULL;
+  fGeometry     = NULL;
+
+}
+
+//_____________________________________________________________________________
+AliTRDpid::AliTRDpid(const char* name, const char* title):TNamed(name,title)
+{
+  //
+  // AliTRDpid constructor
+  //
+
+  fNMom   = 0;
+  fMinMom = 0;
+  fMaxMom = 0;
+  fWidMom = 0;
+
+  fQHist        = NULL;
+  fLQHist       = NULL;
+  fTrackArray   = NULL;
+  fClusterArray = NULL;
+  fGeometry     = NULL;
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDpid::AliTRDpid(const AliTRDpid &p)
+{
+  //
+  // AliTRDpid copy constructor
+  //
+
+  ((AliTRDpid &) p).Copy(*this);
+
+}
+
+//_____________________________________________________________________________
+AliTRDpid::~AliTRDpid()
+{
+  //
+  // AliTRDpid destructor
+  //
+
+  if (fClusterArray) {
+    fClusterArray->Delete();
+    delete fClusterArray;
+  }
+
+  if (fTrackArray) {
+    fTrackArray->Delete();
+    delete fTrackArray;
+  }
+
+  if (fQHist) {
+    fQHist->Delete();
+    delete fQHist;
+  }
+
+  if (fLQHist) {
+    fLQHist->Delete();
+    delete fLQHist;
+  }
+
+}
+
+//_____________________________________________________________________________
+AliTRDpid &AliTRDpid::operator=(const AliTRDpid &p)
+{
+  //
+  // Assignment operator
+  //
+
+  if (this != &p) ((AliTRDpid &) p).Copy(*this);
+  return *this;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpid::Copy(TObject &p)
+{
+  //
+  // Copy function
+  //
+
+  fQHist->Copy(*((AliTRDpid &) p).fQHist);
+  fLQHist->Copy(*((AliTRDpid &) p).fLQHist);
+
+  ((AliTRDpid &) p).fTrackArray   = NULL;    
+  ((AliTRDpid &) p).fClusterArray = NULL;    
+  ((AliTRDpid &) p).fGeometry     = NULL;    
+
+  ((AliTRDpid &) p).fNMom   = fNMom;
+  ((AliTRDpid &) p).fMinMom = fMinMom;
+  ((AliTRDpid &) p).fMaxMom = fMaxMom;
+  ((AliTRDpid &) p).fWidMom = fWidMom;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::Init()
+{
+  //
+  // Initializes the PID object 
+  //
+
+  fClusterArray = new TObjArray();
+  fTrackArray   = new TObjArray();
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::AssignLQ(TObjArray *tarray)
+{
+  //
+  // Assigns the e / pi Q-likelihood to all tracks in the array
+  //
+
+  Bool_t status = kTRUE;
+
+  AliTRDtrack *track;
+
+  TIter nextTrack(tarray);  
+  while ((track = (AliTRDtrack *) nextTrack())) {
+    if (!AssignLQ(track)) {
+      status = kFALSE;
+      break;
+    }
+  }
+
+  return status;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::AssignLQ(AliTRDtrack *t)
+{
+  //
+  // Assigns the e / pi Q-likelihood to a given track
+  //
+
+  const Int_t kNplane = AliTRDgeometry::Nplan();
+  Float_t charge[kNplane];
+
+  // Calculate the total charge in each plane
+  if (!SumCharge(t,charge)) return kFALSE;
+
+  // Assign the likelihoods 
+  t->SetLikelihoodPion(LQPion(charge));
+  t->SetLikelihoodElectron(LQElectron(charge));
+
+  return kTRUE;  
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::CreateHistograms(const Int_t nmom
+                                 , const Float_t minmom
+                                 , const Float_t maxmom)
+{
+  //
+  // Creates the likelihood histograms
+  //
+
+  Int_t imom;
+  Int_t ipid;
+  Int_t ipla;
+
+  const Int_t kNpla = AliTRDgeometry::Nplan();
+
+  fNMom   = nmom;
+  fMinMom = minmom;
+  fMaxMom = maxmom;
+  fWidMom = (maxmom - minmom) / ((Float_t) nmom);
+
+  // The L-Q distributions
+  fLQHist = new TObjArray(kNpid * nmom);
+  for (imom = 0; imom < nmom;  imom++) {
+    for (ipid = 0; ipid < kNpid; ipid++) {
+      Int_t index = GetIndexLQ(imom,ipid);
+      Char_t name[10];
+      Char_t title[80];
+      sprintf(name ,"hLQ%03d",index);
+      if (ipid == kElectron) {
+        sprintf(title,"L-Q electron p-bin %03d",imom);
+      }
+      else {
+        sprintf(title,"L-Q pion p-bin %03d",imom);
+      }
+      TH1F *hTmp  = new TH1F(name,title,416,-0.02,1.02);
+      fLQHist->AddAt(hTmp,index);
+    }
+  }
+
+  // The Q-distributions
+  fQHist = new TObjArray(kNpla * kNpid * nmom);
+  for (imom = 0; imom < nmom;  imom++) {
+    for (ipid = 0; ipid < kNpid; ipid++) {
+      for (ipla = 0; ipla < kNpla; ipla++) {
+        Int_t index = GetIndexQ(imom,ipla,ipid);
+        Char_t name[10];
+        Char_t title[80];
+        sprintf(name ,"hQ%03d",index);
+        if (ipid == kElectron) {
+          sprintf(title,"Q electron p-bin %03d plane %d",imom,ipla);
+        }
+        else {
+          sprintf(title,"Q pion p-bin %03d plane %d",imom,ipla);
+        }
+        TH1F *hTmp  = new TH1F(name,title,100,0.0,5000.0);
+        fQHist->AddAt(hTmp,index);
+      }
+    }
+  }
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::FillQspectra()
+{
+  //
+  // Fills the Q-spectra histograms.
+  //
+  
+  Bool_t status = kTRUE;
+
+  AliTRDtrack *track;
+
+  TIter nextTrack(fTrackArray);
+  while ((track = (AliTRDtrack *) nextTrack())) {
+    if (!FillQspectra(track)) {
+      status = kFALSE;
+      break;
+    }
+  }
+
+  return status;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::FillQspectra(const AliTRDtrack *t)
+{
+  //
+  // Fills the Q-spectra histograms with track <t>.
+  //
+
+  const Int_t kNpla = AliTRDgeometry::Nplan();
+
+  if (isnan(t->GetP())) return kTRUE;
+
+  printf("----------------------------------------------------------\n");
+  Float_t charge[kNpla];
+  Float_t mom  = t->GetP();
+  Int_t   ipid = Pid(t);
+
+  if (!SumCharge(t,charge)) return kFALSE;
+
+  printf(" ipid   = %d\n",ipid);
+  printf(" mom    = %f\n",mom);
+  printf(" charge = %8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f\n"
+        ,charge[0],charge[1],charge[2],charge[3],charge[4],charge[5]);
+
+  for (Int_t ipla = 0; ipla < kNpla; ipla++) {
+    Int_t index = GetIndexQ(mom,ipla,ipid);    
+    if (index > -1) {
+      TH1F *hTmp = (TH1F *) fQHist->UncheckedAt(index);
+      hTmp->Fill(charge[ipla]);
+    }
+  }  
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::Open(const Char_t *name, Int_t event)
+{
+  //
+  // Opens and reads the kine tree, the geometry, the cluster array.
+  // and the track array from the file <name>
+  //
+
+  Bool_t status = kTRUE;
+
+  status = ReadKine(name,event);
+  status = ReadCluster(name);
+  status = ReadTracks(name);
+
+  return status;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::Open(const Char_t *namekine
+                     , const Char_t *namecluster
+                     , const Char_t *nametracks
+                     , Int_t event)
+{
+  //
+  // Opens and reads the kine tree and the geometry from file <namekine>,
+  // the cluster array from file <namecluster>,
+  // and the track array from the file <nametracks>
+  //
+
+  Bool_t status = kTRUE;
+
+  status = ReadKine(namekine,event);
+  status = ReadCluster(namecluster);
+  status = ReadTracks(nametracks);
+
+  return status;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::ReadKine(const Char_t *name, Int_t event)
+{
+  //
+  // Opens and reads the kine tree and the geometry from the file <name>.
+  //
+
+  TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
+  if (!file) {
+    printf("AliTRDpid::ReadKine -- ");
+    printf("Open file %s\n",name);
+    file = new TFile(name);
+    if (!file) {
+      printf("AliTRDpid::ReadKine -- ");
+      printf("Cannot open file %s\n",name);
+      return kFALSE;
+    }
+  }
+
+  gAlice = (AliRun *) file->Get("gAlice");
+  if (!gAlice) {
+    printf("AliTRDpid::ReadKine -- ");
+    printf("No AliRun object found\n");    
+    return kFALSE;
+  }
+  gAlice->GetEvent(event);
+
+  AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
+  if (!trd) {
+    printf("AliTRDpid::ReadKine -- ");
+    printf("No TRD object found\n");    
+    return kFALSE;
+  }
+
+  fGeometry = trd->GetGeometry();
+  if (!fGeometry) {
+    printf("AliTRDpid::ReadKine -- ");
+    printf("No TRD geometry found\n");
+    return kFALSE;
+  }
+
+  file->Close(); 
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::ReadCluster(const Char_t *name)
+{
+  //
+  // Opens and reads the cluster array from the file <name>.
+  //
+
+  fClusterArray->Delete();
+
+  printf("AliTRDpid::ReadCluster -- ");
+  printf("Open file %s\n",name);
+
+  AliTRDtracker *tracker = new AliTRDtracker("dummy","dummy");
+  tracker->ReadClusters(fClusterArray,name);
+
+  if (!fClusterArray) {
+    printf("AliTRDpid::ReadCluster -- ");
+    printf("Error reading the cluster array from file %s\n",name);
+    return kFALSE;
+  }
+
+  delete tracker;
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::ReadTracks(const Char_t *name)
+{
+  //
+  // Opens and reads the track array from the file <name>.
+  //
+
+  fTrackArray->Delete();
+
+  TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
+  if (!file) {
+    printf("AliTRDpid::ReadTracks -- ");
+    printf("Open file %s\n",name);
+    file = new TFile(name);
+    if (!file) {
+      printf("AliTRDpid::ReadTracks -- ");
+      printf("Cannot open file %s\n",name);
+      return kFALSE;
+    }
+  }
+
+  TTree   *trackTree   = (TTree *) file->Get("TreeT");
+  TBranch *trackBranch = trackTree->GetBranch("tracks");
+
+  Int_t nEntry = ((Int_t) trackTree->GetEntries());
+  for (Int_t iEntry = 0; iEntry < nEntry; iEntry++) {
+    AliTRDtrack *track = new AliTRDtrack();
+    trackBranch->SetAddress(&track);
+    trackTree->GetEvent(iEntry);
+    fTrackArray->AddLast(track);
+  }
+
+  file->Close();
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDpid::GetIndexQ(const Float_t mom, const Int_t ipla, const Int_t ipid)
+{
+  //
+  // Returns the Q-spectrum histogram index
+  //
+
+  if ((mom < fMinMom) || (mom > fMaxMom))  return -1;
+  Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
+  return GetIndexQ(imom,ipla,ipid);
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDpid::GetIndexQ(const Int_t imom, const Int_t ipla, const Int_t ipid)
+{
+  //
+  // Returns the Q-spectrum histogram index
+  //
+
+  const Int_t kNpla = AliTRDgeometry::Nplan();
+  if ((ipid < 0) || (ipid >= kNpid)) return -1;
+  return imom * (kNpid * kNpla) + ipla * kNpid + ipid; 
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDpid::GetIndexLQ(const Float_t mom, const Int_t ipid)
+{
+  //
+  // Returns the Q-likelihood histogram index
+  //
+
+  if ((mom < fMinMom) || (mom > fMaxMom))  return -1;
+  Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
+  return GetIndexLQ(imom,ipid);
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDpid::GetIndexLQ(const Int_t imom, const Int_t ipid)
+{
+  //
+  // Returns the Q-likelihood histogram index
+  //
+
+  if ((ipid < 0) || (ipid >= kNpid)) return -1;
+  return imom * kNpid + ipid; 
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDpid::Pid(const AliTRDtrack *t)
+{
+  //
+  // Determines the pid of the track <t>
+  // For a given track to be assigned an electron or pion pid it requires
+  // that more than a fraction of 0.7 of all associated cluster are 'pure'
+  // clusters -- meaning to have only contribution from one particle --
+  // of the specific particle type.
+  //
+
+  // PDG codes
+  const Int_t kPdgEl =  11;
+  const Int_t kPdgPi = 211;
+
+  // Minimum fraction of cluster from one particle
+  const Float_t kRatioMin = 0.7;
+
+  AliTRDcluster *cluster;
+  TParticle     *particle;
+
+  Int_t   nClusterEl = 0;
+  Int_t   nClusterPi = 0;
+  Int_t   nClusterNo = 0;
+
+  Float_t ratioEl    = 0;
+  Float_t ratioPi    = 0;
+
+  Int_t   ipid       = -1;
+
+  if (!fClusterArray) {
+    printf("AliTRDpid::Pid -- ");
+    printf("ClusterArray not defined. Initialize the PID object first\n");
+    return -1;  
+  }
+  
+  // Loop through all clusters associated to this track
+  Int_t nCluster = t->GetNclusters();
+  for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
+
+    // Get a cluster
+    Int_t index = t->GetClusterIndex(iCluster);
+    if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
+      break;
+    } 
+
+    // Get the first two MC track indices
+    Int_t track0 = cluster->GetTrackIndex(0);
+    Int_t track1 = cluster->GetTrackIndex(1);
+    printf(" track0 = %d, track1 = %d\n",track0,track1);
+
+    // Take only 'pure' cluster 
+    if ((track0 > -1) && (track1 == -1)) {
+
+      particle = gAlice->Particle(track0);
+      switch (TMath::Abs(particle->GetPdgCode())) {
+      case kPdgEl:
+        nClusterEl++;
+        break;
+      case kPdgPi:
+        nClusterPi++;
+        break;
+      default:
+        nClusterNo++;
+        break;
+      };
+
+    }
+
+  }
+  
+  if (nCluster) {
+    ratioEl = ((Float_t) nClusterEl) / ((Float_t) nCluster);
+    ratioPi = ((Float_t) nClusterPi) / ((Float_t) nCluster);
+    if      (ratioEl > kRatioMin) {
+      ipid = kElectron;
+    }
+    else if (ratioPi > kRatioMin) {
+      ipid = kPion;
+    }
+  }
+  printf(" nCluster = %d, nClusterEl = %d, nClusterPi = %d, nClusterNo = %d\n"
+        ,nCluster,nClusterEl,nClusterPi,nClusterNo);
+  printf(" ratioEl = %f, ratioPi = %f\n",ratioEl,ratioPi);
+
+  return ipid;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDpid::SumCharge(const AliTRDtrack *t, Float_t *charge)
+{
+  //
+  // Sums up the charge in each plane for track <t>.
+  //
+
+  Bool_t status = kTRUE;
+
+  AliTRDcluster *cluster;
+
+  const Int_t kNplane = AliTRDgeometry::Nplan();
+  for (Int_t iPlane = 0; iPlane < kNplane; iPlane++) {
+    charge[iPlane] = 0.0;
+  }
+
+  if (!fClusterArray) {
+    printf("AliTRDpid::SumCharge -- ");
+    printf("ClusterArray not defined. Initialize the PID object first\n");
+    return kFALSE;  
+  }
+  if (!fGeometry) {
+    printf("AliTRDpid::SumCharge -- ");
+    printf("TRD geometry not defined. Initialize the PID object first\n");
+    return kFALSE;  
+  }
+  
+  // Loop through all clusters associated to this track
+  Int_t nCluster = t->GetNclusters();
+  for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
+
+    // Get a cluster
+    Int_t index = t->GetClusterIndex(iCluster);
+    if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
+      status = kFALSE;
+      break;
+    } 
+
+    // Sum up the charge
+    Int_t plane    = fGeometry->GetPlane(cluster->GetDetector());
+    charge[plane] += cluster->GetQ();
+
+  }
+
+  return status;
+
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDpid::LQElectron(const Float_t *charge)
+{
+  //
+  // Returns the Q-likelihood to be a electron
+  //
+
+  return 1;
+
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDpid::LQPion(const Float_t *charge)
+{
+  //
+  // Returns the Q-likelihood to be a pion
+  //
+
+  return 0;
+
+}
+
diff --git a/TRD/AliTRDpid.h b/TRD/AliTRDpid.h
new file mode 100644 (file)
index 0000000..886d011
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef ALITRDPID_H
+#define ALITRDPID_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <TNamed.h>
+
+class TObjArray;
+class TFile;
+
+class AliTRDgeometry;
+class AliTRDtrack;
+
+class AliTRDpid : public TNamed {
+
+ public:
+
+  AliTRDpid();
+  AliTRDpid(const char* name, const char* title);
+  AliTRDpid(const AliTRDpid &p);
+  virtual ~AliTRDpid();
+  AliTRDpid &operator=(const AliTRDpid &p);
+
+  virtual void          Copy(TObject &p);
+  virtual Bool_t        Init();
+  virtual Bool_t        AssignLQ(TObjArray *tarray);
+  virtual Bool_t        AssignLQ(AliTRDtrack *t);
+  virtual Bool_t        FillQspectra();
+  virtual Bool_t        FillQspectra(const AliTRDtrack *t);
+  virtual Bool_t        CreateHistograms(const Int_t nmom
+                                       , const Float_t minmom
+                                       , const Float_t maxmom);
+  virtual Float_t       LQPion(const Float_t *charge);
+  virtual Float_t       LQElectron(const Float_t *charge);
+  virtual Bool_t        Open(const Char_t *name, Int_t event = 0);
+  virtual Bool_t        Open(const Char_t *namekine
+                           , const Char_t *namecluster
+                           , const Char_t *nametracks, Int_t event = 0); 
+  virtual Int_t         Pid(const AliTRDtrack *t);
+  virtual Bool_t        ReadCluster(const Char_t *name);
+  virtual Bool_t        ReadTracks(const Char_t *name);
+  virtual Bool_t        ReadKine(const Char_t *name, Int_t event);
+  virtual Bool_t        SumCharge(const AliTRDtrack *t, Float_t *charge);
+
+  inline  Int_t         GetIndexLQ(const Int_t imom, const Int_t ipid);
+  inline  Int_t         GetIndexLQ(const Float_t mom, const Int_t ipid);
+  inline  Int_t         GetIndexQ(const Int_t imom, const Int_t ipla, const Int_t ipid);   
+  inline  Int_t         GetIndexQ(const Float_t mom, const Int_t ipla, const Int_t ipid);
+
+          TObjArray*    GetQHist() const                    { return fQHist;  };
+          TObjArray*    GetLQHist() const                   { return fLQHist; };
+
+          void          SetGeometry(AliTRDgeometry *geo)    { fGeometry     = geo;    };
+          void          SetTrackArray(TObjArray *tarray)    { fTrackArray   = tarray; };
+          void          SetClusterArray(TObjArray *carray)  { fClusterArray = carray; };
+
+ protected:
+
+  enum { 
+    kNpid     = 2,                   //  Number of pid types (pion + electron)
+    kElectron = 0,                   //  Electron pid
+    kPion     = 1                    //  Pion pid
+  };
+
+  Int_t           fNMom;             //  Number of momentum bins
+  Float_t         fMinMom;           //  Lower momentum
+  Float_t         fMaxMom;           //  Upper momentum
+  Float_t         fWidMom;           //  Width of the momentum bins
+  TObjArray      *fLQHist;           //  Array of L-Q histograms
+  TObjArray      *fQHist;            //  Array of Q histograms
+  TObjArray      *fTrackArray;       //! Array containing the tracks
+  TObjArray      *fClusterArray;     //! Array containing the cluster
+  AliTRDgeometry *fGeometry;         //! The TRD geometry
+
+  ClassDef(AliTRDpid,1)              //  Assigns e/pi propability to the tracks 
+
+};
+#endif
index 0b45a1064782b7e3e4a9ae65bf1084fbd480dd5c..ef7a76ba5474d302bc06f00c486b1125358e6bde 100644 (file)
@@ -104,7 +104,7 @@ class AliTRDsim : public TObject {
 
   Double_t *fSigma;                //[fSpNBins] Array of sigma values
 
-  TH1D     *fSpectrum;             // TR photon energy spectrum
+  TH1D     *fSpectrum;             //!TR photon energy spectrum
 
   ClassDef(AliTRDsim,1)            // Simulates TR photons
 
index 57ebd5eb32c5a8e88cd8bb712f6a3171ab6af2ed..a2b5f4c795c001cf9c1f163a425c99864c3396e8 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.5  2001/02/05 14:49:32  hristov
+Compare() declared const (R.Brun)
+
 Revision 1.4  2000/12/08 16:07:02  cblume
 Update of the tracking by Sergei
 
@@ -65,6 +68,10 @@ const Double_t cc[15], Double_t xref, Double_t alpha) {
 
   fN=0;
   fIndex[fN++]=index;
+
+  fLhPion     = 0.0;
+  fLhElectron = 0.0;
+
 }                              
            
 //_____________________________________________________________________________
@@ -91,6 +98,10 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) {
 
   fN=t.fN;
   for (Int_t i=0; i<fN; i++) fIndex[i]=t.fIndex[i];
+
+  fLhPion     = t.fLhPion;
+  fLhElectron = t.fLhElectron;
+
 }                                                       
 
 //_____________________________________________________________________________
@@ -399,6 +410,10 @@ void AliTRDtrack::Streamer(TBuffer &R__b)
       R__b >> fCtt;
       R__b >> fN;
       for (Int_t i=0; i<fN; i++) R__b >> fIndex[i];
+      if (R__v > 1) {
+        R__b >> fLhElectron;
+        R__b >> fLhPion;
+      }
    } else {                                
       R__b.WriteVersion(AliTRDtrack::IsA());
       TObject::Streamer(R__b);
@@ -429,6 +444,8 @@ void AliTRDtrack::Streamer(TBuffer &R__b)
       R__b << fCtt;
       R__b << fN;
       for (Int_t i=0; i<fN; i++) R__b << fIndex[i];
+      R__b << fLhElectron;
+      R__b << fLhPion;
    }
 }                                                          
 
index 17ac8f15960734d44a6f874053847c01d10c0566..540bf8f6e48704f1e980a7c062226169b046228c 100644 (file)
@@ -34,7 +34,8 @@ public:
      return TMath::Abs(GetPt())*sqrt(1.+GetTgl()*GetTgl());
    }
    Double_t GetPredictedChi2(const AliTRDcluster*) const ;
-   Double_t GetPt()    const {return 0.3*0.2/GetC()/100;}
+   //Double_t GetPt()    const {return 0.3*0.2/GetC()/100;}
+   Double_t GetPt()    const {return 0.3*0.4/GetC()/100;}
    void     GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const ;
    Double_t GetSigmaC2()   const {return fCcc;}
    Double_t GetSigmaTgl2() const {return fCtt;}
@@ -45,6 +46,9 @@ public:
    Double_t GetY()   const {return fY;} // returns running Y
    Double_t GetZ()   const {return fZ;}
 
+   Float_t  GetLikelihoodPion()     const { return fLhPion;     };
+   Float_t  GetLikelihoodElectron() const { return fLhElectron; };
+
    Bool_t   IsSortable() const {return kTRUE;}
 
    Int_t    PropagateTo(Double_t xr,
@@ -58,6 +62,8 @@ public:
 
    void     Update(const AliTRDcluster* c, Double_t chi2, UInt_t i);
 
+   void     SetLikelihoodPion(Float_t l)     { fLhPion     = l; };
+   void     SetLikelihoodElectron(Float_t l) { fLhElectron = l; };
 
 protected:
 
@@ -82,9 +88,11 @@ protected:
 
    Short_t fN;             // number of clusters associated with the track
    UInt_t  fIndex[200];    // global indexes of these clusters  
-                          
 
-   ClassDef(AliTRDtrack,1)  // TRD reconstructed tracks
+   Float_t fLhElectron;    // Likelihood to be an electron
+   Float_t fLhPion;        // Likelihood to be a pion             
+
+   ClassDef(AliTRDtrack,2)  // TRD reconstructed tracks
 };                     
 
 
index 7b6739588215f4719006bc47465a8dca370b18fd..4cda5e520bf796474fa0df95f01bfd96287bf711 100644 (file)
@@ -15,6 +15,9 @@
                                                       
 /*
 $Log$
+Revision 1.9  2001/02/14 18:22:26  cblume
+Change in the geometry of the padplane
+
 Revision 1.8  2000/12/20 13:00:44  cblume
 Modifications for the HP-compiler
 
@@ -680,7 +683,9 @@ Int_t option)
   TObjArray *ioArray = new TObjArray(400);
 
   if( option < 0 ) {
-    branch = tree->GetBranch("RecPoints");
+    //branch = tree->GetBranch("RecPoints");
+    // changed CBL
+    branch = tree->GetBranch("TRDrecPoints");
 
     for (Int_t i=0; i<nentr; i++) {
       branch->SetAddress(&ioArray);
@@ -696,7 +701,7 @@ Int_t option)
     }
   }
   else {
-    branch = tree->GetBranch("Clusters");
+    branch = tree->GetBranch("TRDcluster");
 
     for (Int_t i=0; i<nentr; i++) {
       branch->SetAddress(&ioArray);
index 6caf38a6d40ff8d67e55af01a963d85f443c3d16..e0b7d8844eef7878b21884689f1850070e2777e1 100644 (file)
@@ -12,7 +12,7 @@ class TParticlePDG;
 class TObjArray;
 
 class AliTRDgeometry;
-// class AliTRDtrackingSector;
+class AliTRDtrackingSector;
 class AliTRDtrack;
 class AliTRDmcTrack;