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