]>
Commit | Line | Data |
---|---|---|
dfd03fc3 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
dfd03fc3 | 16 | /////////////////////////////////////////////////////////////////////////////// |
17 | // // | |
18 | // TRD MCM (Multi Chip Module) simulator // | |
19 | // // | |
20 | /////////////////////////////////////////////////////////////////////////////// | |
21 | ||
ecf39416 | 22 | /* $Id$ */ |
23 | ||
24 | /* | |
25 | ||
26 | New release on 2007/08/17 | |
27 | ||
28 | AliTRDmcmSim is now stably working and zero suppression function seems ok. | |
29 | From now, the default version of raw data is set to 3 in AliTRDfeeParam. | |
30 | ||
31 | The following internal parameters were abolished because it is useless and | |
32 | made trouble: | |
33 | ||
34 | fColOfADCbeg | |
35 | fColOfADCend | |
36 | ||
37 | GetCol member was modified accordingly. | |
38 | ||
39 | New member function DumpData was prepared for diagnostics. | |
40 | ||
41 | ZSMapping member function was debugged. It was causing crash due to | |
42 | wrong indexing in 1 dimensional numbering. Also code was shaped up better. | |
43 | ||
44 | */ | |
45 | ||
46 | /*Semi-final version of TRD raw data simulation code with zero suppression (ZS) | |
47 | similar to TRD FEE. ZS is realized by the class group: | |
48 | ||
49 | AliTRDfeeParam | |
50 | AliTRDmcmSim | |
51 | AliTRDrawData | |
52 | ||
53 | AliTRDfeeParam has been modified to have more parameters like raw data | |
54 | production version and so on. AliTRDmcmSim is new class and this is the core | |
55 | of MCM (PASA+TRAP) simulator. It has still very simple function and it will be | |
56 | another project to improve this to make it closer to the reall FEE. | |
57 | AliTRDrawData has been modified to use new class AliTRDmcmSim. | |
dfd03fc3 | 58 | |
ecf39416 | 59 | These modifications were tested on Aug. 02 HEAD version that code itself |
60 | compiles. I'm sure there must be still bugs and we need testing by as many as | |
61 | possible persons now. Especially it seems HLT part is impacted by problems | |
62 | because some parameters were moved from AliTRDrawData to AliTRDfeeParam (like | |
63 | fRawVersion disappeared from AliTRDrawData). | |
64 | ||
65 | In TRD definition, we have now 4 raw data versions. | |
66 | ||
67 | 0 very old offline version (by Bogdan) | |
68 | 1 test version (no zero suppression) | |
69 | 2 final version (no zero suppression) | |
70 | 3 test version (with zero suppression) | |
71 | ||
72 | The default is still set to 2 in AliTRDfeeParam::fgkRAWversion and it uses | |
73 | previously existing codes. If you set this to 3, AliTRDrawData changes behavior | |
74 | to use AliTRDmcmSim with ZS. | |
75 | ||
76 | Plan is after we make sure it works stably, we delete AliTRDmcm which is obsolete. | |
0c349049 | 77 | However it still take time because tracklet part is not yet touched. |
ecf39416 | 78 | The default raw version is 2. |
79 | ||
80 | Ken Oyama | |
81 | */ | |
82 | ||
1d93b218 | 83 | // if no histo is drawn, these are obsolete |
84 | #include <TH1.h> | |
85 | #include <TCanvas.h> | |
86 | ||
87 | // only needed if I/O of tracklets is activated | |
88 | #include <TObject.h> | |
89 | #include <TFile.h> | |
90 | #include <TTree.h> | |
91 | #include <TSystem.h> | |
92 | ||
ecf39416 | 93 | #include <fstream> |
0c349049 | 94 | |
ecf39416 | 95 | #include <TMath.h> |
0c349049 | 96 | |
dfd03fc3 | 97 | #include "AliLog.h" |
0c349049 | 98 | |
dfd03fc3 | 99 | #include "AliTRDmcmSim.h" |
100 | #include "AliTRDfeeParam.h" | |
ecf39416 | 101 | #include "AliTRDSimParam.h" |
dfd03fc3 | 102 | #include "AliTRDgeometry.h" |
ecf39416 | 103 | #include "AliTRDcalibDB.h" |
dfd03fc3 | 104 | |
1d93b218 | 105 | // additional for new tail filter and/or tracklet |
106 | #include "AliTRDtrapAlu.h" | |
107 | #include "AliTRDpadPlane.h" | |
108 | ||
109 | ||
dfd03fc3 | 110 | ClassImp(AliTRDmcmSim) |
111 | ||
112 | //_____________________________________________________________________________ | |
113 | AliTRDmcmSim::AliTRDmcmSim() :TObject() | |
114 | ,fInitialized(kFALSE) | |
1d93b218 | 115 | ,fNextEvent(-1) //new |
116 | ,fMaxTracklets(-1) //new | |
6e5d4cb2 | 117 | ,fChaId(-1) |
118 | ,fSector(-1) | |
119 | ,fStack(-1) | |
120 | ,fLayer(-1) | |
121 | ,fRobPos(-1) | |
122 | ,fMcmPos(-1) | |
123 | ,fNADC(-1) | |
124 | ,fNTimeBin(-1) | |
125 | ,fRow(-1) | |
126 | ,fADCR(NULL) | |
127 | ,fADCF(NULL) | |
1d93b218 | 128 | ,fADCT(NULL) //new |
129 | ,fPosLUT(NULL) //new | |
130 | ,fMCMT(NULL) //new | |
6e5d4cb2 | 131 | ,fZSM(NULL) |
132 | ,fZSM1Dim(NULL) | |
dfd03fc3 | 133 | ,fFeeParam(NULL) |
ecf39416 | 134 | ,fSimParam(NULL) |
135 | ,fCal(NULL) | |
dfd03fc3 | 136 | ,fGeo(NULL) |
6e5d4cb2 | 137 | { |
138 | // | |
139 | // AliTRDmcmSim default constructor | |
140 | // | |
141 | ||
142 | // By default, nothing is initialized. | |
143 | // It is necessary to issue Init before use. | |
144 | } | |
145 | ||
146 | //_____________________________________________________________________________ | |
147 | AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m) | |
148 | :TObject(m) | |
1d93b218 | 149 | ,fInitialized(kFALSE) |
150 | ,fNextEvent(-1) //new | |
151 | ,fMaxTracklets(-1) //new | |
dfd03fc3 | 152 | ,fChaId(-1) |
153 | ,fSector(-1) | |
154 | ,fStack(-1) | |
155 | ,fLayer(-1) | |
156 | ,fRobPos(-1) | |
157 | ,fMcmPos(-1) | |
158 | ,fNADC(-1) | |
159 | ,fNTimeBin(-1) | |
160 | ,fRow(-1) | |
dfd03fc3 | 161 | ,fADCR(NULL) |
162 | ,fADCF(NULL) | |
1d93b218 | 163 | ,fADCT(NULL) //new |
164 | ,fPosLUT(NULL) //new | |
165 | ,fMCMT(NULL) //new | |
dfd03fc3 | 166 | ,fZSM(NULL) |
167 | ,fZSM1Dim(NULL) | |
6e5d4cb2 | 168 | ,fFeeParam(NULL) |
169 | ,fSimParam(NULL) | |
170 | ,fCal(NULL) | |
171 | ,fGeo(NULL) | |
1d93b218 | 172 | |
dfd03fc3 | 173 | { |
174 | // | |
6e5d4cb2 | 175 | // AliTRDmcmSim copy constructor |
dfd03fc3 | 176 | // |
177 | ||
178 | // By default, nothing is initialized. | |
179 | // It is necessary to issue Init before use. | |
180 | } | |
181 | ||
182 | //_____________________________________________________________________________ | |
183 | AliTRDmcmSim::~AliTRDmcmSim() | |
184 | { | |
185 | // | |
186 | // AliTRDmcmSim destructor | |
187 | // | |
0c349049 | 188 | |
dfd03fc3 | 189 | if( fADCR != NULL ) { |
190 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
16e077d0 | 191 | delete [] fADCR[iadc]; |
192 | delete [] fADCF[iadc]; | |
1d93b218 | 193 | delete [] fADCT[iadc]; |
16e077d0 | 194 | delete [] fZSM [iadc]; |
dfd03fc3 | 195 | } |
16e077d0 | 196 | delete [] fADCR; |
197 | delete [] fADCF; | |
1d93b218 | 198 | delete [] fADCT; |
16e077d0 | 199 | delete [] fZSM; |
200 | delete [] fZSM1Dim; | |
dfd03fc3 | 201 | } |
1d93b218 | 202 | |
203 | if(fInitialized){ | |
204 | delete [] fPosLUT; | |
205 | delete [] fMCMT; | |
206 | } | |
207 | ||
dfd03fc3 | 208 | delete fGeo; |
6e5d4cb2 | 209 | |
210 | } | |
211 | ||
212 | //_____________________________________________________________________________ | |
213 | AliTRDmcmSim &AliTRDmcmSim::operator=(const AliTRDmcmSim &m) | |
214 | { | |
215 | // | |
216 | // Assignment operator | |
217 | // | |
218 | ||
219 | if (this != &m) { | |
220 | ((AliTRDmcmSim &) m).Copy(*this); | |
221 | } | |
222 | return *this; | |
223 | ||
224 | } | |
225 | ||
226 | //_____________________________________________________________________________ | |
227 | void AliTRDmcmSim::Copy(TObject &m) const | |
228 | { | |
229 | // | |
230 | // Copy function | |
231 | // | |
1d93b218 | 232 | ((AliTRDmcmSim &) m).fNextEvent = 0; //new |
233 | ((AliTRDmcmSim &) m).fMaxTracklets = 0; //new | |
6e5d4cb2 | 234 | ((AliTRDmcmSim &) m).fInitialized = 0; |
235 | ((AliTRDmcmSim &) m).fChaId = 0; | |
236 | ((AliTRDmcmSim &) m).fSector = 0; | |
237 | ((AliTRDmcmSim &) m).fStack = 0; | |
238 | ((AliTRDmcmSim &) m).fLayer = 0; | |
239 | ((AliTRDmcmSim &) m).fRobPos = 0; | |
240 | ((AliTRDmcmSim &) m).fMcmPos = 0; | |
241 | ((AliTRDmcmSim &) m).fNADC = 0; | |
242 | ((AliTRDmcmSim &) m).fNTimeBin = 0; | |
243 | ((AliTRDmcmSim &) m).fRow = 0; | |
244 | ((AliTRDmcmSim &) m).fADCR = 0; | |
245 | ((AliTRDmcmSim &) m).fADCF = 0; | |
1d93b218 | 246 | ((AliTRDmcmSim &) m).fADCT = 0; //new |
247 | ((AliTRDmcmSim &) m).fPosLUT = 0; //new | |
248 | ((AliTRDmcmSim &) m).fMCMT = 0; //new | |
6e5d4cb2 | 249 | ((AliTRDmcmSim &) m).fZSM = 0; |
250 | ((AliTRDmcmSim &) m).fZSM1Dim = 0; | |
251 | ((AliTRDmcmSim &) m).fFeeParam = 0; | |
252 | ((AliTRDmcmSim &) m).fSimParam = 0; | |
253 | ((AliTRDmcmSim &) m).fCal = 0; | |
254 | ((AliTRDmcmSim &) m).fGeo = 0; | |
255 | ||
dfd03fc3 | 256 | } |
257 | ||
258 | //_____________________________________________________________________________ | |
1d93b218 | 259 | //void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos, Bool_t newEvent ) // only for readout tree (new event) |
260 | void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos ) // only for readout tree (new event) | |
dfd03fc3 | 261 | { |
0c349049 | 262 | // |
dfd03fc3 | 263 | // Initialize the class with new geometry information |
264 | // fADC array will be reused with filled by zero | |
0c349049 | 265 | // |
dfd03fc3 | 266 | |
1d93b218 | 267 | fNextEvent = 0; //**!!new!!** |
dfd03fc3 | 268 | fFeeParam = AliTRDfeeParam::Instance(); |
ecf39416 | 269 | fSimParam = AliTRDSimParam::Instance(); |
270 | fCal = AliTRDcalibDB::Instance(); | |
dfd03fc3 | 271 | fGeo = new AliTRDgeometry(); |
0c349049 | 272 | fChaId = chaId; |
dfd03fc3 | 273 | fSector = fGeo->GetSector( fChaId ); |
274 | fStack = fGeo->GetChamber( fChaId ); | |
275 | fLayer = fGeo->GetPlane( fChaId ); | |
0c349049 | 276 | fRobPos = robPos; |
277 | fMcmPos = mcmPos; | |
dfd03fc3 | 278 | fNADC = fFeeParam->GetNadcMcm(); |
ecf39416 | 279 | fNTimeBin = fCal->GetNumberOfTimeBins(); |
dfd03fc3 | 280 | fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos ); |
dfd03fc3 | 281 | |
1d93b218 | 282 | fMaxTracklets = fFeeParam->GetMaxNrOfTracklets(); |
283 | ||
284 | ||
285 | ||
286 | //if ((fChaId == 0 && fRobPos == 0 && fMcmPos == 0)) { | |
287 | /* if (newEvent == kTRUE) { | |
288 | fNextEvent = 1; | |
289 | }*/ | |
290 | fNextEvent = 0; | |
291 | ||
292 | ||
dfd03fc3 | 293 | // Allocate ADC data memory if not yet done |
294 | if( fADCR == NULL ) { | |
295 | fADCR = new Int_t *[fNADC]; | |
296 | fADCF = new Int_t *[fNADC]; | |
1d93b218 | 297 | fADCT = new Int_t *[fNADC]; //new |
dfd03fc3 | 298 | fZSM = new Int_t *[fNADC]; |
299 | fZSM1Dim = new Int_t [fNADC]; | |
300 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
301 | fADCR[iadc] = new Int_t[fNTimeBin]; | |
302 | fADCF[iadc] = new Int_t[fNTimeBin]; | |
1d93b218 | 303 | fADCT[iadc] = new Int_t[fNTimeBin]; //new |
dfd03fc3 | 304 | fZSM [iadc] = new Int_t[fNTimeBin]; |
305 | } | |
306 | } | |
307 | ||
308 | // Initialize ADC data | |
309 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
310 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
311 | fADCR[iadc][it] = 0; | |
312 | fADCF[iadc][it] = 0; | |
1d93b218 | 313 | fADCT[iadc][it] = -1; //new |
dfd03fc3 | 314 | fZSM [iadc][it] = 1; // Default unread = 1 |
315 | } | |
316 | fZSM1Dim[iadc] = 1; // Default unread = 1 | |
317 | } | |
ecf39416 | 318 | |
1d93b218 | 319 | //new: |
320 | fPosLUT = new Int_t[128]; | |
321 | for(Int_t i = 0; i<128; i++){ | |
322 | fPosLUT[i] = 0; | |
323 | } | |
324 | ||
325 | fMCMT = new UInt_t[fMaxTracklets]; | |
326 | for(Int_t i = 0; i < fMaxTracklets; i++) { | |
327 | fMCMT[i] = 0; | |
328 | } | |
329 | ||
330 | ||
dfd03fc3 | 331 | fInitialized = kTRUE; |
332 | } | |
333 | ||
ecf39416 | 334 | //_____________________________________________________________________________ |
335 | Bool_t AliTRDmcmSim::CheckInitialized() | |
336 | { | |
0c349049 | 337 | // |
338 | // Check whether object is initialized | |
339 | // | |
340 | ||
ecf39416 | 341 | if( ! fInitialized ) { |
342 | AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called.")); | |
343 | } | |
344 | return fInitialized; | |
345 | } | |
346 | ||
dfd03fc3 | 347 | //_____________________________________________________________________________ |
1d93b218 | 348 | |
349 | ||
350 | void AliTRDmcmSim::SetPosLUT() { | |
351 | Double_t iHi = (Double_t)fCal->GetPRFhi(); | |
352 | Double_t iLo = (Double_t)fCal->GetPRFlo(); | |
353 | Int_t nBin = fCal->GetPRFbin(); | |
354 | Int_t iOff = fLayer * nBin; | |
355 | Int_t kNplan = fGeo->Nplan(); | |
356 | ||
357 | Float_t *sPRFsmp = new Float_t[nBin*kNplan]; | |
358 | Double_t *sPRFlayer = new Double_t[nBin]; | |
359 | ||
360 | ||
361 | for(Int_t i = 0; i<nBin*kNplan; i++){ | |
362 | ||
363 | //printf("%f\n",fCal->GetSampledPRF()[i]); | |
364 | sPRFsmp[i] = fCal->GetSampledPRF()[i]; | |
365 | ||
366 | } | |
367 | ||
368 | Double_t sWidth = (iHi-iLo)/((Double_t) nBin); | |
369 | Int_t sPad = (Int_t) (1.0/sWidth); | |
370 | ||
371 | // get the PRF for actual layer (interpolated to ibin data-points; 61 measured) | |
372 | for(Int_t iBin = 0; iBin < nBin; iBin++){ | |
373 | sPRFlayer[iBin] = (Double_t)sPRFsmp[iOff+iBin]; | |
374 | } | |
375 | ||
376 | Int_t bin0 = (Int_t)(-iLo / sWidth - 0.5); // bin-nr. for pad-position 0 | |
377 | ||
378 | Int_t bin1 = (Int_t)((Double_t)(0.5 - iLo) / sWidth - 0.5); // bin-nr. for pad-position 0.5 | |
379 | bin1 = bin1 + 1; | |
380 | bin0 = bin0 + 1; //avoid negative values in aYest (start right of symmetry center) | |
381 | while (bin0-sPad<0) { | |
382 | bin0 = bin0 + 1; | |
383 | } | |
384 | while (bin1+sPad>=nBin) { | |
385 | bin1 = bin1 - 1; | |
386 | } | |
387 | ||
388 | Double_t* aYest = new Double_t[bin1-bin0+1]; | |
389 | ||
390 | /*TH1F* hist1 = new TH1F("h1","yest(y)",128,0,0.5); | |
391 | TH1F* hist2 = new TH1F("h2","y(yest)",128,0,0.5); | |
392 | TH1F* hist3 = new TH1F("h3","y(yest)-yest",128,0,0.5); | |
393 | TH1F* hist4 = new TH1F("h4","y(yest)-yest,discrete",128,0,0.5); | |
394 | ||
395 | TCanvas *c1 = new TCanvas("c1","c1",800,1000); | |
396 | hist1->Draw(); | |
397 | TCanvas *c2 = new TCanvas("c2","c2",800,1000); | |
398 | hist2->Draw(); | |
399 | TCanvas *c3 = new TCanvas("c3","c3",800,1000); | |
400 | hist3->Draw(); | |
401 | TCanvas *c4 = new TCanvas("c4","c4",800,1000); | |
402 | hist4->Draw();*/ | |
403 | ||
404 | for(Int_t iBin = bin0; iBin <= bin1; iBin++){ | |
405 | aYest[iBin-bin0] = 0.5*(sPRFlayer[iBin-sPad] - sPRFlayer[iBin+sPad])/(sPRFlayer[iBin]); // estimated position from PRF; between 0 and 1 | |
406 | //Double_t position = ((Double_t)(iBin)+0.5)*sWidth+iLo; | |
407 | // hist1->Fill(position,aYest[iBin-bin0]); | |
408 | } | |
409 | ||
410 | ||
411 | ||
412 | Double_t aY[128]; // reversed function | |
413 | ||
414 | AliTRDtrapAlu a; | |
415 | a.Init(1,8,0,31); | |
416 | ||
417 | for(Int_t j = 0; j<128; j++) { // loop over all Yest; LUT has 128 entries; | |
418 | Double_t yest = ((Double_t)j)/256; | |
419 | ||
420 | Int_t iBin = 0; | |
421 | while (yest>aYest[iBin] && iBin<(bin1-bin0)) { | |
422 | iBin = iBin+1; | |
423 | } | |
424 | if((iBin == bin1 - bin0)&&(yest>aYest[iBin])) { | |
425 | aY[j] = 0.5; // yest too big | |
426 | //hist2->Fill(yest,aY[j]); | |
427 | ||
428 | } | |
429 | else { | |
430 | Int_t bin_d = iBin + bin0 - 1; | |
431 | Int_t bin_u = iBin + bin0; | |
432 | Double_t y_d = ((Double_t)bin_d + 0.5)*sWidth + iLo; // lower y | |
433 | Double_t y_u = ((Double_t)bin_u + 0.5)*sWidth + iLo; // upper y | |
434 | Double_t yest_d = aYest[iBin-1]; // lower estimated y | |
435 | Double_t yest_u = aYest[iBin]; // upper estimated y | |
436 | ||
437 | aY[j] = ((yest-yest_d)/(yest_u-yest_d))*(y_u-y_d) + y_d; | |
438 | //hist2->Fill(yest,aY[j]); | |
439 | ||
440 | } | |
441 | aY[j] = aY[j] - yest; | |
442 | //hist3->Fill(yest,aY[j]); | |
443 | // formatting | |
444 | a.AssignDouble(aY[j]); | |
445 | //a.WriteWord(); | |
446 | fPosLUT[j] = a.GetValue(); // 1+8Bit value;128 entries;LUT is steered by abs(Q(i+1)-Q(i-1))/Q(i)=COG and gives the correction to COG/2 | |
447 | //hist4->Fill(yest,fPosLUT[j]); | |
448 | ||
449 | } | |
450 | ||
451 | ||
452 | ||
453 | delete [] sPRFsmp; | |
454 | delete [] sPRFlayer; | |
455 | delete [] aYest; | |
456 | ||
457 | } | |
458 | ||
459 | ||
460 | //_____________________________________________________________________________ | |
461 | Int_t* AliTRDmcmSim::GetPosLUT(){ | |
462 | return fPosLUT; | |
463 | } | |
464 | ||
465 | ||
466 | ||
dfd03fc3 | 467 | void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc ) |
468 | { | |
0c349049 | 469 | // |
dfd03fc3 | 470 | // Store ADC data into array of raw data |
0c349049 | 471 | // |
dfd03fc3 | 472 | |
ecf39416 | 473 | if( !CheckInitialized() ) return; |
dfd03fc3 | 474 | |
475 | if( iadc < 0 || iadc >= fNADC ) { | |
476 | //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1)); | |
477 | return; | |
478 | } | |
479 | ||
480 | for( int it = 0 ; it < fNTimeBin ; it++ ) { | |
481 | fADCR[iadc][it] = (Int_t)(adc[it]); | |
482 | } | |
483 | } | |
484 | ||
485 | //_____________________________________________________________________________ | |
486 | void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc ) | |
487 | { | |
0c349049 | 488 | // |
dfd03fc3 | 489 | // Store ADC data into array of raw data |
0c349049 | 490 | // |
dfd03fc3 | 491 | |
ecf39416 | 492 | if( !CheckInitialized() ) return; |
dfd03fc3 | 493 | |
494 | if( iadc < 0 || iadc >= fNADC ) { | |
495 | //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1)); | |
496 | return; | |
497 | } | |
498 | ||
499 | fADCR[iadc][it] = adc; | |
500 | } | |
501 | ||
502 | //_____________________________________________________________________________ | |
503 | void AliTRDmcmSim::SetDataPedestal( Int_t iadc ) | |
504 | { | |
0c349049 | 505 | // |
dfd03fc3 | 506 | // Store ADC data into array of raw data |
0c349049 | 507 | // |
dfd03fc3 | 508 | |
ecf39416 | 509 | if( !CheckInitialized() ) return; |
dfd03fc3 | 510 | |
511 | if( iadc < 0 || iadc >= fNADC ) { | |
512 | //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1)); | |
513 | return; | |
514 | } | |
515 | ||
516 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
ecf39416 | 517 | fADCR[iadc][it] = fSimParam->GetADCbaseline(); |
dfd03fc3 | 518 | } |
519 | } | |
520 | ||
521 | //_____________________________________________________________________________ | |
522 | Int_t AliTRDmcmSim::GetCol( Int_t iadc ) | |
523 | { | |
0c349049 | 524 | // |
dfd03fc3 | 525 | // Return column id of the pad for the given ADC channel |
0c349049 | 526 | // |
527 | ||
ecf39416 | 528 | if( !CheckInitialized() ) return -1; |
dfd03fc3 | 529 | |
ecf39416 | 530 | return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc); |
dfd03fc3 | 531 | } |
532 | ||
dfd03fc3 | 533 | //_____________________________________________________________________________ |
534 | Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize ) | |
535 | { | |
0c349049 | 536 | // |
dfd03fc3 | 537 | // Produce raw data stream from this MCM and put in buf |
0c349049 | 538 | // Returns number of words filled, or negative value |
539 | // with -1 * number of overflowed words | |
540 | // | |
dfd03fc3 | 541 | |
542 | UInt_t x; | |
543 | UInt_t iEv = 0; | |
544 | Int_t nw = 0; // Number of written words | |
545 | Int_t of = 0; // Number of overflowed words | |
546 | Int_t rawVer = fFeeParam->GetRAWversion(); | |
547 | Int_t **adc; | |
548 | ||
ecf39416 | 549 | if( !CheckInitialized() ) return 0; |
550 | ||
dfd03fc3 | 551 | if( fFeeParam->GetRAWstoreRaw() ) { |
552 | adc = fADCR; | |
553 | } else { | |
554 | adc = fADCF; | |
555 | } | |
556 | ||
557 | // Produce MCM header | |
558 | x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC; | |
559 | if (nw < maxSize) { | |
560 | buf[nw++] = x; | |
561 | } | |
562 | else { | |
563 | of++; | |
564 | } | |
565 | ||
566 | // Produce ADC mask | |
567 | if( rawVer >= 3 ) { | |
568 | x = 0; | |
569 | for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) { | |
570 | if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed | |
571 | x = x | (1 << iAdc); | |
572 | } | |
573 | } | |
574 | if (nw < maxSize) { | |
575 | buf[nw++] = x; | |
576 | } | |
577 | else { | |
578 | of++; | |
579 | } | |
580 | } | |
581 | ||
582 | // Produce ADC data. 3 timebins are packed into one 32 bits word | |
583 | // In this version, different ADC channel will NOT share the same word | |
584 | ||
585 | UInt_t aa=0, a1=0, a2=0, a3=0; | |
586 | ||
587 | for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) { | |
ecf39416 | 588 | if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // suppressed |
dfd03fc3 | 589 | aa = !(iAdc & 1) + 2; |
590 | for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) { | |
591 | a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] : 0; | |
592 | a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0; | |
593 | a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0; | |
ecf39416 | 594 | x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa; |
595 | if (nw < maxSize) { | |
596 | buf[nw++] = x; | |
597 | } | |
598 | else { | |
599 | of++; | |
600 | } | |
dfd03fc3 | 601 | } |
602 | } | |
603 | ||
604 | if( of != 0 ) return -of; else return nw; | |
605 | } | |
606 | ||
1d93b218 | 607 | //_____________________________________________________________________________ |
608 | Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize ) | |
609 | { | |
610 | // | |
611 | // Produce tracklet data stream from this MCM and put in buf | |
612 | // Returns number of words filled, or negative value | |
613 | // with -1 * number of overflowed words | |
614 | // | |
615 | ||
616 | UInt_t x; | |
617 | Int_t nw = 0; // Number of written words | |
618 | Int_t of = 0; // Number of overflowed words | |
619 | ||
620 | if( !CheckInitialized() ) return 0; | |
621 | ||
622 | // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM | |
623 | // fMCMT is filled continuously until no more tracklet words available | |
624 | ||
625 | Int_t wd = 0; | |
626 | while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){ | |
627 | x = fMCMT[wd]; | |
628 | if (nw < maxSize) { | |
629 | buf[nw++] = x; | |
630 | } | |
631 | else { | |
632 | of++; | |
633 | } | |
634 | wd++; | |
635 | } | |
636 | ||
637 | if( of != 0 ) return -of; else return nw; | |
638 | } | |
639 | ||
640 | ||
dfd03fc3 | 641 | //_____________________________________________________________________________ |
642 | void AliTRDmcmSim::Filter() | |
643 | { | |
0c349049 | 644 | // |
dfd03fc3 | 645 | // Apply digital filter |
0c349049 | 646 | // |
dfd03fc3 | 647 | |
ecf39416 | 648 | if( !CheckInitialized() ) return; |
dfd03fc3 | 649 | |
650 | // Initialize filtered data array with raw data | |
651 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
652 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
653 | fADCF[iadc][it] = fADCR[iadc][it]; | |
654 | } | |
655 | } | |
656 | ||
657 | // Then apply fileters one by one to filtered data array | |
acc49af9 | 658 | if( fFeeParam->IsPFon() ) FilterPedestal(); |
659 | if( fFeeParam->IsGFon() ) FilterGain(); | |
660 | if( fFeeParam->IsTFon() ) FilterTail(); | |
dfd03fc3 | 661 | } |
662 | ||
663 | //_____________________________________________________________________________ | |
664 | void AliTRDmcmSim::FilterPedestal() | |
665 | { | |
0c349049 | 666 | // |
dfd03fc3 | 667 | // Apply pedestal filter |
0c349049 | 668 | // |
dfd03fc3 | 669 | |
ecf39416 | 670 | Int_t ap = fSimParam->GetADCbaseline(); // ADC instrinsic pedestal |
dfd03fc3 | 671 | Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal |
6e5d4cb2 | 672 | //Int_t tc = fFeeParam->GetPFtimeConstant(); // this makes no sense yet |
dfd03fc3 | 673 | |
674 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
675 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
676 | fADCF[iadc][it] = fADCF[iadc][it] - ap + ep; | |
677 | } | |
678 | } | |
679 | } | |
680 | ||
681 | //_____________________________________________________________________________ | |
682 | void AliTRDmcmSim::FilterGain() | |
683 | { | |
0c349049 | 684 | // |
dfd03fc3 | 685 | // Apply gain filter (not implemented) |
ecf39416 | 686 | // Later it will be implemented because gain digital filiter will |
687 | // increase noise level. | |
0c349049 | 688 | // |
689 | ||
dfd03fc3 | 690 | } |
691 | ||
692 | //_____________________________________________________________________________ | |
693 | void AliTRDmcmSim::FilterTail() | |
694 | { | |
0c349049 | 695 | // |
dfd03fc3 | 696 | // Apply exponential tail filter (Bogdan's version) |
0c349049 | 697 | // |
dfd03fc3 | 698 | |
699 | Double_t *dtarg = new Double_t[fNTimeBin]; | |
700 | Int_t *itarg = new Int_t[fNTimeBin]; | |
701 | Int_t nexp = fFeeParam->GetTFnExp(); | |
702 | Int_t tftype = fFeeParam->GetTFtype(); | |
703 | ||
704 | switch( tftype ) { | |
705 | ||
706 | case 0: // Exponential Filter Analog Bogdan | |
707 | for (Int_t iCol = 0; iCol < fNADC; iCol++) { | |
708 | FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp); | |
709 | for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) { | |
710 | fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]); | |
711 | } | |
712 | } | |
713 | break; | |
714 | ||
715 | case 1: // Exponential filter digital Bogdan | |
716 | for (Int_t iCol = 0; iCol < fNADC; iCol++) { | |
717 | FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp); | |
718 | for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) { | |
719 | fADCF[iCol][iTime] = itarg[iTime]; | |
720 | } | |
721 | } | |
722 | break; | |
723 | ||
724 | case 2: // Exponential filter Marian special | |
725 | for (Int_t iCol = 0; iCol < fNADC; iCol++) { | |
726 | FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin); | |
727 | for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) { | |
728 | fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]); | |
729 | } | |
730 | } | |
731 | break; | |
1d93b218 | 732 | |
733 | //new | |
734 | case 3: // Exponential filter using AliTRDtrapAlu class | |
735 | for (Int_t iCol = 0; iCol < fNADC; iCol++) { | |
736 | FilterSimDeConvExpEl( fADCF[iCol], itarg, fNTimeBin, nexp); | |
737 | for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) { | |
738 | fADCF[iCol][iTime] = itarg[iTime]>>2; // to be used for raw-data | |
739 | fADCT[iCol][iTime] = itarg[iTime]; // 12bits; to be used for tracklet; tracklet will have own container; | |
740 | } | |
741 | } | |
742 | break; | |
743 | ||
dfd03fc3 | 744 | |
745 | default: | |
746 | AliError(Form("Invalid filter type %d ! \n", tftype )); | |
747 | break; | |
748 | } | |
749 | ||
750 | delete dtarg; | |
751 | delete itarg; | |
752 | } | |
753 | ||
754 | //_____________________________________________________________________________ | |
755 | void AliTRDmcmSim::ZSMapping() | |
756 | { | |
0c349049 | 757 | // |
dfd03fc3 | 758 | // Zero Suppression Mapping implemented in TRAP chip |
759 | // | |
760 | // See detail TRAP manual "Data Indication" section: | |
761 | // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf | |
0c349049 | 762 | // |
dfd03fc3 | 763 | |
0c349049 | 764 | Int_t eBIS = fFeeParam->GetEBsglIndThr(); // TRAP default = 0x4 (Tis=4) |
765 | Int_t eBIT = fFeeParam->GetEBsumIndThr(); // TRAP default = 0x28 (Tit=40) | |
766 | Int_t eBIL = fFeeParam->GetEBindLUT(); // TRAP default = 0xf0 | |
767 | // (lookup table accept (I2,I1,I0)=(111) | |
768 | // or (110) or (101) or (100)) | |
769 | Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity) | |
ecf39416 | 770 | Int_t ep = AliTRDfeeParam::GetPFeffectPedestal(); |
771 | ||
772 | if( !CheckInitialized() ) return; | |
dfd03fc3 | 773 | |
774 | for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) { | |
775 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
776 | ||
ecf39416 | 777 | // Get ADC data currently in filter buffer |
0c349049 | 778 | Int_t ap = fADCF[iadc-1][it] - ep; // previous |
779 | Int_t ac = fADCF[iadc ][it] - ep; // current | |
780 | Int_t an = fADCF[iadc+1][it] - ep; // next | |
ecf39416 | 781 | |
dfd03fc3 | 782 | // evaluate three conditions |
0c349049 | 783 | Int_t i0 = ( ac >= ap && ac >= an ) ? 0 : 1; // peak center detection |
784 | Int_t i1 = ( ap + ac + an > eBIT ) ? 0 : 1; // cluster | |
785 | Int_t i2 = ( ac > eBIS ) ? 0 : 1; // absolute large peak | |
786 | ||
787 | Int_t i = i2 * 4 + i1 * 2 + i0; // Bit position in lookup table | |
788 | Int_t d = (eBIL >> i) & 1; // Looking up (here d=0 means true | |
789 | // and d=1 means false according to TRAP manual) | |
790 | ||
791 | fZSM[iadc][it] &= d; | |
792 | if( eBIN == 0 ) { // turn on neighboring ADCs | |
793 | fZSM[iadc-1][it] &= d; | |
794 | fZSM[iadc+1][it] &= d; | |
dfd03fc3 | 795 | } |
0c349049 | 796 | |
dfd03fc3 | 797 | } |
798 | } | |
799 | ||
800 | // do 1 dim projection | |
801 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
802 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
ecf39416 | 803 | fZSM1Dim[iadc] &= fZSM[iadc][it]; |
804 | } | |
805 | } | |
0c349049 | 806 | |
ecf39416 | 807 | } |
808 | ||
809 | //_____________________________________________________________________________ | |
810 | void AliTRDmcmSim::DumpData( char *f, char *target ) | |
811 | { | |
0c349049 | 812 | // |
ecf39416 | 813 | // Dump data stored (for debugging). |
814 | // target should contain one or multiple of the following characters | |
815 | // R for raw data | |
816 | // F for filtered data | |
817 | // Z for zero suppression map | |
818 | // S Raw dat astream | |
819 | // other characters are simply ignored | |
0c349049 | 820 | // |
821 | ||
ecf39416 | 822 | UInt_t tempbuf[1024]; |
823 | ||
824 | if( !CheckInitialized() ) return; | |
825 | ||
826 | std::ofstream of( f, std::ios::out | std::ios::app ); | |
827 | of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n", | |
828 | fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos ); | |
829 | ||
830 | for( int t=0 ; target[t] != 0 ; t++ ) { | |
831 | switch( target[t] ) { | |
832 | case 'R' : | |
833 | case 'r' : | |
834 | of << Form("fADCR (raw ADC data)\n"); | |
835 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
836 | of << Form(" ADC %02d: ", iadc); | |
837 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
838 | of << Form("% 4d", fADCR[iadc][it]); | |
839 | } | |
840 | of << Form("\n"); | |
841 | } | |
842 | break; | |
843 | case 'F' : | |
844 | case 'f' : | |
845 | of << Form("fADCF (filtered ADC data)\n"); | |
846 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
847 | of << Form(" ADC %02d: ", iadc); | |
848 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
849 | of << Form("% 4d", fADCF[iadc][it]); | |
850 | } | |
851 | of << Form("\n"); | |
852 | } | |
853 | break; | |
854 | case 'Z' : | |
855 | case 'z' : | |
856 | of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n"); | |
857 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
858 | of << Form(" ADC %02d: ", iadc); | |
859 | if( fZSM1Dim[iadc] == 0 ) { of << " R " ; } else { of << " . "; } // R:read .:suppressed | |
860 | for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { | |
861 | if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed | |
862 | } | |
863 | of << Form("\n"); | |
864 | } | |
865 | break; | |
866 | case 'S' : | |
867 | case 's' : | |
868 | Int_t s = ProduceRawStream( tempbuf, 1024 ); | |
869 | of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion()); | |
870 | of << Form(" address data\n"); | |
871 | for( int i = 0 ; i < s ; i++ ) { | |
872 | of << Form(" %04x %08x\n", i, tempbuf[i]); | |
873 | } | |
dfd03fc3 | 874 | } |
875 | } | |
876 | } | |
877 | ||
878 | //_____________________________________________________________________________ | |
0c349049 | 879 | void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target |
880 | , Int_t n, Int_t nexp) | |
dfd03fc3 | 881 | { |
882 | // | |
883 | // Exponential filter "analog" | |
884 | // source will not be changed | |
885 | // | |
886 | ||
887 | Int_t i = 0; | |
888 | Int_t k = 0; | |
889 | Double_t reminder[2]; | |
890 | Double_t correction; | |
891 | Double_t result; | |
892 | Double_t rates[2]; | |
893 | Double_t coefficients[2]; | |
894 | ||
895 | // Initialize (coefficient = alpha, rates = lambda) | |
896 | // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02) | |
897 | ||
898 | Double_t r1 = (Double_t)fFeeParam->GetTFr1(); | |
899 | Double_t r2 = (Double_t)fFeeParam->GetTFr2(); | |
900 | Double_t c1 = (Double_t)fFeeParam->GetTFc1(); | |
901 | Double_t c2 = (Double_t)fFeeParam->GetTFc2(); | |
902 | ||
903 | coefficients[0] = c1; | |
904 | coefficients[1] = c2; | |
905 | ||
906 | Double_t dt = 0.1; | |
907 | rates[0] = TMath::Exp(-dt/(r1)); | |
908 | rates[1] = TMath::Exp(-dt/(r2)); | |
909 | ||
910 | // Attention: computation order is important | |
911 | correction = 0.0; | |
912 | for (k = 0; k < nexp; k++) { | |
913 | reminder[k] = 0.0; | |
914 | } | |
915 | ||
916 | for (i = 0; i < n; i++) { | |
917 | ||
918 | result = ((Double_t)source[i] - correction); // no rescaling | |
919 | target[i] = result; | |
920 | ||
921 | for (k = 0; k < nexp; k++) { | |
922 | reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result); | |
923 | } | |
924 | ||
925 | correction = 0.0; | |
926 | for (k = 0; k < nexp; k++) { | |
927 | correction += reminder[k]; | |
928 | } | |
929 | } | |
930 | } | |
931 | ||
932 | //_____________________________________________________________________________ | |
0c349049 | 933 | void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n |
934 | , Int_t nexp) | |
dfd03fc3 | 935 | { |
936 | // | |
937 | // Exponential filter "digital" | |
938 | // source will not be changed | |
939 | // | |
940 | ||
ecf39416 | 941 | Int_t i = 0; |
dfd03fc3 | 942 | Int_t fAlphaL = 0; |
943 | Int_t fAlphaS = 0; | |
dfd03fc3 | 944 | Int_t fTailPed = 0; |
dfd03fc3 | 945 | Int_t iAlphaL = 0; |
946 | Int_t iAlphaS = 0; | |
dfd03fc3 | 947 | |
948 | // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02) | |
949 | // initialize (coefficient = alpha, rates = lambda) | |
950 | ||
951 | Double_t dt = 0.1; | |
952 | Double_t r1 = (Double_t)fFeeParam->GetTFr1(); | |
953 | Double_t r2 = (Double_t)fFeeParam->GetTFr2(); | |
954 | Double_t c1 = (Double_t)fFeeParam->GetTFc1(); | |
955 | Double_t c2 = (Double_t)fFeeParam->GetTFc2(); | |
956 | ||
ecf39416 | 957 | Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0); |
958 | Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0); | |
959 | Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits | |
960 | Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits | |
dfd03fc3 | 961 | |
962 | if (nexp == 1) { | |
963 | fAlphaL = (Int_t) (c1 * 2048.0); | |
964 | iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter | |
965 | } | |
966 | if (nexp == 2) { | |
967 | fAlphaL = (Int_t) (c1 * 2048.0); | |
968 | fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0); | |
969 | iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter | |
970 | iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits | |
971 | } | |
972 | ||
973 | Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers | |
974 | Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers | |
975 | Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers | |
976 | Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers | |
977 | ||
978 | Int_t h1; | |
979 | Int_t h2; | |
980 | Int_t rem1; | |
981 | Int_t rem2; | |
982 | Int_t correction; | |
983 | Int_t result; | |
984 | Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2; | |
985 | ||
986 | Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the | |
987 | rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects. | |
988 | rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs)); | |
989 | ||
990 | // further initialization | |
991 | if ((rem1 + rem2) > 0x0FFF) { | |
992 | correction = 0x0FFF; | |
993 | } | |
994 | else { | |
995 | correction = (rem1 + rem2) & 0x0FFF; | |
996 | } | |
997 | ||
998 | fTailPed = iFactor - correction; | |
999 | ||
1000 | for (i = 0; i < n; i++) { | |
1001 | ||
ecf39416 | 1002 | result = (source[i] - correction); |
1003 | if (result < 0) { // Too much undershoot | |
dfd03fc3 | 1004 | result = 0; |
1005 | } | |
1006 | ||
1007 | target[i] = result; | |
1008 | ||
1009 | h1 = (rem1 + ((iAlphaL * result) >> 11)); | |
1010 | if (h1 > 0x0FFF) { | |
1011 | h1 = 0x0FFF; | |
1012 | } | |
1013 | else { | |
1014 | h1 &= 0x0FFF; | |
1015 | } | |
1016 | ||
1017 | h2 = (rem2 + ((iAlphaS * result) >> 11)); | |
1018 | if (h2 > 0x0FFF) { | |
1019 | h2 = 0x0FFF; | |
1020 | } | |
1021 | else { | |
1022 | h2 &= 0x0FFF; | |
1023 | } | |
1024 | ||
1025 | rem1 = (iLambdaL * h1 ) >> 11; | |
1026 | rem2 = (iLambdaS * h2 ) >> 11; | |
1027 | ||
1028 | if ((rem1 + rem2) > 0x0FFF) { | |
1029 | correction = 0x0FFF; | |
1030 | } | |
1031 | else { | |
1032 | correction = (rem1 + rem2) & 0x0FFF; | |
1033 | } | |
1034 | ||
1035 | } | |
0c349049 | 1036 | |
dfd03fc3 | 1037 | } |
1038 | ||
1039 | //_____________________________________________________________________________ | |
0c349049 | 1040 | void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target |
1041 | , Int_t n) | |
dfd03fc3 | 1042 | { |
1043 | // | |
1044 | // Exponential filter (M. Ivanov) | |
1045 | // source will not be changed | |
1046 | // | |
1047 | ||
1048 | Int_t i = 0; | |
1049 | Double_t sig1[100]; | |
1050 | Double_t sig2[100]; | |
1051 | Double_t sig3[100]; | |
1052 | ||
1053 | for (i = 0; i < n; i++) { | |
1054 | sig1[i] = (Double_t)source[i]; | |
1055 | } | |
1056 | ||
1057 | Float_t dt = 0.1; | |
1058 | Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt; | |
1059 | Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt; | |
1060 | ||
1061 | FilterSimTailMakerSpline( sig1, sig2, lambda0, n); | |
1062 | FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n); | |
1063 | ||
1064 | for (i = 0; i < n; i++) { | |
1065 | target[i] = sig3[i]; | |
1066 | } | |
1067 | ||
1068 | } | |
1069 | ||
1070 | //______________________________________________________________________________ | |
0c349049 | 1071 | void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout |
1072 | , Double_t lambda, Int_t n) | |
dfd03fc3 | 1073 | { |
1074 | // | |
1075 | // Special filter (M. Ivanov) | |
1076 | // | |
1077 | ||
1078 | Int_t i = 0; | |
1079 | Double_t l = TMath::Exp(-lambda*0.5); | |
1080 | Double_t in[1000]; | |
1081 | Double_t out[1000]; | |
1082 | ||
1083 | // Initialize in[] and out[] goes 0 ... 2*n+19 | |
1084 | for (i = 0; i < n*2+20; i++) { | |
1085 | in[i] = out[i] = 0; | |
1086 | } | |
1087 | ||
1088 | // in[] goes 0, 1 | |
1089 | in[0] = ampin[0]; | |
1090 | in[1] = (ampin[0] + ampin[1]) * 0.5; | |
1091 | ||
1092 | // Add charge to the end | |
1093 | for (i = 0; i < 22; i++) { | |
1094 | in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19 | |
1095 | } | |
1096 | ||
1097 | // Use arithmetic mean | |
1098 | for (i = 1; i < n-1; i++) { | |
1099 | in[2*i] = ampin[i]; // in[] goes 2, 3, ... , 2*n-4, 2*n-3 | |
1100 | in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.; | |
1101 | } | |
1102 | ||
1103 | Double_t temp; | |
1104 | out[2*n] = in[2*n]; | |
1105 | temp = 0; | |
1106 | for (i = 2*n; i >= 0; i--) { | |
1107 | out[i] = in[i] + temp; | |
1108 | temp = l*(temp+in[i]); | |
1109 | } | |
1110 | ||
1111 | for (i = 0; i < n; i++){ | |
1112 | //ampout[i] = out[2*i+1]; // org | |
1113 | ampout[i] = out[2*i]; | |
1114 | } | |
1115 | ||
1116 | } | |
1117 | ||
1118 | //______________________________________________________________________________ | |
0c349049 | 1119 | void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout |
1120 | , Double_t norm, Double_t lambda | |
1121 | , Int_t n) | |
dfd03fc3 | 1122 | { |
1123 | // | |
1124 | // Special filter (M. Ivanov) | |
1125 | // | |
1126 | ||
1127 | Int_t i = 0; | |
1128 | ||
1129 | Double_t l = TMath::Exp(-lambda*0.5); | |
1130 | Double_t k = l*(1.0 - norm*lambda*0.5); | |
1131 | Double_t in[1000]; | |
1132 | Double_t out[1000]; | |
1133 | ||
1134 | // Initialize in[] and out[] goes 0 ... 2*n+19 | |
1135 | for (i = 0; i < n*2+20; i++) { | |
1136 | in[i] = out[i] = 0; | |
1137 | } | |
1138 | ||
1139 | // in[] goes 0, 1 | |
1140 | in[0] = ampin[0]; | |
1141 | in[1] = (ampin[0]+ampin[1])*0.5; | |
1142 | ||
1143 | // Add charge to the end | |
1144 | for (i =-2; i < 22; i++) { | |
1145 | // in[] goes 2*n-4, 2*n-3, ... , 2*n+19 | |
1146 | in[2*(n-1)+i] = ampin[n-1]; | |
1147 | } | |
1148 | ||
1149 | for (i = 1; i < n-2; i++) { | |
1150 | // in[] goes 2, 3, ... , 2*n-6, 2*n-5 | |
1151 | in[2*i] = ampin[i]; | |
1152 | in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0; | |
1153 | //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0; | |
1154 | } | |
1155 | ||
1156 | Double_t temp; | |
1157 | out[0] = in[0]; | |
1158 | temp = in[0]; | |
1159 | for (i = 1; i <= 2*n; i++) { | |
1160 | out[i] = in[i] + (k-l)*temp; | |
1161 | temp = in[i] + k *temp; | |
1162 | } | |
1163 | ||
1164 | for (i = 0; i < n; i++) { | |
1165 | //ampout[i] = out[2*i+1]; // org | |
1166 | //ampout[i] = TMath::Max(out[2*i+1],0.0); // org | |
1167 | ampout[i] = TMath::Max(out[2*i],0.0); | |
1168 | } | |
1169 | } | |
1170 | ||
1d93b218 | 1171 | |
1172 | //_____________________________________________________________________________________ | |
1173 | //the following filter uses AliTRDtrapAlu-class | |
1174 | ||
1175 | void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) { | |
1176 | //static Int_t count = 0; | |
1177 | ||
1178 | Double_t dt = 0.1; | |
1179 | Double_t r1 = (Double_t)fFeeParam->GetTFr1(); | |
1180 | Double_t r2 = (Double_t)fFeeParam->GetTFr2(); | |
1181 | Double_t c1 = (Double_t)fFeeParam->GetTFc1(); | |
1182 | Double_t c2 = (Double_t)fFeeParam->GetTFc2(); | |
1183 | ||
1184 | nexp = 1; | |
1185 | ||
1186 | //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual | |
1187 | //parameters need to be adjusted | |
1188 | AliTRDtrapAlu lambdaL; | |
1189 | AliTRDtrapAlu lambdaS; | |
1190 | AliTRDtrapAlu alphaL; | |
1191 | AliTRDtrapAlu alphaS; | |
1192 | ||
1193 | AliTRDtrapAlu correction; | |
1194 | AliTRDtrapAlu result; | |
1195 | AliTRDtrapAlu bufL; | |
1196 | AliTRDtrapAlu bufS; | |
1197 | ||
1198 | AliTRDtrapAlu bSource; | |
1199 | ||
1200 | lambdaL.Init(1,11); | |
1201 | lambdaS.Init(1,11); | |
1202 | alphaL.Init(1,11); | |
1203 | alphaS.Init(1,11); | |
1204 | ||
1205 | //count=count+1; | |
1206 | ||
1207 | lambdaL.AssignDouble(TMath::Exp(-dt/r1)); | |
1208 | lambdaS.AssignDouble(TMath::Exp(-dt/r2)); | |
1209 | alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants | |
1210 | alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL | |
1211 | ||
1212 | //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly | |
1213 | correction.Init(10,2); | |
1214 | result.Init(10,2); | |
1215 | bufL.Init(10,2); | |
1216 | bufS.Init(10,2); | |
1217 | bSource.Init(10,2); | |
1218 | ||
1219 | for(Int_t i = 0; i < n; i++) { | |
1220 | bSource.AssignInt(source[i]); | |
1221 | result = bSource - correction; // subtraction can produce an underflow | |
1222 | if(result.GetSign() == kTRUE) { | |
1223 | result.AssignInt(0); | |
1224 | } | |
1225 | ||
1226 | //target[i] = result.GetValuePre(); // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly | |
1227 | ||
1228 | target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2 | |
1229 | ||
1230 | //printf("target-Wert zur Zeit %d : %d",i,target[i]); | |
1231 | //printf("\n"); | |
1232 | ||
1233 | bufL = bufL + (result * alphaL); | |
1234 | bufL = bufL * lambdaL; | |
1235 | ||
1236 | bufS = bufS + (result * alphaS); | |
1237 | bufS = bufS * lambdaS; // eventually this should look like: | |
1238 | // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more | |
1239 | ||
1240 | correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF | |
1241 | } | |
1242 | ||
1243 | ||
1244 | } | |
1245 | ||
1246 | ||
1247 | ||
1248 | ||
1249 | ||
1250 | ||
1251 | ||
1252 | ||
1253 | //__________________________________________________________________________________ | |
1254 | ||
1255 | ||
1256 | // in order to use the Tracklets, please first | |
1257 | // -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation | |
1258 | // -- set AliTRDfeeParam::fgkTFtype to 3, in order to use the new tail cancellation filter | |
1259 | // currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3 | |
1260 | ||
1261 | // code is designed such that the less possible calculations with AliTRDtrapAlu class-objects are performed; whenever possible calculations are done with doubles or integers and the results are transformed into the right format | |
1262 | ||
1263 | void AliTRDmcmSim::Tracklet(){ | |
1264 | // tracklet calculation | |
1265 | // if you use this code outside a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam | |
1266 | ||
1267 | if(!CheckInitialized()){ return; } | |
1268 | ||
1269 | Bool_t filtered = kTRUE; | |
1270 | ||
1271 | ||
1272 | ||
1273 | AliTRDtrapAlu data; | |
1274 | data.Init(10,2); | |
1275 | if(fADCT[0][0]==-1){ // check if filter was applied | |
1276 | filtered = kFALSE; | |
1277 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
1278 | for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) { | |
1279 | data.AssignInt(fADCR[iadc][iT]); | |
1280 | fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called this is done correctly | |
1281 | } | |
1282 | } | |
1283 | ||
1284 | } | |
1285 | ||
1286 | // the online ordering of mcm's is reverse to the TRAP-ordering(?)! reverse fADCT (to be consistent to TRAP), then do all calculations | |
1287 | // reverse fADCT: | |
1288 | Int_t** rev0 = new Int_t *[fNADC]; | |
1289 | Int_t** rev1 = new Int_t *[fNADC]; | |
1290 | ||
1291 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
1292 | rev0[iadc] = new Int_t[fNTimeBin]; | |
1293 | rev1[iadc] = new Int_t[fNTimeBin]; | |
1294 | for( Int_t iT = 0; iT < fNTimeBin; iT++) { | |
1295 | if( iadc <= fNADC-iadc-1 ) { | |
1296 | rev0[iadc][iT] = fADCT[fNADC-iadc-1][iT]; | |
1297 | rev1[iadc][iT] = fADCT[iadc][iT]; | |
1298 | fADCT[iadc][iT] = rev0[iadc][iT]; | |
1299 | } | |
1300 | else { | |
1301 | rev0[iadc][iT] = rev1[fNADC-iadc-1][iT]; | |
1302 | fADCT[iadc][iT] = rev0[iadc][iT]; | |
1303 | } | |
1304 | } | |
1305 | } | |
1306 | for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { | |
1307 | delete[] rev0[iadc]; | |
1308 | delete[] rev1[iadc]; | |
1309 | } | |
1310 | ||
1311 | delete[] rev0; | |
1312 | delete[] rev1; | |
1313 | ||
1314 | rev0 = NULL; | |
1315 | rev1 = NULL; | |
1316 | ||
1317 | // get the filtered pedestal; supports only electronic tail-cancellation filter | |
1318 | AliTRDtrapAlu filPed; | |
1319 | Int_t ep = 0; | |
1320 | Int_t *ieffped = new Int_t[fNTimeBin]; | |
1321 | for(Int_t iT = 0; iT < fNTimeBin; iT++){ | |
1322 | ieffped[iT] = ep; | |
1323 | } | |
1324 | ||
1325 | if( filtered == kTRUE ) { | |
1326 | if( fFeeParam->IsPFon() ){ | |
1327 | ep = fFeeParam->GetPFeffectPedestal(); | |
1328 | } | |
1329 | Int_t nexp = fFeeParam->GetTFnExp(); | |
1330 | Int_t *isource = new Int_t[fNTimeBin]; | |
1331 | filPed.Init(10,2); | |
1332 | filPed.AssignInt(ep); | |
1333 | Int_t epf = filPed.GetValue(); | |
1334 | for(Int_t iT = 0; iT < fNTimeBin; iT++){ | |
1335 | isource[iT] = ep; | |
1336 | ieffped[iT] = epf; | |
1337 | } | |
1338 | ||
1339 | if( fFeeParam->IsTFon() ) { | |
1340 | FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp); | |
1341 | } | |
1342 | ||
1343 | delete[] isource; | |
1344 | } | |
1345 | ||
1346 | //the following values should be in AliTRDfeeParam and have to be read in properly | |
1347 | //naming follows conventions in TRAP-manual | |
1348 | ||
1349 | ||
1350 | Bool_t bVBY = kTRUE; // cluster-verification bypass | |
1351 | ||
1352 | Double_t cQTParam = 0; // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10 | |
1353 | AliTRDtrapAlu cQTAlu; | |
1354 | cQTAlu.Init(1,10,0,63); | |
1355 | cQTAlu.AssignDouble(cQTParam); | |
1356 | Int_t cQT = cQTAlu.GetValue(); | |
1357 | ||
1358 | // linear fit | |
1359 | Int_t tFS = fFeeParam->GetLinearFitStart(); // linear fit start | |
1360 | Int_t tFE = fFeeParam->GetLinearFitEnd(); // linear fit stop | |
1361 | ||
1362 | // charge accumulators | |
1363 | Int_t tQS0 = fFeeParam->GetQacc0Start(); // start-time for charge-accumulator 0 | |
1364 | Int_t tQE0 = fFeeParam->GetQacc0End(); // stop-time for charge-accumulator 0 | |
1365 | Int_t tQS1 = fFeeParam->GetQacc1Start(); // start-time for charge-accumulator 1 | |
1366 | Int_t tQE1 = fFeeParam->GetQacc1End(); // stop-time for charge-accumulator 1 | |
1367 | // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N | |
1368 | ||
1369 | Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold | |
1370 | AliTRDtrapAlu cTHAlu; | |
1371 | cTHAlu.Init(12,2); | |
1372 | cTHAlu.AssignDouble(cTHParam); | |
1373 | Int_t cTH = cTHAlu.GetValue(); // cTH used for comparison | |
1374 | ||
1375 | struct List_t { | |
1376 | List_t *next; | |
1377 | Int_t iadc; | |
1378 | Int_t value; | |
1379 | }; | |
1380 | ||
1381 | List_t selection[7]; // list with 7 elements | |
1382 | List_t *list = NULL; | |
1383 | List_t *listLeft = NULL; | |
1384 | ||
1385 | Int_t* qsum = new Int_t[fNADC]; | |
1386 | ||
1387 | // fit sums | |
1388 | AliTRDtrapAlu qsumAlu; | |
1389 | qsumAlu.Init(12,2); // charge sum will be 12+2 bits | |
1390 | AliTRDtrapAlu dCOGAlu; | |
1391 | dCOGAlu.Init(1,7,0,127); // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT | |
1392 | AliTRDtrapAlu yrawAlu; | |
1393 | yrawAlu.Init(1,8,-1,255); | |
1394 | AliTRDtrapAlu yAlu; | |
1395 | yAlu.Init(1,16,-1,0xFF00); // only first 8 past-comma bits filled;additional 8 bits for accuracy;maximum 1 - 2^-8; sign is given by + or - | |
1396 | AliTRDtrapAlu xAlu; | |
1397 | xAlu.Init(5,8); // 8 past-comma bits because value will be added/multiplied to another value with this accuracy | |
1398 | AliTRDtrapAlu xxAlu; | |
1399 | xxAlu.Init(10,0); | |
1400 | AliTRDtrapAlu yyAlu; | |
1401 | yyAlu.Init(1,16,0,0xFFFF); // maximum is 2^16-1; 16Bit for past-commas | |
1402 | AliTRDtrapAlu xyAlu; | |
1403 | xyAlu.Init(6,8); | |
1404 | AliTRDtrapAlu XAlu; | |
1405 | XAlu.Init(9,0); | |
1406 | AliTRDtrapAlu XXAlu; | |
1407 | XXAlu.Init(14,0); | |
1408 | AliTRDtrapAlu YAlu; | |
1409 | YAlu.Init(5,8); // 14 bit, 1 is sign-bit; therefore only 13 bit | |
1410 | AliTRDtrapAlu YYAlu; | |
1411 | YYAlu.Init(5,16); | |
1412 | AliTRDtrapAlu XYAlu; | |
1413 | XYAlu.Init(8,8); // 17 bit, 1 is sign-bit; therefore only 16 bit | |
1414 | AliTRDtrapAlu qtruncAlu; | |
1415 | qtruncAlu.Init(12,0); | |
1416 | AliTRDtrapAlu QT0Alu; | |
1417 | QT0Alu.Init(15,0); | |
1418 | AliTRDtrapAlu QT1Alu; | |
1419 | QT1Alu.Init(16,0); | |
1420 | ||
1421 | ||
1422 | AliTRDtrapAlu inverseNAlu; | |
1423 | inverseNAlu.Init(1,8); // replaces the LUT for 1/N | |
1424 | AliTRDtrapAlu MeanChargeAlu; // mean charge in ADC counts | |
1425 | MeanChargeAlu.Init(8,0); | |
1426 | AliTRDtrapAlu TotalChargeAlu; | |
1427 | TotalChargeAlu.Init(17,8); | |
1428 | //nr of post comma bits should be the same for inverseN and TotalCharge | |
1429 | ||
1430 | ||
1431 | SetPosLUT(); // initialize the position correction LUT for this MCM; | |
1432 | ||
1433 | ||
1434 | // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18 | |
1435 | Int_t *X = new Int_t[fNADC-2]; | |
1436 | Int_t *XX = new Int_t[fNADC-2]; | |
1437 | Int_t *Y = new Int_t[fNADC-2]; | |
1438 | Int_t *YY = new Int_t[fNADC-2]; | |
1439 | Int_t *XY = new Int_t[fNADC-2]; | |
1440 | Int_t *N = new Int_t[fNADC-2]; | |
1441 | Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge | |
1442 | Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge | |
1443 | ||
1444 | for (Int_t iCol = 0; iCol < fNADC-2; iCol++) { | |
1445 | ||
1446 | // initialize fit-sums | |
1447 | X[iCol] = 0; | |
1448 | XX[iCol] = 0; | |
1449 | Y[iCol] = 0; | |
1450 | YY[iCol] = 0; | |
1451 | XY[iCol] = 0; | |
1452 | N[iCol] = 0; | |
1453 | QT0[iCol] = 0; | |
1454 | QT1[iCol] = 0; | |
1455 | } | |
1456 | ||
1457 | ||
1458 | filPed.Init(7,2); // convert filtered pedestal into 7+2Bits | |
1459 | ||
1460 | for(Int_t iT = 0; iT < fNTimeBin; iT++){ | |
1461 | ||
1462 | if(iT<tFS || iT>=tFE) continue; // linear fit yes/no? // !!**enable**!! | |
1463 | ||
1464 | // reset | |
1465 | Int_t portChannel[4] = {-1,-1,-1,-1}; | |
1466 | Int_t clusterCharge[4] = {0,0,0,0}; | |
1467 | Int_t leftCharge[4] = {0,0,0,0}; | |
1468 | Int_t centerCharge[4] = {0,0,0,0}; | |
1469 | Int_t rightCharge[4] = {0,0,0,0}; | |
1470 | ||
1471 | Int_t mark = 0; | |
1472 | ||
1473 | filPed.AssignFormatted(ieffped[iT]); // no size-checking with AssignFormatted; ieffped>=0 | |
1474 | filPed = filPed; // this checks the size | |
1475 | ||
1476 | ieffped[iT] = filPed.GetValue(); | |
1477 | ||
1478 | for(Int_t i = 0; i<7; i++){ | |
1479 | selection[i].next = NULL; | |
1480 | selection[i].iadc = -1; // value of -1: invalid adc | |
1481 | selection[i].value = 0; | |
1482 | ||
1483 | } | |
1484 | // selection[0] is starting list-element; just for pointing | |
1485 | ||
1486 | // loop over inner adc's | |
1487 | for (Int_t iCol = 1; iCol < fNADC-1; iCol++) { | |
1488 | ||
1489 | Int_t left = fADCT[iCol-1][iT]; | |
1490 | Int_t center = fADCT[iCol][iT]; | |
1491 | Int_t right = fADCT[iCol+1][iT]; | |
1492 | ||
1493 | Int_t sum = left + center + right; // cluster charge sum | |
1494 | qsumAlu.AssignFormatted(sum); | |
1495 | qsumAlu = qsumAlu; // size-checking; redundant | |
1496 | ||
1497 | qsum[iCol] = qsumAlu.GetValue(); | |
1498 | ||
1499 | //hit detection and masking | |
1500 | if(center>=left){ | |
1501 | if(center>right){ | |
1502 | if(qsum[iCol]>=(cTH + 3*ieffped[iT])){ // effective pedestal of all three channels must be added to cTH(+20); this is not parallel to TRAP manual; maybe cTH has to be adjusted in fFeeParam; therefore channels are not yet reduced by their pedestal | |
1503 | mark |= 1; // marker | |
1504 | } | |
1505 | } | |
1506 | } | |
1507 | mark = mark<<1; | |
1508 | } | |
1509 | mark = mark>>1; | |
1510 | ||
1511 | ||
1512 | // get selection of 6 adc's and sort,starting with greatest values | |
1513 | ||
1514 | //read three from right side and sort (primitive sorting algorithm) | |
1515 | Int_t i = 0; // adc number | |
1516 | Int_t j = 1; // selection number | |
1517 | while(i<fNADC-2 && j<=3){ | |
1518 | i = i + 1; | |
1519 | if((mark>>(i-1)) & 1 == 1) { | |
1520 | selection[j].iadc = fNADC-1-i; | |
1521 | selection[j].value = qsum[fNADC-1-i]>>6; // for hit-selection only the first 8 out of the 14 Bits are used for comparison | |
1522 | ||
1523 | // insert into sorted list | |
1524 | listLeft = &selection[0]; | |
1525 | list = listLeft->next; | |
1526 | ||
1527 | if(list!=NULL) { | |
1528 | while((list->next != NULL) && (selection[j].value <= list->value)){ | |
1529 | listLeft = list; | |
1530 | list = list->next; | |
1531 | } | |
1532 | ||
1533 | if(selection[j].value<=list->value){ | |
1534 | selection[j].next = list->next; | |
1535 | list->next = &selection[j]; | |
1536 | } | |
1537 | else { | |
1538 | listLeft->next = &selection[j]; | |
1539 | selection[j].next = list; | |
1540 | } | |
1541 | } | |
1542 | else{ | |
1543 | listLeft->next = &selection[j]; | |
1544 | selection[j].next = list; | |
1545 | } | |
1546 | ||
1547 | j = j + 1; | |
1548 | } | |
1549 | } | |
1550 | ||
1551 | ||
1552 | // read three from left side | |
1553 | Int_t k = fNADC-2; | |
1554 | while(k>i && j<=6) { | |
1555 | if((mark>>(k-1)) & 1 == 1) { | |
1556 | selection[j].iadc = fNADC-1-k; | |
1557 | selection[j].value = qsum[fNADC-1-k]>>6; | |
1558 | ||
1559 | listLeft = &selection[0]; | |
1560 | list = listLeft->next; | |
1561 | ||
1562 | if(list!=NULL){ | |
1563 | while((list->next != NULL) && (selection[j].value <= list->value)){ | |
1564 | listLeft = list; | |
1565 | list = list->next; | |
1566 | } | |
1567 | ||
1568 | if(selection[j].value<=list->value){ | |
1569 | selection[j].next = list->next; | |
1570 | list->next = &selection[j]; | |
1571 | } | |
1572 | else { | |
1573 | listLeft->next = &selection[j]; | |
1574 | selection[j].next = list; | |
1575 | } | |
1576 | } | |
1577 | else{ | |
1578 | listLeft->next = &selection[j]; | |
1579 | selection[j].next = list; | |
1580 | } | |
1581 | ||
1582 | j = j + 1; | |
1583 | } | |
1584 | k = k - 1; | |
1585 | } | |
1586 | ||
1587 | // get the four with greatest charge-sum | |
1588 | list = &selection[0]; | |
1589 | for(i = 0; i<4; i++){ | |
1590 | if(list->next == NULL) continue; | |
1591 | list = list->next; | |
1592 | if(list->iadc == -1) continue; | |
1593 | Int_t adc = list->iadc; // channel number with selected hit | |
1594 | ||
1595 | // the following arrays contain the four chosen channels in 1 time-bin | |
1596 | portChannel[i] = adc; | |
1597 | clusterCharge[i] = qsum[adc]; | |
1598 | leftCharge[i] = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal) | |
1599 | centerCharge[i] = fADCT[adc][iT] - ieffped[iT]; | |
1600 | rightCharge[i] = fADCT[adc+1][iT] - ieffped[iT]; | |
1601 | } | |
1602 | ||
1603 | // arithmetic unit | |
1604 | ||
1605 | // cluster verification | |
1606 | if(!bVBY){ | |
1607 | for(i = 0; i<4; i++){ | |
1608 | Int_t lr = leftCharge[i]*rightCharge[i]*1024; | |
1609 | Int_t cc = centerCharge[i]*centerCharge[i]*cQT; | |
1610 | if (lr>=cc){ | |
1611 | portChannel[i] = -1; // set to invalid address | |
1612 | clusterCharge[i] = 0; | |
1613 | } | |
1614 | } | |
1615 | } | |
1616 | ||
1617 | // fit-sums of valid channels | |
1618 | // local hit position | |
1619 | for(i = 0; i<4; i++){ | |
1620 | if (centerCharge[i] == 0) { | |
1621 | portChannel[i] = -1; | |
1622 | }// prevent division by 0 | |
1623 | ||
1624 | if (portChannel[i] == -1) continue; | |
1625 | ||
1626 | Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i]; | |
1627 | ||
1628 | Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE; | |
1629 | dCOG = (sign == kFALSE) ? dCOG : -dCOG; // AssignDouble doesn't allow for signed doubles | |
1630 | dCOGAlu.AssignDouble(dCOG); | |
1631 | Int_t iLUTpos = dCOGAlu.GetValue(); // steers position in LUT | |
1632 | ||
1633 | dCOG = dCOG/2; | |
1634 | yrawAlu.AssignDouble(dCOG); | |
1635 | Int_t iCOG = yrawAlu.GetValue(); | |
1636 | Int_t y = iCOG + fPosLUT[iLUTpos % 128]; // local position in pad-units | |
1637 | yrawAlu.AssignFormatted(y); // 0<y<1 | |
1638 | yAlu = yrawAlu; // convert to 16 past-comma bits | |
1639 | ||
1640 | if(sign == kTRUE) yAlu.SetSign(-1); // buffer width of 9 bits; sign on real (not estimated) position | |
1641 | xAlu.AssignInt(iT); // buffer width of 5 bits | |
1642 | ||
1643 | ||
1644 | xxAlu = xAlu * xAlu; // buffer width of 10 bits -> fulfilled by x*x | |
1645 | ||
1646 | yyAlu = yAlu * yAlu; // buffer width of 16 bits | |
1647 | ||
1648 | xyAlu = xAlu * yAlu; // buffer width of 14 bits | |
1649 | ||
1650 | Int_t adc = portChannel[i]-1; // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18) | |
1651 | ||
1652 | // calculate fit-sums recursively | |
1653 | // interpretation of their bit-length is given as comment | |
1654 | ||
1655 | // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value | |
1656 | ||
1657 | XAlu.AssignFormatted(X[adc]); | |
1658 | XAlu = XAlu + xAlu; // buffer width of 9 bits | |
1659 | X[adc] = XAlu.GetValue(); | |
1660 | ||
1661 | XXAlu.AssignFormatted(XX[adc]); | |
1662 | XXAlu = XXAlu + xxAlu; // buffer width of 14 bits | |
1663 | XX[adc] = XXAlu.GetValue(); | |
1664 | ||
1665 | if (Y[adc] < 0) { | |
1666 | YAlu.AssignFormatted(-Y[adc]); // make sure that only positive values are assigned; sign-setting must be done by hand | |
1667 | YAlu.SetSign(-1); | |
1668 | } | |
1669 | else { | |
1670 | YAlu.AssignFormatted(Y[adc]); | |
1671 | YAlu.SetSign(1); | |
1672 | } | |
1673 | ||
1674 | YAlu = YAlu + yAlu; // buffer width of 14 bits (8 past-comma); | |
1675 | Y[adc] = YAlu.GetSignedValue(); | |
1676 | ||
1677 | YYAlu.AssignFormatted(YY[adc]); | |
1678 | YYAlu = YYAlu + yyAlu; // buffer width of 21 bits (16 past-comma) | |
1679 | YY[adc] = YYAlu.GetValue(); | |
1680 | ||
1681 | if (XY[adc] < 0) { | |
1682 | XYAlu.AssignFormatted(-XY[adc]); | |
1683 | XYAlu.SetSign(-1); | |
1684 | } | |
1685 | else { | |
1686 | XYAlu.AssignFormatted(XY[adc]); | |
1687 | XYAlu.SetSign(1); | |
1688 | } | |
1689 | ||
1690 | XYAlu = XYAlu + xyAlu; // buffer allows 17 bits (8 past-comma) | |
1691 | XY[adc] = XYAlu.GetSignedValue(); | |
1692 | ||
1693 | N[adc] = N[adc] + 1; | |
1694 | ||
1695 | ||
1696 | // accumulated charge | |
1697 | qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped! | |
1698 | qtruncAlu = qsumAlu; | |
1699 | ||
1700 | if(iT>=tQS0 && iT<=tQE0){ | |
1701 | QT0Alu.AssignFormatted(QT0[adc]); | |
1702 | QT0Alu = QT0Alu + qtruncAlu; | |
1703 | QT0[adc] = QT0Alu.GetValue(); | |
1704 | //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit | |
1705 | } | |
1706 | ||
1707 | if(iT>=tQS1 && iT<=tQE1){ | |
1708 | QT1Alu.AssignFormatted(QT1[adc]); | |
1709 | QT1Alu = QT1Alu + qtruncAlu; | |
1710 | QT1[adc] = QT1Alu.GetValue(); | |
1711 | //interpretation of QT1 as 12bit-value; buffer allows 16 Bit | |
1712 | } | |
1713 | }// i | |
1714 | ||
1715 | // remapping is done!! | |
1716 | ||
1717 | }//iT | |
1718 | ||
1719 | ||
1720 | ||
1721 | // tracklet-assembly | |
1722 | ||
1723 | // put into AliTRDfeeParam and take care that values are in proper range | |
1724 | const Int_t cTCL = 1; // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??) | |
1725 | const Int_t cTCT = 8; // joint number of hits; 8<=TCT<=31 | |
1726 | ||
1727 | Int_t mPair = 0; // marker for possible tracklet pairs | |
1728 | Int_t* hitSum = new Int_t[fNADC-3]; | |
1729 | // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18; | |
1730 | ||
1731 | // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs | |
1732 | for (Int_t iCol = 0; iCol < fNADC-3; iCol++) { | |
1733 | hitSum[iCol] = N[iCol] + N[iCol+1]; | |
1734 | if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) { | |
1735 | mPair |= 1; // mark as possible channel-pair | |
1736 | ||
1737 | } | |
1738 | mPair = mPair<<1; | |
1739 | } | |
1740 | mPair = mPair>>1; | |
1741 | ||
1742 | List_t* selectPair = new List_t[fNADC-2]; // list with 18 elements (0..18) containing the left channel-nr and hit sums | |
1743 | // selectPair[18] is starting list-element just for pointing | |
1744 | for(Int_t k = 0; k<fNADC-2; k++){ | |
1745 | selectPair[k].next = NULL; | |
1746 | selectPair[k].iadc = -1; // invalid adc | |
1747 | selectPair[k].value = 0; | |
1748 | ||
1749 | } | |
1750 | ||
1751 | list = NULL; | |
1752 | listLeft = NULL; | |
1753 | ||
1754 | // read marker and sort according to hit-sum | |
1755 | ||
1756 | Int_t adcL = 0; // left adc-channel-number (remapped) | |
1757 | Int_t selNr = 0; // current number in list | |
1758 | ||
1759 | // insert marked channels into list and sort according to hit-sum | |
1760 | while(adcL < fNADC-3 && selNr < fNADC-3){ | |
1761 | ||
1762 | if((mPair>>((fNADC-4)-(adcL))) & 1 == 1) { | |
1763 | selectPair[selNr].iadc = adcL; | |
1764 | selectPair[selNr].value = hitSum[adcL]; | |
1765 | ||
1766 | listLeft = &selectPair[fNADC-3]; | |
1767 | list = listLeft->next; | |
1768 | ||
1769 | if(list!=NULL) { | |
1770 | while((list->next != NULL) && (selectPair[selNr].value <= list->value)){ | |
1771 | listLeft = list; | |
1772 | list = list->next; | |
1773 | } | |
1774 | ||
1775 | if(selectPair[selNr].value <= list->value){ | |
1776 | selectPair[selNr].next = list->next; | |
1777 | list->next = &selectPair[selNr]; | |
1778 | } | |
1779 | else { | |
1780 | listLeft->next = &selectPair[selNr]; | |
1781 | selectPair[selNr].next = list; | |
1782 | } | |
1783 | ||
1784 | } | |
1785 | else{ | |
1786 | listLeft->next = &selectPair[selNr]; | |
1787 | selectPair[selNr].next = list; | |
1788 | } | |
1789 | ||
1790 | selNr = selNr + 1; | |
1791 | } | |
1792 | adcL = adcL + 1; | |
1793 | } | |
1794 | ||
1795 | //select up to 4 channels with maximum number of hits | |
1796 | Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum | |
1797 | Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets | |
1798 | list = &selectPair[fNADC-3]; | |
1799 | ||
1800 | for (Int_t i = 0; i<4; i++) { | |
1801 | if(list->next == NULL) continue; | |
1802 | list = list->next; | |
1803 | if(list->iadc == -1) continue; | |
1804 | lpairChannel[i] = list->iadc; // channel number with selected hit | |
1805 | rpairChannel[i] = lpairChannel[i]+1; | |
1806 | } | |
1807 | ||
1808 | // avoid submission of double tracklets | |
1809 | for (Int_t i = 3; i>0; i--) { | |
1810 | for (Int_t j = i-1; j>-1; j--) { | |
1811 | if(lpairChannel[i] == rpairChannel[j]) { | |
1812 | lpairChannel[i] = -1; | |
1813 | rpairChannel[i] = -1; | |
1814 | break; | |
1815 | } | |
1816 | if(rpairChannel[i] == lpairChannel[j]) { | |
1817 | lpairChannel[i] = -1; | |
1818 | rpairChannel[i] = -1; | |
1819 | break; | |
1820 | } | |
1821 | } | |
1822 | } | |
1823 | ||
1824 | // merging of the fit-sums of the remainig channels | |
1825 | // assume same data-word-width as for fit-sums for 1 channel | |
1826 | // relative scales! | |
1827 | Int_t mADC[4]; | |
1828 | Int_t mN[4]; | |
1829 | Int_t mQT0[4]; | |
1830 | Int_t mQT1[4]; | |
1831 | Int_t mX[4]; | |
1832 | Int_t mXX[4]; | |
1833 | Int_t mY[4]; | |
1834 | Int_t mYY[4]; | |
1835 | Int_t mXY[4]; | |
1836 | Int_t mOffset[4]; | |
1837 | Int_t mSlope[4]; | |
1838 | Int_t mMeanCharge[4]; | |
1839 | Int_t inverseN = 0; | |
1840 | Double_t invN = 0; | |
1841 | ||
1842 | for (Int_t i = 0; i<4; i++){ | |
1843 | mADC[i] = -1; // set to invalid number | |
1844 | mN[i] = 0; | |
1845 | mQT0[i] = 0; | |
1846 | mQT1[i] = 0; | |
1847 | mX[i] = 0; | |
1848 | mXX[i] = 0; | |
1849 | mY[i] = 0; | |
1850 | mYY[i] = 0; | |
1851 | mXY[i] = 0; | |
1852 | mOffset[i] = 0; | |
1853 | mSlope[i] = 0; | |
1854 | mMeanCharge[i] = 0; | |
1855 | } | |
1856 | ||
1857 | ||
1858 | YAlu.AssignInt(1); | |
1859 | Int_t wpad = YAlu.GetValue(); // 1 with 8 past-comma bits | |
1860 | ||
1861 | for (Int_t i = 0; i<4; i++){ | |
1862 | ||
1863 | mADC[i] = lpairChannel[i]; // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17) | |
1864 | // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering | |
1865 | Int_t madc = mADC[i]; | |
1866 | if (madc == -1) continue; | |
1867 | mN[i] = hitSum[madc]; | |
1868 | ||
1869 | // don't merge fit sums in case of a stand-alone tracklet (consisting of only 1 channel); in that case only left channel makes up the fit sums | |
1870 | if (N[madc+1] == 0) { | |
1871 | mQT0[i] = QT0[madc]; | |
1872 | mQT1[i] = QT1[madc]; | |
1873 | ||
1874 | } | |
1875 | else { | |
1876 | ||
1877 | // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums? | |
1878 | ||
1879 | mQT0[i] = QT0[madc] + QT0[madc+1]; | |
1880 | QT0Alu.AssignFormatted(mQT0[i]); | |
1881 | QT0Alu = QT0Alu; // size-check | |
1882 | mQT0[i] = QT0Alu.GetValue(); // write back | |
1883 | ||
1884 | mQT1[i] = QT1[madc] + QT1[madc+1]; | |
1885 | QT1Alu.AssignFormatted(mQT1[i]); | |
1886 | QT1Alu = QT1Alu; | |
1887 | mQT1[i] = QT1Alu.GetValue(); | |
1888 | } | |
1889 | ||
1890 | // calculate the mean charge in adc values; later to be replaced by electron likelihood | |
1891 | mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge | |
1892 | mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge | |
1893 | // simulate LUT for 1/N; LUT is fed with the double-accurate pre-calculated value of 1/N; accuracy of entries has to be adjusted to real TRAP | |
1894 | invN = 1.0/(mN[i]); | |
1895 | inverseNAlu.AssignDouble(invN); | |
1896 | inverseN = inverseNAlu.GetValue(); | |
1897 | mMeanCharge[i] = mMeanCharge[i] * inverseN; // now to be interpreted with 8 past-comma bits | |
1898 | TotalChargeAlu.AssignFormatted(mMeanCharge[i]); | |
1899 | TotalChargeAlu = TotalChargeAlu; | |
1900 | MeanChargeAlu = TotalChargeAlu; | |
1901 | mMeanCharge[i] = MeanChargeAlu.GetValue(); | |
1902 | ||
1903 | if (N[madc+1] == 0) { | |
1904 | mX[i] = X[madc]; | |
1905 | mXX[i] = XX[madc]; | |
1906 | mY[i] = Y[madc]; | |
1907 | mXY[i] = XY[madc]; | |
1908 | mYY[i] = YY[madc]; | |
1909 | } | |
1910 | else { | |
1911 | ||
1912 | mX[i] = X[madc] + X[madc+1]; | |
1913 | XAlu.AssignFormatted(mX[i]); | |
1914 | XAlu = XAlu; | |
1915 | mX[i] = XAlu.GetValue(); | |
1916 | ||
1917 | mXX[i] = XX[madc] + XX[madc+1]; | |
1918 | XXAlu.AssignFormatted(mXX[i]); | |
1919 | XXAlu = XXAlu; | |
1920 | mXX[i] = XXAlu.GetValue(); | |
1921 | ||
1922 | ||
1923 | mY[i] = Y[madc] + Y[madc+1] + wpad; | |
1924 | if (mY[i] < 0) { | |
1925 | YAlu.AssignFormatted(-mY[i]); | |
1926 | YAlu.SetSign(-1); | |
1927 | } | |
1928 | else { | |
1929 | YAlu.AssignFormatted(mY[i]); | |
1930 | YAlu.SetSign(1); | |
1931 | } | |
1932 | YAlu = YAlu; | |
1933 | mY[i] = YAlu.GetSignedValue(); | |
1934 | ||
1935 | mXY[i] = XY[madc] + XY[madc+1] + X[madc+1]*wpad; | |
1936 | ||
1937 | if (mXY[i] < 0) { | |
1938 | XYAlu.AssignFormatted(-mXY[i]); | |
1939 | XYAlu.SetSign(-1); | |
1940 | } | |
1941 | else { | |
1942 | XYAlu.AssignFormatted(mXY[i]); | |
1943 | XYAlu.SetSign(1); | |
1944 | } | |
1945 | XYAlu = XYAlu; | |
1946 | mXY[i] = XYAlu.GetSignedValue(); | |
1947 | ||
1948 | mYY[i] = YY[madc] + YY[madc+1] + 2*Y[madc+1]*wpad + wpad*wpad; | |
1949 | if (mYY[i] < 0) { | |
1950 | YYAlu.AssignFormatted(-mYY[i]); | |
1951 | YYAlu.SetSign(-1); | |
1952 | } | |
1953 | else { | |
1954 | YYAlu.AssignFormatted(mYY[i]); | |
1955 | YYAlu.SetSign(1); | |
1956 | } | |
1957 | ||
1958 | YYAlu = YYAlu; | |
1959 | mYY[i] = YYAlu.GetSignedValue(); | |
1960 | } | |
1961 | ||
1962 | } | |
1963 | ||
1964 | // calculation of offset and slope from the merged fit-sums; | |
1965 | // YY is needed for some error measure only; still to be done | |
1966 | // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale | |
1967 | ||
1968 | // which formats should be chosen? | |
1969 | AliTRDtrapAlu denomAlu; | |
1970 | denomAlu.Init(20,8); | |
1971 | AliTRDtrapAlu numAlu; | |
1972 | numAlu.Init(20,8); | |
1973 | // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset | |
1974 | // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big | |
1975 | ||
1976 | for (Int_t i = 0; i<4; i++) { | |
1977 | if (mADC[i] == -1) continue; | |
1978 | ||
1979 | Int_t num0 = (mN[i]*mXX[i]-mX[i]*mX[i]); | |
1980 | if (num0 < 0) { | |
1981 | denomAlu.AssignInt(-num0); // num0 does not have to be interpreted as having past-comma bits -> AssignInt | |
1982 | denomAlu.SetSign(-1); | |
1983 | } | |
1984 | else { | |
1985 | denomAlu.AssignInt(num0); | |
1986 | denomAlu.SetSign(1); | |
1987 | } | |
1988 | ||
1989 | Int_t num1 = mN[i]*mXY[i] - mX[i]*mY[i]; | |
1990 | if (num1 < 0) { | |
1991 | numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits | |
1992 | numAlu.SetSign(-1); | |
1993 | } | |
1994 | else { | |
1995 | numAlu.AssignFormatted(num1); | |
1996 | numAlu.SetSign(1); | |
1997 | } | |
1998 | numAlu = numAlu/denomAlu; | |
1999 | mSlope[i] = numAlu.GetSignedValue(); | |
2000 | ||
2001 | Int_t num2 = mXX[i]*mY[i] - mX[i]*mXY[i]; | |
2002 | ||
2003 | if (num2 < 0) { | |
2004 | numAlu.AssignFormatted(-num2); | |
2005 | numAlu.SetSign(-1); | |
2006 | } | |
2007 | else { | |
2008 | numAlu.AssignFormatted(num2); | |
2009 | numAlu.SetSign(1); | |
2010 | } | |
2011 | ||
2012 | numAlu = numAlu/denomAlu; | |
2013 | ||
2014 | ||
2015 | mOffset[i] = numAlu.GetSignedValue(); | |
2016 | numAlu.SetSign(1); | |
2017 | denomAlu.SetSign(1); | |
2018 | ||
2019 | ||
2020 | //numAlu.AssignInt(mADC[i]+1); // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away) | |
2021 | numAlu.AssignDouble((Double_t)mADC[i] + 1.5); // numAlu has enough pre-comma place for that; correct trafo, best values | |
2022 | mOffset[i] = mOffset[i] + numAlu.GetValue(); // transform offset to a coord.system relative to chip; +1 to avoid neg. values | |
2023 | ||
2024 | // up to here: adc-mapping according to TRAP and in line with pad-col mapping | |
2025 | // reverse adc-counting to be again in line with the online mapping | |
2026 | mADC[i] = fNADC - 4 - mADC[i]; // fNADC-4-mADC[i]: 0..17; remapping necessary; | |
2027 | mADC[i] = mADC[i] + 2; | |
2028 | // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19 | |
2029 | } | |
2030 | ||
2031 | // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected; | |
2032 | // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction); | |
2033 | ||
2034 | // transform parameters to the local coordinate-system of a stack (used by GTU) | |
2035 | AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack); | |
2036 | ||
2037 | Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm | |
2038 | //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP?? | |
2039 | ||
2040 | // difference between width of inner and outer pads of a row is not accounted for; | |
2041 | ||
2042 | Double_t magField = 0.4; // z-component of magnetic field in Tesla; adjust to current simulation!!; magnetic field can hardly be evaluated for the position of each mcm | |
2043 | Double_t eCharge = 0.3; // unit charge in (GeV/c)/m*T | |
2044 | Double_t ptMin = 2.3; // minimum transverse momentum (GeV/c); to be adjusted(?) | |
2045 | ||
2046 | Double_t granularityOffset = 0.160; // granularity for offset in mm | |
2047 | Double_t granularitySlope = 0.140; // granularity for slope in mm | |
2048 | ||
2049 | // get the coordinates in SM-system; parameters: | |
2050 | ||
2051 | Double_t zPos = (padPlane->GetRowPos(fRow))*10.0; // z-position of the MCM; fRow is counted on a chamber; SM consists of 5 | |
2052 | // zPos is position of pad-borders; | |
2053 | Double_t zOffset = 0.0; | |
2054 | if ( fRow == 0 || fRow == 15 ) { | |
2055 | zOffset = padPlane->GetLengthOPad(); | |
2056 | } | |
2057 | else { | |
2058 | zOffset = padPlane->GetLengthIPad(); | |
2059 | } | |
2060 | zOffset = (-1.0) * zOffset/2.0; | |
2061 | // turn zPos to be z-coordinate at middle of pad-row | |
2062 | zPos = zPos + zOffset; | |
2063 | ||
2064 | ||
2065 | Double_t xPos = 0.0; // x-position of the upper border of the drift-chamber of actual layer | |
2066 | Int_t icol = 0; // column-number of adc-channel | |
2067 | Double_t yPos[4]; // y-position of the pad to which ADC is connected | |
2068 | Double_t dx = 30.0; // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry | |
2069 | Double_t vdrift = fCal->GetVdriftAverage(fChaId); // averaged drift velocity for this detector | |
2070 | Double_t lorTan = fCal->GetOmegaTau(vdrift,magField); // tan of the Lorentz-angle for this detector; could be evaluated and set as a parameter for each mcm | |
2071 | //Double_t lorAngle = 7.0; // Lorentz-angle in degrees | |
2072 | Double_t tiltAngle = padPlane->GetTiltingAngle(); // sign-respecting tilting angle of pads in actual layer | |
2073 | Double_t tiltTan = TMath::Tan(TMath::Pi()/180.0 * tiltAngle); | |
2074 | //Double_t lorTan = TMath::Tan(TMath::Pi()/180.0 * lorAngle); | |
2075 | ||
2076 | Double_t alphaMax[4]; // maximum deflection from the direction to the primary vertex; granularity of hit pads | |
2077 | Double_t slopeMin[4]; // local limits for the deflection | |
2078 | Double_t slopeMax[4]; | |
2079 | Int_t mslopeMin[4]; // in granularity units; to be compared to mSlope[i] | |
2080 | Int_t mslopeMax[4]; | |
2081 | ||
2082 | ||
2083 | //x coord. of upper side of drift chambers in local SM-system (in mm) | |
2084 | //obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm); | |
2085 | //the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane; | |
2086 | switch(fLayer) | |
2087 | { | |
2088 | case 0: | |
2089 | xPos = 3003.0; | |
2090 | break; | |
2091 | case 1: | |
2092 | xPos = 3129.0; | |
2093 | break; | |
2094 | case 2: | |
2095 | xPos = 3255.0; | |
2096 | break; | |
2097 | case 3: | |
2098 | xPos = 3381.0; | |
2099 | break; | |
2100 | case 4: | |
2101 | xPos = 3507.0; | |
2102 | break; | |
2103 | case 5: | |
2104 | xPos = 3633.0; | |
2105 | break; | |
2106 | } | |
2107 | ||
2108 | // calculation of offset-correction n: | |
2109 | ||
2110 | Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4)); | |
2111 | ||
2112 | nCorrectOffset = (nCorrectOffset - 4)*18 - 1; | |
2113 | if (nCorrectOffset < 0) { | |
2114 | numAlu.AssignInt(-nCorrectOffset); | |
2115 | numAlu.SetSign(-1); | |
2116 | } | |
2117 | else { | |
2118 | numAlu.AssignInt(nCorrectOffset); | |
2119 | numAlu.SetSign(1); | |
2120 | } | |
2121 | nCorrectOffset = numAlu.GetSignedValue(); | |
2122 | Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0 | |
2123 | ||
2124 | // calculation of slope-correction | |
2125 | ||
2126 | // this is only true for tracks coming (approx.) from primary vertex | |
2127 | Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope; | |
2128 | // Double_t cCorrectSlope = zPos/xPos*dx*tiltTan/granularitySlope; | |
2129 | // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0 | |
2130 | ||
2131 | if (cCorrectSlope < 0) { | |
2132 | numAlu.AssignDouble(-cCorrectSlope); | |
2133 | numAlu.SetSign(-1); | |
2134 | } | |
2135 | else { | |
2136 | numAlu.AssignDouble(cCorrectSlope); | |
2137 | numAlu.SetSign(1); | |
2138 | } | |
2139 | cCorrectSlope = numAlu.GetSignedValue(); | |
2140 | ||
2141 | // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins) | |
2142 | // different pad-width of outer pads of a pad-plane not taken into account | |
2143 | // tFS: upper plane of drift-volume (not amplification region); this choice is important for offset | |
2144 | // tFE: does !!not!! need to correspond to lower plane of drift-volume; TR hits can be cut; | |
2145 | ||
2146 | Double_t mCorrectSlope = ((Double_t)(fNTimeBin-tFS))*padWidthI/granularitySlope; // >= 0.0 | |
2147 | ||
2148 | AliTRDtrapAlu correctAlu; | |
2149 | correctAlu.Init(20,8); | |
2150 | ||
2151 | AliTRDtrapAlu offsetAlu; | |
2152 | offsetAlu.Init(13,0,-0x1000,0x0FFF); // 13 bit-word; 2-complement (1 sign-bit); asymmetric range | |
2153 | ||
2154 | AliTRDtrapAlu slopeAlu; | |
2155 | slopeAlu.Init(7,0,-0x40,0x3F); // 7 bit-word; 2-complement (1 sign-bit); | |
2156 | ||
2157 | for (Int_t i = 0; i<4; i++) { | |
2158 | ||
2159 | if (mADC[i] == -1) continue; | |
2160 | ||
2161 | icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping | |
2162 | yPos[i] = (padPlane->GetColPos(icol))*10.0; | |
2163 | ||
2164 | ||
2165 | // offset: | |
2166 | ||
2167 | correctAlu.AssignDouble(mCorrectOffset); // done because max. accuracy is 8 bit | |
2168 | mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy | |
2169 | mOffset[i] = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset)); | |
2170 | mOffset[i] = mOffset[i]*(-1); // adjust to direction of y-axes in online simulation | |
2171 | ||
2172 | if (mOffset[i] < 0) { | |
2173 | numAlu.AssignFormatted(-mOffset[i]); | |
2174 | numAlu.SetSign(-1); | |
2175 | } | |
2176 | else { | |
2177 | numAlu.AssignFormatted(mOffset[i]); | |
2178 | numAlu.SetSign(1); | |
2179 | } | |
2180 | ||
2181 | offsetAlu = numAlu; | |
2182 | mOffset[i] = offsetAlu.GetSignedValue(); | |
2183 | ||
2184 | ||
2185 | // slope: | |
2186 | ||
2187 | correctAlu.AssignDouble(mCorrectSlope); | |
2188 | mCorrectSlope = correctAlu.GetValueWhole(); | |
2189 | ||
2190 | mSlope[i] = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope); | |
2191 | ||
2192 | if (mSlope[i] < 0) { | |
2193 | numAlu.AssignFormatted(-mSlope[i]); | |
2194 | numAlu.SetSign(-1); | |
2195 | } | |
2196 | else { | |
2197 | numAlu.AssignFormatted(mSlope[i]); | |
2198 | numAlu.SetSign(1); | |
2199 | } | |
2200 | ||
2201 | slopeAlu = numAlu; // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding) | |
2202 | mSlope[i] = slopeAlu.GetSignedValue(); | |
2203 | ||
2204 | // local (LTU) limits for the deflection | |
2205 | // ATan returns angles in radian | |
2206 | alphaMax[i] = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m | |
2207 | slopeMin[i] = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope; | |
2208 | slopeMax[i] = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope; | |
2209 | ||
2210 | if (slopeMin[i] < 0) { | |
2211 | slopeAlu.AssignDouble(-slopeMin[i]); | |
2212 | slopeAlu.SetSign(-1); | |
2213 | } | |
2214 | else { | |
2215 | slopeAlu.AssignDouble(slopeMin[i]); | |
2216 | slopeAlu.SetSign(1); | |
2217 | } | |
2218 | mslopeMin[i] = slopeAlu.GetSignedValue(); // the borders should lie inside the range of mSlope -> usage of slopeAlu again | |
2219 | ||
2220 | if (slopeMax[i] < 0) { | |
2221 | slopeAlu.AssignDouble(-slopeMax[i]); | |
2222 | slopeAlu.SetSign(-1); | |
2223 | } | |
2224 | else { | |
2225 | slopeAlu.AssignDouble(slopeMax[i]); | |
2226 | slopeAlu.SetSign(1); | |
2227 | } | |
2228 | mslopeMax[i] = slopeAlu.GetSignedValue(); | |
2229 | } | |
2230 | ||
2231 | // suppress submission of tracks with low stiffness | |
2232 | // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber) | |
2233 | ||
2234 | // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible) | |
2235 | // up to now they are sorted according to maximum hit sum | |
2236 | // is the sorting really done in the TRAP-chip? | |
2237 | ||
2238 | Int_t order[4] = {-1,-1,-1,-1}; | |
2239 | Int_t wordnr = 0; // number of tracklet-words | |
2240 | ||
2241 | for(Int_t j = 0; j < fMaxTracklets; j++) { | |
2242 | //if( mADC[j] == -1) continue; | |
2243 | if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut | |
2244 | wordnr++; | |
2245 | if( wordnr-1 == 0) { | |
2246 | order[0] = j; | |
2247 | continue; | |
2248 | } | |
2249 | // wordnr-1>0, wordnr-1<4 | |
2250 | order[wordnr-1] = j; | |
2251 | for( Int_t k = 0; k < wordnr-1; k++) { | |
2252 | if( mOffset[j] < mOffset[order[k]] ) { | |
2253 | for( Int_t l = wordnr-1; l > k; l-- ) { | |
2254 | order[l] = order[l-1]; | |
2255 | } | |
2256 | order[k] = j; | |
2257 | break; | |
2258 | } | |
2259 | ||
2260 | } | |
2261 | } | |
2262 | ||
2263 | // fill the bit-words in ascending order and without gaps | |
2264 | UInt_t bitWord[4] = {0,0,0,0}; // attention: unsigned int to have real 32 bits (no 2-complement) | |
2265 | for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words | |
2266 | //Bool_t rem1 = kTRUE; | |
2267 | ||
2268 | Int_t i = order[j]; | |
2269 | bitWord[j] = 0; // invalid bit-word (bit-word is 2-complement and therefore without sign) | |
2270 | //if( mADC[i] == -1) continue; | |
2271 | ////if( (mADC[i] == -1) || (mSlope[i] < mslopeMin[i]) || (mSlope[i] > mslopeMax[i])) continue; //don't transmit bit word | |
2272 | bitWord[j] = 1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored | |
2273 | //printf("\n"); | |
2274 | ||
2275 | /* | |
2276 | // pad position | |
2277 | if(mOffset[i] < 0) { | |
2278 | rem1 = kFALSE; // don't remove the first 1 | |
2279 | //printf("1"); | |
2280 | for(Int_t iBit = 1; iBit < 13; iBit++) { | |
2281 | bitWord[j] = bitWord[j]<<1; | |
2282 | bitWord[j] |= (1-((-mOffset[i])>>(12-iBit))&1); | |
2283 | //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1)); | |
2284 | } | |
2285 | } | |
2286 | else { | |
2287 | bitWord[j] |= 0; | |
2288 | //printf("0"); | |
2289 | for(Int_t iBit = 1; iBit < 13; iBit++) { | |
2290 | bitWord[j] = bitWord[j]<<1; | |
2291 | bitWord[j] |= (mOffset[i]>>(12-iBit))&1; | |
2292 | //printf("%d",(mOffset[i]>>(12-iBit))&1); | |
2293 | } | |
2294 | } | |
2295 | ||
2296 | // deflection length | |
2297 | bitWord[j] = bitWord[j]<<1; | |
2298 | if(mSlope[i] < 0) { | |
2299 | bitWord[j] |= 1; | |
2300 | //printf("1"); | |
2301 | for(Int_t iBit = 1; iBit < 7; iBit++) { | |
2302 | bitWord[j] = bitWord[j]<<1; | |
2303 | bitWord[j] |= (1-((-mSlope[i])>>(6-iBit))&1); | |
2304 | //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1)); | |
2305 | } | |
2306 | } | |
2307 | else { | |
2308 | bitWord[j] |= 0; | |
2309 | //printf("0"); | |
2310 | for(Int_t iBit = 1; iBit < 7; iBit++) { | |
2311 | bitWord[j] = bitWord[j]<<1; | |
2312 | bitWord[j] |= (mSlope[i]>>(6-iBit))&1; | |
2313 | //printf("%d",(mSlope[i]>>(6-iBit))&1); | |
2314 | } | |
2315 | } | |
2316 | ||
2317 | // pad row | |
2318 | for(Int_t iBit = 0; iBit < 4; iBit++) { | |
2319 | bitWord[j] = bitWord[j]<<1; | |
2320 | bitWord[j] |= (fRow>>(3-iBit))&1; | |
2321 | //printf("%d", (fRow>>(3-iBit))&1); | |
2322 | } | |
2323 | ||
2324 | // electron probability (currently not implemented; the mean charge is just scaled) | |
2325 | for(Int_t iBit = 0; iBit < 8; iBit++) { | |
2326 | bitWord[j] = bitWord[j]<<1; | |
2327 | bitWord[j] |= (mMeanCharge[i]>>(7-iBit))&1; | |
2328 | //printf("0"); | |
2329 | } | |
2330 | ||
2331 | ||
2332 | if (rem1 == kTRUE) { | |
2333 | bitWord[j] = bitWord[j] - (1<<31); | |
2334 | }*/ | |
2335 | ||
2336 | ||
2337 | /*printf("mean charge: %d\n",mMeanCharge[i]); | |
2338 | printf("row: %d\n",fRow); | |
2339 | printf("slope: %d\n",mSlope[i]); | |
2340 | printf("pad position: %d\n",mOffset[i]); | |
2341 | printf("channel: %d\n",mADC[i]);*/ | |
2342 | ||
2343 | // electron probability (currently not implemented; the mean charge is just scaled) | |
2344 | for(Int_t iBit = 0; iBit < 8; iBit++) { | |
2345 | bitWord[j] = bitWord[j]<<1; | |
2346 | bitWord[j] |= (mMeanCharge[i]>>(7-iBit))&1; | |
2347 | //printf("0"); | |
2348 | } | |
2349 | ||
2350 | // pad row | |
2351 | for(Int_t iBit = 0; iBit < 4; iBit++) { | |
2352 | bitWord[j] = bitWord[j]<<1; | |
2353 | bitWord[j] |= (fRow>>(3-iBit))&1; | |
2354 | //printf("%d", (fRow>>(3-iBit))&1); | |
2355 | } | |
2356 | ||
2357 | // deflection length | |
2358 | bitWord[j] = bitWord[j]<<1; | |
2359 | if(mSlope[i] < 0) { | |
2360 | bitWord[j] |= 1; | |
2361 | //printf("1"); | |
2362 | for(Int_t iBit = 1; iBit < 7; iBit++) { | |
2363 | bitWord[j] = bitWord[j]<<1; | |
2364 | bitWord[j] |= (1-((-mSlope[i])>>(6-iBit))&1); | |
2365 | //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1)); | |
2366 | } | |
2367 | } | |
2368 | else { | |
2369 | bitWord[j] |= 0; | |
2370 | //printf("0"); | |
2371 | for(Int_t iBit = 1; iBit < 7; iBit++) { | |
2372 | bitWord[j] = bitWord[j]<<1; | |
2373 | bitWord[j] |= (mSlope[i]>>(6-iBit))&1; | |
2374 | //printf("%d",(mSlope[i]>>(6-iBit))&1); | |
2375 | } | |
2376 | } | |
2377 | ||
2378 | // pad position | |
2379 | bitWord[j] = bitWord[j]<<1; | |
2380 | if(mOffset[i] < 0) { | |
2381 | bitWord[j] |= 1; | |
2382 | //printf("1"); | |
2383 | for(Int_t iBit = 1; iBit < 13; iBit++) { | |
2384 | bitWord[j] = bitWord[j]<<1; | |
2385 | bitWord[j] |= (1-((-mOffset[i])>>(12-iBit))&1); | |
2386 | //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1)); | |
2387 | } | |
2388 | } | |
2389 | else { | |
2390 | bitWord[j] |= 0; | |
2391 | //printf("0"); | |
2392 | for(Int_t iBit = 1; iBit < 13; iBit++) { | |
2393 | bitWord[j] = bitWord[j]<<1; | |
2394 | bitWord[j] |= (mOffset[i]>>(12-iBit))&1; | |
2395 | //printf("%d",(mOffset[i]>>(12-iBit))&1); | |
2396 | } | |
2397 | } | |
2398 | ||
2399 | ||
2400 | ||
2401 | ||
2402 | //printf("bitWord: %u\n",bitWord[j]); | |
2403 | //printf("adc: %d\n",mADC[i]); | |
2404 | fMCMT[j] = bitWord[j]; | |
2405 | } | |
2406 | ||
2407 | //printf("\n"); | |
2408 | ||
2409 | ||
2410 | delete [] qsum; | |
2411 | delete [] ieffped; | |
2412 | ||
2413 | delete [] X; | |
2414 | delete [] XX; | |
2415 | delete [] Y; | |
2416 | delete [] YY; | |
2417 | delete [] XY; | |
2418 | delete [] N; | |
2419 | delete [] QT0; | |
2420 | delete [] QT1; | |
2421 | ||
2422 | delete [] hitSum; | |
2423 | delete [] selectPair; | |
2424 | ||
2425 | delete padPlane; | |
2426 | ||
2427 | ||
2428 | ||
2429 | /* | |
2430 | ||
2431 | ||
2432 | // output-part; creates some dump trees; output should not be organized inside the AliTRDmcmSim-class | |
2433 | ||
2434 | // structure: in system directory "./TRD_Tracklet" a root-file called "TRD_readout_tree.root" is stored with subdirectories SMxx/sx; | |
2435 | // in each of these subdirectories 6 trees according to layers are saved, called lx; | |
2436 | // whenever a mcm of that layer had a bit-word!=0, a branch containing an array with 4 (possibly some valued 0) elements is added; | |
2437 | // branch-name: mcmxxxwd; | |
2438 | // another branch contains the channel-number (mcmxxxch) | |
2439 | ||
2440 | ||
2441 | AliLog::SetClassDebugLevel("AliTRDmcmSim", 10); | |
2442 | AliLog::SetFileOutput("../log/tracklet.log"); | |
2443 | ||
2444 | ||
2445 | UInt_t* trackletWord; | |
2446 | Int_t* adcChannel; | |
2447 | ||
2448 | Int_t u = 0; | |
2449 | ||
2450 | // testing for wordnr in order to speed up the simulation | |
2451 | ||
2452 | if (wordnr == 0 && fNextEvent == 0) { | |
2453 | return; | |
2454 | } | |
2455 | ||
2456 | ||
2457 | Int_t mcmNr = fRobPos * (fGeo->MCMmax()) + fMcmPos; | |
2458 | ||
2459 | Char_t* SMName = new Char_t[4]; | |
2460 | Char_t* stackName = new Char_t[2]; | |
2461 | Char_t* layerName = new Char_t[2]; | |
2462 | ||
2463 | Char_t* treeName = new Char_t[2]; | |
2464 | Char_t* treeTitle = new Char_t[8]; | |
2465 | Char_t* branchNameWd = new Char_t[8]; // mcmxxxwd bit word | |
2466 | Char_t* branchNameCh = new Char_t[8]; // mcmxxxch channel | |
2467 | ||
2468 | Char_t* dirName = NULL; | |
2469 | Char_t* treeFile = NULL; | |
2470 | Char_t* evFile = NULL; | |
2471 | Char_t* curDir = new Char_t[255]; | |
2472 | ||
2473 | for (Int_t i = 0; i<255; i++) { | |
2474 | curDir[i]='n'; | |
2475 | } | |
2476 | sprintf(curDir,"%s",gSystem->BaseName(gSystem->WorkingDirectory())); | |
2477 | Int_t rawEvent = 0; | |
2478 | Int_t nrPos = 3; | |
2479 | ||
2480 | ||
2481 | while(curDir[nrPos]!='n'){ | |
2482 | ||
2483 | ||
2484 | switch(curDir[nrPos]) { | |
2485 | case '0': | |
2486 | rawEvent = rawEvent*10 + 0; | |
2487 | break; | |
2488 | case '1': | |
2489 | rawEvent = rawEvent*10 + 1; | |
2490 | break; | |
2491 | case '2': | |
2492 | rawEvent = rawEvent*10 + 2; | |
2493 | break; | |
2494 | case '3': | |
2495 | rawEvent = rawEvent*10 + 3; | |
2496 | break; | |
2497 | case '4': | |
2498 | rawEvent = rawEvent*10 + 4; | |
2499 | break; | |
2500 | case '5': | |
2501 | rawEvent = rawEvent*10 + 5; | |
2502 | break; | |
2503 | case '6': | |
2504 | rawEvent = rawEvent*10 + 6; | |
2505 | break; | |
2506 | case '7': | |
2507 | rawEvent = rawEvent*10 + 7; | |
2508 | break; | |
2509 | case '8': | |
2510 | rawEvent = rawEvent*10 + 8; | |
2511 | break; | |
2512 | case '9': | |
2513 | rawEvent = rawEvent*10 + 9; | |
2514 | break; | |
2515 | ||
2516 | } | |
2517 | nrPos = nrPos + 1; | |
2518 | } | |
2519 | delete [] curDir; | |
2520 | ||
2521 | if (!gSystem->ChangeDirectory("../TRD_Tracklet")) { | |
2522 | gSystem->MakeDirectory("../TRD_Tracklet"); | |
2523 | gSystem->ChangeDirectory("../TRD_Tracklet"); | |
2524 | } | |
2525 | ||
2526 | TFile *f = new TFile("TRD_readout_tree.root","update"); | |
2527 | TTree *tree = NULL; | |
2528 | TBranch *branch = NULL; | |
2529 | TBranch *branchCh = NULL; | |
2530 | ||
2531 | Int_t iEventNr = 0; | |
2532 | Int_t dignr = 10; // nr of digits of a integer | |
2533 | Int_t space = 1; // additional char-space | |
2534 | ||
2535 | ||
2536 | evFile = new Char_t[2+space]; | |
2537 | sprintf(evFile,"ev%d",iEventNr); | |
2538 | ||
2539 | ||
2540 | while(f->cd(evFile)){ | |
2541 | iEventNr = iEventNr + 1; | |
2542 | if (iEventNr/dignr > 0) { | |
2543 | dignr = dignr * 10; | |
2544 | space = space + 1; | |
2545 | } | |
2546 | delete [] evFile; | |
2547 | evFile = NULL; | |
2548 | evFile = new Char_t[2+space]; | |
2549 | sprintf(evFile,"ev%d",iEventNr); | |
2550 | } | |
2551 | ||
2552 | if(iEventNr == rawEvent) { fNextEvent = 1; } // new event | |
2553 | ||
2554 | if (fNextEvent == 1) { | |
2555 | fNextEvent = 0; | |
2556 | // turn to head directory | |
2557 | f->mkdir(evFile); | |
2558 | f->cd(evFile); | |
2559 | // create all subdirectories and trees in case of new event | |
2560 | ||
2561 | ||
2562 | for (Int_t iSector = 0; iSector < 18; iSector++) { | |
2563 | ||
2564 | if (iSector < 10) { | |
2565 | sprintf(SMName,"SM0%d",iSector); | |
2566 | } | |
2567 | else { | |
2568 | sprintf(SMName,"SM%d",iSector); | |
2569 | } | |
2570 | ||
2571 | for (Int_t iStack = 0; iStack < 5; iStack++) { | |
2572 | sprintf(stackName,"s%d",iStack); | |
2573 | ||
2574 | f->cd(evFile); | |
2575 | if (iStack == 0) { | |
2576 | gDirectory->mkdir(SMName); | |
2577 | } | |
2578 | gDirectory->cd(SMName); | |
2579 | gDirectory->mkdir(stackName); | |
2580 | gDirectory->cd(stackName); | |
2581 | ||
2582 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
2583 | sprintf(layerName,"l%d",iLayer); | |
2584 | sprintf(treeName,"%s",layerName); | |
2585 | sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName); | |
2586 | tree = new TTree(treeName,treeTitle); | |
2587 | tree->Write("",TObject::kOverwrite); | |
2588 | delete tree; | |
2589 | tree = NULL; | |
2590 | } | |
2591 | } | |
2592 | } | |
2593 | ||
2594 | ||
2595 | } | |
2596 | ||
2597 | ||
2598 | else { | |
2599 | iEventNr = iEventNr - 1; | |
2600 | dignr = dignr/10; | |
2601 | if (iEventNr/dignr == 0) space = space - 1; | |
2602 | delete [] evFile; | |
2603 | evFile = NULL; | |
2604 | evFile = new Char_t[2+space]; | |
2605 | sprintf(evFile,"ev%d",iEventNr); | |
2606 | } | |
2607 | ||
2608 | if (wordnr == 0) { | |
2609 | delete [] SMName; | |
2610 | delete [] stackName; | |
2611 | delete [] layerName; | |
2612 | delete [] treeName; | |
2613 | delete [] treeTitle; | |
2614 | delete [] branchNameWd; | |
2615 | delete [] branchNameCh; | |
2616 | delete [] evFile; | |
2617 | f->Close(); | |
2618 | dirName = new Char_t[6+space]; | |
2619 | sprintf(dirName,"../raw%d",iEventNr); | |
2620 | gSystem->ChangeDirectory(dirName); | |
2621 | delete [] dirName; | |
2622 | return; | |
2623 | } | |
2624 | ||
2625 | dirName = new Char_t[6+space]; | |
2626 | sprintf(dirName,"../raw%d",iEventNr); | |
2627 | ||
2628 | ||
2629 | f->cd(evFile); | |
2630 | ||
2631 | if (fSector < 10) { | |
2632 | sprintf(SMName,"SM0%d",fSector); | |
2633 | } | |
2634 | else { | |
2635 | sprintf(SMName,"SM%d",fSector); | |
2636 | } | |
2637 | sprintf(stackName,"s%d",fStack); | |
2638 | sprintf(layerName,"l%d",fLayer); | |
2639 | sprintf(treeName,"%s",layerName); | |
2640 | sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName); | |
2641 | ||
2642 | treeFile = new Char_t[13+space]; | |
2643 | sprintf(treeFile,"%s/%s/%s/%s",evFile,SMName,stackName,treeName); | |
2644 | delete [] evFile; | |
2645 | evFile = NULL; | |
2646 | gDirectory->cd(SMName); | |
2647 | gDirectory->cd(stackName); | |
2648 | tree = (TTree*)f->Get(treeFile); | |
2649 | delete [] treeFile; | |
2650 | treeFile = NULL; | |
2651 | ||
2652 | ||
2653 | //make branch with number of words and fill | |
2654 | ||
2655 | if (mcmNr < 10) { | |
2656 | sprintf(branchNameWd,"mcm00%dwd",mcmNr); | |
2657 | sprintf(branchNameCh,"mcm00%dch",mcmNr); | |
2658 | } | |
2659 | else { | |
2660 | if (mcmNr < 100) { | |
2661 | sprintf(branchNameWd,"mcm0%dwd",mcmNr); | |
2662 | sprintf(branchNameCh,"mcm0%dch",mcmNr); | |
2663 | } | |
2664 | else { | |
2665 | sprintf(branchNameWd,"mcm%dwd",mcmNr); | |
2666 | sprintf(branchNameCh,"mcm%dch",mcmNr); | |
2667 | } | |
2668 | } | |
2669 | ||
2670 | ||
2671 | ||
2672 | // fill the tracklet word; here wordnr > 0 | |
2673 | ||
2674 | trackletWord = new UInt_t[fMaxTracklets]; | |
2675 | adcChannel = new Int_t[fMaxTracklets]; | |
2676 | ||
2677 | for (Int_t j = 0; j < fMaxTracklets; j++) { | |
2678 | Int_t i = order[j]; | |
2679 | trackletWord[j] = 0; | |
2680 | adcChannel[j] = -1; | |
2681 | if (bitWord[j]!=0) { | |
2682 | trackletWord[u] = bitWord[j]; | |
2683 | adcChannel[u] = mADC[i]; // mapping onto the original adc-array to be in line with the digits-adc-ordering (21 channels in total on 1 mcm, 18 belonging to pads); mADC[i] should be >-1 in case bitWord[i]>0 | |
2684 | ||
2685 | //fMCMT[u] = bitWord[j]; | |
2686 | u = u + 1; | |
2687 | } | |
2688 | } | |
2689 | ||
2690 | branch = tree->GetBranch(branchNameWd); | |
2691 | if(!branch) { | |
2692 | //make branch and fill | |
2693 | branch = tree->Branch(branchNameWd,trackletWord,"trackletWord[4]/i"); // 32 bit unsigned integer | |
2694 | branch->Fill(); | |
2695 | } | |
2696 | ||
2697 | branchCh = tree->GetBranch(branchNameCh); | |
2698 | if(!branchCh) { | |
2699 | //make branch and fill | |
2700 | branchCh = tree->Branch(branchNameCh,adcChannel,"adcChannel[4]/i"); // 32 bit unsigned integer | |
2701 | branchCh->Fill(); | |
2702 | } | |
2703 | ||
2704 | ||
2705 | tree->Write("",TObject::kOverwrite); | |
2706 | ||
2707 | ||
2708 | delete [] SMName; | |
2709 | delete [] stackName; | |
2710 | delete [] layerName; | |
2711 | delete [] treeName; | |
2712 | delete [] treeTitle; | |
2713 | delete [] branchNameWd; | |
2714 | delete [] branchNameCh; | |
2715 | delete [] trackletWord; | |
2716 | delete [] adcChannel; | |
2717 | ||
2718 | f->Close(); | |
2719 | gSystem->ChangeDirectory(dirName); | |
2720 | delete [] dirName; | |
2721 | ||
2722 | */ | |
2723 | ||
2724 | // to be done: | |
2725 | // error measure for quality of fit (not necessarily needed for the trigger) | |
2726 | // cluster quality threshold (not yet set) | |
2727 | // electron probability | |
2728 | ||
2729 | ||
2730 | ||
2731 | } | |
2732 | ||
2733 | ||
2734 | ||
2735 | ||
2736 | ||
2737 |