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