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

///////////////////////////////////////////////////////
//                                                   //
//                                                   //
//  Tracklet object created in the local tracking    //
//                                                   //
//                                                   //
///////////////////////////////////////////////////////

#include <TGraph.h>
#include <TMath.h>
#include <TF1.h>

#include "AliLog.h"

#include "AliTRDcalibDB.h"
#include "AliTRDCommonParam.h"
#include "AliTRDpadPlane.h"
#include "AliTRDgeometry.h"
#include "AliTRDmcmTracklet.h"

ClassImp(AliTRDmcmTracklet)

//_____________________________________________________________________________
AliTRDmcmTracklet::AliTRDmcmTracklet() 
  :TObject()
  ,fGeo(0)
  ,fDetector(-1)
  ,fRow(-1)
  ,fTrackLabel(-1)
  ,fNclusters(0)
  ,fN(0)
  ,fGPos(0)
  ,fGAmp(0)
  ,fTime0(0)
  ,fRowz(0)
  ,fSlope(0)
  ,fOffset(0)
  ,fPt(0)
  ,fdQdl(0)
{
  //
  // AliTRDmcmTracklet default constructor
  //

  for (Int_t time = 0; time < kNtimeBins; time++) {
    for (Int_t icl = 0; icl < kNclsPads; icl++) {
      fADC[time][icl] = 0;
    }
    for (Int_t it = 0; it < kNdict; it++) {
      fTrack[time][it] = -1;
    }
    fTime[time]   = 0;
    fCol[time]    = 0;
  }

  fGeo = new AliTRDgeometry();

}

//_____________________________________________________________________________
AliTRDmcmTracklet::AliTRDmcmTracklet(Int_t det, Int_t row, Int_t n) 
  :TObject()
  ,fGeo(0)
  ,fDetector(det)
  ,fRow(row)
  ,fTrackLabel(-1)
  ,fNclusters(0)
  ,fN(n)
  ,fGPos(0)
  ,fGAmp(0)
  ,fTime0(0)
  ,fRowz(0)
  ,fSlope(0)
  ,fOffset(0)
  ,fPt(0)
  ,fdQdl(0)
{
  //
  // AliTRDmcmTracklet default constructor
  //

  for (Int_t time = 0; time < kNtimeBins; time++) {
    for (Int_t icl = 0; icl < kNclsPads; icl++) {
      fADC[time][icl] = 0;
    }
    for (Int_t it = 0; it < kNdict; it++) {
      fTrack[time][it] = -1;
    }
    fTime[time]   = 0;
    fCol[time]    = 0;
  }

  fGPos = new TGraph(0);
  fGAmp = new TGraph(0);

  fGeo  = new AliTRDgeometry();

}

//_____________________________________________________________________________
AliTRDmcmTracklet::AliTRDmcmTracklet(const AliTRDmcmTracklet &t) 
  :TObject(t)
  ,fGeo(0)
  ,fDetector(t.fDetector)
  ,fRow(t.fRow)
  ,fTrackLabel(t.fTrackLabel)
  ,fNclusters(t.fNclusters)
  ,fN(t.fN)
  ,fGPos(NULL)
  ,fGAmp(NULL)
  ,fTime0(t.fTime0)
  ,fRowz(t.fRowz)
  ,fSlope(t.fSlope)
  ,fOffset(t.fOffset)
  ,fPt(t.fPt)
  ,fdQdl(t.fdQdl)
{
  //
  // AliTRDmcmTracklet copy constructor
  //

  for (Int_t time = 0; time < kNtimeBins; time++) {
    for (Int_t icl = 0; icl < kNclsPads; icl++) {
      ((AliTRDmcmTracklet &) t).fADC[time][icl] = 0;
    }
    for (Int_t it = 0; it < kNdict; it++) {
      ((AliTRDmcmTracklet &) t).fTrack[time][it] = -1;
    }
    ((AliTRDmcmTracklet &) t).fTime[time]   = 0;
    ((AliTRDmcmTracklet &) t).fCol[time]    = 0;
  }

  if (fGeo) {
    delete fGeo;
  }
  fGeo = new AliTRDgeometry();

}

//_____________________________________________________________________________
AliTRDmcmTracklet::~AliTRDmcmTracklet() 
{
  //
  // AliTRDmcmTracklet destructor
  //

  if (fGPos != 0) delete fGPos;
  if (fGAmp != 0) delete fGAmp;

  if (fGeo) {
    delete fGeo;
  }

}

//_____________________________________________________________________________
AliTRDmcmTracklet &AliTRDmcmTracklet::operator=(const AliTRDmcmTracklet &t)
{
  //
  // Assignment operator
  //

  if (this != &t) ((AliTRDmcmTracklet &) t).Copy(*this); 
  return *this;

}

