[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.3  2000/06/25 13:23:28  hristov
stdlib.h needed for non-Linux compilation

Revision 1.2  2000/06/15 07:58:48  morsch
Code from MUON-dev joined

Revision 1.1.2.3  2000/06/12 10:11:34  morsch
Dummy copy constructor and assignment operator added

Revision 1.1.2.2  2000/06/09 12:58:05  gosset
Removed comment beginnings in Log sections of .cxx files
Suppressed most violations of coding rules

Revision 1.1.2.1  2000/06/07 14:44:53  gosset
Addition of files for track reconstruction in C++
*/

//__________________________________________________________________________
//
// Reconstructed track in ALICE dimuon spectrometer
//__________________________________________________________________________

#include "AliMUONTrack.h"

#include <iostream.h>

#include <TClonesArray.h>
#include <TMinuit.h>
#include <TMath.h>
#include <TMatrix.h>

#include "AliMUONEventReconstructor.h" 
#include "AliMUONHitForRec.h" 
#include "AliMUONSegment.h" 
#include "AliMUONTrackHit.h"

#include <stdlib.h>

// variables to be known from minimization functions
static AliMUONTrack *trackBeingFitted;
static AliMUONTrackParam *trackParamBeingFitted;

// Functions to be minimized with Minuit
void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);

Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);

ClassImp(AliMUONTrack) // Class implementation in ROOT context

  //__________________________________________________________________________
AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor)
{
  // Constructor from two Segment's
  fEventReconstructor = EventReconstructor; // link back to EventReconstructor
  // memory allocation for the TClonesArray of reconstructed TrackHit's
  fTrackHitsPtr = new  TClonesArray("AliMUONTrackHit", 10);
  fNTrackHits = 0;
  AddSegment(BegSegment); // add hits from BegSegment
  AddSegment(EndSegment); // add hits from EndSegment
  fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
  SetTrackParamAtVertex(); // set track parameters at vertex
  fFitMCS = 0;
  return;
}

  //__________________________________________________________________________
AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor)
{
  // Constructor from one Segment and one HitForRec
  fEventReconstructor = EventReconstructor; // link back to EventReconstructor
  // memory allocation for the TClonesArray of reconstructed TrackHit's
  fTrackHitsPtr = new  TClonesArray("AliMUONTrackHit", 10);
  fNTrackHits = 0;
  AddSegment(Segment); // add hits from Segment
  AddHitForRec(HitForRec); // add HitForRec
  fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
  SetTrackParamAtVertex(); // set track parameters at vertex
  fFitMCS = 0;
  return;
}

AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack)
{
// Dummy copy constructor
}

AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
{
// Dummy assignment operator
    return *this;
}

  //__________________________________________________________________________
void AliMUONTrack::SetFitMCS(Int_t FitMCS)
{
  // Set track fit option with or without multiple Coulomb scattering
  // from "FitMCS" argument: 0 without, 1 with
  fFitMCS = FitMCS;
  // better implementation with enum(with, without) ????
  return;
}

// Inline functions for Get and Set: inline removed because it does not work !!!!
AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex(void) {
  // Get pointer to fTrackParamAtVertex
  return &fTrackParamAtVertex;}
AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) {
  // Get pointer to TrackParamAtFirstHit
  return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
TClonesArray* AliMUONTrack::GetTrackHitsPtr(void) {
  // Get fTrackHitsPtr
  return fTrackHitsPtr;}
Int_t AliMUONTrack::GetNTrackHits(void) {
  // Get fNTrackHits
  return fNTrackHits;}

  //__________________________________________________________________________
void AliMUONTrack::RecursiveDump(void)
{
  // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's
  AliMUONTrackHit *trackHit;
  AliMUONHitForRec *hitForRec;
  cout << "Recursive dump of Track: " << this << endl;
  // Track
  this->Dump();
  for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
    trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[trackHitIndex]);
    // TrackHit
    cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl;
    trackHit->Dump();
    hitForRec = trackHit->GetHitForRecPtr();
    // HitForRec
    cout << "HitForRec: " << hitForRec << endl;
    hitForRec->Dump();
  }
  return;
}

  //__________________________________________________________________________
