correcting compilation error in debug mode
[u/mrichter/AliRoot.git] / HLT / MUON / OnlineAnalysis / AliHLTMUONFullTracker.cxx
1 /**************************************************************************
2  * This file is property of and copyright by the ALICE HLT Project        * 
3  * All rights reserved.                                                   *
4  *                                                                        *
5  * Primary Authors:                                                       *
6  *   Indranil Das <indra.das@saha.ac.in>                                  *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          * 
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 //////////////////////////////////////////////////////////////////////
18 ///Author : Indranil Das, SINP, INDIA
19 ///         
20 ///Email :  indra.das@saha.ac.in || indra.ehep@gmail.com
21 ///
22 /// This class is created to peform the a full track reconstroction
23 /// for online HLT for Dimuon Detector. It is based on the method of 
24 /// straight line tracking for slat detectors and cellular automation
25 /// mehtods for finding the tracks in the quadrant detectos. The track
26 /// segments of the front and back side of the spectromrter is either
27 /// compared using KalmanChi2 method (a slightly modified version of 
28 /// AliMUONTrackReconstructorK) or alternated method of simple 
29 /// extrapolation of tracks through dipole magnet. Finally the Track is
30 /// passed through Absorber to take the MCS effect and Branson Effect 
31 /// into account for correction of track parameters.
32 /////////////////////////////////////////////////////////////////////////
33
34
35 #include "AliHLTMUONFullTracker.h"
36
37 #ifdef PRINT_FULL
38 #define PRINT_POINTS 1
39 #define PRINT_BACK 1
40 #define PRINT_FRONT 1
41 #define PRINT_KALMAN 1
42 #define PRINT_SELECT 1
43 ///#define PRINT_DETAIL_KALMAN 1
44 #endif
45
46 class AliHLTMUONFullTracker;
47
48 const Float_t AliHLTMUONFullTracker::TrackDetCoordinate[3] = {
49   155.179+20.0,  166.234+20.0, 
50   (AliMUONConstants::DefaultChamberZ(4)+ AliMUONConstants::DefaultChamberZ(5))/2.0
51 }; 
52
53 const Double_t AliHLTMUONFullTracker::fgkAbsoedge[4] = {92.0, 317.0,443.0,499.0};
54 const Double_t AliHLTMUONFullTracker::fgkRadLen[3] = {13.875635, 11.273801, 1.765947};
55 const Double_t AliHLTMUONFullTracker::fgkRho[3] = {1.750000, 2.350000, 7.880000};
56 const Double_t AliHLTMUONFullTracker::fgkAtomicZ[3] = {6.000000, 11.098486,25.780000};
57 const Double_t AliHLTMUONFullTracker::fgkAtomicA[3] = {12.010000, 22.334156,55.299670 };
58
59
60 const Int_t AliHLTMUONFullTracker::fgkMaxNofCellsPerCh = 400;
61 const Int_t AliHLTMUONFullTracker::fgkMaxNofPointsPerCh = 300; 
62 const Int_t AliHLTMUONFullTracker::fgkMaxNofCh = 11 ; /// 10tracking + 1 trigger
63 const Int_t AliHLTMUONFullTracker::fgkMaxNofTracks = 200;
64 const Int_t AliHLTMUONFullTracker::fgkMaxNofConnectedTracks = 20;
65
66
67
68 AliHLTMUONFullTracker::AliHLTMUONFullTracker() : 
69   AliHLTLogging(),
70   fChamberGeometryTransformer(0x0),
71   fChPoint(0x0),
72   fChPoint11(0x0),
73   fBackTrackSeg(0x0),
74   fFrontTrackSeg(0x0),
75   fExtrapSt3X(0x0),
76   fExtrapSt3Y(0x0),
77   fInclinationBack(0x0),
78   fNofConnectedfrontTrackSeg(0x0),
79   fBackToFront(0x0),
80   fCharge(0x0),
81   fNofPoints(0x0),
82   fTrackParam(0x0),
83   fTotNofPoints(0),
84   fTotTrackSeg(0),
85   fOverflowed(false),
86   fNofbackTrackSeg(0),
87   fNoffrontTrackSeg(0),
88   fNofConnected(0)
89 {
90
91   /// Constructor of the class
92
93   /// fChPoint = (AliHLTMUONRecHitStruct***)(malloc(10 * sizeof(AliHLTMUONRecHitStruct)));
94   /// for( Int_t ich=0;ich<10;ich++)
95   ///   fChPoint[ich] = (AliHLTMUONRecHitStruct**)(malloc(300 * sizeof(AliHLTMUONRecHitStruct)));
96
97   try{
98     fChPoint = new AliHLTMUONRecHitStruct**[fgkMaxNofCh-1];
99   }catch (const std::bad_alloc&){
100     HLTError("Dynamic memory allocation failed in constructor : fChPoint");
101     throw;
102   }
103
104   for( Int_t ich=0;ich<(fgkMaxNofCh-1);ich++)
105     try{
106       fChPoint[ich] = new AliHLTMUONRecHitStruct*[fgkMaxNofPointsPerCh];
107     }catch (const std::bad_alloc&){
108       HLTError("Dynamic memory allocation failed in constructor : fChPoint[%d]",ich);
109       throw;
110     }
111   
112   try{
113     fChPoint11 = new AliHLTMUONTriggerRecordStruct*[fgkMaxNofPointsPerCh];
114   }catch (const std::bad_alloc&){
115     HLTError("Dynamic memory allocation failed in constructor : fChPoint11");
116     throw;
117   }
118
119   try{
120     fBackTrackSeg = new TrackSeg[fgkMaxNofTracks];
121   }catch (const std::bad_alloc&){
122     HLTError("Dynamic memory allocation failed in constructor : fBackTrackSeg");
123     throw;
124   }
125
126   try{
127     fFrontTrackSeg = new TrackSeg[fgkMaxNofTracks];
128   }catch (const std::bad_alloc&){
129     HLTError("Dynamic memory allocation failed in constructor : fFrontTrackSeg");
130     throw;
131   }
132
133   try{
134     fExtrapSt3X = new float[fgkMaxNofTracks];
135   }catch (const std::bad_alloc&){
136     HLTError("Dynamic memory allocation failed in constructor : fExtrapSt3X");
137     throw;
138   }
139
140   try{
141     fExtrapSt3Y = new float[fgkMaxNofTracks];
142   }catch (const std::bad_alloc&){
143     HLTError("Dynamic memory allocation failed in constructor : fExtrapSt3Y");
144     throw;
145   }
146
147   try{
148     fInclinationBack = new float[fgkMaxNofTracks];
149   }catch (const std::bad_alloc&){
150     HLTError("Dynamic memory allocation failed in constructor : fInclinationBack");
151     throw;
152   }
153   
154   try{
155     fNofConnectedfrontTrackSeg = new int[fgkMaxNofTracks];
156   }catch (const std::bad_alloc&){
157     HLTError("Dynamic memory allocation failed in constructor : fNofConnectedfrontTrackSeg");
158     throw;
159   }
160
161   try{
162     fBackToFront = new int*[fgkMaxNofTracks];
163   }catch (const std::bad_alloc&){
164     HLTError("Dynamic memory allocation failed in constructor : fBackToFront");
165     throw;
166   }
167
168   for( Int_t itr=0;itr<fgkMaxNofTracks;itr++)
169     try{
170       fBackToFront[itr] = new int[fgkMaxNofConnectedTracks];
171     }catch (const std::bad_alloc&){
172       HLTError("Dynamic memory allocation failed in constructor : fBackToFront[%d]",itr);
173       throw;
174     }
175
176   
177   try{
178     fCharge = new float[fgkMaxNofTracks];
179   }catch (const std::bad_alloc&){
180     HLTError("Dynamic memory allocation failed in constructor : fCharge");
181     throw;
182   }
183
184   try{
185     fNofPoints = new int[fgkMaxNofCh];
186   }catch (const std::bad_alloc&){
187     HLTError("Dynamic memory allocation failed in constructor : fNofPoints");
188     throw;
189   }
190
191   try{
192     fTrackParam = new AliMUONTrackParam[fgkMaxNofTracks];
193   }catch (const std::bad_alloc&){
194     HLTError("Dynamic memory allocation failed in constructor : fTrackParam");
195     throw;
196   }
197
198   Clear();
199
200 }
201
202 ///__________________________________________________________________________
203
204 AliHLTMUONFullTracker::~AliHLTMUONFullTracker()
205
206   /// Destructor of the class
207
208   fChamberGeometryTransformer->Delete();
209   ///free(fChPoint);
210   delete []fChPoint;
211   delete []fChPoint11;
212   delete []fBackTrackSeg;
213   delete []fFrontTrackSeg;
214   delete []fExtrapSt3X;
215   delete []fExtrapSt3Y;
216   delete []fInclinationBack;
217   delete []fNofConnectedfrontTrackSeg;
218   delete []fBackToFront;
219   delete []fCharge;
220   delete []fNofPoints;
221   delete []fTrackParam;
222
223   ///delete fChamberGeometryTransformer;
224
225 }
226
227 ///__________________________________________________________________________
228
229 Bool_t AliHLTMUONFullTracker::Clear(){
230
231   /// Clear method to be called after each event
232
233   for( Int_t ich=0;ich<fgkMaxNofCh;ich++)
234     fNofPoints[ich] = 0;
235   
236   fNofbackTrackSeg = 0;
237   fNoffrontTrackSeg = 0;
238   fNofConnected = 0;
239
240   return true;
241 }
242
243 ///__________________________________________________________________________
244
245 void AliHLTMUONFullTracker::Print()
246 {
247   /// Print anything mainly for debugging the codes, will be removed later
248   HLTInfo("\nPrint This--->\n");
249 }
250
251 ///__________________________________________________________________________
252
253 Bool_t AliHLTMUONFullTracker::Init()
254 {
255   /// Initilation to be called once, later can be used to set/load the CDB path/entries
256
257   if (AliGeomManager::GetGeometry() == NULL){
258     AliGeomManager::LoadGeometry();
259     AliGeomManager::ApplyAlignObjsFromCDB("GRP MUON");
260     InitGRP();
261   }
262
263   AliMUONTrackExtrap::SetField();
264   /// see header file for class documentation
265   fChamberGeometryTransformer = new AliMUONGeometryTransformer();
266   if(! fChamberGeometryTransformer->LoadGeometryData()){
267     cerr<<__FILE__<<": Failed to Load Geomerty Data "<<endl;
268   }
269   
270
271   return true;
272 }
273
274 ///__________________________________________________________________________
275
276 Bool_t AliHLTMUONFullTracker::SetInput(AliHLTInt32_t /*ddl*/,const AliHLTMUONRecHitStruct  *data, AliHLTInt32_t size)
277 {
278   /// Set the input for rechit points
279
280   AliHLTUInt16_t detElemID;
281   AliHLTUInt8_t chamber;
282   
283   for( Int_t ipoint=0;ipoint<int(size);ipoint++){
284     AliHLTMUONUtils::UnpackRecHitFlags(data->fFlags,chamber,detElemID);
285     fChPoint[detElemID/100-1][fNofPoints[detElemID/100-1]++]  = (AliHLTMUONRecHitStruct  *)data;
286 #ifdef PRINT_POINTS  
287     HLTInfo("\nch : %02d, detelem : %04d, (X,Y,Z) : (%8.3f,%8.3f,%8.3f)",chamber,detElemID,data->fX,data->fY,data->fZ);
288 #endif
289     data++;
290   }
291
292   return true;
293 }
294
295 ///__________________________________________________________________________
296
297 Bool_t AliHLTMUONFullTracker::SetInput(AliHLTInt32_t /*ddl*/, const AliHLTMUONTriggerRecordStruct  *data, AliHLTInt32_t size)
298 {
299   /// Set the input for trigrecs
300
301   AliHLTUInt16_t detElemID;
302   AliHLTUInt8_t chamber;
303
304   for( Int_t ipoint=0;ipoint<int(size);ipoint++){
305     fChPoint11[fNofPoints[10]++] = (AliHLTMUONTriggerRecordStruct *)data;
306     for( Int_t ich=0;ich<4;ich++){
307       AliHLTMUONUtils::UnpackRecHitFlags((data->fHit[ich]).fFlags,chamber,detElemID);
308 #ifdef PRINT_POINTS  
309       HLTInfo("\nsize : %d, itrig  : %04d, ch : %02d, detelem : %04d, (X,Y,Z) : (%8.3f,%8.3f,%8.3f)",
310              size,ipoint,chamber,detElemID,(data->fHit[ich]).fX,(data->fHit[ich]).fY,(data->fHit[ich]).fZ);
311 #endif
312     }///ich loop
313     data++;
314   }///ipoint
315   return true;
316 }
317
318 ///__________________________________________________________________________
319   
320 Bool_t AliHLTMUONFullTracker::Run( Int_t /*iEvent*/,AliHLTMUONMansoTrackStruct *data, AliHLTUInt32_t& size)
321 {
322   /// Main Run call of the class
323
324   SlatTrackSeg();
325   QuadTrackSeg();
326   KalmanChi2Test();
327   ExtrapolateToOrigin(true);
328   FillOutData(data,size);
329 //   HLTDebug("iEvent: %d, has : %d tracks, triggers : %d, nof slat tracks : %d, quad tracks : %d, connected : %d\n",
330 //       iEvent,size,fNofPoints[10],fNofbackTrackSeg,fNoffrontTrackSeg,fNofConnected);
331   Clear();
332   return true;
333 }
334  
335 ///__________________________________________________________________________
336
337 void AliHLTMUONFullTracker::Sub(AliHLTMUONRecHitStruct *v1, AliHLTMUONRecHitStruct *v2, AliHLTMUONRecHitStruct *v3)
338 {
339   /// Subtraction of position co-odinate of two space points
340
341   v3->fX = v1->fX - v2->fX; 
342   v3->fY = v1->fY - v2->fY; 
343   v3->fZ = v1->fZ - v2->fZ;
344 };
345
346 ///__________________________________________________________________________
347
348 Double_t AliHLTMUONFullTracker::Angle(AliHLTMUONRecHitStruct *v1, AliHLTMUONRecHitStruct *v2)
349 {
350   ///Angle of a straight line formed using v1 and v2
351
352   Float_t ptot2 = ((v1->fX * v1->fX) + (v1->fY * v1->fY) + (v1->fZ * v1->fZ))*
353     ((v2->fX * v2->fX) + (v2->fY * v2->fY) + (v2->fZ * v2->fZ));
354   if(ptot2 <= 0) {
355     return 0.0;
356   } else {
357     Float_t arg = ((v1->fX * v2->fX) + (v1->fY * v2->fY) + (v1->fZ * v2->fZ))/sqrt(ptot2);
358     if(arg >  1.0) arg =  1.0;
359     if(arg < -1.0) arg = -1.0;
360     return TMath::ACos(arg);
361     ///return acos(arg);
362   }
363   
364 }
365
366 ///__________________________________________________________________________
367
368 Bool_t AliHLTMUONFullTracker::FillOutData(AliHLTMUONMansoTrackStruct *track, AliHLTUInt32_t& size)
369 {
370   ///Fill the output data pointers
371
372   size = (AliHLTUInt32_t(fNofbackTrackSeg)<size) ? AliHLTUInt32_t(fNofbackTrackSeg) : size;
373
374   for( Int_t ibackTrackSeg=0;ibackTrackSeg<int(size);ibackTrackSeg++){
375     
376     track->fId = (ibackTrackSeg << 8) | (fBackTrackSeg[ibackTrackSeg].fTrigRec & 0xFF);
377     track->fTrigRec = fBackTrackSeg[ibackTrackSeg].fTrigRec;
378     track->fPx = fTrackParam[ibackTrackSeg].Px();
379     track->fPy = fTrackParam[ibackTrackSeg].Py();
380     track->fPz = fTrackParam[ibackTrackSeg].Pz();  
381     
382     track->fChi2 = 0;
383     Bool_t hitset[4];
384     
385     for( Int_t ipoint=3;ipoint>=0;ipoint--){
386       track->fHit[ipoint] = AliHLTMUONConstants::NilRecHitStruct();
387       hitset[ipoint] = false;
388       
389       if(fBackTrackSeg[ibackTrackSeg].fIndex[ipoint]!=-1){
390
391         track->fHit[ipoint] = *(fChPoint[ipoint+6][fBackTrackSeg[ibackTrackSeg].fIndex[ipoint]]);
392         hitset[ipoint] = true;
393         
394       }
395     }
396     AliHLTMUONParticleSign sign = AliHLTMUONParticleSign(Int_t(TMath::Sign(1.,fTrackParam[ibackTrackSeg].GetInverseBendingMomentum())));
397     track->fFlags = AliHLTMUONUtils::PackMansoTrackFlags(sign,hitset);
398                                                          
399     
400     track++;
401   }
402
403   return true;
404 }
405
406 ///__________________________________________________________________________
407
408 Bool_t AliHLTMUONFullTracker::SlatTrackSeg()
409 {
410
411   ///Find the Slat Track Segments
412
413   Float_t trigX1,trigY1,trigZ1=AliMUONConstants::DefaultChamberZ(10);
414   Float_t trigX2,trigY2,trigZ2=AliMUONConstants::DefaultChamberZ(12);
415   Float_t extrapCh9X,extrapCh9Y,extrapCh9Z=AliMUONConstants::DefaultChamberZ(9);
416   Float_t extrapCh8X,extrapCh8Y,extrapCh8Z=AliMUONConstants::DefaultChamberZ(8);
417   Float_t extrapCh7X,extrapCh7Y,extrapCh7Z=AliMUONConstants::DefaultChamberZ(7);
418   Float_t extrapCh6X,extrapCh6Y,extrapCh6Z=AliMUONConstants::DefaultChamberZ(6);
419
420   Float_t meanX1,meanX2,meanY1,meanY2,meanZ1,meanZ2;
421
422   Double_t distChFront,distChBack;
423   Int_t nofFrontChPoints,nofBackChPoints;
424   Int_t frontIndex[fgkMaxNofTracks], backIndex[fgkMaxNofTracks];
425   
426   AliHLTMUONRecHitStruct p2,p3,pSeg1,pSeg2,pSeg3;
427   Double_t anglediff,anglediff1,anglediff2;
428   Double_t minAngle = 2.0;
429   
430   Bool_t St5TrackletFound = false;
431   Bool_t Ch9PointFound = false;
432   Bool_t Ch8PointFound = false;
433   Bool_t St4TrackletFound = false;
434   Bool_t Ch7PointFound = false;
435   Bool_t Ch6PointFound = false;
436
437   Int_t index1,index2,index3,index4;
438   IntPair Cells[2][fgkMaxNofTracks]; ///cell array  for 5 stn for given trigger
439
440
441   Float_t maxXDeflectionExtrap = 10.0 + 4.0;  ///simulation result 10.0
442   Float_t extrapRatio = 0.2;            ///simulation result 0.2  
443   Float_t circularWindow = 20.0;        ///simulatiuon result 20
444   Float_t minAngleWindow = 2.0 + 1.0;        ///simulation result 2.0
445
446   Float_t tx,ty;
447
448   AliHLTUInt16_t detElemID,prevDetElemID;
449   AliHLTUInt8_t chamber;
450   
451   Int_t minTrgCh,maxTrgCh;
452   for( Int_t itrig=0;itrig<fNofPoints[10];itrig++){
453     
454     
455     St5TrackletFound = false;
456     Ch9PointFound = false;
457     Ch8PointFound = false;
458
459     St4TrackletFound = false;
460     Ch7PointFound = false;
461     Ch6PointFound = false;
462
463     fOverflowed = false;
464
465     minTrgCh = -1;
466     maxTrgCh = -1;
467
468     ///for( Int_t istn=0;istn<5;istn++)
469     fNofCells[0] = 0;
470     fNofCells[1] = 0;
471     for( Int_t ihit=0;ihit<4;ihit++){
472
473       AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[ihit]).fFlags,chamber,detElemID);
474
475       if(ihit==0 and detElemID!=0)
476         minTrgCh = ihit;
477       if(ihit==1 and detElemID!=0 and prevDetElemID==0)
478         minTrgCh = ihit;
479       if(ihit==2 and detElemID!=0)
480         maxTrgCh = ihit;
481       if(ihit==3 and detElemID!=0 and prevDetElemID==0)
482         maxTrgCh = ihit;
483
484       prevDetElemID = detElemID;
485     }
486
487     if(minTrgCh == -1 or maxTrgCh == -1)
488       continue;
489
490     AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[minTrgCh]).fFlags,chamber,detElemID);
491     fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,trigZ1);
492     AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[maxTrgCh]).fFlags,chamber,detElemID);
493     fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,trigZ2);
494     
495     
496     trigX1 = (fChPoint11[itrig]->fHit[minTrgCh]).fX;
497     trigY1 = (fChPoint11[itrig]->fHit[minTrgCh]).fY;
498
499     trigX2 = (fChPoint11[itrig]->fHit[maxTrgCh]).fX;
500     trigY2 = (fChPoint11[itrig]->fHit[maxTrgCh]).fY;
501     
502
503 #ifdef PRINT_BACK
504     Float_t tz=0;    
505     HLTInfo("AliHLTMUONFullTracker::SlatTrackSeg()--Begin\n\n");
506     AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[minTrgCh]).fFlags,chamber,detElemID);
507     fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,tz);
508     HLTInfo("\ttrig 1 : (%f,%f,%f), org : %f\n",trigX1,trigY1,trigZ1,tz);
509     AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[maxTrgCh]).fFlags,chamber,detElemID);
510     fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,trigz);
511     HLTInfo("\ttrig 2 : (%f,%f,%f), org : %f\n",trigX2,trigY2,trigZ2,tz);
512 #endif
513
514     /////////////////////////////////////////////////// Stn 5///////////////////////////////////////////////////////////////
515
516     ///extrapCh9X = trigX1 + (trigX2-trigX1) * (extrapCh9Z-trigZ1)/(trigZ2 - trigZ1) ;
517     ///     extrapCh9Yref = trigY1 * extrapCh9Z/trigZ1 ;
518
519     nofFrontChPoints = 0; nofBackChPoints = 0;
520     for( Int_t ipointch9=0;ipointch9<fNofPoints[9];ipointch9++){
521       AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[9][ipointch9]->fFlags,chamber,detElemID);
522       fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh9Z);
523       
524       extrapCh9X = trigX1 * extrapCh9Z/trigZ1 ;
525       extrapCh9Y = trigY1 + (trigY2-trigY1) * (extrapCh9Z-trigZ1)/(trigZ2 - trigZ1) ;
526       
527 #ifdef PRINT_BACK
528       HLTInfo("\textrap9 : (%f,%f,%f)\n",extrapCh9X,extrapCh9Y,extrapCh9Z);
529 #endif    
530       
531       if(nofBackChPoints < (fgkMaxNofTracks-1) && 
532          fabsf(extrapCh9X-fChPoint[9][ipointch9]->fX) < maxXDeflectionExtrap && 
533          fabsf(extrapCh9Y-fChPoint[9][ipointch9]->fY)/
534          ((fChPoint[9][ipointch9]->fX * fChPoint[9][ipointch9]->fX) + 
535           (fChPoint[9][ipointch9]->fY * fChPoint[9][ipointch9]->fY)) <= extrapRatio ){
536         
537         distChBack = sqrt((extrapCh9X-fChPoint[9][ipointch9]->fX)*(extrapCh9X-fChPoint[9][ipointch9]->fX) 
538                           + (extrapCh9Y-fChPoint[9][ipointch9]->fY)*(extrapCh9Y-fChPoint[9][ipointch9]->fY));
539         if(distChBack>circularWindow) continue;
540
541 #ifdef PRINT_BACK
542         HLTInfo("\t\t\tCh9  :%lf (%f,%f,%f) , circularWindow : %f, distChBack : %f\n",
543                fabsf(extrapCh9X-fChPoint[9][ipointch9]->fX)/distChBack,
544                fChPoint[9][ipointch9]->fX,fChPoint[9][ipointch9]->fY,fChPoint[9][ipointch9]->fZ,
545                circularWindow,distChBack);
546 #endif
547
548         backIndex[nofBackChPoints++] = ipointch9;
549       }///if point found
550       
551     }/// ch10 loop
552     
553     for( Int_t ipointch8=0;ipointch8<fNofPoints[8];ipointch8++){
554       AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[8][ipointch8]->fFlags,chamber,detElemID);
555       fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh8Z);
556
557       ///     extrapCh8Yref = trigY1 * extrapCh8Z/trigZ1 ;
558       extrapCh8X = trigX1 * extrapCh8Z/trigZ1 ;
559       extrapCh8Y = trigY1 + (trigY2-trigY1) * (extrapCh8Z-trigZ1)/(trigZ2 - trigZ1) ;
560
561       
562 #ifdef PRINT_BACK
563       HLTInfo("\textrap8 :  (%f,%f,%f), extrapRatio : %lf, maxXDeflectionExtrap : %lf\n",
564              extrapCh8X,extrapCh8Y,extrapCh8Z,extrapRatio,maxXDeflectionExtrap);
565 #endif    
566       
567       if( nofFrontChPoints < (fgkMaxNofTracks-1) &&
568           fabsf(extrapCh8X-fChPoint[8][ipointch8]->fX) < maxXDeflectionExtrap && 
569           fabsf(extrapCh8Y-fChPoint[8][ipointch8]->fY)/
570           ((fChPoint[8][ipointch8]->fX * fChPoint[8][ipointch8]->fX) + 
571            (fChPoint[8][ipointch8]->fY * fChPoint[8][ipointch8]->fY)) <= extrapRatio ){
572         
573         distChFront = sqrt((extrapCh8X-fChPoint[8][ipointch8]->fX)*(extrapCh8X-fChPoint[8][ipointch8]->fX) 
574                            + (extrapCh8Y-fChPoint[8][ipointch8]->fY)*(extrapCh8Y-fChPoint[8][ipointch8]->fY));
575         
576 #ifdef PRINT_BACK
577         HLTInfo("\t\t\tCh8  : circularWindow : %f, distChFront : %f\n",circularWindow,distChFront);
578 #endif
579         
580         if(distChFront>circularWindow) continue;
581         
582 #ifdef PRINT_BACK
583         HLTInfo("\t\t\tCh8  :%lf (%f,%f,%f)\n",distChFront,
584                fChPoint[8][ipointch8]->fX,fChPoint[8][ipointch8]->fY,fChPoint[8][ipointch8]->fZ);
585 #endif
586         
587         frontIndex[nofFrontChPoints++] = ipointch8;
588       }///if point found
589       
590     }/// ch9 loop
591     
592     minAngle = minAngleWindow;
593     p3.fX = trigX1 ; p3.fY = trigY1 ; p3.fZ = trigZ1 ;
594     for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
595       Sub(&p3,fChPoint[9][backIndex[ibackpoint]],&pSeg2);
596       for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
597         Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[8][frontIndex[ifrontpoint]],&pSeg1);
598         anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2);
599 #ifdef PRINT_BACK
600         HLTInfo("\t\ttracklet-check-St5 : anglediff : %lf, minAngle : %lf\n",anglediff,minAngle);
601 #endif  
602         if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
603           St5TrackletFound = true;
604           Cells[1][fNofCells[1]].fFirst =  frontIndex[ifrontpoint];
605           Cells[1][fNofCells[1]].fSecond =  backIndex[ibackpoint];
606           fNofCells[1]++ ;
607 #ifdef PRINT_BACK
608           HLTInfo("\t\ttracklet-St5 : anglediff : %lf\n",anglediff);
609           HLTInfo("\t\t\tCh9  : (%f,%f,%f)\n",fChPoint[9][backIndex[ibackpoint]]->fX,
610                  fChPoint[9][backIndex[ibackpoint]]->fY,fChPoint[9][backIndex[ibackpoint]]->fZ);
611           HLTInfo("\t\t\tCh8  : (%f,%f,%f)\n",fChPoint[8][frontIndex[ifrontpoint]]->fX,
612                  fChPoint[8][frontIndex[ifrontpoint]]->fY,fChPoint[8][frontIndex[ifrontpoint]]->fZ);
613 #endif
614         }///anglediff condition
615       }///front
616     }///back
617
618     
619
620     
621     /// If tracklet not found, search for the single space point in Ch9 or in Ch8
622     if(!St5TrackletFound){
623       
624       minAngle = minAngleWindow; 
625       p3.fX = trigX2 ; p3.fY = trigY2 ; p3.fZ = trigZ2 ;
626       p2.fX = trigX1 ; p2.fY = trigY1 ; p2.fZ = trigZ1 ;
627       Sub(&p3,&p2,&pSeg2);
628       
629       ///Search in Ch9
630       for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
631         Sub(&p2,fChPoint[9][backIndex[ibackpoint]],&pSeg1);
632         anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
633         if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
634           Ch9PointFound = true;
635           Cells[1][fNofCells[1]].fFirst =  -1;
636           Cells[1][fNofCells[1]].fSecond =  backIndex[ibackpoint];
637           fNofCells[1]++ ;
638         }
639         
640 #ifdef PRINT_BACK
641         HLTInfo("\t\tsppoint-Ch9 : anglediff : %lf\n\n",anglediff);
642 #endif
643       }///back
644       
645       ///Search in Ch8
646       for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
647         Sub(&p2,fChPoint[8][frontIndex[ifrontpoint]],&pSeg1);
648         anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
649         if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
650           Ch8PointFound = true;
651           Cells[1][fNofCells[1]].fFirst = frontIndex[ifrontpoint];
652           Cells[1][fNofCells[1]].fSecond =  -1;
653           fNofCells[1]++ ;
654         }
655         
656 #ifdef PRINT_BACK
657         HLTInfo("\t\tsppoint-Ch8 : anglediff : %lf\n\n",anglediff);
658 #endif
659       }///front
660     }///if no tracklets found condition
661     
662 #ifdef PRINT_BACK
663     HLTInfo("\tnofTracks found after stn 5 : %d\n",fNofCells[1]);
664 #endif
665     
666     if(!St5TrackletFound && !Ch9PointFound && !Ch8PointFound) continue;
667
668     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
669     
670
671     /////////////////////////////////////////////////// Stn 4///////////////////////////////////////////////////////////////
672
673     /// extrapCh7X = trigX1 * extrapCh7Z/trigZ1 ;
674     /// extrapCh7Y = trigY1 + (trigY2-trigY1) * (extrapCh7Z-trigZ1)/(trigZ2 - trigZ1) ;
675     
676     /// extrapCh6X = trigX1 * extrapCh6Z/trigZ1 ;
677     /// extrapCh6Y = trigY1 + (trigY2-trigY1) * (extrapCh6Z-trigZ1)/(trigZ2 - trigZ1) ;
678     
679     
680     nofFrontChPoints = 0; nofBackChPoints = 0;
681     for( Int_t ipointch7=0;ipointch7<fNofPoints[7];ipointch7++){
682
683       ///     distChBack = (fChPoint[7][ipointch7]->fX * fChPoint[7][ipointch7]->fX) 
684       ///       + (fChPoint[7][ipointch7]->fY * fChPoint[7][ipointch7]->fY);
685       AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[7][ipointch7]->fFlags,chamber,detElemID);
686       fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh7Z);
687       extrapCh7X = trigX1 * extrapCh7Z/trigZ1 ;
688       extrapCh7Y = trigY1 + (trigY2-trigY1) * (extrapCh7Z-trigZ1)/(trigZ2 - trigZ1) ;
689       
690 #ifdef PRINT_BACK
691       HLTInfo("\textrap7 : (%f,%f,%f)\n",extrapCh7X,extrapCh7Y,extrapCh7Z);
692 #endif    
693
694       if( nofBackChPoints < (fgkMaxNofTracks-1) &&
695           fabsf(extrapCh7X-fChPoint[7][ipointch7]->fX) < maxXDeflectionExtrap && 
696           fabsf(extrapCh7Y-fChPoint[7][ipointch7]->fY)/
697           ((fChPoint[7][ipointch7]->fX * fChPoint[7][ipointch7]->fX) + 
698            (fChPoint[7][ipointch7]->fY * fChPoint[7][ipointch7]->fY)) <= extrapRatio ){
699         
700         distChBack = sqrt((extrapCh7X-fChPoint[7][ipointch7]->fX)*(extrapCh7X-fChPoint[7][ipointch7]->fX) 
701                           + (extrapCh7Y-fChPoint[7][ipointch7]->fY)*(extrapCh7Y-fChPoint[7][ipointch7]->fY));
702         
703         if(distChBack>circularWindow) continue;
704         
705 #ifdef PRINT_BACK
706         HLTInfo("\t\t\tCh7  :%lf (%f,%f,%f)\n",distChBack,
707                fChPoint[7][ipointch7]->fX,fChPoint[7][ipointch7]->fY,fChPoint[7][ipointch7]->fZ);
708 #endif
709         
710         backIndex[nofBackChPoints++] = ipointch7;
711       }///if point found
712       
713     }///ch8 loop
714     
715     for( Int_t ipointch6=0;ipointch6<fNofPoints[6];ipointch6++){
716
717       AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[6][ipointch6]->fFlags,chamber,detElemID);
718       fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh6Z);        
719       
720       extrapCh6X = trigX1 * extrapCh6Z/trigZ1 ;
721       extrapCh6Y = trigY1 + (trigY2-trigY1) * (extrapCh6Z-trigZ1)/(trigZ2 - trigZ1) ;
722
723 #ifdef PRINT_BACK
724       HLTInfo("\textrap6 :  (%f,%f,%f)\n",extrapCh6X,extrapCh6Y,extrapCh6Z);
725 #endif    
726
727       
728       if(nofFrontChPoints < (fgkMaxNofTracks-1) && 
729          fabsf(extrapCh6X-fChPoint[6][ipointch6]->fX) < maxXDeflectionExtrap && 
730          fabsf(extrapCh6Y-fChPoint[6][ipointch6]->fY)/
731          ((fChPoint[6][ipointch6]->fX * fChPoint[6][ipointch6]->fX) + 
732           (fChPoint[6][ipointch6]->fY * fChPoint[6][ipointch6]->fY)) <= extrapRatio ){
733         
734         distChFront = sqrt((extrapCh6X-fChPoint[6][ipointch6]->fX)*(extrapCh6X-fChPoint[6][ipointch6]->fX) 
735                            + (extrapCh6Y-fChPoint[6][ipointch6]->fY)*(extrapCh6Y-fChPoint[6][ipointch6]->fY));
736         if(distChFront>circularWindow) continue;
737         
738 #ifdef PRINT_BACK
739         HLTInfo("\t\t\tCh6  :%lf (%f,%f,%f)\n",distChFront,
740                fChPoint[6][ipointch6]->fX,fChPoint[6][ipointch6]->fY,fChPoint[6][ipointch6]->fZ);
741 #endif
742
743         frontIndex[nofFrontChPoints++] = ipointch6;
744       }///if point found
745       
746     }/// ch7 loop
747     
748     minAngle = minAngleWindow;
749     p3.fX = trigX1 ; p3.fY = trigY1 ; p3.fZ = trigZ1 ;
750     for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
751       Sub(&p3,fChPoint[7][backIndex[ibackpoint]],&pSeg2);
752       for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
753         Sub(fChPoint[7][backIndex[ibackpoint]],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
754         anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
755 #ifdef PRINT_BACK
756         HLTInfo("\t\ttracklet-check-St4 : anglediff : %lf, minAngle : %lf\n",anglediff,minAngle);
757 #endif  
758         if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
759
760           St4TrackletFound = true;
761
762           Cells[0][fNofCells[0]].fFirst =  frontIndex[ifrontpoint];
763           Cells[0][fNofCells[0]].fSecond =  backIndex[ibackpoint];
764           fNofCells[0]++ ;
765 #ifdef PRINT_BACK
766           HLTInfo("\t\ttracklet-St4 : anglediff : %lf\n",anglediff);
767           HLTInfo("\t\t\tCh7  : (%f,%f,%f)\n",fChPoint[7][backIndex[ibackpoint]]->fX,
768                  fChPoint[7][backIndex[ibackpoint]]->fY,fChPoint[7][backIndex[ibackpoint]]->fZ);
769           HLTInfo("\t\t\tCh6  : (%f,%f,%f)\n",fChPoint[6][frontIndex[ifrontpoint]]->fX,
770                  fChPoint[6][frontIndex[ifrontpoint]]->fY,fChPoint[6][frontIndex[ifrontpoint]]->fZ);
771 #endif
772
773         }///anglediff condn
774       }///front
775     }///back
776
777     
778
779     
780     /// If tracklet not found search for the single space point in Ch7 or in Ch6
781     if(!St4TrackletFound){
782
783       minAngle = minAngleWindow; 
784       p3.fX = trigX2 ; p3.fY = trigY2 ; p3.fZ = trigZ2 ;
785       p2.fX = trigX1 ; p2.fY = trigY1 ; p2.fZ = trigZ1 ;
786       Sub(&p3,&p2,&pSeg2);
787
788       ///Search in Ch7
789       for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
790         Sub(&p2,fChPoint[7][backIndex[ibackpoint]],&pSeg1);
791         
792         anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
793         if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
794           Ch7PointFound = true;
795           Cells[0][fNofCells[0]].fFirst =  -1;
796           Cells[0][fNofCells[0]].fSecond =  backIndex[ibackpoint];
797           fNofCells[0]++ ;
798         }
799         
800 #ifdef PRINT_BACK
801         HLTInfo("\t\tsppoint-Ch7 : anglediff : %lf\n\n",anglediff);
802 #endif
803       }///back
804       
805       ///Search in Ch6
806       for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
807         Sub(&p2,fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
808         anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
809         if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
810           Ch6PointFound = true;
811           Cells[0][fNofCells[0]].fFirst = frontIndex[ifrontpoint];
812           Cells[0][fNofCells[0]].fSecond =  -1;
813           fNofCells[0]++ ;
814         }
815         
816 #ifdef PRINT_BACK
817         HLTInfo("\t\tsppoint-Ch6 : anglediff : %lf\n\n",anglediff);
818 #endif
819       }///front
820     }///if no tracklets found condition
821     
822 #ifdef PRINT_BACK
823     HLTInfo("\tnofTracks found after stn 4 : %d\n",fNofCells[0]);
824 #endif
825     
826     if(!St4TrackletFound && !Ch7PointFound && !Ch6PointFound) continue;
827     
828     /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
829     ;
830       
831     ////////////////////////////////////////////// Analyse and fill trackseg array////////////////////////////////////////
832     ;
833
834     if(St5TrackletFound && St4TrackletFound){
835
836       minAngle = minAngleWindow;
837  
838       for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
839         index1 = Cells[0][itrackletfront].fFirst ;
840         index2 = Cells[0][itrackletfront].fSecond ;
841         Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
842         for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
843           index3 = Cells[1][itrackletback].fFirst ;
844           index4 = Cells[1][itrackletback].fSecond ;
845           Sub(fChPoint[8][index3],fChPoint[7][index2],&pSeg2);
846           Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
847           anglediff = Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3);
848           if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
849             fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
850             fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
851             fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
852             fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
853             fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
854             minAngle = anglediff;
855 #ifdef PRINT_BACK
856             HLTInfo("\t\ttracklet-St4 and St5 : anglediff : %lf\n",anglediff);
857             HLTInfo("\t\t\tCh9  : (%f,%f,%f)\n",fChPoint[9][index4]->fX,
858                    fChPoint[9][index4]->fY,fChPoint[9][index4]->fZ);
859             HLTInfo("\t\t\tCh8  : (%f,%f,%f)\n",fChPoint[8][index3]->fX,
860                    fChPoint[8][index3]->fY,fChPoint[8][index3]->fZ);
861             HLTInfo("\t\t\tCh7  : (%f,%f,%f)\n",fChPoint[7][index2]->fX,
862                    fChPoint[7][index2]->fY,fChPoint[7][index2]->fZ);
863             HLTInfo("\t\t\tCh6  : (%f,%f,%f)\n",fChPoint[6][index1]->fX,
864                    fChPoint[6][index1]->fY,fChPoint[6][index1]->fZ);
865 #endif
866
867           }///if minangle
868         }///for of front ch
869       }///for loop of back ch
870       
871       fNofbackTrackSeg++;
872
873     }else if(St5TrackletFound && (Ch7PointFound || Ch6PointFound)){
874       
875       
876       nofFrontChPoints = 0; nofBackChPoints = 0;
877       for( Int_t ifrontpoint=0;ifrontpoint<fNofCells[0];ifrontpoint++){
878         if(Cells[0][ifrontpoint].fFirst==-1 && nofBackChPoints<(fgkMaxNofTracks-1))
879           backIndex[nofBackChPoints++] = Cells[0][ifrontpoint].fSecond;
880         else if(nofFrontChPoints<(fgkMaxNofTracks-1))
881           frontIndex[nofFrontChPoints++] = Cells[0][ifrontpoint].fFirst; 
882       }
883       
884       minAngle = minAngleWindow;
885       if(nofFrontChPoints>0 && nofBackChPoints>0){
886
887         for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
888           index3 = Cells[1][itrackletback].fFirst ;
889           index4 = Cells[1][itrackletback].fSecond ;
890           Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
891           for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
892             Sub(fChPoint[8][index3],fChPoint[7][backIndex[ibackpoint]],&pSeg2);
893             for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
894               Sub(fChPoint[7][backIndex[ibackpoint]],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
895
896               anglediff = TMath::RadToDeg()*(Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3))/2.0;
897
898               if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
899                 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
900                 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
901                 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
902                 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
903                 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
904                 minAngle = anglediff;
905                 continue;
906               }
907               
908               Sub(fChPoint[8][index3],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
909               anglediff1 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg3);
910               anglediff2 = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
911               
912               if( anglediff1 < anglediff2 && anglediff1<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
913                 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
914                 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = -1;
915                 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
916                 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
917                 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
918                 minAngle = anglediff1;
919                 continue;
920               }
921               
922               if( anglediff2 < anglediff1 && anglediff2<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
923                 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = -1;
924                 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
925                 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
926                 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
927                 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
928                 minAngle = anglediff2;
929               }
930               
931             }///loop of ifrontpoint
932           }///loop of ibackpoint
933         }/// for loop of St5 cells
934       }else if(nofFrontChPoints>0){
935         
936         for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
937           index3 = Cells[1][itrackletback].fFirst ;
938           index4 = Cells[1][itrackletback].fSecond ;
939           Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
940
941           for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
942             Sub(fChPoint[8][index3],fChPoint[6][frontIndex[ifrontpoint]],&pSeg2);
943             
944             anglediff = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
945             if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
946               fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
947               fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = -1;
948               fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
949               fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
950               fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
951               minAngle = anglediff;
952             }///if anglediff
953           }///backch loop
954         }///tracklet loop
955
956       }else{ /// if(nofBackChPoints>0){
957         for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
958           index3 = Cells[1][itrackletback].fFirst ;
959           index4 = Cells[1][itrackletback].fSecond ;
960           Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
961           
962           for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
963             Sub(fChPoint[8][index3],fChPoint[7][backIndex[ibackpoint]],&pSeg2);
964             
965             anglediff = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
966             
967             if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
968               fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = -1;
969               fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
970               fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
971               fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
972               fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
973               minAngle = anglediff;
974             }///if anglediff
975           }///backch loop
976         }///tracklet loop
977
978       }///condn for if(nofFrontChPoints>0)
979
980       fNofbackTrackSeg++;
981
982     }else if((Ch9PointFound || Ch8PointFound) && St4TrackletFound){
983
984       nofFrontChPoints = 0; nofBackChPoints = 0;
985       for( Int_t ibackpoint=0;ibackpoint<fNofCells[1];ibackpoint++){
986         if(Cells[1][ibackpoint].fFirst==-1 && nofBackChPoints<(fgkMaxNofTracks-1))
987           backIndex[nofBackChPoints++] = Cells[1][ibackpoint].fSecond;
988         else if(nofFrontChPoints<(fgkMaxNofTracks-1))
989           frontIndex[nofFrontChPoints++] = Cells[1][ibackpoint].fFirst; 
990       }
991       
992       minAngle = minAngleWindow;
993       if(nofFrontChPoints>0 && nofBackChPoints>0){
994
995         for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
996           index1 = Cells[0][itrackletfront].fFirst ;
997           index2 = Cells[0][itrackletfront].fSecond ;
998           
999           Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1000           
1001           for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1002             
1003
1004             Sub(fChPoint[8][frontIndex[ifrontpoint]],fChPoint[7][index2],&pSeg2);
1005             
1006             for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1007               
1008               Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[8][frontIndex[ifrontpoint]],&pSeg3);
1009               
1010               anglediff = TMath::RadToDeg()*(Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3))/2.0;
1011               if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1012                 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1013                 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1014                 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1015                 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1016                 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1017                 minAngle = anglediff;
1018                 continue;
1019               }
1020               
1021               Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[7][index2],&pSeg3);
1022               
1023               anglediff1 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
1024               anglediff2 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg3);
1025               
1026               if( anglediff1 < anglediff2 && anglediff1<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1027                 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1028                 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1029                 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1030                 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = -1;
1031                 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1032                 minAngle = anglediff1;
1033                 continue;
1034               }
1035               
1036               if( anglediff2 < anglediff1 && anglediff2<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1037                 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1038                 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1039                 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = -1;
1040                 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1041                 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1042                 minAngle = anglediff2;
1043               }
1044               
1045             }///loop of ifrontpoint
1046           }///loop of ibackpoint
1047         }/// for loop of St5 cells
1048       }else if(nofFrontChPoints>0){
1049
1050         for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1051           index1 = Cells[0][itrackletfront].fFirst ;
1052           index2 = Cells[0][itrackletfront].fSecond ;
1053           
1054           Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1055           
1056           for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1057             
1058             Sub(fChPoint[8][frontIndex[ifrontpoint]],fChPoint[7][index2],&pSeg2);
1059             
1060             anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
1061             if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1062               fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1063               fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;;
1064               fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1065               fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = -1;
1066               fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1067               minAngle = anglediff;
1068             }///if anglediff
1069           }///point loop
1070         }///tracklet loop
1071
1072       }else{ /// if(nofBackChPoints>0){
1073
1074         for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1075           index1 = Cells[0][itrackletfront].fFirst ;
1076           index2 = Cells[0][itrackletfront].fSecond ;
1077
1078           Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1079
1080           for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1081             Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[6][index1],&pSeg3);
1082             anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg3);
1083             if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1084               fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1085               fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1086               fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = -1;
1087               fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1088               fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1089               minAngle = anglediff;
1090             }
1091           }///backch loop
1092         }///tracklet loop
1093
1094       }///condn for if(nofFrontChPoints>0)
1095
1096       fNofbackTrackSeg++;
1097       
1098     }else if((Ch9PointFound || Ch8PointFound) && (Ch7PointFound || Ch6PointFound)){
1099       ///To Do :  To be analysed for two points out of four slat chambers
1100     }
1101     
1102     /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1103 #ifdef PRINT_BACK
1104     HLTInfo("\n");
1105 #endif
1106
1107   }///trigger loop
1108
1109   for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
1110     
1111     
1112     if(fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[1]!=-1 ){
1113       meanX1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fX 
1114                 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fX)/2.0 ;
1115       meanY1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fY 
1116                 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fY)/2.0 ;
1117       meanZ1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fZ 
1118                 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fZ)/2.0 ;
1119     }else if(fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[1]==-1 ){
1120       meanX1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fX ;
1121       meanY1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fY ;
1122       meanZ1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fZ ;
1123     }else{
1124       meanX1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fX ;
1125       meanY1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fY ;
1126       meanZ1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fZ ;
1127     }
1128     
1129     if(fBackTrackSeg[ibacktrackseg].fIndex[2]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1 ){
1130       meanX2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fX 
1131                 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fX)/2.0 ;
1132       meanY2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fY 
1133                 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fY)/2.0 ;
1134       meanZ2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fZ 
1135                 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fZ)/2.0 ;
1136     }else if(fBackTrackSeg[ibacktrackseg].fIndex[2]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[3]==-1 ){
1137       meanX2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fX ;
1138       meanY2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fY ;
1139       meanZ2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fZ ;
1140     }else{
1141       meanX2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fX ;
1142       meanY2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fY ;
1143       meanZ2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fZ ;
1144     }
1145     
1146     fExtrapSt3X[ibacktrackseg] = meanX1 + (TrackDetCoordinate[2]-meanZ1)*(meanX2-meanX1)/(meanZ2-meanZ1);
1147     fExtrapSt3Y[ibacktrackseg] = meanY1 + (TrackDetCoordinate[2]-meanZ1)*(meanY2-meanY1)/(meanZ2-meanZ1);  
1148     fInclinationBack[ibacktrackseg] = (meanX2-meanX1)/(meanZ2-meanZ1) ;
1149     fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
1150   }///backtrigseg loop
1151
1152 #ifdef PRINT_BACK
1153   HLTInfo("AliHLTMUONFullTracker::SlatTrackSeg()--End\n");
1154   HLTInfo("\n\n");
1155 #endif
1156   
1157   return true;
1158 }
1159
1160
1161 ///__________________________________________________________________________
1162
1163 Bool_t AliHLTMUONFullTracker::QuadTrackSeg()
1164 {
1165   ///Find the Track Segments in the Quadrant chambers
1166   
1167   Float_t meanX1,meanX2,meanY1,meanY2,meanZ1,meanZ2;
1168   Float_t expectSt3X,expectSt3Y,inclinationFront;
1169
1170   AliHLTMUONRecHitStruct pSeg1,pSeg2,pSeg3;
1171   Double_t anglediff;///,anglediff1,anglediff2;
1172   Double_t minAngle = -1.0;
1173   Int_t  indexMinAngleFront = -1;
1174   Int_t  indexMinAngleBack = -1;
1175   Int_t backIndex = -1;
1176
1177   Int_t ch3CellPoint[fgkMaxNofCellsPerCh],ch2CellPoint[fgkMaxNofCellsPerCh],nofSt2Cells=0;
1178   Int_t ch1CellPoint[fgkMaxNofCellsPerCh],ch0CellPoint[fgkMaxNofCellsPerCh],nofSt1Cells=0;
1179   Bool_t isConnected[fgkMaxNofTracks];
1180
1181   Float_t distDiffX = 4.0;                    ///simulation result 4.0
1182   Float_t distDiffY = 10.0 ;                   ///simulation result 4.0
1183   ///float closeToY0 = 10.0;                    ///simulation result 10.0 
1184   Float_t minAngleWindow = 2.0 ;               ///simulation result 2.0
1185   Float_t diffDistStX = 25.0;                  ///simulation result 25.0
1186   Float_t diffDistStY = 75.0;                  ///simulation result 25.0
1187   Float_t st3WindowX = 40.0   ;                ///simulation result 40.0
1188   Float_t st3WindowY = 10.0;                   ///simulation result 10.0
1189   ///   Float_t inclinationWindow = 0.04;            ///inclination window   
1190   ///   Float_t st3WindowXOp2 = 40.0 ;                 ///simulation result 40.0
1191   ///   Float_t st3WindowYOp2 = 10.0;                ///simulation result 10.0
1192   
1193   
1194 #ifdef PRINT_FRONT
1195   HLTInfo("AliHLTMUONFullTracker::QuadTrackSeg()--Begin\n\n");
1196 #endif
1197   
1198
1199   for( Int_t ibackpoint=0;ibackpoint<fNofPoints[3];ibackpoint++){
1200     for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[2];ifrontpoint++){
1201
1202       if(fabsf(fChPoint[3][ibackpoint]->fX - fChPoint[2][ifrontpoint]->fX)<distDiffX 
1203          && fabsf(fChPoint[3][ibackpoint]->fY - fChPoint[2][ifrontpoint]->fY)<distDiffY){
1204
1205         /// if((fabsf(fChPoint[3][ibackpoint]->fY) > closeToY0 && 
1206         ///     fabsf(fChPoint[3][ibackpoint]->fY) < fabsf(fChPoint[2][ifrontpoint]->fY)) ||
1207         ///    nofSt2Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1208
1209         if(nofSt2Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1210         
1211 #ifdef PRINT_FRONT
1212         HLTInfo("\t\t\tCh3  : %d, (%f,%f,%f)\n",
1213                nofSt2Cells,fChPoint[3][ibackpoint]->fX,fChPoint[3][ibackpoint]->fY,fChPoint[3][ibackpoint]->fZ);
1214         HLTInfo("\t\t\tCh2  :(%f,%f,%f)\n\n",
1215                fChPoint[2][ifrontpoint]->fX,fChPoint[2][ifrontpoint]->fY,fChPoint[2][ifrontpoint]->fZ);
1216 #endif
1217         
1218         ch3CellPoint[nofSt2Cells] = ibackpoint;
1219         ch2CellPoint[nofSt2Cells] = ifrontpoint;
1220         nofSt2Cells++; 
1221         
1222       }///if point found
1223     }///frontch
1224   }///backch
1225   
1226   for( Int_t ibackpoint=0;ibackpoint<fNofPoints[1];ibackpoint++){
1227     for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[0];ifrontpoint++){
1228       if(fabsf(fChPoint[1][ibackpoint]->fX - fChPoint[0][ifrontpoint]->fX)< distDiffX
1229          && fabsf(fChPoint[1][ibackpoint]->fY - fChPoint[0][ifrontpoint]->fY)<distDiffY){
1230         
1231         /// if((fabsf(fChPoint[1][ibackpoint]->fY) > closeToY0 && 
1232         ///     fabsf(fChPoint[1][ibackpoint]->fY) < fabsf(fChPoint[0][ifrontpoint]->fY)) ||
1233         ///    nofSt1Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1234
1235         if(nofSt1Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1236            
1237         
1238 #ifdef PRINT_FRONT
1239         HLTInfo("\t\t\tCh1  : %d, (%f,%f,%f)\n",
1240                nofSt1Cells,fChPoint[1][ibackpoint]->fX,fChPoint[1][ibackpoint]->fY,fChPoint[1][ibackpoint]->fZ);
1241         HLTInfo("\t\t\tCh0  :(%f,%f,%f)\n\n",
1242                fChPoint[0][ifrontpoint]->fX,fChPoint[0][ifrontpoint]->fY,fChPoint[0][ifrontpoint]->fZ);
1243 #endif
1244         ch1CellPoint[nofSt1Cells] = ibackpoint;
1245         ch0CellPoint[nofSt1Cells] = ifrontpoint;
1246         nofSt1Cells++;
1247       }///if point found
1248     }///frontch
1249   }///backch
1250   
1251 #ifdef PRINT_FRONT
1252   HLTInfo("\tnofSt1Cells : %d, nofSt2Cells : %d\n",nofSt1Cells,nofSt2Cells);
1253 #endif
1254
1255   for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++)
1256     isConnected[ibacktrackseg] = false;
1257
1258
1259   ///First Check for Tracklets in two front stations
1260   for( Int_t itrackletfront=0;itrackletfront<nofSt1Cells;itrackletfront++){
1261     
1262     Sub(fChPoint[1][ch1CellPoint[itrackletfront]],fChPoint[0][ch0CellPoint[itrackletfront]],&pSeg1);
1263     
1264     minAngle = minAngleWindow;
1265     indexMinAngleBack = -1;
1266     indexMinAngleFront = -1;
1267
1268     meanX1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fX 
1269               + fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/2.0 ;
1270     meanY1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fY 
1271               + fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/2.0 ;
1272     meanZ1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fZ 
1273               + fChPoint[1][ch1CellPoint[itrackletfront]]->fZ)/2.0 ;
1274     
1275     for( Int_t itrackletback=0;itrackletback<nofSt2Cells;itrackletback++){
1276 #ifdef PRINT_FRONT
1277       cout<<"\tBefore "
1278           <<fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)
1279           <<"\t"
1280           <<fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)
1281           <<endl;
1282 #endif
1283       if(fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fX - 
1284                fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX || 
1285          fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fY - 
1286                fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1287
1288       meanX2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fX 
1289                 + fChPoint[3][ch3CellPoint[itrackletback]]->fX)/2.0 ;
1290       meanY2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fY 
1291                 + fChPoint[3][ch3CellPoint[itrackletback]]->fY)/2.0 ;
1292       meanZ2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fZ 
1293                 + fChPoint[3][ch3CellPoint[itrackletback]]->fZ)/2.0 ;
1294       
1295       expectSt3X = meanX2 + (TrackDetCoordinate[2]-meanZ2)*(meanX2-meanX1)/(meanZ2-meanZ1);
1296       expectSt3Y = meanY2 + (TrackDetCoordinate[2]-meanZ2)*(meanY2-meanY1)/(meanZ2-meanZ1);
1297       inclinationFront = (meanX2-meanX1)/(meanZ2-meanZ1) ;
1298       
1299       for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){  
1300         
1301         if( /// fabsf(inclinationBack[ibacktrackseg]-inclinationFront)<0.04 && 
1302            fabsf((expectSt3X-fExtrapSt3X[ibacktrackseg])) < st3WindowX
1303            && fabsf((expectSt3Y-fExtrapSt3Y[ibacktrackseg])) < st3WindowY){
1304           
1305           Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1306           Sub(fChPoint[3][ch3CellPoint[itrackletback]],fChPoint[2][ch2CellPoint[itrackletback]],&pSeg3);
1307           
1308           anglediff = TMath::RadToDeg()* (Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3));
1309
1310 #ifdef PRINT_FRONT
1311           HLTInfo("\t\t\tanglediff : %lf\n",anglediff);
1312 #endif    
1313           if(anglediff<minAngle){
1314             minAngle = anglediff;
1315             indexMinAngleBack = itrackletback;
1316             indexMinAngleFront = itrackletfront;
1317             backIndex = ibacktrackseg;
1318             isConnected[ibacktrackseg] = true;
1319           }
1320         }///matching tracklet
1321       }///for loop of backtrackseg
1322       
1323     }///for of back ch
1324     
1325     if(minAngle < minAngleWindow && indexMinAngleFront!=-1 
1326        && indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)){
1327       
1328       fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = ch0CellPoint[indexMinAngleFront];
1329       fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = ch1CellPoint[indexMinAngleFront];
1330       fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = ch2CellPoint[indexMinAngleBack];
1331       fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = ch3CellPoint[indexMinAngleBack];
1332       
1333       fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;
1334       fNoffrontTrackSeg++;
1335     }///condition to find valid tracklet
1336     
1337   }///for loop of front ch
1338   
1339   Int_t nofNCfBackTrackSeg = 0;
1340   Int_t ncfBackTrackSeg[fgkMaxNofTracks];
1341   
1342   for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++)
1343     if(!isConnected[ibacktrackseg]) 
1344       ncfBackTrackSeg[nofNCfBackTrackSeg++] = ibacktrackseg;
1345     else
1346       fNofConnected++;
1347   
1348
1349 #ifdef PRINT_FRONT
1350   HLTInfo("\tfNofConnected : %d, nofNCfBackTrackSeg : %d\n",fNofConnected,nofNCfBackTrackSeg);
1351   HLTInfo("\tfNofPoints[3] : %d, fNofPoints[2] : %d\n",fNofPoints[3],fNofPoints[2]);
1352   if(nofNCfBackTrackSeg==0){
1353     HLTInfo("All fBackTrackSegs are connected with fFrontTrackSegs, no need to search further\n");
1354     HLTInfo("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
1355   }
1356 #endif
1357   
1358   if(nofNCfBackTrackSeg==0) return true;
1359
1360
1361   ///Next Check for tracklet in Stn1 and space point in Stn2
1362   Bool_t isbackpoint=false,isfrontpoint=false;
1363   for( Int_t itrackletfront=0;itrackletfront<nofSt1Cells;itrackletfront++){
1364     Sub(fChPoint[1][ch1CellPoint[itrackletfront]],fChPoint[0][ch0CellPoint[itrackletfront]],&pSeg1);
1365     minAngle = minAngleWindow;
1366     indexMinAngleBack = -1;
1367     indexMinAngleFront = -1;
1368
1369     for( Int_t ibackpoint=0;ibackpoint<fNofPoints[3];ibackpoint++){
1370       if(/// hasCh3Cells[ibackpoint] == true && 
1371          fabsf(fChPoint[3][ibackpoint]->fX - 
1372                fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX || 
1373          fabsf(fChPoint[3][ibackpoint]->fY - 
1374                fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1375       
1376       expectSt3X = fChPoint[3][ibackpoint]->fX + (TrackDetCoordinate[2] - fChPoint[3][ibackpoint]->fZ)*
1377         (fChPoint[3][ibackpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1378         (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1379       expectSt3Y = fChPoint[3][ibackpoint]->fY + (TrackDetCoordinate[2] - fChPoint[3][ibackpoint]->fZ)*
1380         (fChPoint[3][ibackpoint]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/
1381         (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1382       inclinationFront = (fChPoint[3][ibackpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1383         (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ) ;
1384       
1385       for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++){  
1386         
1387         if(/// fabsf(inclinationBack[ncfBackTrackSeg[ibacktrackseg]]-inclinationFront)< inclinationWindow && 
1388            fabsf((expectSt3X-fExtrapSt3X[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1389            && fabsf((expectSt3Y-fExtrapSt3Y[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowY){
1390           
1391           Sub(fChPoint[3][ibackpoint],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1392           
1393           anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
1394 #ifdef PRINT_FRONT
1395           HLTInfo("\t\t annglediff(Ch4) : %lf\n",anglediff);
1396 #endif    
1397           if(anglediff<minAngle){
1398             minAngle = anglediff;
1399             indexMinAngleBack = ibackpoint;
1400             indexMinAngleFront = itrackletfront;
1401             backIndex = ibacktrackseg;
1402             isConnected[ncfBackTrackSeg[ibacktrackseg]] = true;
1403             isbackpoint = true;
1404             isfrontpoint = false;
1405           }///angle min condn
1406         }///matching tracklet
1407       }///loop on not found back trackseg
1408     }///backpoint loop
1409     
1410     for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[2];ifrontpoint++){
1411       if(/// hasCh2Cells[ifrontpoint] == true && 
1412          fabsf(fChPoint[2][ifrontpoint]->fX - 
1413                fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX || 
1414          fabsf(fChPoint[2][ifrontpoint]->fY - 
1415                fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1416       
1417       expectSt3X = fChPoint[2][ifrontpoint]->fX + (TrackDetCoordinate[2] - fChPoint[2][ifrontpoint]->fZ)*
1418         (fChPoint[2][ifrontpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1419         (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1420       expectSt3Y = fChPoint[2][ifrontpoint]->fY + (TrackDetCoordinate[2] - fChPoint[2][ifrontpoint]->fZ)*
1421         (fChPoint[2][ifrontpoint]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/
1422         (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1423       inclinationFront = (fChPoint[2][ifrontpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1424         (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ) ;
1425       
1426       for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++){  
1427         
1428         if( /// fabsf(inclinationBack[ncfBackTrackSeg[ibacktrackseg]]-inclinationFront)< inclinationWindow && 
1429            fabsf((expectSt3X-fExtrapSt3X[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1430            && fabsf((expectSt3Y-fExtrapSt3Y[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowY){
1431           
1432           Sub(fChPoint[2][ifrontpoint],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1433
1434           anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
1435 #ifdef PRINT_FRONT
1436           HLTInfo("\t\t annglediff(Ch3) : %lf\n",anglediff);
1437 #endif    
1438           if(anglediff<minAngle){
1439             minAngle = anglediff;
1440             indexMinAngleBack = ifrontpoint;
1441             indexMinAngleFront = itrackletfront;
1442             backIndex = ibacktrackseg;
1443             isConnected[ncfBackTrackSeg[ibacktrackseg]] = true;
1444             isbackpoint = false;
1445             isfrontpoint = true;
1446
1447           }///angle min condn
1448         }///matching tracklet
1449       }///loop on not found back trackseg
1450     }///backpoint loop
1451     
1452     if(minAngle < minAngleWindow && indexMinAngleFront!=-1 && 
1453        indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)){
1454       
1455       fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = ch0CellPoint[indexMinAngleFront];
1456       fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = ch1CellPoint[indexMinAngleFront];
1457       if(isfrontpoint && !isbackpoint){
1458         fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = indexMinAngleBack;
1459         fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = -1;
1460       }else{
1461         fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = -1;
1462         fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = indexMinAngleBack;
1463       }
1464       fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;      
1465       fNoffrontTrackSeg++;
1466     }///condition to find valid tracklet
1467     
1468   }///front ch
1469
1470
1471   Int_t nofSNCfBackTrackSeg = 0;
1472   Int_t sncfBackTrackSeg[fgkMaxNofTracks];
1473   
1474   for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++)
1475     if(!isConnected[ncfBackTrackSeg[ibacktrackseg]]) 
1476       sncfBackTrackSeg[nofSNCfBackTrackSeg++] = ncfBackTrackSeg[ibacktrackseg];
1477     else
1478       fNofConnected++;
1479
1480 #ifdef PRINT_FRONT
1481   HLTInfo("\tfNofConnected : %d, nofSNCfBackTrackSeg : %d\n",fNofConnected,nofSNCfBackTrackSeg);
1482   if(nofSNCfBackTrackSeg==0){
1483     HLTInfo("All fBackTrackSegs are connected with fFrontTrackSegs, no need to search further\n");
1484     HLTInfo("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
1485   }
1486 #endif
1487
1488   if(nofSNCfBackTrackSeg==0) return true;
1489
1490   ///Last Check for tracklet in Stn2 and space point in Stn1
1491   for( Int_t itrackletback=0;itrackletback<nofSt2Cells;itrackletback++){
1492     Sub(fChPoint[3][ch3CellPoint[itrackletback]],fChPoint[2][ch2CellPoint[itrackletback]],&pSeg1);
1493     minAngle = minAngleWindow ;
1494     indexMinAngleBack = -1;
1495     indexMinAngleFront = -1;
1496
1497     for( Int_t ibackpoint=0;ibackpoint<fNofPoints[1];ibackpoint++){
1498       if(/// hasCh1Cells[ibackpoint] == true && 
1499          fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fX - 
1500                fChPoint[1][ibackpoint]->fX) > diffDistStX || 
1501          fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fY - 
1502                fChPoint[1][ibackpoint]->fY) > diffDistStY) continue;
1503       
1504       expectSt3X = fChPoint[2][ch2CellPoint[itrackletback]]->fX + (TrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
1505         (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ibackpoint]->fX)/
1506         (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ);
1507       expectSt3Y = fChPoint[2][ch2CellPoint[itrackletback]]->fY + (TrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
1508         (fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[1][ibackpoint]->fY)/
1509         (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ);
1510       inclinationFront = (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ibackpoint]->fX)/
1511         (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ) ;
1512       
1513       for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++){  
1514         
1515         if(///  fabsf(inclinationBack[sncfBackTrackSeg[ibacktrackseg]]-inclinationFront)<inclinationWindow && 
1516            fabsf((expectSt3X-fExtrapSt3X[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1517            && fabsf((expectSt3Y-fExtrapSt3Y[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowY ){
1518           
1519           Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[1][ibackpoint],&pSeg2);
1520           
1521           anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
1522           if(anglediff<minAngle){
1523             minAngle = anglediff;
1524             indexMinAngleBack = itrackletback;
1525             indexMinAngleFront = ibackpoint;
1526             backIndex = ibacktrackseg;
1527             isConnected[sncfBackTrackSeg[ibacktrackseg]] = true;
1528             isbackpoint = true;
1529             isfrontpoint = false;
1530           }///angle min condn
1531         }///matching tracklet
1532       }///loop on not found back trackseg
1533     }///backpoint loop
1534     
1535     for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[0];ifrontpoint++){
1536       if(/// hasCh0Cells[ifrontpoint] == true &&
1537          fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fX - 
1538                fChPoint[0][ifrontpoint]->fX) > diffDistStX || 
1539          fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fY - 
1540                fChPoint[0][ifrontpoint]->fY) > diffDistStY ) continue;
1541       
1542       expectSt3X = fChPoint[2][ch2CellPoint[itrackletback]]->fX + (TrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
1543         (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[0][ifrontpoint]->fX)/
1544         (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ);
1545       expectSt3Y = fChPoint[2][ch2CellPoint[itrackletback]]->fY + (TrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
1546         (fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[0][ifrontpoint]->fY)/
1547         (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ);
1548       inclinationFront = (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[0][ifrontpoint]->fX)/
1549         (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ) ;
1550       
1551       for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++){  
1552
1553         if(///  fabsf(inclinationBack[sncfBackTrackSeg[ibacktrackseg]]-inclinationFront)<inclinationWindow && 
1554            fabsf((expectSt3X-fExtrapSt3X[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1555            && fabsf((expectSt3Y-fExtrapSt3Y[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowY ){
1556           
1557           Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[0][ifrontpoint],&pSeg2);
1558           
1559           anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2) ;
1560           if(anglediff<minAngle){
1561             minAngle = anglediff;
1562             indexMinAngleBack = itrackletback;
1563             indexMinAngleFront = ifrontpoint;
1564             backIndex = ibacktrackseg;
1565             isConnected[sncfBackTrackSeg[ibacktrackseg]] = true;
1566             isbackpoint = false;
1567             isfrontpoint = true;
1568             
1569           }///angle min condn
1570         }///matching tracklet
1571       }///loop on not found back trackseg
1572     }///backpoint loop
1573     
1574     if(minAngle < minAngleWindow && indexMinAngleFront!=-1 && 
1575        indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)){
1576
1577       if(isfrontpoint){      
1578         fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = indexMinAngleFront;
1579         fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = -1;
1580       }else{
1581         fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = -1;
1582         fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = indexMinAngleFront;
1583       }
1584
1585       fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = ch2CellPoint[indexMinAngleBack];
1586       fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = ch3CellPoint[indexMinAngleBack];
1587
1588       fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;
1589       fNoffrontTrackSeg++;
1590     }///condition to find valid tracklet
1591     
1592   }///front ch
1593   
1594   for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++)
1595     if(isConnected[sncfBackTrackSeg[ibacktrackseg]]) 
1596       fNofConnected++;
1597   
1598 #ifdef PRINT_FRONT
1599   HLTInfo("\tfNofConnected : %d\n",fNofConnected);
1600   HLTInfo("Three spacepoints are found in fFrontTrackSegs\n");
1601   HLTInfo("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
1602 #endif
1603
1604
1605   return true;
1606 }
1607
1608 ///__________________________________________________________________________
1609
1610 Double_t AliHLTMUONFullTracker::KalmanFilter(AliMUONTrackParam &trackParamAtCluster, Cluster *cluster)
1611 {
1612   //// Compute new track parameters and their covariances including new cluster using kalman filter
1613   //// return the additional track chi2
1614
1615 #ifdef PRINT_DETAIL_KALMAN
1616   HLTInfo("AliHLTMUONFullTracker::KalmanFilter()--Begin\n\n");
1617 #endif
1618
1619
1620   /// Get actual track parameters (p)
1621   TMatrixD param(trackParamAtCluster.GetParameters());
1622 #ifdef PRINT_DETAIL_KALMAN
1623   Info("\tKalmanFilter","param.Print() [p]");
1624   param.Print();
1625   HLTInfo("GetZ : %lf\n",trackParamAtCluster.GetZ());
1626 #endif
1627
1628
1629
1630   /// Get new cluster parameters (m)
1631   ///AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
1632   TMatrixD clusterParam(5,1);
1633   clusterParam.Zero();
1634   ///   clusterParam(0,0) = cluster->GetX();
1635   ///   clusterParam(2,0) = cluster->GetY();
1636   clusterParam(0,0) = cluster->fX;
1637   clusterParam(2,0) = cluster->fY;
1638 #ifdef PRINT_DETAIL_KALMAN
1639   Info("\tKalmanFilter","clusterParam.Print() [m]");
1640   clusterParam.Print();
1641 #endif
1642
1643
1644
1645   /// Compute the actual parameter weight (W)
1646   TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
1647 #ifdef PRINT_DETAIL_KALMAN
1648   Info("\tKalmanFilter","Covariance : [C]");
1649   paramWeight.Print();
1650 #endif
1651
1652   if (paramWeight.Determinant() != 0) {
1653     paramWeight.Invert();
1654   } else {
1655     Warning("KalmanFilter"," Determinant = 0");
1656     return 1.e10;
1657   }
1658
1659 #ifdef PRINT_DETAIL_KALMAN
1660   Info("\tKalmanFilter","Weight Matrix inverse of Covariance [W = C^-1]");
1661   paramWeight.Print();
1662 #endif
1663
1664
1665   /// Compute the new cluster weight (U)
1666   TMatrixD clusterWeight(5,5);
1667   clusterWeight.Zero();
1668   clusterWeight(0,0) = 1. / cluster->fErrX2;
1669   clusterWeight(2,2) = 1. / cluster->fErrY2;
1670 #ifdef PRINT_DETAIL_KALMAN
1671   Info("\tKalmanFilter","clusterWeight.Print() [U]");
1672   HLTInfo("\tErrX2 : %lf, ErrY2 : %lf\n",cluster->fErrX2,cluster->fErrY2);
1673   clusterWeight.Print();
1674 #endif
1675
1676
1677
1678
1679   /// Compute the new parameters covariance matrix ( (W+U)^-1 )
1680   TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
1681 #ifdef PRINT_DETAIL_KALMAN
1682   Info("\tKalmanFilter","newParamCov.Print() [(W+U)]");
1683   newParamCov.Print();
1684 #endif
1685   if (newParamCov.Determinant() != 0) {
1686     newParamCov.Invert();
1687   } else {
1688     Warning("RunKalmanFilter"," Determinant = 0");
1689     return 1.e10;
1690   }
1691 #ifdef PRINT_DETAIL_KALMAN
1692   Info("\tKalmanFilter","newParamCov.Print() [(W+U)^-1] (new covariances[W] for trackParamAtCluster)");
1693   newParamCov.Print();
1694 #endif
1695
1696   /// Save the new parameters covariance matrix
1697   trackParamAtCluster.SetCovariances(newParamCov);
1698   
1699   /// Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
1700   TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
1701   TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); /// U(m-p)
1702   TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); /// ((W+U)^-1)U(m-p)
1703   newParam += param; /// ((W+U)^-1)U(m-p) + p
1704 #ifdef PRINT_DETAIL_KALMAN
1705   Info("\tKalmanFilter","newParam.Print() [p' = ((W+U)^-1)U(m-p) + p] (new parameters[p] for trackParamAtCluster)");
1706   newParam.Print();
1707 #endif
1708
1709   /// Save the new parameters
1710   trackParamAtCluster.SetParameters(newParam);
1711   ///   HLTInfo(Form("Pt : %lf\n",TMath::Sqrt(trackParamAtCluster.Px()*trackParamAtCluster.Px() + 
1712   ///                                 trackParamAtCluster.Py()*trackParamAtCluster.Py())));
1713   
1714   
1715   /// Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
1716   tmp = newParam; /// p'
1717   tmp -= param; /// (p'-p)
1718   TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); /// W(p'-p)
1719   TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); /// ((p'-p)^-1)W(p'-p)
1720   tmp = newParam; /// p'
1721   tmp -= clusterParam; /// (p'-m)
1722   TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); /// U(p'-m)
1723   addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); /// ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
1724
1725 #ifdef PRINT_DETAIL_KALMAN
1726   Info("\tKalmanFilter","addChi2Track.Print() [additional chi2 = ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)))]");
1727   addChi2Track.Print();
1728   HLTInfo("AliHLTMUONFullTracker::KalmanFilter()--End\n\n");
1729 #endif
1730
1731
1732   return addChi2Track(0,0);
1733 }
1734
1735 ///__________________________________________________________________________
1736
1737 Double_t AliHLTMUONFullTracker::TryOneCluster(const AliMUONTrackParam &trackParam, Cluster* cluster,
1738                                                      AliMUONTrackParam &trackParamAtCluster, Bool_t updatePropagator)
1739 {
1740   //// Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
1741   //// return the corresponding Chi2
1742   //// return trackParamAtCluster
1743   
1744   /// extrapolate track parameters and covariances at the z position of the tested cluster
1745   trackParamAtCluster = trackParam;
1746   AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->fZ, updatePropagator);
1747   
1748   /// set pointer to cluster into trackParamAtCluster
1749   ///trackParamAtCluster.SetClusterPtr(cluster);
1750   
1751   /// Set differences between trackParam and cluster in the bending and non bending directions
1752   Double_t dX = cluster->fX - trackParamAtCluster.GetNonBendingCoor();
1753   Double_t dY = cluster->fY - trackParamAtCluster.GetBendingCoor();
1754   
1755   /// Calculate errors and covariances
1756   const TMatrixD& kParamCov = trackParamAtCluster.GetCovariances();
1757   Double_t sigmaX2 = kParamCov(0,0) + cluster->fErrX2;
1758   Double_t sigmaY2 = kParamCov(2,2) + cluster->fErrY2;
1759   
1760   /// Compute chi2
1761   return dX * dX / sigmaX2 + dY * dY / sigmaY2;
1762   
1763 }
1764
1765 ///__________________________________________________________________________
1766
1767 Bool_t AliHLTMUONFullTracker::TryOneClusterFast(const AliMUONTrackParam &trackParam, Cluster* cluster)
1768 {
1769   //// Test the compatibility between the track and the cluster
1770   //// given the track resolution + the maximum-distance-to-track value
1771   //// and assuming linear propagation of the track:
1772   //// return kTRUE if they are compatibles
1773   
1774   Float_t SigmaCutForTracking = 6.0;
1775   Float_t MaxNonBendingDistanceToTrack = 1.0;
1776   Float_t MaxBendingDistanceToTrack = 1.0;
1777   
1778   Double_t dZ = cluster->fZ - trackParam.GetZ();
1779   Double_t dX = cluster->fX - (trackParam.GetNonBendingCoor() + trackParam.GetNonBendingSlope() * dZ);
1780   Double_t dY = cluster->fY - (trackParam.GetBendingCoor() + trackParam.GetBendingSlope() * dZ);
1781   const TMatrixD& kParamCov = trackParam.GetCovariances();
1782   Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1);
1783   Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3);
1784
1785   Double_t dXmax = SigmaCutForTracking * TMath::Sqrt(errX2) +
1786     MaxNonBendingDistanceToTrack;
1787   Double_t dYmax = SigmaCutForTracking * TMath::Sqrt(errY2) +
1788     MaxBendingDistanceToTrack;
1789   
1790   if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE;
1791   
1792   return kTRUE;
1793   
1794 }
1795 ///__________________________________________________________________________
1796
1797 void AliHLTMUONFullTracker::PropagateTracks(Double_t charge, Float_t& px, Float_t& py, Float_t& pz, 
1798                                             Float_t& xr, Float_t& yr, Float_t& zr, Float_t zprop)
1799 {
1800   ///
1801   /// propagate in magnetic field between hits of indices i1 and i2
1802   ///
1803
1804   Double_t vect[7], vout[7];
1805   Double_t step = -5.0;
1806   Double_t zMax = zprop+.5;
1807   
1808   vect[0] = xr;
1809   vect[1] = yr;
1810   vect[2] = zr;
1811   vect[6] = sqrt((px)*(px) + (py)*(py) + (pz)*(pz));
1812   vect[3] = px/vect[6];
1813   vect[4] = py/vect[6];
1814   vect[5] = pz/vect[6];
1815   ///   cout<<"vec[2] : "<<vect[2]<<", zMax : "<<zMax<<endl;
1816
1817   Int_t nSteps = 0;
1818   while ((vect[2] < zMax) && (nSteps < 10000)) {
1819     nSteps++;
1820     ///     OneStepRungekutta(charge, step, vect, vout);
1821     OneStepHelix3(charge,step,vect,vout);
1822     ///SetPoint(fCount,vout[0],vout[1],vout[2]);
1823     ///fCount++;
1824     ///     HLTInfo("(x,y,z) : (%f,%f,%f)\n",vout[0],vout[1],vout[2]);
1825     for (Int_t i = 0; i < 7; i++) {
1826       vect[i] = vout[i];
1827     }
1828   }
1829   
1830   xr = vect[0];
1831   yr = vect[1];
1832   zr = vect[2];
1833   
1834   px = vect[3]*vect[6];
1835   py = vect[4]*vect[6];
1836   pz = vect[5]*vect[6];
1837   return;
1838 }
1839
1840 ///__________________________________________________________________________
1841
1842 void AliHLTMUONFullTracker::OneStepHelix3(Double_t field, Double_t step, Double_t *vect, Double_t *vout)
1843 {
1844   //// <pre>
1845   ////  ******************************************************************
1846   ////  *                                                                *
1847   ////  *       Tracking routine in a constant field oriented            *
1848   ////  *       along axis 3                                             *
1849   ////  *       Tracking is performed with a conventional                *
1850   ////  *       helix step method                                        *
1851   ////  *                                                                *
1852   ////  *    ==>Called by : USER, GUSWIM                                 *
1853   ////  *       Authors    R.Brun, M.Hansroul  *********                 *
1854   ////  *       Rewritten  V.Perevoztchikov
1855   ////  *                                                                *
1856   ////  ******************************************************************
1857   //// </pre>
1858
1859   Double_t hxp[3];
1860   Double_t h4, hp, rho, tet;
1861   Double_t sint, sintt, tsint, cos1t;
1862   Double_t f1, f2, f3, f4, f5, f6;
1863
1864   const Int_t kix  = 0;
1865   const Int_t kiy  = 1;
1866   const Int_t kiz  = 2;
1867   const Int_t kipx = 3;
1868   const Int_t kipy = 4;
1869   const Int_t kipz = 5;
1870   const Int_t kipp = 6;
1871
1872   const Double_t kec = 2.9979251e-4;
1873
1874   /// 
1875   ///     ------------------------------------------------------------------
1876   /// 
1877   ///       units are kgauss,centimeters,gev/c
1878   /// 
1879   vout[kipp] = vect[kipp];
1880   h4 = field * kec;
1881
1882   hxp[0] = - vect[kipy];
1883   hxp[1] = + vect[kipx];
1884  
1885   hp = vect[kipz];
1886
1887   rho = -h4/vect[kipp];
1888   tet = rho * step;
1889   if (TMath::Abs(tet) > 0.15) {
1890     sint = TMath::Sin(tet);
1891     sintt = (sint/tet);
1892     tsint = (tet-sint)/tet;
1893     cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet;
1894   } else {
1895     tsint = tet*tet/36.;
1896     sintt = (1. - tsint);
1897     sint = tet*sintt;
1898     cos1t = 0.5*tet;
1899   }
1900
1901   f1 = step * sintt;
1902   f2 = step * cos1t;
1903   f3 = step * tsint * hp;
1904   f4 = -tet*cos1t;
1905   f5 = sint;
1906   f6 = tet * cos1t * hp;
1907  
1908   vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0];
1909   vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1];
1910   vout[kiz] = vect[kiz] + f1*vect[kipz] + f3;
1911  
1912   vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0];
1913   vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1];
1914   vout[kipz] = vect[kipz] + f4*vect[kipz] + f6;
1915
1916   return;
1917 }
1918
1919 ///__________________________________________________________________________
1920
1921 ///______________________________________________________________________________
1922
1923 void AliHLTMUONFullTracker::OneStepRungekutta(Double_t charge, Double_t step,
1924                                               Double_t* vect, Double_t* vout)
1925 {
1926   ////  ******************************************************************
1927   ////  *                                                                *
1928   ////  *  Runge-Kutta method for tracking a particle through a magnetic *
1929   ////  *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
1930   ////  *  Standards, procedure 25.5.20)                                 *
1931   ////  *                                                                *
1932   ////  *  Input parameters                                              *
1933   ////  *       CHARGE    Particle charge                                *
1934   ////  *       STEP      Step size                                      *
1935   ////  *       VECT      Initial co-ords,direction cosines,momentum     *
1936   ////  *  Output parameters                                             *
1937   ////  *       VOUT      Output co-ords,direction cosines,momentum      *
1938   ////  *  User routine called                                           *
1939   ////  *       CALL GUFLD(X,F)                                          *
1940   ////  *                                                                *
1941   ////  *    ==>Called by : <USER>, GUSWIM                               *
1942   ////  *       Authors    R.Brun, M.Hansroul  *********                 *
1943   ////  *                  V.Perevoztchikov (CUT STEP implementation)    *
1944   ////  *                                                                *
1945   ////  *                                                                *
1946   ////  ******************************************************************
1947
1948   Double_t h2, h4, f[4];
1949   Double_t xyzt[3], a, b, c, ph,ph2;
1950   Double_t secxs[4],secys[4],seczs[4],hxp[3];
1951   Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
1952   Double_t est, at, bt, ct, cba;
1953   Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
1954
1955   Double_t x;
1956   Double_t y;
1957   Double_t z;
1958
1959   Double_t xt;
1960   Double_t yt;
1961   Double_t zt;
1962
1963   Double_t maxit = 1992;
1964   Double_t maxcut = 11;
1965
1966   const Double_t kdlt   = 1e-4;
1967   const Double_t kdlt32 = kdlt/32.;
1968   const Double_t kthird = 1./3.;
1969   const Double_t khalf  = 0.5;
1970   const Double_t kec = 2.9979251e-4;
1971
1972   const Double_t kpisqua = 9.86960440109;
1973   const Int_t kix  = 0;
1974   const Int_t kiy  = 1;
1975   const Int_t kiz  = 2;
1976   const Int_t kipx = 3;
1977   const Int_t kipy = 4;
1978   const Int_t kipz = 5;
1979
1980   /// *.
1981   /// *.    ------------------------------------------------------------------
1982   /// *.
1983   /// *             this constant is for units cm,gev/c and kgauss
1984   /// *
1985   Int_t iter = 0;
1986   Int_t ncut = 0;
1987   for(Int_t j = 0; j < 7; j++)
1988     vout[j] = vect[j];
1989
1990   Double_t  pinv   = kec * charge / vect[6];
1991   Double_t tl = 0.;
1992   Double_t h = step;
1993   Double_t rest;
1994
1995
1996   do {
1997     rest  = step - tl;
1998     if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
1999     ///cmodif: call gufld(vout,f) changed into:
2000     TGeoGlobalMagField::Instance()->Field(vout,f);
2001
2002     /// *
2003     /// *             start of integration
2004     /// *
2005     x      = vout[0];
2006     y      = vout[1];
2007     z      = vout[2];
2008     a      = vout[3];
2009     b      = vout[4];
2010     c      = vout[5];
2011
2012     h2     = khalf * h;
2013     h4     = khalf * h2;
2014     ph     = pinv * h;
2015     ph2    = khalf * ph;
2016     secxs[0] = (b * f[2] - c * f[1]) * ph2;
2017     secys[0] = (c * f[0] - a * f[2]) * ph2;
2018     seczs[0] = (a * f[1] - b * f[0]) * ph2;
2019     ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
2020     if (ang2 > kpisqua) break;
2021
2022     dxt    = h2 * a + h4 * secxs[0];
2023     dyt    = h2 * b + h4 * secys[0];
2024     dzt    = h2 * c + h4 * seczs[0];
2025     xt     = x + dxt;
2026     yt     = y + dyt;
2027     zt     = z + dzt;
2028     /// *
2029     /// *              second intermediate point
2030     /// *
2031
2032     est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
2033     if (est > h) {
2034       if (ncut++ > maxcut) break;
2035       h *= khalf;
2036       continue;
2037     }
2038
2039     xyzt[0] = xt;
2040     xyzt[1] = yt;
2041     xyzt[2] = zt;
2042
2043     ///cmodif: call gufld(xyzt,f) changed into:
2044     TGeoGlobalMagField::Instance()->Field(xyzt,f);
2045
2046     at     = a + secxs[0];
2047     bt     = b + secys[0];
2048     ct     = c + seczs[0];
2049
2050     secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
2051     secys[1] = (ct * f[0] - at * f[2]) * ph2;
2052     seczs[1] = (at * f[1] - bt * f[0]) * ph2;
2053     at     = a + secxs[1];
2054     bt     = b + secys[1];
2055     ct     = c + seczs[1];
2056     secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
2057     secys[2] = (ct * f[0] - at * f[2]) * ph2;
2058     seczs[2] = (at * f[1] - bt * f[0]) * ph2;
2059     dxt    = h * (a + secxs[2]);
2060     dyt    = h * (b + secys[2]);
2061     dzt    = h * (c + seczs[2]);
2062     xt     = x + dxt;
2063     yt     = y + dyt;
2064     zt     = z + dzt;
2065     at     = a + 2.*secxs[2];
2066     bt     = b + 2.*secys[2];
2067     ct     = c + 2.*seczs[2];
2068
2069     est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
2070     if (est > 2.*TMath::Abs(h)) {
2071       if (ncut++ > maxcut) break;
2072       h *= khalf;
2073       continue;
2074     }
2075
2076     xyzt[0] = xt;
2077     xyzt[1] = yt;
2078     xyzt[2] = zt;
2079
2080     ///cmodif: call gufld(xyzt,f) changed into:
2081     TGeoGlobalMagField::Instance()->Field(xyzt,f);
2082
2083     z      = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
2084     y      = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
2085     x      = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
2086
2087     secxs[3] = (bt*f[2] - ct*f[1])* ph2;
2088     secys[3] = (ct*f[0] - at*f[2])* ph2;
2089     seczs[3] = (at*f[1] - bt*f[0])* ph2;
2090     a      = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
2091     b      = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
2092     c      = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
2093
2094     est    = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
2095       + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
2096       + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
2097
2098     if (est > kdlt && TMath::Abs(h) > 1.e-4) {
2099       if (ncut++ > maxcut) break;
2100       h *= khalf;
2101       continue;
2102     }
2103
2104     ncut = 0;
2105     /// *               if too many iterations, go to helix
2106     if (iter++ > maxit) break;
2107
2108     tl += h;
2109     if (est < kdlt32)
2110       h *= 2.;
2111     cba    = 1./ TMath::Sqrt(a*a + b*b + c*c);
2112     vout[0] = x;
2113     vout[1] = y;
2114     vout[2] = z;
2115     vout[3] = cba*a;
2116     vout[4] = cba*b;
2117     vout[5] = cba*c;
2118     rest = step - tl;
2119     if (step < 0.) rest = -rest;
2120     if (rest < 1.e-5*TMath::Abs(step)) return;
2121
2122   } while(1);
2123
2124   /// angle too big, use helix
2125
2126   f1  = f[0];
2127   f2  = f[1];
2128   f3  = f[2];
2129   f4  = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
2130   rho = -f4*pinv;
2131   tet = rho * step;
2132
2133   hnorm = 1./f4;
2134   f1 = f1*hnorm;
2135   f2 = f2*hnorm;
2136   f3 = f3*hnorm;
2137
2138   hxp[0] = f2*vect[kipz] - f3*vect[kipy];
2139   hxp[1] = f3*vect[kipx] - f1*vect[kipz];
2140   hxp[2] = f1*vect[kipy] - f2*vect[kipx];
2141
2142   hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
2143
2144   rho1 = 1./rho;
2145   sint = TMath::Sin(tet);
2146   cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
2147
2148   g1 = sint*rho1;
2149   g2 = cost*rho1;
2150   g3 = (tet-sint) * hp*rho1;
2151   g4 = -cost;
2152   g5 = sint;
2153   g6 = cost * hp;
2154
2155   vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
2156   vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
2157   vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
2158
2159   vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
2160   vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
2161   vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
2162
2163   return;
2164 }
2165
2166 ///______________________________________________________________________________
2167
2168
2169 Bool_t AliHLTMUONFullTracker::SelectFront()
2170 {
2171   Cluster clus1,clus2;
2172   Int_t minIndex=0,maxIndex=0;
2173   Int_t minCh=0,maxCh=0;
2174   Int_t ifronttrackseg=0;
2175   
2176   Float_t xf,yf,zf,thetaDev;
2177   Float_t p,spx,spy,spz,px,py,pz,sx,sy,sz,x,y,z;
2178   Double_t charge;
2179   Float_t dist,mindist;
2180   Int_t frontsegIndex ;
2181
2182   for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
2183
2184     if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
2185     
2186     ///     if(fBackTrackSeg[ibacktrackseg].fIndex[2]==-1 || fBackTrackSeg[ibacktrackseg].fIndex[3]==-1) continue;
2187
2188     ifronttrackseg = fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]-1];
2189     
2190 #ifdef PRINT_SELECT
2191     HLTInfo("AliHLTMUONFullTracker::SelectFront()--Begin\n\n");
2192     HLTInfo("\tbacktrack : %d is connected with : %d front tracks\n",
2193            ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg]);
2194 #endif
2195
2196     maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
2197     maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
2198
2199     minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
2200     minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
2201
2202
2203     
2204     clus2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
2205     clus2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
2206     clus2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
2207     
2208     clus1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
2209     clus1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
2210     clus1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
2211
2212     zf = 0.5*(AliMUONConstants::DefaultChamberZ(4) + AliMUONConstants::DefaultChamberZ(5)); 
2213     thetaDev= (1/zf)*(clus1.fY*clus2.fZ - clus2.fY*clus1.fZ)/(clus2.fZ - clus1.fZ);
2214     xf = clus1.fX*zf/clus1.fZ;
2215     yf = clus2.fY - ((clus2.fY - clus1.fY)*(clus2.fZ-zf))/(clus2.fZ - clus1.fZ);
2216     p = 3.0*0.3/thetaDev;
2217     charge = (p>0)?-1:+1;
2218
2219     fCharge[ibacktrackseg] = charge;
2220
2221     /// cout<<"Charge : "<<charge<<endl;
2222     ///fTrackParam[ibacktrackseg].SetCharge(charge);
2223     spx = p*xf/zf;
2224     spy = p*yf/zf;
2225     spz = sqrt((p*p-(spx*spx + spy*spy)));
2226
2227     if(zf<0)spz = -spz;
2228     
2229     sx = clus1.fX; sy = clus1.fY; sz = clus1.fZ;
2230     mindist = 200000.0;
2231     frontsegIndex = -1;
2232     for( Int_t iconnected=0;iconnected<fNofConnectedfrontTrackSeg[ibacktrackseg];iconnected++){
2233       
2234       ifronttrackseg = fBackToFront[ibacktrackseg][iconnected];
2235       
2236       minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2237       minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2238       clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2239       clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2240       clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2241       
2242       x = sx; y = sy; z = sz;
2243       px = spx; py = spy; pz = spz;
2244       PropagateTracks(charge,px,py,pz,x,y,z,clus1.fZ);
2245       
2246       dist = sqrt((clus1.fX-x)*(clus1.fX-x) + 
2247                   (clus1.fY-y)*(clus1.fY-y));
2248       if(dist<mindist){
2249         mindist = dist;
2250         frontsegIndex = ifronttrackseg;
2251       }
2252     }///for loop on all connected front track segs
2253     
2254     fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2255     ///have to check later
2256     if(frontsegIndex==-1) continue;
2257
2258     fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]++] = frontsegIndex; 
2259     
2260     ///     fTrackParam[ibacktrackseg] = trackParam;
2261
2262 #ifdef PRINT_SELECT
2263     HLTInfo("\tbacktrack : %d is connected with : %d front tracks\n",
2264            ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg]);
2265     HLTInfo("AliHLTMUONFullTracker::SelectFront()--End\n\n");
2266
2267 #endif
2268     
2269   }///backtrackSeg loop
2270   
2271   
2272   return true;
2273 }
2274
2275 ///__________________________________________________________________________
2276
2277 Bool_t AliHLTMUONFullTracker::KalmanChi2Test()
2278 {
2279   
2280   /// Kalman Chi2 test for trach segments selection
2281
2282   Cluster clus1,clus2;
2283   Int_t minIndex=0,maxIndex=0;
2284   Int_t minCh=0,maxCh=0;
2285   Int_t ifronttrackseg=0;
2286   for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
2287
2288     if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
2289     
2290     ///     if(fBackTrackSeg[ibacktrackseg].fIndex[2]==-1 || fBackTrackSeg[ibacktrackseg].fIndex[3]==-1) continue;
2291
2292     ifronttrackseg = fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]-1];
2293     
2294 #ifdef PRINT_KALMAN
2295     HLTInfo("AliHLTMUONFullTracker::KalmanChi2Test()--Begin\n\n");
2296     HLTInfo("\tbacktrack : %d is connected with : %d front tracks, front track index %d\n",
2297            ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg],ifronttrackseg);
2298 #endif
2299     maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
2300     maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
2301
2302     minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
2303     minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
2304       
2305     
2306     clus2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
2307     clus2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
2308     clus2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
2309     clus2.fErrX2 =  0.020736;
2310     clus2.fErrY2 =  0.000100;
2311     
2312     clus1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
2313     clus1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
2314     clus1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
2315     clus1.fErrX2 =  0.020736;
2316     clus1.fErrY2 =  0.000100;
2317
2318
2319     AliMUONTrackParam trackParam;
2320     Double_t dZ =  clus2.fZ - clus1.fZ;
2321     trackParam.SetNonBendingCoor(clus2.fX);
2322     trackParam.SetBendingCoor(clus2.fY);
2323     trackParam.SetZ(clus2.fZ);
2324     trackParam.SetNonBendingSlope((clus2.fX - clus1.fX) / dZ);
2325     trackParam.SetBendingSlope((clus2.fY - clus1.fY) / dZ);
2326     Double_t bendingImpact = clus1.fY - clus1.fZ * trackParam.GetBendingSlope();
2327     Double_t inverseBendingMomentum = 1. 
2328       / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
2329     trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
2330     
2331 #ifdef PRINT_KALMAN
2332     HLTInfo("\t\tCh%d  : (%f,%f,%f)\n",maxCh,clus2.fX,clus2.fY,clus2.fZ);
2333     HLTInfo("\t\tCh%d  : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2334     HLTInfo("\t\tFor minCh : %d, GetBeMom : %lf\n",minCh,trackParam.GetInverseBendingMomentum());
2335 #endif      
2336
2337     ///       trackParam->SetClusterPtr(clus[8]);
2338       
2339     /// => Reset track parameter covariances at last clus (as if the other cluss did not exist)
2340     TMatrixD lastParamCov(5,5);
2341     lastParamCov.Zero();
2342     lastParamCov(0,0) = clus1.fErrX2;
2343     lastParamCov(1,1) = 100. * ( clus2.fErrX2 + clus1.fErrX2 ) / dZ / dZ;
2344     lastParamCov(2,2) = clus1.fErrY2;
2345     lastParamCov(3,3) = 100. * ( clus2.fErrY2 + clus1.fErrY2 ) / dZ / dZ;
2346     lastParamCov(4,4) = ((10.0*10.0 + (clus2.fZ * clus2.fZ * clus1.fErrY2 + 
2347                                        clus1.fZ * clus1.fZ * clus2.fErrY2) 
2348                           / dZ / dZ) /bendingImpact / bendingImpact +
2349                          0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum;
2350     trackParam.SetCovariances(lastParamCov);
2351
2352     AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2353     ///     AliMUONTrackExtrap::ExtrapToZCov(trackParam, AliMUONConstants::DefaultChamberZ(minCh),kTRUE);
2354     ///AliMUONTrackExtrap::ExtrapToZCov(&trackParam, clus1.fZ,kTRUE);
2355
2356     LinearExtrapToZ(&trackParam, clus1.fZ);
2357
2358     trackParam.SetTrackChi2(0.);
2359     
2360     Double_t chi2 = 0.0;
2361
2362     chi2 = KalmanFilter(trackParam,&clus1);
2363
2364 #ifdef PRINT_KALMAN
2365     HLTInfo("\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\n",minCh,chi2,trackParam.GetInverseBendingMomentum());
2366     ///       TMatrixD para2(trackParam->GetParameters());
2367     ///       para2.Print();
2368 #endif      
2369
2370     AliMUONTrackParam extrapTrackParamAtCluster2,minChi2Param;
2371     Double_t chi2OfCluster;
2372     Bool_t tryonefast;
2373     Double_t minChi2=1000000.0;
2374     Int_t frontsegIndex = -1;
2375     extrapTrackParamAtCluster2 = trackParam ;
2376     minChi2Param = trackParam ;
2377
2378     for( Int_t iconnected=0;iconnected<fNofConnectedfrontTrackSeg[ibacktrackseg];iconnected++){
2379       
2380       ifronttrackseg = fBackToFront[ibacktrackseg][iconnected];
2381       AliMUONTrackParam extrapTrackParam;
2382       trackParam = extrapTrackParamAtCluster2;
2383       
2384       minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2385       minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2386       clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2387       clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2388       clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2389       clus1.fErrX2 =  0.020736;
2390       clus1.fErrY2 =  0.000100;
2391
2392     
2393       AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2394       trackParam.ResetPropagator();
2395       AliMUONTrackExtrap::ExtrapToZCov(&trackParam,clus1.fZ, kTRUE);
2396       
2397       tryonefast = TryOneClusterFast(trackParam, &clus1);
2398       ///if (!tryonefast) continue;
2399       
2400       /// try to add the current cluster accuratly
2401       chi2OfCluster = TryOneCluster(trackParam, &clus1 , extrapTrackParam,kTRUE);
2402       
2403       extrapTrackParam.SetExtrapParameters(extrapTrackParam.GetParameters());
2404       extrapTrackParam.SetExtrapCovariances(extrapTrackParam.GetCovariances());
2405       
2406       chi2 = KalmanFilter(extrapTrackParam,&clus1);
2407       if(chi2<minChi2){
2408         minChi2 = chi2;
2409         minChi2Param = extrapTrackParam;
2410         frontsegIndex = ifronttrackseg;
2411       }
2412     }///for loop on all connected front track segs
2413     
2414     fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2415     ///have to check later
2416     if(frontsegIndex==-1) continue;
2417
2418     fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]++] = frontsegIndex; 
2419     ifronttrackseg = frontsegIndex;
2420     
2421     trackParam = minChi2Param ;
2422     
2423 #ifdef PRINT_KALMAN
2424     HLTInfo("\t\t\tCh%d  : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2425     HLTInfo("\t\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\t",minCh,chi2,trackParam.GetInverseBendingMomentum());
2426     HLTInfo(Form("Pt : %lf\n",TMath::Sqrt(trackParam.Px()*trackParam.Px() + 
2427                                          trackParam.Py()*trackParam.Py())));
2428 #endif      
2429     
2430     
2431     minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2432     minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2433     clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2434     clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2435     clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2436     clus1.fErrX2 =  0.020736;
2437     clus1.fErrY2 =  0.000100;
2438     
2439     AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2440     trackParam.ResetPropagator();
2441     ///AliMUONTrackExtrap::ExtrapToZCov(&trackParam,clus1.fZ, kTRUE);
2442     LinearExtrapToZ(&trackParam, clus1.fZ);
2443
2444     tryonefast = TryOneClusterFast(trackParam, &clus1);
2445     ///if (!tryonefast) continue;
2446           
2447     chi2OfCluster = TryOneCluster(trackParam, &clus1 , extrapTrackParamAtCluster2,kTRUE);
2448           
2449     extrapTrackParamAtCluster2.SetExtrapParameters(extrapTrackParamAtCluster2.GetParameters());
2450     extrapTrackParamAtCluster2.SetExtrapCovariances(extrapTrackParamAtCluster2.GetCovariances());
2451     
2452     chi2 = KalmanFilter(extrapTrackParamAtCluster2,&clus1);
2453     
2454     trackParam = extrapTrackParamAtCluster2;
2455     
2456 #ifdef PRINT_KALMAN
2457     HLTInfo("\t\tCh%d  : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2458     HLTInfo("\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\t",minCh,chi2,trackParam.GetInverseBendingMomentum());
2459     HLTInfo(Form("Pt : %lf\n",TMath::Sqrt(extrapTrackParamAtCluster2.Px()*extrapTrackParamAtCluster2.Px() + 
2460                                          extrapTrackParamAtCluster2.Py()*extrapTrackParamAtCluster2.Py())));
2461     trackParam.Print();
2462     HLTInfo("AliHLTMUONFullTracker::KalmanChi2Test()--End\n\n");
2463 #endif      
2464     ///AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.);
2465     fTrackParam[ibacktrackseg] = trackParam;
2466     
2467
2468     
2469   }///trig loop
2470   
2471
2472   
2473   return true;
2474 }
2475
2476 ///__________________________________________________________________________
2477
2478
2479 Double_t AliHLTMUONFullTracker::BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ)
2480 {
2481   //// Returns the mean total momentum energy loss of muon with total momentum='pTotal'
2482   //// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
2483   Double_t muMass = 0.105658369; /// GeV
2484   Double_t eMass = 0.510998918e-3; /// GeV
2485   Double_t k = 0.307075e-3; /// GeV.g^-1.cm^2
2486   Double_t i = 9.5e-9; /// mean exitation energy per atomic Z (GeV)
2487   Double_t p2=pTotal*pTotal;
2488   Double_t beta2=p2/(p2 + muMass*muMass);
2489   
2490   Double_t w = k * rho * pathLength * atomicZ / atomicA / beta2;
2491   
2492   if (beta2/(1-beta2)>3.5*3.5)
2493     return w * (log(2.*eMass*3.5/(i*atomicZ)) + 0.5*log(beta2/(1-beta2)) - beta2);
2494   
2495   return w * (log(2.*eMass*beta2/(1-beta2)/(i*atomicZ)) - beta2);
2496 }
2497
2498
2499 ///__________________________________________________________________________
2500
2501 Double_t AliHLTMUONFullTracker::EnergyLossFluctuation2(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ)
2502 {
2503   //// Returns the total momentum energy loss fluctuation of muon with total momentum='pTotal'
2504   //// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
2505   Double_t muMass = 0.105658369; /// GeV
2506   ///Double_t eMass = 0.510998918e-3; /// GeV
2507   Double_t k = 0.307075e-3; /// GeV.g^-1.cm^2
2508   Double_t p2=pTotal*pTotal;
2509   Double_t beta2=p2/(p2 + muMass*muMass);
2510   
2511   Double_t fwhm = 2. * k * rho * pathLength * atomicZ / atomicA / beta2; /// FWHM of the energy loss Landau distribution
2512   Double_t sigma2 = fwhm * fwhm / (8.*log(2.)); /// gaussian: fwmh = 2 * srqt(2*ln(2)) * sigma (i.e. fwmh = 2.35 * sigma)
2513   
2514   ///sigma2 = k * rho * pathLength * atomicZ / atomicA * eMass; /// sigma2 of the energy loss gaussian distribution
2515   
2516   return sigma2;
2517 }
2518
2519
2520 ///__________________________________________________________________________
2521
2522 void AliHLTMUONFullTracker::LinearExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd)
2523 {
2524   //// Track parameters (and their covariances if any) linearly extrapolated to the plane at "zEnd".
2525   //// On return, results from the extrapolation are updated in trackParam.
2526   
2527   if (trackParam->GetZ() == zEnd) return; /// nothing to be done if same z
2528
2529   Double_t dZ = zEnd - trackParam->GetZ();
2530   trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + trackParam->GetNonBendingSlope() * dZ);
2531   trackParam->SetBendingCoor(trackParam->GetBendingCoor() + trackParam->GetBendingSlope() * dZ);
2532   trackParam->SetZ(zEnd);
2533 }
2534
2535 void AliHLTMUONFullTracker::CorrectELossEffectInAbsorber(AliMUONTrackParam* param, Double_t eLoss)
2536 {
2537   Double_t nonBendingSlope = param->GetNonBendingSlope();
2538   Double_t bendingSlope = param->GetBendingSlope();
2539   param->SetInverseBendingMomentum(param->GetCharge() / (param->P() + eLoss) *
2540                                    TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) /
2541                                    TMath::Sqrt(1.0 + bendingSlope*bendingSlope));
2542 }
2543
2544 ///__________________________________________________________________________
2545
2546 void AliHLTMUONFullTracker::Cov2CovP(const TMatrixD &param, TMatrixD &cov)
2547 {
2548   //// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, SlopeX, Y, SlopeY, q*PTot)
2549   //// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
2550   
2551   /// charge * total momentum
2552   Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
2553     TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
2554   
2555   /// Jacobian of the opposite transformation
2556   TMatrixD jacob(5,5);
2557   jacob.UnitMatrix();
2558   jacob(4,1) = qPTot * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
2559   jacob(4,3) = - qPTot * param(1,0) * param(1,0) * param(3,0) /
2560     (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
2561   jacob(4,4) = - qPTot / param(4,0);
2562   
2563   /// compute covariances in new coordinate system
2564   TMatrixD tmp(5,5);
2565   tmp.MultT(cov,jacob);
2566   cov.Mult(jacob,tmp);
2567
2568 }
2569
2570 ///__________________________________________________________________________
2571
2572 void AliHLTMUONFullTracker::CovP2Cov(const TMatrixD &param, TMatrixD &covP)
2573 {
2574   //// change coordinate system: (X, SlopeX, Y, SlopeY, q*PTot) -> (X, SlopeX, Y, SlopeY, q/Pyz)
2575   //// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
2576   
2577   /// charge * total momentum
2578   Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
2579     TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
2580   
2581   /// Jacobian of the transformation
2582   TMatrixD jacob(5,5);
2583   jacob.UnitMatrix();
2584   jacob(4,1) = param(4,0) * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
2585   jacob(4,3) = - param(4,0) * param(1,0) * param(1,0) * param(3,0) /
2586     (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
2587   jacob(4,4) = - param(4,0) / qPTot;
2588   
2589   /// compute covariances in new coordinate system
2590   TMatrixD tmp(5,5);
2591   tmp.MultT(covP,jacob);
2592   covP.Mult(jacob,tmp);
2593
2594 }
2595
2596 ///__________________________________________________________________________
2597
2598
2599 void AliHLTMUONFullTracker::CorrectMCSEffectInAbsorber(AliMUONTrackParam* param,
2600                                                               Double_t xVtx, Double_t yVtx, Double_t zVtx,
2601                                                               Double_t absZBeg, Double_t f1, Double_t f2)
2602 {
2603   //// Correct parameters and corresponding covariances using Branson correction
2604   //// - input param are parameters and covariances at the end of absorber
2605   //// - output param are parameters and covariances at vertex
2606   //// Absorber correction parameters are supposed to be calculated at the current track z-position
2607   
2608   /// Position of the Branson plane (spectro. (z<0))
2609   Double_t zB = (f1>0.) ? absZBeg - f2/f1 : 0.;
2610   
2611   LinearExtrapToZ(param,zB);
2612   
2613   /// compute track parameters at vertex
2614   TMatrixD newParam(5,1);
2615   newParam.Zero();
2616   newParam(0,0) = xVtx;
2617   newParam(1,0) = (param->GetNonBendingCoor() - xVtx) / (zB - zVtx);
2618   newParam(2,0) = yVtx;
2619   newParam(3,0) = (param->GetBendingCoor() - yVtx) / (zB - zVtx);
2620   newParam(4,0) = param->GetCharge() / param->P() *
2621     TMath::Sqrt(1.0 + newParam(1,0)*newParam(1,0) + newParam(3,0)*newParam(3,0)) /
2622     TMath::Sqrt(1.0 + newParam(3,0)*newParam(3,0));
2623   
2624   TMatrixD paramCovP(5,5);
2625   
2626   /// Get covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
2627   paramCovP.Use(param->GetCovariances());
2628   
2629   Cov2CovP(param->GetParameters(),paramCovP);
2630   
2631   /// Get the covariance matrix in the (XVtx, X, YVtx, Y, q*PTot) coordinate system
2632   ///   TMatrixD paramCovVtx(5,5);
2633   Double_t element44 = paramCovP(4,4);
2634   paramCovP.Zero();
2635   paramCovP(4,4) = element44;
2636   
2637   CovP2Cov(newParam,paramCovP);
2638   
2639   /// Set parameters and covariances at vertex
2640   param->SetParameters(newParam);
2641   param->SetZ(zVtx);
2642   param->SetCovariances(paramCovP);
2643 }
2644
2645 ///__________________________________________________________________________
2646
2647 Bool_t AliHLTMUONFullTracker::ExtrapolateToOrigin(Bool_t extrap)
2648 {
2649   /// Extrapolation to origin through absorber
2650
2651   Int_t minIndex=0,maxIndex=0;
2652   Int_t minCh=0,maxCh=0;
2653   Int_t ifronttrackseg = -1;
2654   AliMUONTrackParam trackP;
2655   Double_t BSlope, NBSlope;
2656   AliHLTMUONRecHitStruct p1,p2,pSeg1,pSeg2;
2657   Double_t Pyz = -1.0;
2658   TVector3 v1,v2,v3,v4;
2659   Double_t eLoss1,eLoss2,eLoss3;
2660   Double_t b;
2661   Double_t zE,zB,dzE,dzB; 
2662   Double_t F0,F1,F2;
2663   Double_t F0Sum,F1Sum,F2Sum;
2664   Double_t fXVertex=0.0,fYVertex=0.0,fZVertex=0.0;
2665   
2666   for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
2667
2668     if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
2669     
2670     maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
2671     maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
2672     
2673     minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
2674     minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
2675     
2676     p2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ; 
2677     p2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
2678     p2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
2679
2680     p1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ; 
2681     p1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
2682     p1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
2683     
2684     
2685     Sub(&p2,&p1,&pSeg1);
2686     
2687     ifronttrackseg = fBackToFront[ibacktrackseg][0];
2688
2689     maxIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2690     maxCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2691     
2692     minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2693     minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2694
2695     p2.fX = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fX ;
2696     p2.fY = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fY ;
2697     p2.fZ = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fZ ;
2698
2699     p1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2700     p1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2701     p1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2702
2703     Sub(&p2,&p1,&pSeg2);
2704     
2705     Pyz = -(3.0*0.3/sin(Angle(&pSeg1,&pSeg2)));/// *  sqrt(x3*x3 + y3*y3)/z3 ;
2706     
2707     
2708     NBSlope = (p2.fX - p1.fX)/(p2.fZ - p1.fZ); 
2709     BSlope = (p2.fY - p1.fY)/(p2.fZ - p1.fZ); 
2710     
2711     
2712     trackP.SetZ(p1.fZ);
2713     trackP.SetNonBendingCoor(p1.fX);
2714     trackP.SetNonBendingSlope(NBSlope);
2715     trackP.SetBendingCoor(p1.fY);
2716     trackP.SetBendingSlope(BSlope);
2717     trackP.SetInverseBendingMomentum(1.0/Pyz) ;
2718     
2719     
2720     if(extrap){
2721       
2722       trackP = fTrackParam[ibacktrackseg]    ;
2723     
2724       LinearExtrapToZ(&trackP,fgkAbsoedge[3]);
2725       v4.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
2726       LinearExtrapToZ(&trackP,fgkAbsoedge[2]);
2727       v3.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
2728       LinearExtrapToZ(&trackP,fgkAbsoedge[1]);
2729       v2.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
2730       LinearExtrapToZ(&trackP,fgkAbsoedge[0]);
2731       v1.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
2732     
2733       eLoss1 = BetheBloch(trackP.P(), (v4-v3).Mag(), fgkRho[2], fgkAtomicA[2], fgkAtomicZ[2]);
2734       eLoss2 = BetheBloch(trackP.P(), (v3-v2).Mag(), fgkRho[1], fgkAtomicA[1], fgkAtomicZ[1]);
2735       eLoss3 = BetheBloch(trackP.P(), (v2-v1).Mag(), fgkRho[0], fgkAtomicA[0], fgkAtomicZ[0]);
2736     
2737       ///       sigmaELoss1 = EnergyLossFluctuation2(trackP.P(), (v4-v3).Mag(), rho[2], atomicA[2], atomicZ[2]);
2738       ///       sigmaELoss2 = EnergyLossFluctuation2(trackP.P(), (v3-v2).Mag(), rho[1], atomicA[1], atomicZ[1]);
2739       ///       sigmaELoss3 = EnergyLossFluctuation2(trackP.P(), (v2-v1).Mag(), rho[0], atomicA[0], atomicZ[0]);
2740
2741       ///     eDiff = totELoss-(eLoss1+eLoss2+eLoss3);
2742       ///     sigmaELossDiff = sigmaTotELoss ;///- (sigmaELoss1+sigmaELoss2+sigmaELoss3);
2743
2744       ///       CorrectELossEffectInAbsorber(&trackP, 0.5*(eLoss1+eLoss2+eLoss3), 0.5*(sigmaELoss1+sigmaELoss2+sigmaELoss3));
2745     
2746
2747       ///CorrectELossEffectInAbsorber(&trackP, totELoss,sigmaTotELoss);
2748
2749       CorrectELossEffectInAbsorber(&trackP, 0.7*(eLoss1+eLoss2+eLoss3));
2750
2751       ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2752       F0Sum = 0.0;      F1Sum = 0.0;      F2Sum = 0.0;
2753
2754       b = (v4.Z()-v1.Z())/((v4-v1).Mag());
2755     
2756       zB = v1.Z();
2757       zE = b*((v2-v1).Mag()) + zB;
2758       dzB = zB - v1.Z();
2759       dzE = zE - v1.Z();
2760     
2761       F0 = ((v2-v1).Mag())/fgkRadLen[0];
2762       F1  = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[0] /2.;
2763       F2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[0] / 3.;
2764
2765       F0Sum += F0;
2766       F1Sum += F1;
2767       F2Sum += F2;
2768     
2769       zB = zE;
2770       zE = b*((v3-v2).Mag()) + zB;
2771       dzB = zB - v1.Z();
2772       dzE = zE - v1.Z();
2773     
2774       F0 = ((v3-v2).Mag())/fgkRadLen[1];
2775       F1  = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[1] /2.;
2776       F2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[1] / 3.;
2777
2778       F0Sum += F0;
2779       F1Sum += F1;
2780       F2Sum += F2;
2781     
2782       zB = zE;
2783       zE = b*((v4-v3).Mag()) + zB;
2784       dzB = zB - v1.Z();
2785       dzE = zE - v1.Z();
2786     
2787       F0 = ((v4-v3).Mag())/fgkRadLen[2];
2788       F1  = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[2] /2.;
2789       F2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[2] / 3.;
2790
2791       F0Sum += F0;
2792       F1Sum += F1;
2793       F2Sum += F2;
2794     
2795     
2796       ///AddMCSEffectInAbsorber(&trackP,(v4-v1).Mag(),F0Sum,F1Sum,F2Sum);
2797
2798       CorrectMCSEffectInAbsorber(&trackP,fXVertex,fYVertex, fZVertex,AliMUONConstants::AbsZBeg(),F1Sum,F2Sum);
2799       CorrectELossEffectInAbsorber(&trackP, 0.5*(eLoss1+eLoss2+eLoss3));
2800     }
2801
2802     ///AliMUONTrackExtrap::ExtrapToVertex(&trackP, 0., 0., 0., 0., 0.);
2803     fTrackParam[ibacktrackseg] = trackP;
2804     
2805   }///backtrackseg for loop
2806   
2807   return true;
2808 }
2809
2810 ///__________________________________________________________________________
2811
2812 Bool_t AliHLTMUONFullTracker::InitGRP() 
2813 {
2814   ///GRP handling needed for standalone testing and debug purpose
2815
2816   AliGRPObject*  fGRPData = 0x0;              /// Data from the GRP/GRP/Data CDB folder
2817
2818   ///------------------------------------
2819   /// Initialization of the GRP entry 
2820   ///------------------------------------
2821   AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
2822
2823   if (entry) {
2824
2825     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  /// old GRP entry
2826
2827     if (m) {
2828       HLTInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject\n");
2829       m->Print();
2830       fGRPData = new AliGRPObject();
2831       fGRPData->ReadValuesFromMap(m);
2832     }
2833
2834     else {
2835       HLTInfo("Found an AliGRPObject in GRP/GRP/Data, reading it\n");
2836       fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject());  /// new GRP entry
2837       entry->SetOwner(0);
2838     }
2839
2840     ///    FIX ME: The unloading of GRP entry is temporarily disabled
2841     ///    because ZDC and VZERO are using it in order to initialize
2842     ///    their reconstructor objects. In the future one has to think
2843     ///    of propagating AliRunInfo to the reconstructors.
2844     ///    AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
2845   }
2846
2847   if (!fGRPData) {
2848     HLTInfo("No GRP entry found in OCDB!\n");
2849     return kFALSE;
2850   }
2851
2852   TString lhcState = fGRPData->GetLHCState();
2853   if (lhcState==AliGRPObject::GetInvalidString()) {
2854     HLTInfo("GRP/GRP/Data entry:  missing value for the LHC state ! Using UNKNOWN\n");
2855     lhcState = "UNKNOWN";
2856   }
2857
2858   TString beamType = fGRPData->GetBeamType();
2859   if (beamType==AliGRPObject::GetInvalidString()) {
2860     HLTInfo("GRP/GRP/Data entry:  missing value for the beam type ! Using UNKNOWN\n");
2861     beamType = "UNKNOWN";
2862   }
2863
2864   Float_t beamEnergy = fGRPData->GetBeamEnergy();
2865   if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
2866     HLTInfo("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0\n");
2867     beamEnergy = 0;
2868   }
2869   /// LHC: "multiply by 120 to get the energy in MeV"
2870   beamEnergy *= 0.120;
2871
2872   TString runType = fGRPData->GetRunType();
2873   if (runType==AliGRPObject::GetInvalidString()) {
2874     HLTInfo("GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN\n");
2875     runType = "UNKNOWN";
2876   }
2877
2878   Int_t activeDetectors = fGRPData->GetDetectorMask();
2879   if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
2880     HLTInfo("GRP/GRP/Data entry:  missing value for the detector mask ! Using 1074790399\n");
2881     activeDetectors = 1074790399;
2882   }
2883
2884   AliRunInfo *fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
2885   fRunInfo->Dump();
2886
2887
2888   ///*** Dealing with the magnetic field map
2889   if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
2890     if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
2891       HLTInfo("ExpertMode!!! GRP information will be ignored !\n");
2892       HLTInfo("ExpertMode!!! Running with the externally locked B field !\n");
2893     }
2894     else {
2895       HLTInfo("Destroying existing B field instance!\n");
2896       delete TGeoGlobalMagField::Instance();
2897     }    
2898   }
2899   if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
2900     /// Construct the field map out of the information retrieved from GRP.
2901     Bool_t ok = kTRUE;
2902     /// L3
2903     Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
2904     if (l3Current == AliGRPObject::GetInvalidFloat()) {
2905       HLTInfo("GRP/GRP/Data entry:  missing value for the L3 current !\n");
2906       ok = kFALSE;
2907     }
2908     
2909     Char_t l3Polarity = fGRPData->GetL3Polarity();
2910     if (l3Polarity == AliGRPObject::GetInvalidChar()) {
2911       HLTInfo("GRP/GRP/Data entry:  missing value for the L3 polarity !\n");
2912       ok = kFALSE;
2913     }
2914
2915     /// Dipole
2916     Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
2917     if (diCurrent == AliGRPObject::GetInvalidFloat()) {
2918       HLTInfo("GRP/GRP/Data entry:  missing value for the dipole current !\n");
2919       ok = kFALSE;
2920     }
2921
2922     Char_t diPolarity = fGRPData->GetDipolePolarity();
2923     if (diPolarity == AliGRPObject::GetInvalidChar()) {
2924       HLTInfo("GRP/GRP/Data entry:  missing value for the dipole polarity !\n");
2925       ok = kFALSE;
2926     }
2927
2928     /// read special bits for the polarity convention and map type
2929     Int_t  polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
2930     Bool_t uniformB = fGRPData->IsUniformBMap();
2931
2932     if (ok) { 
2933       AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1), 
2934                                              TMath::Abs(diCurrent) * (diPolarity ? -1:1), 
2935                                              polConvention,uniformB,beamEnergy, beamType.Data());
2936       if (fld) {
2937         TGeoGlobalMagField::Instance()->SetField( fld );
2938         TGeoGlobalMagField::Instance()->Lock();
2939         HLTInfo("Running with the B field constructed out of GRP !\n");
2940       }
2941       else HLTInfo("Failed to create a B field map !\n");
2942     }
2943     else HLTInfo("B field is neither set nor constructed from GRP ! Exitig...\n");
2944   }
2945   
2946
2947   return kTRUE;
2948
2949