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