From 99d5402e0113abf49efc553ced49b525a9143b5d Mon Sep 17 00:00:00 2001 From: fca Date: Wed, 22 Sep 1999 16:59:12 +0000 Subject: [PATCH] Update TRD code from C.Blume --- TRD/AliTRD.cxx | 60 +++++- TRD/AliTRD.h | 28 ++- TRD/AliTRDconst.h | 9 +- TRD/AliTRDmatrix.cxx | 341 ++++++++++++++++++++++++++++++ TRD/AliTRDmatrix.h | 59 ++++++ TRD/AliTRDpixel.cxx | 23 ++ TRD/AliTRDpixel.h | 32 +++ TRD/AliTRDsim.cxx | 310 +++++++++++++++++++++++++++ TRD/AliTRDsim.h | 27 +++ TRD/AliTRDv2.cxx | 489 ++++++++++++++++++++++++++++++++++++++++--- TRD/AliTRDv2.h | 97 +++++++-- TRD/Makefile | 3 +- TRD/TRDLinkDef.h | 3 + TRD/digitsAna.C | 115 ++++++++++ TRD/digitsCreate.C | 57 +++++ 15 files changed, 1589 insertions(+), 64 deletions(-) create mode 100644 TRD/AliTRDmatrix.cxx create mode 100644 TRD/AliTRDmatrix.h create mode 100644 TRD/AliTRDpixel.cxx create mode 100644 TRD/AliTRDpixel.h create mode 100644 TRD/AliTRDsim.cxx create mode 100644 TRD/AliTRDsim.h create mode 100644 TRD/digitsAna.C create mode 100644 TRD/digitsCreate.C diff --git a/TRD/AliTRD.cxx b/TRD/AliTRD.cxx index d895d250603..20ad31d1096 100644 --- a/TRD/AliTRD.cxx +++ b/TRD/AliTRD.cxx @@ -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]; + +} diff --git a/TRD/AliTRD.h b/TRD/AliTRD.h index d299823867c..0bdc0d61764 100644 --- a/TRD/AliTRD.h +++ b/TRD/AliTRD.h @@ -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 diff --git a/TRD/AliTRDconst.h b/TRD/AliTRDconst.h index f7c9c42b68c..5cda8106b9b 100644 --- a/TRD/AliTRDconst.h +++ b/TRD/AliTRDconst.h @@ -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 index 00000000000..b429bade569 --- /dev/null +++ b/TRD/AliTRDmatrix.cxx @@ -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 index 00000000000..fdd4db06b07 --- /dev/null +++ b/TRD/AliTRDmatrix.h @@ -0,0 +1,59 @@ +#ifndef TRDmatrix_h +#define TRDmatrix_h + +#include +#include +#include +#include +#include +#include +#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 index 00000000000..7d121d9fdf4 --- /dev/null +++ b/TRD/AliTRDpixel.cxx @@ -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 index 00000000000..63e132f5645 --- /dev/null +++ b/TRD/AliTRDpixel.h @@ -0,0 +1,32 @@ +#ifndef TRDpixel_h +#define TRDpixel_h + +#include + +////////////////////////////////////////////////////// +// 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 index 00000000000..e81e6f2130e --- /dev/null +++ b/TRD/AliTRDsim.cxx @@ -0,0 +1,310 @@ +#include + +#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; + /* (binsize corr.) */ + om = h100->GetMean(); + /* */ + if (fIrst == 1) { + /* prints the production */ + printf(" Produced TR - , : %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[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 #include +#include #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 ! // diff --git a/TRD/AliTRDv2.h b/TRD/AliTRDv2.h index 22c1022afce..1798d1952f6 100644 --- a/TRD/AliTRDv2.h +++ b/TRD/AliTRDv2.h @@ -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) }; diff --git a/TRD/Makefile b/TRD/Makefile index 7d86e70f210..6877bbe24dd 100644 --- a/TRD/Makefile +++ b/TRD/Makefile @@ -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 diff --git a/TRD/TRDLinkDef.h b/TRD/TRDLinkDef.h index eee4f572dcc..c4c6974b93b 100644 --- a/TRD/TRDLinkDef.h +++ b/TRD/TRDLinkDef.h @@ -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 index 00000000000..9ba651fb1a7 --- /dev/null +++ b/TRD/digitsAna.C @@ -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 index 00000000000..590549a9fb8 --- /dev/null +++ b/TRD/digitsCreate.C @@ -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); + +} -- 2.43.0