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