void AliMUONTrack::Fit(AliMUONTrackParam *TrackParam, Int_t NParam)
{
  // Fit the current track ("this"),
  // starting with track parameters pointed to by "TrackParam",
  // and with 3 or 5 parameters ("NParam"):
  // 3 if one keeps X and Y fixed in "TrackParam",
  // 5 if one lets them vary.
  if ((NParam != 3) && (NParam != 5)) {
    cout << "ERROR in AliMUONTrack::Fit, NParam = " << NParam;
    cout << " , i.e. neither 3 nor 5 ====> EXIT" << endl;
    exit(0); // right instruction for exit ????
  }
  Int_t error = 0;
  Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
  TString parName;
  TMinuit *minuit = new TMinuit(5);
  trackBeingFitted = this; // for the track to be known from the function to minimize
  trackParamBeingFitted = TrackParam; // for the track parameters to be known from the function to minimize; possible to use only Minuit parameters ????
  // try to use TVirtualFitter to get this feature automatically !!!!
  minuit->mninit(5, 10, 7); // sysrd, syswr, syssa: useful ????
  // how to go faster ???? choice of Minuit parameters like EDM ????
  // choice of function to be minimized according to fFitMCS
  if (fFitMCS == 0) minuit->SetFCN(TrackChi2);
  else minuit->SetFCN(TrackChi2MCS);
  minuit->SetPrintLevel(1); // More printing !!!!
  // set first 3 parameters
  // could be tried with no limits for the search (min=max=0) ????
  minuit->mnparm(0, "InvBenP",
		 TrackParam->GetInverseBendingMomentum(),
		 0.003, -0.4, 0.4, error);
  minuit->mnparm(1, "BenS",
		 TrackParam->GetBendingSlope(),
		 0.001, -0.5, 0.5, error);
  minuit->mnparm(2, "NonBenS",
		 TrackParam->GetNonBendingSlope(),
		 0.001, -0.5, 0.5, error);
  if (NParam == 5) {
    // set last 2 parameters (no limits for the search: min=max=0)
    minuit->mnparm(3, "X",
		   TrackParam->GetNonBendingCoor(),
		   0.03, 0.0, 0.0, error);
    minuit->mnparm(4, "Y",
		   TrackParam->GetBendingCoor(),
		   0.10, 0.0, 0.0, error);
  }
  // search without gradient calculation in the function
  minuit->mnexcm("SET NOGRADIENT", arg, 0, error);
  // minimization
  minuit->mnexcm("MINIMIZE", arg, 0, error);
  // exit from Minuit
  minuit->mnexcm("EXIT", arg, 0, error); // necessary ????
  // print results
  minuit->mnpout(0, parName, invBenP, errorParam, lower, upper, error);
  minuit->mnpout(1, parName, benC, errorParam, lower, upper, error);
  minuit->mnpout(2, parName, nonBenC, errorParam, lower, upper, error);
  if (NParam == 5) {
    minuit->mnpout(3, parName, x, errorParam, lower, upper, error);
    minuit->mnpout(4, parName, y, errorParam, lower, upper, error);
  }
  // result of the fit into track parameters
  TrackParam->SetInverseBendingMomentum(invBenP);
  TrackParam->SetBendingSlope(benC);
  TrackParam->SetNonBendingSlope(nonBenC);
  if (NParam == 5) {
    TrackParam->SetNonBendingCoor(x);
    TrackParam->SetBendingCoor(y);
  }
  trackBeingFitted = NULL;
  delete minuit;
}

  //__________________________________________________________________________
void AliMUONTrack::AddSegment(AliMUONSegment* Segment)
{
  // Add Segment
  AddHitForRec(Segment->GetHitForRec1()); // 1st hit
  AddHitForRec(Segment->GetHitForRec2()); // 2nd hit
}

  //__________________________________________________________________________
void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec)
{
  // Add HitForRec
  new ((*fTrackHitsPtr)[fNTrackHits]) AliMUONTrackHit(HitForRec);
  fNTrackHits++;
}

  //__________________________________________________________________________
void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam)
{
  // Set track parameters at TrackHit with index "indexHit"
  // from the track parameters pointed to by "TrackParam".
  AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]);
  trackHit->SetTrackParam(TrackParam);
}

  //__________________________________________________________________________
