Transition to NewIO
[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 <AliITSLoader.h>
12
13 #endif
14
15 void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t NumbofEv=1, const char *filename="galice.root") {
16
17   Int_t evNumber2 = evNumber1+NumbofEv;
18   
19   if (gClassTable->GetID("AliRun") < 0) {
20     gROOT->Macro("loadlibs.C");
21   } else {
22     if(gAlice){
23       delete gAlice->GetRunLoader();
24       delete gAlice;
25       gAlice=0;
26     }
27   }
28
29   
30   AliRunLoader* rl = AliRunLoader::Open("galice.root");
31   if (rl == 0x0){
32     ::Error("AliITSFindPrimaryVertex.C","Can not open session RL=NULL");
33     return;
34   }
35  
36   Int_t retval = rl->LoadHeader();
37   if (retval){
38     cerr<<"AliITSFindPrimaryVertex.C : LoadHeader returned error"<<endl;
39     return;
40   }
41   retval = rl->LoadKinematics();
42   if (retval){
43     cerr<<"AliITSFindPrimaryVertex.C : LoadKinematics returned error"<<endl;
44     return;
45   }
46
47   // Open output file for vertices (default name: ITS.Vertex.root 
48   // and Create vertexer
49   AliITSVertexerIons *vertexer = new AliITSVertexerIons("default");
50   AliITSVertex *V;
51   //   Loop over events 
52   //
53  
54   for (int nev=evNumber1; nev< evNumber2; nev++) {
55     cout<<"=============================================================\n";
56     cout<<" Processing event "<<nev<<endl;
57     cout<<"=============================================================\n";
58     rl->GetEvent(nev);
59     cout << "nev         " << nev <<endl;
60     // The true Z coord. is fetched for comparison
61     AliHeader *header = rl->GetHeader();
62     AliGenEventHeader* genEventHeader = header->GenEventHeader();
63     TArrayF primaryVertex(3);
64     genEventHeader->PrimaryVertex(primaryVertex);
65   
66     TStopwatch timer;
67     timer.Start();
68
69     V=vertexer->FindVertexForCurrentEvent(nev);
70     if(V){
71       Double_t pos[3];
72       for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
73       V->SetTruePos(pos);
74     }
75     timer.Stop();
76     timer.Print();
77
78     cout << endl << "Xv = " << V->GetXv() << " cm" << endl;
79     cout << "X resolution = " << V->GetXRes()*10000 << " microns"  << endl;
80     cout << "Signal/Noise for X = " << V->GetXSNR() << endl;
81     cout << endl << "Yv = " << V->GetYv() << " cm"  << endl;
82     cout << "Y resolution = " << V->GetYRes()*10000 << " microns"  << endl;
83     cout << "Signal/Noise for Y = " << V->GetYSNR() << endl;
84     cout << endl << "Zv = " << V->GetZv() << " cm" << endl;
85     cout << "Z Resolution = " << V->GetZRes()*10000 << " microns" << endl;
86     cout << "Signal/Noise for Z = " << V->GetZSNR() <<endl;
87                  
88     vertexer->WriteCurrentVertex(); 
89   } 
90   
91
92 }
93