1 /**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //-----------------------------------------------------------------
17 // Implementation of the vertexer from tracks
18 // It accepts V2 and ESD tracks
20 // Origin: A.Dainese, Padova,
21 // andrea.dainese@pd.infn.it
23 // massimo.masera@to.infn.it
24 //-----------------------------------------------------------------
26 //---- standard headers ----
27 #include <Riostream.h>
28 //---- Root headers --------
32 //---- AliRoot headers -----
33 #include "AliESDVertex.h"
34 #include "AliITSVertexerTracks.h"
36 #include "AliESDtrack.h"
37 #include "AliVertexerTracks.h"
40 ClassImp(AliITSVertexerTracks)
43 //----------------------------------------------------------------------------
44 AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
46 // Default constructor
56 //----------------------------------------------------------------------------
57 AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
59 Double_t xStart,Double_t yStart) {
61 // Standard constructor
68 SetVtxStart(xStart,yStart);
75 //----------------------------------------------------------------------------
76 AliITSVertexerTracks::AliITSVertexerTracks(TString fn,
77 Double_t xStart,Double_t yStart)
80 // Alternative constructor
84 SetVtxStart(xStart,yStart);
90 //______________________________________________________________________
91 AliITSVertexerTracks::AliITSVertexerTracks(const AliITSVertexerTracks &vtxr) : AliITSVertexer(vtxr) {
93 // Copies are not allowed. The method is protected to avoid misuse.
94 Error("AliITSVertexerTracks","Copy constructor not allowed\n");
97 //______________________________________________________________________
98 AliITSVertexerTracks& AliITSVertexerTracks::operator=(const AliITSVertexerTracks& /* vtxr */){
99 // Assignment operator
100 // Assignment is not allowed. The method is protected to avoid misuse.
101 Error("= operator","Assignment operator not allowed\n");
105 //-----------------------------------------------------------------------------
106 AliITSVertexerTracks::~AliITSVertexerTracks() {
107 // Default Destructor
108 // The objects poited by the following pointers are not owned
109 // by this class and are not deleted
115 //-----------------------------------------------------------------------------
116 Bool_t AliITSVertexerTracks::CheckField() const {
118 // Check if the conv. const. has been set
120 Double_t field = AliTracker::GetBz(); // in kG
122 if(field<1 || field>6) {
123 printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
126 printf("AliITSVertexerTracks::CheckField(): Using B = %3.1f T\n",field);
129 //---------------------------------------------------------------------------
130 void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) {
132 // Max. contr. to the chi2 has been tuned as a function of multiplicity
134 if(nTracks < 7) { fMaxChi2PerTrack = 1.e6;
135 } else { fMaxChi2PerTrack = 100.; }
139 //---------------------------------------------------------------------------
140 void AliITSVertexerTracks::FindVertices() {
142 // Vertices for all events from fFirstEvent to fLastEvent
145 // Check if the conv. const. has been set
146 if(!CheckField()) return;
148 TDirectory *curdir = 0;
151 for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
152 if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
154 FindPrimaryVertexForCurrentEvent(ev);
156 if(!fCurrentVertex) {
157 printf("AliITSVertexerTracks::FindVertices(): no tracks tree for event %d\n",ev);
161 if(fDebug) fCurrentVertex->PrintStatus();
163 // write vertex to file
164 TString vtxName = "Vertex_";
166 fCurrentVertex->SetName(vtxName.Data());
167 fCurrentVertex->SetTitle("VertexerTracks");
168 //WriteCurrentVertex();
171 fCurrentVertex->Write();
174 } // loop over events
178 //---------------------------------------------------------------------------
179 void AliITSVertexerTracks::FindVerticesESD() {
181 // Vertices for all events from fFirstEvent to fLastEvent
184 // Check if the conv. const. has been set
185 if(!CheckField()) return;
187 TDirectory *curdir = 0;
190 TTree *esdTree = (TTree*)fInFile->Get("esdTree");
192 printf("AliITSVertexerTracks::FindVerticesESD(): no tree in file!\n");
195 Int_t nev = (Int_t)esdTree->GetEntries();
197 // loop on events in tree
198 for(Int_t i=0; i<nev; i++) {
199 AliESD *esdEvent = new AliESD;
200 esdTree->SetBranchAddress("ESD",&esdEvent);
201 if(!esdTree->GetEvent(i)) {
202 printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n");
206 ev = (Int_t)esdEvent->GetEventNumber();
207 if(ev<fFirstEvent || ev>fLastEvent) { delete esdEvent; continue; }
208 if(ev % 100 == 0 || fDebug)
209 printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);
211 FindPrimaryVertexForCurrentEvent(esdEvent);
213 if(!fCurrentVertex) {
214 printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
218 if(fDebug) fCurrentVertex->PrintStatus();
220 // write vertex to file
221 TString vtxName = "Vertex_";
223 fCurrentVertex->SetName(vtxName.Data());
224 fCurrentVertex->SetTitle("VertexerTracks");
225 //WriteCurrentVertex();
228 fCurrentVertex->Write();
233 } // end loop over events
237 //----------------------------------------------------------------------------
238 Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
240 // Propagate tracks to initial vertex position and store them in a TObjArray
242 Double_t maxd0rphi = 3.;
243 Double_t alpha,xlStart,d0rphi;
246 Int_t nEntries = (Int_t)trkTree.GetEntries();
248 Double_t field=AliTracker::GetBz();
250 if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
251 fTrkArray.Expand(nEntries);
254 printf(" PrepareTracks()\n");
258 for(Int_t i=0; i<nEntries; i++) {
259 // check tracks to skip
261 for(Int_t j=0; j<fNTrksToSkip; j++) {
262 if(i==fTrksToSkip[j]) {
263 if(fDebug) printf("skipping track: %d\n",i);
267 if(skipThis) continue;
269 AliESDtrack *track = new AliESDtrack;
270 trkTree.SetBranchAddress("tracks",&track);
273 // propagate track to vtxSeed
274 alpha = track->GetAlpha();
275 xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
276 track->PropagateTo(xlStart,field); // to vtxSeed
278 // select tracks with d0rphi < maxd0rphi
279 d0rphi = TMath::Abs(track->GetD(fNominalPos[0],fNominalPos[1],field));
280 if(d0rphi > maxd0rphi) { delete track; continue; }
282 fTrkArray.AddLast(track);
285 if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<" "<<d0rphi<<endl;
288 if(fTrksToSkip) delete [] fTrksToSkip;
292 //----------------------------------------------------------------------------
293 void AliITSVertexerTracks::PrintStatus() const {
296 // printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
297 printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
298 printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
302 //----------------------------------------------------------------------------
303 AliVertex* AliITSVertexerTracks::VertexForSelectedTracks(AliESD *esdEvent,Int_t nofCand, Int_t *trkPos, Int_t opt){
306 // Computes the vertex for selected tracks
307 // trkPos=vector with track positions in ESD
310 esdEvent->GetVertex()->GetXYZ(vtx);
311 TTree *trkTree = new TTree("TreeT","tracks");
312 AliESDtrack *esdTrack = 0;
313 trkTree->Branch("tracks","AliESDtrack",&esdTrack);
314 for(Int_t i=0; i<nofCand;i++){
315 esdTrack = (AliESDtrack*)esdEvent->GetTrack(trkPos[i]);
317 if(!esdTrack->GetStatus()&AliESDtrack::kTPCin) continue;
318 if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
319 if(!esdTrack->GetStatus()&AliESDtrack::kITSrefit) continue;
321 Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
322 if(nclus<6) continue;
326 Int_t nTrks = PrepareTracks(*trkTree);
327 //delete trkTree;// :-))
328 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
329 if(nTrks < fMinTracks) {
330 printf("TooFewTracks\n");
331 AliVertex *theVert=new AliVertex();
332 theVert->SetDispersion(999);
333 theVert->SetNContributors(-5);
337 AliVertexerTracks *vertexer=new AliVertexerTracks(vtx[0],vtx[1]);
338 vertexer->SetFinderAlgorithm(opt);
339 AliVertex *theVert=(AliVertex*)vertexer->VertexForSelectedTracks(&fTrkArray);
340 // beware: newvt object should be deleted by the caller
341 AliVertex *newvt = new AliVertex(*theVert);
345 //----------------------------------------------------------------------------
346 AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(Int_t evnumb) {
348 // Vertex for current event
352 // get tree with tracks from input file
354 TTree *esdTree = (TTree*)fInFile->Get("esdTree");
357 printf("AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(): no tree in file!\n");
358 return fCurrentVertex;
360 AliESD *esdEvent = new AliESD;
361 esdTree->SetBranchAddress("ESD",&esdEvent);
362 esdTree->GetEvent(evnumb);
363 return FindPrimaryVertexForCurrentEvent(esdEvent);
365 //----------------------------------------------------------------------------
366 AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(AliESD *esdEvent)
369 // Vertex for current ESD event
372 Double_t vtx[3],cvtx[6];
374 Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
375 TTree *trkTree = new TTree("TreeT","tracks");
376 AliESDtrack *esdTrack = 0;
377 trkTree->Branch("tracks","AliESDtrack",&esdTrack);
379 for(Int_t i=0; i<entr; i++) {
380 esdTrack = (AliESDtrack*)esdEvent->GetTrack(i);
381 if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
386 // preselect tracks and propagate them to initial vertex position
387 Int_t nTrks = PrepareTracks(*trkTree);
389 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
390 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
392 // Set initial vertex position from ESD
393 esdEvent->GetVertex()->GetXYZ(vtx);
394 SetVtxStart(vtx[0],vtx[1]);
397 ComputeMaxChi2PerTrack(nTrks);
399 if(fDebug) printf(" vertex fit completed\n");
401 // store vertex information in ESD
402 fCurrentVertex->GetXYZ(vtx);
403 fCurrentVertex->GetCovMatrix(cvtx);
406 esdEvent->GetVertex()->GetTruePos(tp);
407 fCurrentVertex->SetTruePos(tp);
409 esdEvent->SetVertex(fCurrentVertex);
411 cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
412 return fCurrentVertex;
414 //---------------------------------------------------------------------------
415 void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
417 // Mark the tracks not ot be used in the vertex finding
420 fTrksToSkip = new Int_t[n];
421 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
424 //---------------------------------------------------------------------------
425 void AliITSVertexerTracks::TooFewTracks() {
427 // When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
428 // and the number of tracks to -1
430 fCurrentVertex = new AliESDVertex(0.,0.,-1);
433 //---------------------------------------------------------------------------
434 void AliITSVertexerTracks::VertexFitter() {
436 // The optimal estimate of the vertex position is given by a "weighted
437 // average of tracks positions"
438 // Original method: CMS Note 97/0051
441 printf(" VertexFitter(): start\n");
444 AliVertexerTracks *vertexer=new AliVertexerTracks(fNominalPos[0],fNominalPos[1]);
445 vertexer->SetFinderAlgorithm(5);
446 AliVertex *thevert=(AliVertex*)vertexer->VertexForSelectedTracks(&fTrkArray);
448 thevert->GetXYZ(initPos);
455 rv(0,0) = initPos[0];
456 rv(1,0) = initPos[1];
458 Double_t xlStart,alpha;
460 Double_t cosRot,sinRot;
464 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
468 Int_t *skipTrack = new Int_t[arrEntries];
469 for(i=0; i<arrEntries; i++) skipTrack[i]=0;
472 // 1st - first estimate of vtx using all tracks
473 // 2nd - apply cut on chi2 max per track
474 // 3rd - estimate of global chi2
475 for(step=0; step<3; step++) {
476 if(fDebug) printf(" step = %d\n",step);
480 TMatrixD sumWiri(3,1);
484 for(j=0; j<3; j++) sumWi(j,i) = 0.;
488 for(k=0; k<arrEntries; k++) {
489 if(skipTrack[k]) continue;
490 // get track from track array
491 t = (AliESDtrack*)fTrkArray.At(k);
492 alpha = t->GetAlpha();
493 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
494 t->PropagateTo(xlStart,AliTracker::GetBz()); // to vtxSeed
496 if(alpha<0.) rotAngle += 2.*TMath::Pi();
497 cosRot = TMath::Cos(rotAngle);
498 sinRot = TMath::Sin(rotAngle);
500 // vector of track global coordinates
502 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
503 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
506 // matrix to go from global (x,y,z) to local (y,z);
515 // covariance matrix of local (y,z) - inverted
517 t->GetExternalCovariance(cc);
523 // weights matrix: wWi = qQiT * uUiInv * qQi
524 if(uUi.Determinant() <= 0.) continue;
525 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
526 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
527 TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
530 TMatrixD deltar = rv; deltar -= ri;
531 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
532 chi2i = deltar(0,0)*wWideltar(0,0)+
533 deltar(1,0)*wWideltar(1,0)+
534 deltar(2,0)*wWideltar(2,0);
537 if(step==1 && chi2i > fMaxChi2PerTrack) {
545 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
551 } // end loop on tracks
553 if(nUsedTrks < fMinTracks) {
558 Double_t determinant = sumWi.Determinant();
559 //cerr<<" determinant: "<<determinant<<endl;
560 if(determinant < 100.) {
561 printf("det(V) = 0\n");
566 // inverted of weights matrix
567 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
570 // position of primary vertex
573 } // end loop on the 3 steps
583 Double_t position[3];
584 position[0] = rv(0,0);
585 position[1] = rv(1,0);
586 position[2] = rv(2,0);
587 Double_t covmatrix[6];
588 covmatrix[0] = vV(0,0);
589 covmatrix[1] = vV(0,1);
590 covmatrix[2] = vV(1,1);
591 covmatrix[3] = vV(0,2);
592 covmatrix[4] = vV(1,2);
593 covmatrix[5] = vV(2,2);
595 // store data in the vertex object
596 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
599 printf(" VertexFitter(): finish\n");
600 printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
601 fCurrentVertex->PrintStatus();
606 //----------------------------------------------------------------------------
607 AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
609 // Return vertex from tracks in trkTree
611 if(fCurrentVertex) fCurrentVertex = 0;
613 // get tracks and propagate them to initial vertex position
614 Int_t nTrks = PrepareTracks(*(&trkTree));
615 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
616 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
619 ComputeMaxChi2PerTrack(nTrks);
621 if(fDebug) printf(" vertex fit completed\n");
623 return fCurrentVertex;
625 //----------------------------------------------------------------------------