void AliMUONTrack::SetTrackParamAtVertex()
{
  // Set track parameters at vertex.
  // TrackHit's are assumed to be only in stations(1..) 4 and 5,
  // and sorted according to increasing Z..
  // Parameters are calculated from information in HitForRec's
  // of first and last TrackHit's.
  AliMUONTrackParam *trackParam =
    &fTrackParamAtVertex; // pointer to track parameters
  // Pointer to HitForRec of first TrackHit
  AliMUONHitForRec *firstHit =
    ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr();
  // Pointer to HitForRec of last TrackHit
  AliMUONHitForRec *lastHit =
    ((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr();
  // Z difference between first and last hits
  Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
  // bending slope in stations(1..) 4 and 5
  Double_t bendingSlope =
    (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
  trackParam->SetBendingSlope(bendingSlope);
  // impact parameter
  Double_t impactParam =
    firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and  lastHit ????
  // signed bending momentum
  Double_t signedBendingMomentum =
    fEventReconstructor->GetBendingMomentumFromImpactParam(impactParam);
  trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum);
  // bending slope at vertex
  trackParam->
    SetBendingSlope(bendingSlope +
		    impactParam / fEventReconstructor->GetSimpleBPosition());
  // non bending slope
  Double_t nonBendingSlope =
    (firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ;
  trackParam->SetNonBendingSlope(nonBendingSlope);
  // vertex coordinates at (0,0,0)
  trackParam->SetZ(0.0);
  trackParam->SetBendingCoor(0.0);
  trackParam->SetNonBendingCoor(0.0);
}

  //__________________________________________________________________________
void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
{
  // Return the "Chi2" to be minimized with Minuit for track fitting,
  // with "NParam" parameters
  // and their current values in array pointed to by "Param".
  // Assumes that the track hits are sorted according to increasing Z.
  // Track parameters at each TrackHit are updated accordingly.
  // Multiple Coulomb scattering is not taken into account
  AliMUONTrackHit* hit;
  AliMUONTrackParam param1;
  Int_t hitNumber;
  Double_t zHit;
  Chi2 = 0.0; // initialize Chi2
  // copy of track parameters to be fitted
  param1 = *trackParamBeingFitted;
  // Minuit parameters to be fitted into this copy
  param1.SetInverseBendingMomentum(Param[0]);
  param1.SetBendingSlope(Param[1]);
  param1.SetNonBendingSlope(Param[2]);
  if (NParam == 5) {
    param1.SetNonBendingCoor(Param[3]);
    param1.SetBendingCoor(Param[4]);
  }
  // Follow track through all planes of track hits
  for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) {
    hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
    zHit = hit->GetHitForRecPtr()->GetZ();
    // do something special if 2 hits with same Z ????
    // security against infinite loop ????
    (&param1)->ExtrapToZ(zHit); // extrapolation
    hit->SetTrackParam(&param1);
    // Increment Chi2
    // done hit per hit, with hit resolution,
    // and not with point and angle like in "reco_muon.F" !!!!
    // Needs to add multiple scattering contribution ????
    Double_t dX =
      hit->GetHitForRecPtr()->GetNonBendingCoor() - (&param1)->GetNonBendingCoor();
    Double_t dY =
      hit->GetHitForRecPtr()->GetBendingCoor() - (&param1)->GetBendingCoor();
    Chi2 =
      Chi2 +
      dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() +
      dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
  }
}

  //__________________________________________________________________________
void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
{
  // Return the "Chi2" to be minimized with Minuit for track fitting,
  // with "NParam" parameters
  // and their current values in array pointed to by "Param".
  // Assumes that the track hits are sorted according to increasing Z.
  // Track parameters at each TrackHit are updated accordingly.
  // Multiple Coulomb scattering is taken into account with covariance matrix.
  AliMUONTrackParam param1;
  Chi2 = 0.0; // initialize Chi2
  // copy of track parameters to be fitted
  param1 = *trackParamBeingFitted;
  // Minuit parameters to be fitted into this copy
  param1.SetInverseBendingMomentum(Param[0]);
  param1.SetBendingSlope(Param[1]);
  param1.SetNonBendingSlope(Param[2]);
  if (NParam == 5) {
    param1.SetNonBendingCoor(Param[3]);
    param1.SetBendingCoor(Param[4]);
  }

  AliMUONTrackHit* hit, hit1, hit2, hit3;
  Bool_t GoodDeterminant;
  Int_t hitNumber, hitNumber1, hitNumber2, hitNumber3;
  Double_t zHit[10], paramBendingCoor[10], paramNonBendingCoor[10], ap[10];
  Double_t hitBendingCoor[10], hitNonBendingCoor[10];
  Double_t hitBendingReso2[10], hitNonBendingReso2[10];
  Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10);
  TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
  TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);

  // Predicted coordinates and  multiple scattering angles are first calculated
  for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
    hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
    zHit[hitNumber] = hit->GetHitForRecPtr()->GetZ();
    // do something special if 2 hits with same Z ????
    // security against infinite loop ????
    (&param1)->ExtrapToZ(zHit[hitNumber]); // extrapolation
    hit->SetTrackParam(&param1);
    paramBendingCoor[hitNumber]= (&param1)->GetBendingCoor();
    paramNonBendingCoor[hitNumber]= (&param1)->GetNonBendingCoor();
    hitBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetBendingCoor();
    hitNonBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingCoor();
    hitBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetBendingReso2();
    hitNonBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingReso2();
    ap[hitNumber] = MultipleScatteringAngle2(hit); // multiple scatt. angle ^2  
  }

  // Calculates the covariance matrix
  for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {    
    for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
      (*covBending)(hitNumber1, hitNumber2) = 0;
      (*covBending)(hitNumber2, hitNumber1) = 0;
      if (hitNumber1 == hitNumber2){ // diagonal elements
	(*covBending)(hitNumber2, hitNumber1) =
	  (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
      }
      // Multiple Scattering...  loop on upstream chambers ??
      for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++){ 	
	(*covBending)(hitNumber2, hitNumber1) =
	  (*covBending)(hitNumber2, hitNumber1) +
	  ((zHit[hitNumber1] - zHit[hitNumber3]) *
	   (zHit[hitNumber2] - zHit[hitNumber3]) * ap[hitNumber3]); 
      }  
      (*covNonBending)(hitNumber1, hitNumber2) = 0;
      (*covNonBending)(hitNumber2, hitNumber1) =
	(*covBending)(hitNumber2, hitNumber1);
      if (hitNumber1 == hitNumber2) {  // diagonal elements
	(*covNonBending)(hitNumber2, hitNumber1) =
	  (*covNonBending)(hitNumber2, hitNumber1) -
	  hitBendingReso2[hitNumber1] + hitNonBendingReso2[hitNumber1] ;
      }      
    }
  }

  // Inverts covariance matrix 
  GoodDeterminant = kTRUE;
  if (covBending->Determinant() != 0) {
    covBending->Invert();
  } else {
    GoodDeterminant = kFALSE;
    cout << "Warning in ChiMCS  Determinant Bending=0: " << endl;  
  }
  if (covNonBending->Determinant() != 0){
    covNonBending->Invert();
  } else {
    GoodDeterminant = kFALSE;
    cout << "Warning in ChiMCS  Determinant non Bending=0: " << endl;  
  }
  
  // Calculates Chi2
  if (GoodDeterminant) { // with Multiple Scattering if inversion correct
    for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++){ 
      for (hitNumber2=0; hitNumber2 < numberOfHit; hitNumber2++){
	Chi2 = Chi2 +
	  ((*covBending)(hitNumber2, hitNumber1) * 
	   (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
	   (hitBendingCoor[hitNumber2] - paramBendingCoor[hitNumber2]));
	Chi2 = Chi2 +
	  ((*covNonBending)(hitNumber2, hitNumber1) *
	   (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
	   (hitNonBendingCoor[hitNumber2] - paramNonBendingCoor[hitNumber2]));
      }
    }
  } else {  // without Multiple Scattering if inversion impossible
    for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++) { 
      Chi2 = Chi2 +
	((hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
	 (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) /
	 hitBendingReso2[hitNumber1]);
      Chi2 = Chi2 +
	((hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
	 (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) /
	 hitNonBendingReso2[hitNumber1]);      
    }
  }
  
  delete covBending;
  delete covNonBending;
}

Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
{
  // Returns square of multiple Coulomb scattering angle
  // at TrackHit pointed to by "TrackHit"
  Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
  Double_t varMultipleScatteringAngle;
  AliMUONTrackParam *param = TrackHit->GetTrackParam();
  slopeBending = param->GetBendingSlope();
  slopeNonBending = param->GetNonBendingSlope();
  // thickness in radiation length for the current track,
  // taking local angle into account
  radiationLength =
    trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() *
    TMath::Sqrt(1.0 +
		slopeBending * slopeBending + slopeNonBending * slopeNonBending);
  inverseBendingMomentum2 = 
    param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
  inverseTotalMomentum2 =
    inverseBendingMomentum2 * (1.0 + slopeBending*slopeBending) /
    (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); 
  varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
  varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
    varMultipleScatteringAngle * varMultipleScatteringAngle;
  return varMultipleScatteringAngle;
}
Commit	Line	Data
	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	/*
	17	$Log$
	18	Revision 1.3 2000/06/25 13:23:28 hristov
	19	stdlib.h needed for non-Linux compilation
	20
	21	Revision 1.2 2000/06/15 07:58:48 morsch
	22	Code from MUON-dev joined
	23
	24	Revision 1.1.2.3 2000/06/12 10:11:34 morsch
	25	Dummy copy constructor and assignment operator added
	26
	27	Revision 1.1.2.2 2000/06/09 12:58:05 gosset
	28	Removed comment beginnings in Log sections of .cxx files
	29	Suppressed most violations of coding rules
	30
	31	Revision 1.1.2.1 2000/06/07 14:44:53 gosset
	32	Addition of files for track reconstruction in C++
	33	*/
	34
	35	//__________________________________________________________________________
	36	//
	37	// Reconstructed track in ALICE dimuon spectrometer
	38	//__________________________________________________________________________
	39
	40	#include "AliMUONTrack.h"
	41
	42	#include <iostream.h>
	43
	44	#include <TClonesArray.h>
	45	#include <TMinuit.h>
	46	#include <TMath.h>
	47	#include <TMatrix.h>
	48
	49	#include "AliMUONEventReconstructor.h"
	50	#include "AliMUONHitForRec.h"
	51	#include "AliMUONSegment.h"
	52	#include "AliMUONTrackHit.h"
	53
	54	#include <stdlib.h>
	55
	56	// variables to be known from minimization functions
	57	static AliMUONTrack *trackBeingFitted;
	58	static AliMUONTrackParam *trackParamBeingFitted;
	59
	60	// Functions to be minimized with Minuit
	61	void TrackChi2(Int_t &NParam, Double_t Gradient, Double_t &Chi2, Double_t Param, Int_t Flag);
	62	void TrackChi2MCS(Int_t &NParam, Double_t Gradient, Double_t &Chi2, Double_t Param, Int_t Flag);
	63
	64	Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
	65
	66	ClassImp(AliMUONTrack) // Class implementation in ROOT context
	67
	68	//__________________________________________________________________________
	69	AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor)
	70	{
	71	// Constructor from two Segment's
	72	fEventReconstructor = EventReconstructor; // link back to EventReconstructor
	73	// memory allocation for the TClonesArray of reconstructed TrackHit's
	74	fTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 10);
	75	fNTrackHits = 0;
	76	AddSegment(BegSegment); // add hits from BegSegment
	77	AddSegment(EndSegment); // add hits from EndSegment
	78	fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
	79	SetTrackParamAtVertex(); // set track parameters at vertex
	80	fFitMCS = 0;
	81	return;
	82	}
	83
	84	//__________________________________________________________________________
	85	AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor)
	86	{
	87	// Constructor from one Segment and one HitForRec
	88	fEventReconstructor = EventReconstructor; // link back to EventReconstructor
	89	// memory allocation for the TClonesArray of reconstructed TrackHit's
	90	fTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 10);
	91	fNTrackHits = 0;
	92	AddSegment(Segment); // add hits from Segment
	93	AddHitForRec(HitForRec); // add HitForRec
	94	fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
	95	SetTrackParamAtVertex(); // set track parameters at vertex
	96	fFitMCS = 0;
	97	return;
	98	}
	99
	100	AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack)
	101	{
	102	// Dummy copy constructor
	103	}
	104
	105	AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
	106	{
	107	// Dummy assignment operator
	108	return *this;
	109	}
	110
	111	//__________________________________________________________________________
	112	void AliMUONTrack::SetFitMCS(Int_t FitMCS)
	113	{
	114	// Set track fit option with or without multiple Coulomb scattering
	115	// from "FitMCS" argument: 0 without, 1 with
	116	fFitMCS = FitMCS;
	117	// better implementation with enum(with, without) ????
	118	return;
	119	}
	120
	121	// Inline functions for Get and Set: inline removed because it does not work !!!!
	122	AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex(void) {
	123	// Get pointer to fTrackParamAtVertex
	124	return &fTrackParamAtVertex;}
	125	AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) {
	126	// Get pointer to TrackParamAtFirstHit
	127	return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
	128	TClonesArray* AliMUONTrack::GetTrackHitsPtr(void) {
	129	// Get fTrackHitsPtr
	130	return fTrackHitsPtr;}
	131	Int_t AliMUONTrack::GetNTrackHits(void) {
	132	// Get fNTrackHits
	133	return fNTrackHits;}
	134
	135	//__________________________________________________________________________
	136	void AliMUONTrack::RecursiveDump(void)
	137	{
	138	// Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's
	139	AliMUONTrackHit *trackHit;
	140	AliMUONHitForRec *hitForRec;
	141	cout << "Recursive dump of Track: " << this << endl;
	142	// Track
	143	this->Dump();
	144	for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
	145	trackHit = (AliMUONTrackHit) ((fTrackHitsPtr)[trackHitIndex]);
	146	// TrackHit
	147	cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl;
	148	trackHit->Dump();
	149	hitForRec = trackHit->GetHitForRecPtr();
	150	// HitForRec
	151	cout << "HitForRec: " << hitForRec << endl;
	152	hitForRec->Dump();
	153	}
	154	return;
	155	}
	156
	157	//__________________________________________________________________________
	158	void AliMUONTrack::Fit(AliMUONTrackParam *TrackParam, Int_t NParam)
	159	{
	160	// Fit the current track ("this"),
	161	// starting with track parameters pointed to by "TrackParam",
	162	// and with 3 or 5 parameters ("NParam"):
	163	// 3 if one keeps X and Y fixed in "TrackParam",
	164	// 5 if one lets them vary.
	165	if ((NParam != 3) && (NParam != 5)) {
	166	cout << "ERROR in AliMUONTrack::Fit, NParam = " << NParam;
	167	cout << " , i.e. neither 3 nor 5 ====> EXIT" << endl;
	168	exit(0); // right instruction for exit ????
	169	}
	170	Int_t error = 0;
	171	Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
	172	TString parName;
	173	TMinuit *minuit = new TMinuit(5);
	174	trackBeingFitted = this; // for the track to be known from the function to minimize
	175	trackParamBeingFitted = TrackParam; // for the track parameters to be known from the function to minimize; possible to use only Minuit parameters ????
	176	// try to use TVirtualFitter to get this feature automatically !!!!
	177	minuit->mninit(5, 10, 7); // sysrd, syswr, syssa: useful ????
	178	// how to go faster ???? choice of Minuit parameters like EDM ????
	179	// choice of function to be minimized according to fFitMCS
	180	if (fFitMCS == 0) minuit->SetFCN(TrackChi2);
	181	else minuit->SetFCN(TrackChi2MCS);
	182	minuit->SetPrintLevel(1); // More printing !!!!
	183	// set first 3 parameters
	184	// could be tried with no limits for the search (min=max=0) ????
	185	minuit->mnparm(0, "InvBenP",
	186	TrackParam->GetInverseBendingMomentum(),
	187	0.003, -0.4, 0.4, error);
	188	minuit->mnparm(1, "BenS",
	189	TrackParam->GetBendingSlope(),
	190	0.001, -0.5, 0.5, error);
	191	minuit->mnparm(2, "NonBenS",
	192	TrackParam->GetNonBendingSlope(),
	193	0.001, -0.5, 0.5, error);
	194	if (NParam == 5) {
	195	// set last 2 parameters (no limits for the search: min=max=0)
	196	minuit->mnparm(3, "X",
	197	TrackParam->GetNonBendingCoor(),
	198	0.03, 0.0, 0.0, error);
	199	minuit->mnparm(4, "Y",
	200	TrackParam->GetBendingCoor(),
	201	0.10, 0.0, 0.0, error);
	202	}
	203	// search without gradient calculation in the function
	204	minuit->mnexcm("SET NOGRADIENT", arg, 0, error);
	205	// minimization
	206	minuit->mnexcm("MINIMIZE", arg, 0, error);
	207	// exit from Minuit
	208	minuit->mnexcm("EXIT", arg, 0, error); // necessary ????
	209	// print results
	210	minuit->mnpout(0, parName, invBenP, errorParam, lower, upper, error);
	211	minuit->mnpout(1, parName, benC, errorParam, lower, upper, error);
	212	minuit->mnpout(2, parName, nonBenC, errorParam, lower, upper, error);
	213	if (NParam == 5) {
	214	minuit->mnpout(3, parName, x, errorParam, lower, upper, error);
	215	minuit->mnpout(4, parName, y, errorParam, lower, upper, error);
	216	}
	217	// result of the fit into track parameters
	218	TrackParam->SetInverseBendingMomentum(invBenP);
	219	TrackParam->SetBendingSlope(benC);
	220	TrackParam->SetNonBendingSlope(nonBenC);
	221	if (NParam == 5) {
	222	TrackParam->SetNonBendingCoor(x);
	223	TrackParam->SetBendingCoor(y);
	224	}
	225	trackBeingFitted = NULL;
	226	delete minuit;
	227	}
	228
	229	//__________________________________________________________________________
	230	void AliMUONTrack::AddSegment(AliMUONSegment* Segment)
	231	{
	232	// Add Segment
	233	AddHitForRec(Segment->GetHitForRec1()); // 1st hit
	234	AddHitForRec(Segment->GetHitForRec2()); // 2nd hit
	235	}
	236
	237	//__________________________________________________________________________
	238	void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec)
	239	{
	240	// Add HitForRec
	241	new ((*fTrackHitsPtr)[fNTrackHits]) AliMUONTrackHit(HitForRec);
	242	fNTrackHits++;
	243	}
	244
	245	//__________________________________________________________________________
	246	void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam)
	247	{
	248	// Set track parameters at TrackHit with index "indexHit"
	249	// from the track parameters pointed to by "TrackParam".
	250	AliMUONTrackHit* trackHit = (AliMUONTrackHit) ((fTrackHitsPtr)[indexHit]);
	251	trackHit->SetTrackParam(TrackParam);
	252	}
	253
	254	//__________________________________________________________________________
	255	void AliMUONTrack::SetTrackParamAtVertex()
	256	{
	257	// Set track parameters at vertex.
	258	// TrackHit's are assumed to be only in stations(1..) 4 and 5,
	259	// and sorted according to increasing Z..
	260	// Parameters are calculated from information in HitForRec's
	261	// of first and last TrackHit's.
	262	AliMUONTrackParam *trackParam =
	263	&fTrackParamAtVertex; // pointer to track parameters
	264	// Pointer to HitForRec of first TrackHit
	265	AliMUONHitForRec *firstHit =
	266	((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr();
	267	// Pointer to HitForRec of last TrackHit
	268	AliMUONHitForRec *lastHit =
	269	((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr();
	270	// Z difference between first and last hits
	271	Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
	272	// bending slope in stations(1..) 4 and 5
	273	Double_t bendingSlope =
	274	(firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
	275	trackParam->SetBendingSlope(bendingSlope);
	276	// impact parameter
	277	Double_t impactParam =
	278	firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and lastHit ????
	279	// signed bending momentum
	280	Double_t signedBendingMomentum =
	281	fEventReconstructor->GetBendingMomentumFromImpactParam(impactParam);
	282	trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum);
	283	// bending slope at vertex
	284	trackParam->
	285	SetBendingSlope(bendingSlope +
	286	impactParam / fEventReconstructor->GetSimpleBPosition());
	287	// non bending slope
	288	Double_t nonBendingSlope =
	289	(firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ;
	290	trackParam->SetNonBendingSlope(nonBendingSlope);
	291	// vertex coordinates at (0,0,0)
	292	trackParam->SetZ(0.0);
	293	trackParam->SetBendingCoor(0.0);
	294	trackParam->SetNonBendingCoor(0.0);
	295	}
	296
	297	//__________________________________________________________________________
	298	void TrackChi2(Int_t &NParam, Double_t Gradient, Double_t &Chi2, Double_t Param, Int_t Flag)
	299	{
	300	// Return the "Chi2" to be minimized with Minuit for track fitting,
	301	// with "NParam" parameters
	302	// and their current values in array pointed to by "Param".
	303	// Assumes that the track hits are sorted according to increasing Z.
	304	// Track parameters at each TrackHit are updated accordingly.
	305	// Multiple Coulomb scattering is not taken into account
	306	AliMUONTrackHit* hit;
	307	AliMUONTrackParam param1;
	308	Int_t hitNumber;
	309	Double_t zHit;
	310	Chi2 = 0.0; // initialize Chi2
	311	// copy of track parameters to be fitted
	312	param1 = *trackParamBeingFitted;
	313	// Minuit parameters to be fitted into this copy
	314	param1.SetInverseBendingMomentum(Param[0]);
	315	param1.SetBendingSlope(Param[1]);
	316	param1.SetNonBendingSlope(Param[2]);
	317	if (NParam == 5) {
	318	param1.SetNonBendingCoor(Param[3]);
	319	param1.SetBendingCoor(Param[4]);
	320	}
	321	// Follow track through all planes of track hits
	322	for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) {
	323	hit = (AliMUONTrackHit) ((trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
	324	zHit = hit->GetHitForRecPtr()->GetZ();
	325	// do something special if 2 hits with same Z ????
	326	// security against infinite loop ????
	327	(&param1)->ExtrapToZ(zHit); // extrapolation
	328	hit->SetTrackParam(&param1);
	329	// Increment Chi2
	330	// done hit per hit, with hit resolution,
	331	// and not with point and angle like in "reco_muon.F" !!!!
	332	// Needs to add multiple scattering contribution ????
	333	Double_t dX =
	334	hit->GetHitForRecPtr()->GetNonBendingCoor() - (&param1)->GetNonBendingCoor();
	335	Double_t dY =
	336	hit->GetHitForRecPtr()->GetBendingCoor() - (&param1)->GetBendingCoor();
	337	Chi2 =
	338	Chi2 +
	339	dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() +
	340	dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
	341	}
	342	}
	343
	344	//__________________________________________________________________________
	345	void TrackChi2MCS(Int_t &NParam, Double_t Gradient, Double_t &Chi2, Double_t Param, Int_t Flag)
	346	{
	347	// Return the "Chi2" to be minimized with Minuit for track fitting,
	348	// with "NParam" parameters
	349	// and their current values in array pointed to by "Param".
	350	// Assumes that the track hits are sorted according to increasing Z.
	351	// Track parameters at each TrackHit are updated accordingly.
	352	// Multiple Coulomb scattering is taken into account with covariance matrix.
	353	AliMUONTrackParam param1;
	354	Chi2 = 0.0; // initialize Chi2
	355	// copy of track parameters to be fitted
	356	param1 = *trackParamBeingFitted;
	357	// Minuit parameters to be fitted into this copy
	358	param1.SetInverseBendingMomentum(Param[0]);
	359	param1.SetBendingSlope(Param[1]);
	360	param1.SetNonBendingSlope(Param[2]);
	361	if (NParam == 5) {
	362	param1.SetNonBendingCoor(Param[3]);
	363	param1.SetBendingCoor(Param[4]);
	364	}
	365
	366	AliMUONTrackHit* hit, hit1, hit2, hit3;
	367	Bool_t GoodDeterminant;
	368	Int_t hitNumber, hitNumber1, hitNumber2, hitNumber3;
	369	Double_t zHit[10], paramBendingCoor[10], paramNonBendingCoor[10], ap[10];
	370	Double_t hitBendingCoor[10], hitNonBendingCoor[10];
	371	Double_t hitBendingReso2[10], hitNonBendingReso2[10];
	372	Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10);
	373	TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
	374	TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
	375
	376	// Predicted coordinates and multiple scattering angles are first calculated
	377	for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
	378	hit = (AliMUONTrackHit) ((trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
	379	zHit[hitNumber] = hit->GetHitForRecPtr()->GetZ();
	380	// do something special if 2 hits with same Z ????
	381	// security against infinite loop ????
	382	(&param1)->ExtrapToZ(zHit[hitNumber]); // extrapolation
	383	hit->SetTrackParam(&param1);
	384	paramBendingCoor[hitNumber]= (&param1)->GetBendingCoor();
	385	paramNonBendingCoor[hitNumber]= (&param1)->GetNonBendingCoor();
	386	hitBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetBendingCoor();
	387	hitNonBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingCoor();
	388	hitBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetBendingReso2();
	389	hitNonBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingReso2();
	390	ap[hitNumber] = MultipleScatteringAngle2(hit); // multiple scatt. angle ^2
	391	}
	392
	393	// Calculates the covariance matrix
	394	for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {
	395	for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
	396	(*covBending)(hitNumber1, hitNumber2) = 0;
	397	(*covBending)(hitNumber2, hitNumber1) = 0;
	398	if (hitNumber1 == hitNumber2){ // diagonal elements
	399	(*covBending)(hitNumber2, hitNumber1) =
	400	(*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
	401	}
	402	// Multiple Scattering... loop on upstream chambers ??
	403	for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++){
	404	(*covBending)(hitNumber2, hitNumber1) =
	405	(*covBending)(hitNumber2, hitNumber1) +
	406	((zHit[hitNumber1] - zHit[hitNumber3]) *
	407	(zHit[hitNumber2] - zHit[hitNumber3]) * ap[hitNumber3]);
	408	}
	409	(*covNonBending)(hitNumber1, hitNumber2) = 0;
	410	(*covNonBending)(hitNumber2, hitNumber1) =
	411	(*covBending)(hitNumber2, hitNumber1);
	412	if (hitNumber1 == hitNumber2) { // diagonal elements
	413	(*covNonBending)(hitNumber2, hitNumber1) =
	414	(*covNonBending)(hitNumber2, hitNumber1) -
	415	hitBendingReso2[hitNumber1] + hitNonBendingReso2[hitNumber1] ;
	416	}
	417	}
	418	}
	419
	420	// Inverts covariance matrix
	421	GoodDeterminant = kTRUE;
	422	if (covBending->Determinant() != 0) {
	423	covBending->Invert();
	424	} else {
	425	GoodDeterminant = kFALSE;
	426	cout << "Warning in ChiMCS Determinant Bending=0: " << endl;
	427	}
	428	if (covNonBending->Determinant() != 0){
	429	covNonBending->Invert();
	430	} else {
	431	GoodDeterminant = kFALSE;
	432	cout << "Warning in ChiMCS Determinant non Bending=0: " << endl;
	433	}
	434
	435	// Calculates Chi2
	436	if (GoodDeterminant) { // with Multiple Scattering if inversion correct
	437	for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++){
	438	for (hitNumber2=0; hitNumber2 < numberOfHit; hitNumber2++){
	439	Chi2 = Chi2 +
	440	((covBending)(hitNumber2, hitNumber1)
	441	(hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
	442	(hitBendingCoor[hitNumber2] - paramBendingCoor[hitNumber2]));
	443	Chi2 = Chi2 +
	444	((covNonBending)(hitNumber2, hitNumber1)
	445	(hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
	446	(hitNonBendingCoor[hitNumber2] - paramNonBendingCoor[hitNumber2]));
	447	}
	448	}
	449	} else { // without Multiple Scattering if inversion impossible
	450	for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++) {
	451	Chi2 = Chi2 +
	452	((hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
	453	(hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) /
	454	hitBendingReso2[hitNumber1]);
	455	Chi2 = Chi2 +
	456	((hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
	457	(hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) /
	458	hitNonBendingReso2[hitNumber1]);
	459	}
	460	}
	461
	462	delete covBending;
	463	delete covNonBending;
	464	}
	465
	466	Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
	467	{
	468	// Returns square of multiple Coulomb scattering angle
	469	// at TrackHit pointed to by "TrackHit"
	470	Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
	471	Double_t varMultipleScatteringAngle;
	472	AliMUONTrackParam *param = TrackHit->GetTrackParam();
	473	slopeBending = param->GetBendingSlope();
	474	slopeNonBending = param->GetNonBendingSlope();
	475	// thickness in radiation length for the current track,
	476	// taking local angle into account
	477	radiationLength =
	478	trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() *
	479	TMath::Sqrt(1.0 +
	480	slopeBending * slopeBending + slopeNonBending * slopeNonBending);
	481	inverseBendingMomentum2 =
	482	param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
	483	inverseTotalMomentum2 =
	484	inverseBendingMomentum2 * (1.0 + slopeBending*slopeBending) /
	485	(1.0 + slopeBending slopeBending + slopeNonBending slopeNonBending);
	486	varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
	487	varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
	488	varMultipleScatteringAngle * varMultipleScatteringAngle;
	489	return varMultipleScatteringAngle;
	490	}