[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"

//_____________________________________________________________
// 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;
  fPadOverlap=0;
  fTimeOverlap=0;
  fRowData=0;
}


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)
{
  fSlice = slice;
  fPatch = patch;
  fPadOverlap=6;
  fTimeOverlap=8;
  
  sprintf(fPath,"%s",path);
  
  AliL3Transform::Init(fPath);
  
  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)
  if(!file->SetBinaryInput(fname))
    {
      cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
      return;
    }
  file->Binary2TrackArray(fTracks);
  file->CloseBinaryInput();
  delete file;
  
  if(houghtracks)
    cout<<"AliL3Modeller is assuming local hough tracksegments!"<<endl;
  else
    cout<<"AliL3Modeller is assuming global tracks!"<<endl;

  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();
  
  if(!houghtracks)
    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 = 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++)
    {
      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,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));
	  
	  Int_t npads=0;
	  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) 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;
		    }
		  
		  /*
		  if(row[index].fUsed == kTRUE)//Only use the digits once....
		    charge = 0;
		  else
		    charge = row[index].fCharge;
		  */
		  
		  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;
		}
	      
	      
	      if(seq_charge)//There was something on this pad.
		{
		  pad += padsign;
		  npads++;
		}
	      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 < AliL3Transform::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,npads);
			  break;
			}
		    }
		  else //Nothing more in this cluster
		    {
		      //cout<<"Filling final cluster"<<endl;
		      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;
    }
  
}

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