Data structures for track and clusters
[u/mrichter/AliRoot.git] / HLT / comp / AliL3Modeller.cxx
CommitLineData
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
19ClassImp(AliL3Modeller)
20
21AliL3Modeller::AliL3Modeller()
22{
23 fMemHandler=0;
24 fTracks=0;
25 fTransform=0;
26}
27
28
29AliL3Modeller::~AliL3Modeller()
30{
31 if(fMemHandler)
32 delete fMemHandler;
33 if(fTracks)
34 delete fTracks;
35 if(fTransform)
36 delete fTransform;
37
38}
39
40void 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
86void 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,ntimes*npads*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 += pad*pad*charge;
166 cluster.fSigmaZ2 += time*time*charge;
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)
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 216void 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
231void 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 = digitcount*sizeof(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 299void 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;
320 }
95a00d93 321 //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
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]);
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)
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
345void 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
377void 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)*padw*padw;
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)*timew*timew;
401 if(s2 != 0)
402 {
403 sigmaZ2 = sigmaZ2*0.169;
404 if(fPatch < 3)
405 sigmaZ2 = sigmaZ2*1.77;
406 }
407
408}