[u/mrichter/AliRoot.git] / EVE / Alieve / MUONTrack.cxx

#include "MUONTrack.h"

#include <AliMagF.h>
#include <AliLog.h>

#include <AliMUONTrack.h>
#include <AliMUONTrackParam.h>
#include <AliMUONConstants.h>

#include <TClonesArray.h>
#include <TMath.h>
#include <TMatrixD.h>
#include <TStyle.h>
#include <TROOT.h>

using namespace Reve;
using namespace Alieve;

//______________________________________________________________________
// MUONTrack
// Produce Reve:Track from AliMUONTrack with dipole field model

ClassImp(MUONTrack)

//______________________________________________________________________
MUONTrack::MUONTrack(Reve::RecTrack* t, TrackRnrStyle* rs) :
  Reve::Track(t,rs),
  fFieldMap(0),
  fTrack(0),
  fCount(0)
{
  //
  // constructor
  //

}

//______________________________________________________________________
MUONTrack::~MUONTrack()
{
  //
  // destructor
  //

}

//______________________________________________________________________
void MUONTrack::MUONTrackInfo()
{
  //
  // MENU function
  //

  Reve::LoadMacro("MUON_track_info.C");
  gROOT->ProcessLine(Form("MUON_track_info(%d);", fLabel));

}

//______________________________________________________________________
void MUONTrack::MUONTriggerInfo()
{
  //
  // MENU function
  //

  Reve::LoadMacro("MUON_trigger_info.C");
  gROOT->ProcessLine(Form("MUON_trigger_info(%d);", fLabel));

}

