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