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