//______________________________________________________________________
void MUONTrack::MakeTrack(AliMUONTrack *mtrack, AliMagF *fmap)
{
  //
  // builds the track with dipole field
  //

  fTrack    = mtrack;
  fFieldMap = fmap;

  Double_t xv, yv;
  Float_t ax, bx, ay, by;
  Float_t xr[20], yr[20], zr[20];
  Float_t xrc[20], yrc[20], zrc[20];
  char form[1000];
    
  TMatrixD smatrix(2,2);
  TMatrixD sums(2,1);
  TMatrixD res(2,1);

  Float_t xRec, xRec0;
  Float_t yRec, yRec0;
  Float_t zRec, zRec0;
  
  // middle z between the two detector planes of the trigger chambers
  Float_t zg[4] = { -1603.5, -1620.5, -1703.5, -1720.5 };

  AliMUONTrackParam *mtp = (AliMUONTrackParam*)mtrack->GetTrackParamAtVertex();
  Float_t pt = TMath::Sqrt(mtp->Px()*mtp->Px()+mtp->Py()*mtp->Py());

  //PH The line below is replaced waiting for a fix in Root
  //PH which permits to use variable siza arguments in CINT
  //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
  if (mtrack->GetMatchTrigger()) {
    //PH      track->SetName(Form("MUONTrack %2d (MT)", fLabel));
    sprintf(form,"MUONTrack %2d (MT)", fLabel);
    SetName(form);
    SetLineStyle(1);
    //SetLineColor(2);
    SetLineColor(ColorIndex(pt));
  } else {
    //PH      track->SetName(Form("MUONTrack %2d     ", fLabel));
    sprintf(form,"MUONTrack %2d     ", fLabel);
    SetName(form);
    SetLineStyle(2);
    //SetLineColor(2);
    SetLineColor(ColorIndex(pt));
  }
  
  AliMUONTrackParam *trackParam = mtrack->GetTrackParamAtVertex(); 
  xRec0  = trackParam->GetNonBendingCoor();
  yRec0  = trackParam->GetBendingCoor();
  zRec0  = trackParam->GetZ();
  
  SetPoint(fCount,xRec0,yRec0,zRec0);
  fCount++;

  for (Int_t i = 0; i < 20; i++) xr[i]=yr[i]=zr[i]=0.0;
  
  Int_t nTrackHits = mtrack->GetNTrackHits();
  
  Bool_t hitChamber[10] = {kFALSE};

  TClonesArray* trackParamAtHit = mtrack->GetTrackParamAtHit();

  for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
    trackParam = (AliMUONTrackParam*) trackParamAtHit->At(iHit); 
    xRec  = trackParam->GetNonBendingCoor();
    yRec  = trackParam->GetBendingCoor();
    zRec  = trackParam->GetZ();
    
    //printf("Hit %d x %f y %f z %f \n",iHit,xRec,yRec,zRec);
    
    xr[iHit] = xRec;
    yr[iHit] = yRec;
    zr[iHit] = zRec;
    
    hitChamber[AliMUONConstants::ChamberNumber(zRec)] = kTRUE;

  }

  Int_t crntCha, lastHitSt12, firstHitSt3, lastHitSt3, firstHitSt45;

  lastHitSt12  = -1;
  firstHitSt3  = -1;
  lastHitSt3   = -1;
  firstHitSt45 = -1;
  for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
    crntCha = AliMUONConstants::ChamberNumber(zr[iHit]);
    if (hitChamber[crntCha] && crntCha >= 0 && crntCha <= 3) {
      lastHitSt12 = iHit;
    }
    if (hitChamber[crntCha] && crntCha >= 4 && crntCha <= 5) {
      if (firstHitSt3 == -1) firstHitSt3 = iHit;
      lastHitSt3 = iHit;
    }
    if (hitChamber[crntCha] && crntCha >= 6 && crntCha <= 9) {
      if (firstHitSt45 == -1) firstHitSt45 = iHit;
    }
  }

  if (lastHitSt12 >= 0) {
    for (Int_t iHit = 0; iHit <= lastHitSt12; iHit++) {
      SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
      fCount++;
    }
    if (firstHitSt3 >= 0) {
      Propagate(xr,yr,zr,lastHitSt12,firstHitSt3);
      SetPoint(fCount,xr[firstHitSt3],yr[firstHitSt3],zr[firstHitSt3]);
      fCount++;
      if (lastHitSt3 >= 0) {
	SetPoint(fCount,xr[lastHitSt3],yr[lastHitSt3],zr[lastHitSt3]);
	fCount++;
	if (firstHitSt45 >= 0) {
	  Propagate(xr,yr,zr,lastHitSt3,firstHitSt45);
	  for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
	    SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
	    fCount++;
	  }
	} else {
	  Propagate(xr,yr,zr,lastHitSt3,9999);
	}
      } else if (firstHitSt45 >= 0) {
	Propagate(xr,yr,zr,firstHitSt3,firstHitSt45);
	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
	  SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
	  fCount++;
	}
      } else {
	Propagate(xr,yr,zr,firstHitSt3,9999);
      }
    } else if (lastHitSt3 >= 0) {
      Propagate(xr,yr,zr,lastHitSt12,lastHitSt3);
      SetPoint(fCount,xr[lastHitSt3],yr[lastHitSt3],zr[lastHitSt3]);
      fCount++;
      if (firstHitSt45 >= 0) {
	Propagate(xr,yr,zr,lastHitSt3,firstHitSt45);
	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
	  SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
	  fCount++;
	}
      } else {
	Propagate(xr,yr,zr,lastHitSt3,9999);
      }
    } else if (firstHitSt45 >= 0){
      Propagate(xr,yr,zr,lastHitSt12,firstHitSt45);
      for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
	fCount++;
      }
    } else {
      Propagate(xr,yr,zr,lastHitSt12,9999);
    }
  } else if (firstHitSt3 >= 0) {
    SetPoint(fCount,xr[firstHitSt3],yr[firstHitSt3],zr[firstHitSt3]);
    fCount++;
    if (lastHitSt3 >= 0) {
      SetPoint(fCount,xr[lastHitSt3],yr[lastHitSt3],zr[lastHitSt3]);
      fCount++;
      if (firstHitSt45) {
	Propagate(xr,yr,zr,lastHitSt3,firstHitSt45);
	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
	  SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
	  fCount++;
	}
      } else {
	Propagate(xr,yr,zr,lastHitSt3,9999);
      }
    } else if (firstHitSt45 >= 0) {
      Propagate(xr,yr,zr,firstHitSt3,firstHitSt45);
      for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
	fCount++;
      }
    } else {
      Propagate(xr,yr,zr,firstHitSt3,9999);
    }
  } else if (lastHitSt3 >= 0) {
    SetPoint(fCount,xr[lastHitSt3],yr[lastHitSt3],zr[lastHitSt3]);
    fCount++;
    if (firstHitSt45 >= 0) {
      Propagate(xr,yr,zr,lastHitSt3,firstHitSt45);
      for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
	fCount++;
      }
    } else {
      Propagate(xr,yr,zr,lastHitSt3,9999);
    }
  } else {
    for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
      SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
      fCount++;
    }
  }

  Int_t nrc = 0;
  if (mtrack->GetMatchTrigger() && 1) {
    
    for (Int_t i = 0; i < nTrackHits; i++) {
      if (TMath::Abs(zr[i]) > 1000.0) {
	//printf("Hit %d x %f y %f z %f \n",iHit,xr[i],yr[i],zr[i]);
	xrc[nrc] = xr[i];
	yrc[nrc] = yr[i];
	zrc[nrc] = zr[i];
	nrc++;
      }
    }
    
    if (nrc < 2) return;
    
    // fit x-z
    smatrix.Zero();
    sums.Zero();
    for (Int_t i = 0; i < nrc; i++) {
      xv = (Double_t)zrc[i];
      yv = (Double_t)xrc[i];
      //printf("x-z: xv %f yv %f \n",xv,yv);
      smatrix(0,0) += 1.0;
      smatrix(1,1) += xv*xv;
      smatrix(0,1) += xv;
      smatrix(1,0) += xv;
      sums(0,0)    += yv;
      sums(1,0)    += xv*yv;
    }
    res = smatrix.Invert() * sums;
    ax = res(0,0);
    bx = res(1,0);
    
    // fit y-z
    smatrix.Zero();
    sums.Zero();
    for (Int_t i = 0; i < nrc; i++) {
      xv = (Double_t)zrc[i];
      yv = (Double_t)yrc[i];
      //printf("y-z: xv %f yv %f \n",xv,yv);
      smatrix(0,0) += 1.0;
      smatrix(1,1) += xv*xv;
      smatrix(0,1) += xv;
      smatrix(1,0) += xv;
      sums(0,0)    += yv;
      sums(1,0)    += xv*yv;
    }
    res = smatrix.Invert() * sums;
    ay = res(0,0);
    by = res(1,0);
    
    Float_t xtc, ytc, ztc;
    for (Int_t ii = 0; ii < 4; ii++) {
      
      ztc = zg[ii];
      ytc = ay+by*zg[ii];
      xtc = ax+bx*zg[ii];
      
      //printf("tc: x %f y %f z %f \n",xtc,ytc,ztc);
      
      SetPoint(fCount,xtc,ytc,ztc);
      fCount++;
      
    }
    
  }  // end match trigger

}

