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 /////////////////////////////////////////////////////////////////////////
36 #include "TGeoGlobalMagField.h"
40 #include "AliCDBManager.h"
41 #include "AliCDBStorage.h"
42 #include "AliGeomManager.h"
45 #include "AliMUONTrackExtrap.h"
46 #include "AliMUONTrackParam.h"
47 #include "AliMUONTrackExtrap.h"
48 #include "AliMUONConstants.h"
49 #include "AliMUONGeometryTransformer.h"
50 #include "AliMUONTrackParam.h"
52 //MUON mappinf includes
53 #include "AliMpDEIterator.h"
56 #include "AliHLTLogging.h"
59 #include "AliHLTMUONConstants.h"
60 #include "AliHLTMUONDataTypes.h"
61 #include "AliHLTMUONUtils.h"
63 //HLT/MUON/OnlineAnalysis includes
64 #include "AliHLTMUONFullTracker.h"
72 class AliMpSegmentation;
78 //#define PRINT_FULL 1
81 #define PRINT_POINTS 1
84 #define PRINT_KALMAN 1
85 #define PRINT_SELECT 1
86 #define PRINT_OUTPUT 1
87 #define PRINT_TRACKSEG 1
88 ///#define PRINT_DETAIL_KALMAN 1
91 class AliHLTMUONFullTracker;
93 const Float_t AliHLTMUONFullTracker::fgkTrackDetCoordinate[3] = {
94 155.179+20.0, 166.234+20.0,
95 (AliMUONConstants::DefaultChamberZ(4)+ AliMUONConstants::DefaultChamberZ(5))/2.0
98 const Double_t AliHLTMUONFullTracker::fgkAbsoedge[4] = {92.0, 317.0,443.0,499.0};
99 const Double_t AliHLTMUONFullTracker::fgkRadLen[3] = {13.875635, 11.273801, 1.765947};
100 const Double_t AliHLTMUONFullTracker::fgkRho[3] = {1.750000, 2.350000, 7.880000};
101 const Double_t AliHLTMUONFullTracker::fgkAtomicZ[3] = {6.000000, 11.098486,25.780000};
102 const Double_t AliHLTMUONFullTracker::fgkAtomicA[3] = {12.010000, 22.334156,55.299670 };
105 const Int_t AliHLTMUONFullTracker::fgkMaxNofCellsPerCh = 1500;
106 const Int_t AliHLTMUONFullTracker::fgkMaxNofPointsPerCh = 600;
107 const Int_t AliHLTMUONFullTracker::fgkMaxNofCh = 11 ; /// 10tracking + 1 trigger
108 const Int_t AliHLTMUONFullTracker::fgkMaxNofTracks = 200;
109 const Int_t AliHLTMUONFullTracker::fgkMaxNofConnectedTracks = 20;
110 const Int_t AliHLTMUONFullTracker::fgkMaxNofTriggers = 20;
114 AliHLTMUONFullTracker::AliHLTMUONFullTracker() :
116 fChamberGeometryTransformer(0x0),
123 fInclinationBack(0x0),
124 fNofConnectedfrontTrackSeg(0x0),
132 fNoffrontTrackSeg(0),
136 fFastTracking(kFALSE),
138 fNofTriggerInputs(0),
139 fNofTrackerInputs(0),
143 /// Constructor of the class
145 /// fChPoint = (AliHLTMUONRecHitStruct***)(malloc(10 * sizeof(AliHLTMUONRecHitStruct)));
146 /// for( Int_t ich=0;ich<10;ich++)
147 /// fChPoint[ich] = (AliHLTMUONRecHitStruct**)(malloc(300 * sizeof(AliHLTMUONRecHitStruct)));
149 fNofCells[0] = fNofCells[1] = 0;
152 fChPoint = new AliHLTMUONRecHitStruct**[fgkMaxNofCh-1];
153 }catch (const std::bad_alloc&){
154 HLTError("Dynamic memory allocation failed in constructor : fChPoint");
158 for( Int_t ich=0;ich<(fgkMaxNofCh-1);ich++)
160 fChPoint[ich] = new AliHLTMUONRecHitStruct*[fgkMaxNofPointsPerCh];
161 }catch (const std::bad_alloc&){
162 HLTError("Dynamic memory allocation failed in constructor : fChPoint[%d]",ich);
167 fChPoint11 = new AliHLTMUONTriggerRecordStruct*[fgkMaxNofPointsPerCh];
168 }catch (const std::bad_alloc&){
169 HLTError("Dynamic memory allocation failed in constructor : fChPoint11");
174 fBackTrackSeg = new TrackSeg[fgkMaxNofTracks];
175 }catch (const std::bad_alloc&){
176 HLTError("Dynamic memory allocation failed in constructor : fBackTrackSeg");
181 fFrontTrackSeg = new TrackSeg[fgkMaxNofTracks];
182 }catch (const std::bad_alloc&){
183 HLTError("Dynamic memory allocation failed in constructor : fFrontTrackSeg");
188 fExtrapSt3X = new float[fgkMaxNofTracks];
189 }catch (const std::bad_alloc&){
190 HLTError("Dynamic memory allocation failed in constructor : fExtrapSt3X");
195 fExtrapSt3Y = new float[fgkMaxNofTracks];
196 }catch (const std::bad_alloc&){
197 HLTError("Dynamic memory allocation failed in constructor : fExtrapSt3Y");
202 fInclinationBack = new float[fgkMaxNofTracks];
203 }catch (const std::bad_alloc&){
204 HLTError("Dynamic memory allocation failed in constructor : fInclinationBack");
209 fNofConnectedfrontTrackSeg = new int[fgkMaxNofTracks];
210 }catch (const std::bad_alloc&){
211 HLTError("Dynamic memory allocation failed in constructor : fNofConnectedfrontTrackSeg");
216 fBackToFront = new int*[fgkMaxNofTracks];
217 }catch (const std::bad_alloc&){
218 HLTError("Dynamic memory allocation failed in constructor : fBackToFront");
222 for( Int_t itr=0;itr<fgkMaxNofTracks;itr++)
224 fBackToFront[itr] = new int[fgkMaxNofConnectedTracks];
225 }catch (const std::bad_alloc&){
226 HLTError("Dynamic memory allocation failed in constructor : fBackToFront[%d]",itr);
232 fNofPoints = new int[fgkMaxNofCh];
233 }catch (const std::bad_alloc&){
234 HLTError("Dynamic memory allocation failed in constructor : fNofPoints");
239 fTrackParam = new AliMUONTrackParam[fgkMaxNofTracks];
240 }catch (const std::bad_alloc&){
241 HLTError("Dynamic memory allocation failed in constructor : fTrackParam");
246 fHalfTrack = new HalfTrack[fgkMaxNofTracks];
247 }catch (const std::bad_alloc&){
248 HLTError("Dynamic memory allocation failed in constructor : fHalfTrack");
256 ///__________________________________________________________________________
258 AliHLTMUONFullTracker::~AliHLTMUONFullTracker()
260 /// Destructor of the class
262 ///delete fChamberGeometryTransformer;
266 delete []fBackTrackSeg;
267 delete []fFrontTrackSeg;
268 delete []fExtrapSt3X;
269 delete []fExtrapSt3Y;
270 delete []fInclinationBack;
271 delete []fNofConnectedfrontTrackSeg;
272 delete []fBackToFront;
274 delete []fTrackParam;
277 fChamberGeometryTransformer->Delete();
279 fDetElemList.clear();
283 ///__________________________________________________________________________
285 Bool_t AliHLTMUONFullTracker::Clear(){
287 /// Clear method to be called after each event
289 #ifdef PRINT_TRACKSEG
290 for(int iFrontTrackSeg=0;iFrontTrackSeg<fNoffrontTrackSeg;iFrontTrackSeg++)
291 for(int ipoint=0;ipoint<4;ipoint++)
292 if(fFrontTrackSeg[iFrontTrackSeg].fIndex[ipoint]!=-1)
293 HLTImportant("FrontTrackSeg : %d\t%f\t%f\t%f\n",iFrontTrackSeg,
294 fChPoint[ipoint][fFrontTrackSeg[iFrontTrackSeg].fIndex[ipoint]]->fX,
295 fChPoint[ipoint][fFrontTrackSeg[iFrontTrackSeg].fIndex[ipoint]]->fY,
296 fChPoint[ipoint][fFrontTrackSeg[iFrontTrackSeg].fIndex[ipoint]]->fZ);
297 HLTImportant("\n\n");
298 for(int iBackTrackSeg=0;iBackTrackSeg<fNofbackTrackSeg;iBackTrackSeg++)
299 for(int ipoint=0;ipoint<4;ipoint++)
300 if(fBackTrackSeg[iBackTrackSeg].fIndex[ipoint]!=-1)
301 HLTImportant("BackTrackSeg : %d\t%f\t%f\t%f, nofFront : %d\n",iBackTrackSeg,
302 fChPoint[ipoint+6][fBackTrackSeg[iBackTrackSeg].fIndex[ipoint]]->fX,
303 fChPoint[ipoint+6][fBackTrackSeg[iBackTrackSeg].fIndex[ipoint]]->fY,
304 fChPoint[ipoint+6][fBackTrackSeg[iBackTrackSeg].fIndex[ipoint]]->fZ,fNofConnectedfrontTrackSeg[iBackTrackSeg]);
309 fNofTriggerInputs = 0;
310 fNofTrackerInputs = 0;
312 for( Int_t ich=0;ich<fgkMaxNofCh;ich++)
315 fNofbackTrackSeg = 0;
316 fNoffrontTrackSeg = 0;
323 ///__________________________________________________________________________
325 void AliHLTMUONFullTracker::Print()
327 /// Print anything mainly for debugging the codes, will be removed later
328 HLTInfo("\nPrint This--->\n");
331 ///__________________________________________________________________________
333 Bool_t AliHLTMUONFullTracker::Init()
335 /// Initilation to be called once, later can be used to set/load the CDB path/entries
336 HLTInfo("path : %s, run number : %d",
337 (AliCDBManager::Instance())->GetDefaultStorage()->GetURI().Data(),(AliCDBManager::Instance())->GetRun());
338 if (AliGeomManager::GetGeometry() == NULL){
339 AliGeomManager::LoadGeometry();
340 AliGeomManager::ApplyAlignObjsFromCDB("GRP MUON");
344 x[0] = 0.0 ; x[1] = 0.0 ; x[2] = -950.0;
346 TGeoGlobalMagField::Instance()->Field(x,b);
347 //The following condition is based on the fact that at the middle of the dipole the field cannot be zero
348 if(TMath::AreEqualAbs(b[0],0.0,1.0e-5) and TMath::AreEqualAbs(b[1],0.0,1.0e-5) and TMath::AreEqualAbs(b[2],0.0,1.0e-5)){
349 HLTWarning("Magnetic field is not set by GRP");
350 if(not TGeoGlobalMagField::Instance()->IsLocked()){
351 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1, AliMagF::k5kG));
352 TGeoGlobalMagField::Instance()->Field(x,b);
355 HLTWarning("Magnetic field is not set and cannot be set since it is locked");
362 HLTImportant("At (X,Y,Z) : (%6.2lf,%6.2lf,%6.2lf) Field (Bx,By,Bz) is (%6.2lf,%6.2lf,%6.2lf)",
363 x[0],x[1],x[2],b[0],b[1],b[2]);
365 AliMUONTrackExtrap::SetField();
366 /// see header file for class documentation
367 fChamberGeometryTransformer = new AliMUONGeometryTransformer();
368 if(! fChamberGeometryTransformer->LoadGeometryData()){
369 HLTError("Failed to Load Geomerty Data ");
372 fDetElemList.clear();
373 for(int ich=0;ich<AliMUONConstants::NCh();ich++){
375 for ( it.First(ich); ! it.IsDone(); it.Next() )
377 Int_t detElemId = it.CurrentDEId();
378 fDetElemList[detElemId] = detElemId;
385 ///__________________________________________________________________________
387 Bool_t AliHLTMUONFullTracker::SetInput(AliHLTInt32_t /*ddl*/, const AliHLTMUONRecHitStruct *data, AliHLTInt32_t size)
389 /// Set the input for rechit points
390 if(int(size)>=fgkMaxNofPointsPerCh/2)
393 AliHLTUInt16_t detElemID;
394 AliHLTUInt8_t chamber;
397 HLTImportant("Received from DDL : %d, nofHits : %d, data : %p",ddl,size,data);
399 for( Int_t ipoint=0;ipoint<int(size);ipoint++){
401 HLTError("Null Data pointer from HitRec");
406 AliHLTMUONUtils::UnpackRecHitFlags(data->fFlags,chamber,detElemID);
408 if((not fDetElemList[detElemID]) or (chamber>=AliMUONConstants::NTrackingCh())){
409 HLTDebug("Invalid tracking detelem : %d or chamber : %d",detElemID,chamber);
414 HLTImportant("ch : %02d, detelem : %04d, (X,Y,Z) : (%8.3f,%8.3f,%8.3f)",
415 chamber,detElemID,data->fX,data->fY,data->fZ);
418 fChPoint[detElemID/100-1][fNofPoints[detElemID/100-1]++] = (AliHLTMUONRecHitStruct *)data;
426 ///__________________________________________________________________________
428 Bool_t AliHLTMUONFullTracker::SetInput(AliHLTInt32_t /*ddl*/, const AliHLTMUONTriggerRecordStruct *data, AliHLTInt32_t size)
430 /// Set the input for trigrecs
432 // if(int(size)>=fgkMaxNofPointsPerCh/2)
435 if(int(size)>fgkMaxNofTriggers){
436 HLTWarning("More triggers (%d) found than max limit : %d (not possible physics events)",int(size),fgkMaxNofTriggers);
441 AliHLTUInt16_t detElemID;
442 AliHLTUInt8_t chamber;
444 for( Int_t ipoint=0;ipoint<int(size);ipoint++){
446 HLTError("Null Data pointer from TrigRec");
450 fChPoint11[fNofPoints[10]++] = (AliHLTMUONTriggerRecordStruct *)data;
451 for( Int_t ich=0;ich<4;ich++){
452 AliHLTMUONUtils::UnpackRecHitFlags((data->fHit[ich]).fFlags,chamber,detElemID);
453 if((not fDetElemList[detElemID]) or (chamber<AliMUONConstants::NTrackingCh()) or (chamber > AliMUONConstants::NCh()) ){
454 HLTDebug("Invalid trigger detelem : %d or chamber : %d",detElemID,chamber);
458 HLTImportant("size : %d, itrig : %04d, ch : %02d, detelem : %04d, (X,Y,Z) : (%8.3f,%8.3f,%8.3f)\n",
459 size,ipoint,chamber,detElemID,(data->fHit[ich]).fX,(data->fHit[ich]).fY,(data->fHit[ich]).fZ);
470 ///__________________________________________________________________________
472 Bool_t AliHLTMUONFullTracker::CheckInput(AliHLTEventID_t iEvent)
474 /// Cross Check all the inputs before the starting of the tracking
476 bool resultOk = true;
478 //if more than expected inputs or no input from one detector, do not do anything
479 if((fNofInputs > 22) or (fNofTrackerInputs > 20) or (fNofTriggerInputs > 2)
480 or (fNofTrackerInputs == 0) or (fNofTriggerInputs == 0)){
485 }else{ // make double sure that no space point pointer is null
487 //HLTImportant("Found fNofInputs : %d less than 22",fNofInputs);
489 //tracker chamber test
490 for(int ich=0;ich<fgkMaxNofCh-1;ich++){
491 for(int ipt=0;ipt<fNofPoints[ich];ipt++){
493 if((not fChPoint[ich][ipt]) or (TMath::AreEqualAbs(fChPoint[ich][ipt]->fX,0.0,1.0e-5) and
494 TMath::AreEqualAbs(fChPoint[ich][ipt]->fY,0.0,1.0e-5) and
495 TMath::AreEqualAbs(fChPoint[ich][ipt]->fZ,0.0,1.0e-5))){
497 HLTError("iEvent : 0x%x, fNofTrackerInputs : %d, Nof tracker point for chamber %d, is not equal to nof valid tracker pointer",
498 iEvent,fNofTrackerInputs,ich);
502 }// // if resultOk not already false
504 //trigger chamber test
505 if(fNofPoints[10] == 0){
510 for(int ipt=0;ipt<fNofPoints[10];ipt++){
511 if(not fChPoint11[ipt]){
513 HLTError("iEvent : 0x%x, fNofTriggerInputs : %d, Nof trigger points, is not equal to nof valid tracker pointer",
514 iEvent,fNofTriggerInputs);
518 }// for loop over points
520 }// if less than expected blocks found
527 ///__________________________________________________________________________
529 Bool_t AliHLTMUONFullTracker::Run( AliHLTEventID_t iEvent,AliHLTMUONTrackStruct *data, AliHLTUInt32_t& size)
531 /// Main Run call of the class
533 bool resultOk = true;
535 HLTDebug("Processing iEvent : 0x%X, nof triggers : %d, fNofInputs : %d",iEvent,fNofPoints[10],fNofInputs);
537 // for(int ich=0;ich<fgkMaxNofCh;ich++)
538 // HLTDebug("\tNof hits in ich [%d] : %d\t",ich,fNofPoints[ich]);
540 resultOk = CheckInput(iEvent);
542 if((not fIsMagfield) and resultOk)
546 resultOk = SlatTrackSeg();
548 HLTDebug("Error happened in tracking through slat chambers, this event will be skipped");
551 HLTDebug("Finishing SlatTrackSeg");
554 resultOk = PrelimMomCalc();
556 HLTDebug("Error happened in calculating preliminary momentum, this event will be skipped");
559 HLTDebug("Finishing PrelimMomCalc");
562 resultOk = QuadTrackSeg();
564 HLTDebug("Error happened in tracking through quadrant chambers, this event will be skipped");
567 HLTDebug("Finishing QuadTrackSeg");
572 resultOk = SelectFront();
574 resultOk = KalmanChi2Test();
577 HLTDebug("Error happened in tracking through in Kalman Chi2 checking, this event will be skipped");
580 HLTDebug("Finishing KalmanChi2Test");
583 resultOk = ExtrapolateToOrigin();
585 HLTDebug("Error happened in tracking extrapolation, this event will be skipped");
588 HLTDebug("Finishing ExtrapolateToOrigin");
591 resultOk = FillOutData(data,size);
593 HLTDebug("Error happened in filling the output data, this event will be skipped");
596 HLTDebug("Finishing FillOutData");
601 HLTDebug("iEvent: %llu, has tracks : %d, triggers : %d, nof slat tracks : %d, quad tracks : %d, connected : %d\n",
602 iEvent,size,fNofPoints[10],fNofbackTrackSeg,fNoffrontTrackSeg,fNofConnected);
609 ///__________________________________________________________________________
611 void AliHLTMUONFullTracker::Sub(const AliHLTMUONRecHitStruct *v1, const AliHLTMUONRecHitStruct *v2, AliHLTMUONRecHitStruct *v3) const
613 /// Subtraction of position co-odinate of two space points
615 v3->fX = v1->fX - v2->fX;
616 v3->fY = v1->fY - v2->fY;
617 v3->fZ = v1->fZ - v2->fZ;
620 ///__________________________________________________________________________
622 Double_t AliHLTMUONFullTracker::Angle(const AliHLTMUONRecHitStruct *v1, const AliHLTMUONRecHitStruct *v2)
624 ///Angle of a straight line formed using v1 and v2
626 Double_t angle = 0.0;
628 Float_t ptot2 = ((v1->fX * v1->fX) + (v1->fY * v1->fY) + (v1->fZ * v1->fZ))*
629 ((v2->fX * v2->fX) + (v2->fY * v2->fY) + (v2->fZ * v2->fZ));
633 Float_t arg = ((v1->fX * v2->fX) + (v1->fY * v2->fY) + (v1->fZ * v2->fZ))/sqrt(ptot2);
634 if(arg > 1.0) arg = 1.0;
635 if(arg < -1.0) arg = -1.0;
636 angle = TMath::ACos(arg);
637 if(TMath::AreEqualAbs(angle,0.0,1.0e-5))
646 ///__________________________________________________________________________
648 Bool_t AliHLTMUONFullTracker::FillOutData(AliHLTMUONTrackStruct *track, AliHLTUInt32_t& size)
650 ///Fill the output data pointers
652 size = (AliHLTUInt32_t(fNofbackTrackSeg)<size) ? AliHLTUInt32_t(fNofbackTrackSeg) : size;
655 for( Int_t ibackTrackSeg=0;ibackTrackSeg<int(size);ibackTrackSeg++){
657 if(fNofConnectedfrontTrackSeg[ibackTrackSeg]>0){
660 if(not TMath::Finite(fTrackParam[ibackTrackSeg].Px())
661 || not TMath::Finite(fTrackParam[ibackTrackSeg].Py())
662 || not TMath::Finite(fTrackParam[ibackTrackSeg].Pz())) continue;
665 HLTImportant("\nsize : %d, itrack : %04d, sign : %2d, Pt : %8.3f, (Px,Py,Pz) : (%8.3f,%8.3f,%8.3f)\n",
666 size,ibackTrackSeg,Int_t(TMath::Sign(1.,fTrackParam[ibackTrackSeg].GetInverseBendingMomentum())),
667 TMath::Sqrt(fTrackParam[ibackTrackSeg].Px()*fTrackParam[ibackTrackSeg].Px() +
668 fTrackParam[ibackTrackSeg].Py()*fTrackParam[ibackTrackSeg].Py()),
669 fTrackParam[ibackTrackSeg].Px(),fTrackParam[ibackTrackSeg].Py(),
670 fTrackParam[ibackTrackSeg].Pz());
674 // Bits 8 and 9 must be kept zero to prevent the track ID from conflicting with other tracker components.
675 track->fId = (ibackTrackSeg << 10) | (fBackTrackSeg[ibackTrackSeg].fTrigRec & 0xFF);
676 track->fTrigRec = fBackTrackSeg[ibackTrackSeg].fTrigRec;
677 track->fPx = fTrackParam[ibackTrackSeg].Px();
678 track->fPy = fTrackParam[ibackTrackSeg].Py();
679 track->fPz = fTrackParam[ibackTrackSeg].Pz();
683 track->fInverseBendingMomentum = fTrackParam[ibackTrackSeg].GetInverseBendingMomentum();
684 track->fThetaY = TMath::Tan(fTrackParam[ibackTrackSeg].GetBendingSlope());
685 track->fThetaX = TMath::Tan(fTrackParam[ibackTrackSeg].GetNonBendingSlope());
687 track->fZ = fTrackParam[ibackTrackSeg].GetZ();
688 track->fY = fTrackParam[ibackTrackSeg].GetBendingCoor();
689 track->fX = fTrackParam[ibackTrackSeg].GetNonBendingCoor();
692 for( Int_t ipoint=15;ipoint>=0;ipoint--){
693 track->fHit[ipoint] = AliHLTMUONConstants::NilRecHitStruct();
694 hitset[ipoint] = false;
695 if(ipoint >= 6 and ipoint <= 9 and fBackTrackSeg[ibackTrackSeg].fIndex[ipoint-6]!=-1 ){
696 track->fHit[ipoint] = *(fChPoint[ipoint][fBackTrackSeg[ibackTrackSeg].fIndex[ipoint-6]]);
697 hitset[ipoint] = true;
699 AliHLTUInt8_t chamber; AliHLTUInt16_t detElemID;
700 AliHLTMUONUtils::UnpackRecHitFlags((track->fHit[ipoint]).fFlags,chamber,detElemID);
701 HLTImportant("(X,Y,Z) : (%lf,%lf,%lf)",(track->fHit[ipoint]).fX,(track->fHit[ipoint]).fY,(track->fHit[ipoint]).fZ);
703 }else if(ipoint <= 3 and fFrontTrackSeg[fBackToFront[ibackTrackSeg][0]].fIndex[ipoint]!=-1 ){
704 track->fHit[ipoint] = *(fChPoint[ipoint][fFrontTrackSeg[fBackToFront[ibackTrackSeg][0]].fIndex[ipoint]]);
705 hitset[ipoint] = true;
707 AliHLTUInt8_t chamber; AliHLTUInt16_t detElemID;
708 AliHLTMUONUtils::UnpackRecHitFlags((track->fHit[ipoint]).fFlags,chamber,detElemID);
709 HLTImportant("(X,Y,Z) : (%lf,%lf,%lf)",(track->fHit[ipoint]).fX,(track->fHit[ipoint]).fY,(track->fHit[ipoint]).fZ);
713 AliHLTMUONParticleSign sign = AliHLTMUONParticleSign(Int_t(TMath::Sign(1.,fTrackParam[ibackTrackSeg].GetInverseBendingMomentum())));
714 track->fFlags = AliHLTMUONUtils::PackTrackFlags(sign,hitset);
722 if(not TMath::Finite(fHalfTrack[ibackTrackSeg].fPx)
723 || not TMath::Finite(fHalfTrack[ibackTrackSeg].fPy)
724 || not TMath::Finite(fHalfTrack[ibackTrackSeg].fPz)) continue;
727 HLTImportant("\nhalftrack : size : %d, itrack : %04d, sign : %2d, Pt : %8.3f, (Px,Py,Pz) : (%8.3f,%8.3f,%8.3f)\n",
728 size,ibackTrackSeg,Int_t(TMath::Sign(1.,fTrackParam[ibackTrackSeg].GetInverseBendingMomentum())),
729 TMath::Sqrt(fHalfTrack[ibackTrackSeg].fPx*fHalfTrack[ibackTrackSeg].fPx +
730 fHalfTrack[ibackTrackSeg].fPy*fHalfTrack[ibackTrackSeg].fPy),
731 fHalfTrack[ibackTrackSeg].fPx,fHalfTrack[ibackTrackSeg].fPy,
732 fHalfTrack[ibackTrackSeg].fPz);
736 // Bits 8 and 9 must be kept zero to prevent the track ID from conflicting with other tracker components.
737 track->fId = (ibackTrackSeg << 10) | (fBackTrackSeg[ibackTrackSeg].fTrigRec & 0xFF);
738 track->fTrigRec = fBackTrackSeg[ibackTrackSeg].fTrigRec;
739 track->fPx = fHalfTrack[ibackTrackSeg].fPx;
740 track->fPy = fHalfTrack[ibackTrackSeg].fPy;
741 track->fPz = fHalfTrack[ibackTrackSeg].fPz;
745 track->fThetaY = TMath::ATan2(track->fPy, track->fPz);
746 track->fThetaX = TMath::ATan2(track->fPx, track->fPz);
752 for( Int_t ipoint=15;ipoint>=0;ipoint--){
753 track->fHit[ipoint] = AliHLTMUONConstants::NilRecHitStruct();
754 hitset[ipoint] = false;
757 for( Int_t ipoint=9;ipoint>=6;ipoint--){
759 if(fBackTrackSeg[ibackTrackSeg].fIndex[ipoint-6]!=-1 ){
761 track->fHit[ipoint] = *(fChPoint[ipoint][fBackTrackSeg[ibackTrackSeg].fIndex[ipoint-6]]);
762 hitset[ipoint] = true;
764 AliHLTUInt8_t chamber; AliHLTUInt16_t detElemID;
765 AliHLTMUONUtils::UnpackRecHitFlags((track->fHit[ipoint]).fFlags,chamber,detElemID);
766 HLTImportant("halftrack : (X,Y,Z) : (%lf,%lf,%lf)",(track->fHit[ipoint]).fX,(track->fHit[ipoint]).fY,(track->fHit[ipoint]).fZ);
770 AliHLTMUONParticleSign sign = AliHLTMUONParticleSign(fHalfTrack[ibackTrackSeg].fCharge);
771 track->fFlags = AliHLTMUONUtils::PackTrackFlags(sign,hitset);
772 TVector3 mom(track->fPx, track->fPy, track->fPz);
776 case kSignMinus: signVal = +1.; break;
777 case kSignUnknown: signVal = 0.; break;
778 case kSignPlus: signVal = -1.; break;
781 track->fInverseBendingMomentum = signVal/mom.Mag();
783 track->fInverseBendingMomentum = 0 ;
788 }//if nof connected is more than zero or not
789 }//back track seg for loop
796 ///__________________________________________________________________________
799 Bool_t AliHLTMUONFullTracker::SlatTrackSeg()
802 ///Find the Slat Track Segments
803 if(fNofPoints[6]==0 and fNofPoints[7]==0){
804 HLTDebug("No Hits found in Stn4");
806 }else if(fNofPoints[8]==0 and fNofPoints[9]==0){
807 HLTDebug("No Hits found in Stn5");
811 Float_t trigX1,trigY1,trigZ1=AliMUONConstants::DefaultChamberZ(10);
812 Float_t trigX2,trigY2,trigZ2=AliMUONConstants::DefaultChamberZ(12);
813 Float_t extrapCh9X,extrapCh9Y,extrapCh9Z=AliMUONConstants::DefaultChamberZ(9);
814 Float_t extrapCh8X,extrapCh8Y,extrapCh8Z=AliMUONConstants::DefaultChamberZ(8);
815 Float_t extrapCh7X,extrapCh7Y,extrapCh7Z=AliMUONConstants::DefaultChamberZ(7);
816 Float_t extrapCh6X,extrapCh6Y,extrapCh6Z=AliMUONConstants::DefaultChamberZ(6);
818 Double_t distChFront,distChBack;
819 Int_t nofFrontChPoints,nofBackChPoints;
820 Int_t frontIndex[fgkMaxNofTracks], backIndex[fgkMaxNofTracks];
822 AliHLTMUONRecHitStruct p2,p3,pSeg1,pSeg2,pSeg3;
823 Double_t anglediff,anglediff1,anglediff2;
824 Double_t minAngle = 2.0;
826 Bool_t st5TrackletFound = false;
827 Bool_t ch9PointFound = false;
828 Bool_t ch8PointFound = false;
829 Bool_t st4TrackletFound = false;
830 Bool_t ch7PointFound = false;
831 Bool_t ch6PointFound = false;
833 Int_t index1,index2,index3,index4;
834 IntPair cells[2][fgkMaxNofTracks]; ///cell array for 5 stn for given trigger
836 Float_t maxXDeflectionExtrap = 10.0 + 4.0; ///simulation result 10.0
837 Float_t extrapRatio = 0.2; ///simulation result 0.2
838 Float_t circularWindow = 20.0 + 5.0 + 25.0; ///simulatiuon result 20
839 Float_t minAngleWindow = 2.0 + 1.0 + 2.0; ///simulation result 2.0
842 maxXDeflectionExtrap = 10.0; ///simulation result 10.0
843 extrapRatio = 0.2; ///simulation result 0.2
844 circularWindow = 20.0 ; ///simulatiuon result 20
845 minAngleWindow = 2.0; ///simulation result 2.0
848 Float_t tx=0.0,ty=0.0;
850 AliHLTUInt16_t detElemID,prevDetElemID=0xFFFF;
851 AliHLTUInt8_t chamber;
852 Int_t minTrgCh,maxTrgCh;
855 HLTImportant("\nAliHLTMUONFullTracker::SlatTrackSeg()--Begin\n\n");
856 HLTImportant("\nAliHLTMUONFullTracker::SlatTrackSeg() : Total number of Triggers : %d\n",fNofPoints[10]);
860 for( Int_t itrig=0;itrig<fNofPoints[10];itrig++){
862 st5TrackletFound = false;
863 ch9PointFound = false;
864 ch8PointFound = false;
866 st4TrackletFound = false;
867 ch7PointFound = false;
868 ch6PointFound = false;
876 for( Int_t ihit=0;ihit<4;ihit++){
878 AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[ihit]).fFlags,chamber,detElemID);
880 if(ihit==0 and detElemID!=0)
882 if(ihit==1 and detElemID!=0 and prevDetElemID==0)
884 if(ihit==2 and detElemID!=0)
886 if(ihit==3 and detElemID!=0 and prevDetElemID==0)
889 prevDetElemID = detElemID;
892 if(minTrgCh == -1 or maxTrgCh == -1){
893 HLTDebug("Trigger hits not found in both trigger station minTrgCh : %d, maxTrgCh : %d, not harmful to HLT chain",minTrgCh,maxTrgCh);
897 if( not fFastTracking){
898 AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[minTrgCh]).fFlags,chamber,detElemID);
899 if(not fDetElemList[detElemID]){
900 HLTDebug("Invalid detection element : %d, not harmful to HLT chain",detElemID);
903 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,trigZ1);
904 AliHLTMUONUtils::UnpackRecHitFlags((fChPoint11[itrig]->fHit[maxTrgCh]).fFlags,chamber,detElemID);
905 if(not fDetElemList[detElemID]){
906 HLTDebug("Invalid detection element : %d, not harmful to HLT chain",detElemID);
909 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,trigZ2);
913 trigX1 = (fChPoint11[itrig]->fHit[minTrgCh]).fX;
914 trigY1 = (fChPoint11[itrig]->fHit[minTrgCh]).fY;
915 trigZ1 = (fChPoint11[itrig]->fHit[minTrgCh]).fZ;
917 trigX2 = (fChPoint11[itrig]->fHit[maxTrgCh]).fX;
918 trigY2 = (fChPoint11[itrig]->fHit[maxTrgCh]).fY;
919 trigZ2 = (fChPoint11[itrig]->fHit[maxTrgCh]).fZ;
922 HLTImportant("itrig : %d trig 1 : (%f,%f,%f) \n",itrig,trigX1,trigY1,trigZ1);
923 HLTImportant("itrig : %d trig 2 : (%f,%f,%f) \n",itrig,trigX2,trigY2,trigZ2);
926 /////////////////////////////////////////////////// Stn 5///////////////////////////////////////////////////////////////
930 // HLTImportant("\textrap9 : (%f,%f,%f)\n",extrapCh9X,extrapCh9Y,extrapCh9Z);
931 // HLTImportant("\textrap8 : (%f,%f,%f)\n",extrapCh8X,extrapCh8Y,extrapCh8Z);
934 nofFrontChPoints = 0; nofBackChPoints = 0;
936 extrapCh9X = trigX1 * extrapCh9Z/trigZ1 ;
937 extrapCh9Y = trigY1 + (trigY2-trigY1) * (extrapCh9Z-trigZ1)/(trigZ2 - trigZ1) ;
938 for( Int_t ipointch9=0;ipointch9<fNofPoints[9];ipointch9++){
940 if(not fFastTracking){
941 AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[9][ipointch9]->fFlags,chamber,detElemID);
942 if(not fDetElemList[detElemID]){
943 HLTDebug("Invalid detection element : %d, not harmful to HLT chain",detElemID);
946 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh9Z);
948 extrapCh9X = trigX1 * extrapCh9Z/trigZ1 ;
949 extrapCh9Y = trigY1 + (trigY2-trigY1) * (extrapCh9Z-trigZ1)/(trigZ2 - trigZ1) ;
952 if(nofBackChPoints < (fgkMaxNofTracks-1) &&
953 TMath::Abs(extrapCh9X-fChPoint[9][ipointch9]->fX) < maxXDeflectionExtrap &&
954 TMath::Abs(extrapCh9Y-fChPoint[9][ipointch9]->fY)/
955 ((fChPoint[9][ipointch9]->fX * fChPoint[9][ipointch9]->fX) +
956 (fChPoint[9][ipointch9]->fY * fChPoint[9][ipointch9]->fY)) <= extrapRatio ){
958 distChBack = sqrt((extrapCh9X-fChPoint[9][ipointch9]->fX)*(extrapCh9X-fChPoint[9][ipointch9]->fX)
959 + (extrapCh9Y-fChPoint[9][ipointch9]->fY)*(extrapCh9Y-fChPoint[9][ipointch9]->fY));
960 if(distChBack>circularWindow) continue;
963 HLTImportant("\t\tpoints selected in Ch9 : (%f,%f,%f)\n",
964 distChBack,fChPoint[9][ipointch9]->fX,fChPoint[9][ipointch9]->fY,fChPoint[9][ipointch9]->fZ);
967 backIndex[nofBackChPoints++] = ipointch9;
972 extrapCh8X = trigX1 * extrapCh8Z/trigZ1 ;
973 extrapCh8Y = trigY1 + (trigY2-trigY1) * (extrapCh8Z-trigZ1)/(trigZ2 - trigZ1) ;
974 for( Int_t ipointch8=0;ipointch8<fNofPoints[8];ipointch8++){
976 if(not fFastTracking){
977 AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[8][ipointch8]->fFlags,chamber,detElemID);
978 if(not fDetElemList[detElemID]){
979 HLTDebug("Invalid detection element : %d, not harmful to HLT chain",detElemID);
982 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh8Z);
984 extrapCh8X = trigX1 * extrapCh8Z/trigZ1 ;
985 extrapCh8Y = trigY1 + (trigY2-trigY1) * (extrapCh8Z-trigZ1)/(trigZ2 - trigZ1) ;
988 if( nofFrontChPoints < (fgkMaxNofTracks-1) &&
989 TMath::Abs(extrapCh8X-fChPoint[8][ipointch8]->fX) < maxXDeflectionExtrap &&
990 TMath::Abs(extrapCh8Y-fChPoint[8][ipointch8]->fY)/
991 ((fChPoint[8][ipointch8]->fX * fChPoint[8][ipointch8]->fX) +
992 (fChPoint[8][ipointch8]->fY * fChPoint[8][ipointch8]->fY)) <= extrapRatio ){
994 distChFront = sqrt((extrapCh8X-fChPoint[8][ipointch8]->fX)*(extrapCh8X-fChPoint[8][ipointch8]->fX)
995 + (extrapCh8Y-fChPoint[8][ipointch8]->fY)*(extrapCh8Y-fChPoint[8][ipointch8]->fY));
997 if(distChFront>circularWindow) continue;
1000 HLTImportant("\t\tpoints selected in Ch8 : (%f,%f,%f)\n",
1001 fChPoint[8][ipointch8]->fX,fChPoint[8][ipointch8]->fY,fChPoint[8][ipointch8]->fZ,distChFront);
1004 frontIndex[nofFrontChPoints++] = ipointch8;
1008 if(nofBackChPoints==0 and nofFrontChPoints==0) continue;
1010 minAngle = minAngleWindow;
1011 p3.fX = trigX1 ; p3.fY = trigY1 ; p3.fZ = trigZ1 ;
1012 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1013 Sub(&p3,fChPoint[9][backIndex[ibackpoint]],&pSeg2);
1014 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1015 Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[8][frontIndex[ifrontpoint]],&pSeg1);
1016 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2);
1017 // #ifdef PRINT_BACK
1018 // HLTImportant("\t\ttracklet-check-St5 : anglediff : %lf, minAngle : %lf\n",anglediff,minAngle);
1020 if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
1021 st5TrackletFound = true;
1022 cells[1][fNofCells[1]].fFirst = frontIndex[ifrontpoint];
1023 cells[1][fNofCells[1]].fSecond = backIndex[ibackpoint];
1026 HLTImportant("\t\ttracklet-St5 : anglediff : %lf\n",anglediff);
1027 HLTImportant("\t\t\tCh9 : (%f,%f,%f)\n",fChPoint[9][backIndex[ibackpoint]]->fX,
1028 fChPoint[9][backIndex[ibackpoint]]->fY,fChPoint[9][backIndex[ibackpoint]]->fZ);
1029 HLTImportant("\t\t\tCh8 : (%f,%f,%f)\n",fChPoint[8][frontIndex[ifrontpoint]]->fX,
1030 fChPoint[8][frontIndex[ifrontpoint]]->fY,fChPoint[8][frontIndex[ifrontpoint]]->fZ);
1032 }///anglediff condition
1039 /// If tracklet not found, search for the single space point in Ch9 or in Ch8
1040 if(!st5TrackletFound){
1042 minAngle = minAngleWindow;
1043 p3.fX = trigX2 ; p3.fY = trigY2 ; p3.fZ = trigZ2 ;
1044 p2.fX = trigX1 ; p2.fY = trigY1 ; p2.fZ = trigZ1 ;
1045 Sub(&p3,&p2,&pSeg2);
1048 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1049 Sub(&p2,fChPoint[9][backIndex[ibackpoint]],&pSeg1);
1050 anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
1051 if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
1052 ch9PointFound = true;
1053 cells[1][fNofCells[1]].fFirst = -1;
1054 cells[1][fNofCells[1]].fSecond = backIndex[ibackpoint];
1057 HLTImportant("\t\tno st tracklet and single point-Ch9 : anglediff : %lf\n",anglediff);
1063 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1064 Sub(&p2,fChPoint[8][frontIndex[ifrontpoint]],&pSeg1);
1065 anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
1066 if(anglediff<minAngle && fNofCells[1]<(fgkMaxNofTracks-1)){
1067 ch8PointFound = true;
1068 cells[1][fNofCells[1]].fFirst = frontIndex[ifrontpoint];
1069 cells[1][fNofCells[1]].fSecond = -1;
1072 HLTImportant("\t\tno st tracklet and single point-Ch8 : anglediff : %lf\n",anglediff);
1076 }///if no tracklets found condition
1079 HLTImportant("\tnofTracks found after stn 5 : %d\n",fNofCells[1]);
1082 if(!st5TrackletFound && !ch9PointFound && !ch8PointFound) continue;
1084 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1087 /////////////////////////////////////////////////// Stn 4///////////////////////////////////////////////////////////////
1089 // #ifdef PRINT_BACK
1090 // HLTImportant("\textrap7 : (%f,%f,%f)\n",extrapCh7X,extrapCh7Y,extrapCh7Z);
1091 // HLTImportant("\textrap6 : (%f,%f,%f)\n",extrapCh6X,extrapCh6Y,extrapCh6Z);
1094 nofFrontChPoints = 0; nofBackChPoints = 0;
1096 extrapCh7X = trigX1 * extrapCh7Z/trigZ1 ;
1097 extrapCh7Y = trigY1 + (trigY2-trigY1) * (extrapCh7Z-trigZ1)/(trigZ2 - trigZ1) ;
1098 for( Int_t ipointch7=0;ipointch7<fNofPoints[7];ipointch7++){
1100 if(not fFastTracking){
1101 AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[7][ipointch7]->fFlags,chamber,detElemID);
1102 if(not fDetElemList[detElemID]){
1103 HLTDebug("Invalid detection element : %d, not harmful to HLT chain",detElemID);
1106 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh7Z);
1108 extrapCh7X = trigX1 * extrapCh7Z/trigZ1 ;
1109 extrapCh7Y = trigY1 + (trigY2-trigY1) * (extrapCh7Z-trigZ1)/(trigZ2 - trigZ1) ;
1112 if( nofBackChPoints < (fgkMaxNofTracks-1) &&
1113 TMath::Abs(extrapCh7X-fChPoint[7][ipointch7]->fX) < maxXDeflectionExtrap &&
1114 TMath::Abs(extrapCh7Y-fChPoint[7][ipointch7]->fY)/
1115 ((fChPoint[7][ipointch7]->fX * fChPoint[7][ipointch7]->fX) +
1116 (fChPoint[7][ipointch7]->fY * fChPoint[7][ipointch7]->fY)) <= extrapRatio ){
1118 distChBack = sqrt((extrapCh7X-fChPoint[7][ipointch7]->fX)*(extrapCh7X-fChPoint[7][ipointch7]->fX)
1119 + (extrapCh7Y-fChPoint[7][ipointch7]->fY)*(extrapCh7Y-fChPoint[7][ipointch7]->fY));
1121 if(distChBack>circularWindow) continue;
1123 HLTImportant("\t\tpoints selected in Ch7 : (%f,%f,%f)\n",
1124 fChPoint[7][ipointch7]->fX,fChPoint[7][ipointch7]->fY,fChPoint[7][ipointch7]->fZ);
1127 backIndex[nofBackChPoints++] = ipointch7;
1131 extrapCh6X = trigX1 * extrapCh6Z/trigZ1 ;
1132 extrapCh6Y = trigY1 + (trigY2-trigY1) * (extrapCh6Z-trigZ1)/(trigZ2 - trigZ1) ;
1133 for( Int_t ipointch6=0;ipointch6<fNofPoints[6];ipointch6++){
1135 if(not fFastTracking){
1136 AliHLTMUONUtils::UnpackRecHitFlags(fChPoint[6][ipointch6]->fFlags,chamber,detElemID);
1137 if(not fDetElemList[detElemID]){
1138 HLTDebug("Invalid detection element : %d, not harmful to HLT chain",detElemID);
1141 fChamberGeometryTransformer->Local2Global(detElemID,0.0,0.0,0.0,tx,ty,extrapCh6Z);
1144 extrapCh6X = trigX1 * extrapCh6Z/trigZ1 ;
1145 extrapCh6Y = trigY1 + (trigY2-trigY1) * (extrapCh6Z-trigZ1)/(trigZ2 - trigZ1) ;
1148 if(nofFrontChPoints < (fgkMaxNofTracks-1) &&
1149 TMath::Abs(extrapCh6X-fChPoint[6][ipointch6]->fX) < maxXDeflectionExtrap &&
1150 TMath::Abs(extrapCh6Y-fChPoint[6][ipointch6]->fY)/
1151 ((fChPoint[6][ipointch6]->fX * fChPoint[6][ipointch6]->fX) +
1152 (fChPoint[6][ipointch6]->fY * fChPoint[6][ipointch6]->fY)) <= extrapRatio ){
1154 distChFront = sqrt((extrapCh6X-fChPoint[6][ipointch6]->fX)*(extrapCh6X-fChPoint[6][ipointch6]->fX)
1155 + (extrapCh6Y-fChPoint[6][ipointch6]->fY)*(extrapCh6Y-fChPoint[6][ipointch6]->fY));
1156 if(distChFront>circularWindow) continue;
1159 HLTImportant("\t\tpoints selected in Ch6 : (%f,%f,%f)\n",
1160 fChPoint[6][ipointch6]->fX,fChPoint[6][ipointch6]->fY,fChPoint[6][ipointch6]->fZ);
1162 frontIndex[nofFrontChPoints++] = ipointch6;
1166 if(nofBackChPoints==0 and nofFrontChPoints==0) continue;
1168 minAngle = minAngleWindow;
1169 p3.fX = trigX1 ; p3.fY = trigY1 ; p3.fZ = trigZ1 ;
1170 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1171 Sub(&p3,fChPoint[7][backIndex[ibackpoint]],&pSeg2);
1172 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1173 Sub(fChPoint[7][backIndex[ibackpoint]],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
1174 anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
1175 if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
1176 st4TrackletFound = true;
1177 cells[0][fNofCells[0]].fFirst = frontIndex[ifrontpoint];
1178 cells[0][fNofCells[0]].fSecond = backIndex[ibackpoint];
1181 HLTImportant("\t\ttracklet-St4 : anglediff : %lf\n",anglediff);
1182 HLTImportant("\t\t\tCh7 : (%f,%f,%f)\n",fChPoint[7][backIndex[ibackpoint]]->fX,
1183 fChPoint[7][backIndex[ibackpoint]]->fY,fChPoint[7][backIndex[ibackpoint]]->fZ);
1184 HLTImportant("\t\t\tCh6 : (%f,%f,%f)\n",fChPoint[6][frontIndex[ifrontpoint]]->fX,
1185 fChPoint[6][frontIndex[ifrontpoint]]->fY,fChPoint[6][frontIndex[ifrontpoint]]->fZ);
1194 /// If tracklet not found search for the single space point in Ch7 or in Ch6
1195 if(!st4TrackletFound){
1197 minAngle = minAngleWindow;
1198 p3.fX = trigX2 ; p3.fY = trigY2 ; p3.fZ = trigZ2 ;
1199 p2.fX = trigX1 ; p2.fY = trigY1 ; p2.fZ = trigZ1 ;
1200 Sub(&p3,&p2,&pSeg2);
1203 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1204 Sub(&p2,fChPoint[7][backIndex[ibackpoint]],&pSeg1);
1205 anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
1206 if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
1207 ch7PointFound = true;
1208 cells[0][fNofCells[0]].fFirst = -1;
1209 cells[0][fNofCells[0]].fSecond = backIndex[ibackpoint];
1212 HLTImportant("\t\tno st tracklet and single point-Ch7 : anglediff : %lf\n",anglediff);
1218 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1219 Sub(&p2,fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
1220 anglediff = TMath::RadToDeg()*Angle(&pSeg1,&pSeg2);
1221 if(anglediff<minAngle && fNofCells[0]<(fgkMaxNofTracks-1)){
1222 ch6PointFound = true;
1223 cells[0][fNofCells[0]].fFirst = frontIndex[ifrontpoint];
1224 cells[0][fNofCells[0]].fSecond = -1;
1227 HLTImportant("\t\tno st tracklet and single point-Ch6 : anglediff : %lf\n",anglediff);
1231 }///if no tracklets found condition
1234 HLTImportant("\tnofTracks found after stn 4 : %d\n",fNofCells[0]);
1237 if(!st4TrackletFound && !ch7PointFound && !ch6PointFound) continue;
1239 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1242 ////////////////////////////////////////////// Analyse and fill trackseg array////////////////////////////////////////
1245 HLTImportant("\tfNofbackTrackSeg : %d, st5TrackletFound : %d, st4TrackletFound : %d\n",fNofbackTrackSeg,st5TrackletFound,st4TrackletFound);
1248 if(st5TrackletFound && st4TrackletFound){
1250 minAngle = minAngleWindow;
1252 for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1253 index1 = cells[0][itrackletfront].fFirst ;
1254 index2 = cells[0][itrackletfront].fSecond ;
1255 Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1256 for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
1257 index3 = cells[1][itrackletback].fFirst ;
1258 index4 = cells[1][itrackletback].fSecond ;
1259 Sub(fChPoint[8][index3],fChPoint[7][index2],&pSeg2);
1260 Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
1261 anglediff = Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3);
1262 if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1263 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1264 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1265 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
1266 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
1267 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1268 minAngle = anglediff;
1270 HLTImportant("\t\ttracklet-St4 and St5 : anglediff : %lf\n",anglediff);
1271 HLTImportant("\t\t\tCh9 : (%f,%f,%f)\n",fChPoint[9][index4]->fX,
1272 fChPoint[9][index4]->fY,fChPoint[9][index4]->fZ);
1273 HLTImportant("\t\t\tCh8 : (%f,%f,%f)\n",fChPoint[8][index3]->fX,
1274 fChPoint[8][index3]->fY,fChPoint[8][index3]->fZ);
1275 HLTImportant("\t\t\tCh7 : (%f,%f,%f)\n",fChPoint[7][index2]->fX,
1276 fChPoint[7][index2]->fY,fChPoint[7][index2]->fZ);
1277 HLTImportant("\t\t\tCh6 : (%f,%f,%f)\n",fChPoint[6][index1]->fX,
1278 fChPoint[6][index1]->fY,fChPoint[6][index1]->fZ);
1282 }///for loop of back ch
1284 if(minAngle<minAngleWindow)
1287 }else if(st5TrackletFound && (ch7PointFound || ch6PointFound)){
1290 nofFrontChPoints = 0; nofBackChPoints = 0;
1291 for( Int_t ifrontpoint=0;ifrontpoint<fNofCells[0];ifrontpoint++){
1292 if(cells[0][ifrontpoint].fFirst==-1 && nofBackChPoints<(fgkMaxNofTracks-1))
1293 backIndex[nofBackChPoints++] = cells[0][ifrontpoint].fSecond;
1294 else if(cells[0][ifrontpoint].fSecond==-1 && nofFrontChPoints<(fgkMaxNofTracks-1))
1295 frontIndex[nofFrontChPoints++] = cells[0][ifrontpoint].fFirst;
1298 minAngle = minAngleWindow;
1299 if(nofFrontChPoints>0 && nofBackChPoints>0){
1300 for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
1301 index3 = cells[1][itrackletback].fFirst ;
1302 index4 = cells[1][itrackletback].fSecond ;
1303 Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
1304 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1305 Sub(fChPoint[8][index3],fChPoint[7][backIndex[ibackpoint]],&pSeg2);
1306 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1307 Sub(fChPoint[7][backIndex[ibackpoint]],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
1308 anglediff = TMath::RadToDeg()*(Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3))/2.0;
1310 if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1311 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
1312 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
1313 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
1314 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
1315 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1316 minAngle = anglediff;
1320 Sub(fChPoint[8][index3],fChPoint[6][frontIndex[ifrontpoint]],&pSeg1);
1321 anglediff1 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg3);
1322 anglediff2 = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
1324 if( anglediff1 < anglediff2 && anglediff1<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1325 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
1326 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = -1;
1327 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
1328 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
1329 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1330 minAngle = anglediff1;
1334 if( anglediff2 < anglediff1 && anglediff2<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1335 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = -1;
1336 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
1337 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
1338 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
1339 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1340 minAngle = anglediff2;
1343 }///loop of ifrontpoint
1344 }///loop of ibackpoint
1345 }/// for loop of St5 cells
1346 }else if(nofFrontChPoints>0){
1348 for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
1349 index3 = cells[1][itrackletback].fFirst ;
1350 index4 = cells[1][itrackletback].fSecond ;
1351 Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
1353 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1354 Sub(fChPoint[8][index3],fChPoint[6][frontIndex[ifrontpoint]],&pSeg2);
1356 anglediff = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
1357 if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1358 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = frontIndex[ifrontpoint];
1359 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = -1;
1360 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
1361 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
1362 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1363 minAngle = anglediff;
1368 }else{ /// if(nofBackChPoints>0){
1369 for( Int_t itrackletback=0;itrackletback<fNofCells[1];itrackletback++){
1370 index3 = cells[1][itrackletback].fFirst ;
1371 index4 = cells[1][itrackletback].fSecond ;
1372 Sub(fChPoint[9][index4],fChPoint[8][index3],&pSeg3);
1374 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1375 Sub(fChPoint[8][index3],fChPoint[7][backIndex[ibackpoint]],&pSeg2);
1377 anglediff = TMath::RadToDeg() * Angle(&pSeg2,&pSeg3);
1379 if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1380 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = -1;
1381 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = backIndex[ibackpoint] ;
1382 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = index3;
1383 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = index4;
1384 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1385 minAngle = anglediff;
1390 }///condn for if(nofFrontChPoints>0)
1392 if(minAngle<minAngleWindow)
1395 }else if((ch9PointFound || ch8PointFound) && st4TrackletFound){
1397 nofFrontChPoints = 0; nofBackChPoints = 0;
1398 for( Int_t ibackpoint=0;ibackpoint<fNofCells[1];ibackpoint++){
1399 if(cells[1][ibackpoint].fFirst==-1 && nofBackChPoints<(fgkMaxNofTracks-1))
1400 backIndex[nofBackChPoints++] = cells[1][ibackpoint].fSecond;
1401 else if(nofFrontChPoints<(fgkMaxNofTracks-1))
1402 frontIndex[nofFrontChPoints++] = cells[1][ibackpoint].fFirst;
1405 minAngle = minAngleWindow;
1406 if(nofFrontChPoints>0 && nofBackChPoints>0){
1408 for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1409 index1 = cells[0][itrackletfront].fFirst ;
1410 index2 = cells[0][itrackletfront].fSecond ;
1412 Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1414 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1416 Sub(fChPoint[8][frontIndex[ifrontpoint]],fChPoint[7][index2],&pSeg2);
1418 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1420 Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[8][frontIndex[ifrontpoint]],&pSeg3);
1422 anglediff = TMath::RadToDeg()*(Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3))/2.0;
1423 if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1424 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1425 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1426 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1427 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1428 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1429 minAngle = anglediff;
1433 Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[7][index2],&pSeg3);
1435 anglediff1 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
1436 anglediff2 = TMath::RadToDeg() * Angle(&pSeg1,&pSeg3);
1438 if( anglediff1 < anglediff2 && anglediff1<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1439 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1440 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1441 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1442 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = -1;
1443 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1444 minAngle = anglediff1;
1448 if( anglediff2 < anglediff1 && anglediff2<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1449 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1450 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1451 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = -1;
1452 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1453 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1454 minAngle = anglediff2;
1457 }///loop of ifrontpoint
1458 }///loop of ibackpoint
1459 }/// for loop of St5 cells
1460 }else if(nofFrontChPoints>0){
1462 for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1463 index1 = cells[0][itrackletfront].fFirst ;
1464 index2 = cells[0][itrackletfront].fSecond ;
1466 Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1468 for( Int_t ifrontpoint=0;ifrontpoint<nofFrontChPoints;ifrontpoint++){
1470 Sub(fChPoint[8][frontIndex[ifrontpoint]],fChPoint[7][index2],&pSeg2);
1472 anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2);
1473 if( anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1474 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1475 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;;
1476 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = frontIndex[ifrontpoint];
1477 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = -1;
1478 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1479 minAngle = anglediff;
1484 }else{ /// if(nofBackChPoints>0){
1486 for( Int_t itrackletfront=0;itrackletfront<fNofCells[0];itrackletfront++){
1487 index1 = cells[0][itrackletfront].fFirst ;
1488 index2 = cells[0][itrackletfront].fSecond ;
1490 Sub(fChPoint[7][index2],fChPoint[6][index1],&pSeg1);
1492 for( Int_t ibackpoint=0;ibackpoint<nofBackChPoints;ibackpoint++){
1493 Sub(fChPoint[9][backIndex[ibackpoint]],fChPoint[6][index1],&pSeg3);
1494 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg3);
1495 if(anglediff<minAngle && fNofbackTrackSeg<(fgkMaxNofTracks-1)){
1496 fBackTrackSeg[fNofbackTrackSeg].fIndex[0] = index1;
1497 fBackTrackSeg[fNofbackTrackSeg].fIndex[1] = index2;
1498 fBackTrackSeg[fNofbackTrackSeg].fIndex[2] = -1;
1499 fBackTrackSeg[fNofbackTrackSeg].fIndex[3] = backIndex[ibackpoint] ;
1500 fBackTrackSeg[fNofbackTrackSeg].fTrigRec = fChPoint11[itrig]->fId;
1501 minAngle = anglediff;
1506 }///condn for if(nofFrontChPoints>0)
1508 if(minAngle<minAngleWindow)
1511 }else if((ch9PointFound || ch8PointFound) && (ch7PointFound || ch6PointFound)){
1512 ///To Do : To be analysed for two points out of four slat chambers
1515 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1522 Float_t meanX1,meanX2,meanY1,meanY2,meanZ1,meanZ2;
1524 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
1527 HLTImportant("Index : (%d,%d,%d,%d) \n",fBackTrackSeg[ibacktrackseg].fIndex[0],
1528 fBackTrackSeg[ibacktrackseg].fIndex[1],fBackTrackSeg[ibacktrackseg].fIndex[2],
1529 fBackTrackSeg[ibacktrackseg].fIndex[3]);
1532 if(fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[1]!=-1 ){
1533 meanX1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fX
1534 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fX)/2.0 ;
1535 meanY1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fY
1536 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fY)/2.0 ;
1537 meanZ1 = (fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fZ
1538 + fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fZ)/2.0 ;
1539 }else if(fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[1]==-1 ){
1540 meanX1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fX ;
1541 meanY1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fY ;
1542 meanZ1 = fChPoint[6][fBackTrackSeg[ibacktrackseg].fIndex[0]]->fZ ;
1544 meanX1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fX ;
1545 meanY1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fY ;
1546 meanZ1 = fChPoint[7][fBackTrackSeg[ibacktrackseg].fIndex[1]]->fZ ;
1549 if(fBackTrackSeg[ibacktrackseg].fIndex[2]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1 ){
1550 meanX2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fX
1551 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fX)/2.0 ;
1552 meanY2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fY
1553 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fY)/2.0 ;
1554 meanZ2 = (fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fZ
1555 + fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fZ)/2.0 ;
1556 }else if(fBackTrackSeg[ibacktrackseg].fIndex[2]!=-1 && fBackTrackSeg[ibacktrackseg].fIndex[3]==-1 ){
1557 meanX2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fX ;
1558 meanY2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fY ;
1559 meanZ2 = fChPoint[8][fBackTrackSeg[ibacktrackseg].fIndex[2]]->fZ ;
1561 meanX2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fX ;
1562 meanY2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fY ;
1563 meanZ2 = fChPoint[9][fBackTrackSeg[ibacktrackseg].fIndex[3]]->fZ ;
1565 fExtrapSt3X[ibacktrackseg] = meanX1 + (fgkTrackDetCoordinate[2]-meanZ1)*(meanX2-meanX1)/(meanZ2-meanZ1);
1566 fExtrapSt3Y[ibacktrackseg] = meanY1 + (fgkTrackDetCoordinate[2]-meanZ1)*(meanY2-meanY1)/(meanZ2-meanZ1);
1567 fInclinationBack[ibacktrackseg] = (meanX2-meanX1)/(meanZ2-meanZ1) ;
1568 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
1569 }///backtrigseg loop
1572 HLTImportant("AliHLTMUONFullTracker::SlatTrackSeg()--End\n");
1573 HLTImportant("\n\n");
1578 ///__________________________________________________________________________
1580 Bool_t AliHLTMUONFullTracker::PrelimMomCalc()
1582 /// momentum calculation for half tracks
1584 Cluster clus1,clus2;
1585 Float_t xf,yf,thetaDev,zf = 0.5*(AliMUONConstants::DefaultChamberZ(4) + AliMUONConstants::DefaultChamberZ(5));
1588 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
1591 Int_t maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
1592 Int_t maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
1594 Int_t minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
1595 Int_t minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
1597 clus2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
1598 clus2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
1599 clus2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
1601 clus1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
1602 clus1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
1603 clus1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
1605 thetaDev= (1/zf)*(clus1.fY*clus2.fZ - clus2.fY*clus1.fZ)/(clus2.fZ - clus1.fZ);
1606 xf = clus1.fX*zf/clus1.fZ;
1607 yf = clus2.fY - ((clus2.fY - clus1.fY)*(clus2.fZ-zf))/(clus2.fZ - clus1.fZ);
1608 p = 3.0*0.3/thetaDev;
1611 pz = sqrt((p*p-(px*px + py*py)));
1614 fHalfTrack[ibacktrackseg].fCharge = (p>0)?-1:+1;
1615 fHalfTrack[ibacktrackseg].fPx = p*xf/zf;
1616 fHalfTrack[ibacktrackseg].fPy = p*yf/zf;
1617 fHalfTrack[ibacktrackseg].fPz = pz;
1624 ///__________________________________________________________________________
1626 Bool_t AliHLTMUONFullTracker::QuadTrackSeg()
1628 ///Find the Track Segments in the Quadrant chambers
1629 if(fNofPoints[0]==0 and fNofPoints[1]==0){
1630 HLTDebug("No Hits found in Stn4");
1632 }else if(fNofPoints[2]==0 and fNofPoints[3]==0){
1633 HLTDebug("No Hits found in Stn5");
1635 }else if(fNofbackTrackSeg==0){
1636 HLTDebug("No Hits found in Stn5 and Stn4 so no tracking is done in quadrants");
1641 Float_t meanX1,meanX2,meanY1,meanY2,meanZ1,meanZ2;
1642 Float_t expectSt3X,expectSt3Y,inclinationFront;
1644 AliHLTMUONRecHitStruct pSeg1,pSeg2,pSeg3;
1645 Double_t anglediff;///,anglediff1,anglediff2;
1646 Double_t minAngle = -1.0;
1647 Int_t indexMinAngleFront = -1;
1648 Int_t indexMinAngleBack = -1;
1649 Int_t backIndex = -1;
1651 Int_t ch3CellPoint[fgkMaxNofCellsPerCh],ch2CellPoint[fgkMaxNofCellsPerCh],nofSt2Cells=0;
1652 Int_t ch1CellPoint[fgkMaxNofCellsPerCh],ch0CellPoint[fgkMaxNofCellsPerCh],nofSt1Cells=0;
1653 Bool_t isConnected[fgkMaxNofTracks];
1655 Float_t distDiffX = 4.0; ///simulation result 4.0
1656 Float_t distDiffY = 10.0 ; ///simulation result 4.0
1657 ///float closeToY0 = 10.0; ///simulation result 10.0
1658 Float_t minAngleWindow = 2.0 ; ///simulation result 2.0
1659 Float_t diffDistStX = 25.0; ///simulation result 25.0
1660 Float_t diffDistStY = 75.0; ///simulation result 25.0
1661 Float_t st3WindowX = 40.0 ; ///simulation result 40.0
1662 Float_t st3WindowY = 10.0; ///simulation result 10.0
1665 distDiffX = 4.0; ///simulation result 4.0
1666 distDiffY = 4.0 ; ///simulation result 4.0
1667 ///float closeToY0 = 10.0; ///simulation result 10.0
1668 minAngleWindow = 2.0 ; ///simulation result 2.0
1669 diffDistStX = 25.0; ///simulation result 25.0
1670 diffDistStY = 25.0; ///simulation result 25.0
1671 st3WindowX = 40.0 ; ///simulation result 40.0
1672 st3WindowY = 10.0; ///simulation result 10.0
1675 /// Float_t inclinationWindow = 0.04; ///inclination window
1676 /// Float_t st3WindowXOp2 = 40.0 ; ///simulation result 40.0
1677 /// Float_t st3WindowYOp2 = 10.0; ///simulation result 10.0
1681 HLTImportant("\nAliHLTMUONFullTracker::QuadTrackSeg()--Begin\n\n");
1682 HLTImportant("Number to back track segment found in St4 and 5 : %d\n",fNofbackTrackSeg);
1685 for( Int_t ibackpoint=0;ibackpoint<fNofPoints[3];ibackpoint++){
1686 for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[2];ifrontpoint++){
1688 if(TMath::Abs(fChPoint[3][ibackpoint]->fX - fChPoint[2][ifrontpoint]->fX)<distDiffX
1689 && TMath::Abs(fChPoint[3][ibackpoint]->fY - fChPoint[2][ifrontpoint]->fY)<distDiffY){
1692 /// if((TMath::Abs(fChPoint[3][ibackpoint]->fY) > closeToY0 &&
1693 /// TMath::Abs(fChPoint[3][ibackpoint]->fY) < TMath::Abs(fChPoint[2][ifrontpoint]->fY)) ||
1694 /// nofSt2Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1696 if(nofSt2Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1698 // #ifdef PRINT_FRONT
1699 // HLTImportant("\t\t\tCh3 : %d, (%f,%f,%f)\n",
1700 // nofSt2Cells,fChPoint[3][ibackpoint]->fX,fChPoint[3][ibackpoint]->fY,fChPoint[3][ibackpoint]->fZ);
1701 // HLTImportant("\t\t\tCh2 :(%f,%f,%f)\n\n",
1702 // fChPoint[2][ifrontpoint]->fX,fChPoint[2][ifrontpoint]->fY,fChPoint[2][ifrontpoint]->fZ);
1705 ch3CellPoint[nofSt2Cells] = ibackpoint;
1706 ch2CellPoint[nofSt2Cells] = ifrontpoint;
1713 // if(nofSt2Cells==0){
1714 // HLTDebug("No Hits found in Stn2");
1718 for( Int_t ibackpoint=0;ibackpoint<fNofPoints[1];ibackpoint++){
1719 for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[0];ifrontpoint++){
1721 if(TMath::Abs(fChPoint[1][ibackpoint]->fX - fChPoint[0][ifrontpoint]->fX)< distDiffX
1722 && TMath::Abs(fChPoint[1][ibackpoint]->fY - fChPoint[0][ifrontpoint]->fY)<distDiffY){
1725 /// if((TMath::Abs(fChPoint[1][ibackpoint]->fY) > closeToY0 &&
1726 /// TMath::Abs(fChPoint[1][ibackpoint]->fY) < TMath::Abs(fChPoint[0][ifrontpoint]->fY)) ||
1727 /// nofSt1Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1729 if(nofSt1Cells >= (fgkMaxNofCellsPerCh-1)) continue;
1732 // #ifdef PRINT_FRONT
1733 // HLTImportant("\t\t\tCh1 : %d, (%f,%f,%f)\n",
1734 // nofSt1Cells,fChPoint[1][ibackpoint]->fX,fChPoint[1][ibackpoint]->fY,fChPoint[1][ibackpoint]->fZ);
1735 // HLTImportant("\t\t\tCh0 :(%f,%f,%f)\n\n",
1736 // fChPoint[0][ifrontpoint]->fX,fChPoint[0][ifrontpoint]->fY,fChPoint[0][ifrontpoint]->fZ);
1738 ch1CellPoint[nofSt1Cells] = ibackpoint;
1739 ch0CellPoint[nofSt1Cells] = ifrontpoint;
1744 if(nofSt1Cells==0 and nofSt2Cells==0){
1745 HLTDebug("No Hit pair found in Stn1 and St2");
1750 HLTImportant("\tnofSt1Cells : %d, nofSt2Cells : %d\n",nofSt1Cells,nofSt2Cells);
1753 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++)
1754 isConnected[ibacktrackseg] = false;
1757 ///First Check for Tracklets in two front stations
1758 for( Int_t itrackletfront=0;itrackletfront<nofSt1Cells;itrackletfront++){
1760 Sub(fChPoint[1][ch1CellPoint[itrackletfront]],fChPoint[0][ch0CellPoint[itrackletfront]],&pSeg1);
1762 minAngle = minAngleWindow;
1763 indexMinAngleBack = -1;
1764 indexMinAngleFront = -1;
1766 meanX1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fX
1767 + fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/2.0 ;
1768 meanY1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fY
1769 + fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/2.0 ;
1770 meanZ1 = (fChPoint[0][ch0CellPoint[itrackletfront]]->fZ
1771 + fChPoint[1][ch1CellPoint[itrackletfront]]->fZ)/2.0 ;
1773 for( Int_t itrackletback=0;itrackletback<nofSt2Cells;itrackletback++){
1774 // #ifdef PRINT_FRONT
1775 // cout<<"\tBefore "
1776 // <<TMath::Abs(fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)
1778 // <<TMath::Abs(fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)
1781 if(TMath::Abs(fChPoint[2][ch2CellPoint[itrackletback]]->fX -
1782 fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX ||
1783 TMath::Abs(fChPoint[2][ch2CellPoint[itrackletback]]->fY -
1784 fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1786 meanX2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fX
1787 + fChPoint[3][ch3CellPoint[itrackletback]]->fX)/2.0 ;
1788 meanY2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fY
1789 + fChPoint[3][ch3CellPoint[itrackletback]]->fY)/2.0 ;
1790 meanZ2 = (fChPoint[2][ch2CellPoint[itrackletback]]->fZ
1791 + fChPoint[3][ch3CellPoint[itrackletback]]->fZ)/2.0 ;
1793 expectSt3X = meanX2 + (fgkTrackDetCoordinate[2]-meanZ2)*(meanX2-meanX1)/(meanZ2-meanZ1);
1794 expectSt3Y = meanY2 + (fgkTrackDetCoordinate[2]-meanZ2)*(meanY2-meanY1)/(meanZ2-meanZ1);
1795 inclinationFront = (meanX2-meanX1)/(meanZ2-meanZ1) ;
1797 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
1799 if( /// TMath::Abs(inclinationBack[ibacktrackseg]-inclinationFront)<0.04 &&
1800 TMath::Abs((expectSt3X-fExtrapSt3X[ibacktrackseg])) < st3WindowX
1801 && TMath::Abs((expectSt3Y-fExtrapSt3Y[ibacktrackseg])) < st3WindowY){
1803 Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1804 Sub(fChPoint[3][ch3CellPoint[itrackletback]],fChPoint[2][ch2CellPoint[itrackletback]],&pSeg3);
1806 anglediff = TMath::RadToDeg()* (Angle(&pSeg1,&pSeg2) + Angle(&pSeg2,&pSeg3));
1808 // #ifdef PRINT_FRONT
1809 // HLTImportant("\t\t\tanglediff : %lf\n",anglediff);
1811 if(anglediff<minAngle){
1812 minAngle = anglediff;
1813 indexMinAngleBack = itrackletback;
1814 indexMinAngleFront = itrackletfront;
1815 backIndex = ibacktrackseg;
1816 isConnected[ibacktrackseg] = true;
1818 }///matching tracklet
1819 }///for loop of backtrackseg
1823 if(minAngle < minAngleWindow && indexMinAngleFront!=-1
1824 && indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)
1825 && fNofConnectedfrontTrackSeg[backIndex]<fgkMaxNofConnectedTracks ){
1827 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = ch0CellPoint[indexMinAngleFront];
1828 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = ch1CellPoint[indexMinAngleFront];
1829 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = ch2CellPoint[indexMinAngleBack];
1830 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = ch3CellPoint[indexMinAngleBack];
1832 fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;
1833 fNoffrontTrackSeg++;
1834 }///condition to find valid tracklet
1836 }///for loop of front ch
1838 Int_t nofNCfBackTrackSeg = 0;
1839 Int_t ncfBackTrackSeg[fgkMaxNofTracks];
1843 HLTImportant("\tfNofConnected : %d, nofNCfBackTrackSeg : %d, fNoffrontTrackSeg : %d\n",
1844 fNofConnected,nofNCfBackTrackSeg,fNoffrontTrackSeg);
1847 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++)
1848 if(!isConnected[ibacktrackseg])
1849 ncfBackTrackSeg[nofNCfBackTrackSeg++] = ibacktrackseg;
1855 HLTImportant("\tfNofConnected : %d, nofNCfBackTrackSeg : %d, fNoffrontTrackSeg : %d\n",fNofConnected,nofNCfBackTrackSeg,fNoffrontTrackSeg);
1856 HLTImportant("\tfNofPoints[3] : %d, fNofPoints[2] : %d\n",fNofPoints[3],fNofPoints[2]);
1857 if(nofNCfBackTrackSeg==0){
1858 HLTImportant("All fBackTrackSegs are connected with fFrontTrackSegs, no need to search further\n");
1859 HLTImportant("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
1863 if(nofNCfBackTrackSeg==0) return true;
1866 ///Next Check for tracklet in Stn1 and space point in Stn2
1867 Bool_t isbackpoint=false,isfrontpoint=false;
1868 for( Int_t itrackletfront=0;itrackletfront<nofSt1Cells;itrackletfront++){
1869 Sub(fChPoint[1][ch1CellPoint[itrackletfront]],fChPoint[0][ch0CellPoint[itrackletfront]],&pSeg1);
1870 minAngle = minAngleWindow;
1871 indexMinAngleBack = -1;
1872 indexMinAngleFront = -1;
1874 for( Int_t ibackpoint=0;ibackpoint<fNofPoints[3];ibackpoint++){
1875 if(/// hasCh3Cells[ibackpoint] == true &&
1876 TMath::Abs(fChPoint[3][ibackpoint]->fX -
1877 fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX ||
1878 TMath::Abs(fChPoint[3][ibackpoint]->fY -
1879 fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1881 expectSt3X = fChPoint[3][ibackpoint]->fX + (fgkTrackDetCoordinate[2] - fChPoint[3][ibackpoint]->fZ)*
1882 (fChPoint[3][ibackpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1883 (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1884 expectSt3Y = fChPoint[3][ibackpoint]->fY + (fgkTrackDetCoordinate[2] - fChPoint[3][ibackpoint]->fZ)*
1885 (fChPoint[3][ibackpoint]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/
1886 (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1887 inclinationFront = (fChPoint[3][ibackpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1888 (fChPoint[3][ibackpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ) ;
1890 for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++){
1892 if(/// TMath::Abs(inclinationBack[ncfBackTrackSeg[ibacktrackseg]]-inclinationFront)< inclinationWindow &&
1893 TMath::Abs((expectSt3X-fExtrapSt3X[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1894 && TMath::Abs((expectSt3Y-fExtrapSt3Y[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowY){
1896 Sub(fChPoint[3][ibackpoint],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1898 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
1899 // #ifdef PRINT_FRONT
1900 // HLTImportant("\t\t annglediff(Ch4) : %lf\n",anglediff);
1902 if(anglediff<minAngle){
1903 minAngle = anglediff;
1904 indexMinAngleBack = ibackpoint;
1905 indexMinAngleFront = itrackletfront;
1906 backIndex = ncfBackTrackSeg[ibacktrackseg];
1907 isConnected[ncfBackTrackSeg[ibacktrackseg]] = true;
1909 isfrontpoint = false;
1911 }///matching tracklet
1912 }///loop on not found back trackseg
1915 for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[2];ifrontpoint++){
1916 if(/// hasCh2Cells[ifrontpoint] == true &&
1917 TMath::Abs(fChPoint[2][ifrontpoint]->fX -
1918 fChPoint[1][ch1CellPoint[itrackletfront]]->fX) > diffDistStX ||
1919 TMath::Abs(fChPoint[2][ifrontpoint]->fY -
1920 fChPoint[1][ch1CellPoint[itrackletfront]]->fY) > diffDistStY ) continue;
1922 expectSt3X = fChPoint[2][ifrontpoint]->fX + (fgkTrackDetCoordinate[2] - fChPoint[2][ifrontpoint]->fZ)*
1923 (fChPoint[2][ifrontpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1924 (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1925 expectSt3Y = fChPoint[2][ifrontpoint]->fY + (fgkTrackDetCoordinate[2] - fChPoint[2][ifrontpoint]->fZ)*
1926 (fChPoint[2][ifrontpoint]->fY - fChPoint[1][ch1CellPoint[itrackletfront]]->fY)/
1927 (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ);
1928 inclinationFront = (fChPoint[2][ifrontpoint]->fX - fChPoint[1][ch1CellPoint[itrackletfront]]->fX)/
1929 (fChPoint[2][ifrontpoint]->fZ - fChPoint[1][ch1CellPoint[itrackletfront]]->fZ) ;
1930 // #ifdef PRINT_FRONT
1931 // HLTImportant("\t\texpectSt3X : %f, expectSt3Y : %f, inclinationFront : %f\n",expectSt3X,expectSt3Y,inclinationFront);
1934 for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++){
1936 if( /// TMath::Abs(inclinationBack[ncfBackTrackSeg[ibacktrackseg]]-inclinationFront)< inclinationWindow &&
1937 TMath::Abs((expectSt3X-fExtrapSt3X[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
1938 && TMath::Abs((expectSt3Y-fExtrapSt3Y[ncfBackTrackSeg[ibacktrackseg]])) < st3WindowY){
1940 Sub(fChPoint[2][ifrontpoint],fChPoint[1][ch1CellPoint[itrackletfront]],&pSeg2);
1942 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
1943 // #ifdef PRINT_FRONT
1944 // HLTImportant("\t\t annglediff(Ch3) : %lf\n",anglediff);
1946 if(anglediff<minAngle){
1947 minAngle = anglediff;
1948 indexMinAngleBack = ifrontpoint;
1949 indexMinAngleFront = itrackletfront;
1950 backIndex = ncfBackTrackSeg[ibacktrackseg];
1951 isConnected[ncfBackTrackSeg[ibacktrackseg]] = true;
1952 isbackpoint = false;
1953 isfrontpoint = true;
1956 }///matching tracklet
1957 }///loop on not found back trackseg
1960 if(minAngle < minAngleWindow && indexMinAngleFront!=-1 &&
1961 indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)
1962 && fNofConnectedfrontTrackSeg[backIndex]<fgkMaxNofConnectedTracks){
1964 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = ch0CellPoint[indexMinAngleFront];
1965 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = ch1CellPoint[indexMinAngleFront];
1966 if(isfrontpoint && !isbackpoint){
1967 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = indexMinAngleBack;
1968 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = -1;
1970 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = -1;
1971 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = indexMinAngleBack;
1973 fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;
1974 fNoffrontTrackSeg++;
1975 }///condition to find valid tracklet
1980 Int_t nofSNCfBackTrackSeg = 0;
1981 Int_t sncfBackTrackSeg[fgkMaxNofTracks];
1983 for( Int_t ibacktrackseg=0;ibacktrackseg<nofNCfBackTrackSeg;ibacktrackseg++)
1984 if(!isConnected[ncfBackTrackSeg[ibacktrackseg]])
1985 sncfBackTrackSeg[nofSNCfBackTrackSeg++] = ncfBackTrackSeg[ibacktrackseg];
1990 HLTImportant("\tfNofConnected : %d, nofSNCfBackTrackSeg : %d, fNoffrontTrackSeg : %d\n",
1991 fNofConnected,nofSNCfBackTrackSeg,fNoffrontTrackSeg);
1992 if(nofSNCfBackTrackSeg==0){
1993 HLTImportant("All fBackTrackSegs are connected with fFrontTrackSegs, no need to search further\n");
1994 HLTImportant("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
1998 if(nofSNCfBackTrackSeg==0) return true;
2000 ///Last Check for tracklet in Stn2 and space point in Stn1
2001 for( Int_t itrackletback=0;itrackletback<nofSt2Cells;itrackletback++){
2002 Sub(fChPoint[3][ch3CellPoint[itrackletback]],fChPoint[2][ch2CellPoint[itrackletback]],&pSeg1);
2003 minAngle = minAngleWindow ;
2004 indexMinAngleBack = -1;
2005 indexMinAngleFront = -1;
2007 for( Int_t ibackpoint=0;ibackpoint<fNofPoints[1];ibackpoint++){
2008 if(/// hasCh1Cells[ibackpoint] == true &&
2009 TMath::Abs(fChPoint[2][ch2CellPoint[itrackletback]]->fX -
2010 fChPoint[1][ibackpoint]->fX) > diffDistStX ||
2011 TMath::Abs(fChPoint[2][ch2CellPoint[itrackletback]]->fY -
2012 fChPoint[1][ibackpoint]->fY) > diffDistStY) continue;
2014 expectSt3X = fChPoint[2][ch2CellPoint[itrackletback]]->fX + (fgkTrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
2015 (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ibackpoint]->fX)/
2016 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ);
2017 expectSt3Y = fChPoint[2][ch2CellPoint[itrackletback]]->fY + (fgkTrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
2018 (fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[1][ibackpoint]->fY)/
2019 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ);
2020 inclinationFront = (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[1][ibackpoint]->fX)/
2021 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[1][ibackpoint]->fZ) ;
2023 for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++){
2025 if(/// TMath::Abs(inclinationBack[sncfBackTrackSeg[ibacktrackseg]]-inclinationFront)<inclinationWindow &&
2026 TMath::Abs((expectSt3X-fExtrapSt3X[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
2027 && TMath::Abs((expectSt3Y-fExtrapSt3Y[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowY ){
2029 Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[1][ibackpoint],&pSeg2);
2031 anglediff = TMath::RadToDeg()* Angle(&pSeg1,&pSeg2) ;
2032 if(anglediff<minAngle){
2033 minAngle = anglediff;
2034 indexMinAngleBack = itrackletback;
2035 indexMinAngleFront = ibackpoint;
2036 backIndex = sncfBackTrackSeg[ibacktrackseg];
2037 isConnected[sncfBackTrackSeg[ibacktrackseg]] = true;
2039 isfrontpoint = false;
2041 }///matching tracklet
2042 }///loop on not found back trackseg
2045 for( Int_t ifrontpoint=0;ifrontpoint<fNofPoints[0];ifrontpoint++){
2046 if(/// hasCh0Cells[ifrontpoint] == true &&
2047 TMath::Abs(fChPoint[2][ch2CellPoint[itrackletback]]->fX -
2048 fChPoint[0][ifrontpoint]->fX) > diffDistStX ||
2049 TMath::Abs(fChPoint[2][ch2CellPoint[itrackletback]]->fY -
2050 fChPoint[0][ifrontpoint]->fY) > diffDistStY ) continue;
2052 expectSt3X = fChPoint[2][ch2CellPoint[itrackletback]]->fX + (fgkTrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
2053 (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[0][ifrontpoint]->fX)/
2054 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ);
2055 expectSt3Y = fChPoint[2][ch2CellPoint[itrackletback]]->fY + (fgkTrackDetCoordinate[2] - fChPoint[2][ch2CellPoint[itrackletback]]->fZ)*
2056 (fChPoint[2][ch2CellPoint[itrackletback]]->fY - fChPoint[0][ifrontpoint]->fY)/
2057 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ);
2058 inclinationFront = (fChPoint[2][ch2CellPoint[itrackletback]]->fX - fChPoint[0][ifrontpoint]->fX)/
2059 (fChPoint[2][ch2CellPoint[itrackletback]]->fZ - fChPoint[0][ifrontpoint]->fZ) ;
2061 for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++){
2063 if(/// TMath::Abs(inclinationBack[sncfBackTrackSeg[ibacktrackseg]]-inclinationFront)<inclinationWindow &&
2064 TMath::Abs((expectSt3X-fExtrapSt3X[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowX
2065 && TMath::Abs((expectSt3Y-fExtrapSt3Y[sncfBackTrackSeg[ibacktrackseg]])) < st3WindowY ){
2067 Sub(fChPoint[2][ch2CellPoint[itrackletback]],fChPoint[0][ifrontpoint],&pSeg2);
2069 anglediff = TMath::RadToDeg() * Angle(&pSeg1,&pSeg2) ;
2070 if(anglediff<minAngle){
2071 minAngle = anglediff;
2072 indexMinAngleBack = itrackletback;
2073 indexMinAngleFront = ifrontpoint;
2074 backIndex = sncfBackTrackSeg[ibacktrackseg];
2075 isConnected[sncfBackTrackSeg[ibacktrackseg]] = true;
2076 isbackpoint = false;
2077 isfrontpoint = true;
2080 }///matching tracklet
2081 }///loop on not found back trackseg
2084 if(minAngle < minAngleWindow && indexMinAngleFront!=-1 &&
2085 indexMinAngleBack!=-1 && fNoffrontTrackSeg<(fgkMaxNofTracks-1)
2086 && fNofConnectedfrontTrackSeg[backIndex]<fgkMaxNofConnectedTracks){
2089 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = indexMinAngleFront;
2090 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = -1;
2092 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[0] = -1;
2093 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[1] = indexMinAngleFront;
2096 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[2] = ch2CellPoint[indexMinAngleBack];
2097 fFrontTrackSeg[fNoffrontTrackSeg].fIndex[3] = ch3CellPoint[indexMinAngleBack];
2099 HLTImportant("backIndex : %d, fNofConnectedfrontTrackSeg[backIndex] : %d, fNoffrontTrackSeg : %d",
2100 backIndex,fNofConnectedfrontTrackSeg[backIndex],fNoffrontTrackSeg);
2103 fBackToFront[backIndex][fNofConnectedfrontTrackSeg[backIndex]++] = fNoffrontTrackSeg;
2104 fNoffrontTrackSeg++;
2105 }///condition to find valid tracklet
2109 for( Int_t ibacktrackseg=0;ibacktrackseg<nofSNCfBackTrackSeg;ibacktrackseg++)
2110 if(isConnected[sncfBackTrackSeg[ibacktrackseg]])
2114 HLTImportant("\tfNofConnected : %d\n",fNofConnected);
2115 HLTImportant("Three spacepoints are found in fFrontTrackSegs\n");
2116 HLTImportant("AliHLTMUONFullTracker::QuadTrackSeg()--End\n\n");
2123 ///__________________________________________________________________________
2125 Double_t AliHLTMUONFullTracker::KalmanFilter(AliMUONTrackParam &trackParamAtCluster, Cluster *cluster)
2127 //// Compute new track parameters and their covariances including new cluster using kalman filter
2128 //// return the additional track chi2
2130 #ifdef PRINT_DETAIL_KALMAN
2131 HLTImportant("AliHLTMUONFullTracker::KalmanFilter()--Begin\n\n");
2135 /// Get actual track parameters (p)
2136 TMatrixD param(trackParamAtCluster.GetParameters());
2137 #ifdef PRINT_DETAIL_KALMAN
2138 Info("\tKalmanFilter","param.Print() [p]");
2140 HLTImportant("GetZ : %lf\n",trackParamAtCluster.GetZ());
2145 /// Get new cluster parameters (m)
2146 ///AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
2147 TMatrixD clusterParam(5,1);
2148 clusterParam.Zero();
2149 /// clusterParam(0,0) = cluster->GetX();
2150 /// clusterParam(2,0) = cluster->GetY();
2151 clusterParam(0,0) = cluster->fX;
2152 clusterParam(2,0) = cluster->fY;
2153 #ifdef PRINT_DETAIL_KALMAN
2154 Info("\tKalmanFilter","clusterParam.Print() [m]");
2155 clusterParam.Print();
2160 /// Compute the actual parameter weight (W)
2161 TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
2162 #ifdef PRINT_DETAIL_KALMAN
2163 Info("\tKalmanFilter","Covariance : [C]");
2164 paramWeight.Print();
2167 if (paramWeight.Determinant() != 0) {
2168 paramWeight.Invert();
2170 Warning("KalmanFilter"," Determinant = 0");
2174 #ifdef PRINT_DETAIL_KALMAN
2175 Info("\tKalmanFilter","Weight Matrix inverse of Covariance [W = C^-1]");
2176 paramWeight.Print();
2180 /// Compute the new cluster weight (U)
2181 TMatrixD clusterWeight(5,5);
2182 clusterWeight.Zero();
2183 clusterWeight(0,0) = 1. / cluster->fErrX2;
2184 clusterWeight(2,2) = 1. / cluster->fErrY2;
2185 #ifdef PRINT_DETAIL_KALMAN
2186 Info("\tKalmanFilter","clusterWeight.Print() [U]");
2187 HLTImportant("\tErrX2 : %lf, ErrY2 : %lf\n",cluster->fErrX2,cluster->fErrY2);
2188 clusterWeight.Print();
2194 /// Compute the new parameters covariance matrix ( (W+U)^-1 )
2195 TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
2196 #ifdef PRINT_DETAIL_KALMAN
2197 Info("\tKalmanFilter","newParamCov.Print() [(W+U)]");
2198 newParamCov.Print();
2200 if (newParamCov.Determinant() != 0) {
2201 newParamCov.Invert();
2203 Warning("RunKalmanFilter"," Determinant = 0");
2206 #ifdef PRINT_DETAIL_KALMAN
2207 Info("\tKalmanFilter","newParamCov.Print() [(W+U)^-1] (new covariances[W] for trackParamAtCluster)");
2208 newParamCov.Print();
2211 /// Save the new parameters covariance matrix
2212 trackParamAtCluster.SetCovariances(newParamCov);
2214 /// Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
2215 TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
2216 TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); /// U(m-p)
2217 TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); /// ((W+U)^-1)U(m-p)
2218 newParam += param; /// ((W+U)^-1)U(m-p) + p
2219 #ifdef PRINT_DETAIL_KALMAN
2220 Info("\tKalmanFilter","newParam.Print() [p' = ((W+U)^-1)U(m-p) + p] (new parameters[p] for trackParamAtCluster)");
2224 /// Save the new parameters
2225 trackParamAtCluster.SetParameters(newParam);
2226 /// HLTImportant(Form("Pt : %lf\n",TMath::Sqrt(trackParamAtCluster.Px()*trackParamAtCluster.Px() +
2227 /// trackParamAtCluster.Py()*trackParamAtCluster.Py())));
2230 /// Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
2231 tmp = newParam; /// p'
2232 tmp -= param; /// (p'-p)
2233 TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); /// W(p'-p)
2234 TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); /// ((p'-p)^-1)W(p'-p)
2235 tmp = newParam; /// p'
2236 tmp -= clusterParam; /// (p'-m)
2237 TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); /// U(p'-m)
2238 addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); /// ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
2240 #ifdef PRINT_DETAIL_KALMAN
2241 Info("\tKalmanFilter","addChi2Track.Print() [additional chi2 = ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)))]");
2242 addChi2Track.Print();
2243 HLTImportant("AliHLTMUONFullTracker::KalmanFilter()--End\n\n");
2247 return addChi2Track(0,0);
2250 ///__________________________________________________________________________
2252 Double_t AliHLTMUONFullTracker::TryOneCluster(const AliMUONTrackParam &trackParam, Cluster* cluster,
2253 AliMUONTrackParam &trackParamAtCluster, Bool_t updatePropagator)
2255 //// Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
2256 //// return the corresponding Chi2
2257 //// return trackParamAtCluster
2259 /// extrapolate track parameters and covariances at the z position of the tested cluster
2260 trackParamAtCluster = trackParam;
2261 AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->fZ, updatePropagator);
2263 /// set pointer to cluster into trackParamAtCluster
2264 ///trackParamAtCluster.SetClusterPtr(cluster);
2266 /// Set differences between trackParam and cluster in the bending and non bending directions
2267 Double_t dX = cluster->fX - trackParamAtCluster.GetNonBendingCoor();
2268 Double_t dY = cluster->fY - trackParamAtCluster.GetBendingCoor();
2270 /// Calculate errors and covariances
2271 const TMatrixD& kParamCov = trackParamAtCluster.GetCovariances();
2272 Double_t sigmaX2 = kParamCov(0,0) + cluster->fErrX2;
2273 Double_t sigmaY2 = kParamCov(2,2) + cluster->fErrY2;
2276 return dX * dX / sigmaX2 + dY * dY / sigmaY2;
2280 ///__________________________________________________________________________
2282 Bool_t AliHLTMUONFullTracker::TryOneClusterFast(const AliMUONTrackParam &trackParam, const Cluster* cluster)
2284 //// Test the compatibility between the track and the cluster
2285 //// given the track resolution + the maximum-distance-to-track value
2286 //// and assuming linear propagation of the track:
2287 //// return kTRUE if they are compatibles
2289 Float_t sigmaCutForTracking = 6.0;
2290 Float_t maxNonBendingDistanceToTrack = 1.0;
2291 Float_t maxBendingDistanceToTrack = 1.0;
2293 Double_t dZ = cluster->fZ - trackParam.GetZ();
2294 Double_t dX = cluster->fX - (trackParam.GetNonBendingCoor() + trackParam.GetNonBendingSlope() * dZ);
2295 Double_t dY = cluster->fY - (trackParam.GetBendingCoor() + trackParam.GetBendingSlope() * dZ);
2296 const TMatrixD& kParamCov = trackParam.GetCovariances();
2297 Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1);
2298 Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3);
2300 Double_t dXmax = sigmaCutForTracking * TMath::Sqrt(errX2) +
2301 maxNonBendingDistanceToTrack;
2302 Double_t dYmax = sigmaCutForTracking * TMath::Sqrt(errY2) +
2303 maxBendingDistanceToTrack;
2305 if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE;
2310 ///__________________________________________________________________________
2312 void AliHLTMUONFullTracker::PropagateTracks(Double_t charge, Float_t& px, Float_t& py, Float_t& pz,
2313 Float_t& xr, Float_t& yr, Float_t& zr, Float_t zprop)
2316 /// propagate in magnetic field between hits of indices i1 and i2
2319 Double_t vect[7], vout[7];
2320 Double_t step = -5.0;
2321 Double_t zMax = zprop+.5;
2326 vect[6] = sqrt((px)*(px) + (py)*(py) + (pz)*(pz));
2327 vect[3] = px/vect[6];
2328 vect[4] = py/vect[6];
2329 vect[5] = pz/vect[6];
2330 /// cout<<"vec[2] : "<<vect[2]<<", zMax : "<<zMax<<endl;
2333 while ((vect[2] < zMax) && (nSteps < 10000)) {
2335 /// OneStepRungekutta(charge, step, vect, vout);
2336 OneStepHelix3(charge,step,vect,vout);
2337 ///SetPoint(fCount,vout[0],vout[1],vout[2]);
2339 /// HLTImportant("(x,y,z) : (%f,%f,%f)\n",vout[0],vout[1],vout[2]);
2340 for (Int_t i = 0; i < 7; i++) {
2349 px = vect[3]*vect[6];
2350 py = vect[4]*vect[6];
2351 pz = vect[5]*vect[6];
2355 ///__________________________________________________________________________
2357 void AliHLTMUONFullTracker::OneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout) const
2360 //// ******************************************************************
2362 //// * Tracking routine in a constant field oriented *
2363 //// * along axis 3 *
2364 //// * Tracking is performed with a conventional *
2365 //// * helix step method *
2367 //// * ==>Called by : USER, GUSWIM *
2368 //// * Authors R.Brun, M.Hansroul ********* *
2369 //// * Rewritten V.Perevoztchikov
2371 //// ******************************************************************
2375 Double_t h4, hp, rho, tet;
2376 Double_t sint, sintt, tsint, cos1t;
2377 Double_t f1, f2, f3, f4, f5, f6;
2379 const Int_t kix = 0;
2380 const Int_t kiy = 1;
2381 const Int_t kiz = 2;
2382 const Int_t kipx = 3;
2383 const Int_t kipy = 4;
2384 const Int_t kipz = 5;
2385 const Int_t kipp = 6;
2387 const Double_t kec = 2.9979251e-4;
2390 /// ------------------------------------------------------------------
2392 /// units are kgauss,centimeters,gev/c
2394 vout[kipp] = vect[kipp];
2397 hxp[0] = - vect[kipy];
2398 hxp[1] = + vect[kipx];
2402 rho = -h4/vect[kipp];
2404 if (TMath::Abs(tet) > 0.15) {
2405 sint = TMath::Sin(tet);
2407 tsint = (tet-sint)/tet;
2408 cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet;
2410 tsint = tet*tet/36.;
2411 sintt = (1. - tsint);
2418 f3 = step * tsint * hp;
2421 f6 = tet * cos1t * hp;
2423 vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0];
2424 vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1];
2425 vout[kiz] = vect[kiz] + f1*vect[kipz] + f3;
2427 vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0];
2428 vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1];
2429 vout[kipz] = vect[kipz] + f4*vect[kipz] + f6;
2434 ///__________________________________________________________________________
2436 ///______________________________________________________________________________
2438 void AliHLTMUONFullTracker::OneStepRungekutta(Double_t charge, Double_t step,
2439 const Double_t* vect, Double_t* vout)
2441 //// ******************************************************************
2443 //// * Runge-Kutta method for tracking a particle through a magnetic *
2444 //// * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
2445 //// * Standards, procedure 25.5.20) *
2447 //// * Input parameters *
2448 //// * CHARGE Particle charge *
2449 //// * STEP Step size *
2450 //// * VECT Initial co-ords,direction cosines,momentum *
2451 //// * Output parameters *
2452 //// * VOUT Output co-ords,direction cosines,momentum *
2453 //// * User routine called *
2454 //// * CALL GUFLD(X,F) *
2456 //// * ==>Called by : <USER>, GUSWIM *
2457 //// * Authors R.Brun, M.Hansroul ********* *
2458 //// * V.Perevoztchikov (CUT STEP implementation) *
2461 //// ******************************************************************
2463 Double_t h2, h4, f[4];
2464 Double_t xyzt[3], a, b, c, ph,ph2;
2465 Double_t secxs[4],secys[4],seczs[4],hxp[3];
2466 Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
2467 Double_t est, at, bt, ct, cba;
2468 Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
2478 Double_t maxit = 1992;
2479 Double_t maxcut = 11;
2481 const Double_t kdlt = 1e-4;
2482 const Double_t kdlt32 = kdlt/32.;
2483 const Double_t kthird = 1./3.;
2484 const Double_t khalf = 0.5;
2485 const Double_t kec = 2.9979251e-4;
2487 const Double_t kpisqua = 9.86960440109;
2488 const Int_t kix = 0;
2489 const Int_t kiy = 1;
2490 const Int_t kiz = 2;
2491 const Int_t kipx = 3;
2492 const Int_t kipy = 4;
2493 const Int_t kipz = 5;
2496 /// *. ------------------------------------------------------------------
2498 /// * this constant is for units cm,gev/c and kgauss
2502 for(Int_t j = 0; j < 7; j++)
2505 Double_t pinv = kec * charge / vect[6];
2513 if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
2514 ///cmodif: call gufld(vout,f) changed into:
2515 TGeoGlobalMagField::Instance()->Field(vout,f);
2518 /// * start of integration
2531 secxs[0] = (b * f[2] - c * f[1]) * ph2;
2532 secys[0] = (c * f[0] - a * f[2]) * ph2;
2533 seczs[0] = (a * f[1] - b * f[0]) * ph2;
2534 ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
2535 if (ang2 > kpisqua) break;
2537 dxt = h2 * a + h4 * secxs[0];
2538 dyt = h2 * b + h4 * secys[0];
2539 dzt = h2 * c + h4 * seczs[0];
2544 /// * second intermediate point
2547 est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
2549 if (ncut++ > maxcut) break;
2558 ///cmodif: call gufld(xyzt,f) changed into:
2559 TGeoGlobalMagField::Instance()->Field(xyzt,f);
2565 secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
2566 secys[1] = (ct * f[0] - at * f[2]) * ph2;
2567 seczs[1] = (at * f[1] - bt * f[0]) * ph2;
2571 secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
2572 secys[2] = (ct * f[0] - at * f[2]) * ph2;
2573 seczs[2] = (at * f[1] - bt * f[0]) * ph2;
2574 dxt = h * (a + secxs[2]);
2575 dyt = h * (b + secys[2]);
2576 dzt = h * (c + seczs[2]);
2580 at = a + 2.*secxs[2];
2581 bt = b + 2.*secys[2];
2582 ct = c + 2.*seczs[2];
2584 est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
2585 if (est > 2.*TMath::Abs(h)) {
2586 if (ncut++ > maxcut) break;
2595 ///cmodif: call gufld(xyzt,f) changed into:
2596 TGeoGlobalMagField::Instance()->Field(xyzt,f);
2598 z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
2599 y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
2600 x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
2602 secxs[3] = (bt*f[2] - ct*f[1])* ph2;
2603 secys[3] = (ct*f[0] - at*f[2])* ph2;
2604 seczs[3] = (at*f[1] - bt*f[0])* ph2;
2605 a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
2606 b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
2607 c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
2609 est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
2610 + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
2611 + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
2613 if (est > kdlt && TMath::Abs(h) > 1.e-4) {
2614 if (ncut++ > maxcut) break;
2620 /// * if too many iterations, go to helix
2621 if (iter++ > maxit) break;
2626 cba = 1./ TMath::Sqrt(a*a + b*b + c*c);
2634 if (step < 0.) rest = -rest;
2635 if (rest < 1.e-5*TMath::Abs(step)) return;
2639 /// angle too big, use helix
2644 f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
2653 hxp[0] = f2*vect[kipz] - f3*vect[kipy];
2654 hxp[1] = f3*vect[kipx] - f1*vect[kipz];
2655 hxp[2] = f1*vect[kipy] - f2*vect[kipx];
2657 hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
2660 sint = TMath::Sin(tet);
2661 cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
2665 g3 = (tet-sint) * hp*rho1;
2670 vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
2671 vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
2672 vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
2674 vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
2675 vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
2676 vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
2681 ///______________________________________________________________________________
2684 Bool_t AliHLTMUONFullTracker::SelectFront()
2686 /// Select the front trackseg connected with back trackseg
2688 Cluster clus1,clus2;
2689 Int_t minIndex=0,maxIndex=0;
2690 Int_t minCh=0,maxCh=0;
2691 Int_t ifronttrackseg=0;
2693 Float_t xf,yf,zf,thetaDev;
2694 Float_t p,spx,spy,spz,px,py,pz,sx,sy,sz,x,y,z;
2696 Float_t dist,mindist;
2697 Int_t frontsegIndex ;
2699 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
2701 if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
2703 /// if(fBackTrackSeg[ibacktrackseg].fIndex[2]==-1 || fBackTrackSeg[ibacktrackseg].fIndex[3]==-1) continue;
2705 ifronttrackseg = fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]-1];
2708 HLTImportant("AliHLTMUONFullTracker::SelectFront()--Begin\n\n");
2709 HLTImportant("\tbacktrack : %d is connected with : %d front tracks\n",
2710 ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg]);
2713 maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
2714 maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
2716 minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
2717 minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
2721 clus2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
2722 clus2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
2723 clus2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
2725 clus1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
2726 clus1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
2727 clus1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
2729 zf = 0.5*(AliMUONConstants::DefaultChamberZ(4) + AliMUONConstants::DefaultChamberZ(5));
2730 thetaDev= (1/zf)*(clus1.fY*clus2.fZ - clus2.fY*clus1.fZ)/(clus2.fZ - clus1.fZ);
2731 xf = clus1.fX*zf/clus1.fZ;
2732 yf = clus2.fY - ((clus2.fY - clus1.fY)*(clus2.fZ-zf))/(clus2.fZ - clus1.fZ);
2733 p = 3.0*0.3/thetaDev;
2734 charge = (p>0)?-1:+1;
2738 spz = sqrt((p*p-(spx*spx + spy*spy)));
2742 sx = clus1.fX; sy = clus1.fY; sz = clus1.fZ;
2745 for( Int_t iconnected=0;iconnected<fNofConnectedfrontTrackSeg[ibacktrackseg];iconnected++){
2747 ifronttrackseg = fBackToFront[ibacktrackseg][iconnected];
2749 minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2750 minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2751 clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2752 clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2753 clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2755 x = sx; y = sy; z = sz;
2756 px = spx; py = spy; pz = spz;
2757 PropagateTracks(charge,px,py,pz,x,y,z,clus1.fZ);
2759 dist = sqrt((clus1.fX-x)*(clus1.fX-x) +
2760 (clus1.fY-y)*(clus1.fY-y));
2763 frontsegIndex = ifronttrackseg;
2765 }///for loop on all connected front track segs
2767 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2768 ///have to check later
2769 if(frontsegIndex==-1) continue;
2771 fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]++] = frontsegIndex;
2773 /// fTrackParam[ibacktrackseg] = trackParam;
2776 HLTImportant("\tbacktrack : %d is connected with : %d front tracks\n",
2777 ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg]);
2778 HLTImportant("AliHLTMUONFullTracker::SelectFront()--End\n\n");
2782 }///backtrackSeg loop
2788 ///__________________________________________________________________________
2790 Bool_t AliHLTMUONFullTracker::KalmanChi2Test()
2793 /// Kalman Chi2 test for trach segments selection
2795 Cluster clus1,clus2;
2796 Int_t minIndex=0,maxIndex=0;
2797 Int_t minCh=0,maxCh=0;
2798 Int_t ifronttrackseg=0;
2799 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
2801 if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
2803 /// if(fBackTrackSeg[ibacktrackseg].fIndex[2]==-1 || fBackTrackSeg[ibacktrackseg].fIndex[3]==-1) continue;
2805 ifronttrackseg = fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]-1];
2808 HLTImportant("AliHLTMUONFullTracker::KalmanChi2Test()--Begin\n\n");
2809 HLTImportant("\tbacktrack : %d is connected with : %d front tracks, front track index %d\n",
2810 ibacktrackseg,fNofConnectedfrontTrackSeg[ibacktrackseg],ifronttrackseg);
2812 maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
2813 maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
2815 minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
2816 minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
2819 clus2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
2820 clus2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
2821 clus2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
2822 clus2.fErrX2 = 0.020736;
2823 clus2.fErrY2 = 0.000100;
2825 clus1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
2826 clus1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
2827 clus1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
2828 clus1.fErrX2 = 0.020736;
2829 clus1.fErrY2 = 0.000100;
2832 AliMUONTrackParam trackParam;
2833 Double_t dZ = clus2.fZ - clus1.fZ;
2834 trackParam.SetNonBendingCoor(clus2.fX);
2835 trackParam.SetBendingCoor(clus2.fY);
2836 trackParam.SetZ(clus2.fZ);
2837 trackParam.SetNonBendingSlope((clus2.fX - clus1.fX) / dZ);
2838 trackParam.SetBendingSlope((clus2.fY - clus1.fY) / dZ);
2839 Double_t bendingImpact = clus1.fY - clus1.fZ * trackParam.GetBendingSlope();
2840 Double_t inverseBendingMomentum = 1.
2841 / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
2842 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
2845 HLTImportant("\t\tCh%d : (%f,%f,%f)\n",maxCh,clus2.fX,clus2.fY,clus2.fZ);
2846 HLTImportant("\t\tCh%d : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2847 HLTImportant("\t\tFor minCh : %d, GetBeMom : %lf\n",minCh,trackParam.GetInverseBendingMomentum());
2850 /// trackParam->SetClusterPtr(clus[8]);
2852 /// => Reset track parameter covariances at last clus (as if the other cluss did not exist)
2853 TMatrixD lastParamCov(5,5);
2854 lastParamCov.Zero();
2855 lastParamCov(0,0) = clus1.fErrX2;
2856 lastParamCov(1,1) = 100. * ( clus2.fErrX2 + clus1.fErrX2 ) / dZ / dZ;
2857 lastParamCov(2,2) = clus1.fErrY2;
2858 lastParamCov(3,3) = 100. * ( clus2.fErrY2 + clus1.fErrY2 ) / dZ / dZ;
2859 lastParamCov(4,4) = ((10.0*10.0 + (clus2.fZ * clus2.fZ * clus1.fErrY2 +
2860 clus1.fZ * clus1.fZ * clus2.fErrY2)
2861 / dZ / dZ) /bendingImpact / bendingImpact +
2862 0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum;
2863 trackParam.SetCovariances(lastParamCov);
2865 AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2866 /// AliMUONTrackExtrap::ExtrapToZCov(trackParam, AliMUONConstants::DefaultChamberZ(minCh),kTRUE);
2867 ///AliMUONTrackExtrap::ExtrapToZCov(&trackParam, clus1.fZ,kTRUE);
2869 LinearExtrapToZ(&trackParam, clus1.fZ);
2871 trackParam.SetTrackChi2(0.);
2873 Double_t chi2 = 0.0;
2875 chi2 = KalmanFilter(trackParam,&clus1);
2877 if( chi2 > 1.0e9 /* is order to check TMath::AreEqualAbs(chi2,,1.0e-5)*/ ) {
2878 HLTWarning("Kalman Chi2 calculation cannot be completed...skipping slat track %d (c.p = 1)",ibacktrackseg);
2879 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2884 HLTImportant("\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\n",minCh,chi2,trackParam.GetInverseBendingMomentum());
2885 /// TMatrixD para2(trackParam->GetParameters());
2889 AliMUONTrackParam extrapTrackParamAtCluster2,minChi2Param;
2890 Double_t chi2OfCluster;
2892 Double_t minChi2=1000000.0;
2893 Int_t frontsegIndex = -1;
2894 extrapTrackParamAtCluster2 = trackParam ;
2895 minChi2Param = trackParam ;
2897 for( Int_t iconnected=0;iconnected<fNofConnectedfrontTrackSeg[ibacktrackseg];iconnected++){
2899 ifronttrackseg = fBackToFront[ibacktrackseg][iconnected];
2900 AliMUONTrackParam extrapTrackParam;
2901 trackParam = extrapTrackParamAtCluster2;
2903 minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2904 minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
2905 clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2906 clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2907 clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2908 clus1.fErrX2 = 0.020736;
2909 clus1.fErrY2 = 0.000100;
2912 AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2913 trackParam.ResetPropagator();
2914 AliMUONTrackExtrap::ExtrapToZCov(&trackParam,clus1.fZ, kTRUE);
2916 tryonefast = TryOneClusterFast(trackParam, &clus1);
2917 ///if (!tryonefast) continue;
2919 /// try to add the current cluster accuratly
2920 chi2OfCluster = TryOneCluster(trackParam, &clus1 , extrapTrackParam,kTRUE);
2922 extrapTrackParam.SetExtrapParameters(extrapTrackParam.GetParameters());
2923 extrapTrackParam.SetExtrapCovariances(extrapTrackParam.GetCovariances());
2925 chi2 = KalmanFilter(extrapTrackParam,&clus1);
2926 if( chi2 > 1.0e9 /* is order to check TMath::AreEqualAbs(chi2,,1.0e-5)*/ ) {
2927 HLTWarning("Kalman Chi2 calculation cannot be completed...skipping slat track %d (c.p = 2)",ibacktrackseg);
2928 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2933 minChi2Param = extrapTrackParam;
2934 frontsegIndex = ifronttrackseg;
2936 }///for loop on all connected front track segs
2938 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2939 ///have to check later
2940 if(frontsegIndex==-1) continue;
2942 fBackToFront[ibacktrackseg][fNofConnectedfrontTrackSeg[ibacktrackseg]++] = frontsegIndex;
2943 ifronttrackseg = frontsegIndex;
2945 trackParam = minChi2Param ;
2948 HLTImportant("\t\t\tCh%d : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2949 HLTImportant("\t\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\t",minCh,chi2,trackParam.GetInverseBendingMomentum());
2950 HLTImportant(Form("Pt : %lf\n",TMath::Sqrt(trackParam.Px()*trackParam.Px() +
2951 trackParam.Py()*trackParam.Py())));
2955 minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2956 minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
2957 clus1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
2958 clus1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
2959 clus1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
2960 clus1.fErrX2 = 0.020736;
2961 clus1.fErrY2 = 0.000100;
2963 AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(minCh),1.);
2964 trackParam.ResetPropagator();
2965 ///AliMUONTrackExtrap::ExtrapToZCov(&trackParam,clus1.fZ, kTRUE);
2966 LinearExtrapToZ(&trackParam, clus1.fZ);
2968 tryonefast = TryOneClusterFast(trackParam, &clus1);
2969 ///if (!tryonefast) continue;
2971 chi2OfCluster = TryOneCluster(trackParam, &clus1 , extrapTrackParamAtCluster2,kTRUE);
2973 extrapTrackParamAtCluster2.SetExtrapParameters(extrapTrackParamAtCluster2.GetParameters());
2974 extrapTrackParamAtCluster2.SetExtrapCovariances(extrapTrackParamAtCluster2.GetCovariances());
2976 chi2 = KalmanFilter(extrapTrackParamAtCluster2,&clus1);
2977 if(chi2 > 1.0e9 /* is order to check TMath::AreEqualAbs(chi2,,1.0e-5)*/ ) {
2978 HLTWarning("Kalman Chi2 calculation cannot be completed...skipping slat track %d (c.p = 3)",ibacktrackseg);
2979 fNofConnectedfrontTrackSeg[ibacktrackseg] = 0;
2983 trackParam = extrapTrackParamAtCluster2;
2986 HLTImportant("\t\tCh%d : (%f,%f,%f)\n",minCh,clus1.fX,clus1.fY,clus1.fZ);
2987 HLTImportant("\t\tFor minCh : %d, Chi2 = %lf, GetBeMom : %lf\t",minCh,chi2,trackParam.GetInverseBendingMomentum());
2988 HLTImportant(Form("Pt : %lf\n",TMath::Sqrt(extrapTrackParamAtCluster2.Px()*extrapTrackParamAtCluster2.Px() +
2989 extrapTrackParamAtCluster2.Py()*extrapTrackParamAtCluster2.Py())));
2991 HLTImportant("AliHLTMUONFullTracker::KalmanChi2Test()--End\n\n");
2993 ///AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.);
2994 //trackParam.SetInverseBendingMomentum(trackParam.GetCharge());
2995 fTrackParam[ibacktrackseg] = trackParam;
3006 ///__________________________________________________________________________
3009 Double_t AliHLTMUONFullTracker::BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ)
3011 //// Returns the mean total momentum energy loss of muon with total momentum='pTotal'
3012 //// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
3013 Double_t muMass = 0.105658369; /// GeV
3014 Double_t eMass = 0.510998918e-3; /// GeV
3015 Double_t k = 0.307075e-3; /// GeV.g^-1.cm^2
3016 Double_t i = 9.5e-9; /// mean exitation energy per atomic Z (GeV)
3017 Double_t p2=pTotal*pTotal;
3018 Double_t beta2=p2/(p2 + muMass*muMass);
3020 Double_t w = k * rho * pathLength * atomicZ / atomicA / beta2;
3022 if (beta2/(1-beta2)>3.5*3.5)
3023 return w * (log(2.*eMass*3.5/(i*atomicZ)) + 0.5*log(beta2/(1-beta2)) - beta2);
3025 return w * (log(2.*eMass*beta2/(1-beta2)/(i*atomicZ)) - beta2);
3029 ///__________________________________________________________________________
3031 Double_t AliHLTMUONFullTracker::EnergyLossFluctuation2(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ)
3033 //// Returns the total momentum energy loss fluctuation of muon with total momentum='pTotal'
3034 //// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
3035 Double_t muMass = 0.105658369; /// GeV
3036 ///Double_t eMass = 0.510998918e-3; /// GeV
3037 Double_t k = 0.307075e-3; /// GeV.g^-1.cm^2
3038 Double_t p2=pTotal*pTotal;
3039 Double_t beta2=p2/(p2 + muMass*muMass);
3041 Double_t fwhm = 2. * k * rho * pathLength * atomicZ / atomicA / beta2; /// FWHM of the energy loss Landau distribution
3042 Double_t sigma2 = fwhm * fwhm / (8.*log(2.)); /// gaussian: fwmh = 2 * srqt(2*ln(2)) * sigma (i.e. fwmh = 2.35 * sigma)
3044 ///sigma2 = k * rho * pathLength * atomicZ / atomicA * eMass; /// sigma2 of the energy loss gaussian distribution
3050 ///__________________________________________________________________________
3052 void AliHLTMUONFullTracker::LinearExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd)
3054 //// Track parameters (and their covariances if any) linearly extrapolated to the plane at "zEnd".
3055 //// On return, results from the extrapolation are updated in trackParam.
3057 if ( TMath::AreEqualAbs(trackParam->GetZ(),zEnd,1.0e-5) ) return; /// nothing to be done if same z
3059 Double_t dZ = zEnd - trackParam->GetZ();
3060 trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + trackParam->GetNonBendingSlope() * dZ);
3061 trackParam->SetBendingCoor(trackParam->GetBendingCoor() + trackParam->GetBendingSlope() * dZ);
3062 trackParam->SetZ(zEnd);
3065 void AliHLTMUONFullTracker::CorrectELossEffectInAbsorber(AliMUONTrackParam* param, Double_t eLoss)
3067 /// Energy loss correction
3068 Double_t nonBendingSlope = param->GetNonBendingSlope();
3069 Double_t bendingSlope = param->GetBendingSlope();
3070 param->SetInverseBendingMomentum(param->GetCharge() / (param->P() + eLoss) *
3071 TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) /
3072 TMath::Sqrt(1.0 + bendingSlope*bendingSlope));
3075 ///__________________________________________________________________________
3077 void AliHLTMUONFullTracker::Cov2CovP(const TMatrixD ¶m, TMatrixD &cov)
3079 //// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, SlopeX, Y, SlopeY, q*PTot)
3080 //// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
3082 /// charge * total momentum
3083 Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
3084 TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
3086 /// Jacobian of the opposite transformation
3087 TMatrixD jacob(5,5);
3089 jacob(4,1) = qPTot * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
3090 jacob(4,3) = - qPTot * param(1,0) * param(1,0) * param(3,0) /
3091 (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
3092 jacob(4,4) = - qPTot / param(4,0);
3094 /// compute covariances in new coordinate system
3096 tmp.MultT(cov,jacob);
3097 cov.Mult(jacob,tmp);
3101 ///__________________________________________________________________________
3103 void AliHLTMUONFullTracker::CovP2Cov(const TMatrixD ¶m, TMatrixD &covP)
3105 //// change coordinate system: (X, SlopeX, Y, SlopeY, q*PTot) -> (X, SlopeX, Y, SlopeY, q/Pyz)
3106 //// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
3108 /// charge * total momentum
3109 Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
3110 TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
3112 /// Jacobian of the transformation
3113 TMatrixD jacob(5,5);
3115 jacob(4,1) = param(4,0) * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
3116 jacob(4,3) = - param(4,0) * param(1,0) * param(1,0) * param(3,0) /
3117 (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
3118 jacob(4,4) = - param(4,0) / qPTot;
3120 /// compute covariances in new coordinate system
3122 tmp.MultT(covP,jacob);
3123 covP.Mult(jacob,tmp);
3127 ///__________________________________________________________________________
3130 void AliHLTMUONFullTracker::CorrectMCSEffectInAbsorber(AliMUONTrackParam* param,
3131 Double_t xVtx, Double_t yVtx, Double_t zVtx,
3132 Double_t absZBeg, Double_t f1, Double_t f2)
3134 //// Correct parameters and corresponding covariances using Branson correction
3135 //// - input param are parameters and covariances at the end of absorber
3136 //// - output param are parameters and covariances at vertex
3137 //// Absorber correction parameters are supposed to be calculated at the current track z-position
3139 /// Position of the Branson plane (spectro. (z<0))
3140 Double_t zB = (f1>0.) ? absZBeg - f2/f1 : 0.;
3142 LinearExtrapToZ(param,zB);
3144 /// compute track parameters at vertex
3145 TMatrixD newParam(5,1);
3147 newParam(0,0) = xVtx;
3148 newParam(1,0) = (param->GetNonBendingCoor() - xVtx) / (zB - zVtx);
3149 newParam(2,0) = yVtx;
3150 newParam(3,0) = (param->GetBendingCoor() - yVtx) / (zB - zVtx);
3151 newParam(4,0) = param->GetCharge() / param->P() *
3152 TMath::Sqrt(1.0 + newParam(1,0)*newParam(1,0) + newParam(3,0)*newParam(3,0)) /
3153 TMath::Sqrt(1.0 + newParam(3,0)*newParam(3,0));
3155 TMatrixD paramCovP(5,5);
3157 /// Get covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
3158 paramCovP.Use(param->GetCovariances());
3160 Cov2CovP(param->GetParameters(),paramCovP);
3162 /// Get the covariance matrix in the (XVtx, X, YVtx, Y, q*PTot) coordinate system
3163 /// TMatrixD paramCovVtx(5,5);
3164 Double_t element44 = paramCovP(4,4);
3166 paramCovP(4,4) = element44;
3168 CovP2Cov(newParam,paramCovP);
3170 /// Set parameters and covariances at vertex
3171 param->SetParameters(newParam);
3173 param->SetCovariances(paramCovP);
3176 ///__________________________________________________________________________
3178 Bool_t AliHLTMUONFullTracker::ExtrapolateToOrigin()
3180 /// Extrapolation to origin through absorber
3182 Int_t minIndex=0,maxIndex=0;
3183 Int_t minCh=0,maxCh=0;
3184 Int_t ifronttrackseg = -1;
3185 AliMUONTrackParam trackP;
3186 Double_t bSlope, nbSlope;
3187 AliHLTMUONRecHitStruct p1,p2,pSeg1,pSeg2;
3188 Double_t pyz = -1.0;
3189 TVector3 v1,v2,v3,v4;
3190 Double_t eLoss1,eLoss2,eLoss3;
3192 Double_t zE,zB,dzE,dzB;
3194 Double_t f0Sum,f1Sum,f2Sum;
3195 Double_t fXVertex=0.0,fYVertex=0.0,fZVertex=0.0;
3198 for( Int_t ibacktrackseg=0;ibacktrackseg<fNofbackTrackSeg;ibacktrackseg++){
3200 if(fNofConnectedfrontTrackSeg[ibacktrackseg]<=0) continue;
3202 maxIndex = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?3:2;
3203 maxCh = (fBackTrackSeg[ibacktrackseg].fIndex[3]!=-1)?9:8;
3205 minIndex = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?0:1;
3206 minCh = (fBackTrackSeg[ibacktrackseg].fIndex[0]!=-1)?6:7;
3208 p2.fX = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fX ;
3209 p2.fY = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fY ;
3210 p2.fZ = fChPoint[maxCh][fBackTrackSeg[ibacktrackseg].fIndex[maxIndex]]->fZ ;
3212 p1.fX = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fX ;
3213 p1.fY = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fY ;
3214 p1.fZ = fChPoint[minCh][fBackTrackSeg[ibacktrackseg].fIndex[minIndex]]->fZ ;
3216 thetaDev= (p1.fY*p2.fZ - p2.fY*p1.fZ)/(p2.fZ - p1.fZ);
3218 Sub(&p2,&p1,&pSeg1);
3220 ifronttrackseg = fBackToFront[ibacktrackseg][0];
3222 maxIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
3223 maxCh = (fFrontTrackSeg[ifronttrackseg].fIndex[3]!=-1)?3:2;
3225 minIndex = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
3226 minCh = (fFrontTrackSeg[ifronttrackseg].fIndex[0]!=-1)?0:1;
3228 p2.fX = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fX ;
3229 p2.fY = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fY ;
3230 p2.fZ = fChPoint[maxCh][fFrontTrackSeg[ifronttrackseg].fIndex[maxIndex]]->fZ ;
3232 p1.fX = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fX ;
3233 p1.fY = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fY ;
3234 p1.fZ = fChPoint[minCh][fFrontTrackSeg[ifronttrackseg].fIndex[minIndex]]->fZ ;
3236 Sub(&p2,&p1,&pSeg2);
3239 pyz = (3.0*0.3/TMath::Sin(Angle(&pSeg1,&pSeg2)));/// * sqrt(x3*x3 + y3*y3)/z3 ;
3241 pyz = -(3.0*0.3/TMath::Sin(Angle(&pSeg1,&pSeg2)));/// * sqrt(x3*x3 + y3*y3)/z3 ;
3243 nbSlope = (p2.fX - p1.fX)/(p2.fZ - p1.fZ);
3244 bSlope = (p2.fY - p1.fY)/(p2.fZ - p1.fZ);
3247 trackP.SetNonBendingCoor(p1.fX);
3248 trackP.SetNonBendingSlope(nbSlope);
3249 trackP.SetBendingCoor(p1.fY);
3250 trackP.SetBendingSlope(bSlope);
3251 trackP.SetInverseBendingMomentum(1.0/pyz) ;
3256 if(not fFastTracking){
3257 trackP = fTrackParam[ibacktrackseg] ;
3260 LinearExtrapToZ(&trackP,fgkAbsoedge[3]);
3261 v4.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
3262 LinearExtrapToZ(&trackP,fgkAbsoedge[2]);
3263 v3.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
3264 LinearExtrapToZ(&trackP,fgkAbsoedge[1]);
3265 v2.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
3266 LinearExtrapToZ(&trackP,fgkAbsoedge[0]);
3267 v1.SetXYZ(trackP.GetNonBendingCoor(),trackP.GetBendingCoor(),trackP.GetZ());
3269 eLoss1 = BetheBloch(trackP.P(), (v4-v3).Mag(), fgkRho[2], fgkAtomicA[2], fgkAtomicZ[2]);
3270 eLoss2 = BetheBloch(trackP.P(), (v3-v2).Mag(), fgkRho[1], fgkAtomicA[1], fgkAtomicZ[1]);
3271 eLoss3 = BetheBloch(trackP.P(), (v2-v1).Mag(), fgkRho[0], fgkAtomicA[0], fgkAtomicZ[0]);
3273 /// sigmaELoss1 = EnergyLossFluctuation2(trackP.P(), (v4-v3).Mag(), rho[2], atomicA[2], atomicZ[2]);
3274 /// sigmaELoss2 = EnergyLossFluctuation2(trackP.P(), (v3-v2).Mag(), rho[1], atomicA[1], atomicZ[1]);
3275 /// sigmaELoss3 = EnergyLossFluctuation2(trackP.P(), (v2-v1).Mag(), rho[0], atomicA[0], atomicZ[0]);
3277 /// eDiff = totELoss-(eLoss1+eLoss2+eLoss3);
3278 /// sigmaELossDiff = sigmaTotELoss ;///- (sigmaELoss1+sigmaELoss2+sigmaELoss3);
3280 /// CorrectELossEffectInAbsorber(&trackP, 0.5*(eLoss1+eLoss2+eLoss3), 0.5*(sigmaELoss1+sigmaELoss2+sigmaELoss3));
3283 ///CorrectELossEffectInAbsorber(&trackP, totELoss,sigmaTotELoss);
3285 CorrectELossEffectInAbsorber(&trackP, 0.7*(eLoss1+eLoss2+eLoss3));
3287 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3288 f0Sum = 0.0; f1Sum = 0.0; f2Sum = 0.0;
3290 b = (v4.Z()-v1.Z())/((v4-v1).Mag());
3293 zE = b*((v2-v1).Mag()) + zB;
3297 f0 = ((v2-v1).Mag())/fgkRadLen[0];
3298 f1 = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[0] /2.;
3299 f2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[0] / 3.;
3306 zE = b*((v3-v2).Mag()) + zB;
3310 f0 = ((v3-v2).Mag())/fgkRadLen[1];
3311 f1 = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[1] /2.;
3312 f2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[1] / 3.;
3319 zE = b*((v4-v3).Mag()) + zB;
3323 f0 = ((v4-v3).Mag())/fgkRadLen[2];
3324 f1 = (dzE*dzE - dzB*dzB) / b / b / fgkRadLen[2] /2.;
3325 f2 = (dzE*dzE*dzE - dzB*dzB*dzB) / b / b / b / fgkRadLen[2] / 3.;
3332 ///AddMCSEffectInAbsorber(&trackP,(v4-v1).Mag(),f0Sum,f1Sum,f2Sum);
3334 CorrectMCSEffectInAbsorber(&trackP,fXVertex,fYVertex, fZVertex,AliMUONConstants::AbsZBeg(),f1Sum,f2Sum);
3335 CorrectELossEffectInAbsorber(&trackP, 0.5*(eLoss1+eLoss2+eLoss3));
3338 ///AliMUONTrackExtrap::ExtrapToVertex(&trackP, 0., 0., 0., 0., 0.);
3340 trackP.SetNonBendingCoor(p1.fX);
3341 trackP.SetBendingCoor(p1.fY);
3342 LinearExtrapToZ(&trackP, 0.0);
3344 fTrackParam[ibacktrackseg] = trackP;
3346 }///backtrackseg for loop