MUON/AliMUONSegment.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 //
  20 // Segment for reconstruction
  21 // in
  22 // ALICE
  23 // dimuon
  24 // spectrometer:
  25 // two hits for reconstruction in the two chambers of one station
  26 //
  27 ///////////////////////////////////////////////////////////
  28
  29 #include "AliMUONSegment.h"
  30 #include "AliMUON.h"
  31 #include "AliMUONChamber.h"
  32 #include "AliMUONHitForRec.h"
  33 #include "AliMUONTrackParam.h"
  34 #include "AliRun.h" // for gAlice
  35
  36 ClassImp(AliMUONSegment) // Class implementation in ROOT context
  37
  38   //__________________________________________________________________________
  39 AliMUONSegment::AliMUONSegment()
  40   : TObject()
  41 {
  42   // Default constructor
  43   fHitForRecPtr1 = 0; // pointer to HitForRec in first chamber
  44   fHitForRecPtr2 = 0; // pointer to HitForRec in second chamber
  45   // Bending plane:
  46   fBendingCoor = 0.0; // Coordinate in bending plane
  47   fBendingSlope = 0.0; // Slope in bending plane
  48   // Covariance in bending plane:
  49   fBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
  50   fBendingSlopeReso2 = 0.0; // Covariance(slope)
  51   fBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
  52   fBendingImpact = 0.0; // Impact parameter in bending plane
  53   // Non Bending plane:
  54   fNonBendingCoor = 0.0; // Coordinate in non bending plane
  55   fNonBendingSlope = 0.0; // Slope in non bending plane
  56   // Covariance in non bending plane:
  57   fNonBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
  58   fNonBendingSlopeReso2 = 0.0; // Covariance(slope)
  59   fNonBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
  60   fNonBendingImpact = 0.0; // Impact parameter in non bending plane
  61   fInTrack = kFALSE; // TRUE if segment belongs to one track
  62 }
  63
  64   //__________________________________________________________________________
  65 AliMUONSegment::AliMUONSegment(AliMUONHitForRec* Hit1, AliMUONHitForRec* Hit2)
  66   : TObject()
  67 {
  68   // Constructor for AliMUONSegment from two HitForRec's,
  69   // one, in the first chamber of the station, pointed to by "Hit1",
  70   // the other one, in the second chamber of the station, pointed to by "Hit1".
  71   // Fills the pointers to both hits,
  72   // the slope, the covariance for (coordinate in first chamber, slope),
  73   // and the impact parameter at vertex (Z=0),
  74   // in bending and non bending planes.
  75   // Puts the "fInTrack" flag to "kFALSE".
  76   Double_t dz;
  77   // pointers to HitForRec's
  78   fHitForRecPtr1 = Hit1;
  79   fHitForRecPtr2 = Hit2;
  80   dz = Hit1->GetZ() - Hit2->GetZ();
  81   // bending plane
  82   fBendingCoor = Hit1->GetBendingCoor();
  83   fBendingSlope = (fBendingCoor - Hit2->GetBendingCoor()) / dz;
  84   fBendingImpact = fBendingCoor - Hit1->GetZ() * fBendingSlope;
  85   fBendingCoorReso2 = Hit1->GetBendingReso2();
  86   fBendingSlopeReso2 = ( Hit1->GetBendingReso2() +
  87                          Hit2->GetBendingReso2() ) / dz / dz;
  88   fBendingCoorSlopeReso2 = Hit1->GetBendingReso2() / dz;
  89   // non bending plane
  90   fNonBendingCoor = Hit1->GetNonBendingCoor();
  91   fNonBendingSlope = (fNonBendingCoor - Hit2->GetNonBendingCoor()) / dz;
  92   fNonBendingImpact = fNonBendingCoor - Hit1->GetZ() * fNonBendingSlope;
  93   fNonBendingCoorReso2 = Hit1->GetNonBendingReso2();
  94   fNonBendingSlopeReso2 = ( Hit1->GetNonBendingReso2() +
  95                             Hit2->GetNonBendingReso2() ) / dz / dz;
  96   fNonBendingCoorSlopeReso2 = Hit1->GetNonBendingReso2() / dz;
  97   // "fInTrack" flag to "kFALSE"
  98   fInTrack = kFALSE;
  99   return;
 100 }
 101
 102 AliMUONSegment::AliMUONSegment (const AliMUONSegment& theMUONSegment)
 103   : TObject(theMUONSegment)
 104 {
 105 // Protected copy constructor
 106
 107   Fatal("AliMUONSegment", "Not implemented.");
 108 }
 109
 110 AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& rhs)
 111 {
 112 // Protected assignement operator
 113
 114   if (this == &rhs) return *this;
 115
 116   Fatal("operator=", "Not implemented.");
 117
 118   return *this;
 119 }
 120
 121   //__________________________________________________________________________
 122 Int_t AliMUONSegment::Compare(const TObject* Segment) const
 123 {
 124   // "Compare" function to sort with increasing absolute value
 125   // of the "impact parameter" in bending plane.
 126   // Returns -1 (0, +1) if |impact parameter| of current Segment
 127   // is smaller than (equal to, larger than) |impact parameter| of Segment
 128   if (TMath::Abs(((AliMUONSegment*)this)->fBendingImpact)
 129       < TMath::Abs(((AliMUONSegment*)Segment)->fBendingImpact))
 130     return(-1);
 131   // continuous parameter, hence no need for testing equal case
 132   else return(+1);
 133 }
 134
 135   //__________________________________________________________________________
 136 Double_t AliMUONSegment::NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut) const
 137 {
 138   // Calculate the normalized Chi2 between the current Segment (this)
 139   // and the Segment pointed to by "Segment",
 140   // i.e. the square deviations between the coordinates and the slopes,
 141   // in both the bending and the non bending plane,
 142   // divided by the variance of the same quantities and by "Sigma2Cut".
 143   // Returns 5 if none of the 4 quantities is OK,
 144   // something smaller than or equal to 4 otherwise.
 145   // Would it be more correct to use a real chi square
 146   // including the non diagonal term ????
 147   Double_t chi2, chi2Max, diff, normDiff;
 148   chi2 = 0.0;
 149   chi2Max = 5.0;
 150   // coordinate in bending plane
 151   diff = this->fBendingCoor - Segment->fBendingCoor;
 152   normDiff = diff * diff /
 153     (this->fBendingCoorReso2 + Segment->fBendingCoorReso2) / Sigma2Cut;
 154   if (normDiff > 1.0) return chi2Max;
 155   chi2 = chi2 + normDiff;
 156   // slope in bending plane
 157   diff = this->fBendingSlope - Segment->fBendingSlope;
 158   normDiff = diff * diff /
 159     (this->fBendingSlopeReso2 + Segment->fBendingSlopeReso2) / Sigma2Cut;
 160   if (normDiff > 1.0) return chi2Max;
 161   chi2 = chi2 + normDiff;
 162   // coordinate in non bending plane
 163   diff = this->fNonBendingCoor - Segment->fNonBendingCoor;
 164   normDiff = diff * diff /
 165     (this->fNonBendingCoorReso2 + Segment->fNonBendingCoorReso2) / Sigma2Cut;
 166   if (normDiff > 1.0) return chi2Max;
 167   chi2 = chi2 + normDiff;
 168   // slope in non bending plane
 169   diff = this->fNonBendingSlope - Segment->fNonBendingSlope;
 170   normDiff = diff * diff /
 171     (this->fNonBendingSlopeReso2 + Segment->fNonBendingSlopeReso2) / Sigma2Cut;
 172   if (normDiff > 1.0) return chi2Max;
 173   chi2 = chi2 + normDiff;
 174   return chi2;
 175 }
 176
 177   //__________________________________________________________________________
 178 AliMUONSegment* AliMUONSegment::CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor) const
 179 {
 180   // Extrapolates linearly the current Segment (this) to station (0..) "Station".
 181   // Multiple Coulomb scattering calculated from "MCSfactor"
 182   // corresponding to one chamber,
 183   // with one chamber for the coordinate, two chambers for the angle,
 184   // due to the arrangement in stations.
 185   // Valid from station(1..) 4 to 5 or vice versa.
 186   // Returns the pointer to the created AliMUONSegment object
 187   // corresponding to this extrapolation.
 188   // The caller has the responsibility to delete this object.
 189   AliMUONSegment* extrapSegment = new AliMUONSegment(); // creates empty new segment
 190   // dZ from first hit of current Segment to first chamber of station "Station"
 191   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 192   Double_t dZ =
 193     (&(pMUON->Chamber(2 * Station)))->Z() - (this->fHitForRecPtr1)->GetZ();
 194   // Data in bending plane
 195   //  coordinate
 196   extrapSegment->fBendingCoor = this->fBendingCoor + this->fBendingSlope * dZ;
 197   //  slope
 198   extrapSegment->fBendingSlope = this->fBendingSlope;
 199   //  covariance, including multiple Coulomb scattering over dZ due to one chamber
 200   extrapSegment->fBendingCoorReso2 = this->fBendingCoorReso2 +
 201     (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ; // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
 202   extrapSegment->fBendingSlopeReso2 = this->fBendingSlopeReso2 + 2.0 * MCSfactor;
 203   extrapSegment->fBendingCoorSlopeReso2 =
 204     this->fBendingCoorSlopeReso2 + this->fBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
 205   // Data in non bending plane
 206   //  coordinate
 207   extrapSegment->fNonBendingCoor =
 208     this->fNonBendingCoor + this->fNonBendingSlope * dZ;
 209   //  slope
 210   extrapSegment->fNonBendingSlope = this->fNonBendingSlope;
 211   //  covariance, including multiple Coulomb scattering over dZ due to one chamber
 212   extrapSegment->fNonBendingCoorReso2 = this->fNonBendingCoorReso2 +
 213     (this->fNonBendingSlopeReso2 + MCSfactor) *dZ * dZ; // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
 214   extrapSegment->fNonBendingSlopeReso2 =
 215     this->fNonBendingSlopeReso2 + 2.0 * MCSfactor;
 216   extrapSegment->fNonBendingCoorSlopeReso2 =
 217     this->fNonBendingCoorSlopeReso2 + this->fNonBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
 218   return extrapSegment;
 219 }
 220
 221   //__________________________________________________________________________
 222 AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor) const
 223 {
 224   // Extrapolates linearly the current Segment (this) to chamber(0..) "Chamber".
 225   // Multiple Coulomb scattering calculated from "MCSfactor"
 226   // corresponding to one chamber.
 227   // Valid from station(1..) 4 to 5 or vice versa.
 228   // Returns the pointer to the created AliMUONHitForRec object
 229   // corresponding to this extrapolation.
 230   // The caller has the responsibility to delete this object.
 231   AliMUONHitForRec* extrapHitForRec = new AliMUONHitForRec(); // creates empty new HitForRec
 232   // dZ from first hit of current Segment to chamber
 233   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 234   Double_t dZ =
 235     (&(pMUON->Chamber(Chamber)))->Z() - (this->fHitForRecPtr1)->GetZ();
 236   // Data in bending plane
 237   //  coordinate
 238   extrapHitForRec->SetBendingCoor(this->fBendingCoor + this->fBendingSlope * dZ);
 239   //  covariance, including multiple Coulomb scattering over dZ due to one chamber
 240   extrapHitForRec->SetBendingReso2(this->fBendingCoorReso2 +
 241                                    (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
 242   // Data in non bending plane
 243   //  coordinate
 244  extrapHitForRec ->SetNonBendingCoor(this->fNonBendingCoor +
 245                                      this->fNonBendingSlope * dZ);
 246   //  covariance, including multiple Coulomb scattering over dZ due to one chamber
 247   extrapHitForRec->
 248     SetNonBendingReso2(this->fNonBendingCoorReso2 +
 249                        (this->fNonBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
 250   return extrapHitForRec;
 251 }
 252
 253   //__________________________________________________________________________
 254 void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t /*MCSfactor*/, Double_t /*Dz1*/, Double_t /*Dz2*/, Double_t /*Dz3*/, Int_t Station, Double_t InverseMomentum)
 255 {
 256   // Fill data members with values calculated from the array of track parameters
 257   // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
 258   // of the station, respectively).
 259   // Multiple Coulomb scattering is taking into account with "MCSfactor"
 260   // corresponding to one chamber,
 261   // with one chamber for the coordinate, two chambers for the angle,
 262   // due to the arrangement in stations.
 263   // Resolution coming from:
 264   // coordinate in closest station at "Dz1" from current "Station",
 265   // slope between closest stations, with "Dz2" interval between them,
 266   // interval "Dz3" between chambers of closest station,
 267   // extrapolation over "Dz1" from closest station,
 268   // "InverseMomentum".
 269   // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
 270   // are assumed to be filled
 271   // with the variance on bending and non bending coordinates.
 272   // The "road" is parametrized from the old reco_muon.F
 273   // with 8 cm between stations.
 274   AliMUONTrackParam *param0;
 275 //   Double_t cReso2, sReso2;
 276   // parameters to define the widths of the searching roads in station 0,1,2
 277   // width = p0 + p1/ (momentum)^2
 278   //                  station number:        0         1          2
 279 //   static Double_t p0BendingCoor[3] =     { 6.43e-2, 1.64e-2,   0.034 };
 280 //   static Double_t p1BendingCoor[3] =     {    986.,    821.,    446. };
 281 //   static Double_t p0BendingSlope[3] =    { 3.54e-6, 3.63e-6,  3.6e-6 };
 282 //   static Double_t p1BendingSlope[3] =    { 4.49e-3,  4.8e-3,   0.011 };
 283 //   static Double_t p0NonBendingCoor[3] =  { 4.66e-2, 4.83e-2,   0.049 };
 284 //   static Double_t p1NonBendingCoor[3] =  {   1444.,    866.,    354. };
 285 //   static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };
 286 //   static Double_t p1NonBendingSlope[3] = {      0.,      0.,      0. };
 287
 288   static Double_t p0BendingCoor[3] =     { 6.43e-2, 6.43e-2,   6.43e-2  };
 289   static Double_t p1BendingCoor[3] =     {    986.,    986.,       986. };
 290   static Double_t p0BendingSlope[3] =    {   3.6e-6,   3.6e-6,     3.6e-6  };
 291   static Double_t p1BendingSlope[3] =    {  1.1e-2,  1.1e-2,    1.1e-2  };
 292   static Double_t p0NonBendingCoor[3] =  {   0.049,   0.049,     0.049  };
 293   static Double_t p1NonBendingCoor[3] =  {   1444.,   1444.,      1444. };
 294   static Double_t p0NonBendingSlope[3] = {   6.8e-4,   6.8e-4,     6.8e-4  };
 295   static Double_t p1NonBendingSlope[3] = {      0.,      0.,         0. };
 296   param0 = &(TrackParam[0]);
 297
 298 // OLD version
 299 //   // Bending plane
 300 //   fBendingCoor = param0->GetBendingCoor(); // coordinate
 301 //   fBendingSlope = param0->GetBendingSlope(); // slope
 302 //   cReso2 = fBendingCoorReso2;
 303 //   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
 304 //   fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
 305 //   fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
 306 //   // Non bending plane
 307 //   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
 308 //   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
 309 //   cReso2 = fNonBendingCoorReso2;
 310 //   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
 311 //   fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
 312 //   fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
 313
 314   // Coordinate and slope
 315   // Bending plane
 316   fBendingCoor = param0->GetBendingCoor(); // coordinate
 317   fBendingSlope = param0->GetBendingSlope(); // slope
 318   // Non bending plane
 319   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
 320   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
 321
 322   // Resolutions
 323   // cReso2 and sReso2 have to be subtracted here from the parametrization
 324   // because they are added in the functions "NormalizedChi2WithSegment"
 325   // and "NormalizedChi2WithHitForRec"
 326   // Bending plane
 327 //   cReso2 = fBendingCoorReso2;
 328 //   sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
 329   fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum ;  // - cReso2
 330   fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum; //  - sReso2;
 331   // Non bending plane
 332 //   cReso2 = fNonBendingCoorReso2;
 333 //   sReso2 =  (2. * cReso2 )/ (Dz3*Dz3) ;
 334   fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum; // - cReso2;
 335   fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum; //  - sReso2;
 336   return;
 337 }
 338
 339 // OLD function, with roads automatically calculated instead from being parametrized
 340 // kept because it would be a better solution,
 341 // if one can really find the right values.
 342 //   //__________________________________________________________________________
 343 // void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
 344 // {
 345 //   // Fill data members with values calculated from the array of track parameters
 346 //   // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
 347 //   // of the station, respectively).
 348 //   // Multiple Coulomb scattering is taking into account with "MCSfactor"
 349 //   // corresponding to one chamber,
 350 //   // with one chamber for the coordinate, two chambers for the angle,
 351 //   // due to the arrangement in stations.
 352 //   // Resolution coming from:
 353 //   // coordinate in closest station at "Dz1",
 354 //   // slope between closest stations, with "Dz2" interval between them,
 355 //   // extrapolation over "Dz" from closest station.
 356 //   // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
 357 //   // are assumed to be filled
 358 //   // with the variance on bending and non bending coordinates.
 359 //   AliMUONTrackParam *param0;
 360 //   Double_t cReso2, sReso2;
 361 //   param0 = &(TrackParam[0]);
 362 //   // Bending plane
 363 //   fBendingCoor = param0->GetBendingCoor(); // coordinate
 364 //   fBendingSlope = param0->GetBendingSlope(); // slope
 365 //   cReso2 = fBendingCoorReso2;
 366 //   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
 367 //   fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
 368 //   fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
 369 //   // Non bending plane
 370 //   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
 371 //   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
 372 //   cReso2 = fNonBendingCoorReso2;
 373 //   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
 374 //   fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
 375 //   fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
 376 //   return;
 377 // }