Update TRD code from C.Blume
authorfca <fca@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Sep 1999 16:59:12 +0000 (16:59 +0000)
committerfca <fca@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Sep 1999 16:59:12 +0000 (16:59 +0000)
15 files changed:
TRD/AliTRD.cxx
TRD/AliTRD.h
TRD/AliTRDconst.h
TRD/AliTRDmatrix.cxx [new file with mode: 0644]
TRD/AliTRDmatrix.h [new file with mode: 0644]
TRD/AliTRDpixel.cxx [new file with mode: 0644]
TRD/AliTRDpixel.h [new file with mode: 0644]
TRD/AliTRDsim.cxx [new file with mode: 0644]
TRD/AliTRDsim.h [new file with mode: 0644]
TRD/AliTRDv2.cxx
TRD/AliTRDv2.h
TRD/Makefile
TRD/TRDLinkDef.h
TRD/digitsAna.C [new file with mode: 0644]
TRD/digitsCreate.C [new file with mode: 0644]

index d895d25..20ad31d 100644 (file)
@@ -36,6 +36,8 @@ AliTRD::AliTRD()
 
   fIshunt      = 0;
   fGasMix      = 0;
+  fHits        = 0;
+  fDigits      = 0;
 
   // The chamber dimensions
   for (Int_t iplan = 0; iplan < kNplan; iplan++) {
@@ -64,7 +66,10 @@ AliTRD::AliTRD(const char *name, const char *title)
 
   // Allocate the hit array
   fHits   = new TClonesArray("AliTRDhit",  405);
-  
+
+  // Allocate the digits array
+  fDigits = new TClonesArray("AliTRDdigit",10000);
+   
   fIshunt = 0;
   fGasMix = 0;
 
@@ -79,7 +84,33 @@ AliTRD::AliTRD(const char *name, const char *title)
   SetMarkerColor(kWhite);   
 
 }
+
+//_____________________________________________________________________________
+AliTRD::~AliTRD()
+{
+  //
+  // TRD destructor
+  //
+
+  fIshunt = 0;
+
+  delete fHits;
+  delete fDigits;
+
+}
+
+//_____________________________________________________________________________
+void AliTRD::AddDigit(Int_t *tracks, Int_t *digits)
+{
+  //
+  // Add a digit for the TRD
+  //
+
+  TClonesArray &ldigits = *fDigits;
+  new(ldigits[fNdigits++]) AliTRDdigit(tracks,digits);
+
+}
+
 //_____________________________________________________________________________
 void AliTRD::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 {
@@ -743,3 +774,28 @@ AliTRDhit::AliTRDhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
   fQ       = hits[3];
 
 }
+
+ClassImp(AliTRDdigit)
+
+//_____________________________________________________________________________
+AliTRDdigit::AliTRDdigit(Int_t *tracks, Int_t *digits)
+            :AliDigit(tracks)
+{
+  //
+  // Create a TRD digit
+  //
+
+  // Store the volume hierarchy
+  fSector  = digits[0];
+  fChamber = digits[1];
+  fPlane   = digits[2];
+
+  // Store the row, pad, and time bucket number
+  fRow     = digits[3];
+  fCol     = digits[4];
+  fTime    = digits[5];
+
+  // Store the signal amplitude
+  fSignal  = digits[6];
+
+}
index d299823..0bdc0d6 100644 (file)
@@ -6,6 +6,7 @@
  
 #include "AliDetector.h"
 #include "AliHit.h" 
+#include "AliDigit.h"
 #include "AliTRDconst.h"
 
 //_____________________________________________________________________________
@@ -22,8 +23,9 @@ protected:
 public:
   AliTRD();
   AliTRD(const char *name, const char *title);
-  virtual      ~AliTRD() {}
+  virtual      ~AliTRD();
   virtual void  AddHit(Int_t, Int_t*, Float_t*);
+  virtual void  AddDigit(Int_t*, Int_t*);    
   virtual void  BuildGeometry();
   virtual void  CreateGeometry();
   virtual void  CreateMaterials();
@@ -54,10 +56,32 @@ public:
 public:
   AliTRDhit() {}
   AliTRDhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits);
-  virtual ~AliTRDhit() {}
+  virtual ~AliTRDhit() {};
  
   ClassDef(AliTRDhit,1)     // Hits for Transition Radiation Detector
 
 };
 
+//_____________________________________________________________________________
+class AliTRDdigit : public AliDigit {
+
+public:
+  Int_t        fSector;     // TRD sector number
+  Int_t        fChamber;    // TRD chamber number
+  Int_t        fPlane;      // TRD plane number
+  Int_t        fRow;        // Pad row number
+  Int_t        fCol;        // Pad col number
+  Int_t        fTime;       // Time bucket
+  Int_t        fSignal;     // Signal amplitude
+
+public:
+  AliTRDdigit() {};
+  AliTRDdigit(Int_t *tracks, Int_t *digits);
+  virtual ~AliTRDdigit() {};
+
+  ClassDef(AliTRDdigit,1)   // Digits for Transition Radiation Detector
+
+};
+
+
 #endif
index f7c9c42..5cda810 100644 (file)
@@ -1,12 +1,10 @@
 //
-// Parameter for the TRD
+// Geometry parameter for the TRD
 //
 
-//_____________________________________________________________________________
-//  Geometry parameter
-
 const Int_t   kNsect   = 18;      // Number of sectors in the full detector
 const Int_t   kNplan   = 6;       // Number of planes of the TRD
+const Int_t   kNcham   = 5;       // Number of chambers in z-direction
 
 const Float_t kRmin    = 294.0;   // r-dimensions of the TRD
 const Float_t kRmax    = 368.0;
@@ -28,9 +26,6 @@ const Float_t kCcthick =   1.0;   // Thickness of the carbon frame
 const Float_t kCwidcha = (kSwidth2 - kSwidth1) / kSheight * (kCheight + kCspace);
 const Float_t kCcframe = kCheight - kCaframe;
 
-// Number of Radiator foils
-const Int_t   kRaFoils = 100;
-
 // Thicknesses of the the material layers
 const Float_t kSeThick = 0.02;                // Radiator seal
 const Float_t kRaThick = 4.8;                 // Radiator
