1 #if !defined(__CINT__) || defined(__MAKECINT__)
5 #include "TStopwatch.h"
8 #include "AliRunLoader.h"
9 #include "AliTPCLoader.h"
10 #include "AliITSLoader.h"
12 #include "AliITSgeom.h"
13 #include "AliITStrackerANN.h"
17 extern AliRun *gAlice;
19 Int_t AliITSFindTracksANN
20 (Int_t nsectors = 10, // number of azimutal sectors
21 Int_t nev = 5) // number of events to process
24 // ==================================
25 // ==== EVENT READING ===============
26 // ==================================
28 // Remove an already existing Run Loader
30 delete AliRunLoader::Instance();
35 // Instance the new Run Loader
36 AliRunLoader* rl = AliRunLoader::Open("galice.root");
38 cerr<<"AliITSFindTracks.C : Can not open session RL=NULL"<< endl;
42 // Instance the ITS Loader
43 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
45 cerr<<"AliITSFindTracksV2.C : Can not get ITS loader"<<endl;
50 if (rl->LoadgAlice()) {
51 cerr<<"AliITSFindTracksV2.C : LoadgAlice returned error"<<endl;
56 // Set NECESSARY conversion constant for magnetic field
57 AliKalmanTrack::SetConvConst
59 1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()
62 // Get the ITS data & geometry & recpoints
63 AliITS *dITS = (AliITS*)rl->GetAliRun()->GetDetector("ITS");
65 cerr<<"AliITSFindClusters.C : Can not find the ITS detector !"<<endl;
68 AliITSgeom *geom = dITS->GetITSgeom();
69 itsl->LoadRecPoints("read");
71 // ==================================
72 // ==== CURVATURE CUT DEFINITION ====
73 // ==================================
75 // These values define the curvature cuts for all steps
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
90 cut = new Double_t[6];
99 // ==========================
100 // ==== OTHER PARAMETERS ====
101 // ==========================
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 };
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 };
115 Double_t temp = 1.0; // temperature parameter
116 Double_t var = 0.00001; // stabilization threshold
118 Double_t exp = 20.0; // straight-line excitator
119 Double_t gtoc = 6.0; // gain/cost contribution ratio
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
125 // Instance the Tracker
126 AliITStrackerANN tracker(geom, 2);
129 tracker.SetVertex(0.0, 0.0, 0.0);
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);
139 tracker.SetPolarInterval(45.0);
141 // Set number of events
142 if (nev > rl->GetNumberOfEvents()) nev = rl->GetNumberOfEvents();
146 TFile *itsf=TFile::Open("AliESDits.root","RECREATE");
147 if ((!itsf)||(!itsf->IsOpen())) {
148 cerr << "Can't AliESDits.root !\n";
151 TFile *tpcf=TFile::Open("AliESDtpc.root");
152 if ((!tpcf)||(!tpcf->IsOpen())) {
153 cerr<<"Can't AliESDtpc.root !\n";
156 AliESD* event = new AliESD;
157 TTree* esdTree = (TTree*) tpcf->Get("esdTree");
159 cerr<<"no ESD tree found !\n";
162 esdTree->SetBranchAddress("ESD", &event);
166 for (Int_t i = 0; i < nev; i++) {
167 cerr << "Processing event number: " << i << endl;
168 esdTree->GetEvent(i);
172 // Get clusters from file
173 TTree *cTree=itsl->TreeR();
175 cerr<<"AliITSFindTracksANN.C : NO clusters tree!" << endl;
179 // Load clusters in tracker
180 tracker.LoadClusters(cTree);
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);
188 // ***************************
189 // NEURAL TRACKING OPERATION
190 // ***************************
192 Bool_t isStable = kFALSE;
193 Int_t i, nUnits = 0, nLinks = 0, nTracks = 0;
194 Int_t sectTracks = 0, totTracks = 0;
196 // tracking through sectors
199 // nUnits = tracker.CreateNetwork(0, ncuts - 1);
200 // cout << endl << nUnits << " NEURONS CREATED" << endl;
203 for (sector = 0; sector < nsectors; sector++) {
205 for(curv = 0; curv < ncuts; curv++) {
207 nUnits = tracker.CreateNetwork(sector, curv);
209 cout << "no neurons --> skipped" << endl;
213 nLinks = tracker.LinkNeurons();
215 cout << "no connections --> skipped" << endl;
220 isStable = tracker.Update();
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;
232 cout << "\n >>> Total tracks found in sector: " << sectTracks << endl;
233 totTracks += sectTracks;
235 tracker.ExportTracks("ITS.Neural.root");
236 cout << "\n\n--- Totally found " << totTracks << " tracks ---\n\n" << endl;
238 // ***************************
240 // End of operations: unload clusters
241 tracker.UnloadClusters();
247 // Close files & delete objects