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