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