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