Corrected destructors (J. Thaeder)
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmSim.cxx
CommitLineData
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 //
b0a41e80 19// which simulated the TRAP processing after the AD-conversion //
20// The relevant parameters (i.e. configuration registers of the TRAP //
21// configuration are taken from AliTRDtrapConfig. //
dfd03fc3 22// //
23///////////////////////////////////////////////////////////////////////////////
24
b0a41e80 25#include <fstream> // needed for raw data dump
ecf39416 26
1d93b218 27#include <TCanvas.h>
b0a41e80 28#include <TH1F.h>
29#include <TH2F.h>
30#include <TGraph.h>
31#include <TLine.h>
ecf39416 32#include <TMath.h>
b0a41e80 33#include <TRandom.h>
34#include <TClonesArray.h>
0c349049 35
dfd03fc3 36#include "AliLog.h"
b0a41e80 37#include "AliRun.h"
38#include "AliRunLoader.h"
39#include "AliLoader.h"
40#include "AliTRDdigit.h"
0c349049 41
dfd03fc3 42#include "AliTRDfeeParam.h"
b0a41e80 43#include "AliTRDtrapConfig.h"
ecf39416 44#include "AliTRDSimParam.h"
dfd03fc3 45#include "AliTRDgeometry.h"
ecf39416 46#include "AliTRDcalibDB.h"
4cc89512 47#include "AliTRDdigitsManager.h"
b65e5048 48#include "AliTRDarrayADC.h"
1d93b218 49#include "AliTRDpadPlane.h"
52c19022 50#include "AliTRDtrackletMCM.h"
b0a41e80 51#include "AliTRDmcmSim.h"
1d93b218 52
dfd03fc3 53ClassImp(AliTRDmcmSim)
54
55//_____________________________________________________________________________
b0a41e80 56AliTRDmcmSim::AliTRDmcmSim() : TObject()
dfd03fc3 57 ,fInitialized(kFALSE)
23200400 58 ,fMaxTracklets(-1)
b0a41e80 59 ,fDetector(-1)
6e5d4cb2 60 ,fRobPos(-1)
61 ,fMcmPos(-1)
23200400 62 ,fRow (-1)
dfd03fc3 63 ,fNADC(-1)
64 ,fNTimeBin(-1)
dfd03fc3 65 ,fADCR(NULL)
66 ,fADCF(NULL)
b0a41e80 67 ,fMCMT(NULL)
68 ,fTrackletArray(NULL)
dfd03fc3 69 ,fZSM(NULL)
70 ,fZSM1Dim(NULL)
6e5d4cb2 71 ,fFeeParam(NULL)
b0a41e80 72 ,fTrapConfig(NULL)
6e5d4cb2 73 ,fSimParam(NULL)
74 ,fCal(NULL)
75 ,fGeo(NULL)
b0a41e80 76 ,fPedAcc(NULL)
77 ,fGainCounterA(NULL)
78 ,fGainCounterB(NULL)
79 ,fTailAmplLong(NULL)
80 ,fTailAmplShort(NULL)
81 ,fNHits(0)
82 ,fFitReg(NULL)
dfd03fc3 83{
84 //
b0a41e80 85 // AliTRDmcmSim default constructor
dfd03fc3 86 // By default, nothing is initialized.
87 // It is necessary to issue Init before use.
88}
89
dfd03fc3 90AliTRDmcmSim::~AliTRDmcmSim()
91{
92 //
93 // AliTRDmcmSim destructor
94 //
0c349049 95
b0a41e80 96 if(fInitialized) {
dfd03fc3 97 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
16e077d0 98 delete [] fADCR[iadc];
99 delete [] fADCF[iadc];
100 delete [] fZSM [iadc];
dfd03fc3 101 }
16e077d0 102 delete [] fADCR;
103 delete [] fADCF;
104 delete [] fZSM;
105 delete [] fZSM1Dim;
1d93b218 106 delete [] fMCMT;
b0a41e80 107
108 delete [] fPedAcc;
109 delete [] fGainCounterA;
110 delete [] fGainCounterB;
111 delete [] fTailAmplLong;
112 delete [] fTailAmplShort;
113 delete [] fFitReg;
114
115 fTrackletArray->Delete();
116 delete fTrackletArray;
117 delete fGeo;
1d93b218 118 }
dfd03fc3 119}
120
b0a41e80 121void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ )
dfd03fc3 122{
0c349049 123 //
dfd03fc3 124 // Initialize the class with new geometry information
125 // fADC array will be reused with filled by zero
0c349049 126 //
96e6312d 127
b0a41e80 128 if (!fInitialized) {
129 fFeeParam = AliTRDfeeParam::Instance();
130 fTrapConfig = AliTRDtrapConfig::Instance();
131 fSimParam = AliTRDSimParam::Instance();
132 fCal = AliTRDcalibDB::Instance();
133 fGeo = new AliTRDgeometry();
134 }
135
136 fDetector = det;
0c349049 137 fRobPos = robPos;
138 fMcmPos = mcmPos;
dfd03fc3 139 fNADC = fFeeParam->GetNadcMcm();
ecf39416 140 fNTimeBin = fCal->GetNumberOfTimeBins();
dfd03fc3 141 fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
b0a41e80 142 fMaxTracklets = fFeeParam->GetMaxNrOfTracklets();
23200400 143
b0a41e80 144 if (!fInitialized) {
dfd03fc3 145 fADCR = new Int_t *[fNADC];
146 fADCF = new Int_t *[fNADC];
147 fZSM = new Int_t *[fNADC];
148 fZSM1Dim = new Int_t [fNADC];
b0a41e80 149 fGainCounterA = new UInt_t[fNADC];
150 fGainCounterB = new UInt_t[fNADC];
dfd03fc3 151 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
152 fADCR[iadc] = new Int_t[fNTimeBin];
153 fADCF[iadc] = new Int_t[fNTimeBin];
154 fZSM [iadc] = new Int_t[fNTimeBin];
155 }
b0a41e80 156
157 // filter registers
158 fPedAcc = new UInt_t[fNADC]; // accumulator for pedestal filter
159 fTailAmplLong = new UShort_t[fNADC];
160 fTailAmplShort = new UShort_t[fNADC];
161
162 // tracklet calculation
163 fFitReg = new FitReg_t[fNADC];
164 fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fMaxTracklets);
165
166 fMCMT = new UInt_t[fMaxTracklets];
dfd03fc3 167 }
168
b0a41e80 169 fInitialized = kTRUE;
170
171 Reset();
172}
173
174void AliTRDmcmSim::Reset()
175{
176 // Resets the data values and internal filter registers
177 // by re-initialising them
178
dfd03fc3 179 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
180 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
181 fADCR[iadc][it] = 0;
182 fADCF[iadc][it] = 0;
183 fZSM [iadc][it] = 1; // Default unread = 1
184 }
185 fZSM1Dim[iadc] = 1; // Default unread = 1
b0a41e80 186 fGainCounterA[iadc] = 0;
187 fGainCounterB[iadc] = 0;
dfd03fc3 188 }
ecf39416 189
1d93b218 190 for(Int_t i = 0; i < fMaxTracklets; i++) {
191 fMCMT[i] = 0;
192 }
b0a41e80 193
194 FilterPedestalInit();
195 FilterGainInit();
196 FilterTailInit(fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP)); //??? not really correct if gain filter is active
197}
1d93b218 198
ab9f7002 199Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm)
b0a41e80 200{
201 // loads the ADC data as obtained from the digitsManager for the specified MCM
202
203 if (!CheckInitialized())
204 Init(det, rob, mcm);
205
206 if (!runloader) {
207 AliError("No Runloader given");
208 return kFALSE;
209 }
210
211 AliLoader *trdLoader = runloader->GetLoader("TRDLoader");
212 if (!trdLoader) {
213 AliError("Could not get TRDLoader");
214 return kFALSE;
215 }
216
217 trdLoader->LoadDigits();
218 AliTRDdigitsManager *digMgr = new AliTRDdigitsManager();
219 digMgr->SetSDigits(0);
220 digMgr->CreateArrays();
221 digMgr->ReadDigits(trdLoader->TreeD());
222 AliTRDarrayADC *digits = (AliTRDarrayADC*) digMgr->GetDigits(det);
223 if (!digits->HasData())
224 return kFALSE;
225 digits->Expand();
226
227 Int_t padrow = fFeeParam->GetPadRowFromMCM(rob, mcm);
228 Int_t padcol = 0;
229 for (Int_t ch = 0; ch < fNADC; ch++) {
230 for (Int_t tb = 0; tb < fNTimeBin; tb++) {
231 padcol = fFeeParam->GetPadColFromADC(rob, mcm, ch);
232 if (padcol < 0) {
233 fADCR[ch][tb] = 0;
234 fADCF[ch][tb] = 0;
235 }
236 else {
237 if (digits->GetData(padrow,padcol, tb) < 0) {
238 fADCR[ch][tb] = 0;
239 fADCF[ch][tb] = 0;
240 }
241 else {
242 fADCR[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits;
243 fADCF[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits;
244 }
245 }
246 }
247 }
248 digMgr->RemoveDigits(det);
249 delete digMgr;
1d93b218 250
b0a41e80 251 return kTRUE;
252}
253
254void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inputGain, Int_t inputTail)
255{
256 // This function can be used to test the filters.
257 // It feeds nsamples of ADC values with a gaussian distribution specified by mean and sigma.
258 // The filter chain implemented here consists of:
259 // Pedestal -> Gain -> Tail
260 // With inputGain and inputTail the input to the gain and tail filter, respectively,
261 // can be chosen where
262 // 0: noise input
263 // 1: pedestal output
264 // 2: gain output
265 // The input has to be chosen from a stage before.
266 // The filter behaviour is controlled by the TRAP parameters from AliTRDtrapConfig in the
267 // same way as in normal simulation.
268 // The functions produces four histograms with the values at the different stages.
269
270 TH1F *h = new TH1F("noise", "Gaussian Noise;sample;ADC count",
271 nsamples, 0, nsamples);
272 TH1F *hfp = new TH1F("pedf", "Noise #rightarrow Pedestal filter;sample;ADC count", nsamples, 0, nsamples);
273 TH1F *hfg = new TH1F("pedg", "Pedestal #rightarrow Gain;sample;ADC count", nsamples, 0, nsamples);
274 TH1F *hft = new TH1F("pedt", "Gain #rightarrow Tail;sample;ADC count", nsamples, 0, nsamples);
275 h->SetStats(kFALSE);
276 hfp->SetStats(kFALSE);
277 hfg->SetStats(kFALSE);
278 hft->SetStats(kFALSE);
279
280 Int_t value; // ADC count with noise (10 bit)
281 Int_t valuep; // pedestal filter output (12 bit)
282 Int_t valueg; // gain filter output (12 bit)
283 Int_t valuet; // tail filter value (12 bit)
284
285 for (Int_t i = 0; i < nsamples; i++) {
286 value = (Int_t) gRandom->Gaus(mean, sigma); // generate noise with gaussian distribution
287 h->SetBinContent(i, value);
288
289 valuep = FilterPedestalNextSample(1, 0, ((Int_t) value) << 2);
290
291 if (inputGain == 0)
292 valueg = FilterGainNextSample(1, ((Int_t) value) << 2);
293 else
294 valueg = FilterGainNextSample(1, valuep);
295
296 if (inputTail == 0)
297 valuet = FilterTailNextSample(1, ((Int_t) value) << 2);
298 else if (inputTail == 1)
299 valuet = FilterTailNextSample(1, valuep);
300 else
301 valuet = FilterTailNextSample(1, valueg);
302
303 hfp->SetBinContent(i, valuep >> 2);
304 hfg->SetBinContent(i, valueg >> 2);
305 hft->SetBinContent(i, valuet >> 2);
306 }
307
308 TCanvas *c = new TCanvas;
309 c->Divide(2,2);
310 c->cd(1);
311 h->Draw();
312 c->cd(2);
313 hfp->Draw();
314 c->cd(3);
315 hfg->Draw();
316 c->cd(4);
317 hft->Draw();
dfd03fc3 318}
319
ecf39416 320Bool_t AliTRDmcmSim::CheckInitialized()
321{
0c349049 322 //
323 // Check whether object is initialized
324 //
325
ecf39416 326 if( ! fInitialized ) {
327 AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
328 }
329 return fInitialized;
330}
331
ab9f7002 332void AliTRDmcmSim::Print(Option_t* const option) const
b0a41e80 333{
334 // Prints the data stored and/or calculated for this MCM.
335 // The output is controlled by option which can be a sequence of any of
336 // the following characters:
337 // R - prints raw ADC data
338 // F - prints filtered data
339 // H - prints detected hits
340 // T - prints found tracklets
341 // The later stages are only useful when the corresponding calculations
342 // have been performed.
343
344 printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector);
345
346 TString opt = option;
347 if (opt.Contains("U")) {
348 printf("Raw ADC data (10 bit):\n");
349 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
350 for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) {
351 printf("%5i", fADCR[iChannel][iTimeBin] >> fgkAddDigits);
352 }
353 printf("\n");
354 }
1d93b218 355 }
356
b0a41e80 357 if (opt.Contains("F")) {
358 printf("Filtered data (12 bit):\n");
359 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
360 for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) {
361 printf("%5i", fADCF[iChannel][iTimeBin]);
362 }
363 printf("\n");
364 }
1d93b218 365 }
366
b0a41e80 367 if (opt.Contains("H")) {
368 printf("Found %i hits:\n", fNHits);
369 for (Int_t iHit = 0; iHit < fNHits; iHit++) {
370 printf("Hit %3i in timebin %2i, ADC %2i has charge %3i and position %3i\n",
ab9f7002 371 iHit, fHits[iHit].fTimebin, fHits[iHit].fChannel, fHits[iHit].fQtot, fHits[iHit].fYpos);
b0a41e80 372 }
1d93b218 373 }
1d93b218 374
b0a41e80 375 if (opt.Contains("T")) {
376 printf("Tracklets:\n");
377 for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntriesFast(); iTrkl++) {
378 printf("tracklet %i: 0x%08x\n", iTrkl, ((AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl])->GetTrackletWord());
379 }
1d93b218 380 }
b0a41e80 381}
1d93b218 382
ab9f7002 383void AliTRDmcmSim::Draw(Option_t* const option)
b0a41e80 384{
385 // Plots the data stored in a 2-dim. timebin vs. ADC channel plot.
386 // The option selects what data is plotted and can be a sequence of
387 // the following characters:
388 // R - plot raw data (default)
389 // F - plot filtered data (meaningless if R is specified)
390 // In addition to the ADC values:
391 // H - plot hits
392 // T - plot tracklets
393
394 TString opt = option;
395
396 TH2F *hist = new TH2F("mcmdata", Form("Data of MCM %i on ROB %i in detector %i", \
397 fMcmPos, fRobPos, fDetector), \
398 fNADC, -0.5, fNADC-.5, fNTimeBin, -.5, fNTimeBin-.5);
399 hist->GetXaxis()->SetTitle("ADC Channel");
400 hist->GetYaxis()->SetTitle("Timebin");
401 hist->SetStats(kFALSE);
402
403 if (opt.Contains("R")) {
404 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
405 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
406 hist->SetBinContent(iAdc+1, iTimeBin+1, fADCR[iAdc][iTimeBin] >> fgkAddDigits);
407 }
1d93b218 408 }
b0a41e80 409 }
410 else {
411 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
412 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
413 hist->SetBinContent(iAdc+1, iTimeBin+1, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
414 }
1d93b218 415 }
1d93b218 416 }
b0a41e80 417 hist->Draw("colz");
1d93b218 418
b0a41e80 419 if (opt.Contains("H")) {
420 TGraph *grHits = new TGraph();
421 for (Int_t iHit = 0; iHit < fNHits; iHit++) {
422 grHits->SetPoint(iHit,
ab9f7002 423 fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256.,
424 fHits[iHit].fTimebin);
b0a41e80 425 }
426 grHits->Draw("*");
427 }
1d93b218 428
b0a41e80 429 if (opt.Contains("T")) {
430 TLine *trklLines = new TLine[4];
431 for (Int_t iTrkl = 0; iTrkl < 4; iTrkl++) {
432 if (fMCMT[iTrkl] == 0x10001000)
433 break;
434 AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
435 AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl];
436 Float_t offset = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19)) + 19 * pp->GetWidthIPad();
437 trklLines[iTrkl].SetX1((offset - trkl->GetY()) / pp->GetWidthIPad());
438 trklLines[iTrkl].SetY1(0);
439 trklLines[iTrkl].SetX2((offset - (trkl->GetY() + ((Float_t) trkl->GetdY())*140e-4)) / pp->GetWidthIPad());
440 trklLines[iTrkl].SetY2(fNTimeBin - 1);
441 trklLines[iTrkl].SetLineColor(2);
442 trklLines[iTrkl].SetLineWidth(2);
443 printf("Tracklet %i: y = %f, dy = %f, offset = %f\n", iTrkl, trkl->GetY(), (trkl->GetdY() * 140e-4), offset);
444 trklLines[iTrkl].Draw();
445 }
446 }
1d93b218 447}
448
ab9f7002 449void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc )
dfd03fc3 450{
0c349049 451 //
dfd03fc3 452 // Store ADC data into array of raw data
0c349049 453 //
dfd03fc3 454
ecf39416 455 if( !CheckInitialized() ) return;
dfd03fc3 456
457 if( iadc < 0 || iadc >= fNADC ) {
b0a41e80 458 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
dfd03fc3 459 return;
460 }
461
462 for( int it = 0 ; it < fNTimeBin ; it++ ) {
b0a41e80 463 fADCR[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
464 fADCF[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
dfd03fc3 465 }
466}
467
dfd03fc3 468void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
469{
0c349049 470 //
dfd03fc3 471 // Store ADC data into array of raw data
0c349049 472 //
dfd03fc3 473
ecf39416 474 if( !CheckInitialized() ) return;
dfd03fc3 475
476 if( iadc < 0 || iadc >= fNADC ) {
477 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
478 return;
479 }
480
b0a41e80 481 fADCR[iadc][it] = adc << fgkAddDigits;
482 fADCF[iadc][it] = adc << fgkAddDigits;
483}
484
ab9f7002 485void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray)
b0a41e80 486{
ab9f7002 487 // Set the ADC data from an AliTRDarrayADC
488
b0a41e80 489 if (!fInitialized) {
ab9f7002 490 AliError("Called uninitialized! Nothing done!");
b0a41e80 491 return;
492 }
493
494 Int_t firstAdc = 0;
495 Int_t lastAdc = fNADC-1;
496
497 if (GetCol(firstAdc) > 143) {
498 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
499 fADCR[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
500 fADCF[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
501 }
502 firstAdc = 1;
503 }
504
505 if (GetCol(lastAdc) < 0) {
506 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
507 fADCR[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
508 fADCF[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
509 }
510 lastAdc = fNADC - 2;
511 }
512
513 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
514 for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
515 Int_t value = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin);
516 if (value < 0) {
517 fADCR[iAdc][iTimeBin] = 0;
518 fADCF[iAdc][iTimeBin] = 0;
519 }
520 else {
521 fADCR[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits;
522 fADCF[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits;
523 }
524 }
525 }
dfd03fc3 526}
527
dfd03fc3 528void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
529{
0c349049 530 //
dfd03fc3 531 // Store ADC data into array of raw data
0c349049 532 //
dfd03fc3 533
ecf39416 534 if( !CheckInitialized() ) return;
dfd03fc3 535
536 if( iadc < 0 || iadc >= fNADC ) {
537 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
538 return;
539 }
540
541 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
b0a41e80 542 fADCR[iadc][it] = fSimParam->GetADCbaseline() << fgkAddDigits;
543 fADCF[iadc][it] = fSimParam->GetADCbaseline() << fgkAddDigits;
dfd03fc3 544 }
545}
546
dfd03fc3 547Int_t AliTRDmcmSim::GetCol( Int_t iadc )
548{
0c349049 549 //
dfd03fc3 550 // Return column id of the pad for the given ADC channel
0c349049 551 //
552
ecf39416 553 if( !CheckInitialized() ) return -1;
dfd03fc3 554
ecf39416 555 return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
dfd03fc3 556}
557
b0a41e80 558Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
987ba9a3 559{
560 //
561 // Produce raw data stream from this MCM and put in buf
562 // Returns number of words filled, or negative value
563 // with -1 * number of overflowed words
564 //
565
566 UInt_t x;
987ba9a3 567 Int_t nw = 0; // Number of written words
568 Int_t of = 0; // Number of overflowed words
569 Int_t rawVer = fFeeParam->GetRAWversion();
570 Int_t **adc;
571 Int_t nActiveADC = 0; // number of activated ADC bits in a word
572
573 if( !CheckInitialized() ) return 0;
574
575 if( fFeeParam->GetRAWstoreRaw() ) {
576 adc = fADCR;
577 } else {
578 adc = fADCF;
579 }
580
b0a41e80 581 // Produce MCM header
6a04e92b 582 x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
b0a41e80 583
987ba9a3 584 if (nw < maxSize) {
585 buf[nw++] = x;
586 //printf("\nMCM header: %X ",x);
587 }
588 else {
589 of++;
590 }
591
592 // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
593 // n : unused , c : ADC count, m : selected ADCs
594 if( rawVer >= 3 ) {
595 x = 0;
596 for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
597 if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed
598 x = x | (1 << (iAdc+4) ); // last 4 digit reserved for 1100=0xc
599 nActiveADC++; // number of 1 in mmm....m
600 }
601 }
602 x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC; // nn = 01, ccccc are inverted, 0xc=1100
603 //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
604
605 if (nw < maxSize) {
606 buf[nw++] = x;
607 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
608 }
609 else {
610 of++;
611 }
612 }
613
614 // Produce ADC data. 3 timebins are packed into one 32 bits word
615 // In this version, different ADC channel will NOT share the same word
616
617 UInt_t aa=0, a1=0, a2=0, a3=0;
618
619 for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
620 if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
b0a41e80 621 aa = !(iAdc & 1) + 2;
987ba9a3 622 for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
b0a41e80 623 a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] >> fgkAddDigits : 0;
624 a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] >> fgkAddDigits : 0;
625 a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] >> fgkAddDigits : 0;
987ba9a3 626 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
627 if (nw < maxSize) {
b0a41e80 628 buf[nw++] = x;
629 //printf("%08X ",x);
987ba9a3 630 }
631 else {
b0a41e80 632 of++;
987ba9a3 633 }
634 }
635 }
636
637 if( of != 0 ) return -of; else return nw;
638}
639
1d93b218 640Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
641{
642 //
643 // Produce tracklet data stream from this MCM and put in buf
644 // Returns number of words filled, or negative value
645 // with -1 * number of overflowed words
646 //
647
648 UInt_t x;
649 Int_t nw = 0; // Number of written words
650 Int_t of = 0; // Number of overflowed words
651
652 if( !CheckInitialized() ) return 0;
653
654 // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM
655 // fMCMT is filled continuously until no more tracklet words available
656
657 Int_t wd = 0;
658 while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){
659 x = fMCMT[wd];
660 if (nw < maxSize) {
661 buf[nw++] = x;
662 }
663 else {
664 of++;
665 }
666 wd++;
667 }
668
669 if( of != 0 ) return -of; else return nw;
670}
671
dfd03fc3 672void AliTRDmcmSim::Filter()
673{
0c349049 674 //
b0a41e80 675 // Filter the raw ADC values. The active filter stages and their
676 // parameters are taken from AliTRDtrapConfig.
677 // The raw data is stored separate from the filtered data. Thus,
678 // it is possible to run the filters on a set of raw values
679 // sequentially for parameter tuning.
0c349049 680 //
dfd03fc3 681
b0a41e80 682 if( !CheckInitialized() ) {
683 AliError("got called before initialization! Nothing done!");
684 return;
dfd03fc3 685 }
686
b0a41e80 687 // Apply filters sequentially. Bypass is handled by filters
688 // since counters and internal registers may be updated even
689 // if the filter is bypassed.
690 // The first filter takes the data from fADCR and
691 // outputs to fADCF.
692
693 // Non-linearity filter not implemented.
694 FilterPedestal();
695 FilterGain();
696 FilterTail();
697 // Crosstalk filter not implemented.
dfd03fc3 698}
699
b0a41e80 700void AliTRDmcmSim::FilterPedestalInit()
701{
702 // Initializes the pedestal filter assuming that the input has
703 // been constant for a long time (compared to the time constant).
704
705// UShort_t fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
706 UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
707 UShort_t shifts[4] = {11, 14, 17, 21}; //??? where to take shifts from?
708
709 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
710 fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
711}
712
713UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
dfd03fc3 714{
b0a41e80 715 // Returns the output of the pedestal filter given the input value.
716 // The output depends on the internal registers and, thus, the
717 // history of the filter.
718
719 UShort_t fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
720 UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
721 UShort_t fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY); // 0..1 the bypass, active low
722 UShort_t shifts[4] = {11, 14, 17, 21}; //??? where to come from
723
724 UShort_t accumulatorShifted;
725 Int_t correction;
726 UShort_t inpAdd;
727
728 inpAdd = value + fpnp;
729
730 if (fpby == 0) //??? before or after update of accumulator
731 return value;
732
733 accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF; // 10 bits
734 if (timebin == 0) // the accumulator is disabled in the drift time
735 {
736 correction = (value & 0x3FF) - accumulatorShifted;
737 fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF; // 31 bits
738 }
739
740 if (inpAdd <= accumulatorShifted)
741 return 0;
742 else
743 {
744 inpAdd = inpAdd - accumulatorShifted;
745 if (inpAdd > 0xFFF)
746 return 0xFFF;
747 else
748 return inpAdd;
749 }
750}
23200400 751
b0a41e80 752void AliTRDmcmSim::FilterPedestal()
753{
0c349049 754 //
dfd03fc3 755 // Apply pedestal filter
0c349049 756 //
b0a41e80 757 // As the first filter in the chain it reads data from fADCR
758 // and outputs to fADCF.
759 // It has only an effect if previous samples have been fed to
760 // find the pedestal. Currently, the simulation assumes that
761 // the input has been stable for a sufficiently long time.
dfd03fc3 762
b0a41e80 763 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
764 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
765 fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]);
dfd03fc3 766 }
767 }
768}
769
b0a41e80 770void AliTRDmcmSim::FilterGainInit()
dfd03fc3 771{
b0a41e80 772 // Initializes the gain filter. In this case, only threshold
773 // counters are reset.
0c349049 774
b0a41e80 775 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
776 // these are counters which in hardware continue
777 // until maximum or reset
778 fGainCounterA[iAdc] = 0;
779 fGainCounterB[iAdc] = 0;
780 }
dfd03fc3 781}
782
b0a41e80 783UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
dfd03fc3 784{
b0a41e80 785 // Apply the gain filter to the given value.
786 // BEGIN_LATEX O_{i}(t) = #gamma_{i} * I_{i}(t) + a_{i} END_LATEX
787 // The output depends on the internal registers and, thus, the
788 // history of the filter.
dfd03fc3 789
b0a41e80 790 UShort_t fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY); // bypass, active low
791 UShort_t fgf = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc)); // 0x700 + (0 & 0x1ff);
792 UShort_t fga = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc)); // 40;
793 UShort_t fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA); // 20;
794 UShort_t fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB); // 2060;
dfd03fc3 795
b0a41e80 796 UInt_t tmp;
dfd03fc3 797
b0a41e80 798 value &= 0xFFF;
799 tmp = (value * fgf) >> 11;
800 if (tmp > 0xFFF) tmp = 0xFFF;
801
802 if (fgby == 1)
803 value = AddUintClipping(tmp, fga, 12);
804
805 // Update threshold counters
806 // not really useful as they are cleared with every new event
807 if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))
808 {
809 if (value >= fgtb)
810 fGainCounterB[adc]++;
811 else if (value >= fgta)
812 fGainCounterA[adc]++;
813 }
814
815 return value;
816}
817
818void AliTRDmcmSim::FilterGain()
819{
820 // Read data from fADCF and apply gain filter.
821
822 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
823 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
824 fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
1d93b218 825 }
b0a41e80 826 }
827}
1d93b218 828
b0a41e80 829void AliTRDmcmSim::FilterTailInit(Int_t baseline)
830{
831 // Initializes the tail filter assuming that the input has
832 // been at the baseline value (configured by FTFP) for a
833 // sufficiently long time.
834
835 // exponents and weight calculated from configuration
836 UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
837 UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
838 UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
839
840 Float_t lambdaL = lambdaLong * 1.0 / (1 << 11);
841 Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
842 Float_t alphaL = alphaLong * 1.0 / (1 << 11);
843 Float_t qup, qdn;
844 qup = (1 - lambdaL) * (1 - lambdaS);
845 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
846 Float_t kdc = qup/qdn;
847
848 Float_t kt, ql, qs;
849 UShort_t aout;
850
851 kt = kdc * baseline;
852 aout = baseline - (UShort_t) kt;
853 ql = lambdaL * (1 - lambdaS) * alphaL;
854 qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
855
856 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
857 fTailAmplLong[iAdc] = (UShort_t) (aout * ql / (ql + qs));
858 fTailAmplShort[iAdc] = (UShort_t) (aout * qs / (ql + qs));
859 }
860}
861
862UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
863{
864 // Returns the output of the tail filter for the given input value.
865 // The output depends on the internal registers and, thus, the
866 // history of the filter.
867
868 // exponents and weight calculated from configuration
869 UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
870 UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
871 UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
872
873 Float_t lambdaL = lambdaLong * 1.0 / (1 << 11);
874 Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
875 Float_t alphaL = alphaLong * 1.0 / (1 << 11);
876 Float_t qup, qdn;
877 qup = (1 - lambdaL) * (1 - lambdaS);
878 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
879// Float_t kdc = qup/qdn;
880
881 UInt_t aDiff;
882 UInt_t alInpv;
883 UShort_t aQ;
884 UInt_t tmp;
885
ab9f7002 886 UShort_t inpVolt = value & 0xFFF; // 12 bits
b0a41e80 887
888 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
889 return value;
890 else
891 {
892 // add the present generator outputs
893 aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
894
895 // calculate the difference between the input the generated signal
ab9f7002 896 if (inpVolt > aQ)
897 aDiff = inpVolt - aQ;
b0a41e80 898 else
899 aDiff = 0;
900
901 // the inputs to the two generators, weighted
902 alInpv = (aDiff * alphaLong) >> 11;
903
904 // the new values of the registers, used next time
905 // long component
906 tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
907 tmp = (tmp * lambdaLong) >> 11;
908 fTailAmplLong[adc] = tmp & 0xFFF;
909 // short component
910 tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
911 tmp = (tmp * lambdaShort) >> 11;
912 fTailAmplShort[adc] = tmp & 0xFFF;
913
914 // the output of the filter
915 return aDiff;
dfd03fc3 916 }
b0a41e80 917}
dfd03fc3 918
b0a41e80 919void AliTRDmcmSim::FilterTail()
920{
921 // Apply tail filter
fabf2e09 922
b0a41e80 923 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
924 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
925 fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]);
926 }
927 }
dfd03fc3 928}
929
dfd03fc3 930void AliTRDmcmSim::ZSMapping()
931{
0c349049 932 //
dfd03fc3 933 // Zero Suppression Mapping implemented in TRAP chip
934 //
935 // See detail TRAP manual "Data Indication" section:
936 // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
0c349049 937 //
dfd03fc3 938
b0a41e80 939 //??? values should come from TRAPconfig
940 Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS); // TRAP default = 0x4 (Tis=4)
941 Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT); // TRAP default = 0x28 (Tit=40)
942 Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL); // TRAP default = 0xf0
943 // (lookup table accept (I2,I1,I0)=(111)
944 // or (110) or (101) or (100))
945 Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); // TRAP default = 1 (no neighbor sensitivity)
ecf39416 946 Int_t ep = AliTRDfeeParam::GetPFeffectPedestal();
947
b0a41e80 948 Int_t **adc = fADCF;
dfd03fc3 949
b0a41e80 950 if( !CheckInitialized() ) {
ab9f7002 951 AliError("got called uninitialized! Nothing done!");
b0a41e80 952 return;
953 }
954
955 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
956 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
dfd03fc3 957
ecf39416 958 // Get ADC data currently in filter buffer
b0a41e80 959 Int_t ap = adc[iadc-1][it] - ep; // previous
960 Int_t ac = adc[iadc ][it] - ep; // current
961 Int_t an = adc[iadc+1][it] - ep; // next
ecf39416 962
dfd03fc3 963 // evaluate three conditions
0c349049 964 Int_t i0 = ( ac >= ap && ac >= an ) ? 0 : 1; // peak center detection
965 Int_t i1 = ( ap + ac + an > eBIT ) ? 0 : 1; // cluster
966 Int_t i2 = ( ac > eBIS ) ? 0 : 1; // absolute large peak
967
968 Int_t i = i2 * 4 + i1 * 2 + i0; // Bit position in lookup table
969 Int_t d = (eBIL >> i) & 1; // Looking up (here d=0 means true
970 // and d=1 means false according to TRAP manual)
971
972 fZSM[iadc][it] &= d;
973 if( eBIN == 0 ) { // turn on neighboring ADCs
974 fZSM[iadc-1][it] &= d;
975 fZSM[iadc+1][it] &= d;
dfd03fc3 976 }
977 }
978 }
979
980 // do 1 dim projection
981 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
982 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
ecf39416 983 fZSM1Dim[iadc] &= fZSM[iadc][it];
984 }
985 }
0c349049 986
ecf39416 987}
988
ecf39416 989void AliTRDmcmSim::DumpData( char *f, char *target )
990{
0c349049 991 //
ecf39416 992 // Dump data stored (for debugging).
993 // target should contain one or multiple of the following characters
994 // R for raw data
995 // F for filtered data
996 // Z for zero suppression map
997 // S Raw dat astream
998 // other characters are simply ignored
0c349049 999 //
1000
ecf39416 1001 UInt_t tempbuf[1024];
1002
1003 if( !CheckInitialized() ) return;
1004
1005 std::ofstream of( f, std::ios::out | std::ios::app );
1006 of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
b0a41e80 1007 fDetector, fGeo->GetSector(fDetector), fGeo->GetStack(fDetector),
1008 fGeo->GetSector(fDetector), fRobPos, fMcmPos );
ecf39416 1009
1010 for( int t=0 ; target[t] != 0 ; t++ ) {
1011 switch( target[t] ) {
1012 case 'R' :
1013 case 'r' :
1014 of << Form("fADCR (raw ADC data)\n");
1015 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1016 of << Form(" ADC %02d: ", iadc);
1017 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1018 of << Form("% 4d", fADCR[iadc][it]);
1019 }
1020 of << Form("\n");
1021 }
1022 break;
1023 case 'F' :
1024 case 'f' :
1025 of << Form("fADCF (filtered ADC data)\n");
1026 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1027 of << Form(" ADC %02d: ", iadc);
1028 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1029 of << Form("% 4d", fADCF[iadc][it]);
1030 }
1031 of << Form("\n");
1032 }
1033 break;
1034 case 'Z' :
1035 case 'z' :
1036 of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
1037 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1038 of << Form(" ADC %02d: ", iadc);
1039 if( fZSM1Dim[iadc] == 0 ) { of << " R " ; } else { of << " . "; } // R:read .:suppressed
1040 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1041 if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
1042 }
1043 of << Form("\n");
1044 }
1045 break;
1046 case 'S' :
1047 case 's' :
1048 Int_t s = ProduceRawStream( tempbuf, 1024 );
1049 of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
1050 of << Form(" address data\n");
1051 for( int i = 0 ; i < s ; i++ ) {
1052 of << Form(" %04x %08x\n", i, tempbuf[i]);
1053 }
dfd03fc3 1054 }
1055 }
1056}
1057
b0a41e80 1058void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label)
dfd03fc3 1059{
b0a41e80 1060 // Add the given hit to the fit register which is lateron used for
1061 // the tracklet calculation.
1062 // In addition to the fit sums in the fit register MC information
1063 // is stored.
1064
1065 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) &&
1066 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
ab9f7002 1067 fFitReg[adc].fQ0 += qtot;
b0a41e80 1068
1069 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) &&
1070 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
ab9f7002 1071 fFitReg[adc].fQ1 += qtot;
b0a41e80 1072
1073 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) &&
1074 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
1075 {
ab9f7002 1076 fFitReg[adc].fSumX += timebin;
1077 fFitReg[adc].fSumX2 += timebin*timebin;
1078 fFitReg[adc].fNhits++;
1079 fFitReg[adc].fSumY += ypos;
1080 fFitReg[adc].fSumY2 += ypos*ypos;
1081 fFitReg[adc].fSumXY += timebin*ypos;
b0a41e80 1082 }
1083
1084 // register hits (MC info)
ab9f7002 1085 fHits[fNHits].fChannel = adc;
1086 fHits[fNHits].fQtot = qtot;
1087 fHits[fNHits].fYpos = ypos;
1088 fHits[fNHits].fTimebin = timebin;
1089 fHits[fNHits].fLabel = label;
b0a41e80 1090 fNHits++;
dfd03fc3 1091}
1092
b0a41e80 1093void AliTRDmcmSim::CalcFitreg()
dfd03fc3 1094{
b0a41e80 1095 // Preprocessing.
1096 // Detect the hits and fill the fit registers.
1097 // Requires 12-bit data from fADCF which means Filter()
1098 // has to be called before even if all filters are bypassed.
1099
1100 //???
1101 // TRAP parameters:
ab9f7002 1102 const uint16_t lutPos[128] = { // move later to some other file
b0a41e80 1103 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15,
1104 16, 16, 16, 17, 17, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26,
1105 27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 27, 27, 27, 27, 26,
1106 26, 26, 26, 25, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 20, 20, 19, 18, 18, 17, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 7};
1107
1108 //??? to be clarified:
ab9f7002 1109 UInt_t adcMask = 0xfffff;
b0a41e80 1110
ab9f7002 1111 UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
b0a41e80 1112 Short_t ypos, fromLeft, fromRight, found;
ab9f7002 1113 UShort_t qTotal[19]; // the last is dummy
1114 UShort_t marked[6], qMarked[6], worse1, worse2;
b0a41e80 1115
1116 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS);
1117 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)
1118 < timebin1)
1119 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1120 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE);
1121 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)
1122 > timebin2)
1123 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1124
1125 // reset the fit registers
1126 fNHits = 0;
1127 for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels
1128 {
ab9f7002 1129 fFitReg[adcch].fNhits = 0;
1130 fFitReg[adcch].fQ0 = 0;
1131 fFitReg[adcch].fQ1 = 0;
1132 fFitReg[adcch].fSumX = 0;
1133 fFitReg[adcch].fSumY = 0;
1134 fFitReg[adcch].fSumX2 = 0;
1135 fFitReg[adcch].fSumY2 = 0;
1136 fFitReg[adcch].fSumXY = 0;
b0a41e80 1137 }
1138
1139 for (timebin = timebin1; timebin < timebin2; timebin++)
1140 {
ab9f7002 1141 // first find the hit candidates and store the total cluster charge in qTotal array
b0a41e80 1142 // in case of not hit store 0 there.
1143 for (adcch = 0; adcch < fNADC-2; adcch++) {
ab9f7002 1144 if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
b0a41e80 1145 {
ab9f7002 1146 adcLeft = fADCF[adcch ][timebin];
1147 adcCentral = fADCF[adcch+1][timebin];
1148 adcRight = fADCF[adcch+2][timebin];
b0a41e80 1149 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1)
ab9f7002 1150 hitQual = ( (adcLeft * adcRight) <
1151 (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) );
b0a41e80 1152 else
ab9f7002 1153 hitQual = 1;
b0a41e80 1154 // The accumulated charge is with the pedestal!!!
ab9f7002 1155 qtotTemp = adcLeft + adcCentral + adcRight;
1156 if ( (hitQual) &&
1157 (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1158 (adcLeft <= adcCentral) &&
1159 (adcCentral > adcRight) )
1160 qTotal[adcch] = qtotTemp;
b0a41e80 1161 else
ab9f7002 1162 qTotal[adcch] = 0;
1163 //printf("ch %2d qTotal %5d\n",adcch, qTotal[adcch]);
1d93b218 1164 }
b0a41e80 1165 else
ab9f7002 1166 qTotal[adcch] = 0; //jkl
1d93b218 1167 }
b0a41e80 1168
1169 fromLeft = -1;
1170 adcch = 0;
1171 found = 0;
1172 marked[4] = 19; // invalid channel
1173 marked[5] = 19; // invalid channel
ab9f7002 1174 qTotal[19] = 0;
b0a41e80 1175 while ((adcch < 16) && (found < 3))
1176 {
ab9f7002 1177 if (qTotal[adcch] > 0)
b0a41e80 1178 {
1179 fromLeft = adcch;
1180 marked[2*found+1]=adcch;
1181 found++;
1d93b218 1182 }
b0a41e80 1183 adcch++;
1d93b218 1184 }
1d93b218 1185
b0a41e80 1186 fromRight = -1;
1187 adcch = 18;
1188 found = 0;
1189 while ((adcch > 2) && (found < 3))
1190 {
ab9f7002 1191 if (qTotal[adcch] > 0)
b0a41e80 1192 {
1193 marked[2*found]=adcch;
1194 found++;
1195 fromRight = adcch;
1196 }
1197 adcch--;
1d93b218 1198 }
1d93b218 1199
b0a41e80 1200 //printf("Fromleft=%d, Fromright=%d\n",fromLeft, fromRight);
1201 // here mask the hit candidates in the middle, if any
1202 if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight))
1203 for (adcch = fromLeft+1; adcch < fromRight; adcch++)
ab9f7002 1204 qTotal[adcch] = 0;
1d93b218 1205
b0a41e80 1206 found = 0;
1207 for (adcch = 0; adcch < 19; adcch++)
ab9f7002 1208 if (qTotal[adcch] > 0) found++;
b0a41e80 1209 // NOT READY
1d93b218 1210
b0a41e80 1211 if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1212 {
1213 if (marked[4] == marked[5]) marked[5] = 19;
1214 for (found=0; found<6; found++)
1215 {
ab9f7002 1216 qMarked[found] = qTotal[marked[found]] >> 4;
1217 //printf("ch_%d qTotal %d qTotals %d |",marked[found],qTotal[marked[found]],qMarked[found]);
1d93b218 1218 }
b0a41e80 1219 //printf("\n");
1220
1221 Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
ab9f7002 1222 qMarked[0],
1223 qMarked[3],
1224 qMarked[4],
1225 qMarked[1],
1226 qMarked[2],
1227 qMarked[5],
b0a41e80 1228 &worse1, &worse2);
1229 // Now mask the two channels with the smallest charge
1230 if (worse1 < 19)
1231 {
ab9f7002 1232 qTotal[worse1] = 0;
b0a41e80 1233 //printf("Kill ch %d\n",worse1);
1d93b218 1234 }
b0a41e80 1235 if (worse2 < 19)
1236 {
ab9f7002 1237 qTotal[worse2] = 0;
b0a41e80 1238 //printf("Kill ch %d\n",worse2);
1d93b218 1239 }
1d93b218 1240 }
1d93b218 1241
b0a41e80 1242 for (adcch = 0; adcch < 19; adcch++) {
ab9f7002 1243 if (qTotal[adcch] > 0) // the channel is marked for processing
b0a41e80 1244 {
ab9f7002 1245 adcLeft = fADCF[adcch ][timebin];
1246 adcCentral = fADCF[adcch+1][timebin];
1247 adcRight = fADCF[adcch+2][timebin];
b0a41e80 1248 // hit detected, in TRAP we have 4 units and a hit-selection, here we proceed all channels!
1249 // subtract the pedestal TPFP, clipping instead of wrapping
1250
ab9f7002 1251 Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP);
1252// printf("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n",
1253// timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP,
b0a41e80 1254// fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT));
1255
ab9f7002 1256 if (adcLeft < regTPFP) adcLeft = 0; else adcLeft -= regTPFP;
1257 if (adcCentral < regTPFP) adcCentral = 0; else adcCentral -= regTPFP;
1258 if (adcRight < regTPFP) adcRight = 0; else adcRight -= regTPFP;
b0a41e80 1259 // Calculate the center of gravity
ab9f7002 1260 ypos = 128*(adcLeft - adcRight) / adcCentral;
b0a41e80 1261 if (ypos < 0) ypos = -ypos;
1262 // make the correction using the LUT
ab9f7002 1263 ypos = ypos + lutPos[ypos & 0x7F];
1264 if (adcLeft > adcRight) ypos = -ypos;
1265 AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, -1);
1d93b218 1266 }
1d93b218 1267 }
1d93b218 1268 }
b0a41e80 1269}
1d93b218 1270
b0a41e80 1271void AliTRDmcmSim::TrackletSelection()
1272{
1273 // Select up to 4 tracklet candidates from the fit registers
1274 // and assign them to the CPUs.
1275
ab9f7002 1276 UShort_t adcIdx, i, j, ntracks, tmp;
1277 UShort_t trackletCand[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates
b0a41e80 1278
1279 ntracks = 0;
ab9f7002 1280 for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs
1281 if ( (fFitReg[adcIdx].fNhits
b0a41e80 1282 >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) &&
ab9f7002 1283 (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits
b0a41e80 1284 >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT)))
1285 {
ab9f7002 1286 trackletCand[ntracks][0] = adcIdx;
1287 trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits;
1288 //printf("%d %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]);
b0a41e80 1289 ntracks++;
1290 };
1291
ab9f7002 1292 // for (i=0; i<ntracks;i++) printf("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]);
b0a41e80 1293
1294 if (ntracks > 4)
1295 {
1296 // primitive sorting according to the number of hits
1297 for (j = 0; j < (ntracks-1); j++)
1298 {
1299 for (i = j+1; i < ntracks; i++)
1300 {
ab9f7002 1301 if ( (trackletCand[j][1] < trackletCand[i][1]) ||
1302 ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
b0a41e80 1303 {
1304 // swap j & i
ab9f7002 1305 tmp = trackletCand[j][1];
1306 trackletCand[j][1] = trackletCand[i][1];
1307 trackletCand[i][1] = tmp;
1308 tmp = trackletCand[j][0];
1309 trackletCand[j][0] = trackletCand[i][0];
1310 trackletCand[i][0] = tmp;
b0a41e80 1311 }
1d93b218 1312 }
1d93b218 1313 }
b0a41e80 1314 ntracks = 4; // cut the rest, 4 is the max
1d93b218 1315 }
b0a41e80 1316 // else is not necessary to sort
1d93b218 1317
b0a41e80 1318 // now sort, so that the first tracklet going to CPU0 corresponds to the highest adc channel - as in the TRAP
1319 for (j = 0; j < (ntracks-1); j++)
1320 {
1321 for (i = j+1; i < ntracks; i++)
1322 {
ab9f7002 1323 if (trackletCand[j][0] < trackletCand[i][0])
b0a41e80 1324 {
1325 // swap j & i
ab9f7002 1326 tmp = trackletCand[j][1];
1327 trackletCand[j][1] = trackletCand[i][1];
1328 trackletCand[i][1] = tmp;
1329 tmp = trackletCand[j][0];
1330 trackletCand[j][0] = trackletCand[i][0];
1331 trackletCand[i][0] = tmp;
1d93b218 1332 }
1d93b218 1333 }
1334 }
b0a41e80 1335 for (i = 0; i < ntracks; i++) // CPUs with tracklets.
ab9f7002 1336 fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i]
b0a41e80 1337 for (i = ntracks; i < 4; i++) // CPUs without tracklets
1338 fFitPtr[i] = 31; // pointer to the left channel with tracklet for CPU[i] = 31 (invalid)
1339// printf("found %i tracklet candidates\n", ntracks);
1340}
1d93b218 1341
b0a41e80 1342void AliTRDmcmSim::FitTracklet()
1343{
1344 // Perform the actual tracklet fit based on the fit sums
1345 // which have been filled in the fit registers.
1346
1347 // parameters in fitred.asm (fit program)
1348 Int_t decPlaces = 5;
1349 Int_t rndAdd = 0;
1350 if (decPlaces > 1)
1351 rndAdd = (1 << (decPlaces-1)) + 1;
1352 else if (decPlaces == 1)
1353 rndAdd = 1;
1354
1355 // should come from trapConfig (DMEM)
1356 AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
1357 Long64_t shift = ((Long64_t) 1 << 32);
ab9f7002 1358 UInt_t scaleY = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4)));
1359 UInt_t scaleD = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4)));
b0a41e80 1360 int padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos);
1361 int yoffs = (fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19) - fFeeParam->GetNcol()/2) << (8 + decPlaces);
1362 int ndrift = 20; //??? value in simulation?
ab9f7002 1363 int deflCorr = 0; // -370;
b0a41e80 1364 int minslope = -10000; // no pt-cut so far
1365 int maxslope = 10000; // no pt-cut so far
1366
1367 // local variables for calculation
1368 Long64_t mult, temp, denom; //???
1369 UInt_t q0, q1, qTotal; // charges in the two windows and total charge
1370 UShort_t nHits; // number of hits
1371 Int_t slope, offset; // slope and offset of the tracklet
1372 Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers
1373 //int32_t SumY2; // not used in the current TRAP program
1374 FitReg_t *fit0, *fit1; // pointers to relevant fit registers
1375
1376// const uint32_t OneDivN[32] = { // 2**31/N : exactly like in the TRAP, the simple division here gives the same result!
1377// 0x00000000, 0x80000000, 0x40000000, 0x2AAAAAA0, 0x20000000, 0x19999990, 0x15555550, 0x12492490,
1378// 0x10000000, 0x0E38E380, 0x0CCCCCC0, 0x0BA2E8B0, 0x0AAAAAA0, 0x09D89D80, 0x09249240, 0x08888880,
1379// 0x08000000, 0x07878780, 0x071C71C0, 0x06BCA1A0, 0x06666660, 0x06186180, 0x05D17450, 0x0590B210,
1380// 0x05555550, 0x051EB850, 0x04EC4EC0, 0x04BDA120, 0x04924920, 0x0469EE50, 0x04444440, 0x04210840};
1381
1382 for (Int_t cpu = 0; cpu < 4; cpu++) {
1383 if (fFitPtr[cpu] == 31)
1384 {
1385 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1d93b218 1386 }
b0a41e80 1387 else
1388 {
1389 fit0 = &fFitReg[fFitPtr[cpu] ];
1390 fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1391
1392 mult = 1;
1393 mult = mult << (32 + decPlaces);
1394 mult = -mult;
1395
1396 // Merging
ab9f7002 1397 nHits = fit0->fNhits + fit1->fNhits; // number of hits
1398 sumX = fit0->fSumX + fit1->fSumX;
1399 sumX2 = fit0->fSumX2 + fit1->fSumX2;
b0a41e80 1400 denom = nHits*sumX2 - sumX*sumX;
1401
1402 mult = mult / denom; // exactly like in the TRAP program
ab9f7002 1403 q0 = fit0->fQ0 + fit1->fQ0;
1404 q1 = fit0->fQ1 + fit1->fQ1;
1405 sumY = fit0->fSumY + fit1->fSumY + 256*fit1->fNhits;
1406 sumXY = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
b0a41e80 1407
1408 slope = nHits*sumXY - sumX * sumY;
1409 offset = sumX2*sumY - sumX * sumXY;
1410 temp = mult * slope;
1411 slope = temp >> 32; // take the upper 32 bits
1412 temp = mult * offset;
1413 offset = temp >> 32; // take the upper 32 bits
1414
1415 offset = offset + yoffs + (18 << (8 + decPlaces));
ab9f7002 1416 slope = slope * ndrift + deflCorr;
b0a41e80 1417 offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
1d93b218 1418
b0a41e80 1419 if ((slope < minslope) || (slope > maxslope))
1420 {
1421 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1422 }
1423 else
1424 {
1425 temp = slope;
ab9f7002 1426 temp = temp * scaleD;
b0a41e80 1427 slope = (temp >> 32);
1428
1429 temp = offset;
ab9f7002 1430 temp = temp * scaleY;
b0a41e80 1431 offset = (temp >> 32);
1432
1433 // rounding, like in the TRAP
1434 slope = (slope + rndAdd) >> decPlaces;
1435 offset = (offset + rndAdd) >> decPlaces;
1436
1437 if (slope > 0x3f || slope < -0x3f)
1438 AliWarning("Overflow in slope");
1439 slope = slope & 0x7F; // 7 bit
1440
1441 if (offset > 0xfff || offset < 0xfff)
1442 AliWarning("Overflow in offset");
1443 offset = offset & 0x1FFF; // 13 bit
1444
1445 qTotal = (q1 / nHits) >> 1;
1446 if (qTotal > 0xff)
1447 AliWarning("Overflow in charge");
1448 qTotal = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1449
1450 // assemble and store the tracklet word
1451 fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1452 new ((*fTrackletArray)[cpu]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1453 }
1d93b218 1454 }
1d93b218 1455 }
b0a41e80 1456}
1d93b218 1457
b0a41e80 1458void AliTRDmcmSim::Tracklet()
1459{
ab9f7002 1460 // Run the tracklet calculation by calling sequentially:
1461 // CalcFitreg(); TrackletSelection(); FitTracklet()
1462 // and store the tracklets
1463
b0a41e80 1464 if (!fInitialized) {
ab9f7002 1465 AliError("Called uninitialized! Nothing done!");
b0a41e80 1466 return;
1d93b218 1467 }
96e6312d 1468
b0a41e80 1469 fTrackletArray->Delete();
96e6312d 1470
b0a41e80 1471 CalcFitreg();
1472 TrackletSelection();
1473 FitTracklet();
1d93b218 1474
b0a41e80 1475 AliRunLoader *rl = AliRunLoader::Instance();
1476 AliDataLoader *dl = 0x0;
1477 if (rl)
1478 dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1479 if (!dl) {
1480 AliError("Could not get the tracklets data loader!");
1d93b218 1481 }
1482 else {
b0a41e80 1483 TTree *trackletTree = dl->Tree();
1484 if (!trackletTree)
1485 dl->MakeTree();
1486 trackletTree = dl->Tree();
1d93b218 1487
b0a41e80 1488 AliTRDtrackletMCM *trkl = 0x0;
1489 TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1490 if (!trkbranch)
1491 trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
1492// trkbranch = trackletTree->Branch("mcmtrklbranch", &fTrackletArray, 32000, 2);
1d93b218 1493
b0a41e80 1494 for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
1495 trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
1496 trkbranch->SetAddress(&trkl);
1497// printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord());
1498 trkbranch->Fill();
1d93b218 1499 }
b0a41e80 1500 dl->WriteData("OVERWRITE");
1d93b218 1501 }
b0a41e80 1502}
1d93b218 1503
b0a41e80 1504void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
1505{
1506 // write back the processed data configured by EBSF
1507 // EBSF = 1: unfiltered data; EBSF = 0: filtered data
1508 // zero-suppressed valued are written as -1 to digits
1d93b218 1509
b0a41e80 1510 if (!fInitialized) {
ab9f7002 1511 AliError("Called uninitialized! Nothing done!");
b0a41e80 1512 return;
1d93b218 1513 }
1d93b218 1514
b0a41e80 1515 Int_t firstAdc = 0;
1516 Int_t lastAdc = fNADC - 1;
1d93b218 1517
b0a41e80 1518 if (GetCol(firstAdc) > 143)
1519 firstAdc = 1;
1d93b218 1520
b0a41e80 1521 if (GetCol(lastAdc) < 0)
1522 lastAdc = fNADC - 2;
1d93b218 1523
b0a41e80 1524 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1525 {
1526 for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
1527 if (fZSM1Dim[iAdc] == 1) {
1528 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1529 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1);
1530// printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin,
1531// digits->GetData(GetRow(), GetCol(iAdc), iTimeBin));
1532 }
1d93b218 1533 }
b0a41e80 1534 }
1d93b218 1535 }
52c19022 1536 else {
b0a41e80 1537 for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
1538 if (fZSM1Dim[iAdc] == 0) {
1539 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1540 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
1541 }
1542 }
1543 else {
1544 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1545 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1);
1546// printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin);
1547 }
1548 }
52c19022 1549 }
52c19022 1550 }
b0a41e80 1551}
23200400 1552
b0a41e80 1553// help functions, to be cleaned up
1d93b218 1554
ab9f7002 1555UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
b0a41e80 1556{
1557 //
1558 // This function adds a and b (unsigned) and clips to
1559 // the specified number of bits.
1560 //
1561
1562 UInt_t sum = a + b;
1563 if (nbits < 32)
1564 {
1565 UInt_t maxv = (1 << nbits) - 1;;
1566 if (sum > maxv)
1567 sum = maxv;
1568 }
1569 else
1570 {
1571 if ((sum < a) || (sum < b))
1572 sum = 0xFFFFFFFF;
1573 }
1574 return sum;
1d93b218 1575}
b0a41e80 1576
1577void AliTRDmcmSim::Sort2(uint16_t idx1i, uint16_t idx2i, \
1578 uint16_t val1i, uint16_t val2i, \
1579 uint16_t *idx1o, uint16_t *idx2o, \
ab9f7002 1580 uint16_t *val1o, uint16_t *val2o) const
b65e5048 1581{
ab9f7002 1582 // sorting for tracklet selection
b65e5048 1583
b0a41e80 1584 if (val1i > val2i)
b65e5048 1585 {
b0a41e80 1586 *idx1o = idx1i;
1587 *idx2o = idx2i;
1588 *val1o = val1i;
1589 *val2o = val2i;
b65e5048 1590 }
b0a41e80 1591 else
b65e5048 1592 {
b0a41e80 1593 *idx1o = idx2i;
1594 *idx2o = idx1i;
1595 *val1o = val2i;
1596 *val2o = val1i;
b65e5048 1597 }
1598}
b65e5048 1599
b0a41e80 1600void AliTRDmcmSim::Sort3(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, \
1601 uint16_t val1i, uint16_t val2i, uint16_t val3i, \
1602 uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, \
1603 uint16_t *val1o, uint16_t *val2o, uint16_t *val3o)
b65e5048 1604{
ab9f7002 1605 // sorting for tracklet selection
1606
b0a41e80 1607 int sel;
b65e5048 1608
b65e5048 1609
b0a41e80 1610 if (val1i > val2i) sel=4; else sel=0;
1611 if (val2i > val3i) sel=sel + 2;
1612 if (val3i > val1i) sel=sel + 1;
1613 //printf("input channels %d %d %d, charges %d %d %d sel=%d\n",idx1i, idx2i, idx3i, val1i, val2i, val3i, sel);
1614 switch(sel)
b65e5048 1615 {
b0a41e80 1616 case 6 : // 1 > 2 > 3 => 1 2 3
1617 case 0 : // 1 = 2 = 3 => 1 2 3 : in this case doesn't matter, but so is in hardware!
1618 *idx1o = idx1i;
1619 *idx2o = idx2i;
1620 *idx3o = idx3i;
1621 *val1o = val1i;
1622 *val2o = val2i;
1623 *val3o = val3i;
1624 break;
1625
1626 case 4 : // 1 > 2, 2 <= 3, 3 <= 1 => 1 3 2
1627 *idx1o = idx1i;
1628 *idx2o = idx3i;
1629 *idx3o = idx2i;
1630 *val1o = val1i;
1631 *val2o = val3i;
1632 *val3o = val2i;
1633 break;
1634
1635 case 2 : // 1 <= 2, 2 > 3, 3 <= 1 => 2 1 3
1636 *idx1o = idx2i;
1637 *idx2o = idx1i;
1638 *idx3o = idx3i;
1639 *val1o = val2i;
1640 *val2o = val1i;
1641 *val3o = val3i;
1642 break;
1643
1644 case 3 : // 1 <= 2, 2 > 3, 3 > 1 => 2 3 1
1645 *idx1o = idx2i;
1646 *idx2o = idx3i;
1647 *idx3o = idx1i;
1648 *val1o = val2i;
1649 *val2o = val3i;
1650 *val3o = val1i;
1651 break;
1652
1653 case 1 : // 1 <= 2, 2 <= 3, 3 > 1 => 3 2 1
1654 *idx1o = idx3i;
1655 *idx2o = idx2i;
1656 *idx3o = idx1i;
1657 *val1o = val3i;
1658 *val2o = val2i;
1659 *val3o = val1i;
1660 break;
1661
1662 case 5 : // 1 > 2, 2 <= 3, 3 > 1 => 3 1 2
1663 *idx1o = idx3i;
1664 *idx2o = idx1i;
1665 *idx3o = idx2i;
1666 *val1o = val3i;
1667 *val2o = val1i;
1668 *val3o = val2i;
1669 break;
1670
1671 default: // the rest should NEVER happen!
1672 printf("ERROR in Sort3!!!\n");
1673 break;
1674 }
1675// printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
b65e5048 1676}
b0a41e80 1677
1678void AliTRDmcmSim::Sort6To4(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \
1679 uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \
1680 uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, uint16_t *idx4o, \
1681 uint16_t *val1o, uint16_t *val2o, uint16_t *val3o, uint16_t *val4o)
b65e5048 1682{
ab9f7002 1683 // sorting for tracklet selection
b65e5048 1684
b0a41e80 1685 uint16_t idx21s, idx22s, idx23s, dummy;
1686 uint16_t val21s, val22s, val23s;
1687 uint16_t idx23as, idx23bs;
1688 uint16_t val23as, val23bs;
1689
1690 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1691 idx1o, &idx21s, &idx23as,
1692 val1o, &val21s, &val23as);
1693
1694 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1695 idx2o, &idx22s, &idx23bs,
1696 val2o, &val22s, &val23bs);
1697
1698 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
1699
1700 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1701 idx3o, idx4o, &dummy,
1702 val3o, val4o, &dummy);
1703
b65e5048 1704}
b0a41e80 1705
1706void AliTRDmcmSim::Sort6To2Worst(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \
1707 uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \
1708 uint16_t *idx5o, uint16_t *idx6o)
0d64b05f 1709{
ab9f7002 1710 // sorting for tracklet selection
0d64b05f 1711
b0a41e80 1712 uint16_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
1713 uint16_t val21s, val22s, val23s;
1714 uint16_t idx23as, idx23bs;
1715 uint16_t val23as, val23bs;
1716
1717 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1718 &dummy1, &idx21s, &idx23as,
1719 &dummy2, &val21s, &val23as);
1720
1721 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1722 &dummy1, &idx22s, &idx23bs,
1723 &dummy2, &val22s, &val23bs);
1724
1725 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
b65e5048 1726
b0a41e80 1727 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1728 &dummy1, &dummy2, idx6o,
1729 &dummy3, &dummy4, &dummy5);
1730// printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n",
1731// idx21s, idx23as, idx22s, idx23bs, *idx5o, *idx6o);
0d64b05f 1732}
1d93b218 1733