Add class and function docs.
[u/mrichter/AliRoot.git] / ITS / testITSMultReco.C
CommitLineData
ac903f1b 1#if !defined(__CINT__) || defined(__MAKECINT__)
2
3#include <Riostream.h>
4#include <TFile.h>
5#include <TTree.h>
6#include <TBranch.h>
968e8539 7#include <TGeoManager.h>
ac903f1b 8
9#include "AliRunLoader.h"
10#include "AliESD.h"
11#include "AliRun.h"
12
ac903f1b 13#include "AliITSgeom.h"
14#include "AliITSLoader.h"
15#include "AliITSMultReconstructor.h"
16
17#endif
18
19void testITSMultReco(Char_t* dir = ".") {
20
21 Char_t str[256];
22
23 // ########################################################
24 // defining pointers
25 AliRunLoader* runLoader;
26 TFile* esdFile = 0;
27 TTree* esdTree = 0;
28 AliESD* esd = 0;
29
30 // #########################################################
968e8539 31 // get runloader
ac903f1b 32
33 if (gAlice) {
34 delete gAlice->GetRunLoader();
35 delete gAlice;
36 gAlice=0;
37 }
38
39 sprintf(str,"%s/galice.root",dir);
40 runLoader = AliRunLoader::Open(str);
41 if (runLoader == 0x0) {
42 cout << "Can not open session"<<endl;
43 return;
44 }
968e8539 45 // get geometry (here geometry.root is used, change it if needed)
46 if (!gGeoManager) {
47 sprintf(str,"%s/geometry.root",dir);
48 TGeoManager::Import(str);
49 if(!gGeoManager) {
50 cout << "Can not access the geometry file"<<endl;
51 return;
52 }
53 }
ac903f1b 54
ac903f1b 55
56 // #########################################################
57 // open esd file and get the tree
58
59 // close it first to avoid memory leak
60 if (esdFile)
61 if (esdFile->IsOpen())
62 esdFile->Close();
63
64 sprintf(str,"%s/AliESDs.root",dir);
65 esdFile = TFile::Open(str);
66 esdTree = (TTree*)esdFile->Get("esdTree");
67 TBranch * esdBranch = esdTree->GetBranch("ESD");
68 esdBranch->SetAddress(&esd);
69
70
71 // #########################################################
72 // setup its stuff
73
ac903f1b 74 AliITSLoader* itsLoader = (AliITSLoader*)runLoader->GetLoader("ITSLoader");
75 if (!itsLoader) {
76 cout << " Can't get the ITS loader!" << endl;
77 return ;
78 }
968e8539 79 AliITSgeom* itsGeo=itsLoader->GetITSgeom();
ac903f1b 80 itsLoader->LoadRecPoints("read");
81
82 // #########################################################
83 AliITSMultReconstructor* multReco = new AliITSMultReconstructor();
84 multReco->SetGeometry(itsGeo);
85
86 // #########################################################
87 // getting number of events
88
89 Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
90 Int_t nESDEvents = esdBranch->GetEntries();
91
92 if (nEvents!=nESDEvents) {
93 cout << " Different number of events from runloader and esdtree!!!"
94 << nEvents << " / " << nESDEvents << endl;
95 return;
96 }
97
98 // ########################################################
99 // loop over number of events
100 cout << nEvents << " event(s) found in the file set" << endl;
101 for(Int_t i=0; i<nEvents; i++) {
102
103 cout << "-------------------------" << endl << " event# " << i << endl;
104
105 runLoader->GetEvent(i);
106 esdBranch->GetEntry(i);
107
108 // ########################################################
109 // get the EDS vertex
110 const AliESDVertex* vtxESD = esd->GetVertex();
111 Double_t vtx[3];
112 vtxESD->GetXYZ(vtx);
113 Float_t esdVtx[3];
114 esdVtx[0] = vtx[0];
115 esdVtx[1] = vtx[1];
116 esdVtx[2] = vtx[2];
117
118 ///#########################################################
119 // get ITS clusters
120 TTree* itsClusterTree = itsLoader->TreeR();
121 if (!itsClusterTree) {
122 cerr<< " Can't get the ITS cluster tree !\n";
123 return;
124 }
125 multReco->SetHistOn(kTRUE);
126 multReco->Reconstruct(itsClusterTree, esdVtx, esdVtx);
127
968e8539 128 cout <<" >>>> Number of tracklets: "<<multReco->GetNTracklets()<<endl;
ac903f1b 129 for (Int_t t=0; t<multReco->GetNTracklets(); t++) {
130
131 cout << " tracklet " << t
132 << " , theta = " << multReco->GetTracklet(t)[0]
968e8539 133 << " , phi = " << multReco->GetTracklet(t)[1]
134 << " , DeltaPhi = " << multReco->GetTracklet(t)[2]<< endl;
135 }
136 cout <<" >>>> Number of single layer 1 clusters: "<<multReco->GetNSingleClusters()<<endl;
137 for (Int_t t=0; t<multReco->GetNSingleClusters(); t++) {
138
139 cout << " cluster " << t
140 << " , theta = " << multReco->GetCluster(t)[0]
141 << " , phi = " << multReco->GetCluster(t)[1] << endl;
ac903f1b 142 }
143
144 }
145
146 TFile* fout = new TFile("out.root","RECREATE");
147
148 multReco->SaveHists();
149 fout->Write();
150 fout->Close();
151
152
153}