//_____________________________________________________________________________
void AliTRDmcmTracklet::Copy(TObject &t) const
{
  //
  // Copy function
  //

  ((AliTRDmcmTracklet &) t).fDetector    = fDetector;
  ((AliTRDmcmTracklet &) t).fRow         = fRow;
  ((AliTRDmcmTracklet &) t).fTrackLabel  = fTrackLabel;
  ((AliTRDmcmTracklet &) t).fNclusters   = fNclusters;
  ((AliTRDmcmTracklet &) t).fN           = fN;
  ((AliTRDmcmTracklet &) t).fGPos        = NULL;
  ((AliTRDmcmTracklet &) t).fGAmp        = NULL;
  ((AliTRDmcmTracklet &) t).fTime0       = fTime0;
  ((AliTRDmcmTracklet &) t).fRowz        = fRowz;
  ((AliTRDmcmTracklet &) t).fSlope       = fSlope;
  ((AliTRDmcmTracklet &) t).fOffset      = fOffset;
  ((AliTRDmcmTracklet &) t).fPt          = fPt;
  ((AliTRDmcmTracklet &) t).fdQdl        = fdQdl;

  for (Int_t time = 0; time < kNtimeBins; time++) {
    for (Int_t icl = 0; icl < kNclsPads; icl++) {
      ((AliTRDmcmTracklet &) t).fADC[time][icl] = 0;
    }
    for (Int_t it = 0; it < kNdict; it++) {
      ((AliTRDmcmTracklet &) t).fTrack[time][it] = -1;
    }
    ((AliTRDmcmTracklet &) t).fTime[time]   = 0;
    ((AliTRDmcmTracklet &) t).fCol[time]    = 0;
  }

}

//_____________________________________________________________________________
void AliTRDmcmTracklet::Reset() 
{
  //
  // Reset the tracklet information
  //

  fDetector = -1;
  fRow      = -1;

  for (Int_t time = 0; time < kNtimeBins; time++) {
    for (Int_t icl = 0; icl < kNclsPads; icl++) {
      fADC[time][icl] = 0;
    }
    for (Int_t it = 0; it < kNdict; it++) {
      fTrack[time][it] = -1;
    }
    fTime[time]   = 0;
    fCol[time]    = 0;
  }

  fNclusters  =  0;
  fN          =  0;
  fTrackLabel = -1;

  fGPos->Set(0);
  fGAmp->Set(0);

  fSlope  = 0.0;
  fOffset = 0.0;
  fTime0  = 0.0;
  fRowz   = 0.0;
  fPt     = 0.0;
  fdQdl   = 0.0;

}

//_____________________________________________________________________________
void AliTRDmcmTracklet::AddCluster(Int_t icol, Int_t itb, Float_t *adc, Int_t *track) 
{
  //
  // Add a cluster to the tracklet
  //
 
  if (fNclusters >= kNtimeBins) return;

  for (Int_t icl = 0; icl < kNclsPads; icl++) {
    fADC[fNclusters][icl] = adc[icl]; 
  }

  fTrack[fNclusters][0] = track[0];
  fTrack[fNclusters][1] = track[1];
  fTrack[fNclusters][2] = track[2];
  fTime[fNclusters] = itb;
  fCol[fNclusters] = icol;

  fNclusters++;

}

