Updated version of the mutiplicity reconstruction with tracklets (Massimo)
[u/mrichter/AliRoot.git] / ITS / testITSMultReco.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2
3 #include <Riostream.h>
4 #include <TFile.h>
5 #include <TTree.h>
6 #include <TBranch.h>
7 #include <TGeoManager.h>
8
9 #include "AliRunLoader.h"
10 #include "AliESD.h"
11 #include "AliRun.h"
12
13 #include "AliITSgeom.h"
14 #include "AliITSLoader.h"
15 #include "AliITSMultReconstructor.h"
16
17 #endif
18
19 void 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   // get 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   // 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   }
54
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
74   AliITSLoader* itsLoader = (AliITSLoader*)runLoader->GetLoader("ITSLoader");
75   if (!itsLoader) {
76     cout << " Can't get the ITS loader!" << endl;
77     return ;
78   }
79   AliITSgeom* itsGeo=itsLoader->GetITSgeom();
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
128     cout <<" >>>> Number of tracklets: "<<multReco->GetNTracklets()<<endl;     
129     for (Int_t t=0; t<multReco->GetNTracklets(); t++) {
130       
131       cout << "  tracklet " << t 
132            << " , theta = " << multReco->GetTracklet(t)[0]
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; 
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 }