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