//_____________________________________________________________________________
void AliTRDmcmTracklet::MakeTrackletGraph(AliTRDgeometry *geo, Float_t field) 
{
  //
  // Tracklet graph of positions (global coordinates, rotated [cm])
  //
  
  if (!geo) {
    AliError("No geometry.");
    return;
  }
  
  if (!AliTRDCommonParam::Instance()) {
    AliError("No common parameters.");
    return;
  }

  AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
  if (!calibration) {
    AliError("No instance of AliTRDcalibDB.");
    return;
  }

  Int_t ilayer, istack;

  ilayer = geo->GetLayer(fDetector);
  istack = geo->GetStack(fDetector);

  AliTRDpadPlane *padPlane = fGeo->GetPadPlane(ilayer,istack);

  Float_t samplFreq = AliTRDCommonParam::Instance()->GetSamplingFrequency();

  Int_t time, col;
  Float_t amp[3];
  Float_t xpos, ypos, xzero, colSize, timeBinSize;
  Float_t vDrift, omegaTau, lorentzAngle, thetaSlope;
  Float_t tiltingAngle;
  Int_t npg = 0;

  for (Int_t icl = 0; icl < fNclusters; icl++) {

    time   = GetClusterTime(icl);

    amp[0] = GetClusterADC(icl)[0];
    amp[1] = GetClusterADC(icl)[1];
    amp[2] = GetClusterADC(icl)[2];

    col    = GetClusterCol(icl);

    if (amp[0] < 0.0 || amp[1] < 0.0 || amp[2] < 0.0) continue;

    ypos = GetClusY(amp,ilayer);

    colSize = padPlane->GetColSize(col);
    vDrift = calibration->GetVdrift(fDetector,col,fRow);
    timeBinSize = vDrift/samplFreq;
    
    // From v4-03-Release to HEAD28Mar06 the sign has changed from "-" to "+" 
    // due to a change in the digitizer
    omegaTau     = TMath::Sign(1.0,(Double_t)field)*GetOmegaTau(vDrift,TMath::Abs(field));
    lorentzAngle = TMath::ATan(omegaTau)*180.0/TMath::Pi();
    
    xpos = (time+0.5) * timeBinSize;
    xpos = geo->GetTime0(ilayer) - xpos;

    ypos = padPlane->GetColPos(col) - (ypos + 0.5) * colSize;

    // ExB correction
    xzero = geo->GetTime0(ilayer);
    ypos = ypos + (xpos-xzero) * omegaTau;

    // tilted pads correction
    thetaSlope = - padPlane->GetRowPos(fRow)/geo->GetTime0(ilayer);
    tiltingAngle = padPlane->GetTiltingAngle()/180.0*TMath::Pi();
    ypos = ypos - (xpos-xzero) * thetaSlope * TMath::Sin(tiltingAngle);

    fGPos->SetPoint(npg,(Double_t)xpos,(Double_t)ypos);
    npg++;

  }

  fGPos->Set(npg);
  
  fTime0 = geo->GetTime0(ilayer) - AliTRDgeometry::CdrHght() - 0.5*AliTRDgeometry::CamHght();
  fRowz = padPlane->GetRowPos(fRow) - padPlane->GetRowSize(fRow)/2.0;

  Double_t xMin = 0;
  Double_t xMax = 0;
  Double_t x, y;
  fGPos->GetPoint(0    ,x,y); 
  xMax = x + 0.1;
  fGPos->GetPoint(npg-1,x,y); 
  xMin = x - 0.1;
  
  TF1 *line = new TF1("line","[0]+x*[1]",xMin,xMax);
  fGPos->Fit(line,"WRQ0");

  fOffset = line->Eval(fTime0);
  fSlope  = TMath::ATan(line->GetParameter(1))*180.0/TMath::Pi();

  line->Delete();
  
  Float_t fx = fTime0;
  Float_t fy = fOffset;
  
  Float_t infSlope = TMath::ATan(fy/fx)/TMath::Pi()*180.0;    
  Float_t alpha    = fSlope - infSlope;
  Float_t r        = TMath::Sqrt(fx*fx + fy*fy)/(2.0*TMath::Sin(alpha/180.0*TMath::Pi()));

  fPt = 0.3 * field * 0.01 * r;
  
  return;

}

//_____________________________________________________________________________
void AliTRDmcmTracklet::MakeClusAmpGraph() 
{
  //
  // Tracklet graph of cluster charges
  //

  Int_t   time;
  Float_t amp[3];
  Int_t   npg = 0;

  fdQdl = 0.0;
  for (Int_t icl = 0; icl < fNclusters; icl++) {

    time   = GetClusterTime(icl);

    amp[0] = GetClusterADC(icl)[0];
    amp[1] = GetClusterADC(icl)[1];
    amp[2] = GetClusterADC(icl)[2];

    fGAmp->SetPoint(npg,(Double_t)(time+0.5),(Double_t)(amp[0]+amp[1]+amp[2]));
    npg++;

    fdQdl += amp[0]+amp[1]+amp[2];

  }

  fGAmp->Set(npg);

  fdQdl /= (Float_t)npg;

  return;

}

//_____________________________________________________________________________
Float_t AliTRDmcmTracklet::GetClusY(Float_t *adc, Int_t layer) const
{
  //
  // Cluster position in the phi direction in pad units (relative to the pad border)
  //

  Float_t ypos = 0.0;

  Float_t a0 = adc[0];
  Float_t a1 = adc[1];
  Float_t a2 = adc[2];

  Float_t t1, t2, w1 ,w2;

  Float_t w = 1.0;  // pad units

  Float_t sigma = 0.0;

  switch(layer) {
  case 0:
    sigma = 0.515; break;
  case 1:
    sigma = 0.501; break;
  case 2:
    sigma = 0.491; break;
  case 3:
    sigma = 0.481; break;
  case 4:
    sigma = 0.471; break;
  case 5:
    sigma = 0.463; break;
  default:
    AliError("Wrong layer number.");
    return 0.0;
  }

  sigma *= w;

  t1 = 0.0;
  w1 = 0.0;
  if( a0 > 0 ) {
    w1 = a0*a0;
    //w1 = a0;
    t1 = w1*((sigma*sigma)/w*TMath::Log(a1/a0)-0.5*w);
  }
  t2 = 0.0;
  w2 = 0.0;
  if( a2 > 0 ) {
    w2 = a2*a2;
    //w2 = a2;
    t2 = w2*((sigma*sigma)/w*TMath::Log(a2/a1)+0.5*w);
  }

  ypos = w*(t1+t2)/(w1+w2);  // range: -0.5*w ... +0.5*w

  return ypos;

}

