The present commit corresponds to an important change in the way the
[u/mrichter/AliRoot.git] / ITS / AliITSFindTracksANN.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2   #include "Riostream.h"
3   #include "TFile.h"
4   #include "TTree.h"
5   #include "TStopwatch.h"
6
7   #include "AliRun.h"
8   #include "AliRunLoader.h"
9   #include "AliTPCLoader.h"
10   #include "AliITSLoader.h"
11   #include "AliITS.h"
12   #include "AliITSgeom.h"
13   #include "AliITStrackerANN.h"
14   #include "AliESD.h"
15 #endif
16
17 extern AliRun *gAlice;
18
19 Int_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) {
30                 delete AliRunLoader::GetRunLoader();
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         }
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);
163         
164         // Loop on events
165         TStopwatch timer; 
166         for (Int_t i = 0; i < nev; i++) {
167                 cerr << "Processing event number: " << i << endl;
168                 esdTree->GetEvent(i);
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                 
243         }
244         timer.Stop(); 
245         timer.Print();
246         
247         // Close files & delete objects
248         delete event;
249         tpcf->Close();
250         itsf->Close();
251         delete rl;
252         return rc;
253 }
254
255