[u/mrichter/AliRoot.git] / MUON / AliMUONSegment.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////
//
// Segment for reconstruction
// in 
// ALICE 
// dimuon 
// spectrometer:
// two hits for reconstruction in the two chambers of one station
//
///////////////////////////////////////////////////////////

#include "AliMUONSegment.h" 
#include "AliMUON.h"
#include "AliMUONChamber.h" 
#include "AliMUONHitForRec.h" 
#include "AliMUONTrackParam.h" 
#include "AliRun.h" // for gAlice

ClassImp(AliMUONSegment) // Class implementation in ROOT context

  //__________________________________________________________________________
AliMUONSegment::AliMUONSegment()
  : TObject()
{
  // Default constructor
  fHitForRecPtr1 = 0; // pointer to HitForRec in first chamber
  fHitForRecPtr2 = 0; // pointer to HitForRec in second chamber
  // Bending plane:
  fBendingCoor = 0.0; // Coordinate in bending plane
  fBendingSlope = 0.0; // Slope in bending plane
  // Covariance in bending plane:
  fBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
  fBendingSlopeReso2 = 0.0; // Covariance(slope)
  fBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
  fBendingImpact = 0.0; // Impact parameter in bending plane
  // Non Bending plane:
  fNonBendingCoor = 0.0; // Coordinate in non bending plane
  fNonBendingSlope = 0.0; // Slope in non bending plane
  // Covariance in non bending plane:
  fNonBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
  fNonBendingSlopeReso2 = 0.0; // Covariance(slope)
  fNonBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
  fNonBendingImpact = 0.0; // Impact parameter in non bending plane
  fInTrack = kFALSE; // TRUE if segment belongs to one track
}

  //__________________________________________________________________________
AliMUONSegment::AliMUONSegment(AliMUONHitForRec* Hit1, AliMUONHitForRec* Hit2)
  : TObject()
{
  // Constructor for AliMUONSegment from two HitForRec's,
  // one, in the first chamber of the station, pointed to by "Hit1",
  // the other one, in the second chamber of the station, pointed to by "Hit1".
  // Fills the pointers to both hits,
  // the slope, the covariance for (coordinate in first chamber, slope),
  // and the impact parameter at vertex (Z=0),
  // in bending and non bending planes.
  // Puts the "fInTrack" flag to "kFALSE".
  Double_t dz;
  // pointers to HitForRec's
  fHitForRecPtr1 = Hit1;
  fHitForRecPtr2 = Hit2;
  dz = Hit1->GetZ() - Hit2->GetZ();
  // bending plane
  fBendingCoor = Hit1->GetBendingCoor();
  fBendingSlope = (fBendingCoor - Hit2->GetBendingCoor()) / dz;
  fBendingImpact = fBendingCoor - Hit1->GetZ() * fBendingSlope;
  fBendingCoorReso2 = Hit1->GetBendingReso2();
  fBendingSlopeReso2 = ( Hit1->GetBendingReso2() +
			 Hit2->GetBendingReso2() ) / dz / dz;
  fBendingCoorSlopeReso2 = Hit1->GetBendingReso2() / dz;
  // non bending plane
  fNonBendingCoor = Hit1->GetNonBendingCoor();
  fNonBendingSlope = (fNonBendingCoor - Hit2->GetNonBendingCoor()) / dz;
  fNonBendingImpact = fNonBendingCoor - Hit1->GetZ() * fNonBendingSlope;
  fNonBendingCoorReso2 = Hit1->GetNonBendingReso2();
  fNonBendingSlopeReso2 = ( Hit1->GetNonBendingReso2() +
			    Hit2->GetNonBendingReso2() ) / dz / dz;
  fNonBendingCoorSlopeReso2 = Hit1->GetNonBendingReso2() / dz;
  // "fInTrack" flag to "kFALSE"
  fInTrack = kFALSE;
  return;
}

AliMUONSegment::AliMUONSegment (const AliMUONSegment& theMUONSegment)
  : TObject(theMUONSegment)
{
// Protected copy constructor

  Fatal("AliMUONSegment", "Not implemented.");
}

AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& rhs)
{
// Protected assignement operator

  if (this == &rhs) return *this;

  Fatal("operator=", "Not implemented.");
    
  return *this;  
}

  //__________________________________________________________________________
Int_t AliMUONSegment::Compare(const TObject* Segment) const
{
  // "Compare" function to sort with increasing absolute value
  // of the "impact parameter" in bending plane.
  // Returns -1 (0, +1) if |impact parameter| of current Segment
  // is smaller than (equal to, larger than) |impact parameter| of Segment
  if (TMath::Abs(((AliMUONSegment*)this)->fBendingImpact)
      < TMath::Abs(((AliMUONSegment*)Segment)->fBendingImpact))
    return(-1);
  // continuous parameter, hence no need for testing equal case
  else return(+1);
}

  //__________________________________________________________________________
Double_t AliMUONSegment::NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut) const
{
  // Calculate the normalized Chi2 between the current Segment (this)
  // and the Segment pointed to by "Segment",
  // i.e. the square deviations between the coordinates and the slopes,
  // in both the bending and the non bending plane,
  // divided by the variance of the same quantities and by "Sigma2Cut".
  // Returns 5 if none of the 4 quantities is OK,
  // something smaller than or equal to 4 otherwise.
  // Would it be more correct to use a real chi square
  // including the non diagonal term ????
  Double_t chi2, chi2Max, diff, normDiff;
  chi2 = 0.0;
  chi2Max = 5.0;
  // coordinate in bending plane
  diff = this->fBendingCoor - Segment->fBendingCoor;
  normDiff = diff * diff /
    (this->fBendingCoorReso2 + Segment->fBendingCoorReso2) / Sigma2Cut;
  if (normDiff > 1.0) return chi2Max;
  chi2 = chi2 + normDiff;
  // slope in bending plane
  diff = this->fBendingSlope - Segment->fBendingSlope;
  normDiff = diff * diff /
    (this->fBendingSlopeReso2 + Segment->fBendingSlopeReso2) / Sigma2Cut;
  if (normDiff > 1.0) return chi2Max;
  chi2 = chi2 + normDiff;
  // coordinate in non bending plane
  diff = this->fNonBendingCoor - Segment->fNonBendingCoor;
  normDiff = diff * diff /
    (this->fNonBendingCoorReso2 + Segment->fNonBendingCoorReso2) / Sigma2Cut;
  if (normDiff > 1.0) return chi2Max;
  chi2 = chi2 + normDiff;
  // slope in non bending plane
  diff = this->fNonBendingSlope - Segment->fNonBendingSlope;
  normDiff = diff * diff /
    (this->fNonBendingSlopeReso2 + Segment->fNonBendingSlopeReso2) / Sigma2Cut;
  if (normDiff > 1.0) return chi2Max;
  chi2 = chi2 + normDiff;
  return chi2;
}

  //__________________________________________________________________________
AliMUONSegment* AliMUONSegment::CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor) const
{
  // Extrapolates linearly the current Segment (this) to station (0..) "Station".
  // Multiple Coulomb scattering calculated from "MCSfactor"
  // corresponding to one chamber,
  // with one chamber for the coordinate, two chambers for the angle,
  // due to the arrangement in stations.
  // Valid from station(1..) 4 to 5 or vice versa.
  // Returns the pointer to the created AliMUONSegment object
  // corresponding to this extrapolation.
  // The caller has the responsibility to delete this object.
  AliMUONSegment* extrapSegment = new AliMUONSegment(); // creates empty new segment
  // dZ from first hit of current Segment to first chamber of station "Station"
  AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
  Double_t dZ =
    (&(pMUON->Chamber(2 * Station)))->Z() - (this->fHitForRecPtr1)->GetZ();
  // Data in bending plane
  //  coordinate
  extrapSegment->fBendingCoor = this->fBendingCoor + this->fBendingSlope * dZ;
  //  slope
  extrapSegment->fBendingSlope = this->fBendingSlope;
  //  covariance, including multiple Coulomb scattering over dZ due to one chamber
  extrapSegment->fBendingCoorReso2 = this->fBendingCoorReso2 +
    (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ; // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
  extrapSegment->fBendingSlopeReso2 = this->fBendingSlopeReso2 + 2.0 * MCSfactor;
  extrapSegment->fBendingCoorSlopeReso2 =
    this->fBendingCoorSlopeReso2 + this->fBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
  // Data in non bending plane
  //  coordinate
  extrapSegment->fNonBendingCoor =
    this->fNonBendingCoor + this->fNonBendingSlope * dZ;
  //  slope
  extrapSegment->fNonBendingSlope = this->fNonBendingSlope;
  //  covariance, including multiple Coulomb scattering over dZ due to one chamber
  extrapSegment->fNonBendingCoorReso2 = this->fNonBendingCoorReso2 +
    (this->fNonBendingSlopeReso2 + MCSfactor) *dZ * dZ; // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
  extrapSegment->fNonBendingSlopeReso2 =
    this->fNonBendingSlopeReso2 + 2.0 * MCSfactor;
  extrapSegment->fNonBendingCoorSlopeReso2 =
    this->fNonBendingCoorSlopeReso2 + this->fNonBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
  return extrapSegment;
}

  //__________________________________________________________________________
AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor) const
{
  // Extrapolates linearly the current Segment (this) to chamber(0..) "Chamber".
  // Multiple Coulomb scattering calculated from "MCSfactor"
  // corresponding to one chamber.
  // Valid from station(1..) 4 to 5 or vice versa.
  // Returns the pointer to the created AliMUONHitForRec object
  // corresponding to this extrapolation.
  // The caller has the responsibility to delete this object.
  AliMUONHitForRec* extrapHitForRec = new AliMUONHitForRec(); // creates empty new HitForRec
  // dZ from first hit of current Segment to chamber
  AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
  Double_t dZ =
    (&(pMUON->Chamber(Chamber)))->Z() - (this->fHitForRecPtr1)->GetZ();
  // Data in bending plane
  //  coordinate
  extrapHitForRec->SetBendingCoor(this->fBendingCoor + this->fBendingSlope * dZ);
  //  covariance, including multiple Coulomb scattering over dZ due to one chamber
  extrapHitForRec->SetBendingReso2(this->fBendingCoorReso2 +
				   (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
  // Data in non bending plane
  //  coordinate
 extrapHitForRec ->SetNonBendingCoor(this->fNonBendingCoor +
				     this->fNonBendingSlope * dZ);
  //  covariance, including multiple Coulomb scattering over dZ due to one chamber
  extrapHitForRec->
    SetNonBendingReso2(this->fNonBendingCoorReso2 +
		       (this->fNonBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
  return extrapHitForRec;
}

  //__________________________________________________________________________
void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t /*MCSfactor*/, Double_t /*Dz1*/, Double_t /*Dz2*/, Double_t /*Dz3*/, Int_t Station, Double_t InverseMomentum)
{
  // Fill data members with values calculated from the array of track parameters
  // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
  // of the station, respectively).
  // Multiple Coulomb scattering is taking into account with "MCSfactor"
  // corresponding to one chamber,
  // with one chamber for the coordinate, two chambers for the angle,
  // due to the arrangement in stations.
  // Resolution coming from:
  // coordinate in closest station at "Dz1" from current "Station",
  // slope between closest stations, with "Dz2" interval between them,
  // interval "Dz3" between chambers of closest station,
  // extrapolation over "Dz1" from closest station,
  // "InverseMomentum".
  // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
  // are assumed to be filled
  // with the variance on bending and non bending coordinates.
  // The "road" is parametrized from the old reco_muon.F
  // with 8 cm between stations.
  AliMUONTrackParam *param0;
//   Double_t cReso2, sReso2;
  // parameters to define the widths of the searching roads in station 0,1,2
  // width = p0 + p1/ (momentum)^2
  //                  station number:        0         1          2
//   static Double_t p0BendingCoor[3] =     { 6.43e-2, 1.64e-2,   0.034 };   
//   static Double_t p1BendingCoor[3] =     {    986.,    821.,    446. };  
//   static Double_t p0BendingSlope[3] =    { 3.54e-6, 3.63e-6,  3.6e-6 };  
//   static Double_t p1BendingSlope[3] =    { 4.49e-3,  4.8e-3,   0.011 };  
//   static Double_t p0NonBendingCoor[3] =  { 4.66e-2, 4.83e-2,   0.049 };   
//   static Double_t p1NonBendingCoor[3] =  {   1444.,    866.,    354. };  
//   static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };  
//   static Double_t p1NonBendingSlope[3] = {      0.,      0.,      0. };
  
  static Double_t p0BendingCoor[3] =     { 6.43e-2, 6.43e-2,   6.43e-2  };   
  static Double_t p1BendingCoor[3] =     {    986.,    986.,       986. };  
  static Double_t p0BendingSlope[3] =    {   3.6e-6,   3.6e-6,     3.6e-6  };  
  static Double_t p1BendingSlope[3] =    {  1.1e-2,  1.1e-2,    1.1e-2  };  
  static Double_t p0NonBendingCoor[3] =  {   0.049,   0.049,     0.049  };   
  static Double_t p1NonBendingCoor[3] =  {   1444.,   1444.,      1444. };  
  static Double_t p0NonBendingSlope[3] = {   6.8e-4,   6.8e-4,     6.8e-4  };  
  static Double_t p1NonBendingSlope[3] = {      0.,      0.,         0. };  
  param0 = &(TrackParam[0]);

// OLD version
//   // Bending plane
//   fBendingCoor = param0->GetBendingCoor(); // coordinate
//   fBendingSlope = param0->GetBendingSlope(); // slope
//   cReso2 = fBendingCoorReso2;
//   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
//   fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
//   fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
//   // Non bending plane
//   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
//   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
//   cReso2 = fNonBendingCoorReso2;
//   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
//   fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
//   fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;

  // Coordinate and slope
  // Bending plane
  fBendingCoor = param0->GetBendingCoor(); // coordinate
  fBendingSlope = param0->GetBendingSlope(); // slope
  // Non bending plane
  fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
  fNonBendingSlope = param0->GetNonBendingSlope(); // slope

  // Resolutions
  // cReso2 and sReso2 have to be subtracted here from the parametrization
  // because they are added in the functions "NormalizedChi2WithSegment"
  // and "NormalizedChi2WithHitForRec"
  // Bending plane
//   cReso2 = fBendingCoorReso2;
//   sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
  fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum ;  // - cReso2
  fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum; //  - sReso2;
  // Non bending plane
//   cReso2 = fNonBendingCoorReso2;
//   sReso2 =  (2. * cReso2 )/ (Dz3*Dz3) ;
  fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum; // - cReso2;
  fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum; //  - sReso2;
  return;
}

// OLD function, with roads automatically calculated instead from being parametrized
// kept because it would be a better solution,
// if one can really find the right values.
//   //__________________________________________________________________________
// void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
// {
//   // Fill data members with values calculated from the array of track parameters
//   // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
//   // of the station, respectively).
//   // Multiple Coulomb scattering is taking into account with "MCSfactor"
//   // corresponding to one chamber,
//   // with one chamber for the coordinate, two chambers for the angle,
//   // due to the arrangement in stations.
//   // Resolution coming from:
//   // coordinate in closest station at "Dz1",
//   // slope between closest stations, with "Dz2" interval between them,
//   // extrapolation over "Dz" from closest station.
//   // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
//   // are assumed to be filled
//   // with the variance on bending and non bending coordinates.
//   AliMUONTrackParam *param0;
//   Double_t cReso2, sReso2;
//   param0 = &(TrackParam[0]);
//   // Bending plane
//   fBendingCoor = param0->GetBendingCoor(); // coordinate
//   fBendingSlope = param0->GetBendingSlope(); // slope
//   cReso2 = fBendingCoorReso2;
//   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
//   fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
//   fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
//   // Non bending plane
//   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
//   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
//   cReso2 = fNonBendingCoorReso2;
//   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
//   fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
//   fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
//   return;
// }
Commit	Line	Data
a9e2aefa	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
88cb7938	16	/* $Id$ */
a9e2aefa	17
3831f268	18	///////////////////////////////////////////////////////////
a9e2aefa	19	//
3831f268	20	// Segment for reconstruction
	21	// in
	22	// ALICE
	23	// dimuon
	24	// spectrometer:
a9e2aefa	25	// two hits for reconstruction in the two chambers of one station
3831f268	26	//
3831f268	27	///////////////////////////////////////////////////////////
a9e2aefa	28
30178c30	29	#include "AliMUONSegment.h"
a9e2aefa	30	#include "AliMUON.h"
3831f268	31	#include "AliMUONChamber.h"
a9e2aefa	32	#include "AliMUONHitForRec.h"
a9e2aefa	33	#include "AliMUONTrackParam.h"
3831f268	34	#include "AliRun.h" // for gAlice
a9e2aefa	35
	36	ClassImp(AliMUONSegment) // Class implementation in ROOT context
	37
762647ff	38	//__________________________________________________________________________
762647ff	39	AliMUONSegment::AliMUONSegment()
30178c30	40	: TObject()
762647ff	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
a9e2aefa	64	//__________________________________________________________________________
a9e2aefa	65	AliMUONSegment::AliMUONSegment(AliMUONHitForRec* Hit1, AliMUONHitForRec* Hit2)
30178c30	66	: TObject()
a9e2aefa	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
30178c30	102	AliMUONSegment::AliMUONSegment (const AliMUONSegment& theMUONSegment)
30178c30	103	: TObject(theMUONSegment)
a9e2aefa	104	{
30178c30	105	// Protected copy constructor
	106
	107	Fatal("AliMUONSegment", "Not implemented.");
a9e2aefa	108	}
a9e2aefa	109
30178c30	110	AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& rhs)
a9e2aefa	111	{
30178c30	112	// Protected assignement operator
	113
	114	if (this == &rhs) return *this;
	115
	116	Fatal("operator=", "Not implemented.");
	117
	118	return *this;
a9e2aefa	119	}
a9e2aefa	120
a9e2aefa	121	//__________________________________________________________________________
2a941f4e	122	Int_t AliMUONSegment::Compare(const TObject* Segment) const
a9e2aefa	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
31817925	128	if (TMath::Abs(((AliMUONSegment*)this)->fBendingImpact)
31817925	129	< TMath::Abs(((AliMUONSegment*)Segment)->fBendingImpact))
a9e2aefa	130	return(-1);
	131	// continuous parameter, hence no need for testing equal case
	132	else return(+1);
	133	}
	134
	135	//__________________________________________________________________________
30178c30	136	Double_t AliMUONSegment::NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut) const
a9e2aefa	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	//__________________________________________________________________________
30178c30	178	AliMUONSegment* AliMUONSegment::CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor) const
a9e2aefa	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	//__________________________________________________________________________
30178c30	222	AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor) const
a9e2aefa	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	//__________________________________________________________________________
d9a3473d	254	void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam TrackParam, Double_t /MCSfactor/, Double_t /Dz1/, Double_t /Dz2/, Double_t /Dz3*/, Int_t Station, Double_t InverseMomentum)
a9e2aefa	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:
04b5ea16	264	// coordinate in closest station at "Dz1" from current "Station",
a9e2aefa	265	// slope between closest stations, with "Dz2" interval between them,
04b5ea16	266	// interval "Dz3" between chambers of closest station,
	267	// extrapolation over "Dz1" from closest station,
	268	// "InverseMomentum".
