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