Made Getters const
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 6 Oct 2000 16:50:31 +0000 (16:50 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 6 Oct 2000 16:50:31 +0000 (16:50 +0000)
16 files changed:
TRD/AliTRDcluster.cxx [new file with mode: 0644]
TRD/AliTRDcluster.h [new file with mode: 0644]
TRD/AliTRDfindTracks.C [new file with mode: 0644]
TRD/AliTRDpoints.cxx [new file with mode: 0644]
TRD/AliTRDpoints.h [new file with mode: 0644]
TRD/AliTRDsim.cxx [new file with mode: 0644]
TRD/AliTRDsim.h [new file with mode: 0644]
TRD/AliTRDtimeBin.cxx [new file with mode: 0644]
TRD/AliTRDtimeBin.h [new file with mode: 0644]
TRD/AliTRDtrack.cxx [new file with mode: 0644]
TRD/AliTRDtrack.h [new file with mode: 0644]
TRD/AliTRDtracker.cxx [new file with mode: 0644]
TRD/AliTRDtracker.h [new file with mode: 0644]
TRD/AliTRDtrackingSector.cxx [new file with mode: 0644]
TRD/AliTRDtrackingSector.h [new file with mode: 0644]
TRD/TrackDisplay.C [new file with mode: 0644]

diff --git a/TRD/AliTRDcluster.cxx b/TRD/AliTRDcluster.cxx
new file mode 100644 (file)
index 0000000..bd1954a
--- /dev/null
@@ -0,0 +1,64 @@
+/**************************************************************************
+ * 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$
+Revision 1.1.2.1  2000/09/22 14:47:52  cblume
+Add the tracking code
+
+*/
+
+#include "AliTRDcluster.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDrecPoint.h"
+
+
+ClassImp(AliTRDcluster)
+//_____________________________________________________________________________
+AliTRDcluster::AliTRDcluster() {
+  //default constructor
+
+  fDetector = fTimeBin = 0;
+  fTracks[0]=fTracks[1]=fTracks[2]=0; 
+  fY=fZ=fQ=fSigmaY2=fSigmaZ2=0.;
+}
+
+//_____________________________________________________________________________
+AliTRDcluster::AliTRDcluster(AliTRDrecPoint *rp)
+{
+  //
+  // constructor from AliTRDrecPoint
+  //
+
+  fDetector   = rp->GetDetector();
+  fTimeBin    = rp->GetLocalTimeBin();
+
+  fTracks[0]  = rp->GetTrackIndex(0);
+  fTracks[1]  = rp->GetTrackIndex(1);
+  fTracks[2]  = rp->GetTrackIndex(2);
+
+  fQ          = rp->GetEnergy();
+
+  fY          = rp->GetY();
+  fZ          = rp->GetZ();
+  fSigmaY2    = rp->GetSigmaY2();
+  fSigmaZ2    = rp->GetSigmaZ2();  
+
+  fSigmaY2    = 1;
+  fSigmaZ2    = 5;  
+
+}
+
diff --git a/TRD/AliTRDcluster.h b/TRD/AliTRDcluster.h
new file mode 100644 (file)
index 0000000..4756148
--- /dev/null
@@ -0,0 +1,50 @@
+#ifndef ALITRDCLUSTER_H
+#define ALITRDCLUSTER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <TObject.h>
+
+class AliTRDrecPoint;
+
+class AliTRDcluster : public TObject {
+
+ public:
+
+  AliTRDcluster();
+  AliTRDcluster(AliTRDrecPoint *rp);
+  
+  virtual Int_t   GetDetector() const           { return fDetector; };
+  virtual Int_t   GetLocalTimeBin() const       { return fTimeBin; }
+
+  virtual Float_t GetSigmaY2() const            { return fSigmaY2; }
+  virtual Float_t GetSigmaZ2() const            { return fSigmaZ2; }
+  virtual Float_t GetY() const                  { return fY; }
+  virtual Float_t GetZ() const                  { return fZ; }
+  virtual Float_t GetQ() const                  { return fQ; }
+
+  Int_t   IsUsed() const                        { return (fQ<0) ? 1 : 0; }
+  void    Use()                                 { fQ=-fQ; }
+  Int_t   GetTrackIndex(Int_t i) const          { return fTracks[i]; }
+
+
+ protected:
+
+  Int_t    fDetector;    // TRD detector number
+  Int_t    fTimeBin;     // Time bin number within the detector
+
+  Int_t    fTracks[3];   // labels of overlapped tracks
+  Float_t  fQ;           // amplitude 
+  Float_t  fY;           // local Rphi coordinate (cm) within tracking sector
+  Float_t  fZ;           // local Z coordinate (cm) within tracking sector
+  Float_t  fSigmaY2;     // Y variance (cm)
+  Float_t  fSigmaZ2;     // Z variance (cm)  
+
+
+  ClassDef(AliTRDcluster,1) // Reconstructed point for the TRD
+};
+
+#endif
diff --git a/TRD/AliTRDfindTracks.C b/TRD/AliTRDfindTracks.C
new file mode 100644 (file)
index 0000000..1153c02
--- /dev/null
@@ -0,0 +1,38 @@
+void AliTRDfindTracks() {
+
+//////////////////////////////////////////////////////////////////////////
+//
+// Reads TRD clusters from file, finds tracks, and stores them in the file 
+//
+//////////////////////////////////////////////////////////////////////////
+
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+    cout << "Loaded shared libraries" << endl;
+  }       
+
+  Char_t *alifile = "AliTRDclusters.root"; 
+
+  cerr<<"got this far"<<endl;
+
+  AliTRDtracker *Tracker =
+    new AliTRDtracker("TheOnlyTRDtrackerSoFar","UniqueVariationOfIt");
+
+
+  Tracker->GetEvent(alifile);
+
+  Int_t inner, outer, delta=60;
+  for(Int_t i=0; i<1; i++) {
+    outer=179-i; inner=outer-delta;
+    Tracker->MakeSeeds(inner,outer);
+  }
+
+  Tracker->FindTracks();
+
+  Tracker->WriteTracks();
+
+
+}
diff --git a/TRD/AliTRDpoints.cxx b/TRD/AliTRDpoints.cxx
new file mode 100644 (file)
index 0000000..459586b
--- /dev/null
@@ -0,0 +1,146 @@
+/**************************************************************************
+ * 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$
+Revision 1.1.2.1  2000/09/18 13:44:21  cblume
+New class AliTRDpoints to display the TR photon hits
+
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  This class contains the TRD points for the ALICE event display           //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDpoints.h"
+#include "AliRun.h"
+#include "AliDetector.h"
+#include "TPad.h"
+#include "TView.h"
+ClassImp(AliTRDpoints)
+
+//_____________________________________________________________________________
+AliTRDpoints::AliTRDpoints():AliPoints()
+{
+  //
+  // Default constructor
+  //
+
+  fNTRpoints    = 0;
+  fTRpolyMarker = 0;
+
+}
+
+//_____________________________________________________________________________
+AliTRDpoints::AliTRDpoints(Int_t nhitsE, Int_t nhitsT):AliPoints(nhitsE)
+{
+  //
+  // Standard constructor
+  //
+
+  fNTRpoints    = nhitsT;
+  fTRpolyMarker = 0;
+
+}
+        
+//_____________________________________________________________________________
+AliTRDpoints::AliTRDpoints(const AliTRDpoints &p)
+{
+  //
+  // Copy contructor
+  //
+  ((AliTRDpoints &) p).Copy(*this);
+
+}
+
+//_____________________________________________________________________________
+AliTRDpoints::~AliTRDpoints()
+{
+  //
+  // Default destructor
+  //
+
+  if (fTRpolyMarker) delete fTRpolyMarker;
+
+}
+
+//_____________________________________________________________________________
+AliTRDpoints &AliTRDpoints::operator=(const AliTRDpoints &p)
+{
+  //
+  // Assignment operator 
+  //
+
+  if (this != &p) ((AliTRDpoints &) p).Copy(*this);
+  return *this;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpoints::Copy(TObject &p)
+{
+  //
+  // Copy function
+  //
+
+  ((AliTRDpoints &) p).fNTRpoints = fNTRpoints; 
+  for (Int_t i = 0; i < 3*fNTRpoints; i++) {
+    ((AliTRDpoints &) p).fTRpoints[i] = fTRpoints[i];
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpoints::Draw(Option_t *option) 
+{
+  //
+  // Draws a TRD point
+  //
+
+  AliPoints::Draw(option);
+
+  //if (fTRpolyMarker) delete fTRpolyMarker;
+  if (fNTRpoints) {
+    fTRpolyMarker = new TPolyMarker3D(fNTRpoints,fTRpoints,29);
+    fTRpolyMarker->SetMarkerColor(2); 
+    fTRpolyMarker->SetMarkerSize(0.8);
+    fTRpolyMarker->Draw(option);
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpoints::SetTRpoints(Int_t n, Float_t *coor) 
+{
+  //
+  // Sets the number and the coordinates of the photon hits
+  //
+
+  if (kNTRpoints >= 3 * n) {
+    fNTRpoints = n;
+    for (Int_t i = 0; i < 3*n; i++) {
+      fTRpoints[i] = coor[i];
+    } 
+  }
+  else {
+    printf("AliTRDpoints::SetTRpoints -- ");
+    printf("Boundary error: %d/%d\n",3*n,kNTRpoints);
+  }
+
+}
diff --git a/TRD/AliTRDpoints.h b/TRD/AliTRDpoints.h
new file mode 100644 (file)
index 0000000..4fd36d8
--- /dev/null
@@ -0,0 +1,36 @@
+#ifndef ALITRDPOINTS_H
+#define ALITRDPOINTS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliPoints.h"
+
+const Int_t kNTRpoints = 75;
+
+class AliTRDpoints : public AliPoints {
+
+ public:
+
+  AliTRDpoints();
+  AliTRDpoints(const AliTRDpoints &p);
+  AliTRDpoints(Int_t nhitsE, Int_t nhitsT);
+  virtual ~AliTRDpoints();
+  AliTRDpoints &operator=(const AliTRDpoints &p);
+
+  virtual void           Copy(TObject &p);   
+  virtual void           Draw(Option_t *option);
+
+  virtual void           SetTRpoints(Int_t n, Float_t *coor);
+
+ protected:
+
+  Float_t          fTRpoints[kNTRpoints];       //  The hits from TR photons
+  Int_t            fNTRpoints;                  //  The number of TR photon hits
+  TPolyMarker3D   *fTRpolyMarker;               //! Polymarker to draw the photon hits
+  
+  ClassDef(AliTRDpoints,1)                      // Class to draw TRD hits 
+
+};
+#endif
diff --git a/TRD/AliTRDsim.cxx b/TRD/AliTRDsim.cxx
new file mode 100644 (file)
index 0000000..1696637
--- /dev/null
@@ -0,0 +1,781 @@
+/**************************************************************************
+ * 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$
+Revision 1.3.2.1  2000/09/18 13:45:30  cblume
+New class AliTRDsim that simulates TR photons
+
+Revision 1.2  1999/09/29 09:24:35  fca
+Introduction of the Copyright and cvs Log
+
+*/
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  TRD simulation - multimodule (regular rad.)                              //
+//  after: M. CASTELLANO et al., COMP. PHYS. COMM. 51 (1988) 431             //
+//                             + COMP. PHYS. COMM. 61 (1990) 395             //
+//                                                                           //
+//   17.07.1998 - A.Andronic                                                 //
+//   08.12.1998 - simplified version                                         //
+//   11.07.2000 - Adapted code to aliroot environment (C.Blume)              //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <stdlib.h>
+
+#include "TH1.h"
+#include "TRandom.h"
+#include "TMath.h"
+#include "TParticle.h"
+
+#include "AliTRDsim.h"
+#include "AliTRDconst.h"
+#include "AliModule.h"
+
+ClassImp(AliTRDsim)
+
+//_____________________________________________________________________________
+AliTRDsim::AliTRDsim():TObject()
+{
+  //
+  // AliTRDsim default constructor
+  // 
+
+  Init();
+
+}
+
+//_____________________________________________________________________________
+AliTRDsim::AliTRDsim(AliModule *mod, Int_t foil, Int_t gap)
+{
+  //
+  // AliTRDsim constructor. Takes the material properties of the radiator
+  // foils and the gas in the gaps from AliModule <mod>.
+  // The default number of foils is 100 with a thickness of 20 mu. The 
+  // thickness of the gaps is 500 mu.
+  //
+
+  Float_t aFoil, zFoil, rhoFoil;
+  Float_t aGap,  zGap,  rhoGap;
+  Float_t rad, abs;
+  Char_t  name[21];
+
+  Init();
+
+  mod->AliGetMaterial(foil,name,aFoil,zFoil,rhoFoil,rad,abs);
+  mod->AliGetMaterial(gap ,name,aGap ,zGap ,rhoGap ,rad,abs);
+
+  fFoilDens  = rhoFoil;
+  fFoilA     = aFoil;
+  fFoilZ     = zFoil;
+  fFoilOmega = Omega(fFoilDens,fFoilZ,fFoilA);
+
+  fGapDens   = rhoGap;
+  fGapA      = aGap;
+  fGapZ      = zGap;
+  fGapOmega  = Omega(fGapDens ,fGapZ ,fGapA );
+
+}
+
+//_____________________________________________________________________________
+AliTRDsim::AliTRDsim(const AliTRDsim &s)
+{
+  //
+  // AliTRDsim copy constructor
+  //
+
+  ((AliTRDsim &) s).Copy(*this);
+
+}
+
+//_____________________________________________________________________________
+AliTRDsim::~AliTRDsim() 
+{
+  //
+  // AliTRDsim destructor
+  //
+
+  if (fSpectrum) delete fSpectrum;
+  if (fSigma)    delete fSigma;
+
+}
+
+//_____________________________________________________________________________
+AliTRDsim &AliTRDsim::operator=(const AliTRDsim &s)
+{
+  //
+  // Assignment operator
+  //
+
+  if (this != &s) ((AliTRDsim &) s).Copy(*this);
+  return *this;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDsim::Copy(TObject &s)
+{
+  //
+  // Copy function
+  //
+
+  ((AliTRDsim &) s).fNFoils     = fNFoils;
+  ((AliTRDsim &) s).fFoilThick  = fFoilThick;
+  ((AliTRDsim &) s).fFoilDens   = fFoilDens;
+  ((AliTRDsim &) s).fFoilOmega  = fFoilOmega;
+  ((AliTRDsim &) s).fFoilZ      = fFoilZ;
+  ((AliTRDsim &) s).fFoilA      = fFoilA;
+  ((AliTRDsim &) s).fGapThick   = fGapThick;
+  ((AliTRDsim &) s).fGapDens    = fGapDens;
+  ((AliTRDsim &) s).fGapOmega   = fGapOmega;
+  ((AliTRDsim &) s).fGapZ       = fGapZ;
+  ((AliTRDsim &) s).fGapA       = fGapA;
+  ((AliTRDsim &) s).fTemp       = fTemp;
+  ((AliTRDsim &) s).fSpNBins    = fSpNBins;
+  ((AliTRDsim &) s).fSpRange    = fSpRange;
+  ((AliTRDsim &) s).fSpBinWidth = fSpBinWidth;
+  ((AliTRDsim &) s).fSpLower    = fSpLower;
+  ((AliTRDsim &) s).fSpUpper    = fSpUpper;
+
+  if (((AliTRDsim &) s).fSigma) delete ((AliTRDsim &) s).fSigma;
+  ((AliTRDsim &) s).fSigma = new Double_t[fSpNBins];
+  for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
+    ((AliTRDsim &) s).fSigma[iBin] = fSigma[iBin];
+  }  
+
+  fSpectrum->Copy(*((AliTRDsim &) s).fSpectrum);
+
+}
+
+//_____________________________________________________________________________
+void AliTRDsim::Init()
+{
+  //
+  // Initialization 
+  // The default radiator are 100 prolypropilene foils of 20 mu thickness
+  // with gaps of 500 mu filled with CO2.
+  //      
+  // 
+
+  fNFoils     = 100;
+
+  fFoilThick  = 0.0020;
+  fFoilDens   = 0.92;   
+  fFoilZ      = 5.28571;
+  fFoilA      = 10.4286;
+  fFoilOmega  = Omega(fFoilDens,fFoilZ,fFoilA);
+
+  fGapThick   = 0.0500;
+  fGapDens    = 0.001977;  
+  fGapZ       = 7.45455;
+  fGapA       = 14.9091;
+  fGapOmega   = Omega(fGapDens ,fGapZ ,fGapA );
+
+  fTemp       = 293.16;
+
+  fSpNBins    = 200;
+  fSpRange    = 100;
+  fSpBinWidth = fSpRange / fSpNBins;
+  fSpLower    = 1.0 - 0.5 * fSpBinWidth;
+  fSpUpper    = fSpLower + fSpRange;
+
+  if (fSpectrum) delete fSpectrum;
+  fSpectrum   = new TH1D("TRspectrum","TR spectrum",fSpNBins,fSpLower,fSpUpper);
+
+  // Set the sigma values 
+  SetSigma();
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDsim::CreatePhotons(Int_t pdg, Float_t p
+                             , Int_t &nPhoton, Float_t *ePhoton)
+{
+  //
+  // Create TRD photons for a charged particle of type <pdg> with the total 
+  // momentum <p>. 
+  // Number of produced TR photons:       <nPhoton>
+  // Energies of the produced TR photons: <ePhoton>
+  //
+
+  // PDG codes
+  const Int_t kPdgEle  =  11;
+  const Int_t kPdgMuon =  13;
+  const Int_t kPdgPion = 211;
+  const Int_t kPdgKaon = 321;
+
+  Float_t  mass        = 0;
+  switch (TMath::Abs(pdg)) {
+  case kPdgEle:
+    mass      =  5.11e-4;
+    break;
+  case kPdgMuon:
+    mass      =  0.10566;
+    break;
+  case kPdgPion:
+    mass      =  0.13957;
+    break;
+  case kPdgKaon:
+    mass      =  0.4937;
+    break;
+  default:
+    return 0;
+    break;
+  };
+
+  // Calculate gamma
+  Double_t gamma = TMath::Sqrt(p*p + mass*mass) / mass;
+
+  // Calculate the TR photons
+  return TrPhotons(gamma, nPhoton, ePhoton);
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDsim::TrPhotons(Double_t gamma, Int_t &nPhoton, Float_t *ePhoton)
+{
+  //
+  // Produces TR photons.
+  //
+
+  const Double_t kAlpha  = 0.0072973;
+  const Int_t    kSumMax = 10;
+
+  Double_t kappa = fGapThick / fFoilThick;
+
+  fSpectrum->Reset();
+
+  // The TR spectrum
+  Double_t stemp = 0;
+  for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
+
+    // keV -> eV
+    Double_t energyeV = (fSpBinWidth * iBin + 1.0) * 1e3;
+
+    Double_t csFoil   = fFoilOmega / energyeV;
+    Double_t csGap    = fGapOmega  / energyeV;
+
+    Double_t rho1     = energyeV * fFoilThick * 1e4 * 2.5 
+                                 * (1.0 / (gamma*gamma) + csFoil*csFoil);
+    Double_t rho2     = energyeV * fFoilThick * 1e4 * 2.5 
+                                 * (1.0 / (gamma*gamma) + csGap *csGap);
+
+    // Calculate the sum
+    Double_t sum = 0;
+    for (Int_t iSum = 0; iSum < kSumMax; iSum++) {
+      Double_t tetan = (TMath::Pi() * 2.0 * (iSum+1) - (rho1 + kappa * rho2)) 
+                     / (kappa + 1.0);
+      if (tetan < 0.0) tetan = 0.0;
+      Double_t aux   = 1.0 / (rho1 + tetan) - 1.0 / (rho2 + tetan);
+               sum  += tetan * (aux*aux) * (1.0 - TMath::Cos(rho1 + tetan));
+    }
+
+    // Absorbtion
+    Double_t conv      = 1.0 - TMath::Exp(-fNFoils * fSigma[iBin]);
+
+    // eV -> keV
+    Float_t  energykeV = energyeV * 0.001;
+
+    // dN / domega
+    Double_t wn        = kAlpha * 4.0 / (fSigma[iBin] * (kappa + 1.0)) 
+                                * conv * sum / energykeV;
+    fSpectrum->SetBinContent(iBin,wn);
+
+    stemp += wn;
+
+  }
+
+  // <nTR> (binsize corr.)
+  Float_t ntr = stemp * fSpBinWidth;
+  // Number of TR photons from Poisson distribution with mean <ntr>
+  nPhoton = gRandom->Poisson(ntr);
+  // Energy of the TR photons
+  for (Int_t iPhoton = 0; iPhoton < nPhoton; iPhoton++) {
+    ePhoton[iPhoton] = fSpectrum->GetRandom();
+  }
+
+  return 1;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDsim::SetSigma() 
+{
+  //
+  // Sets the absorbtion crosssection for the energies of the TR spectrum
+  //
+
+  if (fSigma) delete fSigma;
+  fSigma = new Double_t[fSpNBins];
+  for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
+    Double_t energykeV = iBin * fSpBinWidth + 1.0;
+    fSigma[iBin]       = Sigma(energykeV);
+    //printf("SetSigma(): iBin = %d fSigma %g\n",iBin,fSigma[iBin]);
+  }
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::Sigma(Double_t energykeV)
+{
+  //
+  // Calculates the absorbtion crosssection for a one-foil-one-gap-radiator
+  //
+
+  // Gas at 0 C
+  const Double_t kTemp0 = 273.16;
+
+  // keV -> MeV
+  Double_t energyMeV = energykeV * 0.001;
+  if (energyMeV >= 0.001) {
+    return(GetMuPo(energyMeV) * fFoilDens * fFoilThick + 
+           GetMuCO(energyMeV) * fGapDens  * fGapThick  * fTemp/kTemp0);
+  }
+  else {
+    return 1e6;
+  }
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuPo(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for polypropylene
+  //
+
+  const Int_t kN = 36;
+
+  Double_t mu[kN] = { 1.894E+03, 5.999E+02, 2.593E+02
+                    , 7.743E+01, 3.242E+01, 1.643E+01
+                    , 9.432E+00, 3.975E+00, 2.088E+00
+                    , 7.452E-01, 4.315E-01, 2.706E-01
+                    , 2.275E-01, 2.084E-01, 1.970E-01
+                    , 1.823E-01, 1.719E-01, 1.534E-01
+                    , 1.402E-01, 1.217E-01, 1.089E-01
+                    , 9.947E-02, 9.198E-02, 8.078E-02
+                    , 7.262E-02, 6.495E-02, 5.910E-02   
+                    , 5.064E-02, 4.045E-02, 3.444E-02
+                    , 3.045E-02, 2.760E-02, 2.383E-02
+                   , 2.145E-02, 1.819E-02, 1.658E-02 };
+
+  Double_t en[kN] = { 1.000E-03, 1.500E-03, 2.000E-03
+                    , 3.000E-03, 4.000E-03, 5.000E-03
+                    , 6.000E-03, 8.000E-03, 1.000E-02
+                    , 1.500E-02, 2.000E-02, 3.000E-02
+                    , 4.000E-02, 5.000E-02, 6.000E-02
+                    , 8.000E-02, 1.000E-01, 1.500E-01
+                    , 2.000E-01, 3.000E-01, 4.000E-01
+                    , 5.000E-01, 6.000E-01, 8.000E-01
+                    , 1.000E+00, 1.250E+00, 1.500E+00
+                    , 2.000E+00, 3.000E+00, 4.000E+00
+                    , 5.000E+00, 6.000E+00, 8.000E+00
+                   , 1.000E+01, 1.500E+01, 2.000E+01 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuCO(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for CO2
+  //
+
+  const Int_t kN = 36;
+
+  Double_t mu[kN] = { 0.39383E+04, 0.13166E+04, 0.58750E+03
+                    , 0.18240E+03, 0.77996E+02, 0.40024E+02
+                    , 0.23116E+02, 0.96997E+01, 0.49726E+01
+                    , 0.15543E+01, 0.74915E+00, 0.34442E+00
+                    , 0.24440E+00, 0.20589E+00, 0.18632E+00
+                    , 0.16578E+00, 0.15394E+00, 0.13558E+00
+                    , 0.12336E+00, 0.10678E+00, 0.95510E-01
+                    , 0.87165E-01, 0.80587E-01, 0.70769E-01
+                    , 0.63626E-01, 0.56894E-01, 0.51782E-01
+                    , 0.44499E-01, 0.35839E-01, 0.30825E-01
+                    , 0.27555E-01, 0.25269E-01, 0.22311E-01
+                   , 0.20516E-01, 0.18184E-01, 0.17152E-01 };
+
+  Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
+                    , 0.30000E-02, 0.40000E-02, 0.50000E-02
+                    , 0.60000E-02, 0.80000E-02, 0.10000E-01
+                    , 0.15000E-01, 0.20000E-01, 0.30000E-01
+                    , 0.40000E-01, 0.50000E-01, 0.60000E-01
+                    , 0.80000E-01, 0.10000E+00, 0.15000E+00
+                    , 0.20000E+00, 0.30000E+00, 0.40000E+00
+                    , 0.50000E+00, 0.60000E+00, 0.80000E+00
+                    , 0.10000E+01, 0.12500E+01, 0.15000E+01
+                    , 0.20000E+01, 0.30000E+01, 0.40000E+01
+                    , 0.50000E+01, 0.60000E+01, 0.80000E+01
+                   , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuXe(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for xenon
+  //
+
+  const Int_t kN = 48;
+
+  Double_t mu[kN] = { 9.413E+03, 8.151E+03, 7.035E+03
+                    , 7.338E+03, 4.085E+03, 2.088E+03
+                    , 7.780E+02, 3.787E+02, 2.408E+02
+                    , 6.941E+02, 6.392E+02, 6.044E+02
+                    , 8.181E+02, 7.579E+02, 6.991E+02
+                    , 8.064E+02, 6.376E+02, 3.032E+02
+                    , 1.690E+02, 5.743E+01, 2.652E+01
+                    , 8.930E+00, 6.129E+00, 3.316E+01
+                    , 2.270E+01, 1.272E+01, 7.825E+00
+                    , 3.633E+00, 2.011E+00, 7.202E-01
+                    , 3.760E-01, 1.797E-01, 1.223E-01
+                    , 9.699E-02, 8.281E-02, 6.696E-02
+                    , 5.785E-02, 5.054E-02, 4.594E-02
+                    , 4.078E-02, 3.681E-02, 3.577E-02
+                    , 3.583E-02, 3.634E-02, 3.797E-02
+                   , 3.987E-02, 4.445E-02, 4.815E-02 };
+
+  Double_t en[kN] = { 1.00000E-03, 1.07191E-03, 1.14900E-03
+                    , 1.14900E-03, 1.50000E-03, 2.00000E-03
+                    , 3.00000E-03, 4.00000E-03, 4.78220E-03
+                    , 4.78220E-03, 5.00000E-03, 5.10370E-03
+                    , 5.10370E-03, 5.27536E-03, 5.45280E-03
+                    , 5.45280E-03, 6.00000E-03, 8.00000E-03
+                    , 1.00000E-02, 1.50000E-02, 2.00000E-02
+                    , 3.00000E-02, 3.45614E-02, 3.45614E-02
+                    , 4.00000E-02, 5.00000E-02, 6.00000E-02
+                    , 8.00000E-02, 1.00000E-01, 1.50000E-01
+                    , 2.00000E-01, 3.00000E-01, 4.00000E-01
+                    , 5.00000E-01, 6.00000E-01, 8.00000E-01
+                    , 1.00000E+00, 1.25000E+00, 1.50000E+00
+                    , 2.00000E+00, 3.00000E+00, 4.00000E+00
+                    , 5.00000E+00, 6.00000E+00, 8.00000E+00
+                   , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuBu(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for isobutane
+  //
+
+  const Int_t kN = 36;
+
+  Double_t mu[kN] = { 0.38846E+03, 0.12291E+03, 0.53225E+02
+                    , 0.16091E+02, 0.69114E+01, 0.36541E+01
+                    , 0.22282E+01, 0.11149E+01, 0.72887E+00
+                    , 0.45053E+00, 0.38167E+00, 0.33920E+00
+                    , 0.32155E+00, 0.30949E+00, 0.29960E+00
+                    , 0.28317E+00, 0.26937E+00, 0.24228E+00
+                    , 0.22190E+00, 0.19289E+00, 0.17288E+00
+                    , 0.15789E+00, 0.14602E+00, 0.12829E+00
+                    , 0.11533E+00, 0.10310E+00, 0.93790E-01
+                    , 0.80117E-01, 0.63330E-01, 0.53229E-01
+                    , 0.46390E-01, 0.41425E-01, 0.34668E-01
+                   , 0.30267E-01, 0.23910E-01, 0.20509E-01 };
+
+  Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
+                    , 0.30000E-02, 0.40000E-02, 0.50000E-02
+                    , 0.60000E-02, 0.80000E-02, 0.10000E-01
+                    , 0.15000E-01, 0.20000E-01, 0.30000E-01
+                    , 0.40000E-01, 0.50000E-01, 0.60000E-01
+                    , 0.80000E-01, 0.10000E+00, 0.15000E+00
+                    , 0.20000E+00, 0.30000E+00, 0.40000E+00
+                    , 0.50000E+00, 0.60000E+00, 0.80000E+00
+                    , 0.10000E+01, 0.12500E+01, 0.15000E+01
+                    , 0.20000E+01, 0.30000E+01, 0.40000E+01
+                    , 0.50000E+01, 0.60000E+01, 0.80000E+01
+                   , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuMy(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for mylar
+  //
+
+  const Int_t kN = 36;
+
+  Double_t mu[kN] = { 2.911E+03, 9.536E+02, 4.206E+02
+                    , 1.288E+02, 5.466E+01, 2.792E+01
+                    , 1.608E+01, 6.750E+00, 3.481E+00
+                    , 1.132E+00, 5.798E-01, 3.009E-01
+                    , 2.304E-01, 2.020E-01, 1.868E-01
+                    , 1.695E-01, 1.586E-01, 1.406E-01
+                    , 1.282E-01, 1.111E-01, 9.947E-02
+                    , 9.079E-02, 8.395E-02, 7.372E-02
+                    , 6.628E-02, 5.927E-02, 5.395E-02
+                    , 4.630E-02, 3.715E-02, 3.181E-02
+                    , 2.829E-02, 2.582E-02, 2.257E-02
+                    , 2.057E-02, 1.789E-02, 1.664E-02 };
+
+  Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
+                    , 3.00000E-03, 4.00000E-03, 5.00000E-03
+                    , 6.00000E-03, 8.00000E-03, 1.00000E-02
+                    , 1.50000E-02, 2.00000E-02, 3.00000E-02
+                    , 4.00000E-02, 5.00000E-02, 6.00000E-02
+                    , 8.00000E-02, 1.00000E-01, 1.50000E-01
+                    , 2.00000E-01, 3.00000E-01, 4.00000E-01
+                    , 5.00000E-01, 6.00000E-01, 8.00000E-01
+                    , 1.00000E+00, 1.25000E+00, 1.50000E+00
+                    , 2.00000E+00, 3.00000E+00, 4.00000E+00
+                    , 5.00000E+00, 6.00000E+00, 8.00000E+00
+                    , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuN2(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for nitrogen
+  //
+
+  const Int_t kN = 36;
+
+  Double_t mu[kN] = { 3.311E+03, 1.083E+03, 4.769E+02
+                    , 1.456E+02, 6.166E+01, 3.144E+01
+                    , 1.809E+01, 7.562E+00, 3.879E+00
+                    , 1.236E+00, 6.178E-01, 3.066E-01
+                    , 2.288E-01, 1.980E-01, 1.817E-01
+                    , 1.639E-01, 1.529E-01, 1.353E-01
+                    , 1.233E-01, 1.068E-01, 9.557E-02
+                    , 8.719E-02, 8.063E-02, 7.081E-02
+                    , 6.364E-02, 5.693E-02, 5.180E-02
+                    , 4.450E-02, 3.579E-02, 3.073E-02
+                    , 2.742E-02, 2.511E-02, 2.209E-02
+                    , 2.024E-02, 1.782E-02, 1.673E-02 };
+
+  Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
+                    , 3.00000E-03, 4.00000E-03, 5.00000E-03
+                    , 6.00000E-03, 8.00000E-03, 1.00000E-02
+                    , 1.50000E-02, 2.00000E-02, 3.00000E-02
+                    , 4.00000E-02, 5.00000E-02, 6.00000E-02
+                    , 8.00000E-02, 1.00000E-01, 1.50000E-01
+                    , 2.00000E-01, 3.00000E-01, 4.00000E-01
+                    , 5.00000E-01, 6.00000E-01, 8.00000E-01
+                    , 1.00000E+00, 1.25000E+00, 1.50000E+00
+                    , 2.00000E+00, 3.00000E+00, 4.00000E+00
+                    , 5.00000E+00, 6.00000E+00, 8.00000E+00
+                    , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuO2(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for oxygen
+  //
+
+  const Int_t kN = 36;
+
+  Double_t mu[kN] = { 4.590E+03, 1.549E+03, 6.949E+02
+                    , 2.171E+02, 9.315E+01, 4.790E+01
+                    , 2.770E+01, 1.163E+01, 5.952E+00
+                    , 1.836E+00, 8.651E-01, 3.779E-01
+                    , 2.585E-01, 2.132E-01, 1.907E-01
+                    , 1.678E-01, 1.551E-01, 1.361E-01
+                    , 1.237E-01, 1.070E-01, 9.566E-02
+                    , 8.729E-02, 8.070E-02, 7.087E-02
+                    , 6.372E-02, 5.697E-02, 5.185E-02
+                    , 4.459E-02, 3.597E-02, 3.100E-02
+                    , 2.777E-02, 2.552E-02, 2.263E-02
+                    , 2.089E-02, 1.866E-02, 1.770E-02 };
+
+  Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
+                    , 3.00000E-03, 4.00000E-03, 5.00000E-03
+                    , 6.00000E-03, 8.00000E-03, 1.00000E-02
+                    , 1.50000E-02, 2.00000E-02, 3.00000E-02
+                    , 4.00000E-02, 5.00000E-02, 6.00000E-02
+                    , 8.00000E-02, 1.00000E-01, 1.50000E-01
+                    , 2.00000E-01, 3.00000E-01, 4.00000E-01
+                    , 5.00000E-01, 6.00000E-01, 8.00000E-01
+                    , 1.00000E+00, 1.25000E+00, 1.50000E+00
+                    , 2.00000E+00, 3.00000E+00, 4.00000E+00
+                    , 5.00000E+00, 6.00000E+00, 8.00000E+00
+                    , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuHe(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for helium
+  //
+
+  const Int_t kN = 36;
+
+  Double_t mu[kN] = { 6.084E+01, 1.676E+01, 6.863E+00
+                    , 2.007E+00, 9.329E-01, 5.766E-01
+                    , 4.195E-01, 2.933E-01, 2.476E-01
+                    , 2.092E-01, 1.960E-01, 1.838E-01
+                    , 1.763E-01, 1.703E-01, 1.651E-01
+                    , 1.562E-01, 1.486E-01, 1.336E-01
+                    , 1.224E-01, 1.064E-01, 9.535E-02
+                    , 8.707E-02, 8.054E-02, 7.076E-02
+                    , 6.362E-02, 5.688E-02, 5.173E-02
+                    , 4.422E-02, 3.503E-02, 2.949E-02
+                    , 2.577E-02, 2.307E-02, 1.940E-02
+                    , 1.703E-02, 1.363E-02, 1.183E-02 };
+
+  Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
+                    , 3.00000E-03, 4.00000E-03, 5.00000E-03
+                    , 6.00000E-03, 8.00000E-03, 1.00000E-02
+                    , 1.50000E-02, 2.00000E-02, 3.00000E-02
+                    , 4.00000E-02, 5.00000E-02, 6.00000E-02
+                    , 8.00000E-02, 1.00000E-01, 1.50000E-01
+                    , 2.00000E-01, 3.00000E-01, 4.00000E-01
+                    , 5.00000E-01, 6.00000E-01, 8.00000E-01
+                    , 1.00000E+00, 1.25000E+00, 1.50000E+00
+                    , 2.00000E+00, 3.00000E+00, 4.00000E+00
+                    , 5.00000E+00, 6.00000E+00, 8.00000E+00
+                    , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDsim::Interpolate(Double_t energyMeV
+                              , Double_t *en, Double_t *mu, Int_t n)
+{
+  //
+  // Interpolates the photon absorbtion cross section 
+  // for a given energy <energyMeV>.
+  //
+
+  Double_t de    = 0;
+  Int_t    index = 0;
+  Int_t    istat = Locate(en,n,energyMeV,index,de);
+  if (istat == 0) {
+    return (mu[index] - de * (mu[index]   - mu[index+1]) 
+                           / (en[index+1] - en[index]  ));
+  }
+  else {
+    return 0.0; 
+  }
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDsim::Locate(Double_t *xv, Int_t n, Double_t xval
+                      , Int_t &kl, Double_t &dx) 
+{
+  //
+  // Locates a point (xval) in a 1-dim grid (xv(n))
+  //
+
+  if (xval >= xv[n-1]) return  1;
+  if (xval <  xv[0])   return -1;
+
+  Int_t km;
+  Int_t kh = n - 1;
+
+  kl = 0;
+  while (kh - kl > 1) {
+    if (xval < xv[km = (kl+kh)/2]) kh = km; 
+    else                           kl = km;
+  }
+  if (xval < xv[kl] || xval > xv[kl+1] || kl >= n-1) {
+    printf("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n"
+          ,kl,xv[kl],xval,kl+1,xv[kl+1]);
+    exit(1);
+  }
+
+  dx = xval - xv[kl];
+
+  return 0;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDsim::Streamer(TBuffer &R__b)
+{
+  //
+  // Stream an object of class AliTRDsim.
+  //
+
+  if (R__b.IsReading()) {
+    Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+    TObject::Streamer(R__b);
+    R__b >> fNFoils;
+    R__b >> fFoilThick;
+    R__b >> fGapThick;
+    R__b >> fFoilDens;
+    R__b >> fGapDens;
+    R__b >> fFoilOmega;
+    R__b >> fGapOmega;
+    R__b >> fFoilZ;
+    R__b >> fGapZ;
+    R__b >> fFoilA;
+    R__b >> fGapA;
+    R__b >> fTemp;
+    R__b >> fSpNBins;
+    R__b >> fSpRange;
+    R__b >> fSpBinWidth;
+    R__b >> fSpLower;
+    R__b >> fSpUpper;
+    R__b.ReadArray(fSigma);
+    R__b >> fSpectrum;
+  } 
+  else {
+    R__b.WriteVersion(AliTRDsim::IsA());
+    TObject::Streamer(R__b);
+    R__b << fNFoils;
+    R__b << fFoilThick;
+    R__b << fGapThick;
+    R__b << fFoilDens;
+    R__b << fGapDens;
+    R__b << fFoilOmega;
+    R__b << fGapOmega;
+    R__b << fFoilZ;
+    R__b << fGapZ;
+    R__b << fFoilA;
+    R__b << fGapA;
+    R__b << fTemp;
+    R__b << fSpNBins;
+    R__b << fSpRange;
+    R__b << fSpBinWidth;
+    R__b << fSpLower;
+    R__b << fSpUpper;
+    R__b.WriteArray(fSigma, fSpNBins);
+    R__b << (TObject*) fSpectrum;
+  }
+
+}
diff --git a/TRD/AliTRDsim.h b/TRD/AliTRDsim.h
new file mode 100644 (file)
index 0000000..0d1f1ff
--- /dev/null
@@ -0,0 +1,112 @@
+#ifndef ALITRDSIM_H
+#define ALITRDSIM_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <TObject.h>
+
+class TH1D;
+
+class AliModule;
+
+class AliTRDsim : public TObject {
+
+ public:
+
+  AliTRDsim();
+  AliTRDsim(const AliTRDsim &s);
+  AliTRDsim(AliModule *mod, Int_t foil, Int_t gap);
+  virtual ~AliTRDsim();
+  AliTRDsim &operator=(const AliTRDsim &s);
+
+  virtual void          Copy(TObject &s);
+  virtual void          Init();
+  virtual Int_t         CreatePhotons(Int_t pdg, Float_t p
+                                    , Int_t &nPhoton, Float_t *ePhoton);
+  virtual Int_t         TrPhotons(Double_t gamma, Int_t &nPhoton, Float_t *ePhoton);
+  virtual Double_t      Sigma(Double_t energykeV);
+  virtual Double_t      Interpolate(Double_t energyMeV
+                                  , Double_t *en, Double_t *mu, Int_t n);
+  virtual Int_t         Locate(Double_t *xv, Int_t n, Double_t xval
+                             , Int_t &kl, Double_t &dx);
+  virtual Double_t      Omega(Float_t rho, Float_t z, Float_t a) 
+                             { return (28.8 * TMath::Sqrt(rho * z / a)); };
+
+          void          SetNFoils(Int_t n)      { fNFoils    = n;       };
+          void          SetFoilThick(Float_t t) { fFoilThick = t;
+                           SetSigma();                                  };
+          void          SetGapThick(Float_t t)  { fGapThick  = t;
+                           SetSigma();                                  };
+          void          SetFoilDens(Float_t d)  { fFoilDens  = d; 
+                           fFoilOmega = Omega(fFoilDens,fFoilZ,fFoilA);
+                           SetSigma();                                  };
+          void          SetFoilZ(Float_t z)     { fFoilZ     = z; 
+                           fFoilOmega = Omega(fFoilDens,fFoilZ,fFoilA); };
+          void          SetFoilA(Float_t a)     { fFoilA     = a; 
+                           fFoilOmega = Omega(fFoilDens,fFoilZ,fFoilA); };
+          void          SetGapDens(Float_t d)   { fGapDens   = d;
+                           fGapOmega  = Omega(fGapDens ,fGapZ ,fGapA );
+                           SetSigma();                                  };
+          void          SetGapZ(Float_t z)      { fGapZ      = z;
+                           fGapOmega  = Omega(fGapDens ,fGapZ ,fGapA ); };
+          void          SetGapA(Float_t a)      { fGapA      = a;
+                           fGapOmega  = Omega(fGapDens ,fGapZ ,fGapA ); };
+          void          SetTemp(Float_t t)      { fTemp      = t; 
+                           SetSigma();                                  };
+          void          SetSigma();
+
+  virtual Double_t      GetMuPo(Double_t energyMeV);
+  virtual Double_t      GetMuCO(Double_t energyMeV);
+  virtual Double_t      GetMuXe(Double_t energyMeV);
+  virtual Double_t      GetMuBu(Double_t energyMeV);
+  virtual Double_t      GetMuMy(Double_t energyMeV);
+  virtual Double_t      GetMuN2(Double_t energyMeV);
+  virtual Double_t      GetMuO2(Double_t energyMeV);
+  virtual Double_t      GetMuHe(Double_t energyMeV);
+
+          Int_t         GetNFoils() const       { return fNFoils;    };
+          Float_t       GetFoilThick() const    { return fFoilThick; };
+          Float_t       GetGapThick() const     { return fGapThick;  };
+          Float_t       GetFoilDens() const     { return fFoilDens;  };
+          Float_t       GetGapDens() const      { return fGapDens;   };
+          Double_t      GetFoilOmega() const    { return fFoilOmega; };
+          Double_t      GetGapOmega() const     { return fGapOmega;  };
+          Float_t       GetTemp() const         { return fTemp;      };
+          TH1D         *GetSpectrum() const     { return fSpectrum;  };
+
+ protected:
+
+  Int_t     fNFoils;               // Number of foils in the radiator stack
+  Float_t   fFoilThick;            // Thickness of the foils (cm)
+  Float_t   fGapThick;             // Thickness of the gaps between the foils (cm)
+
+  Float_t   fFoilDens;             // Density of the radiator foils (g/cm^3) 
+  Float_t   fGapDens;              // Density of the gas in the radiator gaps (g/cm^3)
+
+  Double_t  fFoilOmega;            // Plasma frequency of the radiator foils
+  Double_t  fGapOmega;             // Plasma frequency of the gas in the radiator gaps
+
+  Float_t   fFoilZ;                // Z of the foil material
+  Float_t   fGapZ;                 // Z of the gas in the gaps
+
+  Float_t   fFoilA;                // A of the foil material
+  Float_t   fGapA;                 // A of the gas in the gaps
+
+  Float_t   fTemp;                 // Temperature of the radiator gas (Kelvin)
+
+  Int_t     fSpNBins;              // Number of bins of the TR spectrum
+  Float_t   fSpRange;              // Range of the TR spectrum
+  Float_t   fSpBinWidth;           // Bin width of the TR spectrum
+  Float_t   fSpLower;              // Lower border of the TR spectrum
+  Float_t   fSpUpper;              // Upper border of the TR spectrum
+
+  Double_t *fSigma;                // Array of sigma values
+
+  TH1D     *fSpectrum;             // TR photon energy spectrum
+
+  ClassDef(AliTRDsim,1)            // Simulates TR photons
+
+};
+#endif
diff --git a/TRD/AliTRDtimeBin.cxx b/TRD/AliTRDtimeBin.cxx
new file mode 100644 (file)
index 0000000..c6dd571
--- /dev/null
@@ -0,0 +1,75 @@
+/**************************************************************************
+ * 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$
+Revision 1.1.2.2  2000/10/04 16:34:58  cblume
+Replace include files by forward declarations
+
+Revision 1.1.2.1  2000/09/22 14:47:52  cblume
+Add the tracking code
+
+*/                        
+                                
+#include "AliTRDcluster.h" 
+#include "AliTRDconst.h"
+#include "AliTRDtimeBin.h" 
+
+ClassImp(AliTRDtimeBin)
+
+//______________________________________________________
+
+void AliTRDtimeBin::InsertCluster(AliTRDcluster* c, UInt_t index) {
+
+// Insert cluster in TimeBin cluster array.
+// Clusters are sorted according to Y coordinate.  
+
+  if (fN==kMAX_CLUSTER_PER_TIME_BIN) {
+    printf("AliTRDtimeBin::InsertCluster(): Too many clusters !\n"); 
+    return;
+  }
+  if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
+  Int_t i=Find(c->GetY());
+  memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
+  memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t)); 
+  fIndex[i]=index; fClusters[i]=c; fN++;
+}  
+
+//______________________________________________________
+
+Int_t AliTRDtimeBin::Find(Double_t y) const {
+
+// Returns index of the cluster nearest in Y    
+
+  if (y <= fClusters[0]->GetY()) return 0;
+  if (y > fClusters[fN-1]->GetY()) return fN;
+  Int_t b=0, e=fN-1, m=(b+e)/2;
+  for (; b<e; m=(b+e)/2) {
+    if (y > fClusters[m]->GetY()) b=m+1;
+    else e=m;
+  }
+  return m;
+}    
+
+//______________________________________________________
+AliTRDcluster *AliTRDtimeBin::operator[](Int_t i)
+{
+  //
+  // Index operator
+  //
+  return fClusters[i];
+
+}
diff --git a/TRD/AliTRDtimeBin.h b/TRD/AliTRDtimeBin.h
new file mode 100644 (file)
index 0000000..a02c3d8
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALITRDTIMEBIN_H
+#define ALITRDTIMEBIN_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDtimeBin.h,v */
+
+#include <TObject.h>
+
+class AliTRDcluster;
+
+const unsigned kMAX_CLUSTER_PER_TIME_BIN=3500; 
+
+//----------------------------------------------------------------- 
+class AliTRDtimeBin : public TObject {
+
+// Provides tools to address clusters which lay within one time bin
+public: 
+
+  AliTRDtimeBin() {fN=0;}
+  void InsertCluster(AliTRDcluster*,UInt_t);
+  operator Int_t() const {return fN;}
+  AliTRDcluster* operator[](Int_t i);
+  UInt_t GetIndex(Int_t i) const {return fIndex[i];} 
+
+  Int_t Find(Double_t y) const; 
+
+ protected:
+   unsigned fN;
+   AliTRDcluster *fClusters[kMAX_CLUSTER_PER_TIME_BIN];
+   UInt_t fIndex[kMAX_CLUSTER_PER_TIME_BIN]; 
+
+  ClassDef(AliTRDtimeBin,1) // Provides tools to address clusters which lay within one time bin
+
+}; 
+
+#endif 
+
diff --git a/TRD/AliTRDtrack.cxx b/TRD/AliTRDtrack.cxx
new file mode 100644 (file)
index 0000000..0d0b0af
--- /dev/null
@@ -0,0 +1,451 @@
+/**************************************************************************
+ * 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$
+Revision 1.1.2.1  2000/09/22 14:47:52  cblume
+Add the tracking code
+
+*/                                                        
+
+#include <iostream.h>
+
+#include <TObject.h>
+
+#include "AliTRD.h" 
+#include "AliTRDconst.h"
+#include "AliTRDgeometry.h" 
+#include "AliTRDcluster.h" 
+#include "AliTRDtrack.h"
+
+ClassImp(AliTRDtrack)
+
+
+//_____________________________________________________________________________
+
+AliTRDtrack::AliTRDtrack(UInt_t index, const Double_t xx[5],
+const Double_t cc[15], Double_t xref, Double_t alpha) {
+  //-----------------------------------------------------------------
+  // This is the main track constructor.
+  //-----------------------------------------------------------------
+  fLab=-1;
+  fChi2=0.;
+  fdEdx=0.;
+
+  fAlpha=alpha;
+  fX=xref;
+
+  fY=xx[0]; fZ=xx[1]; fC=xx[2]; fE=xx[3]; fT=xx[4];
+
+  fCyy=cc[0];
+  fCzy=cc[1];  fCzz=cc[2];
+  fCcy=cc[3];  fCcz=cc[4];  fCcc=cc[5];
+  fCey=cc[6];  fCez=cc[7];  fCec=cc[8];  fCee=cc[9];
+  fCty=cc[10]; fCtz=cc[11]; fCtc=cc[12]; fCte=cc[13]; fCtt=cc[14];
+
+  fN=0;
+  fIndex[fN++]=index;
+}                              
+           
+//_____________________________________________________________________________
+AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) {
+  //
+  // Copy constructor.
+  //
+
+  fLab=t.fLab;
+
+  fChi2=t.fChi2;
+  fdEdx=t.fdEdx;
+
+  fAlpha=t.fAlpha;
+  fX=t.fX;
+
+  fY=t.fY; fZ=t.fZ; fC=t.fC; fE=t.fE; fT=t.fT;
+
+  fCyy=t.fCyy;
+  fCzy=t.fCzy;  fCzz=t.fCzz;
+  fCcy=t.fCcy;  fCcz=t.fCcz;  fCcc=t.fCcc;
+  fCey=t.fCey;  fCez=t.fCez;  fCec=t.fCec;  fCee=t.fCee;
+  fCty=t.fCty;  fCtz=t.fCtz;  fCtc=t.fCtc;  fCte=t.fCte;  fCtt=t.fCtt;
+
+  fN=t.fN;
+  for (Int_t i=0; i<fN; i++) fIndex[i]=t.fIndex[i];
+}                                                       
+
+//_____________________________________________________________________________
+void AliTRDtrack::GetCovariance(Double_t cc[15]) const {
+  cc[0]=fCyy;
+  cc[1]=fCzy;  cc[2]=fCzz;
+  cc[3]=fCcy;  cc[4]=fCcz;  cc[5]=fCcc;
+  cc[6]=fCey;  cc[7]=fCez;  cc[8]=fCec;  cc[9]=fCee;
+  cc[10]=fCty; cc[11]=fCtz; cc[12]=fCtc; cc[13]=fCte; cc[14]=fCtt;
+}    
+
+//_____________________________________________________________________________
+Int_t AliTRDtrack::Compare(TObject *o) {
+
+// Compares tracks according to their Y2
+
+  AliTRDtrack *t=(AliTRDtrack*)o;
+  //  Double_t co=t->GetSigmaY2();
+  //  Double_t c =GetSigmaY2();
+
+  Double_t co=TMath::Abs(t->GetC());
+  Double_t c =TMath::Abs(GetC());  
+
+  if (c>co) return 1;
+  else if (c<co) return -1;
+  return 0;
+}                
+
+//_____________________________________________________________________________
+Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
+{
+  // Propagates a track of particle with mass=pm to a reference plane 
+  // defined by x=xk through media of density=rho and radiationLength=x0
+
+  if (TMath::Abs(fC*xk - fE) >= 0.99999) {
+    if (fN>4) cerr<<fN<<" AliTRDtrack warning: Propagation failed !\n";
+    return 0;
+  }
+
+  Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=fY, z1=fZ;
+  Double_t c1=fC*x1 - fE, r1=sqrt(1.- c1*c1);
+  Double_t c2=fC*x2 - fE, r2=sqrt(1.- c2*c2);
+
+  fY += dx*(c1+c2)/(r1+r2);
+  fZ += dx*(c1+c2)/(c1*r2 + c2*r1)*fT;
+
+  //f = F - 1
+  Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
+  Double_t f02= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
+  Double_t f03=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
+  Double_t cr=c1*r2+c2*r1;
+  Double_t f12= dx*fT*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
+  Double_t f13=-dx*fT*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
+  Double_t f14= dx*cc/cr;
+
+  //b = C*ft
+  Double_t b00=f02*fCcy + f03*fCey, b01=f12*fCcy + f13*fCey + f14*fCty;
+  Double_t b10=f02*fCcz + f03*fCez, b11=f12*fCcz + f13*fCez + f14*fCtz;
+  Double_t b20=f02*fCcc + f03*fCec, b21=f12*fCcc + f13*fCec + f14*fCtc;
+  Double_t b30=f02*fCec + f03*fCee, b31=f12*fCec + f13*fCee + f14*fCte;
+  Double_t b40=f02*fCtc + f03*fCte, b41=f12*fCtc + f13*fCte + f14*fCtt;
+
+  //a = f*b = f*C*ft
+  Double_t a00=f02*b20+f03*b30,a01=f02*b21+f03*b31,a11=f12*b21+f13*b31+f14*b41;
+
+  //F*C*Ft = C + (a + b + bt)
+  fCyy += a00 + 2*b00;
+  fCzy += a01 + b01 + b10;
+  fCcy += b20;
+  fCey += b30;
+  fCty += b40;
+  fCzz += a11 + 2*b11;
+  fCcz += b21;
+  fCez += b31;
+  fCtz += b41;                  
+
+  fX=x2;
+
+
+  //Multiple scattering  ******************
+
+  Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-fY)*(y1-fY)+(z1-fZ)*(z1-fZ));
+  Double_t p2=GetPt()*GetPt()*(1.+fT*fT);
+  Double_t beta2=p2/(p2 + pm*pm);
+
+  Double_t ey=fC*fX - fE, ez=fT;
+  Double_t xz=fC*ez, zz1=ez*ez+1, xy=fE+ey;
+
+  Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
+  fCcc += xz*xz*theta2;
+  fCec += xz*ez*xy*theta2;
+  fCtc += xz*zz1*theta2;
+  fCee += (2*ey*ez*ez*fE+1-ey*ey+ez*ez+fE*fE*ez*ez)*theta2;
+  fCte += ez*zz1*xy*theta2;
+  fCtt += zz1*zz1*theta2;
+
+
+  //Energy losses************************
+
+  Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
+  if (x1 < x2) dE=-dE;
+  fC*=(1.- sqrt(p2+pm*pm)/p2*dE);
+  //fE*=(1.- sqrt(p2+pm*pm)/p2*dE);
+
+  return 1;        
+
+}     
+
+
+//_____________________________________________________________________________
+void AliTRDtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
+{
+  // This function propagates tracks to the "vertex".
+
+  Double_t c=fC*fX - fE;
+  Double_t tgf=-fE/(fC*fY + sqrt(1-c*c));
+  Double_t snf=tgf/sqrt(1.+ tgf*tgf);
+  Double_t xv=(fE+snf)/fC;
+  PropagateTo(xv,x0,rho,pm); 
+}          
+
+
+//_____________________________________________________________________________
+void AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index)
+{
+  // Assignes found cluster to the track and updates track information
+
+  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2()*12;
+  r00+=fCyy; r01+=fCzy; r11+=fCzz;
+  Double_t det=r00*r11 - r01*r01;
+  Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
+
+  Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
+  Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
+  Double_t k20=fCcy*r00+fCcz*r01, k21=fCcy*r01+fCcz*r11;
+  Double_t k30=fCey*r00+fCez*r01, k31=fCey*r01+fCez*r11;
+  Double_t k40=fCty*r00+fCtz*r01, k41=fCty*r01+fCtz*r11;
+
+  Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
+  Double_t cur=fC + k20*dy + k21*dz, eta=fE + k30*dy + k31*dz;
+  if (TMath::Abs(cur*fX-eta) >= 0.99999) {
+    if (fN>4) cerr<<fN<<" AliTRDtrack warning: Filtering failed !\n";
+    return;
+  }
+
+  fY += k00*dy + k01*dz;
+  fZ += k10*dy + k11*dz;
+  fC  = cur;
+  fE  = eta;
+  fT += k40*dy + k41*dz;
+
+  Double_t c01=fCzy, c02=fCcy, c03=fCey, c04=fCty;
+  Double_t c12=fCcz, c13=fCez, c14=fCtz;
+
+  fCyy-=k00*fCyy+k01*fCzy; fCzy-=k00*c01+k01*fCzz;
+  fCcy-=k00*c02+k01*c12; fCey-=k00*c03+k01*c13;
+  fCty-=k00*c04+k01*c14;
+
+  fCzz-=k10*c01+k11*fCzz;
+  fCcz-=k10*c02+k11*c12; fCez-=k10*c03+k11*c13;
+  fCtz-=k10*c04+k11*c14;
+
+  fCcc-=k20*c02+k21*c12; fCec-=k20*c03+k21*c13;
+  fCtc-=k20*c04+k21*c14;
+
+  fCee-=k30*c03+k31*c13;
+  fCte-=k30*c04+k31*c14;        
+
+  fCtt-=k40*c04+k41*c14;
+
+  fIndex[fN++]=index;
+  fChi2 += chisq;   
+
+  //  cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
+}                     
+
+//_____________________________________________________________________________
+Int_t AliTRDtrack::Rotate(Double_t alpha)
+{
+  // Rotates track parameters in R*phi plane
+
+  fAlpha += alpha;
+
+  Double_t x1=fX, y1=fY;
+  Double_t ca=cos(alpha), sa=sin(alpha);
+  Double_t r1=fC*fX - fE;
+
+  fX = x1*ca + y1*sa;
+  fY=-x1*sa + y1*ca;
+  fE=fE*ca + (fC*y1 + sqrt(1.- r1*r1))*sa;
+
+  Double_t r2=fC*fX - fE;
+  if (TMath::Abs(r2) >= 0.99999) {
+    if (fN>4) cerr<<fN<<" AliTRDtrack warning: Rotation failed !\n";
+    return 0;
+  }
+
+  Double_t y0=fY + sqrt(1.- r2*r2)/fC;
+  if ((fY-y0)*fC >= 0.) {
+    if (fN>4) cerr<<fN<<" AliTRDtrack warning: Rotation failed !!!\n";
+    return 0;
+  }
+
+  //f = F - 1
+  Double_t f00=ca-1,    f32=(y1 - r1*x1/sqrt(1.- r1*r1))*sa,
+           f30=fC*sa, f33=(ca + sa*r1/sqrt(1.- r1*r1))-1;
+
+  //b = C*ft
+  Double_t b00=fCyy*f00, b03=fCyy*f30+fCcy*f32+fCey*f33;
+  Double_t b10=fCzy*f00, b13=fCzy*f30+fCcz*f32+fCez*f33;
+  Double_t b20=fCcy*f00, b23=fCcy*f30+fCcc*f32+fCec*f33;
+  Double_t b30=fCey*f00, b33=fCey*f30+fCec*f32+fCee*f33;
+  Double_t b40=fCty*f00, b43=fCty*f30+fCtc*f32+fCte*f33;
+
+  //a = f*b = f*C*ft
+  Double_t a00=f00*b00, a03=f00*b03, a33=f30*b03+f32*b23+f33*b33;
+
+  // *** Double_t dy2=fCyy;  
+          
+  //F*C*Ft = C + (a + b + bt)
+  fCyy += a00 + 2*b00;
+  fCzy += b10;
+  fCcy += b20;
+  fCey += a03+b30+b03;
+  fCty += b40;
+  fCez += b13;
+  fCec += b23;
+  fCee += a33 + 2*b33;
+  fCte += b43;
+
+  // *** fCyy+=dy2*sa*sa*r1*r1/(1.- r1*r1);
+  // *** fCzz+=d2y*sa*sa*fT*fT/(1.- r1*r1);   
+
+  return 1;
+}                         
+
+
+
+
+//_____________________________________________________________________________
+Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c) const
+{
+  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2()*12;
+  r00+=fCyy; r01+=fCzy; r11+=fCzz;
+
+  Double_t det=r00*r11 - r01*r01;
+  if (TMath::Abs(det) < 1.e-10) {
+    if (fN>4) cerr<<fN<<" AliTRDtrack warning: Singular matrix !\n";
+    return 1e10;
+  }
+  Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
+
+  Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
+
+  return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;  
+}            
+
+
+//_________________________________________________________________________
+void AliTRDtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
+{
+  // Returns reconstructed track momentum in the global system.
+
+  Double_t pt=TMath::Abs(GetPt()); // GeV/c
+  Double_t r=fC*fX-fE;
+  Double_t y0=fY + sqrt(1.- r*r)/fC;
+  px=-pt*(fY-y0)*fC;    //cos(phi);
+  py=-pt*(fE-fX*fC);   //sin(phi);
+  pz=pt*fT;
+  Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha);
+  py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha);
+  px=tmp;            
+
+}                                
+
+//____________________________________________________________________________
+void AliTRDtrack::Streamer(TBuffer &R__b)
+{
+   if (R__b.IsReading()) {
+      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+      TObject::Streamer(R__b);
+      R__b >> fLab;
+      R__b >> fChi2;
+      R__b >> fdEdx;
+      R__b >> fAlpha;
+      R__b >> fX;
+      R__b >> fY;
+      R__b >> fZ;
+      R__b >> fC;
+      R__b >> fE;
+      R__b >> fT;
+      R__b >> fCyy;
+      R__b >> fCzy;
+      R__b >> fCzz;
+      R__b >> fCcy;
+      R__b >> fCcz;
+      R__b >> fCcc;
+      R__b >> fCey;
+      R__b >> fCez;
+      R__b >> fCec;
+      R__b >> fCee;
+      R__b >> fCty;
+      R__b >> fCtz;
+      R__b >> fCtc;
+      R__b >> fCte;
+      R__b >> fCtt;
+      R__b >> fN;
+      for (Int_t i=0; i<fN; i++) R__b >> fIndex[i];
+   } else {                                
+      R__b.WriteVersion(AliTRDtrack::IsA());
+      TObject::Streamer(R__b);
+      R__b << fLab;
+      R__b << fChi2;
+      R__b << fdEdx;
+      R__b << fAlpha;
+      R__b << fX;
+      R__b << fY;
+      R__b << fZ;
+      R__b << fC;
+      R__b << fE;
+      R__b << fT;
+      R__b << fCyy;
+      R__b << fCzy;
+      R__b << fCzz;
+      R__b << fCcy;
+      R__b << fCcz;
+      R__b << fCcc;
+      R__b << fCey;
+      R__b << fCez;
+      R__b << fCec;
+      R__b << fCee;
+      R__b << fCty;
+      R__b << fCtz;
+      R__b << fCtc;
+      R__b << fCte;
+      R__b << fCtt;
+      R__b << fN;
+      for (Int_t i=0; i<fN; i++) R__b << fIndex[i];
+   }
+}                                                          
+
+//_____________________________________________________________________________
+void AliTRDseed::CookdEdx(Double_t low, Double_t up) {
+
+  // Calculates dE/dX within the "low" and "up" cuts.
+
+  Int_t i;
+  Int_t nc=this->GetNclusters();
+
+  Int_t swap;//stupid sorting
+  do {
+    swap=0;
+    for (i=0; i<nc-1; i++) {
+      if (fdEdx[i]<=fdEdx[i+1]) continue;
+      Float_t tmp=fdEdx[i]; fdEdx[i]=fdEdx[i+1]; fdEdx[i+1]=tmp;
+      swap++;
+    }
+  } while (swap);
+
+  Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
+  Float_t dedx=0;
+  for (i=nl; i<=nu; i++) dedx += fdEdx[i];
+  dedx /= (nu-nl+1);
+  SetdEdx(dedx);
+}
+
diff --git a/TRD/AliTRDtrack.h b/TRD/AliTRDtrack.h
new file mode 100644 (file)
index 0000000..834863e
--- /dev/null
@@ -0,0 +1,104 @@
+#ifndef ALITRDTRACK_H
+#define ALITRDTRACK_H  
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */ 
+
+#include <TObject.h> 
+
+class AliTRDcluster;
+
+class AliTRDtrack : public TObject {
+
+// Represents reconstructed TRD track
+
+   Int_t    fLab;         // track label  
+   Double_t fChi2;        // total chi2 value for the track  
+   Float_t  fdEdx;        // dE/dx 
+
+   Double_t fAlpha;       // rotation angle
+   Double_t fX;           // running local X-coordinate of the track (time bin)
+
+   Double_t fY;           // Y-coordinate of the track
+   Double_t fZ;           // Z-coordinate of the track
+   Double_t fC;           // track curvature
+   Double_t fE;           // C*x0
+   Double_t fT;           // tangent of the track dip angle   
+
+   Double_t fCyy;                         // covariance
+   Double_t fCzy, fCzz;                   // matrix
+   Double_t fCcy, fCcz, fCcc;             // of the
+   Double_t fCey, fCez, fCec, fCee;       // track
+   Double_t fCty, fCtz, fCtc, fCte, fCtt; // parameters   
+
+   Short_t fN;             // number of clusters associated with the track
+   UInt_t  fIndex[200];    // global indexes of these clusters  
+                          
+
+public:
+
+   AliTRDtrack() { fN=0;}
+   AliTRDtrack(UInt_t index, const Double_t xx[5],
+               const Double_t cc[15], Double_t xr, Double_t alpha);  
+   AliTRDtrack(const AliTRDtrack& t);    
+
+   Int_t  Compare(TObject *o);
+
+   Double_t GetAlpha() const {return fAlpha;}
+   Double_t GetC()     const {return fC;}
+   Int_t    GetClusterIndex(Int_t i) const {return fIndex[i];}    
+   Double_t GetChi2()  const {return fChi2;}
+   void     GetCovariance(Double_t cov[15]) const;  
+   Double_t GetdEdX()  const {return fdEdx;}
+   Double_t GetEta()   const {return fE;}
+   Int_t    GetLabel() const {return fLab;}
+   Int_t    GetNclusters() const {return fN;}
+   Double_t GetP()     const {  
+     return TMath::Abs(GetPt())*sqrt(1.+GetTgl()*GetTgl());
+   }
+   Double_t GetPredictedChi2(const AliTRDcluster*) const ;
+   Double_t GetPt()    const {return 0.3*0.2/GetC()/100;}
+   void     GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const ;
+   Double_t GetSigmaC2()   const {return fCcc;}
+   Double_t GetSigmaTgl2() const {return fCtt;}
+   Double_t GetSigmaY2()   const {return fCyy;}
+   Double_t GetSigmaZ2()   const {return fCzz;}
+   Double_t GetTgl() const {return fT;}
+   Double_t GetX()   const {return fX;}
+   Double_t GetY()   const {return fY;} // returns running Y
+   Double_t GetZ()   const {return fZ;}
+
+   Bool_t   IsSortable() const {return kTRUE;}
+
+   Int_t    PropagateTo(Double_t xr,
+                   Double_t x0=8.72,Double_t rho=5.86e-3,Double_t pm=0.139);
+   void     PropagateToVertex(
+                   Double_t x0=36.66,Double_t rho=1.2e-3,Double_t pm=0.139);
+   Int_t    Rotate(Double_t angle);
+
+   void     SetLabel(Int_t lab=0) {fLab=lab;}
+   void     SetdEdx(Float_t dedx) {fdEdx=dedx;}  
+
+   void     Update(const AliTRDcluster* c, Double_t chi2, UInt_t i);
+
+   ClassDef(AliTRDtrack,1)  // TRD reconstructed tracks
+};                     
+
+
+//-----------------------------------------------------------------
+class AliTRDseed : public AliTRDtrack {
+   Float_t fdEdx[200];
+public:
+   AliTRDseed():AliTRDtrack(){}
+   AliTRDseed(UInt_t index, const Double_t xx[5],
+              const Double_t cc[15], Double_t xr, Double_t alpha):
+              AliTRDtrack(index, xx, cc, xr, alpha) {}
+   void SetSampledEdx(Float_t q, Int_t i) {
+     Double_t c=GetC(), e=GetEta(), t=GetTgl(), x=GetX();
+     q *= TMath::Sqrt((1-(c*x-e)*(c*x-e))/(1+t*t));
+     fdEdx[i]=q;
+   }
+   void CookdEdx(Double_t low=0.05, Double_t up=0.70);
+};            
+
+#endif   
diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx
new file mode 100644 (file)
index 0000000..a5bd09d
--- /dev/null
@@ -0,0 +1,680 @@
+/**************************************************************************
+ * 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$
+Revision 1.1.2.2  2000/10/04 16:34:58  cblume
+Replace include files by forward declarations
+
+Revision 1.1.2.1  2000/09/22 14:47:52  cblume
+Add the tracking code
+
+*/   
+
+#include <TFile.h>
+#include <TROOT.h>
+#include <TBranch.h>
+#include <TTree.h>
+
+#include "AliRun.h"
+
+#include "AliTRD.h"
+#include "AliTRDconst.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDrecPoint.h" 
+#include "AliTRDcluster.h" 
+#include "AliTRDtrack.h"
+#include "AliTRDtrackingSector.h"
+#include "AliTRDtimeBin.h"
+
+#include "AliTRDtracker.h"
+
+ClassImp(AliTRDtracker) 
+
+//____________________________________________________________________
+AliTRDtracker::AliTRDtracker()
+{
+  //
+  // Default constructor
+  //   
+
+  fInputFile = NULL;
+  fEvent     = 0;
+
+  fGeom      = NULL;
+  fNclusters = 0;
+  fClusters  = NULL; 
+  fNseeds    = 0;
+  fSeeds     = NULL;
+  fNtracks   = 0;
+  fTracks    = NULL;
+
+}   
+
+//____________________________________________________________________
+AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
+                  :TNamed(name, title)
+{
+  fInputFile = NULL;
+  fEvent     = 0;
+
+  fGeom      = NULL;
+  fNclusters = 0;
+  fClusters  = new TObjArray(2000); 
+  fNseeds    = 0;
+  fSeeds     = new TObjArray(20000);
+  fNtracks   = 0;
+  fTracks    = new TObjArray(10000);
+
+}   
+
+
+//___________________________________________________________________
+AliTRDtracker::~AliTRDtracker()
+{
+  if (fInputFile) {
+    fInputFile->Close();
+    delete fInputFile;
+  }
+  delete fClusters;
+  delete fTracks;
+  delete fSeeds;
+  delete fGeom;
+}   
+
+
+//_____________________________________________________________________
+static Double_t SigmaY2trd(Double_t r, Double_t tgl, Double_t pt)
+{
+  // Parametrised "expected" error of the cluster reconstruction in Y 
+
+  Double_t s = 1.;    
+  return s;
+}
+
+//_____________________________________________________________________
+static Double_t SigmaZ2trd(Double_t r, Double_t tgl)
+{
+  // Parametrised "expected" error of the cluster reconstruction in Z 
+
+  AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
+  AliTRDgeometry   *fGeom;   
+  fGeom = TRD->GetGeometry();
+  Double_t s, pad = fGeom->GetRowPadSize();
+  s = pad * pad /12.;  
+  return s;
+}                  
+
+//_____________________________________________________________________
+inline Double_t f1trd(Double_t x1,Double_t y1,
+                      Double_t x2,Double_t y2,
+                      Double_t x3,Double_t y3)
+{
+  // Initial approximation of the track curvature
+  // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+
+  Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
+  Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
+                  (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
+  Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
+                  (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
+
+  Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
+
+  return -xr*yr/sqrt(xr*xr+yr*yr);
+}          
+
+//_____________________________________________________________________
+inline Double_t f2trd(Double_t x1,Double_t y1,
+                      Double_t x2,Double_t y2,
+                      Double_t x3,Double_t y3)
+{
+  // Initial approximation of the track curvature times center of curvature
+  // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+
+  Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
+  Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
+                  (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
+  Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
+                  (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
+
+  Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
+
+  return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
+}          
+
+//_____________________________________________________________________
+inline Double_t f3trd(Double_t x1,Double_t y1,
+                      Double_t x2,Double_t y2,
+                      Double_t z1,Double_t z2)
+{
+  // Initial approximation of the tangent of the track dip angle
+  // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+
+  return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
+}            
+
+
+//___________________________________________________________________
+
+static Int_t FindProlongation(AliTRDseed& t, AliTRDtrackingSector *sec,
+                            Int_t s, Int_t rf=0)
+{
+  // Starting from current position on track=t this function tries 
+  // to extrapolate the track up to timeBin=rf and to confirm prolongation
+  // if a close cluster is found. *sec is a pointer to allocated
+  // array of sectors, in which the initial sector has index=s. 
+
+  const Int_t TIME_BINS_TO_SKIP=Int_t(0.2*sec->GetNtimeBins());
+  const Double_t MAX_CHI2=12.;
+  Int_t try_again=TIME_BINS_TO_SKIP;
+
+  Double_t alpha=sec->GetAlpha();
+
+  Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
+
+  for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
+    Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr);
+    if (!t.PropagateTo(x)) {
+      cerr<<"Can't propagate to x = "<<x<<endl;
+      return 0;
+    }
+
+    AliTRDcluster *cl=0;
+    UInt_t index=0;
+
+    Double_t max_chi2=MAX_CHI2;
+    //    const AliTRDtimeBin& time_bin=sec[s][nr];
+    AliTRDtimeBin& time_bin=sec[s][nr];
+    Double_t sy2=SigmaY2trd(t.GetX(),t.GetTgl(),t.GetPt());
+    Double_t sz2=SigmaZ2trd(t.GetX(),t.GetTgl());
+    Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
+
+    if (road>30) {
+      if (t.GetNclusters() > 4) {
+       cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
+      }
+      return 0;
+    }       
+
+    if (time_bin) {
+      for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
+        AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
+        if (c->GetY() > y+road) break;
+        if (c->IsUsed() > 0) continue;
+
+       if((c->GetZ()-z)*(c->GetZ()-z) > 25. + sz2) continue;
+
+        Double_t chi2=t.GetPredictedChi2(c);
+
+        if (chi2 > max_chi2) continue;
+        max_chi2=chi2;
+        cl=c;
+        index=time_bin.GetIndex(i);
+      }   
+    }
+    if (cl) {
+
+      //      Float_t l=sec->GetPitch();
+      //      t.SetSampledEdx(cl->fQ/l,Int_t(t));  
+     
+      t.Update(cl,max_chi2,index);
+
+      try_again=TIME_BINS_TO_SKIP;
+    } else {
+      if (try_again==0) break;
+      if (y > ymax) {
+       cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
+         s = (s+1) % ns;
+         if (!t.Rotate(alpha)) {
+          cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
+          return 0;
+        }
+      } else if (y <-ymax) {
+       cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
+         s = (s-1+ns) % ns;
+         if (!t.Rotate(-alpha)) {
+          cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
+          return 0;
+        }
+      }
+      try_again--;
+    }
+  }
+
+  return 1;
+}          
+
+
+
+//_____________________________________________________________________________
+void AliTRDtracker::GetEvent(const Char_t *name, Int_t nEvent)
+{
+  // Opens a ROOT-file with TRD-clusters and reads in the cluster-tree
+
+  // Connect the AliRoot file containing Geometry, Kine, and Hits
+  fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
+  if (!fInputFile) {
+    printf("AliTRDtracker::Open -- ");
+    printf("Open the ALIROOT-file %s.\n",name);
+    fInputFile = new TFile(name,"UPDATE");
+  }
+  else {
+    printf("AliTRDtracker::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("AliTRDtracker::GetEvent -- ");
+      printf("AliRun object found on file.\n");
+    }
+    else {
+      printf("AliTRDtracker::GetEvent -- ");
+      printf("Could not find AliRun object.\n");
+    }
+  //}        
+
+  fEvent = nEvent;
+
+  // Import the Trees for the event nEvent in the file
+  Int_t nparticles = gAlice->GetEvent(fEvent);
+  cerr<<"nparticles = "<<nparticles<<endl;
+
+  if (nparticles <= 0) {
+    printf("AliTRDtracker::GetEvent -- ");
+    printf("No entries in the trees for event %d.\n",fEvent);
+  }
+
+  AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
+  fGeom = TRD->GetGeometry();
+
+  Char_t treeName[14];
+  sprintf(treeName,"TRDrecPoints%d", nEvent);
+
+  TTree *tree=(TTree*)fInputFile->Get(treeName);
+  TBranch *branch=tree->GetBranch("TRDrecPoints");
+  Int_t nentr = (Int_t) tree->GetEntries();
+  printf("found %d entries in %s tree.\n",nentr,treeName);
+
+  for (Int_t i=0; i<nentr; i++) {
+    TObjArray *ioArray = new TObjArray(400);
+    branch->SetAddress(&ioArray);
+    tree->GetEvent(i);
+    Int_t npoints = ioArray->GetEntriesFast();
+    printf("Read %d rec. points from entry %d \n", npoints, i);
+    for(Int_t j=0; j<npoints; j++) {
+      AliTRDrecPoint *p=(AliTRDrecPoint*)ioArray->UncheckedAt(j);
+      AliTRDcluster  *c = new AliTRDcluster(p); 
+      fClusters->AddLast(c);
+    }
+  }
+
+}     
+
+
+//_____________________________________________________________________________
+void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
+{
+  // Fills clusters into TRD tracking_sectors 
+  // Note that the numbering scheme for the TRD tracking_sectors 
+  // differs from that of TRD sectors
+
+  for (Int_t i=0; i<kNsect; i++) sec[i].SetUp();
+
+  //  Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's 
+
+  cerr<<"MakeSeeds: sorting clusters"<<endl;
+              
+  Int_t ncl=fClusters->GetEntriesFast();
+  UInt_t index;
+  while (ncl--) {
+    printf("\r %d left",ncl); 
+    AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
+    Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
+    Int_t sector=fGeom->GetSector(detector);
+
+    Int_t tracking_sector=sector;
+    if(sector > 0) tracking_sector = kNsect - sector;
+
+    Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin); 
+    index=ncl;
+    sec[tracking_sector][tb].InsertCluster(c,index);
+  }    
+  printf("\r\n");
+}
+
+
+//_____________________________________________________________________________
+void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
+{
+  // Creates track seeds using clusters in timeBins=i1,i2
+
+  Int_t i2 = inner, i1=outer; 
+
+  if (!fClusters) return; 
+
+  AliTRDtrackingSector fTrSec[kNsect];
+  SetUpSectors(fTrSec);
+
+  // find seeds
+
+  Double_t x[5], c[15];
+  Int_t max_sec=kNsect;
+
+  Double_t alpha=fTrSec[0].GetAlpha();
+  Double_t shift=fTrSec[0].GetAlpha()/2.;
+  Double_t cs=cos(alpha), sn=sin(alpha);  
+
+  Double_t x1 =fTrSec[0].GetX(i1);
+  Double_t xx2=fTrSec[0].GetX(i2);  
+
+  for (Int_t ns=0; ns<max_sec; ns++) {
+
+    Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2];
+    Int_t nm=fTrSec[ns][i2];
+    Int_t nu=fTrSec[(ns+1)%max_sec][i2];
+
+    AliTRDtimeBin& r1=fTrSec[ns][i1];
+
+    for (Int_t is=0; is < r1; is++) {
+      Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
+
+      for (Int_t js=0; js < nl+nm+nu; js++) {
+       const AliTRDcluster *cl;
+        Double_t x2,   y2,   z2;
+        Double_t x3=0.,y3=0.;  
+
+        if (js<nl) {
+         AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2];
+         cl=r2[js]; 
+         y2=cl->GetY(); z2=cl->GetZ();
+
+          x2= xx2*cs+y2*sn;
+          y2=-xx2*sn+y2*cs;          
+
+        } else
+          if (js<nl+nm) {
+           AliTRDtimeBin& r2=fTrSec[ns][i2];
+           cl=r2[js-nl];
+           x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
+          } else {
+           AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2];
+            cl=r2[js-nl-nm];
+           y2=cl->GetY(); z2=cl->GetZ();
+
+            x2=xx2*cs-y2*sn;
+            y2=xx2*sn+y2*cs;   
+
+          }
+
+        Double_t zz=z1 - z1/x1*(x1-x2);
+       
+        if (TMath::Abs(zz-z2)>30.) continue;   
+
+        Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
+        if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
+
+        x[0]=y1;
+        x[1]=z1;
+        x[2]=f1trd(x1,y1,x2,y2,x3,y3);
+
+        if (TMath::Abs(x[2]) >= 0.01) continue;
+
+        x[3]=f2trd(x1,y1,x2,y2,x3,y3);
+
+        if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
+
+        x[4]=f3trd(x1,y1,x2,y2,z1,z2);
+
+        if (TMath::Abs(x[4]) > 1.2) continue;
+
+        Double_t a=asin(x[3]);
+        Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
+        if (TMath::Abs(zv)>200.) continue;    
+
+        Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2()*12;
+        Double_t sy2=cl->GetSigmaY2(),     sz2=cl->GetSigmaZ2()*12;
+        Double_t sy3=100*0.025, sy=0.1, sz=0.1;
+
+        Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
+        Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
+        Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
+        Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
+        Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
+        Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
+        Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
+        Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
+        Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
+        Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
+
+        c[0]=sy1;
+        c[1]=0.;       c[2]=sz1;
+        c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
+        c[6]=f30*sy1;  c[7]=0.;       c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
+                                      c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
+        c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
+        c[13]=f40*sy1*f30+f42*sy2*f32;
+        c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;   
+
+        UInt_t index=r1.GetIndex(is);
+        AliTRDseed *track=new AliTRDseed(index, x, c, x1, ns*alpha+shift); 
+
+       //        Float_t l=fTrSec->GetPitch();
+       //        track->SetSampledEdx(r1[is]->fQ/l,0); 
+
+        Int_t rc=FindProlongation(*track,fTrSec,ns,i2);
+       
+        if ((rc < 0) || (track->GetNclusters() < (i1-i2)/2) ) delete track;
+        else { 
+         fSeeds->AddLast(track); fNseeds++; 
+         cerr<<"found seed "<<fNseeds<<endl;
+       }
+      }
+    }
+  }
+
+  //  delete[] fTrSec;
+
+  fSeeds->Sort();
+}          
+
+//___________________________________________________________________
+void AliTRDtracker::FindTracks() 
+{
+  if (!fClusters) return; 
+
+  AliTRDtrackingSector fTrSec[kNsect];
+  SetUpSectors(fTrSec);
+
+  // find tracks
+
+  Int_t num_of_time_bins = fTrSec[0].GetNtimeBins(); 
+  Int_t nseed=fSeeds->GetEntriesFast();
+
+  Int_t nSeedClusters;
+  for (Int_t i=0; i<nseed; i++) {
+    cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
+
+    AliTRDseed& t=*((AliTRDseed*)fSeeds->UncheckedAt(i));
+
+    nSeedClusters = t.GetNclusters();
+    Double_t alpha=t.GetAlpha();
+
+    if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
+    if (alpha < 0.            ) alpha += 2.*TMath::Pi();
+    Int_t ns=Int_t(alpha/fTrSec[0].GetAlpha())%kNsect;
+
+    if (FindProlongation(t,fTrSec,ns)) {
+      cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl; 
+      if ((t.GetNclusters() >= Int_t(0.3*num_of_time_bins)) && 
+         ((t.GetNclusters()-nSeedClusters)>60)) {
+       t.CookdEdx();
+       Int_t label = GetTrackLabel(t);
+       t.SetLabel(label);
+       AliTRDtrack* pt=&t;
+        fTracks->AddLast(pt); fNtracks++;
+       UseClusters(t);
+       cerr<<"found track "<<fNtracks<<endl;
+      }                         
+    }     
+  }            
+}
+
+//__________________________________________________________________
+void AliTRDtracker::UseClusters(AliTRDseed t) {
+  Int_t ncl=t.GetNclusters();
+  for (Int_t i=0; i<ncl; i++) {
+    Int_t index = t.GetClusterIndex(i);
+    AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
+    c->Use();
+  }
+}
+
+//__________________________________________________________________
+Int_t AliTRDtracker::GetTrackLabel(AliTRDseed t) {
+
+  Int_t label=123456789, index, i, j;
+  Int_t ncl=t.GetNclusters();
+  const Int_t range=300;
+  Bool_t label_added;
+
+  Int_t s[range][2];
+  for (i=0; i<range; i++) {
+    s[i][0]=-1;
+    s[i][1]=0;
+  }
+
+  Int_t t0,t1,t2;
+  for (i=0; i<ncl; i++) {
+    index=t.GetClusterIndex(i);
+    AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
+    t0=c->GetTrackIndex(0);
+    t1=c->GetTrackIndex(1);
+    t2=c->GetTrackIndex(2);
+  }
+
+  for (i=0; i<ncl; i++) {
+    index=t.GetClusterIndex(i);
+    AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
+    for (Int_t k=0; k<3; k++) { 
+      label=c->GetTrackIndex(k);
+      label_added=false; j=0;
+      if (label >= 0) {
+       while ( (!label_added) && ( j < range ) ) {
+         if (s[j][0]==label || s[j][1]==0) {
+           s[j][0]=label; 
+           s[j][1]=s[j][1]+1; 
+           label_added=true;
+         }
+         j++;
+       }
+      }
+    }
+  }
+
+
+  Int_t max=0;
+  label = -123456789;
+
+  for (i=0; i<range; i++) {
+    if (s[i][1]>max) {
+      max=s[i][1]; label=s[i][0];
+    }
+  }
+  if(max > ncl/2) return label;
+  else return -1;
+}
+
+//___________________________________________________________________
+
+Int_t AliTRDtracker::WriteTracks() {
+
+  TDirectory *savedir=gDirectory;   
+
+  TFile *out=TFile::Open("AliTRDtracks.root","RECREATE");
+
+  TTree tracktree("TreeT","Tree with TRD tracks");
+
+  AliTRDtrack *iotrack=0;
+  tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);  
+
+  Int_t ntracks=fTracks->GetEntriesFast();
+
+  for (Int_t i=0; i<ntracks; i++) {
+    AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
+    iotrack=pt;
+    tracktree.Fill(); 
+    cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
+  }
+
+  tracktree.Write(); 
+  out->Close(); 
+
+  savedir->cd();  
+
+  cerr<<"WriteTracks: done"<<endl;
+  return 0;
+}
+
+
+
+//_____________________________________________________________________________
+void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename, Int_t nEvent = 0, Int_t option = 1)
+{
+  //
+  // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) 
+  // from the file. The names of the cluster tree and branches 
+  // should match the ones used in AliTRDclusterizer::WriteClusters()
+  //
+
+  TDirectory *savedir=gDirectory; 
+
+  TFile *file = TFile::Open(filename);
+  if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;} 
+
+  Char_t treeName[14];
+  sprintf(treeName,"TRDrecPoints%d", nEvent);
+
+  TTree *tree=(TTree*)file->Get(treeName);
+  TBranch *branch=tree->GetBranch("TRDrecPoints");
+  Int_t nentr = (Int_t) tree->GetEntries();
+  printf("found %d entries in %s tree.\n",nentr,treeName);
+
+  for (Int_t i=0; i<nentr; i++) {
+    TObjArray *ioArray = new TObjArray(400);
+    branch->SetAddress(&ioArray);
+    tree->GetEvent(i);
+    Int_t npoints = ioArray->GetEntriesFast();
+    printf("Read %d rec. points from entry %d \n", npoints, i);
+    for(Int_t j=0; j<npoints; j++) {
+      AliTRDrecPoint *p=(AliTRDrecPoint*)ioArray->UncheckedAt(j);
+      AliTRDcluster  *c = new AliTRDcluster(p); 
+      if( option >= 0) array->AddLast(c);
+      else array->AddLast(p);
+    }
+  }
+
+  file->Close();                   
+  savedir->cd(); 
+
+}
+
+
+
+
diff --git a/TRD/AliTRDtracker.h b/TRD/AliTRDtracker.h
new file mode 100644 (file)
index 0000000..2d726ac
--- /dev/null
@@ -0,0 +1,54 @@
+#ifndef ALITRDTRACKER_H
+#define ALITRDTRACKER_H   
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */ 
+
+#include <TNamed.h>
+
+class TFile;
+class TObjArray;
+
+class AliTRDgeometry;
+// class AliTRDtrackingSector;
+class AliTRDtrack;
+class AliTRDseed;
+
+
+class AliTRDtracker : public TNamed { 
+
+ public:
+
+  AliTRDtracker();
+  AliTRDtracker(const Text_t* name, const Text_t* title);
+  ~AliTRDtracker(); 
+
+  virtual void  GetEvent(const Char_t *name, Int_t nEvent = 0);
+  virtual void  SetUpSectors(AliTRDtrackingSector *sec);
+  virtual void  MakeSeeds(Int_t inner, Int_t outer);
+  virtual void  FindTracks();
+  virtual void  UseClusters(AliTRDseed t);
+  virtual Int_t GetTrackLabel(AliTRDseed t);
+  virtual Int_t WriteTracks(); 
+  virtual void  ReadClusters(TObjArray *array, const Char_t *filename, Int_t nEvent = 0, Int_t option = 1);
+
+ protected:
+
+  TFile            *fInputFile;       // AliROOT input file
+  AliTRDgeometry   *fGeom;            // Pointer to TRD geometry
+  Int_t            fEvent;            // Event number
+
+  Int_t            fNclusters;        // Number of clusters in TRD 
+  TObjArray        *fClusters;        // List of clusters for all sectors
+
+  Int_t            fNseeds;           // Number of track seeds  
+  TObjArray        *fSeeds;           // List of track seeds
+   
+  Int_t            fNtracks;          // Number of reconstructed tracks 
+  TObjArray        *fTracks;          // List of reconstructed tracks   
+
+  ClassDef(AliTRDtracker,1)           // manager base class  
+
+};
+
+#endif 
diff --git a/TRD/AliTRDtrackingSector.cxx b/TRD/AliTRDtrackingSector.cxx
new file mode 100644 (file)
index 0000000..ae58d2c
--- /dev/null
@@ -0,0 +1,157 @@
+/**************************************************************************
+ * 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$
+Revision 1.1.2.2  2000/10/04 16:34:58  cblume
+Replace include files by forward declarations
+
+Revision 1.1.2.1  2000/09/22 14:47:52  cblume
+Add the tracking code
+
+*/                        
+                                
+#include <TObject.h>
+
+#include "AliRun.h"
+
+#include "AliTRD.h" 
+#include "AliTRDconst.h"
+#include "AliTRDgeometry.h" 
+#include "AliTRDcluster.h" 
+#include "AliTRDtimeBin.h" 
+
+#include "AliTRDtrackingSector.h" 
+
+
+ClassImp(AliTRDtrackingSector) 
+
+//_______________________________________________________
+
+AliTRDtrackingSector::~AliTRDtrackingSector()
+{
+  //
+  // Destructor
+  //
+
+  delete[] fTimeBin;
+
+}
+
+//_______________________________________________________
+AliTRDtimeBin &AliTRDtrackingSector::operator[](Int_t i)
+{
+  //
+  // Index operator 
+  //
+
+  return *(fTimeBin+i);
+
+}
+
+//_______________________________________________________
+
+void AliTRDtrackingSector::SetUp()
+{ 
+
+  AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
+  fGeom = TRD->GetGeometry();
+
+  fTimeBinSize = fGeom->GetTimeBinSize();
+
+  fN = kNplan * (Int_t(kDrThick/fTimeBinSize) + 1);
+
+  fTimeBin = new AliTRDtimeBin[fN]; 
+
+}
+
+//______________________________________________________
+
+Double_t AliTRDtrackingSector::GetX(Int_t l) const
+{
+  if( (l<0) || (l>fN-1)) {
+    fprintf(stderr,"AliTRDtrackingSector::GetX: TimeBin index is out of range !\n");
+    return -99999.;
+  }
+  else { 
+    
+    Int_t tb_per_plane = fN/kNplan;
+    Int_t plane = l/tb_per_plane;
+    Int_t time_slice = l%(Int_t(kDrThick/fTimeBinSize) + 1);
+
+    Float_t t0 = fGeom->GetTime0(plane);
+    Double_t x = t0 + time_slice * fTimeBinSize;
+
+    //    cerr<<"plane, tb, x = "<<plane<<","<<time_slice<<","<<x<<endl;
+
+    return x;
+  }
+}
+
+//______________________________________________________
+
+Double_t AliTRDtrackingSector::GetMaxY(Int_t l) const 
+{ 
+
+  if((l<(fN-1)) && (l>-1)) { 
+    Int_t tb_per_plane = fN/kNplan;
+    Int_t plane = l/tb_per_plane;
+    return fGeom->GetChamberWidth(plane); 
+  }
+  else {
+    fprintf(stderr,
+    "AliTRDtrackingSector::GetMaxY: TimeBin index is out of range !\n");
+    if(l<0) return fGeom->GetChamberWidth(0);
+    else return fGeom->GetChamberWidth(kNplan-1);
+  }
+}
+
+//______________________________________________________
+
+Int_t AliTRDtrackingSector::GetTimeBinNumber(Double_t x) const
+{
+  Float_t r_out = fGeom->GetTime0(kNplan-1) + kDrThick; 
+  Float_t r_in = fGeom->GetTime0(0);
+  //  cerr<<"GetTimeBinNumber: r_in,r_out = "<<r_in<<","<<r_out<<endl;
+
+  if(x >= r_out) return fN-1;
+  if(x <= r_in) return 0;
+
+  Float_t gap = fGeom->GetTime0(1) - fGeom->GetTime0(0);
+  //  cerr<<"GetTimeBinNumber: gap = "<<gap<<endl;
+
+  Int_t plane = Int_t((x - r_in + fTimeBinSize/2)/gap);
+  //  cerr<<"GetTimeBinNumber: plane="<<plane<<endl;
+
+  Int_t local_tb = Int_t((x-fGeom->GetTime0(plane))/fTimeBinSize + 0.5);
+  //  cerr<<"GetTimeBinNumber: local_tb="<<local_tb<<endl;
+
+  Int_t time_bin = plane * (Int_t(kDrThick/fTimeBinSize) + 1) + local_tb;
+  //  cerr<<"GetTimeBinNumber: time_bin = "<<time_bin<<endl;
+  
+
+  return time_bin;
+}
+
+//______________________________________________________
+
+Int_t AliTRDtrackingSector::GetTimeBin(Int_t det, Int_t local_tb) const 
+{
+  Int_t plane = fGeom->GetPlane(det);
+
+  Int_t time_bin = plane * (Int_t(kDrThick/fTimeBinSize) + 1) + local_tb;  
+  return time_bin;
+}
+
diff --git a/TRD/AliTRDtrackingSector.h b/TRD/AliTRDtrackingSector.h
new file mode 100644 (file)
index 0000000..1317ec6
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALITRDTRACKINGSECTOR_H
+#define ALITRDTRACKINGSECTOR_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDtrackingSector.h,v */
+
+#include <TObject.h>
+
+class AliTRDtimeBin;
+class AliTRDgeometry;
+
+class AliTRDtrackingSector : public TObject {
+
+// Provides tools to address clusters which lay within one sector
+
+protected:
+  Int_t fN;
+  AliTRDtimeBin           *fTimeBin;    // Pointer to array of AliTRDtimeBin
+  AliTRDgeometry          *fGeom;       // Pointer to TRD geometry
+  Float_t                  fTimeBinSize;  // Time bin size in cm  
+
+public:
+  AliTRDtrackingSector() {fN=0; fTimeBin=0; fGeom=0; fTimeBinSize=0;}
+  virtual ~AliTRDtrackingSector();
+  virtual void SetUp();
+  AliTRDtimeBin& operator[](Int_t i);
+  Int_t GetNtimeBins() const { return fN; }
+  Double_t GetX(Int_t l) const;
+  Double_t GetMaxY(Int_t l) const; 
+  Double_t GetAlpha() const { return 2*TMath::Pi()/kNsect; } 
+  Int_t GetTimeBinNumber(Double_t x) const;
+  Int_t GetTimeBin(Int_t det, Int_t local_tb) const;
+  Float_t GetPitch() const {return fTimeBinSize;}   
+
+  ClassDef(AliTRDtrackingSector,1)  // Provides tools to address clusters which lay within one sector
+
+}; 
+
+#endif 
diff --git a/TRD/TrackDisplay.C b/TRD/TrackDisplay.C
new file mode 100644 (file)
index 0000000..83a031a
--- /dev/null
@@ -0,0 +1,141 @@
+void TrackDisplay(Int_t track) {
+
+  c1 = new TCanvas( "RecPoints", "RecPoints Display", 10, 10, 710, 740);
+  c1->SetFillColor(1);
+  TView *v=new TView(1);
+
+      v->SetRange(-350.,-350.,-400.,350.,350.,400.); // full
+  //  v->SetRange(0.,0.,0.,350.,350.,400.);  // top right
+  //  v->SetRange(-350.,0.,0.,0.,350.,400.); // top left
+  //  v->SetRange(0.,-350.,0.,350.,0.,400.); // bottom right 
+  //  v->SetRange(-350.,-350.,0.,0.,0.,400.); // bottom left
+  //  v->Side();
+  v->Top();
+  cerr<<"psi = "<<v->GetPsi()<<endl;
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+    cout << "Loaded shared libraries" << endl;
+  }       
+
+
+// Load tracks
+
+   TFile *tf=TFile::Open("AliTRDtracks.root");
+   if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return 3;}
+   TObjArray tarray(2000);
+   TTree *tracktree=(TTree*)tf->Get("TreeT");
+   TBranch *tbranch=tracktree->GetBranch("tracks");
+   Int_t nentr=tracktree->GetEntries();
+   cerr<<"Found "<<nentr<<" entries in the track tree"<<endl;
+   for (Int_t i=0; i<nentr; i++) {
+       AliTRDtrack *iotrack=new AliTRDtrack;
+       tbranch->SetAddress(&iotrack);
+       tracktree->GetEvent(i);
+       tarray.AddLast(iotrack);
+       cerr<<"got track "<<i<<": index ="<<iotrack->GetLabel()<<endl;
+   }
+   tf->Close();                 
+
+
+  // Load clusters
+  Char_t *alifile = "AliTRDclusters.root";
+  Int_t   nEvent  = 0;
+  TObjArray rparray(2000);
+  TObjArray *RecPointsArray = &rparray;
+  AliTRDtracker *Tracker = new AliTRDtracker("dummy","dummy");
+  Tracker->ReadClusters(RecPointsArray,alifile,nEvent,-1); 
+  Int_t nRecPoints = RecPointsArray->GetEntriesFast();
+  cerr<<"Found "<<nRecPoints<<" rec. points"<<endl;
+
+
+  // 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");
+
+
+  AliTRDv1       *TRD           = (AliTRDv1*) gAlice->GetDetector("TRD");
+  AliTRDgeometry *fGeom   = TRD->GetGeometry();   
+
+  Int_t i,j,index,det,sector, ti[3];
+  Double_t x,y,z, cs,sn,tmp;
+  Float_t global[3], local[3];
+
+  // Display all TRD RecPoints
+  TPolyMarker3D *pm = new TPolyMarker3D(nRecPoints);
+
+  for (Int_t i = 0; i < nRecPoints; i++) {
+    printf("\r point %d out of %d",i,nRecPoints);
+    AliTRDrecPoint *rp = (AliTRDrecPoint *) RecPointsArray->UncheckedAt(i);    
+
+    local[0]=rp->GetLocalRow(); 
+    local[1]=rp->GetLocalCol(); 
+    local[2]=rp->GetLocalTime();
+    det=rp->GetDetector();
+
+    for (j = 0; j < 3; j++) { ti[j] = rp->GetTrackIndex(j); }
+
+    if((track < 0) ||
+       ((ti[0]==track)||(ti[1]==track)||(ti[2]==track))) {
+
+      if(fGeom->Local2Global(det,local,global)) {     
+       x=global[0]; y=global[1]; z=global[2];
+       pm->SetPoint(i,x,y,z);
+      }
+    }
+  }
+
+  pm->SetMarkerSize(1); pm->SetMarkerColor(10); pm->SetMarkerStyle(1);
+  pm->Draw();
+  
+
+  Int_t ntracks = tarray.GetEntriesFast();
+
+  for (i = 0; i < ntracks; i++) {
+    AliTRDtrack *pt = (AliTRDtrack *) tarray.UncheckedAt(i), &t=*pt;
+    Int_t nclusters = t.GetNclusters();
+    cerr<<"in track "<<i<<" found "<<nclusters<<" clusters"<<endl;
+   
+    TPolyMarker3D *pm = new TPolyMarker3D(nclusters);
+    for(j = 0; j < nclusters; j++) {
+      index = t.GetClusterIndex(j);
+      AliTRDrecPoint *rp = (AliTRDrecPoint *) RecPointsArray->UncheckedAt(index);    
+      local[0]=rp->GetLocalRow(); 
+      local[1]=rp->GetLocalCol(); 
+      local[2]=rp->GetLocalTime();
+      det=rp->GetDetector();
+
+      if(fGeom->Local2Global(det,local,global)) {     
+       x=global[0]; y=global[1]; z=global[2];
+       pm->SetPoint(j,x,y,z);
+      }
+    } 
+    pm->SetMarkerSize(1); pm->SetMarkerColor(i%6+3); pm->SetMarkerStyle(1);
+    pm->Draw();
+  }
+  
+  TGeometry *geom=(TGeometry*)gafl->Get("AliceGeom");
+  geom->Draw("same");
+
+  c1->Modified(); 
+  
+  c1->Update();
+  
+  gafl->Close();
+
+}