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