1 /**************************************************************************
2 * This file is property of and copyright by the ALICE HLT Project *
3 * All rights reserved. *
6 * Indranil Das <indra.das@saha.ac.in> *
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 **************************************************************************/
17 //////////////////////////////////////////////////////////////////////
18 ///Author : Indranil Das, SINP, INDIA
20 ///Email : indra.das@saha.ac.in || indra.ehep@gmail.com
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 /////////////////////////////////////////////////////////////////////////
35 #include "AliHLTMUONFullTracker.h"
43 #include "TGeoGlobalMagField.h"
46 #include "AliCDBManager.h"
47 #include "AliGeomManager.h"
48 #include "AliGRPObject.h"
50 #include "AliMUONTrackExtrap.h"
51 #include "AliMUONTrackParam.h"
52 #include "AliMUONTrackExtrap.h"
53 #include "AliMUONConstants.h"
54 #include "AliMUONGeometryTransformer.h"
55 #include "AliMUONTrackParam.h"
57 #include "AliMpDDLStore.h"
58 #include "AliMpSegmentation.h"
61 #include "AliCDBEntry.h"
64 #include "AliRunInfo.h"
66 #include "AliHLTMUONConstants.h"
67 #include "AliHLTMUONDataTypes.h"
68 #include "AliHLTMUONUtils.h"
69 #include "AliHLTMUONRecHitsBlockStruct.h"
70 #include "AliHLTMUONTriggerRecordsBlockStruct.h"
71 #include "AliHLTMUONMansoTracksBlockStruct.h"
76 #define PRINT_POINTS 1
79 #define PRINT_KALMAN 1
80 #define PRINT_SELECT 1
81 ///#define PRINT_DETAIL_KALMAN 1
84 class AliHLTMUONFullTracker;
86 const Float_t AliHLTMUONFullTracker::fgkTrackDetCoordinate[3] = {
87 155.179+20.0, 166.234+20.0,
88 (AliMUONConstants::DefaultChamberZ(4)+ AliMUONConstants::DefaultChamberZ(5))/2.0
91 const Double_t AliHLTMUONFullTracker::fgkAbsoedge[4] = {92.0, 317.0,443.0,499.0};
92 const Double_t AliHLTMUONFullTracker::fgkRadLen[3] = {13.875635, 11.273801, 1.765947};
93 const Double_t AliHLTMUONFullTracker::fgkRho[3] = {1.750000, 2.350000, 7.880000};
94 const Double_t AliHLTMUONFullTracker::fgkAtomicZ[3] = {6.000000, 11.098486,25.780000};
95 const Double_t AliHLTMUONFullTracker::fgkAtomicA[3] = {12.010000, 22.334156,55.299670 };
98 const Int_t AliHLTMUONFullTracker::fgkMaxNofCellsPerCh = 400;
99 const Int_t AliHLTMUONFullTracker::fgkMaxNofPointsPerCh = 300;
100 const Int_t AliHLTMUONFullTracker::fgkMaxNofCh = 11 ; /// 10tracking + 1 trigger
101 const Int_t AliHLTMUONFullTracker::fgkMaxNofTracks = 200;
102 const Int_t AliHLTMUONFullTracker::fgkMaxNofConnectedTracks = 20;
106 AliHLTMUONFullTracker::AliHLTMUONFullTracker() :
108 fChamberGeometryTransformer(0x0),
115 fInclinationBack(0x0),
116 fNofConnectedfrontTrackSeg(0x0),
125 fNoffrontTrackSeg(0),
129 /// Constructor of the class
131 /// fChPoint = (AliHLTMUONRecHitStruct***)(malloc(10 * sizeof(AliHLTMUONRecHitStruct)));
132 /// for( Int_t ich=0;ich<10;ich++)
133 /// fChPoint[ich] = (AliHLTMUONRecHitStruct**)(malloc(300 * sizeof(AliHLTMUONRecHitStruct)));
136 fChPoint = new AliHLTMUONRecHitStruct**[fgkMaxNofCh-1];
137 }catch (const std::bad_alloc&){
138 HLTError("Dynamic memory allocation failed in constructor : fChPoint");
142 for( Int_t ich=0;ich<(fgkMaxNofCh-1);ich++)
144 fChPoint[ich] = new AliHLTMUONRecHitStruct*[fgkMaxNofPointsPerCh];
145 }catch (const std::bad_alloc&){
146 HLTError("Dynamic memory allocation failed in constructor : fChPoint[%d]",ich);
151 fChPoint11 = new AliHLTMUONTriggerRecordStruct*[fgkMaxNofPointsPerCh];
152 }catch (const std::bad_alloc&){
153 HLTError("Dynamic memory allocation failed in constructor : fChPoint11");
158 fBackTrackSeg = new TrackSeg[fgkMaxNofTracks];
159 }catch (const std::bad_alloc&){
160 HLTError("Dynamic memory allocation failed in constructor : fBackTrackSeg");
165 fFrontTrackSeg = new TrackSeg[fgkMaxNofTracks];
166 }catch (const std::bad_alloc&){
167 HLTError("Dynamic memory allocation failed in constructor : fFrontTrackSeg");
172 fExtrapSt3X = new float[fgkMaxNofTracks];
173 }catch (const std::bad_alloc&){
174 HLTError("Dynamic memory allocation failed in constructor : fExtrapSt3X");
179 fExtrapSt3Y = new float[fgkMaxNofTracks];
180 }catch (const std::bad_alloc&){
181 HLTError("Dynamic memory allocation failed in constructor : fExtrapSt3Y");
186 fInclinationBack = new float[fgkMaxNofTracks];
187 }catch (const std::bad_alloc&){
188 HLTError("Dynamic memory allocation failed in constructor : fInclinationBack");
193 fNofConnectedfrontTrackSeg = new int[fgkMaxNofTracks];
194 }catch (const std::bad_alloc&){
195 HLTError("Dynamic memory allocation failed in constructor : fNofConnectedfrontTrackSeg");
200 fBackToFront = new int*[fgkMaxNofTracks];
201 }catch (const std::bad_alloc&){
202 HLTError("Dynamic memory allocation failed in constructor : fBackToFront");
206 for( Int_t itr=0;itr<fgkMaxNofTracks;itr++)
208 fBackToFront[itr] = new int[fgkMaxNofConnectedTracks];
209 }catch (const std::bad_alloc&){
210 HLTError("Dynamic memory allocation failed in constructor : fBackToFront[%d]",itr);
216 fCharge = new float[fgkMaxNofTracks];
217 }catch (const std::bad_alloc&){
218 HLTError("Dynamic memory allocation failed in constructor : fCharge");
223 fNofPoints = new int[fgkMaxNofCh];
224 }catch (const std::bad_alloc&){
225 HLTError("Dynamic memory allocation failed in constructor : fNofPoints");
230 fTrackParam = new AliMUONTrackParam[fgkMaxNofTracks];
231 }catch (const std::bad_alloc&){
232 HLTError("Dynamic memory allocation failed in constructor : fTrackParam");
240 ///__________________________________________________________________________
242 AliHLTMUONFullTracker::~AliHLTMUONFullTracker()
244 /// Destructor of the class
246 fChamberGeometryTransformer->Delete();
250 delete []fBackTrackSeg;
251 delete []fFrontTrackSeg;
252 delete []fExtrapSt3X;
253 delete []fExtrapSt3Y;
254 delete []fInclinationBack;
255 delete []fNofConnectedfrontTrackSeg;
256 delete []fBackToFront;
259 delete []fTrackParam;
261 ///delete fChamberGeometryTransformer;
265 ///__________________________________________________________________________
267 Bool_t AliHLTMUONFullTracker::Clear(){
269 /// Clear method to be called after each event
271 for( Int_t ich=0;ich<fgkMaxNofCh;ich++)
274 fNofbackTrackSeg = 0;
275 fNoffrontTrackSeg = 0;
281 ///__________________________________________________________________________
283 void AliHLTMUONFullTracker::Print()
285 /// Print anything mainly for debugging the codes, will be removed later
286 HLTInfo("\nPrint This--->\n");
289 ///__________________________________________________________________________
291 Bool_t AliHLTMUONFullTracker::Init()
293 /// Initilation to be called once, later can be used to set/load the CDB path/entries
295 //FIXME: remove the code for setting the magnetic field. This functionality is
296 // handled by the HLT framework. Doing so here can cause problems.
297 if (AliGeomManager::GetGeometry() == NULL){
298 AliGeomManager::LoadGeometry();
299 AliGeomManager::ApplyAlignObjsFromCDB("GRP MUON");
303 AliMUONTrackExtrap::SetField();
304 /// see header file for class documentation
305 fChamberGeometryTransformer = new AliMUONGeometryTransformer();
306 if(! fChamberGeometryTransformer->LoadGeometryData()){
307 cerr<<__FILE__<<": Failed to Load Geomerty Data "<<endl;
314 ///__________________________________________________________________________
316 Bool_t AliHLTMUONFullTracker::SetInput(AliHLTInt32_t /*ddl*/,const AliHLTMUONRecHitStruct *data, AliHLTInt32_t size)
318 /// Set the input for rechit points
320 AliHLTUInt16_t detElemID;
321 AliHLTUInt8_t chamber;
323 for( Int_t ipoint=0;ipoint<int(size);ipoint++){
324 AliHLTMUONUtils::UnpackRecHitFlags(data->fFlags,chamber,detElemID);
325 fChPoint[detElemID/100-1][fNofPoints[detElemID/100-1]++] = (AliHLTMUONRecHitStruct *)data;
327 HLTInfo("\nch : %02d, detelem : %04d, (X,Y,Z) : (%8.3f,%8.3f,%8.3f)",chamber,detElemID,data->fX,data->fY,data->fZ);
335 ///__________________________________________________________________________
337 Bool_t AliHLTMUONFullTracker::SetInput(AliHLTInt32_t /*ddl*/, const AliHLTMUONTriggerRecordStruct *data, AliHLTInt32_t size)
339 /// Set the input for trigrecs
341 AliHLTUInt16_t detElemID;
342 AliHLTUInt8_t chamber;
344 for( Int_t ipoint=0;ipoint<int(size);ipoint++){
345 fChPoint11[fNofPoints[10]++] = (AliHLTMUONTriggerRecordStruct *)data;
346 for( Int_t ich=0;ich<4;ich++){
347 AliHLTMUONUtils::UnpackRecHitFlags((data->fHit[ich]).fFlags,chamber,detElemID);
349 HLTInfo("\nsize : %d, itrig : %04d, ch : %02d, detelem : %04d, (X,Y,Z) : (%8.3f,%8.3f,%8.3f)",
350 size,ipoint,chamber,detElemID,(data->fHit[ich]).fX,(data->fHit[ich]).fY,(data->fHit[ich]).fZ);
358 ///__________________________________________________________________________
360 Bool_t AliHLTMUONFullTracker::Run( Int_t /*iEvent*/,AliHLTMUONMansoTrackStruct *data, AliHLTUInt32_t& size)
362 /// Main Run call of the class
367 ExtrapolateToOrigin(true);
368 FillOutData(data,size);
369 // HLTDebug("iEvent: %d, has : %d tracks, triggers : %d, nof slat tracks : %d, quad tracks : %d, connected : %d\n",
370 // iEvent,size,fNofPoints[10],fNofbackTrackSeg,fNoffrontTrackSeg,fNofConnected);
375 ///__________________________________________________________________________
377 void AliHLTMUONFullTracker::Sub(const AliHLTMUONRecHitStruct *v1, const AliHLTMUONRecHitStruct *v2, AliHLTMUONRecHitStruct *v3) const
379 /// Subtraction of position co-odinate of two space points
381 v3->fX = v1->fX - v2->fX;
382 v3->fY = v1->fY - v2->fY;
383 v3->fZ = v1->fZ - v2->fZ;
386 ///__________________________________________________________________________
388 Double_t AliHLTMUONFullTracker::Angle(const AliHLTMUONRecHitStruct *v1, const AliHLTMUONRecHitStruct *v2) const
390 ///Angle of a straight line formed using v1 and v2
392 Float_t ptot2 = ((v1->fX * v1->fX) + (v1->fY * v1->fY) + (v1->fZ * v1->fZ))*
393 ((v2->fX * v2->fX) + (v2->fY * v2->fY) + (v2->fZ * v2->fZ));
397 Float_t arg = ((v1->fX * v2->fX) + (v1->fY * v2->fY) + (v1->fZ * v2->fZ))/sqrt(ptot2);
398 if(arg > 1.0) arg = 1.0;
399 if(arg < -1.0) arg = -1.0;
400 return TMath::ACos(arg);
406 ///__________________________________________________________________________
408 Bool_t AliHLTMUONFullTracker::FillOutData(AliHLTMUONMansoTrackStruct *track, AliHLTUInt32_t& size) const
410 ///Fill the output data pointers
412 size = (AliHLTUInt32_t(fNofbackTrackSeg)<size) ? AliHLTUInt32_t(fNofbackTrackSeg) : size;
414 for( Int_t ibackTrackSeg=0;ibackTrackSeg<int(size);ibackTrackSeg++){
416 track->fId = (ibackTrackSeg << 8) | (fBackTrackSeg[ibackTrackSeg].fTrigRec & 0xFF);
417 track->fTrigRec = fBackTrackSeg[ibackTrackSeg].fTrigRec;
418 track->fPx = fTrackParam[ibackTrackSeg].Px();
419 track->fPy = fTrackParam[ibackTrackSeg].Py();
420 track->fPz = fTrackParam[ibackTrackSeg].Pz();
425 for( Int_t ipoint=3;ipoint>=0;ipoint--){
426 track->fHit[ipoint] = AliHLTMUONConstants::NilRecHitStruct();
427 hitset[ipoint] = false;
429 if(fBackTrackSeg[ibackTrackSeg].fIndex[ipoint]!=-1){
431 track->fHit[ipoint] = *(fChPoint[ipoint+6][fBackTrackSeg[ibackTrackSeg].fIndex[ipoint]]);
432 hitset[ipoint] = true;
436 AliHLTMUONParticleSign sign = AliHLTMUONParticleSign(Int_t(TMath::Sign(1.,fTrackParam[ibackTrackSeg].GetInverseBendingMomentum())));
437 track->fFlags = AliHLTMUONUtils::PackMansoTrackFlags(sign,hitset);
446 ///__________________________________________________________________________
448 Bool_t AliHLTMUONFullTracker::SlatTrackSeg()
451 ///Find the Slat Track Segments
453 Float_t trigX1,trigY1,trigZ1=AliMUONConstants::DefaultChamberZ(10);
454 Float_t trigX2,trigY2,trigZ2=AliMUONConstants::DefaultChamberZ(12);
455 Float_t extrapCh9X,extrapCh9Y,extrapCh9Z=AliMUONConstants::DefaultChamberZ(9);
456 Float_t extrapCh8X,extrapCh8Y,extrapCh8Z=AliMUONConstants::DefaultChamberZ(8);
457 Float_t extrapCh7X,extrapCh7Y,extrapCh7Z=AliMUONConstants::DefaultChamberZ(7);
458 Float_t extrapCh6X,extrapCh6Y,extrapCh6Z=AliMUONConstants::DefaultChamberZ(6);
460 Float_t meanX1,meanX2,meanY1,meanY2,meanZ1,meanZ2;
462 Double_t distChFront,distChBack;
463 Int_t nofFrontChPoints,nofBackChPoints;
464 Int_t frontIndex[fgkMaxNofTracks], backIndex[fgkMaxNofTracks];
466 AliHLTMUONRecHitStruct p2,p3,pSeg1,pSeg2,pSeg3;
467 Double_t anglediff,anglediff1,anglediff2;
468 Double_t minAngle = 2.0;
470 Bool_t st5TrackletFound = false;
471 Bool_t ch9PointFound = false;
472 Bool_t ch8PointFound = false;
473 Bool_t st4TrackletFound = false;
474 Bool_t ch7PointFound = false;
475 Bool_t ch6PointFound = false;
477 Int_t index1,index2,index3,index4;
478 IntPair cells[2][fgkMaxNofTracks]; ///cell array for 5 stn for given trigger
481 Float_t maxXDeflectionExtrap = 10.0 + 4.0; ///simulation result 10.0
482 Float_t extrapRatio = 0.2; ///simulation result 0.2
483 Float_t circularWindow = 20.0; ///simulatiuon result 20
484 Float_t minAngleWindow = 2.0 + 1.0; ///simulation result 2.0
488 AliHLTUInt16_t detElemID,prevDetElemID;
489 AliHLTUInt8_t chamber;
491 Int_t minTrgCh,maxTrgCh;
492 for( Int_t itrig=0;itrig<fNofPoints[10];itrig++){
495 st5TrackletFound = false;
496 ch9PointFound = false;
497 ch8PointFound = false;
499 st4TrackletFound = false;
500 ch7PointFound = false;
501 ch6PointFound = false;
508 ///for( Int_t istn=0;istn<5;istn++)
511 for( Int_t ihit=0;ihit<4;ihit++){
513 AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[ihit]).fFlags,chamber,detElemID);
515 if(ihit==0 and detElemID!=0)
517 if(ihit==1 and detElemID!=0 and prevDetElemID==0)
519 if(ihit==2 and detElemID!=0)
521 if(ihit==3 and detElemID!=0 and prevDetElemID==0)
524 prevDetElemID = detElemID;
527 if(minTrgCh == -1 or maxTrgCh == -1)
530 AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[minTrgCh]).fFlags,chamber,detElemID);
531 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,trigZ1);
532 AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[maxTrgCh]).fFlags,chamber,detElemID);
533 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,trigZ2);
536 trigX1 = (fChPoint11[itrig]->fHit[minTrgCh]).fX;
537 trigY1 = (fChPoint11[itrig]->fHit[minTrgCh]).fY;
539 trigX2 = (fChPoint11[itrig]->fHit[maxTrgCh]).fX;
540 trigY2 = (fChPoint11[itrig]->fHit[maxTrgCh]).fY;
545 HLTInfo("AliHLTMUONFullTracker::SlatTrackSeg()--Begin\n\n");
546 AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[minTrgCh]).fFlags,chamber,detElemID);
547 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,tz);
548 HLTInfo("\ttrig 1 : (%f,%f,%f), org : %f\n",trigX1,trigY1,trigZ1,tz);
549 AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[maxTrgCh]).fFlags,chamber,detElemID);
550 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,trigz);
551 HLTInfo("\ttrig 2 : (%f,%f,%f), org : %f\n",trigX2,trigY2,trigZ2,tz);
554 /////////////////////////////////////////////////// Stn 5///////////////////////////////////////////////////////////////
556 ///extrapCh9X = trigX1 + (trigX2-trigX1) * (extrapCh9Z-trigZ1)/(trigZ2 - trigZ1) ;
557 /// extrapCh9Yref = trigY1 * extrapCh9Z/trigZ1 ;
559 nofFrontChPoints = 0; nofBackChPoints = 0;
560 for( Int_t ipointch9=0;ipointch9<fNofPoints[9];ipointch9++){
561 AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[9][ipointch9]->fFlags,chamber,detElemID);
562 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh9Z);
564 extrapCh9X = trigX1 * extrapCh9Z/trigZ1 ;
565 extrapCh9Y = trigY1 + (trigY2-trigY1) * (extrapCh9Z-trigZ1)/(trigZ2 - trigZ1) ;
568 HLTInfo("\textrap9 : (%f,%f,%f)\n",extrapCh9X,extrapCh9Y,extrapCh9Z);
571 if(nofBackChPoints < (fgkMaxNofTracks-1) &&
572 fabsf(extrapCh9X-fChPoint[9][ipointch9]->fX) < maxXDeflectionExtrap &&
573 fabsf(extrapCh9Y-fChPoint[9][ipointch9]->fY)/
574 ((fChPoint[9][ipointch9]->fX * fChPoint[9][ipointch9]->fX) +
575 (fChPoint[9][ipointch9]->fY * fChPoint[9][ipointch9]->fY)) <= extrapRatio ){
577 distChBack = sqrt((extrapCh9X-fChPoint[9][ipointch9]->fX)*(extrapCh9X-fChPoint[9][ipointch9]->fX)
578 + (extrapCh9Y-fChPoint[9][ipointch9]->fY)*(extrapCh9Y-fChPoint[9][ipointch9]->fY));
579 if(distChBack>circularWindow) continue;
582 HLTInfo("\t\t\tCh9 :%lf (%f,%f,%f) , circularWindow : %f, distChBack : %f\n",
583 fabsf(extrapCh9X-fChPoint[9][ipointch9]->fX)/distChBack,
584 fChPoint[9][ipointch9]->fX,fChPoint[9][ipointch9]->fY,fChPoint[9][ipointch9]->fZ,
585 circularWindow,distChBack);
588 backIndex[nofBackChPoints++] = ipointch9;
593 for( Int_t ipointch8=0;ipointch8<fNofPoints[8];ipointch8++){
594 AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[8][ipointch8]->fFlags,chamber,detElemID);
595 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh8Z);
597 /// extrapCh8Yref = trigY1 * extrapCh8Z/trigZ1 ;
598 extrapCh8X = trigX1 * extrapCh8Z/trigZ1 ;
599 extrapCh8Y = trigY1 + (trigY2-trigY1) * (extrapCh8Z-trigZ1)/(trigZ2 - trigZ1) ;
603 HLTInfo("\textrap8 : (%f,%f,%f), extrapRatio : %lf, maxXDeflectionExtrap : %lf\n",
604 extrapCh8X,extrapCh8Y,extrapCh8Z,extrapRatio,maxXDeflectionExtrap);
607 if( nofFrontChPoints < (fgkMaxNofTracks-1) &&
608 fabsf(extrapCh8X-fChPoint[8][ipointch8]->fX) < maxXDeflectionExtrap &&
609 fabsf(extrapCh8Y-fChPoint[8][ipointch8]->fY)/
610 ((fChPoint[8][ipointch8]->fX * fChPoint[8][ipointch8]->fX) +
611 (fChPoint[8][ipointch8]->fY * fChPoint[8][ipointch8]->fY)) <= extrapRatio ){
613 distChFront = sqrt((extrapCh8X-fChPoint[8][ipointch8]->fX)*(extrapCh8X-fChPoint[8][ipointch8]->fX)
614 + (extrapCh8Y-fChPoint[8][ipointch8]->fY)*(extrapCh8Y-fChPoint[8][ipointch8]->fY));
617 HLTInfo("\t\t\tCh8 : circularWindow : %f, distChFront : %f\n",circularWindow,distChFront);
620 if(distChFront>circularWindow) continue;
623 HLTInfo("\t\t\tCh8 :%lf (%f,%f,%f)\n",distChFront,
624 fChPoint[8][ipointch8]->fX,fChPoint[8][ipointch8]->fY,fChPoint[8][ipointch8]->fZ);
627 frontIndex[nofFrontChPoints++] = ipointch8;
632 minAngle = minAngleWindow;
633 p3.fX = trigX1 ; p3.fY = trigY1 ; p3.fZ = trigZ1 ;
634 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
635 Sub(&p3,fChPoint[9][backIndex[ibackpoint]],&pSeg2);
636 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
637 Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[8][frontIndex[ifrontpoint]],&pSeg1);
638 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2);
640 HLTInfo("\t\ttracklet-check-St5 : anglediff : %lf, minAngle : %lf\n",anglediff,minAngle);
642 if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
643 st5TrackletFound = true;
644 cells[1][fNofCells[1]].fFirst = frontIndex[ifrontpoint];
645 cells[1][fNofCells[1]].fSecond = backIndex[ibackpoint];
648 HLTInfo("\t\ttracklet-St5 : anglediff : %lf\n",anglediff);
649 HLTInfo("\t\t\tCh9 : (%f,%f,%f)\n",fChPoint[9][backIndex[ibackpoint]]->fX,
650 fChPoint[9][backIndex[ibackpoint]]->fY,fChPoint[9][backIndex[ibackpoint]]->fZ);
651 HLTInfo("\t\t\tCh8 : (%f,%f,%f)\n",fChPoint[8][frontIndex[ifrontpoint]]->fX,
652 fChPoint[8][frontIndex[ifrontpoint]]->fY,fChPoint[8][frontIndex[ifrontpoint]]->fZ);
654 }///anglediff condition
661 /// If tracklet not found, search for the single space point in Ch9 or in Ch8
662 if(!st5TrackletFound){
664 minAngle = minAngleWindow;
665 p3.fX = trigX2 ; p3.fY = trigY2 ; p3.fZ = trigZ2 ;
666 p2.fX = trigX1 ; p2.fY = trigY1 ; p2.fZ = trigZ1 ;
670 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
671 Sub(&p2,fChPoint[9][backIndex[ibackpoint]],&pSeg1);
672 anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
673 if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
674 ch9PointFound = true;
675 cells[1][fNofCells[1]].fFirst = -1;
676 cells[1][fNofCells[1]].fSecond = backIndex[ibackpoint];
681 HLTInfo("\t\tsppoint-Ch9 : anglediff : %lf\n\n",anglediff);
686 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
687 Sub(&p2,fChPoint[8][frontIndex[ifrontpoint]],&pSeg1);
688 anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
689 if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
690 ch8PointFound = true;
691 cells[1][fNofCells[1]].fFirst = frontIndex[ifrontpoint];
692 cells[1][fNofCells[1]].fSecond = -1;
697 HLTInfo("\t\tsppoint-Ch8 : anglediff : %lf\n\n",anglediff);
700 }///if no tracklets found condition
703 HLTInfo("\tnofTracks found after stn 5 : %d\n",fNofCells[1]);
706 if(!st5TrackletFound && !ch9PointFound && !ch8PointFound) continue;
708 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
711 /////////////////////////////////////////////////// Stn 4///////////////////////////////////////////////////////////////
713 /// extrapCh7X = trigX1 * extrapCh7Z/trigZ1 ;
714 /// extrapCh7Y = trigY1 + (trigY2-trigY1) * (extrapCh7Z-trigZ1)/(trigZ2 - trigZ1) ;
716 /// extrapCh6X = trigX1 * extrapCh6Z/trigZ1 ;
717 /// extrapCh6Y = trigY1 + (trigY2-trigY1) * (extrapCh6Z-trigZ1)/(trigZ2 - trigZ1) ;
720 nofFrontChPoints = 0; nofBackChPoints = 0;
721 for( Int_t ipointch7=0;ipointch7<fNofPoints[7];ipointch7++){
723 /// distChBack = (fChPoint[7][ipointch7]->fX * fChPoint[7][ipointch7]->fX)
724 /// + (fChPoint[7][ipointch7]->fY * fChPoint[7][ipointch7]->fY);
725 AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[7][ipointch7]->fFlags,chamber,detElemID);
726 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh7Z);
727 extrapCh7X = trigX1 * extrapCh7Z/trigZ1 ;
728 extrapCh7Y = trigY1 + (trigY2-trigY1) * (extrapCh7Z-trigZ1)/(trigZ2 - trigZ1) ;
731 HLTInfo("\textrap7 : (%f,%f,%f)\n",extrapCh7X,extrapCh7Y,extrapCh7Z);
734 if( nofBackChPoints < (fgkMaxNofTracks-1) &&
735 fabsf(extrapCh7X-fChPoint[7][ipointch7]->fX) < maxXDeflectionExtrap &&
736 fabsf(extrapCh7Y-fChPoint[7][ipointch7]->fY)/
737 ((fChPoint[7][ipointch7]->fX * fChPoint[7][ipointch7]->fX) +
738 (fChPoint[7][ipointch7]->fY * fChPoint[7][ipointch7]->fY)) <= extrapRatio ){
740 distChBack = sqrt((extrapCh7X-fChPoint[7][ipointch7]->fX)*(extrapCh7X-fChPoint[7][ipointch7]->fX)
741 + (extrapCh7Y-fChPoint[7][ipointch7]->fY)*(extrapCh7Y-fChPoint[7][ipointch7]->fY));
743 if(distChBack>circularWindow) continue;
746 HLTInfo("\t\t\tCh7 :%lf (%f,%f,%f)\n",distChBack,
747 fChPoint[7][ipointch7]->fX,fChPoint[7][ipointch7]->fY,fChPoint[7][ipointch7]->fZ);
750 backIndex[nofBackChPoints++] = ipointch7;
755 for( Int_t ipointch6=0;ipointch6<fNofPoints[6];ipointch6++){
757 AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[6][ipointch6]->fFlags,chamber,detElemID);
758 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh6Z);
760 extrapCh6X = trigX1 * extrapCh6Z/trigZ1 ;
761 extrapCh6Y = trigY1 + (trigY2-trigY1) * (extrapCh6Z-trigZ1)/(trigZ2 - trigZ1) ;
764 HLTInfo("\textrap6 : (%f,%f,%f)\n",extrapCh6X,extrapCh6Y,extrapCh6Z);
768 if(nofFrontChPoints < (fgkMaxNofTracks-1) &&
769 fabsf(extrapCh6X-fChPoint[6][ipointch6]->fX) < maxXDeflectionExtrap &&
770 fabsf(extrapCh6Y-fChPoint[6][ipointch6]->fY)/
771 ((fChPoint[6][ipointch6]->fX * fChPoint[6][ipointch6]->fX) +
772 (fChPoint[6][ipointch6]->fY * fChPoint[6][ipointch6]->fY)) <= extrapRatio ){
774 distChFront = sqrt((extrapCh6X-fChPoint[6][ipointch6]->fX)*(extrapCh6X-fChPoint[6][ipointch6]->fX)
775 + (extrapCh6Y-fChPoint[6][ipointch6]->fY)*(extrapCh6Y-fChPoint[6][ipointch6]->fY));
776 if(distChFront>circularWindow) continue;
779 HLTInfo("\t\t\tCh6 :%lf (%f,%f,%f)\n",distChFront,
780 fChPoint[6][ipointch6]->fX,fChPoint[6][ipointch6]->fY,fChPoint[6][ipointch6]->fZ);
783 frontIndex[nofFrontChPoints++] = ipointch6;
788 minAngle = minAngleWindow;
789 p3.fX = trigX1 ; p3.fY = trigY1 ; p3.fZ = trigZ1 ;
790 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
791 Sub(&p3,fChPoint[7][backIndex[ibackpoint]],&pSeg2);
792 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
793 Sub(fChPoint[7][backIndex[ibackpoint]],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
794 anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
796 HLTInfo("\t\ttracklet-check-St4 : anglediff : %lf, minAngle : %lf\n",anglediff,minAngle);
798 if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
800 st4TrackletFound = true;
802 cells[0][fNofCells[0]].fFirst = frontIndex[ifrontpoint];
803 cells[0][fNofCells[0]].fSecond = backIndex[ibackpoint];
806 HLTInfo("\t\ttracklet-St4 : anglediff : %lf\n",anglediff);
807 HLTInfo("\t\t\tCh7 : (%f,%f,%f)\n",fChPoint[7][backIndex[ibackpoint]]->fX,
808 fChPoint[7][backIndex[ibackpoint]]->fY,fChPoint[7][backIndex[ibackpoint]]->fZ);
809 HLTInfo("\t\t\tCh6 : (%f,%f,%f)\n",fChPoint[6][frontIndex[ifrontpoint]]->fX,
810 fChPoint[6][frontIndex[ifrontpoint]]->fY,fChPoint[6][frontIndex[ifrontpoint]]->fZ);
820 /// If tracklet not found search for the single space point in Ch7 or in Ch6
821 if(!st4TrackletFound){
823 minAngle = minAngleWindow;
824 p3.fX = trigX2 ; p3.fY = trigY2 ; p3.fZ = trigZ2 ;
825 p2.fX = trigX1 ; p2.fY = trigY1 ; p2.fZ = trigZ1 ;
829 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
830 Sub(&p2,fChPoint[7][backIndex[ibackpoint]],&pSeg1);
832 anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
833 if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
834 ch7PointFound = true;
835 cells[0][fNofCells[0]].fFirst = -1;
836 cells[0][fNofCells[0]].fSecond = backIndex[ibackpoint];
841 HLTInfo("\t\tsppoint-Ch7 : anglediff : %lf\n\n",anglediff);
846 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
847 Sub(&p2,fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
848 anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
849 if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
850 ch6PointFound = true;
851 cells[0][fNofCells[0]].fFirst = frontIndex[ifrontpoint];
852 cells[0][fNofCells[0]].fSecond = -1;
857 HLTInfo("\t\tsppoint-Ch6 : anglediff : %lf\n\n",anglediff);
860 }///if no tracklets found condition
863 HLTInfo("\tnofTracks found after stn 4 : %d\n",fNofCells[0]);
866 if(!st4TrackletFound && !ch7PointFound && !ch6PointFound) continue;
868 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
871 ////////////////////////////////////////////// Analyse and fill trackseg array////////////////////////////////////////
874 if(st5TrackletFound && st4TrackletFound){
876 minAngle = minAngleWindow;
878 for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
879 index1 = cells[0][itrackletfront].fFirst ;
880 index2 = cells[0][itrackletfront].fSecond ;
881 Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
882 for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
883 index3 = cells[1][itrackletback].fFirst ;
884 index4 = cells[1][itrackletback].fSecond ;
885 Sub(fChPoint[8][index3],fChPoint[7][index2],&pSeg2);
886 Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
887 anglediff = Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3);
888 if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
889 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
890 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
891 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
892 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
893 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
894 minAngle = anglediff;
896 HLTInfo("\t\ttracklet-St4 and St5 : anglediff : %lf\n",anglediff);
897 HLTInfo("\t\t\tCh9 : (%f,%f,%f)\n",fChPoint[9][index4]->fX,
898 fChPoint[9][index4]->fY,fChPoint[9][index4]->fZ);
899 HLTInfo("\t\t\tCh8 : (%f,%f,%f)\n",fChPoint[8][index3]->fX,
900 fChPoint[8][index3]->fY,fChPoint[8][index3]->fZ);
901 HLTInfo("\t\t\tCh7 : (%f,%f,%f)\n",fChPoint[7][index2]->fX,
902 fChPoint[7][index2]->fY,fChPoint[7][index2]->fZ);
903 HLTInfo("\t\t\tCh6 : (%f,%f,%f)\n",fChPoint[6][index1]->fX,
904 fChPoint[6][index1]->fY,fChPoint[6][index1]->fZ);
909 }///for loop of back ch
913 }else if(st5TrackletFound && (ch7PointFound || ch6PointFound)){
916 nofFrontChPoints = 0; nofBackChPoints = 0;
917 for( Int_t ifrontpoint=0;ifrontpoint<fNofCells[0];ifrontpoint++){
918 if(cells[0][ifrontpoint].fFirst==-1 && nofBackChPoints<(fgkMaxNofTracks-1))
919 backIndex[nofBackChPoints++] = cells[0][ifrontpoint].fSecond;
920 else if(nofFrontChPoints<(fgkMaxNofTracks-1))
921 frontIndex[nofFrontChPoints++] = cells[0][ifrontpoint].fFirst;
924 minAngle = minAngleWindow;
925 if(nofFrontChPoints>0 && nofBackChPoints>0){
927 for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
928 index3 = cells[1][itrackletback].fFirst ;
929 index4 = cells[1][itrackletback].fSecond ;
930 Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
931 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
932 Sub(fChPoint[8][index3],fChPoint[7][backIndex[ibackpoint]],&pSeg2);
933 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
934 Sub(fChPoint[7][backIndex[ibackpoint]],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
936 anglediff = TMath::RadToDeg()*(Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3))/2.0;
938 if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
939 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
940 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
941 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
942 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
943 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
944 minAngle = anglediff;
948 Sub(fChPoint[8][index3],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
949 anglediff1 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg3);
950 anglediff2 = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
952 if( anglediff1 < anglediff2 && anglediff1<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
953 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
954 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = -1;
955 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
956 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
957 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
958 minAngle = anglediff1;
962 if( anglediff2 < anglediff1 && anglediff2<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
963 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = -1;
964 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
965 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
966 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
967 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
968 minAngle = anglediff2;
971 }///loop of ifrontpoint
972 }///loop of ibackpoint
973 }/// for loop of St5 cells
974 }else if(nofFrontChPoints>0){
976 for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
977 index3 = cells[1][itrackletback].fFirst ;
978 index4 = cells[1][itrackletback].fSecond ;
979 Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
981 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
982 Sub(fChPoint[8][index3],fChPoint[6][frontIndex[ifrontpoint]],&pSeg2);
984 anglediff = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
985 if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
986 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
987 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = -1;
988 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
989 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
990 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
991 minAngle = anglediff;
996 }else{ /// if(nofBackChPoints>0){
997 for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
998 index3 = cells[1][itrackletback].fFirst ;
999 index4 = cells[1][itrackletback].fSecond ;
1000 Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
1002 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1003 Sub(fChPoint[8][index3],fChPoint[7][backIndex[ibackpoint]],&pSeg2);
1005 anglediff = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
1007 if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1008 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = -1;
1009 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
1010 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
1011 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
1012 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1013 minAngle = anglediff;
1018 }///condn for if(nofFrontChPoints>0)
1022 }else if((ch9PointFound || ch8PointFound) && st4TrackletFound){
1024 nofFrontChPoints = 0; nofBackChPoints = 0;
1025 for( Int_t ibackpoint=0;ibackpoint<fNofCells[1];ibackpoint++){
1026 if(cells[1][ibackpoint].fFirst==-1 && nofBackChPoints<(fgkMaxNofTracks-1))
1027 backIndex[nofBackChPoints++] = cells[1][ibackpoint].fSecond;
1028 else if(nofFrontChPoints<(fgkMaxNofTracks-1))
1029 frontIndex[nofFrontChPoints++] = cells[1][ibackpoint].fFirst;
1032 minAngle = minAngleWindow;
1033 if(nofFrontChPoints>0 && nofBackChPoints>0){
1035 for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1036 index1 = cells[0][itrackletfront].fFirst ;
1037 index2 = cells[0][itrackletfront].fSecond ;
1039 Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1041 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1044 Sub(fChPoint[8][frontIndex[ifrontpoint]],fChPoint[7][index2],&pSeg2);
1046 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1048 Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[8][frontIndex[ifrontpoint]],&pSeg3);
1050 anglediff = TMath::RadToDeg()*(Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3))/2.0;
1051 if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1052 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1053 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1054 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1055 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1056 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1057 minAngle = anglediff;
1061 Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[7][index2],&pSeg3);
1063 anglediff1 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
1064 anglediff2 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg3);
1066 if( anglediff1 < anglediff2 && anglediff1<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1067 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1068 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1069 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1070 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = -1;
1071 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1072 minAngle = anglediff1;
1076 if( anglediff2 < anglediff1 && anglediff2<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1077 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1078 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1079 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = -1;
1080 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1081 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1082 minAngle = anglediff2;
1085 }///loop of ifrontpoint
1086 }///loop of ibackpoint
1087 }/// for loop of St5 cells
1088 }else if(nofFrontChPoints>0){
1090 for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1091 index1 = cells[0][itrackletfront].fFirst ;
1092 index2 = cells[0][itrackletfront].fSecond ;
1094 Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1096 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1098 Sub(fChPoint[8][frontIndex[ifrontpoint]],fChPoint[7][index2],&pSeg2);
1100 anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
1101 if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1102 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1103 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;;
1104 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1105 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = -1;
1106 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1107 minAngle = anglediff;
1112 }else{ /// if(nofBackChPoints>0){
1114 for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1115 index1 = cells[0][itrackletfront].fFirst ;
1116 index2 = cells[0][itrackletfront].fSecond ;
1118 Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1120 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1121 Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[6][index1],&pSeg3);
1122 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg3);
1123 if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1124 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1125 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1126 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = -1;
1127 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1128 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1129 minAngle = anglediff;
1134 }///condn for if(nofFrontChPoints>0)
1138 }else if((ch9PointFound || ch8PointFound) && (ch7PointFound || ch6PointFound)){
1139 ///To Do : To be analysed for two points out of four slat chambers
1142 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1149 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
1152 if(fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[1]!=-1 ){
1153 meanX1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fX
1154 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fX)/2.0 ;
1155 meanY1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fY
1156 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fY)/2.0 ;
1157 meanZ1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fZ
1158 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fZ)/2.0 ;
1159 }else if(fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[1]==-1 ){
1160 meanX1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fX ;
1161 meanY1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fY ;
1162 meanZ1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fZ ;
1164 meanX1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fX ;
1165 meanY1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fY ;
1166 meanZ1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fZ ;
1169 if(fBackTrackSeg[ibacktrackseg].fIndex[2]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1 ){
1170 meanX2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fX
1171 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fX)/2.0 ;
1172 meanY2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fY
1173 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fY)/2.0 ;
1174 meanZ2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fZ
1175 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fZ)/2.0 ;
1176 }else if(fBackTrackSeg[ibacktrackseg].fIndex[2]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[3]==-1 ){
1177 meanX2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fX ;
1178 meanY2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fY ;
1179 meanZ2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fZ ;
1181 meanX2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fX ;
1182 meanY2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fY ;
1183 meanZ2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fZ ;
1186 fExtrapSt3X[ibacktrackseg] = meanX1 + (fgkTrackDetCoordinate[2]-meanZ1)*(meanX2-meanX1)/(meanZ2-meanZ1);
1187 fExtrapSt3Y[ibacktrackseg] = meanY1 + (fgkTrackDetCoordinate[2]-meanZ1)*(meanY2-meanY1)/(meanZ2-meanZ1);
1188 fInclinationBack[ibacktrackseg] = (meanX2-meanX1)/(meanZ2-meanZ1) ;
1189 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
1190 }///backtrigseg loop
1193 HLTInfo("AliHLTMUONFullTracker::SlatTrackSeg()--End\n");
1201 ///__________________________________________________________________________
1203 Bool_t AliHLTMUONFullTracker::QuadTrackSeg()
1205 ///Find the Track Segments in the Quadrant chambers
1207 Float_t meanX1,meanX2,meanY1,meanY2,meanZ1,meanZ2;
1208 Float_t expectSt3X,expectSt3Y,inclinationFront;
1210 AliHLTMUONRecHitStruct pSeg1,pSeg2,pSeg3;
1211 Double_t anglediff;///,anglediff1,anglediff2;
1212 Double_t minAngle = -1.0;
1213 Int_t indexMinAngleFront = -1;
1214 Int_t indexMinAngleBack = -1;
1215 Int_t backIndex = -1;
1217 Int_t ch3CellPoint[fgkMaxNofCellsPerCh],ch2CellPoint[fgkMaxNofCellsPerCh],nofSt2Cells=0;
1218 Int_t ch1CellPoint[fgkMaxNofCellsPerCh],ch0CellPoint[fgkMaxNofCellsPerCh],nofSt1Cells=0;
1219 Bool_t isConnected[fgkMaxNofTracks];
1221 Float_t distDiffX = 4.0; ///simulation result 4.0
1222 Float_t distDiffY = 10.0 ; ///simulation result 4.0
1223 ///float closeToY0 = 10.0; ///simulation result 10.0
1224 Float_t minAngleWindow = 2.0 ; ///simulation result 2.0
1225 Float_t diffDistStX = 25.0; ///simulation result 25.0
1226 Float_t diffDistStY = 75.0; ///simulation result 25.0
1227 Float_t st3WindowX = 40.0 ; ///simulation result 40.0
1228 Float_t st3WindowY = 10.0; ///simulation result 10.0
1229 /// Float_t inclinationWindow = 0.04; ///inclination window
1230 /// Float_t st3WindowXOp2 = 40.0 ; ///simulation result 40.0
1231 /// Float_t st3WindowYOp2 = 10.0; ///simulation result 10.0
1235 HLTInfo("AliHLTMUONFullTracker::QuadTrackSeg()--Begin\n\n");
1239 for( Int_t ibackpoint=0;ibackpoint<fNofPoints[3];ibackpoint++){
1240 for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[2];ifrontpoint++){
1242 if(fabsf(fChPoint[3][ibackpoint]->fX - fChPoint[2][ifrontpoint]->fX)<distDiffX
1243 && fabsf(fChPoint[3][ibackpoint]->fY - fChPoint[2][ifrontpoint]->fY)<distDiffY){
1245 /// if((fabsf(fChPoint[3][ibackpoint]->fY) > closeToY0 &&
1246 /// fabsf(fChPoint[3][ibackpoint]->fY) < fabsf(fChPoint[2][ifrontpoint]->fY)) ||
1247 /// nofSt2Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1249 if(nofSt2Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1252 HLTInfo("\t\t\tCh3 : %d, (%f,%f,%f)\n",
1253 nofSt2Cells,fChPoint[3][ibackpoint]->fX,fChPoint[3][ibackpoint]->fY,fChPoint[3][ibackpoint]->fZ);
1254 HLTInfo("\t\t\tCh2 :(%f,%f,%f)\n\n",
1255 fChPoint[2][ifrontpoint]->fX,fChPoint[2][ifrontpoint]->fY,fChPoint[2][ifrontpoint]->fZ);
1258 ch3CellPoint[nofSt2Cells] = ibackpoint;
1259 ch2CellPoint[nofSt2Cells] = ifrontpoint;
1266 for( Int_t ibackpoint=0;ibackpoint<fNofPoints[1];ibackpoint++){
1267 for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[0];ifrontpoint++){
1268 if(fabsf(fChPoint[1][ibackpoint]->fX - fChPoint[0][ifrontpoint]->fX)< distDiffX
1269 && fabsf(fChPoint[1][ibackpoint]->fY - fChPoint[0][ifrontpoint]->fY)<distDiffY){
1271 /// if((fabsf(fChPoint[1][ibackpoint]->fY) > closeToY0 &&
1272 /// fabsf(fChPoint[1][ibackpoint]->fY) < fabsf(fChPoint[0][ifrontpoint]->fY)) ||
1273 /// nofSt1Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1275 if(nofSt1Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1279 HLTInfo("\t\t\tCh1 : %d, (%f,%f,%f)\n",
1280 nofSt1Cells,fChPoint[1][ibackpoint]->fX,fChPoint[1][ibackpoint]->fY,fChPoint[1][ibackpoint]->fZ);
1281 HLTInfo("\t\t\tCh0 :(%f,%f,%f)\n\n",
1282 fChPoint[0][ifrontpoint]->fX,fChPoint[0][ifrontpoint]->fY,fChPoint[0][ifrontpoint]->fZ);
1284 ch1CellPoint[nofSt1Cells] = ibackpoint;
1285 ch0CellPoint[nofSt1Cells] = ifrontpoint;
1292 HLTInfo("\tnofSt1Cells : %d, nofSt2Cells : %d\n",nofSt1Cells,nofSt2Cells);
1295 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++)
1296 isConnected[ibacktrackseg] = false;
1299 ///First Check for Tracklets in two front stations
1300 for( Int_t itrackletfront=0;itrackletfront<nofSt1Cells;itrackletfront++){
1302 Sub(fChPoint[1][ch1CellPoint[itrackletfront]],fChPoint[0][ch0CellPoint[itrackletfront]],&pSeg1);
1304 minAngle = minAngleWindow;
1305 indexMinAngleBack = -1;
1306 indexMinAngleFront = -1;
1308 meanX1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fX
1309 + fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/2.0 ;
1310 meanY1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fY
1311 + fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/2.0 ;
1312 meanZ1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fZ
1313 + fChPoint[1][ch1CellPoint[itrackletfront]]->fZ)/2.0 ;
1315 for( Int_t itrackletback=0;itrackletback<nofSt2Cells;itrackletback++){
1318 <<fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)
1320 <<fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)
1323 if(fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fX -
1324 fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX ||
1325 fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fY -
1326 fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1328 meanX2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fX
1329 + fChPoint[3][ch3CellPoint[itrackletback]]->fX)/2.0 ;
1330 meanY2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fY
1331 + fChPoint[3][ch3CellPoint[itrackletback]]->fY)/2.0 ;
1332 meanZ2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fZ
1333 + fChPoint[3][ch3CellPoint[itrackletback]]->fZ)/2.0 ;
1335 expectSt3X = meanX2 + (fgkTrackDetCoordinate[2]-meanZ2)*(meanX2-meanX1)/(meanZ2-meanZ1);
1336 expectSt3Y = meanY2 + (fgkTrackDetCoordinate[2]-meanZ2)*(meanY2-meanY1)/(meanZ2-meanZ1);
1337 inclinationFront = (meanX2-meanX1)/(meanZ2-meanZ1) ;
1339 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
1341 if( /// fabsf(inclinationBack[ibacktrackseg]-inclinationFront)<0.04 &&
1342 fabsf((expectSt3X-fExtrapSt3X[ibacktrackseg])) < st3WindowX
1343 && fabsf((expectSt3Y-fExtrapSt3Y[ibacktrackseg])) < st3WindowY){
1345 Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1346 Sub(fChPoint[3][ch3CellPoint[itrackletback]],fChPoint[2][ch2CellPoint[itrackletback]],&pSeg3);
1348 anglediff = TMath::RadToDeg()* (Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3));
1351 HLTInfo("\t\t\tanglediff : %lf\n",anglediff);
1353 if(anglediff<minAngle){
1354 minAngle = anglediff;
1355 indexMinAngleBack = itrackletback;
1356 indexMinAngleFront = itrackletfront;
1357 backIndex = ibacktrackseg;
1358 isConnected[ibacktrackseg] = true;
1360 }///matching tracklet
1361 }///for loop of backtrackseg
1365 if(minAngle < minAngleWindow && indexMinAngleFront!=-1
1366 && indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)){
1368 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = ch0CellPoint[indexMinAngleFront];
1369 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = ch1CellPoint[indexMinAngleFront];
1370 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = ch2CellPoint[indexMinAngleBack];
1371 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = ch3CellPoint[indexMinAngleBack];
1373 fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;
1374 fNoffrontTrackSeg++;
1375 }///condition to find valid tracklet
1377 }///for loop of front ch
1379 Int_t nofNCfBackTrackSeg = 0;
1380 Int_t ncfBackTrackSeg[fgkMaxNofTracks];
1382 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++)
1383 if(!isConnected[ibacktrackseg])
1384 ncfBackTrackSeg[nofNCfBackTrackSeg++] = ibacktrackseg;
1390 HLTInfo("\tfNofConnected : %d, nofNCfBackTrackSeg : %d\n",fNofConnected,nofNCfBackTrackSeg);
1391 HLTInfo("\tfNofPoints[3] : %d, fNofPoints[2] : %d\n",fNofPoints[3],fNofPoints[2]);
1392 if(nofNCfBackTrackSeg==0){
1393 HLTInfo("All fBackTrackSegs are connected with fFrontTrackSegs, no need to search further\n");
1394 HLTInfo("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
1398 if(nofNCfBackTrackSeg==0) return true;
1401 ///Next Check for tracklet in Stn1 and space point in Stn2
1402 Bool_t isbackpoint=false,isfrontpoint=false;
1403 for( Int_t itrackletfront=0;itrackletfront<nofSt1Cells;itrackletfront++){
1404 Sub(fChPoint[1][ch1CellPoint[itrackletfront]],fChPoint[0][ch0CellPoint[itrackletfront]],&pSeg1);
1405 minAngle = minAngleWindow;
1406 indexMinAngleBack = -1;
1407 indexMinAngleFront = -1;
1409 for( Int_t ibackpoint=0;ibackpoint<fNofPoints[3];ibackpoint++){
1410 if(/// hasCh3Cells[ibackpoint] == true &&
1411 fabsf(fChPoint[3][ibackpoint]->fX -
1412 fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX ||
1413 fabsf(fChPoint[3][ibackpoint]->fY -
1414 fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1416 expectSt3X = fChPoint[3][ibackpoint]->fX + (fgkTrackDetCoordinate[2] - fChPoint[3][ibackpoint]->fZ)*
1417 (fChPoint[3][ibackpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1418 (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1419 expectSt3Y = fChPoint[3][ibackpoint]->fY + (fgkTrackDetCoordinate[2] - fChPoint[3][ibackpoint]->fZ)*
1420 (fChPoint[3][ibackpoint]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/
1421 (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1422 inclinationFront = (fChPoint[3][ibackpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1423 (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ) ;
1425 for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++){
1427 if(/// fabsf(inclinationBack[ncfBackTrackSeg[ibacktrackseg]]-inclinationFront)< inclinationWindow &&
1428 fabsf((expectSt3X-fExtrapSt3X[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1429 && fabsf((expectSt3Y-fExtrapSt3Y[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowY){
1431 Sub(fChPoint[3][ibackpoint],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1433 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
1435 HLTInfo("\t\t annglediff(Ch4) : %lf\n",anglediff);
1437 if(anglediff<minAngle){
1438 minAngle = anglediff;
1439 indexMinAngleBack = ibackpoint;
1440 indexMinAngleFront = itrackletfront;
1441 backIndex = ibacktrackseg;
1442 isConnected[ncfBackTrackSeg[ibacktrackseg]] = true;
1444 isfrontpoint = false;
1446 }///matching tracklet
1447 }///loop on not found back trackseg
1450 for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[2];ifrontpoint++){
1451 if(/// hasCh2Cells[ifrontpoint] == true &&
1452 fabsf(fChPoint[2][ifrontpoint]->fX -
1453 fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX ||
1454 fabsf(fChPoint[2][ifrontpoint]->fY -
1455 fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1457 expectSt3X = fChPoint[2][ifrontpoint]->fX + (fgkTrackDetCoordinate[2] - fChPoint[2][ifrontpoint]->fZ)*
1458 (fChPoint[2][ifrontpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1459 (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1460 expectSt3Y = fChPoint[2][ifrontpoint]->fY + (fgkTrackDetCoordinate[2] - fChPoint[2][ifrontpoint]->fZ)*
1461 (fChPoint[2][ifrontpoint]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/
1462 (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1463 inclinationFront = (fChPoint[2][ifrontpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1464 (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ) ;
1466 for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++){
1468 if( /// fabsf(inclinationBack[ncfBackTrackSeg[ibacktrackseg]]-inclinationFront)< inclinationWindow &&
1469 fabsf((expectSt3X-fExtrapSt3X[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1470 && fabsf((expectSt3Y-fExtrapSt3Y[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowY){
1472 Sub(fChPoint[2][ifrontpoint],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1474 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
1476 HLTInfo("\t\t annglediff(Ch3) : %lf\n",anglediff);
1478 if(anglediff<minAngle){
1479 minAngle = anglediff;
1480 indexMinAngleBack = ifrontpoint;
1481 indexMinAngleFront = itrackletfront;
1482 backIndex = ibacktrackseg;
1483 isConnected[ncfBackTrackSeg[ibacktrackseg]] = true;
1484 isbackpoint = false;
1485 isfrontpoint = true;
1488 }///matching tracklet
1489 }///loop on not found back trackseg
1492 if(minAngle < minAngleWindow && indexMinAngleFront!=-1 &&
1493 indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)){
1495 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = ch0CellPoint[indexMinAngleFront];
1496 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = ch1CellPoint[indexMinAngleFront];
1497 if(isfrontpoint && !isbackpoint){
1498 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = indexMinAngleBack;
1499 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = -1;
1501 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = -1;
1502 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = indexMinAngleBack;
1504 fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;
1505 fNoffrontTrackSeg++;
1506 }///condition to find valid tracklet
1511 Int_t nofSNCfBackTrackSeg = 0;
1512 Int_t sncfBackTrackSeg[fgkMaxNofTracks];
1514 for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++)
1515 if(!isConnected[ncfBackTrackSeg[ibacktrackseg]])
1516 sncfBackTrackSeg[nofSNCfBackTrackSeg++] = ncfBackTrackSeg[ibacktrackseg];
1521 HLTInfo("\tfNofConnected : %d, nofSNCfBackTrackSeg : %d\n",fNofConnected,nofSNCfBackTrackSeg);
1522 if(nofSNCfBackTrackSeg==0){
1523 HLTInfo("All fBackTrackSegs are connected with fFrontTrackSegs, no need to search further\n");
1524 HLTInfo("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
1528 if(nofSNCfBackTrackSeg==0) return true;
1530 ///Last Check for tracklet in Stn2 and space point in Stn1
1531 for( Int_t itrackletback=0;itrackletback<nofSt2Cells;itrackletback++){
1532 Sub(fChPoint[3][ch3CellPoint[itrackletback]],fChPoint[2][ch2CellPoint[itrackletback]],&pSeg1);
1533 minAngle = minAngleWindow ;
1534 indexMinAngleBack = -1;
1535 indexMinAngleFront = -1;
1537 for( Int_t ibackpoint=0;ibackpoint<fNofPoints[1];ibackpoint++){
1538 if(/// hasCh1Cells[ibackpoint] == true &&
1539 fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fX -
1540 fChPoint[1][ibackpoint]->fX) > diffDistStX ||
1541 fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fY -
1542 fChPoint[1][ibackpoint]->fY) > diffDistStY) continue;
1544 expectSt3X = fChPoint[2][ch2CellPoint[itrackletback]]->fX + (fgkTrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
1545 (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ibackpoint]->fX)/
1546 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ);
1547 expectSt3Y = fChPoint[2][ch2CellPoint[itrackletback]]->fY + (fgkTrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
1548 (fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[1][ibackpoint]->fY)/
1549 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ);
1550 inclinationFront = (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ibackpoint]->fX)/
1551 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ) ;
1553 for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++){
1555 if(/// fabsf(inclinationBack[sncfBackTrackSeg[ibacktrackseg]]-inclinationFront)<inclinationWindow &&
1556 fabsf((expectSt3X-fExtrapSt3X[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1557 && fabsf((expectSt3Y-fExtrapSt3Y[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowY ){
1559 Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[1][ibackpoint],&pSeg2);
1561 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
1562 if(anglediff<minAngle){
1563 minAngle = anglediff;
1564 indexMinAngleBack = itrackletback;
1565 indexMinAngleFront = ibackpoint;
1566 backIndex = ibacktrackseg;
1567 isConnected[sncfBackTrackSeg[ibacktrackseg]] = true;
1569 isfrontpoint = false;
1571 }///matching tracklet
1572 }///loop on not found back trackseg
1575 for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[0];ifrontpoint++){
1576 if(/// hasCh0Cells[ifrontpoint] == true &&
1577 fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fX -
1578 fChPoint[0][ifrontpoint]->fX) > diffDistStX ||
1579 fabsf(fChPoint[2][ch2CellPoint[itrackletback]]->fY -
1580 fChPoint[0][ifrontpoint]->fY) > diffDistStY ) continue;
1582 expectSt3X = fChPoint[2][ch2CellPoint[itrackletback]]->fX + (fgkTrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
1583 (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[0][ifrontpoint]->fX)/
1584 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ);
1585 expectSt3Y = fChPoint[2][ch2CellPoint[itrackletback]]->fY + (fgkTrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
1586 (fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[0][ifrontpoint]->fY)/
1587 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ);
1588 inclinationFront = (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[0][ifrontpoint]->fX)/
1589 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ) ;
1591 for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++){
1593 if(/// fabsf(inclinationBack[sncfBackTrackSeg[ibacktrackseg]]-inclinationFront)<inclinationWindow &&
1594 fabsf((expectSt3X-fExtrapSt3X[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1595 && fabsf((expectSt3Y-fExtrapSt3Y[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowY ){
1597 Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[0][ifrontpoint],&pSeg2);
1599 anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2) ;
1600 if(anglediff<minAngle){
1601 minAngle = anglediff;
1602 indexMinAngleBack = itrackletback;
1603 indexMinAngleFront = ifrontpoint;
1604 backIndex = ibacktrackseg;
1605 isConnected[sncfBackTrackSeg[ibacktrackseg]] = true;
1606 isbackpoint = false;
1607 isfrontpoint = true;
1610 }///matching tracklet
1611 }///loop on not found back trackseg
1614 if(minAngle < minAngleWindow && indexMinAngleFront!=-1 &&
1615 indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)){
1618 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = indexMinAngleFront;
1619 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = -1;
1621 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = -1;
1622 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = indexMinAngleFront;
1625 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = ch2CellPoint[indexMinAngleBack];
1626 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = ch3CellPoint[indexMinAngleBack];
1628 fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;
1629 fNoffrontTrackSeg++;
1630 }///condition to find valid tracklet
1634 for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++)
1635 if(isConnected[sncfBackTrackSeg[ibacktrackseg]])
1639 HLTInfo("\tfNofConnected : %d\n",fNofConnected);
1640 HLTInfo("Three spacepoints are found in fFrontTrackSegs\n");
1641 HLTInfo("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
1648 ///__________________________________________________________________________
1650 Double_t AliHLTMUONFullTracker::KalmanFilter(AliMUONTrackParam &trackParamAtCluster, Cluster *cluster)
1652 //// Compute new track parameters and their covariances including new cluster using kalman filter
1653 //// return the additional track chi2
1655 #ifdef PRINT_DETAIL_KALMAN
1656 HLTInfo("AliHLTMUONFullTracker::KalmanFilter()--Begin\n\n");
1660 /// Get actual track parameters (p)
1661 TMatrixD param(trackParamAtCluster.GetParameters());
1662 #ifdef PRINT_DETAIL_KALMAN
1663 Info("\tKalmanFilter","param.Print() [p]");
1665 HLTInfo("GetZ : %lf\n",trackParamAtCluster.GetZ());
1670 /// Get new cluster parameters (m)
1671 ///AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
1672 TMatrixD clusterParam(5,1);
1673 clusterParam.Zero();
1674 /// clusterParam(0,0) = cluster->GetX();
1675 /// clusterParam(2,0) = cluster->GetY();
1676 clusterParam(0,0) = cluster->fX;
1677 clusterParam(2,0) = cluster->fY;
1678 #ifdef PRINT_DETAIL_KALMAN
1679 Info("\tKalmanFilter","clusterParam.Print() [m]");
1680 clusterParam.Print();
1685 /// Compute the actual parameter weight (W)
1686 TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
1687 #ifdef PRINT_DETAIL_KALMAN
1688 Info("\tKalmanFilter","Covariance : [C]");
1689 paramWeight.Print();
1692 if (paramWeight.Determinant() != 0) {
1693 paramWeight.Invert();
1695 Warning("KalmanFilter"," Determinant = 0");
1699 #ifdef PRINT_DETAIL_KALMAN
1700 Info("\tKalmanFilter","Weight Matrix inverse of Covariance [W = C^-1]");
1701 paramWeight.Print();
1705 /// Compute the new cluster weight (U)
1706 TMatrixD clusterWeight(5,5);
1707 clusterWeight.Zero();
1708 clusterWeight(0,0) = 1. / cluster->fErrX2;
1709 clusterWeight(2,2) = 1. / cluster->fErrY2;
1710 #ifdef PRINT_DETAIL_KALMAN
1711 Info("\tKalmanFilter","clusterWeight.Print() [U]");
1712 HLTInfo("\tErrX2 : %lf, ErrY2 : %lf\n",cluster->fErrX2,cluster->fErrY2);
1713 clusterWeight.Print();
1719 /// Compute the new parameters covariance matrix ( (W+U)^-1 )
1720 TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
1721 #ifdef PRINT_DETAIL_KALMAN
1722 Info("\tKalmanFilter","newParamCov.Print() [(W+U)]");
1723 newParamCov.Print();
1725 if (newParamCov.Determinant() != 0) {
1726 newParamCov.Invert();
1728 Warning("RunKalmanFilter"," Determinant = 0");
1731 #ifdef PRINT_DETAIL_KALMAN
1732 Info("\tKalmanFilter","newParamCov.Print() [(W+U)^-1] (new covariances[W] for trackParamAtCluster)");
1733 newParamCov.Print();
1736 /// Save the new parameters covariance matrix
1737 trackParamAtCluster.SetCovariances(newParamCov);
1739 /// Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
1740 TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
1741 TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); /// U(m-p)
1742 TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); /// ((W+U)^-1)U(m-p)
1743 newParam += param; /// ((W+U)^-1)U(m-p) + p
1744 #ifdef PRINT_DETAIL_KALMAN
1745 Info("\tKalmanFilter","newParam.Print() [p' = ((W+U)^-1)U(m-p) + p] (new parameters[p] for trackParamAtCluster)");
1749 /// Save the new parameters
1750 trackParamAtCluster.SetParameters(newParam);
1751 /// HLTInfo(Form("Pt : %lf\n",TMath::Sqrt(trackParamAtCluster.Px()*trackParamAtCluster.Px() +
1752 /// trackParamAtCluster.Py()*trackParamAtCluster.Py())));
1755 /// Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
1756 tmp = newParam; /// p'
1757 tmp -= param; /// (p'-p)
1758 TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); /// W(p'-p)
1759 TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); /// ((p'-p)^-1)W(p'-p)
1760 tmp = newParam; /// p'
1761 tmp -= clusterParam; /// (p'-m)
1762 TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); /// U(p'-m)
1763 addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); /// ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
1765 #ifdef PRINT_DETAIL_KALMAN
1766 Info("\tKalmanFilter","addChi2Track.Print() [additional chi2 = ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)))]");
1767 addChi2Track.Print();
1768 HLTInfo("AliHLTMUONFullTracker::KalmanFilter()--End\n\n");
1772 return addChi2Track(0,0);
1775 ///__________________________________________________________________________
1777 Double_t AliHLTMUONFullTracker::TryOneCluster(const AliMUONTrackParam &trackParam, Cluster* cluster,
1778 AliMUONTrackParam &trackParamAtCluster, Bool_t updatePropagator)
1780 //// Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
1781 //// return the corresponding Chi2
1782 //// return trackParamAtCluster
1784 /// extrapolate track parameters and covariances at the z position of the tested cluster
1785 trackParamAtCluster = trackParam;
1786 AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->fZ, updatePropagator);
1788 /// set pointer to cluster into trackParamAtCluster
1789 ///trackParamAtCluster.SetClusterPtr(cluster);
1791 /// Set differences between trackParam and cluster in the bending and non bending directions
1792 Double_t dX = cluster->fX - trackParamAtCluster.GetNonBendingCoor();
1793 Double_t dY = cluster->fY - trackParamAtCluster.GetBendingCoor();
1795 /// Calculate errors and covariances
1796 const TMatrixD& kParamCov = trackParamAtCluster.GetCovariances();
1797 Double_t sigmaX2 = kParamCov(0,0) + cluster->fErrX2;
1798 Double_t sigmaY2 = kParamCov(2,2) + cluster->fErrY2;
1801 return dX * dX / sigmaX2 + dY * dY / sigmaY2;
1805 ///__________________________________________________________________________
1807 Bool_t AliHLTMUONFullTracker::TryOneClusterFast(const AliMUONTrackParam &trackParam, const Cluster* cluster) const
1809 //// Test the compatibility between the track and the cluster
1810 //// given the track resolution + the maximum-distance-to-track value
1811 //// and assuming linear propagation of the track:
1812 //// return kTRUE if they are compatibles
1814 Float_t sigmaCutForTracking = 6.0;
1815 Float_t maxNonBendingDistanceToTrack = 1.0;
1816 Float_t maxBendingDistanceToTrack = 1.0;
1818 Double_t dZ = cluster->fZ - trackParam.GetZ();
1819 Double_t dX = cluster->fX - (trackParam.GetNonBendingCoor() + trackParam.GetNonBendingSlope() * dZ);
1820 Double_t dY = cluster->fY - (trackParam.GetBendingCoor() + trackParam.GetBendingSlope() * dZ);
1821 const TMatrixD& kParamCov = trackParam.GetCovariances();
1822 Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1);
1823 Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3);
1825 Double_t dXmax = sigmaCutForTracking * TMath::Sqrt(errX2) +
1826 maxNonBendingDistanceToTrack;
1827 Double_t dYmax = sigmaCutForTracking * TMath::Sqrt(errY2) +
1828 maxBendingDistanceToTrack;
1830 if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE;
1835 ///__________________________________________________________________________
1837 void AliHLTMUONFullTracker::PropagateTracks(Double_t charge, Float_t& px, Float_t& py, Float_t& pz,
1838 Float_t& xr, Float_t& yr, Float_t& zr, Float_t zprop) const
1841 /// propagate in magnetic field between hits of indices i1 and i2
1844 Double_t vect[7], vout[7];
1845 Double_t step = -5.0;
1846 Double_t zMax = zprop+.5;
1851 vect[6] = sqrt((px)*(px) + (py)*(py) + (pz)*(pz));
1852 vect[3] = px/vect[6];
1853 vect[4] = py/vect[6];
1854 vect[5] = pz/vect[6];
1855 /// cout<<"vec[2] : "<<vect[2]<<", zMax : "<<zMax<<endl;
1858 while ((vect[2] < zMax) && (nSteps < 10000)) {
1860 /// OneStepRungekutta(charge, step, vect, vout);
1861 OneStepHelix3(charge,step,vect,vout);
1862 ///SetPoint(fCount,vout[0],vout[1],vout[2]);
1864 /// HLTInfo("(x,y,z) : (%f,%f,%f)\n",vout[0],vout[1],vout[2]);
1865 for (Int_t i = 0; i < 7; i++) {
1874 px = vect[3]*vect[6];
1875 py = vect[4]*vect[6];
1876 pz = vect[5]*vect[6];
1880 ///__________________________________________________________________________
1882 void AliHLTMUONFullTracker::OneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout) const
1885 //// ******************************************************************
1887 //// * Tracking routine in a constant field oriented *
1888 //// * along axis 3 *
1889 //// * Tracking is performed with a conventional *
1890 //// * helix step method *
1892 //// * ==>Called by : USER, GUSWIM *
1893 //// * Authors R.Brun, M.Hansroul ********* *
1894 //// * Rewritten V.Perevoztchikov
1896 //// ******************************************************************
1900 Double_t h4, hp, rho, tet;
1901 Double_t sint, sintt, tsint, cos1t;
1902 Double_t f1, f2, f3, f4, f5, f6;
1904 const Int_t kix = 0;
1905 const Int_t kiy = 1;
1906 const Int_t kiz = 2;
1907 const Int_t kipx = 3;
1908 const Int_t kipy = 4;
1909 const Int_t kipz = 5;
1910 const Int_t kipp = 6;
1912 const Double_t kec = 2.9979251e-4;
1915 /// ------------------------------------------------------------------
1917 /// units are kgauss,centimeters,gev/c
1919 vout[kipp] = vect[kipp];
1922 hxp[0] = - vect[kipy];
1923 hxp[1] = + vect[kipx];
1927 rho = -h4/vect[kipp];
1929 if (TMath::Abs(tet) > 0.15) {
1930 sint = TMath::Sin(tet);
1932 tsint = (tet-sint)/tet;
1933 cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet;
1935 tsint = tet*tet/36.;
1936 sintt = (1. - tsint);
1943 f3 = step * tsint * hp;
1946 f6 = tet * cos1t * hp;
1948 vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0];
1949 vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1];
1950 vout[kiz] = vect[kiz] + f1*vect[kipz] + f3;
1952 vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0];
1953 vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1];
1954 vout[kipz] = vect[kipz] + f4*vect[kipz] + f6;
1959 ///__________________________________________________________________________
1961 ///______________________________________________________________________________
1963 void AliHLTMUONFullTracker::OneStepRungekutta(Double_t charge, Double_t step,
1964 const Double_t* vect, Double_t* vout) const
1966 //// ******************************************************************
1968 //// * Runge-Kutta method for tracking a particle through a magnetic *
1969 //// * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
1970 //// * Standards, procedure 25.5.20) *
1972 //// * Input parameters *
1973 //// * CHARGE Particle charge *
1974 //// * STEP Step size *
1975 //// * VECT Initial co-ords,direction cosines,momentum *
1976 //// * Output parameters *
1977 //// * VOUT Output co-ords,direction cosines,momentum *
1978 //// * User routine called *
1979 //// * CALL GUFLD(X,F) *
1981 //// * ==>Called by : <USER>, GUSWIM *
1982 //// * Authors R.Brun, M.Hansroul ********* *
1983 //// * V.Perevoztchikov (CUT STEP implementation) *
1986 //// ******************************************************************
1988 Double_t h2, h4, f[4];
1989 Double_t xyzt[3], a, b, c, ph,ph2;
1990 Double_t secxs[4],secys[4],seczs[4],hxp[3];
1991 Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
1992 Double_t est, at, bt, ct, cba;
1993 Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
2003 Double_t maxit = 1992;
2004 Double_t maxcut = 11;
2006 const Double_t kdlt = 1e-4;
2007 const Double_t kdlt32 = kdlt/32.;
2008 const Double_t kthird = 1./3.;
2009 const Double_t khalf = 0.5;
2010 const Double_t kec = 2.9979251e-4;
2012 const Double_t kpisqua = 9.86960440109;
2013 const Int_t kix = 0;
2014 const Int_t kiy = 1;
2015 const Int_t kiz = 2;
2016 const Int_t kipx = 3;
2017 const Int_t kipy = 4;
2018 const Int_t kipz = 5;
2021 /// *. ------------------------------------------------------------------
2023 /// * this constant is for units cm,gev/c and kgauss
2027 for(Int_t j = 0; j < 7; j++)
2030 Double_t pinv = kec * charge / vect[6];
2038 if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
2039 ///cmodif: call gufld(vout,f) changed into:
2040 TGeoGlobalMagField::Instance()->Field(vout,f);
2043 /// * start of integration
2056 secxs[0] = (b * f[2] - c * f[1]) * ph2;
2057 secys[0] = (c * f[0] - a * f[2]) * ph2;
2058 seczs[0] = (a * f[1] - b * f[0]) * ph2;
2059 ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
2060 if (ang2 > kpisqua) break;
2062 dxt = h2 * a + h4 * secxs[0];
2063 dyt = h2 * b + h4 * secys[0];
2064 dzt = h2 * c + h4 * seczs[0];
2069 /// * second intermediate point
2072 est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
2074 if (ncut++ > maxcut) break;
2083 ///cmodif: call gufld(xyzt,f) changed into:
2084 TGeoGlobalMagField::Instance()->Field(xyzt,f);
2090 secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
2091 secys[1] = (ct * f[0] - at * f[2]) * ph2;
2092 seczs[1] = (at * f[1] - bt * f[0]) * ph2;
2096 secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
2097 secys[2] = (ct * f[0] - at * f[2]) * ph2;
2098 seczs[2] = (at * f[1] - bt * f[0]) * ph2;
2099 dxt = h * (a + secxs[2]);
2100 dyt = h * (b + secys[2]);
2101 dzt = h * (c + seczs[2]);
2105 at = a + 2.*secxs[2];
2106 bt = b + 2.*secys[2];
2107 ct = c + 2.*seczs[2];
2109 est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
2110 if (est > 2.*TMath::Abs(h)) {
2111 if (ncut++ > maxcut) break;
2120 ///cmodif: call gufld(xyzt,f) changed into:
2121 TGeoGlobalMagField::Instance()->Field(xyzt,f);
2123 z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
2124 y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
2125 x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
2127 secxs[3] = (bt*f[2] - ct*f[1])* ph2;
2128 secys[3] = (ct*f[0] - at*f[2])* ph2;
2129 seczs[3] = (at*f[1] - bt*f[0])* ph2;
2130 a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
2131 b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
2132 c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
2134 est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
2135 + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
2136 + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
2138 if (est > kdlt && TMath::Abs(h) > 1.e-4) {
2139 if (ncut++ > maxcut) break;
2145 /// * if too many iterations, go to helix
2146 if (iter++ > maxit) break;
2151 cba = 1./ TMath::Sqrt(a*a + b*b + c*c);
2159 if (step < 0.) rest = -rest;
2160 if (rest < 1.e-5*TMath::Abs(step)) return;
2164 /// angle too big, use helix
2169 f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
2178 hxp[0] = f2*vect[kipz] - f3*vect[kipy];
2179 hxp[1] = f3*vect[kipx] - f1*vect[kipz];
2180 hxp[2] = f1*vect[kipy] - f2*vect[kipx];
2182 hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
2185 sint = TMath::Sin(tet);
2186 cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
2190 g3 = (tet-sint) * hp*rho1;
2195 vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
2196 vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
2197 vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
2199 vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
2200 vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
2201 vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
2206 ///______________________________________________________________________________
2209 Bool_t AliHLTMUONFullTracker::SelectFront()
2211 // Track extrapolation through dipole magnet to connect front and back track seg.
2213 Cluster clus1,clus2;
2214 Int_t minIndex=0,maxIndex=0;
2215 Int_t minCh=0,maxCh=0;
2216 Int_t ifronttrackseg=0;
2218 Float_t xf,yf,zf,thetaDev;
2219 Float_t p,spx,spy,spz,px,py,pz,sx,sy,sz,x,y,z;
2221 Float_t dist,mindist;
2222 Int_t frontsegIndex ;
2224 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
2226 if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
2228 /// if(fBackTrackSeg[ibacktrackseg].fIndex[2]==-1 || fBackTrackSeg[ibacktrackseg].fIndex[3]==-1) continue;
2230 ifronttrackseg = fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]-1];
2233 HLTInfo("AliHLTMUONFullTracker::SelectFront()--Begin\n\n");
2234 HLTInfo("\tbacktrack : %d is connected with : %d front tracks\n",
2235 ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg]);
2238 maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
2239 maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
2241 minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
2242 minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
2246 clus2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
2247 clus2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
2248 clus2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
2250 clus1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
2251 clus1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
2252 clus1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
2254 zf = 0.5*(AliMUONConstants::DefaultChamberZ(4) + AliMUONConstants::DefaultChamberZ(5));
2255 thetaDev= (1/zf)*(clus1.fY*clus2.fZ - clus2.fY*clus1.fZ)/(clus2.fZ - clus1.fZ);
2256 xf = clus1.fX*zf/clus1.fZ;
2257 yf = clus2.fY - ((clus2.fY - clus1.fY)*(clus2.fZ-zf))/(clus2.fZ - clus1.fZ);
2258 p = 3.0*0.3/thetaDev;
2259 charge = (p>0)?-1:+1;
2261 fCharge[ibacktrackseg] = charge;
2263 /// cout<<"Charge : "<<charge<<endl;
2264 ///fTrackParam[ibacktrackseg].SetCharge(charge);
2267 spz = sqrt((p*p-(spx*spx + spy*spy)));
2271 sx = clus1.fX; sy = clus1.fY; sz = clus1.fZ;
2274 for( Int_t iconnected=0;iconnected<fNofConnectedfrontTrackSeg[ibacktrackseg];iconnected++){
2276 ifronttrackseg = fBackToFront[ibacktrackseg][iconnected];
2278 minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2279 minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2280 clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2281 clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2282 clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2284 x = sx; y = sy; z = sz;
2285 px = spx; py = spy; pz = spz;
2286 PropagateTracks(charge,px,py,pz,x,y,z,clus1.fZ);
2288 dist = sqrt((clus1.fX-x)*(clus1.fX-x) +
2289 (clus1.fY-y)*(clus1.fY-y));
2292 frontsegIndex = ifronttrackseg;
2294 }///for loop on all connected front track segs
2296 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2297 ///have to check later
2298 if(frontsegIndex==-1) continue;
2300 fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]++] = frontsegIndex;
2302 /// fTrackParam[ibacktrackseg] = trackParam;
2305 HLTInfo("\tbacktrack : %d is connected with : %d front tracks\n",
2306 ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg]);
2307 HLTInfo("AliHLTMUONFullTracker::SelectFront()--End\n\n");
2311 }///backtrackSeg loop
2317 ///__________________________________________________________________________
2319 Bool_t AliHLTMUONFullTracker::KalmanChi2Test()
2322 /// Kalman Chi2 test for trach segments selection
2324 Cluster clus1,clus2;
2325 Int_t minIndex=0,maxIndex=0;
2326 Int_t minCh=0,maxCh=0;
2327 Int_t ifronttrackseg=0;
2328 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
2330 if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
2332 /// if(fBackTrackSeg[ibacktrackseg].fIndex[2]==-1 || fBackTrackSeg[ibacktrackseg].fIndex[3]==-1) continue;
2334 ifronttrackseg = fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]-1];
2337 HLTInfo("AliHLTMUONFullTracker::KalmanChi2Test()--Begin\n\n");
2338 HLTInfo("\tbacktrack : %d is connected with : %d front tracks, front track index %d\n",
2339 ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg],ifronttrackseg);
2341 maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
2342 maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
2344 minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
2345 minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
2348 clus2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
2349 clus2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
2350 clus2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
2351 clus2.fErrX2 = 0.020736;
2352 clus2.fErrY2 = 0.000100;
2354 clus1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
2355 clus1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
2356 clus1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
2357 clus1.fErrX2 = 0.020736;
2358 clus1.fErrY2 = 0.000100;
2361 AliMUONTrackParam trackParam;
2362 Double_t dZ = clus2.fZ - clus1.fZ;
2363 trackParam.SetNonBendingCoor(clus2.fX);
2364 trackParam.SetBendingCoor(clus2.fY);
2365 trackParam.SetZ(clus2.fZ);
2366 trackParam.SetNonBendingSlope((clus2.fX - clus1.fX) / dZ);
2367 trackParam.SetBendingSlope((clus2.fY - clus1.fY) / dZ);
2368 Double_t bendingImpact = clus1.fY - clus1.fZ * trackParam.GetBendingSlope();
2369 Double_t inverseBendingMomentum = 1.
2370 / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
2371 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
2374 HLTInfo("\t\tCh%d : (%f,%f,%f)\n",maxCh,clus2.fX,clus2.fY,clus2.fZ);
2375 HLTInfo("\t\tCh%d : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2376 HLTInfo("\t\tFor minCh : %d, GetBeMom : %lf\n",minCh,trackParam.GetInverseBendingMomentum());
2379 /// trackParam->SetClusterPtr(clus[8]);
2381 /// => Reset track parameter covariances at last clus (as if the other cluss did not exist)
2382 TMatrixD lastParamCov(5,5);
2383 lastParamCov.Zero();
2384 lastParamCov(0,0) = clus1.fErrX2;
2385 lastParamCov(1,1) = 100. * ( clus2.fErrX2 + clus1.fErrX2 ) / dZ / dZ;
2386 lastParamCov(2,2) = clus1.fErrY2;
2387 lastParamCov(3,3) = 100. * ( clus2.fErrY2 + clus1.fErrY2 ) / dZ / dZ;
2388 lastParamCov(4,4) = ((10.0*10.0 + (clus2.fZ * clus2.fZ * clus1.fErrY2 +
2389 clus1.fZ * clus1.fZ * clus2.fErrY2)
2390 / dZ / dZ) /bendingImpact / bendingImpact +
2391 0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum;
2392 trackParam.SetCovariances(lastParamCov);
2394 AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2395 /// AliMUONTrackExtrap::ExtrapToZCov(trackParam, AliMUONConstants::DefaultChamberZ(minCh),kTRUE);
2396 ///AliMUONTrackExtrap::ExtrapToZCov(&trackParam, clus1.fZ,kTRUE);
2398 LinearExtrapToZ(&trackParam, clus1.fZ);
2400 trackParam.SetTrackChi2(0.);
2402 Double_t chi2 = 0.0;
2404 chi2 = KalmanFilter(trackParam,&clus1);
2407 HLTInfo("\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\n",minCh,chi2,trackParam.GetInverseBendingMomentum());
2408 /// TMatrixD para2(trackParam->GetParameters());
2412 AliMUONTrackParam extrapTrackParamAtCluster2,minChi2Param;
2413 Double_t chi2OfCluster;
2415 Double_t minChi2=1000000.0;
2416 Int_t frontsegIndex = -1;
2417 extrapTrackParamAtCluster2 = trackParam ;
2418 minChi2Param = trackParam ;
2420 for( Int_t iconnected=0;iconnected<fNofConnectedfrontTrackSeg[ibacktrackseg];iconnected++){
2422 ifronttrackseg = fBackToFront[ibacktrackseg][iconnected];
2423 AliMUONTrackParam extrapTrackParam;
2424 trackParam = extrapTrackParamAtCluster2;
2426 minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2427 minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2428 clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2429 clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2430 clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2431 clus1.fErrX2 = 0.020736;
2432 clus1.fErrY2 = 0.000100;
2435 AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2436 trackParam.ResetPropagator();
2437 AliMUONTrackExtrap::ExtrapToZCov(&trackParam,clus1.fZ, kTRUE);
2439 tryonefast = TryOneClusterFast(trackParam, &clus1);
2440 ///if (!tryonefast) continue;
2442 /// try to add the current cluster accuratly
2443 chi2OfCluster = TryOneCluster(trackParam, &clus1 , extrapTrackParam,kTRUE);
2445 extrapTrackParam.SetExtrapParameters(extrapTrackParam.GetParameters());
2446 extrapTrackParam.SetExtrapCovariances(extrapTrackParam.GetCovariances());
2448 chi2 = KalmanFilter(extrapTrackParam,&clus1);
2451 minChi2Param = extrapTrackParam;
2452 frontsegIndex = ifronttrackseg;
2454 }///for loop on all connected front track segs
2456 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2457 ///have to check later
2458 if(frontsegIndex==-1) continue;
2460 fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]++] = frontsegIndex;
2461 ifronttrackseg = frontsegIndex;
2463 trackParam = minChi2Param ;
2466 HLTInfo("\t\t\tCh%d : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2467 HLTInfo("\t\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\t",minCh,chi2,trackParam.GetInverseBendingMomentum());
2468 HLTInfo(Form("Pt : %lf\n",TMath::Sqrt(trackParam.Px()*trackParam.Px() +
2469 trackParam.Py()*trackParam.Py())));
2473 minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2474 minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2475 clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2476 clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2477 clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2478 clus1.fErrX2 = 0.020736;
2479 clus1.fErrY2 = 0.000100;
2481 AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2482 trackParam.ResetPropagator();
2483 ///AliMUONTrackExtrap::ExtrapToZCov(&trackParam,clus1.fZ, kTRUE);
2484 LinearExtrapToZ(&trackParam, clus1.fZ);
2486 tryonefast = TryOneClusterFast(trackParam, &clus1);
2487 ///if (!tryonefast) continue;
2489 chi2OfCluster = TryOneCluster(trackParam, &clus1 , extrapTrackParamAtCluster2,kTRUE);
2491 extrapTrackParamAtCluster2.SetExtrapParameters(extrapTrackParamAtCluster2.GetParameters());
2492 extrapTrackParamAtCluster2.SetExtrapCovariances(extrapTrackParamAtCluster2.GetCovariances());
2494 chi2 = KalmanFilter(extrapTrackParamAtCluster2,&clus1);
2496 trackParam = extrapTrackParamAtCluster2;
2499 HLTInfo("\t\tCh%d : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2500 HLTInfo("\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\t",minCh,chi2,trackParam.GetInverseBendingMomentum());
2501 HLTInfo(Form("Pt : %lf\n",TMath::Sqrt(extrapTrackParamAtCluster2.Px()*extrapTrackParamAtCluster2.Px() +
2502 extrapTrackParamAtCluster2.Py()*extrapTrackParamAtCluster2.Py())));
2504 HLTInfo("AliHLTMUONFullTracker::KalmanChi2Test()--End\n\n");
2506 ///AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.);
2507 fTrackParam[ibacktrackseg] = trackParam;
2518 ///__________________________________________________________________________
2521 Double_t AliHLTMUONFullTracker::BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ)
2523 //// Returns the mean total momentum energy loss of muon with total momentum='pTotal'
2524 //// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
2525 Double_t muMass = 0.105658369; /// GeV
2526 Double_t eMass = 0.510998918e-3; /// GeV
2527 Double_t k = 0.307075e-3; /// GeV.g^-1.cm^2
2528 Double_t i = 9.5e-9; /// mean exitation energy per atomic Z (GeV)
2529 Double_t p2=pTotal*pTotal;
2530 Double_t beta2=p2/(p2 + muMass*muMass);
2532 Double_t w = k * rho * pathLength * atomicZ / atomicA / beta2;
2534 if (beta2/(1-beta2)>3.5*3.5)
2535 return w * (log(2.*eMass*3.5/(i*atomicZ)) + 0.5*log(beta2/(1-beta2)) - beta2);
2537 return w * (log(2.*eMass*beta2/(1-beta2)/(i*atomicZ)) - beta2);
2541 ///__________________________________________________________________________
2543 Double_t AliHLTMUONFullTracker::EnergyLossFluctuation2(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ)
2545 //// Returns the total momentum energy loss fluctuation of muon with total momentum='pTotal'
2546 //// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
2547 Double_t muMass = 0.105658369; /// GeV
2548 ///Double_t eMass = 0.510998918e-3; /// GeV
2549 Double_t k = 0.307075e-3; /// GeV.g^-1.cm^2
2550 Double_t p2=pTotal*pTotal;
2551 Double_t beta2=p2/(p2 + muMass*muMass);
2553 Double_t fwhm = 2. * k * rho * pathLength * atomicZ / atomicA / beta2; /// FWHM of the energy loss Landau distribution
2554 Double_t sigma2 = fwhm * fwhm / (8.*log(2.)); /// gaussian: fwmh = 2 * srqt(2*ln(2)) * sigma (i.e. fwmh = 2.35 * sigma)
2556 ///sigma2 = k * rho * pathLength * atomicZ / atomicA * eMass; /// sigma2 of the energy loss gaussian distribution
2562 ///__________________________________________________________________________
2564 void AliHLTMUONFullTracker::LinearExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd) const
2566 //// Track parameters (and their covariances if any) linearly extrapolated to the plane at "zEnd".
2567 //// On return, results from the extrapolation are updated in trackParam.
2569 //if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z
2571 Double_t dZ = zEnd - trackParam->GetZ();
2572 if (dZ == 0) return; // nothing to be done if same z
2573 trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + trackParam->GetNonBendingSlope() * dZ);
2574 trackParam->SetBendingCoor(trackParam->GetBendingCoor() + trackParam->GetBendingSlope() * dZ);
2575 trackParam->SetZ(zEnd);
2578 void AliHLTMUONFullTracker::CorrectELossEffectInAbsorber(AliMUONTrackParam* param, Double_t eLoss)
2580 // Energy loss coreection in front absorber.
2582 Double_t nonBendingSlope = param->GetNonBendingSlope();
2583 Double_t bendingSlope = param->GetBendingSlope();
2584 param->SetInverseBendingMomentum(param->GetCharge() / (param->P() + eLoss) *
2585 TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) /
2586 TMath::Sqrt(1.0 + bendingSlope*bendingSlope));
2589 ///__________________________________________________________________________
2591 void AliHLTMUONFullTracker::Cov2CovP(const TMatrixD ¶m, TMatrixD &cov)
2593 //// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, SlopeX, Y, SlopeY, q*PTot)
2594 //// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
2596 /// charge * total momentum
2597 Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
2598 TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
2600 /// Jacobian of the opposite transformation
2601 TMatrixD jacob(5,5);
2603 jacob(4,1) = qPTot * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
2604 jacob(4,3) = - qPTot * param(1,0) * param(1,0) * param(3,0) /
2605 (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
2606 jacob(4,4) = - qPTot / param(4,0);
2608 /// compute covariances in new coordinate system
2610 tmp.MultT(cov,jacob);
2611 cov.Mult(jacob,tmp);
2615 ///__________________________________________________________________________
2617 void AliHLTMUONFullTracker::CovP2Cov(const TMatrixD ¶m, TMatrixD &covP)
2619 //// change coordinate system: (X, SlopeX, Y, SlopeY, q*PTot) -> (X, SlopeX, Y, SlopeY, q/Pyz)
2620 //// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
2622 /// charge * total momentum
2623 Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
2624 TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
2626 /// Jacobian of the transformation
2627 TMatrixD jacob(5,5);
2629 jacob(4,1) = param(4,0) * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
2630 jacob(4,3) = - param(4,0) * param(1,0) * param(1,0) * param(3,0) /
2631 (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
2632 jacob(4,4) = - param(4,0) / qPTot;
2634 /// compute covariances in new coordinate system
2636 tmp.MultT(covP,jacob);
2637 covP.Mult(jacob,tmp);
2641 ///__________________________________________________________________________
2644 void AliHLTMUONFullTracker::CorrectMCSEffectInAbsorber(AliMUONTrackParam* param,
2645 Double_t xVtx, Double_t yVtx, Double_t zVtx,
2646 Double_t absZBeg, Double_t f1, Double_t f2)
2648 //// Correct parameters and corresponding covariances using Branson correction
2649 //// - input param are parameters and covariances at the end of absorber
2650 //// - output param are parameters and covariances at vertex
2651 //// Absorber correction parameters are supposed to be calculated at the current track z-position
2653 /// Position of the Branson plane (spectro. (z<0))
2654 Double_t zB = (f1>0.) ? absZBeg - f2/f1 : 0.;
2656 LinearExtrapToZ(param,zB);
2658 /// compute track parameters at vertex
2659 TMatrixD newParam(5,1);
2661 newParam(0,0) = xVtx;
2662 newParam(1,0) = (param->GetNonBendingCoor() - xVtx) / (zB - zVtx);
2663 newParam(2,0) = yVtx;
2664 newParam(3,0) = (param->GetBendingCoor() - yVtx) / (zB - zVtx);
2665 newParam(4,0) = param->GetCharge() / param->P() *
2666 TMath::Sqrt(1.0 + newParam(1,0)*newParam(1,0) + newParam(3,0)*newParam(3,0)) /
2667 TMath::Sqrt(1.0 + newParam(3,0)*newParam(3,0));
2669 TMatrixD paramCovP(5,5);
2671 /// Get covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
2672 paramCovP.Use(param->GetCovariances());
2674 Cov2CovP(param->GetParameters(),paramCovP);
2676 /// Get the covariance matrix in the (XVtx, X, YVtx, Y, q*PTot) coordinate system
2677 /// TMatrixD paramCovVtx(5,5);
2678 Double_t element44 = paramCovP(4,4);
2680 paramCovP(4,4) = element44;
2682 CovP2Cov(newParam,paramCovP);
2684 /// Set parameters and covariances at vertex
2685 param->SetParameters(newParam);
2687 param->SetCovariances(paramCovP);
2690 ///__________________________________________________________________________
2692 Bool_t AliHLTMUONFullTracker::ExtrapolateToOrigin(Bool_t extrap)
2694 /// Extrapolation to origin through absorber
2696 Int_t minIndex=0,maxIndex=0;
2697 Int_t minCh=0,maxCh=0;
2698 Int_t ifronttrackseg = -1;
2699 AliMUONTrackParam trackP;
2700 Double_t slopeB, slopeNB;
2701 AliHLTMUONRecHitStruct p1,p2,pSeg1,pSeg2;
2702 Double_t pyz = -1.0;
2703 TVector3 v1,v2,v3,v4;
2704 Double_t eLoss1,eLoss2,eLoss3;
2706 Double_t zE,zB,dzE,dzB;
2708 Double_t f0Sum,f1Sum,f2Sum;
2709 Double_t fXVertex=0.0,fYVertex=0.0,fZVertex=0.0;
2711 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
2713 if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
2715 maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
2716 maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
2718 minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
2719 minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
2721 p2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
2722 p2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
2723 p2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
2725 p1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
2726 p1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
2727 p1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
2730 Sub(&p2,&p1,&pSeg1);
2732 ifronttrackseg = fBackToFront[ibacktrackseg][0];
2734 maxIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2735 maxCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2737 minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2738 minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2740 p2.fX = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fX ;
2741 p2.fY = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fY ;
2742 p2.fZ = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fZ ;
2744 p1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2745 p1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2746 p1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2748 Sub(&p2,&p1,&pSeg2);
2750 pyz = -(3.0*0.3/sin(Angle(&pSeg1,&pSeg2)));/// * sqrt(x3*x3 + y3*y3)/z3 ;
2753 slopeNB = (p2.fX - p1.fX)/(p2.fZ - p1.fZ);
2754 slopeB = (p2.fY - p1.fY)/(p2.fZ - p1.fZ);
2758 trackP.SetNonBendingCoor(p1.fX);
2759 trackP.SetNonBendingSlope(slopeNB);
2760 trackP.SetBendingCoor(p1.fY);
2761 trackP.SetBendingSlope(slopeB);
2762 trackP.SetInverseBendingMomentum(1.0/pyz) ;
2767 trackP = fTrackParam[ibacktrackseg] ;
2769 LinearExtrapToZ(&trackP,fgkAbsoedge[3]);
2770 v4.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
2771 LinearExtrapToZ(&trackP,fgkAbsoedge[2]);
2772 v3.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
2773 LinearExtrapToZ(&trackP,fgkAbsoedge[1]);
2774 v2.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
2775 LinearExtrapToZ(&trackP,fgkAbsoedge[0]);
2776 v1.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
2778 eLoss1 = BetheBloch(trackP.P(), (v4-v3).Mag(), fgkRho[2], fgkAtomicA[2], fgkAtomicZ[2]);
2779 eLoss2 = BetheBloch(trackP.P(), (v3-v2).Mag(), fgkRho[1], fgkAtomicA[1], fgkAtomicZ[1]);
2780 eLoss3 = BetheBloch(trackP.P(), (v2-v1).Mag(), fgkRho[0], fgkAtomicA[0], fgkAtomicZ[0]);
2782 /// sigmaELoss1 = EnergyLossFluctuation2(trackP.P(), (v4-v3).Mag(), rho[2], atomicA[2], atomicZ[2]);
2783 /// sigmaELoss2 = EnergyLossFluctuation2(trackP.P(), (v3-v2).Mag(), rho[1], atomicA[1], atomicZ[1]);
2784 /// sigmaELoss3 = EnergyLossFluctuation2(trackP.P(), (v2-v1).Mag(), rho[0], atomicA[0], atomicZ[0]);
2786 /// eDiff = totELoss-(eLoss1+eLoss2+eLoss3);
2787 /// sigmaELossDiff = sigmaTotELoss ;///- (sigmaELoss1+sigmaELoss2+sigmaELoss3);
2789 /// CorrectELossEffectInAbsorber(&trackP, 0.5*(eLoss1+eLoss2+eLoss3), 0.5*(sigmaELoss1+sigmaELoss2+sigmaELoss3));
2792 ///CorrectELossEffectInAbsorber(&trackP, totELoss,sigmaTotELoss);
2794 CorrectELossEffectInAbsorber(&trackP, 0.7*(eLoss1+eLoss2+eLoss3));
2796 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2797 f0Sum = 0.0; f1Sum = 0.0; f2Sum = 0.0;
2799 b = (v4.Z()-v1.Z())/((v4-v1).Mag());
2802 zE = b*((v2-v1).Mag()) + zB;
2806 f0 = ((v2-v1).Mag())/fgkRadLen[0];
2807 f1 = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[0] /2.;
2808 f2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[0] / 3.;
2815 zE = b*((v3-v2).Mag()) + zB;
2819 f0 = ((v3-v2).Mag())/fgkRadLen[1];
2820 f1 = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[1] /2.;
2821 f2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[1] / 3.;
2828 zE = b*((v4-v3).Mag()) + zB;
2832 f0 = ((v4-v3).Mag())/fgkRadLen[2];
2833 f1 = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[2] /2.;
2834 f2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[2] / 3.;
2841 ///AddMCSEffectInAbsorber(&trackP,(v4-v1).Mag(),f0Sum,f1Sum,f2Sum);
2843 CorrectMCSEffectInAbsorber(&trackP,fXVertex,fYVertex, fZVertex,AliMUONConstants::AbsZBeg(),f1Sum,f2Sum);
2844 CorrectELossEffectInAbsorber(&trackP, 0.5*(eLoss1+eLoss2+eLoss3));
2847 ///AliMUONTrackExtrap::ExtrapToVertex(&trackP, 0., 0., 0., 0., 0.);
2848 fTrackParam[ibacktrackseg] = trackP;
2850 }///backtrackseg for loop
2855 ///__________________________________________________________________________
2857 //FIXME: remove the InitGRP method. This is handled by the HLT framework and should not be redone in the component.
2858 Bool_t AliHLTMUONFullTracker::InitGRP()
2860 ///GRP handling needed for standalone testing and debug purpose
2862 AliGRPObject* fGRPData = 0x0; /// Data from the GRP/GRP/Data CDB folder
2864 ///------------------------------------
2865 /// Initialization of the GRP entry
2866 ///------------------------------------
2867 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
2871 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); /// old GRP entry
2874 HLTInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject\n");
2876 fGRPData = new AliGRPObject();
2877 fGRPData->ReadValuesFromMap(m);
2881 HLTInfo("Found an AliGRPObject in GRP/GRP/Data, reading it\n");
2882 fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject()); /// new GRP entry
2886 /// FIX ME: The unloading of GRP entry is temporarily disabled
2887 /// because ZDC and VZERO are using it in order to initialize
2888 /// their reconstructor objects. In the future one has to think
2889 /// of propagating AliRunInfo to the reconstructors.
2890 /// AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
2894 HLTInfo("No GRP entry found in OCDB!\n");
2898 TString lhcState = fGRPData->GetLHCState();
2899 if (lhcState==AliGRPObject::GetInvalidString()) {
2900 HLTInfo("GRP/GRP/Data entry: missing value for the LHC state ! Using UNKNOWN\n");
2901 lhcState = "UNKNOWN";
2904 TString beamType = fGRPData->GetBeamType();
2905 if (beamType==AliGRPObject::GetInvalidString()) {
2906 HLTInfo("GRP/GRP/Data entry: missing value for the beam type ! Using UNKNOWN\n");
2907 beamType = "UNKNOWN";
2910 Float_t beamEnergy = fGRPData->GetBeamEnergy();
2911 if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
2912 HLTInfo("GRP/GRP/Data entry: missing value for the beam energy ! Using 0\n");
2915 /// LHC: "multiply by 120 to get the energy in MeV"
2916 beamEnergy *= 0.120;
2918 TString runType = fGRPData->GetRunType();
2919 if (runType==AliGRPObject::GetInvalidString()) {
2920 HLTInfo("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN\n");
2921 runType = "UNKNOWN";
2924 Int_t activeDetectors = fGRPData->GetDetectorMask();
2925 if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
2926 HLTInfo("GRP/GRP/Data entry: missing value for the detector mask ! Using 1074790399\n");
2927 activeDetectors = 1074790399;
2930 AliRunInfo *fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
2934 // Dealing with the magnetic field map
2935 if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
2936 if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
2937 HLTInfo("ExpertMode!!! GRP information will be ignored !\n");
2938 HLTInfo("ExpertMode!!! Running with the externally locked B field !\n");
2941 HLTInfo("Destroying existing B field instance!\n");
2942 delete TGeoGlobalMagField::Instance();
2945 if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
2946 /// Construct the field map out of the information retrieved from GRP.
2949 Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
2950 if (l3Current == AliGRPObject::GetInvalidFloat()) {
2951 HLTInfo("GRP/GRP/Data entry: missing value for the L3 current !\n");
2955 Char_t l3Polarity = fGRPData->GetL3Polarity();
2956 if (l3Polarity == AliGRPObject::GetInvalidChar()) {
2957 HLTInfo("GRP/GRP/Data entry: missing value for the L3 polarity !\n");
2962 Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
2963 if (diCurrent == AliGRPObject::GetInvalidFloat()) {
2964 HLTInfo("GRP/GRP/Data entry: missing value for the dipole current !\n");
2968 Char_t diPolarity = fGRPData->GetDipolePolarity();
2969 if (diPolarity == AliGRPObject::GetInvalidChar()) {
2970 HLTInfo("GRP/GRP/Data entry: missing value for the dipole polarity !\n");
2974 /// read special bits for the polarity convention and map type
2975 Int_t polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
2976 Bool_t uniformB = fGRPData->IsUniformBMap();
2979 AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1),
2980 TMath::Abs(diCurrent) * (diPolarity ? -1:1),
2981 polConvention,uniformB,beamEnergy, beamType.Data());
2983 TGeoGlobalMagField::Instance()->SetField( fld );
2984 TGeoGlobalMagField::Instance()->Lock();
2985 HLTInfo("Running with the B field constructed out of GRP !\n");
2987 else HLTInfo("Failed to create a B field map !\n");
2989 else HLTInfo("B field is neither set nor constructed from GRP ! Exitig...\n");