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