Reconstruction of particle multiplicity (T.Virgili, C.Jorgensen)
[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>
7
8#include "AliRunLoader.h"
9#include "AliESD.h"
10#include "AliRun.h"
11
12#include "AliITS.h"
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 // #########################################################
31 // setup galice and runloader
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 }
45 runLoader->LoadgAlice();
46
47 gAlice = runLoader->GetAliRun();
48 runLoader->LoadKinematics();
49 runLoader->LoadHeader();
50
51 // #########################################################
52 // open esd file and get the tree
53
54 // close it first to avoid memory leak
55 if (esdFile)
56 if (esdFile->IsOpen())
57 esdFile->Close();
58
59 sprintf(str,"%s/AliESDs.root",dir);
60 esdFile = TFile::Open(str);
61 esdTree = (TTree*)esdFile->Get("esdTree");
62 TBranch * esdBranch = esdTree->GetBranch("ESD");
63 esdBranch->SetAddress(&esd);
64
65
66 // #########################################################
67 // setup its stuff
68
69 AliITS* its=(AliITS*)runLoader->GetAliRun()->GetDetector("ITS");
70 if (!its) {
71 cout << " Can't get the ITS!" << endl;
72 return ;
73 }
74 AliITSgeom* itsGeo=its->GetITSgeom();
75 if (!itsGeo) {
76 cout << " Can't get the ITS geometry!" << endl;
77 return ;
78 }
79 AliITSLoader* itsLoader = (AliITSLoader*)runLoader->GetLoader("ITSLoader");
80 if (!itsLoader) {
81 cout << " Can't get the ITS loader!" << endl;
82 return ;
83 }
84 itsLoader->LoadRecPoints("read");
85
86 // #########################################################
87 AliITSMultReconstructor* multReco = new AliITSMultReconstructor();
88 multReco->SetGeometry(itsGeo);
89
90 // #########################################################
91 // getting number of events
92
93 Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
94 Int_t nESDEvents = esdBranch->GetEntries();
95
96 if (nEvents!=nESDEvents) {
97 cout << " Different number of events from runloader and esdtree!!!"
98 << nEvents << " / " << nESDEvents << endl;
99 return;
100 }
101
102 // ########################################################
103 // loop over number of events
104 cout << nEvents << " event(s) found in the file set" << endl;
105 for(Int_t i=0; i<nEvents; i++) {
106
107 cout << "-------------------------" << endl << " event# " << i << endl;
108
109 runLoader->GetEvent(i);
110 esdBranch->GetEntry(i);
111
112 // ########################################################
113 // get the EDS vertex
114 const AliESDVertex* vtxESD = esd->GetVertex();
115 Double_t vtx[3];
116 vtxESD->GetXYZ(vtx);
117 Float_t esdVtx[3];
118 esdVtx[0] = vtx[0];
119 esdVtx[1] = vtx[1];
120 esdVtx[2] = vtx[2];
121
122 ///#########################################################
123 // get ITS clusters
124 TTree* itsClusterTree = itsLoader->TreeR();
125 if (!itsClusterTree) {
126 cerr<< " Can't get the ITS cluster tree !\n";
127 return;
128 }
129 multReco->SetHistOn(kTRUE);
130 multReco->Reconstruct(itsClusterTree, esdVtx, esdVtx);
131
132
133 for (Int_t t=0; t<multReco->GetNTracklets(); t++) {
134
135 cout << " tracklet " << t
136 << " , theta = " << multReco->GetTracklet(t)[0]
137 << " , phi = " << multReco->GetTracklet(t)[1] << endl;
138 }
139
140 }
141
142 TFile* fout = new TFile("out.root","RECREATE");
143
144 multReco->SaveHists();
145 fout->Write();
146 fout->Close();
147
148
149}