a9e2aefa	269	// When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
	270	// are assumed to be filled
	271	// with the variance on bending and non bending coordinates.
04b5ea16	272	// The "road" is parametrized from the old reco_muon.F
04b5ea16	273	// with 8 cm between stations.
a9e2aefa	274	AliMUONTrackParam *param0;
d0bfce8d	275	// Double_t cReso2, sReso2;
04b5ea16	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
d0bfce8d	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. };
a9e2aefa	296	param0 = &(TrackParam[0]);
04b5ea16	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
a9e2aefa	315	// Bending plane
	316	fBendingCoor = param0->GetBendingCoor(); // coordinate
	317	fBendingSlope = param0->GetBendingSlope(); // slope
a9e2aefa	318	// Non bending plane
	319	fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
	320	fNonBendingSlope = param0->GetNonBendingSlope(); // slope
04b5ea16	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
d0bfce8d	327	// cReso2 = fBendingCoorReso2;
	328	// sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
	329	fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]InverseMomentumInverseMomentum ; // - cReso2
	330	fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]InverseMomentumInverseMomentum; // - sReso2;
04b5ea16	331	// Non bending plane
d0bfce8d	332	// cReso2 = fNonBendingCoorReso2;
	333	// sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
	334	fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]InverseMomentumInverseMomentum; // - cReso2;
	335	fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]InverseMomentumInverseMomentum; // - sReso2;
a9e2aefa	336	return;
a9e2aefa	337	}
04b5ea16	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	// }