Changes related to the extraction of the V0 finder into a separate class (A. Dainese...
[u/mrichter/AliRoot.git] / ITS / AliITSFindPrimaryVertex.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TROOT.h>
3 #include <TClassTable.h>
4 #include <TStopwatch.h>
5 #include <TFile.h>
6 #include <TTree.h>
7 #include <Riostream.h>
8 #include <AliRun.h>
9 #include <AliHeader.h>
10 #include <AliGenEventHeader.h>
11 #include <AliITSVertexerIons.h>
12 #include <AliRunLoader.h>
13 #include <AliITSVertexerIons.h>
14 #include <AliITSLoader.h>
15 #include <AliGenHijingEventHeader.h>
16 #include <unistd.h>
17
18 #endif
19
20 void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t NumbofEv=1, const char *filename="galice.root") {
21
22   Int_t evNumber2 = evNumber1+NumbofEv;
23   
24   if (gClassTable->GetID("AliRun") < 0) {
25     gROOT->Macro("loadlibs.C");
26   } else {
27     if(gAlice){
28       delete AliRunLoader::Instance();
29       delete gAlice;
30       gAlice=0;
31     }
32   }
33   
34   AliRunLoader* rl = AliRunLoader::Open("galice.root");
35   if (rl == 0x0){
36     ::Error("AliITSFindPrimaryVertex.C","Can not open session RL=NULL");
37     return;
38   }
39   
40   // Open output file for vertices (default name: ITS.Vertex.root 
41   // and Create vertexer
42
43   AliITSVertexerIons *vertexer = new AliITSVertexerIons();
44   vertexer->Init("default");
45   //vertexer->SetDebug(1);
46   
47   AliESDVertex *V;
48   //   Loop over events 
49    
50   AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
51   itsloader->LoadRecPoints("read");
52
53   for (int nev=evNumber1; nev< evNumber2; nev++) {
54     cout<<"=============================================================\n";
55     cout<<" Processing event "<<nev<<endl;
56     cout<<"=============================================================\n";
57     rl->GetEvent(nev);
58     AliHeader *header = rl->GetHeader();
59     AliGenEventHeader* genEventHeader = header->GenEventHeader();
60     TArrayF primaryVertex(3);
61     genEventHeader->PrimaryVertex(primaryVertex);
62     
63     AliGenHijingEventHeader* hijingHeader = (AliGenHijingEventHeader*)  genEventHeader;
64     Float_t b = hijingHeader->ImpactParameter();   
65     cout << "Impact parameter = " << b << " fm" << endl;
66
67     TStopwatch timer;
68     timer.Start();
69
70     TTree* cltree = itsloader->TreeR();
71     V=vertexer->FindVertexForCurrentEvent(cltree);
72
73     TVector3 vtrue(primaryVertex[0],primaryVertex[1],primaryVertex[2]);
74     TVector3 vfound(V->GetXv(),V->GetYv(),V->GetZv());
75     TVector3 dif=vtrue-vfound;
76     cout << "True vertex coordinates (cm) = " << vtrue.X() << " " << vtrue.Y() << " " << vtrue.Z() << endl;
77     cout << "Found vertex coordinates  (cm) = " << vfound.X() << " " << vfound.Y() << " " << vfound.Z() << endl;    cout << "Difference true - found (cm) = " << dif.Mag() << " " << dif.X() << " " << dif.Y() << " " << dif.Z() << endl;
78     
79     timer.Stop();
80     timer.Print();
81     
82     vertexer->WriteCurrentVertex(); 
83   } 
84   
85
86 }
87