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

//$Id$

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

#include <stream.h>
#include <iostream.h>
#include <math.h>

#include "AliL3Modeller.h"
#include "AliL3MemHandler.h"
#include "AliL3TrackArray.h"
#include "AliL3ModelTrack.h"
#include "AliL3DigitData.h"
#include "AliL3Transform.h"

#include "AliL3Defs.h"

ClassImp(AliL3Modeller)

AliL3Modeller::AliL3Modeller()
{
  fMemHandler=0;
  fTracks=0;
  fTransform=0;
}


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

}

void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path)
{
  fSlice = slice;
  fPatch = patch;
  fPadOverlap=6;
  fTimeOverlap=8;
  fTrackThreshold=10;
  sprintf(fPath,"%s",path);
  
  fTransform = new AliL3Transform();
  fTracks = new AliL3TrackArray("AliL3ModelTrack");
  
  Char_t fname[100];
  AliL3MemHandler *file = new AliL3MemHandler();
  sprintf(fname,"%s/tracks.raw",trackdata);
  if(!file->SetBinaryInput(fname))
    {
      cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
      return;
    }
  file->Binary2TrackArray(fTracks);
  file->CloseBinaryInput();
  delete file;
  
  fTracks->QSort();
  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
    {
      AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
      if(!track) continue;
      track->Init(fSlice,fPatch);
      track->Rotate(fSlice,kTRUE);
      track->CalculateHelix();
    }    
  
  CalculateCrossingPoints();
  CheckForOverlaps();

  fMemHandler = new AliL3MemHandler();
  sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
  if(!fMemHandler->SetBinaryInput(fname))
    {
      cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
      return;
    }
  UInt_t ndigits;
  AliL3DigitRowData *digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
  
  SetInputData(digits);
}

void AliL3Modeller::FindClusters()
{
  
  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 = fTransform->GetNTimeBins()+1;
  Int_t npads = fTransform->GetNPads(NRows[fPatch][1])+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=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
    {
      memset((void*)row,0,ntimes*npads*sizeof(Digit));
      digPt = (AliL3DigitData*)rowPt->fDigitData;
      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->GetOverlap(i)>=0) continue;//Track is overlapping
	  
	  if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetOverlap(i)>=0)
	    {
	      track->SetCluster(i,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));
	  //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));
	  
	  while(1)//Process this padrow
	    {
	      if(pad < 0 || pad >= fTransform->GetNPads(i)) 
		{
		  //cout<<"Pad = "<<pad<<" on row "<<i<<endl;
		  FillCluster(track,&cluster,i);
		  break;
		}
	      seq_charge=0;
	      timesign=-1;
	      time = hittime;
	      while(1) //Process sequence on this pad:
		{
		  if(time < 0) 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 && timesign==-1)
		    {
		      if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
			{
			  time--;
			  continue;
			}
		      else
			{
			  time = hittime+1;
			  timesign=1;
			  continue;
			}
		    }
		  else if(charge==0 && timesign==1)
		    {
		      if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
			{
			  time++;
			  continue;
			}
		      else
			{
			  //cerr<<"Breaking off at pad "<<pad<<" and time "<<time<<endl;
			  break;
			}
		    }
		  
		  //cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
		  
		  seq_charge += charge;
		  		  		  
		  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;
		}
	      //cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
	      if(seq_charge)
		pad += padsign;
	      else //Nothing more on this pad, goto next pad
		{
		  if(padsign==-1) 
		    {
		      if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2 && pad > 0)
			{
			  pad--; //In this case, we haven't found anything yet, 
			}        //so we will try to expand our search within the natural boundaries.
		      else
			{
			  pad=hitpad+1; 
			  padsign=1; 
			}
		      continue;
		    }
		  
		  else if(padsign==1)
		    {
		      if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2 && pad < fTransform->GetNPads(i)-2)
			{
			  pad++;     //In this case, we haven't found anything yet, 
			  continue;  //so we will try to expand our search within the natural boundaries.
			}
		      else //We are out of range, or cluster if finished.
			{
			  //cout<<"Outof range; charge "<<cluster.fCharge<<" paddiff "<<abs(pad-hitpad)<<endl;
			  FillCluster(track,&cluster,i);
			  break;
			}
		    }
		  else //Nothing more in this cluster
		    {
		      //cout<<"Filling final cluster"<<endl;
		      FillCluster(track,&cluster,i);
		      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() != NumRows[fPatch])
	cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<NumRows[fPatch]<<endl<<endl;
    }
  
}

void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row)
{
  if(cluster->fCharge==0)
    {
      track->SetCluster(row,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);
  
}

void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Digit *row)
{
  //Fill zero where data has been used.
  
  Int_t ntimes = fTransform->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[(NumRows[fPatch])];
  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
    {
      AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
      ndigits[(i-NRows[fPatch][0])]=0;
      for(UInt_t j=0; j<rowPt->fNDigit; j++)
	{
	  if(digPt[j].fCharge==0) continue;
	  
	  digitcount++;
	  ndigits[(i-NRows[fPatch][0])]++;
	}
      //cout<<"Difference "<<(int)ndigits[(i-NRows[fPatch][0])]<<" "<<(int)rowPt->fNDigit<<endl;
      fMemHandler->UpdateRowPointer(rowPt);
    }
  
  Int_t size = digitcount*sizeof(AliL3DigitData) + NumRows[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=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
    {
      Int_t localcount=0;
      tempPt->fRow = i;
      tempPt->fNDigit = ndigits[(i-NRows[fPatch][0])];
      AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
      for(UInt_t j=0; j<rowPt->fNDigit; j++)
	{
	  if(digPt[j].fCharge==0) continue;
	  if(localcount >= ndigits[(i-NRows[fPatch][0])])
	    {
	      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-NRows[fPatch][0])] != 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-NRows[fPatch][0])]*sizeof(AliL3DigitData);
      tmp += size;
      tempPt = (AliL3DigitRowData*)tmp;
    }

  Char_t fname[100];
  AliL3MemHandler *mem = new AliL3MemHandler();
  sprintf(fname,"%s/remains_%d_%d.raw",fPath,fSlice,fPatch);
  mem->SetBinaryOutput(fname);
  mem->Memory2CompBinary((UInt_t)NumRows[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];
  
  
  //Remove tracks which are no good:
  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
    {
      AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
      if(!track) continue;
      if(track->GetFirstPointX() > fTransform->Row2X(NRows[fPatch][1]) || track->GetPt()<0.1)
	fTracks->Remove(i);
    }
  
  Int_t sector,row;
  for(Int_t i=NRows[fPatch][1]; i>=NRows[fPatch][0]; i--)
    {
      for(Int_t j=0; j<fTracks->GetNTracks(); j++)
	{
	  AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
	  if(!track) continue;
	  	  
	  if(track->GetNHits() < fTrackThreshold)
	    {
	      fTracks->Remove(j);
	      continue;
	    }
	  
	  if(!track->GetCrossingPoint(i,hit)) 
	    {
	      cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
		"First point "<<track->GetFirstPointX()<<
		" nhits "<<track->GetNHits()<<endl;//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
		//"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
		//"--------"<<endl;
	      fTracks->Remove(j);
	      continue;
	    }
	  //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
	  
	  fTransform->Slice2Sector(fSlice,i,sector,row);
	  fTransform->Local2Raw(hit,sector,row);
	  if(hit[1]<0 || hit[1]>fTransform->GetNPads(i) ||
	     hit[2]<0 || hit[2]>fTransform->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...";
  
  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=NRows[fPatch][0]; k<=NRows[fPatch][1]; k++)
	    {
	      if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
		 track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
		continue;
	      if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) <= fPadOverlap &&
		 fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) <= fTimeOverlap)
		{
		  track2->SetOverlap(k,i);
		  track1->SetOverlap(k,j);
		}
	    }
	}
    }
  cout<<"done"<<endl;
}


void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
{
  
  Float_t padw,timew;
  if(fPatch < 3)
    padw = fTransform->GetPadPitchWidthLow();
  else
    padw = fTransform->GetPadPitchWidthUp();
  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;
  sigmaY2 = (s2 + 1./12)*padw*padw;

  if(s2 != 0)
    {
      sigmaY2 = sigmaY2*0.108;
      if(fPatch<3)
	sigmaY2 = sigmaY2*2.07;
    }
  
  s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
  timew = fTransform->GetZWidth();
  sigmaZ2 = (s2 +1./12)*timew*timew;
  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
	5
	6	#include <stream.h>
	7	#include <iostream.h>
	8	#include <math.h>
	9
	10	#include "AliL3Modeller.h"
	11	#include "AliL3MemHandler.h"
	12	#include "AliL3TrackArray.h"
	13	#include "AliL3ModelTrack.h"
	14	#include "AliL3DigitData.h"
	15	#include "AliL3Transform.h"
	16
	17	#include "AliL3Defs.h"
	18
	19	ClassImp(AliL3Modeller)
	20
	21	AliL3Modeller::AliL3Modeller()
	22	{
	23	fMemHandler=0;
	24	fTracks=0;
	25	fTransform=0;
	26	}
	27
	28
	29	AliL3Modeller::~AliL3Modeller()
	30	{
	31	if(fMemHandler)
	32	delete fMemHandler;
	33	if(fTracks)
	34	delete fTracks;
	35	if(fTransform)
	36	delete fTransform;
	37
	38	}
	39
92a876e2	40	void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t trackdata,Char_t path)
735e167e	41	{
	42	fSlice = slice;
	43	fPatch = patch;
92a876e2	44	fPadOverlap=6;
92a876e2	45	fTimeOverlap=8;
029912b7	46	fTrackThreshold=10;
be6ddb10	47	sprintf(fPath,"%s",path);
029912b7	48
735e167e	49	fTransform = new AliL3Transform();
735e167e	50	fTracks = new AliL3TrackArray("AliL3ModelTrack");
92a876e2	51
92a876e2	52	Char_t fname[100];
735e167e	53	AliL3MemHandler *file = new AliL3MemHandler();
92a876e2	54	sprintf(fname,"%s/tracks.raw",trackdata);
92a876e2	55	if(!file->SetBinaryInput(fname))
735e167e	56	{
92a876e2	57	cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
735e167e	58	return;
	59	}
	60	file->Binary2TrackArray(fTracks);
	61	file->CloseBinaryInput();
	62	delete file;
	63
	64	fTracks->QSort();
	65	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	66	{
	67	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	68	if(!track) continue;
	69	track->Init(fSlice,fPatch);
	70	track->Rotate(fSlice,kTRUE);
	71	track->CalculateHelix();
	72	}
	73
	74	CalculateCrossingPoints();
	75	CheckForOverlaps();
	76
	77	fMemHandler = new AliL3MemHandler();
be6ddb10	78	sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
735e167e	79	if(!fMemHandler->SetBinaryInput(fname))
	80	{
	81	cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
	82	return;
	83	}
	84	UInt_t ndigits;
	85	AliL3DigitRowData digits=(AliL3DigitRowData)fMemHandler->CompBinary2Memory(ndigits);
	86
	87	SetInputData(digits);
	88	}
	89
029912b7	90	void AliL3Modeller::FindClusters()
735e167e	91	{
	92
	93	if(!fTracks)
	94	{
	95	cerr<<"AliL3Modeller::Process : No tracks"<<endl;
	96	return;
	97	}
	98	if(!fRowData)
	99	{
	100	cerr<<"AliL3Modeller::Process : No data "<<endl;
	101	return;
	102	}
	103
	104	AliL3DigitRowData *rowPt = fRowData;
	105	AliL3DigitData *digPt=0;
	106
	107	Int_t ntimes = fTransform->GetNTimeBins()+1;
	108	Int_t npads = fTransform->GetNPads(NRows[fPatch][1])+1;//Max num of pads.
029912b7	109	Int_t bounds = ntimes*npads;
029912b7	110	Digit *row = new Digit[bounds];
735e167e	111
735e167e	112	Int_t seq_charge;
029912b7	113	Int_t pad,time,index;
735e167e	114	Short_t charge;
	115	Cluster cluster;
	116
	117	for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
	118	{
	119	memset((void)row,0,ntimesnpads*sizeof(Digit));
	120	digPt = (AliL3DigitData*)rowPt->fDigitData;
	121	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	122	{
	123	pad = digPt[j].fPad;
	124	time = digPt[j].fTime;
	125	charge = digPt[j].fCharge;
	126	row[ntimes*pad+time].fCharge = charge;
	127	row[ntimes*pad+time].fUsed = kFALSE;
029912b7	128	//cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
735e167e	129	}
	130
	131	for(Int_t k=0; k<fTracks->GetNTracks(); k++)
	132	{
	133	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(k);
	134	if(!track) continue;
029912b7	135	//if(track->GetOverlap(i)>=0) continue;//Track is overlapping
	136
	137	if(track->GetPadHit(i)<0 \|\| track->GetTimeHit(i)<0 \|\| track->GetOverlap(i)>=0)
95a00d93	138	{
029912b7	139	track->SetCluster(i,0,0,0,0,0); //The track has left the patch.
95a00d93	140	continue;
95a00d93	141	}
735e167e	142
	143	Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
	144	Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
029912b7	145	//cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
735e167e	146	pad = hitpad;
	147	time = hittime;
	148	Int_t padsign=-1;
	149	Int_t timesign=-1;
	150
	151	memset(&cluster,0,sizeof(Cluster));
	152
	153	while(1)//Process this padrow
	154	{
029912b7	155	if(pad < 0 \|\| pad >= fTransform->GetNPads(i))
029912b7	156	{
92a876e2	157	//cout<<"Pad = "<<pad<<" on row "<<i<<endl;
029912b7	158	FillCluster(track,&cluster,i);
	159	break;
	160	}
735e167e	161	seq_charge=0;
	162	timesign=-1;
	163	time = hittime;
	164	while(1) //Process sequence on this pad:
	165	{
029912b7	166	if(time < 0) break;
	167	index = ntimes*pad + time;
	168	if(index < 0 \|\| index >= bounds)
	169	{
	170	cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
	171	<<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
	172	break;
	173	}
	174	charge = row[index].fCharge;
735e167e	175	if(charge==0 && timesign==-1)
029912b7	176	{
	177	if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
	178	{
	179	time--;
	180	continue;
	181	}
	182	else
	183	{
	184	time = hittime+1;
	185	timesign=1;
	186	continue;
	187	}
	188	}
735e167e	189	else if(charge==0 && timesign==1)
029912b7	190	{
	191	if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
	192	{
	193	time++;
	194	continue;
	195	}
	196	else
92a876e2	197	{
	198	//cerr<<"Breaking off at pad "<<pad<<" and time "<<time<<endl;
	199	break;
	200	}
029912b7	201	}
029912b7	202
735e167e	203	//cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
	204
	205	seq_charge += charge;
	206
	207	cluster.fTime += time*charge;
	208	cluster.fPad += pad*charge;
	209	cluster.fCharge += charge;
	210	cluster.fSigmaY2 += padpadcharge;
	211	cluster.fSigmaZ2 += timetimecharge;
	212
	213	row[ntimes*pad+time].fUsed = kTRUE;
	214	time += timesign;
	215	}
95a00d93	216	//cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
735e167e	217	if(seq_charge)
735e167e	218	pad += padsign;
735e167e	219	else //Nothing more on this pad, goto next pad
	220	{
	221	if(padsign==-1)
95a00d93	222	{
92a876e2	223	if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2 && pad > 0)
95a00d93	224	{
	225	pad--; //In this case, we haven't found anything yet,
	226	} //so we will try to expand our search within the natural boundaries.
	227	else
	228	{
	229	pad=hitpad+1;
	230	padsign=1;
	231	}
	232	continue;
	233	}
029912b7	234
029912b7	235	else if(padsign==1)
95a00d93	236	{
92a876e2	237	if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2 && pad < fTransform->GetNPads(i)-2)
029912b7	238	{
	239	pad++; //In this case, we haven't found anything yet,
	240	continue; //so we will try to expand our search within the natural boundaries.
	241	}
	242	else //We are out of range, or cluster if finished.
	243	{
	244	//cout<<"Outof range; charge "<<cluster.fCharge<<" paddiff "<<abs(pad-hitpad)<<endl;
	245	FillCluster(track,&cluster,i);
	246	break;
	247	}
95a00d93	248	}
735e167e	249	else //Nothing more in this cluster
735e167e	250	{
029912b7	251	//cout<<"Filling final cluster"<<endl;
029912b7	252	FillCluster(track,&cluster,i);
735e167e	253	break;
	254	}
	255	}
735e167e	256	}
029912b7	257	//cout<<"done"<<endl;
735e167e	258	}
029912b7	259
95a00d93	260	FillZeros(rowPt,row);
735e167e	261	fMemHandler->UpdateRowPointer(rowPt);
735e167e	262	}
735e167e	263	delete [] row;
029912b7	264	cout<<"done processing"<<endl;
	265
	266
	267	//Debug:
	268	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	269	{
	270	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	271	if(!track) continue;
	272	if(track->GetNClusters() != NumRows[fPatch])
	273	cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<NumRows[fPatch]<<endl<<endl;
	274	}
	275
	276	}
	277
	278	void AliL3Modeller::FillCluster(AliL3ModelTrack track,Cluster cluster,Int_t row)
	279	{
	280	if(cluster->fCharge==0)
	281	{
	282	track->SetCluster(row,0,0,0,0,0);
	283	return;
	284	}
	285	Float_t fcharge = (Float_t)cluster->fCharge;
	286	Float_t fpad = ((Float_t)cluster->fPad/fcharge);
	287	Float_t ftime = ((Float_t)cluster->fTime/fcharge);
	288	Float_t sigmaY2,sigmaZ2;
	289	CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
	290	track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2);
735e167e	291
	292	}
	293
95a00d93	294	void AliL3Modeller::FillZeros(AliL3DigitRowData rowPt,Digit row)
	295	{
	296	//Fill zero where data has been used.
	297
	298	Int_t ntimes = fTransform->GetNTimeBins()+1;
	299	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	300	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	301	{
	302	Int_t pad = digPt[j].fPad;
	303	Int_t time = digPt[j].fTime;
	304	if(row[ntimes*pad+time].fUsed==kTRUE)
	305	digPt[j].fCharge = 0;
	306	}
	307	}
	308
be6ddb10	309	void AliL3Modeller::WriteRemaining()
95a00d93	310	{
	311	//Write remaining (nonzero) digits to file.
	312
95a00d93	313	AliL3DigitRowData *rowPt;
	314	rowPt = (AliL3DigitRowData*)fRowData;
	315	Int_t digitcount=0;
	316	Int_t ndigits[(NumRows[fPatch])];
029912b7	317	for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
95a00d93	318	{
	319	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	320	ndigits[(i-NRows[fPatch][0])]=0;
	321	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	322	{
	323	if(digPt[j].fCharge==0) continue;
92a876e2	324
95a00d93	325	digitcount++;
	326	ndigits[(i-NRows[fPatch][0])]++;
	327	}
	328	//cout<<"Difference "<<(int)ndigits[(i-NRows[fPatch][0])]<<" "<<(int)rowPt->fNDigit<<endl;
	329	fMemHandler->UpdateRowPointer(rowPt);
	330	}
	331
	332	Int_t size = digitcountsizeof(AliL3DigitData) + NumRows[fPatch]sizeof(AliL3DigitRowData);
	333	Byte_t *data = new Byte_t[size];
	334	memset(data,0,size);
	335	AliL3DigitRowData tempPt = (AliL3DigitRowData)data;
	336	rowPt = (AliL3DigitRowData*)fRowData;
	337
029912b7	338	for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
95a00d93	339	{
	340	Int_t localcount=0;
	341	tempPt->fRow = i;
	342	tempPt->fNDigit = ndigits[(i-NRows[fPatch][0])];
	343	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	344	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	345	{
	346	if(digPt[j].fCharge==0) continue;
	347	if(localcount >= ndigits[(i-NRows[fPatch][0])])
	348	{
	349	cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
	350	return;
	351	}
	352	tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
	353	tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
	354	tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
	355	localcount++;
	356	}
	357	if(ndigits[(i-NRows[fPatch][0])] != localcount)
	358	{
	359	cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
	360	return;
	361	}
	362	fMemHandler->UpdateRowPointer(rowPt);
	363	Byte_t tmp = (Byte_t)tempPt;
	364	Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-NRows[fPatch][0])]*sizeof(AliL3DigitData);
	365	tmp += size;
	366	tempPt = (AliL3DigitRowData*)tmp;
	367	}
	368
be6ddb10	369	Char_t fname[100];
95a00d93	370	AliL3MemHandler *mem = new AliL3MemHandler();
be6ddb10	371	sprintf(fname,"%s/remains_%d_%d.raw",fPath,fSlice,fPatch);
be6ddb10	372	mem->SetBinaryOutput(fname);
95a00d93	373	mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
	374	mem->CloseBinaryOutput();
	375	delete mem;
029912b7	376	delete [] data;
95a00d93	377	}
	378
	379
735e167e	380	void AliL3Modeller::CalculateCrossingPoints()
	381	{
	382	cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
	383	if(!fTracks)
	384	{
	385	cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
	386	return;
	387	}
	388	Float_t hit[3];
029912b7	389
	390
	391	//Remove tracks which are no good:
	392	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	393	{
	394	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	395	if(!track) continue;
	396	if(track->GetFirstPointX() > fTransform->Row2X(NRows[fPatch][1]) \|\| track->GetPt()<0.1)
	397	fTracks->Remove(i);
	398	}
	399
	400	Int_t sector,row;
95a00d93	401	for(Int_t i=NRows[fPatch][1]; i>=NRows[fPatch][0]; i--)
735e167e	402	{
	403	for(Int_t j=0; j<fTracks->GetNTracks(); j++)
	404	{
	405	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(j);
	406	if(!track) continue;
029912b7	407
	408	if(track->GetNHits() < fTrackThreshold)
	409	{
	410	fTracks->Remove(j);
	411	continue;
	412	}
	413
95a00d93	414	if(!track->GetCrossingPoint(i,hit))
735e167e	415	{
029912b7	416	cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
	417	"First point "<<track->GetFirstPointX()<<
	418	" nhits "<<track->GetNHits()<<endl;//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
	419	//"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
	420	//"--------"<<endl;
	421	fTracks->Remove(j);
735e167e	422	continue;
735e167e	423	}
95a00d93	424	//cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
95a00d93	425
029912b7	426	fTransform->Slice2Sector(fSlice,i,sector,row);
029912b7	427	fTransform->Local2Raw(hit,sector,row);
95a00d93	428	if(hit[1]<0 \|\| hit[1]>fTransform->GetNPads(i) \|\|
	429	hit[2]<0 \|\| hit[2]>fTransform->GetNTimeBins())
	430	{//Track is leaving the patch, so flag the track hits (<0)
	431	track->SetPadHit(i,-1);
	432	track->SetTimeHit(i,-1);
	433	continue;
	434	}
	435
735e167e	436	track->SetPadHit(i,hit[1]);
735e167e	437	track->SetTimeHit(i,hit[2]);
95a00d93	438
735e167e	439	//if(hit[1]<0 \|\| hit[2]>445)
95a00d93	440	//if(hit[2]<0 \|\| hit[2]>445)
95a00d93	441	//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	442	//cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
	443	}
	444	}
	445	fTracks->Compress();
029912b7	446	cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
735e167e	447	}
	448
	449	void AliL3Modeller::CheckForOverlaps()
	450	{
	451	//Flag the tracks that overlap
	452
029912b7	453	cout<<"Checking for overlaps...";
735e167e	454
	455	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	456	{
	457	AliL3ModelTrack track1 = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	458	if(!track1) continue;
	459	for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
	460	{
	461	AliL3ModelTrack track2 = (AliL3ModelTrack)fTracks->GetCheckedTrack(j);
	462	if(!track2) continue;
029912b7	463	for(Int_t k=NRows[fPatch][0]; k<=NRows[fPatch][1]; k++)
735e167e	464	{
95a00d93	465	if(track1->GetPadHit(k)<0 \|\| track1->GetTimeHit(k)<0 \|\|
	466	track2->GetPadHit(k)<0 \|\| track2->GetTimeHit(k)<0)
	467	continue;
029912b7	468	if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) <= fPadOverlap &&
029912b7	469	fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) <= fTimeOverlap)
735e167e	470	{
029912b7	471	track2->SetOverlap(k,i);
029912b7	472	track1->SetOverlap(k,j);
735e167e	473	}
	474	}
	475	}
	476	}
029912b7	477	cout<<"done"<<endl;
735e167e	478	}
	479
	480
	481	void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
	482	{
	483
	484	Float_t padw,timew;
	485	if(fPatch < 3)
	486	padw = fTransform->GetPadPitchWidthLow();
	487	else
	488	padw = fTransform->GetPadPitchWidthUp();
	489	Float_t charge = (Float_t)cl->fCharge;
	490	Float_t pad = (Float_t)cl->fPad/charge;
	491	Float_t time = (Float_t)cl->fTime/charge;
	492	Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
	493	sigmaY2 = (s2 + 1./12)padwpadw;
	494
	495	if(s2 != 0)
	496	{
	497	sigmaY2 = sigmaY2*0.108;
	498	if(fPatch<3)
	499	sigmaY2 = sigmaY2*2.07;
	500	}
	501
	502	s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
	503	timew = fTransform->GetZWidth();
	504	sigmaZ2 = (s2 +1./12)timewtimew;
	505	if(s2 != 0)
	506	{
	507	sigmaZ2 = sigmaZ2*0.169;
	508	if(fPatch < 3)
	509	sigmaZ2 = sigmaZ2*1.77;
	510	}
	511
	512	}