//______________________________________________________________________
void MUONTrack::Propagate(Float_t *xr, Float_t *yr, Float_t *zr, Int_t i1, Int_t i2)
{
  //
  // propagate in magnetic field between hits of indices i1 and i2
  //

  Double_t vect[7], vout[7];
  Double_t step = 1.0;
  Double_t zMax;

  if (i2 == 9999) {
    //zMax = -1750.0;
    zMax = zr[i1]+2.0;
  } else {
    zMax = zr[i2]+2.0;
  }

  AliMUONTrackParam *trackParam = fTrack->GetTrackParamAtVertex(); 
  Int_t charge = (Int_t)TMath::Sign(1.0,trackParam->GetInverseBendingMomentum());
  
  TClonesArray* trackParamAtHit = fTrack->GetTrackParamAtHit();
  trackParam = (AliMUONTrackParam*)trackParamAtHit->At(i1); 

  vect[0] = xr[i1];
  vect[1] = yr[i1];
  vect[2] = zr[i1];
  vect[3] = trackParam->Px()/trackParam->P();
  vect[4] = trackParam->Py()/trackParam->P();
  vect[5] = trackParam->Pz()/trackParam->P();
  vect[6] = trackParam->P();

  Int_t nSteps = 0;
  while ((vect[2] > zMax) && (nSteps < 10000)) {
    nSteps++;
    OneStepRungekutta(charge, step, vect, vout);
    //printf("%f %f %f \n",vect[0],vect[1],vect[2]);
    SetPoint(fCount,vout[0],vout[1],vout[2]);
    fCount++;
    for (Int_t i = 0; i < 7; i++) {
      vect[i] = vout[i];
    }
  }

}

//______________________________________________________________________
void MUONTrack::GetField(Double_t *position, Double_t *field)
{
  // 
  // returns field components at position, for a give field map
  //

  /// interface for arguments in double precision (Why ? ChF)
  Float_t x[3], b[3];

  x[0] = position[0]; x[1] = position[1]; x[2] = position[2];

  if (fFieldMap) fFieldMap->Field(x,b);
  else {
    //AliWarning("No field map");
    field[0] = field[1] = field[2] = 0.0;
    return;
  }
  
  // force components
  //b[1] = 0.0;
  //b[2] = 0.0;

  field[0] = b[0]; field[1] = b[1]; field[2] = b[2];

  return;

}

