MUON/AliMUONTrack.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------------
  19 // Class AliMUONTrack
  20 //-------------------
  21 // Reconstructed track in ALICE dimuon spectrometer
  22 //-----------------------------------------------------------------------------
  23
  24 #include "AliMUONTrack.h"
  25
  26 #include "AliMUONTrackParam.h"
  27 #include "AliMUONHitForRec.h"
  28 #include "AliMUONObjectPair.h"
  29 #include "AliMUONConstants.h"
  30 #include "AliMUONTrackExtrap.h"
  31
  32 #include "AliLog.h"
  33
  34 #include <TMath.h>
  35 #include <TMatrixD.h>
  36
  37 #include <Riostream.h>
  38
  39 /// \cond CLASSIMP
  40 ClassImp(AliMUONTrack) // Class implementation in ROOT context
  41 /// \endcond
  42
  43 //__________________________________________________________________________
  44 AliMUONTrack::AliMUONTrack()
  45   : TObject(),
  46     fTrackParamAtVertex(),
  47     fTrackParamAtHit(0x0),
  48     fHitForRecAtHit(0x0),
  49     fNTrackHits(0),
  50     fFitWithVertex(kFALSE),
  51     fVertex(0x0),
  52     fFitWithMCS(kFALSE),
  53     fHitWeightsNonBending(0x0),
  54     fHitWeightsBending(0x0),
  55     fGlobalChi2(-1.),
  56     fImproved(kFALSE),
  57     fMatchTrigger(-1),
  58     floTrgNum(-1),
  59     fChi2MatchTrigger(0.),
  60     fTrackID(0),
  61     fHitsPatternInTrigCh(0),
  62     fLocalTrigger(234)
  63 {
  64   /// Default constructor
  65 }
  66
  67   //__________________________________________________________________________
  68 AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment)
  69   : TObject(),
  70     fTrackParamAtVertex(),
  71     fTrackParamAtHit(0x0),
  72     fHitForRecAtHit(0x0),
  73     fNTrackHits(0),
  74     fFitWithVertex(kFALSE),
  75     fVertex(0x0),
  76     fFitWithMCS(kFALSE),
  77     fHitWeightsNonBending(0x0),
  78     fHitWeightsBending(0x0),
  79     fGlobalChi2(0.),
  80     fImproved(kFALSE),
  81     fMatchTrigger(-1),
  82     floTrgNum(-1),
  83     fChi2MatchTrigger(0.),
  84     fTrackID(0),
  85     fHitsPatternInTrigCh(0),
  86     fLocalTrigger(234)
  87 {
  88   /// Constructor from thw hitForRec's
  89
  90   fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
  91   fTrackParamAtHit->SetOwner(kTRUE);
  92   fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
  93   fHitForRecAtHit->SetOwner(kTRUE);
  94
  95   if (!segment) return; //AZ
  96
  97   // Pointers to hits from the segment
  98   AliMUONHitForRec* hit1 = (AliMUONHitForRec*) segment->First();
  99   AliMUONHitForRec* hit2 = (AliMUONHitForRec*) segment->Second();
 100
 101   // check sorting in -Z (spectro z<0)
 102   if (hit1->GetZ() < hit2->GetZ()) {
 103     hit1 = hit2;
 104     hit2 = (AliMUONHitForRec*) segment->First();
 105   }
 106
 107   // order the hits into the track according to the station the segment belong to
 108   //(the hit first attached is the one from which we will start the tracking procedure)
 109   if (hit1->GetChamberNumber() == 8) {
 110     AddTrackParamAtHit(0,hit1);
 111     AddTrackParamAtHit(0,hit2);
 112   } else {
 113     AddTrackParamAtHit(0,hit2);
 114     AddTrackParamAtHit(0,hit1);
 115   }
 116
 117   AliMUONTrackParam* trackParamAtFirstHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
 118   AliMUONHitForRec* firstHit = trackParamAtFirstHit->GetHitForRecPtr();
 119   AliMUONTrackParam* trackParamAtLastHit = (AliMUONTrackParam*) fTrackParamAtHit->Last();
 120   AliMUONHitForRec* lastHit = trackParamAtLastHit->GetHitForRecPtr();
 121
 122
 123   // Compute track parameters
 124   Double_t dZ = firstHit->GetZ() - lastHit->GetZ();
 125   // Non bending plane
 126   Double_t nonBendingCoor1 = firstHit->GetNonBendingCoor();
 127   Double_t nonBendingCoor2 = lastHit->GetNonBendingCoor();
 128   Double_t nonBendingSlope = (nonBendingCoor1 - nonBendingCoor2) / dZ;
 129   // Bending plane
 130   Double_t bendingCoor1 = firstHit->GetBendingCoor();
 131   Double_t bendingCoor2 = lastHit->GetBendingCoor();
 132   Double_t bendingSlope = (bendingCoor1 - bendingCoor2) / dZ;
 133   // Inverse bending momentum
 134   Double_t bendingImpact = bendingCoor1 - firstHit->GetZ() * bendingSlope;
 135   Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
 136
 137
 138   // Set track parameters at first hit
 139   trackParamAtFirstHit->SetNonBendingCoor(nonBendingCoor1);
 140   trackParamAtFirstHit->SetNonBendingSlope(nonBendingSlope);
 141   trackParamAtFirstHit->SetBendingCoor(bendingCoor1);
 142   trackParamAtFirstHit->SetBendingSlope(bendingSlope);
 143   trackParamAtFirstHit->SetInverseBendingMomentum(inverseBendingMomentum);
 144
 145
 146   // Set track parameters at last hit
 147   trackParamAtLastHit->SetNonBendingCoor(nonBendingCoor2);
 148   trackParamAtLastHit->SetNonBendingSlope(nonBendingSlope);
 149   trackParamAtLastHit->SetBendingCoor(bendingCoor2);
 150   trackParamAtLastHit->SetBendingSlope(bendingSlope);
 151   trackParamAtLastHit->SetInverseBendingMomentum(inverseBendingMomentum);
 152
 153
 154   // Compute and set track parameters covariances at first hit
 155   TMatrixD paramCov1(5,5);
 156   paramCov1.Zero();
 157   // Non bending plane
 158   paramCov1(0,0) = firstHit->GetNonBendingReso2();
 159   paramCov1(0,1) = firstHit->GetNonBendingReso2() / dZ;
 160   paramCov1(1,0) = paramCov1(0,1);
 161   paramCov1(1,1) = ( firstHit->GetNonBendingReso2() + lastHit->GetNonBendingReso2() ) / dZ / dZ;
 162   // Bending plane
 163   paramCov1(2,2) = firstHit->GetBendingReso2();
 164   paramCov1(2,3) = firstHit->GetBendingReso2() / dZ;
 165   paramCov1(3,2) = paramCov1(2,3);
 166   paramCov1(3,3) = ( firstHit->GetBendingReso2() + lastHit->GetBendingReso2() ) / dZ / dZ;
 167   // Inverse bending momentum (50% error)
 168   paramCov1(4,4) = 0.5*inverseBendingMomentum * 0.5*inverseBendingMomentum;
 169   // Set covariances
 170   trackParamAtFirstHit->SetCovariances(paramCov1);
 171
 172
 173   // Compute and set track parameters covariances at last hit (as if the first hit did not exist)
 174   TMatrixD paramCov2(5,5);
 175   paramCov2.Zero();
 176   // Non bending plane
 177   paramCov2(0,0) = paramCov1(0,0);
 178   paramCov2(1,1) = 100.*paramCov1(1,1);
 179   // Bending plane
 180   paramCov2(2,2) = paramCov1(2,2);
 181   paramCov2(3,3) = 100.*paramCov1(3,3);
 182   // Inverse bending momentum
 183   paramCov2(4,4) = paramCov1(4,4);
 184   // Set covariances
 185   trackParamAtLastHit->SetCovariances(paramCov2);
 186
 187
 188   // Flag first hit as being removable
 189   trackParamAtFirstHit->SetRemovable(kTRUE);
 190
 191   // Flag last hit as being removable
 192   trackParamAtLastHit->SetRemovable(kTRUE);
 193
 194 }
 195
 196   //__________________________________________________________________________
 197 AliMUONTrack::AliMUONTrack (const AliMUONTrack& track)
 198   : TObject(track),
 199     fTrackParamAtVertex(track.fTrackParamAtVertex),
 200     fTrackParamAtHit(0x0),
 201     fHitForRecAtHit(0x0),
 202     fNTrackHits(track.fNTrackHits),
 203     fFitWithVertex(track.fFitWithVertex),
 204     fVertex(0x0),
 205     fFitWithMCS(track.fFitWithMCS),
 206     fHitWeightsNonBending(0x0),
 207     fHitWeightsBending(0x0),
 208     fGlobalChi2(track.fGlobalChi2),
 209     fImproved(track.fImproved),
 210     fMatchTrigger(track.fMatchTrigger),
 211     floTrgNum(track.floTrgNum),
 212     fChi2MatchTrigger(track.fChi2MatchTrigger),
 213     fTrackID(track.fTrackID),
 214     fHitsPatternInTrigCh(track.fHitsPatternInTrigCh),
 215     fLocalTrigger(track.fLocalTrigger)
 216 {
 217   ///copy constructor
 218   Int_t maxIndex = 0;
 219
 220   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
 221   if (track.fTrackParamAtHit) {
 222     maxIndex = (track.fTrackParamAtHit)->GetEntriesFast();
 223     fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",maxIndex);
 224     for (Int_t index = 0; index < maxIndex; index++) {
 225       new ((*fTrackParamAtHit)[index]) AliMUONTrackParam(*(AliMUONTrackParam*)track.fTrackParamAtHit->At(index));
 226     }
 227   }
 228
 229   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
 230   if (track.fHitForRecAtHit) {
 231     maxIndex = (track.fHitForRecAtHit)->GetEntriesFast();
 232     fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",maxIndex);
 233     for (Int_t index = 0; index < maxIndex; index++) {
 234       new ((*fHitForRecAtHit)[index]) AliMUONHitForRec(*(AliMUONHitForRec*)track.fHitForRecAtHit->At(index));
 235     }
 236   }
 237
 238   // copy vertex used during the tracking procedure if any
 239   if (track.fVertex) fVertex = new AliMUONHitForRec(*(track.fVertex));
 240
 241   // copy hit weights matrices if any
 242   if (track.fHitWeightsNonBending) fHitWeightsNonBending = new TMatrixD(*(track.fHitWeightsNonBending));
 243   if (track.fHitWeightsBending) fHitWeightsBending = new TMatrixD(*(track.fHitWeightsBending));
 244
 245 }
 246
 247   //__________________________________________________________________________
 248 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track)
 249 {
 250   /// Asignment operator
 251   // check assignement to self
 252   if (this == &track)
 253     return *this;
 254
 255   // base class assignement
 256   TObject::operator=(track);
 257
 258   fTrackParamAtVertex = track.fTrackParamAtVertex;
 259
 260   Int_t maxIndex = 0;
 261
 262   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
 263   if (track.fTrackParamAtHit) {
 264     if (fTrackParamAtHit) fTrackParamAtHit->Clear();
 265     else fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
 266     maxIndex = (track.fTrackParamAtHit)->GetEntriesFast();
 267     for (Int_t index = 0; index < maxIndex; index++) {
 268       new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()])
 269         AliMUONTrackParam(*(AliMUONTrackParam*)(track.fTrackParamAtHit)->At(index));
 270     }
 271   } else if (fTrackParamAtHit) {
 272     delete fTrackParamAtHit;
 273     fTrackParamAtHit = 0x0;
 274   }
 275
 276   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
 277   if (track.fHitForRecAtHit) {
 278     if (fHitForRecAtHit) fHitForRecAtHit->Clear();
 279     else fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
 280     maxIndex = (track.fHitForRecAtHit)->GetEntriesFast();
 281     for (Int_t index = 0; index < maxIndex; index++) {
 282       new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()])
 283         AliMUONHitForRec(*(AliMUONHitForRec*)(track.fHitForRecAtHit)->At(index));
 284     }
 285   } else if (fHitForRecAtHit) {
 286     delete fHitForRecAtHit;
 287     fHitForRecAtHit = 0x0;
 288   }
 289
 290   // copy vertex used during the tracking procedure if any.
 291   if (track.fVertex) {
 292     if (fVertex) *fVertex = *(track.fVertex);
 293     else fVertex = new AliMUONHitForRec(*(track.fVertex));
 294   } else if (fVertex) {
 295     delete fVertex;
 296     fVertex = 0x0;
 297   }
 298
 299   // copy hit weights matrix if any
 300   if (track.fHitWeightsNonBending) {
 301     if (fHitWeightsNonBending) {
 302       fHitWeightsNonBending->ResizeTo(*(track.fHitWeightsNonBending));
 303       *fHitWeightsNonBending = *(track.fHitWeightsNonBending);
 304     } else fHitWeightsNonBending = new TMatrixD(*(track.fHitWeightsNonBending));
 305   } else if (fHitWeightsNonBending) {
 306     delete fHitWeightsNonBending;
 307     fHitWeightsNonBending = 0x0;
 308   }
 309
 310   // copy hit weights matrix if any
 311   if (track.fHitWeightsBending) {
 312     if (fHitWeightsBending) {
 313       fHitWeightsBending->ResizeTo(*(track.fHitWeightsBending));
 314       *fHitWeightsBending = *(track.fHitWeightsBending);
 315     } else fHitWeightsBending = new TMatrixD(*(track.fHitWeightsBending));
 316   } else if (fHitWeightsBending) {
 317     delete fHitWeightsBending;
 318     fHitWeightsBending = 0x0;
 319   }
 320
 321   fNTrackHits         =  track.fNTrackHits;
 322   fFitWithVertex      =  track.fFitWithVertex;
 323   fFitWithMCS         =  track.fFitWithMCS;
 324   fGlobalChi2         =  track.fGlobalChi2;
 325   fImproved           =  track.fImproved;
 326   fMatchTrigger       =  track.fMatchTrigger;
 327   floTrgNum           =  track.floTrgNum;
 328   fChi2MatchTrigger   =  track.fChi2MatchTrigger;
 329   fTrackID            =  track.fTrackID;
 330   fHitsPatternInTrigCh = track.fHitsPatternInTrigCh;
 331   fLocalTrigger        = track.fLocalTrigger;
 332
 333   return *this;
 334 }
 335
 336   //__________________________________________________________________________
 337 AliMUONTrack::~AliMUONTrack()
 338 {
 339   /// Destructor
 340   delete fTrackParamAtHit;
 341   delete fHitForRecAtHit;
 342   delete fVertex;
 343   delete fHitWeightsNonBending;
 344   delete fHitWeightsBending;
 345 }
 346
 347   //__________________________________________________________________________
 348 void AliMUONTrack::Clear(Option_t* opt)
 349 {
 350   /// Clear arrays
 351   if ( fTrackParamAtHit ) fTrackParamAtHit->Clear(opt);
 352   if ( fHitForRecAtHit ) fHitForRecAtHit->Clear(opt);
 353   delete fVertex; fVertex = 0x0;
 354   delete fHitWeightsNonBending; fHitWeightsNonBending = 0x0;
 355   delete fHitWeightsBending; fHitWeightsBending = 0x0;
 356 }
 357
 358   //__________________________________________________________________________
 359 void AliMUONTrack::AddTrackParamAtHit(const AliMUONTrackParam *trackParam, AliMUONHitForRec *hitForRec)
 360 {
 361   /// Add TrackParamAtHit if "trackParam" != NULL
 362   /// else create empty TrackParamAtHit and set the z position to the one of "hitForRec" if any
 363   /// Update link to HitForRec if "hitForRec" != NULL
 364   if (!fTrackParamAtHit) {
 365     fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
 366     fNTrackHits = 0;
 367   }
 368   AliMUONTrackParam* trackParamAtHit;
 369   if (trackParam) {
 370     trackParamAtHit = new ((*fTrackParamAtHit)[fNTrackHits]) AliMUONTrackParam(*trackParam);
 371     if (hitForRec) {
 372       if (hitForRec->GetZ() != trackParam->GetZ())
 373         AliWarning("Added track parameters at a different z position than the one of the attached hit");
 374     }
 375   } else {
 376     trackParamAtHit = new ((*fTrackParamAtHit)[fNTrackHits]) AliMUONTrackParam();
 377     if (hitForRec) trackParamAtHit->SetZ(hitForRec->GetZ());
 378   }
 379   if (hitForRec) trackParamAtHit->SetHitForRecPtr(hitForRec);
 380   fNTrackHits++;
 381 }
 382
 383   //__________________________________________________________________________
 384 void AliMUONTrack::RemoveTrackParamAtHit(AliMUONTrackParam *trackParam)
 385 {
 386   /// Remove trackParam from the array of TrackParamAtHit
 387   if (!fTrackParamAtHit) {
 388     AliWarning("array fTrackParamAtHit does not exist");
 389     return;
 390   }
 391
 392   if (!fTrackParamAtHit->Remove(trackParam)) {
 393     AliWarning("object to remove does not exist in array fTrackParamAtHit");
 394     return;
 395   }
 396
 397   fTrackParamAtHit->Compress();
 398   fNTrackHits--;
 399 }
 400
 401   //__________________________________________________________________________
 402 void AliMUONTrack::AddHitForRecAtHit(const AliMUONHitForRec *hitForRec)
 403 {
 404   /// Add hitForRec to the array of hitForRec at hit
 405   if (!fHitForRecAtHit)
 406     fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
 407
 408   if (!hitForRec)
 409     AliFatal("AliMUONTrack::AddHitForRecAtHit: hitForRec == NULL");
 410
 411   new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) AliMUONHitForRec(*hitForRec);
 412 }
 413
 414   //__________________________________________________________________________
 415 void AliMUONTrack::UpdateTrackParamAtHit()
 416 {
 417   /// Update track parameters at each attached hit
 418
 419   if (fNTrackHits == 0) {
 420     AliWarning("no hit attached to the track");
 421     return;
 422   }
 423
 424   Double_t z;
 425   AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtHit->First();
 426   AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(startingTrackParam);
 427   while (trackParamAtHit) {
 428
 429     // save current z
 430     z = trackParamAtHit->GetZ();
 431
 432     // reset track parameters and their covariances
 433     trackParamAtHit->SetParameters(startingTrackParam->GetParameters());
 434     trackParamAtHit->SetZ(startingTrackParam->GetZ());
 435
 436     // extrapolation to the given z
 437     AliMUONTrackExtrap::ExtrapToZ(trackParamAtHit, z);
 438
 439     // prepare next step
 440     startingTrackParam = trackParamAtHit;
 441     trackParamAtHit = (AliMUONTrackParam*) (fTrackParamAtHit->After(trackParamAtHit));
 442   }
 443
 444 }
 445
 446   //__________________________________________________________________________
 447 void AliMUONTrack::UpdateCovTrackParamAtHit()
 448 {
 449   /// Update track parameters and their covariances at each attached hit
 450
 451   if (fNTrackHits == 0) {
 452     AliWarning("no hit attached to the track");
 453     return;
 454   }
 455
 456   Double_t z;
 457   AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtHit->First();
 458   AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(startingTrackParam);
 459   while (trackParamAtHit) {
 460
 461     // save current z
 462     z = trackParamAtHit->GetZ();
 463
 464     // reset track parameters and their covariances
 465     trackParamAtHit->SetParameters(startingTrackParam->GetParameters());
 466     trackParamAtHit->SetZ(startingTrackParam->GetZ());
 467     trackParamAtHit->SetCovariances(startingTrackParam->GetCovariances());
 468
 469     // extrapolation to the given z
 470     AliMUONTrackExtrap::ExtrapToZCov(trackParamAtHit, z);
 471
 472     // prepare next step
 473     startingTrackParam = trackParamAtHit;
 474     trackParamAtHit = (AliMUONTrackParam*) (fTrackParamAtHit->After(trackParamAtHit));
 475   }
 476
 477 }
 478
 479   //__________________________________________________________________________
 480 void AliMUONTrack::SetVertex(const AliMUONHitForRec* vertex)
 481 {
 482   /// Set the vertex used during the tracking procedure
 483   if (!fVertex) fVertex = new AliMUONHitForRec(*vertex);
 484   else *fVertex = *vertex;
 485 }
 486
 487
 488   //__________________________________________________________________________
 489 Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
 490 {
 491   /// Compute the removable hit contribution to the chi2 of the track
 492   /// accounting for multiple scattering or not according to the flag
 493   /// - Also recompute the weight matrices of the attached hits if accountForMCS=kTRUE
 494   /// - Assume that track parameters at each hit are corrects
 495   /// - Return kFALSE if computation failed
 496
 497   // Check hits (if the first one exist, assume that the other ones exit too!)
 498   AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
 499   if (!trackParamAtHit->GetHitForRecPtr()) {
 500     AliWarning("hit is missing");
 501     return kFALSE;
 502   }
 503
 504   if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects
 505
 506     // Compute MCS covariance matrix only once
 507     TMatrixD mcsCovariances(fNTrackHits,fNTrackHits);
 508     ComputeMCSCovariances(mcsCovariances);
 509
 510     // Make sure hit weights are consistent with following calculations
 511     if (!ComputeHitWeights(&mcsCovariances)) {
 512       AliWarning("cannot take into account the multiple scattering effects");
 513       return ComputeLocalChi2(kFALSE);
 514     }
 515
 516     // Compute chi2 of the track
 517     Double_t globalChi2 = ComputeGlobalChi2(kTRUE);
 518     if (globalChi2 < 0.) return kFALSE;
 519
 520     // Loop over removable hits and compute their local chi2
 521     AliMUONTrackParam* trackParamAtHit1;
 522     AliMUONHitForRec *hitForRec, *discardedHit;
 523     Int_t hitNumber1, hitNumber2, currentHitNumber1, currentHitNumber2;
 524     TMatrixD hitWeightsNB(fNTrackHits-1,fNTrackHits-1);
 525     TMatrixD hitWeightsB(fNTrackHits-1,fNTrackHits-1);
 526     Double_t *dX = new Double_t[fNTrackHits-1];
 527     Double_t *dY = new Double_t[fNTrackHits-1];
 528     Double_t globalChi2b;
 529     while (trackParamAtHit) {
 530
 531       discardedHit = trackParamAtHit->GetHitForRecPtr();
 532
 533       // Recompute hit weights without the current hit
 534       if (!ComputeHitWeights(hitWeightsNB, hitWeightsB, &mcsCovariances, discardedHit)) {
 535         AliWarning("cannot take into account the multiple scattering effects");
 536         ComputeLocalChi2(kFALSE);
 537       }
 538
 539       // Compute track chi2 without the current hit
 540       globalChi2b = 0.;
 541       currentHitNumber1 = 0;
 542       for (hitNumber1 = 0; hitNumber1 < fNTrackHits ; hitNumber1++) {
 543         trackParamAtHit1 = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1);
 544         hitForRec = trackParamAtHit1->GetHitForRecPtr();
 545
 546         if (hitForRec == discardedHit) continue;
 547
 548         // Compute and save residuals
 549         dX[currentHitNumber1] = hitForRec->GetNonBendingCoor() - trackParamAtHit1->GetNonBendingCoor();
 550         dY[currentHitNumber1] = hitForRec->GetBendingCoor() - trackParamAtHit1->GetBendingCoor();
 551
 552         currentHitNumber2 = 0;
 553         for (hitNumber2 = 0; hitNumber2 < hitNumber1; hitNumber2++) {
 554           hitForRec = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber2))->GetHitForRecPtr();
 555
 556           if (hitForRec == discardedHit) continue;
 557
 558           // Add contribution from covariances
 559           globalChi2b += (hitWeightsNB(currentHitNumber1, currentHitNumber2) +
 560                           hitWeightsNB(currentHitNumber2, currentHitNumber1)) * dX[currentHitNumber1] * dX[currentHitNumber2] +
 561                          (hitWeightsB(currentHitNumber1, currentHitNumber2) +
 562                           hitWeightsB(currentHitNumber2, currentHitNumber1)) * dY[currentHitNumber1] * dY[currentHitNumber2];
 563
 564           currentHitNumber2++;
 565         }
 566
 567         // Add contribution from variances
 568         globalChi2b += hitWeightsNB(currentHitNumber1, currentHitNumber1) * dX[currentHitNumber1] * dX[currentHitNumber1] +
 569                        hitWeightsB(currentHitNumber1, currentHitNumber1) * dY[currentHitNumber1] * dY[currentHitNumber1];
 570
 571         currentHitNumber1++;
 572       }
 573
 574       // Set local chi2
 575       trackParamAtHit->SetLocalChi2(globalChi2 - globalChi2b);
 576
 577       trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(trackParamAtHit);
 578     }
 579
 580     delete [] dX;
 581     delete [] dY;
 582
 583   } else { // without multiple scattering effects
 584
 585     AliMUONHitForRec *discardedHit;
 586     Double_t dX, dY;
 587     while (trackParamAtHit) {
 588
 589       discardedHit = trackParamAtHit->GetHitForRecPtr();
 590
 591       // Compute residuals
 592       dX = discardedHit->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
 593       dY = discardedHit->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
 594
 595       // Set local chi2
 596       trackParamAtHit->SetLocalChi2(dX * dX / discardedHit->GetNonBendingReso2() + dY * dY / discardedHit->GetBendingReso2());
 597
 598     trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(trackParamAtHit);
 599     }
 600
 601   }
 602
 603   return kTRUE;
 604
 605 }
 606
 607   //__________________________________________________________________________
 608 Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
 609 {
 610   /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
 611   /// - Assume that track parameters at each hit are corrects
 612   /// - Assume the hits weights matrices are corrects
 613   /// - Return negative value if chi2 computation failed
 614
 615   // Check hits (if the first one exist, assume that the other ones exit too!)
 616   AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
 617   if (!trackParamAtHit->GetHitForRecPtr()) {
 618     AliWarning("hit is missing");
 619     return -1.;
 620   }
 621
 622   Double_t chi2 = 0.;
 623
 624   if (accountForMCS) {
 625
 626     // Check the weight matrices. If weight matrices are not available compute chi2 without MCS
 627     if (!fHitWeightsNonBending || !fHitWeightsBending) {
 628       AliWarning("hit weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
 629       return ComputeGlobalChi2(kFALSE);
 630     }
 631     if (fHitWeightsNonBending->GetNrows() != fNTrackHits || fHitWeightsBending->GetNcols() != fNTrackHits) {
 632       AliWarning("hit weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
 633       return ComputeGlobalChi2(kFALSE);
 634     }
 635
 636     // Compute chi2
 637     AliMUONHitForRec *hitForRec;
 638     Double_t *dX = new Double_t[fNTrackHits];
 639     Double_t *dY = new Double_t[fNTrackHits];
 640     Int_t hitNumber1, hitNumber2;
 641     for (hitNumber1 = 0; hitNumber1 < fNTrackHits ; hitNumber1++) {
 642       trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1);
 643       hitForRec = trackParamAtHit->GetHitForRecPtr();
 644       dX[hitNumber1] = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
 645       dY[hitNumber1] = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
 646       for (hitNumber2 = 0; hitNumber2 < hitNumber1; hitNumber2++) {
 647         chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber2) + (*fHitWeightsNonBending)(hitNumber2, hitNumber1)) * dX[hitNumber1] * dX[hitNumber2] +
 648                 ((*fHitWeightsBending)(hitNumber1, hitNumber2) + (*fHitWeightsBending)(hitNumber2, hitNumber1)) * dY[hitNumber1] * dY[hitNumber2];
 649       }
 650       chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber1) * dX[hitNumber1] * dX[hitNumber1]) +
 651               ((*fHitWeightsBending)(hitNumber1, hitNumber1) * dY[hitNumber1] * dY[hitNumber1]);
 652     }
 653     delete [] dX;
 654     delete [] dY;
 655
 656   } else {
 657
 658     AliMUONHitForRec *hitForRec;
 659     Double_t dX, dY;
 660     for (Int_t hitNumber = 0; hitNumber < fNTrackHits ; hitNumber++) {
 661       trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
 662       hitForRec = trackParamAtHit->GetHitForRecPtr();
 663       dX = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
 664       dY = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
 665       chi2 += dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
 666     }
 667
 668   }
 669
 670   return chi2;
 671
 672 }
 673
 674   //__________________________________________________________________________
 675 Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD* mcsCovariances)
 676 {
 677   /// Compute the weight matrices of the attached hits, in non bending and bending direction,
 678   /// accounting for multiple scattering correlations and hits resolution
 679   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
 680   /// - Assume that track parameters at each hit are corrects
 681   /// - Return kFALSE if computation failed
 682
 683   // Alocate memory
 684   if (!fHitWeightsNonBending) fHitWeightsNonBending = new TMatrixD(fNTrackHits,fNTrackHits);
 685   if (!fHitWeightsBending) fHitWeightsBending = new TMatrixD(fNTrackHits,fNTrackHits);
 686
 687   // Check hits (if the first one exist, assume that the other ones exit too!)
 688   if (!((AliMUONTrackParam*) fTrackParamAtHit->First())->GetHitForRecPtr()) {
 689     AliWarning("hit is missing");
 690     fHitWeightsNonBending->ResizeTo(0,0);
 691     fHitWeightsBending->ResizeTo(0,0);
 692     return kFALSE;
 693   }
 694
 695   // Compute weights matrices
 696   if (!ComputeHitWeights(*fHitWeightsNonBending, *fHitWeightsBending, mcsCovariances)) return kFALSE;
 697
 698   return kTRUE;
 699
 700 }
 701
 702   //__________________________________________________________________________
 703 Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD& hitWeightsNB, TMatrixD& hitWeightsB, TMatrixD* mcsCovariances, AliMUONHitForRec* discardedHit) const
 704 {
 705   /// Compute the weight matrices, in non bending and bending direction,
 706   /// of the other attached hits assuming the discarded one does not exist
 707   /// accounting for multiple scattering correlations and hits resolution
 708   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
 709   /// - Return kFALSE if computation failed
 710
 711   // Check MCS covariance matrix and recompute it if need
 712   Bool_t deleteMCSCov = kFALSE;
 713   if (!mcsCovariances) {
 714     mcsCovariances = new TMatrixD(fNTrackHits,fNTrackHits);
 715     deleteMCSCov = kTRUE;
 716     ComputeMCSCovariances(*mcsCovariances);
 717   }
 718
 719   // Resize the weights matrices; alocate memory
 720   if (discardedHit) {
 721     hitWeightsNB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
 722     hitWeightsB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
 723   } else {
 724     hitWeightsNB.ResizeTo(fNTrackHits,fNTrackHits);
 725     hitWeightsB.ResizeTo(fNTrackHits,fNTrackHits);
 726   }
 727
 728   // Define variables
 729   AliMUONHitForRec *hitForRec1, *hitForRec2;
 730   Int_t currentHitNumber1, currentHitNumber2;
 731
 732   // Compute the covariance matrices
 733   currentHitNumber1 = 0;
 734   for (Int_t hitNumber1 = 0; hitNumber1 < fNTrackHits; hitNumber1++) {
 735     hitForRec1 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1))->GetHitForRecPtr();
 736
 737     if (hitForRec1 == discardedHit) continue;
 738
 739     // Loop over next hits
 740     currentHitNumber2 = currentHitNumber1;
 741     for (Int_t hitNumber2 = hitNumber1; hitNumber2 < fNTrackHits; hitNumber2++) {
 742       hitForRec2 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber2))->GetHitForRecPtr();
 743
 744       if (hitForRec2 == discardedHit) continue;
 745
 746       // Fill with MCS covariances
 747       hitWeightsNB(currentHitNumber1, currentHitNumber2) = (*mcsCovariances)(hitNumber1,hitNumber2);
 748
 749       // Equal contribution from multiple scattering in non bending and bending directions
 750       hitWeightsB(currentHitNumber1, currentHitNumber2) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
 751
 752       // Add contribution from hit resolution to diagonal element and symmetrize the matrix
 753       if (currentHitNumber1 == currentHitNumber2) {
 754
 755         // In non bending plane
 756         hitWeightsNB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetNonBendingReso2();
 757         // In bending plane
 758         hitWeightsB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetBendingReso2();
 759
 760       } else {
 761
 762         // In non bending plane
 763         hitWeightsNB(currentHitNumber2, currentHitNumber1) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
 764         // In bending plane
 765         hitWeightsB(currentHitNumber2, currentHitNumber1) = hitWeightsB(currentHitNumber1, currentHitNumber2);
 766
 767       }
 768
 769       currentHitNumber2++;
 770     }
 771
 772     currentHitNumber1++;
 773   }
 774
 775   // Inversion of covariance matrices to get the weights
 776   if (hitWeightsNB.Determinant() != 0 && hitWeightsB.Determinant() != 0) {
 777     hitWeightsNB.Invert();
 778     hitWeightsB.Invert();
 779   } else {
 780     AliWarning(" Determinant = 0");
 781     hitWeightsNB.ResizeTo(0,0);
 782     hitWeightsB.ResizeTo(0,0);
 783     if(deleteMCSCov) delete mcsCovariances;
 784     return kFALSE;
 785   }
 786
 787   if(deleteMCSCov) delete mcsCovariances;
 788
 789   return kTRUE;
 790
 791 }
 792
 793   //__________________________________________________________________________
 794 void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
 795 {
 796   /// Compute the multiple scattering covariance matrix
 797   /// (assume that track parameters at each hit are corrects)
 798
 799   // Reset the size of the covariance matrix if needed
 800   if (mcsCovariances.GetNrows() != fNTrackHits) mcsCovariances.ResizeTo(fNTrackHits,fNTrackHits);
 801
 802   // Define variables
 803   Int_t nChambers = AliMUONConstants::NTrackingCh();
 804   AliMUONTrackParam* trackParamAtHit;
 805   AliMUONHitForRec *hitForRec;
 806   AliMUONTrackParam extrapTrackParam;
 807   Int_t currentChamber = 0, expectedChamber = 0, size = 0;
 808   Double_t *mcsAngle2 = new Double_t[2*nChambers];
 809   Double_t *zMCS = new Double_t[2*nChambers];
 810   Int_t *indices = new Int_t[2*fNTrackHits];
 811
 812   // Compute multiple scattering dispersion angle at each chamber
 813   // and save the z position where it is calculated
 814   for (Int_t hitNumber = 0; hitNumber < fNTrackHits; hitNumber++) {
 815     trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
 816     hitForRec = trackParamAtHit->GetHitForRecPtr();
 817
 818     // look for missing chambers if any
 819     currentChamber = hitForRec->GetChamberNumber();
 820     while (currentChamber > expectedChamber) {
 821
 822       // Save the z position where MCS dispersion is calculated
 823       zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
 824
 825       // Do not take into account MCS in chambers prior the first hit
 826       if (hitNumber > 0) {
 827
 828         // Get track parameters at missing chamber z
 829         extrapTrackParam = *trackParamAtHit;
 830         AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[expectedChamber]);
 831
 832         // Save multiple scattering dispersion angle in missing chamber
 833         mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
 834
 835       } else mcsAngle2[size] = 0.;
 836
 837       expectedChamber++;
 838       size++;
 839     }
 840
 841     // Save z position where MCS dispersion is calculated
 842     zMCS[size] = trackParamAtHit->GetZ();
 843
 844     // Save multiple scattering dispersion angle in current chamber
 845     mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.);
 846
 847     // Save indice in zMCS array corresponding to the current cluster
 848     indices[hitNumber] = size;
 849
 850     expectedChamber = currentChamber + 1;
 851     size++;
 852   }
 853
 854   // complete array of z if last hit is on the last but one chamber
 855   if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
 856
 857   // Compute the covariance matrix
 858   for (Int_t hitNumber1 = 0; hitNumber1 < fNTrackHits; hitNumber1++) {
 859
 860     for (Int_t hitNumber2 = hitNumber1; hitNumber2 < fNTrackHits; hitNumber2++) {
 861
 862       // Initialization to 0 (diagonal plus upper triangular part)
 863       mcsCovariances(hitNumber1,hitNumber2) = 0.;
 864
 865       // Compute contribution from multiple scattering in upstream chambers
 866       for (Int_t k = 0; k < indices[hitNumber1]; k++) {
 867         mcsCovariances(hitNumber1,hitNumber2) += (zMCS[indices[hitNumber1]] - zMCS[k]) * (zMCS[indices[hitNumber2]] - zMCS[k]) * mcsAngle2[k];
 868       }
 869
 870       // Symetrize the matrix
 871       mcsCovariances(hitNumber2,hitNumber1) = mcsCovariances(hitNumber1,hitNumber2);
 872     }
 873
 874   }
 875
 876   delete [] mcsAngle2;
 877   delete [] zMCS;
 878   delete [] indices;
 879
 880 }
 881
 882   //__________________________________________________________________________
 883 Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* track) const
 884 {
 885   /// Returns the number of hits in common between the current track ("this")
 886   /// and the track pointed to by "track".
 887   Int_t hitsInCommon = 0;
 888   AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2;
 889   // Loop over hits of first track
 890   trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->First();
 891   while (trackParamAtHit1) {
 892     // Loop over hits of second track
 893     trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->First();
 894     while (trackParamAtHit2) {
 895       // Increment "hitsInCommon" if both TrackParamAtHits point to the same HitForRec
 896       if ((trackParamAtHit1->GetHitForRecPtr()) == (trackParamAtHit2->GetHitForRecPtr())) {
 897         hitsInCommon++;
 898         break;
 899       }
 900       trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->After(trackParamAtHit2);
 901     } // trackParamAtHit2
 902     trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->After(trackParamAtHit1);
 903   } // trackParamAtHit1
 904   return hitsInCommon;
 905 }
 906
 907   //__________________________________________________________________________
 908 Double_t AliMUONTrack::GetNormalizedChi2() const
 909 {
 910   /// return the chi2 value divided by the number of degrees of freedom (or 1.e10 if ndf < 0)
 911
 912   Double_t numberOfDegFree = (2. * fNTrackHits - 5.);
 913   if (numberOfDegFree > 0.) return fGlobalChi2 / numberOfDegFree;
 914   else return 1.e10;
 915 }
 916
 917   //__________________________________________________________________________
 918 Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * track, Double_t sigma2Cut) const
 919 {
 920   /// Return kTRUE/kFALSE for each chamber if hit is compatible or not
 921   TClonesArray *hitArray, *thisHitArray;
 922   AliMUONHitForRec *hit, *thisHit;
 923   Int_t chamberNumber;
 924   Float_t deltaZ;
 925   Float_t deltaZMax = 1.; // 1 cm
 926   Float_t chi2 = 0;
 927   Bool_t *nCompHit = new Bool_t[AliMUONConstants::NTrackingCh()];
 928
 929   for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 930     nCompHit[ch] = kFALSE;
 931   }
 932
 933   thisHitArray = this->GetHitForRecAtHit();
 934
 935   hitArray =  track->GetHitForRecAtHit();
 936
 937   for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) {
 938     thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis);
 939     chamberNumber = thisHit->GetChamberNumber();
 940     if (chamberNumber < 0 || chamberNumber > AliMUONConstants::NTrackingCh()) continue;
 941     nCompHit[chamberNumber] = kFALSE;
 942     for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) {
 943       hit = (AliMUONHitForRec*) hitArray->At(iH);
 944       deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ());
 945       chi2 = thisHit->NormalizedChi2WithHitForRec(hit,sigma2Cut); // set cut to 4 sigmas
 946       if (chi2 < 3. && deltaZ < deltaZMax) {
 947         nCompHit[chamberNumber] = kTRUE;
 948         break;
 949       }
 950     }
 951   }
 952
 953   return nCompHit;
 954 }
 955
 956   //__________________________________________________________________________
 957 void AliMUONTrack::RecursiveDump(void) const
 958 {
 959   /// Recursive dump of AliMUONTrack, i.e. with dump of TrackParamAtHit's and attached HitForRec's
 960   AliMUONTrackParam *trackParamAtHit;
 961   AliMUONHitForRec *hitForRec;
 962   cout << "Recursive dump of Track: " << this << endl;
 963   // Track
 964   this->Dump();
 965   for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
 966     trackParamAtHit = (AliMUONTrackParam*) ((*fTrackParamAtHit)[trackHitIndex]);
 967     // TrackHit
 968     cout << "TrackParamAtHit: " << trackParamAtHit << " (index: " << trackHitIndex << ")" << endl;
 969     trackParamAtHit->Dump();
 970     hitForRec = trackParamAtHit->GetHitForRecPtr();
 971     // HitForRec
 972     cout << "HitForRec: " << hitForRec << endl;
 973     hitForRec->Dump();
 974   }
 975   return;
 976 }
 977
 978 //_____________________________________________-
 979 void AliMUONTrack::Print(Option_t*) const
 980 {
 981   /// Printing Track information
 982
 983   cout << "<AliMUONTrack> No.Clusters=" << setw(2)   << GetNTrackHits() <<
 984       ", Match2Trig=" << setw(1) << GetMatchTrigger()  <<
 985       ", LoTrgNum=" << setw(3) << GetLoTrgNum()  <<
 986     ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) <<  GetChi2MatchTrigger();
 987   cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh) << endl;
 988   GetTrackParamAtHit()->First()->Print("FULL");
 989 }
 990
 991 //__________________________________________________________________________
 992 void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt)
 993 {
 994   /// pack the local trigger information and store
 995
 996   if (loCirc < 0 || loCirc > 233) return;
 997
 998   fLocalTrigger = 0;
 999   fLocalTrigger += loCirc;
1000   fLocalTrigger += loStripX << 8;
1001   fLocalTrigger += loStripY << 13;
1002   fLocalTrigger += loDev    << 17;
1003   fLocalTrigger += loLpt    << 22;
1004   fLocalTrigger += loHpt    << 24;
1005
1006 }
1007