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