Add new TRD classes
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Feb 2000 19:05:17 +0000 (19:05 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Feb 2000 19:05:17 +0000 (19:05 +0000)
29 files changed:
TRD/AliTRDarrayI.cxx [new file with mode: 0644]
TRD/AliTRDarrayI.h [new file with mode: 0644]
TRD/AliTRDclusterizer.cxx [new file with mode: 0644]
TRD/AliTRDclusterizer.h [new file with mode: 0644]
TRD/AliTRDclusterizerV0.cxx [new file with mode: 0644]
TRD/AliTRDclusterizerV0.h [new file with mode: 0644]
TRD/AliTRDclusterizerV1.cxx [new file with mode: 0644]
TRD/AliTRDclusterizerV1.h [new file with mode: 0644]
TRD/AliTRDdataArray.cxx [new file with mode: 0644]
TRD/AliTRDdataArray.h [new file with mode: 0644]
TRD/AliTRDdetDigits.cxx [new file with mode: 0644]
TRD/AliTRDdetDigits.h [new file with mode: 0644]
TRD/AliTRDdigitizer.cxx [new file with mode: 0644]
TRD/AliTRDdigitizer.h [new file with mode: 0644]
TRD/AliTRDgeometry.cxx [new file with mode: 0644]
TRD/AliTRDgeometry.h [new file with mode: 0644]
TRD/AliTRDgeometryFull.cxx [new file with mode: 0644]
TRD/AliTRDgeometryFull.h [new file with mode: 0644]
TRD/AliTRDgeometryHole.cxx [new file with mode: 0644]
TRD/AliTRDgeometryHole.h [new file with mode: 0644]
TRD/AliTRDrecPoint.cxx [new file with mode: 0644]
TRD/AliTRDrecPoint.h [new file with mode: 0644]
TRD/AliTRDsegmentArray.cxx [new file with mode: 0644]
TRD/AliTRDsegmentArray.h [new file with mode: 0644]
TRD/AliTRDsegmentArrayBase.cxx [new file with mode: 0644]
TRD/AliTRDsegmentArrayBase.h [new file with mode: 0644]
TRD/AliTRDsegmentID.cxx [new file with mode: 0644]
TRD/AliTRDsegmentID.h [new file with mode: 0644]
TRD/slowClusterAna.C [new file with mode: 0644]

diff --git a/TRD/AliTRDarrayI.cxx b/TRD/AliTRDarrayI.cxx
new file mode 100644 (file)
index 0000000..609e671
--- /dev/null
@@ -0,0 +1,52 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////
+//   Added additional functionality  to original TArrayI             //
+//   multiple inheritance from TObject to be possible use automatic  //
+//   branch mechanism for tree                                       //
+//   function Expand to be possible expand array without deleting    //
+//   array contents                                                  //
+//                                                                   //
+//  Origin:  Marian Ivanov, Uni. of Bratislava, ivanov@fmph.uniba.sk // 
+//                                                                   //  
+///////////////////////////////////////////////////////////////////////
+
+#include "AliTRDarrayI.h"
+
+ClassImp(AliTRDarrayI)
+
+AliTRDarrayI::~AliTRDarrayI()
+{
+  //
+  //default destructor
+}
+
+void AliTRDarrayI::Expand(Int_t n)
+{
+  //
+  // Set array size of TArrayI object to n integers and copy old array
+  // If n<0 leave array unchanged.
+  // user are responsible for appopiate size of array
+  //
+  if (n < 0) return;  
+  fArray = (Int_t*)  TStorage::ReAlloc(fArray, n * sizeof(Int_t),fN * sizeof(Int_t));
+  if (fArray!=0) fN= n; 
+}
+
diff --git a/TRD/AliTRDarrayI.h b/TRD/AliTRDarrayI.h
new file mode 100644 (file)
index 0000000..fdae3ec
--- /dev/null
@@ -0,0 +1,23 @@
+#ifndef AliTRDArrayI_H
+#define AliTRDArrayI_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDarrayI.h,v */
+
+#include "TObject.h"
+#include "TArrayI.h"
+
+class AliTRDarrayI: public TObject ,public TArrayI {
+
+public:
+
+  ~AliTRDarrayI();
+  void Expand(Int_t n);  
+
+  ClassDef(AliTRDarrayI,1)  
+
+};
+
+#endif 
+
diff --git a/TRD/AliTRDclusterizer.cxx b/TRD/AliTRDclusterizer.cxx
new file mode 100644 (file)
index 0000000..67b8c12
--- /dev/null
@@ -0,0 +1,169 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  TRD cluster finder base class                                            //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliRun.h"
+
+#include "AliTRD.h"
+#include "AliTRDclusterizer.h"
+
+ClassImp(AliTRDclusterizer)
+
+//_____________________________________________________________________________
+AliTRDclusterizer::AliTRDclusterizer():TNamed()
+{
+  //
+  // AliTRDclusterizer default constructor
+  //
+
+  fInputFile = NULL;
+  fEvent     = 0;
+
+}
+
+//_____________________________________________________________________________
+AliTRDclusterizer::AliTRDclusterizer(const Text_t* name, const Text_t* title)
+                  :TNamed(name, title)
+{
+  //
+  // AliTRDclusterizer default constructor
+  //
+
+  fInputFile = NULL;
+  fEvent     = 0;
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDclusterizer::~AliTRDclusterizer()
+{
+
+  if (fInputFile) {
+    fInputFile->Close();
+    delete fInputFile;
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDclusterizer::Init()
+{
+  //
+  // Initializes the cluster finder
+  //
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
+{
+  //
+  // Opens a ROOT-file with TRD-hits and reads in the digits-tree
+  //
+
+  // Connect the AliRoot file containing Geometry, Kine, and Hits
+  fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
+  if (!fInputFile) {
+    printf("AliTRDclusterizer::Open -- ");
+    printf("Open the ALIROOT-file %s.\n",name);
+    fInputFile = new TFile(name,"UPDATE");
+  }
+  else {
+    printf("AliTRDclusterizer::Open -- ");
+    printf("%s is already open.\n",name);
+  }
+
+  // Get AliRun object from file or create it if not on file
+  //if (!gAlice) {
+    gAlice = (AliRun*) fInputFile->Get("gAlice");
+    if (gAlice) {
+      printf("AliTRDclusterizer::Open -- ");
+      printf("AliRun object found on file.\n");
+    }
+    else {
+      printf("AliTRDclusterizer::Open -- ");
+      printf("Could not find AliRun object.\n");
+      return kFALSE;
+    }
+  //}
+
+  fEvent = nEvent;
+
+  // Import the Trees for the event nEvent in the file
+  Int_t nparticles = gAlice->GetEvent(fEvent);
+  if (nparticles <= 0) {
+    printf("AliTRDclusterizer::Open -- ");
+    printf("No entries in the trees for event %d.\n",fEvent);
+    return kFALSE;
+  }
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDclusterizer::WriteCluster()
+{
+  //
+  // Writes out the TRD-cluster
+  //
+
+  // Write the new tree into the input file (use overwrite option)
+  Char_t treeName[7];
+  sprintf(treeName,"TreeR%d",fEvent);
+  printf("AliTRDclusterizer::WriteCluster -- ");
+  printf("Write the cluster tree %s for event %d.\n"
+        ,treeName,fEvent);
+  gAlice->TreeR()->Write(treeName,2);
+
+  return kTRUE;
+
+}
+
+ClassImp(AliTRDcluster)
+
+//_____________________________________________________________________________
+AliTRDcluster::AliTRDcluster(Int_t *tracks, Int_t *cluster, Float_t energy, Float_t* position)
+              :TObject()
+{
+  //
+  // Create a TRD cluster
+  //
+
+  fDetector  = cluster[0];
+
+  fTimeSlice = cluster[1];
+  fEnergy    = energy;
+
+  fX         = position[0];
+  fY         = position[1];
+  fZ         = position[2];
+
+  fTracks[0] = tracks[0];
+  fTracks[1] = tracks[1];
+  fTracks[2] = tracks[2];
+
+}
diff --git a/TRD/AliTRDclusterizer.h b/TRD/AliTRDclusterizer.h
new file mode 100644 (file)
index 0000000..5638104
--- /dev/null
@@ -0,0 +1,66 @@
+#ifndef TRDclusterizer_h
+#define TRDclusterizer_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <TNamed.h>
+#include <TFile.h>
+
+///////////////////////////////////////////////////////
+//  Finds and handles cluster                        //
+///////////////////////////////////////////////////////
+
+class AliTRDclusterizer : public TNamed {
+
+ public:
+
+  AliTRDclusterizer();
+  AliTRDclusterizer(const Text_t* name, const Text_t* title);
+  ~AliTRDclusterizer();
+
+  virtual void    Init();
+  virtual Bool_t  Open(const Char_t *name, Int_t nEvent = 0);
+  virtual Bool_t  MakeCluster() = 0;
+  virtual Bool_t  WriteCluster();
+
+ protected:
+
+  TFile   *fInputFile;             //! AliROOT input file
+  
+  Int_t    fEvent;                 //! Event number
+
+  ClassDef(AliTRDclusterizer,1)    // TRD-Cluster manager base class
+
+};
+
+//_____________________________________________________________________________
+class AliTRDcluster : public TObject {
+
+public:
+
+  Int_t    fDetector;       // TRD detector number
+
+  Int_t    fTimeSlice;      // Timeslice in chamber where cluster has been found
+  Float_t  fEnergy;         // Charge sum of this cluster
+
+  Float_t  fX;              // X coord in ALICE reference frame
+  Float_t  fY;              // Y coord in ALICE reference frame
+  Float_t  fZ;              // Z coord in ALICE reference frame
+
+  Int_t    fTracks[3];      // Track information
+
+public:
+
+  AliTRDcluster() {};
+  AliTRDcluster(Int_t *tracks, Int_t *cluster, Float_t energy, Float_t *pos);
+  virtual ~AliTRDcluster() {};
+
+  inline virtual Int_t *GetTracks() { return &fTracks[0]; }
+
+  ClassDef(AliTRDcluster,1) // Cluster for Transition Radiation Detector
+
+};
+
+#endif
diff --git a/TRD/AliTRDclusterizerV0.cxx b/TRD/AliTRDclusterizerV0.cxx
new file mode 100644 (file)
index 0000000..c55b93a
--- /dev/null
@@ -0,0 +1,273 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// TRD cluster finder for the fast simulator. It takes the hits from the     //
+// fast simulator (one hit per plane) and transforms them                    //
+// into cluster, by applying position smearing and merging                   //
+// of nearby cluster. The searing is done uniformly in z-direction           //
+// over the length of a readout pad. In rphi-direction a Gaussian            //
+// smearing is applied with a sigma given by fRphiSigma.                     //
+// Clusters are considered as overlapping when they are closer in            //
+// rphi-direction than the value defined in fRphiDist.                       //
+// Use the macro fastClusterCreate.C to create the cluster.                  //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TRandom.h>
+
+#include "AliTRDclusterizerV0.h"
+#include "AliTRDconst.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDrecPoint.h"
+
+ClassImp(AliTRDclusterizerV0)
+
+//_____________________________________________________________________________
+AliTRDclusterizerV0::AliTRDclusterizerV0():AliTRDclusterizer()
+{
+  //
+  // AliTRDclusterizerV0 default constructor
+  //
+
+}
+
+//_____________________________________________________________________________
+AliTRDclusterizerV0::AliTRDclusterizerV0(const Text_t* name, const Text_t* title)
+                    :AliTRDclusterizer(name,title)
+{
+  //
+  // AliTRDclusterizerV0 default constructor
+  //
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDclusterizerV0::~AliTRDclusterizerV0()
+{
+
+}
+
+//_____________________________________________________________________________
+void AliTRDclusterizerV0::Init()
+{
+  //
+  // Initializes the cluster finder
+  //
+
+  // Position resolution in rphi-direction
+  fRphiSigma  = 0.02;
+  // Minimum distance of non-overlapping cluster
+  fRphiDist   = 1.0;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDclusterizerV0::MakeCluster()
+{
+  //
+  // Generates the cluster
+  //
+
+  // Get the pointer to the detector class and check for version 1
+  AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
+  if (TRD->IsVersion() != 0) {
+    printf("AliTRDclusterizerV0::MakeCluster -- ");
+    printf("TRD must be version 0 (fast simulator).\n");
+    return kFALSE; 
+  }
+
+  // Get the geometry
+  AliTRDgeometry *Geo = TRD->GetGeometry();
+  
+  printf("AliTRDclusterizerV0::MakeCluster -- ");
+  printf("Start creating cluster.\n");
+
+  Int_t nBytes = 0;
+
+  AliTRDhit      *Hit;
+  
+  // Get the pointer to the hit tree
+  TTree *HitTree     = gAlice->TreeH();
+  // Get the pointer to the reconstruction tree
+  TTree *ClusterTree = gAlice->TreeR();
+
+  TObjArray *Chamber = new TObjArray();
+
+  // Get the number of entries in the hit tree
+  // (Number of primary particles creating a hit somewhere)
+  Int_t nTrack = (Int_t) HitTree->GetEntries();
+
+  // Loop through all the chambers
+  for (Int_t icham = 0; icham < kNcham; icham++) {
+    for (Int_t iplan = 0; iplan < kNplan; iplan++) {
+      for (Int_t isect = 0; isect < kNsect; isect++) {
+
+        Int_t   nColMax    = Geo->GetColMax(iplan);
+        Float_t row0       = Geo->GetRow0(iplan,icham,isect);
+        Float_t col0       = Geo->GetCol0(iplan);
+
+        Float_t rowPadSize = Geo->GetRowPadSize();
+        Float_t colPadSize = Geo->GetColPadSize();
+
+        // Loop through all entries in the tree
+        for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
+
+          gAlice->ResetHits();
+          nBytes += HitTree->GetEvent(iTrack);
+
+          // Get the number of hits in the TRD created by this particle
+          Int_t nHit = TRD->Hits()->GetEntriesFast();
+
+          // Loop through the TRD hits  
+          for (Int_t iHit = 0; iHit < nHit; iHit++) {
+
+            if (!(Hit = (AliTRDhit *) TRD->Hits()->UncheckedAt(iHit))) 
+              continue;
+
+            Float_t pos[3];
+                    pos[0]   = Hit->fX;
+                    pos[1]   = Hit->fY;
+                    pos[2]   = Hit->fZ;
+            Int_t   track    = Hit->fTrack;
+            Int_t   detector = Hit->fDetector;
+            Int_t   plane    = Geo->GetPlane(detector);
+            Int_t   sector   = Geo->GetSector(detector);
+            Int_t   chamber  = Geo->GetChamber(detector);        
+
+            if ((sector  != isect) ||
+                (plane   != iplan) ||
+                (chamber != icham)) 
+              continue;
+
+            // Rotate the sectors on top of each other
+            Float_t rot[3];
+            Geo->Rotate(detector,pos,rot);
+
+            // Add this cluster to the temporary cluster-array for this chamber
+            Int_t   tracks[3];
+            tracks[0] = track;
+            Int_t   clusters[2];
+            clusters[0] = detector;
+            clusters[1] = 0;
+            Float_t position[3];
+            position[0] = rot[2];
+            position[1] = rot[1];
+            position[2] = rot[0];
+           AliTRDcluster *Cluster = new AliTRDcluster(tracks,clusters,0,position);
+            Chamber->Add(Cluster);
+
+         }
+
+       }
+  
+        // Loop through the temporary cluster-array
+        for (Int_t iClus1 = 0; iClus1 < Chamber->GetEntries(); iClus1++) {
+
+          AliTRDcluster *Cluster1 = (AliTRDcluster *) Chamber->UncheckedAt(iClus1);
+          Float_t x1 = Cluster1->fX;
+          Float_t y1 = Cluster1->fY;
+          Float_t z1 = Cluster1->fZ;
+
+          if (!(z1)) continue;             // Skip marked cluster  
+
+          const Int_t nSave = 2;
+          Int_t idxSave[nSave];
+          Int_t iSave = 0;
+
+          Int_t tracks[3];
+          tracks[0] = Cluster1->fTracks[0];
+
+          // Check the other cluster to see, whether there are close ones
+          for (Int_t iClus2 = iClus1 + 1; iClus2 < Chamber->GetEntries(); iClus2++) {
+            AliTRDcluster *Cluster2 = (AliTRDcluster *) Chamber->UncheckedAt(iClus2);
+            Float_t x2 = Cluster2->fX;
+            Float_t y2 = Cluster2->fY;
+            if ((TMath::Abs(x1 - x2) < rowPadSize) ||
+                (TMath::Abs(y1 - y2) <  fRphiDist)) {
+              if (iSave == nSave) {
+                printf("AliTRDclusterizerV0::MakeCluster -- ");
+                printf("Boundary error: iSave = %d, nSave = %d.\n"
+                      ,iSave,nSave);
+             }
+              else {                
+                idxSave[iSave]  = iClus2;
+                tracks[iSave+1] = Cluster2->fTracks[0];
+             }
+              iSave++;
+           }
+         }
+     
+          // Merge close cluster
+          Float_t yMerge = y1;
+          Float_t xMerge = x1;
+          if (iSave) {
+            for (Int_t iMerge = 0; iMerge < iSave; iMerge++) {
+              AliTRDcluster *Cluster2 = 
+                (AliTRDcluster *) Chamber->UncheckedAt(idxSave[iMerge]);
+              xMerge += Cluster2->fX;
+              yMerge += Cluster2->fY;
+              Cluster2->fZ = 0;            // Mark merged cluster
+           }
+            xMerge /= (iSave + 1);
+            yMerge /= (iSave + 1);
+          }
+
+          Float_t smear[3];
+
+          // The position smearing in z-direction (uniform over pad width)            
+          Int_t row = (Int_t) ((xMerge - row0) / rowPadSize);
+          smear[0]  = (row + gRandom->Rndm()) * rowPadSize + row0;
+
+          // The position smearing in rphi-direction (Gaussian)
+          smear[1] = 0;
+          do
+            smear[1] = gRandom->Gaus(yMerge,fRphiSigma);
+          while ((smear[1] < col0                        ) ||
+                 (smear[1] > col0 + nColMax * colPadSize));
+
+          // Time direction stays unchanged
+          smear[2] = z1;
+
+          // Add the smeared cluster to the output array 
+          Int_t detector  = Cluster1->fDetector;
+          Int_t digits[3] = {0};
+          TRD->AddRecPoint(smear,digits,detector,0.0);
+
+       }
+
+        // Clear the temporary cluster-array and delete the cluster
+        Chamber->Delete();
+
+      }
+    }
+  }
+
+  printf("AliTRDclusterizerV0::MakeCluster -- ");
+  printf("Found %d points.\n",TRD->RecPoints()->GetEntries());
+  printf("AliTRDclusterizerV0::MakeCluster -- ");
+  printf("Fill the cluster tree.\n");
+  ClusterTree->Fill();
+
+  return kTRUE;
+
+}
diff --git a/TRD/AliTRDclusterizerV0.h b/TRD/AliTRDclusterizerV0.h
new file mode 100644 (file)
index 0000000..b5d84e7
--- /dev/null
@@ -0,0 +1,41 @@
+#ifndef TRDclusterizerV0_h
+#define TRDclusterizerV0_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliTRD.h"
+#include "AliTRDclusterizer.h"
+
+///////////////////////////////////////////////////////
+//  Finds and handles cluster (fast simulation)      //
+///////////////////////////////////////////////////////
+
+class AliTRDclusterizerV0 : public AliTRDclusterizer {
+
+ public:
+
+  AliTRDclusterizerV0();
+  AliTRDclusterizerV0(const Text_t* name, const Text_t* title);
+  ~AliTRDclusterizerV0();
+
+  virtual void    Init();
+  virtual Bool_t  MakeCluster();
+  
+  virtual void    SetRphiSigma(Float_t sigma) { fRphiSigma = sigma; };
+  virtual void    SetRphiDist(Float_t dist)   { fRphiDist  = dist;  };
+
+  virtual Float_t GetRphiSigma()              { return fRphiSigma;  };
+  virtual Float_t GetRphiDist()               { return fRphiDist;   };
+
+ protected:
+
+  Float_t      fRphiSigma;           // Gaussian position smearing in rphi-direction
+  Float_t      fRphiDist;            // Maximum distance for non-overlapping cluster
+
+  ClassDef(AliTRDclusterizerV0,1)    // TRD-Cluster manager, fast simulator
+
+};
+
+#endif
diff --git a/TRD/AliTRDclusterizerV1.cxx b/TRD/AliTRDclusterizerV1.cxx
new file mode 100644 (file)
index 0000000..afb19bb
--- /dev/null
@@ -0,0 +1,416 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// TRD cluster finder for the slow simulator. 
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TF1.h>
+
+#include "AliTRDclusterizerV1.h"
+#include "AliTRDmatrix.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDdigitizer.h"
+#include "AliTRDrecPoint.h"
+#include "AliTRDdataArray.h"
+
+ClassImp(AliTRDclusterizerV1)
+
+//_____________________________________________________________________________
+AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
+{
+  //
+  // AliTRDclusterizerV1 default constructor
+  //
+
+  fDigitsArray = NULL;
+
+}
+
+//_____________________________________________________________________________
+AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
+                    :AliTRDclusterizer(name,title)
+{
+  //
+  // AliTRDclusterizerV1 default constructor
+  //
+
+  fDigitsArray = NULL;
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDclusterizerV1::~AliTRDclusterizerV1()
+{
+
+  if (fDigitsArray) {
+    fDigitsArray->Delete();
+    delete fDigitsArray;
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDclusterizerV1::Init()
+{
+  //
+  // Initializes the cluster finder
+  //
+
+  // The default parameter for the clustering
+  fClusMaxThresh = 5.0;
+  fClusSigThresh = 2.0;
+  fClusMethod    = 1;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDclusterizerV1::ReadDigits()
+{
+  //
+  // Reads the digits arrays from the input aliroot file
+  //
+
+  if (!fInputFile) {
+    printf("AliTRDclusterizerV1::ReadDigits -- ");
+    printf("No input file open\n");
+    return kFALSE;
+  }
+
+  // Create a new segment array for the digits 
+  fDigitsArray = new AliTRDsegmentArray(kNsect*kNplan*kNcham);
+
+  // Read in the digit arrays
+  return (fDigitsArray->LoadArray("TRDdigits"));  
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDclusterizerV1::MakeCluster()
+{
+  //
+  // Generates the cluster.
+  //
+
+  Int_t row, col, time;
+
+  // Get the pointer to the detector class and check for version 1
+  AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
+  if (TRD->IsVersion() != 1) {
+    printf("AliTRDclusterizerV1::MakeCluster -- ");
+    printf("TRD must be version 1 (slow simulator).\n");
+    return kFALSE; 
+  }
+
+  // Get the geometry
+  AliTRDgeometry *Geo = TRD->GetGeometry();
+
+  printf("AliTRDclusterizerV1::MakeCluster -- ");
+  printf("Start creating clusters.\n");
+
+  AliTRDdataArray *Digits;
+
+  // Parameters
+  Float_t maxThresh        = fClusMaxThresh;   // threshold value for maximum
+  Float_t signalThresh     = fClusSigThresh;   // threshold value for digit signal
+  Int_t   clusteringMethod = fClusMethod;      // clustering method option (for testing)
+
+  // Iteration limit for unfolding procedure
+  const Float_t epsilon = 0.01;             
+
+  const Int_t   nClus   = 3;  
+  const Int_t   nSig    = 5;
+
+  Int_t chamBeg = 0;
+  Int_t chamEnd = kNcham;
+  if (TRD->GetSensChamber() >= 0) {
+    chamBeg = TRD->GetSensChamber();
+    chamEnd = chamEnd + 1;
+  }
+  Int_t planBeg = 0;
+  Int_t planEnd = kNplan;
+  if (TRD->GetSensPlane()   >= 0) {
+    planBeg = TRD->GetSensPlane();
+    planEnd = planBeg + 1;
+  }
+  Int_t sectBeg = 0;
+  Int_t sectEnd = kNsect;
+  if (TRD->GetSensSector()  >= 0) {
+    sectBeg = TRD->GetSensSector();
+    sectEnd = sectBeg + 1;
+  }
+
+  // *** Start clustering *** in every chamber
+  for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
+    for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
+      for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
+
+        Int_t idet = Geo->GetDetector(iplan,icham,isect);
+
+        Int_t nClusters = 0;
+        printf("AliTRDclusterizerV1::MakeCluster -- ");
+        printf("Analyzing chamber %d, plane %d, sector %d.\n"
+               ,icham,iplan,isect);
+
+        Int_t   nRowMax  = Geo->GetRowMax(iplan,icham,isect);
+        Int_t   nColMax  = Geo->GetColMax(iplan);
+        Int_t   nTimeMax = Geo->GetTimeMax();
+
+        // Create a detector matrix to keep maxima
+        AliTRDmatrix *digitMatrix  = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
+                                                     ,isect,icham,iplan);
+        // Create a matrix to contain maximum flags
+        AliTRDmatrix *maximaMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
+                                                     ,isect,icham,iplan);
+
+        // Read in the digits
+        Digits = (AliTRDdataArray *) fDigitsArray->At(idet);
+
+        // Loop through the detector pixel
+        for (time = 0; time < nTimeMax; time++) {
+          for ( col = 0;  col <  nColMax;  col++) {
+            for ( row = 0;  row <  nRowMax;  row++) {
+
+              Int_t signal = Digits->GetData(row,col,time);
+              Int_t index  = Digits->GetIndex(row,col,time);
+
+              // Fill the detector matrix
+              if (signal > signalThresh) {
+               // Store the signal amplitude
+                digitMatrix->SetSignal(row,col,time,signal);
+               // Store the digits number
+                digitMatrix->AddTrack(row,col,time,index);
+              }
+
+           }
+         }
+       }
+
+        // Loop chamber and find maxima in digitMatrix
+        for ( row = 0;  row <  nRowMax;  row++) {
+          for ( col = 1;  col <  nColMax;  col++) {
+            for (time = 0; time < nTimeMax; time++) {
+
+              if (digitMatrix->GetSignal(row,col,time) 
+                  < digitMatrix->GetSignal(row,col - 1,time)) {
+                // really maximum?
+                if (col > 1) {
+                  if (digitMatrix->GetSignal(row,col - 2,time)
+                      < digitMatrix->GetSignal(row,col - 1,time)) {
+                    // yes, so set maximum flag
+                    maximaMatrix->SetSignal(row,col - 1,time,1);
+                  }
+                  else maximaMatrix->SetSignal(row,col - 1,time,0);
+                }
+              }
+
+            }   // time
+          }     // col
+        }       // row
+
+        // now check maxima and calculate cluster position
+        for ( row = 0;  row <  nRowMax;  row++) {
+          for ( col = 1;  col <  nColMax;  col++) {
+            for (time = 0; time < nTimeMax; time++) {
+
+              if ((maximaMatrix->GetSignal(row,col,time) > 0)
+                  && (digitMatrix->GetSignal(row,col,time) > maxThresh)) {
+
+                // Ratio resulting from unfolding
+                Float_t ratio                =  0;    
+                // Signals on max and neighbouring pads
+                Float_t padSignal[nSig]      = {0};   
+                // Signals from cluster
+                Float_t clusterSignal[nClus] = {0};
+                // Cluster pad info
+                Float_t clusterPads[nClus]   = {0};   
+                // Cluster digit info
+                Int_t   clusterDigit[nClus]  = {0};
+
+                for (Int_t iPad = 0; iPad < nClus; iPad++) {
+                  clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
+                  clusterDigit[iPad]  = digitMatrix->GetTrack(row,col-1+iPad,time,0);
+                }
+
+                // neighbouring maximum on right side?
+                if (col < nColMax - 2) {
+                  if (maximaMatrix->GetSignal(row,col + 2,time) > 0) {
+
+                    for (Int_t iPad = 0; iPad < 5; iPad++) {
+                      padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
+                    }
+
+                    // unfold:
+                    ratio = Unfold(epsilon, padSignal);
+
+                    // set signal on overlapping pad to ratio
+                    clusterSignal[2] *= ratio;
+
+                  }
+                }
+                
+               // Calculate the position of the cluster
+                switch (clusteringMethod) {
+                case 1:
+                  // method 1: simply center of mass
+                  clusterPads[0] = row + 0.5;
+                  clusterPads[1] = col - 0.5 + (clusterSignal[2] - clusterSignal[0]) /
+                                   (clusterSignal[1] + clusterSignal[2] + clusterSignal[3]);
+                  clusterPads[2] = time + 0.5;
+
+                  nClusters++;
+                  break;
+                case 2:
+                  // method 2: integral gauss fit on 3 pads
+                  TH1F *hPadCharges = new TH1F("hPadCharges", "Charges on center 3 pads"
+                                                           , 5, -1.5, 3.5);
+                  for (Int_t iCol = -1; iCol <= 3; iCol++) {
+                    if (clusterSignal[iCol] < 1) clusterSignal[iCol] = 1;
+                    hPadCharges->Fill(iCol, clusterSignal[iCol]);
+                  }
+                  hPadCharges->Fit("gaus", "IQ", "SAME", -0.5, 2.5);
+                  TF1     *fPadChargeFit = hPadCharges->GetFunction("gaus");
+                  Double_t  colMean = fPadChargeFit->GetParameter(1);
+
+                  clusterPads[0] = row + 0.5;
+                  clusterPads[1] = col - 1.5 + colMean;
+                  clusterPads[2] = time + 0.5;
+
+                  delete hPadCharges;
+
+                  nClusters++;
+                  break;
+                }
+
+                Float_t clusterCharge =   clusterSignal[0]
+                                        + clusterSignal[1]
+                                        + clusterSignal[2];
+
+                // Add the cluster to the output array 
+                TRD->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge);
+
+              }
+            }  // time
+          }    // col
+        }      // row
+
+        printf("AliTRDclusterizerV1::MakeCluster -- ");
+        printf("Number of clusters found: %d\n",nClusters);
+
+        delete digitMatrix;
+        delete maximaMatrix;
+
+      }          // isect
+    }            // iplan
+  }              // icham
+
+  printf("AliTRDclusterizerV1::MakeCluster -- ");
+  printf("Total number of points found: %d\n"
+        ,TRD->RecPoints()->GetEntries());
+
+  // Get the pointer to the cluster branch
+  TTree *ClusterTree = gAlice->TreeR(); 
+
+  // Fill the cluster-branch
+  printf("AliTRDclusterizerV1::MakeCluster -- ");
+  printf("Fill the cluster tree.\n");
+  ClusterTree->Fill();
+  printf("AliTRDclusterizerV1::MakeCluster -- ");
+  printf("Done.\n");
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
+{
+  //
+  // Method to unfold neighbouring maxima.
+  // The charge ratio on the overlapping pad is calculated
+  // until there is no more change within the range given by eps.
+  // The resulting ratio is then returned to the calling method.
+  //
+
+  Int_t   itStep            = 0;      // count iteration steps
+
+  Float_t ratio             = 0.5;    // start value for ratio
+  Float_t prevRatio         = 0;      // store previous ratio
+
+  Float_t newLeftSignal[3]  = {0};    // array to store left cluster signal
+  Float_t newRightSignal[3] = {0};    // array to store right cluster signal
+
+  // start iteration:
+  while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
+
+    itStep++;
+    prevRatio = ratio;
+
+    // cluster position according to charge ratio
+    Float_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) /
+                       (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
+    Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) /
+                       ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
+
+    // set cluster charge ratio
+    Float_t ampLeft  = padSignal[1];
+    Float_t ampRight = padSignal[3];
+
+    // apply pad response to parameters
+    newLeftSignal[0] = ampLeft*PadResponse(-1 - maxLeft);
+    newLeftSignal[1] = ampLeft*PadResponse( 0 - maxLeft);
+    newLeftSignal[2] = ampLeft*PadResponse( 1 - maxLeft);
+
+    newRightSignal[0] = ampRight*PadResponse(-1 - maxRight);
+    newRightSignal[1] = ampRight*PadResponse( 0 - maxRight);
+    newRightSignal[2] = ampRight*PadResponse( 1 - maxRight);
+
+    // calculate new overlapping ratio
+    ratio = newLeftSignal[2]/(newLeftSignal[2] + newRightSignal[0]);
+
+  }
+
+  return ratio;
+
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
+{
+  //
+  // The pad response for the chevron pads. 
+  // We use a simple Gaussian approximation which should be good
+  // enough for our purpose.
+  //
+
+  // The parameters for the response function
+  const Float_t aa  =  0.8872;
+  const Float_t bb  = -0.00573;
+  const Float_t cc  =  0.454;
+  const Float_t cc2 =  cc*cc;
+
+  Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2)));
+
+  return (pr);
+
+}
diff --git a/TRD/AliTRDclusterizerV1.h b/TRD/AliTRDclusterizerV1.h
new file mode 100644 (file)
index 0000000..84af1af
--- /dev/null
@@ -0,0 +1,55 @@
+#ifndef TRDclusterizerV1_h
+#define TRDclusterizerV1_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliTRD.h"
+#include "AliTRDsegmentArray.h"
+#include "AliTRDclusterizer.h"
+
+///////////////////////////////////////////////////////
+//  Finds and handles cluster (slow simulation)      //
+///////////////////////////////////////////////////////
+
+class AliTRDclusterizerV1 : public AliTRDclusterizer {
+
+ public:
+
+  AliTRDclusterizerV1();
+  AliTRDclusterizerV1(const Text_t* name, const Text_t* title);
+  ~AliTRDclusterizerV1();
+
+  virtual void    Init();
+  virtual Bool_t  MakeCluster();
+  virtual Bool_t  ReadDigits();
+
+  virtual void    SetDigitsArray(AliTRDsegmentArray *Array) { fDigitsArray   = Array;  };
+
+  virtual void    SetClusMaxThresh(Float_t thresh)          { fClusMaxThresh = thresh; };
+  virtual void    SetClusSigThresh(Float_t thresh)          { fClusSigThresh = thresh; };
+  virtual void    SetClusMethod(Int_t meth)                 { fClusMethod    = meth;   };
+
+  virtual Float_t GetClusMaxThresh()                        { return fClusMaxThresh; };
+  virtual Float_t GetClusSigThresh()                        { return fClusSigThresh; };
+  virtual Int_t   GetClusMethod()                           { return fClusMethod;    };
+
+ protected:
+
+  AliTRDsegmentArray *fDigitsArray;    //! Array of detector segments containing the digits  
+
+  Float_t             fClusMaxThresh;  // Threshold value for cluster maximum
+  Float_t             fClusSigThresh;  // Threshold value for cluster signal
+  Int_t               fClusMethod;     // Clustering method
+
+ private:
+
+  virtual Float_t  Unfold(Float_t eps, Float_t *padSignal);
+  virtual Float_t  PadResponse(Float_t x);
+
+  ClassDef(AliTRDclusterizerV1,1)      // TRD-Cluster manager, slow simulator
+
+};
+
+#endif
diff --git a/TRD/AliTRDdataArray.cxx b/TRD/AliTRDdataArray.cxx
new file mode 100644 (file)
index 0000000..35b8df3
--- /dev/null
@@ -0,0 +1,657 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  General container for data of a TRD detector segment.                    //
+//  Adapted from AliDigits (origin: M.Ivanov).                               //
+//  The main difference is that we used 4 byte integer, so that this class   //
+//  can also be used as a dictionary between digits and MC particles.        //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "TClass.h"
+#include "TError.h"
+#include "AliTRDsegmentID.h"
+#include "AliTRDarrayI.h"
+#include "AliTRDdataArray.h"
+
+ClassImp(AliTRDdataArray)
+
+//_____________________________________________________________________________
+AliTRDdataArray::AliTRDdataArray()
+{
+  //
+  // Default constructor
+  //
+
+  fIndex     = 0;
+  fElements  = 0;
+  fThreshold = 0;
+  Reset();
+
+}
+
+//_____________________________________________________________________________
+AliTRDdataArray::~AliTRDdataArray()
+{
+  //
+  // Destructor
+  //
+
+  if (fIndex)    fIndex->Delete();;
+  if (fElements) fElements->Delete();
+  
+}
+
+
+//_____________________________________________________________________________
+void AliTRDdataArray::Reset() 
+{ 
+  //
+  // Reset the array (old content gets deleted)
+  //
+
+  if (fIndex)    delete fIndex;
+  fIndex    = new AliTRDarrayI;
+  
+  if (fElements) delete fElements;
+  fElements = new AliTRDarrayI;
+
+  fNdim1 = fNdim2 = fNelems = -1; 
+  fElements->Set(0); 
+  fIndex->Set(0); 
+  fBufType = -1;
+
+  fNrow  = 0;
+  fNcol  = 0;
+  fNtime = 0;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdataArray::Allocate(Int_t nrow, Int_t ncol, Int_t ntime)
+{
+  //
+  // Allocate an empty buffer of the size <nrow> x <ncol> x <ntime>
+  //
+
+  Reset();
+
+  if (nrow  <= 0) {
+    Error("AliTRDdataArray::Allocate","The number of rows has to be positive");
+    return;
+  }
+  if (ncol  <= 0) {
+    Error("AliTRDdataArray::Allocate","The number of columns has to be positive");
+    return;
+  }
+  if (ntime <= 0) {
+    Error("AliTRDdataArray::Allocate","The number of timebins has to be positive");
+    return;
+  }
+
+  // The two-dimensional array row/column gets mapped into the first 
+  // dimension of the array. The second array dimension, which is not compressible,
+  // corresponds to the time direction
+  fNdim1  = nrow * ncol;
+  fNdim2  = ntime;
+  fNelems = fNdim1 * fNdim2;
+
+  fNrow   = nrow;
+  fNcol   = ncol;
+  fNtime  = ntime;
+
+  fElements->Set(fNelems);
+  fIndex->Set(fNdim2);
+
+  for (Int_t i = 0, k = 0; i < fNdim2; i++, k += fNdim1) { 
+    (*fIndex)[i] = k;
+  }
+
+  fBufType = 0;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDdataArray::GetSize()
+{
+  //
+  // Returns the size of the complete object
+  //
+
+  Int_t size = sizeof(this);
+
+  if (fIndex)    size += sizeof(fIndex)    
+                         + fIndex->GetSize()    * sizeof(Int_t);
+  if (fElements) size += sizeof(fElements) 
+                         + fElements->GetSize() * sizeof(Int_t);
+
+  return size;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDdataArray::GetDataSize()
+{
+  //
+  // Returns the size of only the data part
+  //
+
+  if (fElements == 0) 
+    return 0;
+  else 
+    return sizeof(fElements) + fElements->GetSize() * sizeof(Int_t);
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDdataArray::GetOverThreshold(Float_t threshold)
+{
+  //
+  // Returns the number of entries over threshold
+  //
+  if ((fElements == 0) || (fElements->GetSize() <= 0))
+    return 0;
+  Int_t  over = 0;
+
+  for (Bool_t cont = First(); cont == kTRUE; cont = Next()) {
+    if ((fCurrentIdx1 < 0) || (fCurrentIdx1 > fNdim1)) continue;
+    if ((fCurrentIdx2 < 0) || (fCurrentIdx2 > fNdim2)) continue;
+    if (fElements->At(fCurrentIndex) > threshold) over++;
+  }
+
+  return over;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDdataArray::GetData(Int_t row, Int_t col, Int_t time)
+{
+  //
+  // Returns the data value at a given position of the array
+  //
+
+  if (fBufType == 0) return GetDataFast(GetIdx1(row,col),time);
+  if (fBufType == 1) return GetData1(GetIdx1(row,col),time);
+
+  return 0;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdataArray::SetData(Int_t row, Int_t col, Int_t time, Int_t value)
+{
+  //
+  // Sets the data value at a given position of the array
+  //
+
+  SetDataFast(GetIdx1(row,col),time,value);
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdataArray::Expand()
+{  
+  //
+  // Expands the compressed buffer
+  //
+
+  if (fBufType  < 0) {
+    Error("AliTRDdataArray::Expand","Buffer does not exist");
+    return;
+  }
+  if (fBufType == 0) {  
+    return;
+  } 
+  // Expand a buffer of type 1
+  if (fBufType == 1) Expand1();
+  
+  fBufType = 0;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdataArray::Compress(Int_t bufferType)
+{
+  //
+  // Compresses the buffer
+  //
+
+  if (fBufType  < 0) {
+    Error("AliTRDdataArray::Compress","Buffer does not exist");
+    return;
+  }
+  if (fBufType == bufferType) {
+    return;
+  }  
+  if (fBufType > 0) {
+    Expand();
+  }
+  if (fBufType !=0)  {
+    Error("AliTRDdataArray::Compress","Buffer does not exist");
+    return;
+  }
+
+  // Compress a buffer of type 1
+  if (bufferType == 1) {
+    Compress1();
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdataArray::Compress(Int_t bufferType, Int_t threshold)
+{
+  //
+  // Compresses the buffer
+  //
+
+  fThreshold = threshold;
+  Compress(bufferType);
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdataArray::First()
+{
+  //
+  // Returns the position of the first valid data value
+  //
+
+  if (fBufType == 0) return First0();
+  if (fBufType == 1) return First1();
+  return kFALSE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t  AliTRDdataArray::Next()
+{
+  //
+  // Returns the position of the next valid data value
+  //
+
+  if (fBufType == 0) return Next0();
+  if (fBufType == 1) return Next1();
+  return kFALSE;
+
+}
+//_____________________________________________________________________________
+void AliTRDdataArray::Expand1()
+{
+  //
+  // Expands a buffer of type 1
+  //
+
+  Int_t i, k;
+
+  fNelems = fNdim1 * fNdim2;
+
+  Int_t *buf = new Int_t[fNelems];
+
+  fIndex->Set(fNdim2);
+
+  for (i = 0, k = 0; i < fNdim2; i++, k += fNdim1) (*fIndex)[i] = k;
+
+  Int_t idx1 = 0;
+  Int_t idx2 = 0;
+  Int_t N    = fElements->fN;
+
+  for (i = 0; i < N; i++){
+
+    // Negative sign counts the unwritten values (under threshold)
+    if ((*fElements)[i] < 0) {
+      idx1 -= fElements->At(i);
+    } 
+    else {
+      buf[(*fIndex)[idx2] + idx1] = fElements->At(i);
+      idx1++;
+    }
+    if (idx1 == fNdim1) {
+      idx1 = 0;
+      idx2++;
+    }
+    else { 
+      if (idx1 > fNdim1){
+       Reset();
+       return;
+      }      
+    }
+  }
+
+  fElements->Adopt(fNelems,buf); 
+   
+}
+
+//_____________________________________________________________________________
+void AliTRDdataArray::Compress1()
+{
+  //
+  // Compress a buffer of type 1
+  //
+
+  AliTRDarrayI  buf;  
+  buf.Set(fNelems);
+  AliTRDarrayI  index;
+  index.Set(fNdim2);
+
+  Int_t icurrent = -1;
+  Int_t izero;
+  for (Int_t idx2 = 0; idx2 < fNdim2; idx2++){      
+
+    // Set the idx2 pointer
+    index[idx2] = icurrent + 1;
+
+    // Reset the zero counter 
+    izero = 0;  
+
+    for (Int_t idx1 = 0; idx1 < fNdim1; idx1++){
+      // If below threshold
+      if (GetDataFast(idx1,idx2) <= fThreshold) {
+        izero++;
+      }
+      else {
+       if (izero > 0) {
+         // If we have currently izero counts under threshold
+         icurrent++;     
+         if (icurrent >= buf.fN) buf.Expand(icurrent*2);
+          // Store the number of entries below zero
+         buf[icurrent] = -izero;  
+         izero = 0;
+       } 
+       icurrent++;
+       if (icurrent >= buf.fN) buf.Expand(icurrent*2);
+       buf[icurrent] = GetDataFast(idx1,idx2);     
+      } // If signal larger than threshold             
+    } // End of loop over idx1
+
+    if (izero > 0) {
+      icurrent++;        
+      if (icurrent >= buf.fN) buf.Expand(icurrent*2);
+      // Store the number of entries below zero
+      buf[icurrent] = -izero; 
+    }
+
+  }
+
+  buf.Expand(icurrent+1);
+  (*fElements) = buf;
+  fNelems   = fElements->fN;
+  fBufType  = 1;
+  (*fIndex) = index;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdataArray::Expand2()
+{
+  //
+  // Expands a buffer of type 2 
+  //
+
+  Int_t i, k;
+  Int_t *buf = new Int_t[fNelems];
+
+  fNelems = fNdim1 * fNdim2;
+  fIndex->Set(fNdim2);
+
+  for (i = 0, k = 0; i < fNdim2; i++, k += fNdim1) (*fIndex)[i] = k;
+
+  Int_t idx1 = 0;
+  Int_t idx2 = 0;
+  Int_t N    = fElements->fN;
+  for (i = 0; i < N; i++){
+    // Negative sign counts the unwritten values (under threshold)
+    if ((*fElements)[i] < 0) {
+      idx1 -= fElements->At(i); 
+    }
+    else {
+      buf[(*fIndex)[idx2]+idx1] = fElements->At(i);
+      idx1++;
+    }
+    if (idx1 == fNdim1) {
+      idx1 = 0;
+      idx2++;
+    }
+    else { 
+      if (idx1 > fNdim1){
+       Reset();
+       return;
+      }      
+    }
+  }
+
+  fElements->Adopt(fNelems,buf);    
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdataArray::Compress2()
+{
+  /*
+
+  //
+  // Compress a buffer of type 2
+  //
+
+  AliArrayS  buf;  //lets have the nearly the "worst case"
+  buf.Set(fNelems);
+  AliTRDarrayI  index;
+  index.Set(fNdim2);
+  Int_t icurrent=-1;
+  Int_t izero;
+
+  for (Int_t col = 0; col<fNdim2; col++){      
+    index[col]=icurrent+1;//set collumn pointer
+    izero = 0;  //reset zer counter at the begining of the column
+    Int_t lastrow=0;
+    Int_t lastrowval=GetDigitFast(row,0);
+
+    for (Int_t row = 1; row< fNdim1;row++){
+      //if under threshold
+      Int_t val = GetDigitFast(row,col);
+      Int_t dif = val -lastrowval;
+      
+      if (TMath::Abs(dif)<fThreshold)  izero++;
+      else{
+       if (izero>0) {
+         //if we have currently izero count under threshold
+         icurrent++;     
+         if (icurrent>=buf.fN) buf.Expand(icurrent*2);
+         buf[icurrent]= -izero;  //write how many under zero
+         izero = 0;
+       } //end of reseting izero
+       icurrent++;
+       if (icurrent>=buf.fN) buf.Expand(icurrent*2);
+       buf[icurrent] = GetDigitFast(row,col);      
+      }//if signal bigger then threshold               
+    } //end of loop over rows
+    
+    if (izero>0) {
+      icurrent++;        
+      if (icurrent>=buf.fN) buf.Expand(icurrent*2);
+      buf[icurrent]= -izero;  //write how many under zero
+    }
+  }//end of lopping over digits
+  buf.Expand(icurrent+1);
+  (*fElements)=buf;
+  fNelems = fElements->fN;
+  fBufType = 1;
+  (*fIndex) =index;
+  //end of compresing bufer of type 1 
+
+  */
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdataArray::First0()
+{
+  //
+  // Returns the first entry for a buffer of type 0
+  //
+
+  fCurrentIdx1  = -1;
+  fCurrentIdx2  = -1;
+  fCurrentIndex = -1;
+
+  Int_t i;
+  for (i = 0; ((i < fNelems) && (fElements->At(i) <= fThreshold)); i++)
+  if (i == fNelems) return kFALSE;
+
+  fCurrentIdx1  = i % fNdim1;
+  fCurrentIdx2  = i / fNdim1;
+  fCurrentIndex = i;
+  return kTRUE;        
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdataArray::Next0()
+{
+  //
+  // Returns the next entry for a buffer of type 0
+  //
+
+  if (fCurrentIndex < 0) return kFALSE; 
+
+  Int_t i;
+  for (i = fCurrentIndex + 1; 
+       ((i < fNelems) && (fElements->At(i) <= fThreshold)); 
+       i++);
+  if (i >= fNelems)  {
+    fCurrentIndex = -1;
+    return kFALSE;
+  }
+
+  fCurrentIdx1  = i % fNdim1;
+  fCurrentIdx2  = i / fNdim1;
+  fCurrentIndex = i;
+  return kTRUE;        
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdataArray::First1()
+{
+  //
+  // Returns the first entry for a buffer of type 1
+  //
+
+  fCurrentIdx1  = -1;
+  fCurrentIdx2  =  0;
+  fCurrentIndex = -1;
+
+  Int_t i;
+  for (i = 0; i < fNelems; i++){
+    if (fElements->At(i) < 0) {
+      fCurrentIdx1-=fElements->At(i);
+    }
+    else {     
+      fCurrentIdx1++;
+    }
+    if (fCurrentIdx1 >= fNdim1) {
+      fCurrentIdx2++;
+      fCurrentIdx1 -= fNdim1;
+    }
+    if (fElements->At(i) > fThreshold) break;
+  }
+
+  fCurrentIndex = i;
+  if (fCurrentIndex >= 0) return kTRUE;
+  fCurrentIdx1  = -1;
+  fCurrentIdx2  = -1;
+  return kFALSE;       
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdataArray::Next1()
+{
+  //
+  // Returns the next entry for a buffer of type 1
+  //
+
+  if (fCurrentIndex < 0) return kFALSE;
+
+  Int_t i;
+  for (i = fCurrentIndex + 1; i < fNelems; i++){
+    if (fElements->At(i) < 0) {
+      fCurrentIdx1 -= fElements->At(i);
+    }
+    else {      
+      fCurrentIdx1++;
+    }
+    if (fCurrentIdx1 >= fNdim1) {
+      fCurrentIdx2++;
+      fCurrentIdx1 -= fNdim1;
+    }
+    if (fElements->At(i) > fThreshold) break;
+  }
+
+  fCurrentIndex =  i;
+  if ((i >= 0) && (i < fNelems)) return kTRUE;
+  fCurrentIdx1  = -1;
+  fCurrentIdx2  = -1;
+  return kFALSE;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDdataArray::GetData1(Int_t idx1, Int_t idx2)
+{
+  //
+  // Returns the value at a given position of the array
+  //
+  
+  Int_t i, n2;
+
+  if ((idx2 + 1) >= fNdim2) {
+    n2 = fNelems;
+  }
+  else {
+    n2 = fIndex->At(idx2 + 1);
+  }
+
+  // Current idx1    
+  Int_t curidx1 = 0; 
+  for (i = fIndex->At(idx2); ((i < n2) && (curidx1 < idx1)); i++){
+    if (fElements->At(i) < 0) {
+      curidx1 -= fElements->At(i);
+    }
+    else {      
+      curidx1++;
+    }
+  }
+
+  if ((curidx1 == idx1) && (fElements->At(i) > 0)) {
+    return fElements->At(i);
+  }
+  else {
+    return 0;
+  }
+
+}
+
diff --git a/TRD/AliTRDdataArray.h b/TRD/AliTRDdataArray.h
new file mode 100644 (file)
index 0000000..8cce169
--- /dev/null
@@ -0,0 +1,174 @@
+#ifndef TRDdataArray_H
+#define TRDdataArray_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+#include   "AliTRDarrayI.h"
+#include   "AliTRDsegmentID.h"
+
+/////////////////////////////////////////////////////////////
+//  General container for data from TRD detector segments  //
+//  Adapted from AliDigits, origin M.Ivanov                //
+/////////////////////////////////////////////////////////////
+
+class AliTRDdataArray : public AliTRDsegmentID {
+
+ public:
+
+  AliTRDdataArray();
+  ~AliTRDdataArray();
+
+  virtual void   Allocate(Int_t nrow, Int_t ncol,Int_t ntime);  
+  virtual void   Compress(Int_t bufferType, Int_t threshold);
+  virtual void   Compress(Int_t bufferType); 
+  virtual void   Expand();
+  virtual void   Reset();  
+  virtual Bool_t First();
+  virtual Bool_t Next(); 
+
+  virtual void   SetData(Int_t row, Int_t col, Int_t time, Int_t value);
+
+  virtual void   SetThreshold(Int_t threshold) { fThreshold = threshold; };
+
+  virtual Int_t  GetThreshold()                { return fThreshold;  };
+  virtual Int_t  GetNRow()                     { return fNrow;       };
+  virtual Int_t  GetNCol()                     { return fNcol;       };
+  virtual Int_t  GetNtime()                    { return fNtime;      };
+
+  virtual Int_t  GetSize();
+  virtual Int_t  GetDataSize(); 
+  virtual Int_t  GetOverThreshold(Float_t threshold);  
+
+  virtual Int_t  GetData(Int_t row, Int_t col, Int_t time);
+  inline  Int_t  GetIndex(Int_t row, Int_t col, Int_t time);
+
+ protected:
+
+  inline  Int_t  GetIdx1(Int_t row, Int_t col);
+  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);
+  inline  void   SetDataFast(Int_t idx1, Int_t idx2, Int_t value); 
+  inline  Int_t  GetDataFast(Int_t idx1, Int_t idx2); 
+
+  Int_t          GetData1(Int_t idx1, Int_t idx2); 
+  void           Expand1(); 
+  void           Compress1(); 
+  void           Expand2();
+  void           Compress2();
+  Bool_t         First0();
+  Bool_t         Next0(); 
+  Bool_t         First1();
+  Bool_t         Next1();
+  Int_t          fNrow;            // Number of rows of the detector segement
+  Int_t          fNcol;            // Number of columns of the detector segment
+  Int_t          fNtime;           // Number of timebins of the detector segment
+
+  Int_t          fNdim1;           // First dimension of the array (row * column)
+  Int_t          fNdim2;           // Second dimension of the array (time, not compressed) 
+ private:
+
+  AliTRDarrayI  *fElements;        // Buffer of 4 bytes integers for the array content
+  AliTRDarrayI  *fIndex;           // Index position of column
+  Int_t          fBufType;         // Type of the buffer - defines the compression algorithm  
+  Int_t          fThreshold;       // Treshold for zero suppression
+  Int_t          fNelems;          // Total number of elements 
+  Int_t          fCurrentIdx1;     // !Current index 1
+  Int_t          fCurrentIdx2;     // !Current index 2
+  Int_t          fCurrentIndex;    // !Current index in field
+  ClassDef(AliTRDdataArray,1)      // Data container for one TRD detector segment
+
+};
+//_____________________________________________________________________________
+inline Int_t AliTRDdataArray::GetIndex(Int_t row, Int_t col, Int_t time)
+{
+  //
+  // Maps the row/column/time into one number
+  // 
+
+  return time * fNrow*fNcol + GetIdx1(row,col);
+
+}
+
+//_____________________________________________________________________________
+inline Int_t AliTRDdataArray::GetIdx1(Int_t row, Int_t col)
+{
+  //
+  // Maps the two-dimensional row/column plane into an one-dimensional array.
+  //
+
+  return row * fNrow + col;
+
+}
+
+//_____________________________________________________________________________
+inline Bool_t AliTRDdataArray::CheckBounds(const char *where
+                                          , Int_t idx1, Int_t idx2) 
+{
+  //
+  // Does the boundary checking
+  //
+
+  if ((idx2 >= fNdim2) || (idx2 < 0)) 
+    return OutOfBoundsError(where,idx1,idx2);
+
+  Int_t index = (*fIndex).At(idx2) + idx1;
+  if ((index < 0) || (index > fNelems)) 
+    return OutOfBoundsError(where,idx1,idx2);
+
+  return kTRUE;  
+
+}
+
+//_____________________________________________________________________________
+inline Bool_t AliTRDdataArray::OutOfBoundsError(const char *where
+                                               , Int_t idx1, Int_t idx2) 
+{
+  //
+  // Generate an out-of-bounds error. Always returns false.
+  //
+
+  ::Error(where, "idx1 %d  idx2 %d out of bounds (size: %d x %d, this: 0x%08x)"
+        ,idx1,idx2,fNdim1,fNdim2,this);
+
+  return kFALSE;
+
+}
+
+//_____________________________________________________________________________
+inline Int_t AliTRDdataArray::GetDataFast(Int_t idx1, Int_t idx2)
+{
+  //
+  // Returns the value at a given position in the array
+  //
+
+  return fElements->At(fIndex->At(idx2) + idx1); 
+
+}
+
+//_____________________________________________________________________________
+inline void  AliTRDdataArray::SetDataFast(Int_t idx1, Int_t idx2, Int_t value)
+{
+  //
+  // Set the value at a given position in the array
+  //
+
+  if ((idx1 < 0) || (idx1 >= fNdim1) || 
+      (idx2 < 0) || (idx2 >= fNdim2)) { 
+    ::Error("AliTRDdataArray::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; 
+
+}
+
+#endif
+
diff --git a/TRD/AliTRDdetDigits.cxx b/TRD/AliTRDdetDigits.cxx
new file mode 100644 (file)
index 0000000..6e13b29
--- /dev/null
@@ -0,0 +1,70 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Container class for the TRD digits of one detector segment (chamber).    //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDdetDigits.h"
+
+ClassImp(AliTRDdetDigits)
+
+//_____________________________________________________________________________
+AliTRDdetDigits::AliTRDdetDigits():AliDigits()
+{
+  //
+  // Default constructor
+  //
+
+  fNrowTRD = 0;
+  fNcolTRD = 0;
+
+}
+
+//_____________________________________________________________________________
+AliTRDdetDigits::Allocate(Int_t nrow, Int_t ncol, Int_t ntime)
+{
+  //
+  // Allocates an empty buffer for the digits data of the size
+  // <nrow> * <ncol> * <ntime>
+  //
+
+  fNrowTRD = nrow;
+  fNcolTRD = ncol;
+  
+  // The two-dimensional row/column structure of the TRD gets mapped into 
+  // an one-dimensional row array which can be used inside AliDigits.
+  // The TRD timebins correspond to the columns in AliDigits.
+  AliDigits::Allocate(fNrowTRD*fNcolTRD,ntime);
+
+}
+
+//_____________________________________________________________________________
+AliTRDdetDigits::SetDigits(Int_t row, Int_t col, Int_t time, Short_t value)
+{
+  //
+  // Sets the value of one given digit
+  //
+
+  Int_t index = row * fNrowTRD + col;
+  AliDigits::SetDigitFast(value,index,time);
+
+}
diff --git a/TRD/AliTRDdetDigits.h b/TRD/AliTRDdetDigits.h
new file mode 100644 (file)
index 0000000..6c32a59
--- /dev/null
@@ -0,0 +1,34 @@
+#ifndef TRDdetDigits_H
+#define TRDdetDigits_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+/////////////////////////////////////////////
+//  Digits container for one TRD detector  //
+/////////////////////////////////////////////
+
+#include "AliDigits.h"
+
+//_____________________________________________________________________________
+class AliTRDdetDigits : public AliDigits {
+
+ public:
+
+  AliTRDdetDigits();
+  ~AliTRDdetDigits() { };
+
+  virtual void Allocate(Int_t nrow, Int_t ncol, Int_t ntime);
+  virtual void SetDigit(Int_t row, Int_t col, Int_t time, Short_t value);
+
+ protected:
+
+  Int_t   fNrowTRD;                  // Number of rows in the TRD
+  Int_t   fNcolTRD;                  // Number of colums in the TRD
+
+  ClassDef(AliTRDdetDigits,1)        // Digits container for one TRD detector
+
+};
+
+#endif
diff --git a/TRD/AliTRDdigitizer.cxx b/TRD/AliTRDdigitizer.cxx
new file mode 100644 (file)
index 0000000..0c30826
--- /dev/null
@@ -0,0 +1,696 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Creates and handles digits from TRD hits                                 //
+//                                                                           //
+//  The following effects are included:                                      //
+//      - Diffusion                                                          //
+//      - ExB effects                                                        //
+//      - Gas gain including fluctuations                                    //
+//      - Pad-response (simple Gaussian approximation)                       //
+//      - Electronics noise                                                  //
+//      - Electronics gain                                                   //
+//      - Digitization                                                       //
+//      - ADC threshold                                                      //
+//  The corresponding parameter can be adjusted via the various              //
+//  Set-functions. If these parameters are not explicitly set, default       //
+//  values are used (see Init-function).                                     //
+//  To produce digits from a root-file with TRD-hits use the                 //
+//  slowDigitsCreate.C macro.                                                //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TMath.h>
+#include <TVector.h>
+#include <TRandom.h>
+
+#include "AliTRD.h"
+#include "AliTRDdigitizer.h"
+#include "AliTRDmatrix.h"
+
+ClassImp(AliTRDdigitizer)
+
+//_____________________________________________________________________________
+AliTRDdigitizer::AliTRDdigitizer():TNamed()
+{
+  //
+  // AliTRDdigitizer default constructor
+  //
+
+  fInputFile = NULL;
+  fEvent     = 0;
+
+}
+
+//_____________________________________________________________________________
+AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
+                :TNamed(name,title)
+{
+  //
+  // AliTRDdigitizer default constructor
+  //
+
+  fInputFile   = NULL;
+  fEvent       = 0;
+
+  fDigitsArray   = new AliTRDsegmentArray(kNsect*kNplan*kNcham);
+  for (Int_t iDict = 0; iDict < kNDict; iDict++) {
+    fDictionary[iDict] = new AliTRDsegmentArray(kNsect*kNplan*kNcham);
+  }
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDdigitizer::~AliTRDdigitizer()
+{
+
+  if (fInputFile) {
+    fInputFile->Close();
+    delete fInputFile;
+  }
+
+  if (fDigitsArray) {
+    fDigitsArray->Delete();
+    delete fDigitsArray;
+  }
+
+  for (Int_t iDict = 0; iDict < kNDict; iDict++) {
+    fDictionary[iDict]->Delete();
+    delete fDictionary[iDict];
+  }
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDdigitizer::Diffusion(Float_t driftlength, Float_t *xyz)
+{
+  //
+  // Applies the diffusion smearing to the position of a single electron
+  //
+
+  Float_t driftSqrt = TMath::Sqrt(driftlength);
+  Float_t sigmaT = driftSqrt * fDiffusionT;
+  Float_t sigmaL = driftSqrt * fDiffusionL;
+  xyz[0] = gRandom->Gaus(xyz[0], sigmaL * fLorentzFactor);
+  xyz[1] = gRandom->Gaus(xyz[1], sigmaT * fLorentzFactor);
+  xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
+  return 1;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDdigitizer::ExB(Float_t driftlength, Float_t *xyz)
+{
+  //
+  // Applies E x B effects to the position of a single electron
+  //
+
+  xyz[0] = xyz[0];
+  xyz[1] = xyz[1] + fLorentzAngle * driftlength;
+  xyz[2] = xyz[2];
+
+  return 1;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdigitizer::Init()
+{
+  //
+  // Initializes the digitization procedure with standard values
+  //
+
+  // The default parameter for the digitization
+  fGasGain       = 2.0E3;
+  fNoise         = 3000.;
+  fChipGain      = 10.;
+  fADCoutRange   = 255.;
+  fADCinRange    = 2000.;
+  fADCthreshold  = 1;
+
+  // Transverse and longitudinal diffusion coefficients (Xe/Isobutane)
+  fDiffusionOn   = 1;
+  fDiffusionT    = 0.060;
+  fDiffusionL    = 0.017;
+
+  // Propability for electron attachment
+  fElAttachOn    = 0;
+  fElAttachProp  = 0.0;
+
+  // E x B effects
+  fExBOn         = 0;
+  // omega * tau. (tau ~ 12 * 10^-12, B = 0.2T)
+  fLorentzAngle  = 17.6 * 12.0 * 0.2 * 0.01;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::Open(const Char_t *name, Int_t nEvent)
+{
+  //
+  // Opens a ROOT-file with TRD-hits and reads in the hit-tree
+  //
+
+  // Connect the AliRoot file containing Geometry, Kine, and Hits
+  fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
+  if (!fInputFile) {
+    printf("AliTRDdigitizer::Open -- ");
+    printf("Open the ALIROOT-file %s.\n",name);
+    fInputFile = new TFile(name,"UPDATE");
+  }
+  else {
+    printf("AliTRDdigitizer::Open -- ");
+    printf("%s is already open.\n",name);
+  }
+
+  // Get AliRun object from file or create it if not on file
+  //if (!gAlice) {
+    gAlice = (AliRun*) fInputFile->Get("gAlice");
+    if (gAlice) {
+      printf("AliTRDdigitizer::Open -- ");
+      printf("AliRun object found on file.\n");
+    }
+    else {
+      printf("AliTRDdigitizer::Open -- ");
+      printf("Could not find AliRun object.\n");
+      return kFALSE;
+    }
+  //}
+
+  fEvent = nEvent;
+
+  // Import the Trees for the event nEvent in the file
+  Int_t nparticles = gAlice->GetEvent(fEvent);
+  if (nparticles <= 0) {
+    printf("AliTRDdigitizer::Open -- ");
+    printf("No entries in the trees for event %d.\n",fEvent);
+    return kFALSE;
+  }
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDdigitizer::PadResponse(Float_t x)
+{
+  //
+  // The pad response for the chevron pads. 
+  // We use a simple Gaussian approximation which should be good
+  // enough for our purpose.
+  //
+
+  // The parameters for the response function
+  const Float_t aa  =  0.8872;
+  const Float_t bb  = -0.00573;
+  const Float_t cc  =  0.454;
+  const Float_t cc2 =  cc*cc;
+
+  Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2)));
+
+  return (pr);
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::MakeDigits()
+{
+  //
+  // Loops through the TRD-hits and creates the digits.
+  //
+
+  // Get the pointer to the detector class and check for version 1
+  AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
+  if (TRD->IsVersion() != 1) {
+    printf("AliTRDdigitizer::MakeDigits -- ");
+    printf("TRD must be version 1 (slow simulator).\n");
+    return kFALSE; 
+  }
+
+  // Get the geometry
+  AliTRDgeometry *Geo = TRD->GetGeometry();
+  printf("AliTRDdigitizer::MakeDigits -- ");
+  printf("Geometry version %d\n",Geo->IsVersion());
+
+  printf("AliTRDdigitizer::MakeDigits -- ");
+  printf("Start creating digits.\n");
+
+  ///////////////////////////////////////////////////////////////
+  // Parameter 
+  ///////////////////////////////////////////////////////////////
+
+  // Converts number of electrons to fC
+  const Float_t el2fC  = 1.602E-19 * 1.0E15; 
+
+  ///////////////////////////////////////////////////////////////
+
+  Int_t   iRow, iCol, iTime;
+  Int_t   nBytes = 0;
+
+  Int_t   totalSizeDigits = 0;
+  Int_t   totalSizeDict0  = 0;
+  Int_t   totalSizeDict1  = 0;
+  Int_t   totalSizeDict2  = 0;
+
+  AliTRDhit       *Hit;
+  AliTRDdataArray *Digits;
+  AliTRDdataArray *Dictionary[kNDict];
+
+  // Get the pointer to the hit tree
+  TTree *HitTree    = gAlice->TreeH();
+
+  // The Lorentz factor
+  if (fExBOn) {
+    fLorentzFactor = 1.0 / (1.0 + fLorentzAngle*fLorentzAngle);
+  }
+  else {
+    fLorentzFactor = 1.0;
+  }
+
+  // Get the number of entries in the hit tree
+  // (Number of primary particles creating a hit somewhere)
+  Int_t nTrack = (Int_t) HitTree->GetEntries();
+
+  Int_t chamBeg = 0;
+  Int_t chamEnd = kNcham;
+  if (TRD->GetSensChamber() >= 0) {
+    chamBeg = TRD->GetSensChamber();
+    chamEnd = chamEnd + 1;
+  }
+  Int_t planBeg = 0;
+  Int_t planEnd = kNplan;
+  if (TRD->GetSensPlane()   >= 0) {
+    planBeg = TRD->GetSensPlane();
+    planEnd = planBeg + 1;
+  }
+  Int_t sectBeg = 0;
+  Int_t sectEnd = kNsect;
+  if (TRD->GetSensSector()  >= 0) {
+    sectBeg = TRD->GetSensSector();
+    sectEnd = sectBeg + 1;
+  }
+
+  // Loop through all the chambers
+  for (Int_t iCham = chamBeg; iCham < chamEnd; iCham++) {
+    for (Int_t iPlan = planBeg; iPlan < planEnd; iPlan++) {
+      for (Int_t iSect = sectBeg; iSect < sectEnd; iSect++) {
+
+        Int_t nDigits = 0;
+
+        Int_t iDet = Geo->GetDetector(iPlan,iCham,iSect);
+
+        printf("AliTRDdigitizer::MakeDigits -- ");
+        printf("Digitizing chamber %d, plane %d, sector %d.\n"
+              ,iCham,iPlan,iSect);
+
+        Int_t   nRowMax     = Geo->GetRowMax(iPlan,iCham,iSect);
+        Int_t   nColMax     = Geo->GetColMax(iPlan);
+        Int_t   nTimeMax    = Geo->GetTimeMax();
+        Float_t row0        = Geo->GetRow0(iPlan,iCham,iSect);
+        Float_t col0        = Geo->GetCol0(iPlan);
+        Float_t time0       = Geo->GetTime0(iPlan);
+        Float_t rowPadSize  = Geo->GetRowPadSize();
+        Float_t colPadSize  = Geo->GetColPadSize();
+        Float_t timeBinSize = Geo->GetTimeBinSize();
+
+        // Create a detector matrix to keep the signal and track numbers
+        AliTRDmatrix *Matrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
+                                               ,iSect,iCham,iPlan);
+
+        // Loop through all entries in the tree
+        for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
+
+          gAlice->ResetHits();
+          nBytes += HitTree->GetEvent(iTrack);
+
+          // Get the number of hits in the TRD created by this particle
+          Int_t nHit = TRD->Hits()->GetEntriesFast();
+
+          // Loop through the TRD hits  
+          for (Int_t iHit = 0; iHit < nHit; iHit++) {
+
+            if (!(Hit = (AliTRDhit *) TRD->Hits()->UncheckedAt(iHit))) 
+              continue;
+
+            Float_t pos[3];
+                    pos[0]   = Hit->fX;
+                    pos[1]   = Hit->fY;
+                    pos[2]   = Hit->fZ;
+            Float_t q        = Hit->fQ;
+            Int_t   track    = Hit->fTrack;
+            Int_t   detector = Hit->fDetector;
+            Int_t   plane    = Geo->GetPlane(detector);
+            Int_t   sector   = Geo->GetSector(detector);
+            Int_t   chamber  = Geo->GetChamber(detector);
+
+            if ((sector  != iSect) ||
+                (plane   != iPlan) ||
+                (chamber != iCham)) 
+              continue;
+
+            // Rotate the sectors on top of each other
+            Float_t rot[3];
+            Geo->Rotate(detector,pos,rot);
+
+            // The hit position in pad coordinates (center pad)
+            // The pad row (z-direction)
+            Int_t  rowH = (Int_t) ((rot[2] -  row0) /  rowPadSize);
+            // The pad column (rphi-direction)  
+            Int_t  colH = (Int_t) ((rot[1] -  col0) /  colPadSize);
+            // The time bucket
+            Int_t timeH = (Int_t) ((rot[0] - time0) / timeBinSize);
+
+            // Array to sum up the signal in a box surrounding the
+            // hit postition
+            const Int_t timeBox = 7;
+            const Int_t  colBox = 9;
+            const Int_t  rowBox = 7;
+            Float_t signalSum[rowBox][colBox][timeBox];
+            for (iRow  = 0;  iRow <  rowBox; iRow++ ) {
+              for (iCol  = 0;  iCol <  colBox; iCol++ ) {
+                for (iTime = 0; iTime < timeBox; iTime++) {
+                  signalSum[iRow][iCol][iTime] = 0;
+               }
+             }
+           }
+
+            // Loop over all electrons of this hit
+            Int_t nEl = (Int_t) q;
+            for (Int_t iEl = 0; iEl < nEl; iEl++) {
+
+             // The driftlength
+              Float_t driftlength = rot[0] - time0;
+              if ((driftlength <        0) || 
+                  (driftlength > kDrThick)) break;
+              Float_t driftlengthL = driftlength;
+              if (fExBOn) driftlengthL /= TMath::Sqrt(fLorentzFactor);
+              Float_t xyz[3];
+              xyz[0] = rot[0];
+              xyz[1] = rot[1];
+              xyz[2] = rot[2];
+
+              // Electron attachment
+              if (fElAttachOn) {
+                if (gRandom->Rndm() < (driftlengthL * fElAttachProp / 100.)) continue;
+             }
+
+              // Apply the diffusion smearing
+              if (fDiffusionOn) {
+                if (!(Diffusion(driftlengthL,xyz))) continue;
+             }
+
+              // Apply E x B effects
+              if (fExBOn) { 
+                if (!(ExB(driftlength,xyz))) continue;   
+             }
+
+              // The electron position and the distance to the hit position
+             // in pad units
+              // The pad row (z-direction)
+              Int_t  rowE = (Int_t) ((xyz[2] -  row0) /  rowPadSize);
+              Int_t  rowD =  rowH -  rowE;
+              // The pad column (rphi-direction)
+              Int_t  colE = (Int_t) ((xyz[1] -  col0) /  colPadSize);
+              Int_t  colD =  colH -  colE;
+              // The time bucket
+              Int_t timeE = (Int_t) ((xyz[0] - time0) / timeBinSize);
+              Int_t timeD = timeH - timeE;
+
+              // Apply the gas gain including fluctuations
+              Int_t signal = (Int_t) (-fGasGain * TMath::Log(gRandom->Rndm()));
+
+             // The distance of the electron to the center of the pad 
+             // in units of pad width
+              Float_t dist = (xyz[1] - col0 - (colE + 0.5) * colPadSize) 
+                           / colPadSize;
+
+              // Sum up the signal in the different pixels
+              // and apply the pad response
+              Int_t  rowIdx =  rowD + (Int_t) ( rowBox / 2);
+              Int_t  colIdx =  colD + (Int_t) ( colBox / 2);
+              Int_t timeIdx = timeD + (Int_t) (timeBox / 2);
+
+              if (( rowIdx < 0) || ( rowIdx >  rowBox)) {
+                printf("AliTRDdigitizer::MakeDigits -- ");
+                printf("Boundary error. rowIdx = %d (%d)\n", rowIdx, rowBox);
+                continue;
+             }
+              if (( colIdx < 0) || ( colIdx >  colBox)) {
+                printf("AliTRDdigitizer::MakeDigits -- ");
+                printf("Boundary error. colIdx = %d (%d)\n", colIdx, colBox);
+                continue;
+             }
+              if ((timeIdx < 0) || (timeIdx > timeBox)) {
+                printf("AliTRDdigitizer::MakeDigits -- ");
+                printf("Boundary error. timeIdx = %d (%d)\n",timeIdx,timeBox);
+                continue;
+             }
+              signalSum[rowIdx][colIdx-1][timeIdx] += PadResponse(dist-1.) * signal;
+              signalSum[rowIdx][colIdx  ][timeIdx] += PadResponse(dist   ) * signal;
+              signalSum[rowIdx][colIdx+1][timeIdx] += PadResponse(dist+1.) * signal;
+
+            }
+
+            // Add the padcluster to the detector matrix
+            for (iRow  = 0;  iRow <  rowBox; iRow++ ) {
+              for (iCol  = 0;  iCol <  colBox; iCol++ ) {
+                for (iTime = 0; iTime < timeBox; iTime++) {
+
+                  Int_t  rowB =  rowH + iRow  - (Int_t) ( rowBox / 2); 
+                  Int_t  colB =  colH + iCol  - (Int_t) ( colBox / 2);
+                  Int_t timeB = timeH + iTime - (Int_t) (timeBox / 2);
+
+                  Float_t signalB = signalSum[iRow][iCol][iTime];
+                  if (signalB > 0.0) {
+                    Matrix->AddSignal(rowB,colB,timeB,signalB);
+                    if (!(Matrix->AddTrack(rowB,colB,timeB,track))) { 
+                      printf("AliTRDdigitizer::MakeDigits -- ");
+                      printf("More than three tracks in a pixel!\n");
+                   }
+                 }
+
+               }
+             }
+           }
+
+          }
+
+       }
+
+        // Add a container for the digits of this detector
+        Digits = (AliTRDdataArray *) fDigitsArray->At(iDet);        
+        // Allocate memory space for the digits buffer
+        Digits->Allocate(nRowMax,nColMax,nTimeMax);
+
+        for (Int_t iDict = 0; iDict < kNDict; iDict++) {
+          Dictionary[iDict] = (AliTRDdataArray *) fDictionary[iDict]->At(iDet);
+          Dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeMax);
+       }
+
+        // Create the hits for this chamber
+        for (iRow  = 0; iRow  <  nRowMax; iRow++ ) {
+          for (iCol  = 0; iCol  <  nColMax; iCol++ ) {
+            for (iTime = 0; iTime < nTimeMax; iTime++) {         
+
+              Float_t signalAmp = Matrix->GetSignal(iRow,iCol,iTime);
+
+              // Add the noise
+              signalAmp  = TMath::Max(gRandom->Gaus(signalAmp,fNoise),(Float_t) 0.0);
+             // Convert to fC
+              signalAmp *= el2fC;
+              // Convert to mV
+              signalAmp *= fChipGain;
+             // Convert to ADC counts
+              Int_t adc  = (Int_t) (signalAmp * (fADCoutRange / fADCinRange));
+
+              // Store the amplitude of the digit
+              Digits->SetData(iRow,iCol,iTime,adc);
+
+              // Store the track index in the dictionary
+              // Note: We store index+1 in order to allow the array to be compressed
+              for (Int_t iDict = 0; iDict < kNDict; iDict++) {
+                Dictionary[iDict]->SetData(iRow,iCol,iTime
+                                          ,Matrix->GetTrack(iRow,iCol,iTime,iDict)+1);
+             }
+
+              if (adc > fADCthreshold) nDigits++;
+
+           }
+         }
+       }
+
+        // Compress the arrays
+        Digits->Compress(1,fADCthreshold);
+        for (Int_t iDict = 0; iDict < kNDict; iDict++) {
+          Dictionary[iDict]->Compress(1,0);
+       }
+
+        totalSizeDigits += Digits->GetSize();
+        totalSizeDict0  += Dictionary[0]->GetSize();
+        totalSizeDict1  += Dictionary[1]->GetSize();
+        totalSizeDict2  += Dictionary[2]->GetSize();
+
+        printf("AliTRDdigitizer::MakeDigits -- ");
+        printf("Number of digits found: %d.\n",nDigits);
+
+       // Clean up
+        if (Matrix) delete Matrix;
+
+      }
+    }
+  }
+
+  printf("AliTRDdigitizer::MakeDigits -- ");
+  printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
+                                                    ,totalSizeDict0
+                                                    ,totalSizeDict1
+                                                    ,totalSizeDict2);        
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::MakeBranch()
+{
+  //
+  // Creates the branches for the digits and the dictionary
+  //
+
+  Int_t buffersize = 64000;
+
+  Bool_t status = kTRUE;
+
+  if (gAlice->TreeD()) {
+
+    // Make the branch for the digits
+    if (fDigitsArray) {
+      const AliTRDdataArray *Digits = 
+           (AliTRDdataArray *) fDigitsArray->At(0);
+      if (Digits) {
+        gAlice->TreeD()->Branch("TRDdigits",Digits->IsA()->GetName()
+                                           ,&Digits,buffersize,1);
+        printf("AliTRDdigitizer::MakeBranch -- ");
+        printf("Making branch TRDdigits\n");
+      }
+      else {
+        status = kFALSE;
+      }
+    }
+    else {
+      status = kFALSE;
+    }
+
+    // Make the branches for the dictionaries
+    for (Int_t iDict = 0; iDict < kNDict; iDict++) {
+
+      Char_t branchname[15];
+      sprintf(branchname,"TRDdictionary%d",iDict);
+      if (fDictionary[iDict]) {
+        const AliTRDdataArray *Dictionary = 
+             (AliTRDdataArray *) fDictionary[iDict]->At(0);
+        if (Dictionary) {
+          gAlice->TreeD()->Branch(branchname,Dictionary->IsA()->GetName()
+                                            ,&Dictionary,buffersize,1);
+          printf("AliTRDdigitizer::MakeBranch -- ");
+          printf("Making branch %s\n",branchname);
+       }
+        else {
+          status = kFALSE;
+       }
+      }
+      else {
+        status = kFALSE;
+      }
+    }
+
+  }
+  else {
+    status = kFALSE;
+  }
+
+  return status;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::WriteDigits()
+{
+  //
+  // Writes out the TRD-digits and the dictionaries
+  //
+
+  // Create the branches
+  if (!(gAlice->TreeD()->GetBranch("TRDdigits"))) { 
+    if (!MakeBranch()) return kFALSE;
+  }
+
+  // Store the contents of the segment array in the tree
+  if (!fDigitsArray->StoreArray("TRDdigits")) {
+    printf("AliTRDdigitizer::WriteDigits -- ");
+    printf("Error while storing digits in branch TRDdigits\n");
+    return kFALSE;
+  }
+  for (Int_t iDict = 0; iDict < kNDict; iDict++) {
+    Char_t branchname[15];
+    sprintf(branchname,"TRDdictionary%d",iDict);
+    if (!fDictionary[iDict]->StoreArray(branchname)) {
+      printf("AliTRDdigitizer::WriteDigits -- ");
+      printf("Error while storing dictionary in branch %s\n",branchname);
+      return kFALSE;
+    }
+  }
+
+  // Write the new tree into the input file (use overwrite option)
+  Char_t treeName[7];
+  sprintf(treeName,"TreeD%d",fEvent);
+  printf("AliTRDdigitizer::WriteDigits -- ");
+  printf("Write the digits tree %s for event %d.\n"
+        ,treeName,fEvent);
+  gAlice->TreeD()->Write(treeName,2);
+  return kTRUE;
+
+}
+
+ClassImp(AliTRDdigit)
+
+//_____________________________________________________________________________
+AliTRDdigit::AliTRDdigit(Int_t *digits):AliDigitNew()
+{
+  //
+  // Create a TRD digit
+  //
+
+  // Store the volume hierarchy
+  fDetector  = digits[0];
+
+  // Store the row, pad, and time bucket number
+  fRow       = digits[1];
+  fCol       = digits[2];
+  fTime      = digits[3];
+
+  // Store the signal amplitude
+  fAmplitude = digits[4];
+
+}
diff --git a/TRD/AliTRDdigitizer.h b/TRD/AliTRDdigitizer.h
new file mode 100644 (file)
index 0000000..6289a54
--- /dev/null
@@ -0,0 +1,126 @@
+#ifndef TRDdigitizer_h
+#define TRDdigitizer_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <TNamed.h>
+#include <TFile.h>
+
+#include "AliHit.h" 
+#include "AliDigit.h"
+#include "AliTRDconst.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDsegmentArray.h"
+
+///////////////////////////////////////////////////////
+//  Produces digits from the hits information        //
+///////////////////////////////////////////////////////
+
+const Int_t kNDict = 3;
+
+class AliTRDdigitizer : public TNamed {
+
+ public:
+
+  AliTRDdigitizer();
+  AliTRDdigitizer(const Text_t* name, const Text_t* title);
+  ~AliTRDdigitizer();
+
+  virtual void        Init();
+  virtual Bool_t      Open(const Char_t *name, Int_t nEvent = 0);
+  virtual Bool_t      MakeBranch();
+  virtual Bool_t      MakeDigits();
+  virtual Bool_t      WriteDigits();
+
+  virtual void        SetGasGain(Float_t gasgain)     { fGasGain       = gasgain;  };
+  virtual void        SetNoise(Float_t noise)         { fNoise         = noise;    };
+  virtual void        SetChipGain(Float_t chipgain)   { fChipGain      = chipgain; };
+  virtual void        SetADCoutRange(Float_t range)   { fADCoutRange   = range;    };
+  virtual void        SetADCinRange(Float_t range)    { fADCinRange    = range;    };
+  virtual void        SetADCthreshold(Int_t thresh)   { fADCthreshold  = thresh;   };
+  virtual void        SetDiffusion(Int_t diff_on = 1) { fDiffusionOn   = diff_on;  };
+  virtual void        SetDiffusionT(Float_t diff)     { fDiffusionT    = diff;     };
+  virtual void        SetDiffusionL(Float_t diff)     { fDiffusionL    = diff;     };
+  virtual void        SetElAttach(Int_t el_on = 1)    { fElAttachOn    = el_on;    };
+  virtual void        SetElAttachProp(Float_t prop)   { fElAttachProp  = prop;     };
+  virtual void        SetExB(Int_t exb_on = 1)        { fExBOn         = exb_on;   };
+  virtual void        SetLorentzAngle(Float_t angle)  { fLorentzAngle  = angle;    };
+
+  AliTRDsegmentArray *DigitsArray()                   { return fDigitsArray;   };
+  AliTRDsegmentArray *Dictionary(Int_t i)             { return fDictionary[i]; };
+
+  virtual Float_t     GetGasGain()                    { return fGasGain;       };
+  virtual Float_t     GetNoise()                      { return fNoise;         };
+  virtual Float_t     GetChipGain()                   { return fChipGain;      };
+  virtual Float_t     GetADCoutRange()                { return fADCoutRange;   };
+  virtual Float_t     GetADCinRange()                 { return fADCinRange;    };
+  virtual Int_t       GetADCthreshold()               { return fADCthreshold;  };
+  virtual Float_t     GetDiffusionT()                 { return fDiffusionT;    };
+  virtual Float_t     GetDiffusionL()                 { return fDiffusionL;    };
+  virtual Float_t     GetElAttachProp()               { return fElAttachProp;  };
+  virtual Float_t     GetLorentzAngle()               { return fLorentzAngle;  };
+
+ protected:
+
+  TFile              *fInputFile;       //! ALIROOT-filename
+
+  AliTRDsegmentArray *fDigitsArray;     //! Array of detector segments containing the digits
+  AliTRDsegmentArray *fDictionary[3];   //! Dictionary array, connecting MC tracks to the digits
+  
+  Int_t               fEvent;           //! Event number
+
+  Float_t             fGasGain;         // Gas gain
+  Float_t             fNoise;           // Electronics noise
+  Float_t             fChipGain;        // Electronics gain
+  Float_t             fADCoutRange;     // ADC output range (number of channels)
+  Float_t             fADCinRange;      // ADC input range (input charge)
+  Int_t               fADCthreshold;    // ADC threshold in ADC channel
+  Int_t               fDiffusionOn;     // Switch for the diffusion
+  Float_t             fDiffusionT;      // Diffusion in transverse direction
+  Float_t             fDiffusionL;      // Diffusion in longitudinal direction
+  Int_t               fElAttachOn;      // Switch for the electron attachment
+  Float_t             fElAttachProp;    // Propability for electron attachment (for 1m)
+  Int_t               fExBOn;           // Switch for the ExB effects
+  Float_t             fLorentzAngle;    // Lorentz angle 
+  Float_t             fLorentzFactor;   // Factor due to Lorentz force
+
+ private:
+
+  virtual Int_t       Diffusion(Float_t driftlength, Float_t *xyz);
+  virtual Int_t       ExB(Float_t driftlength, Float_t *xyz);  
+  virtual Float_t     PadResponse(Float_t x);
+  
+  ClassDef(AliTRDdigitizer,1)           // TRD-Digits manager
+
+};
+
+//_____________________________________________________________________________
+class AliTRDdigit : public AliDigitNew {
+
+ public:
+
+  AliTRDdigit() {};
+  AliTRDdigit(Int_t *digits);
+  virtual ~AliTRDdigit() {};
+
+  virtual Int_t GetAmp()      { return fAmplitude; };
+  virtual Int_t GetDetector() { return fDetector;  };
+  virtual Int_t GetRow()      { return fRow;       };
+  virtual Int_t GetCol()      { return fCol;       };
+  virtual Int_t GetTime()     { return fTime;      };
+
+ protected:
+
+  Int_t        fDetector;   // TRD detector number
+  Int_t        fRow;        // Pad row number
+  Int_t        fCol;        // Pad col number
+  Int_t        fTime;       // Time bucket
+  Int_t        fAmplitude;  // Signal amplitude
+
+  ClassDef(AliTRDdigit,1)   // Digits for Transition Radiation Detector
+
+};
+
+#endif
diff --git a/TRD/AliTRDgeometry.cxx b/TRD/AliTRDgeometry.cxx
new file mode 100644 (file)
index 0000000..6240059
--- /dev/null
@@ -0,0 +1,464 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  TRD geometry class                                                       //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDgeometry.h"
+#include "AliTRDrecPoint.h"
+
+ClassImp(AliTRDgeometry)
+
+//_____________________________________________________________________________
+AliTRDgeometry::AliTRDgeometry():AliGeometry()
+{
+  //
+  // AliTRDgeometry default constructor
+  //
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDgeometry::~AliTRDgeometry()
+{
+
+}
+
+//_____________________________________________________________________________
+void AliTRDgeometry::Init()
+{
+  //
+  // Initializes the geometry parameter
+  //
+
+  Int_t iplan;
+
+  // The width of the chambers
+  fCwidth[0]    =  99.6;
+  fCwidth[1]    = 104.1;
+  fCwidth[2]    = 108.5;
+  fCwidth[3]    = 112.9;
+  fCwidth[4]    = 117.4;
+  fCwidth[5]    = 121.8;
+
+  // The default pad dimensions
+  fRowPadSize  = 4.5;
+  fColPadSize  = 1.0;
+  fTimeBinSize = 0.1;
+
+  // The maximum number of pads
+  // and the position of pad 0,0,0 
+  // 
+  // chambers seen from the top:
+  //     +----------------------------+
+  //     |                            |
+  //     |                            |     ^
+  //     |                            | rphi|
+  //     |                            |     |
+  //     |0                           |     | 
+  //     +----------------------------+     +------>
+  //                                             z 
+  // chambers seen from the side:           ^
+  //     +----------------------------+ time|
+  //     |                            |     |
+  //     |0                           |     |
+  //     +----------------------------+     +------>
+  //                                             z
+  //                                             
+
+  // The pad column (rphi-direction)  
+  for (iplan = 0; iplan < kNplan; iplan++) {
+    fColMax[iplan] = 1 + TMath::Nint((fCwidth[iplan] - 2. * kCcthick) 
+                                                     / fColPadSize - 0.5);
+    fCol0[iplan]   = -fCwidth[iplan]/2. + kCcthick;
+  }
+
+  // The time bucket
+  fTimeMax = 1 + TMath::Nint(kDrThick / fTimeBinSize - 0.5);
+  for (iplan = 0; iplan < kNplan; iplan++) {
+    fTime0[iplan]  = kRmin + kCcframe/2. + kDrZpos - 0.5 * kDrThick
+                           + iplan * (kCheight + kCspace);
+  } 
+
+}
+
+//_____________________________________________________________________________
+void AliTRDgeometry::CreateGeometry(Int_t *idtmed)
+{
+  //
+  // Create the TRD geometry
+  //
+  // Author: Christoph Blume (C.Blume@gsi.de) 20/07/99
+  //
+  // The volumes:
+  //    TRD1-3     (Air)   --- The TRD mother volumes for one sector. 
+  //                           To be placed into the spaceframe.
+  //
+  //    UAFI(/M/O) (Al)    --- The aluminum frame of the inner(/middle/outer) chambers (readout)
+  //    UCFI(/M/O) (C)     --- The carbon frame of the inner(/middle/outer) chambers 
+  //                           (driftchamber + radiator)
+  //    UAII(/M/O) (Air)   --- The inner part of the readout of the inner(/middle/outer) chambers
+  //    UFII(/M/O) (Air)   --- The inner part of the chamner and radiator of the 
+  //                           inner(/middle/outer) chambers
+  //
+  // The material layers in one chamber:
+  //    UL01       (G10)   --- The gas seal of the radiator
+  //    UL02       (CO2)   --- The gas in the radiator
+  //    UL03       (PE)    --- The foil stack
+  //    UL04       (Mylar) --- Entrance window to the driftvolume and HV-cathode
+  //    UL05       (Xe)    --- The driftvolume
+  //    UL06       (Xe)    --- The amplification region
+  //    
+  //    UL07       (Cu)    --- The pad plane
+  //    UL08       (G10)   --- The Nomex honeycomb support structure
+  //    UL09       (Cu)    --- FEE and signal lines
+  //    UL10       (PE)    --- The cooling devices
+  //    UL11       (Water) --- The cooling water
+
+  const Int_t npar_cha = 3;
+
+  Float_t par_dum[3];
+  Float_t par_cha[npar_cha];
+
+  Float_t xpos, ypos, zpos;
+
+  // The aluminum frames - readout + electronics (Al)
+  // The inner chambers
+  gMC->Gsvolu("UAFI","BOX ",idtmed[1301-1],par_dum,0);
+  // The middle chambers
+  gMC->Gsvolu("UAFM","BOX ",idtmed[1301-1],par_dum,0);
+  // The outer chambers
+  gMC->Gsvolu("UAFO","BOX ",idtmed[1301-1],par_dum,0);
+
+  // The inner part of the aluminum frames (Air)
+  // The inner chambers
+  gMC->Gsvolu("UAII","BOX ",idtmed[1302-1],par_dum,0);
+  // The middle chambers
+  gMC->Gsvolu("UAIM","BOX ",idtmed[1302-1],par_dum,0);
+  // The outer chambers
+  gMC->Gsvolu("UAIO","BOX ",idtmed[1302-1],par_dum,0);
+
+  // The carbon frames - radiator + driftchamber (C)
+  // The inner chambers
+  gMC->Gsvolu("UCFI","BOX ",idtmed[1307-1],par_dum,0);
+  // The middle chambers
+  gMC->Gsvolu("UCFM","BOX ",idtmed[1307-1],par_dum,0);
+  // The outer chambers
+  gMC->Gsvolu("UCFO","BOX ",idtmed[1307-1],par_dum,0);
+
+  // The inner part of the carbon frames (Air)
+  // The inner chambers
+  gMC->Gsvolu("UCII","BOX ",idtmed[1302-1],par_dum,0);
+  // The middle chambers
+  gMC->Gsvolu("UCIM","BOX ",idtmed[1302-1],par_dum,0);
+  // The outer chambers
+  gMC->Gsvolu("UCIO","BOX ",idtmed[1302-1],par_dum,0);
+
+  // The material layers inside the chambers
+  par_cha[0] = -1.;
+  par_cha[1] = -1.;
+  // G10 layer (radiator seal)
+  par_cha[2] = kSeThick/2;
+  gMC->Gsvolu("UL01","BOX ",idtmed[1313-1],par_cha,npar_cha);
+  // CO2 layer (radiator)
+  par_cha[2] = kRaThick/2;
+  gMC->Gsvolu("UL02","BOX ",idtmed[1312-1],par_cha,npar_cha);
+  // PE layer (radiator)
+  par_cha[2] = kPeThick/2;
+  gMC->Gsvolu("UL03","BOX ",idtmed[1303-1],par_cha,npar_cha);
+  // Mylar layer (entrance window + HV cathode) 
+  par_cha[2] = kMyThick/2;
+  gMC->Gsvolu("UL04","BOX ",idtmed[1308-1],par_cha,npar_cha);
+  // Xe/Isobutane layer (drift volume, sensitive) 
+  par_cha[2] = kDrThick/2.;
+  gMC->Gsvolu("UL05","BOX ",idtmed[1309-1],par_cha,npar_cha);
+  // Xe/Isobutane layer (amplification volume, not sensitive)
+  par_cha[2] = kAmThick/2.;
+  gMC->Gsvolu("UL06","BOX ",idtmed[1309-1],par_cha,npar_cha);
+  
+  // Cu layer (pad plane)
+  par_cha[2] = kCuThick/2;
+  gMC->Gsvolu("UL07","BOX ",idtmed[1305-1],par_cha,npar_cha);
+  // G10 layer (support structure)
+  par_cha[2] = kSuThick/2;
+  gMC->Gsvolu("UL08","BOX ",idtmed[1313-1],par_cha,npar_cha);
+  // Cu layer (FEE + signal lines)
+  par_cha[2] = kFeThick/2;
+  gMC->Gsvolu("UL09","BOX ",idtmed[1305-1],par_cha,npar_cha);
+  // PE layer (cooling devices)
+  par_cha[2] = kCoThick/2;
+  gMC->Gsvolu("UL10","BOX ",idtmed[1303-1],par_cha,npar_cha);
+  // Water layer (cooling)
+  par_cha[2] = kWaThick/2;
+  gMC->Gsvolu("UL11","BOX ",idtmed[1314-1],par_cha,npar_cha);
+
+  // Position the layers in the chambers
+  xpos = 0;
+  ypos = 0;
+
+  // G10 layer (radiator seal)
+  zpos = kSeZpos;
+  gMC->Gspos("UL01",1,"UCII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL01",2,"UCIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL01",3,"UCIO",xpos,ypos,zpos,0,"ONLY");
+  // CO2 layer (radiator)
+  zpos = kRaZpos;
+  gMC->Gspos("UL02",1,"UCII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL02",2,"UCIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL02",3,"UCIO",xpos,ypos,zpos,0,"ONLY");
+  // PE layer (radiator)
+  zpos = 0;
+  gMC->Gspos("UL03",1,"UL02",xpos,ypos,zpos,0,"ONLY");
+  // Mylar layer (entrance window + HV cathode)   
+  zpos = kMyZpos;
+  gMC->Gspos("UL04",1,"UCII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL04",2,"UCIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL04",3,"UCIO",xpos,ypos,zpos,0,"ONLY");
+  // Xe/Isobutane layer (drift volume) 
+  zpos = kDrZpos;
+  gMC->Gspos("UL05",1,"UCII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL05",2,"UCIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL05",3,"UCIO",xpos,ypos,zpos,0,"ONLY");
+  // Xe/Isobutane layer (amplification volume)
+  zpos = kAmZpos;
+  gMC->Gspos("UL06",1,"UCII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL06",2,"UCIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL06",3,"UCIO",xpos,ypos,zpos,0,"ONLY");
+
+  // Cu layer (pad plane)
+  zpos = kCuZpos;
+  gMC->Gspos("UL07",1,"UAII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL07",2,"UAIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL07",3,"UAIO",xpos,ypos,zpos,0,"ONLY");
+  // G10 layer (support structure)
+  zpos = kSuZpos;
+  gMC->Gspos("UL08",1,"UAII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL08",2,"UAIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL08",3,"UAIO",xpos,ypos,zpos,0,"ONLY");
+  // Cu layer (FEE + signal lines)
+  zpos = kFeZpos; 
+  gMC->Gspos("UL09",1,"UAII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL09",2,"UAIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL09",3,"UAIO",xpos,ypos,zpos,0,"ONLY");
+  // PE layer (cooling devices)
+  zpos = kCoZpos;
+  gMC->Gspos("UL10",1,"UAII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL10",2,"UAIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL10",3,"UAIO",xpos,ypos,zpos,0,"ONLY");
+  // Water layer (cooling)
+  zpos = kWaZpos;
+  gMC->Gspos("UL11",1,"UAII",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL11",1,"UAIM",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("UL11",1,"UAIO",xpos,ypos,zpos,0,"ONLY");
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDgeometry::Local2Global(Int_t idet, Float_t *local, Float_t *global)
+{
+  //
+  // Converts local pad-coordinates (row,col,time) into 
+  // global ALICE reference frame coordinates (x,y,z)
+  //
+
+  Int_t        icham     = GetChamber(idet);    // Chamber info (0-4)
+  Int_t        isect     = GetSector(idet);     // Sector info  (0-17)
+  Int_t        iplan     = GetPlane(idet);      // Plane info   (0-5)
+
+  Float_t      padRow    = local[0];            // Pad Row position
+  Float_t      padCol    = local[1];            // Pad Column position
+  Float_t      timeSlice = local[2];            // Time "position"
+
+  Float_t      row0      = GetRow0(iplan,icham,isect);
+  Float_t      col0      = GetCol0(iplan);
+  Float_t      time0     = GetTime0(iplan);
+
+  Float_t      rot[3];
+
+  // calculate (x,y) position in rotated chamber
+  rot[1] = col0  + padCol    * fColPadSize;
+  rot[0] = time0 + timeSlice * fTimeBinSize;
+  // calculate z-position:
+  rot[2] = row0  + padRow    * fRowPadSize;
+
+  // Rotate back to original position
+  return RotateBack(idet,rot,global);
+
+}
+//_____________________________________________________________________________
+Bool_t AliTRDgeometry::Local2Global(Int_t iplan, Int_t icham, Int_t isect
+                                  , Float_t *local, Float_t *global)
+{
+  //
+  // Converts local pad-coordinates (row,col,time) into 
+  // global ALICE reference frame coordinates (x,y,z)
+  //
+
+  Int_t        idet      = GetDetector(iplan,icham,isect); // Detector number
+
+  Float_t      padRow    = local[0];                       // Pad Row position
+  Float_t      padCol    = local[1];                       // Pad Column position
+  Float_t      timeSlice = local[2];                       // Time "position"
+
+  Float_t      row0      = GetRow0(iplan,icham,isect);
+  Float_t      col0      = GetCol0(iplan);
+  Float_t      time0     = GetTime0(iplan);
+
+  Float_t      rot[3];
+
+  // calculate (x,y,z) position in rotated chamber
+  rot[1] = col0  + padCol    * fColPadSize;
+  rot[0] = time0 + timeSlice * fTimeBinSize;
+  rot[2] = row0  + padRow    * fRowPadSize;
+
+  // Rotate back to original position
+  return RotateBack(idet,rot,global);
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDgeometry::Rotate(Int_t d, Float_t *pos, Float_t *rot)
+{
+  //
+  // Rotates all chambers in the position of sector 0 and transforms
+  // the coordinates in the ALICE restframe <pos> into the 
+  // corresponding local frame <rot>.
+  //
+
+  Int_t   sector = GetSector(d);
+
+  Float_t phi    = -2.0 * kPI /  (Float_t) kNsect * ((Float_t) sector + 0.5);
+
+  rot[0] =  pos[0] * TMath::Cos(phi) + pos[1] * TMath::Sin(phi);
+  rot[1] = -pos[0] * TMath::Sin(phi) + pos[1] * TMath::Cos(phi);
+  rot[2] =  pos[2];
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDgeometry::RotateBack(Int_t d, Float_t *rot, Float_t *pos)
+{
+  //
+  // Rotates a chambers from the position of sector 0 into its
+  // original position and transforms the corresponding local frame 
+  // coordinates <rot> into the coordinates of the ALICE restframe <pos>.
+  //
+
+  Int_t   sector = GetSector(d);
+
+  Float_t phi    =  2.0 * kPI /  (Float_t) kNsect * ((Float_t) sector + 0.5);
+
+  rot[0] =  pos[0] * TMath::Cos(phi) + pos[1] * TMath::Sin(phi);
+  rot[1] = -pos[0] * TMath::Sin(phi) + pos[1] * TMath::Cos(phi);
+  rot[2] =  pos[2];
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDgeometry::GetDetector(Int_t p, Int_t c, Int_t s)
+{
+  //
+  // Convert plane / chamber / sector into detector number
+  //
+
+  return (p + c * kNplan + s * kNplan * kNcham);
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDgeometry::GetPlane(Int_t d)
+{
+  //
+  // Reconstruct the plane number from the detector number
+  //
+
+  return ((Int_t) (d % kNplan));
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDgeometry::GetChamber(Int_t d)
+{
+  //
+  // Reconstruct the chamber number from the detector number
+  //
+
+  return ((Int_t) (d % (kNplan * kNcham)) / kNplan);
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDgeometry::GetSector(Int_t d)
+{
+  //
+  // Reconstruct the sector number from the detector number
+  //
+
+  return ((Int_t) (d / (kNplan * kNcham)));
+
+}
+
+//_____________________________________________________________________________
+void AliTRDgeometry::GetGlobal(const AliRecPoint *p, TVector3 &pos, TMatrix &mat)
+{
+  // 
+  // Returns the global coordinate and error matrix of a AliTRDrecPoint
+  //
+
+  GetGlobal(p,pos);
+
+}
+
+//_____________________________________________________________________________
+void AliTRDgeometry::GetGlobal(const AliRecPoint *p, TVector3 &pos)
+{
+  // 
+  // Returns the global coordinate and error matrix of a AliTRDrecPoint
+  //
+
+  Int_t detector = ((AliTRDrecPoint *) p)->GetDetector();
+
+  Float_t global[3];
+  Float_t local[3];
+  local[0] = pos.X();
+  local[1] = pos.Y();
+  local[2] = pos.Z();
+
+  if (Local2Global(detector,local,global)) {
+    pos.SetX(global[0]);
+    pos.SetY(global[1]);
+    pos.SetZ(global[2]);
+  }
+  else {
+    pos.SetX(0.0);
+    pos.SetY(0.0);
+    pos.SetZ(0.0);
+  }
+
+}
diff --git a/TRD/AliTRDgeometry.h b/TRD/AliTRDgeometry.h
new file mode 100644 (file)
index 0000000..0cb221c
--- /dev/null
@@ -0,0 +1,75 @@
+#ifndef TRDgeometry_h
+#define TRDgeometry_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <TObject.h>
+#include <TMath.h>
+
+#include "AliRun.h"
+#include "AliRecPoint.h"
+
+#include "AliTRDconst.h"
+
+class AliTRDgeometry : public AliGeometry {
+
+ public:
+
+  AliTRDgeometry();
+  ~AliTRDgeometry();
+
+  virtual void    CreateGeometry(Int_t *);
+  virtual Int_t   IsVersion() const = 0;
+  virtual void    Init();
+  virtual Bool_t  Local2Global(Int_t d, Float_t *local, Float_t *global);
+  virtual Bool_t  Local2Global(Int_t p, Int_t c, Int_t s, Float_t *local, Float_t *global);
+  virtual Bool_t  Rotate(Int_t d, Float_t *pos, Float_t *rot);
+  virtual Bool_t  RotateBack(Int_t d, Float_t *rot, Float_t *pos);
+
+  virtual void    SetRowPadSize(Float_t size)          { fRowPadSize  = size; };
+  virtual void    SetColPadSize(Float_t size)          { fColPadSize  = size; };
+  virtual void    SetTimeBinSize(Float_t size)         { fTimeBinSize = size; };
+
+  virtual Int_t   GetDetector(Int_t p, Int_t c, Int_t s);
+  virtual Int_t   GetPlane(Int_t d);
+  virtual Int_t   GetChamber(Int_t d);
+  virtual Int_t   GetSector(Int_t d);
+
+  virtual Int_t   GetRowMax(Int_t p, Int_t c, Int_t s) { return fRowMax[p][c][s]; };
+  virtual Int_t   GetColMax(Int_t p)                   { return fColMax[p];       };
+  virtual Int_t   GetTimeMax()                         { return fTimeMax;         };
+  virtual Float_t GetRow0(Int_t p, Int_t c, Int_t s)   { return fRow0[p][c][s]; };
+  virtual Float_t GetCol0(Int_t p)                     { return fCol0[p];       };
+  virtual Float_t GetTime0(Int_t p)                    { return fTime0[p];      };
+
+  virtual Float_t GetRowPadSize()                      { return fRowPadSize;  };
+  virtual Float_t GetColPadSize()                      { return fColPadSize;  };
+  virtual Float_t GetTimeBinSize()                     { return fTimeBinSize; };
+
+  virtual void    GetGlobal(const AliRecPoint * p, TVector3 & pos, TMatrix & mat); 
+  virtual void    GetGlobal(const AliRecPoint * p, TVector3 & pos);   
+
+ protected:
+
+  Float_t         fCwidth[kNplan];                 // Width of the chambers
+
+  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
+
+  Float_t         fRow0[kNplan][kNcham][kNsect];   // Row-position of pad 0
+  Float_t         fCol0[kNplan];                   // Column-position of pad 0
+  Float_t         fTime0[kNplan];                  // Time-position of pad 0
+
+  Float_t         fRowPadSize;                     // Pad size in z-direction
+  Float_t         fColPadSize;                     // Pad size in rphi-direction
+  Float_t         fTimeBinSize;                    // Size of the time buckets
+
+  ClassDef(AliTRDgeometry,1)                       // TRD geometry base class
+
+};
+
+#endif
diff --git a/TRD/AliTRDgeometryFull.cxx b/TRD/AliTRDgeometryFull.cxx
new file mode 100644 (file)
index 0000000..c3d4fe2
--- /dev/null
@@ -0,0 +1,284 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  TRD geometry without holes                                               //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDgeometryFull.h"
+
+ClassImp(AliTRDgeometryFull)
+
+//_____________________________________________________________________________
+AliTRDgeometryFull::AliTRDgeometryFull():AliTRDgeometry()
+{
+  //
+  // AliTRDgeometryFull default constructor
+  //
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDgeometryFull::~AliTRDgeometryFull()
+{
+
+}
+
+//_____________________________________________________________________________
+void AliTRDgeometryFull::Init()
+{
+  //
+  // Initializes the geometry parameter
+  //
+
+  Int_t iplan;
+
+  // The length of the inner chambers
+  for (iplan = 0; iplan < kNplan; iplan++) 
+    fClengthI[iplan] = 110.0;
+
+  // The length of the middle chambers
+  fClengthM[0] = 123.5;
+  fClengthM[1] = 131.0;
+  fClengthM[2] = 138.5;
+  fClengthM[3] = 146.0;
+  fClengthM[4] = 153.0;
+  fClengthM[5] = 160.5;
+
+  // The length of the outer chambers
+  fClengthO[0] = 123.5;
+  fClengthO[1] = 131.0;
+  fClengthO[2] = 134.5;
+  fClengthO[3] = 142.0;
+  fClengthO[4] = 142.0;
+  fClengthO[5] = 134.5;
+
+  // The maximum number of pads
+  // and the position of pad 0,0,0 
+  // 
+  // chambers seen from the top:
+  //     +----------------------------+
+  //     |                            |
+  //     |                            |     ^
+  //     |                            | rphi|
+  //     |                            |     |
+  //     |0                           |     | 
+  //     +----------------------------+     +------>
+  //                                             z 
+  // chambers seen from the side:           ^
+  //     +----------------------------+ time|
+  //     |                            |     |
+  //     |0                           |     |
+  //     +----------------------------+     +------>
+  //                                             z
+  //                                             
+
+  // The pad row (z-direction)
+  for (iplan = 0; iplan < kNplan; iplan++) {
+
+    for (Int_t isect = 0; isect < kNsect; isect++) {
+      Float_t clengthI = fClengthI[iplan];
+      Float_t clengthM = fClengthM[iplan];
+      Float_t clengthO = fClengthO[iplan];
+      fRowMax[iplan][0][isect] = 1 + TMath::Nint((clengthO - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRowMax[iplan][1][isect] = 1 + TMath::Nint((clengthM - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRowMax[iplan][2][isect] = 1 + TMath::Nint((clengthI - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRowMax[iplan][3][isect] = 1 + TMath::Nint((clengthM - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRowMax[iplan][4][isect] = 1 + TMath::Nint((clengthO - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRow0[iplan][0][isect]   = -clengthI/2. - clengthM - clengthO + kCcthick; 
+      fRow0[iplan][1][isect]   = -clengthI/2. - clengthM            + kCcthick;
+      fRow0[iplan][2][isect]   = -clengthI/2.                       + kCcthick;
+      fRow0[iplan][3][isect]   =  clengthI/2.                       + kCcthick; 
+      fRow0[iplan][4][isect]   =  clengthI/2. + clengthM            + kCcthick; 
+    }
+
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDgeometryFull::CreateGeometry(Int_t *idtmed)
+{
+  //
+  // Create the TRD geometry without hole
+  //
+
+  Int_t iplan;
+
+  const Int_t npar_trd = 4;
+  const Int_t npar_cha = 3;
+
+  Float_t par_trd[npar_trd];
+  Float_t par_cha[npar_cha];
+
+  Float_t xpos, ypos, zpos;
+
+  AliTRDgeometry::CreateGeometry(idtmed);
+
+  // The TRD mother volume for one sector (Air) (dimensions identical to BTR1)
+  par_trd[0] = kSwidth1/2.;
+  par_trd[1] = kSwidth2/2.;
+  par_trd[2] = kSlenTR1/2.;
+  par_trd[3] = kSheight/2.;
+  gMC->Gsvolu("TRD1","TRD1",idtmed[1302-1],par_trd,npar_trd);
+  
+  // Position the chambers in the TRD mother volume
+  for (iplan = 1; iplan <= kNplan; iplan++) {
+
+    // The inner chambers ---------------------------------------------------------------
+
+    // the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthI[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = 0.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFI",iplan       ,"TRD1",xpos,ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.   - kCathick;
+    par_cha[1] = fClengthI[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = 0.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAII",iplan       ,"TRD1",xpos,ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthI[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = 0.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFI",iplan       ,"TRD1",xpos,ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.   - kCcthick;
+    par_cha[1] = fClengthI[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = 0.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCII",iplan       ,"TRD1",xpos,ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // The middle chambers --------------------------------------------------------------
+
+    // the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthM[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2.  + fClengthM[iplan-1]/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFM",iplan         ,"TRD1",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    gMC->Gsposp("UAFM",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCathick;
+    par_cha[1] = fClengthM[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2.  + fClengthM[iplan-1]/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAIM",iplan         ,"TRD1",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    gMC->Gsposp("UAIM",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthM[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM[iplan-1]/2.;
+    zpos       = kCcframe/2.           - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFM",iplan         ,"TRD1",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    gMC->Gsposp("UCFM",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCcthick;
+    par_cha[1] = fClengthM[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM[iplan-1]/2.;
+    zpos       = kCcframe/2.           - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCIM",iplan         ,"TRD1",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    gMC->Gsposp("UCIM",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // The outer chambers ---------------------------------------------------------------
+
+    // the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthO[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM[iplan-1]    + fClengthO[iplan-1]/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFO",iplan         ,"TRD1",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    gMC->Gsposp("UAFO",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCathick;
+    par_cha[1] = fClengthO[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM[iplan-1]    + fClengthO[iplan-1]/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAIO",iplan         ,"TRD1",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    gMC->Gsposp("UAIO",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthO[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM[iplan-1]    + fClengthO[iplan-1]/2.;
+    zpos       = kCcframe/2.           - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFO",iplan,         "TRD1",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    gMC->Gsposp("UCFO",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCcthick;
+    par_cha[1] = fClengthO[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM[iplan-1]    + fClengthO[iplan-1]/2.;
+    zpos       = kCcframe/2.           - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCIO",iplan         ,"TRD1",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    gMC->Gsposp("UCIO",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+  }
+
+  xpos     = 0.;
+  ypos     = 0.;
+  zpos     = 0.;
+  gMC->Gspos("TRD1",1,"BTR1",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("TRD1",2,"BTR2",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("TRD1",3,"BTR3",xpos,ypos,zpos,0,"ONLY");
+
+}
diff --git a/TRD/AliTRDgeometryFull.h b/TRD/AliTRDgeometryFull.h
new file mode 100644 (file)
index 0000000..96cea6a
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef TRDgeometryFull_h
+#define TRDgeometryFull_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliTRDgeometry.h"
+
+class AliTRDgeometryFull : public AliTRDgeometry {
+
+ public:
+
+  AliTRDgeometryFull();
+  ~AliTRDgeometryFull();
+
+  virtual void    CreateGeometry(Int_t *);
+  virtual Int_t   IsVersion() const { return 1; };
+  virtual void    Init();
+
+ protected:
+
+  Float_t         fClengthI[kNplan];               // Length of the inner chambers
+  Float_t         fClengthM[kNplan];               // Length of the middle chambers
+  Float_t         fClengthO[kNplan];               // Length of the outer chambers
+  Float_t         fCwidth[kNplan];                 // Width of the chambers
+
+  ClassDef(AliTRDgeometryFull,1)                   // TRD geometry without hole
+
+};
+
+#endif
diff --git a/TRD/AliTRDgeometryHole.cxx b/TRD/AliTRDgeometryHole.cxx
new file mode 100644 (file)
index 0000000..464bd08
--- /dev/null
@@ -0,0 +1,418 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  TRD geometry with holes                                                  //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDgeometryHole.h"
+
+ClassImp(AliTRDgeometryHole)
+
+//_____________________________________________________________________________
+AliTRDgeometryHole::AliTRDgeometryHole():AliTRDgeometry()
+{
+  //
+  // AliTRDgeometryHole default constructor
+  //
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDgeometryHole::~AliTRDgeometryHole()
+{
+
+}
+
+//_____________________________________________________________________________
+void AliTRDgeometryHole::Init()
+{
+  //
+  // Initializes the geometry parameter
+  //
+
+  Int_t iplan;
+
+  // The length of the inner chambers
+  for (iplan = 0; iplan < kNplan; iplan++) 
+    fClengthI[iplan] = 110.0;
+
+  // The length of the middle chambers
+  fClengthM1[0] = 123.5;
+  fClengthM1[1] = 131.0;
+  fClengthM1[2] = 138.5;
+  fClengthM1[3] = 146.0;
+  fClengthM1[4] = 153.0;
+  fClengthM1[5] = 160.5;
+
+  fClengthM2[0] = 123.5 - 7.0;
+  fClengthM2[1] = 131.0 - 7.0;
+  fClengthM2[2] = 138.5 - 7.0;
+  fClengthM2[3] = 146.0 - 7.0;
+  fClengthM2[4] = 153.0 - 7.0;
+  fClengthM2[5] = 160.4 - 7.0;
+
+  // The length of the outer chambers
+  fClengthO1[0] = 123.5;
+  fClengthO1[1] = 131.0;
+  fClengthO1[2] = 134.5;
+  fClengthO1[3] = 142.0;
+  fClengthO1[4] = 142.0;
+  fClengthO1[5] = 134.5;
+
+  fClengthO2[0] = 123.5;
+  fClengthO2[1] = 131.0;
+  fClengthO2[2] = 134.5;
+  fClengthO2[3] = 142.0;
+  fClengthO2[4] = 142.0;
+  fClengthO2[5] = 134.5;
+
+  fClengthO3[0] =  86.5;
+  fClengthO3[1] = 101.5;
+  fClengthO3[2] = 112.5;
+  fClengthO3[3] = 127.5;
+  fClengthO3[4] = 134.5;
+  fClengthO3[5] = 134.5;
+
+  // The maximum number of pads
+  // and the position of pad 0,0,0 
+  // 
+  // chambers seen from the top:
+  //     +----------------------------+
+  //     |                            |
+  //     |                            |     ^
+  //     |                            | rphi|
+  //     |                            |     |
+  //     |0                           |     | 
+  //     +----------------------------+     +------>
+  //                                             z 
+  // chambers seen from the side:           ^
+  //     +----------------------------+ time|
+  //     |                            |     |
+  //     |0                           |     |
+  //     +----------------------------+     +------>
+  //                                             z
+  //                                             
+
+  // The pad row (z-direction)
+  for (iplan = 0; iplan < kNplan; iplan++) {
+
+    for (Int_t isect = 0; isect < kNsect; isect++) {
+      Float_t clengthI = fClengthI[iplan];
+      Float_t clengthM = fClengthM1[iplan];
+      Float_t clengthO = fClengthO1[iplan];
+      switch (isect) {
+      case 12:
+      case 13:
+      case 14:
+      case 15:
+      case 16:
+        clengthM = fClengthM2[iplan];
+        clengthO = fClengthO2[iplan];
+        break;
+      case 4:
+      case 5:
+      case 6:
+        clengthO = fClengthO3[iplan];
+        break;
+      };
+      fRowMax[iplan][0][isect] = 1 + TMath::Nint((clengthO - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRowMax[iplan][1][isect] = 1 + TMath::Nint((clengthM - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRowMax[iplan][2][isect] = 1 + TMath::Nint((clengthI - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRowMax[iplan][3][isect] = 1 + TMath::Nint((clengthM - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRowMax[iplan][4][isect] = 1 + TMath::Nint((clengthO - 2. * kCcthick) 
+                                                           / fRowPadSize - 0.5);
+      fRow0[iplan][0][isect]   = -clengthI/2. - clengthM - clengthO + kCcthick; 
+      fRow0[iplan][1][isect]   = -clengthI/2. - clengthM            + kCcthick;
+      fRow0[iplan][2][isect]   = -clengthI/2.                       + kCcthick;
+      fRow0[iplan][3][isect]   =  clengthI/2.                       + kCcthick; 
+      fRow0[iplan][4][isect]   =  clengthI/2. + clengthM            + kCcthick; 
+    }
+
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDgeometryHole::CreateGeometry(Int_t *idtmed)
+{
+  //
+  // Create the TRD geometry with hole
+  //
+
+  Int_t iplan;
+
+  const Int_t npar_trd = 4;
+  const Int_t npar_cha = 3;
+
+  Float_t par_trd[npar_trd];
+  Float_t par_cha[npar_cha];
+
+  Float_t xpos, ypos, zpos;
+
+  AliTRDgeometry::CreateGeometry(idtmed);
+
+  // The TRD mother volume for one sector (Air) (dimensions identical to BTR1)
+  par_trd[0] = kSwidth1/2.;
+  par_trd[1] = kSwidth2/2.;
+  par_trd[2] = kSlenTR1/2.;
+  par_trd[3] = kSheight/2.;
+  gMC->Gsvolu("TRD1","TRD1",idtmed[1302-1],par_trd,npar_trd);
+  
+  // The TRD mother volume for one sector (Air) (dimensions identical to BTR2)
+  par_trd[0] = kSwidth1/2.;
+  par_trd[1] = kSwidth2/2.;
+  par_trd[2] = kSlenTR2/2.;
+  par_trd[3] = kSheight/2.;
+  gMC->Gsvolu("TRD2","TRD1",idtmed[1302-1],par_trd,npar_trd);
+
+  // The TRD mother volume for one sector (Air) (dimensions identical to BTR3)
+  par_trd[0] = kSwidth1/2.;
+  par_trd[1] = kSwidth2/2.;
+  par_trd[2] = kSlenTR3/2.;
+  par_trd[3] = kSheight/2.;
+  gMC->Gsvolu("TRD3","TRD1",idtmed[1302-1],par_trd,npar_trd);
+
+  // Position the chambers in the TRD mother volume
+  for (iplan = 1; iplan <= kNplan; iplan++) {
+
+    // The inner chambers ---------------------------------------------------------------
+
+    // the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthI[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = 0.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFI",iplan       ,"TRD1",xpos,ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.   - kCathick;
+    par_cha[1] = fClengthI[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = 0.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAII",iplan       ,"TRD1",xpos,ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthI[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = 0.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFI",iplan       ,"TRD1",xpos,ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.   - kCcthick;
+    par_cha[1] = fClengthI[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = 0.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCII",iplan       ,"TRD1",xpos,ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // The middle chambers --------------------------------------------------------------
+
+    // the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthM1[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2.  + fClengthM1[iplan-1]/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFM",iplan         ,"TRD1",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    gMC->Gsposp("UAFM",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"MANY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthM2[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthM2[iplan-1]/2. - kSlenTR2/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFM",iplan+2*kNplan,"TRD2",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCathick;
+    par_cha[1] = fClengthM1[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2.  + fClengthM1[iplan-1]/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAIM",iplan         ,"TRD1",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    gMC->Gsposp("UAIM",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCathick;
+    par_cha[1] = fClengthM2[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthM2[iplan-1]/2. - kSlenTR2/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAIM",iplan+2*kNplan,"TRD2",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthM1[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM1[iplan-1]/2.;
+    zpos       = kCcframe/2.           - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFM",iplan         ,"TRD1",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    gMC->Gsposp("UCFM",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"MANY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthM2[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthM2[iplan-1]/2. - kSlenTR2/2.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFM",iplan+2*kNplan,"TRD2",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCcthick;
+    par_cha[1] = fClengthM1[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM1[iplan-1]/2.;
+    zpos       = kCcframe/2.           - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCIM",iplan         ,"TRD1",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    gMC->Gsposp("UCIM",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCcthick;
+    par_cha[1] = fClengthM2[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthM2[iplan-1]/2. - kSlenTR2/2.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCIM",iplan+2*kNplan,"TRD2",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // The outer chambers ---------------------------------------------------------------
+
+    // the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthO1[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM1[iplan-1]    + fClengthO1[iplan-1]/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFO",iplan         ,"TRD1",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    gMC->Gsposp("UAFO",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"MANY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthO2[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthM2[iplan-1]    + fClengthO2[iplan-1]/2. - kSlenTR2/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFO",iplan+2*kNplan,"TRD2",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthO3[iplan-1]/2.;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthO3[iplan-1]/2. - kSlenTR3/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAFO",iplan+4*kNplan,"TRD3",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the aluminum frame
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCathick;
+    par_cha[1] = fClengthO1[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM1[iplan-1]    + fClengthO1[iplan-1]/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAIO",iplan         ,"TRD1",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    gMC->Gsposp("UAIO",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCathick;
+    par_cha[1] = fClengthO2[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthM2[iplan-1]    + fClengthO2[iplan-1]/2. - kSlenTR2/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAIO",iplan+2*kNplan,"TRD2",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCathick;
+    par_cha[1] = fClengthO3[iplan-1]/2. - kCathick;
+    par_cha[2] = kCaframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthO3[iplan-1]/2. - kSlenTR3/2.;
+    zpos       = kCheight - kCaframe/2. - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UAIO",iplan+4*kNplan,"TRD3",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+    // the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthO1[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM1[iplan-1]    + fClengthO1[iplan-1]/2.;
+    zpos       = kCcframe/2.           - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFO",iplan,         "TRD1",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    gMC->Gsposp("UCFO",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"MANY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthO2[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthM2[iplan-1]    + fClengthO2[iplan-1]/2. - kSlenTR2/2.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFO",iplan+2*kNplan,"TRD2",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.;
+    par_cha[1] = fClengthO3[iplan-1]/2.;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthO3[iplan-1]/2. - kSlenTR3/2.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCFO",iplan+4*kNplan,"TRD3",xpos, ypos,zpos,0,"MANY",par_cha,npar_cha);
+
+    // the inner part of the carbon frame
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCcthick;
+    par_cha[1] = fClengthO1[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthI[iplan-1]/2. + fClengthM1[iplan-1]    + fClengthO1[iplan-1]/2.;
+    zpos       = kCcframe/2.           - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCIO",iplan         ,"TRD1",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    gMC->Gsposp("UCIO",iplan+  kNplan,"TRD1",xpos,-ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCcthick;
+    par_cha[1] = fClengthO2[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthM2[iplan-1]    + fClengthO2[iplan-1]/2. - kSlenTR2/2.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCIO",iplan+2*kNplan,"TRD2",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+    par_cha[0] = fCwidth[iplan-1]/2.    - kCcthick;
+    par_cha[1] = fClengthO3[iplan-1]/2. - kCcthick;
+    par_cha[2] = kCcframe/2.;
+    xpos       = 0.;
+    ypos       = fClengthO3[iplan-1]/2. - kSlenTR3/2.;
+    zpos       = kCcframe/2.            - kSheight/2. + (iplan-1) * (kCheight + kCspace);
+    gMC->Gsposp("UCIO",iplan+4*kNplan,"TRD3",xpos, ypos,zpos,0,"ONLY",par_cha,npar_cha);
+
+  }
+
+  xpos     = 0.;
+  ypos     = 0.;
+  zpos     = 0.;
+  gMC->Gspos("TRD1",1,"BTR1",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("TRD2",1,"BTR2",xpos,ypos,zpos,0,"ONLY");
+  gMC->Gspos("TRD3",1,"BTR3",xpos,ypos,zpos,0,"ONLY");
+
+}
diff --git a/TRD/AliTRDgeometryHole.h b/TRD/AliTRDgeometryHole.h
new file mode 100644 (file)
index 0000000..b263469
--- /dev/null
@@ -0,0 +1,34 @@
+#ifndef TRDgeometryHole_h
+#define TRDgeometryHole_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliTRDgeometry.h"
+
+class AliTRDgeometryHole : public AliTRDgeometry {
+
+ public:
+
+  AliTRDgeometryHole();
+  ~AliTRDgeometryHole();
+
+  virtual void    CreateGeometry(Int_t *); 
+  virtual Int_t   IsVersion() const { return 0; };
+  virtual void    Init();
+
+ protected:
+
+  Float_t         fClengthI[kNplan];               // Length of the inner chambers
+  Float_t         fClengthM1[kNplan];              // Length of the middle chambers
+  Float_t         fClengthM2[kNplan];              // 
+  Float_t         fClengthO1[kNplan];              // Length of the outer chambers
+  Float_t         fClengthO2[kNplan];              // 
+  Float_t         fClengthO3[kNplan];              // 
+
+  ClassDef(AliTRDgeometryHole,1)                   // TRD geometry with hole
+
+};
+
+#endif
diff --git a/TRD/AliTRDrecPoint.cxx b/TRD/AliTRDrecPoint.cxx
new file mode 100644 (file)
index 0000000..58173cf
--- /dev/null
@@ -0,0 +1,105 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  TRD reconstructed point                                                  //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDgeometry.h"
+#include "AliTRDrecPoint.h"
+#include "AliTRD.h"
+
+ClassImp(AliTRDrecPoint)
+
+//_____________________________________________________________________________
+AliTRDrecPoint::AliTRDrecPoint():AliRecPoint()
+{
+  //
+  // Standard constructor
+  //
+
+  fDetector = 0;
+
+  AliTRD *TRD;
+  if ((gAlice) &&
+      (TRD = ((AliTRD*) gAlice->GetDetector("TRD")))) {
+    fGeom = TRD->GetGeometry();
+  }
+  else {
+    fGeom = NULL;
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDrecPoint::AddDigit(Int_t digit)
+{
+  //
+  // Adds the index of a digit to the digits list
+  //
+
+  // First resize the list 
+  // (no clusters with more than 3 digits for the TRD
+  if ((fMulDigit == 0) && (fMaxDigit >= 5)) {
+    fMaxDigit = 5;
+    delete fDigitsList;
+    fDigitsList = new int[fMaxDigit];
+  }
+
+  // Increase the size of the list if necessary
+  if (fMulDigit >= fMaxDigit) { 
+    int *tempo = new (int[fMaxDigit*=2]); 
+    Int_t index; 
+    for (index = 0; index < fMulDigit; index++)
+      tempo[index] = fDigitsList[index]; 
+    delete fDigitsList; 
+    fDigitsList = tempo; 
+  }
+  
+  fDigitsList[fMulDigit++] = digit;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDrecPoint::SetLocalPosition(TVector3 &pos)
+{
+  //
+  // Sets the position of the point in the local coordinate system
+  // (row,col,time) and calculates the error matrix in the same
+  // system.
+  //
+
+  const Float_t sq12 = 3.464101615;
+
+  // Set the position
+  fLocPos = pos;
+
+  // Set the error matrix
+  // row:  pad-size / sqrt(12)
+  // col:  not defined yet
+  // time: bin-size / sqrt(12)
+  fLocPosM->operator()(0,0) = ((AliTRDgeometry *) fGeom)->GetRowPadSize()  
+                            / sq12;
+  fLocPosM->operator()(1,1) = 0.0;
+  fLocPosM->operator()(2,2) = ((AliTRDgeometry *) fGeom)->GetTimeBinSize() 
+                            / sq12;
+
+}
diff --git a/TRD/AliTRDrecPoint.h b/TRD/AliTRDrecPoint.h
new file mode 100644 (file)
index 0000000..e407647
--- /dev/null
@@ -0,0 +1,33 @@
+#ifndef TRDrecPoint_h
+#define TRDrecPoint_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliRecPoint.h"
+
+class AliTRDrecPoint : public AliRecPoint {
+
+ public:
+
+  AliTRDrecPoint();
+  virtual ~AliTRDrecPoint() {};
+  virtual void  Print(Option_t * opt = "void") {};
+  virtual void  AddDigit(Int_t digit);
+
+  virtual void  SetAmplitude(Float_t amp)       { fAmp      = amp; };
+  virtual void  SetDetector(Int_t det)          { fDetector = det; };
+  virtual void  SetLocalPosition(TVector3 &pos);
+
+  virtual Int_t GetDetector()                   { return fDetector; };
+
+ protected:
+
+  Int_t    fDetector;        // TRD detector number
+
+  ClassDef(AliTRDrecPoint,1) // Reconstructed point for the TRD
+
+};
+
+#endif
diff --git a/TRD/AliTRDsegmentArray.cxx b/TRD/AliTRDsegmentArray.cxx
new file mode 100644 (file)
index 0000000..15306cc
--- /dev/null
@@ -0,0 +1,166 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRD.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDsegmentArray.h"
+
+ClassImp(AliTRDsegmentArray)
+
+//_____________________________________________________________________________
+AliTRDsegmentArray::AliTRDsegmentArray():AliTRDsegmentArrayBase()
+{
+  //
+  // Default constructor
+  //
+
+}
+
+//_____________________________________________________________________________
+AliTRDsegmentArray::AliTRDsegmentArray(Int_t n)
+                   :AliTRDsegmentArrayBase("AliTRDdataArray",n)
+{
+  //
+  // Constructor creating an array of AliTRDdataArray of size <n>
+  //
+
+  AliTRDdataArray *DataArray;  
+
+  for (Int_t i = 0; i < n; i++) {
+    DataArray = (AliTRDdataArray *) AddSegment(i);
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDsegmentArray::Delete()
+{
+  //
+  // Deletes all detector segments from the array
+  //
+
+  for (Int_t iDet = 0; iDet < fNSegment; iDet++) {
+    ClearSegment(iDet);
+  }
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDsegmentArray::LoadArray(const Char_t *branchname)
+{
+  //
+  // Loads all segments of the array from the branch <branchname> of
+  // the digits tree
+  //
+
+  // Connect the digits tree
+  fTree = gAlice->TreeD();
+  if (!fTree) return kFALSE;
+
+  // Get the branch
+  fBranch = fTree->GetBranch(branchname);
+  if (!fBranch) return kFALSE;
+
+  // Loop through all segments and read them from the tree
+  Bool_t status = kTRUE;
+  for (Int_t iSegment = 0; iSegment < fNSegment; iSegment++) {
+    AliTRDdataArray *DataArray = (AliTRDdataArray *) fSegment->At(iSegment);
+    if (!DataArray) {
+      status = kFALSE;
+      break;    
+    }
+    fBranch->SetAddress(&DataArray);
+    fBranch->GetEntry(iSegment);
+  }
+
+  return status;
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDsegmentArray::StoreArray(const Char_t *branchname)
+{
+  //
+  // Stores all segments of the array in the branch <branchname> of 
+  // the digits tree
+  //
+
+  // Connect the digits tree
+  fTree = gAlice->TreeD();
+  if (!fTree) return kFALSE;
+
+  // Get the branch
+  fBranch = fTree->GetBranch(branchname);
+  if (!fBranch) return kFALSE;
+
+  // Loop through all segments and fill them into the tree
+  Bool_t status = kTRUE;
+  for (Int_t iSegment = 0; iSegment < fNSegment; iSegment++) {
+    const AliTRDdataArray *DataArray = 
+         (AliTRDdataArray *) AliTRDsegmentArrayBase::At(iSegment);
+    if (!DataArray) {
+      status = kFALSE;
+      break;
+    }
+    fBranch->SetAddress(&DataArray);
+    fBranch->Fill();
+  }
+
+  return status;
+
+}
+
+//_____________________________________________________________________________
+AliTRDdataArray *AliTRDsegmentArray::GetDataArray(Int_t det)
+{
+  //
+  // Returns the data array for a given detector
+  //
+
+  return ((AliTRDdataArray *) AliTRDsegmentArrayBase::At(det));
+
+}
+
+//_____________________________________________________________________________
+AliTRDdataArray *AliTRDsegmentArray::GetDataArray(Int_t pla, Int_t cha, Int_t sec)
+{
+  //
+  // Returns the data array for a given detector
+  //
+
+  if (gAlice) {
+
+    AliTRDgeometry *Geo = ((AliTRD*) gAlice->GetDetector("TRD"))->GetGeometry();  
+    Int_t det = Geo->GetDetector(pla,cha,sec);
+    return GetDataArray(det);
+
+  }
+  else {
+
+    printf("AliTRDsegmentArray::GetDigits -- ");
+    printf("gAlice is not defined\n");
+    return NULL;
+
+  }
+
+}
diff --git a/TRD/AliTRDsegmentArray.h b/TRD/AliTRDsegmentArray.h
new file mode 100644 (file)
index 0000000..2d296ea
--- /dev/null
@@ -0,0 +1,38 @@
+#ifndef TRDsegmentArray_H
+#define TRDsegmentArray_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////
+//  Array for TRD detector segments containing digits //
+////////////////////////////////////////////////////////
+
+#include "AliTRDsegmentArrayBase.h"
+#include "AliTRDdataArray.h"
+
+//_____________________________________________________________________________
+class AliTRDsegmentArray : public AliTRDsegmentArrayBase {
+
+ public:
+
+  AliTRDsegmentArray();
+  AliTRDsegmentArray(Int_t n);
+  ~AliTRDsegmentArray() { };
+
+  virtual void             Delete();
+
+  virtual Bool_t           LoadArray(const Char_t *branchname);
+  virtual Bool_t           StoreArray(const Char_t *branchname);
+
+  virtual AliTRDdataArray *GetDataArray(Int_t det);
+  virtual AliTRDdataArray *GetDataArray(Int_t sec, Int_t cha, Int_t pla);
+
+ protected:
+
+  ClassDef(AliTRDsegmentArray,1)        // TRD detector segment array 
+
+};
+
+#endif
diff --git a/TRD/AliTRDsegmentArrayBase.cxx b/TRD/AliTRDsegmentArrayBase.cxx
new file mode 100644 (file)
index 0000000..37b1f04
--- /dev/null
@@ -0,0 +1,325 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Alice segment manager object                                             //
+//                                                                           //
+//  AliTRDsegmentIDArray object  is array of pointers to object derived from //
+//  AliTRDsegmentID object                                                   //
+//  AliTRDsegmentID - object in comparison with TObject enhalt               //
+//  additional information fSegmentID                                        //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include  <TROOT.h>
+#include <TTree.h>
+#include "TClonesArray.h"
+#include "TDirectory.h"
+#include "AliTRDarrayI.h"
+#include "TError.h"
+#include "TClass.h"
+
+#include "AliTRDsegmentID.h"
+#include "AliTRDsegmentArrayBase.h"
+
+//_____________________________________________________________________________
+ClassImp(AliTRDsegmentArrayBase)
+  
+AliTRDsegmentArrayBase::AliTRDsegmentArrayBase()
+{
+  //
+  //
+  //
+  fNSegment=0;
+  fSegment =0; 
+  fTreeIndex = 0;
+  fTree  = 0;
+  fClass = 0;
+}
+
+AliTRDsegmentArrayBase::AliTRDsegmentArrayBase(Text_t *classname, Int_t n)
+{
+  //
+  //constructor which 
+  // 
+  //  Create an array of objects of classname. The class must inherit from
+  //  AliTRDsegmentID .  The second argument adjust number of entries in 
+  //  the array.
+  fNSegment=0;
+  fSegment =0; 
+  fTreeIndex = 0;
+  fTree  = 0;
+  fClass = 0;
+  SetClass(classname);
+  if (MakeArray(n)==kFALSE){
+     Error("AliTRDsegmentArrayBase", "can't allocate %d segments in memory",n);
+     return;
+   }
+}
+
+Bool_t AliTRDsegmentArrayBase:: SetClass(Text_t *classname)
+{
+  //
+  //set class of stored object
+  if ( fClass !=0 ) {
+    delete fClass;
+    fClass = 0;
+  }
+  if (fTree !=0) {
+    delete fTree;
+    fTree = 0;
+    fBranch = 0;
+    delete fTreeIndex;
+    fTreeIndex = 0;
+  } 
+  if (fSegment != 0) {
+    fSegment->Delete();
+    delete fSegment;
+    fSegment = 0;
+  }
+  if (!gROOT)
+      ::Fatal("AliTRDsegmentArrayBase::AliTRDsegmentArrayBase", "ROOT system not initialized");
+   
+   fClass = gROOT->GetClass(classname);
+   if (!fClass) {
+      Error("AliTRDsegmentArrayBase", "%s is not a valid class name", classname);
+      return kFALSE;
+   }
+   if (!fClass->InheritsFrom(AliTRDsegmentID::Class())) {
+      Error("AliTRDsegmentArrayBase", "%s does not inherit from AliTRDsegmentID", classname);
+      return kFALSE;
+   }  
+   return kTRUE;
+}
+
+//Bool_t AliTRDsegmentArrayBase::ClassError( )
+//{
+  //signalize class error 
+  //  if (!fClass) {
+  //    Error("AliTRDsegmentArrayBase", "%s is not a valid class name", classname);
+  //    return kFALSE;
+  // }
+////  return kFALSE;
+//}
+
+AliTRDsegmentArrayBase::~AliTRDsegmentArrayBase()
+{
+  if (fNSegment>0){
+    fSegment->Delete();
+    delete fSegment;
+  }
+  if (fTree) delete fTree;
+  if (fTreeIndex) delete fTreeIndex;
+  if (fClass!=0) delete fClass;
+}
+
+AliTRDsegmentID * AliTRDsegmentArrayBase::NewSegment()
+{
+  //
+  //create object according class information
+  if (fClass==0) return 0;
+  AliTRDsegmentID * segment = (AliTRDsegmentID * )fClass->New();
+  if (segment == 0) return 0;
+  return segment;
+}
+
+
+Bool_t AliTRDsegmentArrayBase::AddSegment(AliTRDsegmentID *segment)
+{
+  //
+  // add segment to array
+  //
+  if (segment==0) return kFALSE;
+  if (fSegment==0) return kFALSE;
+  if (fClass==0) return kFALSE;
+  if (!(segment->IsA()->InheritsFrom(fClass))){
+    Error("AliTRDsegmentArrayBase", "added class %s  is not of proper type ",
+         segment->IsA()->GetName());
+      return kFALSE;
+  }
+  fSegment->AddAt(segment,segment->GetID());
+  fNSegment = fSegment->GetLast()+1;
+  return kTRUE;
+}
+
+AliTRDsegmentID * AliTRDsegmentArrayBase::AddSegment(Int_t index)
+{
+  //
+  // add segment to array
+  //
+  if (fSegment==0) return 0;
+  if (fClass==0) return 0;
+  //  AliTRDsegmentID * segment = (AliTRDsegmentID * )fClass->New();
+  AliTRDsegmentID * segment = NewSegment();
+  if (segment == 0) return 0;
+  fSegment->AddAt(segment,index);
+  segment->SetID(index);
+  fNSegment = fSegment->GetLast()+1;
+  return segment;
+}
+
+
+
+Bool_t AliTRDsegmentArrayBase::MakeArray(Int_t n)
+{
+  //
+  //make array of pointers to Segments
+  //
+  if (fSegment) {
+    fSegment->Delete();
+    delete fSegment;
+  }
+  if (fTreeIndex) delete   fTreeIndex;  
+  fSegment = new TObjArray(n);
+  fTreeIndex = new AliTRDarrayI;
+  fTreeIndex->Set(n);
+  fNSegment=n;
+  if ( (fSegment) && (fTreeIndex)) return kTRUE;
+  else return kFALSE;            
+}
+
+
+void AliTRDsegmentArrayBase::ClearSegment(Int_t index)
+{
+  //
+  //remove segment from active memory    
+  //
+  if ((*fSegment)[index]){
+    //    (*fSegment)[index]->Delete(); //not working for TClonesArray
+    delete (*fSegment)[index]; //because problem with deleting TClonesArray
+    fSegment->RemoveAt(index);
+  }
+}
+
+
+void AliTRDsegmentArrayBase::MakeTree()
+{
+  //  AliTRDsegmentID  segment;
+  AliTRDsegmentID * psegment = NewSegment();  
+  if (fTree) delete fTree;
+  fTree = new TTree("Segment Tree","Tree with segments");
+  fBranch = fTree->Branch("Segment",psegment->IsA()->GetName(),&psegment,64000,1);
+  delete psegment;
+}              
+
+Bool_t AliTRDsegmentArrayBase::ConnectTree(const char * treeName)
+{
+  //connect tree from current directory  
+  if (fTree){
+    delete fTree;
+    fTree = 0;
+    fBranch = 0;
+  }
+  fTree =(TTree*)gDirectory->Get(treeName);
+  if (fTree == 0)    return kFALSE;
+  fBranch = fTree->GetBranch("Segment");
+  if (fBranch==0) return kFALSE;
+  MakeDictionary(TMath::Max(fNSegment,Int_t(fTree->GetEntries())));
+  return kTRUE;
+}
+
+AliTRDsegmentID *AliTRDsegmentArrayBase::LoadSegment(Int_t index)
+{
+  //
+  //load segment with index to the memory
+  //
+  //
+  if (fTreeIndex ==0 ) MakeDictionary(3000);
+  //firstly try to load dictionary 
+  if (fTreeIndex ==0 ) return 0;
+  if (fBranch==0) return 0;
+  if (index>fTreeIndex->fN) return 0;
+  AliTRDsegmentID *s = (AliTRDsegmentID*)(*fSegment)[index];
+  if (s==0)  s=  NewSegment();
+  s->SetID(index);
+  //  new AliTRDsegmentID(index);
+  
+  if (s!=0) {
+    Int_t treeIndex =(*fTreeIndex)[index];
+    if (treeIndex<1) return 0;
+    else treeIndex--;   //I don't like it Int table I have index shifted by 1                 
+    fBranch->SetAddress(&s);
+    fTree->GetEvent(treeIndex);
+    (*fSegment)[index] = (TObject*) s;
+  }
+  else 
+    return 0;
+  return s;
+  //  AbstractMethod("LoadSegment");
+}
+AliTRDsegmentID *AliTRDsegmentArrayBase::LoadEntry(Int_t index)
+{
+  //
+  //load segment at position inex in tree  to the memory
+  //
+  //
+  if (fBranch==0) return 0;
+  if (index>fTree->GetEntries()) return 0;
+  AliTRDsegmentID * s =  NewSegment();
+  
+  if (s) {
+    fBranch->SetAddress(&s);
+    fTree->GetEvent(index);
+  }
+  else 
+    return 0;
+  Int_t nindex = s->GetID();
+  ClearSegment(nindex);
+  (*fSegment)[nindex] = (TObject*) s;
+  return s;
+  //  AbstractMethod("LoadSegment");
+}
+
+void AliTRDsegmentArrayBase::StoreSegment(Int_t index)
+{
+  //
+  //make segment persistent 
+  //
+  const AliTRDsegmentID *  segment = (*this)[index];
+  if (segment == 0 ) return;
+  if (fTree==0) MakeTree();
+  fBranch->SetAddress(&segment);
+  fTree->Fill();
+}
+
+Bool_t  AliTRDsegmentArrayBase::MakeDictionary(Int_t size)
+{
+  //
+  //create index table for tree
+  //  
+  if (size<1) return kFALSE;
+  if (fTreeIndex) delete fTreeIndex;
+  fTreeIndex = new AliTRDarrayI(); 
+  fTreeIndex->Set(size);
+  
+  AliTRDsegmentID  segment;
+  AliTRDsegmentID * psegment = &segment;
+  fBranch->SetAddress(&psegment);
+  TBranch * brindix = fTree->GetBranch("fSegmentID");
+  Int_t nevent = (Int_t)fTree->GetEntries();  
+  for (Int_t i = 0; i<nevent; i++){
+    brindix->GetEvent(i);
+    Int_t treeIndex=segment.GetID();
+    if (fTreeIndex->fN<treeIndex) fTreeIndex->Expand(Int_t(Float_t(treeIndex)*1.5)+1);
+    //    Int_t index = segment.GetID(); 
+    (*fTreeIndex)[treeIndex]=i+1; // MI 19.5. I'm sorry  -index 0 couldn't be use in AliTRDarrayI   
+  }
+  return kTRUE;
+}
diff --git a/TRD/AliTRDsegmentArrayBase.h b/TRD/AliTRDsegmentArrayBase.h
new file mode 100644 (file)
index 0000000..cd6ceae
--- /dev/null
@@ -0,0 +1,79 @@
+#ifndef ALISEGARRAYBASE_H
+#define ALISEGARRAYBASE_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDsegmentArrayBase.h,v */
+
+////////////////////////////////////////////////
+//  Manager class generaol Alice segment 
+//  segment is for example one pad row in TPC //
+////////////////////////////////////////////////
+
+#include "TNamed.h"
+#include "TError.h"
+#include "TObjArray.h"
+
+//#include "AliTRDsegmentID.h"
+
+class TTree;
+class TBranch;
+class AliTRDarrayI;
+class AliTRDsegmentID;
+class TObjArray;
+class AliTRDsegmentArrayBase: public TNamed{
+public:
+  AliTRDsegmentArrayBase();
+  AliTRDsegmentArrayBase(Text_t *classname, Int_t n);  //
+  Bool_t  SetClass(Text_t *classname);  //set class of stored object
+  ~AliTRDsegmentArrayBase();
+  inline const AliTRDsegmentID * At(Int_t i); //return pointer to segment with index i 
+  inline const AliTRDsegmentID * operator[](Int_t i); //return pointer to segment with index i
+
+  Bool_t AddSegment(AliTRDsegmentID *segment); // add segment to array
+  AliTRDsegmentID * AddSegment(Int_t index);   //create objet and set index
+  Bool_t   MakeArray(Int_t n);       //make array of pointers to Segments
+  void ClearSegment(Int_t index); //remove segment from active   
+  virtual AliTRDsegmentID * NewSegment(); //dynamicaly create new segment 
+  //input output functions
+  TTree * GetTree(){return fTree;}      //return pointer to connected tree
+  
+  virtual void MakeTree();              //Make tree with the name
+  virtual Bool_t ConnectTree(const char * treeName); //connect tree from current directory 
+  virtual AliTRDsegmentID * LoadSegment(Int_t index);//load segment with index to the memory
+  virtual AliTRDsegmentID * LoadEntry(Int_t index); //load segment entry from position index in tree
+  virtual void StoreSegment(Int_t index);//write segmen persistent  
+  Bool_t  MakeDictionary(Int_t size);//create index table for tree
+  TClass * GetClass() {return fClass;}
+public:
+  TObjArray  * fSegment;  //!pointer to array of pointers to segment
+  AliTRDarrayI    * fTreeIndex; //!pointers(index) table in tree
+  Int_t      fNSegment;   
+  TTree    * fTree;   //!tree with segment objects
+  TBranch  * fBranch; //!total branch 
+private: 
+  TClass  *   fClass;    //!class type of included objects 
+  ClassDef(AliTRDsegmentArrayBase,1) 
+};
+
+
+
+const AliTRDsegmentID*  AliTRDsegmentArrayBase::operator[](Int_t i)
+{
+  //
+  //return segment with given index
+  //
+  if ( (i<0) || (i>=fNSegment)) return 0; 
+  return (AliTRDsegmentID *)fSegment->At(i);
+}
+const AliTRDsegmentID*  AliTRDsegmentArrayBase::At(Int_t i)
+{
+  //
+  //return segment with given index
+  //
+  if ( (i<0) || (i>=fNSegment)) return 0; 
+  return (AliTRDsegmentID *)((*fSegment)[i]);
+}
+
+#endif //ALISEGARRAY_H
diff --git a/TRD/AliTRDsegmentID.cxx b/TRD/AliTRDsegmentID.cxx
new file mode 100644 (file)
index 0000000..d13a982
--- /dev/null
@@ -0,0 +1,34 @@
+/**************************************************************************
+ * 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$
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Alice  AliSementID   object                                             //
+//                                
+//                                                                           //
+//                                                                          //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDsegmentID.h"
+
+ClassImp(AliTRDsegmentID)
+
+AliTRDsegmentID::AliTRDsegmentID()
+{
+}
diff --git a/TRD/AliTRDsegmentID.h b/TRD/AliTRDsegmentID.h
new file mode 100644 (file)
index 0000000..6a98df3
--- /dev/null
@@ -0,0 +1,29 @@
+#ifndef ALITRDSEGMENTID_H
+#define ALITRDSEGMENTID_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDsegmentID.h,v */
+
+////////////////////////////////////////////////
+//  Manager class generaol Alice segment 
+//  segment is for example one pad row in TPC //
+////////////////////////////////////////////////
+
+#include "TObject.h"
+
+class AliTRDsegmentID: public TObject{
+public:
+  AliTRDsegmentID();
+  AliTRDsegmentID(Int_t index){fSegmentID = index;}
+  Int_t GetID() {return fSegmentID;}
+  void  SetID(Int_t index){fSegmentID = index;} 
+  virtual Int_t  GetSize() {return sizeof(*this);} //function which return size of object 
+protected:
+  Int_t fSegmentID;   //identification number of Segment
+  ClassDef(AliTRDsegmentID,1) 
+};
+   
+#endif //ALISEGMENTID_H
+
diff --git a/TRD/slowClusterAna.C b/TRD/slowClusterAna.C
new file mode 100644 (file)
index 0000000..bb646c2
--- /dev/null
@@ -0,0 +1,94 @@
+void slowClusterAna() {
+
+/////////////////////////////////////////////////////////////////////////
+//
+// Example macro for the analysis of the TRD cluster
+//
+/////////////////////////////////////////////////////////////////////////
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+  }
+
+  // Input file name
+  Char_t *alifile = "galice_c_v1.root"; 
+
+  // Event number
+  Int_t   nEvent  = 0;
+
+  // Define the histograms
+  TH1F *hEnergy = new TH1F("hEnergy","Cluster energy",100,0.0,1000.0);
+
+  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
+  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
+  if (!gafl) {
+    cout << "Open the ALIROOT-file " << alifile << endl;
+    gafl = new TFile(alifile);
+  }
+  else {
+    cout << alifile << " is already open" << endl;
+  }
+
+  // Get AliRun object from file or create it if not on file
+  gAlice = (AliRun*) gafl->Get("gAlice");
+  if (gAlice)  
+    cout << "AliRun object found on file" << endl;
+  else
+    gAlice = new AliRun("gAlice","Alice test program");
+
+  // Import the Trees for the event nEvent in the file
+  Int_t nparticles = gAlice->GetEvent(nEvent);
+  if (nparticles <= 0) break;
+  
+  // Get the pointer to the hit-tree
+  TTree          *RecTree       = gAlice->TreeR();
+  RecTree->Print();
+  // Get the pointer to the detector classes
+  AliTRDv1       *TRD           = (AliTRDv1*) gAlice->GetDetector("TRD");
+  // Get the geometry
+  AliTRDgeometry *TRDgeometry   = TRD->GetGeometry();
+  // Get the pointer to the hit container
+  TClonesArray   *RecPointArray = TRD->RecPoints();
+  // Set the branch address
+  RecTree->GetBranch("TRDrecPoints")->SetAddress(&RecPointArray);
+
+  Int_t nEntries = RecTree->GetEntries();
+
+  // Loop through all entries in the tree
+  Int_t nbytes;
+  for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
+
+    cout << "iEntry = " << iEntry << endl;
+
+    // Import the tree
+    nbytes += RecTree->GetEvent(iEntry);
+
+    // Get the number of digits in the detector 
+    Int_t nRecPoint = RecPointArray->GetEntriesFast();
+    cout << " nRecPoint = " << nRecPoint << endl;    
+
+    // Loop through all TRD digits
+    for (Int_t iRecPoint = 0; iRecPoint < nRecPoint; iRecPoint++) {
+
+      // Get the information for this digit
+      AliTRDrecPoint *RecPoint = (AliTRDrecPoint *) RecPointArray->UncheckedAt(iRecPoint);
+      Int_t    detector = RecPoint->GetDetector();      
+      Int_t    sector   = TRDgeometry->GetSector(detector);
+      Int_t    plane    = TRDgeometry->GetPlane(detector);
+      Int_t    chamber  = TRDgeometry->GetChamber(detector);
+      Int_t    energy   = RecPoint->GetEnergy();
+      TVector3 pos;
+      RecPoint->GetLocalPosition(pos);
+
+      hEnergy->Fill((Float_t) energy);
+
+    }
+
+  }
+
+  TCanvas *c1 = new TCanvas("c1","Cluster",50,50,600,400);
+  hEnergy->Draw();
+
+}