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