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 "AliTRDarrayDictionary.h"
50 #include "AliTRDpadPlane.h"
51 #include "AliTRDtrackletMCM.h"
52 #include "AliTRDmcmSim.h"
55 #include "TGeoGlobalMagField.h"
57 ClassImp(AliTRDmcmSim)
59 Bool_t AliTRDmcmSim::fgApplyCut = kTRUE;
61 //_____________________________________________________________________________
62 AliTRDmcmSim::AliTRDmcmSim() : TObject()
93 // AliTRDmcmSim default constructor
94 // By default, nothing is initialized.
95 // It is necessary to issue Init before use.
98 AliTRDmcmSim::~AliTRDmcmSim()
101 // AliTRDmcmSim destructor
105 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
106 delete [] fADCR[iadc];
107 delete [] fADCF[iadc];
108 delete [] fZSM [iadc];
117 delete [] fGainCounterA;
118 delete [] fGainCounterB;
119 delete [] fTailAmplLong;
120 delete [] fTailAmplShort;
123 fTrackletArray->Delete();
124 delete fTrackletArray;
129 void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ )
132 // Initialize the class with new geometry information
133 // fADC array will be reused with filled by zero
137 fFeeParam = AliTRDfeeParam::Instance();
138 fTrapConfig = AliTRDtrapConfig::Instance();
139 fSimParam = AliTRDSimParam::Instance();
140 fCommonParam = AliTRDCommonParam::Instance();
141 fCal = AliTRDcalibDB::Instance();
142 fGeo = new AliTRDgeometry();
148 fNADC = fFeeParam->GetNadcMcm();
149 fNTimeBin = fCal->GetNumberOfTimeBins();
150 fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
151 fMaxTracklets = fFeeParam->GetMaxNrOfTracklets();
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];
167 fPedAcc = new UInt_t[fNADC]; // accumulator for pedestal filter
168 fTailAmplLong = new UShort_t[fNADC];
169 fTailAmplShort = new UShort_t[fNADC];
171 // tracklet calculation
172 fFitReg = new FitReg_t[fNADC];
173 fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fMaxTracklets);
175 fMCMT = new UInt_t[fMaxTracklets];
178 fInitialized = kTRUE;
183 void AliTRDmcmSim::Reset()
185 // Resets the data values and internal filter registers
186 // by re-initialising them
188 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
189 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
192 fZSM [iadc][it] = 1; // Default unread = 1
194 fZSM1Dim[iadc] = 1; // Default unread = 1
195 fGainCounterA[iadc] = 0;
196 fGainCounterB[iadc] = 0;
199 for(Int_t i = 0; i < fMaxTracklets; i++) {
203 FilterPedestalInit();
205 FilterTailInit(fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP)); //??? not really correct if gain filter is active
208 Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm)
210 // loads the ADC data as obtained from the digitsManager for the specified MCM
215 AliError("No Runloader given");
219 AliLoader *trdLoader = runloader->GetLoader("TRDLoader");
221 AliError("Could not get TRDLoader");
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())
236 Int_t padrow = fFeeParam->GetPadRowFromMCM(rob, mcm);
238 for (Int_t ch = 0; ch < fNADC; ch++) {
240 for (Int_t tb = 0; tb < fNTimeBin; tb++) {
246 if (digits->GetData(padrow,padcol, tb) < 0) {
251 fADCR[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits;
252 fADCF[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits;
262 void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inputGain, Int_t inputTail)
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
271 // 1: pedestal 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.
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);
284 hfp->SetStats(kFALSE);
285 hfg->SetStats(kFALSE);
286 hft->SetStats(kFALSE);
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)
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);
297 valuep = FilterPedestalNextSample(1, 0, ((Int_t) value) << 2);
300 valueg = FilterGainNextSample(1, ((Int_t) value) << 2);
302 valueg = FilterGainNextSample(1, valuep);
305 valuet = FilterTailNextSample(1, ((Int_t) value) << 2);
306 else if (inputTail == 1)
307 valuet = FilterTailNextSample(1, valuep);
309 valuet = FilterTailNextSample(1, valueg);
311 hfp->SetBinContent(i, valuep >> 2);
312 hfg->SetBinContent(i, valueg >> 2);
313 hft->SetBinContent(i, valuet >> 2);
316 TCanvas *c = new TCanvas;
328 Bool_t AliTRDmcmSim::CheckInitialized()
331 // Check whether object is initialized
334 if( ! fInitialized ) {
335 AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
340 void AliTRDmcmSim::Print(Option_t* const option) const
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.
352 printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector);
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);
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]);
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);
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());
391 void AliTRDmcmSim::Draw(Option_t* const option)
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:
400 // T - plot tracklets
402 TString opt = option;
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);
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);
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);
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);
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();
455 void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc )
458 // Store ADC data into array of raw data
461 if( !CheckInitialized() ) return;
463 if( iadc < 0 || iadc >= fNADC ) {
464 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
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;
474 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
477 // Store ADC data into array of raw data
480 if( !CheckInitialized() ) return;
482 if( iadc < 0 || iadc >= fNADC ) {
483 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
487 fADCR[iadc][it] = adc << fgkAddDigits;
488 fADCF[iadc][it] = adc << fgkAddDigits;
491 void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager)
493 // Set the ADC data from an AliTRDarrayADC
496 AliError("Called uninitialized! Nothing done!");
500 fDigitsManager = digitsManager;
503 Int_t lastAdc = fNADC-1;
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;
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;
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);
525 fADCR[iAdc][iTimeBin] = 0;
526 fADCF[iAdc][iTimeBin] = 0;
529 fADCR[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits;
530 fADCF[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits;
536 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
539 // Store ADC data into array of raw data
542 if( !CheckInitialized() ) return;
544 if( iadc < 0 || iadc >= fNADC ) {
545 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
549 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
550 fADCR[iadc][it] = fSimParam->GetADCbaseline() << fgkAddDigits;
551 fADCF[iadc][it] = fSimParam->GetADCbaseline() << fgkAddDigits;
555 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
558 // Return column id of the pad for the given ADC channel
561 if( !CheckInitialized() )
564 Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
565 if (col < 0 || col >= fFeeParam->GetNcol())
571 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
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
580 Int_t nw = 0; // Number of written words
581 Int_t of = 0; // Number of overflowed words
582 Int_t rawVer = fFeeParam->GetRAWversion();
584 Int_t nActiveADC = 0; // number of activated ADC bits in a word
586 if( !CheckInitialized() ) return 0;
588 if( fFeeParam->GetRAWstoreRaw() ) {
594 // Produce MCM header
595 x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
599 //printf("\nMCM header: %X ",x);
605 // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
606 // n : unused , c : ADC count, m : selected ADCs
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
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 );
620 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
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
630 UInt_t aa=0, a1=0, a2=0, a3=0;
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;
650 if( of != 0 ) return -of; else return nw;
653 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
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
661 Int_t nw = 0; // Number of written words
662 Int_t of = 0; // Number of overflowed words
664 if( !CheckInitialized() ) return 0;
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
669 for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
671 buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord();
676 if( of != 0 ) return -of; else return nw;
679 void AliTRDmcmSim::Filter()
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.
689 if( !CheckInitialized() ) {
690 AliError("got called before initialization! Nothing done!");
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
700 // Non-linearity filter not implemented.
704 // Crosstalk filter not implemented.
707 void AliTRDmcmSim::FilterPedestalInit()
709 // Initializes the pedestal filter assuming that the input has
710 // been constant for a long time (compared to the time constant).
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?
716 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
717 fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
720 UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
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.
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
731 UShort_t accumulatorShifted;
735 inpAdd = value + fpnp;
737 if (fpby == 0) //??? before or after update of accumulator
740 accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF; // 10 bits
741 if (timebin == 0) // the accumulator is disabled in the drift time
743 correction = (value & 0x3FF) - accumulatorShifted;
744 fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF; // 31 bits
747 if (inpAdd <= accumulatorShifted)
751 inpAdd = inpAdd - accumulatorShifted;
759 void AliTRDmcmSim::FilterPedestal()
762 // Apply pedestal filter
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.
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]);
777 void AliTRDmcmSim::FilterGainInit()
779 // Initializes the gain filter. In this case, only threshold
780 // counters are reset.
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;
790 UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
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.
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;
806 tmp = (value * fgf) >> 11;
807 if (tmp > 0xFFF) tmp = 0xFFF;
810 value = AddUintClipping(tmp, fga, 12);
812 // Update threshold counters
813 // not really useful as they are cleared with every new event
814 if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))
817 fGainCounterB[adc]++;
818 else if (value >= fgta)
819 fGainCounterA[adc]++;
825 void AliTRDmcmSim::FilterGain()
827 // Read data from fADCF and apply gain filter.
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]);
836 void AliTRDmcmSim::FilterTailInit(Int_t baseline)
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.
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
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);
851 qup = (1 - lambdaL) * (1 - lambdaS);
852 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
853 Float_t kdc = qup/qdn;
859 aout = baseline - (UShort_t) kt;
860 ql = lambdaL * (1 - lambdaS) * alphaL;
861 qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
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));
869 UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
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.
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
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);
884 qup = (1 - lambdaL) * (1 - lambdaS);
885 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
886 // Float_t kdc = qup/qdn;
893 UShort_t inpVolt = value & 0xFFF; // 12 bits
895 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
899 // add the present generator outputs
900 aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
902 // calculate the difference between the input the generated signal
904 aDiff = inpVolt - aQ;
908 // the inputs to the two generators, weighted
909 alInpv = (aDiff * alphaLong) >> 11;
911 // the new values of the registers, used next time
913 tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
914 tmp = (tmp * lambdaLong) >> 11;
915 fTailAmplLong[adc] = tmp & 0xFFF;
917 tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
918 tmp = (tmp * lambdaShort) >> 11;
919 fTailAmplShort[adc] = tmp & 0xFFF;
921 // the output of the filter
926 void AliTRDmcmSim::FilterTail()
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]);
937 void AliTRDmcmSim::ZSMapping()
940 // Zero Suppression Mapping implemented in TRAP chip
942 // See detail TRAP manual "Data Indication" section:
943 // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
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
957 if( !CheckInitialized() ) {
958 AliError("got called uninitialized! Nothing done!");
962 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
963 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
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
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
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)
980 if( eBIN == 0 ) { // turn on neighboring ADCs
981 fZSM[iadc-1][it] &= d;
982 fZSM[iadc+1][it] &= d;
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];
995 void AliTRDmcmSim::DumpData( char *f, char *target )
998 // Dump data stored (for debugging).
999 // target should contain one or multiple of the following characters
1001 // F for filtered data
1002 // Z for zero suppression map
1003 // S Raw dat astream
1004 // other characters are simply ignored
1007 UInt_t tempbuf[1024];
1009 if( !CheckInitialized() ) return;
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 );
1016 for( int t=0 ; target[t] != 0 ; t++ ) {
1017 switch( target[t] ) {
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]);
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]);
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
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]);
1064 void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label)
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
1071 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) &&
1072 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
1073 fFitReg[adc].fQ0 += qtot;
1075 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) &&
1076 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
1077 fFitReg[adc].fQ1 += qtot;
1079 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) &&
1080 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
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;
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;
1099 void AliTRDmcmSim::CalcFitreg()
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.
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};
1114 //??? to be clarified:
1115 UInt_t adcMask = 0xffffffff;
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;
1122 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS);
1123 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)
1125 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1126 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE);
1127 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)
1129 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1131 // reset the fit registers
1133 for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels
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;
1145 for (timebin = timebin1; timebin < timebin2; timebin++)
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
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) );
1160 // The accumulated charge is with the pedestal!!!
1161 qtotTemp = adcLeft + adcCentral + adcRight;
1163 (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1164 (adcLeft <= adcCentral) &&
1165 (adcCentral > adcRight) )
1166 qTotal[adcch] = qtotTemp;
1169 //printf("ch %2d qTotal %5d\n",adcch, qTotal[adcch]);
1172 qTotal[adcch] = 0; //jkl
1178 marked[4] = 19; // invalid channel
1179 marked[5] = 19; // invalid channel
1181 while ((adcch < 16) && (found < 3))
1183 if (qTotal[adcch] > 0)
1186 marked[2*found+1]=adcch;
1195 while ((adcch > 2) && (found < 3))
1197 if (qTotal[adcch] > 0)
1199 marked[2*found]=adcch;
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++)
1213 for (adcch = 0; adcch < 19; adcch++)
1214 if (qTotal[adcch] > 0) found++;
1217 if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1219 if (marked[4] == marked[5]) marked[5] = 19;
1220 for (found=0; found<6; found++)
1222 qMarked[found] = qTotal[marked[found]] >> 4;
1223 //printf("ch_%d qTotal %d qTotals %d |",marked[found],qTotal[marked[found]],qMarked[found]);
1227 Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
1235 // Now mask the two channels with the smallest charge
1239 //printf("Kill ch %d\n",worse1);
1244 //printf("Kill ch %d\n",worse2);
1248 for (adcch = 0; adcch < 19; adcch++) {
1249 if (qTotal[adcch] > 0) // the channel is marked for processing
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
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));
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;
1266 // Calculate the center of gravity
1267 // checking for adcCentral != 0 (in case of "bad" configuration)
1268 if (adcCentral == 0)
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;
1276 // label calculation
1278 if (fDigitsManager) {
1279 Int_t label[9] = { 0 }; // up to 9 different labels possible
1280 Int_t count[9] = { 0 };
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");
1294 AliTRDarrayDictionary *dict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict);
1295 if (dict->GetDim() == 0) {
1296 AliError("Dictionary has dim. 0");
1300 for (Int_t iPad = 0; iPad < 3; iPad++) {
1301 if (padcol[iPad] < 0)
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]) {
1308 if (count[iLabel] > maxCount) {
1309 maxCount = count[iLabel];
1316 if (currLabel > 0) {
1317 label[nLabels++] = currLabel;
1322 mcLabel = label[maxIdx];
1325 // add the hit to the fitregister
1326 AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, mcLabel);
1332 void AliTRDmcmSim::TrackletSelection()
1334 // Select up to 4 tracklet candidates from the fit registers
1335 // and assign them to the CPUs.
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
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)))
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]);
1353 // for (i=0; i<ntracks;i++) printf("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]);
1357 // primitive sorting according to the number of hits
1358 for (j = 0; j < (ntracks-1); j++)
1360 for (i = j+1; i < ntracks; i++)
1362 if ( (trackletCand[j][1] < trackletCand[i][1]) ||
1363 ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
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;
1375 ntracks = 4; // cut the rest, 4 is the max
1377 // else is not necessary to sort
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++)
1382 for (i = j+1; i < ntracks; i++)
1384 if (trackletCand[j][0] < trackletCand[i][0])
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;
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]);
1405 void AliTRDmcmSim::FitTracklet()
1407 // Perform the actual tracklet fit based on the fit sums
1408 // which have been filled in the fit registers.
1410 // parameters in fitred.asm (fit program)
1411 Int_t decPlaces = 5;
1414 rndAdd = (1 << (decPlaces-1)) + 1;
1415 else if (decPlaces == 1)
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();
1439 bz = 0.1 * fld->SolenoidField(); // kGauss -> Tesla
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);
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
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};
1467 for (Int_t cpu = 0; cpu < 4; cpu++) {
1468 if (fFitPtr[cpu] == 31)
1470 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1474 fit0 = &fFitReg[fFitPtr[cpu] ];
1475 fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1478 mult = mult << (32 + decPlaces);
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;
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;
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
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));
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)))
1512 // printf("rejected\n");
1513 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1517 // printf("accepted\n");
1519 temp = temp * scaleD;
1520 slope = (temp >> 32);
1521 // printf("slope after scaling: %i\n", slope);
1524 temp = temp * scaleY;
1525 offset = (temp >> 32);
1527 // rounding, like in the TRAP
1528 slope = (slope + rndAdd) >> decPlaces;
1529 // printf("slope after shifting: %i\n", slope);
1530 offset = (offset + rndAdd) >> decPlaces;
1532 if (slope > 63) { // wrapping in TRAP!
1533 AliError(Form("Overflow in slope: %i, tracklet discarded!", slope));
1534 fMCMT[cpu] = 0x10001000;
1537 else if (slope < -64) {
1538 AliError(Form("Underflow in slope: %i, tracklet discarded!", slope));
1539 fMCMT[cpu] = 0x10001000;
1543 slope = slope & 0x7F; // 7 bit
1545 // printf("slope after clipping: 0x%02x\n", slope);
1547 if (offset > 0xfff || offset < -0xfff)
1548 AliWarning("Overflow in offset");
1549 offset = offset & 0x1FFF; // 13 bit
1551 qTotal = (q1 / nHits) >> 1;
1553 AliWarning("Overflow in charge");
1554 qTotal = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1556 // assemble and store the tracklet word
1557 fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1559 // calculate MC label
1561 if (fDigitsManager) {
1562 Int_t label[30] = {0}; // up to 30 different labels possible
1563 Int_t count[30] = {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))
1571 Int_t currLabel = fHits[iHit].fLabel;
1572 for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1573 if (currLabel == label[iLabel]) {
1575 if (count[iLabel] > maxCount) {
1576 maxCount = count[iLabel];
1583 if (currLabel > 0) {
1584 label[nLabels++] = currLabel;
1588 mcLabel = label[maxIdx];
1590 new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1591 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
1597 void AliTRDmcmSim::Tracklet()
1599 // Run the tracklet calculation by calling sequentially:
1600 // CalcFitreg(); TrackletSelection(); FitTracklet()
1601 // and store the tracklets
1603 if (!fInitialized) {
1604 AliError("Called uninitialized! Nothing done!");
1608 fTrackletArray->Delete();
1613 TrackletSelection();
1617 Bool_t AliTRDmcmSim::StoreTracklets()
1619 if (fTrackletArray->GetEntriesFast() == 0)
1622 AliRunLoader *rl = AliRunLoader::Instance();
1623 AliDataLoader *dl = 0x0;
1625 dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1627 AliError("Could not get the tracklets data loader!");
1631 TTree *trackletTree = dl->Tree();
1632 if (!trackletTree) {
1634 trackletTree = dl->Tree();
1637 AliTRDtrackletMCM *trkl = 0x0;
1638 TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1640 trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
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());
1648 dl->WriteData("OVERWRITE");
1653 void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
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
1659 if (!fInitialized) {
1660 AliError("Called uninitialized! Nothing done!");
1665 Int_t lastAdc = fNADC - 1;
1667 while (GetCol(firstAdc) < 0)
1670 while (GetCol(lastAdc) < 0)
1673 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
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));
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);
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);
1702 // help functions, to be cleaned up
1704 UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
1707 // This function adds a and b (unsigned) and clips to
1708 // the specified number of bits.
1714 UInt_t maxv = (1 << nbits) - 1;;
1720 if ((sum < a) || (sum < b))
1726 void 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
1731 // sorting for tracklet selection
1749 void 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)
1754 // sorting for tracklet selection
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);
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!
1775 case 4 : // 1 > 2, 2 <= 3, 3 <= 1 => 1 3 2
1784 case 2 : // 1 <= 2, 2 > 3, 3 <= 1 => 2 1 3
1793 case 3 : // 1 <= 2, 2 > 3, 3 > 1 => 2 3 1
1802 case 1 : // 1 <= 2, 2 <= 3, 3 > 1 => 3 2 1
1811 case 5 : // 1 > 2, 2 <= 3, 3 > 1 => 3 1 2
1820 default: // the rest should NEVER happen!
1821 AliError("ERROR in Sort3!!!\n");
1824 // printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
1827 void 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)
1832 // sorting for tracklet selection
1834 UShort_t idx21s, idx22s, idx23s, dummy;
1835 UShort_t val21s, val22s, val23s;
1836 UShort_t idx23as, idx23bs;
1837 UShort_t val23as, val23bs;
1839 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1840 idx1o, &idx21s, &idx23as,
1841 val1o, &val21s, &val23as);
1843 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1844 idx2o, &idx22s, &idx23bs,
1845 val2o, &val22s, &val23bs);
1847 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
1849 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1850 idx3o, idx4o, &dummy,
1851 val3o, val4o, &dummy);
1855 void 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)
1859 // sorting for tracklet selection
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;
1866 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1867 &dummy1, &idx21s, &idx23as,
1868 &dummy2, &val21s, &val23as);
1870 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1871 &dummy1, &idx22s, &idx23bs,
1872 &dummy2, &val22s, &val23bs);
1874 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
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);