1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////////
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. //
23 ///////////////////////////////////////////////////////////////////////////////
25 #include <fstream> // needed for raw data dump
34 #include <TClonesArray.h>
38 #include "AliRunLoader.h"
39 #include "AliLoader.h"
40 #include "AliTRDdigit.h"
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 "AliTRDpadPlane.h"
50 #include "AliTRDtrackletMCM.h"
51 #include "AliTRDmcmSim.h"
53 ClassImp(AliTRDmcmSim)
55 //_____________________________________________________________________________
56 AliTRDmcmSim::AliTRDmcmSim() : TObject()
85 // AliTRDmcmSim default constructor
86 // By default, nothing is initialized.
87 // It is necessary to issue Init before use.
90 AliTRDmcmSim::~AliTRDmcmSim()
93 // AliTRDmcmSim destructor
97 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
98 delete [] fADCR[iadc];
99 delete [] fADCF[iadc];
100 delete [] fZSM [iadc];
109 delete [] fGainCounterA;
110 delete [] fGainCounterB;
111 delete [] fTailAmplLong;
112 delete [] fTailAmplShort;
115 fTrackletArray->Delete();
116 delete fTrackletArray;
121 void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ )
124 // Initialize the class with new geometry information
125 // fADC array will be reused with filled by zero
129 fFeeParam = AliTRDfeeParam::Instance();
130 fTrapConfig = AliTRDtrapConfig::Instance();
131 fSimParam = AliTRDSimParam::Instance();
132 fCal = AliTRDcalibDB::Instance();
133 fGeo = new AliTRDgeometry();
139 fNADC = fFeeParam->GetNadcMcm();
140 fNTimeBin = fCal->GetNumberOfTimeBins();
141 fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
142 fMaxTracklets = fFeeParam->GetMaxNrOfTracklets();
145 fADCR = new Int_t *[fNADC];
146 fADCF = new Int_t *[fNADC];
147 fZSM = new Int_t *[fNADC];
148 fZSM1Dim = new Int_t [fNADC];
149 fGainCounterA = new UInt_t[fNADC];
150 fGainCounterB = new UInt_t[fNADC];
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];
158 fPedAcc = new UInt_t[fNADC]; // accumulator for pedestal filter
159 fTailAmplLong = new UShort_t[fNADC];
160 fTailAmplShort = new UShort_t[fNADC];
162 // tracklet calculation
163 fFitReg = new FitReg_t[fNADC];
164 fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fMaxTracklets);
166 fMCMT = new UInt_t[fMaxTracklets];
169 fInitialized = kTRUE;
174 void AliTRDmcmSim::Reset()
176 // Resets the data values and internal filter registers
177 // by re-initialising them
179 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
180 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
183 fZSM [iadc][it] = 1; // Default unread = 1
185 fZSM1Dim[iadc] = 1; // Default unread = 1
186 fGainCounterA[iadc] = 0;
187 fGainCounterB[iadc] = 0;
190 for(Int_t i = 0; i < fMaxTracklets; i++) {
194 FilterPedestalInit();
196 FilterTailInit(fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP)); //??? not really correct if gain filter is active
199 Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm)
201 // loads the ADC data as obtained from the digitsManager for the specified MCM
206 AliError("No Runloader given");
210 AliLoader *trdLoader = runloader->GetLoader("TRDLoader");
212 AliError("Could not get TRDLoader");
216 trdLoader->LoadDigits();
217 AliTRDdigitsManager *digMgr = new AliTRDdigitsManager();
218 digMgr->SetSDigits(0);
219 digMgr->CreateArrays();
220 digMgr->ReadDigits(trdLoader->TreeD());
221 AliTRDarrayADC *digits = (AliTRDarrayADC*) digMgr->GetDigits(det);
222 if (!digits->HasData())
226 Int_t padrow = fFeeParam->GetPadRowFromMCM(rob, mcm);
228 for (Int_t ch = 0; ch < fNADC; ch++) {
230 for (Int_t tb = 0; tb < fNTimeBin; tb++) {
236 if (digits->GetData(padrow,padcol, tb) < 0) {
241 fADCR[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits;
242 fADCF[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits;
247 digMgr->RemoveDigits(det);
253 void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inputGain, Int_t inputTail)
255 // This function can be used to test the filters.
256 // It feeds nsamples of ADC values with a gaussian distribution specified by mean and sigma.
257 // The filter chain implemented here consists of:
258 // Pedestal -> Gain -> Tail
259 // With inputGain and inputTail the input to the gain and tail filter, respectively,
260 // can be chosen where
262 // 1: pedestal output
264 // The input has to be chosen from a stage before.
265 // The filter behaviour is controlled by the TRAP parameters from AliTRDtrapConfig in the
266 // same way as in normal simulation.
267 // The functions produces four histograms with the values at the different stages.
269 TH1F *h = new TH1F("noise", "Gaussian Noise;sample;ADC count",
270 nsamples, 0, nsamples);
271 TH1F *hfp = new TH1F("pedf", "Noise #rightarrow Pedestal filter;sample;ADC count", nsamples, 0, nsamples);
272 TH1F *hfg = new TH1F("pedg", "Pedestal #rightarrow Gain;sample;ADC count", nsamples, 0, nsamples);
273 TH1F *hft = new TH1F("pedt", "Gain #rightarrow Tail;sample;ADC count", nsamples, 0, nsamples);
275 hfp->SetStats(kFALSE);
276 hfg->SetStats(kFALSE);
277 hft->SetStats(kFALSE);
279 Int_t value; // ADC count with noise (10 bit)
280 Int_t valuep; // pedestal filter output (12 bit)
281 Int_t valueg; // gain filter output (12 bit)
282 Int_t valuet; // tail filter value (12 bit)
284 for (Int_t i = 0; i < nsamples; i++) {
285 value = (Int_t) gRandom->Gaus(mean, sigma); // generate noise with gaussian distribution
286 h->SetBinContent(i, value);
288 valuep = FilterPedestalNextSample(1, 0, ((Int_t) value) << 2);
291 valueg = FilterGainNextSample(1, ((Int_t) value) << 2);
293 valueg = FilterGainNextSample(1, valuep);
296 valuet = FilterTailNextSample(1, ((Int_t) value) << 2);
297 else if (inputTail == 1)
298 valuet = FilterTailNextSample(1, valuep);
300 valuet = FilterTailNextSample(1, valueg);
302 hfp->SetBinContent(i, valuep >> 2);
303 hfg->SetBinContent(i, valueg >> 2);
304 hft->SetBinContent(i, valuet >> 2);
307 TCanvas *c = new TCanvas;
319 Bool_t AliTRDmcmSim::CheckInitialized()
322 // Check whether object is initialized
325 if( ! fInitialized ) {
326 AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
331 void AliTRDmcmSim::Print(Option_t* const option) const
333 // Prints the data stored and/or calculated for this MCM.
334 // The output is controlled by option which can be a sequence of any of
335 // the following characters:
336 // R - prints raw ADC data
337 // F - prints filtered data
338 // H - prints detected hits
339 // T - prints found tracklets
340 // The later stages are only useful when the corresponding calculations
341 // have been performed.
343 printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector);
345 TString opt = option;
346 if (opt.Contains("U")) {
347 printf("Raw ADC data (10 bit):\n");
348 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
349 for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) {
350 printf("%5i", fADCR[iChannel][iTimeBin] >> fgkAddDigits);
356 if (opt.Contains("F")) {
357 printf("Filtered data (12 bit):\n");
358 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
359 for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) {
360 printf("%5i", fADCF[iChannel][iTimeBin]);
366 if (opt.Contains("H")) {
367 printf("Found %i hits:\n", fNHits);
368 for (Int_t iHit = 0; iHit < fNHits; iHit++) {
369 printf("Hit %3i in timebin %2i, ADC %2i has charge %3i and position %3i\n",
370 iHit, fHits[iHit].fTimebin, fHits[iHit].fChannel, fHits[iHit].fQtot, fHits[iHit].fYpos);
374 if (opt.Contains("T")) {
375 printf("Tracklets:\n");
376 for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntriesFast(); iTrkl++) {
377 printf("tracklet %i: 0x%08x\n", iTrkl, ((AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl])->GetTrackletWord());
382 void AliTRDmcmSim::Draw(Option_t* const option)
384 // Plots the data stored in a 2-dim. timebin vs. ADC channel plot.
385 // The option selects what data is plotted and can be a sequence of
386 // the following characters:
387 // R - plot raw data (default)
388 // F - plot filtered data (meaningless if R is specified)
389 // In addition to the ADC values:
391 // T - plot tracklets
393 TString opt = option;
395 TH2F *hist = new TH2F("mcmdata", Form("Data of MCM %i on ROB %i in detector %i", \
396 fMcmPos, fRobPos, fDetector), \
397 fNADC, -0.5, fNADC-.5, fNTimeBin, -.5, fNTimeBin-.5);
398 hist->GetXaxis()->SetTitle("ADC Channel");
399 hist->GetYaxis()->SetTitle("Timebin");
400 hist->SetStats(kFALSE);
402 if (opt.Contains("R")) {
403 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
404 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
405 hist->SetBinContent(iAdc+1, iTimeBin+1, fADCR[iAdc][iTimeBin] >> fgkAddDigits);
410 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
411 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
412 hist->SetBinContent(iAdc+1, iTimeBin+1, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
418 if (opt.Contains("H")) {
419 TGraph *grHits = new TGraph();
420 for (Int_t iHit = 0; iHit < fNHits; iHit++) {
421 grHits->SetPoint(iHit,
422 fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256.,
423 fHits[iHit].fTimebin);
428 if (opt.Contains("T")) {
429 TLine *trklLines = new TLine[4];
430 for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) {
431 AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
432 AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl];
433 Float_t offset = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19)) + 19 * pp->GetWidthIPad();
434 trklLines[iTrkl].SetX1((offset - trkl->GetY()) / pp->GetWidthIPad());
435 trklLines[iTrkl].SetY1(0);
436 trklLines[iTrkl].SetX2((offset - (trkl->GetY() + ((Float_t) trkl->GetdY())*140e-4)) / pp->GetWidthIPad());
437 trklLines[iTrkl].SetY2(fNTimeBin - 1);
438 trklLines[iTrkl].SetLineColor(2);
439 trklLines[iTrkl].SetLineWidth(2);
440 printf("Tracklet %i: y = %f, dy = %f, offset = %f\n", iTrkl, trkl->GetY(), (trkl->GetdY() * 140e-4), offset);
441 trklLines[iTrkl].Draw();
446 void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc )
449 // Store ADC data into array of raw data
452 if( !CheckInitialized() ) return;
454 if( iadc < 0 || iadc >= fNADC ) {
455 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
459 for( int it = 0 ; it < fNTimeBin ; it++ ) {
460 fADCR[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
461 fADCF[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
465 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
468 // Store ADC data into array of raw data
471 if( !CheckInitialized() ) return;
473 if( iadc < 0 || iadc >= fNADC ) {
474 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
478 fADCR[iadc][it] = adc << fgkAddDigits;
479 fADCF[iadc][it] = adc << fgkAddDigits;
482 void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray)
484 // Set the ADC data from an AliTRDarrayADC
487 AliError("Called uninitialized! Nothing done!");
492 Int_t lastAdc = fNADC-1;
494 while (GetCol(firstAdc) < 0) {
495 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
496 fADCR[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
497 fADCF[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
502 while (GetCol(lastAdc) < 0) {
503 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
504 fADCR[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
505 fADCF[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits;
510 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
511 for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
512 Int_t value = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin);
514 fADCR[iAdc][iTimeBin] = 0;
515 fADCF[iAdc][iTimeBin] = 0;
518 fADCR[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits;
519 fADCF[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits;
525 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
528 // Store ADC data into array of raw data
531 if( !CheckInitialized() ) return;
533 if( iadc < 0 || iadc >= fNADC ) {
534 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
538 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
539 fADCR[iadc][it] = fSimParam->GetADCbaseline() << fgkAddDigits;
540 fADCF[iadc][it] = fSimParam->GetADCbaseline() << fgkAddDigits;
544 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
547 // Return column id of the pad for the given ADC channel
550 if( !CheckInitialized() )
553 Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
554 if (col < 0 || col >= fFeeParam->GetNcol())
560 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
563 // Produce raw data stream from this MCM and put in buf
564 // Returns number of words filled, or negative value
565 // with -1 * number of overflowed words
569 Int_t nw = 0; // Number of written words
570 Int_t of = 0; // Number of overflowed words
571 Int_t rawVer = fFeeParam->GetRAWversion();
573 Int_t nActiveADC = 0; // number of activated ADC bits in a word
575 if( !CheckInitialized() ) return 0;
577 if( fFeeParam->GetRAWstoreRaw() ) {
583 // Produce MCM header
584 x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
588 //printf("\nMCM header: %X ",x);
594 // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
595 // n : unused , c : ADC count, m : selected ADCs
598 for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
599 if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed
600 x = x | (1 << (iAdc+4) ); // last 4 digit reserved for 1100=0xc
601 nActiveADC++; // number of 1 in mmm....m
604 x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC; // nn = 01, ccccc are inverted, 0xc=1100
605 //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
609 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
616 // Produce ADC data. 3 timebins are packed into one 32 bits word
617 // In this version, different ADC channel will NOT share the same word
619 UInt_t aa=0, a1=0, a2=0, a3=0;
621 for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
622 if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
623 aa = !(iAdc & 1) + 2;
624 for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
625 a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] >> fgkAddDigits : 0;
626 a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] >> fgkAddDigits : 0;
627 a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] >> fgkAddDigits : 0;
628 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
639 if( of != 0 ) return -of; else return nw;
642 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
645 // Produce tracklet data stream from this MCM and put in buf
646 // Returns number of words filled, or negative value
647 // with -1 * number of overflowed words
650 Int_t nw = 0; // Number of written words
651 Int_t of = 0; // Number of overflowed words
653 if( !CheckInitialized() ) return 0;
655 // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM
656 // fMCMT is filled continuously until no more tracklet words available
658 for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
660 buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord();
665 if( of != 0 ) return -of; else return nw;
668 void AliTRDmcmSim::Filter()
671 // Filter the raw ADC values. The active filter stages and their
672 // parameters are taken from AliTRDtrapConfig.
673 // The raw data is stored separate from the filtered data. Thus,
674 // it is possible to run the filters on a set of raw values
675 // sequentially for parameter tuning.
678 if( !CheckInitialized() ) {
679 AliError("got called before initialization! Nothing done!");
683 // Apply filters sequentially. Bypass is handled by filters
684 // since counters and internal registers may be updated even
685 // if the filter is bypassed.
686 // The first filter takes the data from fADCR and
689 // Non-linearity filter not implemented.
693 // Crosstalk filter not implemented.
696 void AliTRDmcmSim::FilterPedestalInit()
698 // Initializes the pedestal filter assuming that the input has
699 // been constant for a long time (compared to the time constant).
701 // UShort_t fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
702 UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
703 UShort_t shifts[4] = {11, 14, 17, 21}; //??? where to take shifts from?
705 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
706 fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
709 UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
711 // Returns the output of the pedestal filter given the input value.
712 // The output depends on the internal registers and, thus, the
713 // history of the filter.
715 UShort_t fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
716 UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
717 UShort_t fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY); // 0..1 the bypass, active low
718 UShort_t shifts[4] = {11, 14, 17, 21}; //??? where to come from
720 UShort_t accumulatorShifted;
724 inpAdd = value + fpnp;
726 if (fpby == 0) //??? before or after update of accumulator
729 accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF; // 10 bits
730 if (timebin == 0) // the accumulator is disabled in the drift time
732 correction = (value & 0x3FF) - accumulatorShifted;
733 fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF; // 31 bits
736 if (inpAdd <= accumulatorShifted)
740 inpAdd = inpAdd - accumulatorShifted;
748 void AliTRDmcmSim::FilterPedestal()
751 // Apply pedestal filter
753 // As the first filter in the chain it reads data from fADCR
754 // and outputs to fADCF.
755 // It has only an effect if previous samples have been fed to
756 // find the pedestal. Currently, the simulation assumes that
757 // the input has been stable for a sufficiently long time.
759 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
760 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
761 fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]);
766 void AliTRDmcmSim::FilterGainInit()
768 // Initializes the gain filter. In this case, only threshold
769 // counters are reset.
771 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
772 // these are counters which in hardware continue
773 // until maximum or reset
774 fGainCounterA[iAdc] = 0;
775 fGainCounterB[iAdc] = 0;
779 UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
781 // Apply the gain filter to the given value.
782 // BEGIN_LATEX O_{i}(t) = #gamma_{i} * I_{i}(t) + a_{i} END_LATEX
783 // The output depends on the internal registers and, thus, the
784 // history of the filter.
786 UShort_t fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY); // bypass, active low
787 UShort_t fgf = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc)); // 0x700 + (0 & 0x1ff);
788 UShort_t fga = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc)); // 40;
789 UShort_t fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA); // 20;
790 UShort_t fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB); // 2060;
795 tmp = (value * fgf) >> 11;
796 if (tmp > 0xFFF) tmp = 0xFFF;
799 value = AddUintClipping(tmp, fga, 12);
801 // Update threshold counters
802 // not really useful as they are cleared with every new event
803 if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))
806 fGainCounterB[adc]++;
807 else if (value >= fgta)
808 fGainCounterA[adc]++;
814 void AliTRDmcmSim::FilterGain()
816 // Read data from fADCF and apply gain filter.
818 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
819 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
820 fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
825 void AliTRDmcmSim::FilterTailInit(Int_t baseline)
827 // Initializes the tail filter assuming that the input has
828 // been at the baseline value (configured by FTFP) for a
829 // sufficiently long time.
831 // exponents and weight calculated from configuration
832 UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
833 UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
834 UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
836 Float_t lambdaL = lambdaLong * 1.0 / (1 << 11);
837 Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
838 Float_t alphaL = alphaLong * 1.0 / (1 << 11);
840 qup = (1 - lambdaL) * (1 - lambdaS);
841 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
842 Float_t kdc = qup/qdn;
848 aout = baseline - (UShort_t) kt;
849 ql = lambdaL * (1 - lambdaS) * alphaL;
850 qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
852 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
853 fTailAmplLong[iAdc] = (UShort_t) (aout * ql / (ql + qs));
854 fTailAmplShort[iAdc] = (UShort_t) (aout * qs / (ql + qs));
858 UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
860 // Returns the output of the tail filter for the given input value.
861 // The output depends on the internal registers and, thus, the
862 // history of the filter.
864 // exponents and weight calculated from configuration
865 UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
866 UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
867 UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
869 Float_t lambdaL = lambdaLong * 1.0 / (1 << 11);
870 Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
871 Float_t alphaL = alphaLong * 1.0 / (1 << 11);
873 qup = (1 - lambdaL) * (1 - lambdaS);
874 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
875 // Float_t kdc = qup/qdn;
882 UShort_t inpVolt = value & 0xFFF; // 12 bits
884 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
888 // add the present generator outputs
889 aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
891 // calculate the difference between the input the generated signal
893 aDiff = inpVolt - aQ;
897 // the inputs to the two generators, weighted
898 alInpv = (aDiff * alphaLong) >> 11;
900 // the new values of the registers, used next time
902 tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
903 tmp = (tmp * lambdaLong) >> 11;
904 fTailAmplLong[adc] = tmp & 0xFFF;
906 tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
907 tmp = (tmp * lambdaShort) >> 11;
908 fTailAmplShort[adc] = tmp & 0xFFF;
910 // the output of the filter
915 void AliTRDmcmSim::FilterTail()
919 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
920 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
921 fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]);
926 void AliTRDmcmSim::ZSMapping()
929 // Zero Suppression Mapping implemented in TRAP chip
931 // See detail TRAP manual "Data Indication" section:
932 // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
935 //??? values should come from TRAPconfig
936 Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS); // TRAP default = 0x4 (Tis=4)
937 Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT); // TRAP default = 0x28 (Tit=40)
938 Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL); // TRAP default = 0xf0
939 // (lookup table accept (I2,I1,I0)=(111)
940 // or (110) or (101) or (100))
941 Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); // TRAP default = 1 (no neighbor sensitivity)
942 Int_t ep = 0; // fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); //??? really subtracted here
946 if( !CheckInitialized() ) {
947 AliError("got called uninitialized! Nothing done!");
951 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
952 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
954 // Get ADC data currently in filter buffer
955 Int_t ap = adc[iadc-1][it] - ep; // previous
956 Int_t ac = adc[iadc ][it] - ep; // current
957 Int_t an = adc[iadc+1][it] - ep; // next
959 // evaluate three conditions
960 Int_t i0 = ( ac >= ap && ac >= an ) ? 0 : 1; // peak center detection
961 Int_t i1 = ( ap + ac + an > eBIT ) ? 0 : 1; // cluster
962 Int_t i2 = ( ac > eBIS ) ? 0 : 1; // absolute large peak
964 Int_t i = i2 * 4 + i1 * 2 + i0; // Bit position in lookup table
965 Int_t d = (eBIL >> i) & 1; // Looking up (here d=0 means true
966 // and d=1 means false according to TRAP manual)
969 if( eBIN == 0 ) { // turn on neighboring ADCs
970 fZSM[iadc-1][it] &= d;
971 fZSM[iadc+1][it] &= d;
976 // do 1 dim projection
977 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
978 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
979 fZSM1Dim[iadc] &= fZSM[iadc][it];
984 void AliTRDmcmSim::DumpData( char *f, char *target )
987 // Dump data stored (for debugging).
988 // target should contain one or multiple of the following characters
990 // F for filtered data
991 // Z for zero suppression map
993 // other characters are simply ignored
996 UInt_t tempbuf[1024];
998 if( !CheckInitialized() ) return;
1000 std::ofstream of( f, std::ios::out | std::ios::app );
1001 of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
1002 fDetector, fGeo->GetSector(fDetector), fGeo->GetStack(fDetector),
1003 fGeo->GetSector(fDetector), fRobPos, fMcmPos );
1005 for( int t=0 ; target[t] != 0 ; t++ ) {
1006 switch( target[t] ) {
1009 of << Form("fADCR (raw ADC data)\n");
1010 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1011 of << Form(" ADC %02d: ", iadc);
1012 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1013 of << Form("% 4d", fADCR[iadc][it]);
1020 of << Form("fADCF (filtered ADC data)\n");
1021 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1022 of << Form(" ADC %02d: ", iadc);
1023 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1024 of << Form("% 4d", fADCF[iadc][it]);
1031 of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
1032 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1033 of << Form(" ADC %02d: ", iadc);
1034 if( fZSM1Dim[iadc] == 0 ) { of << " R " ; } else { of << " . "; } // R:read .:suppressed
1035 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1036 if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
1043 Int_t s = ProduceRawStream( tempbuf, 1024 );
1044 of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
1045 of << Form(" address data\n");
1046 for( int i = 0 ; i < s ; i++ ) {
1047 of << Form(" %04x %08x\n", i, tempbuf[i]);
1053 void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label)
1055 // Add the given hit to the fit register which is lateron used for
1056 // the tracklet calculation.
1057 // In addition to the fit sums in the fit register MC information
1060 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) &&
1061 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
1062 fFitReg[adc].fQ0 += qtot;
1064 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) &&
1065 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
1066 fFitReg[adc].fQ1 += qtot;
1068 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) &&
1069 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
1071 fFitReg[adc].fSumX += timebin;
1072 fFitReg[adc].fSumX2 += timebin*timebin;
1073 fFitReg[adc].fNhits++;
1074 fFitReg[adc].fSumY += ypos;
1075 fFitReg[adc].fSumY2 += ypos*ypos;
1076 fFitReg[adc].fSumXY += timebin*ypos;
1079 // register hits (MC info)
1080 fHits[fNHits].fChannel = adc;
1081 fHits[fNHits].fQtot = qtot;
1082 fHits[fNHits].fYpos = ypos;
1083 fHits[fNHits].fTimebin = timebin;
1084 fHits[fNHits].fLabel = label;
1088 void AliTRDmcmSim::CalcFitreg()
1091 // Detect the hits and fill the fit registers.
1092 // Requires 12-bit data from fADCF which means Filter()
1093 // has to be called before even if all filters are bypassed.
1097 const uint16_t lutPos[128] = { // move later to some other file
1098 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15,
1099 16, 16, 16, 17, 17, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26,
1100 27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 27, 27, 27, 27, 26,
1101 26, 26, 26, 25, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 20, 20, 19, 18, 18, 17, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 7};
1103 //??? to be clarified:
1104 UInt_t adcMask = 0xffffffff;
1106 UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
1107 Short_t ypos, fromLeft, fromRight, found;
1108 UShort_t qTotal[19]; // the last is dummy
1109 UShort_t marked[6], qMarked[6], worse1, worse2;
1111 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS);
1112 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)
1114 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1115 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE);
1116 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)
1118 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1120 // reset the fit registers
1122 for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels
1124 fFitReg[adcch].fNhits = 0;
1125 fFitReg[adcch].fQ0 = 0;
1126 fFitReg[adcch].fQ1 = 0;
1127 fFitReg[adcch].fSumX = 0;
1128 fFitReg[adcch].fSumY = 0;
1129 fFitReg[adcch].fSumX2 = 0;
1130 fFitReg[adcch].fSumY2 = 0;
1131 fFitReg[adcch].fSumXY = 0;
1134 for (timebin = timebin1; timebin < timebin2; timebin++)
1136 // first find the hit candidates and store the total cluster charge in qTotal array
1137 // in case of not hit store 0 there.
1138 for (adcch = 0; adcch < fNADC-2; adcch++) {
1139 if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
1141 adcLeft = fADCF[adcch ][timebin];
1142 adcCentral = fADCF[adcch+1][timebin];
1143 adcRight = fADCF[adcch+2][timebin];
1144 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1)
1145 hitQual = ( (adcLeft * adcRight) <
1146 (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) );
1149 // The accumulated charge is with the pedestal!!!
1150 qtotTemp = adcLeft + adcCentral + adcRight;
1152 (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1153 (adcLeft <= adcCentral) &&
1154 (adcCentral > adcRight) )
1155 qTotal[adcch] = qtotTemp;
1158 //printf("ch %2d qTotal %5d\n",adcch, qTotal[adcch]);
1161 qTotal[adcch] = 0; //jkl
1167 marked[4] = 19; // invalid channel
1168 marked[5] = 19; // invalid channel
1170 while ((adcch < 16) && (found < 3))
1172 if (qTotal[adcch] > 0)
1175 marked[2*found+1]=adcch;
1184 while ((adcch > 2) && (found < 3))
1186 if (qTotal[adcch] > 0)
1188 marked[2*found]=adcch;
1195 //printf("Fromleft=%d, Fromright=%d\n",fromLeft, fromRight);
1196 // here mask the hit candidates in the middle, if any
1197 if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight))
1198 for (adcch = fromLeft+1; adcch < fromRight; adcch++)
1202 for (adcch = 0; adcch < 19; adcch++)
1203 if (qTotal[adcch] > 0) found++;
1206 if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1208 if (marked[4] == marked[5]) marked[5] = 19;
1209 for (found=0; found<6; found++)
1211 qMarked[found] = qTotal[marked[found]] >> 4;
1212 //printf("ch_%d qTotal %d qTotals %d |",marked[found],qTotal[marked[found]],qMarked[found]);
1216 Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
1224 // Now mask the two channels with the smallest charge
1228 //printf("Kill ch %d\n",worse1);
1233 //printf("Kill ch %d\n",worse2);
1237 for (adcch = 0; adcch < 19; adcch++) {
1238 if (qTotal[adcch] > 0) // the channel is marked for processing
1240 adcLeft = fADCF[adcch ][timebin];
1241 adcCentral = fADCF[adcch+1][timebin];
1242 adcRight = fADCF[adcch+2][timebin];
1243 // hit detected, in TRAP we have 4 units and a hit-selection, here we proceed all channels!
1244 // subtract the pedestal TPFP, clipping instead of wrapping
1246 Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP);
1247 // printf("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n",
1248 // timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP,
1249 // fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT));
1251 if (adcLeft < regTPFP) adcLeft = 0; else adcLeft -= regTPFP;
1252 if (adcCentral < regTPFP) adcCentral = 0; else adcCentral -= regTPFP;
1253 if (adcRight < regTPFP) adcRight = 0; else adcRight -= regTPFP;
1255 // Calculate the center of gravity
1256 // checking for adcCentral != 0 (in case of "bad" configuration)
1257 if (adcCentral == 0)
1259 ypos = 128*(adcLeft - adcRight) / adcCentral;
1260 if (ypos < 0) ypos = -ypos;
1261 // make the correction using the LUT
1262 ypos = ypos + lutPos[ypos & 0x7F];
1263 if (adcLeft > adcRight) ypos = -ypos;
1264 AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, -1);
1270 void AliTRDmcmSim::TrackletSelection()
1272 // Select up to 4 tracklet candidates from the fit registers
1273 // and assign them to the CPUs.
1275 UShort_t adcIdx, i, j, ntracks, tmp;
1276 UShort_t trackletCand[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates
1279 for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs
1280 if ( (fFitReg[adcIdx].fNhits
1281 >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) &&
1282 (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits
1283 >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT)))
1285 trackletCand[ntracks][0] = adcIdx;
1286 trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits;
1287 //printf("%d %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]);
1291 // for (i=0; i<ntracks;i++) printf("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]);
1295 // primitive sorting according to the number of hits
1296 for (j = 0; j < (ntracks-1); j++)
1298 for (i = j+1; i < ntracks; i++)
1300 if ( (trackletCand[j][1] < trackletCand[i][1]) ||
1301 ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
1304 tmp = trackletCand[j][1];
1305 trackletCand[j][1] = trackletCand[i][1];
1306 trackletCand[i][1] = tmp;
1307 tmp = trackletCand[j][0];
1308 trackletCand[j][0] = trackletCand[i][0];
1309 trackletCand[i][0] = tmp;
1313 ntracks = 4; // cut the rest, 4 is the max
1315 // else is not necessary to sort
1317 // now sort, so that the first tracklet going to CPU0 corresponds to the highest adc channel - as in the TRAP
1318 for (j = 0; j < (ntracks-1); j++)
1320 for (i = j+1; i < ntracks; i++)
1322 if (trackletCand[j][0] < trackletCand[i][0])
1325 tmp = trackletCand[j][1];
1326 trackletCand[j][1] = trackletCand[i][1];
1327 trackletCand[i][1] = tmp;
1328 tmp = trackletCand[j][0];
1329 trackletCand[j][0] = trackletCand[i][0];
1330 trackletCand[i][0] = tmp;
1334 for (i = 0; i < ntracks; i++) // CPUs with tracklets.
1335 fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i]
1336 for (i = ntracks; i < 4; i++) // CPUs without tracklets
1337 fFitPtr[i] = 31; // pointer to the left channel with tracklet for CPU[i] = 31 (invalid)
1338 // printf("found %i tracklet candidates\n", ntracks);
1341 void AliTRDmcmSim::FitTracklet()
1343 // Perform the actual tracklet fit based on the fit sums
1344 // which have been filled in the fit registers.
1346 // parameters in fitred.asm (fit program)
1347 Int_t decPlaces = 5;
1350 rndAdd = (1 << (decPlaces-1)) + 1;
1351 else if (decPlaces == 1)
1354 // should come from trapConfig (DMEM)
1355 AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
1356 Long64_t shift = ((Long64_t) 1 << 32);
1357 UInt_t scaleY = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4)));
1358 UInt_t scaleD = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4)));
1359 int padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos);
1360 int yoffs = (fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19) - fFeeParam->GetNcol()/2) << (8 + decPlaces);
1361 int ndrift = 20; //??? value in simulation?
1362 int deflCorr = 0; // -370;
1363 int minslope = -10000; // no pt-cut so far
1364 int maxslope = 10000; // no pt-cut so far
1366 // local variables for calculation
1367 Long64_t mult, temp, denom; //???
1368 UInt_t q0, q1, qTotal; // charges in the two windows and total charge
1369 UShort_t nHits; // number of hits
1370 Int_t slope, offset; // slope and offset of the tracklet
1371 Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers
1372 //int32_t SumY2; // not used in the current TRAP program
1373 FitReg_t *fit0, *fit1; // pointers to relevant fit registers
1375 // const uint32_t OneDivN[32] = { // 2**31/N : exactly like in the TRAP, the simple division here gives the same result!
1376 // 0x00000000, 0x80000000, 0x40000000, 0x2AAAAAA0, 0x20000000, 0x19999990, 0x15555550, 0x12492490,
1377 // 0x10000000, 0x0E38E380, 0x0CCCCCC0, 0x0BA2E8B0, 0x0AAAAAA0, 0x09D89D80, 0x09249240, 0x08888880,
1378 // 0x08000000, 0x07878780, 0x071C71C0, 0x06BCA1A0, 0x06666660, 0x06186180, 0x05D17450, 0x0590B210,
1379 // 0x05555550, 0x051EB850, 0x04EC4EC0, 0x04BDA120, 0x04924920, 0x0469EE50, 0x04444440, 0x04210840};
1381 for (Int_t cpu = 0; cpu < 4; cpu++) {
1382 if (fFitPtr[cpu] == 31)
1384 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1388 fit0 = &fFitReg[fFitPtr[cpu] ];
1389 fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1392 mult = mult << (32 + decPlaces);
1396 nHits = fit0->fNhits + fit1->fNhits; // number of hits
1397 sumX = fit0->fSumX + fit1->fSumX;
1398 sumX2 = fit0->fSumX2 + fit1->fSumX2;
1399 denom = nHits*sumX2 - sumX*sumX;
1401 mult = mult / denom; // exactly like in the TRAP program
1402 q0 = fit0->fQ0 + fit1->fQ0;
1403 q1 = fit0->fQ1 + fit1->fQ1;
1404 sumY = fit0->fSumY + fit1->fSumY + 256*fit1->fNhits;
1405 sumXY = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
1407 slope = nHits*sumXY - sumX * sumY;
1408 offset = sumX2*sumY - sumX * sumXY;
1409 temp = mult * slope;
1410 slope = temp >> 32; // take the upper 32 bits
1411 temp = mult * offset;
1412 offset = temp >> 32; // take the upper 32 bits
1414 offset = offset + yoffs + (18 << (8 + decPlaces));
1415 slope = slope * ndrift + deflCorr;
1416 offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
1418 if ((slope < minslope) || (slope > maxslope))
1420 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1425 temp = temp * scaleD;
1426 slope = (temp >> 32);
1429 temp = temp * scaleY;
1430 offset = (temp >> 32);
1432 // rounding, like in the TRAP
1433 slope = (slope + rndAdd) >> decPlaces;
1434 offset = (offset + rndAdd) >> decPlaces;
1436 if (slope > 0x3f || slope < -0x3f)
1437 AliWarning("Overflow in slope");
1438 slope = slope & 0x7F; // 7 bit
1440 if (offset > 0xfff || offset < -0xfff)
1441 AliWarning("Overflow in offset");
1442 offset = offset & 0x1FFF; // 13 bit
1444 qTotal = (q1 / nHits) >> 1;
1446 AliWarning("Overflow in charge");
1447 qTotal = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1449 // assemble and store the tracklet word
1450 fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1451 new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1457 void AliTRDmcmSim::Tracklet()
1459 // Run the tracklet calculation by calling sequentially:
1460 // CalcFitreg(); TrackletSelection(); FitTracklet()
1461 // and store the tracklets
1463 if (!fInitialized) {
1464 AliError("Called uninitialized! Nothing done!");
1468 fTrackletArray->Delete();
1471 TrackletSelection();
1474 AliRunLoader *rl = AliRunLoader::Instance();
1475 AliDataLoader *dl = 0x0;
1477 dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1479 AliError("Could not get the tracklets data loader!");
1482 TTree *trackletTree = dl->Tree();
1483 if (!trackletTree) {
1485 trackletTree = dl->Tree();
1488 AliTRDtrackletMCM *trkl = 0x0;
1489 TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1491 trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
1493 for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
1494 trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
1495 trkbranch->SetAddress(&trkl);
1496 // printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord());
1499 dl->WriteData("OVERWRITE");
1503 void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
1505 // write back the processed data configured by EBSF
1506 // EBSF = 1: unfiltered data; EBSF = 0: filtered data
1507 // zero-suppressed valued are written as -1 to digits
1509 if (!fInitialized) {
1510 AliError("Called uninitialized! Nothing done!");
1515 Int_t lastAdc = fNADC - 1;
1517 while (GetCol(firstAdc) < 0)
1520 while (GetCol(lastAdc) < 0)
1523 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1525 for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
1526 if (fZSM1Dim[iAdc] == 1) {
1527 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1528 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1);
1529 // printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin,
1530 // digits->GetData(GetRow(), GetCol(iAdc), iTimeBin));
1536 for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
1537 if (fZSM1Dim[iAdc] == 0) {
1538 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1539 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
1543 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1544 digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1);
1545 // printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin);
1552 // help functions, to be cleaned up
1554 UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
1557 // This function adds a and b (unsigned) and clips to
1558 // the specified number of bits.
1564 UInt_t maxv = (1 << nbits) - 1;;
1570 if ((sum < a) || (sum < b))
1576 void AliTRDmcmSim::Sort2(uint16_t idx1i, uint16_t idx2i, \
1577 uint16_t val1i, uint16_t val2i, \
1578 uint16_t *idx1o, uint16_t *idx2o, \
1579 uint16_t *val1o, uint16_t *val2o) const
1581 // sorting for tracklet selection
1599 void AliTRDmcmSim::Sort3(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, \
1600 uint16_t val1i, uint16_t val2i, uint16_t val3i, \
1601 uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, \
1602 uint16_t *val1o, uint16_t *val2o, uint16_t *val3o)
1604 // sorting for tracklet selection
1609 if (val1i > val2i) sel=4; else sel=0;
1610 if (val2i > val3i) sel=sel + 2;
1611 if (val3i > val1i) sel=sel + 1;
1612 //printf("input channels %d %d %d, charges %d %d %d sel=%d\n",idx1i, idx2i, idx3i, val1i, val2i, val3i, sel);
1615 case 6 : // 1 > 2 > 3 => 1 2 3
1616 case 0 : // 1 = 2 = 3 => 1 2 3 : in this case doesn't matter, but so is in hardware!
1625 case 4 : // 1 > 2, 2 <= 3, 3 <= 1 => 1 3 2
1634 case 2 : // 1 <= 2, 2 > 3, 3 <= 1 => 2 1 3
1643 case 3 : // 1 <= 2, 2 > 3, 3 > 1 => 2 3 1
1652 case 1 : // 1 <= 2, 2 <= 3, 3 > 1 => 3 2 1
1661 case 5 : // 1 > 2, 2 <= 3, 3 > 1 => 3 1 2
1670 default: // the rest should NEVER happen!
1671 printf("ERROR in Sort3!!!\n");
1674 // printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
1677 void AliTRDmcmSim::Sort6To4(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \
1678 uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \
1679 uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, uint16_t *idx4o, \
1680 uint16_t *val1o, uint16_t *val2o, uint16_t *val3o, uint16_t *val4o)
1682 // sorting for tracklet selection
1684 uint16_t idx21s, idx22s, idx23s, dummy;
1685 uint16_t val21s, val22s, val23s;
1686 uint16_t idx23as, idx23bs;
1687 uint16_t val23as, val23bs;
1689 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1690 idx1o, &idx21s, &idx23as,
1691 val1o, &val21s, &val23as);
1693 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1694 idx2o, &idx22s, &idx23bs,
1695 val2o, &val22s, &val23bs);
1697 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
1699 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1700 idx3o, idx4o, &dummy,
1701 val3o, val4o, &dummy);
1705 void AliTRDmcmSim::Sort6To2Worst(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \
1706 uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \
1707 uint16_t *idx5o, uint16_t *idx6o)
1709 // sorting for tracklet selection
1711 uint16_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
1712 uint16_t val21s, val22s, val23s;
1713 uint16_t idx23as, idx23bs;
1714 uint16_t val23as, val23bs;
1716 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1717 &dummy1, &idx21s, &idx23as,
1718 &dummy2, &val21s, &val23as);
1720 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1721 &dummy1, &idx22s, &idx23bs,
1722 &dummy2, &val22s, &val23bs);
1724 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
1726 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1727 &dummy1, &dummy2, idx6o,
1728 &dummy3, &dummy4, &dummy5);
1729 // printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n",
1730 // idx21s, idx23as, idx22s, idx23bs, *idx5o, *idx6o);