//_____________________________________________________________________________
void AliTRDmcmTracklet::CookLabel(Float_t frac) 
{
  //
  // Cook the track label from cluster labels
  //

  const Int_t kMaxTracks = 10;
  Int_t trackLabel[kMaxTracks];
  Int_t trackCount[kMaxTracks];
  for (Int_t it = 0; it < kMaxTracks; it++) {
    trackLabel[it] = -1;
    trackCount[it] = 0;
  }

  Bool_t counted;
  Int_t label, nTracks = 0;
  for (Int_t icl = 0; icl < fNclusters; icl++) {

    for (Int_t id = 0; id < kNdict; id++) {

      if (fTrack[icl][id] == -1) continue;

      label = fTrack[icl][id];

      counted = kFALSE;
      for (Int_t it = 0; it < nTracks; it++) {
	if (label == trackLabel[it]) {
	  trackCount[it]++;
	  counted = kTRUE;
	  break;
	}
      }
      if (!counted) {
	trackLabel[nTracks] = label;
	trackCount[nTracks]++;
	nTracks++;
	if (nTracks == kMaxTracks) {
	  AliWarning("Too many tracks for this tracklet.");
	  nTracks--;
	  break;
	}
      }

    }

  }

  for (Int_t it = 0; it < kMaxTracks; it++) {
    if (trackCount[it] >= (Int_t)(frac*fNclusters)) {
      fTrackLabel = trackLabel[it];
      break;
    }
  }

}

