Changed AliRunLoader::GetRunLoader() into AliRunLoader::Instance()
[u/mrichter/AliRoot.git] / ITS / AliITSFindTracksANN.C
CommitLineData
cec46807 1#if !defined(__CINT__) || defined(__MAKECINT__)
2 #include "Riostream.h"
8b462fd8 3 #include "TFile.h"
4 #include "TTree.h"
cec46807 5 #include "TStopwatch.h"
6
7 #include "AliRun.h"
cec46807 8 #include "AliRunLoader.h"
9 #include "AliTPCLoader.h"
10 #include "AliITSLoader.h"
11 #include "AliITS.h"
12 #include "AliITSgeom.h"
8b462fd8 13 #include "AliITStrackerANN.h"
cec46807 14 #include "AliESD.h"
15#endif
16
17extern AliRun *gAlice;
18
19Int_t AliITSFindTracksANN
20(Int_t nsectors = 10, // number of azimutal sectors
21 Int_t nev = 5) // number of events to process
22{
23
24 // ==================================
25 // ==== EVENT READING ===============
26 // ==================================
27
28 // Remove an already existing Run Loader
29 if (gAlice) {
33c3c91a 30 delete AliRunLoader::Instance();
cec46807 31 delete gAlice;
32 gAlice=0;
33 }
34
35 // Instance the new Run Loader
36 AliRunLoader* rl = AliRunLoader::Open("galice.root");
37 if (rl == 0x0) {
38 cerr<<"AliITSFindTracks.C : Can not open session RL=NULL"<< endl;
39 return 3;
40 }
41
42 // Instance the ITS Loader
43 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
44 if (itsl == 0x0) {
45 cerr<<"AliITSFindTracksV2.C : Can not get ITS loader"<<endl;
46 return 4;
47 }
48
49 // Load the gAlice
50 if (rl->LoadgAlice()) {
51 cerr<<"AliITSFindTracksV2.C : LoadgAlice returned error"<<endl;
52 delete rl;
53 return 3;
54 }
55
56 // Set NECESSARY conversion constant for magnetic field
57 AliKalmanTrack::SetConvConst
58 (
59 1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()
60 );
61
62 // Get the ITS data & geometry & recpoints
63 AliITS *dITS = (AliITS*)rl->GetAliRun()->GetDetector("ITS");
64 if (!dITS) {
65 cerr<<"AliITSFindClusters.C : Can not find the ITS detector !"<<endl;
66 return 6;
67 }
68 AliITSgeom *geom = dITS->GetITSgeom();
69 itsl->LoadRecPoints("read");
70
71 // ==================================
72 // ==== CURVATURE CUT DEFINITION ====
73 // ==================================
74
75 // These values define the curvature cuts for all steps
76 // within a sector.
77 // For a greater clarity, the cuts are indicated in units
78 // of transverse momentum (GeV/c) but these value have no
79 // exact physical meaning, but are useful to understand
80 // well what means a choice in the value of a certain
81 // curvature constraint
82 // NOTE: be careful to make sure that the 'ncuts' variable
83 // have the same value of the dimension of the allocated arrays
84
85 Int_t ncuts;
86 Double_t *p, *cut;
87
88 ncuts = 6;
89 p = new Double_t[6];
90 cut = new Double_t[6];
91 cut[ 0] = 0.0006;
92 cut[ 1] = 0.0010;
93 cut[ 2] = 0.0013;
94 cut[ 3] = 0.0018;
95 cut[ 4] = 0.0020;
96 cut[ 5] = 0.0030;
97
98
99 // ==========================
100 // ==== OTHER PARAMETERS ====
101 // ==========================
102
103
104 Double_t helix[5] = { 0.03, 0.02, 0.01, 0.03, 0.2 };
105 Double_t theta2D[5] = { 1.5, 1.0, 1.0, 3.0, 10.0 };
106 Double_t theta3D[5] = { 200., 200., 200., 200.0, 200.0 };
107
108
109 /*
110 Double_t helix[5] = { 0.3, 0.3, 0.3, 0.5, 0.5 };
111 Double_t theta2D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
112 Double_t theta3D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
113 */
114
115 Double_t temp = 1.0; // temperature parameter
116 Double_t var = 0.00001; // stabilization threshold
117
118 Double_t exp = 20.0; // straight-line excitator
119 Double_t gtoc = 6.0; // gain/cost contribution ratio
120
121 Double_t min = 0.4; // minimum in random activation initialization
122 Double_t max = 0.6; // maximum in random activation initialization
123 Double_t actmin = 0.55; // activation threshold for binary map conversion
124
125 // Instance the Tracker
126 AliITStrackerANN tracker(geom, 2);
127
128 // Set cuts
129 tracker.SetVertex(0.0, 0.0, 0.0);
130
131 tracker.SetCuts(ncuts, cut, theta2D, theta3D, helix);
132 tracker.SetTemperature(temp);
133 tracker.SetVariationLimit(var);
134 tracker.SetGainToCostRatio(gtoc);
135 tracker.SetWeightExponent(exp);
136 tracker.SetInitInterval(min, max);
137 tracker.SetActThreshold(actmin);
138
139 tracker.SetPolarInterval(45.0);
140
141 // Set number of events
142 if (nev > rl->GetNumberOfEvents()) nev = rl->GetNumberOfEvents();
143 Int_t rc = 0;
144
145 // Get ESD files
146 TFile *itsf=TFile::Open("AliESDits.root","RECREATE");
147 if ((!itsf)||(!itsf->IsOpen())) {
148 cerr << "Can't AliESDits.root !\n";
149 return 1;
150 }
151 TFile *tpcf=TFile::Open("AliESDtpc.root");
152 if ((!tpcf)||(!tpcf->IsOpen())) {
153 cerr<<"Can't AliESDtpc.root !\n";
154 return 1;
155 }
8b462fd8 156 AliESD* event = new AliESD;
157 TTree* esdTree = (TTree*) tpcf->Get("esdTree");
158 if (!esdTree) {
159 cerr<<"no ESD tree found !\n";
160 return 1;
161 }
162 esdTree->SetBranchAddress("ESD", &event);
cec46807 163
164 // Loop on events
cec46807 165 TStopwatch timer;
166 for (Int_t i = 0; i < nev; i++) {
cec46807 167 cerr << "Processing event number: " << i << endl;
8b462fd8 168 esdTree->GetEvent(i);
cec46807 169
170 rl->GetEvent(i);
171
172 // Get clusters from file
173 TTree *cTree=itsl->TreeR();
174 if (!cTree) {
175 cerr<<"AliITSFindTracksANN.C : NO clusters tree!" << endl;
176 return 4;
177 }
178
179 // Load clusters in tracker
180 tracker.LoadClusters(cTree);
181
182 // Create array structure and arrange points in it
183 tracker.CreateArrayStructure(nsectors);
184 tracker.ArrangePoints("its.recpoints.txt");
185 tracker.StoreOverallMatches();
186 //tracker.PrintMatches(kFALSE);
187
188 // ***************************
189 // NEURAL TRACKING OPERATION
190 // ***************************
191
192 Bool_t isStable = kFALSE;
193 Int_t i, nUnits = 0, nLinks = 0, nTracks = 0;
194 Int_t sectTracks = 0, totTracks = 0;
195
196 // tracking through sectors
197 cout << endl;
198 Int_t sector, curv;
199// nUnits = tracker.CreateNetwork(0, ncuts - 1);
200// cout << endl << nUnits << " NEURONS CREATED" << endl;
201// return;
202
203 for (sector = 0; sector < nsectors; sector++) {
204 sectTracks = 0;
205 for(curv = 0; curv < ncuts; curv++) {
206 // units creation
207 nUnits = tracker.CreateNetwork(sector, curv);
208 if (!nUnits) {
209 cout << "no neurons --> skipped" << endl;
210 continue;
211 }
212 // units connection
213 nLinks = tracker.LinkNeurons();
214 if (!nLinks) {
215 cout << "no connections --> skipped" << endl;
216 continue;
217 }
218 // neural updating
219 for (i = 0;; i++) {
220 isStable = tracker.Update();
221 if (isStable) break;
222 }
223 // tracks saving
224 tracker.CleanNetwork();
225 nTracks = tracker.StoreTracks();
226 cout << nUnits << " units, ";
227 cout << nLinks << " connections, ";
228 cout << i << " cycles --> ";
229 cout << nTracks << " tracks" << endl;
230 sectTracks += nTracks;
231 }
232 cout << "\n >>> Total tracks found in sector: " << sectTracks << endl;
233 totTracks += sectTracks;
234 }
235 tracker.ExportTracks("ITS.Neural.root");
236 cout << "\n\n--- Totally found " << totTracks << " tracks ---\n\n" << endl;
237
238 // ***************************
239
240 // End of operations: unload clusters
241 tracker.UnloadClusters();
242
cec46807 243 }
244 timer.Stop();
245 timer.Print();
246
247 // Close files & delete objects
8b462fd8 248 delete event;
cec46807 249 tpcf->Close();
250 itsf->Close();
251 delete rl;
252 return rc;
253}
254
255