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