//_____________________________________________________________________________
Float_t AliTRDmcmTracklet::GetOmegaTau(Float_t vdrift, Float_t field) const
{
  //
  // Returns omega*tau (tan(Lorentz-angle)) for a given drift velocity <vd> 
  // and a B-field <b> for Xe/CO2 (15%).
  // The values are according to a GARFIELD simulation.
  //

  //
  // Copy of the "AliTRDcalibDB" function, taking as argument the magnetic field too
  //
  
  const Int_t kNb = 5;
  Float_t p0[kNb] = {  0.004810,  0.007412,  0.010252,  0.013409,  0.016888 };
  Float_t p1[kNb] = {  0.054875,  0.081534,  0.107333,  0.131983,  0.155455 };
  Float_t p2[kNb] = { -0.008682, -0.012896, -0.016987, -0.020880, -0.024623 };
  Float_t p3[kNb] = {  0.000155,  0.000238,  0.000330,  0.000428,  0.000541 };

  Int_t ib = ((Int_t) (10 * (field - 0.15)));
  ib       = TMath::Max(  0,ib);
  ib       = TMath::Min(kNb,ib);

  Float_t alphaL = p0[ib] 
                 + p1[ib] * vdrift
                 + p2[ib] * vdrift*vdrift
                 + p3[ib] * vdrift*vdrift*vdrift;

  return TMath::Tan(alphaL);

}
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	// //
	18	// //
	19	// Tracklet object created in the local tracking //
	20	// //
	21	// //
	22	///////////////////////////////////////////////////////
	23
	24	#include <TGraph.h>
	25	#include <TMath.h>
	26	#include <TF1.h>
	27
	28	#include "AliLog.h"
	29
	30	#include "AliTRDcalibDB.h"
	31	#include "AliTRDCommonParam.h"
	32	#include "AliTRDpadPlane.h"
	33	#include "AliTRDgeometry.h"
	34	#include "AliTRDmcmTracklet.h"
	35
	36	ClassImp(AliTRDmcmTracklet)
	37
	38	//_____________________________________________________________________________
	39	AliTRDmcmTracklet::AliTRDmcmTracklet()
	40	:TObject()
	41	,fGeo(0)
	42	,fDetector(-1)
	43	,fRow(-1)
	44	,fTrackLabel(-1)
	45	,fNclusters(0)
	46	,fN(0)
	47	,fGPos(0)
	48	,fGAmp(0)
	49	,fTime0(0)
	50	,fRowz(0)
	51	,fSlope(0)
	52	,fOffset(0)
	53	,fPt(0)
	54	,fdQdl(0)
	55	{
	56	//
	57	// AliTRDmcmTracklet default constructor
	58	//
	59
	60	for (Int_t time = 0; time < kNtimeBins; time++) {
	61	for (Int_t icl = 0; icl < kNclsPads; icl++) {
	62	fADC[time][icl] = 0;
	63	}
	64	for (Int_t it = 0; it < kNdict; it++) {
	65	fTrack[time][it] = -1;
	66	}
	67	fTime[time] = 0;
	68	fCol[time] = 0;
	69	}
	70
	71	fGeo = new AliTRDgeometry();
	72
	73	}
	74
	75	//_____________________________________________________________________________
	76	AliTRDmcmTracklet::AliTRDmcmTracklet(Int_t det, Int_t row, Int_t n)
	77	:TObject()
	78	,fGeo(0)
	79	,fDetector(det)
	80	,fRow(row)
	81	,fTrackLabel(-1)
	82	,fNclusters(0)
	83	,fN(n)
	84	,fGPos(0)
	85	,fGAmp(0)
	86	,fTime0(0)
	87	,fRowz(0)
	88	,fSlope(0)
	89	,fOffset(0)
	90	,fPt(0)
	91	,fdQdl(0)
	92	{
	93	//
	94	// AliTRDmcmTracklet default constructor
	95	//
	96
	97	for (Int_t time = 0; time < kNtimeBins; time++) {
	98	for (Int_t icl = 0; icl < kNclsPads; icl++) {
	99	fADC[time][icl] = 0;
	100	}
	101	for (Int_t it = 0; it < kNdict; it++) {
	102	fTrack[time][it] = -1;
	103	}
	104	fTime[time] = 0;
	105	fCol[time] = 0;
	106	}
	107
	108	fGPos = new TGraph(0);
	109	fGAmp = new TGraph(0);
	110
	111	fGeo = new AliTRDgeometry();
	112
	113	}
	114
	115	//_____________________________________________________________________________
	116	AliTRDmcmTracklet::AliTRDmcmTracklet(const AliTRDmcmTracklet &t)
	117	:TObject(t)
	118	,fGeo(0)
	119	,fDetector(t.fDetector)
	120	,fRow(t.fRow)
	121	,fTrackLabel(t.fTrackLabel)
	122	,fNclusters(t.fNclusters)
	123	,fN(t.fN)
	124	,fGPos(NULL)
	125	,fGAmp(NULL)
	126	,fTime0(t.fTime0)
	127	,fRowz(t.fRowz)
	128	,fSlope(t.fSlope)
	129	,fOffset(t.fOffset)
	130	,fPt(t.fPt)
	131	,fdQdl(t.fdQdl)
	132	{
	133	//
	134	// AliTRDmcmTracklet copy constructor
	135	//
	136
	137	for (Int_t time = 0; time < kNtimeBins; time++) {
	138	for (Int_t icl = 0; icl < kNclsPads; icl++) {
	139	((AliTRDmcmTracklet &) t).fADC[time][icl] = 0;
	140	}
	141	for (Int_t it = 0; it < kNdict; it++) {
	142	((AliTRDmcmTracklet &) t).fTrack[time][it] = -1;
	143	}
	144	((AliTRDmcmTracklet &) t).fTime[time] = 0;
	145	((AliTRDmcmTracklet &) t).fCol[time] = 0;
	146	}
	147
	148	if (fGeo) {
	149	delete fGeo;
	150	}
	151	fGeo = new AliTRDgeometry();
	152
	153	}
	154
	155	//_____________________________________________________________________________
	156	AliTRDmcmTracklet::~AliTRDmcmTracklet()
	157	{
	158	//
	159	// AliTRDmcmTracklet destructor
	160	//
	161
	162	if (fGPos != 0) delete fGPos;
	163	if (fGAmp != 0) delete fGAmp;
	164
	165	if (fGeo) {
	166	delete fGeo;
	167	}
	168
	169	}
	170
	171	//_____________________________________________________________________________
	172	AliTRDmcmTracklet &AliTRDmcmTracklet::operator=(const AliTRDmcmTracklet &t)
	173	{
	174	//
	175	// Assignment operator
	176	//
	177
	178	if (this != &t) ((AliTRDmcmTracklet &) t).Copy(*this);
	179	return *this;
	180
	181	}
	182
	183	//_____________________________________________________________________________
	184	void AliTRDmcmTracklet::Copy(TObject &t) const
	185	{
	186	//
	187	// Copy function
	188	//
	189
	190	((AliTRDmcmTracklet &) t).fDetector = fDetector;
	191	((AliTRDmcmTracklet &) t).fRow = fRow;
	192	((AliTRDmcmTracklet &) t).fTrackLabel = fTrackLabel;
	193	((AliTRDmcmTracklet &) t).fNclusters = fNclusters;
	194	((AliTRDmcmTracklet &) t).fN = fN;
	195	((AliTRDmcmTracklet &) t).fGPos = NULL;
	196	((AliTRDmcmTracklet &) t).fGAmp = NULL;
	197	((AliTRDmcmTracklet &) t).fTime0 = fTime0;
	198	((AliTRDmcmTracklet &) t).fRowz = fRowz;
	199	((AliTRDmcmTracklet &) t).fSlope = fSlope;
	200	((AliTRDmcmTracklet &) t).fOffset = fOffset;
	201	((AliTRDmcmTracklet &) t).fPt = fPt;
	202	((AliTRDmcmTracklet &) t).fdQdl = fdQdl;
	203
	204	for (Int_t time = 0; time < kNtimeBins; time++) {
	205	for (Int_t icl = 0; icl < kNclsPads; icl++) {
	206	((AliTRDmcmTracklet &) t).fADC[time][icl] = 0;
	207	}
	208	for (Int_t it = 0; it < kNdict; it++) {
	209	((AliTRDmcmTracklet &) t).fTrack[time][it] = -1;
	210	}
	211	((AliTRDmcmTracklet &) t).fTime[time] = 0;
	212	((AliTRDmcmTracklet &) t).fCol[time] = 0;
	213	}
	214
	215	}
	216
	217	//_____________________________________________________________________________
	218	void AliTRDmcmTracklet::Reset()
	219	{
	220	//
	221	// Reset the tracklet information
	222	//
	223
	224	fDetector = -1;
	225	fRow = -1;
	226
	227	for (Int_t time = 0; time < kNtimeBins; time++) {
	228	for (Int_t icl = 0; icl < kNclsPads; icl++) {
	229	fADC[time][icl] = 0;
	230	}
	231	for (Int_t it = 0; it < kNdict; it++) {
	232	fTrack[time][it] = -1;
	233	}
	234	fTime[time] = 0;
	235	fCol[time] = 0;
	236	}
	237
	238	fNclusters = 0;
	239	fN = 0;
	240	fTrackLabel = -1;
	241
	242	fGPos->Set(0);
	243	fGAmp->Set(0);
	244
	245	fSlope = 0.0;
	246	fOffset = 0.0;
	247	fTime0 = 0.0;
	248	fRowz = 0.0;
	249	fPt = 0.0;
	250	fdQdl = 0.0;
	251
	252	}
	253
	254	//_____________________________________________________________________________
	255	void AliTRDmcmTracklet::AddCluster(Int_t icol, Int_t itb, Float_t adc, Int_t track)
	256	{
	257	//
	258	// Add a cluster to the tracklet
	259	//
	260
	261	if (fNclusters >= kNtimeBins) return;
	262
	263	for (Int_t icl = 0; icl < kNclsPads; icl++) {
	264	fADC[fNclusters][icl] = adc[icl];
	265	}
	266
	267	fTrack[fNclusters][0] = track[0];
	268	fTrack[fNclusters][1] = track[1];
	269	fTrack[fNclusters][2] = track[2];
	270	fTime[fNclusters] = itb;
	271	fCol[fNclusters] = icol;
	272
	273	fNclusters++;
	274
	275	}
	276
	277	//_____________________________________________________________________________
	278	void AliTRDmcmTracklet::MakeTrackletGraph(AliTRDgeometry *geo, Float_t field)
	279	{
	280	//
	281	// Tracklet graph of positions (global coordinates, rotated [cm])
	282	//
	283
	284	if (!geo) {
	285	AliError("No geometry.");
	286	return;
	287	}
	288
	289	if (!AliTRDCommonParam::Instance()) {
	290	AliError("No common parameters.");
	291	return;
	292	}
	293
	294	AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
	295	if (!calibration) {
	296	AliError("No instance of AliTRDcalibDB.");
	297	return;
	298	}
	299
	300	Int_t ilayer, istack;
	301
	302	ilayer = geo->GetLayer(fDetector);
	303	istack = geo->GetStack(fDetector);
	304
	305	AliTRDpadPlane *padPlane = fGeo->GetPadPlane(ilayer,istack);
	306
	307	Float_t samplFreq = AliTRDCommonParam::Instance()->GetSamplingFrequency();
	308
	309	Int_t time, col;
	310	Float_t amp[3];
	311	Float_t xpos, ypos, xzero, colSize, timeBinSize;
	312	Float_t vDrift, omegaTau, lorentzAngle, thetaSlope;
	313	Float_t tiltingAngle;
	314	Int_t npg = 0;
	315
	316	for (Int_t icl = 0; icl < fNclusters; icl++) {
	317
	318	time = GetClusterTime(icl);
	319
	320	amp[0] = GetClusterADC(icl)[0];
	321	amp[1] = GetClusterADC(icl)[1];
	322	amp[2] = GetClusterADC(icl)[2];
	323
	324	col = GetClusterCol(icl);
	325
	326	if (amp[0] < 0.0 \|\| amp[1] < 0.0 \|\| amp[2] < 0.0) continue;
	327
	328	ypos = GetClusY(amp,ilayer);
	329
	330	colSize = padPlane->GetColSize(col);
	331	vDrift = calibration->GetVdrift(fDetector,col,fRow);
	332	timeBinSize = vDrift/samplFreq;
	333
	334	// From v4-03-Release to HEAD28Mar06 the sign has changed from "-" to "+"
	335	// due to a change in the digitizer
	336	omegaTau = TMath::Sign(1.0,(Double_t)field)*GetOmegaTau(vDrift,TMath::Abs(field));
	337	lorentzAngle = TMath::ATan(omegaTau)*180.0/TMath::Pi();
	338
	339	xpos = (time+0.5) * timeBinSize;
	340	xpos = geo->GetTime0(ilayer) - xpos;
	341
	342	ypos = padPlane->GetColPos(col) - (ypos + 0.5) * colSize;
	343
	344	// ExB correction
	345	xzero = geo->GetTime0(ilayer);
	346	ypos = ypos + (xpos-xzero) * omegaTau;
	347
	348	// tilted pads correction
	349	thetaSlope = - padPlane->GetRowPos(fRow)/geo->GetTime0(ilayer);
	350	tiltingAngle = padPlane->GetTiltingAngle()/180.0*TMath::Pi();
	351	ypos = ypos - (xpos-xzero) * thetaSlope * TMath::Sin(tiltingAngle);
	352
	353	fGPos->SetPoint(npg,(Double_t)xpos,(Double_t)ypos);
	354	npg++;
	355
	356	}
	357
	358	fGPos->Set(npg);
	359
	360	fTime0 = geo->GetTime0(ilayer) - AliTRDgeometry::CdrHght() - 0.5*AliTRDgeometry::CamHght();
	361	fRowz = padPlane->GetRowPos(fRow) - padPlane->GetRowSize(fRow)/2.0;
	362
	363	Double_t xMin = 0;
	364	Double_t xMax = 0;
	365	Double_t x, y;
	366	fGPos->GetPoint(0 ,x,y);
	367	xMax = x + 0.1;
	368	fGPos->GetPoint(npg-1,x,y);
	369	xMin = x - 0.1;
	370
	371	TF1 line = new TF1("line","[0]+x[1]",xMin,xMax);
	372	fGPos->Fit(line,"WRQ0");
	373
	374	fOffset = line->Eval(fTime0);
	375	fSlope = TMath::ATan(line->GetParameter(1))*180.0/TMath::Pi();
	376
	377	line->Delete();
	378
	379	Float_t fx = fTime0;
	380	Float_t fy = fOffset;
	381
	382	Float_t infSlope = TMath::ATan(fy/fx)/TMath::Pi()*180.0;
	383	Float_t alpha = fSlope - infSlope;
	384	Float_t r = TMath::Sqrt(fxfx + fyfy)/(2.0TMath::Sin(alpha/180.0TMath::Pi()));
	385
	386	fPt = 0.3 * field * 0.01 * r;
	387
	388	return;
	389
	390	}
	391
	392	//_____________________________________________________________________________
	393	void AliTRDmcmTracklet::MakeClusAmpGraph()
	394	{
	395	//
	396	// Tracklet graph of cluster charges
	397	//
	398
	399	Int_t time;
	400	Float_t amp[3];
	401	Int_t npg = 0;
	402
	403	fdQdl = 0.0;
	404	for (Int_t icl = 0; icl < fNclusters; icl++) {
	405
	406	time = GetClusterTime(icl);
	407
	408	amp[0] = GetClusterADC(icl)[0];
	409	amp[1] = GetClusterADC(icl)[1];
	410	amp[2] = GetClusterADC(icl)[2];
	411
	412	fGAmp->SetPoint(npg,(Double_t)(time+0.5),(Double_t)(amp[0]+amp[1]+amp[2]));
	413	npg++;
	414
	415	fdQdl += amp[0]+amp[1]+amp[2];
	416
	417	}
	418
	419	fGAmp->Set(npg);
	420
	421	fdQdl /= (Float_t)npg;
	422
	423	return;
	424
	425	}
	426
	427	//_____________________________________________________________________________
	428	Float_t AliTRDmcmTracklet::GetClusY(Float_t *adc, Int_t layer) const
	429	{
	430	//
	431	// Cluster position in the phi direction in pad units (relative to the pad border)
	432	//
	433
	434	Float_t ypos = 0.0;
	435
	436	Float_t a0 = adc[0];
	437	Float_t a1 = adc[1];
	438	Float_t a2 = adc[2];
	439
	440	Float_t t1, t2, w1 ,w2;
	441
	442	Float_t w = 1.0; // pad units
	443
	444	Float_t sigma = 0.0;
	445
	446	switch(layer) {
	447	case 0:
	448	sigma = 0.515; break;
	449	case 1:
	450	sigma = 0.501; break;
	451	case 2:
	452	sigma = 0.491; break;
	453	case 3:
	454	sigma = 0.481; break;
	455	case 4:
	456	sigma = 0.471; break;
	457	case 5:
	458	sigma = 0.463; break;
	459	default:
	460	AliError("Wrong layer number.");
	461	return 0.0;
	462	}
	463
	464	sigma *= w;
	465
	466	t1 = 0.0;
	467	w1 = 0.0;
	468	if( a0 > 0 ) {
	469	w1 = a0*a0;
	470	//w1 = a0;
	471	t1 = w1((sigmasigma)/wTMath::Log(a1/a0)-0.5w);
	472	}
	473	t2 = 0.0;
	474	w2 = 0.0;
	475	if( a2 > 0 ) {
	476	w2 = a2*a2;
	477	//w2 = a2;
	478	t2 = w2((sigmasigma)/wTMath::Log(a2/a1)+0.5w);
	479	}
	480
	481	ypos = w(t1+t2)/(w1+w2); // range: -0.5w ... +0.5*w
	482
	483	return ypos;
	484
	485	}
	486
	487	//_____________________________________________________________________________
	488	void AliTRDmcmTracklet::CookLabel(Float_t frac)
	489	{
	490	//
	491	// Cook the track label from cluster labels
	492	//
	493
	494	const Int_t kMaxTracks = 10;
	495	Int_t trackLabel[kMaxTracks];
	496	Int_t trackCount[kMaxTracks];
	497	for (Int_t it = 0; it < kMaxTracks; it++) {
	498	trackLabel[it] = -1;
	499	trackCount[it] = 0;
	500	}
	501
	502	Bool_t counted;
	503	Int_t label, nTracks = 0;
	504	for (Int_t icl = 0; icl < fNclusters; icl++) {
	505
	506	for (Int_t id = 0; id < kNdict; id++) {
	507
	508	if (fTrack[icl][id] == -1) continue;
	509
	510	label = fTrack[icl][id];
	511
	512	counted = kFALSE;
	513	for (Int_t it = 0; it < nTracks; it++) {
	514	if (label == trackLabel[it]) {
	515	trackCount[it]++;
	516	counted = kTRUE;
	517	break;
	518	}
	519	}
	520	if (!counted) {
	521	trackLabel[nTracks] = label;
	522	trackCount[nTracks]++;
	523	nTracks++;
	524	if (nTracks == kMaxTracks) {
	525	AliWarning("Too many tracks for this tracklet.");
	526	nTracks--;
	527	break;
	528	}
	529	}
	530
	531	}
	532
	533	}
	534
	535	for (Int_t it = 0; it < kMaxTracks; it++) {
	536	if (trackCount[it] >= (Int_t)(frac*fNclusters)) {
	537	fTrackLabel = trackLabel[it];
	538	break;
	539	}
	540	}
	541
	542	}
	543
	544	//_____________________________________________________________________________
	545	Float_t AliTRDmcmTracklet::GetOmegaTau(Float_t vdrift, Float_t field) const
	546	{
	547	//
	548	// Returns omega*tau (tan(Lorentz-angle)) for a given drift velocity <vd>
	549	// and a B-field <b> for Xe/CO2 (15%).
	550	// The values are according to a GARFIELD simulation.
	551	//
	552
	553	//
	554	// Copy of the "AliTRDcalibDB" function, taking as argument the magnetic field too
	555	//
	556
	557	const Int_t kNb = 5;
	558	Float_t p0[kNb] = { 0.004810, 0.007412, 0.010252, 0.013409, 0.016888 };
	559	Float_t p1[kNb] = { 0.054875, 0.081534, 0.107333, 0.131983, 0.155455 };
	560	Float_t p2[kNb] = { -0.008682, -0.012896, -0.016987, -0.020880, -0.024623 };
	561	Float_t p3[kNb] = { 0.000155, 0.000238, 0.000330, 0.000428, 0.000541 };
	562
	563	Int_t ib = ((Int_t) (10 * (field - 0.15)));
	564	ib = TMath::Max( 0,ib);
	565	ib = TMath::Min(kNb,ib);
	566
	567	Float_t alphaL = p0[ib]
	568	+ p1[ib] * vdrift
	569	+ p2[ib] * vdrift*vdrift
	570	+ p3[ib] * vdriftvdriftvdrift;
	571
	572	return TMath::Tan(alphaL);
	573
	574	}