Including latest changes
[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");
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,"%s/digits_%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
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
130
131 Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
132 Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
133 cout<<"Checking track with pad "<<hitpad<<" time "<<hittime<<endl;
134 pad = hitpad;
135 time = hittime;
136 Int_t padsign=-1;
137 Int_t timesign=-1;
138
139 memset(&cluster,0,sizeof(Cluster));
140
141 while(1)//Process this padrow
142 {
143 seq_charge=0;
144 timesign=-1;
145 time = hittime;
146 while(1) //Process sequence on this pad:
147 {
148 charge = row[ntimes*pad+time].fCharge;
149 if(charge==0 && timesign==-1)
150 {time=hittime+1; timesign=1; continue;}
151 else if(charge==0 && timesign==1)
152 break;
153 //cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
154
155 seq_charge += charge;
156
157 cluster.fTime += time*charge;
158 cluster.fPad += pad*charge;
159 cluster.fCharge += charge;
160 cluster.fSigmaY2 += pad*pad*charge;
161 cluster.fSigmaZ2 += time*time*charge;
162
163 row[ntimes*pad+time].fUsed = kTRUE;
164 time += timesign;
165 }
166 cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
167 if(seq_charge)
168 pad += padsign;
169
170 else //Nothing more on this pad, goto next pad
171 {
172 if(padsign==-1)
173 {pad=hitpad+1; padsign=1; continue;}
174 else //Nothing more in this cluster
175 {
176 Float_t fcharge = (Float_t)cluster.fCharge;
177 Float_t fpad = ((Float_t)cluster.fPad/fcharge);
178 Float_t ftime = ((Float_t)cluster.fTime/fcharge);
179 Float_t sigmaY2,sigmaZ2;
180 CalcClusterWidth(&cluster,sigmaY2,sigmaZ2);
181 //cout<<"row "<<i<<" pad "<<dpad<<endl;
182 track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
183 break;
184 }
185 }
186 // pad += padsign;
187 }
188 }
189 fMemHandler->UpdateRowPointer(rowPt);
190
191 }
192 delete [] row;
193
194}
195
196void AliL3Modeller::CalculateCrossingPoints()
197{
198 cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
199 if(!fTracks)
200 {
201 cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
202 return;
203 }
204 Float_t hit[3];
205 for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
206 {
207 for(Int_t j=0; j<fTracks->GetNTracks(); j++)
208 {
209 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
210 if(!track) continue;
211 //if(!track->GetCrossingPoint(i,hit))
212 // fTracks->Remove(j); //Track is bending out.
213 track->CalculatePoint(fTransform->Row2X(i));
214 if(!track->IsPoint())
215 {
216 fTracks->Remove(j);
217 continue;
218 }
219 hit[1]=track->GetPointY();
220 hit[2]=track->GetPointZ();
221 fTransform->Local2Raw(hit,fSlice,i);
222 track->SetPadHit(i,hit[1]);
223 track->SetTimeHit(i,hit[2]);
224 //if(hit[1]<0 || hit[2]>445)
225 cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
226 //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
227 }
228 }
229 fTracks->Compress();
230 cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
231}
232
233void AliL3Modeller::CheckForOverlaps()
234{
235 //Flag the tracks that overlap
236
237 cout<<"Checking for overlaps"<<endl;
238
239 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
240 {
241 AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
242 if(!track1) continue;
243 for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
244 {
245 AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
246 if(!track2) continue;
247 for(Int_t k=NRows[fPatch][0]; k<NRows[fPatch][1]; k++)
248 {
249 if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) < fPadOverlap &&
250 fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) < fTimeOverlap)
251 {
252 track1->SetOverlap(j);
253 track2->SetOverlap(i);
254 }
255 }
256 }
257 }
258
259}
260
261
262void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
263{
264
265 Float_t padw,timew;
266 if(fPatch < 3)
267 padw = fTransform->GetPadPitchWidthLow();
268 else
269 padw = fTransform->GetPadPitchWidthUp();
270 Float_t charge = (Float_t)cl->fCharge;
271 Float_t pad = (Float_t)cl->fPad/charge;
272 Float_t time = (Float_t)cl->fTime/charge;
273 Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
274 sigmaY2 = (s2 + 1./12)*padw*padw;
275
276 if(s2 != 0)
277 {
278 sigmaY2 = sigmaY2*0.108;
279 if(fPatch<3)
280 sigmaY2 = sigmaY2*2.07;
281 }
282
283 s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
284 timew = fTransform->GetZWidth();
285 sigmaZ2 = (s2 +1./12)*timew*timew;
286 if(s2 != 0)
287 {
288 sigmaZ2 = sigmaZ2*0.169;
289 if(fPatch < 3)
290 sigmaZ2 = sigmaZ2*1.77;
291 }
292
293}