--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/*
+$Log$
+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;
+
+}
+
--- /dev/null
+#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
--- /dev/null
+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();
+
+
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/*
+$Log$
+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);
+ }
+
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/*
+$Log$
+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;
+ }
+
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/*
+$Log$
+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];
+
+}
--- /dev/null
+#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
+
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/*
+$Log$
+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);
+}
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/*
+$Log$
+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();
+
+}
+
+
+
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/*
+$Log$
+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;
+}
+
--- /dev/null
+#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
--- /dev/null
+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();
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff