/*
$Log$
+Revision 1.14 2002/10/14 14:57:43 hristov
+Merging the VirtualMC branch to the main development branch (HEAD)
+
Revision 1.11.6.2 2002/10/11 07:26:37 hristov
Updating VirtualMC to v3-09-02
Revision 1.13 2002/09/18 09:20:53 cblume
Write the parameter class into the cluster file
+Revision 1.12 2002/06/12 09:54:35 cblume
+Update of tracking code provided by Sergei
+
Revision 1.11 2001/11/27 08:50:33 hristov
BranchOld replaced by Branch
AliTRDgeometry *geo = fTRD->GetGeometry();
geo->SetName("TRDgeometry");
geo->Write();
- fPar->Write();
+ fPar->Write();
return kTRUE;
AliTRDclusterizerV1 *clusterizer =
new AliTRDclusterizerV1("clusterizer","Clusterizer class");
- // Define the parameter object
- // If no external parameter object is defined,
- // default parameter will be used
- AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
- ,"TRD parameter class");
- clusterizer->SetParameter(parameter);
+ // Read the parameter
+ TFile *parfile = TFile::Open(infile);
+ AliTRDparameter *par = (AliTRDparameter *) parfile->Get("TRDparameter");
+ par->ReInit();
+ clusterizer->SetParameter(par);
+
+ // Set the parameter
+ clusterizer->SetVerbose(1);
// Open the AliRoot file
+ // clusterizer->Open(infile,0);
clusterizer->Open(infile,outfile,0);
+
// Load the digits
clusterizer->ReadDigits();
+ clusterizer->Dump();
// Find the cluster
clusterizer->MakeClusters();
* provided "as is" without express or implied warranty. *
**************************************************************************/
-///////////////////////////////////////////////////////////////////////////////
-// //
-// TRD MC track //
-// Used for efficiency estimates and matching of reconstructed tracks //
-// to MC particles //
-// //
-///////////////////////////////////////////////////////////////////////////////
-
#include "AliTRDmcTrack.h"
#include "AliTRDgeometry.h"
+
ClassImp(AliTRDmcTrack)
//_____________________________________________________________________________
fCharge = 0;
fN = 0;
- for(Int_t i=0; i<200; i++) fIndex[i]=-1;
+ for(Int_t i=0; i<kMAX_CLUSTERS_PER_MC_TRACK; i++) fIndex[i]=-1;
for(Int_t i=0; i<6; i++) {
for(Int_t j=0; j<3; j++) {
- fPin[i][j]=0.;
- fPout[i][j] = 0.;
+ Pin[i][j]=0.;
+ Pout[i][j] = 0.;
+ XYZin[i][j]=0.;
+ XYZout[i][j] = 0.;
}
}
fPDG = pdg;
fN = 0;
- for(Int_t i=0; i<200; i++) fIndex[i]=-1;
+ for(Int_t i=0; i<kMAX_CLUSTERS_PER_MC_TRACK; i++) fIndex[i]=-1;
for(Int_t i=0; i<6; i++) {
for(Int_t j=0; j<3; j++) {
- fPin[i][j]=0.;
- fPout[i][j] = 0.;
+ Pin[i][j]=0.;
+ Pout[i][j] = 0.;
+ XYZin[i][j]=0.;
+ XYZout[i][j] = 0.;
}
}
}
//_____________________________________________________________________________
-void AliTRDmcTrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz,
- Int_t opt) const
+void AliTRDmcTrack::GetPxPyPzXYZ(Double_t& px, Double_t& py, Double_t& pz,
+ Double_t& x, Double_t& y, Double_t& z,
+ Int_t opt) const
{
//
- // Returns momentum components at the entrance (opt >= 0), or
- // exit (opt < 0) of TRD.
+ // Returns track momentum components and coordinates at the entrance
+ // (opt >= 0), or exit (opt < 0) of TRD.
//
Int_t i;
if(opt >= 0) {
for(i = 0; i < AliTRDgeometry::Nplan(); i++) {
- if( fPin[i][0] * fPin[i][0]
- + fPin[i][1] * fPin[i][1]
- + fPin[i][2] * fPin[i][2] > 0.0005) break;
+ if( Pin[i][0] * Pin[i][0]
+ + Pin[i][1] * Pin[i][1]
+ + Pin[i][2] * Pin[i][2] > 0.0005) break;
}
- px = fPin[0][0]; py = fPin[0][1]; pz = fPin[0][2];
+ px = Pin[i][0]; py = Pin[i][1]; pz = Pin[i][2];
+ x = XYZin[i][0]; y = XYZin[i][1]; z = XYZin[i][2];
}
else {
for(i = AliTRDgeometry::Nplan() - 1; i >= 0; i--) {
- if( fPout[i][0] * fPout[i][0]
- + fPout[i][1] * fPout[i][1]
- + fPout[i][2] * fPout[i][2] > 0.0005) break;
+ if( Pout[i][0] * Pout[i][0]
+ + Pout[i][1] * Pout[i][1]
+ + Pout[i][2] * Pout[i][2] > 0.0005) break;
}
- px = fPout[i][0]; py = fPout[i][1]; pz = fPout[i][2];
+ px = Pout[i][0]; py = Pout[i][1]; pz = Pout[i][2];
+ x = XYZout[i][0]; y = XYZout[i][1]; z = XYZout[i][2];
}
-
return;
-
}
//_____________________________________________________________________________
//
if(opt >= 0) {
- px = fPin[plane][0]; py = fPin[plane][1]; pz = fPin[plane][2];
+ px = Pin[plane][0]; py = Pin[plane][1]; pz = Pin[plane][2];
}
else {
- px = fPout[plane][0]; py = fPout[plane][1]; pz = fPout[plane][2];
+ px = Pout[plane][0]; py = Pout[plane][1]; pz = Pout[plane][2];
}
-
return;
-
}
+
+
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-///////////////////////////////////////////////////////////////////////////////
-// //
-// TRD MC track //
-// Used for efficiency estimates and matching of reconstructed tracks //
-// to MC particles //
-// //
-///////////////////////////////////////////////////////////////////////////////
-
#include <TObject.h>
class AliTRDgeometry;
+const Int_t kMAX_CLUSTERS_PER_MC_TRACK=210;
+
class AliTRDmcTrack : public TObject {
// Represents TRD related info about generated track
AliTRDmcTrack(Int_t label, Bool_t primary, Float_t mass, Int_t charge, Int_t pdg);
void SetPin(Int_t plane, Double_t px, Double_t py, Double_t pz)
- { fPin[plane][0] = px; fPin[plane][1] = py; fPin[plane][2] = pz; }
+ { Pin[plane][0] = px; Pin[plane][1] = py; Pin[plane][2] = pz; }
void SetPout(Int_t plane, Double_t px, Double_t py, Double_t pz)
- { fPout[plane][0] = px; fPout[plane][1] = py; fPout[plane][2] = pz; }
+ {Pout[plane][0] = px; Pout[plane][1] = py; Pout[plane][2] = pz;}
+
+ void SetXYZin(Int_t plane, Double_t x, Double_t y, Double_t z)
+ { XYZin[plane][0] = x; XYZin[plane][1] = y; XYZin[plane][2] = z; }
+
+ void SetXYZout(Int_t plane, Double_t x, Double_t y, Double_t z)
+ { XYZout[plane][0] = x; XYZout[plane][1] = y; XYZout[plane][2] = z; }
+
+ void GetPxPyPzXYZ(Double_t &px, Double_t &py, Double_t &pz,
+ Double_t &x, Double_t &y, Double_t &z,
+ Int_t opt = 0) const;
- void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz, Int_t opt = 0) const;
void GetPlanePxPyPz(Double_t &px, Double_t &py, Double_t &pz
,Int_t plane, Int_t opt = 0) const;
- void Update(Int_t index) { if (fN < 199) fIndex[fN++] = index; }
+ void GetXYZin(Int_t plane, Double_t &x, Double_t &y, Double_t &z) const
+ { x = XYZin[plane][0]; y = XYZin[plane][1]; z = XYZin[plane][2]; return; }
+
+ void GetXYZout(Int_t plane, Double_t &x, Double_t &y, Double_t &z) const
+ {x = XYZout[plane][0]; y = XYZout[plane][1]; z = XYZout[plane][2]; return;}
+
+ void Update(Int_t index) {
+ if (fN < kMAX_CLUSTERS_PER_MC_TRACK-1) fIndex[fN++] = index; }
Int_t GetTrackIndex() const { return fLab; }
Bool_t IsPrimary() const { return fPrimary; }
Int_t fCharge; // Charge of the MC track
Int_t fPDG; // PDG code of the MC track
- Int_t fN; // Number of TRD clusters associated with the track
- Int_t fIndex[200]; // Indices of these clusters
+ Int_t fN; // Number of TRD clusters associated with the track
+ Int_t fIndex[kMAX_CLUSTERS_PER_MC_TRACK]; // Indices of these clusters
- Double_t fPin[6][3]; // Px,Py,Pz at the entrance of each TRD plane
- Double_t fPout[6][3]; // Px,Py,Pz at the exit of each TRD plane
+ Double_t Pin[6][3]; // Px,Py,Pz at the entrance of each TRD plane
+ Double_t Pout[6][3]; // Px,Py,Pz at the exit of each TRD plane
+
+ Double_t XYZin[6][3]; // x,y,z at the entrance of the TRD
+ Double_t XYZout[6][3]; // x,y,z at the exit of the TRD
ClassDef(AliTRDmcTrack,1) // TRD MC track
virtual void BuildPhysics() { };
virtual void ProcessEvent();
virtual void ProcessRun(Int_t nevent) { };
- //virtual AliMCGeomType GetMCGeomType() const { return kGeant3; }
+ //virtual TMCGeomType GetMCGeomType() const { return kGeant3; }
// External Decayer
virtual void SetExternalDecayer(AliDecayer* decayer) { };
/*
$Log$
+Revision 1.13 2002/10/22 15:53:08 alibrary
+Introducing Riostream.h
+
Revision 1.12 2002/10/14 14:57:44 hristov
Merging the VirtualMC branch to the main development branch (HEAD)
Revision 1.8.10.2 2002/07/24 10:09:31 alibrary
Updating VirtualMC
-Revision 1.11 2002/06/13 12:09:58 hristov
+RRevision 1.11 2002/06/13 12:09:58 hristov
Minor corrections
Revision 1.10 2002/06/12 09:54:35 cblume
+
//_____________________________________________________________________________
Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
{
}
+
//_____________________________________________________________________________
-Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index)
+Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index, Double_t h01)
{
// Assignes found cluster to the track and updates track information
- Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
- r00+=fCyy; r01+=fCzy; r11+=fCzz;
+
+ Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2()*100.0;
+
+ r00+=(fCyy+2.0*h01*fCzy+h01*h01*fCzz);
+ r01+=(fCzy+h01*fCzz);
+ 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 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 k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
+ Double_t k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
+ Double_t k20=fCcy*r00+fCcz*(r01+h01*r00),k21=fCcy*r01+fCcz*(r11+h01*r01);
+ Double_t k30=fCey*r00+fCez*(r01+h01*r00),k31=fCey*r01+fCez*(r11+h01*r01);
+ Double_t k40=fCty*r00+fCtz*(r01+h01*r00),k41=fCty*r01+fCtz*(r11+h01*r01);
Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
+
+ dy=dy+h01*dz;
+
Double_t cur=fC + k20*dy + k21*dz, eta=fE + k30*dy + k31*dz;
if (TMath::Abs(cur*fX-eta) >= 0.99999) {
- Int_t n=GetNumberOfClusters();
- if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
+ // if (fN>4) cerr<<fN<<" AliTRDtrack warning: Filtering failed !\n";
return 0;
}
fE = eta;
fT += k40*dy + k41*dz;
+
+ k01+=h01*k00;
+ k11+=h01*k10;
+ k21+=h01*k20;
+ k31+=h01*k30;
+ k41+=h01*k40;
+
Double_t c01=fCzy, c02=fCcy, c03=fCey, c04=fCty;
Double_t c12=fCcz, c13=fCez, c14=fCtz;
SetNumberOfClusters(n+1);
SetChi2(GetChi2()+chisq);
+
// cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
return 1;
+
}
+
//_____________________________________________________________________________
Int_t AliTRDtrack::Rotate(Double_t alpha)
{
//_____________________________________________________________________________
-Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c) const
+Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
{
- /*
+
Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
+
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";
+ Int_t n=GetNumberOfClusters();
+ if (n>4) cerr<<n<<" 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;
- */
-
- Double_t dy=c->GetY() - fY;
- Double_t r00=c->GetSigmaY2();
+ dy=dy+h01*dz;
- return (dy*dy)/r00;
+ return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
-}
+}
//_________________________________________________________________________
Double_t GetP() const {
return TMath::Abs(GetPt())*sqrt(1.+GetTgl()*GetTgl());
}
- Double_t GetPredictedChi2(const AliTRDcluster*) const ;
+ Double_t GetPredictedChi2(const AliTRDcluster*, Double_t h01) const ;
Double_t GetPt() const {return 1./Get1Pt();}
void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const ;
void GetGlobalXYZ(Double_t &x, Double_t &y, Double_t &z) const ;
void SetSeedLabel(Int_t lab) { fSeedLab=lab; }
- Int_t Update(const AliTRDcluster* c, Double_t chi2, UInt_t i);
+ Int_t Update(const AliTRDcluster* c, Double_t chi2, UInt_t i,
+ Double_t h01);
+
+
protected:
/*
$Log$
+Revision 1.19 2002/10/22 15:53:08 alibrary
+Introducing Riostream.h
+
Revision 1.18 2002/10/14 14:57:44 hristov
Merging the VirtualMC branch to the main development branch (HEAD)
for (Int_t i=0; i<n; i++) {
seedTree->GetEvent(i);
seed->ResetCovariance();
- fSeeds->AddLast(new AliTRDtrack(*seed));
+ AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
+ fSeeds->AddLast(tr);
fNseeds++;
}
delete seed;
+ delete seedTree;
}
}
out->cd();
+
+
// find tracks from loaded seeds
Int_t nseed=fSeeds->GetEntriesFast();
}
iotrack_trd = pt;
trd_tree.Fill();
- cerr<<found++<<'\r';
+ cout<<found++<<'\r';
if(PropagateToTPC(t)) {
AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
Int_t inner=outer-gap;
nseed=fSeeds->GetEntriesFast();
- printf("\n initial number of seeds %d\n", nseed);
MakeSeeds(inner, outer, turn);
nseed=fSeeds->GetEntriesFast();
- printf("\n number of seeds after MakeSeeds %d\n", nseed);
-
-
+ printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
+ turn, i, nseed);
+
for (Int_t i=0; i<nseed; i++) {
AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
FollowProlongation(t,innerTB);
UseClusters(&t);
CookLabel(pt, 1-fLabelFraction);
t.CookdEdx();
- cerr<<found++<<'\r';
+ cout<<found++<<'\r';
iotrack_trd = pt;
trd_tree.Fill();
if(PropagateToTPC(t)) {
Int_t foundClr = s.GetNumberOfClusters();
Int_t last_tb = fTrSec[0]->GetLayerNumber(s.GetX());
- printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
- i, foundClr, expectedClr, foundMin);
+ // printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
+ // i, foundClr, expectedClr, foundMin);
if (foundClr >= foundMin) {
s.CookdEdx();
minDY = TMath::Abs(c->GetY() - y);
wYcorrect = c->GetY();
wZcorrect = c->GetZ();
- wChi2 = t.GetPredictedChi2(c);
+
+ Double_t h01 = GetTiltFactor(c);
+ wChi2 = t.GetPredictedChi2(c, h01);
}
}
if (c->GetY() > y+road) break;
if (c->IsUsed() > 0) continue;
if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
- Double_t chi2=t.GetPredictedChi2(c);
+
+ Double_t h01 = GetTiltFactor(c);
+ Double_t chi2=t.GetPredictedChi2(c,h01);
if (chi2 > max_chi2) continue;
max_chi2=chi2;
if (c->IsUsed() > 0) continue;
if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
- Double_t chi2=t.GetPredictedChi2(c);
+ Double_t h01 = GetTiltFactor(c);
+ Double_t chi2=t.GetPredictedChi2(c, h01);
if (chi2 > max_chi2) continue;
max_chi2=chi2;
if (cl) {
wYclosest = cl->GetY();
wZclosest = cl->GetZ();
+ Double_t h01 = GetTiltFactor(cl);
t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
- if(!t.Update(cl,max_chi2,index)) {
+ if(!t.Update(cl,max_chi2,index,h01)) {
if(!try_again--) return 0;
}
else try_again=fMaxGap;
minDY = TMath::Abs(c->GetY() - y);
wYcorrect = c->GetY();
wZcorrect = c->GetZ();
- wChi2 = t.GetPredictedChi2(c);
+
+ Double_t h01 = GetTiltFactor(c);
+ wChi2 = t.GetPredictedChi2(c, h01);
}
}
if (c->GetY() > y+road) break;
if (c->IsUsed() > 0) continue;
if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
- Double_t chi2=t.GetPredictedChi2(c);
+
+ Double_t h01 = GetTiltFactor(c);
+ Double_t chi2=t.GetPredictedChi2(c,h01);
if (chi2 > max_chi2) continue;
max_chi2=chi2;
if (c->IsUsed() > 0) continue;
if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
- Double_t chi2=t.GetPredictedChi2(c);
+ Double_t h01 = GetTiltFactor(c);
+ Double_t chi2=t.GetPredictedChi2(c,h01);
if (chi2 > max_chi2) continue;
max_chi2=chi2;
wZclosest = cl->GetZ();
t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
- if(!t.Update(cl,max_chi2,index)) {
+ Double_t h01 = GetTiltFactor(cl);
+ if(!t.Update(cl,max_chi2,index,h01)) {
if(!try_again--) return 0;
}
else try_again=fMaxGap;
Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;
+
+ // Tilt changes
+ Double_t h01 = GetTiltFactor(r1[is]);
+ sy1=sy1+sz1*h01*h01;
+ Double_t syz=sz1*h01;
+ // end of tilt changes
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 f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
c[0]=sy1;
- c[1]=0.; c[2]=sz1;
+ // c[1]=0.; c[2]=sz1;
+ c[1]=syz; c[2]=sz1*100;
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;
return m;
}
+//---------------------------------------------------------
+Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
+//
+// Returns correction factor for tilted pads geometry
+//
+ Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
+ Int_t det = c->GetDetector();
+ Int_t plane = fGeom->GetPlane(det);
+ if((plane == 0) || (plane == 2) || (plane == 4)) h01=-h01;
+
+ return h01;
+}
void SetEventNumber(Int_t event) { fEvent = event; }
void SetAddTRDseeds() { fAddTRDseeds = kTRUE; }
+ Double_t GetTiltFactor(const AliTRDcluster* c);
+
void ReadClusters(TObjArray *array, const Char_t *filename);
Int_t CookSectorIndex(Int_t gs) { return kTRACKING_SECTORS - 1 - gs; }