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