//---- standard headers ----
#include <Riostream.h>
//---- Root headers --------
-#include <TKey.h>
#include <TFile.h>
#include <TTree.h>
#include <TMatrixD.h>
//---- AliRoot headers -----
-#include "AliITSStrLine.h"
-#include "AliITStrackV2.h"
#include "AliESDVertex.h"
#include "AliITSVertexerTracks.h"
#include "AliESD.h"
#include "AliESDtrack.h"
+#include "AliVertexerTracks.h"
ClassImp(AliITSVertexerTracks)
SetMinTracks();
fTrksToSkip = 0;
fNTrksToSkip = 0;
- for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
+ fDCAcut=0;
}
//----------------------------------------------------------------------------
AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
- Double_t field,
Int_t fEv,Int_t lEv,
Double_t xStart,Double_t yStart) {
//
fOutFile = outFile;
SetFirstEvent(fEv);
SetLastEvent(lEv);
- SetField(field);
SetVtxStart(xStart,yStart);
SetMinTracks();
fTrksToSkip = 0;
fNTrksToSkip = 0;
- for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
+ fDCAcut=0;
SetDebug();
}
//----------------------------------------------------------------------------
-AliITSVertexerTracks::AliITSVertexerTracks(Double_t field, TString fn,
+AliITSVertexerTracks::AliITSVertexerTracks(TString fn,
Double_t xStart,Double_t yStart)
:AliITSVertexer(fn) {
//
//
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.;
+ fDCAcut=0;
}
//______________________________________________________________________
AliITSVertexerTracks::AliITSVertexerTracks(const AliITSVertexerTracks &vtxr) : AliITSVertexer(vtxr) {
//
// Check if the conv. const. has been set
//
- AliITStrackV2 t;
- Double_t cc = t.GetConvConst();
- Double_t field = 100./0.299792458/cc;
+ Double_t field = AliTracker::GetBz(); // in kG
- if(field<0.1 || field>0.6) {
+ if(field<1 || field>6) {
printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
return kFALSE;
}
for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
- FindVertexForCurrentEvent(ev);
+ FindPrimaryVertexForCurrentEvent(ev);
if(!fCurrentVertex) {
- printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
+ printf("AliITSVertexerTracks::FindVertices(): no tracks tree for event %d\n",ev);
continue;
}
// 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;
+ // Check if the conv. const. has been set
+ if(!CheckField()) return;
+
+ TDirectory *curdir = 0;
+
+ fInFile->cd();
+ TTree *esdTree = (TTree*)fInFile->Get("esdTree");
+ if(!esdTree) {
+ printf("AliITSVertexerTracks::FindVerticesESD(): no tree in file!\n");
+ return;
+ }
+ Int_t nev = (Int_t)esdTree->GetEntries();
+ Int_t ev;
+ // loop on events in tree
+ for(Int_t i=0; i<nev; i++) {
+ AliESD *esdEvent = new AliESD;
+ esdTree->SetBranchAddress("ESD",&esdEvent);
+ if(!esdTree->GetEvent(i)) {
+ printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n");
+ delete esdEvent;
+ return;
+ }
+ 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);
+
+ FindPrimaryVertexForCurrentEvent(esdEvent);
+
+ if(!fCurrentVertex) {
+ printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
+ continue;
+ }
+
+ if(fDebug) fCurrentVertex->PrintStatus();
+
+ // write vertex to file
+ TString vtxName = "Vertex_";
+ vtxName += ev;
+ fCurrentVertex->SetName(vtxName.Data());
+ fCurrentVertex->SetTitle("VertexerTracks");
+ //WriteCurrentVertex();
+ curdir = gDirectory;
+ fOutFile->cd();
+ fCurrentVertex->Write();
+ curdir->cd();
+ fCurrentVertex = 0;
+ esdEvent = 0;
+ delete esdEvent;
+ } // end loop over events
+
+ return;
}
//----------------------------------------------------------------------------
Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
-//
-// Propagate tracks to initial vertex position and store them in a TObjArray
-//
- Double_t maxd0rphi = 3.;
+ //
+ // Propagate tracks to initial vertex position and store them in a TObjArray
+ //
+ Double_t maxd0rphi = 3.;
Double_t alpha,xlStart,d0rphi;
Int_t nTrks = 0;
Bool_t skipThis;
-
Int_t nEntries = (Int_t)trkTree.GetEntries();
+ Double_t field=AliTracker::GetBz();
+
if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
fTrkArray.Expand(nEntries);
}
if(skipThis) continue;
- AliITStrackV2 *itstrack = new AliITStrackV2;
- trkTree.SetBranchAddress("tracks",&itstrack);
+ AliESDtrack *track = new AliESDtrack;
+ trkTree.SetBranchAddress("tracks",&track);
trkTree.GetEvent(i);
-
// propagate track to vtxSeed
- alpha = itstrack->GetAlpha();
+ alpha = track->GetAlpha();
xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
- itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
- itstrack->PropagateTo(xlStart,0.,0.); // to vtxSeed
+ track->PropagateTo(xlStart,field); // to vtxSeed
// select tracks with d0rphi < maxd0rphi
- d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
- if(d0rphi > maxd0rphi) { delete itstrack; continue; }
-
- fTrkArray.AddLast(itstrack);
+ d0rphi = TMath::Abs(track->GetD(fNominalPos[0],fNominalPos[1],field));
+ if(d0rphi > maxd0rphi) { delete track; continue; }
+ fTrkArray.AddLast(track);
+
nTrks++;
+ if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<" "<<d0rphi<<endl;
+
}
-
if(fTrksToSkip) delete [] fTrksToSkip;
return nTrks;
void AliITSVertexerTracks::PrintStatus() const {
//
// Print status
-//
- printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
- printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
+// printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
return;
}
//----------------------------------------------------------------------------
-AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
+AliVertex* AliITSVertexerTracks::VertexForSelectedTracks(AliESD *esdEvent,Int_t nofCand, Int_t *trkPos, Int_t opt){
+
+ //
+ // Computes the vertex for selected tracks
+ // trkPos=vector with track positions in ESD
+ //
+ Double_t vtx[3];
+ esdEvent->GetVertex()->GetXYZ(vtx);
+ TTree *trkTree = new TTree("TreeT","tracks");
+ AliESDtrack *esdTrack = 0;
+ trkTree->Branch("tracks","AliESDtrack",&esdTrack);
+ for(Int_t i=0; i<nofCand;i++){
+ esdTrack = (AliESDtrack*)esdEvent->GetTrack(trkPos[i]);
+
+ if(!esdTrack->GetStatus()&AliESDtrack::kTPCin) continue;
+ if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
+ if(!esdTrack->GetStatus()&AliESDtrack::kITSrefit) continue;
+
+ Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
+ if(nclus<6) continue;
+ trkTree->Fill();
+ }
+ delete esdTrack;
+ Int_t nTrks = PrepareTracks(*trkTree);
+ //delete trkTree;// :-))
+ if(fDebug) printf(" tracks prepared: %d\n",nTrks);
+ if(nTrks < fMinTracks) {
+ printf("TooFewTracks\n");
+ AliVertex *theVert=new AliVertex();
+ theVert->SetDispersion(999);
+ theVert->SetNContributors(-5);
+ return theVert;
+ }
+
+ AliVertexerTracks *vertexer=new AliVertexerTracks(vtx[0],vtx[1]);
+ vertexer->SetFinderAlgorithm(opt);
+ AliVertex *theVert=(AliVertex*)vertexer->VertexForSelectedTracks(&fTrkArray);
+// beware: newvt object should be deleted by the caller
+ AliVertex *newvt = new AliVertex(*theVert);
+ delete vertexer;
+ return newvt;
+}
+//----------------------------------------------------------------------------
+AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(Int_t evnumb) {
//
// Vertex for current event
//
fCurrentVertex = 0;
// get tree with tracks from input file
- TString treeName = "TreeT_ITS_";
- treeName += evnumb;
- TTree *trkTree=(TTree*)fInFile->Get(treeName.Data());
- if(!trkTree) return fCurrentVertex;
-
-
- // get 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; }
-
- // VERTEX FINDER
- VertexFinder();
-
- // VERTEX FITTER
- ComputeMaxChi2PerTrack(nTrks);
- VertexFitter();
- if(fDebug) printf(" vertex fit completed\n");
-
- return fCurrentVertex;
+ fInFile->cd();
+ TTree *esdTree = (TTree*)fInFile->Get("esdTree");
+
+ if(!esdTree) {
+ printf("AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(): no tree in file!\n");
+ return fCurrentVertex;
+ }
+ AliESD *esdEvent = new AliESD;
+ esdTree->SetBranchAddress("ESD",&esdEvent);
+ esdTree->GetEvent(evnumb);
+ return FindPrimaryVertexForCurrentEvent(esdEvent);
}
//----------------------------------------------------------------------------
-AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
+AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(AliESD *esdEvent)
{
//
// Vertex for current ESD event
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);
+ TTree *trkTree = new TTree("TreeT","tracks");
+ AliESDtrack *esdTrack = 0;
+ trkTree->Branch("tracks","AliESDtrack",&esdTrack);
for(Int_t i=0; i<entr; i++) {
- AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(i);
- if(!esdTrack->GetStatus()&AliESDtrack::kITSin)
- { delete esdTrack; continue; }
- itstrack = new AliITStrackV2(*esdTrack);
+ esdTrack = (AliESDtrack*)esdEvent->GetTrack(i);
+ if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
trkTree->Fill();
- itstrack = 0;
- delete esdTrack;
}
- delete itstrack;
+ delete esdTrack;
// preselect tracks and propagate them to initial vertex position
Int_t nTrks = PrepareTracks(*trkTree);
esdEvent->GetVertex()->GetXYZ(vtx);
SetVtxStart(vtx[0],vtx[1]);
- // VERTEX FINDER
- VertexFinder();
-
// VERTEX FITTER
ComputeMaxChi2PerTrack(nTrks);
VertexFitter();
// store vertex information in ESD
fCurrentVertex->GetXYZ(vtx);
fCurrentVertex->GetCovMatrix(cvtx);
+
+ Double_t tp[3];
+ esdEvent->GetVertex()->GetTruePos(tp);
+ fCurrentVertex->SetTruePos(tp);
+
esdEvent->SetVertex(fCurrentVertex);
cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
return;
}
//---------------------------------------------------------------------------
-void AliITSVertexerTracks::VertexFinder() {
-
- // Get estimate of vertex position in (x,y) from tracks DCA
- // Then this estimate is stored to the data member fInitPos
- // (previous values are overwritten)
-
-
-
- /*
-******* TEMPORARY!!! FOR TEST ONLY!!! **********************************
-
-fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
-fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
- */
-
- fInitPos[2] = 0.;
- for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
-
- Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
-
- Double_t aver[3]={0.,0.,0.};
- Int_t ncombi = 0;
- AliITStrackV2 *track1;
- AliITStrackV2 *track2;
- for(Int_t i=0; i<nacc; i++){
- track1 = (AliITStrackV2*)fTrkArray.At(i);
- if(fDebug>5){
- Double_t xv,par[5];
- track1->GetExternalParameters(xv,par);
- cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
- for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
- cout<<endl;
- }
- Double_t mom1[3];
- Double_t alpha = track1->GetAlpha();
- Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
- Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
- mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
- mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
- mom1[2] = TMath::Cos(theta);
-
- 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);
- for(Int_t j=i+1; j<nacc; j++){
- track2 = (AliITStrackV2*)fTrkArray.At(j);
- Double_t mom2[3];
- alpha = track2->GetAlpha();
- azim = TMath::ASin(track2->GetSnp())+alpha;
- theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
- mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
- mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
- mom2[2] = TMath::Cos(theta);
- 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);
- Double_t crosspoint[3];
- 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();
- }
- else {
- ncombi++;
- for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
- if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
- 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;
- }
- else {
- Warning("VertexFinder","Finder did not succed");
- }
-
-
- //************************************************************************
- return;
-
-}
-//---------------------------------------------------------------------------
void AliITSVertexerTracks::VertexFitter() {
//
// The optimal estimate of the vertex position is given by a "weighted
printf(" VertexFitter(): start\n");
PrintStatus();
}
+ AliVertexerTracks *vertexer=new AliVertexerTracks(fNominalPos[0],fNominalPos[1]);
+ vertexer->SetFinderAlgorithm(5);
+ AliVertex *thevert=(AliVertex*)vertexer->VertexForSelectedTracks(&fTrkArray);
+ Double_t initPos[3];
+ thevert->GetXYZ(initPos);
+ delete vertexer;
Int_t i,j,k,step=0;
TMatrixD rv(3,1);
TMatrixD vV(3,3);
- rv(0,0) = fInitPos[0];
- rv(1,0) = fInitPos[1];
+ rv(0,0) = initPos[0];
+ rv(1,0) = initPos[1];
rv(2,0) = 0.;
Double_t xlStart,alpha;
Double_t rotAngle;
Int_t nUsedTrks;
Double_t chi2,chi2i;
Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
- AliITStrackV2 *t = 0;
+ AliESDtrack *t = 0;
Int_t failed = 0;
Int_t *skipTrack = new Int_t[arrEntries];
for(k=0; k<arrEntries; k++) {
if(skipTrack[k]) continue;
// get track from track array
- t = (AliITStrackV2*)fTrkArray.At(k);
+ t = (AliESDtrack*)fTrkArray.At(k);
alpha = t->GetAlpha();
- xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
- t->PropagateTo(xlStart,0.,0.); // to vtxSeed
+ xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
+ t->PropagateTo(xlStart,AliTracker::GetBz()); // to vtxSeed
rotAngle = alpha;
if(alpha<0.) rotAngle += 2.*TMath::Pi();
cosRot = TMath::Cos(rotAngle);
if(fDebug) printf(" tracks prepared: %d\n",nTrks);
if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
- // VERTEX FINDER
- VertexFinder();
-
// VERTEX FITTER
ComputeMaxChi2PerTrack(nTrks);
VertexFitter();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff