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