Always export to FES the regional configuration file
[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++) {
a6d08b7f 230 padcol = GetCol(ch);
b0a41e80 231 for (Int_t tb = 0; tb < fNTimeBin; tb++) {
b0a41e80 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
a6d08b7f 497 while (GetCol(firstAdc) < 0) {
b0a41e80 498 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
499 fADCR[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
500 fADCF[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
501 }
a6d08b7f 502 firstAdc++;
b0a41e80 503 }
504
a6d08b7f 505 while (GetCol(lastAdc) < 0) {
b0a41e80 506 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
507 fADCR[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
508 fADCF[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
509 }
a6d08b7f 510 lastAdc--;
b0a41e80 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
a6d08b7f 555 Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
556 if (col < 0 || col >= fFeeParam->GetNcol())
557 return -1;
558 else
559 return col;
dfd03fc3 560}
561
b0a41e80 562Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
987ba9a3 563{
564 //
565 // Produce raw data stream from this MCM and put in buf
566 // Returns number of words filled, or negative value
567 // with -1 * number of overflowed words
568 //
569
570 UInt_t x;
987ba9a3 571 Int_t nw = 0; // Number of written words
572 Int_t of = 0; // Number of overflowed words
573 Int_t rawVer = fFeeParam->GetRAWversion();
574 Int_t **adc;
575 Int_t nActiveADC = 0; // number of activated ADC bits in a word
576
577 if( !CheckInitialized() ) return 0;
578
579 if( fFeeParam->GetRAWstoreRaw() ) {
580 adc = fADCR;
581 } else {
582 adc = fADCF;
583 }
584
b0a41e80 585 // Produce MCM header
6a04e92b 586 x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
b0a41e80 587
987ba9a3 588 if (nw < maxSize) {
589 buf[nw++] = x;
590 //printf("\nMCM header: %X ",x);
591 }
592 else {
593 of++;
594 }
595
596 // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
597 // n : unused , c : ADC count, m : selected ADCs
598 if( rawVer >= 3 ) {
599 x = 0;
600 for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
601 if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed
602 x = x | (1 << (iAdc+4) ); // last 4 digit reserved for 1100=0xc
603 nActiveADC++; // number of 1 in mmm....m
604 }
605 }
606 x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC; // nn = 01, ccccc are inverted, 0xc=1100
607 //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
608
609 if (nw < maxSize) {
610 buf[nw++] = x;
611 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
612 }
613 else {
614 of++;
615 }
616 }
617
618 // Produce ADC data. 3 timebins are packed into one 32 bits word
619 // In this version, different ADC channel will NOT share the same word
620
621 UInt_t aa=0, a1=0, a2=0, a3=0;
622
623 for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
624 if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
b0a41e80 625 aa = !(iAdc & 1) + 2;
987ba9a3 626 for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
b0a41e80 627 a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] >> fgkAddDigits : 0;
628 a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] >> fgkAddDigits : 0;
629 a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] >> fgkAddDigits : 0;
987ba9a3 630 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
631 if (nw < maxSize) {
b0a41e80 632 buf[nw++] = x;
633 //printf("%08X ",x);
987ba9a3 634 }
635 else {
b0a41e80 636 of++;
987ba9a3 637 }
638 }
639 }
640
641 if( of != 0 ) return -of; else return nw;
642}
643
1d93b218 644Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
645{
646 //
647 // Produce tracklet data stream from this MCM and put in buf
648 // Returns number of words filled, or negative value
649 // with -1 * number of overflowed words
650 //
651
652 UInt_t x;
653 Int_t nw = 0; // Number of written words
654 Int_t of = 0; // Number of overflowed words
655
656 if( !CheckInitialized() ) return 0;
657
658 // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM
659 // fMCMT is filled continuously until no more tracklet words available
660
661 Int_t wd = 0;
662 while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){
663 x = fMCMT[wd];
664 if (nw < maxSize) {
665 buf[nw++] = x;
666 }
667 else {
668 of++;
669 }
670 wd++;
671 }
672
673 if( of != 0 ) return -of; else return nw;
674}
675
dfd03fc3 676void AliTRDmcmSim::Filter()
677{
0c349049 678 //
b0a41e80 679 // Filter the raw ADC values. The active filter stages and their
680 // parameters are taken from AliTRDtrapConfig.
681 // The raw data is stored separate from the filtered data. Thus,
682 // it is possible to run the filters on a set of raw values
683 // sequentially for parameter tuning.
0c349049 684 //
dfd03fc3 685
b0a41e80 686 if( !CheckInitialized() ) {
687 AliError("got called before initialization! Nothing done!");
688 return;
dfd03fc3 689 }
690
b0a41e80 691 // Apply filters sequentially. Bypass is handled by filters
692 // since counters and internal registers may be updated even
693 // if the filter is bypassed.
694 // The first filter takes the data from fADCR and
695 // outputs to fADCF.
696
697 // Non-linearity filter not implemented.
698 FilterPedestal();
699 FilterGain();
700 FilterTail();
701 // Crosstalk filter not implemented.
dfd03fc3 702}
703
b0a41e80 704void AliTRDmcmSim::FilterPedestalInit()
705{
706 // Initializes the pedestal filter assuming that the input has
707 // been constant for a long time (compared to the time constant).
708
709// UShort_t fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
710 UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
711 UShort_t shifts[4] = {11, 14, 17, 21}; //??? where to take shifts from?
712
713 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
714 fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
715}
716
717UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
dfd03fc3 718{
b0a41e80 719 // Returns the output of the pedestal filter given the input value.
720 // The output depends on the internal registers and, thus, the
721 // history of the filter.
722
723 UShort_t fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
724 UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
725 UShort_t fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY); // 0..1 the bypass, active low
726 UShort_t shifts[4] = {11, 14, 17, 21}; //??? where to come from
727
728 UShort_t accumulatorShifted;
729 Int_t correction;
730 UShort_t inpAdd;
731
732 inpAdd = value + fpnp;
733
734 if (fpby == 0) //??? before or after update of accumulator
735 return value;
736
737 accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF; // 10 bits
738 if (timebin == 0) // the accumulator is disabled in the drift time
739 {
740 correction = (value & 0x3FF) - accumulatorShifted;
741 fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF; // 31 bits
742 }
743
744 if (inpAdd <= accumulatorShifted)
745 return 0;
746 else
747 {
748 inpAdd = inpAdd - accumulatorShifted;
749 if (inpAdd > 0xFFF)
750 return 0xFFF;
751 else
752 return inpAdd;
753 }
754}
23200400 755
b0a41e80 756void AliTRDmcmSim::FilterPedestal()
757{
0c349049 758 //
dfd03fc3 759 // Apply pedestal filter
0c349049 760 //
b0a41e80 761 // As the first filter in the chain it reads data from fADCR
762 // and outputs to fADCF.
763 // It has only an effect if previous samples have been fed to
764 // find the pedestal. Currently, the simulation assumes that
765 // the input has been stable for a sufficiently long time.
dfd03fc3 766
b0a41e80 767 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
768 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
769 fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]);
dfd03fc3 770 }
771 }
772}
773
b0a41e80 774void AliTRDmcmSim::FilterGainInit()
dfd03fc3 775{
b0a41e80 776 // Initializes the gain filter. In this case, only threshold
777 // counters are reset.
0c349049 778
b0a41e80 779 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
780 // these are counters which in hardware continue
781 // until maximum or reset
782 fGainCounterA[iAdc] = 0;
783 fGainCounterB[iAdc] = 0;
784 }
dfd03fc3 785}
786
b0a41e80 787UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
dfd03fc3 788{
b0a41e80 789 // Apply the gain filter to the given value.
790 // BEGIN_LATEX O_{i}(t) = #gamma_{i} * I_{i}(t) + a_{i} END_LATEX
791 // The output depends on the internal registers and, thus, the
792 // history of the filter.
dfd03fc3 793
b0a41e80 794 UShort_t fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY); // bypass, active low
795 UShort_t fgf = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc)); // 0x700 + (0 & 0x1ff);
796 UShort_t fga = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc)); // 40;
797 UShort_t fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA); // 20;
798 UShort_t fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB); // 2060;
dfd03fc3 799
b0a41e80 800 UInt_t tmp;
dfd03fc3 801
b0a41e80 802 value &= 0xFFF;
803 tmp = (value * fgf) >> 11;
804 if (tmp > 0xFFF) tmp = 0xFFF;
805
806 if (fgby == 1)
807 value = AddUintClipping(tmp, fga, 12);
808
809 // Update threshold counters
810 // not really useful as they are cleared with every new event
811 if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))
812 {
813 if (value >= fgtb)
814 fGainCounterB[adc]++;
815 else if (value >= fgta)
816 fGainCounterA[adc]++;
817 }
818
819 return value;
820}
821
822void AliTRDmcmSim::FilterGain()
823{
824 // Read data from fADCF and apply gain filter.
825
826 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
827 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
828 fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
1d93b218 829 }
b0a41e80 830 }
831}
1d93b218 832
b0a41e80 833void AliTRDmcmSim::FilterTailInit(Int_t baseline)
834{
835 // Initializes the tail filter assuming that the input has
836 // been at the baseline value (configured by FTFP) for a
837 // sufficiently long time.
838
839 // exponents and weight calculated from configuration
840 UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
841 UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
842 UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
843
844 Float_t lambdaL = lambdaLong * 1.0 / (1 << 11);
845 Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
846 Float_t alphaL = alphaLong * 1.0 / (1 << 11);
847 Float_t qup, qdn;
848 qup = (1 - lambdaL) * (1 - lambdaS);
849 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
850 Float_t kdc = qup/qdn;
851
852 Float_t kt, ql, qs;
853 UShort_t aout;
854
855 kt = kdc * baseline;
856 aout = baseline - (UShort_t) kt;
857 ql = lambdaL * (1 - lambdaS) * alphaL;
858 qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
859
860 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
861 fTailAmplLong[iAdc] = (UShort_t) (aout * ql / (ql + qs));
862 fTailAmplShort[iAdc] = (UShort_t) (aout * qs / (ql + qs));
863 }
864}
865
866UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
867{
868 // Returns the output of the tail filter for the given input value.
869 // The output depends on the internal registers and, thus, the
870 // history of the filter.
871
872 // exponents and weight calculated from configuration
873 UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
874 UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
875 UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
876
877 Float_t lambdaL = lambdaLong * 1.0 / (1 << 11);
878 Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
879 Float_t alphaL = alphaLong * 1.0 / (1 << 11);
880 Float_t qup, qdn;
881 qup = (1 - lambdaL) * (1 - lambdaS);
882 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
883// Float_t kdc = qup/qdn;
884
885 UInt_t aDiff;
886 UInt_t alInpv;
887 UShort_t aQ;
888 UInt_t tmp;
889
ab9f7002 890 UShort_t inpVolt = value & 0xFFF; // 12 bits
b0a41e80 891
892 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
893 return value;
894 else
895 {
896 // add the present generator outputs
897 aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
898
899 // calculate the difference between the input the generated signal
ab9f7002 900 if (inpVolt > aQ)
901 aDiff = inpVolt - aQ;
b0a41e80 902 else
903 aDiff = 0;
904
905 // the inputs to the two generators, weighted
906 alInpv = (aDiff * alphaLong) >> 11;
907
908 // the new values of the registers, used next time
909 // long component
910 tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
911 tmp = (tmp * lambdaLong) >> 11;
912 fTailAmplLong[adc] = tmp & 0xFFF;
913 // short component
914 tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
915 tmp = (tmp * lambdaShort) >> 11;
916 fTailAmplShort[adc] = tmp & 0xFFF;
917
918 // the output of the filter
919 return aDiff;
dfd03fc3 920 }
b0a41e80 921}
dfd03fc3 922
b0a41e80 923void AliTRDmcmSim::FilterTail()
924{
925 // Apply tail filter
fabf2e09 926
b0a41e80 927 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
928 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
929 fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]);
930 }
931 }
dfd03fc3 932}
933
dfd03fc3 934void AliTRDmcmSim::ZSMapping()
935{
0c349049 936 //
dfd03fc3 937 // Zero Suppression Mapping implemented in TRAP chip
938 //
939 // See detail TRAP manual "Data Indication" section:
940 // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
0c349049 941 //
dfd03fc3 942
b0a41e80 943 //??? values should come from TRAPconfig
944 Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS); // TRAP default = 0x4 (Tis=4)
945 Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT); // TRAP default = 0x28 (Tit=40)
946 Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL); // TRAP default = 0xf0
947 // (lookup table accept (I2,I1,I0)=(111)
948 // or (110) or (101) or (100))
949 Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); // TRAP default = 1 (no neighbor sensitivity)
ecf39416 950 Int_t ep = AliTRDfeeParam::GetPFeffectPedestal();
951
b0a41e80 952 Int_t **adc = fADCF;
dfd03fc3 953
b0a41e80 954 if( !CheckInitialized() ) {
ab9f7002 955 AliError("got called uninitialized! Nothing done!");
b0a41e80 956 return;
957 }
958
959 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
960 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
dfd03fc3 961
ecf39416 962 // Get ADC data currently in filter buffer
b0a41e80 963 Int_t ap = adc[iadc-1][it] - ep; // previous
964 Int_t ac = adc[iadc ][it] - ep; // current
965 Int_t an = adc[iadc+1][it] - ep; // next
ecf39416 966
dfd03fc3 967 // evaluate three conditions
0c349049 968 Int_t i0 = ( ac >= ap && ac >= an ) ? 0 : 1; // peak center detection
969 Int_t i1 = ( ap + ac + an > eBIT ) ? 0 : 1; // cluster
970 Int_t i2 = ( ac > eBIS ) ? 0 : 1; // absolute large peak
971
972 Int_t i = i2 * 4 + i1 * 2 + i0; // Bit position in lookup table
973 Int_t d = (eBIL >> i) & 1; // Looking up (here d=0 means true
974 // and d=1 means false according to TRAP manual)
975
976 fZSM[iadc][it] &= d;
977 if( eBIN == 0 ) { // turn on neighboring ADCs
978 fZSM[iadc-1][it] &= d;
979 fZSM[iadc+1][it] &= d;
dfd03fc3 980 }
981 }
982 }
983
984 // do 1 dim projection
985 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
986 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
ecf39416 987 fZSM1Dim[iadc] &= fZSM[iadc][it];
988 }
989 }
0c349049 990
ecf39416 991}
992
ecf39416 993void AliTRDmcmSim::DumpData( char *f, char *target )
994{
0c349049 995 //
ecf39416 996 // Dump data stored (for debugging).
997 // target should contain one or multiple of the following characters
998 // R for raw data
999 // F for filtered data
1000 // Z for zero suppression map
1001 // S Raw dat astream
1002 // other characters are simply ignored
0c349049 1003 //
1004
ecf39416 1005 UInt_t tempbuf[1024];
1006
1007 if( !CheckInitialized() ) return;
1008
1009 std::ofstream of( f, std::ios::out | std::ios::app );
1010 of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
b0a41e80 1011 fDetector, fGeo->GetSector(fDetector), fGeo->GetStack(fDetector),
1012 fGeo->GetSector(fDetector), fRobPos, fMcmPos );
ecf39416 1013
1014 for( int t=0 ; target[t] != 0 ; t++ ) {
1015 switch( target[t] ) {
1016 case 'R' :
1017 case 'r' :
1018 of << Form("fADCR (raw ADC data)\n");
1019 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1020 of << Form(" ADC %02d: ", iadc);
1021 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1022 of << Form("% 4d", fADCR[iadc][it]);
1023 }
1024 of << Form("\n");
1025 }
1026 break;
1027 case 'F' :
1028 case 'f' :
1029 of << Form("fADCF (filtered ADC data)\n");
1030 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1031 of << Form(" ADC %02d: ", iadc);
1032 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1033 of << Form("% 4d", fADCF[iadc][it]);
1034 }
1035 of << Form("\n");
1036 }
1037 break;
1038 case 'Z' :
1039 case 'z' :
1040 of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
1041 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1042 of << Form(" ADC %02d: ", iadc);
1043 if( fZSM1Dim[iadc] == 0 ) { of << " R " ; } else { of << " . "; } // R:read .:suppressed
1044 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1045 if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
1046 }
1047 of << Form("\n");
1048 }
1049 break;
1050 case 'S' :
1051 case 's' :
1052 Int_t s = ProduceRawStream( tempbuf, 1024 );
1053 of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
1054 of << Form(" address data\n");
1055 for( int i = 0 ; i < s ; i++ ) {
1056 of << Form(" %04x %08x\n", i, tempbuf[i]);
1057 }
dfd03fc3 1058 }
1059 }
1060}
1061
b0a41e80 1062void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label)
dfd03fc3 1063{
b0a41e80 1064 // Add the given hit to the fit register which is lateron used for
1065 // the tracklet calculation.
1066 // In addition to the fit sums in the fit register MC information
1067 // is stored.
1068
1069 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) &&
1070 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
ab9f7002 1071 fFitReg[adc].fQ0 += qtot;
b0a41e80 1072
1073 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) &&
1074 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
ab9f7002 1075 fFitReg[adc].fQ1 += qtot;
b0a41e80 1076
1077 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) &&
1078 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
1079 {
ab9f7002 1080 fFitReg[adc].fSumX += timebin;
1081 fFitReg[adc].fSumX2 += timebin*timebin;
1082 fFitReg[adc].fNhits++;
1083 fFitReg[adc].fSumY += ypos;
1084 fFitReg[adc].fSumY2 += ypos*ypos;
1085 fFitReg[adc].fSumXY += timebin*ypos;
b0a41e80 1086 }
1087
1088 // register hits (MC info)
ab9f7002 1089 fHits[fNHits].fChannel = adc;
1090 fHits[fNHits].fQtot = qtot;
1091 fHits[fNHits].fYpos = ypos;
1092 fHits[fNHits].fTimebin = timebin;
1093 fHits[fNHits].fLabel = label;
b0a41e80 1094 fNHits++;
dfd03fc3 1095}
1096
b0a41e80 1097void AliTRDmcmSim::CalcFitreg()
dfd03fc3 1098{
b0a41e80 1099 // Preprocessing.
1100 // Detect the hits and fill the fit registers.
1101 // Requires 12-bit data from fADCF which means Filter()
1102 // has to be called before even if all filters are bypassed.
1103
1104 //???
1105 // TRAP parameters:
ab9f7002 1106 const uint16_t lutPos[128] = { // move later to some other file
b0a41e80 1107 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,
1108 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,
1109 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,
1110 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};
1111
1112 //??? to be clarified:
ab9f7002 1113 UInt_t adcMask = 0xfffff;
b0a41e80 1114
ab9f7002 1115 UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
b0a41e80 1116 Short_t ypos, fromLeft, fromRight, found;
ab9f7002 1117 UShort_t qTotal[19]; // the last is dummy
1118 UShort_t marked[6], qMarked[6], worse1, worse2;
b0a41e80 1119
1120 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS);
1121 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)
1122 < timebin1)
1123 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1124 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE);
1125 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)
1126 > timebin2)
1127 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1128
1129 // reset the fit registers
1130 fNHits = 0;
1131 for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels
1132 {
ab9f7002 1133 fFitReg[adcch].fNhits = 0;
1134 fFitReg[adcch].fQ0 = 0;
1135 fFitReg[adcch].fQ1 = 0;
1136 fFitReg[adcch].fSumX = 0;
1137 fFitReg[adcch].fSumY = 0;
1138 fFitReg[adcch].fSumX2 = 0;
1139 fFitReg[adcch].fSumY2 = 0;
1140 fFitReg[adcch].fSumXY = 0;
b0a41e80 1141 }
1142
1143 for (timebin = timebin1; timebin < timebin2; timebin++)
1144 {
ab9f7002 1145 // first find the hit candidates and store the total cluster charge in qTotal array
b0a41e80 1146 // in case of not hit store 0 there.
1147 for (adcch = 0; adcch < fNADC-2; adcch++) {
ab9f7002 1148 if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
b0a41e80 1149 {
ab9f7002 1150 adcLeft = fADCF[adcch ][timebin];
1151 adcCentral = fADCF[adcch+1][timebin];
1152 adcRight = fADCF[adcch+2][timebin];
b0a41e80 1153 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1)
ab9f7002 1154 hitQual = ( (adcLeft * adcRight) <
1155 (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) );
b0a41e80 1156 else
ab9f7002 1157 hitQual = 1;
b0a41e80 1158 // The accumulated charge is with the pedestal!!!
ab9f7002 1159 qtotTemp = adcLeft + adcCentral + adcRight;
1160 if ( (hitQual) &&
1161 (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1162 (adcLeft <= adcCentral) &&
1163 (adcCentral > adcRight) )
1164 qTotal[adcch] = qtotTemp;
b0a41e80 1165 else
ab9f7002 1166 qTotal[adcch] = 0;
1167 //printf("ch %2d qTotal %5d\n",adcch, qTotal[adcch]);
1d93b218 1168 }
b0a41e80 1169 else
ab9f7002 1170 qTotal[adcch] = 0; //jkl
1d93b218 1171 }
b0a41e80 1172
1173 fromLeft = -1;
1174 adcch = 0;
1175 found = 0;
1176 marked[4] = 19; // invalid channel
1177 marked[5] = 19; // invalid channel
ab9f7002 1178 qTotal[19] = 0;
b0a41e80 1179 while ((adcch < 16) && (found < 3))
1180 {
ab9f7002 1181 if (qTotal[adcch] > 0)
b0a41e80 1182 {
1183 fromLeft = adcch;
1184 marked[2*found+1]=adcch;
1185 found++;
1d93b218 1186 }
b0a41e80 1187 adcch++;
1d93b218 1188 }
1d93b218 1189
b0a41e80 1190 fromRight = -1;
1191 adcch = 18;
1192 found = 0;
1193 while ((adcch > 2) && (found < 3))
1194 {
ab9f7002 1195 if (qTotal[adcch] > 0)
b0a41e80 1196 {
1197 marked[2*found]=adcch;
1198 found++;
1199 fromRight = adcch;
1200 }
1201 adcch--;
1d93b218 1202 }
1d93b218 1203
b0a41e80 1204 //printf("Fromleft=%d, Fromright=%d\n",fromLeft, fromRight);
1205 // here mask the hit candidates in the middle, if any
1206 if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight))
1207 for (adcch = fromLeft+1; adcch < fromRight; adcch++)
ab9f7002 1208 qTotal[adcch] = 0;
1d93b218 1209
b0a41e80 1210 found = 0;
1211 for (adcch = 0; adcch < 19; adcch++)
ab9f7002 1212 if (qTotal[adcch] > 0) found++;
b0a41e80 1213 // NOT READY
1d93b218 1214
b0a41e80 1215 if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1216 {
1217 if (marked[4] == marked[5]) marked[5] = 19;
1218 for (found=0; found<6; found++)
1219 {
ab9f7002 1220 qMarked[found] = qTotal[marked[found]] >> 4;
1221 //printf("ch_%d qTotal %d qTotals %d |",marked[found],qTotal[marked[found]],qMarked[found]);
1d93b218 1222 }
b0a41e80 1223 //printf("\n");
1224
1225 Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
ab9f7002 1226 qMarked[0],
1227 qMarked[3],
1228 qMarked[4],
1229 qMarked[1],
1230 qMarked[2],
1231 qMarked[5],
b0a41e80 1232 &worse1, &worse2);
1233 // Now mask the two channels with the smallest charge
1234 if (worse1 < 19)
1235 {
ab9f7002 1236 qTotal[worse1] = 0;
b0a41e80 1237 //printf("Kill ch %d\n",worse1);
1d93b218 1238 }
b0a41e80 1239 if (worse2 < 19)
1240 {
ab9f7002 1241 qTotal[worse2] = 0;
b0a41e80 1242 //printf("Kill ch %d\n",worse2);
1d93b218 1243 }
1d93b218 1244 }
1d93b218 1245
b0a41e80 1246 for (adcch = 0; adcch < 19; adcch++) {
ab9f7002 1247 if (qTotal[adcch] > 0) // the channel is marked for processing
b0a41e80 1248 {
ab9f7002 1249 adcLeft = fADCF[adcch ][timebin];
1250 adcCentral = fADCF[adcch+1][timebin];
1251 adcRight = fADCF[adcch+2][timebin];
b0a41e80 1252 // hit detected, in TRAP we have 4 units and a hit-selection, here we proceed all channels!
1253 // subtract the pedestal TPFP, clipping instead of wrapping
1254
ab9f7002 1255 Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP);
1256// printf("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n",
1257// timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP,
b0a41e80 1258// fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT));
1259
ab9f7002 1260 if (adcLeft < regTPFP) adcLeft = 0; else adcLeft -= regTPFP;
1261 if (adcCentral < regTPFP) adcCentral = 0; else adcCentral -= regTPFP;
1262 if (adcRight < regTPFP) adcRight = 0; else adcRight -= regTPFP;
b0a41e80 1263 // Calculate the center of gravity
ab9f7002 1264 ypos = 128*(adcLeft - adcRight) / adcCentral;
b0a41e80 1265 if (ypos < 0) ypos = -ypos;
1266 // make the correction using the LUT
ab9f7002 1267 ypos = ypos + lutPos[ypos & 0x7F];
1268 if (adcLeft > adcRight) ypos = -ypos;
1269 AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, -1);
1d93b218 1270 }
1d93b218 1271 }
1d93b218 1272 }
b0a41e80 1273}
1d93b218 1274
b0a41e80 1275void AliTRDmcmSim::TrackletSelection()
1276{
1277 // Select up to 4 tracklet candidates from the fit registers
1278 // and assign them to the CPUs.
1279
ab9f7002 1280 UShort_t adcIdx, i, j, ntracks, tmp;
1281 UShort_t trackletCand[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates
b0a41e80 1282
1283 ntracks = 0;
ab9f7002 1284 for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs
1285 if ( (fFitReg[adcIdx].fNhits
b0a41e80 1286 >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) &&
ab9f7002 1287 (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits
b0a41e80 1288 >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT)))
1289 {
ab9f7002 1290 trackletCand[ntracks][0] = adcIdx;
1291 trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits;
1292 //printf("%d %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]);
b0a41e80 1293 ntracks++;
1294 };
1295
ab9f7002 1296 // for (i=0; i<ntracks;i++) printf("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]);
b0a41e80 1297
1298 if (ntracks > 4)
1299 {
1300 // primitive sorting according to the number of hits
1301 for (j = 0; j < (ntracks-1); j++)
1302 {
1303 for (i = j+1; i < ntracks; i++)
1304 {
ab9f7002 1305 if ( (trackletCand[j][1] < trackletCand[i][1]) ||
1306 ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
b0a41e80 1307 {
1308 // swap j & i
ab9f7002 1309 tmp = trackletCand[j][1];
1310 trackletCand[j][1] = trackletCand[i][1];
1311 trackletCand[i][1] = tmp;
1312 tmp = trackletCand[j][0];
1313 trackletCand[j][0] = trackletCand[i][0];
1314 trackletCand[i][0] = tmp;
b0a41e80 1315 }
1d93b218 1316 }
1d93b218 1317 }
b0a41e80 1318 ntracks = 4; // cut the rest, 4 is the max
1d93b218 1319 }
b0a41e80 1320 // else is not necessary to sort
1d93b218 1321
b0a41e80 1322 // now sort, so that the first tracklet going to CPU0 corresponds to the highest adc channel - as in the TRAP
1323 for (j = 0; j < (ntracks-1); j++)
1324 {
1325 for (i = j+1; i < ntracks; i++)
1326 {
ab9f7002 1327 if (trackletCand[j][0] < trackletCand[i][0])
b0a41e80 1328 {
1329 // swap j & i
ab9f7002 1330 tmp = trackletCand[j][1];
1331 trackletCand[j][1] = trackletCand[i][1];
1332 trackletCand[i][1] = tmp;
1333 tmp = trackletCand[j][0];
1334 trackletCand[j][0] = trackletCand[i][0];
1335 trackletCand[i][0] = tmp;
1d93b218 1336 }
1d93b218 1337 }
1338 }
b0a41e80 1339 for (i = 0; i < ntracks; i++) // CPUs with tracklets.
ab9f7002 1340 fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i]
b0a41e80 1341 for (i = ntracks; i < 4; i++) // CPUs without tracklets
1342 fFitPtr[i] = 31; // pointer to the left channel with tracklet for CPU[i] = 31 (invalid)
1343// printf("found %i tracklet candidates\n", ntracks);
1344}
1d93b218 1345
b0a41e80 1346void AliTRDmcmSim::FitTracklet()
1347{
1348 // Perform the actual tracklet fit based on the fit sums
1349 // which have been filled in the fit registers.
1350
1351 // parameters in fitred.asm (fit program)
1352 Int_t decPlaces = 5;
1353 Int_t rndAdd = 0;
1354 if (decPlaces > 1)
1355 rndAdd = (1 << (decPlaces-1)) + 1;
1356 else if (decPlaces == 1)
1357 rndAdd = 1;
1358
1359 // should come from trapConfig (DMEM)
1360 AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
1361 Long64_t shift = ((Long64_t) 1 << 32);
ab9f7002 1362 UInt_t scaleY = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4)));
1363 UInt_t scaleD = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4)));
b0a41e80 1364 int padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos);
1365 int yoffs = (fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19) - fFeeParam->GetNcol()/2) << (8 + decPlaces);
1366 int ndrift = 20; //??? value in simulation?
ab9f7002 1367 int deflCorr = 0; // -370;
b0a41e80 1368 int minslope = -10000; // no pt-cut so far
1369 int maxslope = 10000; // no pt-cut so far
1370
1371 // local variables for calculation
1372 Long64_t mult, temp, denom; //???
1373 UInt_t q0, q1, qTotal; // charges in the two windows and total charge
1374 UShort_t nHits; // number of hits
1375 Int_t slope, offset; // slope and offset of the tracklet
1376 Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers
1377 //int32_t SumY2; // not used in the current TRAP program
1378 FitReg_t *fit0, *fit1; // pointers to relevant fit registers
1379
1380// const uint32_t OneDivN[32] = { // 2**31/N : exactly like in the TRAP, the simple division here gives the same result!
1381// 0x00000000, 0x80000000, 0x40000000, 0x2AAAAAA0, 0x20000000, 0x19999990, 0x15555550, 0x12492490,
1382// 0x10000000, 0x0E38E380, 0x0CCCCCC0, 0x0BA2E8B0, 0x0AAAAAA0, 0x09D89D80, 0x09249240, 0x08888880,
1383// 0x08000000, 0x07878780, 0x071C71C0, 0x06BCA1A0, 0x06666660, 0x06186180, 0x05D17450, 0x0590B210,
1384// 0x05555550, 0x051EB850, 0x04EC4EC0, 0x04BDA120, 0x04924920, 0x0469EE50, 0x04444440, 0x04210840};
1385
1386 for (Int_t cpu = 0; cpu < 4; cpu++) {
1387 if (fFitPtr[cpu] == 31)
1388 {
1389 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1d93b218 1390 }
b0a41e80 1391 else
1392 {
1393 fit0 = &fFitReg[fFitPtr[cpu] ];
1394 fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1395
1396 mult = 1;
1397 mult = mult << (32 + decPlaces);
1398 mult = -mult;
1399
1400 // Merging
ab9f7002 1401 nHits = fit0->fNhits + fit1->fNhits; // number of hits
1402 sumX = fit0->fSumX + fit1->fSumX;
1403 sumX2 = fit0->fSumX2 + fit1->fSumX2;
b0a41e80 1404 denom = nHits*sumX2 - sumX*sumX;
1405
1406 mult = mult / denom; // exactly like in the TRAP program
ab9f7002 1407 q0 = fit0->fQ0 + fit1->fQ0;
1408 q1 = fit0->fQ1 + fit1->fQ1;
1409 sumY = fit0->fSumY + fit1->fSumY + 256*fit1->fNhits;
1410 sumXY = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
b0a41e80 1411
1412 slope = nHits*sumXY - sumX * sumY;
1413 offset = sumX2*sumY - sumX * sumXY;
1414 temp = mult * slope;
1415 slope = temp >> 32; // take the upper 32 bits
1416 temp = mult * offset;
1417 offset = temp >> 32; // take the upper 32 bits
1418
1419 offset = offset + yoffs + (18 << (8 + decPlaces));
ab9f7002 1420 slope = slope * ndrift + deflCorr;
b0a41e80 1421 offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
1d93b218 1422
b0a41e80 1423 if ((slope < minslope) || (slope > maxslope))
1424 {
1425 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1426 }
1427 else
1428 {
1429 temp = slope;
ab9f7002 1430 temp = temp * scaleD;
b0a41e80 1431 slope = (temp >> 32);
1432
1433 temp = offset;
ab9f7002 1434 temp = temp * scaleY;
b0a41e80 1435 offset = (temp >> 32);
1436
1437 // rounding, like in the TRAP
1438 slope = (slope + rndAdd) >> decPlaces;
1439 offset = (offset + rndAdd) >> decPlaces;
1440
1441 if (slope > 0x3f || slope < -0x3f)
1442 AliWarning("Overflow in slope");
1443 slope = slope & 0x7F; // 7 bit
1444
1445 if (offset > 0xfff || offset < 0xfff)
1446 AliWarning("Overflow in offset");
1447 offset = offset & 0x1FFF; // 13 bit
1448
1449 qTotal = (q1 / nHits) >> 1;
1450 if (qTotal > 0xff)
1451 AliWarning("Overflow in charge");
1452 qTotal = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1453
1454 // assemble and store the tracklet word
1455 fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1456 new ((*fTrackletArray)[cpu]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1457 }
1d93b218 1458 }
1d93b218 1459 }
b0a41e80 1460}
1d93b218 1461
b0a41e80 1462void AliTRDmcmSim::Tracklet()
1463{
ab9f7002 1464 // Run the tracklet calculation by calling sequentially:
1465 // CalcFitreg(); TrackletSelection(); FitTracklet()
1466 // and store the tracklets
1467
b0a41e80 1468 if (!fInitialized) {
ab9f7002 1469 AliError("Called uninitialized! Nothing done!");
b0a41e80 1470 return;
1d93b218 1471 }
96e6312d 1472
b0a41e80 1473 fTrackletArray->Delete();
96e6312d 1474
b0a41e80 1475 CalcFitreg();
1476 TrackletSelection();
1477 FitTracklet();
1d93b218 1478
b0a41e80 1479 AliRunLoader *rl = AliRunLoader::Instance();
1480 AliDataLoader *dl = 0x0;
1481 if (rl)
1482 dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1483 if (!dl) {
1484 AliError("Could not get the tracklets data loader!");
1d93b218 1485 }
1486 else {
b0a41e80 1487 TTree *trackletTree = dl->Tree();
1488 if (!trackletTree)
1489 dl->MakeTree();
1490 trackletTree = dl->Tree();
1d93b218 1491
b0a41e80 1492 AliTRDtrackletMCM *trkl = 0x0;
1493 TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1494 if (!trkbranch)
1495 trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
1496// trkbranch = trackletTree->Branch("mcmtrklbranch", &fTrackletArray, 32000, 2);
1d93b218 1497
b0a41e80 1498 for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
1499 trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
1500 trkbranch->SetAddress(&trkl);
1501// printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord());
1502 trkbranch->Fill();
1d93b218 1503 }
b0a41e80 1504 dl->WriteData("OVERWRITE");
1d93b218 1505 }
b0a41e80 1506}
1d93b218 1507
b0a41e80 1508void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
1509{
1510 // write back the processed data configured by EBSF
1511 // EBSF = 1: unfiltered data; EBSF = 0: filtered data
1512 // zero-suppressed valued are written as -1 to digits
1d93b218 1513
b0a41e80 1514 if (!fInitialized) {
ab9f7002 1515 AliError("Called uninitialized! Nothing done!");
b0a41e80 1516 return;
1d93b218 1517 }
1d93b218 1518
b0a41e80 1519 Int_t firstAdc = 0;
1520 Int_t lastAdc = fNADC - 1;
1d93b218 1521
a6d08b7f 1522 while (GetCol(firstAdc) < 0)
1523 firstAdc++;
1d93b218 1524
a6d08b7f 1525 while (GetCol(lastAdc) < 0)
1526 lastAdc--;
1d93b218 1527
b0a41e80 1528 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1529 {
1530 for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
1531 if (fZSM1Dim[iAdc] == 1) {
1532 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1533 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1);
1534// printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin,
1535// digits->GetData(GetRow(), GetCol(iAdc), iTimeBin));
1536 }
1d93b218 1537 }
b0a41e80 1538 }
1d93b218 1539 }
52c19022 1540 else {
b0a41e80 1541 for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
1542 if (fZSM1Dim[iAdc] == 0) {
1543 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1544 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
1545 }
1546 }
1547 else {
1548 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1549 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1);
1550// printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin);
1551 }
1552 }
52c19022 1553 }
52c19022 1554 }
b0a41e80 1555}
23200400 1556
b0a41e80 1557// help functions, to be cleaned up
1d93b218 1558
ab9f7002 1559UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
b0a41e80 1560{
1561 //
1562 // This function adds a and b (unsigned) and clips to
1563 // the specified number of bits.
1564 //
1565
1566 UInt_t sum = a + b;
1567 if (nbits < 32)
1568 {
1569 UInt_t maxv = (1 << nbits) - 1;;
1570 if (sum > maxv)
1571 sum = maxv;
1572 }
1573 else
1574 {
1575 if ((sum < a) || (sum < b))
1576 sum = 0xFFFFFFFF;
1577 }
1578 return sum;
1d93b218 1579}
b0a41e80 1580
1581void AliTRDmcmSim::Sort2(uint16_t idx1i, uint16_t idx2i, \
1582 uint16_t val1i, uint16_t val2i, \
1583 uint16_t *idx1o, uint16_t *idx2o, \
ab9f7002 1584 uint16_t *val1o, uint16_t *val2o) const
b65e5048 1585{
ab9f7002 1586 // sorting for tracklet selection
b65e5048 1587
b0a41e80 1588 if (val1i > val2i)
b65e5048 1589 {
b0a41e80 1590 *idx1o = idx1i;
1591 *idx2o = idx2i;
1592 *val1o = val1i;
1593 *val2o = val2i;
b65e5048 1594 }
b0a41e80 1595 else
b65e5048 1596 {
b0a41e80 1597 *idx1o = idx2i;
1598 *idx2o = idx1i;
1599 *val1o = val2i;
1600 *val2o = val1i;
b65e5048 1601 }
1602}
b65e5048 1603
b0a41e80 1604void AliTRDmcmSim::Sort3(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, \
1605 uint16_t val1i, uint16_t val2i, uint16_t val3i, \
1606 uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, \
1607 uint16_t *val1o, uint16_t *val2o, uint16_t *val3o)
b65e5048 1608{
ab9f7002 1609 // sorting for tracklet selection
1610
b0a41e80 1611 int sel;
b65e5048 1612
b65e5048 1613
b0a41e80 1614 if (val1i > val2i) sel=4; else sel=0;
1615 if (val2i > val3i) sel=sel + 2;
1616 if (val3i > val1i) sel=sel + 1;
1617 //printf("input channels %d %d %d, charges %d %d %d sel=%d\n",idx1i, idx2i, idx3i, val1i, val2i, val3i, sel);
1618 switch(sel)
b65e5048 1619 {
b0a41e80 1620 case 6 : // 1 > 2 > 3 => 1 2 3
1621 case 0 : // 1 = 2 = 3 => 1 2 3 : in this case doesn't matter, but so is in hardware!
1622 *idx1o = idx1i;
1623 *idx2o = idx2i;
1624 *idx3o = idx3i;
1625 *val1o = val1i;
1626 *val2o = val2i;
1627 *val3o = val3i;
1628 break;
1629
1630 case 4 : // 1 > 2, 2 <= 3, 3 <= 1 => 1 3 2
1631 *idx1o = idx1i;
1632 *idx2o = idx3i;
1633 *idx3o = idx2i;
1634 *val1o = val1i;
1635 *val2o = val3i;
1636 *val3o = val2i;
1637 break;
1638
1639 case 2 : // 1 <= 2, 2 > 3, 3 <= 1 => 2 1 3
1640 *idx1o = idx2i;
1641 *idx2o = idx1i;
1642 *idx3o = idx3i;
1643 *val1o = val2i;
1644 *val2o = val1i;
1645 *val3o = val3i;
1646 break;
1647
1648 case 3 : // 1 <= 2, 2 > 3, 3 > 1 => 2 3 1
1649 *idx1o = idx2i;
1650 *idx2o = idx3i;
1651 *idx3o = idx1i;
1652 *val1o = val2i;
1653 *val2o = val3i;
1654 *val3o = val1i;
1655 break;
1656
1657 case 1 : // 1 <= 2, 2 <= 3, 3 > 1 => 3 2 1
1658 *idx1o = idx3i;
1659 *idx2o = idx2i;
1660 *idx3o = idx1i;
1661 *val1o = val3i;
1662 *val2o = val2i;
1663 *val3o = val1i;
1664 break;
1665
1666 case 5 : // 1 > 2, 2 <= 3, 3 > 1 => 3 1 2
1667 *idx1o = idx3i;
1668 *idx2o = idx1i;
1669 *idx3o = idx2i;
1670 *val1o = val3i;
1671 *val2o = val1i;
1672 *val3o = val2i;
1673 break;
1674
1675 default: // the rest should NEVER happen!
1676 printf("ERROR in Sort3!!!\n");
1677 break;
1678 }
1679// printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
b65e5048 1680}
b0a41e80 1681
1682void AliTRDmcmSim::Sort6To4(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \
1683 uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \
1684 uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, uint16_t *idx4o, \
1685 uint16_t *val1o, uint16_t *val2o, uint16_t *val3o, uint16_t *val4o)
b65e5048 1686{
ab9f7002 1687 // sorting for tracklet selection
b65e5048 1688
b0a41e80 1689 uint16_t idx21s, idx22s, idx23s, dummy;
1690 uint16_t val21s, val22s, val23s;
1691 uint16_t idx23as, idx23bs;
1692 uint16_t val23as, val23bs;
1693
1694 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1695 idx1o, &idx21s, &idx23as,
1696 val1o, &val21s, &val23as);
1697
1698 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1699 idx2o, &idx22s, &idx23bs,
1700 val2o, &val22s, &val23bs);
1701
1702 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
1703
1704 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1705 idx3o, idx4o, &dummy,
1706 val3o, val4o, &dummy);
1707
b65e5048 1708}
b0a41e80 1709
1710void AliTRDmcmSim::Sort6To2Worst(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \
1711 uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \
1712 uint16_t *idx5o, uint16_t *idx6o)
0d64b05f 1713{
ab9f7002 1714 // sorting for tracklet selection
0d64b05f 1715
b0a41e80 1716 uint16_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
1717 uint16_t val21s, val22s, val23s;
1718 uint16_t idx23as, idx23bs;
1719 uint16_t val23as, val23bs;
1720
1721 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1722 &dummy1, &idx21s, &idx23as,
1723 &dummy2, &val21s, &val23as);
1724
1725 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1726 &dummy1, &idx22s, &idx23bs,
1727 &dummy2, &val22s, &val23bs);
1728
1729 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
b65e5048 1730
b0a41e80 1731 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1732 &dummy1, &dummy2, idx6o,
1733 &dummy3, &dummy4, &dummy5);
1734// printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n",
1735// idx21s, idx23as, idx22s, idx23bs, *idx5o, *idx6o);
0d64b05f 1736}