Data structures for track and clusters
[u/mrichter/AliRoot.git] / HLT / comp / AliL3Modeller.cxx
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,"%sdigits_%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,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
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             }
135           
136           Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
137           Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
138           //cout<<"Checking track with pad "<<hitpad<<" time "<<hittime<<endl;
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                 }
171               //cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
172               if(seq_charge)
173                 pad += padsign;
174               else //Nothing more on this pad, goto next pad
175                 {
176                   if(padsign==-1) 
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                     }
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);
201                       //cout<<"row "<<i<<" charge "<<fcharge<<endl;
202                       track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
203                       break;
204                     } 
205                 }
206               // pad += padsign;
207             }
208         }
209       FillZeros(rowPt,row);
210       fMemHandler->UpdateRowPointer(rowPt);
211     }
212   delete [] row;
213   
214 }
215
216 void 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
231 void 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
299 void 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];
308   for(Int_t i=NRows[fPatch][1]; i>=NRows[fPatch][0]; i--)
309     {
310       for(Int_t j=0; j<fTracks->GetNTracks(); j++)
311         {
312           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
313           if(!track) continue;
314           if(track->GetNHits() < 100)
315             fTracks->Remove(j);
316           if(!track->GetCrossingPoint(i,hit)) 
317             {
318               cerr<<"AliL3Modeller::CalculateCrossingPoints : Track does not intersect line "<<endl;
319               continue;
320             }
321           //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
322           
323           fTransform->Local2Raw(hit,fSlice,i);
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             
332           track->SetPadHit(i,hit[1]);
333           track->SetTimeHit(i,hit[2]);
334           
335           //if(hit[1]<0 || hit[2]>445)
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;
338           //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
339         }
340     }
341   fTracks->Compress();
342   //cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
343 }
344
345 void 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             {
361               if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
362                  track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
363                 continue;
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
377 void 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 }