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 //                                                                           //
  20 ///////////////////////////////////////////////////////////////////////////////
  21
  22 /* $Id$ */
  23
  24 /*
  25
  26   New release on 2007/08/17
  27
  28 AliTRDmcmSim is now stably working and zero suppression function seems ok.
  29 From now, the default version of raw data is set to 3 in AliTRDfeeParam.
  30
  31 The following internal parameters were abolished because it is useless and
  32 made trouble:
  33
  34    fColOfADCbeg
  35    fColOfADCend
  36
  37 GetCol member was modified accordingly.
  38
  39 New member function DumpData was prepared for diagnostics.
  40
  41 ZSMapping member function was debugged. It was causing crash due to
  42 wrong indexing in 1 dimensional numbering. Also code was shaped up better.
  43
  44 */
  45
  46 /*Semi-final version of TRD raw data simulation code with zero suppression (ZS)
  47 similar to TRD FEE. ZS is realized by the class group:
  48
  49   AliTRDfeeParam
  50   AliTRDmcmSim
  51   AliTRDrawData
  52
  53 AliTRDfeeParam has been modified to have more parameters like raw data
  54 production version and so on. AliTRDmcmSim is new class and this is the core
  55 of MCM (PASA+TRAP) simulator. It has still very simple function and it will be
  56 another project to improve this to make it closer to the reall FEE.
  57 AliTRDrawData has been modified to use new class AliTRDmcmSim.
  58
  59 These modifications were tested on Aug. 02 HEAD version that code itself
  60 compiles. I'm sure there must be still bugs and we need testing by as many as
  61 possible persons now. Especially it seems HLT part is impacted by problems
  62 because some parameters were moved from AliTRDrawData to AliTRDfeeParam (like
  63 fRawVersion disappeared from AliTRDrawData).
  64
  65 In TRD definition, we have now 4 raw data versions.
  66
  67   0 very old offline version (by Bogdan)
  68   1 test version (no zero suppression)
  69   2 final version (no zero suppression)
  70   3 test version (with zero suppression)
  71
  72 The default is still set to 2 in AliTRDfeeParam::fgkRAWversion and it uses
  73 previously existing codes. If you set this to 3, AliTRDrawData changes behavior
  74 to use AliTRDmcmSim with ZS.
  75
  76 Plan is after we make sure it works stably, we delete AliTRDmcm which is obsolete.
  77 However it still take time because tracklet part is not yet touched.
  78 The default raw version is 2.
  79
  80                                                                  Ken Oyama
  81 */
  82
  83 // if no histo is drawn, these are obsolete
  84 #include <TH1.h>
  85 #include <TCanvas.h>
  86
  87 // only needed if I/O of tracklets is activated
  88 #include <TObject.h>
  89 #include <TFile.h>
  90 #include <TTree.h>
  91 #include <TSystem.h>
  92
  93 #include <fstream>
  94
  95 #include <TMath.h>
  96
  97 #include "AliLog.h"
  98
  99 #include "AliTRDmcmSim.h"
 100 #include "AliTRDfeeParam.h"
 101 #include "AliTRDSimParam.h"
 102 #include "AliTRDgeometry.h"
 103 #include "AliTRDcalibDB.h"
 104 #include "AliTRDdigitsManager.h"
 105
 106 // additional for new tail filter and/or tracklet
 107 #include "AliTRDtrapAlu.h"
 108 #include "AliTRDpadPlane.h"
 109 #include "AliTRDtrackletMCM.h"
 110
 111 #include "AliRun.h"
 112 #include "AliLoader.h"
 113
 114 ClassImp(AliTRDmcmSim)
 115
 116 //_____________________________________________________________________________
 117 AliTRDmcmSim::AliTRDmcmSim() :TObject()
 118   ,fInitialized(kFALSE)
 119   ,fNextEvent(-1)
 120   ,fMaxTracklets(-1)
 121   ,fChaId(-1)
 122   ,fSector(-1)
 123   ,fStack(-1)
 124   ,fLayer(-1)
 125   ,fRobPos(-1)
 126   ,fMcmPos(-1)
 127   ,fNADC(-1)
 128   ,fNTimeBin(-1)
 129   ,fRow (-1)
 130   ,fADCR(NULL)
 131   ,fADCF(NULL)
 132   ,fADCT(NULL)
 133   ,fPosLUT(NULL)
 134   ,fMCMT(NULL)
 135   ,fZSM(NULL)
 136   ,fZSM1Dim(NULL)
 137   ,fFeeParam(NULL)
 138   ,fSimParam(NULL)
 139   ,fCal(NULL)
 140   ,fGeo(NULL)
 141 {
 142   //
 143   // AliTRDmcmSim default constructor
 144   //
 145
 146   // By default, nothing is initialized.
 147   // It is necessary to issue Init before use.
 148 }
 149
 150 //_____________________________________________________________________________
 151 AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m)
 152   :TObject(m)
 153   ,fInitialized(kFALSE)
 154   ,fNextEvent(-1)
 155   ,fMaxTracklets(-1)
 156   ,fChaId(-1)
 157   ,fSector(-1)
 158   ,fStack(-1)
 159   ,fLayer(-1)
 160   ,fRobPos(-1)
 161   ,fMcmPos(-1)
 162   ,fNADC(-1)
 163   ,fNTimeBin(-1)
 164   ,fRow(-1)
 165   ,fADCR(NULL)
 166   ,fADCF(NULL)
 167   ,fADCT(NULL)
 168   ,fPosLUT(NULL)
 169   ,fMCMT(NULL)
 170   ,fZSM(NULL)
 171   ,fZSM1Dim(NULL)
 172   ,fFeeParam(NULL)
 173   ,fSimParam(NULL)
 174   ,fCal(NULL)
 175   ,fGeo(NULL)
 176
 177 {
 178   //
 179   // AliTRDmcmSim copy constructor
 180   //
 181
 182   // By default, nothing is initialized.
 183   // It is necessary to issue Init before use.
 184 }
 185
 186 //_____________________________________________________________________________
 187 AliTRDmcmSim::~AliTRDmcmSim()
 188 {
 189   //
 190   // AliTRDmcmSim destructor
 191   //
 192
 193   if( fADCR != NULL ) {
 194     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 195       delete [] fADCR[iadc];
 196       delete [] fADCF[iadc];
 197       delete [] fADCT[iadc];
 198       delete [] fZSM [iadc];
 199     }
 200     delete [] fADCR;
 201     delete [] fADCF;
 202     delete [] fADCT;
 203     delete [] fZSM;
 204     delete [] fZSM1Dim;
 205   }
 206
 207   if(fInitialized){
 208     delete [] fPosLUT;
 209     delete [] fMCMT;
 210   }
 211
 212   delete fGeo;
 213
 214 }
 215
 216 //_____________________________________________________________________________
 217 AliTRDmcmSim &AliTRDmcmSim::operator=(const AliTRDmcmSim &m)
 218 {
 219   //
 220   // Assignment operator
 221   //
 222
 223   if (this != &m) {
 224     ((AliTRDmcmSim &) m).Copy(*this);
 225   }
 226   return *this;
 227
 228 }
 229
 230 //_____________________________________________________________________________
 231 void AliTRDmcmSim::Copy(TObject &m) const
 232 {
 233   //
 234   // Copy function
 235   //
 236   ((AliTRDmcmSim &) m).fNextEvent     = 0; //new
 237   ((AliTRDmcmSim &) m).fMaxTracklets  = 0; //new
 238   ((AliTRDmcmSim &) m).fInitialized   = 0;
 239   ((AliTRDmcmSim &) m).fChaId         = 0;
 240   ((AliTRDmcmSim &) m).fSector        = 0;
 241   ((AliTRDmcmSim &) m).fStack         = 0;
 242   ((AliTRDmcmSim &) m).fLayer         = 0;
 243   ((AliTRDmcmSim &) m).fRobPos        = 0;
 244   ((AliTRDmcmSim &) m).fMcmPos        = 0;
 245   ((AliTRDmcmSim &) m).fNADC          = 0;
 246   ((AliTRDmcmSim &) m).fNTimeBin      = 0;
 247   ((AliTRDmcmSim &) m).fRow           = 0;
 248   ((AliTRDmcmSim &) m).fADCR          = 0;
 249   ((AliTRDmcmSim &) m).fADCF          = 0;
 250   ((AliTRDmcmSim &) m).fADCT          = 0; //new
 251   ((AliTRDmcmSim &) m).fPosLUT        = 0; //new
 252   ((AliTRDmcmSim &) m).fMCMT          = 0; //new
 253   ((AliTRDmcmSim &) m).fZSM           = 0;
 254   ((AliTRDmcmSim &) m).fZSM1Dim       = 0;
 255   ((AliTRDmcmSim &) m).fFeeParam      = 0;
 256   ((AliTRDmcmSim &) m).fSimParam      = 0;
 257   ((AliTRDmcmSim &) m).fCal           = 0;
 258   ((AliTRDmcmSim &) m).fGeo           = 0;
 259
 260 }
 261
 262 //_____________________________________________________________________________
 263
 264 //void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
 265 void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos, Bool_t newEvent = kFALSE ) // only for readout tree (new event)
 266 {
 267   //
 268   // Initialize the class with new geometry information
 269   // fADC array will be reused with filled by zero
 270   //
 271
 272   fNextEvent     = 0;
 273   fFeeParam      = AliTRDfeeParam::Instance();
 274   fSimParam      = AliTRDSimParam::Instance();
 275   fCal           = AliTRDcalibDB::Instance();
 276   fGeo           = new AliTRDgeometry();
 277   fChaId         = chaId;
 278   fSector        = fGeo->GetSector( fChaId );
 279   fStack         = fGeo->GetStack( fChaId );
 280   fLayer         = fGeo->GetLayer( fChaId );
 281   fRobPos        = robPos;
 282   fMcmPos        = mcmPos;
 283   fNADC          = fFeeParam->GetNadcMcm();
 284   fNTimeBin      = fCal->GetNumberOfTimeBins();
 285   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
 286
 287   fMaxTracklets  = fFeeParam->GetMaxNrOfTracklets();
 288
 289
 290
 291
 292   if (newEvent == kTRUE) {
 293       fNextEvent = 1;
 294   }
 295
 296
 297
 298   // Allocate ADC data memory if not yet done
 299   if( fADCR == NULL ) {
 300     fADCR    = new Int_t *[fNADC];
 301     fADCF    = new Int_t *[fNADC];
 302     fADCT    = new Int_t *[fNADC]; //new
 303     fZSM     = new Int_t *[fNADC];
 304     fZSM1Dim = new Int_t  [fNADC];
 305     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 306       fADCR[iadc] = new Int_t[fNTimeBin];
 307       fADCF[iadc] = new Int_t[fNTimeBin];
 308       fADCT[iadc] = new Int_t[fNTimeBin]; //new
 309       fZSM [iadc] = new Int_t[fNTimeBin];
 310     }
 311   }
 312
 313   // Initialize ADC data
 314   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 315     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 316       fADCR[iadc][it] = 0;
 317       fADCF[iadc][it] = 0;
 318       fADCT[iadc][it] = -1;  //new
 319       fZSM [iadc][it] = 1;   // Default unread = 1
 320     }
 321     fZSM1Dim[iadc] = 1;      // Default unread = 1
 322   }
 323
 324   //new:
 325   fPosLUT = new Int_t[128];
 326   for(Int_t i = 0; i<128; i++){
 327     fPosLUT[i] = 0;
 328   }
 329
 330   fMCMT = new UInt_t[fMaxTracklets];
 331   for(Int_t i = 0; i < fMaxTracklets; i++) {
 332     fMCMT[i] = 0;
 333   }
 334
 335
 336   fInitialized = kTRUE;
 337 }
 338
 339 //_____________________________________________________________________________
 340 Bool_t AliTRDmcmSim::CheckInitialized()
 341 {
 342   //
 343   // Check whether object is initialized
 344   //
 345
 346   if( ! fInitialized ) {
 347     AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
 348   }
 349   return fInitialized;
 350 }
 351
 352 //_____________________________________________________________________________
 353
 354
 355 void AliTRDmcmSim::SetPosLUT() {
 356   Double_t iHi  = (Double_t)fCal->GetPRFhi();
 357   Double_t iLo  = (Double_t)fCal->GetPRFlo();
 358   Int_t   nBin  = fCal->GetPRFbin();
 359   Int_t   iOff  = fLayer * nBin;
 360   Int_t kNlayer = fGeo->Nlayer();
 361
 362   Float_t  *sPRFsmp   = new Float_t[nBin*kNlayer];
 363   Double_t *sPRFlayer = new Double_t[nBin];
 364
 365
 366   for(Int_t i = 0; i<nBin*kNlayer; i++){
 367
 368     //printf("%f\n",fCal->GetSampledPRF()[i]);
 369     sPRFsmp[i] = fCal->GetSampledPRF()[i];
 370
 371   }
 372
 373   Double_t sWidth = (iHi-iLo)/((Double_t) nBin);
 374   Int_t   sPad    = (Int_t) (1.0/sWidth);
 375
 376   // get the PRF for actual layer (interpolated to ibin data-points; 61 measured)
 377   for(Int_t iBin = 0; iBin < nBin; iBin++){
 378     sPRFlayer[iBin] = (Double_t)sPRFsmp[iOff+iBin];
 379   }
 380
 381   Int_t bin0 = (Int_t)(-iLo / sWidth - 0.5);                           // bin-nr. for pad-position 0
 382
 383   Int_t bin1 = (Int_t)((Double_t)(0.5 - iLo) / sWidth - 0.5);          // bin-nr. for pad-position 0.5
 384   bin1 = bin1 + 1;
 385   bin0 = bin0 + 1;  //avoid negative values in aYest (start right of symmetry center)
 386   while (bin0-sPad<0) {
 387     bin0 = bin0 + 1;
 388   }
 389   while (bin1+sPad>=nBin) {
 390     bin1 = bin1 - 1;
 391   }
 392
 393   Double_t* aYest = new Double_t[bin1-bin0+1];
 394
 395   /*TH1F* hist1 = new TH1F("h1","yest(y)",128,0,0.5);
 396   TH1F* hist2 = new TH1F("h2","y(yest)",128,0,0.5);
 397   TH1F* hist3 = new TH1F("h3","y(yest)-yest",128,0,0.5);
 398   TH1F* hist4 = new TH1F("h4","y(yest)-yest,discrete",128,0,0.5);
 399
 400   TCanvas *c1 = new TCanvas("c1","c1",800,1000);
 401   hist1->Draw();
 402   TCanvas *c2 = new TCanvas("c2","c2",800,1000);
 403   hist2->Draw();
 404   TCanvas *c3 = new TCanvas("c3","c3",800,1000);
 405   hist3->Draw();
 406   TCanvas *c4 = new TCanvas("c4","c4",800,1000);
 407   hist4->Draw();*/
 408
 409   for(Int_t iBin = bin0; iBin <= bin1; iBin++){
 410     aYest[iBin-bin0] = 0.5*(sPRFlayer[iBin-sPad] - sPRFlayer[iBin+sPad])/(sPRFlayer[iBin]); // estimated position from PRF; between 0 and 1
 411     //Double_t position = ((Double_t)(iBin)+0.5)*sWidth+iLo;
 412     //  hist1->Fill(position,aYest[iBin-bin0]);
 413   }
 414
 415
 416
 417   Double_t aY[128]; // reversed function
 418
 419   AliTRDtrapAlu a;
 420   a.Init(1,8,0,31);
 421
 422   for(Int_t j = 0; j<128; j++) { // loop over all Yest; LUT has 128 entries;
 423     Double_t yest = ((Double_t)j)/256;
 424
 425     Int_t iBin = 0;
 426     while (yest>aYest[iBin] && iBin<(bin1-bin0)) {
 427       iBin = iBin+1;
 428     }
 429     if((iBin == bin1 - bin0)&&(yest>aYest[iBin])) {
 430       aY[j] = 0.5;                      // yest too big
 431       //hist2->Fill(yest,aY[j]);
 432
 433     }
 434     else {
 435       Int_t bin_d = iBin + bin0 - 1;
 436       Int_t bin_u = iBin + bin0;
 437       Double_t y_d = ((Double_t)bin_d + 0.5)*sWidth + iLo; // lower y
 438       Double_t y_u = ((Double_t)bin_u + 0.5)*sWidth + iLo; // upper y
 439       Double_t yest_d = aYest[iBin-1];                     // lower estimated y
 440       Double_t yest_u = aYest[iBin];                       // upper estimated y
 441
 442       aY[j] = ((yest-yest_d)/(yest_u-yest_d))*(y_u-y_d) + y_d;
 443       //hist2->Fill(yest,aY[j]);
 444
 445     }
 446     aY[j] = aY[j] - yest;
 447     //hist3->Fill(yest,aY[j]);
 448     // formatting
 449     a.AssignDouble(aY[j]);
 450     //a.WriteWord();
 451     fPosLUT[j] = a.GetValue(); // 1+8Bit value;128 entries;LUT is steered by abs(Q(i+1)-Q(i-1))/Q(i)=COG and gives the correction to COG/2
 452     //hist4->Fill(yest,fPosLUT[j]);
 453
 454   }
 455
 456
 457
 458   delete [] sPRFsmp;
 459   delete [] sPRFlayer;
 460   delete [] aYest;
 461
 462 }
 463
 464
 465 //_____________________________________________________________________________
 466 Int_t* AliTRDmcmSim::GetPosLUT(){
 467   return fPosLUT;
 468 }
 469
 470
 471
 472 void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
 473 {
 474   //
 475   // Store ADC data into array of raw data
 476   //
 477
 478   if( !CheckInitialized() ) return;
 479
 480   if( iadc < 0 || iadc >= fNADC ) {
 481     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 482     return;
 483   }
 484
 485   for( int it = 0 ;  it < fNTimeBin ; it++ ) {
 486     fADCR[iadc][it] = (Int_t)(adc[it]);
 487   }
 488 }
 489
 490 //_____________________________________________________________________________
 491 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
 492 {
 493   //
 494   // Store ADC data into array of raw data
 495   //
 496
 497   if( !CheckInitialized() ) return;
 498
 499   if( iadc < 0 || iadc >= fNADC ) {
 500     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 501     return;
 502   }
 503
 504   fADCR[iadc][it] = adc;
 505 }
 506
 507 //_____________________________________________________________________________
 508 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
 509 {
 510   //
 511   // Store ADC data into array of raw data
 512   //
 513
 514   if( !CheckInitialized() ) return;
 515
 516   if( iadc < 0 || iadc >= fNADC ) {
 517     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 518     return;
 519   }
 520
 521   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 522     fADCR[iadc][it] = fSimParam->GetADCbaseline();
 523   }
 524 }
 525
 526 //_____________________________________________________________________________
 527 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
 528 {
 529   //
 530   // Return column id of the pad for the given ADC channel
 531   //
 532
 533   if( !CheckInitialized() ) return -1;
 534
 535   return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
 536 }
 537
 538 //_____________________________________________________________________________
 539 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
 540 {
 541   //
 542   // Produce raw data stream from this MCM and put in buf
 543   // Returns number of words filled, or negative value
 544   // with -1 * number of overflowed words
 545   //
 546
 547   UInt_t  x;
 548   UInt_t  iEv = 0;
 549   Int_t   nw  = 0;  // Number of written words
 550   Int_t   of  = 0;  // Number of overflowed words
 551   Int_t   rawVer   = fFeeParam->GetRAWversion();
 552   Int_t **adc;
 553
 554   if( !CheckInitialized() ) return 0;
 555
 556   if( fFeeParam->GetRAWstoreRaw() ) {
 557     adc = fADCR;
 558   } else {
 559     adc = fADCF;
 560   }
 561
 562   // Produce MCM header
 563   x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
 564   if (nw < maxSize) {
 565     buf[nw++] = x;
 566   }
 567   else {
 568     of++;
 569   }
 570
 571   // Produce ADC mask
 572   if( rawVer >= 3 ) {
 573     x = 0;
 574     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
 575       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
 576         x = x | (1 << iAdc);
 577       }
 578     }
 579     if (nw < maxSize) {
 580       buf[nw++] = x;
 581     }
 582     else {
 583       of++;
 584     }
 585   }
 586
 587   // Produce ADC data. 3 timebins are packed into one 32 bits word
 588   // In this version, different ADC channel will NOT share the same word
 589
 590   UInt_t aa=0, a1=0, a2=0, a3=0;
 591
 592   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
 593     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // suppressed
 594     aa = !(iAdc & 1) + 2;
 595     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
 596       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] : 0;
 597       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
 598       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
 599       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
 600       if (nw < maxSize) {
 601         buf[nw++] = x;
 602       }
 603       else {
 604         of++;
 605       }
 606     }
 607   }
 608
 609   if( of != 0 ) return -of; else return nw;
 610 }
 611
 612 //_____________________________________________________________________________
 613 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
 614 {
 615   //
 616   // Produce tracklet data stream from this MCM and put in buf
 617   // Returns number of words filled, or negative value
 618   // with -1 * number of overflowed words
 619   //
 620
 621   UInt_t  x;
 622   Int_t   nw  = 0;  // Number of written words
 623   Int_t   of  = 0;  // Number of overflowed words
 624
 625   if( !CheckInitialized() ) return 0;
 626
 627   // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM
 628   // fMCMT is filled continuously until no more tracklet words available
 629
 630   Int_t wd = 0;
 631   while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){
 632       x = fMCMT[wd];
 633       if (nw < maxSize) {
 634         buf[nw++] = x;
 635       }
 636       else {
 637         of++;
 638       }
 639       wd++;
 640   }
 641
 642   if( of != 0 ) return -of; else return nw;
 643 }
 644
 645
 646 //_____________________________________________________________________________
 647 void AliTRDmcmSim::Filter()
 648 {
 649   //
 650   // Apply digital filter
 651   //
 652
 653   if( !CheckInitialized() ) return;
 654
 655   // Initialize filtered data array with raw data
 656   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 657     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 658       fADCF[iadc][it] = fADCR[iadc][it];
 659     }
 660   }
 661
 662   // Then apply fileters one by one to filtered data array
 663   if( fFeeParam->IsPFon() ) FilterPedestal();
 664   if( fFeeParam->IsGFon() ) FilterGain();
 665   if( fFeeParam->IsTFon() ) FilterTail();
 666 }
 667
 668 //_____________________________________________________________________________
 669 void AliTRDmcmSim::FilterPedestal()
 670 {
 671
 672
 673   //
 674   // Apply pedestal filter
 675   //
 676
 677   Int_t ap = fSimParam->GetADCbaseline();      // ADC instrinsic pedestal
 678   Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
 679   //Int_t tc = fFeeParam->GetPFtimeConstant();   // this makes no sense yet
 680
 681   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 682     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 683       fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
 684     }
 685   }
 686 }
 687
 688 //_____________________________________________________________________________
 689 void AliTRDmcmSim::FilterGain()
 690 {
 691   //
 692   // Apply gain filter (not implemented)
 693   // Later it will be implemented because gain digital filiter will
 694   // increase noise level.
 695   //
 696
 697 }
 698
 699 //_____________________________________________________________________________
 700 void AliTRDmcmSim::FilterTail()
 701 {
 702   //
 703   // Apply exponential tail filter (Bogdan's version)
 704   //
 705
 706   Double_t *dtarg  = new Double_t[fNTimeBin];
 707   Int_t    *itarg  = new Int_t[fNTimeBin];
 708   Int_t     nexp   = fFeeParam->GetTFnExp();
 709   Int_t     tftype = fFeeParam->GetTFtype();
 710
 711   switch( tftype ) {
 712
 713   case 0: // Exponential Filter Analog Bogdan
 714     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 715       FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
 716       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 717         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
 718       }
 719     }
 720     break;
 721
 722   case 1: // Exponential filter digital Bogdan
 723     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 724       FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
 725       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 726         fADCF[iCol][iTime] = itarg[iTime];
 727       }
 728     }
 729     break;
 730
 731   case 2: // Exponential filter Marian special
 732     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 733       FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
 734       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 735         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
 736       }
 737     }
 738     break;
 739
 740     //new
 741   case 3: // Exponential filter using AliTRDtrapAlu class
 742     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 743       FilterSimDeConvExpEl( fADCF[iCol], itarg, fNTimeBin, nexp);
 744       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 745         fADCF[iCol][iTime] = itarg[iTime]>>2; // to be used for raw-data
 746         fADCT[iCol][iTime] = itarg[iTime];    // 12bits; to be used for tracklet; tracklet will have own container;
 747       }
 748     }
 749     break;
 750
 751
 752   default:
 753     AliError(Form("Invalid filter type %d ! \n", tftype ));
 754     break;
 755   }
 756
 757   delete dtarg;
 758   delete itarg;
 759 }
 760
 761 //_____________________________________________________________________________
 762 void AliTRDmcmSim::ZSMapping()
 763 {
 764   //
 765   // Zero Suppression Mapping implemented in TRAP chip
 766   //
 767   // See detail TRAP manual "Data Indication" section:
 768   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
 769   //
 770
 771   Int_t eBIS = fFeeParam->GetEBsglIndThr();       // TRAP default = 0x4  (Tis=4)
 772   Int_t eBIT = fFeeParam->GetEBsumIndThr();       // TRAP default = 0x28 (Tit=40)
 773   Int_t eBIL = fFeeParam->GetEBindLUT();          // TRAP default = 0xf0
 774                                                   // (lookup table accept (I2,I1,I0)=(111)
 775                                                   // or (110) or (101) or (100))
 776   Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
 777   Int_t ep   = AliTRDfeeParam::GetPFeffectPedestal();
 778
 779   if( !CheckInitialized() ) return;
 780
 781   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
 782     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 783
 784       // Get ADC data currently in filter buffer
 785       Int_t ap = fADCF[iadc-1][it] - ep; // previous
 786       Int_t ac = fADCF[iadc  ][it] - ep; // current
 787       Int_t an = fADCF[iadc+1][it] - ep; // next
 788
 789       // evaluate three conditions
 790       Int_t i0 = ( ac >=  ap && ac >=  an ) ? 0 : 1; // peak center detection
 791       Int_t i1 = ( ap + ac + an > eBIT )    ? 0 : 1; // cluster
 792       Int_t i2 = ( ac > eBIS )              ? 0 : 1; // absolute large peak
 793
 794       Int_t i = i2 * 4 + i1 * 2 + i0;    // Bit position in lookup table
 795       Int_t d = (eBIL >> i) & 1;         // Looking up  (here d=0 means true
 796                                          // and d=1 means false according to TRAP manual)
 797
 798       fZSM[iadc][it] &= d;
 799       if( eBIN == 0 ) {  // turn on neighboring ADCs
 800         fZSM[iadc-1][it] &= d;
 801         fZSM[iadc+1][it] &= d;
 802       }
 803
 804     }
 805   }
 806
 807   // do 1 dim projection
 808   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 809     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 810       fZSM1Dim[iadc] &= fZSM[iadc][it];
 811     }
 812   }
 813
 814 }
 815
 816 //_____________________________________________________________________________
 817 void AliTRDmcmSim::DumpData( char *f, char *target )
 818 {
 819   //
 820   // Dump data stored (for debugging).
 821   // target should contain one or multiple of the following characters
 822   //   R   for raw data
 823   //   F   for filtered data
 824   //   Z   for zero suppression map
 825   //   S   Raw dat astream
 826   // other characters are simply ignored
 827   //
 828
 829   UInt_t tempbuf[1024];
 830
 831   if( !CheckInitialized() ) return;
 832
 833   std::ofstream of( f, std::ios::out | std::ios::app );
 834   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
 835              fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos );
 836
 837   for( int t=0 ; target[t] != 0 ; t++ ) {
 838     switch( target[t] ) {
 839     case 'R' :
 840     case 'r' :
 841       of << Form("fADCR (raw ADC data)\n");
 842       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 843         of << Form("  ADC %02d: ", iadc);
 844         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 845           of << Form("% 4d",  fADCR[iadc][it]);
 846         }
 847         of << Form("\n");
 848       }
 849       break;
 850     case 'F' :
 851     case 'f' :
 852       of << Form("fADCF (filtered ADC data)\n");
 853       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 854         of << Form("  ADC %02d: ", iadc);
 855         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 856           of << Form("% 4d",  fADCF[iadc][it]);
 857         }
 858         of << Form("\n");
 859       }
 860       break;
 861     case 'Z' :
 862     case 'z' :
 863       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
 864       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 865         of << Form("  ADC %02d: ", iadc);
 866         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
 867         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 868           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
 869         }
 870         of << Form("\n");
 871       }
 872       break;
 873     case 'S' :
 874     case 's' :
 875       Int_t s = ProduceRawStream( tempbuf, 1024 );
 876       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
 877       of << Form("  address  data\n");
 878       for( int i = 0 ; i < s ; i++ ) {
 879         of << Form("  %04x     %08x\n", i, tempbuf[i]);
 880       }
 881     }
 882   }
 883 }
 884
 885 //_____________________________________________________________________________
 886 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
 887                                      , Int_t n, Int_t nexp)
 888 {
 889   //
 890   // Exponential filter "analog"
 891   // source will not be changed
 892   //
 893
 894   Int_t    i = 0;
 895   Int_t    k = 0;
 896   Double_t reminder[2];
 897   Double_t correction;
 898   Double_t result;
 899   Double_t rates[2];
 900   Double_t coefficients[2];
 901
 902   // Initialize (coefficient = alpha, rates = lambda)
 903   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 904
 905   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
 906   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
 907   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
 908   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
 909
 910   coefficients[0] = c1;
 911   coefficients[1] = c2;
 912
 913   Double_t dt = 0.1;
 914   rates[0] = TMath::Exp(-dt/(r1));
 915   rates[1] = TMath::Exp(-dt/(r2));
 916
 917   // Attention: computation order is important
 918   correction = 0.0;
 919   for (k = 0; k < nexp; k++) {
 920     reminder[k] = 0.0;
 921   }
 922
 923   for (i = 0; i < n; i++) {
 924
 925     result    = ((Double_t)source[i] - correction);    // no rescaling
 926     target[i] = result;
 927
 928     for (k = 0; k < nexp; k++) {
 929       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
 930     }
 931
 932     correction = 0.0;
 933     for (k = 0; k < nexp; k++) {
 934       correction += reminder[k];
 935     }
 936   }
 937 }
 938
 939 //_____________________________________________________________________________
 940 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
 941                                      , Int_t nexp)
 942 {
 943   //
 944   // Exponential filter "digital"
 945   // source will not be changed
 946   //
 947
 948   Int_t i        = 0;
 949   Int_t fAlphaL  = 0;
 950   Int_t fAlphaS  = 0;
 951   Int_t fTailPed = 0;
 952   Int_t iAlphaL  = 0;
 953   Int_t iAlphaS  = 0;
 954
 955   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 956   // initialize (coefficient = alpha, rates = lambda)
 957
 958   Double_t dt = 0.1;
 959   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
 960   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
 961   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
 962   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
 963
 964   Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
 965   Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
 966   Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; //  9 bit paramter + fixed bits
 967   Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; //  9 bit paramter + fixed bits
 968
 969   if (nexp == 1) {
 970     fAlphaL = (Int_t) (c1 * 2048.0);
 971     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 972   }
 973   if (nexp == 2) {
 974     fAlphaL = (Int_t) (c1 * 2048.0);
 975     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
 976     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 977     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
 978   }
 979
 980   Double_t iAl = iAlphaL  / 2048.0;            // alpha  L: correspondence to floating point numbers
 981   Double_t iAs = iAlphaS  / 2048.0;            // alpha  S: correspondence to floating point numbers
 982   Double_t iLl = iLambdaL / 2048.0;            // lambda L: correspondence to floating point numbers
 983   Double_t iLs = iLambdaS / 2048.0;            // lambda S: correspondence to floating point numbers
 984
 985   Int_t h1;
 986   Int_t h2;
 987   Int_t rem1;
 988   Int_t rem2;
 989   Int_t correction;
 990   Int_t result;
 991   Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
 992
 993   Double_t xi = 1 - (iLl*iAs + iLs*iAl);             // Calculation of equilibrium values of the
 994   rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
 995   rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
 996
 997   // further initialization
 998   if ((rem1 + rem2) > 0x0FFF) {
 999     correction = 0x0FFF;
1000   }
1001   else {
1002     correction = (rem1 + rem2) & 0x0FFF;
1003   }
1004
1005   fTailPed = iFactor - correction;
1006
1007   for (i = 0; i < n; i++) {
1008
1009     result = (source[i]  - correction);
1010     if (result < 0) { // Too much undershoot
1011       result = 0;
1012     }
1013
1014     target[i] = result;
1015
1016     h1 = (rem1 + ((iAlphaL * result) >> 11));
1017     if (h1 > 0x0FFF) {
1018       h1 = 0x0FFF;
1019     }
1020     else {
1021       h1 &= 0x0FFF;
1022     }
1023
1024     h2 = (rem2 + ((iAlphaS * result) >> 11));
1025     if (h2 > 0x0FFF) {
1026       h2 = 0x0FFF;
1027     }
1028     else {
1029       h2 &= 0x0FFF;
1030     }
1031
1032     rem1 = (iLambdaL * h1 ) >> 11;
1033     rem2 = (iLambdaS * h2 ) >> 11;
1034
1035     if ((rem1 + rem2) > 0x0FFF) {
1036       correction = 0x0FFF;
1037     }
1038     else {
1039       correction = (rem1 + rem2) & 0x0FFF;
1040     }
1041
1042   }
1043
1044 }
1045
1046 //_____________________________________________________________________________
1047 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
1048                                       , Int_t n)
1049 {
1050   //
1051   // Exponential filter (M. Ivanov)
1052   // source will not be changed
1053   //
1054
1055   Int_t i = 0;
1056   Double_t sig1[100];
1057   Double_t sig2[100];
1058   Double_t sig3[100];
1059
1060   for (i = 0; i < n; i++) {
1061     sig1[i] = (Double_t)source[i];
1062   }
1063
1064   Float_t dt      = 0.1;
1065   Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
1066   Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
1067
1068   FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
1069   FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
1070
1071   for (i = 0; i < n; i++) {
1072     target[i] = sig3[i];
1073   }
1074
1075 }
1076
1077 //______________________________________________________________________________
1078 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
1079                                           , Double_t lambda, Int_t n)
1080 {
1081   //
1082   // Special filter (M. Ivanov)
1083   //
1084
1085   Int_t    i = 0;
1086   Double_t l = TMath::Exp(-lambda*0.5);
1087   Double_t in[1000];
1088   Double_t out[1000];
1089
1090   // Initialize in[] and out[] goes 0 ... 2*n+19
1091   for (i = 0; i < n*2+20; i++) {
1092     in[i] = out[i] = 0;
1093   }
1094
1095   // in[] goes 0, 1
1096   in[0] = ampin[0];
1097   in[1] = (ampin[0] + ampin[1]) * 0.5;
1098
1099   // Add charge to the end
1100   for (i = 0; i < 22; i++) {
1101     in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
1102   }
1103
1104   // Use arithmetic mean
1105   for (i = 1; i < n-1; i++) {
1106     in[2*i]   = ampin[i];    // in[] goes 2, 3, ... , 2*n-4, 2*n-3
1107     in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
1108   }
1109
1110   Double_t temp;
1111   out[2*n]    = in[2*n];
1112   temp        = 0;
1113   for (i = 2*n; i >= 0; i--) {
1114     out[i]    = in[i] + temp;
1115     temp      = l*(temp+in[i]);
1116   }
1117
1118   for (i = 0; i < n; i++){
1119     //ampout[i] = out[2*i+1];  // org
1120     ampout[i] = out[2*i];
1121   }
1122
1123 }
1124
1125 //______________________________________________________________________________
1126 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
1127                                             , Double_t norm, Double_t lambda
1128                                             , Int_t n)
1129 {
1130   //
1131   // Special filter (M. Ivanov)
1132   //
1133
1134   Int_t    i = 0;
1135
1136   Double_t l = TMath::Exp(-lambda*0.5);
1137   Double_t k = l*(1.0 - norm*lambda*0.5);
1138   Double_t in[1000];
1139   Double_t out[1000];
1140
1141   // Initialize in[] and out[] goes 0 ... 2*n+19
1142   for (i = 0; i < n*2+20; i++) {
1143     in[i] = out[i] = 0;
1144   }
1145
1146   // in[] goes 0, 1
1147   in[0] = ampin[0];
1148   in[1] = (ampin[0]+ampin[1])*0.5;
1149
1150   // Add charge to the end
1151   for (i =-2; i < 22; i++) {
1152     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
1153     in[2*(n-1)+i] = ampin[n-1];
1154   }
1155
1156   for (i = 1; i < n-2; i++) {
1157     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
1158     in[2*i]    = ampin[i];
1159     in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
1160     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
1161   }
1162
1163   Double_t temp;
1164   out[0] = in[0];
1165   temp   = in[0];
1166   for (i = 1; i <= 2*n; i++) {
1167     out[i] = in[i] + (k-l)*temp;
1168     temp   = in[i] +  k   *temp;
1169   }
1170
1171   for (i = 0; i < n; i++) {
1172     //ampout[i] = out[2*i+1];  // org
1173     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
1174     ampout[i] = TMath::Max(out[2*i],0.0);
1175   }
1176 }
1177
1178
1179 //_____________________________________________________________________________________
1180 //the following filter uses AliTRDtrapAlu-class
1181
1182 void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) {
1183   //static Int_t count = 0;
1184
1185   Double_t dt = 0.1;
1186   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1187   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1188   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1189   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1190
1191   nexp = 1;
1192
1193   //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual
1194   //parameters need to be adjusted
1195   AliTRDtrapAlu lambdaL;
1196   AliTRDtrapAlu lambdaS;
1197   AliTRDtrapAlu alphaL;
1198   AliTRDtrapAlu alphaS;
1199
1200   AliTRDtrapAlu correction;
1201   AliTRDtrapAlu result;
1202   AliTRDtrapAlu bufL;
1203   AliTRDtrapAlu bufS;
1204
1205   AliTRDtrapAlu bSource;
1206
1207   lambdaL.Init(1,11);
1208   lambdaS.Init(1,11);
1209   alphaL.Init(1,11);
1210   alphaS.Init(1,11);
1211
1212   //count=count+1;
1213
1214   lambdaL.AssignDouble(TMath::Exp(-dt/r1));
1215   lambdaS.AssignDouble(TMath::Exp(-dt/r2));
1216   alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants
1217   alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL
1218
1219   //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly
1220   correction.Init(10,2);
1221   result.Init(10,2);
1222   bufL.Init(10,2);
1223   bufS.Init(10,2);
1224   bSource.Init(10,2);
1225
1226   for(Int_t i = 0; i < n; i++) {
1227     bSource.AssignInt(source[i]);
1228     result = bSource - correction; // subtraction can produce an underflow
1229     if(result.GetSign() == kTRUE) {
1230       result.AssignInt(0);
1231     }
1232
1233     //target[i] = result.GetValuePre();  // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly
1234
1235     target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2
1236
1237     //printf("target-Wert zur Zeit %d : %d",i,target[i]);
1238     //printf("\n");
1239
1240     bufL  =  bufL + (result * alphaL);
1241     bufL  =  bufL * lambdaL;
1242
1243     bufS  =  bufS + (result * alphaS);
1244     bufS  =  bufS * lambdaS;  // eventually this should look like:
1245     // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more
1246
1247     correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF
1248   }
1249
1250
1251 }
1252
1253
1254
1255
1256
1257
1258
1259 //__________________________________________________________________________________
1260
1261
1262 // in order to use the Tracklets, please first
1263 // -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
1264 // -- set AliTRDfeeParam::fgkTFtype   to 3, in order to use the new tail cancellation filter
1265 //    currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
1266 // -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with           information about the Tracklet position (MCM, channel number)
1267
1268 // The code is designed such that the less possible calculations with AliTRDtrapAlu class-objects are performed; whenever possible calculations are done with doubles or integers and the results are transformed into the right format
1269
1270 void AliTRDmcmSim::Tracklet(){
1271     // tracklet calculation
1272     // if you use this code after a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
1273
1274   if(!CheckInitialized()){ return; }
1275
1276   Bool_t filtered = kTRUE;
1277
1278
1279
1280   AliTRDtrapAlu data;
1281   data.Init(10,2);
1282   if(fADCT[0][0]==-1){                      // check if filter was applied
1283     filtered = kFALSE;
1284     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1285       for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) {
1286         data.AssignInt(fADCR[iadc][iT]);
1287         fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called, this is done correctly
1288       }
1289     }
1290
1291   }
1292
1293   // the online ordering of mcm's is reverse to the TRAP-manual-ordering! reverse fADCT (to be consistent to TRAP), then do all calculations
1294   // reverse fADCT:
1295   Int_t** rev0 = new Int_t *[fNADC];
1296   Int_t** rev1 = new Int_t *[fNADC];
1297
1298   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1299     rev0[iadc] = new Int_t[fNTimeBin];
1300     rev1[iadc] = new Int_t[fNTimeBin];
1301     for( Int_t iT = 0; iT < fNTimeBin; iT++) {
1302       if( iadc <= fNADC-iadc-1 ) {
1303         rev0[iadc][iT]  = fADCT[fNADC-iadc-1][iT];
1304         rev1[iadc][iT]  = fADCT[iadc][iT];
1305         fADCT[iadc][iT] = rev0[iadc][iT];
1306       }
1307       else {
1308         rev0[iadc][iT]  = rev1[fNADC-iadc-1][iT];
1309         fADCT[iadc][iT] = rev0[iadc][iT];
1310       }
1311     }
1312   }
1313   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1314     delete[] rev0[iadc];
1315     delete[] rev1[iadc];
1316   }
1317
1318   delete[] rev0;
1319   delete[] rev1;
1320
1321   rev0 = NULL;
1322   rev1 = NULL;
1323
1324   // get the filtered pedestal; supports only electronic tail-cancellation filter
1325   AliTRDtrapAlu filPed;
1326   Int_t ep = 0;
1327   Int_t *ieffped = new Int_t[fNTimeBin];
1328   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1329     ieffped[iT] = ep;
1330   }
1331
1332   if( filtered == kTRUE ) {
1333     if( fFeeParam->IsPFon() ){
1334       ep = fFeeParam->GetPFeffectPedestal();
1335     }
1336     Int_t      nexp  = fFeeParam->GetTFnExp();
1337     Int_t  *isource  = new Int_t[fNTimeBin];
1338     filPed.Init(10,2);
1339     filPed.AssignInt(ep);
1340     Int_t epf = filPed.GetValue();
1341     for(Int_t iT = 0; iT < fNTimeBin; iT++){
1342       isource[iT] = ep;
1343       ieffped[iT] = epf;
1344     }
1345
1346     if( fFeeParam->IsTFon() ) {
1347       FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp);
1348     }
1349
1350     delete[] isource;
1351   }
1352
1353   //the following values should go to AliTRDfeeParam once they are defined; then they have to be read in properly
1354   //naming follows conventions in TRAP-manual
1355
1356
1357   Bool_t bVBY = kTRUE;                         // cluster-verification bypass
1358
1359   Double_t cQTParam = 0;                      // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10
1360   AliTRDtrapAlu cQTAlu;
1361   cQTAlu.Init(1,10,0,63);
1362   cQTAlu.AssignDouble(cQTParam);
1363   Int_t cQT = cQTAlu.GetValue();
1364
1365   // linear fit
1366   Int_t tFS = fFeeParam->GetLinearFitStart();  // linear fit start
1367   Int_t tFE = fFeeParam->GetLinearFitEnd();    // linear fit stop
1368
1369   // charge accumulators
1370   Int_t tQS0 = fFeeParam->GetQacc0Start();     // start-time for charge-accumulator 0
1371   Int_t tQE0 = fFeeParam->GetQacc0End();       // stop-time for charge-accumulator 0
1372   Int_t tQS1 = fFeeParam->GetQacc1Start();     // start-time for charge-accumulator 1
1373   Int_t tQE1 = fFeeParam->GetQacc1End();       // stop-time for charge-accumulator 1
1374   // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N
1375
1376   Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold
1377   AliTRDtrapAlu cTHAlu;
1378   cTHAlu.Init(12,2);
1379   cTHAlu.AssignDouble(cTHParam);
1380   Int_t cTH = cTHAlu.GetValue();                                 // cTH used for comparison
1381
1382   struct List_t {
1383     List_t *next;
1384     Int_t iadc;
1385     Int_t value;
1386   };
1387
1388   List_t selection[7];            // list with 7 elements
1389   List_t *list = NULL;
1390   List_t *listLeft = NULL;
1391
1392   Int_t* qsum = new Int_t[fNADC];
1393
1394   // fit sums
1395   AliTRDtrapAlu qsumAlu;
1396   qsumAlu.Init(12,2);           // charge sum will be 12+2 bits
1397   AliTRDtrapAlu dCOGAlu;
1398   dCOGAlu.Init(1,7,0,127);      // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT
1399   AliTRDtrapAlu yrawAlu;
1400   yrawAlu.Init(1,8,-1,255);
1401   AliTRDtrapAlu yAlu;
1402   yAlu.Init(1,16,-1,0xFF00);    // only first 8 past-comma bits filled;additional 8 bits for accuracy;maximum 1 - 2^-8; sign is given by + or -
1403   AliTRDtrapAlu xAlu;
1404   xAlu.Init(5,8);               // 8 past-comma bits because value will be added/multiplied to another value with this accuracy
1405   AliTRDtrapAlu xxAlu;
1406   xxAlu.Init(10,0);
1407   AliTRDtrapAlu yyAlu;
1408   yyAlu.Init(1,16,0,0xFFFF);    // maximum is 2^16-1; 16Bit for past-commas
1409   AliTRDtrapAlu xyAlu;
1410   xyAlu.Init(6,8);
1411   AliTRDtrapAlu XAlu;
1412   XAlu.Init(9,0);
1413   AliTRDtrapAlu XXAlu;
1414   XXAlu.Init(14,0);
1415   AliTRDtrapAlu YAlu;
1416   YAlu.Init(5,8);               // 14 bit, 1 is sign-bit; therefore only 13 bit
1417   AliTRDtrapAlu YYAlu;
1418   YYAlu.Init(5,16);
1419   AliTRDtrapAlu XYAlu;
1420   XYAlu.Init(8,8);              // 17 bit, 1 is sign-bit; therefore only 16 bit
1421   AliTRDtrapAlu qtruncAlu;
1422   qtruncAlu.Init(12,0);
1423   AliTRDtrapAlu QT0Alu;
1424   QT0Alu.Init(15,0);
1425   AliTRDtrapAlu QT1Alu;
1426   QT1Alu.Init(16,0);
1427
1428   AliTRDtrapAlu oneAlu;
1429   oneAlu.Init(1,8);
1430
1431
1432   AliTRDtrapAlu inverseNAlu;
1433   inverseNAlu.Init(1,8);        // simulates the LUT for 1/N
1434   AliTRDtrapAlu MeanChargeAlu;  // mean charge in ADC counts
1435   MeanChargeAlu.Init(8,0);
1436   AliTRDtrapAlu TotalChargeAlu;
1437   TotalChargeAlu.Init(17,8);
1438   //nr of post comma bits should be the same for inverseN and TotalCharge
1439
1440
1441   SetPosLUT();                    // initialize the position correction LUT for this MCM;
1442
1443
1444   // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
1445   Int_t *X   = new Int_t[fNADC-2];
1446   Int_t *XX  = new Int_t[fNADC-2];
1447   Int_t *Y   = new Int_t[fNADC-2];
1448   Int_t *YY  = new Int_t[fNADC-2];
1449   Int_t *XY  = new Int_t[fNADC-2];
1450   Int_t *N   = new Int_t[fNADC-2];
1451   Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
1452   Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
1453
1454   for (Int_t iCol = 0; iCol < fNADC-2; iCol++) {
1455
1456       // initialize fit-sums
1457       X[iCol]   = 0;
1458       XX[iCol]  = 0;
1459       Y[iCol]   = 0;
1460       YY[iCol]  = 0;
1461       XY[iCol]  = 0;
1462       N[iCol]   = 0;
1463       QT0[iCol] = 0;
1464       QT1[iCol] = 0;
1465   }
1466
1467
1468   filPed.Init(7,2);                         // convert filtered pedestal into 7+2Bits
1469
1470   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1471
1472     if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no?
1473
1474     // reset
1475     Int_t portChannel[4]   = {-1,-1,-1,-1};
1476     Int_t clusterCharge[4] = {0,0,0,0};
1477     Int_t leftCharge[4]    = {0,0,0,0};
1478     Int_t centerCharge[4]  = {0,0,0,0};
1479     Int_t rightCharge[4]   = {0,0,0,0};
1480
1481     Int_t mark = 0;
1482
1483     filPed.AssignFormatted(ieffped[iT]);   // no size-checking when using AssignFormatted; ieffped>=0
1484     filPed = filPed;                       // this checks the size
1485
1486     ieffped[iT] = filPed.GetValue();
1487
1488     for(Int_t i = 0; i<7; i++){
1489       selection[i].next       = NULL;
1490       selection[i].iadc       =   -1;     // value of -1: invalid adc
1491       selection[i].value      =    0;
1492
1493     }
1494     // selection[0] is starting list-element; just for pointing
1495
1496     // loop over inner adc's
1497     for (Int_t iCol = 1; iCol < fNADC-1; iCol++) {
1498
1499       Int_t left   = fADCT[iCol-1][iT];
1500       Int_t center = fADCT[iCol][iT];
1501       Int_t right  = fADCT[iCol+1][iT];
1502
1503       Int_t sum = left + center + right;            // cluster charge sum
1504       qsumAlu.AssignFormatted(sum);
1505       qsumAlu = qsumAlu;                        // size-checking; redundant
1506
1507       qsum[iCol] = qsumAlu.GetValue();
1508
1509       //hit detection and masking
1510       if(center>=left){
1511         if(center>right){
1512           if(qsum[iCol]>=(cTH + 3*ieffped[iT])){    // effective pedestal of all three channels must be added to cTH(+20); this is not parallel to TRAP manual; maybe cTH has to be adjusted in fFeeParam; therefore channels are not yet reduced by their pedestal
1513             mark |= 1;                              // marker
1514           }
1515         }
1516       }
1517       mark = mark<<1;
1518     }
1519     mark = mark>>1;
1520
1521
1522     // get selection of 6 adc's and sort,starting with greatest values
1523
1524     //read three from right side and sort (primitive sorting algorithm)
1525     Int_t i = 0; // adc number
1526     Int_t j = 1; // selection number
1527     while(i<fNADC-2 && j<=3){
1528       i = i + 1;
1529       if( ((mark>>(i-1)) & 1) == 1) {
1530         selection[j].iadc  = fNADC-1-i;
1531         selection[j].value = qsum[fNADC-1-i]>>6;   // for hit-selection only the first 8 out of the 14 Bits are used for comparison
1532
1533         // insert into sorted list
1534         listLeft = &selection[0];
1535         list = listLeft->next;
1536
1537         if(list!=NULL) {
1538           while((list->next != NULL) && (selection[j].value <= list->value)){
1539             listLeft = list;
1540             list = list->next;
1541           }
1542
1543           if(selection[j].value<=list->value){
1544             selection[j].next = list->next;
1545             list->next = &selection[j];
1546           }
1547           else {
1548             listLeft->next = &selection[j];
1549             selection[j].next = list;
1550           }
1551         }
1552         else{
1553           listLeft->next = &selection[j];
1554           selection[j].next = list;
1555         }
1556
1557         j = j + 1;
1558       }
1559     }
1560
1561
1562     // read three from left side
1563     Int_t k = fNADC-2;
1564     while(k>i && j<=6) {
1565       if( ((mark>>(k-1)) & 1) == 1) {
1566         selection[j].iadc  = fNADC-1-k;
1567         selection[j].value = qsum[fNADC-1-k]>>6;
1568
1569         listLeft = &selection[0];
1570         list = listLeft->next;
1571
1572         if(list!=NULL){
1573           while((list->next != NULL) && (selection[j].value <= list->value)){
1574             listLeft = list;
1575             list = list->next;
1576           }
1577
1578           if(selection[j].value<=list->value){
1579             selection[j].next = list->next;
1580             list->next = &selection[j];
1581           }
1582           else {
1583             listLeft->next = &selection[j];
1584             selection[j].next = list;
1585           }
1586         }
1587         else{
1588           listLeft->next = &selection[j];
1589           selection[j].next = list;
1590         }
1591
1592         j = j + 1;
1593       }
1594       k = k - 1;
1595     }
1596
1597     // get the four with greatest charge-sum
1598     list = &selection[0];
1599     for(i = 0; i<4; i++){
1600       if(list->next == NULL) continue;
1601       list = list->next;
1602       if(list->iadc == -1) continue;
1603       Int_t adc = list->iadc;                              // channel number with selected hit
1604
1605       // the following arrays contain the four chosen channels in 1 time-bin
1606       portChannel[i]   = adc;
1607       clusterCharge[i] = qsum[adc];
1608       leftCharge[i]    = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
1609       centerCharge[i]  = fADCT[adc][iT] - ieffped[iT];
1610       rightCharge[i]   = fADCT[adc+1][iT] - ieffped[iT];
1611     }
1612
1613     // arithmetic unit
1614
1615     // cluster verification
1616     if(!bVBY){
1617       for(i = 0; i<4; i++){
1618         Int_t lr = leftCharge[i]*rightCharge[i]*1024;
1619         Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
1620         if (lr>=cc){
1621           portChannel[i]   = -1;                                 // set to invalid address
1622           clusterCharge[i] = 0;
1623         }
1624       }
1625     }
1626
1627     // fit-sums of valid channels
1628     // local hit position
1629     for(i = 0; i<4; i++){
1630       if (centerCharge[i] ==  0) {
1631         portChannel[i] = -1;
1632       }// prevent division by 0
1633
1634       if (portChannel[i]  == -1) continue;
1635
1636       Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
1637
1638       Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
1639       dCOG = (sign == kFALSE) ? dCOG : -dCOG;     // AssignDouble doesn't allow for signed doubles
1640       dCOGAlu.AssignDouble(dCOG);
1641       Int_t iLUTpos = dCOGAlu.GetValue();       // steers position in LUT
1642
1643       dCOG = dCOG/2;
1644       yrawAlu.AssignDouble(dCOG);
1645       Int_t iCOG = yrawAlu.GetValue();
1646       Int_t y = iCOG + fPosLUT[iLUTpos % 128];    // local position in pad-units
1647       yrawAlu.AssignFormatted(y);               // 0<y<1
1648       yAlu  = yrawAlu;                        // convert to 16 past-comma bits
1649
1650       if(sign == kTRUE) yAlu.SetSign(-1);       // buffer width of 9 bits; sign on real (not estimated) position
1651       xAlu.AssignInt(iT);                       // buffer width of 5 bits
1652
1653
1654       xxAlu = xAlu * xAlu;                  // buffer width of 10 bits -> fulfilled by x*x
1655
1656       yyAlu = yAlu * yAlu;                  // buffer width of 16 bits
1657
1658       xyAlu = xAlu * yAlu;                  // buffer width of 14 bits
1659
1660       Int_t adc = portChannel[i]-1;              // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
1661
1662       // calculate fit-sums recursively
1663       // interpretation of their bit-length is given as comment
1664
1665       // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
1666
1667       XAlu.AssignFormatted(X[adc]);
1668       XAlu = XAlu + xAlu;                   // buffer width of 9 bits
1669       X[adc] = XAlu.GetValue();
1670
1671       XXAlu.AssignFormatted(XX[adc]);
1672       XXAlu = XXAlu + xxAlu;                // buffer width of 14 bits
1673       XX[adc] = XXAlu.GetValue();
1674
1675       if (Y[adc] < 0) {
1676         YAlu.AssignFormatted(-Y[adc]);          // make sure that only positive values are assigned; sign-setting must be done by hand
1677         YAlu.SetSign(-1);
1678       }
1679       else {
1680         YAlu.AssignFormatted(Y[adc]);
1681         YAlu.SetSign(1);
1682       }
1683
1684       YAlu = YAlu + yAlu;                   // buffer width of 14 bits (8 past-comma);
1685       Y[adc] = YAlu.GetSignedValue();
1686
1687       YYAlu.AssignFormatted(YY[adc]);
1688       YYAlu = YYAlu + yyAlu;                // buffer width of 21 bits (16 past-comma)
1689       YY[adc] = YYAlu.GetValue();
1690
1691       if (XY[adc] < 0) {
1692         XYAlu.AssignFormatted(-XY[adc]);
1693         XYAlu.SetSign(-1);
1694       }
1695       else {
1696         XYAlu.AssignFormatted(XY[adc]);
1697         XYAlu.SetSign(1);
1698       }
1699
1700       XYAlu = XYAlu + xyAlu;                // buffer allows 17 bits (8 past-comma)
1701       XY[adc] = XYAlu.GetSignedValue();
1702
1703       N[adc]  = N[adc] + 1;
1704
1705
1706       // accumulated charge
1707       qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
1708       qtruncAlu = qsumAlu;
1709
1710       if(iT>=tQS0 && iT<=tQE0){
1711         QT0Alu.AssignFormatted(QT0[adc]);
1712         QT0Alu = QT0Alu + qtruncAlu;
1713         QT0[adc] = QT0Alu.GetValue();
1714         //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
1715       }
1716
1717       if(iT>=tQS1 && iT<=tQE1){
1718         QT1Alu.AssignFormatted(QT1[adc]);
1719         QT1Alu = QT1Alu + qtruncAlu;
1720         QT1[adc] = QT1Alu.GetValue();
1721         //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
1722       }
1723     }// i
1724
1725     // remapping is done!!
1726
1727   }//iT
1728
1729
1730
1731   // tracklet-assembly
1732
1733   // put into AliTRDfeeParam and take care that values are in proper range
1734   const Int_t cTCL = 1;      // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??)
1735   const Int_t cTCT = 8;      // joint number of hits;     8<=TCT<=31; note that according to TRAP manual this number cannot be lower than 8; however it should be adjustable to the number of hits in the fit time range (40%)
1736
1737   Int_t mPair   = 0;         // marker for possible tracklet pairs
1738   Int_t* hitSum = new Int_t[fNADC-3];
1739   // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18;
1740
1741   // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
1742   for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
1743     hitSum[iCol] = N[iCol] + N[iCol+1];
1744     if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
1745         mPair |= 1;         // mark as possible channel-pair
1746
1747     }
1748     mPair = mPair<<1;
1749   }
1750   mPair = mPair>>1;
1751
1752   List_t* selectPair = new List_t[fNADC-2];      // list with 18 elements (0..18) containing the left channel-nr and hit sums
1753                                                  // selectPair[18] is starting list-element just for pointing
1754   for(Int_t k = 0; k<fNADC-2; k++){
1755       selectPair[k].next       = NULL;
1756       selectPair[k].iadc       =   -1;           // invalid adc
1757       selectPair[k].value      =    0;
1758
1759     }
1760
1761  list = NULL;
1762  listLeft = NULL;
1763
1764   // read marker and sort according to hit-sum
1765
1766   Int_t adcL  = 0;            // left adc-channel-number (remapped)
1767   Int_t selNr = 0;            // current number in list
1768
1769   // insert marked channels into list and sort according to hit-sum
1770   while(adcL < fNADC-3 && selNr < fNADC-3){
1771
1772     if( ((mPair>>((fNADC-4)-(adcL))) & 1) == 1) {
1773       selectPair[selNr].iadc  = adcL;
1774       selectPair[selNr].value = hitSum[adcL];
1775
1776       listLeft = &selectPair[fNADC-3];
1777       list = listLeft->next;
1778
1779       if(list!=NULL) {
1780         while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
1781           listLeft = list;
1782           list = list->next;
1783         }
1784
1785         if(selectPair[selNr].value <= list->value){
1786           selectPair[selNr].next = list->next;
1787           list->next = &selectPair[selNr];
1788         }
1789         else {
1790           listLeft->next = &selectPair[selNr];
1791           selectPair[selNr].next = list;
1792         }
1793
1794       }
1795       else{
1796         listLeft->next = &selectPair[selNr];
1797         selectPair[selNr].next = list;
1798       }
1799
1800       selNr = selNr + 1;
1801     }
1802     adcL = adcL + 1;
1803   }
1804
1805   //select up to 4 channels with maximum number of hits
1806   Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
1807   Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
1808   list = &selectPair[fNADC-3];
1809
1810   for (Int_t i = 0; i<4; i++) {
1811     if(list->next == NULL) continue;
1812     list = list->next;
1813     if(list->iadc == -1) continue;
1814     lpairChannel[i] = list->iadc;        // channel number with selected hit
1815     rpairChannel[i] = lpairChannel[i]+1;
1816   }
1817
1818   // avoid submission of double tracklets
1819   for (Int_t i = 3; i>0; i--) {
1820     for (Int_t j = i-1; j>-1; j--) {
1821       if(lpairChannel[i] == rpairChannel[j]) {
1822         lpairChannel[i] = -1;
1823         rpairChannel[i] = -1;
1824         break;
1825       }
1826       /* if(rpairChannel[i] == lpairChannel[j]) {
1827         lpairChannel[i] = -1;
1828         rpairChannel[i] = -1;
1829         break;
1830         }*/
1831     }
1832   }
1833
1834   // merging of the fit-sums of the remainig channels
1835   // assume same data-word-width as for fit-sums for 1 channel
1836   // relative scales!
1837   Int_t mADC[4];
1838   Int_t mN[4];
1839   Int_t mQT0[4];
1840   Int_t mQT1[4];
1841   Int_t mX[4];
1842   Int_t mXX[4];
1843   Int_t mY[4];
1844   Int_t mYY[4];
1845   Int_t mXY[4];
1846   Int_t mOffset[4];
1847   Int_t mSlope[4];
1848   Int_t mMeanCharge[4];
1849   Int_t inverseN = 0;
1850   Double_t invN = 0;
1851   Int_t one = 0;
1852
1853   for (Int_t i = 0; i<4; i++){
1854     mADC[i] = -1;                        // set to invalid number
1855     mN[i]   =  0;
1856     mQT0[i] =  0;
1857     mQT1[i] =  0;
1858     mX[i]   =  0;
1859     mXX[i]  =  0;
1860     mY[i]   =  0;
1861     mYY[i]  =  0;
1862     mXY[i]  =  0;
1863     mOffset[i] = 0;
1864     mSlope[i]  = 0;
1865     mMeanCharge[i] = 0;
1866   }
1867
1868   oneAlu.AssignInt(1);
1869   one = oneAlu.GetValue();              // one with 8 past comma bits
1870
1871   for (Int_t i = 0; i<4; i++){
1872
1873
1874     mADC[i] = lpairChannel[i];          // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
1875                                         // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
1876     Int_t madc = mADC[i];
1877     if (madc == -1) continue;
1878
1879     YAlu.AssignInt(N[rpairChannel[i]]);
1880     Int_t wpad  = YAlu.GetValue();       // enlarge hit counter of right channel by 8 past-comma bits; YAlu can have 5 pre-comma bits (values up to 63); hit counter<=nr of time bins (24)
1881
1882     mN[i]    = hitSum[madc];
1883
1884     // don't merge fit sums in case of a stand-alone tracklet (consisting of only 1 channel); in that case only left channel makes up the fit sums
1885     if (N[madc+1] == 0) {
1886         mQT0[i] = QT0[madc];
1887         mQT1[i] = QT1[madc];
1888
1889     }
1890     else {
1891
1892         // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
1893
1894         mQT0[i]   = QT0[madc] + QT0[madc+1];
1895         QT0Alu.AssignFormatted(mQT0[i]);
1896         QT0Alu  = QT0Alu;                // size-check
1897         mQT0[i]   = QT0Alu.GetValue();     // write back
1898
1899         mQT1[i]   = QT1[madc] + QT1[madc+1];
1900         QT1Alu.AssignFormatted(mQT1[i]);
1901         QT1Alu  = QT1Alu;
1902         mQT1[i]   = QT1Alu.GetValue();
1903     }
1904
1905     // calculate the mean charge in adc values; later to be replaced by electron likelihood
1906     mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
1907     mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
1908     // simulate LUT for 1/N; LUT is fed with the double-accurate pre-calculated value of 1/N; accuracy of entries has to be adjusted to real TRAP
1909     invN = 1.0/(mN[i]);
1910     inverseNAlu.AssignDouble(invN);
1911     inverseN = inverseNAlu.GetValue();
1912     mMeanCharge[i] = mMeanCharge[i] * inverseN;  // now to be interpreted with 8 past-comma bits
1913     TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
1914     TotalChargeAlu = TotalChargeAlu;
1915     MeanChargeAlu = TotalChargeAlu;
1916     mMeanCharge[i] = MeanChargeAlu.GetValue();
1917
1918     // this check is not necessary; it is just for efficiency reasons
1919     if (N[madc+1] == 0) {
1920         mX[i]     =   X[madc];
1921         mXX[i]    =  XX[madc];
1922         mY[i]     =   Y[madc];
1923         mXY[i]    =  XY[madc];
1924         mYY[i]    =  YY[madc];
1925     }
1926     else {
1927
1928         mX[i]     =   X[madc] +  X[madc+1];
1929         XAlu.AssignFormatted(mX[i]);
1930         XAlu      = XAlu;
1931         mX[i]     = XAlu.GetValue();
1932
1933         mXX[i]    =  XX[madc] + XX[madc+1];
1934         XXAlu.AssignFormatted(mXX[i]);
1935         XXAlu     = XXAlu;
1936         mXX[i]    = XXAlu.GetValue();
1937
1938
1939         mY[i]     =   Y[madc] + Y[madc+1] + wpad;
1940         if (mY[i] < 0) {
1941             YAlu.AssignFormatted(-mY[i]);
1942             YAlu.SetSign(-1);
1943         }
1944         else {
1945             YAlu.AssignFormatted(mY[i]);
1946             YAlu.SetSign(1);
1947         }
1948         YAlu    = YAlu;
1949         mY[i]     = YAlu.GetSignedValue();
1950
1951         mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*one;    // multiplication by one to maintain the data format
1952
1953         if (mXY[i] < 0) {
1954             XYAlu.AssignFormatted(-mXY[i]);
1955             XYAlu.SetSign(-1);
1956         }
1957         else {
1958             XYAlu.AssignFormatted(mXY[i]);
1959             XYAlu.SetSign(1);
1960         }
1961         XYAlu   = XYAlu;
1962         mXY[i]    = XYAlu.GetSignedValue();
1963
1964         mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
1965         if (mYY[i] < 0) {
1966             YYAlu.AssignFormatted(-mYY[i]);
1967             YYAlu.SetSign(-1);
1968         }
1969         else {
1970             YYAlu.AssignFormatted(mYY[i]);
1971             YYAlu.SetSign(1);
1972         }
1973
1974         YYAlu   = YYAlu;
1975         mYY[i]    = YYAlu.GetSignedValue();
1976     }
1977
1978   }
1979
1980   // calculation of offset and slope from the merged fit-sums;
1981   // YY is needed for some error measure only; still to be done
1982   // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
1983
1984   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1985   // !!important note: the offset is calculated from hits in the time bin range between tFS and tFE; it corresponds to the value at the height of the time bin tFS which does NOT need to correspond to the upper side of the drift   !!
1986   // !!volume (cathode wire plane). The offset cannot be rescaled as long as it is unknown which is the first time bin that contains hits from the drift region and thus to which distance from the cathode plane tFS corresponds.    !!
1987   // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below).                                                                                             !!
1988   // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used.                                                                                          !!
1989   // !!The corrections to the offset (e.g. no ExB correction applied as offset is supposed to be on top of drift region; however not at anode wire, so some inclination of drifting clusters due to Lorentz angle exists) are only    !!
1990   // !!valid (in approximation) if tFS is close to the beginning of the drift region.                                                                                                                                                 !!
1991   // !!The slope however can be converted to a deflection length between electrode and cathode wire plane as it is clear that the drift region is sampled 20 times                                                                    !!
1992   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1993
1994   // which formats should be chosen?
1995   AliTRDtrapAlu denomAlu;
1996   denomAlu.Init(20,8);
1997   AliTRDtrapAlu numAlu;
1998   numAlu.Init(20,8);
1999   // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
2000   // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
2001
2002   for (Int_t i = 0; i<4; i++) {
2003     if (mADC[i] == -1) continue;
2004
2005     Int_t num0  = (mN[i]*mXX[i]-mX[i]*mX[i]);
2006     if (num0 < 0) {
2007       denomAlu.AssignInt(-num0);    // num0 does not have to be interpreted as having past-comma bits -> AssignInt
2008       denomAlu.SetSign(-1);
2009     }
2010     else {
2011       denomAlu.AssignInt(num0);
2012       denomAlu.SetSign(1);
2013     }
2014
2015     Int_t num1  = mN[i]*mXY[i] - mX[i]*mY[i];
2016     if (num1 < 0) {
2017       numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
2018       numAlu.SetSign(-1);
2019     }
2020     else {
2021       numAlu.AssignFormatted(num1);
2022       numAlu.SetSign(1);
2023     }
2024     numAlu    = numAlu/denomAlu;
2025     mSlope[i]   = numAlu.GetSignedValue();
2026
2027     Int_t num2  = mXX[i]*mY[i] - mX[i]*mXY[i];
2028
2029     if (num2 < 0) {
2030       numAlu.AssignFormatted(-num2);
2031       numAlu.SetSign(-1);
2032     }
2033     else {
2034       numAlu.AssignFormatted(num2);
2035       numAlu.SetSign(1);
2036     }
2037
2038     numAlu    = numAlu/denomAlu;
2039
2040
2041     mOffset[i]  = numAlu.GetSignedValue();
2042     numAlu.SetSign(1);
2043     denomAlu.SetSign(1);
2044
2045
2046     //numAlu.AssignInt(mADC[i]+1);   // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)
2047     numAlu.AssignDouble((Double_t)mADC[i] + 1.5);      // numAlu has enough pre-comma place for that; correct trafo, best values
2048     mOffset[i]  = mOffset[i] + numAlu.GetValue();      // transform offset to a coord.system relative to chip; +1 to avoid neg. values
2049
2050     // up to here: adc-mapping according to TRAP manual and in line with pad-col mapping
2051     // reverse adc-counting to be again in line with the online mapping
2052     mADC[i]     = fNADC - 4 - mADC[i];                 // fNADC-4-mADC[i]: 0..17; remapping necessary;
2053     mADC[i]     = mADC[i] + 2;
2054     // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
2055   }
2056
2057   // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected;
2058   // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
2059
2060   // transform parameters to the local coordinate-system of a stack (used by GTU)
2061   AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
2062
2063   Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
2064   //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
2065
2066   // difference between width of inner and outer pads of a row is not accounted for;
2067
2068   Double_t magField = 0.4;                           // z-component of magnetic field in Tesla; adjust to current simulation!!; magnetic field can hardly be evaluated for the position of each mcm
2069   Double_t eCharge  = 0.3;                           // unit charge in (GeV/c)/m*T
2070   Double_t ptMin   = 2.3;                            // minimum transverse momentum (GeV/c); to be adjusted(?)
2071
2072   Double_t granularityOffset = 0.160;                // granularity for offset in mm
2073   Double_t granularitySlope  = 0.140;                // granularity for slope  in mm
2074
2075   // get the coordinates in SM-system; parameters:
2076
2077   Double_t zPos       =  (padPlane->GetRowPos(fRow))*10.0;  // z-position of the MCM; fRow is counted on a chamber; SM consists of 5
2078   // zPos is position of pad-borders;
2079   Double_t zOffset = 0.0;
2080   if ( fRow == 0 || fRow == 15 ) {
2081       zOffset = padPlane->GetLengthOPad();
2082   }
2083   else {
2084       zOffset = padPlane->GetLengthIPad();
2085   }
2086   zOffset = (-1.0) * zOffset/2.0;
2087   // turn zPos to be z-coordinate at middle of pad-row
2088   zPos = zPos + zOffset;
2089
2090
2091   Double_t xPos       =  0.0;                               // x-position of the upper border of the drift-chamber of actual layer
2092   Int_t    icol       =  0;                                 // column-number of adc-channel
2093   Double_t yPos[4];                                         // y-position of the pad to which ADC is connected
2094   Double_t dx         = 30.0;                               // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
2095   Double_t freqSample = fFeeParam->GetSamplingFrequency();  // retrieve the sampling frequency (10.019750 MHz)
2096   Double_t vdrift     = fCal->GetVdriftAverage(fChaId);     // averaged drift velocity for this detector (1.500000 cm/us)
2097   Int_t    nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
2098   Int_t    nrOfAmplTimeBins  = 2;                           // the number of time bins between anode wire and cathode wires in ampl.region (3.5mm)(guess)(suppose v_drift+3.5cm/us there=>all clusters arrive at anode wire within one time bin (100ns))
2099   Int_t    nrOfOffsetCorrTimeBins = tFS - nrOfAmplTimeBins - 1; // -1 is  to be conservative; offset correction will not remove the shift but is supposed to improve it; if tFS = 5, 2 drift time bins before tFS are assumed
2100   if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed
2101   Double_t lorTan     = fCal->GetOmegaTau(vdrift,magField); // tan of the Lorentz-angle for this detector; could be evaluated and set as a parameter for each mcm
2102   //Double_t lorAngle   =  7.0;                             // Lorentz-angle in degrees
2103   Double_t tiltAngle  = padPlane->GetTiltingAngle();        // sign-respecting tilting angle of pads in actual layer
2104   Double_t tiltTan    = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
2105   //Double_t lorTan     = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
2106
2107   Double_t alphaMax[4];                            // maximum deflection from the direction to the primary vertex; granularity of hit pads
2108   Double_t slopeMin[4];                            // local limits for the deflection
2109   Double_t slopeMax[4];
2110   Int_t   mslopeMin[4];                            // in granularity units; to be compared to mSlope[i]
2111   Int_t   mslopeMax[4];
2112
2113
2114   // x coord. of upper side of drift chambers in local SM-system (in mm)
2115   // obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
2116   // the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
2117   switch(fLayer)
2118     {
2119     case 0:
2120       xPos = 3003.0;
2121       break;
2122     case 1:
2123       xPos = 3129.0;
2124       break;
2125     case 2:
2126       xPos = 3255.0;
2127       break;
2128     case 3:
2129       xPos = 3381.0;
2130       break;
2131     case 4:
2132       xPos = 3507.0;
2133       break;
2134     case 5:
2135       xPos = 3633.0;
2136       break;
2137     }
2138
2139   // calculation of offset-correction n:
2140
2141   Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));
2142
2143   nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
2144   if (nCorrectOffset < 0) {
2145     numAlu.AssignInt(-nCorrectOffset);
2146     numAlu.SetSign(-1);
2147   }
2148   else {
2149     numAlu.AssignInt(nCorrectOffset);
2150     numAlu.SetSign(1);
2151   }
2152   nCorrectOffset = numAlu.GetSignedValue();
2153
2154   // the Lorentz correction to the offset
2155   Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
2156
2157
2158   lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
2159
2160   if(lorCorrectOffset < 0) {
2161       numAlu.AssignDouble(-lorCorrectOffset);
2162       numAlu.SetSign(-1);
2163   }
2164   else{
2165       numAlu.AssignDouble(lorCorrectOffset);
2166       numAlu.SetSign(1);
2167   }
2168
2169   Int_t mlorCorrectOffset = numAlu.GetSignedValue();
2170
2171
2172   Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
2173
2174   // calculation of slope-correction
2175
2176   // this is only true for tracks coming (approx.) from primary vertex
2177   // everything is evaluated for a tracklet covering the whole drift chamber
2178   Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
2179   // Double_t cCorrectSlope =  zPos/xPos*dx*tiltTan/granularitySlope;
2180   // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
2181
2182   if (cCorrectSlope < 0) {
2183       numAlu.AssignDouble(-cCorrectSlope);
2184       numAlu.SetSign(-1);
2185   }
2186   else {
2187       numAlu.AssignDouble(cCorrectSlope);
2188       numAlu.SetSign(1);
2189   }
2190   cCorrectSlope = numAlu.GetSignedValue();
2191
2192   // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
2193   // different pad-width of outer pads of a pad-plane not taken into account
2194   // note that the fit was only done in the range tFS to tFE, however this range does not need to cover the whole drift region (neither start nor end of it)
2195   // however the tracklets are supposed to be a fit in the drift region thus the linear function is stretched to fit the drift region of 30 mm
2196
2197
2198   Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope;  // >= 0.0
2199
2200   AliTRDtrapAlu correctAlu;
2201   correctAlu.Init(20,8);
2202
2203   AliTRDtrapAlu offsetAlu;
2204   offsetAlu.Init(13,0,-0x1000,0x0FFF);          // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
2205
2206   AliTRDtrapAlu slopeAlu;
2207   slopeAlu.Init(7,0,-0x40,0x3F);                // 7 bit-word;  2-complement (1 sign-bit);
2208
2209   for (Int_t i = 0; i<4; i++) {
2210
2211     if (mADC[i] == -1) continue;
2212
2213     icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
2214     yPos[i]   = (padPlane->GetColPos(icol))*10.0;
2215
2216
2217     // offset:
2218
2219     correctAlu.AssignDouble(mCorrectOffset);     // done because max. accuracy is 8 bit
2220     mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
2221     mOffset[i]  = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset));
2222     //mOffset[i]  = mOffset[i]*(-1);                   // adjust to direction of y-axes in online simulation
2223
2224     if (mOffset[i] < 0) {
2225       numAlu.AssignFormatted(-mOffset[i]);
2226       numAlu.SetSign(-1);
2227     }
2228     else {
2229       numAlu.AssignFormatted(mOffset[i]);
2230       numAlu.SetSign(1);
2231     }
2232
2233     offsetAlu = numAlu;
2234     mOffset[i]  = offsetAlu.GetSignedValue();
2235
2236
2237     // slope:
2238
2239     correctAlu.AssignDouble(mCorrectSlope);
2240     mCorrectSlope = correctAlu.GetValueWhole();
2241
2242     mSlope[i]   = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
2243
2244     if (mSlope[i] < 0) {
2245       numAlu.AssignFormatted(-mSlope[i]);
2246       numAlu.SetSign(-1);
2247     }
2248     else {
2249       numAlu.AssignFormatted(mSlope[i]);
2250       numAlu.SetSign(1);
2251     }
2252
2253     slopeAlu  = numAlu;     // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
2254     mSlope[i]   = slopeAlu.GetSignedValue();
2255
2256     // local (LTU) limits for the deflection
2257     // ATan returns angles in radian
2258     alphaMax[i]  = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
2259     slopeMin[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
2260     slopeMax[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
2261
2262     if (slopeMin[i] < 0) {
2263       slopeAlu.AssignDouble(-slopeMin[i]);
2264       slopeAlu.SetSign(-1);
2265     }
2266     else {
2267       slopeAlu.AssignDouble(slopeMin[i]);
2268       slopeAlu.SetSign(1);
2269     }
2270     mslopeMin[i] = slopeAlu.GetSignedValue();  // the borders should lie inside the range of mSlope -> usage of slopeAlu again
2271
2272     if (slopeMax[i] < 0) {
2273       slopeAlu.AssignDouble(-slopeMax[i]);
2274       slopeAlu.SetSign(-1);
2275     }
2276     else {
2277       slopeAlu.AssignDouble(slopeMax[i]);
2278       slopeAlu.SetSign(1);
2279     }
2280     mslopeMax[i] = slopeAlu.GetSignedValue();
2281   }
2282
2283   // suppress submission of tracks with low stiffness
2284   // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber)
2285
2286   // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
2287   // up to now they are sorted according to maximum hit sum
2288   // is the sorting really done in the TRAP-chip?
2289
2290   Int_t order[4] = {-1,-1,-1,-1};
2291   Int_t wordnr = 0;   // number of tracklet-words
2292
2293   for(Int_t j = 0; j < fMaxTracklets; j++) {
2294       //if( mADC[j] == -1) continue;
2295       if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
2296       wordnr++;
2297       if( wordnr-1 == 0) {
2298           order[0] = j;
2299           continue;
2300       }
2301       // wordnr-1>0, wordnr-1<4
2302       order[wordnr-1] = j;
2303       for( Int_t k = 0; k < wordnr-1; k++) {
2304           if( mOffset[j] < mOffset[order[k]] ) {
2305               for( Int_t l = wordnr-1; l > k; l-- ) {
2306                   order[l] = order[l-1];
2307               }
2308               order[k] = j;
2309               break;
2310           }
2311
2312       }
2313   }
2314
2315   // fill the bit-words in ascending order and without gaps
2316   UInt_t bitWord[4] = {0,0,0,0};                 // attention: unsigned int to have real 32 bits (no 2-complement)
2317   for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
2318       //Bool_t rem1 = kTRUE;
2319
2320     Int_t i = order[j];
2321     //bit-word is 2-complement and therefore without sign
2322     bitWord[j] =   1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
2323     //printf("\n");
2324     UInt_t shift  = 0;
2325     UInt_t shift2 = 0;
2326
2327
2328
2329
2330     /*printf("mean charge: %d\n",mMeanCharge[i]);
2331     printf("row: %d\n",fRow);
2332     printf("slope: %d\n",mSlope[i]);
2333     printf("pad position: %d\n",mOffset[i]);
2334     printf("channel: %d\n",mADC[i]);*/
2335
2336     // electron probability (currently not implemented; the mean charge is just scaled)
2337     shift = (UInt_t)mMeanCharge[i];
2338     for(Int_t iBit = 0; iBit < 8; iBit++) {
2339       bitWord[j]  = bitWord[j]<<1;
2340       bitWord[j] |= (shift>>(7-iBit))&1;
2341       //printf("0");
2342     }
2343
2344     // pad row
2345     shift = (UInt_t)fRow;
2346     for(Int_t iBit = 0; iBit < 4; iBit++) {
2347       bitWord[j]  = bitWord[j]<<1;
2348       bitWord[j] |= (shift>>(3-iBit))&1;
2349       //printf("%d", (fRow>>(3-iBit))&1);
2350     }
2351
2352     // deflection length
2353     if(mSlope[i] < 0) {
2354         shift = (UInt_t)(-mSlope[i]);
2355         // shift2 is 2-complement of shift
2356         shift2 = 1;
2357         for(Int_t iBit = 1; iBit < 7; iBit++) {
2358             shift2  = shift2<<1;
2359             shift2 |= (1- (((shift)>>(6-iBit))&1) );
2360             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2361         }
2362         shift2 = shift2 + 1;
2363         //printf("1");
2364         for(Int_t iBit = 0; iBit < 7; iBit++) {
2365             bitWord[j]  = bitWord[j]<<1;
2366             bitWord[j] |= (shift2>>(6-iBit))&1;
2367             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2368         }
2369     }
2370     else {
2371         shift = (UInt_t)(mSlope[i]);
2372         bitWord[j]  = bitWord[j]<<1;
2373         bitWord[j]   |= 0;
2374         //printf("0");
2375         for(Int_t iBit = 1; iBit < 7; iBit++) {
2376             bitWord[j]  = bitWord[j]<<1;
2377             bitWord[j] |= (shift>>(6-iBit))&1;
2378             //printf("%d",(mSlope[i]>>(6-iBit))&1);
2379         }
2380     }
2381
2382     // pad position
2383     if(mOffset[i] < 0) {
2384         shift = (UInt_t)(-mOffset[i]);
2385         shift2 = 1;
2386         for(Int_t iBit = 1; iBit < 13; iBit++) {
2387             shift2  = shift2<<1;
2388             shift2 |= (1-(((shift)>>(12-iBit))&1));
2389             //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
2390         }
2391         shift2 = shift2 + 1;
2392         //printf("1");
2393         for(Int_t iBit = 0; iBit < 13; iBit++) {
2394             bitWord[j]  = bitWord[j]<<1;
2395             bitWord[j] |= (shift2>>(12-iBit))&1;
2396             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2397         }
2398     }
2399     else {
2400         shift = (UInt_t)mOffset[i];
2401         bitWord[j] = bitWord[j]<<1;
2402         bitWord[j]   |= 0;
2403         //printf("0");
2404         for(Int_t iBit = 1; iBit < 13; iBit++) {
2405             bitWord[j]  = bitWord[j]<<1;
2406             bitWord[j] |= (shift>>(12-iBit))&1;
2407             //printf("%d",(mOffset[i]>>(12-iBit))&1);
2408         }
2409     }
2410
2411
2412
2413     //printf("bitWord: %u\n",bitWord[j]);
2414     //printf("adc: %d\n",mADC[i]);
2415     fMCMT[j] = bitWord[j];
2416   }
2417
2418   //printf("\n");
2419
2420
2421   delete [] qsum;
2422   delete [] ieffped;
2423
2424   delete [] X;
2425   delete [] XX;
2426   delete [] Y;
2427   delete [] YY;
2428   delete [] XY;
2429   delete [] N;
2430   delete [] QT0;
2431   delete [] QT1;
2432
2433   delete [] hitSum;
2434   delete [] selectPair;
2435
2436   delete padPlane;
2437
2438 //if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
2439
2440   if (!fFeeParam->GetMCTrackletOutput())
2441       return;
2442
2443   AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
2444   AliLog::SetFileOutput("../log/tracklet.log");
2445
2446   // testing for wordnr in order to speed up the simulation
2447   if (wordnr == 0)
2448     return;
2449
2450   UInt_t        *trackletWord = new UInt_t[fMaxTracklets];
2451   Int_t         *adcChannel   = new Int_t[fMaxTracklets];
2452   Int_t         *trackRef     = new Int_t[fMaxTracklets];
2453
2454   Int_t u = 0;
2455
2456   AliTRDdigitsManager *digman = new AliTRDdigitsManager();
2457   digman->ReadDigits(gAlice->GetRunLoader()->GetLoader("TRDLoader")->TreeD());
2458   digman->SetUseDictionaries(kTRUE);
2459   AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
2460
2461   for (Int_t j = 0; j < fMaxTracklets; j++) {
2462       Int_t i = order[j];
2463       trackletWord[j] = 0;
2464       adcChannel[j] = -1;
2465       if (bitWord[j]!=0) {
2466           trackletWord[u] = bitWord[j];
2467           adcChannel[u]   = mADC[i];   // mapping onto the original adc-array to be in line with the digits-adc-ordering (21 channels in total on 1 mcm, 18 belonging to pads); mADC[i] should be >-1 in case bitWord[i]>0
2468
2469 // Finding label of MC track
2470           TH1F *hTrkRef = new TH1F("trackref", "trackref", 100000, 0, 100000);
2471           Int_t track[3];
2472           Int_t padcol = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u]);
2473           Int_t padcol_ngb = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u] - 1);
2474           Int_t padrow = 4 * (fRobPos / 2) + fMcmPos / 4;
2475           Int_t det = 30 * fSector + 6 * fStack + fLayer;
2476           for(Int_t iTimebin = feeParam->GetLinearFitStart(); iTimebin < feeParam->GetLinearFitEnd(); iTimebin++) {
2477               track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2478               track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2479               track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2480               hTrkRef->Fill(track[0]);
2481               if (track[1] != track[0] && track[1] != -1)
2482                   hTrkRef->Fill(track[1]);
2483               if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2484                   hTrkRef->Fill(track[2]);
2485               if (padcol_ngb >= 0) {
2486                   track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2487                   track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2488                   track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2489                   hTrkRef->Fill(track[0]);
2490                   if (track[1] != track[0] && track[1] != -1)
2491                       hTrkRef->Fill(track[1]);
2492                   if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2493                       hTrkRef->Fill(track[2]);
2494               }
2495           }
2496           trackRef[u] = hTrkRef->GetMaximumBin() - 1;
2497           delete hTrkRef;
2498           u = u + 1;
2499       }
2500   }
2501
2502  AliDataLoader *dl = gAlice->GetRunLoader()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
2503   if (!dl) {
2504     AliError("Could not get the tracklets data loader!");
2505   }
2506   else {
2507     TTree *trackletTree = dl->Tree();
2508     if (!trackletTree)
2509       dl->MakeTree();
2510     trackletTree = dl->Tree();
2511
2512    AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM();
2513    TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
2514    if (!trkbranch)
2515        trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
2516     trkbranch->SetAddress(&trkl);
2517
2518     for (Int_t iTracklet = 0; iTracklet < fMaxTracklets; iTracklet++) {
2519         if (trackletWord[iTracklet] == 0)
2520             continue;
2521         trkl->SetTrackletWord(trackletWord[iTracklet]);
2522         trkl->SetDetector(30*fSector + 6*fStack + fLayer);
2523         trkl->SetROB(fRobPos);
2524         trkl->SetMCM(fMcmPos);
2525         trkl->SetLabel(trackRef[iTracklet]);
2526         trackletTree->Fill();
2527     }
2528     delete trkl;
2529     dl->WriteData("OVERWRITE");
2530   }
2531
2532   delete [] trackletWord;
2533   delete [] adcChannel;
2534   delete [] trackRef;
2535   delete digman;
2536
2537   // to be done:
2538   // error measure for quality of fit (not necessarily needed for the trigger)
2539   // cluster quality threshold (not yet set)
2540   // electron probability
2541 }
2542