]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSVertexerTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Consistent usage of the magnetic field interface (Yu.Belikov)
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerTracks.cxx
CommitLineData
cab6ff9b 1/**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16//-----------------------------------------------------------------
17// Implementation of the vertexer from tracks
a48d1ed3 18// It accepts V2 and ESD tracks
cab6ff9b 19//
a48d1ed3 20// Origin: A.Dainese, Padova,
21// andrea.dainese@pd.infn.it
22// M.Masera, Torino,
23// massimo.masera@to.infn.it
cab6ff9b 24//-----------------------------------------------------------------
25
26//---- standard headers ----
27#include <Riostream.h>
28//---- Root headers --------
11ba84a4 29#include <TKey.h>
cab6ff9b 30#include <TFile.h>
31#include <TTree.h>
cab6ff9b 32#include <TMatrixD.h>
cab6ff9b 33//---- AliRoot headers -----
cab6ff9b 34#include "AliITSStrLine.h"
35#include "AliITStrackV2.h"
d681bb2d 36#include "AliESDVertex.h"
cab6ff9b 37#include "AliITSVertexerTracks.h"
11ba84a4 38#include "AliESD.h"
39#include "AliESDtrack.h"
cab6ff9b 40
41
42ClassImp(AliITSVertexerTracks)
43
44
45//----------------------------------------------------------------------------
46AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
47//
48// Default constructor
49//
11ba84a4 50 fInFile = 0;
51 fOutFile = 0;
cab6ff9b 52 SetVtxStart();
53 SetMinTracks();
cab6ff9b 54 fTrksToSkip = 0;
55 fNTrksToSkip = 0;
66e811e2 56 fDCAcut=0;
cab6ff9b 57}
cab6ff9b 58//----------------------------------------------------------------------------
11ba84a4 59AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
11ba84a4 60 Int_t fEv,Int_t lEv,
61 Double_t xStart,Double_t yStart) {
cab6ff9b 62//
63// Standard constructor
64//
11ba84a4 65 fCurrentVertex = 0;
66 fInFile = inFile;
67 fOutFile = outFile;
68 SetFirstEvent(fEv);
69 SetLastEvent(lEv);
cab6ff9b 70 SetVtxStart(xStart,yStart);
71 SetMinTracks();
cab6ff9b 72 fTrksToSkip = 0;
73 fNTrksToSkip = 0;
66e811e2 74 fDCAcut=0;
11ba84a4 75 SetDebug();
cab6ff9b 76}
77//----------------------------------------------------------------------------
7d0f8548 78AliITSVertexerTracks::AliITSVertexerTracks(TString fn,
11ba84a4 79 Double_t xStart,Double_t yStart)
80 :AliITSVertexer(fn) {
81//
82// Alternative constructor
83//
84 fInFile = 0;
85 fOutFile = 0;
11ba84a4 86 SetVtxStart(xStart,yStart);
87 SetMinTracks();
88 fTrksToSkip = 0;
89 fNTrksToSkip = 0;
66e811e2 90 fDCAcut=0;
11ba84a4 91}
a48d1ed3 92//______________________________________________________________________
93AliITSVertexerTracks::AliITSVertexerTracks(const AliITSVertexerTracks &vtxr) : AliITSVertexer(vtxr) {
94 // Copy constructor
95 // Copies are not allowed. The method is protected to avoid misuse.
96 Error("AliITSVertexerTracks","Copy constructor not allowed\n");
97}
98
99//______________________________________________________________________
100AliITSVertexerTracks& AliITSVertexerTracks::operator=(const AliITSVertexerTracks& /* vtxr */){
101 // Assignment operator
102 // Assignment is not allowed. The method is protected to avoid misuse.
103 Error("= operator","Assignment operator not allowed\n");
104 return *this;
105}
106
11ba84a4 107//-----------------------------------------------------------------------------
108AliITSVertexerTracks::~AliITSVertexerTracks() {
109 // Default Destructor
110 // The objects poited by the following pointers are not owned
111 // by this class and are not deleted
112
113 fCurrentVertex = 0;
114 fInFile = 0;
115 fOutFile = 0;
116}
117//-----------------------------------------------------------------------------
cab6ff9b 118Bool_t AliITSVertexerTracks::CheckField() const {
119//
120// Check if the conv. const. has been set
121//
122 AliITStrackV2 t;
123 Double_t cc = t.GetConvConst();
124 Double_t field = 100./0.299792458/cc;
125
126 if(field<0.1 || field>0.6) {
127 printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
128 return kFALSE;
129 }
130 printf("AliITSVertexerTracks::CheckField(): Using B = %3.1f T\n",field);
131 return kTRUE;
132}
133//---------------------------------------------------------------------------
134void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) {
135//
136// Max. contr. to the chi2 has been tuned as a function of multiplicity
137//
138 if(nTracks < 7) { fMaxChi2PerTrack = 1.e6;
139 } else { fMaxChi2PerTrack = 100.; }
140
141 return;
142}
143//---------------------------------------------------------------------------
144void AliITSVertexerTracks::FindVertices() {
145//
146// Vertices for all events from fFirstEvent to fLastEvent
147//
148
149 // Check if the conv. const. has been set
150 if(!CheckField()) return;
151
11ba84a4 152 TDirectory *curdir = 0;
cab6ff9b 153
154 // loop over events
155 for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
156 if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
157
66e811e2 158 FindPrimaryVertexForCurrentEvent(ev);
cab6ff9b 159
160 if(!fCurrentVertex) {
161 printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
162 continue;
163 }
164
165 if(fDebug) fCurrentVertex->PrintStatus();
11ba84a4 166
167 // write vertex to file
d4221964 168 TString vtxName = "Vertex_";
169 vtxName += ev;
11ba84a4 170 fCurrentVertex->SetName(vtxName.Data());
d4221964 171 fCurrentVertex->SetTitle("VertexerTracks");
11ba84a4 172 //WriteCurrentVertex();
173 curdir = gDirectory;
174 fOutFile->cd();
175 fCurrentVertex->Write();
176 curdir->cd();
177 fCurrentVertex = 0;
178 } // loop over events
179
180 return;
181}
182//---------------------------------------------------------------------------
183void AliITSVertexerTracks::FindVerticesESD() {
184//
185// Vertices for all events from fFirstEvent to fLastEvent
186//
187
188 // Check if the conv. const. has been set
189 if(!CheckField()) return;
190
191 TDirectory *curdir = 0;
192
193 fInFile->cd();
194 TKey *key=0;
195 TIter next(fInFile->GetListOfKeys());
196 // loop on events in file
197 while ((key=(TKey*)next())!=0) {
198 AliESD *esdEvent=(AliESD*)key->ReadObj();
199 if(!esdEvent) {
200 printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n");
201 return;
202 }
203 Int_t ev = (Int_t)esdEvent->GetEventNumber();
204 if(ev<fFirstEvent || ev>fLastEvent) { delete esdEvent; continue; }
205 if(ev % 100 == 0 || fDebug)
206 printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);
207
66e811e2 208 FindPrimaryVertexForCurrentEvent(esdEvent);
11ba84a4 209
210 if(!fCurrentVertex) {
211 printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
212 continue;
213 }
214
66e811e2 215
11ba84a4 216 if(fDebug) fCurrentVertex->PrintStatus();
217
218 // write the ESD to file
219 curdir = gDirectory;
220 Char_t ename[100];
221 sprintf(ename,"%d",ev);
222 fOutFile->cd();
223 esdEvent->Dump();
224 esdEvent->Write(ename,TObject::kOverwrite);
225 curdir->cd();
226 fCurrentVertex = 0;
cab6ff9b 227 } // loop over events
228
229 return;
230}
231//----------------------------------------------------------------------------
232Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
66e811e2 233 //
234 // Propagate tracks to initial vertex position and store them in a TObjArray
235 //
236 Double_t maxd0rphi = 3.;
cab6ff9b 237 Int_t nTrks = 0;
238 Bool_t skipThis;
66e811e2 239 Double_t d0rphi;
cab6ff9b 240 Int_t nEntries = (Int_t)trkTree.GetEntries();
241
242 if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
243 fTrkArray.Expand(nEntries);
244
245 if(fDebug) {
246 printf(" PrepareTracks()\n");
247 trkTree.Print();
248 }
66e811e2 249 cout<<" entr tree its tracks = "<<nEntries<<endl;
cab6ff9b 250 for(Int_t i=0; i<nEntries; i++) {
251 // check tracks to skip
252 skipThis = kFALSE;
253 for(Int_t j=0; j<fNTrksToSkip; j++) {
254 if(i==fTrksToSkip[j]) {
255 if(fDebug) printf("skipping track: %d\n",i);
256 skipThis = kTRUE;
257 }
258 }
259 if(skipThis) continue;
260
261 AliITStrackV2 *itstrack = new AliITStrackV2;
262 trkTree.SetBranchAddress("tracks",&itstrack);
263 trkTree.GetEvent(i);
66e811e2 264 d0rphi=Prepare(itstrack);
265 if(d0rphi> maxd0rphi) { if(fDebug)cout<<" !!!! d0rphi "<<d0rphi<<endl;continue; }
266 fTrkArray.AddLast(itstrack);
267
268 nTrks++;
269 if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<" "<<d0rphi<<endl;
270
271 }
272 if(fTrksToSkip) delete [] fTrksToSkip;
cab6ff9b 273
66e811e2 274 return nTrks;
275}
276//----------------------------------------------------------------------------
277Int_t AliITSVertexerTracks::PrepareTracks(AliESD* esdEvent,Int_t nofCand, Int_t *trkPos) {
278 //
279 // Propagate tracks to initial vertex position and store them in a TObjArray
280 //
281 Int_t nTrks = 0;
282 Double_t maxd0rphi = 3.;
283 Double_t d0rphi;
cab6ff9b 284
66e811e2 285 if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
286 fTrkArray.Expand(100);
cab6ff9b 287
66e811e2 288 if(fDebug) {
289 printf(" PrepareTracks()\n");
290 }
291 AliITStrackV2* itstrack;
292
293 for(Int_t i=0; i<nofCand;i++){
294 AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(trkPos[i]);
295 UInt_t status=esdTrack->GetStatus();
296 if ((status&AliESDtrack::kTPCin)==0)continue;
297 if ((status&AliESDtrack::kITSin)==0)continue;
298 if ((status&AliESDtrack::kITSrefit)==0) continue;
299
300 itstrack = new AliITStrackV2(*esdTrack);
301 d0rphi=Prepare(itstrack);
302 if(d0rphi> maxd0rphi) { if(fDebug)cout<<" !!!! d0rphi "<<d0rphi<<endl;continue; }
303 Int_t nclus=itstrack->GetNumberOfClusters();
304
305 if(nclus<6){delete itstrack;continue;}
cab6ff9b 306 fTrkArray.AddLast(itstrack);
66e811e2 307
cab6ff9b 308 nTrks++;
66e811e2 309 if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<" "<<d0rphi<<endl;
310 //delete itstrack;
cab6ff9b 311 }
66e811e2 312
cab6ff9b 313
314 if(fTrksToSkip) delete [] fTrksToSkip;
cab6ff9b 315 return nTrks;
66e811e2 316}
317//----------------------------------------------------------------------------
318Double_t AliITSVertexerTracks::Prepare(AliITStrackV2* itstrack){
319
320 //
321 Double_t alpha,xlStart,d0rphi;
322 // propagate track to vtxSeed
323 alpha = itstrack->GetAlpha();
324 xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
325 if(itstrack->GetX()>3.)itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
326 itstrack->PropagateTo(xlStart,0.,0.);
327 // select tracks with d0rphi < maxd0rphi
328 d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
329 return d0rphi;
330
331}
cab6ff9b 332//----------------------------------------------------------------------------
333void AliITSVertexerTracks::PrintStatus() const {
334//
335// Print status
336//
337 printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
66e811e2 338 printf(" Vertex position after vertex finder:\n");
339 fSimpVert.Print();
cab6ff9b 340 printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
341 printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
cab6ff9b 342
343 return;
344}
345//----------------------------------------------------------------------------
66e811e2 346AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(Int_t evnumb) {
cab6ff9b 347//
348// Vertex for current event
349//
350 fCurrentVertex = 0;
351
352 // get tree with tracks from input file
353 TString treeName = "TreeT_ITS_";
354 treeName += evnumb;
11ba84a4 355 TTree *trkTree=(TTree*)fInFile->Get(treeName.Data());
cab6ff9b 356 if(!trkTree) return fCurrentVertex;
357
358
359 // get tracks and propagate them to initial vertex position
360 Int_t nTrks = PrepareTracks(*trkTree);
361 delete trkTree;
362 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
363 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
364
365 // VERTEX FINDER
366 VertexFinder();
367
368 // VERTEX FITTER
369 ComputeMaxChi2PerTrack(nTrks);
cab6ff9b 370 VertexFitter();
371 if(fDebug) printf(" vertex fit completed\n");
372
cab6ff9b 373 return fCurrentVertex;
374}
cab6ff9b 375//----------------------------------------------------------------------------
66e811e2 376AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(AliESD *esdEvent)
11ba84a4 377{
cab6ff9b 378//
11ba84a4 379// Vertex for current ESD event
cab6ff9b 380//
11ba84a4 381 fCurrentVertex = 0;
382 Double_t vtx[3],cvtx[6];
383
384 // put tracks reco in ITS in a tree
385 Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
386 TTree *trkTree = new TTree("TreeT_ITS","its tracks");
387 AliITStrackV2 *itstrack = 0;
388 trkTree->Branch("tracks","AliITStrackV2",&itstrack,entr,0);
389
390 for(Int_t i=0; i<entr; i++) {
391 AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(i);
c84a5e9e 392 if(!esdTrack->GetStatus()&AliESDtrack::kITSin)
11ba84a4 393 { delete esdTrack; continue; }
70d3d424 394 try {
395 itstrack = new AliITStrackV2(*esdTrack);
396 }
397 catch (const Char_t *msg) {
66e811e2 398 Warning("FindPrimaryVertexForCurrentEvent",msg);
70d3d424 399 delete esdTrack;
400 continue;
401 }
402
11ba84a4 403 trkTree->Fill();
404 itstrack = 0;
405 delete esdTrack;
406 }
407 delete itstrack;
408
409 // preselect tracks and propagate them to initial vertex position
410 Int_t nTrks = PrepareTracks(*trkTree);
411 delete trkTree;
412 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
413 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
cab6ff9b 414
11ba84a4 415 // Set initial vertex position from ESD
2257f27e 416 esdEvent->GetVertex()->GetXYZ(vtx);
11ba84a4 417 SetVtxStart(vtx[0],vtx[1]);
11ba84a4 418 // VERTEX FINDER
419 VertexFinder();
cab6ff9b 420
11ba84a4 421 // VERTEX FITTER
422 ComputeMaxChi2PerTrack(nTrks);
423 VertexFitter();
424 if(fDebug) printf(" vertex fit completed\n");
cab6ff9b 425
11ba84a4 426 // store vertex information in ESD
427 fCurrentVertex->GetXYZ(vtx);
428 fCurrentVertex->GetCovMatrix(cvtx);
70d3d424 429
430 Double_t tp[3];
431 esdEvent->GetVertex()->GetTruePos(tp);
432 fCurrentVertex->SetTruePos(tp);
433
2257f27e 434 esdEvent->SetVertex(fCurrentVertex);
cab6ff9b 435
11ba84a4 436 cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
437 return fCurrentVertex;
cab6ff9b 438}
66e811e2 439//----------------------------------------------------------------------------
440AliITSSimpleVertex* AliITSVertexerTracks::VertexForSelectedTracks(AliESD *esdEvent,Int_t nofCand, Int_t *trkPos, Int_t opt){
441
442 //
443 // Computes the vertex for selected tracks
444 // trkPos=vector with track positions in ESD
445 // values of opt -> see AliITSVertexerTracks.h
446 //
447 Double_t vtx[3]={0,0,0};
448
449 Int_t nTrks = PrepareTracks(esdEvent,nofCand, trkPos);
450 //delete trkTree;// :-))
451 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
452 if(nTrks < fMinTracks) {
453 fSimpVert.SetXYZ(vtx);
454 fSimpVert.SetDispersion(999);
455 fSimpVert.SetNContributors(-5);
456 return &fSimpVert;
457 }
458
459 // Set initial vertex position from ESD
460 esdEvent->GetVertex()->GetXYZ(vtx);
461 SetVtxStart(vtx[0],vtx[1]);
462 if(opt==1) StrLinVertexFinderMinDist(1);
463 if(opt==2) StrLinVertexFinderMinDist(0);
464 if(opt==3) HelixVertexFinder();
465 if(opt==4) VertexFinder(1);
466 if(opt==5) VertexFinder(0);
467 return &fSimpVert;
468}
469
cab6ff9b 470//---------------------------------------------------------------------------
11ba84a4 471void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
cab6ff9b 472//
11ba84a4 473// Mark the tracks not ot be used in the vertex finding
cab6ff9b 474//
11ba84a4 475 fNTrksToSkip = n;
476 fTrksToSkip = new Int_t[n];
477 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
478 return;
cab6ff9b 479}
480//---------------------------------------------------------------------------
481void AliITSVertexerTracks::TooFewTracks() {
482//
483// When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
484// and the number of tracks to -1
485//
d681bb2d 486 fCurrentVertex = new AliESDVertex(0.,0.,-1);
cab6ff9b 487 return;
488}
489//---------------------------------------------------------------------------
66e811e2 490void AliITSVertexerTracks::VertexFinder(Int_t OptUseWeights) {
cab6ff9b 491
492 // Get estimate of vertex position in (x,y) from tracks DCA
66e811e2 493 // Then this estimate is stored to the data member fSimpVert
cab6ff9b 494 // (previous values are overwritten)
495
cab6ff9b 496
66e811e2 497 Double_t initPos[3];
498 initPos[2] = 0.;
499 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
cab6ff9b 500 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
cab6ff9b 501 Double_t aver[3]={0.,0.,0.};
66e811e2 502 Double_t aversq[3]={0.,0.,0.};
503 Double_t sigmasq[3]={0.,0.,0.};
504 Double_t sigma=0;
cab6ff9b 505 Int_t ncombi = 0;
506 AliITStrackV2 *track1;
507 AliITStrackV2 *track2;
66e811e2 508 Double_t alpha,mindist;
509
cab6ff9b 510 for(Int_t i=0; i<nacc; i++){
511 track1 = (AliITStrackV2*)fTrkArray.At(i);
66e811e2 512 alpha=track1->GetAlpha();
513 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
514 AliITSStrLine *line1 = new AliITSStrLine();
515 track1->ApproximateHelixWithLine(mindist,line1);
516
cab6ff9b 517 if(fDebug>5){
518 Double_t xv,par[5];
519 track1->GetExternalParameters(xv,par);
520 cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
521 for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
522 cout<<endl;
523 }
66e811e2 524
cab6ff9b 525 for(Int_t j=i+1; j<nacc; j++){
526 track2 = (AliITStrackV2*)fTrkArray.At(j);
66e811e2 527 alpha=track2->GetAlpha();
cab6ff9b 528 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
66e811e2 529 AliITSStrLine *line2 = new AliITSStrLine();
530 track2->ApproximateHelixWithLine(mindist,line2);
531 Double_t distCA=line2->GetDCA(line1);
532 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
533 Double_t pnt1[3],pnt2[3],crosspoint[3];
534
535 if(OptUseWeights<=0){
536 Int_t retcode = line2->Cross(line1,crosspoint);
537 if(retcode<0){
538 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
539 if(fDebug>10)cout<<"bad intersection\n";
540 line1->PrintStatus();
541 line2->PrintStatus();
542 }
543 else {
544 ncombi++;
545 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
546 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
547 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
548 if(fDebug>10)cout<<"\n Cross point: ";
549 if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
550 }
551 }
552 if(OptUseWeights>0){
553 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
554 if(retcode>=0){
555 Double_t alpha, cs, sn;
556 alpha=track1->GetAlpha();
557 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
558 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
559 alpha=track2->GetAlpha();
560 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
561 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
562 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
563 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
564 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
565 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
566 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
567 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
568 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
569
570 ncombi++;
571 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
572 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
573 }
574 }
cab6ff9b 575 }
66e811e2 576 delete line2;
577 }
578 delete line1;
579 }
580 if(ncombi>0){
581 for(Int_t jj=0;jj<3;jj++){
582 initPos[jj] = aver[jj]/ncombi;
583 aversq[jj]/=ncombi;
584 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
585 sigma+=sigmasq[jj];
586 }
587 sigma=TMath::Sqrt(TMath::Abs(sigma));
588 }
589 else {
590 Warning("VertexFinder","Finder did not succed");
591 sigma=999;
592 }
593 fSimpVert.SetXYZ(initPos);
594 fSimpVert.SetDispersion(sigma);
595 fSimpVert.SetNContributors(ncombi);
596}
597//---------------------------------------------------------------------------
598void AliITSVertexerTracks::HelixVertexFinder() {
599
600 // Get estimate of vertex position in (x,y) from tracks DCA
601 // Then this estimate is stored to the data member fSimpVert
602 // (previous values are overwritten)
603
604
605 Double_t initPos[3];
606 initPos[2] = 0.;
607 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
608
609 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
610
611 Double_t aver[3]={0.,0.,0.};
612 Double_t averquad[3]={0.,0.,0.};
613 Double_t sigmaquad[3]={0.,0.,0.};
614 Double_t sigma=0;
615 Int_t ncombi = 0;
616 AliITStrackV2 *track1;
617 AliITStrackV2 *track2;
618 Double_t distCA;
619 Double_t x, par[5];
620 Double_t alpha, cs, sn;
621 Double_t crosspoint[3];
622 for(Int_t i=0; i<nacc; i++){
623 track1 = (AliITStrackV2*)fTrkArray.At(i);
624
625
626 for(Int_t j=i+1; j<nacc; j++){
627 track2 = (AliITStrackV2*)fTrkArray.At(j);
628
629
630 distCA=track2->PropagateToDCA(track1);
631
632 if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
633 track1->GetExternalParameters(x,par);
634 alpha=track1->GetAlpha();
635 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
636 Double_t x1=x*cs - par[0]*sn;
637 Double_t y1=x*sn + par[0]*cs;
638 Double_t z1=par[1];
639 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
640
641 track2->GetExternalParameters(x,par);
642 alpha=track2->GetAlpha();
643 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
644 Double_t x2=x*cs - par[0]*sn;
645 Double_t y2=x*sn + par[0]*cs;
646 Double_t z2=par[1];
647 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
648 // printf("Track %d pos=(%f,%f,%f) - dca=%f\n",i,x1,y1,z1,distCA);
649 //printf("Track %d pos=(%f,%f,%f)\n",j,x2,y2,z2);
650
651 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
652 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
653 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
654 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
655 crosspoint[0]=wx1*x1 + wx2*x2;
656 crosspoint[1]=wy1*y1 + wy2*y2;
657 crosspoint[2]=wz1*z1 + wz2*z2;
658
cab6ff9b 659 ncombi++;
660 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
66e811e2 661 for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
cab6ff9b 662 }
cab6ff9b 663 }
66e811e2 664
cab6ff9b 665 }
666 if(ncombi>0){
66e811e2 667 for(Int_t jj=0;jj<3;jj++){
668 initPos[jj] = aver[jj]/ncombi;
669 averquad[jj]/=ncombi;
670 sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
671 sigma+=sigmaquad[jj];
672 }
673 sigma=TMath::Sqrt(TMath::Abs(sigma));
cab6ff9b 674 }
675 else {
676 Warning("VertexFinder","Finder did not succed");
66e811e2 677 sigma=999;
678 }
679 fSimpVert.SetXYZ(initPos);
680 fSimpVert.SetDispersion(sigma);
681 fSimpVert.SetNContributors(ncombi);
682}
683//---------------------------------------------------------------------------
684void AliITSVertexerTracks::StrLinVertexFinderMinDist(Int_t OptUseWeights){
685
686 // Calculate the point at minimum distance to prepared tracks
687 // Then this estimate is stored to the data member fSimpVert
688 // (previous values are overwritten)
689
690 Double_t initPos[3];
691 initPos[2] = 0.;
692 Double_t sigma=0;
693 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
694 const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
695
696 AliITStrackV2 *track1;
697 Double_t (*vectP0)[3]=new Double_t [knacc][3];
698 Double_t (*vectP1)[3]=new Double_t [knacc][3];
699
700 Double_t sum[3][3];
701 Double_t dsum[3]={0,0,0};
702 for(Int_t i=0;i<3;i++)
703 for(Int_t j=0;j<3;j++)sum[i][j]=0;
704 for(Int_t i=0; i<knacc; i++){
705 track1 = (AliITStrackV2*)fTrkArray.At(i);
706 Double_t alpha=track1->GetAlpha();
707 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
708 AliITSStrLine *line1 = new AliITSStrLine();
709 track1->ApproximateHelixWithLine(mindist,line1);
710
711 Double_t p0[3],cd[3];
712 line1->GetP0(p0);
713 line1->GetCd(cd);
714 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
715 vectP0[i][0]=p0[0];
716 vectP0[i][1]=p0[1];
717 vectP0[i][2]=p0[2];
718 vectP1[i][0]=p1[0];
719 vectP1[i][1]=p1[1];
720 vectP1[i][2]=p1[2];
721
722 Double_t matr[3][3];
723 Double_t dknow[3];
724 if(OptUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
725 if(OptUseWeights==1){
726 Double_t sigmasq[3];
727 sigmasq[0]=track1->GetSigmaY2();
728 sigmasq[1]=track1->GetSigmaY2();
729 sigmasq[2]=track1->GetSigmaZ2();
730 GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
731 }
732
733 for(Int_t iii=0;iii<3;iii++){
734 dsum[iii]+=dknow[iii];
735 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
736 }
737 delete line1;
738 }
739
740 Double_t vett[3][3];
741 Double_t det=GetDeterminant3X3(sum);
742
743 if(det!=0){
744 for(Int_t zz=0;zz<3;zz++){
745 for(Int_t ww=0;ww<3;ww++){
746 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
747 }
748 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
749 initPos[zz]=GetDeterminant3X3(vett)/det;
750 }
751
752
753 for(Int_t i=0; i<knacc; i++){
754 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
755 for(Int_t ii=0;ii<3;ii++){
756 p0[ii]=vectP0[i][ii];
757 p1[ii]=vectP1[i][ii];
758 }
759 sigma+=GetStrLinMinDist(p0,p1,initPos);
760 }
761
762 sigma=TMath::Sqrt(sigma);
763 }else{
764 Warning("VertexFinder","Finder did not succed");
765 sigma=999;
cab6ff9b 766 }
66e811e2 767 delete vectP0;
768 delete vectP1;
769 fSimpVert.SetXYZ(initPos);
770 fSimpVert.SetDispersion(sigma);
771 fSimpVert.SetNContributors(knacc);
772}
773//_______________________________________________________________________
774Double_t AliITSVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
775 //
776 Double_t det=matr[0][0]*matr[1][1]*matr[2][2]-matr[0][0]*matr[1][2]*matr[2][1]-matr[0][1]*matr[1][0]*matr[2][2]+matr[0][1]*matr[1][2]*matr[2][0]+matr[0][2]*matr[1][0]*matr[2][1]-matr[0][2]*matr[1][1]*matr[2][0];
777 return det;
778}
779//____________________________________________________________________________
780void AliITSVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t m[][3],Double_t *d){
781
782 //
783 Double_t x12=p0[0]-p1[0];
784 Double_t y12=p0[1]-p1[1];
785 Double_t z12=p0[2]-p1[2];
786 Double_t kk=x12*x12+y12*y12+z12*z12;
787 m[0][0]=2-2/kk*x12*x12;
788 m[0][1]=-2/kk*x12*y12;
789 m[0][2]=-2/kk*x12*z12;
790 m[1][0]=-2/kk*x12*y12;
791 m[1][1]=2-2/kk*y12*y12;
792 m[1][2]=-2/kk*y12*z12;
793 m[2][0]=-2/kk*x12*z12;
794 m[2][1]=-2*y12*z12;
795 m[2][2]=2-2/kk*z12*z12;
796 d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
797 d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
798 d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
cab6ff9b 799
66e811e2 800}
801//____________________________________________________________________________
802void AliITSVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t m[][3],Double_t *d){
803 //
804 Double_t x12=p1[0]-p0[0];
805 Double_t y12=p1[1]-p0[1];
806 Double_t z12=p1[2]-p0[2];
cab6ff9b 807
66e811e2 808 Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
809
810 Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
811
812 Double_t cc[3];
813 cc[0]=-x12/sigmasq[0];
814 cc[1]=-y12/sigmasq[1];
815 cc[2]=-z12/sigmasq[2];
816
817 Double_t ww=(-p0[0]*x12*sigmasq[1]*sigmasq[2]-p0[1]*y12*sigmasq[0]*sigmasq[2]-p0[2]*z12*sigmasq[0]*sigmasq[1])/den;
818
819 Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
820
821 Double_t aa[3];
822 aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
823 aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
824 aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
825
826 m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
827 m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
828 m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
829
830 m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
831 m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
832 m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
833
834 m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
835 m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
836 m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
837
838 d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
839 d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
840 d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
841
842}
843//_____________________________________________________________________________
844Double_t AliITSVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
845 //
846 Double_t x12=p0[0]-p1[0];
847 Double_t y12=p0[1]-p1[1];
848 Double_t z12=p0[2]-p1[2];
849 Double_t x10=p0[0]-x0[0];
850 Double_t y10=p0[1]-x0[1];
851 Double_t z10=p0[2]-x0[2];
852 return ((x10*x10+y10*y10+z10*z10)*(x12*x12+y12*y12+z12*z12)-(x10*x12+y10*y12+z10*z12)*(x10*x12+y10*y12+z10*z12))/(x12*x12+y12*y12+z12*z12);
cab6ff9b 853
854}
cab6ff9b 855//---------------------------------------------------------------------------
856void AliITSVertexerTracks::VertexFitter() {
857//
858// The optimal estimate of the vertex position is given by a "weighted
859// average of tracks positions"
860// Original method: CMS Note 97/0051
861//
862 if(fDebug) {
863 printf(" VertexFitter(): start\n");
864 PrintStatus();
865 }
866
867
868 Int_t i,j,k,step=0;
869 TMatrixD rv(3,1);
a48d1ed3 870 TMatrixD vV(3,3);
66e811e2 871 Double_t initPos[3];
872 fSimpVert.GetXYZ(initPos);
873 rv(0,0) = initPos[0];
874 rv(1,0) = initPos[1];
cab6ff9b 875 rv(2,0) = 0.;
876 Double_t xlStart,alpha;
877 Double_t rotAngle;
878 Double_t cosRot,sinRot;
879 Double_t cc[15];
880 Int_t nUsedTrks;
881 Double_t chi2,chi2i;
882 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
883 AliITStrackV2 *t = 0;
884 Int_t failed = 0;
885
886 Int_t *skipTrack = new Int_t[arrEntries];
887 for(i=0; i<arrEntries; i++) skipTrack[i]=0;
888
889 // 3 steps:
890 // 1st - first estimate of vtx using all tracks
891 // 2nd - apply cut on chi2 max per track
892 // 3rd - estimate of global chi2
893 for(step=0; step<3; step++) {
894 if(fDebug) printf(" step = %d\n",step);
895 chi2 = 0.;
896 nUsedTrks = 0;
897
a48d1ed3 898 TMatrixD sumWiri(3,1);
899 TMatrixD sumWi(3,3);
cab6ff9b 900 for(i=0; i<3; i++) {
a48d1ed3 901 sumWiri(i,0) = 0.;
902 for(j=0; j<3; j++) sumWi(j,i) = 0.;
cab6ff9b 903 }
904
905 // loop on tracks
906 for(k=0; k<arrEntries; k++) {
907 if(skipTrack[k]) continue;
908 // get track from track array
909 t = (AliITStrackV2*)fTrkArray.At(k);
910 alpha = t->GetAlpha();
66e811e2 911 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
cab6ff9b 912 t->PropagateTo(xlStart,0.,0.); // to vtxSeed
11ba84a4 913 rotAngle = alpha;
cab6ff9b 914 if(alpha<0.) rotAngle += 2.*TMath::Pi();
915 cosRot = TMath::Cos(rotAngle);
916 sinRot = TMath::Sin(rotAngle);
917
918 // vector of track global coordinates
919 TMatrixD ri(3,1);
920 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
921 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
922 ri(2,0) = t->GetZ();
923
924 // matrix to go from global (x,y,z) to local (y,z);
a48d1ed3 925 TMatrixD qQi(2,3);
926 qQi(0,0) = -sinRot;
927 qQi(0,1) = cosRot;
928 qQi(0,2) = 0.;
929 qQi(1,0) = 0.;
930 qQi(1,1) = 0.;
931 qQi(1,2) = 1.;
cab6ff9b 932
933 // covariance matrix of local (y,z) - inverted
a48d1ed3 934 TMatrixD uUi(2,2);
cab6ff9b 935 t->GetExternalCovariance(cc);
a48d1ed3 936 uUi(0,0) = cc[0];
937 uUi(0,1) = cc[1];
938 uUi(1,0) = cc[1];
939 uUi(1,1) = cc[2];
cab6ff9b 940
a48d1ed3 941 // weights matrix: wWi = qQiT * uUiInv * qQi
942 if(uUi.Determinant() <= 0.) continue;
943 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
944 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
945 TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
cab6ff9b 946
947 // track chi2
948 TMatrixD deltar = rv; deltar -= ri;
a48d1ed3 949 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
950 chi2i = deltar(0,0)*wWideltar(0,0)+
951 deltar(1,0)*wWideltar(1,0)+
952 deltar(2,0)*wWideltar(2,0);
cab6ff9b 953
954
955 if(step==1 && chi2i > fMaxChi2PerTrack) {
956 skipTrack[k] = 1;
957 continue;
958 }
959
960 // add to total chi2
961 chi2 += chi2i;
962
a48d1ed3 963 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
cab6ff9b 964
a48d1ed3 965 sumWiri += wWiri;
966 sumWi += wWi;
cab6ff9b 967
968 nUsedTrks++;
969 } // end loop on tracks
970
971 if(nUsedTrks < fMinTracks) {
972 failed=1;
973 continue;
974 }
975
a48d1ed3 976 Double_t determinant = sumWi.Determinant();
cab6ff9b 977 //cerr<<" determinant: "<<determinant<<endl;
978 if(determinant < 100.) {
979 printf("det(V) = 0\n");
980 failed=1;
981 continue;
982 }
983
984 // inverted of weights matrix
a48d1ed3 985 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
986 vV = invsumWi;
cab6ff9b 987
988 // position of primary vertex
a48d1ed3 989 rv.Mult(vV,sumWiri);
cab6ff9b 990
991 } // end loop on the 3 steps
992
993 delete [] skipTrack;
994 delete t;
995
996 if(failed) {
997 TooFewTracks();
998 return;
999 }
1000
1001 Double_t position[3];
1002 position[0] = rv(0,0);
1003 position[1] = rv(1,0);
1004 position[2] = rv(2,0);
1005 Double_t covmatrix[6];
a48d1ed3 1006 covmatrix[0] = vV(0,0);
1007 covmatrix[1] = vV(0,1);
1008 covmatrix[2] = vV(1,1);
1009 covmatrix[3] = vV(0,2);
1010 covmatrix[4] = vV(1,2);
1011 covmatrix[5] = vV(2,2);
cab6ff9b 1012
1013 // store data in the vertex object
d681bb2d 1014 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
cab6ff9b 1015
1016 if(fDebug) {
1017 printf(" VertexFitter(): finish\n");
1018 printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
1019 fCurrentVertex->PrintStatus();
1020 }
1021
1022 return;
1023}
1024//----------------------------------------------------------------------------
d681bb2d 1025AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
cab6ff9b 1026//
1027// Return vertex from tracks in trkTree
1028//
1029 if(fCurrentVertex) fCurrentVertex = 0;
1030
1031 // get tracks and propagate them to initial vertex position
1032 Int_t nTrks = PrepareTracks(*(&trkTree));
1033 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
1034 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
1035
1036 // VERTEX FINDER
1037 VertexFinder();
1038
1039 // VERTEX FITTER
1040 ComputeMaxChi2PerTrack(nTrks);
cab6ff9b 1041 VertexFitter();
1042 if(fDebug) printf(" vertex fit completed\n");
1043
1044 return fCurrentVertex;
1045}
1046//----------------------------------------------------------------------------
1047
1048
1049