//-----------------------------------------------------------------
// Implementation of the vertexer from tracks
+// It accepts V2 and ESD tracks
//
-// Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it
-// M.Masera, Torino, massimo.masera@to.infn.it
+// Origin: A.Dainese, Padova,
+// andrea.dainese@pd.infn.it
+// M.Masera, Torino,
+// massimo.masera@to.infn.it
//-----------------------------------------------------------------
//---- standard headers ----
#include <Riostream.h>
//---- Root headers --------
+#include <TKey.h>
#include <TFile.h>
#include <TTree.h>
-#include <TVector3.h>
#include <TMatrixD.h>
-#include <TObjArray.h>
-#include <TRandom.h> // TEMPORARY !!!!!!!
//---- AliRoot headers -----
-#include <AliRun.h>
-#include "AliKalmanTrack.h"
#include "AliITSStrLine.h"
#include "AliITStrackV2.h"
-#include "AliITSVertex.h"
+#include "AliESDVertex.h"
#include "AliITSVertexerTracks.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
ClassImp(AliITSVertexerTracks)
//
// Default constructor
//
-
+ fInFile = 0;
+ fOutFile = 0;
SetVtxStart();
SetMinTracks();
- SetUseThrustFrame();
- SetPhiThrust();
fTrksToSkip = 0;
fNTrksToSkip = 0;
for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
}
-
//----------------------------------------------------------------------------
AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
- Double_t field,
- Double_t xStart,Double_t yStart,
- Int_t useThFr)
- :AliITSVertexer(inFile,outFile) {
+ Double_t field,
+ Int_t fEv,Int_t lEv,
+ Double_t xStart,Double_t yStart) {
//
// Standard constructor
//
-
+ fCurrentVertex = 0;
+ fInFile = inFile;
+ fOutFile = outFile;
+ SetFirstEvent(fEv);
+ SetLastEvent(lEv);
SetField(field);
SetVtxStart(xStart,yStart);
SetMinTracks();
- SetUseThrustFrame(useThFr);
- SetPhiThrust();
fTrksToSkip = 0;
fNTrksToSkip = 0;
- for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
+ for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
+ SetDebug();
}
//----------------------------------------------------------------------------
+AliITSVertexerTracks::AliITSVertexerTracks(Double_t field, TString fn,
+ Double_t xStart,Double_t yStart)
+ :AliITSVertexer(fn) {
+//
+// Alternative constructor
+//
+ fInFile = 0;
+ fOutFile = 0;
+ SetField(field);
+ SetVtxStart(xStart,yStart);
+ SetMinTracks();
+ fTrksToSkip = 0;
+ fNTrksToSkip = 0;
+ for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
+}
+//______________________________________________________________________
+AliITSVertexerTracks::AliITSVertexerTracks(const AliITSVertexerTracks &vtxr) : AliITSVertexer(vtxr) {
+ // Copy constructor
+ // Copies are not allowed. The method is protected to avoid misuse.
+ Error("AliITSVertexerTracks","Copy constructor not allowed\n");
+}
+
+//______________________________________________________________________
+AliITSVertexerTracks& AliITSVertexerTracks::operator=(const AliITSVertexerTracks& /* vtxr */){
+ // Assignment operator
+ // Assignment is not allowed. The method is protected to avoid misuse.
+ Error("= operator","Assignment operator not allowed\n");
+ return *this;
+}
+
+//-----------------------------------------------------------------------------
+AliITSVertexerTracks::~AliITSVertexerTracks() {
+ // Default Destructor
+ // The objects poited by the following pointers are not owned
+ // by this class and are not deleted
+
+ fCurrentVertex = 0;
+ fInFile = 0;
+ fOutFile = 0;
+}
+//-----------------------------------------------------------------------------
Bool_t AliITSVertexerTracks::CheckField() const {
//
// Check if the conv. const. has been set
// Check if the conv. const. has been set
if(!CheckField()) return;
+ TDirectory *curdir = 0;
// loop over events
for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
}
if(fDebug) fCurrentVertex->PrintStatus();
+
+ // write vertex to file
TString vtxName = "Vertex_";
vtxName += ev;
fCurrentVertex->SetName(vtxName.Data());
fCurrentVertex->SetTitle("VertexerTracks");
- WriteCurrentVertex();
+ //WriteCurrentVertex();
+ curdir = gDirectory;
+ fOutFile->cd();
+ fCurrentVertex->Write();
+ curdir->cd();
+ fCurrentVertex = 0;
+ } // loop over events
+
+ return;
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::FindVerticesESD() {
+//
+// Vertices for all events from fFirstEvent to fLastEvent
+//
+
+ // Check if the conv. const. has been set
+ if(!CheckField()) return;
+
+ TDirectory *curdir = 0;
+
+ fInFile->cd();
+ TKey *key=0;
+ TIter next(fInFile->GetListOfKeys());
+ // loop on events in file
+ while ((key=(TKey*)next())!=0) {
+ AliESD *esdEvent=(AliESD*)key->ReadObj();
+ if(!esdEvent) {
+ printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n");
+ return;
+ }
+ Int_t ev = (Int_t)esdEvent->GetEventNumber();
+ if(ev<fFirstEvent || ev>fLastEvent) { delete esdEvent; continue; }
+ if(ev % 100 == 0 || fDebug)
+ printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);
+
+ FindVertexForCurrentEvent(esdEvent);
+
+ if(!fCurrentVertex) {
+ printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
+ continue;
+ }
+
+ cout<<"VERTICE TROVATO\n";
+ if(fDebug) fCurrentVertex->PrintStatus();
+
+ // write the ESD to file
+ curdir = gDirectory;
+ Char_t ename[100];
+ sprintf(ename,"%d",ev);
+ fOutFile->cd();
+ esdEvent->Dump();
+ esdEvent->Write(ename,TObject::kOverwrite);
+ curdir->cd();
+ fCurrentVertex = 0;
} // loop over events
return;
printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
- printf(" Using Thrust Frame: %d fPhiThrust = %f\n",fUseThrustFrame,fPhiThrust);
return;
}
//----------------------------------------------------------------------------
-AliITSVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
+AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
//
// Vertex for current event
//
// VERTEX FITTER
ComputeMaxChi2PerTrack(nTrks);
- if(fUseThrustFrame) ThrustFinderXY();
- if(fDebug) printf(" thrust found: phi = %f\n",fPhiThrust);
VertexFitter();
if(fDebug) printf(" vertex fit completed\n");
- TString vtxName;
- vtxName = "Vertex_";
- vtxName += evnumb;
- fCurrentVertex->SetName(vtxName.Data());
return fCurrentVertex;
}
-//---------------------------------------------------------------------------
-void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
-//
-// Mark the tracks not ot be used in the vertex finding
-//
- fNTrksToSkip = n;
- fTrksToSkip = new Int_t[n];
- for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
- return;
-}
//----------------------------------------------------------------------------
-Double_t AliITSVertexerTracks::SumPl(TTree &momTree,Double_t phi) const {
+AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
+{
//
-// Function to be maximized for thrust determination
+// Vertex for current ESD event
//
- TVector3 u(1.,1.,0);
- u.SetMag(1.);
- u.SetPhi(phi);
- Double_t pl;
-
- Double_t sum = 0.;
-
- TVector3* mom=0;
- momTree.SetBranchAddress("momenta",&mom);
- Int_t entries = (Int_t)momTree.GetEntries();
-
- for(Int_t i=0; i<entries; i++) {
- momTree.GetEvent(i);
- pl = mom->Dot(u);
- if(pl>0.) sum += pl;
- }
-
- delete mom;
+ fCurrentVertex = 0;
+ Double_t vtx[3],cvtx[6];
+
+ // put tracks reco in ITS in a tree
+ Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
+ TTree *trkTree = new TTree("TreeT_ITS","its tracks");
+ AliITStrackV2 *itstrack = 0;
+ trkTree->Branch("tracks","AliITStrackV2",&itstrack,entr,0);
+
+ for(Int_t i=0; i<entr; i++) {
+ AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(i);
+ if(!esdTrack->GetStatus()&AliESDtrack::kITSin)
+ { delete esdTrack; continue; }
+ try {
+ itstrack = new AliITStrackV2(*esdTrack);
+ }
+ catch (const Char_t *msg) {
+ Warning("FindVertexForCurrentEvent",msg);
+ delete esdTrack;
+ continue;
+ }
- return sum;
-}
-//---------------------------------------------------------------------------
-void AliITSVertexerTracks::ThrustFinderXY() {
-//
-// This function looks for the thrust direction, \vec{u}, in the (x,y) plane.
-// The following function is maximized:
-// \Sum_{\vec{p}\cdot\vec{u}} \vec{p}\cdot\vec{u} / \Sum |\vec{p}|
-// where \vec{p} = (p_x,p_y)
-//
- Double_t pt,alpha,phi;
- Double_t totPt;
+ trkTree->Fill();
+ itstrack = 0;
+ delete esdTrack;
+ }
+ delete itstrack;
- // tree for thrust determination
- TVector3 *ioMom = new TVector3;
- TTree *t = new TTree("Tree_Momenta","Tree with momenta");
- t->Branch("momenta","TVector3",&ioMom);
- totPt = 0.;
+ // preselect tracks and propagate them to initial vertex position
+ Int_t nTrks = PrepareTracks(*trkTree);
+ delete trkTree;
+ if(fDebug) printf(" tracks prepared: %d\n",nTrks);
+ if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
- AliITStrackV2 *itstrack = 0;
- Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
+ // Set initial vertex position from ESD
+ esdEvent->GetVertex()->GetXYZ(vtx);
+ SetVtxStart(vtx[0],vtx[1]);
- // loop on tracks
- for(Int_t i=0; i<arrEntries; i++) {
- itstrack = (AliITStrackV2*)fTrkArray.At(i);
- // momentum of the track at the vertex
- pt = 1./TMath::Abs(itstrack->Get1Pt());
- alpha = itstrack->GetAlpha();
- phi = alpha+TMath::ASin(itstrack->GetSnp());
- ioMom->SetX(pt*TMath::Cos(phi));
- ioMom->SetY(pt*TMath::Sin(phi));
- ioMom->SetZ(0.);
-
- totPt += ioMom->Pt();
- t->Fill();
- } // end loop on tracks
-
- Double_t tValue=0.,tPhi=0.;
- Double_t maxSumPl = 0.;
- Double_t thisSumPl;
- Double_t dPhi;
- Int_t nSteps,iStep;
-
- phi = 0.;
- nSteps = 100;
- dPhi = 2.*TMath::Pi()/(Double_t)nSteps;
-
- for(iStep=0; iStep<nSteps; iStep++) {
- phi += dPhi;
- thisSumPl = SumPl(*t,phi);
- if(thisSumPl > maxSumPl) {
- maxSumPl = thisSumPl;
- tPhi = phi;
- }
- }
-
- phi = tPhi-dPhi;
- nSteps = 10;
- dPhi /= 5.;
+ // VERTEX FINDER
+ VertexFinder();
- for(iStep=0; iStep<nSteps; iStep++) {
- phi += dPhi;
- thisSumPl = SumPl(*t,phi);
- if(thisSumPl > maxSumPl) {
- maxSumPl = thisSumPl;
- tPhi = phi;
- }
- }
+ // VERTEX FITTER
+ ComputeMaxChi2PerTrack(nTrks);
+ VertexFitter();
+ if(fDebug) printf(" vertex fit completed\n");
- tValue = 2.*maxSumPl/totPt;
- if(tPhi<0.) tPhi += 2.*TMath::Pi();
- if(tPhi>2.*TMath::Pi()) tPhi -= 2.*TMath::Pi();
+ // store vertex information in ESD
+ fCurrentVertex->GetXYZ(vtx);
+ fCurrentVertex->GetCovMatrix(cvtx);
- SetPhiThrust(tPhi);
+ Double_t tp[3];
+ esdEvent->GetVertex()->GetTruePos(tp);
+ fCurrentVertex->SetTruePos(tp);
- delete t;
- delete ioMom;
+ esdEvent->SetVertex(fCurrentVertex);
- return;
+ cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
+ return fCurrentVertex;
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
+//
+// Mark the tracks not ot be used in the vertex finding
+//
+ fNTrksToSkip = n;
+ fTrksToSkip = new Int_t[n];
+ for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
+ return;
}
//---------------------------------------------------------------------------
void AliITSVertexerTracks::TooFewTracks() {
// When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
// and the number of tracks to -1
//
- fCurrentVertex = new AliITSVertex(0.,0.,-1);
+ fCurrentVertex = new AliESDVertex(0.,0.,-1);
return;
}
//---------------------------------------------------------------------------
Int_t ncombi = 0;
AliITStrackV2 *track1;
AliITStrackV2 *track2;
+ AliITSStrLine line1;
+ AliITSStrLine line2;
for(Int_t i=0; i<nacc; i++){
track1 = (AliITStrackV2*)fTrkArray.At(i);
if(fDebug>5){
Double_t pos1[3];
Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
- AliITSStrLine *line1 = new AliITSStrLine(pos1,mom1);
+ line1.SetP0(pos1);
+ line1.SetCd(mom1);
for(Int_t j=i+1; j<nacc; j++){
track2 = (AliITStrackV2*)fTrkArray.At(j);
Double_t mom2[3];
Double_t pos2[3];
mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
- AliITSStrLine *line2 = new AliITSStrLine(pos2,mom2);
+ line2.SetP0(pos2);
+ line2.SetCd(mom2);
Double_t crosspoint[3];
- Int_t retcode = line2->Cross(line1,crosspoint);
+ Int_t retcode = line2.Cross(&line1,crosspoint);
if(retcode<0){
if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
if(fDebug>10)cout<<"bad intersection\n";
- line1->PrintStatus();
- line2->PrintStatus();
+ line1.PrintStatus();
+ line2.PrintStatus();
}
else {
ncombi++;
if(fDebug>10)cout<<"\n Cross point: ";
if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
}
- delete line2;
}
- delete line1;
}
if(ncombi>0){
for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
return;
}
-
//---------------------------------------------------------------------------
void AliITSVertexerTracks::VertexFitter() {
//
Int_t i,j,k,step=0;
TMatrixD rv(3,1);
- TMatrixD V(3,3);
+ TMatrixD vV(3,3);
rv(0,0) = fInitPos[0];
rv(1,0) = fInitPos[1];
rv(2,0) = 0.;
chi2 = 0.;
nUsedTrks = 0;
- TMatrixD SumWiri(3,1);
- TMatrixD SumWi(3,3);
+ TMatrixD sumWiri(3,1);
+ TMatrixD sumWi(3,3);
for(i=0; i<3; i++) {
- SumWiri(i,0) = 0.;
- for(j=0; j<3; j++) SumWi(j,i) = 0.;
+ sumWiri(i,0) = 0.;
+ for(j=0; j<3; j++) sumWi(j,i) = 0.;
}
// loop on tracks
alpha = t->GetAlpha();
xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
t->PropagateTo(xlStart,0.,0.); // to vtxSeed
- rotAngle = alpha-fPhiThrust;
+ rotAngle = alpha;
if(alpha<0.) rotAngle += 2.*TMath::Pi();
cosRot = TMath::Cos(rotAngle);
sinRot = TMath::Sin(rotAngle);
ri(2,0) = t->GetZ();
// matrix to go from global (x,y,z) to local (y,z);
- TMatrixD Qi(2,3);
- Qi(0,0) = -sinRot;
- Qi(0,1) = cosRot;
- Qi(0,2) = 0.;
- Qi(1,0) = 0.;
- Qi(1,1) = 0.;
- Qi(1,2) = 1.;
+ TMatrixD qQi(2,3);
+ qQi(0,0) = -sinRot;
+ qQi(0,1) = cosRot;
+ qQi(0,2) = 0.;
+ qQi(1,0) = 0.;
+ qQi(1,1) = 0.;
+ qQi(1,2) = 1.;
// covariance matrix of local (y,z) - inverted
- TMatrixD Ui(2,2);
+ TMatrixD uUi(2,2);
t->GetExternalCovariance(cc);
- Ui(0,0) = cc[0];
- Ui(0,1) = cc[1];
- Ui(1,0) = cc[1];
- Ui(1,1) = cc[2];
+ uUi(0,0) = cc[0];
+ uUi(0,1) = cc[1];
+ uUi(1,0) = cc[1];
+ uUi(1,1) = cc[2];
- // weights matrix: Wi = QiT * UiInv * Qi
- if(Ui.Determinant() <= 0.) continue;
- TMatrixD UiInv(TMatrixD::kInverted,Ui);
- TMatrixD UiInvQi(UiInv,TMatrixD::kMult,Qi);
- TMatrixD Wi(Qi,TMatrixD::kTransposeMult,UiInvQi);
+ // weights matrix: wWi = qQiT * uUiInv * qQi
+ if(uUi.Determinant() <= 0.) continue;
+ TMatrixD uUiInv(TMatrixD::kInverted,uUi);
+ TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
+ TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
// track chi2
TMatrixD deltar = rv; deltar -= ri;
- TMatrixD Wideltar(Wi,TMatrixD::kMult,deltar);
- chi2i = deltar(0,0)*Wideltar(0,0)+
- deltar(1,0)*Wideltar(1,0)+
- deltar(2,0)*Wideltar(2,0);
+ TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
+ chi2i = deltar(0,0)*wWideltar(0,0)+
+ deltar(1,0)*wWideltar(1,0)+
+ deltar(2,0)*wWideltar(2,0);
if(step==1 && chi2i > fMaxChi2PerTrack) {
// add to total chi2
chi2 += chi2i;
- TMatrixD Wiri(Wi,TMatrixD::kMult,ri);
+ TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
- SumWiri += Wiri;
- SumWi += Wi;
+ sumWiri += wWiri;
+ sumWi += wWi;
nUsedTrks++;
} // end loop on tracks
continue;
}
- Double_t determinant = SumWi.Determinant();
+ Double_t determinant = sumWi.Determinant();
//cerr<<" determinant: "<<determinant<<endl;
if(determinant < 100.) {
printf("det(V) = 0\n");
}
// inverted of weights matrix
- TMatrixD InvSumWi(TMatrixD::kInverted,SumWi);
- V = InvSumWi;
+ TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
+ vV = invsumWi;
// position of primary vertex
- rv.Mult(V,SumWiri);
+ rv.Mult(vV,sumWiri);
} // end loop on the 3 steps
position[1] = rv(1,0);
position[2] = rv(2,0);
Double_t covmatrix[6];
- covmatrix[0] = V(0,0);
- covmatrix[1] = V(0,1);
- covmatrix[2] = V(1,1);
- covmatrix[3] = V(0,2);
- covmatrix[4] = V(1,2);
- covmatrix[5] = V(2,2);
+ covmatrix[0] = vV(0,0);
+ covmatrix[1] = vV(0,1);
+ covmatrix[2] = vV(1,1);
+ covmatrix[3] = vV(0,2);
+ covmatrix[4] = vV(1,2);
+ covmatrix[5] = vV(2,2);
// store data in the vertex object
- fCurrentVertex = new AliITSVertex(fPhiThrust,position,covmatrix,chi2,nUsedTrks);
+ fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
if(fDebug) {
printf(" VertexFitter(): finish\n");
return;
}
//----------------------------------------------------------------------------
-AliITSVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
+AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
//
// Return vertex from tracks in trkTree
//
// VERTEX FITTER
ComputeMaxChi2PerTrack(nTrks);
- if(fUseThrustFrame) ThrustFinderXY();
- if(fDebug) printf(" thrust found: phi = %f\n",fPhiThrust);
VertexFitter();
if(fDebug) printf(" vertex fit completed\n");
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff