[u/mrichter/AliRoot.git] / HLT / comp / AliL3Modeller.cxx

//$Id$

// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
//*-- Copyright &copy ASV

#include "AliL3StandardIncludes.h"

#include "AliL3Modeller.h"
#include "AliL3MemHandler.h"
#ifdef use_aliroot
#include "AliL3FileHandler.h"
#endif
#include "AliL3TrackArray.h"
#include "AliL3ModelTrack.h"
#include "AliL3DigitData.h"
#include "AliL3Transform.h"

#if GCCVERSION == 3
using namespace std;
#endif

//_____________________________________________________________
// AliL3Modeller
//
// Class for modeling TPC data.
// 
// This performs the cluster finding, based on track parameters.
// Basically it propagates the tracks to all padrows, and looks 
// for a corresponding cluster. For the moment only cog is calculated,
// and no deconvolution is done. 

ClassImp(AliL3Modeller)

AliL3Modeller::AliL3Modeller()
{
  fMemHandler=0;
  fTracks=0;
  fTrackThreshold=0;
  SetOverlap();
  SetTrackThreshold();
  SetSearchRange();
}


AliL3Modeller::~AliL3Modeller()
{
  if(fMemHandler)
    delete fMemHandler;
  if(fTracks)
    delete fTracks;
}

void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,Bool_t houghtracks,Bool_t binary)
{
  fSlice = slice;
  fPatch = patch;

  sprintf(fPath,"%s",path);
  
  fTracks = new AliL3TrackArray("AliL3ModelTrack");
  
  Char_t fname[100];
  AliL3MemHandler *file = new AliL3MemHandler();
  //sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
  sprintf(fname,"%s/tracks_ho_%d.raw",trackdata,fSlice);
  if(!file->SetBinaryInput(fname))
    {
      cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
      return;
    }
  file->Binary2TrackArray(fTracks);
  file->CloseBinaryInput();
  delete file;
  

  if(!houghtracks)
    fTracks->QSort();
  
  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
    {
      AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
      if(!track) continue;
      track->Init(fSlice,fPatch);

      //Only if the tracks has been merged across sector boundaries:
      //if(!houghtracks)
      //track->Rotate(fSlice,kTRUE); //!!!!!!!!!!!!!!!!!!!
      
      track->CalculateHelix();
    }    
  
  CalculateCrossingPoints();
  
  CheckForOverlaps();
  
  UInt_t ndigits=0;
  AliL3DigitRowData *digits=0;
#ifdef use_aliroot
  fMemHandler = new AliL3FileHandler();
  fMemHandler->Init(slice,patch);
  if(binary == kFALSE)
    {
      sprintf(fname,"%s/digitfile.root",fPath);
      fMemHandler->SetAliInput(fname);
      digits = fMemHandler->AliDigits2Memory(ndigits);
    }
  else
    {
      sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
      if(!fMemHandler->SetBinaryInput(fname))
	{
	  cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
	  return;
	}
      digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
    }
#else
  fMemHandler = new AliL3MemHandler();
  fMemHandler->Init(slice,patch);
  if(binary == kFALSE)
    {
      cerr<<"AliL3Modeller::Init : Compile with AliROOT if you want rootfile as input"<<endl;
      return;
    }
  else
    {
      sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
      if(!fMemHandler->SetBinaryInput(fname))
	{
	  cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
	  return;
	}
    }
  digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
#endif
  
  SetInputData(digits);
}

