New vertex finders optimized for 3 prongs secondary vertices (AliITSVertexerTracks...
[u/mrichter/AliRoot.git] / ITS / AliITSSecondaryVertices.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 //-- --- standard headers------------- 
3 #include <Riostream.h>
4 //--------Root headers ---------------
5 #include <TSystem.h>
6 #include <TFile.h>
7 #include <TStopwatch.h>
8 #include <TObject.h>
9 #include <TTree.h>
10 #include <TChain.h>
11 #include <TF1.h>
12 #include <TLegend.h>
13 #include <TLegendEntry.h>
14 #include <TH1F.h>
15 #include <TH2F.h>
16 #include <TCut.h>
17 #include <TParticle.h>
18 #include <TCanvas.h>
19 #include <TStyle.h>
20 #include <TGraph.h>
21 //----- AliRoot headers ---------------
22 //#include "alles.h"
23 #include "AliRun.h"
24 #include "AliMagF.h"
25 #include "AliKalmanTrack.h"
26 #include "AliESDVertex.h"
27 #include "AliITSVertexer.h"
28 #include "AliITSVertexerTracks.h"
29 #include <AliHeader.h>
30 #include <AliGenEventHeader.h>
31 #include "AliRun.h"
32 #include "AliRunLoader.h"
33 #include "AliITStrackV2.h"
34 #include "AliITSSimpleVertex.h"
35 #include "AliITSStrLine.h"
36 #include "AliITSLoader.h"
37 #include <AliESD.h>
38 #include <AliStack.h>
39 //-------------------------------------
40 #endif  
41
42
43 void AliITSSecondaryVertices(Int_t ntracks, Int_t *pos) {
44  
45   // example of usige AliITSVertexerTracks to get a (secodnary) vertex from a set of tracks.
46
47   Double_t field = 0.4;
48  
49   AliRunLoader *rl = AliRunLoader::Open("galice.root");
50   if (!rl) {
51     cerr<<"Can't start session !\n";
52   }
53
54   rl->LoadgAlice();
55   if (rl->GetAliRun()){
56     AliMagF * magf = gAlice->Field();
57     AliKalmanTrack:: SetFieldMap(magf);
58     AliKalmanTrack::SetUniformFieldTracking();
59   }else {
60     cerr<<"Can't get AliRun !\n";
61     
62   }
63
64   field=rl->GetAliRun()->Field()->SolenoidField()/10.;
65   rl->LoadHeader();
66  
67   
68   TFile *ef=TFile::Open("AliESDs.root");
69   if (!ef || !ef->IsOpen()) {cerr<<"Can't AliESDs.root !n";}
70   AliESD* event = new AliESD;
71   TTree* tree = (TTree*) ef->Get("esdTree");
72   if (!tree) {cerr<<"no ESD tree foundn";};
73   tree->SetBranchAddress("ESD", &event);
74
75
76   Int_t e=0;
77   tree->GetEvent(e);
78   Int_t ntrk=event->GetNumberOfTracks();
79   cout<<"Number of ESD tracks : "<<ntrk<<endl; 
80
81   // Open input and output files
82   TFile *outFile = TFile::Open("AliITSVertices.root","recreate");
83     
84
85   TTree *trkTree = new TTree("TreeT_ITS","its tracks");
86   AliITStrackV2 *itstrack = 0;
87   trkTree->Branch("tracks","AliITStrackV2",&itstrack,ntrk,0);
88   Int_t pipe=0;
89   for(Int_t i=0; i<ntrk; i++) {
90     AliESDtrack *esdTrack = (AliESDtrack*)event->GetTrack(i);
91     UInt_t status=esdTrack->GetStatus();
92     UInt_t flags=AliESDtrack::kTPCin|AliESDtrack::kITSin;
93     if ((status&AliESDtrack::kTPCin)==0)continue;
94     if ((status&AliESDtrack::kITSin)==0)continue;    
95     if ((status&AliESDtrack::kITSrefit)==0) continue;
96     if ((status&flags)==status) pipe=1;
97
98     itstrack = new AliITStrackV2(*esdTrack);
99     if (pipe!=0) {
100       itstrack->PropagateTo(3.,0.0028,65.19);
101       itstrack->PropagateToVertex();  // This is not "exactly" the vertex 
102     }
103     Int_t nclus=itstrack->GetNumberOfClusters();
104
105     if(nclus<6)continue;
106     trkTree->Fill();
107     itstrack = 0;
108     // delete esdTrack; 
109   }
110   // Primary vertex
111   Double_t vtx[3];
112   event->GetVertex()->GetXYZ(vtx);
113   cout<<"Primary Vertex:"<<endl;
114   event->GetVertex()->PrintStatus();
115   cout<<"\n\n\n"<<endl;
116   
117
118   // Create vertexer
119   AliITSVertexerTracks *vertexer = new AliITSVertexerTracks();
120   vertexer->SetDCAcut(0);
121   vertexer->SetDebug(0);
122
123
124   Double_t vertF1[3],vertF2[3],vertF3[3],vertF4[3],vertF5[3];
125   Int_t ncombi1,ncombi2,ncombi3,ncombi4,ncombi5;
126   Double_t sig1,sig2,sig3,sig4,sig5; 
127   AliITSSimpleVertex *vert=0;
128
129   
130   // Calculate vertex from tracks in pos
131   // BEWARE: the method VertexForSelectedTracks returns the address of the data member fSimpVert 
132   // which is always the same for all calls to VertexForSelectedTracks.
133
134   cout<<"\nStraight Line Min Dist finder with weights"<<endl;
135   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,1);
136   vert->Print();
137   vert->GetXYZ(vertF1);
138   ncombi1=vert->GetNContributors();
139   sig1=vert->GetDispersion();
140
141
142   cout<<"Straight Line Min Dist finder"<<endl;
143   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,2);   
144   vert->Print();
145   vert->GetXYZ(vertF2);
146   ncombi2=vert->GetNContributors();
147   sig2=vert->GetDispersion();
148   
149   cout<<"Helix finder"<<endl;
150   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,3);   
151   vert->Print();
152   vert->GetXYZ(vertF3);
153   ncombi3=vert->GetNContributors();
154   sig3=vert->GetDispersion();
155  
156   cout<<"Straight Line finder with weights"<<endl;
157   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,4);   
158   vert->Print();
159   vert->GetXYZ(vertF4);
160   ncombi4=vert->GetNContributors();
161   sig4=vert->GetDispersion();
162   
163   cout<<"Straight Line finder with weights"<<endl;
164   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,5);
165   vert->Print();
166   vert->GetXYZ(vertF5);
167   ncombi5=vert->GetNContributors();
168   sig5=vert->GetDispersion();
169   
170
171   delete vertexer;
172
173   outFile->Close();
174   return;
175 }