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