1 #if !defined(__CINT__) || defined(__MAKECINT__)
5 #include "TStopwatch.h"
9 #include "AliRunLoader.h"
10 #include "AliTPCLoader.h"
11 #include "AliITSLoader.h"
13 #include "AliITSgeom.h"
14 #include "AliITStrackerANN.h"
18 extern AliRun *gAlice;
20 Int_t AliITSFindTracksANN
21 (Int_t nsectors = 10, // number of azimutal sectors
22 Int_t nev = 5) // number of events to process
25 // ==================================
26 // ==== EVENT READING ===============
27 // ==================================
29 // Remove an already existing Run Loader
31 delete gAlice->GetRunLoader();
36 // Instance the new Run Loader
37 AliRunLoader* rl = AliRunLoader::Open("galice.root");
39 cerr<<"AliITSFindTracks.C : Can not open session RL=NULL"<< endl;
43 // Instance the ITS Loader
44 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
46 cerr<<"AliITSFindTracksV2.C : Can not get ITS loader"<<endl;
51 if (rl->LoadgAlice()) {
52 cerr<<"AliITSFindTracksV2.C : LoadgAlice returned error"<<endl;
57 // Set NECESSARY conversion constant for magnetic field
58 AliKalmanTrack::SetConvConst
60 1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()
63 // Get the ITS data & geometry & recpoints
64 AliITS *dITS = (AliITS*)rl->GetAliRun()->GetDetector("ITS");
66 cerr<<"AliITSFindClusters.C : Can not find the ITS detector !"<<endl;
69 AliITSgeom *geom = dITS->GetITSgeom();
70 itsl->LoadRecPoints("read");
72 // ==================================
73 // ==== CURVATURE CUT DEFINITION ====
74 // ==================================
76 // These values define the curvature cuts for all steps
78 // For a greater clarity, the cuts are indicated in units
79 // of transverse momentum (GeV/c) but these value have no
80 // exact physical meaning, but are useful to understand
81 // well what means a choice in the value of a certain
82 // curvature constraint
83 // NOTE: be careful to make sure that the 'ncuts' variable
84 // have the same value of the dimension of the allocated arrays
91 cut = new Double_t[6];
100 // ==========================
101 // ==== OTHER PARAMETERS ====
102 // ==========================
105 Double_t helix[5] = { 0.03, 0.02, 0.01, 0.03, 0.2 };
106 Double_t theta2D[5] = { 1.5, 1.0, 1.0, 3.0, 10.0 };
107 Double_t theta3D[5] = { 200., 200., 200., 200.0, 200.0 };
111 Double_t helix[5] = { 0.3, 0.3, 0.3, 0.5, 0.5 };
112 Double_t theta2D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
113 Double_t theta3D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
116 Double_t temp = 1.0; // temperature parameter
117 Double_t var = 0.00001; // stabilization threshold
119 Double_t exp = 20.0; // straight-line excitator
120 Double_t gtoc = 6.0; // gain/cost contribution ratio
122 Double_t min = 0.4; // minimum in random activation initialization
123 Double_t max = 0.6; // maximum in random activation initialization
124 Double_t actmin = 0.55; // activation threshold for binary map conversion
126 // Instance the Tracker
127 AliITStrackerANN tracker(geom, 2);
130 tracker.SetVertex(0.0, 0.0, 0.0);
132 tracker.SetCuts(ncuts, cut, theta2D, theta3D, helix);
133 tracker.SetTemperature(temp);
134 tracker.SetVariationLimit(var);
135 tracker.SetGainToCostRatio(gtoc);
136 tracker.SetWeightExponent(exp);
137 tracker.SetInitInterval(min, max);
138 tracker.SetActThreshold(actmin);
140 tracker.SetPolarInterval(45.0);
142 // Set number of events
143 if (nev > rl->GetNumberOfEvents()) nev = rl->GetNumberOfEvents();
147 TFile *itsf=TFile::Open("AliESDits.root","RECREATE");
148 if ((!itsf)||(!itsf->IsOpen())) {
149 cerr << "Can't AliESDits.root !\n";
152 TFile *tpcf=TFile::Open("AliESDtpc.root");
153 if ((!tpcf)||(!tpcf->IsOpen())) {
154 cerr<<"Can't AliESDtpc.root !\n";
157 AliESD* event = new AliESD;
158 TTree* esdTree = (TTree*) tpcf->Get("esdTree");
160 cerr<<"no ESD tree found !\n";
163 esdTree->SetBranchAddress("ESD", &event);
167 for (Int_t i = 0; i < nev; i++) {
168 cerr << "Processing event number: " << i << endl;
169 esdTree->GetEvent(i);
173 // Get clusters from file
174 TTree *cTree=itsl->TreeR();
176 cerr<<"AliITSFindTracksANN.C : NO clusters tree!" << endl;
180 // Load clusters in tracker
181 tracker.LoadClusters(cTree);
183 // Create array structure and arrange points in it
184 tracker.CreateArrayStructure(nsectors);
185 tracker.ArrangePoints("its.recpoints.txt");
186 tracker.StoreOverallMatches();
187 //tracker.PrintMatches(kFALSE);
189 // ***************************
190 // NEURAL TRACKING OPERATION
191 // ***************************
193 Bool_t isStable = kFALSE;
194 Int_t i, nUnits = 0, nLinks = 0, nTracks = 0;
195 Int_t sectTracks = 0, totTracks = 0;
197 // tracking through sectors
200 // nUnits = tracker.CreateNetwork(0, ncuts - 1);
201 // cout << endl << nUnits << " NEURONS CREATED" << endl;
204 for (sector = 0; sector < nsectors; sector++) {
206 for(curv = 0; curv < ncuts; curv++) {
208 nUnits = tracker.CreateNetwork(sector, curv);
210 cout << "no neurons --> skipped" << endl;
214 nLinks = tracker.LinkNeurons();
216 cout << "no connections --> skipped" << endl;
221 isStable = tracker.Update();
225 tracker.CleanNetwork();
226 nTracks = tracker.StoreTracks();
227 cout << nUnits << " units, ";
228 cout << nLinks << " connections, ";
229 cout << i << " cycles --> ";
230 cout << nTracks << " tracks" << endl;
231 sectTracks += nTracks;
233 cout << "\n >>> Total tracks found in sector: " << sectTracks << endl;
234 totTracks += sectTracks;
236 tracker.ExportTracks("ITS.Neural.root");
237 cout << "\n\n--- Totally found " << totTracks << " tracks ---\n\n" << endl;
239 // ***************************
241 // End of operations: unload clusters
242 tracker.UnloadClusters();
248 // Close files & delete objects