[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,"%s/digits_%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
	  
	  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) 
		    {pad=hitpad+1; padsign=1; 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<<" pad "<<dpad<<endl;
		      track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
		      break;
		    } 
		}
	      // pad += padsign;
	    }
	}
      fMemHandler->UpdateRowPointer(rowPt);
      
    }
  delete [] row;
  
}

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][0]; i<NRows[fPatch][1]; i++)
    {
      for(Int_t j=0; j<fTracks->GetNTracks(); j++)
	{
	  AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
	  if(!track) continue;
	  //if(!track->GetCrossingPoint(i,hit)) 
	  // fTracks->Remove(j); //Track is bending out.
	  track->CalculatePoint(fTransform->Row2X(i));
	  if(!track->IsPoint())
	    {
	      fTracks->Remove(j);
	      continue;
	    }
	  hit[1]=track->GetPointY();
	  hit[2]=track->GetPointZ();
	  fTransform->Local2Raw(hit,fSlice,i);
	  track->SetPadHit(i,hit[1]);
	  track->SetTimeHit(i,hit[2]);
	  //if(hit[1]<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(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");
	48
	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];
74	sprintf(fname,"%s/digits_%d_%d.raw",path,fSlice,fPatch);
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
130
131	Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
132	Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
133	cout<<"Checking track with pad "<<hitpad<<" time "<<hittime<<endl;
134	pad = hitpad;
135	time = hittime;
136	Int_t padsign=-1;
137	Int_t timesign=-1;
138
139	memset(&cluster,0,sizeof(Cluster));
140
141	while(1)//Process this padrow
142	{
143	seq_charge=0;
144	timesign=-1;
145	time = hittime;
146	while(1) //Process sequence on this pad:
147	{
148	charge = row[ntimes*pad+time].fCharge;
149	if(charge==0 && timesign==-1)
150	{time=hittime+1; timesign=1; continue;}
151	else if(charge==0 && timesign==1)
152	break;
153	//cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
154
155	seq_charge += charge;
156
157	cluster.fTime += time*charge;
158	cluster.fPad += pad*charge;
159	cluster.fCharge += charge;
160	cluster.fSigmaY2 += padpadcharge;
161	cluster.fSigmaZ2 += timetimecharge;
162
163	row[ntimes*pad+time].fUsed = kTRUE;
164	time += timesign;
165	}
166	cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
167	if(seq_charge)
168	pad += padsign;
169
170	else //Nothing more on this pad, goto next pad
171	{
172	if(padsign==-1)
173	{pad=hitpad+1; padsign=1; continue;}
174	else //Nothing more in this cluster
175	{
176	Float_t fcharge = (Float_t)cluster.fCharge;
177	Float_t fpad = ((Float_t)cluster.fPad/fcharge);
178	Float_t ftime = ((Float_t)cluster.fTime/fcharge);
179	Float_t sigmaY2,sigmaZ2;
180	CalcClusterWidth(&cluster,sigmaY2,sigmaZ2);
181	//cout<<"row "<<i<<" pad "<<dpad<<endl;
182	track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
183	break;
184	}
185	}
186	// pad += padsign;
187	}
188	}
189	fMemHandler->UpdateRowPointer(rowPt);
190
191	}
192	delete [] row;
193
194	}
195
196	void AliL3Modeller::CalculateCrossingPoints()
197	{
198	cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
199	if(!fTracks)
200	{
201	cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
202	return;
203	}
204	Float_t hit[3];
205	for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
206	{
207	for(Int_t j=0; j<fTracks->GetNTracks(); j++)
208	{
209	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(j);
210	if(!track) continue;
211	//if(!track->GetCrossingPoint(i,hit))
212	// fTracks->Remove(j); //Track is bending out.
213	track->CalculatePoint(fTransform->Row2X(i));
214	if(!track->IsPoint())
215	{
216	fTracks->Remove(j);
217	continue;
218	}
219	hit[1]=track->GetPointY();
220	hit[2]=track->GetPointZ();
221	fTransform->Local2Raw(hit,fSlice,i);
222	track->SetPadHit(i,hit[1]);
223	track->SetTimeHit(i,hit[2]);
224	//if(hit[1]<0 \|\| hit[2]>445)
225	cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
226	//cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
227	}
228	}
229	fTracks->Compress();
230	cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
231	}
232
233	void AliL3Modeller::CheckForOverlaps()
234	{
235	//Flag the tracks that overlap
236
237	cout<<"Checking for overlaps"<<endl;
238
239	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
240	{
241	AliL3ModelTrack track1 = (AliL3ModelTrack)fTracks->GetCheckedTrack(i);
242	if(!track1) continue;
243	for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
244	{
245	AliL3ModelTrack track2 = (AliL3ModelTrack)fTracks->GetCheckedTrack(j);
246	if(!track2) continue;
247	for(Int_t k=NRows[fPatch][0]; k<NRows[fPatch][1]; k++)
248	{
249	if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) < fPadOverlap &&
250	fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) < fTimeOverlap)
251	{
252	track1->SetOverlap(j);
253	track2->SetOverlap(i);
254	}
255	}
256	}
257	}
258
259	}
260
261
262	void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
263	{
264
265	Float_t padw,timew;
266	if(fPatch < 3)
267	padw = fTransform->GetPadPitchWidthLow();
268	else
269	padw = fTransform->GetPadPitchWidthUp();
270	Float_t charge = (Float_t)cl->fCharge;
271	Float_t pad = (Float_t)cl->fPad/charge;
272	Float_t time = (Float_t)cl->fTime/charge;
273	Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
274	sigmaY2 = (s2 + 1./12)padwpadw;
275
276	if(s2 != 0)
277	{
278	sigmaY2 = sigmaY2*0.108;
279	if(fPatch<3)
280	sigmaY2 = sigmaY2*2.07;
281	}
282
283	s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
284	timew = fTransform->GetZWidth();
285	sigmaZ2 = (s2 +1./12)timewtimew;
286	if(s2 != 0)
287	{
288	sigmaZ2 = sigmaZ2*0.169;
289	if(fPatch < 3)
290	sigmaZ2 = sigmaZ2*1.77;
291	}
292
293	}