//______________________________________________________________________
void MUONTrack::OneStepRungekutta(Double_t charge, Double_t step, 
				  Double_t* vect, Double_t* vout)
{
///	******************************************************************
///	*								 *
///	*  Runge-Kutta method for tracking a particle through a magnetic *
///	*  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of	 *
///	*  Standards, procedure 25.5.20)				 *
///	*								 *
///	*  Input parameters						 *
///	*	CHARGE    Particle charge				 *
///	*	STEP	  Step size					 *
///	*	VECT	  Initial co-ords,direction cosines,momentum	 *
///	*  Output parameters						 *
///	*	VOUT	  Output co-ords,direction cosines,momentum	 *
///	*  User routine called  					 *
///	*	CALL GUFLD(X,F) 					 *
///	*								 *
///	*    ==>Called by : <USER>, GUSWIM				 *
///	*	Authors    R.Brun, M.Hansroul  *********		 *
///	*		   V.Perevoztchikov (CUT STEP implementation)	 *
///	*								 *
///	*								 *
///	******************************************************************

    Double_t h2, h4, f[4];
    Double_t xyzt[3], a, b, c, ph,ph2;
    Double_t secxs[4],secys[4],seczs[4],hxp[3];
    Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
    Double_t est, at, bt, ct, cba;
    Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
    
    Double_t x;
    Double_t y;
    Double_t z;
    
    Double_t xt;
    Double_t yt;
    Double_t zt;

    Double_t maxit = 1992;
    Double_t maxcut = 11;

    const Double_t kdlt   = 1e-4;
    const Double_t kdlt32 = kdlt/32.;
    const Double_t kthird = 1./3.;
    const Double_t khalf  = 0.5;
    const Double_t kec = 2.9979251e-4;

    const Double_t kpisqua = 9.86960440109;
    const Int_t kix  = 0;
    const Int_t kiy  = 1;
    const Int_t kiz  = 2;
    const Int_t kipx = 3;
    const Int_t kipy = 4;
    const Int_t kipz = 5;
  
    // *.
    // *.    ------------------------------------------------------------------
    // *.
    // *             this constant is for units cm,gev/c and kgauss
    // *
    Int_t iter = 0;
    Int_t ncut = 0;
    for(Int_t j = 0; j < 7; j++)
      vout[j] = vect[j];

    Double_t  pinv   = kec * charge / vect[6];
    Double_t tl = 0.;
    Double_t h = step;
    Double_t rest;

 
    do {
      rest  = step - tl;
      if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
      //cmodif: call gufld(vout,f) changed into:

      GetField(vout,f);

      // *
      // *             start of integration
      // *
      x      = vout[0];
      y      = vout[1];
      z      = vout[2];
      a      = vout[3];
      b      = vout[4];
      c      = vout[5];

      h2     = khalf * h;
      h4     = khalf * h2;
      ph     = pinv * h;
      ph2    = khalf * ph;
      secxs[0] = (b * f[2] - c * f[1]) * ph2;
      secys[0] = (c * f[0] - a * f[2]) * ph2;
      seczs[0] = (a * f[1] - b * f[0]) * ph2;
      ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
      if (ang2 > kpisqua) break;

      dxt    = h2 * a + h4 * secxs[0];
      dyt    = h2 * b + h4 * secys[0];
      dzt    = h2 * c + h4 * seczs[0];
      xt     = x + dxt;
      yt     = y + dyt;
      zt     = z + dzt;
      // *
      // *              second intermediate point
      // *

      est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
      if (est > h) {
	if (ncut++ > maxcut) break;
	h *= khalf;
	continue;
      }
 
      xyzt[0] = xt;
      xyzt[1] = yt;
      xyzt[2] = zt;

      //cmodif: call gufld(xyzt,f) changed into:
      GetField(xyzt,f);

      at     = a + secxs[0];
      bt     = b + secys[0];
      ct     = c + seczs[0];

      secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
      secys[1] = (ct * f[0] - at * f[2]) * ph2;
      seczs[1] = (at * f[1] - bt * f[0]) * ph2;
      at     = a + secxs[1];
      bt     = b + secys[1];
      ct     = c + seczs[1];
      secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
      secys[2] = (ct * f[0] - at * f[2]) * ph2;
      seczs[2] = (at * f[1] - bt * f[0]) * ph2;
      dxt    = h * (a + secxs[2]);
      dyt    = h * (b + secys[2]);
      dzt    = h * (c + seczs[2]);
      xt     = x + dxt;
      yt     = y + dyt;
      zt     = z + dzt;
      at     = a + 2.*secxs[2];
      bt     = b + 2.*secys[2];
      ct     = c + 2.*seczs[2];

      est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
      if (est > 2.*TMath::Abs(h)) {
	if (ncut++ > maxcut) break;
	h *= khalf;
	continue;
      }
 
      xyzt[0] = xt;
      xyzt[1] = yt;
      xyzt[2] = zt;

      //cmodif: call gufld(xyzt,f) changed into:
      GetField(xyzt,f);

      z      = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
      y      = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
      x      = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;

      secxs[3] = (bt*f[2] - ct*f[1])* ph2;
      secys[3] = (ct*f[0] - at*f[2])* ph2;
      seczs[3] = (at*f[1] - bt*f[0])* ph2;
      a      = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
      b      = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
      c      = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;

      est    = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
	+ TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
	+ TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));

      if (est > kdlt && TMath::Abs(h) > 1.e-4) {
	if (ncut++ > maxcut) break;
	h *= khalf;
	continue;
      }

      ncut = 0;
      // *               if too many iterations, go to helix
      if (iter++ > maxit) break;

      tl += h;
      if (est < kdlt32) 
	h *= 2.;
      cba    = 1./ TMath::Sqrt(a*a + b*b + c*c);
      vout[0] = x;
      vout[1] = y;
      vout[2] = z;
      vout[3] = cba*a;
      vout[4] = cba*b;
      vout[5] = cba*c;
      rest = step - tl;
      if (step < 0.) rest = -rest;
      if (rest < 1.e-5*TMath::Abs(step)) return;

    } while(1);

    // angle too big, use helix

    f1  = f[0];
    f2  = f[1];
    f3  = f[2];
    f4  = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
    rho = -f4*pinv;
    tet = rho * step;
 
    hnorm = 1./f4;
    f1 = f1*hnorm;
    f2 = f2*hnorm;
    f3 = f3*hnorm;

    hxp[0] = f2*vect[kipz] - f3*vect[kipy];
    hxp[1] = f3*vect[kipx] - f1*vect[kipz];
    hxp[2] = f1*vect[kipy] - f2*vect[kipx];
 
    hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];

    rho1 = 1./rho;
    sint = TMath::Sin(tet);
    cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);

    g1 = sint*rho1;
    g2 = cost*rho1;
    g3 = (tet-sint) * hp*rho1;
    g4 = -cost;
    g5 = sint;
    g6 = cost * hp;
 
    vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
    vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
    vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
 
    vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
    vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
    vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;

    return;
}

//______________________________________________________________________
Int_t MUONTrack::ColorIndex(Float_t val)
{
  //
  // returns color index in the palette for a give value
  //

  Float_t threshold =  0.0;
  Float_t maxVal    =  2.0;

  Float_t div  = TMath::Max(1, (Int_t)(maxVal - threshold));
  Int_t   nCol = gStyle->GetNumberOfColors();
  Int_t   cBin = (Int_t) TMath::Nint(nCol*(val - threshold)/div);

  return gStyle->GetColorPalette(TMath::Min(nCol - 1, cBin));

}
Commit	Line	Data
e19f3684	1	#include "MUONTrack.h"
	2
	3	#include <AliMagF.h>
	4	#include <AliLog.h>
	5
	6	#include <AliMUONTrack.h>
	7	#include <AliMUONTrackParam.h>
	8	#include <AliMUONConstants.h>
	9
	10	#include <TClonesArray.h>
	11	#include <TMath.h>
	12	#include <TMatrixD.h>
	13	#include <TStyle.h>
