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 --------
30 #include <TObjArray.h>
31 #include <TRandom.h> // TEMPORARY !!!!!!!
32 //---- AliRoot headers -----
34 #include "AliKalmanTrack.h"
35 #include "AliITSStrLine.h"
36 #include "AliITStrackV2.h"
37 #include "AliITSVertex.h"
38 #include "AliITSVertexerTracks.h"
41 ClassImp(AliITSVertexerTracks)
44 //----------------------------------------------------------------------------
45 AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
47 // Default constructor
56 for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
59 //----------------------------------------------------------------------------
60 AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
62 Double_t xStart,Double_t yStart,
64 :AliITSVertexer(inFile,outFile) {
66 // Standard constructor
70 SetVtxStart(xStart,yStart);
72 SetUseThrustFrame(useThFr);
76 for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
78 //----------------------------------------------------------------------------
79 Bool_t AliITSVertexerTracks::CheckField() const {
81 // Check if the conv. const. has been set
84 Double_t cc = t.GetConvConst();
85 Double_t field = 100./0.299792458/cc;
87 if(field<0.1 || field>0.6) {
88 printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
91 printf("AliITSVertexerTracks::CheckField(): Using B = %3.1f T\n",field);
94 //---------------------------------------------------------------------------
95 void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) {
97 // Max. contr. to the chi2 has been tuned as a function of multiplicity
99 if(nTracks < 7) { fMaxChi2PerTrack = 1.e6;
100 } else { fMaxChi2PerTrack = 100.; }
104 //---------------------------------------------------------------------------
105 void AliITSVertexerTracks::FindVertices() {
107 // Vertices for all events from fFirstEvent to fLastEvent
110 // Check if the conv. const. has been set
111 if(!CheckField()) return;
115 for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
116 if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
118 FindVertexForCurrentEvent(ev);
120 if(!fCurrentVertex) {
121 printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
125 if(fDebug) fCurrentVertex->PrintStatus();
127 WriteCurrentVertex();
128 } // loop over events
132 //----------------------------------------------------------------------------
133 Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
135 // Propagate tracks to initial vertex position and store them in a TObjArray
137 Double_t maxd0rphi = 3.;
138 Double_t alpha,xlStart,d0rphi;
142 Int_t nEntries = (Int_t)trkTree.GetEntries();
144 if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
145 fTrkArray.Expand(nEntries);
148 printf(" PrepareTracks()\n");
152 for(Int_t i=0; i<nEntries; i++) {
153 // check tracks to skip
155 for(Int_t j=0; j<fNTrksToSkip; j++) {
156 if(i==fTrksToSkip[j]) {
157 if(fDebug) printf("skipping track: %d\n",i);
161 if(skipThis) continue;
163 AliITStrackV2 *itstrack = new AliITStrackV2;
164 trkTree.SetBranchAddress("tracks",&itstrack);
168 // propagate track to vtxSeed
169 alpha = itstrack->GetAlpha();
170 xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
171 itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
172 itstrack->PropagateTo(xlStart,0.,0.); // to vtxSeed
174 // select tracks with d0rphi < maxd0rphi
175 d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
176 if(d0rphi > maxd0rphi) { delete itstrack; continue; }
178 fTrkArray.AddLast(itstrack);
183 if(fTrksToSkip) delete [] fTrksToSkip;
187 //----------------------------------------------------------------------------
188 void AliITSVertexerTracks::PrintStatus() const {
192 printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
193 printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
194 printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
195 printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
196 printf(" Using Thrust Frame: %d fPhiThrust = %f\n",fUseThrustFrame,fPhiThrust);
200 //----------------------------------------------------------------------------
201 AliITSVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
203 // Vertex for current event
207 // get tree with tracks from input file
208 TString treeName = "TreeT_ITS_";
210 TTree *trkTree=(TTree*)fInFile->Get(treeName.Data());
211 if(!trkTree) return fCurrentVertex;
214 // get tracks and propagate them to initial vertex position
215 Int_t nTrks = PrepareTracks(*trkTree);
217 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
218 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
224 ComputeMaxChi2PerTrack(nTrks);
225 if(fUseThrustFrame) ThrustFinderXY();
226 if(fDebug) printf(" thrust found: phi = %f\n",fPhiThrust);
228 if(fDebug) printf(" vertex fit completed\n");
231 vtxName = "VertexTracks_";
233 fCurrentVertex->SetName(vtxName.Data());
234 return fCurrentVertex;
236 //---------------------------------------------------------------------------
237 void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
239 // Mark the tracks not ot be used in the vertex finding
242 fTrksToSkip = new Int_t[n];
243 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
246 //----------------------------------------------------------------------------
247 Double_t AliITSVertexerTracks::SumPl(TTree &momTree,Double_t phi) const {
249 // Function to be maximized for thrust determination
259 momTree.SetBranchAddress("momenta",&mom);
260 Int_t entries = (Int_t)momTree.GetEntries();
262 for(Int_t i=0; i<entries; i++) {
272 //---------------------------------------------------------------------------
273 void AliITSVertexerTracks::ThrustFinderXY() {
275 // This function looks for the thrust direction, \vec{u}, in the (x,y) plane.
276 // The following function is maximized:
277 // \Sum_{\vec{p}\cdot\vec{u}} \vec{p}\cdot\vec{u} / \Sum |\vec{p}|
278 // where \vec{p} = (p_x,p_y)
280 Double_t pt,alpha,phi;
283 // tree for thrust determination
284 TVector3 *ioMom = new TVector3;
285 TTree *t = new TTree("Tree_Momenta","Tree with momenta");
286 t->Branch("momenta","TVector3",&ioMom);
289 AliITStrackV2 *itstrack = 0;
290 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
293 for(Int_t i=0; i<arrEntries; i++) {
294 itstrack = (AliITStrackV2*)fTrkArray.At(i);
295 // momentum of the track at the vertex
296 pt = 1./TMath::Abs(itstrack->Get1Pt());
297 alpha = itstrack->GetAlpha();
298 phi = alpha+TMath::ASin(itstrack->GetSnp());
299 ioMom->SetX(pt*TMath::Cos(phi));
300 ioMom->SetY(pt*TMath::Sin(phi));
303 totPt += ioMom->Pt();
305 } // end loop on tracks
307 Double_t tValue=0.,tPhi=0.;
308 Double_t maxSumPl = 0.;
315 dPhi = 2.*TMath::Pi()/(Double_t)nSteps;
317 for(iStep=0; iStep<nSteps; iStep++) {
319 thisSumPl = SumPl(*t,phi);
320 if(thisSumPl > maxSumPl) {
321 maxSumPl = thisSumPl;
330 for(iStep=0; iStep<nSteps; iStep++) {
332 thisSumPl = SumPl(*t,phi);
333 if(thisSumPl > maxSumPl) {
334 maxSumPl = thisSumPl;
339 tValue = 2.*maxSumPl/totPt;
340 if(tPhi<0.) tPhi += 2.*TMath::Pi();
341 if(tPhi>2.*TMath::Pi()) tPhi -= 2.*TMath::Pi();
350 //---------------------------------------------------------------------------
351 void AliITSVertexerTracks::TooFewTracks() {
353 // When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
354 // and the number of tracks to -1
356 fCurrentVertex = new AliITSVertex(0.,0.,-1);
359 //---------------------------------------------------------------------------
360 void AliITSVertexerTracks::VertexFinder() {
362 // Get estimate of vertex position in (x,y) from tracks DCA
363 // Then this estimate is stored to the data member fInitPos
364 // (previous values are overwritten)
369 ******* TEMPORARY!!! FOR TEST ONLY!!! **********************************
371 fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
372 fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
376 for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
378 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
380 Double_t aver[3]={0.,0.,0.};
382 AliITStrackV2 *track1;
383 AliITStrackV2 *track2;
384 for(Int_t i=0; i<nacc; i++){
385 track1 = (AliITStrackV2*)fTrkArray.At(i);
388 track1->GetExternalParameters(xv,par);
389 cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
390 for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
394 Double_t alpha = track1->GetAlpha();
395 Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
396 Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
397 mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
398 mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
399 mom1[2] = TMath::Cos(theta);
402 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
403 track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
404 AliITSStrLine *line1 = new AliITSStrLine(pos1,mom1);
405 for(Int_t j=i+1; j<nacc; j++){
406 track2 = (AliITStrackV2*)fTrkArray.At(j);
408 alpha = track2->GetAlpha();
409 azim = TMath::ASin(track2->GetSnp())+alpha;
410 theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
411 mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
412 mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
413 mom2[2] = TMath::Cos(theta);
415 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
416 track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
417 AliITSStrLine *line2 = new AliITSStrLine(pos2,mom2);
418 Double_t crosspoint[3];
419 Int_t retcode = line2->Cross(line1,crosspoint);
421 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
422 if(fDebug>10)cout<<"bad intersection\n";
423 line1->PrintStatus();
424 line2->PrintStatus();
428 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
429 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
430 if(fDebug>10)cout<<"\n Cross point: ";
431 if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
438 for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
441 Warning("VertexFinder","Finder did not succed");
445 //************************************************************************
450 //---------------------------------------------------------------------------
451 void AliITSVertexerTracks::VertexFitter() {
453 // The optimal estimate of the vertex position is given by a "weighted
454 // average of tracks positions"
455 // Original method: CMS Note 97/0051
458 printf(" VertexFitter(): start\n");
466 rv(0,0) = fInitPos[0];
467 rv(1,0) = fInitPos[1];
469 Double_t xlStart,alpha;
471 Double_t cosRot,sinRot;
475 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
476 AliITStrackV2 *t = 0;
479 Int_t *skipTrack = new Int_t[arrEntries];
480 for(i=0; i<arrEntries; i++) skipTrack[i]=0;
483 // 1st - first estimate of vtx using all tracks
484 // 2nd - apply cut on chi2 max per track
485 // 3rd - estimate of global chi2
486 for(step=0; step<3; step++) {
487 if(fDebug) printf(" step = %d\n",step);
491 TMatrixD SumWiri(3,1);
495 for(j=0; j<3; j++) SumWi(j,i) = 0.;
499 for(k=0; k<arrEntries; k++) {
500 if(skipTrack[k]) continue;
501 // get track from track array
502 t = (AliITStrackV2*)fTrkArray.At(k);
503 alpha = t->GetAlpha();
504 xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
505 t->PropagateTo(xlStart,0.,0.); // to vtxSeed
506 rotAngle = alpha-fPhiThrust;
507 if(alpha<0.) rotAngle += 2.*TMath::Pi();
508 cosRot = TMath::Cos(rotAngle);
509 sinRot = TMath::Sin(rotAngle);
511 // vector of track global coordinates
513 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
514 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
517 // matrix to go from global (x,y,z) to local (y,z);
526 // covariance matrix of local (y,z) - inverted
528 t->GetExternalCovariance(cc);
534 // weights matrix: Wi = QiT * UiInv * Qi
535 if(Ui.Determinant() <= 0.) continue;
536 TMatrixD UiInv(TMatrixD::kInverted,Ui);
537 TMatrixD UiInvQi(UiInv,TMatrixD::kMult,Qi);
538 TMatrixD Wi(Qi,TMatrixD::kTransposeMult,UiInvQi);
541 TMatrixD deltar = rv; deltar -= ri;
542 TMatrixD Wideltar(Wi,TMatrixD::kMult,deltar);
543 chi2i = deltar(0,0)*Wideltar(0,0)+
544 deltar(1,0)*Wideltar(1,0)+
545 deltar(2,0)*Wideltar(2,0);
548 if(step==1 && chi2i > fMaxChi2PerTrack) {
556 TMatrixD Wiri(Wi,TMatrixD::kMult,ri);
562 } // end loop on tracks
564 if(nUsedTrks < fMinTracks) {
569 Double_t determinant = SumWi.Determinant();
570 //cerr<<" determinant: "<<determinant<<endl;
571 if(determinant < 100.) {
572 printf("det(V) = 0\n");
577 // inverted of weights matrix
578 TMatrixD InvSumWi(TMatrixD::kInverted,SumWi);
581 // position of primary vertex
584 } // end loop on the 3 steps
594 Double_t position[3];
595 position[0] = rv(0,0);
596 position[1] = rv(1,0);
597 position[2] = rv(2,0);
598 Double_t covmatrix[6];
599 covmatrix[0] = V(0,0);
600 covmatrix[1] = V(0,1);
601 covmatrix[2] = V(1,1);
602 covmatrix[3] = V(0,2);
603 covmatrix[4] = V(1,2);
604 covmatrix[5] = V(2,2);
606 // store data in the vertex object
607 fCurrentVertex = new AliITSVertex(fPhiThrust,position,covmatrix,chi2,nUsedTrks);
610 printf(" VertexFitter(): finish\n");
611 printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
612 fCurrentVertex->PrintStatus();
617 //----------------------------------------------------------------------------
618 AliITSVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
620 // Return vertex from tracks in trkTree
622 if(fCurrentVertex) fCurrentVertex = 0;
624 // get tracks and propagate them to initial vertex position
625 Int_t nTrks = PrepareTracks(*(&trkTree));
626 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
627 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
633 ComputeMaxChi2PerTrack(nTrks);
634 if(fUseThrustFrame) ThrustFinderXY();
635 if(fDebug) printf(" thrust found: phi = %f\n",fPhiThrust);
637 if(fDebug) printf(" vertex fit completed\n");
639 return fCurrentVertex;
641 //----------------------------------------------------------------------------