GetPythiaParticleName removed.
[u/mrichter/AliRoot.git] / ITS / AliITSneuralTracker.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: Alberto Pulvirenti.                                            *
5  *                                                                        *
6  * Permission to use, copy, modify and distribute this software and its   *
7  * documentation strictly for non-commercial purposes is hereby granted   *
8  * without fee, provided that the above copyright notice appears in all   *
9  * copies and that both the copyright notice and this permission notice   *
10  * appear in the supporting documentation. The authors make no claims     *
11  * about the suitability of this software for any purpose. It is          *
12  * provided "as is" without express or implied warranty.                  *
13  *                                                                        *
14  * track finder for ITS stand-alone with neural network algorithm         *
15  * This class defines a Hopfield MFT neural network simulation            *
16  * which reads all recpoints from an event and produces a tree with       *
17  * the points belonging to recognized tracks                              *
18  * the TTree obtained as the output file will be saved in a .root file    * 
19  * and it must be read by groups of six entries at a time                 *
20  **************************************************************************/
21
22 #include <Riostream.h>
23 #include <Riostream.h>
24 #include <stdlib.h>
25
26 #include <TROOT.h>
27 #include <TFile.h>
28 #include <TObjArray.h>
29 #include <TTree.h>
30 #include <TMath.h>
31 #include <TRandom.h>
32 #include <TVector3.h>
33
34 #include "AliRun.h"
35 #include "AliITS.h"
36 #include "AliITSgeom.h"
37 #include "AliITSgeomMatrix.h"
38 #include "AliITSRecPoint.h"
39 #include "AliITSglobalRecPoint.h"
40 #include "AliITSneuralTrack.h"
41 #include "AliITSneuralTracker.h"
42
43 ClassImp(AliITSneuralTracker)
44
45
46 // ====================================================================================================
47 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuralTracker <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
48 // ====================================================================================================
49
50 AliITSneuralTracker::AliITSneuralTracker() 
51 {
52         //
53         // Default constructor (not for use)
54         //
55         Int_t i = 0;
56         fCurvCut = 0;
57         for (i = 0; i < 6; i++) {
58                 fPoints[i] = 0;
59                 if (i < 5) fNeurons[i] = 0;
60         }
61         fTracks = 0;
62 }
63 //---------------------------------------------------------------------------------------------------------
64 AliITSneuralTracker::AliITSneuralTracker (Int_t nsecs, Int_t ncurv, Double_t *curv, Double_t theta)
65
66         // Argumented constructor
67         // ----------------------
68         // gets the number of azymuthal sectors,
69         // the number of curvature cut steps and their cut values, 
70         // and then the polar angle cut.
71         // Be careful not to put 'ncurv' greater than the dimension of the 
72         // 'curv' array (--> segmentation fault)
73         
74         fSectorNum   = nsecs;
75         fSectorWidth = 2.0 * TMath::Pi()/(Double_t)nsecs;
76                 
77         fCurvNum = ncurv;
78         fCurvCut = new Double_t[ncurv];
79         Int_t i;
80         for (i = 0; i < ncurv; i++) fCurvCut[i] = curv[i];
81         
82         fThetaCut = theta;
83         fThetaNum = (Int_t)(TMath::Pi() / theta) + 1;
84         if (fThetaNum < 1) fThetaNum = 1;
85                 
86         cout << "\nClass AliITSneuralTracker invoked with " << fSectorNum << " phi secs.";
87         cout << " and " << fThetaNum << " theta secs.\n" << endl;
88         
89         Int_t k, j;
90         for (k = 0; k < 6; k++) {
91                 fPoints[k] = new TObjArray**[fSectorNum];
92                 for (i = 0; i < fSectorNum; i++) {
93                         fPoints[k][i] = new TObjArray*[fThetaNum];
94                         for (j = 0; j < fThetaNum; j++) fPoints[k][i][j] = new TObjArray(0);
95                 }
96                 if (k < 6) fNeurons[k] = new TObjArray(0);
97         }
98         
99         fTracks = new TObjArray(1);
100 }
101 //---------------------------------------------------------------------------------------------------------
102 AliITSneuralTracker::~AliITSneuralTracker() 
103 {
104         // Destructor
105         // ----------
106         // clears the TObjArrays and all heap objects
107         
108         Int_t lay, sec, k;
109         delete [] fCurvCut;
110         
111         for (lay = 0; lay < 6; lay++) {
112                 for (sec = 0; sec < fSectorNum; sec++) {
113                         for (k = 0; k < fThetaNum; k++) {
114                                 fPoints[lay][sec][k]->Delete();
115                         }
116                 }
117                 delete [] fPoints[lay];
118                 if (lay < 5) {
119                         fNeurons[lay]->Delete();
120                         delete fNeurons[lay];
121                 }
122         }
123         
124         fTracks->Delete();
125 }
126 //---------------------------------------------------------------------------------------------------------
127 Int_t AliITSneuralTracker::ReadFile(const Text_t *fname, Int_t evnum)
128 {
129         // File reader
130         // -----------
131         // Reads a ROOT file in order to retrieve recpoints.
132         // the 'evnum' argument can have two meanings:
133         // if it is >= 0 it represents the event number to retrieve the correct gAlice
134         // if it is < 0, instead, this means that the open file contains 
135         // a 'TreeP' tree, with all points which remained unused 
136         // after the Kalman tracking procedure (produced via the ITSafterKalman.C macro)
137         // the return value is the number of collected points 
138         // (then, if errors occur, the method returns 0)
139         
140         TFile *file = new TFile(fname, "read");
141         TTree *tree = 0;
142         Int_t i, nentries, total, sector, mesh;
143         
144         if (evnum < 0) {
145                 tree = (TTree*)file->Get("TreeP0");
146                 if (!tree) {
147                         Error("", "Specifying evnum < 0, I expect to find a 'TreeP' tree in the file, but there isn't any");
148                         return 0;
149                 }
150                 AliITSglobalRecPoint *p = 0;
151                 tree->SetBranchAddress("Points", &p);
152                 nentries = (Int_t)tree->GetEntries();
153                 total = 0; 
154                 for (i = 0; i < nentries; i++) {
155                         tree->GetEntry(i);
156                         AliITSglobalRecPoint *q = new AliITSglobalRecPoint(*p);
157                         sector = (Int_t)(q->fPhi / fSectorWidth);
158                         mesh = (Int_t)(q->fTheta / fThetaCut);
159                         q->fUsed = 0;
160                         fPoints[q->fLayer][sector][mesh]->AddLast(q);
161                         total++;
162                 }
163                 return total;
164         }
165         
166         if (gAlice) { delete gAlice; gAlice = 0; }
167         gAlice = (AliRun*)file->Get("gAlice");
168         if (!gAlice) {
169                 Error("", "gAlice is NULL, maybe wrong filename...");
170                 return 0;
171         }
172
173         Int_t nparticles = (Int_t)gAlice->GetEvent(evnum);
174         if (!nparticles) {
175                 Error("", "No particles found");
176                 return 0;
177         }
178
179         AliITS *its = (AliITS*)gAlice->GetModule("ITS");
180         if (!its) {
181                 Error("", "No ITS found");
182                 return 0;
183         }
184
185         AliITSgeom *geom = (AliITSgeom*)its->GetITSgeom();
186         if (!geom) {
187                 Error("", "AliITSgeom is NULL");
188                 return 0;
189         }
190
191         tree = gAlice->TreeR();
192         if (!tree) {
193                 Error("", "TreeR() returned NULL");
194                 return 0;
195         }
196
197         nentries = (Int_t)tree->GetEntries();
198         if (!nentries) {
199                 Error("", "TreeR is empty");
200                 return 0;
201         }
202
203         // Reading objects initialization
204         Int_t k, entry, lay, lad, det;
205         Double_t l[3], g[3], ll[3][3], gg[3][3];
206         
207         TObjArray *recsArray = its->RecPoints();
208         AliITSRecPoint *recp;
209         AliITSglobalRecPoint *pnt = 0;
210         
211         total = 0;
212         for(entry = 0; entry < nentries; entry++) {
213                 if (entry > geom->GetLastSSD()) continue;
214                 its->ResetRecPoints();
215                 tree->GetEvent(entry);
216                 TObjArrayIter recs(recsArray);
217                 AliITSgeomMatrix *gm = geom->GetGeomMatrix(entry);
218                 geom->GetModuleId(entry, lay, lad, det);
219                 lay--;
220                 while ((recp = (AliITSRecPoint*)recs.Next())) {
221                         // conversion to global coordinate system
222                         for (i = 0; i < 3; i++) {
223                                 l[i] = g[i] = 0.0;
224                                 for (k = 0; k < 3; k++) {
225                                         ll[i][k] = 0.0;
226                                         gg[i][k] = 0.0;
227                                 }
228                         }
229                         // local to global conversions of coords
230                         l[0] = recp->fX;
231                         l[1] = 0.0;
232                         l[2] = recp->fZ;
233                         gm->LtoGPosition(l, g);
234                         // local to global conversions of sigmas
235                         ll[0][0] = recp->fSigmaX2;
236                         ll[2][2] = recp->fSigmaZ2;
237                         gm->LtoGPositionError(ll, gg);
238                         // instantiation of a new global rec-point
239                         pnt = new AliITSglobalRecPoint(g[0], g[1], g[2], gg[0][0], gg[1][1], gg[2][2], lay);
240                         // copy of owner track labels
241                         for (i = 0; i < 3; i++) pnt->fLabel[i] = recp->fTracks[i];
242                         sector = (Int_t)(pnt->fPhi / fSectorWidth);
243                         mesh = (Int_t)(pnt->fTheta / fThetaCut);
244                         pnt->fUsed = 0;
245                         fPoints[lay][sector][mesh]->AddLast(pnt);
246                         total++;
247                 }
248         }
249         
250         return total;
251 }
252 //---------------------------------------------------------------------------------------------------------
253 void AliITSneuralTracker::Go(const char* filesave, Bool_t flag)
254 {
255         // Global working method
256         // ---------------------
257         // It's the method which calls all other ones,
258         // in order to analyze the whole ITS points ensemble
259         // It's the only method to use, besides the parameter setters,
260         // so, all other methods are defined as private
261         // (see header file)
262         // the flag meaning is explained in the 'Initialize' method
263         
264         Int_t i; 
265         for (i = 0; i < fSectorNum; i++) DoSector(i, flag);
266         cout << endl << "Saving results in file " << filesave << "..." << flush;
267         TFile *f = new TFile(filesave, "recreate");
268         fTracks->Write();
269         cout << "done" << endl;
270         f->Close();
271 }
272
273 //=========================================================================================================
274 // * * *   P R I V A T E   S E C T I O N   * * *
275 //=========================================================================================================
276
277 Int_t AliITSneuralTracker::Initialize(Int_t secnum, Int_t curvnum, Bool_t flag)
278 {
279         // Network initializer
280         // -------------------
281         // Fills the neuron arrays, following the cut criteria for the selected step
282         // (secnum = sector to analyze, curvnum = curvature cut step to use)
283         // It also sets the initial random activation.
284         // In order to avoid a large number of 'new' operations, all existing neurons
285         // are reset and initialized with the new values, and are created new unit only if
286         // it turns out to be necessary
287         // the 'flag' argument is used to decide if the lower edge in the intevral
288         // of the accepted curvatures is given by zero (kFALSE) or by the preceding used cut (kTRUE)
289         // (if this is the first step, anyway, the minimum is always zero)
290         
291         Int_t l, m, k, neurons = 0;
292         
293         for (l = 0; l < 5; l++) fNeurons[l]->Delete();
294         
295         AliITSneuron *unit = 0;
296         Double_t abscurv, max = fCurvCut[curvnum], min = 0.0;
297         if (flag && curvnum > 0) min = fCurvCut[curvnum - 1];
298         AliITSglobalRecPoint *inner = 0, *outer = 0;
299         for (l = 0; l < 5; l++) {
300                 for (m = 0; m < fThetaNum; m++) {
301                         TObjArrayIter inners(fPoints[l][secnum][m]);
302                         while ((inner = (AliITSglobalRecPoint*)inners.Next())) {
303                                 if (inner->fUsed > 0) continue; // points can't be recycled
304                                 for (k = m-1; k <= m+1; k++) {
305                                         if (k < 0 || k >= fThetaNum) continue; // to avoid seg faults 
306                                         TObjArrayIter outers(fPoints[l+1][secnum][k]);
307                                         while ((outer = (AliITSglobalRecPoint*)outers.Next())) {
308                                                 if (outer->fUsed > 0) continue;  // points can't be recycled
309                                                 if (inner->DTheta(outer) > fThetaCut) continue;
310                                                 unit = new AliITSneuron;
311                                                 unit->fInner = inner;
312                                                 unit->fOuter = outer;
313                                                 CalcParams(unit);
314                                                 abscurv = TMath::Abs(unit->fCurv);
315                                                 if (unit->fDiff > fDiff || abscurv < min || abscurv > max) {
316                                                         delete unit;
317                                                         continue;
318                                                 }
319                                                 unit->fActivation = gRandom->Rndm() * (fMax - fMin) + fMin;
320                                                 fNeurons[l]->AddLast(unit);
321                                                 neurons++;
322                                         } // end loop on candidates for outer point for neurons
323                                 }
324                         } // end loop on candidates inner points for neuron
325                 }
326         } // end loop on layers
327         return neurons;
328 }
329 //---------------------------------------------------------------------------------------------------------
330 Bool_t AliITSneuralTracker::Update()
331 {
332         // Updating cycle
333         // --------------
334         // Performs an updating cycle, by summing all the
335         // gain and cost contribution for each neuron and
336         // updating its activation.
337         // An asynchronous method is followed, with the order
338         // according that the neurons have been created
339         // Time wasting is avoided by dividing the neurons 
340         // into different arrays, depending on the layer of their inner point.
341         // The return value answers the question: "The network has stabilized?"
342         
343         Int_t j, l;
344         Double_t actOld, actNew, argNew, totDiff = 0.0;
345         Double_t sumInh, sumExc, units = 0.0;
346         AliITSneuron *me = 0, *it = 0;
347         for (l = 0; l < 5; l++) {
348                 TObjArrayIter meIter(fNeurons[l]);
349                 while ((me = (AliITSneuron*)meIter.Next())) {
350                         units++;
351                         TObjArrayIter inhIter(fNeurons[l]);
352                         
353                         // inhibition (search only for neurons starting in the same layer)
354                         sumInh = 0.0;
355                         while((it = (AliITSneuron*)inhIter.Next())) {
356                                 if (it->fOuter == me->fOuter || it->fInner == me->fInner)       
357                                         sumInh += it->fActivation;
358                         }
359                         
360                         // excitation (search only for neurons 
361                         // which start in the previous or next layers)
362                         sumExc = 0.0;
363                         for (j = l - 1; j <= l + 1; j += 2) {
364                                 if (j < 0 || j >= 5) continue;
365                                 TObjArrayIter itIter(fNeurons[j]);
366                                 while ((it = (AliITSneuron*)itIter.Next())) {
367                                         if (it->fInner == me->fOuter || it->fOuter == me->fInner && it->fCurv * me->fCurv > 0.0) {
368                                                 sumExc += Weight(me, it) * it->fActivation;
369                                         }
370                                 }
371                         } // end search for excitories
372                         
373                         actOld = me->fActivation;
374                         argNew = fRatio * sumExc - sumInh;
375                         actNew = 1.0 / (1.0 + TMath::Exp(-argNew / fTemp));
376                         me->fActivation = actNew;
377                         // calculation the relative activation variation 
378                         // (for stabilization check)
379                         totDiff += TMath::Abs((actNew - actOld) / actOld);
380                 } // end loop over 'me' (updated neuron)
381         } // end loop over layers
382         totDiff /= units;
383         
384         return (totDiff < fVar);
385 }
386 //---------------------------------------------------------------------------------------------------------
387 Int_t AliITSneuralTracker::Save()
388 {
389         // Tracki saving method
390         // --------------------
391         // Stores the recognized tracks into the TObjArray 'fTracks'
392         // but doesn't save it into a file until the whole work has ended
393         
394         Int_t l;
395         Double_t test = 0.5;
396         
397         // Before all, for each neuron is found the best active sequences
398         // among the incoming and outgoing other units
399         AliITSneuron *main, *fwd, *back, *start;
400         for (l = 0; l < 5; l++) {
401                 TObjArrayIter mainIter(fNeurons[l]);
402                 while((main = (AliITSneuron*)mainIter.Next())) {
403                         if (l < 4) {
404                                 TObjArrayIter fwdIter(fNeurons[l+1]);
405                                 test = 0.5;
406                                 while ((fwd = (AliITSneuron*)fwdIter.Next())) {
407                                         if (main->fOuter == fwd->fInner && fwd->fActivation > test) {
408                                                 test = fwd->fActivation;
409                                                 main->fNext = fwd;
410                                         }
411                                 };
412                         }
413                         if (l > 0) {
414                                 TObjArrayIter backIter(fNeurons[l-1]);
415                                 test = 0.5;
416                                 while ((back = (AliITSneuron*)backIter.Next())) {
417                                         if (main->fInner == back->fOuter && back->fActivation > test) {
418                                                 test = back->fActivation;
419                                                 main->fPrev = back;
420                                         }
421                                 };
422                         }
423                 }
424         }
425         
426         // Then, are resolved the mismatches in the chains found above 
427         // (when unit->next->prev != unit or unit->prev->next != unit)
428         for (l = 0; l < 5; l++) {
429                 TObjArrayIter mainIter(fNeurons[l]);
430                 while((main = (AliITSneuron*)mainIter.Next())) {
431                         fwd  = main->fNext;
432                         back = main->fPrev;
433                         if (fwd  != 0 && fwd->fPrev  != main) main->fNext = 0; 
434                         if (back != 0 && back->fNext != main) main->fPrev = 0; 
435                 }
436         }
437         
438         Int_t pos, neg, incr, decr, ntracks;
439         AliITSneuralTrack *trk = 0;
440                 
441         TObjArrayIter startIter(fNeurons[0]);
442         for (;;) {
443                 start = (AliITSneuron*)startIter.Next();
444                 if (!start) break;
445                 Int_t pts = 0;
446                 // with the chain above defined, a track can be followed
447                 // like a linked list structure
448                 for (main = start; main; main = main->fNext) pts++;
449                 // the count will be  > 5 only if the track contains 6 points
450                 // (what we want)
451                 if (pts < 5) continue;
452                 // track storage
453                 trk = new AliITSneuralTrack;
454                 trk->CopyPoint(start->fInner);
455                 for (main = start; main; main = main->fNext) {
456                         //main->fInner->fUsed = kTRUE;
457                         //main->fOuter->fUsed = kTRUE;
458                         trk->CopyPoint(main->fOuter);
459                 }
460                 trk->Kinks(pos, neg, incr, decr);
461                 if (pos != 0 && neg != 0) {
462                         cout << "pos, neg kinks = " << pos << " " << neg << endl;
463                         continue;
464                 }
465                 if (incr != 0 && decr != 0) {
466                         cout << "increment, decrements in delta phi = " << incr << " " << decr << endl;
467                         continue;
468                 }
469                 for (main = start; main; main = main->fNext) {
470                         main->fInner->fUsed++;
471                         main->fOuter->fUsed++;
472                 }
473                 fTracks->AddLast(trk);
474         }
475         ntracks = (Int_t)fTracks->GetEntriesFast();
476         return ntracks;
477 }
478 //---------------------------------------------------------------------------------------------------------
479 void AliITSneuralTracker::DoSector(Int_t sect, Bool_t flag)
480 {
481         // Sector recognition
482         // ------------------
483         // This is a private method which works on a single sector
484         // just for an easier readability of the code.
485         // The 'flag' argument is needed to be given to the 'Initialize' method (see it)
486         
487         Int_t curvnum, nTracks = 0, nUnits;
488         Bool_t isStable = kFALSE;
489         for (curvnum = 0; curvnum < fCurvNum; curvnum++) {
490                 cout << "\rSector " << sect << ":" << ((curvnum+1)*100)/fCurvNum << "%" << flush;
491                 nUnits = Initialize(sect, curvnum, flag);
492                 if (!nUnits) continue;
493                 do {
494                         isStable = Update();
495                 } while (!isStable);
496                 nTracks = Save();
497         }
498         cout << "\rTracks stored after sector " << sect << ": " << nTracks << endl;
499 }
500 //---------------------------------------------------------------------------------------------------------
501 Int_t AliITSneuralTracker::CountExits(AliITSneuron *n)
502 {
503         // Counter for neurons which enter in 'n' tail
504         
505         Int_t count = 0, l = n->fOuter->fLayer;
506         if (l == 5) return 0;
507         
508         AliITSneuron *test = 0;
509         TObjArrayIter iter(fNeurons[l]);
510         while ((test = (AliITSneuron*)iter.Next())) {
511                 if (test->fInner == n->fOuter) count++;
512         }
513         return count;
514 }
515 //---------------------------------------------------------------------------------------------------------
516 Int_t AliITSneuralTracker::CountEnters(AliITSneuron *n)
517 {
518         // Counter for neurons which exit from 'n' head
519         
520         Int_t count = 0, l = n->fInner->fLayer;
521         if (l == 0) return 0;
522         
523         AliITSneuron *test = 0;
524         TObjArrayIter iter(fNeurons[l-1]);
525         while ((test = (AliITSneuron*)iter.Next())) {
526                 if (test->fOuter == n->fInner) count++;
527         }
528         return count;
529 }
530 //---------------------------------------------------------------------------------------------------------
531 Bool_t AliITSneuralTracker::ResetNeuron(AliITSneuron *n)
532 {
533         // A method which resets the neuron's pointers to zero
534         
535         if (!n) return kFALSE;
536         n->fNext = 0;
537         n->fPrev = 0;
538         n->fInner = 0;
539         n->fOuter = 0;
540         // the other datamembers don't need to be reinitialized
541         return kTRUE;
542 }
543 //---------------------------------------------------------------------------------------------------------
544 Bool_t AliITSneuralTracker::SetEdge(AliITSneuron *n, AliITSglobalRecPoint *p, const char what)
545 {
546         // Sets the indicated edge point to the neuron
547         // (if the arguments are suitable ones...)
548         
549         if (!n || !p || (what != 'i' && what != 'I' && what != 'o' && what != 'O')) return kFALSE;
550         if (what == 'i' || what == 'I') 
551                 n->fInner = p;
552         else
553                 n->fOuter = p;
554         return kTRUE;
555 }
556 //---------------------------------------------------------------------------------------------------------
557 Bool_t AliITSneuralTracker::CalcParams(AliITSneuron *n)
558 {
559         // This method evaluates the helix parameters from the edged of a neuron
560         // (curvature, phase parameter, tangent of lambda angle)
561         // if the p[arameters assume too strange (and unphysical) values, the 
562         // method return kFALSE, else it return kTRUE;
563         
564         if (!n) return kFALSE;
565         
566         Double_t sign(1.0), det, rc = 0.0;
567         Double_t tanl1, tanl2, l1, l2;
568         // changing the reference frame into the one centered in the vertex
569         Double_t x1 = n->fInner->fGX - fVPos[0];
570         Double_t y1 = n->fInner->fGY - fVPos[1];
571         Double_t z1 = n->fInner->fGZ - fVPos[2];
572         Double_t r1 = x1 * x1 + y1 * y1;
573         Double_t x2 = n->fOuter->fGX - fVPos[0];
574         Double_t y2 = n->fOuter->fGY - fVPos[1];
575         Double_t z2 = n->fOuter->fGZ - fVPos[2];
576         Double_t r2 = x2 * x2 + y2 * y2;
577         // initialization of the XY-plane data (curvature and centre)
578         // (will remain these values if we encounter errors)
579         n->fCurv = n->fCX = n->fCY = n->fTanL = n->fPhase = 0.0;
580         n->fLength = (n->fOuter->fGX - n->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX);
581         n->fLength += (n->fOuter->fGY - n->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY);
582         n->fLength += (n->fOuter->fGZ - n->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ);
583         n->fLength = TMath::Sqrt(n->fLength);
584         // this determinant which is zero when the points are in a line
585         det = 2.0 * (x1 * y2 - x2 * y1);
586         if (det == 0.0)
587                 return kFALSE; // it is difficult having perfectly aligned points --- it's more likely to be an error
588         else {
589                 n->fCX = (r1*y2 - r2*y1) / det;
590                 n->fCY = (r2*x1 - r1*x2) / det;
591                 rc = TMath::Sqrt(n->fCX * n->fCX + n->fCY * n->fCY);
592                 r1 = TMath::Sqrt(r1);
593                 r2 = TMath::Sqrt(r2);
594                 sign = (n->fOuter->fPhi >= n->fInner->fPhi) ? 1.0 : -1.0;
595                 if (rc > 0.0) {
596                         n->fCurv = sign / rc;
597                         n->fPhase = TMath::ATan2(-n->fCX, -n->fCY);
598                         if (n->fPhase < 0.0) n->fPhase += 2.0 * TMath::Pi(); // angles in 0 --> 2pi
599                         if (r1 <= 2.0 * rc && r2 <= 2.0 * rc) {
600                                 l1 = 2.0 * rc * TMath::ASin(r1 / (2.0 * rc));
601                                 l2 = 2.0 * rc * TMath::ASin(r2 / (2.0 * rc));
602                                 tanl1 = z1 / l1;
603                                 tanl2 = z2 / l2;
604                                 n->fDiff = TMath::Abs(tanl1 - tanl2);
605                                 n->fTanL = 0.5 * (tanl1 + tanl2);
606                         }
607                         else {
608                                 n->fTanL = 0.0;
609                                 n->fDiff = 100.0;
610                                 return kFALSE; // it' even more difficult that the arguments above aren't acceptable for arcsine
611                         }
612                 }
613                 else
614                         return kFALSE;
615         }
616         return kTRUE;
617 }
618 //---------------------------------------------------------------------------------------------------------
619 Double_t AliITSneuralTracker::Angle(AliITSneuron *n, AliITSneuron *m)
620 {
621         // calculates the angle between two segments
622         // but doesn't check if they have a common point.
623         // The calculation is made by the ratio of their scalar product
624         // to the product of their magnitudes
625         
626         Double_t x = (m->fOuter->fGX - m->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX);
627         Double_t y = (m->fOuter->fGY - m->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY);
628         Double_t z = (m->fOuter->fGZ - m->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ);
629                 
630         Double_t cosine = (x + y + z) / (n->fLength * m->fLength);
631         
632         if (cosine > 1.0 || cosine < -1.0) {
633                 Warning("AliITSneuronV1::Angle", "Strange value of cosine");
634                 return 0.0;
635         }
636         
637         return TMath::ACos(cosine);
638 }
639 //---------------------------------------------------------------------------------------------------------
640 Double_t AliITSneuralTracker::Weight(AliITSneuron *n, AliITSneuron *m)
641 {
642         // calculation of the excitoy weight
643         
644 //      Double_t v1 = TMath::Abs(fTanL - n->fTanL);
645 //      Double_t v2 = TMath::Abs(TMath::Abs(fCurv) - TMath::Abs(n->fCurv));
646 //      Double_t v3 = TMath::Abs(fPhase - n->fPhase);
647 //      Double_t q = (v1 + v2 + v3) / 3.0;
648 //      Double_t a = 1.0; // 0.02 / (fDiff + n->fDiff);
649         Double_t b = 1.0 - TMath::Sin(Angle(n, m));
650         return TMath::Power(b, fExp);
651 //      return TMath::Exp(q / par1);
652
653 }
654 //---------------------------------------------------------------------------------------------------------
655 // ====================================================================================================
656 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuron <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
657 // ====================================================================================================
658
659 AliITSneuron::AliITSneuron()
660 {
661         // Default constructor
662         // -------------------
663         // sets all parameters to default (but unuseful) 
664         // (in particular all pointers are set to zero)
665         // There is no alternative constructor, because 
666         // this class should not be used by itself, but 
667         // it is only necessary to allow the neural tracker to work
668         // In fact, all methods which operate on neuron's datamembers
669         // are stored in the AliITSneuralTracker class anyway
670         
671         fActivation = 0.0;
672         
673         fCurv = 0.0;
674         fCX = 0.0;
675         fCY = 0.0;
676         fTanL = 0.0;
677         fPhase = 0.0;
678         fDiff = 0.0;
679         fLength = 0.0;
680         
681         fNext = fPrev = 0;
682         fInner = fOuter = 0;
683 }
684 //---------------------------------------------------------------------------------------------------------
685 AliITSneuron::~AliITSneuron() 
686 {
687         // Destructor
688         // ----------
689         // does nothing, because there is no need to affect parameters,
690         // while the AliITSglobalRecPoint pointers can't be deleted,
691         // and the AliITSneuron pointers belong to a TObjArray which will
692         // be cleared when necessary...
693 }
694