TRD/AliTRDmcmSim.cxx

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