--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+//-- --- standard headers-------------
+#include <Riostream.h>
+//--------Root headers ---------------
+#include <TSystem.h>
+#include <TFile.h>
+#include <TStopwatch.h>
+#include <TObject.h>
+#include <TTree.h>
+#include <TChain.h>
+#include <TF1.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TCut.h>
+#include <TParticle.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TGraph.h>
+//----- AliRoot headers ---------------
+//#include "alles.h"
+#include "AliRun.h"
+#include "AliMagF.h"
+#include "AliKalmanTrack.h"
+#include "AliESDVertex.h"
+#include "AliITSVertexer.h"
+#include "AliITSVertexerTracks.h"
+#include <AliHeader.h>
+#include <AliGenEventHeader.h>
+#include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliITStrackV2.h"
+#include "AliITSSimpleVertex.h"
+#include "AliITSStrLine.h"
+#include "AliITSLoader.h"
+#include <AliESD.h>
+#include <AliStack.h>
+//-------------------------------------
+#endif
+
+
+void AliITSSecondaryVertices(Int_t ntracks, Int_t *pos) {
+
+ // example of usige AliITSVertexerTracks to get a (secodnary) vertex from a set of tracks.
+
+ Double_t field = 0.4;
+
+ AliRunLoader *rl = AliRunLoader::Open("galice.root");
+ if (!rl) {
+ cerr<<"Can't start session !\n";
+ }
+
+ rl->LoadgAlice();
+ if (rl->GetAliRun()){
+ AliMagF * magf = gAlice->Field();
+ AliKalmanTrack:: SetFieldMap(magf);
+ AliKalmanTrack::SetUniformFieldTracking();
+ }else {
+ cerr<<"Can't get AliRun !\n";
+
+ }
+
+ field=rl->GetAliRun()->Field()->SolenoidField()/10.;
+ rl->LoadHeader();
+
+
+ TFile *ef=TFile::Open("AliESDs.root");
+ if (!ef || !ef->IsOpen()) {cerr<<"Can't AliESDs.root !n";}
+ AliESD* event = new AliESD;
+ TTree* tree = (TTree*) ef->Get("esdTree");
+ if (!tree) {cerr<<"no ESD tree foundn";};
+ tree->SetBranchAddress("ESD", &event);
+
+
+ Int_t e=0;
+ tree->GetEvent(e);
+ Int_t ntrk=event->GetNumberOfTracks();
+ cout<<"Number of ESD tracks : "<<ntrk<<endl;
+
+ // Open input and output files
+ TFile *outFile = TFile::Open("AliITSVertices.root","recreate");
+
+
+ TTree *trkTree = new TTree("TreeT_ITS","its tracks");
+ AliITStrackV2 *itstrack = 0;
+ trkTree->Branch("tracks","AliITStrackV2",&itstrack,ntrk,0);
+ Int_t pipe=0;
+ for(Int_t i=0; i<ntrk; i++) {
+ AliESDtrack *esdTrack = (AliESDtrack*)event->GetTrack(i);
+ UInt_t status=esdTrack->GetStatus();
+ UInt_t flags=AliESDtrack::kTPCin|AliESDtrack::kITSin;
+ if ((status&AliESDtrack::kTPCin)==0)continue;
+ if ((status&AliESDtrack::kITSin)==0)continue;
+ if ((status&AliESDtrack::kITSrefit)==0) continue;
+ if ((status&flags)==status) pipe=1;
+
+ itstrack = new AliITStrackV2(*esdTrack);
+ if (pipe!=0) {
+ itstrack->PropagateTo(3.,0.0028,65.19);
+ itstrack->PropagateToVertex(); // This is not "exactly" the vertex
+ }
+ Int_t nclus=itstrack->GetNumberOfClusters();
+
+ if(nclus<6)continue;
+ trkTree->Fill();
+ itstrack = 0;
+ // delete esdTrack;
+ }
+ // Primary vertex
+ Double_t vtx[3];
+ event->GetVertex()->GetXYZ(vtx);
+ cout<<"Primary Vertex:"<<endl;
+ event->GetVertex()->PrintStatus();
+ cout<<"\n\n\n"<<endl;
+
+
+ // Create vertexer
+ AliITSVertexerTracks *vertexer = new AliITSVertexerTracks();
+ vertexer->SetDCAcut(0);
+ vertexer->SetDebug(0);
+
+
+ Double_t vertF1[3],vertF2[3],vertF3[3],vertF4[3],vertF5[3];
+ Int_t ncombi1,ncombi2,ncombi3,ncombi4,ncombi5;
+ Double_t sig1,sig2,sig3,sig4,sig5;
+ AliITSSimpleVertex *vert=0;
+
+
+ // Calculate vertex from tracks in pos
+ // BEWARE: the method VertexForSelectedTracks returns the address of the data member fSimpVert
+ // which is always the same for all calls to VertexForSelectedTracks.
+
+ cout<<"\nStraight Line Min Dist finder with weights"<<endl;
+ vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,1);
+ vert->Print();
+ vert->GetXYZ(vertF1);
+ ncombi1=vert->GetNContributors();
+ sig1=vert->GetDispersion();
+
+
+ cout<<"Straight Line Min Dist finder"<<endl;
+ vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,2);
+ vert->Print();
+ vert->GetXYZ(vertF2);
+ ncombi2=vert->GetNContributors();
+ sig2=vert->GetDispersion();
+
+ cout<<"Helix finder"<<endl;
+ vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,3);
+ vert->Print();
+ vert->GetXYZ(vertF3);
+ ncombi3=vert->GetNContributors();
+ sig3=vert->GetDispersion();
+
+ cout<<"Straight Line finder with weights"<<endl;
+ vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,4);
+ vert->Print();
+ vert->GetXYZ(vertF4);
+ ncombi4=vert->GetNContributors();
+ sig4=vert->GetDispersion();
+
+ cout<<"Straight Line finder with weights"<<endl;
+ vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,5);
+ vert->Print();
+ vert->GetXYZ(vertF5);
+ ncombi5=vert->GetNContributors();
+ sig5=vert->GetDispersion();
+
+
+ delete vertexer;
+
+ outFile->Close();
+ return;
+}
SetMinTracks();
fTrksToSkip = 0;
fNTrksToSkip = 0;
- for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
+ fDCAcut=0;
}
//----------------------------------------------------------------------------
AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
SetMinTracks();
fTrksToSkip = 0;
fNTrksToSkip = 0;
- for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
+ fDCAcut=0;
SetDebug();
}
//----------------------------------------------------------------------------
SetMinTracks();
fTrksToSkip = 0;
fNTrksToSkip = 0;
- for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
+ fDCAcut=0;
}
//______________________________________________________________________
AliITSVertexerTracks::AliITSVertexerTracks(const AliITSVertexerTracks &vtxr) : AliITSVertexer(vtxr) {
for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
- FindVertexForCurrentEvent(ev);
+ FindPrimaryVertexForCurrentEvent(ev);
if(!fCurrentVertex) {
printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
if(ev % 100 == 0 || fDebug)
printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);
- FindVertexForCurrentEvent(esdEvent);
+ FindPrimaryVertexForCurrentEvent(esdEvent);
if(!fCurrentVertex) {
printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
continue;
}
- cout<<"VERTICE TROVATO\n";
+
if(fDebug) fCurrentVertex->PrintStatus();
// write the ESD to file
}
//----------------------------------------------------------------------------
Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
-//
-// Propagate tracks to initial vertex position and store them in a TObjArray
-//
- Double_t maxd0rphi = 3.;
- Double_t alpha,xlStart,d0rphi;
+ //
+ // Propagate tracks to initial vertex position and store them in a TObjArray
+ //
+ Double_t maxd0rphi = 3.;
Int_t nTrks = 0;
Bool_t skipThis;
-
+ Double_t d0rphi;
Int_t nEntries = (Int_t)trkTree.GetEntries();
if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
printf(" PrepareTracks()\n");
trkTree.Print();
}
-
+ cout<<" entr tree its tracks = "<<nEntries<<endl;
for(Int_t i=0; i<nEntries; i++) {
// check tracks to skip
skipThis = kFALSE;
AliITStrackV2 *itstrack = new AliITStrackV2;
trkTree.SetBranchAddress("tracks",&itstrack);
trkTree.GetEvent(i);
+ d0rphi=Prepare(itstrack);
+ if(d0rphi> maxd0rphi) { if(fDebug)cout<<" !!!! d0rphi "<<d0rphi<<endl;continue; }
+ fTrkArray.AddLast(itstrack);
+
+ nTrks++;
+ if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<" "<<d0rphi<<endl;
+
+ }
+ if(fTrksToSkip) delete [] fTrksToSkip;
+ return nTrks;
+}
+//----------------------------------------------------------------------------
+Int_t AliITSVertexerTracks::PrepareTracks(AliESD* esdEvent,Int_t nofCand, Int_t *trkPos) {
+ //
+ // Propagate tracks to initial vertex position and store them in a TObjArray
+ //
+ Int_t nTrks = 0;
+ Double_t maxd0rphi = 3.;
+ Double_t d0rphi;
- // propagate track to vtxSeed
- alpha = itstrack->GetAlpha();
- xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
- itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
- itstrack->PropagateTo(xlStart,0.,0.); // to vtxSeed
+ if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
+ fTrkArray.Expand(100);
- // select tracks with d0rphi < maxd0rphi
- d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
- if(d0rphi > maxd0rphi) { delete itstrack; continue; }
-
+ if(fDebug) {
+ printf(" PrepareTracks()\n");
+ }
+ AliITStrackV2* itstrack;
+
+ for(Int_t i=0; i<nofCand;i++){
+ AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(trkPos[i]);
+ UInt_t status=esdTrack->GetStatus();
+ if ((status&AliESDtrack::kTPCin)==0)continue;
+ if ((status&AliESDtrack::kITSin)==0)continue;
+ if ((status&AliESDtrack::kITSrefit)==0) continue;
+
+ itstrack = new AliITStrackV2(*esdTrack);
+ d0rphi=Prepare(itstrack);
+ if(d0rphi> maxd0rphi) { if(fDebug)cout<<" !!!! d0rphi "<<d0rphi<<endl;continue; }
+ Int_t nclus=itstrack->GetNumberOfClusters();
+
+ if(nclus<6){delete itstrack;continue;}
fTrkArray.AddLast(itstrack);
-
+
nTrks++;
+ if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<" "<<d0rphi<<endl;
+ //delete itstrack;
}
+
if(fTrksToSkip) delete [] fTrksToSkip;
-
return nTrks;
-}
+}
+//----------------------------------------------------------------------------
+Double_t AliITSVertexerTracks::Prepare(AliITStrackV2* itstrack){
+
+ //
+ Double_t alpha,xlStart,d0rphi;
+ // propagate track to vtxSeed
+ alpha = itstrack->GetAlpha();
+ xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
+ if(itstrack->GetX()>3.)itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
+ itstrack->PropagateTo(xlStart,0.,0.);
+ // select tracks with d0rphi < maxd0rphi
+ d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
+ return d0rphi;
+
+}
//----------------------------------------------------------------------------
void AliITSVertexerTracks::PrintStatus() const {
//
// Print status
//
printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
- printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
+ printf(" Vertex position after vertex finder:\n");
+ fSimpVert.Print();
printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
return;
}
//----------------------------------------------------------------------------
-AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
+AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(Int_t evnumb) {
//
// Vertex for current event
//
return fCurrentVertex;
}
//----------------------------------------------------------------------------
-AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
+AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(AliESD *esdEvent)
{
//
// Vertex for current ESD event
itstrack = new AliITStrackV2(*esdTrack);
}
catch (const Char_t *msg) {
- Warning("FindVertexForCurrentEvent",msg);
+ Warning("FindPrimaryVertexForCurrentEvent",msg);
delete esdTrack;
continue;
}
// Set initial vertex position from ESD
esdEvent->GetVertex()->GetXYZ(vtx);
SetVtxStart(vtx[0],vtx[1]);
-
// VERTEX FINDER
VertexFinder();
cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
return fCurrentVertex;
}
+//----------------------------------------------------------------------------
+AliITSSimpleVertex* AliITSVertexerTracks::VertexForSelectedTracks(AliESD *esdEvent,Int_t nofCand, Int_t *trkPos, Int_t opt){
+
+ //
+ // Computes the vertex for selected tracks
+ // trkPos=vector with track positions in ESD
+ // values of opt -> see AliITSVertexerTracks.h
+ //
+ Double_t vtx[3]={0,0,0};
+
+ Int_t nTrks = PrepareTracks(esdEvent,nofCand, trkPos);
+ //delete trkTree;// :-))
+ if(fDebug) printf(" tracks prepared: %d\n",nTrks);
+ if(nTrks < fMinTracks) {
+ fSimpVert.SetXYZ(vtx);
+ fSimpVert.SetDispersion(999);
+ fSimpVert.SetNContributors(-5);
+ return &fSimpVert;
+ }
+
+ // Set initial vertex position from ESD
+ esdEvent->GetVertex()->GetXYZ(vtx);
+ SetVtxStart(vtx[0],vtx[1]);
+ if(opt==1) StrLinVertexFinderMinDist(1);
+ if(opt==2) StrLinVertexFinderMinDist(0);
+ if(opt==3) HelixVertexFinder();
+ if(opt==4) VertexFinder(1);
+ if(opt==5) VertexFinder(0);
+ return &fSimpVert;
+}
+
//---------------------------------------------------------------------------
void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
//
return;
}
//---------------------------------------------------------------------------
-void AliITSVertexerTracks::VertexFinder() {
+void AliITSVertexerTracks::VertexFinder(Int_t OptUseWeights) {
// Get estimate of vertex position in (x,y) from tracks DCA
- // Then this estimate is stored to the data member fInitPos
+ // Then this estimate is stored to the data member fSimpVert
// (previous values are overwritten)
-
- /*
-******* TEMPORARY!!! FOR TEST ONLY!!! **********************************
-
-fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
-fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
- */
-
- fInitPos[2] = 0.;
- for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
-
+ Double_t initPos[3];
+ initPos[2] = 0.;
+ for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
-
Double_t aver[3]={0.,0.,0.};
+ Double_t aversq[3]={0.,0.,0.};
+ Double_t sigmasq[3]={0.,0.,0.};
+ Double_t sigma=0;
Int_t ncombi = 0;
AliITStrackV2 *track1;
AliITStrackV2 *track2;
- AliITSStrLine line1;
- AliITSStrLine line2;
+ Double_t alpha,mindist;
+
for(Int_t i=0; i<nacc; i++){
track1 = (AliITStrackV2*)fTrkArray.At(i);
+ alpha=track1->GetAlpha();
+ mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
+ AliITSStrLine *line1 = new AliITSStrLine();
+ track1->ApproximateHelixWithLine(mindist,line1);
+
if(fDebug>5){
Double_t xv,par[5];
track1->GetExternalParameters(xv,par);
for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
cout<<endl;
}
- Double_t mom1[3];
- Double_t alpha = track1->GetAlpha();
- Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
- Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
- mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
- mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
- mom1[2] = TMath::Cos(theta);
-
- Double_t pos1[3];
- Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
- track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
- line1.SetP0(pos1);
- line1.SetCd(mom1);
+
for(Int_t j=i+1; j<nacc; j++){
track2 = (AliITStrackV2*)fTrkArray.At(j);
- Double_t mom2[3];
- alpha = track2->GetAlpha();
- azim = TMath::ASin(track2->GetSnp())+alpha;
- theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
- mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
- mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
- mom2[2] = TMath::Cos(theta);
- Double_t pos2[3];
+ alpha=track2->GetAlpha();
mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
- track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
- line2.SetP0(pos2);
- line2.SetCd(mom2);
- Double_t crosspoint[3];
- Int_t retcode = line2.Cross(&line1,crosspoint);
- if(retcode<0){
- if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
- if(fDebug>10)cout<<"bad intersection\n";
- line1.PrintStatus();
- line2.PrintStatus();
+ AliITSStrLine *line2 = new AliITSStrLine();
+ track2->ApproximateHelixWithLine(mindist,line2);
+ Double_t distCA=line2->GetDCA(line1);
+ if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
+ Double_t pnt1[3],pnt2[3],crosspoint[3];
+
+ if(OptUseWeights<=0){
+ Int_t retcode = line2->Cross(line1,crosspoint);
+ if(retcode<0){
+ if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
+ if(fDebug>10)cout<<"bad intersection\n";
+ line1->PrintStatus();
+ line2->PrintStatus();
+ }
+ else {
+ ncombi++;
+ for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
+ for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
+ if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
+ if(fDebug>10)cout<<"\n Cross point: ";
+ if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
+ }
+ }
+ if(OptUseWeights>0){
+ Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
+ if(retcode>=0){
+ Double_t alpha, cs, sn;
+ alpha=track1->GetAlpha();
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
+ alpha=track2->GetAlpha();
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
+ Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
+ Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
+ Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
+ Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
+ crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
+ crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
+ crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
+
+ ncombi++;
+ for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
+ for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
+ }
+ }
}
- else {
+ delete line2;
+ }
+ delete line1;
+ }
+ if(ncombi>0){
+ for(Int_t jj=0;jj<3;jj++){
+ initPos[jj] = aver[jj]/ncombi;
+ aversq[jj]/=ncombi;
+ sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
+ sigma+=sigmasq[jj];
+ }
+ sigma=TMath::Sqrt(TMath::Abs(sigma));
+ }
+ else {
+ Warning("VertexFinder","Finder did not succed");
+ sigma=999;
+ }
+ fSimpVert.SetXYZ(initPos);
+ fSimpVert.SetDispersion(sigma);
+ fSimpVert.SetNContributors(ncombi);
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::HelixVertexFinder() {
+
+ // Get estimate of vertex position in (x,y) from tracks DCA
+ // Then this estimate is stored to the data member fSimpVert
+ // (previous values are overwritten)
+
+
+ Double_t initPos[3];
+ initPos[2] = 0.;
+ for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
+
+ Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
+
+ Double_t aver[3]={0.,0.,0.};
+ Double_t averquad[3]={0.,0.,0.};
+ Double_t sigmaquad[3]={0.,0.,0.};
+ Double_t sigma=0;
+ Int_t ncombi = 0;
+ AliITStrackV2 *track1;
+ AliITStrackV2 *track2;
+ Double_t distCA;
+ Double_t x, par[5];
+ Double_t alpha, cs, sn;
+ Double_t crosspoint[3];
+ for(Int_t i=0; i<nacc; i++){
+ track1 = (AliITStrackV2*)fTrkArray.At(i);
+
+
+ for(Int_t j=i+1; j<nacc; j++){
+ track2 = (AliITStrackV2*)fTrkArray.At(j);
+
+
+ distCA=track2->PropagateToDCA(track1);
+
+ if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
+ track1->GetExternalParameters(x,par);
+ alpha=track1->GetAlpha();
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ Double_t x1=x*cs - par[0]*sn;
+ Double_t y1=x*sn + par[0]*cs;
+ Double_t z1=par[1];
+ Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
+
+ track2->GetExternalParameters(x,par);
+ alpha=track2->GetAlpha();
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ Double_t x2=x*cs - par[0]*sn;
+ Double_t y2=x*sn + par[0]*cs;
+ Double_t z2=par[1];
+ Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
+ // printf("Track %d pos=(%f,%f,%f) - dca=%f\n",i,x1,y1,z1,distCA);
+ //printf("Track %d pos=(%f,%f,%f)\n",j,x2,y2,z2);
+
+ Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
+ Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
+ Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
+ Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
+ crosspoint[0]=wx1*x1 + wx2*x2;
+ crosspoint[1]=wy1*y1 + wy2*y2;
+ crosspoint[2]=wz1*z1 + wz2*z2;
+
ncombi++;
for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
- if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
- if(fDebug>10)cout<<"\n Cross point: ";
- if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
+ for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
}
}
+
}
if(ncombi>0){
- for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
+ for(Int_t jj=0;jj<3;jj++){
+ initPos[jj] = aver[jj]/ncombi;
+ averquad[jj]/=ncombi;
+ sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
+ sigma+=sigmaquad[jj];
+ }
+ sigma=TMath::Sqrt(TMath::Abs(sigma));
}
else {
Warning("VertexFinder","Finder did not succed");
+ sigma=999;
}
+ fSimpVert.SetXYZ(initPos);
+ fSimpVert.SetDispersion(sigma);
+ fSimpVert.SetNContributors(ncombi);
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::StrLinVertexFinderMinDist(Int_t OptUseWeights){
+ // Calculate the point at minimum distance to prepared tracks
+ // Then this estimate is stored to the data member fSimpVert
+ // (previous values are overwritten)
+
+ Double_t initPos[3];
+ initPos[2] = 0.;
+ Double_t sigma=0;
+ for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
+ const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
- //************************************************************************
- return;
+ AliITStrackV2 *track1;
+ Double_t (*vectP0)[3]=new Double_t [knacc][3];
+ Double_t (*vectP1)[3]=new Double_t [knacc][3];
+
+ Double_t sum[3][3];
+ Double_t dsum[3]={0,0,0};
+ for(Int_t i=0;i<3;i++)
+ for(Int_t j=0;j<3;j++)sum[i][j]=0;
+ for(Int_t i=0; i<knacc; i++){
+ track1 = (AliITStrackV2*)fTrkArray.At(i);
+ Double_t alpha=track1->GetAlpha();
+ Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
+ AliITSStrLine *line1 = new AliITSStrLine();
+ track1->ApproximateHelixWithLine(mindist,line1);
+
+ Double_t p0[3],cd[3];
+ line1->GetP0(p0);
+ line1->GetCd(cd);
+ Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
+ vectP0[i][0]=p0[0];
+ vectP0[i][1]=p0[1];
+ vectP0[i][2]=p0[2];
+ vectP1[i][0]=p1[0];
+ vectP1[i][1]=p1[1];
+ vectP1[i][2]=p1[2];
+
+ Double_t matr[3][3];
+ Double_t dknow[3];
+ if(OptUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
+ if(OptUseWeights==1){
+ Double_t sigmasq[3];
+ sigmasq[0]=track1->GetSigmaY2();
+ sigmasq[1]=track1->GetSigmaY2();
+ sigmasq[2]=track1->GetSigmaZ2();
+ GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
+ }
+
+ for(Int_t iii=0;iii<3;iii++){
+ dsum[iii]+=dknow[iii];
+ for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
+ }
+ delete line1;
+ }
+
+ Double_t vett[3][3];
+ Double_t det=GetDeterminant3X3(sum);
+
+ if(det!=0){
+ for(Int_t zz=0;zz<3;zz++){
+ for(Int_t ww=0;ww<3;ww++){
+ for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
+ }
+ for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
+ initPos[zz]=GetDeterminant3X3(vett)/det;
+ }
+
+
+ for(Int_t i=0; i<knacc; i++){
+ Double_t p0[3]={0,0,0},p1[3]={0,0,0};
+ for(Int_t ii=0;ii<3;ii++){
+ p0[ii]=vectP0[i][ii];
+ p1[ii]=vectP1[i][ii];
+ }
+ sigma+=GetStrLinMinDist(p0,p1,initPos);
+ }
+
+ sigma=TMath::Sqrt(sigma);
+ }else{
+ Warning("VertexFinder","Finder did not succed");
+ sigma=999;
+ }
+ delete vectP0;
+ delete vectP1;
+ fSimpVert.SetXYZ(initPos);
+ fSimpVert.SetDispersion(sigma);
+ fSimpVert.SetNContributors(knacc);
+}
+//_______________________________________________________________________
+Double_t AliITSVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
+ //
+ 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];
+ return det;
+}
+//____________________________________________________________________________
+void AliITSVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t m[][3],Double_t *d){
+
+ //
+ Double_t x12=p0[0]-p1[0];
+ Double_t y12=p0[1]-p1[1];
+ Double_t z12=p0[2]-p1[2];
+ Double_t kk=x12*x12+y12*y12+z12*z12;
+ m[0][0]=2-2/kk*x12*x12;
+ m[0][1]=-2/kk*x12*y12;
+ m[0][2]=-2/kk*x12*z12;
+ m[1][0]=-2/kk*x12*y12;
+ m[1][1]=2-2/kk*y12*y12;
+ m[1][2]=-2/kk*y12*z12;
+ m[2][0]=-2/kk*x12*z12;
+ m[2][1]=-2*y12*z12;
+ m[2][2]=2-2/kk*z12*z12;
+ d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
+ d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
+ d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
+
+}
+//____________________________________________________________________________
+void AliITSVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t m[][3],Double_t *d){
+ //
+ Double_t x12=p1[0]-p0[0];
+ Double_t y12=p1[1]-p0[1];
+ Double_t z12=p1[2]-p0[2];
+
+ Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
+
+ Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
+
+ Double_t cc[3];
+ cc[0]=-x12/sigmasq[0];
+ cc[1]=-y12/sigmasq[1];
+ cc[2]=-z12/sigmasq[2];
+
+ 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;
+
+ Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
+
+ Double_t aa[3];
+ aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
+ aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
+ aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
+
+ m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
+ m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
+ m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
+
+ m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
+ m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
+ m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
+
+ m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
+ m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
+ m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
+
+ d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
+ d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
+ d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
+
+}
+//_____________________________________________________________________________
+Double_t AliITSVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
+ //
+ Double_t x12=p0[0]-p1[0];
+ Double_t y12=p0[1]-p1[1];
+ Double_t z12=p0[2]-p1[2];
+ Double_t x10=p0[0]-x0[0];
+ Double_t y10=p0[1]-x0[1];
+ Double_t z10=p0[2]-x0[2];
+ 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);
}
//---------------------------------------------------------------------------
Int_t i,j,k,step=0;
TMatrixD rv(3,1);
TMatrixD vV(3,3);
- rv(0,0) = fInitPos[0];
- rv(1,0) = fInitPos[1];
+ Double_t initPos[3];
+ fSimpVert.GetXYZ(initPos);
+ rv(0,0) = initPos[0];
+ rv(1,0) = initPos[1];
rv(2,0) = 0.;
Double_t xlStart,alpha;
Double_t rotAngle;
// get track from track array
t = (AliITStrackV2*)fTrkArray.At(k);
alpha = t->GetAlpha();
- xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
+ xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
t->PropagateTo(xlStart,0.,0.); // to vtxSeed
rotAngle = alpha;
if(alpha<0.) rotAngle += 2.*TMath::Pi();