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