diff --git a/TRD/AliTRDmatrix.cxx b/TRD/AliTRDmatrix.cxx
new file mode 100644 (file)
index 0000000..b429bad
--- /dev/null
@@ -0,0 +1,341 @@
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Contains the pixel information for one TRD chamber                       //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDmatrix.h"
+
+ClassImp(AliTRDmatrix)
+
+//_____________________________________________________________________________
+AliTRDmatrix::AliTRDmatrix():TObject()
+{
+  //
+  // Create a TRD detector matrix
+  // 
+  fRow        = 0;
+  fCol        = 0;
+  fTime       = 0;
+  fPixel      = 0;
+  fSector     = 0;
+  fChamber    = 0;
+  fPlane      = 0;
+  fPixelArray = NULL;
+
+}
+
+//_____________________________________________________________________________
+AliTRDmatrix::AliTRDmatrix(Int_t nRow, Int_t nCol, Int_t nTime
+                         , Int_t iSec, Int_t iCha, Int_t iPla)
+             :TObject()
+{
+  //
+  // Create a TRD detector matrix with a given size
+  // 
+
+  fRow        = nRow;
+  fCol        = nCol;
+  fTime       = nTime;
+  fPixel      = nRow * nCol * nTime;
+  fSector     = iSec;
+  fChamber    = iCha;
+  fPlane      = iPla;
+  fPixelArray = new TObjArray(fPixel);
+  for (Int_t iPixel = 0; iPixel < fPixel; iPixel++) {
+    AliTRDpixel *pixel = new AliTRDpixel();
+    fPixelArray->Add(pixel);
+  }
+
+}
+
+//_____________________________________________________________________________
+AliTRDmatrix::~AliTRDmatrix()
+{
+
+  if (fPixelArray) {
+    fPixelArray->Delete();
+    delete fPixelArray;
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmatrix::AddSignal(Int_t iRow, Int_t iCol, Int_t iTime, Float_t signal)
+{
+  //
+  // Add a value to the amplitude of the signal for one specific pixel
+  //
+
+  AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+  if (pixel) {
+    signal += pixel->GetSignal();
+    pixel->SetSignal(signal);
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmatrix::Draw()
+{
+  //
+  // Draws a 3D view of the detector matrix
+  //
+
+  Char_t ctitle[50];
+  sprintf(ctitle,"Matrix (Sector:%d Chamber:%d Plane:%d)"
+                ,fSector,fChamber,fPlane);
+  TH3F *hMatrix = new TH3F("hMatrix",ctitle,fRow ,-0.5,fRow +0.5
+                                           ,fCol ,-0.5,fCol +0.5
+                                           ,fTime,-0.5,fTime+0.5);
+
+  for (Int_t iRow  = 0; iRow  < fRow;  iRow++ ) {
+    for (Int_t iCol  = 0; iCol  < fCol;  iCol++ ) {
+      for (Int_t iTime = 0; iTime < fTime; iTime++) {
+        AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+        if (pixel) hMatrix->Fill3(iRow,iCol,iTime,pixel->GetSignal());
+      }
+    }
+  }
+
+  gStyle->SetOptStat(0);
+  TCanvas *cMatrix = new TCanvas("cMatrix","Detector matrix 3D-view"
+                                          ,50,50,600,400);
+  cMatrix->ToggleEventStatus();
+  hMatrix->SetXTitle("Pad-row (z)");
+  hMatrix->SetYTitle("Pad-column (rphi)");
+  hMatrix->SetZTitle("Timebucket");
+  hMatrix->Draw("BOX");
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmatrix::DrawRow(Int_t iRow)
+{
+  //
+  // Draws a 2D slice of the detector matrix along one row
+  //
+
+  if ((iRow < 0) || (iRow >= fRow)) {
+    printf("Index out of bounds (%d/%d)\n",iRow,fRow);
+    return;
+  }
+
+  Char_t ctitle[50];
+  sprintf(ctitle,"Pad-row %d (Sector:%d Chamber:%d Plane:%d)"
+                ,iRow,fSector,fChamber,fPlane);
+  TH2F *hSliceRow = new TH2F("hSliceRow",ctitle,fCol ,-0.5,fCol +0.5
+                                               ,fTime,-0.5,fTime+0.5);
+
+  for (Int_t iCol  = 0; iCol  < fCol;  iCol++ ) {
+    for (Int_t iTime = 0; iTime < fTime; iTime++) {
+      AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+      if (pixel) hSliceRow->Fill(iCol,iTime,pixel->GetSignal());
+    }
+  }
+
+  gStyle->SetOptStat(0);
+  TCanvas *cSliceRow = new TCanvas("cSliceRow","Detector matrix 2D-slice"
+                                              ,50,50,600,400);
+  cSliceRow->ToggleEventStatus();
+  hSliceRow->SetXTitle("Pad-column (rphi)");
+  hSliceRow->SetYTitle("Timebucket");
+  hSliceRow->Draw("COLZ");
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmatrix::DrawCol(Int_t iCol)
+{
+  //
+  // Draws a 2D slice of the detector matrix along one column
+  //
+
+  if ((iCol < 0) || (iCol >= fCol)) {
+    printf("Index out of bounds (%d/%d)\n",iCol,fCol);
+    return;
+  }
+
+  Char_t ctitle[50];
+  sprintf(ctitle,"Pad-column %d (Sector:%d Chamber:%d Plane:%d)"
+                ,iCol,fSector,fChamber,fPlane);
+  TH2F *hSliceCol = new TH2F("hSliceCol",ctitle,fRow ,-0.5,fRow +0.5
+                                               ,fTime,-0.5,fTime+0.5);
+
+  for (Int_t iRow  = 0; iRow  < fRow;  iRow++ ) {
+    for (Int_t iTime = 0; iTime < fTime; iTime++) {
+      AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+      if (pixel) hSliceCol->Fill(iRow,iTime,pixel->GetSignal());
+    }
+  }
+
+  gStyle->SetOptStat(0);
+  TCanvas *cSliceCol = new TCanvas("cSliceCol","Detector matrix 2D-slice"
+                                              ,50,50,600,400);
+  cSliceCol->ToggleEventStatus();
+  hSliceCol->SetXTitle("Pad-row (z)");
+  hSliceCol->SetYTitle("Timebucket");
+  hSliceCol->Draw("COLZ");
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmatrix::DrawTime(Int_t iTime)
+{
+  //
+  // Draws a 2D slice of the detector matrix along one time slice
+  //
+
+  if ((iTime < 0) || (iTime >= fTime)) {
+    printf("Index out of bounds (%d/%d)\n",iTime,fTime);
+    return;
+  }
+
+  Char_t ctitle[50];
+  sprintf(ctitle,"Time-slice %d (Sector:%d Chamber:%d Plane:%d)"
+                ,iTime,fSector,fChamber,fPlane);
+  TH2F *hSliceTime = new TH2F("hSliceTime",ctitle,fRow,-0.5,fRow+0.5
+                                                 ,fCol,-0.5,fCol+0.5);
+
+  for (Int_t iRow = 0; iRow < fRow; iRow++) {
+    for (Int_t iCol = 0; iCol < fCol; iCol++) {
+      AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+      if (pixel) hSliceTime->Fill(iRow,iCol,pixel->GetSignal());
+    }
+  }
+
+  gStyle->SetOptStat(0);
+  TCanvas *cSliceTime = new TCanvas("cSliceTime","Detector matrix 2D-slice"
+                                                ,50,50,600,400);
+  cSliceTime->ToggleEventStatus();
+  hSliceTime->SetXTitle("Pad-row (z)");
+  hSliceTime->SetYTitle("Pad-column (rphi)");
+  hSliceTime->Draw("COLZ");
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmatrix::SetSignal(Int_t iRow, Int_t iCol, Int_t iTime, Float_t signal)
+{
+  //
+  // Set the amplitude of the signal for one specific pixel
+  //
+
+  AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+  if (pixel) {
+    pixel->SetSignal(signal);
+  }
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDmatrix::AddTrack(Int_t iRow, Int_t iCol, Int_t iTime, Int_t track)
+{
+  // 
+  // Add this track number to the stored tracks passing through this pixel. 
+  // If there are already three stored the return status is FALSE.  
+  //
+
+  AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+  if (!(pixel)) return kTRUE;
+
+  Bool_t trackSet = kFALSE;
+  for (Int_t i = 0; i < kTrackPixel; i++) {
+    if (pixel->GetTrack(i) == track) {
+      trackSet = kTRUE;
+      break;
+    }
+    if (pixel->GetTrack(i) ==     0) {
+      pixel->SetTrack(i,track);
+      trackSet = kTRUE;
+      break;
+    }
+  }    
+
+  return trackSet;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmatrix::SetTrack(Int_t iRow, Int_t iCol, Int_t iTime
+                          , Int_t iTrack, Int_t track)
+{
+  //
+  // Store the number of a track which is passing through this pixel
+  //
+
+  AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+  if (pixel) {
+    pixel->SetTrack(iTrack,track);
+  }
+
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDmatrix::GetSignal(Int_t iRow, Int_t iCol, Int_t iTime)
+{
+  //
+  // Returns the amplitude of the signal for one specific pixel
+  //
+
+  AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+  if (pixel) {
+    return (pixel->GetSignal());
+  }
+  else {
+    return 0;
+  }
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDmatrix::GetTrack(Int_t iRow, Int_t iCol, Int_t iTime, Int_t iTrack)
+{
+  //
+  // Returns the numbers of the tracks passing through one specific pixel
+  //
+
+  if ((iTrack < 0) || (iTrack >= kTrackPixel)) {
+    printf("GetTrack: Index out of bounds (%d)\n",iTrack);
+    return 0;
+  }
+
+  AliTRDpixel *pixel = GetPixel(iRow,iCol,iTime);
+  if (pixel) {
+    return (pixel->GetTrack(iTrack));
+  }
+  else {
+    return 0;
+  }
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDmatrix::GetIndex(Int_t iRow, Int_t iCol, Int_t iTime)
+{
+
+  if ((iRow  >= 0) && (iRow  < fRow ) &&
+      (iCol  >= 0) && (iCol  < fCol ) &&
+      (iTime >= 0) && (iTime < fTime)) {
+    return (iTime + iCol * fTime + iRow * fTime * fCol);
+  }
+  else {
+    return -1;
+  }
+
+}
+
+//_____________________________________________________________________________
+AliTRDpixel *AliTRDmatrix::GetPixel(Int_t iRow, Int_t iCol, Int_t iTime)
+{
+
+  Int_t iPixel = GetIndex(iRow,iCol,iTime);
+  if (iPixel < 0) {
+    return NULL;
+  }
+  else {
+    return ((AliTRDpixel *) fPixelArray->At(iPixel));
+  }
+
+}
diff --git a/TRD/AliTRDmatrix.h b/TRD/AliTRDmatrix.h
new file mode 100644 (file)
index 0000000..fdd4db0
--- /dev/null
@@ -0,0 +1,59 @@
+#ifndef TRDmatrix_h
+#define TRDmatrix_h
+
+#include <TObject.h>
+#include <TObjArray.h>
+#include <TH2.h>
+#include <TH3.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include "AliTRDpixel.h"
+
+///////////////////////////////////////////////////////
+//  Stores the pixel-information of one TRD chamber  //
+///////////////////////////////////////////////////////
+
+class AliTRDmatrix : public TObject {
+
+protected:
+  Int_t         fRow;            // Number of pad-rows
+  Int_t         fCol;            // Number of pad-columns
+  Int_t         fTime;           // Number of time buckets
+  Int_t         fPixel;          // Number of pixels
+  Int_t         fSector;         // Sector number
+  Int_t         fChamber;        // Chamber number
+  Int_t         fPlane;          // Plane number
+  TObjArray    *fPixelArray;     // Array of pixels
+
+  virtual Int_t        GetIndex(Int_t iRow, Int_t iCol, Int_t iTime);
+  virtual AliTRDpixel *GetPixel(Int_t iRow, Int_t iCol, Int_t iTime);
+
+public:
+  AliTRDmatrix();
+  AliTRDmatrix(Int_t nRow, Int_t nCol, Int_t nTime, Int_t iSec, Int_t iCha, Int_t iPla);
+  virtual ~AliTRDmatrix();
+
+  virtual void         AddSignal(Int_t iRow, Int_t iCol, Int_t iTime, Float_t signal);
+  virtual Bool_t       AddTrack(Int_t iRow, Int_t iCol, Int_t iTime, Int_t track);
+
+  virtual void         Draw();
+  virtual void         DrawRow(Int_t iRow);
+  virtual void         DrawCol(Int_t iCol);
+  virtual void         DrawTime(Int_t iTime);
+
+  virtual void         SetSignal(Int_t iRow, Int_t iCol, Int_t iTime, Float_t signal);
+  virtual void         SetTrack(Int_t iRow, Int_t iCol, Int_t iTime
+                              , Int_t iTrack, Int_t track);
+
+  virtual Float_t      GetSignal(Int_t iRow, Int_t iCol, Int_t iTime);
+  virtual Int_t        GetTrack(Int_t iRow, Int_t iCol, Int_t iTime, Int_t iTrack);
+
+  virtual Int_t        GetSector()  { return fSector;  };
+  virtual Int_t        GetChamber() { return fChamber; };
+  virtual Int_t        GetPlane()   { return fPlane;   };
+
+  ClassDef(AliTRDmatrix,1)
+
+};
+
+#endif
diff --git a/TRD/AliTRDpixel.cxx b/TRD/AliTRDpixel.cxx
new file mode 100644 (file)
index 0000000..7d121d9
--- /dev/null
@@ -0,0 +1,23 @@
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Contains the information for one TRD pixel                               //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDpixel.h"
+
+ClassImp(AliTRDpixel)
+
+//_____________________________________________________________________________
+AliTRDpixel::AliTRDpixel():TObject()
+{
+  //
+  // Create a TRD pixel
+  // 
+
+  fSignal   = 0;
+  fTrack[0] = 0;
+  fTrack[1] = 0;
+  fTrack[2] = 0;
+
+}
diff --git a/TRD/AliTRDpixel.h b/TRD/AliTRDpixel.h
new file mode 100644 (file)
index 0000000..63e132f
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef TRDpixel_h
+#define TRDpixel_h
+
+#include <TObject.h>
+
+//////////////////////////////////////////////////////
+//  Stores the information for one detector pixel   //
+//////////////////////////////////////////////////////
+
+const Int_t kTrackPixel = 3;
+
+class AliTRDpixel : public TObject {
+
+protected:
+  Float_t      fSignal;                    // Signal sum
+  Int_t        fTrack[kTrackPixel];        // Tracks contributing to this pixel
+
+public:
+  AliTRDpixel();
+  virtual ~AliTRDpixel() {};
+
+  virtual void    SetSignal(Float_t signal)      { fSignal   = signal; };
+  virtual void    SetTrack(Int_t i, Int_t track) { fTrack[i] = track;  };
+
+  virtual Float_t GetSignal()                    { return fSignal;     };
+  virtual Int_t   GetTrack(Int_t i)              { return fTrack[i];   };
+
+  ClassDef(AliTRDpixel,1)
+
+};
+
+#endif
diff --git a/TRD/AliTRDsim.cxx b/TRD/AliTRDsim.cxx
new file mode 100644 (file)
index 0000000..e81e6f2
--- /dev/null
@@ -0,0 +1,310 @@
+#include <stdlib.h>
+
+#include "TH1.h"
+#include "TRandom.h"
+#include "TMath.h"
+
+#include "AliTRDsim.h"
+#include "AliTRDconst.h"
+
+ClassImp(AliTRDsim)
+
+
+const Float_t kD1 = kPeThick / kRaFoils;
+const Float_t kD2 = kRaThick / (kRaFoils + 1);
+
+//Root specials, to be removed
+
+static TH1F *h100, *h101, *h102;
+
+AliTRDsim::AliTRDsim()
+{
+  fNj=200;
+  fIrst=1;
+}
+
+AliTRDsim::AliTRDsim(AliModule* mod, Int_t foil, Int_t gas)
+{
+  Float_t a1, z1, ro1, rad, abs;
+  Float_t a2, z2, ro2;
+  char * name[21];
+  mod->AliGetMaterial(foil, name, a1, z1, ro1, rad, abs);
+  mod->AliGetMaterial(gas, name, a2, z2, ro2, rad, abs);
+  fOmega1 = 28.8*TMath::Sqrt(ro1*z1/a1);
+  fOmega2 = 28.8*TMath::Sqrt(ro2*z2/a2);
+}
+
+void AliTRDsim::trd_sim()
+{
+
+  const Float_t amass[4] = { 5.11e-4,.13957,.4937,.10566 };
+  const Double_t of[4] = { 20.9,24.4,14.27,26.9 };
+  const Double_t og[4] = { .28,.7,.74,.74 };
+  Int_t ifl = 0;
+  Int_t ig = 1;
+  Int_t nev = 1000;
+  Double_t gamma = -10.;
+  
+  /* Local variables */
+  static Float_t temp, pres;
+  static Int_t i, j;
+  static Float_t o, sigma[200];
+  static Float_t trEn[10];
+  static Double_t omega1, omega2;
+  static Float_t am;
+  static Int_t np;
+  static Int_t ipa;
+  
+  /* ***********************************************************************
+   */
+  /*  TRD simulation - multimodule (regular rad.) */
+  /*     after: M. CASTELLANO et al., */
+  /*   COMP. PHYS. COMM. 51 (1988) 431 + COMP. PHYS. COMM. 61 (1990) 395 */
+  
+  /*   17.07.1998 - A.Andronic */
+  /*   08.12.1998 - simplified version */
+  
+  ipa = 0;
+  /* that's electron */
+  am = amass[ipa];
+  omega1 = of[ifl];
+  /* plasma frequency: foil and gap */
+  omega2 = og[ig - 1];
+  if (gamma < -1e5) printf("*** Momentum steps !!! ***\n");
+  if (gamma < 0. && gamma >= -1e5) {
+    gamma = sqrt(gamma * gamma + am * am) / am;
+    printf("*** Gamma (electron) = %f\n",gamma);
+  }
+  temp = 20.;
+  pres = 1.;
+  fBin = 100. /  fNj;
+  /* binsize */
+  fL = 1. - fBin / 2.;
+  fU = fL + 100.;
+  /*  setting the stage ................................... */
+  for (j = 0; j < fNj; ++j) {
+    /* getting the sigma values - for fixed energy values */
+    o = fBin *  j + 1.;
+    /* omega in keV */
+    /* abs. in rad. (1 foi */
+    sigma[j] = fsigmaRad(ifl, ig, o);
+  }
+  printf(" Working...\n");
+  /*  sampling over some events ........................... */
+  for (i = 0; i < nev; ++i) {
+    xtr(gamma, omega1, omega2, ro1, ro2, sigma, np, trEn);
+    /* TR: n, E */
+    h101->Fill(np);
+    /* sample nTR distr. */
+    for (j = 0; j < np; ++j) {
+      h102->Fill(trEn[j], 1. / fBin);
+      /* sample the TR en. distr. */
+    }
+  }
+  /* ------------------------------------------------------------------- */
+  /*      else  !ns steps */
+  /*      enddo  !imod */
+  /* events */
+  h100->Draw();
+  h101->Draw();
+  h102->Draw();
+} /* trd_sim__ */
+
+void AliTRDsim::xtr(Double_t gamma, Double_t omega1, Double_t omega2, Double_t ro1,
+                   Double_t ro2,
+                    Float_t *sigmaRad, Int_t &np, Float_t *trEn)
+{
+    /* Initialized data */
+
+  static Double_t alfa = .0072973;
+  static Double_t pi = 3.14159265;
+  
+  /* Local variables */
+  static Double_t conv, a;
+  static Int_t i, j;
+  static Float_t o, w[200], omega;
+  static Double_t tetan, stemp;
+  static Float_t om;
+  static Double_t sk;
+  static Float_t wn[200];
+  static Double_t cs1, cs2;
+  static Double_t ro11, ro22, aux;
+  static Float_t ntr;
+  static Double_t sum;
+  
+  /************************************************************************
+   ******/
+  /*   TR: number and energy distr. */
+  
+  /* Function Body */
+  sk = kD2 / kD1;
+  /* -------------- starts with the TR spectrum ------------- */
+  
+  stemp = 0.;
+  for (j = 0; j < fNj; ++j) {
+    /* TR spectrum */
+    omega = (fBin *  j + 1.) * 1e3;
+    /* keV->eV */
+    cs1 = omega1 / omega;
+    cs2 = omega2 / omega;
+    ro11 = omega * kD1 * 2.5 * (1. / (gamma * gamma) + cs1*cs1);
+    ro22 = omega * kD1 * 2.5 * (1. / (gamma * gamma) + cs2*cs2);
+    sum = 0.;
+    for (i = 0; i < 10; ++i) {
+/* 30 - it matters a bit */
+      tetan = (pi * 2. *  (i+1) - (ro11 + sk * ro22)) / (sk + 1.);
+      if (tetan < 0.) {
+       tetan = 0.;
+      }
+      aux = 1. / (ro11 + tetan) - 1. / (ro22 + tetan);
+      a = tetan * (aux * aux) * (1. - cos(ro11 + tetan));
+      sum += a;
+    }
+    o = omega * .001;
+    /* eV->keV */
+    conv = 1. - exp(-kRaFoils * sigmaRad[j]);
+    w[j] = alfa * 4. / (sigmaRad[j] * (sk + 1.)) * conv * sum;
+    /* dW/domega */
+    wn[j] = w[j] / o;
+    /* dN/domega */
+    stemp += wn[j];
+    if (fIrst == 1) {
+      h100->Fill(o, wn[j]);
+      /* double precision not accepted */
+    }
+  }
+  /* -------------- done with the spectrum ------------- */
+  /* j (omega spectrum) */
+  ntr = stemp * fBin;
+  /* <nTR> (binsize corr.) */
+  om = h100->GetMean();
+  /* <Etr> */
+  if (fIrst == 1) {
+    /* prints the production */
+    printf(" Produced TR - <n>, <E>: %5.2f  %6.2f  KeV\n",ntr,om);
+    fIrst = 0;
+  }
+  /* prob. distr. */
+  np = gRandom->Poisson(ntr);
+  /* Np TR photons Poiss distr. from mean */
+  for (j = 0; j < np; ++j) {
+    /* TR energy (binsize corr.) */
+    trEn[j] = hisran(wn, fNj, fL, fBin);
+    /* their energy */
+  }
+}
+
+Float_t AliTRDsim::fsigmaRad(Float_t ro1, Float_t ro2, Int_t ig, Float_t o)
+{
+
+  /* Local variables */
+  static Float_t pres;
+  static Double_t mumu;
+  static Int_t j;
+  static Double_t t;
+  static Int_t i1, i2;
+  static Double_t x1;
+  static Double_t mu1, mu2, deo, omf[36], omg[36], muf[36], mug[36];
+  
+  static Bool_t first = kTRUE;
+  
+  /* cccccccccccccccccccccccccccccccccccccccccccc */
+  /*    calculates sigma for radiator - one foil+one gap */
+  
+  if(first) {
+    FILE* inp = fopen("po.x","r");
+    for (j=0;j<36;++j) {
+      fscanf(inp,"%lf %lf %lf",&omf[j],&muf[j],&mumu);
+    }
+    fclose(inp);
+    inp = fopen("he.x","r");
+    for (j=0;j<36;++j) {
+      fscanf(inp,"%lf %lf %lf",&omg[j],&mug[j],&mumu);
+    }
+    fclose(inp);
+    first=kFALSE;
+  }
+  /* first */
+  x1 = o * .001;
+  /* keV->MeV */
+  if (x1 >= .001) {
+    locate(omf, 36, x1, i1, deo);
+    mu1 = muf[i1] - deo * (muf[i1] - muf[i1+1]) / (omf[i1+1] - omf[i1]);
+    locate(omg, 36, x1, i2, deo);
+    mu2 = mug[i2] - deo * (mug[i2] - mug[i2+1]) / (omg[i2+1] - omg[i2]);
+    t = 273.16;
+    /* gases at 0 C */
+    return (mu1*ro1*kD1+mu2*293.16/t * ro2*kD2)/1e4;
+    /* mu */
+  } else {
+    return 1e6;
+  }
+} 
+
+Int_t AliTRDsim::locate(Double_t *xv, Int_t n, Double_t xval, 
+                      Int_t &kl, Double_t &dx)
+{
+  /* -------------------------------------------------------------- */
+  /*  locates a point (xval) in a 1-dim grid (xv(n)) --> iloc,dx,ier */
+  /* -------------------------------------------------------------- */
+  /* Function Body */
+  if (xval >= xv[n-1]) return 1;
+  if (xval < xv[0]) return -1;
+  Int_t km,kh=n-1;
+  kl=0;
+  while(kh-kl>1) if(xval<xv[km=(kl+kh)/2]) kh=km; else kl=km;
+  if(xval<xv[kl] || xval > xv[kl+1] || kl >= n-1) {
+    printf("locate failed xv[%d] %f xval %f xv[%d] %f!!!\n",
+          kl,xv[kl],xval,kl+1,xv[kl+1]);
+    exit(1);
+  }
+  dx=xval-xv[kl];
+  return 0;
+}
+
+Float_t AliTRDsim::hisran(Float_t *y, Int_t n, Float_t xlo, Float_t xwid)
+{
+    /* Local variables */
+    Float_t yinv, ytot=0;
+    Int_t i;
+    Float_t yr;
+
+/*         SUBROUTINE TO GENERATE RANDOM NUMBERS */
+/*         ACCORDING TO AN EMPIRICAL DISTRIBUTION */
+/*         SUPPLIED BY THE USER IN THE FORM OF A HISTOGRAM */
+/*         F. JAMES,    MAY, 1976 */
+
+    if (y[n-1] != 1.) {
+
+/*         INITIALIZE HISTOGRAM TO FORM CUMULATIVE DISTRIBUTION */
+
+      ytot = 0.;
+      for (i = 0; i < n; ++i) {
+       if (y[i] < 0.) {
+         printf("hisran: found value y[%d] = %f\n",i,y[i]);
+         exit(1);
+       }
+       ytot += y[i];
+       y[i] = ytot;
+      }
+      if (ytot <= 0.) {
+       printf("hisran: total probability %f < 0\n",ytot);
+       exit(1);
+      }
+      yinv = 1. / ytot;
+      for (i = 0; i < n-1; ++i) {
+       y[i] *= yinv;
+      }
+      y[n-1] = 1.;
+    }
+/*         NOW GENERATE RANDOM NUMBER BETWEEN 0 AND ONE */
+    yr = gRandom->Rndm();
+/*         AND TRANSFORM IT INTO THE CORRESPONDING X-VALUE */
+    if(yr<=y[0]) return xlo + xwid * (yr / y[0]);
+    else {
+      Int_t km,kl=0,kh=n-1;
+      while(kh-kl>1) if(yr<y[km=(kl+kh)/2]) kh=km; else kl=km;
+      return xlo + xwid * (kl + (yr - y[kl]) / (y[kl + 1] - y[kl]));
+    }
+}
+
diff --git a/TRD/AliTRDsim.h b/TRD/AliTRDsim.h
new file mode 100644 (file)
index 0000000..1c1e6d8
--- /dev/null
@@ -0,0 +1,27 @@
+#ifndef AliTRDsim_H
+#define AliTRDsim_H
+
+#include "TObject.h"
+
+class AliTRDsim : public TObject {
+private:
+  Int_t fNj;            // Number of channel of the histogram
+  Float_t fU, fL, fBin;
+  Double_t fRo1, fRo2, fOmega1, fOmega2;
+  Int_t fIrst;
+
+public:
+  AliTRDsim();
+  virtual ~AliTRDsim() {}
+  virtual void trd_sim();
+  virtual void xtr(Double_t gamma, Double_t omega1, Double_t omega2,
+                   Float_t *sigmaRad, Int_t &np, Float_t *trEn);
+  virtual Float_t fsigmaRad(Int_t ifl, Int_t ig, Float_t o);
+  virtual Int_t locate(Double_t *xv, Int_t n, Double_t xval, 
+                      Int_t &kl, Double_t &dx);
+  virtual Float_t hisran(Float_t *y, Int_t n, Float_t xlo, Float_t xwid);
+
+  
+  ClassDef(AliTRDsim,1)  //Base class for all Alice hits
+};
+#endif
index 3f76921..cb23f49 100644 (file)
 
 #include <TMath.h>
 #include <TVector.h>
+#include <TRandom.h>
 
 #include "AliTRDv2.h"
+#include "AliTRDmatrix.h"
 #include "AliRun.h"
 #include "AliMC.h"
 #include "AliConst.h"
@@ -29,31 +31,54 @@ AliTRDv2::AliTRDv2(const char *name, const char *title)
   // Standard constructor for Transition Radiation Detector version 2
   //
 
-  fIdSens      = 0;
+  fIdSens       = 0;
 
-  fIdSpace1    = 0;
-  fIdSpace2    = 0;
-  fIdSpace3    = 0;
+  fIdSpace1     = 0;
+  fIdSpace2     = 0;
+  fIdSpace3     = 0;
 
-  fIdChamber1  = 0;
-  fIdChamber2  = 0;
-  fIdChamber3  = 0;
+  fIdChamber1   = 0;
+  fIdChamber2   = 0;
+  fIdChamber3   = 0;
 
-  fSensSelect  = 0;
-  fSensPlane   = 0;
-  fSensChamber = 0;
-  fSensSector  = 0;
+  fSensSelect   = 0;
+  fSensPlane    = 0;
+  fSensChamber  = 0;
+  fSensSector   = 0;
 
-  fDeltaE      = NULL;
+  for (Int_t iplan = 0; iplan < kNplan; iplan++) {
+    for (Int_t icham = 0; icham < kNcham; icham++) {
+      fRowMax[iplan][icham] = 0;
+    }
+    fColMax[iplan] = 0;
+  }
+  fTimeMax      = 0;
+
+  fRowPadSize   = 0;
+  fColPadSize   = 0;
+  fTimeBinSize  = 0;
+
+  fGasGain      = 0;
+  fNoise        = 0;
+  fChipGain     = 0;
+  fADCoutRange  = 0;
+  fADCinRange   = 0;
+  fADCthreshold = 0;
+
+  fDiffusionT   = 0;
+  fDiffusionL   = 0;
+
+  fDeltaE       = NULL;
 
   SetBufferSize(128000);
 
 }
 
+//_____________________________________________________________________________
 AliTRDv2::~AliTRDv2()
 {
 
-  if (fDeltaE)  delete fDeltaE;
+  if (fDeltaE) delete fDeltaE;
 
 }
  
@@ -98,16 +123,329 @@ void AliTRDv2::CreateMaterials()
 }
 
 //_____________________________________________________________________________
-void AliTRDv2::Init() 
+void AliTRDv2::Diffusion(Float_t driftlength, Float_t *xyz)
 {
   //
-  // Initialise Transition Radiation Detector after geometry has been built
+  // Applies the diffusion smearing to the position of a single electron
   //
 
-  // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
-  const Float_t kPoti = 12.1;
-  // Maximum energy (50 keV);
-  const Float_t kEend = 50000.0;
+  if ((driftlength >        0) && 
+      (driftlength < kDrThick)) {
+    Float_t driftSqrt = TMath::Sqrt(driftlength);
+    Float_t sigmaT = driftSqrt * fDiffusionT;
+    Float_t sigmaL = driftSqrt * fDiffusionL;
+    xyz[0] = gRandom->Gaus(xyz[0], sigmaL);
+    xyz[1] = gRandom->Gaus(xyz[1], sigmaT);
+    xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
+  }
+  else {
+    xyz[0] = 0.0;
+    xyz[1] = 0.0;
+    xyz[2] = 0.0;
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDv2::Hits2Digits()
+{
+  //
+  // Creates TRD digits from hits. This procedure includes the following:
+  //      - Diffusion
+  //      - 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 galice.root file with TRD-hits use the
+  // digitsCreate.C macro.
+  //
+
+  printf(" Start creating digits\n");
+
+  ///////////////////////////////////////////////////////////////
+  // Parameter 
+  ///////////////////////////////////////////////////////////////
+
+  // Converts number of electrons to fC
+  const Float_t el2fC = 1.602E-19 * 1.0E15; 
+
+  ///////////////////////////////////////////////////////////////
+
+  Int_t nBytes = 0;
+
+  AliTRDhit *TRDhit;
+
+  // 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)
+  Float_t row0[kNplan][kNcham];
+  for (Int_t iplan = 0; iplan < kNplan; iplan++) {
+    row0[iplan][0] = -fClengthI[iplan]/2. - fClengthM[iplan] - fClengthO[iplan] 
+                   + kCcthick; 
+    row0[iplan][1] = -fClengthI[iplan]/2. - fClengthM[iplan]                    
+                   + kCcthick;
+    row0[iplan][2] = -fClengthI[iplan]/2.                                       
+                   + kCcthick;
+    row0[iplan][3] =  fClengthI[iplan]/2.                                       
+                   + kCcthick; 
+    row0[iplan][4] =  fClengthI[iplan]/2. + fClengthM[iplan]                    
+                   + kCcthick; 
+  }
+  // The pad column (rphi-direction)  
+  Float_t col0[kNplan];
+  for (Int_t iplan = 0; iplan < kNplan; iplan++) {
+    col0[iplan]    = -fCwidth[iplan]/2. + kCcthick;
+  }
+  // The time bucket
+  Float_t time0[kNplan];
+  for (Int_t iplan = 0; iplan < kNplan; iplan++) {
+    time0[iplan]   = kRmin + kCcframe/2. + kDrZpos - 0.5 * kDrThick
+                           + iplan * (kCheight + kCspace);
+  } 
+
+  // Get the pointer to the hit tree
+  TTree *HitTree    = gAlice->TreeH();
+  // Get the pointer to the digits tree
+  TTree *DigitsTree = gAlice->TreeD();
+
+  // 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 (fSensChamber) chamEnd = chamBeg = fSensChamber;
+  Int_t planBeg = 0;
+  Int_t planEnd = kNplan;
+  if (fSensPlane)   planEnd = planBeg = fSensPlane;
+  Int_t sectBeg = 0;
+  Int_t sectEnd = kNsect;
+  if (fSensSector)  sectEnd = sectBeg = fSensSector;
+
+  // 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++) {
+
+        printf(" Digitizing chamber %d, plane %d, sector %d\n"
+              ,icham+1,iplan+1,isect+1);
+
+        // Create a detector matrix to keep the signal and track numbers
+        AliTRDmatrix *matrix = new AliTRDmatrix(fRowMax[iplan][icham]
+                                               ,fColMax[iplan]
+                                               ,fTimeMax
+                                               ,isect+1,icham+1,iplan+1);
+
+        // 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 = fHits->GetEntriesFast();
+
+          // Loop through the TRD hits  
+          for (Int_t iHit = 0; iHit < nHit; iHit++) {
+
+            if (!(TRDhit = (AliTRDhit *) fHits->UncheckedAt(iHit))) 
+              continue;
+
+            Float_t x       = TRDhit->fX;
+            Float_t y       = TRDhit->fY;
+            Float_t z       = TRDhit->fZ;
+            Float_t q       = TRDhit->fQ;
+            Int_t   track   = TRDhit->fTrack;
+            Int_t   plane   = TRDhit->fPlane;
+            Int_t   sector  = TRDhit->fSector;
+            Int_t   chamber = TRDhit->fChamber;        
+
+            if ((sector  != isect+1) ||
+                (plane   != iplan+1) ||
+                (chamber != icham+1)) 
+              continue;
+
+            // Rotate the sectors on top of each other
+            Float_t phi  = 2.0 * kPI /  (Float_t) kNsect 
+                               * ((Float_t) sector - 0.5);
+            Float_t xRot = -x * TMath::Cos(phi) + y * TMath::Sin(phi);
+            Float_t yRot =  x * TMath::Sin(phi) + y * TMath::Cos(phi);
+            Float_t zRot =  z;
+
+            // The hit position in pad coordinates (center pad)
+            // The pad row (z-direction)
+            Int_t rowH  = (Int_t) ((zRot -  row0[iplan][icham]) / fRowPadSize);
+            // The pad column (rphi-direction)  
+            Int_t colH  = (Int_t) ((yRot -  col0[iplan]       ) / fColPadSize);
+            // The time bucket
+            Int_t timeH = (Int_t) ((xRot - time0[iplan]       ) / fTimeBinSize);
+
+            // Array to sum up the signal in a box surrounding the
+            // hit postition
+            const Int_t timeBox = 5;
+            const Int_t  colBox = 7;
+            const Int_t  rowBox = 5;
+            Float_t signalSum[rowBox][colBox][timeBox];
+            for (Int_t iRow  = 0;  iRow <  rowBox; iRow++ ) {
+              for (Int_t iCol  = 0;  iCol <  colBox; iCol++ ) {
+                for (Int_t 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++) {
+
+              // Apply the diffusion smearing
+              Float_t driftlength = xRot - time0[iplan];
+              Float_t xyz[3];
+              xyz[0] = xRot;
+              xyz[1] = yRot;
+              xyz[2] = zRot;
+              Diffusion(driftlength,xyz);
+
+              // At this point absorption effects that depend on the 
+             // driftlength could be taken into account.              
+
+              // 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[iplan][icham]) / fRowPadSize);
+              Int_t  rowD =  rowH -  rowE;
+              // The pad column (rphi-direction)
+              Int_t  colE = (Int_t) ((xyz[1] -  col0[iplan]       ) / fColPadSize);
+              Int_t  colD =  colH -  colE;
+              // The time bucket
+              Int_t timeE = (Int_t) ((xyz[0] - time0[iplan]       ) / fTimeBinSize);
+              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[iplan] - (colE + 0.5) * fColPadSize) 
+                           / fColPadSize;
+
+              // 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);
+              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 (Int_t iRow  = 0;  iRow <  rowBox; iRow++ ) {
+              for (Int_t iCol  = 0;  iCol <  colBox; iCol++ ) {
+                for (Int_t 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("More than three tracks in a pixel!\n");
+                 }
+
+               }
+             }
+           }
+
+          }
+
+       }
+
+        // Create the hits for this chamber
+        for (Int_t iRow  = 0; iRow  <  fRowMax[iplan][icham]; iRow++ ) {
+          for (Int_t iCol  = 0; iCol  <  fColMax[iplan]         ; iCol++ ) {
+            for (Int_t iTime = 0; iTime < fTimeMax                ; iTime++) {         
+
+              Float_t signalAmp = matrix->GetSignal(iRow,iCol,iTime);
+
+              // Add the noise
+              signalAmp  = TMath::Max(gRandom->Gaus(signalAmp,fNoise),0.0);
+             // Convert to fC
+              signalAmp *= el2fC;
+              // Convert to mV
+              signalAmp *= fChipGain;
+             // Convert to ADC counts
+              Int_t adc  = (Int_t) (signalAmp * (fADCoutRange / fADCinRange));
+
+             // Apply threshold on ADC value
+              if (adc > fADCthreshold) {
+
+                Int_t trackSave[3];
+                for (Int_t ii = 0; ii < 3; ii++) {
+                  trackSave[ii] = matrix->GetTrack(iRow,iCol,iTime,ii);
+               }
+
+                Int_t digits[7];
+                digits[0] = matrix->GetSector();
+                digits[1] = matrix->GetChamber();
+                digits[2] = matrix->GetPlane();
+                digits[3] = iRow;
+                digits[4] = iCol;
+                digits[5] = iTime;
+                digits[6] = adc;
+
+               // Add this digit to the TClonesArray
+                AddDigit(trackSave,digits);
+
+             }
+
+           }
+         }
+       }
+
+       // Clean up
+        delete matrix;
+
+      }
+    }
+  }
+
+  // Fill the digits-tree
+  DigitsTree->Fill();
+
+}
+
+//_____________________________________________________________________________
+void AliTRDv2::Init() 
+{
+  //
+  // Initialise Transition Radiation Detector after geometry has been built.
+  // Includes the default settings of all parameter for the digitization.
+  //
 
   AliTRD::Init();
 
@@ -118,9 +456,13 @@ void AliTRDv2::Init()
   if (fSensSector)
     printf("          Only sector %d is sensitive\n",fSensSector);
 
-  for(Int_t i = 0; i < 80; i++) printf("*");
+  for (Int_t i = 0; i < 80; i++) printf("*");
   printf("\n");
 
+  // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
+  const Float_t kPoti = 12.1;
+  // Maximum energy (50 keV);
+  const Float_t kEend = 50000.0;
   // Ermilova distribution for the delta-ray spectrum
   Float_t Poti = TMath::Log(kPoti);
   Float_t Eend = TMath::Log(kEend);
@@ -139,6 +481,95 @@ void AliTRDv2::Init()
   fIdChamber2 = gMC->VolId("UCIM");
   fIdChamber3 = gMC->VolId("UCII");
 
+  // The default pad dimensions
+  if (!(fRowPadSize))  fRowPadSize  = 4.5;
+  if (!(fColPadSize))  fColPadSize  = 1.0;
+  if (!(fTimeBinSize)) fTimeBinSize = 0.1;
+
+  // The maximum number of pads
+  for (Int_t iplan = 0; iplan < kNplan; iplan++) {
+    // Rows 
+    fRowMax[iplan][0] = 1 + TMath::Nint((fClengthO[iplan] - 2. * kCcthick) 
+                                                          / fRowPadSize - 0.5);
+    fRowMax[iplan][1] = 1 + TMath::Nint((fClengthM[iplan] - 2. * kCcthick) 
+                                                          / fRowPadSize - 0.5);
+    fRowMax[iplan][2] = 1 + TMath::Nint((fClengthI[iplan] - 2. * kCcthick) 
+                                                          / fRowPadSize - 0.5);
+    fRowMax[iplan][3] = 1 + TMath::Nint((fClengthM[iplan] - 2. * kCcthick) 
+                                                          / fRowPadSize - 0.5);
+    fRowMax[iplan][4] = 1 + TMath::Nint((fClengthO[iplan] - 2. * kCcthick) 
+                                                          / fRowPadSize - 0.5);
+    // Columns
+    fColMax[iplan]    = 1 + TMath::Nint((fCwidth[iplan]   - 2. * kCcthick) 
+                                                          / fColPadSize - 0.5);
+  }
+  // Time buckets
+  fTimeMax = 1 + TMath::Nint(kDrThick / fTimeBinSize - 0.5);
+
+  // The default parameter for the digitization
+  if (!(fGasGain))      fGasGain      = 2.0E3;
+  if (!(fNoise))        fNoise        = 3000.;
+  if (!(fChipGain))     fChipGain     = 10.;
+  if (!(fADCoutRange))  fADCoutRange  = 255.;
+  if (!(fADCinRange))   fADCinRange   = 2000.;
+  if (!(fADCthreshold)) fADCthreshold = 0;
+
+  // Transverse and longitudinal diffusion coefficients (Xe/Isobutane)
+  if (!(fDiffusionT))   fDiffusionT   = 0.060;
+  if (!(fDiffusionL))   fDiffusionL   = 0.017;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDv2::MakeBranch(Option_t* option)
+{
+  //
+  // Create Tree branches for the TRD digits.
+  //
+
+  Int_t  buffersize = 4000;
+  Char_t branchname[10];
+
+  sprintf(branchname,"%s",GetName());
+
+  AliDetector::MakeBranch(option); 
+
+  Char_t *D = strstr(option,"D");
+  if (fDigits && gAlice->TreeD() && D) {
+    gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
+    printf("Making Branch %s for digits\n",branchname);
+  }
+
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDv2::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)));
+
+  //TF1 *funPR = new TF1("funPR","[0]*([1]+exp(-x*x /(2.*[2])))",-3,3);
+  //funPR->SetParameter(0,aa );
+  //funPR->SetParameter(1,bb );
+  //funPR->SetParameter(2,cc2);
+  //
+  //Float_t pr = funPR->Eval(distance);
+  //
+  //delete funPR;
+
+  return (pr);
+
 }
 
 //_____________________________________________________________________________
@@ -364,14 +795,6 @@ Double_t AliTRDv2::BetheBloch(Double_t bg)
   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
   //
 
-  // The parameters have been adjusted to Xe-data found in:
-  // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
-  //const Double_t kP1 = 0.76176E-1;
-  //const Double_t kP2 = 10.632;
-  //const Double_t kP3 = 3.17983E-6;
-  //const Double_t kP4 = 1.8631;
-  //const Double_t kP5 = 1.9479;
-
   // This parameters have been adjusted to averaged values from GEANT
   const Double_t kP1 = 7.17960e-02;
   const Double_t kP2 = 8.54196;
@@ -379,6 +802,14 @@ Double_t AliTRDv2::BetheBloch(Double_t bg)
   const Double_t kP4 = 5.30972;
   const Double_t kP5 = 2.83798;
 
+  // This parameters have been adjusted to Xe-data found in:
+  // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
+  //const Double_t kP1 = 0.76176E-1;
+  //const Double_t kP2 = 10.632;
+  //const Double_t kP3 = 3.17983E-6;
+  //const Double_t kP4 = 1.8631;
+  //const Double_t kP5 = 1.9479;
+
   if (bg > 0) {
     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
     Double_t aa = TMath::Power(yy,kP4);
@@ -395,7 +826,7 @@ Double_t AliTRDv2::BetheBloch(Double_t bg)
 Double_t Ermilova(Double_t *x, Double_t *)
 {
   //
-  // Calculates the delta-ray energy distribution according to Ermilova
+  // Calculates the delta-ray energy distribution according to Ermilova.
   // Logarithmic scale !
   //
 
index 22c1022..1798d19 100644 (file)
@@ -15,38 +15,89 @@ class AliTRDv2 : public AliTRD {
 public:
   AliTRDv2() {}
   AliTRDv2(const char *name, const char *title);
-  virtual      ~AliTRDv2();
-  virtual void  CreateGeometry();
-  virtual void  CreateMaterials();
-  virtual Int_t IsVersion() const {return 2;}
-  virtual void  StepManager();
-  virtual void  SetSensPlane(Int_t iplane = 0);
-  virtual void  SetSensChamber(Int_t ichamber = 0);
-  virtual void  SetSensSector(Int_t isector = 0);
-  virtual void  Init();
+  virtual        ~AliTRDv2();
+  virtual void    CreateGeometry();
+  virtual void    CreateMaterials();
+  virtual Int_t   IsVersion() const {return 2;}
+  virtual void    MakeBranch(Option_t* option);
+  virtual void    StepManager();
+  virtual void    SetSensPlane(Int_t iplane = 0);
+  virtual void    SetSensChamber(Int_t ichamber = 0);
+  virtual void    SetSensSector(Int_t isector = 0);
+  virtual void    Init();
+  virtual void    Hits2Digits(); 
+
+  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 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    SetDiffusionT(Float_t diff)         { fDiffusionT   = diff;     };
+  virtual void    SetDiffusionL(Float_t diff)         { fDiffusionL   = diff;     };
+
+  virtual Float_t GetRowPadSize()                     { return fRowPadSize;   };
+  virtual Float_t GetColPadSize()                     { return fColPadSize;   };
+  virtual Float_t GetTimeBinSize()                    { return fTimeBinSize;  };
+
+  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 Int_t   GetRowMax(Int_t iplan, Int_t icham) { return fRowMax[iplan-1][icham-1]; };
+  virtual Int_t   GetColMax(Int_t iplan)              { return fColMax[iplan-1];          };
+  virtual Int_t   GetTimeMax()                        { return fTimeMax;                  };
 
 protected:
-  Int_t        fIdSens;          // Sensitive volume identifier
+  Int_t        fIdSens;                 // Sensitive volume identifier
+
+  Int_t        fIdSpace1;               // Spaceframe volume identifier
+  Int_t        fIdSpace2;               // 
+  Int_t        fIdSpace3;               // 
+
+  Int_t        fIdChamber1;             // Driftchamber volume identifier
+  Int_t        fIdChamber2;             // 
+  Int_t        fIdChamber3;             // 
+
+  Int_t        fSensSelect;             // Switch to select only parts of the detector
+  Int_t        fSensPlane;              // Sensitive detector plane
+  Int_t        fSensChamber;            // Sensitive detector chamber
+  Int_t        fSensSector;             // Sensitive detector sector
 
-  Int_t        fIdSpace1;        // Spaceframe volume identifier
-  Int_t        fIdSpace2;        // 
-  Int_t        fIdSpace3;        // 
+  Int_t        fRowMax[kNplan][kNcham]; // Number of pad-rows
+  Int_t        fColMax[kNplan];         // Number of pad-columns
+  Int_t        fTimeMax;                // Number of time buckets
 
-  Int_t        fIdChamber1;      // Driftchamber volume identifier
-  Int_t        fIdChamber2;      // 
-  Int_t        fIdChamber3;      // 
+  Float_t      fRowPadSize;             // Pad size in z-direction
+  Float_t      fColPadSize;             // Pad size in rphi-direction
+  Float_t      fTimeBinSize;            // Size of the time buckets
 
-  Int_t        fSensSelect;      // Switch to select only parts of the detector
-  Int_t        fSensPlane;       // Sensitive detector plane
-  Int_t        fSensChamber;     // Sensitive detector chamber
-  Int_t        fSensSector;      // Sensitive detector sector
+  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
+  Float_t      fDiffusionT;             // Diffusion in transverse direction
+  Float_t      fDiffusionL;             // Diffusion in logitudinal direction
 
 private:
   virtual Double_t BetheBloch(Double_t bg);
+  virtual void     Diffusion(Float_t driftlength, Float_t *xyz);
+  virtual Float_t  PadResponse(Float_t x);
 
-  TF1         *fDeltaE;          // Energy distribution of the delta-electrons
-  
-  ClassDef(AliTRDv2,1)           // Transition Radiation Detector version 2
+  TF1         *fDeltaE;                 // Energy distribution of the delta-electrons
+   
+  ClassDef(AliTRDv2,1)                  // Transition Radiation Detector version 2 (slow simulator)
 
 };
 
index 7d86e70..6877bbe 100644 (file)
@@ -11,7 +11,8 @@ PACKAGE             = TRD
 
 # C++ sources
 
-SRCS          = AliTRD.cxx     AliTRDv0.cxx      AliTRDv1.cxx    AliTRDv2.cxx 
+SRCS          = AliTRD.cxx AliTRDv0.cxx AliTRDv1.cxx AliTRDv2.cxx \
+               AliTRDmatrix.cxx AliTRDpixel.cxx
 
 # C++ Headers
 
index eee4f57..c4c6974 100644 (file)
@@ -9,5 +9,8 @@
 #pragma link C++ class  AliTRDv1;
 #pragma link C++ class  AliTRDv2;
 #pragma link C++ class  AliTRDhit;
+#pragma link C++ class  AliTRDdigit;
+#pragma link C++ class  AliTRDpixel;
+#pragma link C++ class  AliTRDmatrix;
 
 #endif
diff --git a/TRD/digitsAna.C b/TRD/digitsAna.C
new file mode 100644 (file)
index 0000000..9ba651f
--- /dev/null
@@ -0,0 +1,115 @@
+{
+
+/////////////////////////////////////////////////////////////////////////
+//
+// Example macro for the analysis of the TRD digits and the use
+// of the AliTRDmatrix class.
+//
+/////////////////////////////////////////////////////////////////////////
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+  }
+
+  // Input file name
+  Char_t *alifile = "galice_v2.root"; 
+
+  // Event number
+  Int_t   nEvent  = 0;
+
+  // Define the objects
+  AliTRDv2     *TRD;
+  TClonesArray *TRDDigits;
+  AliTRDdigit  *OneTRDDigit;
+
+  // 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
+  if (!gAlice) {
+    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 *DigitsTree = gAlice->TreeD();
+
+  // Get the pointer to the detector classes
+  TRD = (AliTRDv2*) gAlice->GetDetector("TRD");
+  // Get the pointer to the hit container
+  if (TRD) TRDDigits = TRD->Digits();
+
+  // Define the detector matrix for one chamber (Sector 6, Chamber 3, Plane 1)
+  const Int_t iSec = 6;
+  const Int_t iCha = 3;
+  const Int_t iPla = 1;
+  Int_t  rowMax = TRD->GetRowMax(iPla,iCha);
+  Int_t  colMax = TRD->GetColMax(iPla);
+  Int_t timeMax = TRD->GetTimeMax();
+  cout << " rowMax = "  << rowMax
+       << " colMax = "  << colMax
+       << " timeMax = " << timeMax << endl;
+  AliTRDmatrix *TRDMatrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
+
+  Int_t nEntries = DigitsTree->GetEntries();
+  cout << "Number of entries in digits tree = " << nEntries << endl; 
+
+  // 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
+    gAlice->ResetDigits();
+    nbytes += DigitsTree->GetEvent(iEntry);
+
+    // Get the number of digits in the detector 
+    Int_t nTRDDigits = TRDDigits->GetEntriesFast();
+    cout << " nTRDDigits = " << nTRDDigits << endl;    
+
+    // Loop through all TRD digits
+    for (Int_t iTRDDigits = 0; iTRDDigits < nTRDDigits; iTRDDigits++) {
+
+      // Get the information for this digit
+      OneTRDDigit = (AliTRDdigit*) TRDDigits->UncheckedAt(iTRDDigits);
+      Int_t signal    = OneTRDDigit->fSignal;
+      Int_t   sector  = OneTRDDigit->fSector;
+      Int_t   chamber = OneTRDDigit->fChamber;
+      Int_t   plane   = OneTRDDigit->fPlane;
+      Int_t   row     = OneTRDDigit->fRow;
+      Int_t   col     = OneTRDDigit->fCol;
+      Int_t   time    = OneTRDDigit->fTime;
+
+      // Fill the detector matrix
+      if (signal > 1) {
+        TRDMatrix->SetSignal(row,col,time,signal);
+      }
+
+    }
+
+  }
+
+  // Display the detector matrix
+  TRDMatrix->Draw();
+  TRDMatrix->DrawRow(18);
+  TRDMatrix->DrawCol(58);
+  TRDMatrix->DrawTime(20);
+
+}
diff --git a/TRD/digitsCreate.C b/TRD/digitsCreate.C
new file mode 100644 (file)
index 0000000..590549a
--- /dev/null
@@ -0,0 +1,57 @@
+{
+
+/////////////////////////////////////////////////////////////////////////
+//
+// Creates the digits from the hit information. An additional hit-tree
+// is added to the input file.
+//
+/////////////////////////////////////////////////////////////////////////
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+  }
+
+  // Input (and output) file name
+  Char_t *alifile = "galice_v2.root"; 
+
+  // Event number
+  Int_t   nEvent  = 0;
+
+  // Connect the AliRoot file containing Geometry, Kine, and Hits
+  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
+  if (!gafl) {
+    cout << "Open the ALIROOT-file " << alifile << endl;
+    gafl = new TFile(alifile,"UPDATE");
+  }
+  else {
+    cout << alifile << " is already open" << endl;
+  }
+
+  // Get AliRun object from file or create it if not on file
+  if (!gAlice) {
+    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 detector class
+  AliTRDv2 *TRD = (AliTRDv2*) gAlice->GetDetector("TRD");
+
+  // Create the digitd and fill the digits-tree
+  TRD->Hits2Digits();
+
+  // Write the new tree into the input file
+  cout << "Entries in hit tree = " << gAlice->TreeD()->GetEntries()) << endl;
+  Char_t treeName[7];
+  sprintf(treeName,"TreeD%d",nEvent);
+  gAlice->TreeD()->Write(treeName);
+
+}