TRD/AliTRDmcmSim.cxx

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