void AliL3Modeller::FindClusters()
{
  cout<<"AliL3Modeller::FindClusters : Processing slice "<<fSlice<<" patch "<<fPatch<<endl;
  if(!fTracks)
    {
      cerr<<"AliL3Modeller::Process : No tracks"<<endl;
      return;
    }
  if(!fRowData)
    {
      cerr<<"AliL3Modeller::Process : No data "<<endl;
      return;
    }
  
  AliL3DigitRowData *rowPt = fRowData;
  AliL3DigitData *digPt=0;

  Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
  Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
  Int_t bounds = ntimes*npads;
  Digit *row = new Digit[bounds];
  
  Int_t seq_charge;
  Int_t pad,time,index;
  Short_t charge;
  Cluster cluster;

  for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
    {
      fCurrentPadRow = i;
      memset((void*)row,0,ntimes*npads*sizeof(Digit));
      digPt = (AliL3DigitData*)rowPt->fDigitData;
      //cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
      for(UInt_t j=0; j<rowPt->fNDigit; j++)
	{
	  pad = digPt[j].fPad;
	  time = digPt[j].fTime;
	  charge = digPt[j].fCharge;
	  row[ntimes*pad+time].fCharge = charge;
	  row[ntimes*pad+time].fUsed = kFALSE;
	  //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
	}
      
      for(Int_t k=0; k<fTracks->GetNTracks(); k++)
	{
	  AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
	  if(!track) continue;
	  
	  if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetOverlap(i)>=0)
	    {
	      track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch.
	      continue;
	    }
	  
	  Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
	  Int_t hittime = (Int_t)rint(track->GetTimeHit(i));

	  if(!CheckCluster(row,hitpad,hittime)) //Not a good cluster.
	    {
	      track->SetCluster(i,0,0,0,0,0,0); 
	      continue;
	    }

	  //cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
	  pad = hitpad;
	  time = hittime;

	  Int_t padsign=-1;
	  Int_t timesign=-1;
	  
	  memset(&cluster,0,sizeof(Cluster));
	  //cout<<"Processing padrow "<<i<<" with hittime "<<hittime<<" hitpad "<<hitpad<<" charge "<<row[ntimes*pad + time].fCharge<<" index "<<ntimes*pad + time<<endl;
	  Int_t npads=0;
	  Int_t last_mean = hittime;
	  while(1)//Process this padrow
	    {
	      if(pad < 0 || pad >= AliL3Transform::GetNPads(i)) 
		{
		  //cout<<"Pad = "<<pad<<" on row "<<i<<endl;
		  FillCluster(track,&cluster,i,npads);
		  break;
		}
	      seq_charge=0;
	      timesign=-1;
	      time = hittime;
	      
	      while(1) //Process sequence on this pad:
		{
		  if(time < 0 || time >= AliL3Transform::GetNTimeBins()) 
		    break;
		
		  index = ntimes*pad + time;
		  if(index < 0 || index >= bounds)
		    {
		      cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
			<<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
		      break;
		    }
		  
		  charge = row[index].fCharge;
		  
		  if(charge==0) //Empty timebin, so try and expand the search within the limits on this pad
		    {
		      if(seq_charge==0 && abs(time-hittime) <= fTimeSearch)//No sequence yet, keep looking
			{
			  if(timesign==-1)
			    time--;
			  else
			    time++;
			  continue;
			}
		      else if(timesign==-1) //Switch search direction
			{
			  time = hittime + 1;
			  timesign = 1;
			  continue;
			}
		      else //ok, boundary reached, this pad is done.
			break;
		    }
		  
		  if(row[ntimes*pad+time].fUsed==kTRUE) //Don't use digits several times. This leads to mult. rec.tracks.
		    {
		      time += timesign;
		      continue;
		    }
		  
		  seq_charge += charge;
		  
		  //Update the cluster parameters with this timebin
		  cluster.fTime += time*charge;
		  cluster.fPad += pad*charge;
		  cluster.fCharge += charge;
		  cluster.fSigmaY2 += pad*pad*charge;
		  cluster.fSigmaZ2 += time*time*charge;
		  
		  row[ntimes*pad+time].fUsed = kTRUE;
		  time += timesign;
		}
	      
	      //Always compare with the current mean in time of the cluster under construction.
	      if(cluster.fCharge)
		last_mean = (Int_t)rint((Float_t)(cluster.fTime/cluster.fCharge));
	      
	      if(seq_charge)//There was something on this pad, so keep looking on the neighbouring pad
		{
		  pad += padsign;
		  npads++;
		}
	      else
		{
		  if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadSearch && pad > 0)
		    {
		      if(padsign==-1)
			pad--;
		      else
			pad++;
		      continue;
		    }
		  else if(padsign==-1)
		    {
		      pad = hitpad + 1;
		      padsign = 1;
		      continue;
		    }
		  else
		    {
		      FillCluster(track,&cluster,i,npads);
		      break;
		    }
		}
	    }
	  //cout<<"done"<<endl;
	}
      FillZeros(rowPt,row);
      fMemHandler->UpdateRowPointer(rowPt);
    }
  delete [] row;
  //cout<<"done processing"<<endl;
  
  
  //Debug:
  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
    {
      AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
      if(!track) continue;
      if(track->GetNClusters() != AliL3Transform::GetNRows(fPatch))
	cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<AliL3Transform::GetNRows(fPatch)<<endl<<endl;
    }
  
}

Bool_t AliL3Modeller::CheckCluster(Digit *row,Int_t hitpad,Int_t hittime)
{
  Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
  Bool_t seq_was_falling=kFALSE;
  Int_t last_seq_charge=0;
  Int_t cluster_charge=0;
  Int_t time,pad,charge,padsign=-1,timesign=-1;
  time = hittime;
  pad = hitpad;
  //cout<<"Checking cluster "<<hitpad<<" hittime "<<hittime<<endl;
  while(1)
    {
      
      Bool_t last_was_falling=kFALSE;
      Int_t last_charge=0;
      Int_t seq_charge=0;
      while(1)
	{
	  
	  charge = row[(ntimes*pad+time)].fCharge;
	  //cout<<"Charge "<<charge<<" on pad "<<pad<<" time "<<time<<endl;
	  if(charge==0)
	    {
	      if(seq_charge==0 && abs(time-hittime) <= fTimeSearch)
		{
		  if(timesign==-1)
		    time--;
		  else
		    time++;
		  continue;
		}
	      else if(timesign==-1)
		{
		  last_charge = row[(ntimes*hitpad+hittime)].fCharge;
		  time = hittime+1;
		  timesign=1;
		  if(last_charge < row[(ntimes*hitpad+(hittime-1))].fCharge)
		    last_was_falling=kTRUE;
		  else
		    last_was_falling=kFALSE;
		  continue;
		}
	      else
		break;
	    }

	  if(charge > last_charge)
	    {
	      if(last_was_falling) //This is a overlapping cluster.
		return kFALSE;
	    }
	  else
	    last_was_falling = kTRUE;

	  cluster_charge+=charge;
	  seq_charge += charge;
	  last_charge = charge;
	  time += timesign;
	}
      
      if(seq_charge)
	pad += padsign;
      else
	{
	  if(cluster_charge==0 && abs(pad-hitpad)<=fPadSearch && pad > 0)
	    {
	      if(padsign==-1)
		pad--;
	      else
		pad++;
	      continue;
	    }
	  else if(padsign==-1)
	    {
	      pad = hitpad+1;
	      padsign = 1;
	      seq_was_falling=kFALSE;
	      continue;
	    }
	  else
	    break;
	}
      
      if(seq_charge > last_seq_charge)
	{
	  if(seq_was_falling)
	    return kFALSE;
	  else
	    seq_was_falling=kTRUE;
	}
      last_seq_charge = seq_charge;
    }
  
  return kTRUE;
}

void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row,Int_t npads)
{
  if(cluster->fCharge==0)
    {
      track->SetCluster(row,0,0,0,0,0,0);
      return;
    }
  Float_t fcharge = (Float_t)cluster->fCharge;
  Float_t fpad = ((Float_t)cluster->fPad/fcharge);
  Float_t ftime = ((Float_t)cluster->fTime/fcharge);
  Float_t sigmaY2,sigmaZ2;
  CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
  track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2,npads);
}

void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Digit *row)
{
  //Fill zero where data has been used.
  
  Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
  AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
  for(UInt_t j=0; j<rowPt->fNDigit; j++)
    {
      Int_t pad = digPt[j].fPad;
      Int_t time = digPt[j].fTime;
      if(row[ntimes*pad+time].fUsed==kTRUE)
	digPt[j].fCharge = 0;
    }
}

void AliL3Modeller::WriteRemaining()
{
  //Write remaining (nonzero) digits to file.
  
  AliL3DigitRowData *rowPt;
  rowPt = (AliL3DigitRowData*)fRowData;
  Int_t digitcount=0;
  Int_t ndigits[(AliL3Transform::GetNRows(fPatch))];
  for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
    {
      AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
      ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]=0;
      for(UInt_t j=0; j<rowPt->fNDigit; j++)
	{
	  if(digPt[j].fCharge==0) continue;
	  digitcount++;
	  ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]++;
	}
      //cout<<"Difference "<<(int)ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]<<" "<<(int)rowPt->fNDigit<<endl;
      fMemHandler->UpdateRowPointer(rowPt);
    }
  
  Int_t size = digitcount*sizeof(AliL3DigitData) + AliL3Transform::GetNRows(fPatch)*sizeof(AliL3DigitRowData);
  Byte_t *data = new Byte_t[size];
  memset(data,0,size);
  AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
  rowPt = (AliL3DigitRowData*)fRowData;
  
  for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
    {
      Int_t localcount=0;
      tempPt->fRow = i;
      tempPt->fNDigit = ndigits[(i-AliL3Transform::GetFirstRow(fPatch))];
      AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
      for(UInt_t j=0; j<rowPt->fNDigit; j++)
	{
	  if(digPt[j].fCharge==0) continue;
	  if(localcount >= ndigits[(i-AliL3Transform::GetFirstRow(fPatch))])
	    {
	      cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
	      return;
	    }
	  tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
	  tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
	  tempPt->fDigitData[localcount].fTime = digPt[j].fTime;

	  localcount++;
	}
      if(ndigits[(i-AliL3Transform::GetFirstRow(fPatch))] != localcount)
	{
	  cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
	  return;
	}
      fMemHandler->UpdateRowPointer(rowPt);
      Byte_t *tmp = (Byte_t*)tempPt;
      Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]*sizeof(AliL3DigitData);
      tmp += size;
      tempPt = (AliL3DigitRowData*)tmp;
    }

  Char_t fname[100];
  AliL3MemHandler *mem = new AliL3MemHandler();
  sprintf(fname,"%s/comp/remains_%d_%d.raw",fPath,fSlice,fPatch);
  mem->SetBinaryOutput(fname);
  mem->Memory2CompBinary((UInt_t)AliL3Transform::GetNRows(fPatch),(AliL3DigitRowData*)data);
  mem->CloseBinaryOutput();
  delete mem;
  delete [] data;
}


void AliL3Modeller::CalculateCrossingPoints()
{
  //cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
  if(!fTracks)
    {
      cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
      return;
    }
  Float_t hit[3];
  
  Int_t sector,row;
  for(Int_t i=AliL3Transform::GetLastRow(fPatch); i>=AliL3Transform::GetFirstRow(fPatch); i--)
    {
      for(Int_t j=0; j<fTracks->GetNTracks(); j++)
	{
	  AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
	  if(!track) continue;

	  if(!track->GetCrossingPoint(i,hit)) 
	    {
	      //cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
	      //" pt "<<track->GetPt()<<
	      //" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
		//"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
		//"--------"<<endl;
	      fTracks->Remove(j);
	      continue;
	    }
	  //cout<<"X "<<hit[0]<<" Y "<<hit[1]<<" Z "<<hit[2]<<" tgl "<<track->GetTgl()<<endl;
	  
	  AliL3Transform::Slice2Sector(fSlice,i,sector,row);
	  AliL3Transform::Local2Raw(hit,sector,row);
	  //cout<<"Pad "<<hit[1]<<" time "<<hit[2]<<" in sector "<<sector<<" row "<<row<<endl;
	  if(hit[1]<0 || hit[1]>AliL3Transform::GetNPads(i) ||
	     hit[2]<0 || hit[2]>AliL3Transform::GetNTimeBins())
	    {//Track is leaving the patch, so flag the track hits (<0)
	      track->SetPadHit(i,-1);
	      track->SetTimeHit(i,-1);
	      continue;
	    }
	  
	  
	  track->SetPadHit(i,hit[1]);
	  track->SetTimeHit(i,hit[2]);
	  
	  //if(hit[1]<0 || hit[2]>445)
	  //if(hit[2]<0 || hit[2]>445)
	  //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
	  //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
	}
    }
  fTracks->Compress();
  //cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
}

void AliL3Modeller::CheckForOverlaps()
{
  //Flag the tracks that overlap
  
  //cout<<"Checking for overlaps...";
  Int_t counter=0;
  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
    {
      AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
      if(!track1) continue;
      for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
	{
	  AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
	  if(!track2) continue;
	  for(Int_t k=AliL3Transform::GetFirstRow(fPatch); k<=AliL3Transform::GetLastRow(fPatch); k++)
	    {
	      if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
		 track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
		continue;
	      
	      if(track1->GetOverlap(k)>=0 || track2->GetOverlap(k)>=0) continue;
	      
	      if(abs((Int_t)rint(track1->GetPadHit(k))-(Int_t)rint(track2->GetPadHit(k))) <= fPadOverlap &&
		 abs((Int_t)rint(track1->GetTimeHit(k))-(Int_t)rint(track2->GetTimeHit(k))) <= fTimeOverlap)
		{
		  track2->SetOverlap(k,i);
		  //track1->SetOverlap(k,j);
		  counter++;
		}
	    }
	}
    }
  //cout<<"found "<<counter<<" done"<<endl;
}


void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
{
  
  Float_t padw,timew;
  
  padw = AliL3Transform::GetPadPitchWidth(fPatch);
  
  Float_t charge = (Float_t)cl->fCharge;
  Float_t pad = (Float_t)cl->fPad/charge;
  Float_t time = (Float_t)cl->fTime/charge;
  Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
  
  //Save the sigmas in pad and time:
  
  sigmaY2 = (s2);// + 1./12);//*padw*padw;
  
  /*Constants added by offline
    if(s2 != 0)
    {
    sigmaY2 = sigmaY2*0.108;
    if(fPatch<3)
    sigmaY2 = sigmaY2*2.07;
    }
  */

  s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
  timew = AliL3Transform::GetZWidth();
  sigmaZ2 = (s2);// +1./12);//*timew*timew;
  

  /*Constants added by offline
    if(s2 != 0)
    {
    sigmaZ2 = sigmaZ2*0.169;
    if(fPatch < 3)
    sigmaZ2 = sigmaZ2*1.77;
    }
  */
}
Commit	Line	Data
735e167e	1	//$Id$
735e167e	2
029912b7	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
735e167e	4	//*-- Copyright &copy ASV
735e167e	5
4994a25d	6	#include "AliL3StandardIncludes.h"
735e167e	7
	8	#include "AliL3Modeller.h"
	9	#include "AliL3MemHandler.h"
8487f697	10	#ifdef use_aliroot
	11	#include "AliL3FileHandler.h"
	12	#endif
735e167e	13	#include "AliL3TrackArray.h"
	14	#include "AliL3ModelTrack.h"
	15	#include "AliL3DigitData.h"
	16	#include "AliL3Transform.h"
	17
4994a25d	18	#if GCCVERSION == 3
	19	using namespace std;
	20	#endif
	21
4a838220	22	//_____________________________________________________________
	23	// AliL3Modeller
	24	//
	25	// Class for modeling TPC data.
	26	//
	27	// This performs the cluster finding, based on track parameters.
	28	// Basically it propagates the tracks to all padrows, and looks
	29	// for a corresponding cluster. For the moment only cog is calculated,
	30	// and no deconvolution is done.
735e167e	31
	32	ClassImp(AliL3Modeller)
	33
	34	AliL3Modeller::AliL3Modeller()
	35	{
	36	fMemHandler=0;
	37	fTracks=0;
4a838220	38	fTrackThreshold=0;
2357bb38	39	SetOverlap();
2357bb38	40	SetTrackThreshold();
f64c4364	41	SetSearchRange();
735e167e	42	}
	43
	44
	45	AliL3Modeller::~AliL3Modeller()
	46	{
	47	if(fMemHandler)
	48	delete fMemHandler;
	49	if(fTracks)
	50	delete fTracks;
735e167e	51	}
735e167e	52
4994a25d	53	void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t trackdata,Char_t path,Bool_t houghtracks,Bool_t binary)
735e167e	54	{
	55	fSlice = slice;
	56	fPatch = patch;
2357bb38	57
be6ddb10	58	sprintf(fPath,"%s",path);
029912b7	59
735e167e	60	fTracks = new AliL3TrackArray("AliL3ModelTrack");
92a876e2	61
92a876e2	62	Char_t fname[100];
735e167e	63	AliL3MemHandler *file = new AliL3MemHandler();
8487f697	64	//sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
8487f697	65	sprintf(fname,"%s/tracks_ho_%d.raw",trackdata,fSlice);
92a876e2	66	if(!file->SetBinaryInput(fname))
735e167e	67	{
92a876e2	68	cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
735e167e	69	return;
	70	}
	71	file->Binary2TrackArray(fTracks);
	72	file->CloseBinaryInput();
	73	delete file;
	74
4a838220	75
	76	if(!houghtracks)
	77	fTracks->QSort();
	78
735e167e	79	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	80	{
	81	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	82	if(!track) continue;
	83	track->Init(fSlice,fPatch);
4a838220	84
	85	//Only if the tracks has been merged across sector boundaries:
	86	//if(!houghtracks)
	87	//track->Rotate(fSlice,kTRUE); //!!!!!!!!!!!!!!!!!!!
	88
735e167e	89	track->CalculateHelix();
	90	}
	91
	92	CalculateCrossingPoints();
4a838220	93
2357bb38	94	CheckForOverlaps();
2357bb38	95
8487f697	96	UInt_t ndigits=0;
	97	AliL3DigitRowData *digits=0;
	98	#ifdef use_aliroot
	99	fMemHandler = new AliL3FileHandler();
2357bb38	100	fMemHandler->Init(slice,patch);
8487f697	101	if(binary == kFALSE)
8487f697	102	{
2357bb38	103	sprintf(fname,"%s/digitfile.root",fPath);
8487f697	104	fMemHandler->SetAliInput(fname);
	105	digits = fMemHandler->AliDigits2Memory(ndigits);
	106	}
2357bb38	107	else
	108	{
	109	sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
	110	if(!fMemHandler->SetBinaryInput(fname))
	111	{
	112	cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
	113	return;
	114	}
	115	digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
	116	}
8487f697	117	#else
735e167e	118	fMemHandler = new AliL3MemHandler();
2357bb38	119	fMemHandler->Init(slice,patch);
8487f697	120	if(binary == kFALSE)
735e167e	121	{
8487f697	122	cerr<<"AliL3Modeller::Init : Compile with AliROOT if you want rootfile as input"<<endl;
735e167e	123	return;
735e167e	124	}
8487f697	125	else
	126	{
	127	sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
	128	if(!fMemHandler->SetBinaryInput(fname))
	129	{
	130	cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
	131	return;
	132	}
	133	}
	134	digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
	135	#endif
735e167e	136
	137	SetInputData(digits);
	138	}
	139
029912b7	140	void AliL3Modeller::FindClusters()
735e167e	141	{
f64c4364	142	cout<<"AliL3Modeller::FindClusters : Processing slice "<<fSlice<<" patch "<<fPatch<<endl;
735e167e	143	if(!fTracks)
	144	{
	145	cerr<<"AliL3Modeller::Process : No tracks"<<endl;
	146	return;
	147	}
	148	if(!fRowData)
	149	{
	150	cerr<<"AliL3Modeller::Process : No data "<<endl;
	151	return;
	152	}
	153
	154	AliL3DigitRowData *rowPt = fRowData;
	155	AliL3DigitData *digPt=0;
	156
4a838220	157	Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
4a838220	158	Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
029912b7	159	Int_t bounds = ntimes*npads;
029912b7	160	Digit *row = new Digit[bounds];
735e167e	161
735e167e	162	Int_t seq_charge;
029912b7	163	Int_t pad,time,index;
735e167e	164	Short_t charge;
	165	Cluster cluster;
	166
4a838220	167	for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
735e167e	168	{
8487f697	169	fCurrentPadRow = i;
735e167e	170	memset((void)row,0,ntimesnpads*sizeof(Digit));
735e167e	171	digPt = (AliL3DigitData*)rowPt->fDigitData;
2357bb38	172	//cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
735e167e	173	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	174	{
	175	pad = digPt[j].fPad;
	176	time = digPt[j].fTime;
	177	charge = digPt[j].fCharge;
	178	row[ntimes*pad+time].fCharge = charge;
	179	row[ntimes*pad+time].fUsed = kFALSE;
029912b7	180	//cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
735e167e	181	}
	182
	183	for(Int_t k=0; k<fTracks->GetNTracks(); k++)
	184	{
	185	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(k);
	186	if(!track) continue;
029912b7	187
029912b7	188	if(track->GetPadHit(i)<0 \|\| track->GetTimeHit(i)<0 \|\| track->GetOverlap(i)>=0)
95a00d93	189	{
4a838220	190	track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch.
95a00d93	191	continue;
95a00d93	192	}
735e167e	193
	194	Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
	195	Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
fe1c53fa	196
	197	if(!CheckCluster(row,hitpad,hittime)) //Not a good cluster.
	198	{
	199	track->SetCluster(i,0,0,0,0,0,0);
	200	continue;
	201	}
	202
029912b7	203	//cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
735e167e	204	pad = hitpad;
735e167e	205	time = hittime;
fe1c53fa	206
735e167e	207	Int_t padsign=-1;
	208	Int_t timesign=-1;
	209
	210	memset(&cluster,0,sizeof(Cluster));
f64c4364	211	//cout<<"Processing padrow "<<i<<" with hittime "<<hittime<<" hitpad "<<hitpad<<" charge "<<row[ntimespad + time].fCharge<<" index "<<ntimespad + time<<endl;
4a838220	212	Int_t npads=0;
f64c4364	213	Int_t last_mean = hittime;
735e167e	214	while(1)//Process this padrow
735e167e	215	{
4a838220	216	if(pad < 0 \|\| pad >= AliL3Transform::GetNPads(i))
029912b7	217	{
92a876e2	218	//cout<<"Pad = "<<pad<<" on row "<<i<<endl;
4a838220	219	FillCluster(track,&cluster,i,npads);
029912b7	220	break;
029912b7	221	}
735e167e	222	seq_charge=0;
	223	timesign=-1;
	224	time = hittime;
4a838220	225
735e167e	226	while(1) //Process sequence on this pad:
735e167e	227	{
f64c4364	228	if(time < 0 \|\| time >= AliL3Transform::GetNTimeBins())
	229	break;
	230
029912b7	231	index = ntimes*pad + time;
	232	if(index < 0 \|\| index >= bounds)
	233	{
	234	cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
	235	<<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
	236	break;
	237	}
4a838220	238
029912b7	239	charge = row[index].fCharge;
fe1c53fa	240
fe1c53fa	241	if(charge==0) //Empty timebin, so try and expand the search within the limits on this pad
029912b7	242	{
fe1c53fa	243	if(seq_charge==0 && abs(time-hittime) <= fTimeSearch)//No sequence yet, keep looking
029912b7	244	{
fe1c53fa	245	if(timesign==-1)
	246	time--;
	247	else
	248	time++;
029912b7	249	continue;
029912b7	250	}
fe1c53fa	251	else if(timesign==-1) //Switch search direction
029912b7	252	{
fe1c53fa	253	time = hittime + 1;
fe1c53fa	254	timesign = 1;
029912b7	255	continue;
029912b7	256	}
fe1c53fa	257	else //ok, boundary reached, this pad is done.
fe1c53fa	258	break;
029912b7	259	}
029912b7	260
2357bb38	261	if(row[ntimes*pad+time].fUsed==kTRUE) //Don't use digits several times. This leads to mult. rec.tracks.
	262	{
	263	time += timesign;
	264	continue;
	265	}
735e167e	266
735e167e	267	seq_charge += charge;
2357bb38	268
2357bb38	269	//Update the cluster parameters with this timebin
735e167e	270	cluster.fTime += time*charge;
	271	cluster.fPad += pad*charge;
	272	cluster.fCharge += charge;
	273	cluster.fSigmaY2 += padpadcharge;
	274	cluster.fSigmaZ2 += timetimecharge;
	275
	276	row[ntimes*pad+time].fUsed = kTRUE;
	277	time += timesign;
	278	}
4a838220	279
f64c4364	280	//Always compare with the current mean in time of the cluster under construction.
	281	if(cluster.fCharge)
	282	last_mean = (Int_t)rint((Float_t)(cluster.fTime/cluster.fCharge));
4a838220	283
2357bb38	284	if(seq_charge)//There was something on this pad, so keep looking on the neighbouring pad
4a838220	285	{
	286	pad += padsign;
	287	npads++;
	288	}
fe1c53fa	289	else
735e167e	290	{
fe1c53fa	291	if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadSearch && pad > 0)
95a00d93	292	{
fe1c53fa	293	if(padsign==-1)
	294	pad--;
	295	else
	296	pad++;
	297	continue;
95a00d93	298	}
fe1c53fa	299	else if(padsign==-1)
95a00d93	300	{
fe1c53fa	301	pad = hitpad + 1;
	302	padsign = 1;
	303	continue;
95a00d93	304	}
fe1c53fa	305	else
735e167e	306	{
4a838220	307	FillCluster(track,&cluster,i,npads);
735e167e	308	break;
fe1c53fa	309	}
735e167e	310	}
735e167e	311	}
029912b7	312	//cout<<"done"<<endl;
735e167e	313	}
95a00d93	314	FillZeros(rowPt,row);
735e167e	315	fMemHandler->UpdateRowPointer(rowPt);
735e167e	316	}
735e167e	317	delete [] row;
f64c4364	318	//cout<<"done processing"<<endl;
029912b7	319
	320
	321	//Debug:
	322	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	323	{
	324	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	325	if(!track) continue;
4a838220	326	if(track->GetNClusters() != AliL3Transform::GetNRows(fPatch))
4a838220	327	cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<AliL3Transform::GetNRows(fPatch)<<endl<<endl;
029912b7	328	}
	329
	330	}
	331
fe1c53fa	332	Bool_t AliL3Modeller::CheckCluster(Digit *row,Int_t hitpad,Int_t hittime)
	333	{
	334	Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
	335	Bool_t seq_was_falling=kFALSE;
	336	Int_t last_seq_charge=0;
	337	Int_t cluster_charge=0;
	338	Int_t time,pad,charge,padsign=-1,timesign=-1;
	339	time = hittime;
	340	pad = hitpad;
	341	//cout<<"Checking cluster "<<hitpad<<" hittime "<<hittime<<endl;
	342	while(1)
	343	{
	344
	345	Bool_t last_was_falling=kFALSE;
	346	Int_t last_charge=0;
	347	Int_t seq_charge=0;
	348	while(1)
	349	{
	350
	351	charge = row[(ntimes*pad+time)].fCharge;
	352	//cout<<"Charge "<<charge<<" on pad "<<pad<<" time "<<time<<endl;
	353	if(charge==0)
	354	{
	355	if(seq_charge==0 && abs(time-hittime) <= fTimeSearch)
	356	{
	357	if(timesign==-1)
	358	time--;
	359	else
	360	time++;
	361	continue;
	362	}
	363	else if(timesign==-1)
	364	{
	365	last_charge = row[(ntimes*hitpad+hittime)].fCharge;
	366	time = hittime+1;
	367	timesign=1;
	368	if(last_charge < row[(ntimes*hitpad+(hittime-1))].fCharge)
	369	last_was_falling=kTRUE;
	370	else
	371	last_was_falling=kFALSE;
	372	continue;
	373	}
	374	else
	375	break;
	376	}
	377
	378	if(charge > last_charge)
	379	{
	380	if(last_was_falling) //This is a overlapping cluster.
	381	return kFALSE;
	382	}
	383	else
	384	last_was_falling = kTRUE;
	385
	386	cluster_charge+=charge;
	387	seq_charge += charge;
	388	last_charge = charge;
	389	time += timesign;
	390	}
	391
	392	if(seq_charge)
	393	pad += padsign;
	394	else
	395	{
396	if(cluster_charge==0 && abs(pad-hitpad)<=fPadSearch && pad > 0)
397	{
398	if(padsign==-1)
399	pad--;
400	else
401	pad++;
402	continue;
403	}
404	else if(padsign==-1)
405	{
406	pad = hitpad+1;
407	padsign = 1;
408	seq_was_falling=kFALSE;
409	continue;
410	}
411	else
412	break;
413	}
414
415	if(seq_charge > last_seq_charge)
416	{
417	if(seq_was_falling)
418	return kFALSE;
419	else
420	seq_was_falling=kTRUE;
421	}
422	last_seq_charge = seq_charge;
423	}
424
425	return kTRUE;
426	}
427
4a838220	428	void AliL3Modeller::FillCluster(AliL3ModelTrack track,Cluster cluster,Int_t row,Int_t npads)
029912b7	429	{
	430	if(cluster->fCharge==0)
	431	{
4a838220	432	track->SetCluster(row,0,0,0,0,0,0);
029912b7	433	return;
	434	}
	435	Float_t fcharge = (Float_t)cluster->fCharge;
	436	Float_t fpad = ((Float_t)cluster->fPad/fcharge);
	437	Float_t ftime = ((Float_t)cluster->fTime/fcharge);
	438	Float_t sigmaY2,sigmaZ2;
	439	CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
4a838220	440	track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2,npads);
735e167e	441	}
735e167e	442
95a00d93	443	void AliL3Modeller::FillZeros(AliL3DigitRowData rowPt,Digit row)
	444	{
	445	//Fill zero where data has been used.
	446
4a838220	447	Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
95a00d93	448	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	449	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	450	{
	451	Int_t pad = digPt[j].fPad;
	452	Int_t time = digPt[j].fTime;
	453	if(row[ntimes*pad+time].fUsed==kTRUE)
	454	digPt[j].fCharge = 0;
	455	}
	456	}
	457
be6ddb10	458	void AliL3Modeller::WriteRemaining()
95a00d93	459	{
	460	//Write remaining (nonzero) digits to file.
	461
95a00d93	462	AliL3DigitRowData *rowPt;
	463	rowPt = (AliL3DigitRowData*)fRowData;
	464	Int_t digitcount=0;
4a838220	465	Int_t ndigits[(AliL3Transform::GetNRows(fPatch))];
4a838220	466	for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
95a00d93	467	{
95a00d93	468	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
4a838220	469	ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]=0;
95a00d93	470	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	471	{
	472	if(digPt[j].fCharge==0) continue;
	473	digitcount++;
4a838220	474	ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]++;
95a00d93	475	}
4a838220	476	//cout<<"Difference "<<(int)ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]<<" "<<(int)rowPt->fNDigit<<endl;
95a00d93	477	fMemHandler->UpdateRowPointer(rowPt);
	478	}
	479
4a838220	480	Int_t size = digitcountsizeof(AliL3DigitData) + AliL3Transform::GetNRows(fPatch)sizeof(AliL3DigitRowData);
95a00d93	481	Byte_t *data = new Byte_t[size];
	482	memset(data,0,size);
	483	AliL3DigitRowData tempPt = (AliL3DigitRowData)data;
	484	rowPt = (AliL3DigitRowData*)fRowData;
	485
4a838220	486	for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
95a00d93	487	{
	488	Int_t localcount=0;
	489	tempPt->fRow = i;
4a838220	490	tempPt->fNDigit = ndigits[(i-AliL3Transform::GetFirstRow(fPatch))];
95a00d93	491	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	492	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	493	{
	494	if(digPt[j].fCharge==0) continue;
4a838220	495	if(localcount >= ndigits[(i-AliL3Transform::GetFirstRow(fPatch))])
95a00d93	496	{
	497	cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
	498	return;
	499	}
	500	tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
	501	tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
	502	tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
4a838220	503
95a00d93	504	localcount++;
95a00d93	505	}
4a838220	506	if(ndigits[(i-AliL3Transform::GetFirstRow(fPatch))] != localcount)
95a00d93	507	{
	508	cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
	509	return;
	510	}
	511	fMemHandler->UpdateRowPointer(rowPt);
	512	Byte_t tmp = (Byte_t)tempPt;
4a838220	513	Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]*sizeof(AliL3DigitData);
95a00d93	514	tmp += size;
	515	tempPt = (AliL3DigitRowData*)tmp;
	516	}
	517
be6ddb10	518	Char_t fname[100];
95a00d93	519	AliL3MemHandler *mem = new AliL3MemHandler();
2357bb38	520	sprintf(fname,"%s/comp/remains_%d_%d.raw",fPath,fSlice,fPatch);
be6ddb10	521	mem->SetBinaryOutput(fname);
4a838220	522	mem->Memory2CompBinary((UInt_t)AliL3Transform::GetNRows(fPatch),(AliL3DigitRowData*)data);
95a00d93	523	mem->CloseBinaryOutput();
95a00d93	524	delete mem;
029912b7	525	delete [] data;
95a00d93	526	}
	527
	528
735e167e	529	void AliL3Modeller::CalculateCrossingPoints()
735e167e	530	{
f64c4364	531	//cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
735e167e	532	if(!fTracks)
	533	{
	534	cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
	535	return;
	536	}
	537	Float_t hit[3];
029912b7	538
029912b7	539	Int_t sector,row;
4a838220	540	for(Int_t i=AliL3Transform::GetLastRow(fPatch); i>=AliL3Transform::GetFirstRow(fPatch); i--)
735e167e	541	{
	542	for(Int_t j=0; j<fTracks->GetNTracks(); j++)
	543	{
	544	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(j);
	545	if(!track) continue;
4a838220	546
95a00d93	547	if(!track->GetCrossingPoint(i,hit))
735e167e	548	{
f64c4364	549	//cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
	550	//" pt "<<track->GetPt()<<
	551	//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
029912b7	552	//"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
	553	//"--------"<<endl;
	554	fTracks->Remove(j);
735e167e	555	continue;
735e167e	556	}
2357bb38	557	//cout<<"X "<<hit[0]<<" Y "<<hit[1]<<" Z "<<hit[2]<<" tgl "<<track->GetTgl()<<endl;
95a00d93	558
4a838220	559	AliL3Transform::Slice2Sector(fSlice,i,sector,row);
4a838220	560	AliL3Transform::Local2Raw(hit,sector,row);
2357bb38	561	//cout<<"Pad "<<hit[1]<<" time "<<hit[2]<<" in sector "<<sector<<" row "<<row<<endl;
4a838220	562	if(hit[1]<0 \|\| hit[1]>AliL3Transform::GetNPads(i) \|\|
4a838220	563	hit[2]<0 \|\| hit[2]>AliL3Transform::GetNTimeBins())
95a00d93	564	{//Track is leaving the patch, so flag the track hits (<0)
	565	track->SetPadHit(i,-1);
	566	track->SetTimeHit(i,-1);
	567	continue;
	568	}
4a838220	569
4a838220	570
735e167e	571	track->SetPadHit(i,hit[1]);
735e167e	572	track->SetTimeHit(i,hit[2]);
95a00d93	573
735e167e	574	//if(hit[1]<0 \|\| hit[2]>445)
95a00d93	575	//if(hit[2]<0 \|\| hit[2]>445)
95a00d93	576	//cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
735e167e	577	//cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
	578	}
	579	}
	580	fTracks->Compress();
f64c4364	581	//cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
735e167e	582	}
	583
	584	void AliL3Modeller::CheckForOverlaps()
	585	{
	586	//Flag the tracks that overlap
	587
f64c4364	588	//cout<<"Checking for overlaps...";
2357bb38	589	Int_t counter=0;
735e167e	590	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	591	{
	592	AliL3ModelTrack track1 = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	593	if(!track1) continue;
	594	for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
	595	{
	596	AliL3ModelTrack track2 = (AliL3ModelTrack)fTracks->GetCheckedTrack(j);
	597	if(!track2) continue;
4a838220	598	for(Int_t k=AliL3Transform::GetFirstRow(fPatch); k<=AliL3Transform::GetLastRow(fPatch); k++)
735e167e	599	{
95a00d93	600	if(track1->GetPadHit(k)<0 \|\| track1->GetTimeHit(k)<0 \|\|
	601	track2->GetPadHit(k)<0 \|\| track2->GetTimeHit(k)<0)
	602	continue;
2357bb38	603
	604	if(track1->GetOverlap(k)>=0 \|\| track2->GetOverlap(k)>=0) continue;
	605
	606	if(abs((Int_t)rint(track1->GetPadHit(k))-(Int_t)rint(track2->GetPadHit(k))) <= fPadOverlap &&
	607	abs((Int_t)rint(track1->GetTimeHit(k))-(Int_t)rint(track2->GetTimeHit(k))) <= fTimeOverlap)
735e167e	608	{
029912b7	609	track2->SetOverlap(k,i);
2357bb38	610	//track1->SetOverlap(k,j);
2357bb38	611	counter++;
735e167e	612	}
	613	}
	614	}
	615	}
f64c4364	616	//cout<<"found "<<counter<<" done"<<endl;
735e167e	617	}
	618
	619
	620	void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
	621	{
	622
	623	Float_t padw,timew;
4a838220	624
	625	padw = AliL3Transform::GetPadPitchWidth(fPatch);
	626
735e167e	627	Float_t charge = (Float_t)cl->fCharge;
	628	Float_t pad = (Float_t)cl->fPad/charge;
	629	Float_t time = (Float_t)cl->fTime/charge;
	630	Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
4a838220	631
	632	//Save the sigmas in pad and time:
	633
	634	sigmaY2 = (s2);// + 1./12);//padwpadw;
	635
	636	/*Constants added by offline
	637	if(s2 != 0)
735e167e	638	{
4a838220	639	sigmaY2 = sigmaY2*0.108;
	640	if(fPatch<3)
	641	sigmaY2 = sigmaY2*2.07;
735e167e	642	}
4a838220	643	*/
4a838220	644
735e167e	645	s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
4a838220	646	timew = AliL3Transform::GetZWidth();
	647	sigmaZ2 = (s2);// +1./12);//timewtimew;
	648
	649
	650
	651	/*Constants added by offline
	652	if(s2 != 0)
735e167e	653	{
4a838220	654	sigmaZ2 = sigmaZ2*0.169;
	655	if(fPatch < 3)
	656	sigmaZ2 = sigmaZ2*1.77;
735e167e	657	}
4a838220	658	*/
735e167e	659	}
f64c4364	660