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