New topdir Makefile for compiling all libraries in the HLT tree.
[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
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
27ClassImp(AliL3Modeller)
28
29AliL3Modeller::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
40AliL3Modeller::~AliL3Modeller()
41{
42 if(fMemHandler)
43 delete fMemHandler;
44 if(fTracks)
45 delete fTracks;
735e167e 46}
47
4a838220 48void 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;
53 fTimeOverlap=8;
4a838220 54
be6ddb10 55 sprintf(fPath,"%s",path);
029912b7 56
4a838220 57 AliL3Transform::Init(fPath);
58
735e167e 59 fTracks = new AliL3TrackArray("AliL3ModelTrack");
92a876e2 60
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
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 112void 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;
130 Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
029912b7 131 Int_t bounds = ntimes*npads;
132 Digit *row = new Digit[bounds];
735e167e 133
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,ntimes*npads*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;
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
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;
183 }
735e167e 184 seq_charge=0;
185 timesign=-1;
186 time = hittime;
4a838220 187
735e167e 188 while(1) //Process sequence on this pad:
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 }
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 += pad*pad*charge;
243 cluster.fSigmaZ2 += time*time*charge;
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
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;
283 }
95a00d93 284 }
735e167e 285 else //Nothing more in this cluster
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 }
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))
309 cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<AliL3Transform::GetNRows(fPatch)<<endl<<endl;
029912b7 310 }
311
312}
313
4a838220 314void 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 330void 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 345void 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))];
353 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
95a00d93 354 {
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 = digitcount*sizeof(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++;
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);
408 mem->SetBinaryOutput(fname);
4a838220 409 mem->Memory2CompBinary((UInt_t)AliL3Transform::GetNRows(fPatch),(AliL3DigitRowData*)data);
95a00d93 410 mem->CloseBinaryOutput();
411 delete mem;
029912b7 412 delete [] data;
95a00d93 413}
414
415
735e167e 416void 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;
458 }
95a00d93 459 //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
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
471
735e167e 472 track->SetPadHit(i,hit[1]);
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)
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
485void 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 &&
505 fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) <= fTimeOverlap)
735e167e 506 {
029912b7 507 track2->SetOverlap(k,i);
508 track1->SetOverlap(k,j);
735e167e 509 }
510 }
511 }
512 }
029912b7 513 cout<<"done"<<endl;
735e167e 514}
515
516
517void 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);//*padw*padw;
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 */
541
735e167e 542 s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
4a838220 543 timew = AliL3Transform::GetZWidth();
544 sigmaZ2 = (s2);// +1./12);//*timew*timew;
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}