f7e76528	14	#include <TROOT.h>
e19f3684	15
	16	using namespace Reve;
	17	using namespace Alieve;
	18
	19	//______________________________________________________________________
	20	// MUONTrack
	21	// Produce Reve:Track from AliMUONTrack with dipole field model
	22
	23	ClassImp(MUONTrack)
	24
	25	//______________________________________________________________________
	26	MUONTrack::MUONTrack(Reve::RecTrack* t, TrackRnrStyle* rs) :
	27	Reve::Track(t,rs),
	28	fFieldMap(0),
	29	fTrack(0),
	30	fCount(0)
	31	{
	32	//
	33	// constructor
	34	//
	35
	36	}
	37
	38	//______________________________________________________________________
	39	MUONTrack::~MUONTrack()
	40	{
	41	//
	42	// destructor
	43	//
	44
	45	}
	46
	47	//______________________________________________________________________
	48	void MUONTrack::MUONTrackInfo()
	49	{
	50	//
	51	// MENU function
	52	//
	53
	54	Reve::LoadMacro("MUON_track_info.C");
	55	gROOT->ProcessLine(Form("MUON_track_info(%d);", fLabel));
	56
	57	}
	58
	59	//______________________________________________________________________
	60	void MUONTrack::MUONTriggerInfo()
	61	{
	62	//
	63	// MENU function
	64	//
	65
	66	Reve::LoadMacro("MUON_trigger_info.C");
	67	gROOT->ProcessLine(Form("MUON_trigger_info(%d);", fLabel));
	68
	69	}
	70
	71	//______________________________________________________________________
	72	void MUONTrack::MakeTrack(AliMUONTrack mtrack, AliMagF fmap)
	73	{
	74	//
	75	// builds the track with dipole field
	76	//
	77
	78	fTrack = mtrack;
79	fFieldMap = fmap;
80
81	Double_t xv, yv;
82	Float_t ax, bx, ay, by;
83	Float_t xr[20], yr[20], zr[20];
84	Float_t xrc[20], yrc[20], zrc[20];
85	char form[1000];
86
87	TMatrixD smatrix(2,2);
88	TMatrixD sums(2,1);
89	TMatrixD res(2,1);
90
91	Float_t xRec, xRec0;
92	Float_t yRec, yRec0;
93	Float_t zRec, zRec0;
94
95	// middle z between the two detector planes of the trigger chambers
96	Float_t zg[4] = { -1603.5, -1620.5, -1703.5, -1720.5 };
97
98	AliMUONTrackParam mtp = (AliMUONTrackParam)mtrack->GetTrackParamAtVertex();
99	Float_t pt = TMath::Sqrt(mtp->Px()mtp->Px()+mtp->Py()mtp->Py());
100
101	//PH The line below is replaced waiting for a fix in Root
102	//PH which permits to use variable siza arguments in CINT
103	//PH on some platforms (alphalinuxgcc, solariscc5, etc.)
104	if (mtrack->GetMatchTrigger()) {
105	//PH track->SetName(Form("MUONTrack %2d (MT)", fLabel));
106	sprintf(form,"MUONTrack %2d (MT)", fLabel);
107	SetName(form);
108	SetLineStyle(1);
109	//SetLineColor(2);
110	SetLineColor(ColorIndex(pt));
111	} else {
112	//PH track->SetName(Form("MUONTrack %2d ", fLabel));
113	sprintf(form,"MUONTrack %2d ", fLabel);
114	SetName(form);
115	SetLineStyle(2);
116	//SetLineColor(2);
117	SetLineColor(ColorIndex(pt));
118	}
119
120	AliMUONTrackParam *trackParam = mtrack->GetTrackParamAtVertex();
121	xRec0 = trackParam->GetNonBendingCoor();
122	yRec0 = trackParam->GetBendingCoor();
123	zRec0 = trackParam->GetZ();
124
125	SetPoint(fCount,xRec0,yRec0,zRec0);
126	fCount++;
127
128	for (Int_t i = 0; i < 20; i++) xr[i]=yr[i]=zr[i]=0.0;
129
130	Int_t nTrackHits = mtrack->GetNTrackHits();
131
132	Bool_t hitChamber[10] = {kFALSE};
133
134	TClonesArray* trackParamAtHit = mtrack->GetTrackParamAtHit();
135
136	for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
137	trackParam = (AliMUONTrackParam*) trackParamAtHit->At(iHit);
138	xRec = trackParam->GetNonBendingCoor();
139	yRec = trackParam->GetBendingCoor();
140	zRec = trackParam->GetZ();
141
142	//printf("Hit %d x %f y %f z %f \n",iHit,xRec,yRec,zRec);
143
144	xr[iHit] = xRec;
145	yr[iHit] = yRec;
146	zr[iHit] = zRec;
147
148	hitChamber[AliMUONConstants::ChamberNumber(zRec)] = kTRUE;
149
150	}
151
152	Int_t crntCha, lastHitSt12, firstHitSt3, lastHitSt3, firstHitSt45;
153
154	lastHitSt12 = -1;
155	firstHitSt3 = -1;
156	lastHitSt3 = -1;
157	firstHitSt45 = -1;
158	for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
159	crntCha = AliMUONConstants::ChamberNumber(zr[iHit]);
160	if (hitChamber[crntCha] && crntCha >= 0 && crntCha <= 3) {
161	lastHitSt12 = iHit;
162	}
163	if (hitChamber[crntCha] && crntCha >= 4 && crntCha <= 5) {
164	if (firstHitSt3 == -1) firstHitSt3 = iHit;
165	lastHitSt3 = iHit;
166	}
167	if (hitChamber[crntCha] && crntCha >= 6 && crntCha <= 9) {
168	if (firstHitSt45 == -1) firstHitSt45 = iHit;
169	}
170	}
171
172	if (lastHitSt12 >= 0) {
173	for (Int_t iHit = 0; iHit <= lastHitSt12; iHit++) {
174	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
175	fCount++;
176	}
177	if (firstHitSt3 >= 0) {
178	Propagate(xr,yr,zr,lastHitSt12,firstHitSt3);
179	SetPoint(fCount,xr[firstHitSt3],yr[firstHitSt3],zr[firstHitSt3]);
180	fCount++;
181	if (lastHitSt3 >= 0) {
182	SetPoint(fCount,xr[lastHitSt3],yr[lastHitSt3],zr[lastHitSt3]);
183	fCount++;
184	if (firstHitSt45 >= 0) {
185	Propagate(xr,yr,zr,lastHitSt3,firstHitSt45);
186	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
187	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
188	fCount++;
189	}
190	} else {
191	Propagate(xr,yr,zr,lastHitSt3,9999);
192	}
193	} else if (firstHitSt45 >= 0) {
194	Propagate(xr,yr,zr,firstHitSt3,firstHitSt45);
195	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
196	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
197	fCount++;
198	}
199	} else {
200	Propagate(xr,yr,zr,firstHitSt3,9999);
201	}
202	} else if (lastHitSt3 >= 0) {
203	Propagate(xr,yr,zr,lastHitSt12,lastHitSt3);
204	SetPoint(fCount,xr[lastHitSt3],yr[lastHitSt3],zr[lastHitSt3]);
205	fCount++;
206	if (firstHitSt45 >= 0) {
207	Propagate(xr,yr,zr,lastHitSt3,firstHitSt45);
208	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
209	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
210	fCount++;
211	}
212	} else {
213	Propagate(xr,yr,zr,lastHitSt3,9999);
214	}
215	} else if (firstHitSt45 >= 0){
216	Propagate(xr,yr,zr,lastHitSt12,firstHitSt45);
217	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
218	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
219	fCount++;
220	}
221	} else {
222	Propagate(xr,yr,zr,lastHitSt12,9999);
223	}
224	} else if (firstHitSt3 >= 0) {
225	SetPoint(fCount,xr[firstHitSt3],yr[firstHitSt3],zr[firstHitSt3]);
226	fCount++;
227	if (lastHitSt3 >= 0) {
228	SetPoint(fCount,xr[lastHitSt3],yr[lastHitSt3],zr[lastHitSt3]);
229	fCount++;
230	if (firstHitSt45) {
231	Propagate(xr,yr,zr,lastHitSt3,firstHitSt45);
232	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
233	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
234	fCount++;
235	}
236	} else {
237	Propagate(xr,yr,zr,lastHitSt3,9999);
238	}
239	} else if (firstHitSt45 >= 0) {
240	Propagate(xr,yr,zr,firstHitSt3,firstHitSt45);
241	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
242	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
243	fCount++;
244	}
245	} else {
246	Propagate(xr,yr,zr,firstHitSt3,9999);
247	}
248	} else if (lastHitSt3 >= 0) {
249	SetPoint(fCount,xr[lastHitSt3],yr[lastHitSt3],zr[lastHitSt3]);
250	fCount++;
251	if (firstHitSt45 >= 0) {
252	Propagate(xr,yr,zr,lastHitSt3,firstHitSt45);
253	for (Int_t iHit = firstHitSt45; iHit < nTrackHits; iHit++) {
254	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
255	fCount++;
256	}
257	} else {
258	Propagate(xr,yr,zr,lastHitSt3,9999);
259	}
260	} else {
261	for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
262	SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
263	fCount++;
264	}
265	}
266
267	Int_t nrc = 0;
268	if (mtrack->GetMatchTrigger() && 1) {
269
270	for (Int_t i = 0; i < nTrackHits; i++) {
271	if (TMath::Abs(zr[i]) > 1000.0) {
272	//printf("Hit %d x %f y %f z %f \n",iHit,xr[i],yr[i],zr[i]);
273	xrc[nrc] = xr[i];
274	yrc[nrc] = yr[i];
275	zrc[nrc] = zr[i];
276	nrc++;
277	}
278	}
279
280	if (nrc < 2) return;
281
282	// fit x-z
283	smatrix.Zero();
284	sums.Zero();
285	for (Int_t i = 0; i < nrc; i++) {
286	xv = (Double_t)zrc[i];
287	yv = (Double_t)xrc[i];
288	//printf("x-z: xv %f yv %f \n",xv,yv);
289	smatrix(0,0) += 1.0;
290	smatrix(1,1) += xv*xv;
291	smatrix(0,1) += xv;
292	smatrix(1,0) += xv;
293	sums(0,0) += yv;
294	sums(1,0) += xv*yv;
295	}
296	res = smatrix.Invert() * sums;
297	ax = res(0,0);
298	bx = res(1,0);
299
300	// fit y-z
301	smatrix.Zero();
302	sums.Zero();
303	for (Int_t i = 0; i < nrc; i++) {
304	xv = (Double_t)zrc[i];
305	yv = (Double_t)yrc[i];
306	//printf("y-z: xv %f yv %f \n",xv,yv);
307	smatrix(0,0) += 1.0;
308	smatrix(1,1) += xv*xv;
309	smatrix(0,1) += xv;
310	smatrix(1,0) += xv;
311	sums(0,0) += yv;
312	sums(1,0) += xv*yv;
313	}
314	res = smatrix.Invert() * sums;
315	ay = res(0,0);
316	by = res(1,0);
317
318	Float_t xtc, ytc, ztc;
319	for (Int_t ii = 0; ii < 4; ii++) {
320
321	ztc = zg[ii];
322	ytc = ay+by*zg[ii];
323	xtc = ax+bx*zg[ii];
324
325	//printf("tc: x %f y %f z %f \n",xtc,ytc,ztc);
326
327	SetPoint(fCount,xtc,ytc,ztc);
328	fCount++;
329
330	}
331
332	} // end match trigger
333
334	}
335
336	//______________________________________________________________________
337	void MUONTrack::Propagate(Float_t xr, Float_t yr, Float_t *zr, Int_t i1, Int_t i2)
338	{
339	//
340	// propagate in magnetic field between hits of indices i1 and i2
341	//
342
343	Double_t vect[7], vout[7];
344	Double_t step = 1.0;
345	Double_t zMax;
346
347	if (i2 == 9999) {
348	//zMax = -1750.0;
349	zMax = zr[i1]+2.0;
350	} else {
351	zMax = zr[i2]+2.0;
352	}
353
354	AliMUONTrackParam *trackParam = fTrack->GetTrackParamAtVertex();
355	Int_t charge = (Int_t)TMath::Sign(1.0,trackParam->GetInverseBendingMomentum());
356
357	TClonesArray* trackParamAtHit = fTrack->GetTrackParamAtHit();
358	trackParam = (AliMUONTrackParam*)trackParamAtHit->At(i1);
359
360	vect[0] = xr[i1];
361	vect[1] = yr[i1];
362	vect[2] = zr[i1];
363	vect[3] = trackParam->Px()/trackParam->P();
364	vect[4] = trackParam->Py()/trackParam->P();
365	vect[5] = trackParam->Pz()/trackParam->P();
366	vect[6] = trackParam->P();
367
368	Int_t nSteps = 0;
369	while ((vect[2] > zMax) && (nSteps < 10000)) {
370	nSteps++;
371	OneStepRungekutta(charge, step, vect, vout);
372	//printf("%f %f %f \n",vect[0],vect[1],vect[2]);
373	SetPoint(fCount,vout[0],vout[1],vout[2]);
374	fCount++;
375	for (Int_t i = 0; i < 7; i++) {
376	vect[i] = vout[i];
377	}
378	}
379
380	}
381
382	//______________________________________________________________________
383	void MUONTrack::GetField(Double_t position, Double_t field)
384	{
385	//
386	// returns field components at position, for a give field map
387	//
388
389	/// interface for arguments in double precision (Why ? ChF)
390	Float_t x[3], b[3];
391
392	x[0] = position[0]; x[1] = position[1]; x[2] = position[2];
393
394	if (fFieldMap) fFieldMap->Field(x,b);
395	else {
396	//AliWarning("No field map");
397	field[0] = field[1] = field[2] = 0.0;
398	return;
399	}
400
401	// force components
402	//b[1] = 0.0;
403	//b[2] = 0.0;
404
405	field[0] = b[0]; field[1] = b[1]; field[2] = b[2];
406
407	return;
408
409	}
410
411	//______________________________________________________________________
412	void MUONTrack::OneStepRungekutta(Double_t charge, Double_t step,
413	Double_t* vect, Double_t* vout)
414	{
415	/// ******************************************************************
416	/// * *
417	/// * Runge-Kutta method for tracking a particle through a magnetic *
418	/// * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
419	/// * Standards, procedure 25.5.20) *
420	/// * *
421	/// * Input parameters *
422	/// * CHARGE Particle charge *
423	/// * STEP Step size *
424	/// * VECT Initial co-ords,direction cosines,momentum *
425	/// * Output parameters *
426	/// * VOUT Output co-ords,direction cosines,momentum *
427	/// * User routine called *
428	/// * CALL GUFLD(X,F) *
429	/// * *
430	/// * ==>Called by : <USER>, GUSWIM *
431	/// * Authors R.Brun, M.Hansroul ********* *
432	/// * V.Perevoztchikov (CUT STEP implementation) *
433	/// * *
434	/// * *
435	/// ******************************************************************
436
437	Double_t h2, h4, f[4];
438	Double_t xyzt[3], a, b, c, ph,ph2;
439	Double_t secxs[4],secys[4],seczs[4],hxp[3];
440	Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
441	Double_t est, at, bt, ct, cba;
442	Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
443
444	Double_t x;
445	Double_t y;
446	Double_t z;
447
448	Double_t xt;
449	Double_t yt;
450	Double_t zt;
451
452	Double_t maxit = 1992;
453	Double_t maxcut = 11;
454
455	const Double_t kdlt = 1e-4;
456	const Double_t kdlt32 = kdlt/32.;
457	const Double_t kthird = 1./3.;
458	const Double_t khalf = 0.5;
459	const Double_t kec = 2.9979251e-4;
460
461	const Double_t kpisqua = 9.86960440109;
462	const Int_t kix = 0;
463	const Int_t kiy = 1;
464	const Int_t kiz = 2;
465	const Int_t kipx = 3;
466	const Int_t kipy = 4;
467	const Int_t kipz = 5;
468
469	// *.
470	// *. ------------------------------------------------------------------
471	// *.
472	// * this constant is for units cm,gev/c and kgauss
473	// *
474	Int_t iter = 0;
475	Int_t ncut = 0;
476	for(Int_t j = 0; j < 7; j++)
477	vout[j] = vect[j];
478
479	Double_t pinv = kec * charge / vect[6];
480	Double_t tl = 0.;
481	Double_t h = step;
482	Double_t rest;
483
484
485	do {
486	rest = step - tl;
487	if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
488	//cmodif: call gufld(vout,f) changed into:
489
490	GetField(vout,f);
491
492	// *
493	// * start of integration
494	// *
495	x = vout[0];
496	y = vout[1];
497	z = vout[2];
498	a = vout[3];
499	b = vout[4];
500	c = vout[5];
501
502	h2 = khalf * h;
503	h4 = khalf * h2;
504	ph = pinv * h;
505	ph2 = khalf * ph;
506	secxs[0] = (b * f[2] - c * f[1]) * ph2;
507	secys[0] = (c * f[0] - a * f[2]) * ph2;
508	seczs[0] = (a * f[1] - b * f[0]) * ph2;
509	ang2 = (secxs[0]secxs[0] + secys[0]secys[0] + seczs[0]*seczs[0]);
510	if (ang2 > kpisqua) break;
511
512	dxt = h2 * a + h4 * secxs[0];
513	dyt = h2 * b + h4 * secys[0];
514	dzt = h2 * c + h4 * seczs[0];
515	xt = x + dxt;
516	yt = y + dyt;
517	zt = z + dzt;
518	// *
519	// * second intermediate point
520	// *
521
522	est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
523	if (est > h) {
524	if (ncut++ > maxcut) break;
525	h *= khalf;
526	continue;
527	}
528
529	xyzt[0] = xt;
530	xyzt[1] = yt;
531	xyzt[2] = zt;
532
533	//cmodif: call gufld(xyzt,f) changed into:
534	GetField(xyzt,f);
535
536	at = a + secxs[0];
537	bt = b + secys[0];
538	ct = c + seczs[0];
539
540	secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
541	secys[1] = (ct * f[0] - at * f[2]) * ph2;
542	seczs[1] = (at * f[1] - bt * f[0]) * ph2;
543	at = a + secxs[1];
544	bt = b + secys[1];
545	ct = c + seczs[1];
546	secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
547	secys[2] = (ct * f[0] - at * f[2]) * ph2;
548	seczs[2] = (at * f[1] - bt * f[0]) * ph2;
549	dxt = h * (a + secxs[2]);
550	dyt = h * (b + secys[2]);
551	dzt = h * (c + seczs[2]);
552	xt = x + dxt;
553	yt = y + dyt;
554	zt = z + dzt;
555	at = a + 2.*secxs[2];
556	bt = b + 2.*secys[2];
557	ct = c + 2.*seczs[2];
558
559	est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
560	if (est > 2.*TMath::Abs(h)) {
561	if (ncut++ > maxcut) break;
562	h *= khalf;
563	continue;
564	}
565
566	xyzt[0] = xt;
567	xyzt[1] = yt;
568	xyzt[2] = zt;
569
570	//cmodif: call gufld(xyzt,f) changed into:
571	GetField(xyzt,f);
572
573	z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
574	y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
575	x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
576
577	secxs[3] = (btf[2] - ctf[1])* ph2;
578	secys[3] = (ctf[0] - atf[2])* ph2;
579	seczs[3] = (atf[1] - btf[0])* ph2;
580	a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
581	b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
582	c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
583
584	est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
585	+ TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
586	+ TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
587
588	if (est > kdlt && TMath::Abs(h) > 1.e-4) {
589	if (ncut++ > maxcut) break;
590	h *= khalf;
591	continue;
592	}
593
594	ncut = 0;
595	// * if too many iterations, go to helix
596	if (iter++ > maxit) break;
597
598	tl += h;
599	if (est < kdlt32)
600	h *= 2.;
601	cba = 1./ TMath::Sqrt(aa + bb + c*c);
602	vout[0] = x;
603	vout[1] = y;
604	vout[2] = z;
605	vout[3] = cba*a;
606	vout[4] = cba*b;
607	vout[5] = cba*c;
608	rest = step - tl;
609	if (step < 0.) rest = -rest;
610	if (rest < 1.e-5*TMath::Abs(step)) return;
611
612	} while(1);
613
614	// angle too big, use helix
615
616	f1 = f[0];
617	f2 = f[1];
618	f3 = f[2];
619	f4 = TMath::Sqrt(f1f1+f2f2+f3*f3);
620	rho = -f4*pinv;
621	tet = rho * step;
622
623	hnorm = 1./f4;
624	f1 = f1*hnorm;
625	f2 = f2*hnorm;
626	f3 = f3*hnorm;
627
628	hxp[0] = f2vect[kipz] - f3vect[kipy];
629	hxp[1] = f3vect[kipx] - f1vect[kipz];
630	hxp[2] = f1vect[kipy] - f2vect[kipx];
631
632	hp = f1vect[kipx] + f2vect[kipy] + f3*vect[kipz];
633
634	rho1 = 1./rho;
635	sint = TMath::Sin(tet);
636	cost = 2.TMath::Sin(khalftet)TMath::Sin(khalftet);
637
638	g1 = sint*rho1;
639	g2 = cost*rho1;
640	g3 = (tet-sint) * hp*rho1;
641	g4 = -cost;
642	g5 = sint;
643	g6 = cost * hp;
644
645	vout[kix] = vect[kix] + g1vect[kipx] + g2hxp[0] + g3*f1;
646	vout[kiy] = vect[kiy] + g1vect[kipy] + g2hxp[1] + g3*f2;
647	vout[kiz] = vect[kiz] + g1vect[kipz] + g2hxp[2] + g3*f3;
648
649	vout[kipx] = vect[kipx] + g4vect[kipx] + g5hxp[0] + g6*f1;
650	vout[kipy] = vect[kipy] + g4vect[kipy] + g5hxp[1] + g6*f2;
651	vout[kipz] = vect[kipz] + g4vect[kipz] + g5hxp[2] + g6*f3;
652
653	return;
654	}
655
656	//______________________________________________________________________
657	Int_t MUONTrack::ColorIndex(Float_t val)
658	{
659	//
660	// returns color index in the palette for a give value
661	//
662
663	Float_t threshold = 0.0;
664	Float_t maxVal = 2.0;
665
666	Float_t div = TMath::Max(1, (Int_t)(maxVal - threshold));
667	Int_t nCol = gStyle->GetNumberOfColors();
668	Int_t cBin = (Int_t) TMath::Nint(nCol*(val - threshold)/div);
669
670	return gStyle->GetColorPalette(TMath::Min(nCol - 1, cBin));
671
672	}