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