fIshunt = 0;
fGasMix = 0;
+ fHits = 0;
+ fDigits = 0;
// The chamber dimensions
for (Int_t iplan = 0; iplan < kNplan; iplan++) {
// Allocate the hit array
fHits = new TClonesArray("AliTRDhit", 405);
-
+
+ // Allocate the digits array
+ fDigits = new TClonesArray("AliTRDdigit",10000);
+
fIshunt = 0;
fGasMix = 0;
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)
{
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];
+
+}
#include "AliDetector.h"
#include "AliHit.h"
+#include "AliDigit.h"
#include "AliTRDconst.h"
//_____________________________________________________________________________
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();
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
//
-// 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;
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
--- /dev/null
+///////////////////////////////////////////////////////////////////////////////
+// //
+// 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));
+ }
+
+}
--- /dev/null
+#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
--- /dev/null
+///////////////////////////////////////////////////////////////////////////////
+// //
+// 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;
+
+}
--- /dev/null
+#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
--- /dev/null
+#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]));
+ }
+}
+
--- /dev/null
+#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
#include <TMath.h>
#include <TVector.h>
+#include <TRandom.h>
#include "AliTRDv2.h"
+#include "AliTRDmatrix.h"
#include "AliRun.h"
#include "AliMC.h"
#include "AliConst.h"
// 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;
}
}
//_____________________________________________________________________________
-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();
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);
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);
+
}
//_____________________________________________________________________________
// 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;
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);
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 !
//
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)
};
# 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
#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
--- /dev/null
+{
+
+/////////////////////////////////////////////////////////////////////////
+//
+// 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);
+
+}
--- /dev/null
+{
+
+/////////////////////////////////////////////////////////////////////////
+//
+// 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);
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff