[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 *path)
{
  fSlice = slice;
  fPatch = patch;
  fPadOverlap=4;
  fTimeOverlap=4;
  fTransform = new AliL3Transform();
  fTracks = new AliL3TrackArray("AliL3ModelTrack");

  AliL3MemHandler *file = new AliL3MemHandler();
  if(!file->SetBinaryInput("tracks.raw"))
    {
      cerr<<"AliL3Modeller::Init : Error opening trackfile"<<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();
  Char_t fname[100];
  sprintf(fname,"%sdigits_%d_%d.raw",path,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::Process()
{
  
  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.
  Digit *row = new Digit[(ntimes)*(npads)];
  
  Int_t seq_charge;
  Int_t pad,time;
  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;
	}
      
      for(Int_t k=0; k<fTracks->GetNTracks(); k++)
	{
	  AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
	  if(!track) continue;
	  if(track->GetOverlap()>=0) continue;//Track is overlapping
	  if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0)
	    {
	      track->SetCluster(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 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
	    {
	      seq_charge=0;
	      timesign=-1;
	      time = hittime;
	      while(1) //Process sequence on this pad:
		{
		  charge = row[ntimes*pad+time].fCharge;
		  if(charge==0 && timesign==-1)
		    {time=hittime+1; timesign=1; continue;}
		  else if(charge==0 && timesign==1)
		    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--; //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 && cluster.fCharge==0 && abs(pad-hitpad) < fPadOverlap/2)
		    {
		      pad++;
		      continue;
		    }
		  else //Nothing more in this cluster
		    {
		      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);
		      //cout<<"row "<<i<<" charge "<<fcharge<<endl;
		      track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
		      break;
		    } 
		}
	      // pad += padsign;
	    }
	}
      FillZeros(rowPt,row);
      fMemHandler->UpdateRowPointer(rowPt);
    }
  delete [] row;
  
}

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(Char_t *output)
{
  //Write remaining (nonzero) digits to file.
  
  cout<<"Writing remaining data to file "<<output<<endl;
  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;
    }

  AliL3MemHandler *mem = new AliL3MemHandler();
  mem->SetBinaryOutput(output);
  mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
  mem->CloseBinaryOutput();
  delete mem;
}


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];
  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() < 100)
	    fTracks->Remove(j);
	  if(!track->GetCrossingPoint(i,hit)) 
	    {
	      cerr<<"AliL3Modeller::CalculateCrossingPoints : Track does not intersect line "<<endl;
	      continue;
	    }
	  //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
	  
	  fTransform->Local2Raw(hit,fSlice,i);
	  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"<<endl;
  
  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)
		{
		  track1->SetOverlap(j);
		  track2->SetOverlap(i);
		}
	    }
	}
    }
  
}


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$
	2
	3	// Author: Anders Vestbo <mailto:vestbo$fi.uib.no>
	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
	40	void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *path)
	41	{
	42	fSlice = slice;
	43	fPatch = patch;
	44	fPadOverlap=4;
	45	fTimeOverlap=4;
	46	fTransform = new AliL3Transform();
	47	fTracks = new AliL3TrackArray("AliL3ModelTrack");
95a00d93	48
735e167e	49	AliL3MemHandler *file = new AliL3MemHandler();
	50	if(!file->SetBinaryInput("tracks.raw"))
	51	{
	52	cerr<<"AliL3Modeller::Init : Error opening trackfile"<<endl;
	53	return;
	54	}
	55	file->Binary2TrackArray(fTracks);
	56	file->CloseBinaryInput();
	57	delete file;
	58
	59	fTracks->QSort();
	60	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	61	{
	62	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	63	if(!track) continue;
	64	track->Init(fSlice,fPatch);
	65	track->Rotate(fSlice,kTRUE);
	66	track->CalculateHelix();
	67	}
	68
	69	CalculateCrossingPoints();
	70	CheckForOverlaps();
	71
	72	fMemHandler = new AliL3MemHandler();
	73	Char_t fname[100];
95a00d93	74	sprintf(fname,"%sdigits_%d_%d.raw",path,fSlice,fPatch);
735e167e	75	if(!fMemHandler->SetBinaryInput(fname))
	76	{
	77	cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
	78	return;
	79	}
	80	UInt_t ndigits;
	81	AliL3DigitRowData digits=(AliL3DigitRowData)fMemHandler->CompBinary2Memory(ndigits);
	82
	83	SetInputData(digits);
	84	}
	85
	86	void AliL3Modeller::Process()
	87	{
	88
	89	if(!fTracks)
	90	{
	91	cerr<<"AliL3Modeller::Process : No tracks"<<endl;
	92	return;
	93	}
	94	if(!fRowData)
	95	{
	96	cerr<<"AliL3Modeller::Process : No data "<<endl;
	97	return;
	98	}
	99
	100	AliL3DigitRowData *rowPt = fRowData;
	101	AliL3DigitData *digPt=0;
	102
	103	Int_t ntimes = fTransform->GetNTimeBins()+1;
	104	Int_t npads = fTransform->GetNPads(NRows[fPatch][1])+1;//Max num of pads.
	105	Digit row = new Digit[(ntimes)(npads)];
	106
	107	Int_t seq_charge;
	108	Int_t pad,time;
	109	Short_t charge;
	110	Cluster cluster;
	111
	112	for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
	113	{
	114	memset((void)row,0,ntimesnpads*sizeof(Digit));
	115	digPt = (AliL3DigitData*)rowPt->fDigitData;
	116	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	117	{
	118	pad = digPt[j].fPad;
	119	time = digPt[j].fTime;
	120	charge = digPt[j].fCharge;
	121	row[ntimes*pad+time].fCharge = charge;
	122	row[ntimes*pad+time].fUsed = kFALSE;
	123	}
	124
	125	for(Int_t k=0; k<fTracks->GetNTracks(); k++)
	126	{
	127	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(k);
	128	if(!track) continue;
	129	if(track->GetOverlap()>=0) continue;//Track is overlapping
95a00d93	130	if(track->GetPadHit(i)<0 \|\| track->GetTimeHit(i)<0)
	131	{
	132	track->SetCluster(0,0,0,0,0); //The track has left the patch.
	133	continue;
	134	}
735e167e	135
	136	Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
	137	Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
95a00d93	138	//cout<<"Checking track with pad "<<hitpad<<" time "<<hittime<<endl;
735e167e	139	pad = hitpad;
	140	time = hittime;
	141	Int_t padsign=-1;
	142	Int_t timesign=-1;
	143
	144	memset(&cluster,0,sizeof(Cluster));
	145
	146	while(1)//Process this padrow
	147	{
	148	seq_charge=0;
	149	timesign=-1;
	150	time = hittime;
	151	while(1) //Process sequence on this pad:
	152	{
	153	charge = row[ntimes*pad+time].fCharge;
	154	if(charge==0 && timesign==-1)
	155	{time=hittime+1; timesign=1; continue;}
	156	else if(charge==0 && timesign==1)
	157	break;
	158	//cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
	159
	160	seq_charge += charge;
	161
	162	cluster.fTime += time*charge;
	163	cluster.fPad += pad*charge;
	164	cluster.fCharge += charge;
	165	cluster.fSigmaY2 += padpadcharge;
	166	cluster.fSigmaZ2 += timetimecharge;
	167
	168	row[ntimes*pad+time].fUsed = kTRUE;
	169	time += timesign;
	170	}
95a00d93	171	//cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
735e167e	172	if(seq_charge)
735e167e	173	pad += padsign;
735e167e	174	else //Nothing more on this pad, goto next pad
	175	{
	176	if(padsign==-1)
95a00d93	177	{
	178	if(cluster.fCharge==0 && abs(pad-hitpad) < fPadOverlap/2)
	179	{
	180	pad--; //In this case, we haven't found anything yet,
	181	} //so we will try to expand our search within the natural boundaries.
	182	else
	183	{
	184	pad=hitpad+1;
	185	padsign=1;
	186	}
	187	continue;
	188	}
	189	else if(padsign==1 && cluster.fCharge==0 && abs(pad-hitpad) < fPadOverlap/2)
	190	{
	191	pad++;
	192	continue;
	193	}
735e167e	194	else //Nothing more in this cluster
	195	{
	196	Float_t fcharge = (Float_t)cluster.fCharge;
	197	Float_t fpad = ((Float_t)cluster.fPad/fcharge);
	198	Float_t ftime = ((Float_t)cluster.fTime/fcharge);
	199	Float_t sigmaY2,sigmaZ2;
	200	CalcClusterWidth(&cluster,sigmaY2,sigmaZ2);
95a00d93	201	//cout<<"row "<<i<<" charge "<<fcharge<<endl;
735e167e	202	track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
	203	break;
	204	}
	205	}
	206	// pad += padsign;
	207	}
	208	}
95a00d93	209	FillZeros(rowPt,row);
735e167e	210	fMemHandler->UpdateRowPointer(rowPt);
735e167e	211	}
	212	delete [] row;
	213
	214	}
	215
95a00d93	216	void AliL3Modeller::FillZeros(AliL3DigitRowData rowPt,Digit row)
	217	{
	218	//Fill zero where data has been used.
	219
	220	Int_t ntimes = fTransform->GetNTimeBins()+1;
	221	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	222	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	223	{
	224	Int_t pad = digPt[j].fPad;
	225	Int_t time = digPt[j].fTime;
	226	if(row[ntimes*pad+time].fUsed==kTRUE)
	227	digPt[j].fCharge = 0;
	228	}
	229	}
	230
	231	void AliL3Modeller::WriteRemaining(Char_t *output)
	232	{
	233	//Write remaining (nonzero) digits to file.
	234
	235	cout<<"Writing remaining data to file "<<output<<endl;
	236	AliL3DigitRowData *rowPt;
	237	rowPt = (AliL3DigitRowData*)fRowData;
	238	Int_t digitcount=0;
	239	Int_t ndigits[(NumRows[fPatch])];
	240	for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
	241	{
	242	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	243	ndigits[(i-NRows[fPatch][0])]=0;
	244	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	245	{
	246	if(digPt[j].fCharge==0) continue;
	247	digitcount++;
	248	ndigits[(i-NRows[fPatch][0])]++;
	249	}
	250	//cout<<"Difference "<<(int)ndigits[(i-NRows[fPatch][0])]<<" "<<(int)rowPt->fNDigit<<endl;
	251	fMemHandler->UpdateRowPointer(rowPt);
	252	}
	253
	254	Int_t size = digitcountsizeof(AliL3DigitData) + NumRows[fPatch]sizeof(AliL3DigitRowData);
	255	Byte_t *data = new Byte_t[size];
	256	memset(data,0,size);
	257	AliL3DigitRowData tempPt = (AliL3DigitRowData)data;
	258	rowPt = (AliL3DigitRowData*)fRowData;
	259
	260	for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
	261	{
	262	Int_t localcount=0;
	263	tempPt->fRow = i;
	264	tempPt->fNDigit = ndigits[(i-NRows[fPatch][0])];
	265	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	266	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	267	{
	268	if(digPt[j].fCharge==0) continue;
	269	if(localcount >= ndigits[(i-NRows[fPatch][0])])
	270	{
	271	cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
	272	return;
	273	}
	274	tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
	275	tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
	276	tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
	277	localcount++;
	278	}
	279	if(ndigits[(i-NRows[fPatch][0])] != localcount)
280	{
281	cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
282	return;
283	}
284	fMemHandler->UpdateRowPointer(rowPt);
285	Byte_t tmp = (Byte_t)tempPt;
286	Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-NRows[fPatch][0])]*sizeof(AliL3DigitData);
287	tmp += size;
288	tempPt = (AliL3DigitRowData*)tmp;
289	}
290
291	AliL3MemHandler *mem = new AliL3MemHandler();
292	mem->SetBinaryOutput(output);
293	mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
294	mem->CloseBinaryOutput();
295	delete mem;
296	}
297
298
735e167e	299	void AliL3Modeller::CalculateCrossingPoints()
	300	{
	301	cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
	302	if(!fTracks)
	303	{
	304	cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
	305	return;
	306	}
	307	Float_t hit[3];
95a00d93	308	for(Int_t i=NRows[fPatch][1]; i>=NRows[fPatch][0]; i--)
735e167e	309	{
	310	for(Int_t j=0; j<fTracks->GetNTracks(); j++)
	311	{
	312	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(j);
	313	if(!track) continue;
95a00d93	314	if(track->GetNHits() < 100)
	315	fTracks->Remove(j);
	316	if(!track->GetCrossingPoint(i,hit))
735e167e	317	{
95a00d93	318	cerr<<"AliL3Modeller::CalculateCrossingPoints : Track does not intersect line "<<endl;
735e167e	319	continue;
735e167e	320	}
95a00d93	321	//cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
95a00d93	322
735e167e	323	fTransform->Local2Raw(hit,fSlice,i);
95a00d93	324	if(hit[1]<0 \|\| hit[1]>fTransform->GetNPads(i) \|\|
	325	hit[2]<0 \|\| hit[2]>fTransform->GetNTimeBins())
	326	{//Track is leaving the patch, so flag the track hits (<0)
	327	track->SetPadHit(i,-1);
	328	track->SetTimeHit(i,-1);
	329	continue;
	330	}
	331
735e167e	332	track->SetPadHit(i,hit[1]);
735e167e	333	track->SetTimeHit(i,hit[2]);
95a00d93	334
735e167e	335	//if(hit[1]<0 \|\| hit[2]>445)
95a00d93	336	//if(hit[2]<0 \|\| hit[2]>445)
95a00d93	337	//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	338	//cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
	339	}
	340	}
	341	fTracks->Compress();
95a00d93	342	//cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
735e167e	343	}
	344
	345	void AliL3Modeller::CheckForOverlaps()
	346	{
	347	//Flag the tracks that overlap
	348
	349	cout<<"Checking for overlaps"<<endl;
	350
	351	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	352	{
	353	AliL3ModelTrack track1 = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
	354	if(!track1) continue;
	355	for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
	356	{
	357	AliL3ModelTrack track2 = (AliL3ModelTrack)fTracks->GetCheckedTrack(j);
	358	if(!track2) continue;
	359	for(Int_t k=NRows[fPatch][0]; k<NRows[fPatch][1]; k++)
	360	{
95a00d93	361	if(track1->GetPadHit(k)<0 \|\| track1->GetTimeHit(k)<0 \|\|
	362	track2->GetPadHit(k)<0 \|\| track2->GetTimeHit(k)<0)
	363	continue;
735e167e	364	if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) < fPadOverlap &&
	365	fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) < fTimeOverlap)
	366	{
	367	track1->SetOverlap(j);
	368	track2->SetOverlap(i);
	369	}
	370	}
	371	}
	372	}
	373
	374	}
	375
	376
	377	void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
	378	{
	379
	380	Float_t padw,timew;
	381	if(fPatch < 3)
	382	padw = fTransform->GetPadPitchWidthLow();
	383	else
	384	padw = fTransform->GetPadPitchWidthUp();
	385	Float_t charge = (Float_t)cl->fCharge;
	386	Float_t pad = (Float_t)cl->fPad/charge;
	387	Float_t time = (Float_t)cl->fTime/charge;
	388	Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
	389	sigmaY2 = (s2 + 1./12)padwpadw;
	390
	391	if(s2 != 0)
	392	{
	393	sigmaY2 = sigmaY2*0.108;
	394	if(fPatch<3)
	395	sigmaY2 = sigmaY2*2.07;
	396	}
	397
	398	s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
	399	timew = fTransform->GetZWidth();
	400	sigmaZ2 = (s2 +1./12)timewtimew;
	401	if(s2 != 0)
	402	{
	403	sigmaZ2 = sigmaZ2*0.169;
	404	if(fPatch < 3)
	405	sigmaZ2 = sigmaZ2*1.77;
	406	}
	407
	408	}