[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
0ee00e25	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
e3b2b5e5	16	///////////////////////////////////////////////////////
	17	// //
	18	// //
	19	// Tracklet object created in the local tracking //
	20	// //
	21	// //
	22	///////////////////////////////////////////////////////
	23
0ee00e25	24	#include <TGraph.h>
	25	#include <TMath.h>
	26	#include <TF1.h>
	27
6d50f529	28	#include "AliLog.h"
6d50f529	29
0ee00e25	30	#include "AliTRDcalibDB.h"
	31	#include "AliTRDCommonParam.h"
	32	#include "AliTRDpadPlane.h"
	33	#include "AliTRDgeometry.h"
0ee00e25	34	#include "AliTRDmcmTracklet.h"
	35
	36	ClassImp(AliTRDmcmTracklet)
	37
	38	//_____________________________________________________________________________
	39	AliTRDmcmTracklet::AliTRDmcmTracklet()
6d50f529	40	:TObject()
f162af62	41	,fGeo(0)
6d50f529	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)
0ee00e25	55	{
0ee00e25	56	//
	57	// AliTRDmcmTracklet default constructor
	58	//
	59
0ee00e25	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
f162af62	71	fGeo = new AliTRDgeometry();
f162af62	72
0ee00e25	73	}
	74
	75	//_____________________________________________________________________________
	76	AliTRDmcmTracklet::AliTRDmcmTracklet(Int_t det, Int_t row, Int_t n)
6d50f529	77	:TObject()
f162af62	78	,fGeo(0)
6d50f529	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)
0ee00e25	92	{
0ee00e25	93	//
	94	// AliTRDmcmTracklet default constructor
	95	//
	96
0ee00e25	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
0ee00e25	108	fGPos = new TGraph(0);
	109	fGAmp = new TGraph(0);
	110
f162af62	111	fGeo = new AliTRDgeometry();
f162af62	112
6d50f529	113	}
	114
	115	//_____________________________________________________________________________
	116	AliTRDmcmTracklet::AliTRDmcmTracklet(const AliTRDmcmTracklet &t)
	117	:TObject(t)
f162af62	118	,fGeo(0)
6d50f529	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	}
0ee00e25	147
f162af62	148	if (fGeo) {
	149	delete fGeo;
	150	}
	151	fGeo = new AliTRDgeometry();
	152
0ee00e25	153	}
	154
	155	//_____________________________________________________________________________
	156	AliTRDmcmTracklet::~AliTRDmcmTracklet()
	157	{
0ee00e25	158	//
	159	// AliTRDmcmTracklet destructor
	160	//
3fcb9908	161
	162	if (fGPos != 0) delete fGPos;
	163	if (fGAmp != 0) delete fGAmp;
	164
f162af62	165	if (fGeo) {
	166	delete fGeo;
	167	}
	168
e3b2b5e5	169	}
	170
	171	//_____________________________________________________________________________
	172	AliTRDmcmTracklet &AliTRDmcmTracklet::operator=(const AliTRDmcmTracklet &t)
	173	{
	174	//
6d50f529	175	// Assignment operator
e3b2b5e5	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	//
6d50f529	187	// Copy function
e3b2b5e5	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	//_____________________________________________________________________________
3fcb9908	218	void AliTRDmcmTracklet::Reset()
e3b2b5e5	219	{
e3b2b5e5	220	//
3fcb9908	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	{
e3b2b5e5	257	//
3fcb9908	258	// Add a cluster to the tracklet
e3b2b5e5	259	//
	260
	261	if (fNclusters >= kNtimeBins) return;
	262
	263	for (Int_t icl = 0; icl < kNclsPads; icl++) {
e3b2b5e5	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++;
0ee00e25	274
	275	}
	276
	277	//_____________________________________________________________________________
	278	void AliTRDmcmTracklet::MakeTrackletGraph(AliTRDgeometry *geo, Float_t field)
	279	{
0ee00e25	280	//
	281	// Tracklet graph of positions (global coordinates, rotated [cm])
	282	//
	283
	284	if (!geo) {
6d50f529	285	AliError("No geometry.");
0ee00e25	286	return;
	287	}
	288
f162af62	289	if (!AliTRDCommonParam::Instance()) {
6d50f529	290	AliError("No common parameters.");
0ee00e25	291	return;
	292	}
	293
	294	AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
6d50f529	295	if (!calibration) {
6d50f529	296	AliError("No instance of AliTRDcalibDB.");
0ee00e25	297	return;
	298	}
	299
053767a4	300	Int_t ilayer, istack;
0ee00e25	301
053767a4	302	ilayer = geo->GetLayer(fDetector);
053767a4	303	istack = geo->GetStack(fDetector);
0ee00e25	304
053767a4	305	AliTRDpadPlane *padPlane = fGeo->GetPadPlane(ilayer,istack);
0ee00e25	306
f162af62	307	Float_t samplFreq = AliTRDCommonParam::Instance()->GetSamplingFrequency();
0ee00e25	308
	309	Int_t time, col;
	310	Float_t amp[3];
e3b2b5e5	311	Float_t xpos, ypos, xzero, colSize, timeBinSize;
	312	Float_t vDrift, omegaTau, lorentzAngle, thetaSlope;
	313	Float_t tiltingAngle;
0ee00e25	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
053767a4	328	ypos = GetClusY(amp,ilayer);
0ee00e25	329
e3b2b5e5	330	colSize = padPlane->GetColSize(col);
0ee00e25	331	vDrift = calibration->GetVdrift(fDetector,col,fRow);
e3b2b5e5	332	timeBinSize = vDrift/samplFreq;
0ee00e25	333
	334	// From v4-03-Release to HEAD28Mar06 the sign has changed from "-" to "+"
	335	// due to a change in the digitizer
6d50f529	336	omegaTau = TMath::Sign(1.0,(Double_t)field)*GetOmegaTau(vDrift,TMath::Abs(field));
e3b2b5e5	337	lorentzAngle = TMath::ATan(omegaTau)*180.0/TMath::Pi();
0ee00e25	338
e3b2b5e5	339	xpos = (time+0.5) * timeBinSize;
053767a4	340	xpos = geo->GetTime0(ilayer) - xpos;
0ee00e25	341
e3b2b5e5	342	ypos = padPlane->GetColPos(col) - (ypos + 0.5) * colSize;
0ee00e25	343
0ee00e25	344	// ExB correction
053767a4	345	xzero = geo->GetTime0(ilayer);
e3b2b5e5	346	ypos = ypos + (xpos-xzero) * omegaTau;
0ee00e25	347
0ee00e25	348	// tilted pads correction
053767a4	349	thetaSlope = - padPlane->GetRowPos(fRow)/geo->GetTime0(ilayer);
e3b2b5e5	350	tiltingAngle = padPlane->GetTiltingAngle()/180.0*TMath::Pi();
e3b2b5e5	351	ypos = ypos - (xpos-xzero) * thetaSlope * TMath::Sin(tiltingAngle);
0ee00e25	352
e3b2b5e5	353	fGPos->SetPoint(npg,(Double_t)xpos,(Double_t)ypos);
0ee00e25	354	npg++;
	355
	356	}
	357
	358	fGPos->Set(npg);
	359
053767a4	360	fTime0 = geo->GetTime0(ilayer) - AliTRDgeometry::CdrHght() - 0.5*AliTRDgeometry::CamHght();
9e79a757	361	fRowz = padPlane->GetRowPos(fRow) - padPlane->GetRowSize(fRow)/2.0;
0ee00e25	362
6d50f529	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;
0ee00e25	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
e3b2b5e5	379	Float_t fx = fTime0;
e3b2b5e5	380	Float_t fy = fOffset;
0ee00e25	381
e3b2b5e5	382	Float_t infSlope = TMath::ATan(fy/fx)/TMath::Pi()*180.0;
6d50f529	383	Float_t alpha = fSlope - infSlope;
6d50f529	384	Float_t r = TMath::Sqrt(fxfx + fyfy)/(2.0TMath::Sin(alpha/180.0TMath::Pi()));
0ee00e25	385
e3b2b5e5	386	fPt = 0.3 * field * 0.01 * r;
0ee00e25	387
	388	return;
	389
	390	}
	391
	392	//_____________________________________________________________________________
	393	void AliTRDmcmTracklet::MakeClusAmpGraph()
	394	{
	395	//
	396	// Tracklet graph of cluster charges
	397	//
	398
6d50f529	399	Int_t time;
0ee00e25	400	Float_t amp[3];
6d50f529	401	Int_t npg = 0;
6d50f529	402
0ee00e25	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	//_____________________________________________________________________________
053767a4	428	Float_t AliTRDmcmTracklet::GetClusY(Float_t *adc, Int_t layer) const
0ee00e25	429	{
	430	//
	431	// Cluster position in the phi direction in pad units (relative to the pad border)
	432	//
	433
e3b2b5e5	434	Float_t ypos = 0.0;
0ee00e25	435
e3b2b5e5	436	Float_t a0 = adc[0];
	437	Float_t a1 = adc[1];
	438	Float_t a2 = adc[2];
0ee00e25	439
e3b2b5e5	440	Float_t t1, t2, w1 ,w2;
0ee00e25	441
e3b2b5e5	442	Float_t w = 1.0; // pad units
0ee00e25	443
e3b2b5e5	444	Float_t sigma = 0.0;
0ee00e25	445
053767a4	446	switch(layer) {
0ee00e25	447	case 0:
e3b2b5e5	448	sigma = 0.515; break;
0ee00e25	449	case 1:
e3b2b5e5	450	sigma = 0.501; break;
0ee00e25	451	case 2:
e3b2b5e5	452	sigma = 0.491; break;
0ee00e25	453	case 3:
e3b2b5e5	454	sigma = 0.481; break;
0ee00e25	455	case 4:
e3b2b5e5	456	sigma = 0.471; break;
0ee00e25	457	case 5:
e3b2b5e5	458	sigma = 0.463; break;
0ee00e25	459	default:
053767a4	460	AliError("Wrong layer number.");
0ee00e25	461	return 0.0;
	462	}
	463
e3b2b5e5	464	sigma *= w;
0ee00e25	465
e3b2b5e5	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);
0ee00e25	472	}
e3b2b5e5	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);
0ee00e25	479	}
0ee00e25	480
e3b2b5e5	481	ypos = w(t1+t2)/(w1+w2); // range: -0.5w ... +0.5*w
0ee00e25	482
e3b2b5e5	483	return ypos;
0ee00e25	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) {
6d50f529	525	AliWarning("Too many tracks for this tracklet.");
0ee00e25	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	//_____________________________________________________________________________
e3b2b5e5	545	Float_t AliTRDmcmTracklet::GetOmegaTau(Float_t vdrift, Float_t field) const
0ee00e25	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]
6d50f529	568	+ p1[ib] * vdrift
	569	+ p2[ib] * vdrift*vdrift
	570	+ p3[ib] * vdriftvdriftvdrift;
0ee00e25	571
	572	return TMath::Tan(alphaL);
	573
	574	}