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
19 // Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it
20 // M.Masera, Torino, massimo.masera@to.infn.it
21 //-----------------------------------------------------------------
23 //---- standard headers ----
24 #include <Riostream.h>
25 //---- Root headers --------
32 //---- AliRoot headers -----
34 #include "AliKalmanTrack.h"
35 #include "AliITSStrLine.h"
36 #include "AliITStrackV2.h"
37 #include "AliESDVertex.h"
38 #include "AliITSVertexerTracks.h"
40 #include "AliESDtrack.h"
43 ClassImp(AliITSVertexerTracks)
46 //----------------------------------------------------------------------------
47 AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
49 // Default constructor
57 for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
59 //----------------------------------------------------------------------------
60 AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
63 Double_t xStart,Double_t yStart) {
65 // Standard constructor
73 SetVtxStart(xStart,yStart);
77 for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
80 //----------------------------------------------------------------------------
81 AliITSVertexerTracks::AliITSVertexerTracks(Double_t field, TString fn,
82 Double_t xStart,Double_t yStart)
85 // Alternative constructor
90 SetVtxStart(xStart,yStart);
94 for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
96 //-----------------------------------------------------------------------------
97 AliITSVertexerTracks::~AliITSVertexerTracks() {
99 // The objects poited by the following pointers are not owned
100 // by this class and are not deleted
106 //-----------------------------------------------------------------------------
107 Bool_t AliITSVertexerTracks::CheckField() const {
109 // Check if the conv. const. has been set
112 Double_t cc = t.GetConvConst();
113 Double_t field = 100./0.299792458/cc;
115 if(field<0.1 || field>0.6) {
116 printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
119 printf("AliITSVertexerTracks::CheckField(): Using B = %3.1f T\n",field);
122 //---------------------------------------------------------------------------
123 void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) {
125 // Max. contr. to the chi2 has been tuned as a function of multiplicity
127 if(nTracks < 7) { fMaxChi2PerTrack = 1.e6;
128 } else { fMaxChi2PerTrack = 100.; }
132 //---------------------------------------------------------------------------
133 void AliITSVertexerTracks::FindVertices() {
135 // Vertices for all events from fFirstEvent to fLastEvent
138 // Check if the conv. const. has been set
139 if(!CheckField()) return;
141 TDirectory *curdir = 0;
144 for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
145 if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
147 FindVertexForCurrentEvent(ev);
149 if(!fCurrentVertex) {
150 printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
154 if(fDebug) fCurrentVertex->PrintStatus();
156 // write vertex to file
157 TString vtxName = "Vertex_";
159 fCurrentVertex->SetName(vtxName.Data());
160 fCurrentVertex->SetTitle("VertexerTracks");
161 //WriteCurrentVertex();
164 fCurrentVertex->Write();
167 } // loop over events
171 //---------------------------------------------------------------------------
172 void AliITSVertexerTracks::FindVerticesESD() {
174 // Vertices for all events from fFirstEvent to fLastEvent
177 // Check if the conv. const. has been set
178 if(!CheckField()) return;
180 TDirectory *curdir = 0;
184 TIter next(fInFile->GetListOfKeys());
185 // loop on events in file
186 while ((key=(TKey*)next())!=0) {
187 AliESD *esdEvent=(AliESD*)key->ReadObj();
189 printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n");
192 Int_t ev = (Int_t)esdEvent->GetEventNumber();
193 if(ev<fFirstEvent || ev>fLastEvent) { delete esdEvent; continue; }
194 if(ev % 100 == 0 || fDebug)
195 printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);
197 FindVertexForCurrentEvent(esdEvent);
199 if(!fCurrentVertex) {
200 printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
204 cout<<"VERTICE TROVATO\n";
205 if(fDebug) fCurrentVertex->PrintStatus();
207 // write the ESD to file
210 sprintf(ename,"%d",ev);
213 esdEvent->Write(ename,TObject::kOverwrite);
216 } // loop over events
220 //----------------------------------------------------------------------------
221 Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
223 // Propagate tracks to initial vertex position and store them in a TObjArray
225 Double_t maxd0rphi = 3.;
226 Double_t alpha,xlStart,d0rphi;
230 Int_t nEntries = (Int_t)trkTree.GetEntries();
232 if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
233 fTrkArray.Expand(nEntries);
236 printf(" PrepareTracks()\n");
240 for(Int_t i=0; i<nEntries; i++) {
241 // check tracks to skip
243 for(Int_t j=0; j<fNTrksToSkip; j++) {
244 if(i==fTrksToSkip[j]) {
245 if(fDebug) printf("skipping track: %d\n",i);
249 if(skipThis) continue;
251 AliITStrackV2 *itstrack = new AliITStrackV2;
252 trkTree.SetBranchAddress("tracks",&itstrack);
256 // propagate track to vtxSeed
257 alpha = itstrack->GetAlpha();
258 xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
259 itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
260 itstrack->PropagateTo(xlStart,0.,0.); // to vtxSeed
262 // select tracks with d0rphi < maxd0rphi
263 d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
264 if(d0rphi > maxd0rphi) { delete itstrack; continue; }
266 fTrkArray.AddLast(itstrack);
271 if(fTrksToSkip) delete [] fTrksToSkip;
275 //----------------------------------------------------------------------------
276 void AliITSVertexerTracks::PrintStatus() const {
280 printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
281 printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
282 printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
283 printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
287 //----------------------------------------------------------------------------
288 AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
290 // Vertex for current event
294 // get tree with tracks from input file
295 TString treeName = "TreeT_ITS_";
297 TTree *trkTree=(TTree*)fInFile->Get(treeName.Data());
298 if(!trkTree) return fCurrentVertex;
301 // get tracks and propagate them to initial vertex position
302 Int_t nTrks = PrepareTracks(*trkTree);
304 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
305 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
311 ComputeMaxChi2PerTrack(nTrks);
313 if(fDebug) printf(" vertex fit completed\n");
315 return fCurrentVertex;
317 //----------------------------------------------------------------------------
318 AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
321 // Vertex for current ESD event
324 Double_t vtx[3],cvtx[6];
326 // put tracks reco in ITS in a tree
327 Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
328 TTree *trkTree = new TTree("TreeT_ITS","its tracks");
329 AliITStrackV2 *itstrack = 0;
330 trkTree->Branch("tracks","AliITStrackV2",&itstrack,entr,0);
332 for(Int_t i=0; i<entr; i++) {
333 AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(i);
334 if(!esdTrack->GetStatus()&AliESDtrack::kITSin)
335 { delete esdTrack; continue; }
336 itstrack = new AliITStrackV2(*esdTrack);
343 // preselect tracks and propagate them to initial vertex position
344 Int_t nTrks = PrepareTracks(*trkTree);
346 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
347 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
349 // Set initial vertex position from ESD
350 esdEvent->GetVertex()->GetXYZ(vtx);
351 SetVtxStart(vtx[0],vtx[1]);
357 ComputeMaxChi2PerTrack(nTrks);
359 if(fDebug) printf(" vertex fit completed\n");
361 // store vertex information in ESD
362 fCurrentVertex->GetXYZ(vtx);
363 fCurrentVertex->GetCovMatrix(cvtx);
364 esdEvent->SetVertex(fCurrentVertex);
366 cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
367 return fCurrentVertex;
369 //---------------------------------------------------------------------------
370 void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
372 // Mark the tracks not ot be used in the vertex finding
375 fTrksToSkip = new Int_t[n];
376 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
379 //---------------------------------------------------------------------------
380 void AliITSVertexerTracks::TooFewTracks() {
382 // When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
383 // and the number of tracks to -1
385 fCurrentVertex = new AliESDVertex(0.,0.,-1);
388 //---------------------------------------------------------------------------
389 void AliITSVertexerTracks::VertexFinder() {
391 // Get estimate of vertex position in (x,y) from tracks DCA
392 // Then this estimate is stored to the data member fInitPos
393 // (previous values are overwritten)
398 ******* TEMPORARY!!! FOR TEST ONLY!!! **********************************
400 fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
401 fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
405 for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
407 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
409 Double_t aver[3]={0.,0.,0.};
411 AliITStrackV2 *track1;
412 AliITStrackV2 *track2;
413 for(Int_t i=0; i<nacc; i++){
414 track1 = (AliITStrackV2*)fTrkArray.At(i);
417 track1->GetExternalParameters(xv,par);
418 cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
419 for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
423 Double_t alpha = track1->GetAlpha();
424 Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
425 Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
426 mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
427 mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
428 mom1[2] = TMath::Cos(theta);
431 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
432 track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
433 AliITSStrLine *line1 = new AliITSStrLine(pos1,mom1);
434 for(Int_t j=i+1; j<nacc; j++){
435 track2 = (AliITStrackV2*)fTrkArray.At(j);
437 alpha = track2->GetAlpha();
438 azim = TMath::ASin(track2->GetSnp())+alpha;
439 theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
440 mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
441 mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
442 mom2[2] = TMath::Cos(theta);
444 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
445 track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
446 AliITSStrLine *line2 = new AliITSStrLine(pos2,mom2);
447 Double_t crosspoint[3];
448 Int_t retcode = line2->Cross(line1,crosspoint);
450 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
451 if(fDebug>10)cout<<"bad intersection\n";
452 line1->PrintStatus();
453 line2->PrintStatus();
457 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
458 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
459 if(fDebug>10)cout<<"\n Cross point: ";
460 if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
467 for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
470 Warning("VertexFinder","Finder did not succed");
474 //************************************************************************
478 //---------------------------------------------------------------------------
479 void AliITSVertexerTracks::VertexFitter() {
481 // The optimal estimate of the vertex position is given by a "weighted
482 // average of tracks positions"
483 // Original method: CMS Note 97/0051
486 printf(" VertexFitter(): start\n");
494 rv(0,0) = fInitPos[0];
495 rv(1,0) = fInitPos[1];
497 Double_t xlStart,alpha;
499 Double_t cosRot,sinRot;
503 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
504 AliITStrackV2 *t = 0;
507 Int_t *skipTrack = new Int_t[arrEntries];
508 for(i=0; i<arrEntries; i++) skipTrack[i]=0;
511 // 1st - first estimate of vtx using all tracks
512 // 2nd - apply cut on chi2 max per track
513 // 3rd - estimate of global chi2
514 for(step=0; step<3; step++) {
515 if(fDebug) printf(" step = %d\n",step);
519 TMatrixD SumWiri(3,1);
523 for(j=0; j<3; j++) SumWi(j,i) = 0.;
527 for(k=0; k<arrEntries; k++) {
528 if(skipTrack[k]) continue;
529 // get track from track array
530 t = (AliITStrackV2*)fTrkArray.At(k);
531 alpha = t->GetAlpha();
532 xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
533 t->PropagateTo(xlStart,0.,0.); // to vtxSeed
535 if(alpha<0.) rotAngle += 2.*TMath::Pi();
536 cosRot = TMath::Cos(rotAngle);
537 sinRot = TMath::Sin(rotAngle);
539 // vector of track global coordinates
541 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
542 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
545 // matrix to go from global (x,y,z) to local (y,z);
554 // covariance matrix of local (y,z) - inverted
556 t->GetExternalCovariance(cc);
562 // weights matrix: Wi = QiT * UiInv * Qi
563 if(Ui.Determinant() <= 0.) continue;
564 TMatrixD UiInv(TMatrixD::kInverted,Ui);
565 TMatrixD UiInvQi(UiInv,TMatrixD::kMult,Qi);
566 TMatrixD Wi(Qi,TMatrixD::kTransposeMult,UiInvQi);
569 TMatrixD deltar = rv; deltar -= ri;
570 TMatrixD Wideltar(Wi,TMatrixD::kMult,deltar);
571 chi2i = deltar(0,0)*Wideltar(0,0)+
572 deltar(1,0)*Wideltar(1,0)+
573 deltar(2,0)*Wideltar(2,0);
576 if(step==1 && chi2i > fMaxChi2PerTrack) {
584 TMatrixD Wiri(Wi,TMatrixD::kMult,ri);
590 } // end loop on tracks
592 if(nUsedTrks < fMinTracks) {
597 Double_t determinant = SumWi.Determinant();
598 //cerr<<" determinant: "<<determinant<<endl;
599 if(determinant < 100.) {
600 printf("det(V) = 0\n");
605 // inverted of weights matrix
606 TMatrixD InvSumWi(TMatrixD::kInverted,SumWi);
609 // position of primary vertex
612 } // end loop on the 3 steps
622 Double_t position[3];
623 position[0] = rv(0,0);
624 position[1] = rv(1,0);
625 position[2] = rv(2,0);
626 Double_t covmatrix[6];
627 covmatrix[0] = V(0,0);
628 covmatrix[1] = V(0,1);
629 covmatrix[2] = V(1,1);
630 covmatrix[3] = V(0,2);
631 covmatrix[4] = V(1,2);
632 covmatrix[5] = V(2,2);
634 // store data in the vertex object
635 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
638 printf(" VertexFitter(): finish\n");
639 printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
640 fCurrentVertex->PrintStatus();
645 //----------------------------------------------------------------------------
646 AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
648 // Return vertex from tracks in trkTree
650 if(fCurrentVertex) fCurrentVertex = 0;
652 // get tracks and propagate them to initial vertex position
653 Int_t nTrks = PrepareTracks(*(&trkTree));
654 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
655 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
661 ComputeMaxChi2PerTrack(nTrks);
663 if(fDebug) printf(" vertex fit completed\n");
665 return fCurrentVertex;
667 //----------------------------------------------------------------------------