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 #include "AliTRDarrayADC.h"
 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::ProduceRawStreamV2( UInt_t *buf, Int_t maxSize )
 614 {
 615   //
 616   // Produce raw 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   UInt_t  iEv = 0;
 623   Int_t   nw  = 0;  // Number of written words
 624   Int_t   of  = 0;  // Number of overflowed words
 625   Int_t   rawVer   = fFeeParam->GetRAWversion();
 626   Int_t **adc;
 627   Int_t   nActiveADC = 0;       // number of activated ADC bits in a word
 628
 629   if( !CheckInitialized() ) return 0;
 630
 631   if( fFeeParam->GetRAWstoreRaw() ) {
 632     adc = fADCR;
 633   } else {
 634     adc = fADCF;
 635   }
 636
 637   // Produce MCM header
 638   x = (1<<31) | ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
 639   if (nw < maxSize) {
 640     buf[nw++] = x;
 641         //printf("\nMCM header: %X ",x);
 642   }
 643   else {
 644     of++;
 645   }
 646
 647   // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
 648   //                            n : unused , c : ADC count, m : selected ADCs
 649   if( rawVer >= 3 ) {
 650     x = 0;
 651     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
 652       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
 653                 x = x | (1 << (iAdc+4) );       // last 4 digit reserved for 1100=0xc
 654                 nActiveADC++;           // number of 1 in mmm....m
 655       }
 656     }
 657         x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC;   // nn = 01, ccccc are inverted, 0xc=1100
 658         //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
 659
 660     if (nw < maxSize) {
 661       buf[nw++] = x;
 662           //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
 663     }
 664     else {
 665       of++;
 666     }
 667   }
 668
 669   // Produce ADC data. 3 timebins are packed into one 32 bits word
 670   // In this version, different ADC channel will NOT share the same word
 671
 672   UInt_t aa=0, a1=0, a2=0, a3=0;
 673
 674   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
 675     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
 676     aa = !(iAdc & 1) + 2;
 677     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
 678       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] : 0;
 679       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
 680       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
 681       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
 682       if (nw < maxSize) {
 683         buf[nw++] = x;
 684         //printf("%08X ",x);
 685       }
 686       else {
 687         of++;
 688       }
 689     }
 690   }
 691
 692   if( of != 0 ) return -of; else return nw;
 693 }
 694
 695 //_____________________________________________________________________________
 696 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
 697 {
 698   //
 699   // Produce tracklet data stream from this MCM and put in buf
 700   // Returns number of words filled, or negative value
 701   // with -1 * number of overflowed words
 702   //
 703
 704   UInt_t  x;
 705   Int_t   nw  = 0;  // Number of written words
 706   Int_t   of  = 0;  // Number of overflowed words
 707
 708   if( !CheckInitialized() ) return 0;
 709
 710   // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM
 711   // fMCMT is filled continuously until no more tracklet words available
 712
 713   Int_t wd = 0;
 714   while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){
 715       x = fMCMT[wd];
 716       if (nw < maxSize) {
 717         buf[nw++] = x;
 718       }
 719       else {
 720         of++;
 721       }
 722       wd++;
 723   }
 724
 725   if( of != 0 ) return -of; else return nw;
 726 }
 727
 728
 729 //_____________________________________________________________________________
 730 void AliTRDmcmSim::Filter()
 731 {
 732   //
 733   // Apply digital filter
 734   //
 735
 736   if( !CheckInitialized() ) return;
 737
 738   // Initialize filtered data array with raw data
 739   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 740     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 741       fADCF[iadc][it] = fADCR[iadc][it];
 742     }
 743   }
 744
 745   // Then apply fileters one by one to filtered data array
 746   if( fFeeParam->IsPFon() ) FilterPedestal();
 747   if( fFeeParam->IsGFon() ) FilterGain();
 748   if( fFeeParam->IsTFon() ) FilterTail();
 749 }
 750
 751 //_____________________________________________________________________________
 752 void AliTRDmcmSim::FilterPedestal()
 753 {
 754
 755
 756   //
 757   // Apply pedestal filter
 758   //
 759
 760   Int_t ap = fSimParam->GetADCbaseline();      // ADC instrinsic pedestal
 761   Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
 762   //Int_t tc = fFeeParam->GetPFtimeConstant();   // this makes no sense yet
 763
 764   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 765     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 766       fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
 767     }
 768   }
 769 }
 770
 771 //_____________________________________________________________________________
 772 void AliTRDmcmSim::FilterGain()
 773 {
 774   //
 775   // Apply gain filter (not implemented)
 776   // Later it will be implemented because gain digital filiter will
 777   // increase noise level.
 778   //
 779
 780 }
 781
 782 //_____________________________________________________________________________
 783 void AliTRDmcmSim::FilterTail()
 784 {
 785   //
 786   // Apply exponential tail filter (Bogdan's version)
 787   //
 788
 789   Double_t *dtarg  = new Double_t[fNTimeBin];
 790   Int_t    *itarg  = new Int_t[fNTimeBin];
 791   Int_t     nexp   = fFeeParam->GetTFnExp();
 792   Int_t     tftype = fFeeParam->GetTFtype();
 793
 794   switch( tftype ) {
 795
 796   case 0: // Exponential Filter Analog Bogdan
 797     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 798       FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
 799       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 800         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
 801       }
 802     }
 803     break;
 804
 805   case 1: // Exponential filter digital Bogdan
 806     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 807       FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
 808       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 809         fADCF[iCol][iTime] = itarg[iTime];
 810       }
 811     }
 812     break;
 813
 814   case 2: // Exponential filter Marian special
 815     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 816       FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
 817       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 818         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
 819       }
 820     }
 821     break;
 822
 823     //new
 824   case 3: // Exponential filter using AliTRDtrapAlu class
 825     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 826       FilterSimDeConvExpEl( fADCF[iCol], itarg, fNTimeBin, nexp);
 827       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 828         fADCF[iCol][iTime] = itarg[iTime]>>2; // to be used for raw-data
 829         fADCT[iCol][iTime] = itarg[iTime];    // 12bits; to be used for tracklet; tracklet will have own container;
 830       }
 831     }
 832     break;
 833
 834
 835   default:
 836     AliError(Form("Invalid filter type %d ! \n", tftype ));
 837     break;
 838   }
 839
 840   delete [] dtarg;
 841   delete [] itarg;
 842
 843 }
 844
 845 //_____________________________________________________________________________
 846 void AliTRDmcmSim::ZSMapping()
 847 {
 848   //
 849   // Zero Suppression Mapping implemented in TRAP chip
 850   //
 851   // See detail TRAP manual "Data Indication" section:
 852   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
 853   //
 854
 855   Int_t eBIS = fFeeParam->GetEBsglIndThr();       // TRAP default = 0x4  (Tis=4)
 856   Int_t eBIT = fFeeParam->GetEBsumIndThr();       // TRAP default = 0x28 (Tit=40)
 857   Int_t eBIL = fFeeParam->GetEBindLUT();          // TRAP default = 0xf0
 858                                                   // (lookup table accept (I2,I1,I0)=(111)
 859                                                   // or (110) or (101) or (100))
 860   Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
 861   Int_t ep   = AliTRDfeeParam::GetPFeffectPedestal();
 862
 863   if( !CheckInitialized() ) return;
 864
 865   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
 866     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 867
 868       // Get ADC data currently in filter buffer
 869       Int_t ap = fADCF[iadc-1][it] - ep; // previous
 870       Int_t ac = fADCF[iadc  ][it] - ep; // current
 871       Int_t an = fADCF[iadc+1][it] - ep; // next
 872
 873       // evaluate three conditions
 874       Int_t i0 = ( ac >=  ap && ac >=  an ) ? 0 : 1; // peak center detection
 875       Int_t i1 = ( ap + ac + an > eBIT )    ? 0 : 1; // cluster
 876       Int_t i2 = ( ac > eBIS )              ? 0 : 1; // absolute large peak
 877
 878       Int_t i = i2 * 4 + i1 * 2 + i0;    // Bit position in lookup table
 879       Int_t d = (eBIL >> i) & 1;         // Looking up  (here d=0 means true
 880                                          // and d=1 means false according to TRAP manual)
 881
 882       fZSM[iadc][it] &= d;
 883       if( eBIN == 0 ) {  // turn on neighboring ADCs
 884         fZSM[iadc-1][it] &= d;
 885         fZSM[iadc+1][it] &= d;
 886       }
 887
 888     }
 889   }
 890
 891   // do 1 dim projection
 892   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 893     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 894       fZSM1Dim[iadc] &= fZSM[iadc][it];
 895     }
 896   }
 897
 898 }
 899
 900 //_____________________________________________________________________________
 901 void AliTRDmcmSim::DumpData( char *f, char *target )
 902 {
 903   //
 904   // Dump data stored (for debugging).
 905   // target should contain one or multiple of the following characters
 906   //   R   for raw data
 907   //   F   for filtered data
 908   //   Z   for zero suppression map
 909   //   S   Raw dat astream
 910   // other characters are simply ignored
 911   //
 912
 913   UInt_t tempbuf[1024];
 914
 915   if( !CheckInitialized() ) return;
 916
 917   std::ofstream of( f, std::ios::out | std::ios::app );
 918   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
 919              fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos );
 920
 921   for( int t=0 ; target[t] != 0 ; t++ ) {
 922     switch( target[t] ) {
 923     case 'R' :
 924     case 'r' :
 925       of << Form("fADCR (raw ADC data)\n");
 926       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 927         of << Form("  ADC %02d: ", iadc);
 928         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 929           of << Form("% 4d",  fADCR[iadc][it]);
 930         }
 931         of << Form("\n");
 932       }
 933       break;
 934     case 'F' :
 935     case 'f' :
 936       of << Form("fADCF (filtered ADC data)\n");
 937       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 938         of << Form("  ADC %02d: ", iadc);
 939         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 940           of << Form("% 4d",  fADCF[iadc][it]);
 941         }
 942         of << Form("\n");
 943       }
 944       break;
 945     case 'Z' :
 946     case 'z' :
 947       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
 948       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 949         of << Form("  ADC %02d: ", iadc);
 950         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
 951         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 952           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
 953         }
 954         of << Form("\n");
 955       }
 956       break;
 957     case 'S' :
 958     case 's' :
 959       Int_t s = ProduceRawStream( tempbuf, 1024 );
 960       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
 961       of << Form("  address  data\n");
 962       for( int i = 0 ; i < s ; i++ ) {
 963         of << Form("  %04x     %08x\n", i, tempbuf[i]);
 964       }
 965     }
 966   }
 967 }
 968
 969 //_____________________________________________________________________________
 970 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
 971                                      , Int_t n, Int_t nexp)
 972 {
 973   //
 974   // Exponential filter "analog"
 975   // source will not be changed
 976   //
 977
 978   Int_t    i = 0;
 979   Int_t    k = 0;
 980   Double_t reminder[2];
 981   Double_t correction;
 982   Double_t result;
 983   Double_t rates[2];
 984   Double_t coefficients[2];
 985
 986   // Initialize (coefficient = alpha, rates = lambda)
 987   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 988
 989   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
 990   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
 991   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
 992   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
 993
 994   coefficients[0] = c1;
 995   coefficients[1] = c2;
 996
 997   Double_t dt = 0.1;
 998   rates[0] = TMath::Exp(-dt/(r1));
 999   rates[1] = TMath::Exp(-dt/(r2));
1000
1001   // Attention: computation order is important
1002   correction = 0.0;
1003   for (k = 0; k < nexp; k++) {
1004     reminder[k] = 0.0;
1005   }
1006
1007   for (i = 0; i < n; i++) {
1008
1009     result    = ((Double_t)source[i] - correction);    // no rescaling
1010     target[i] = result;
1011
1012     for (k = 0; k < nexp; k++) {
1013       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1014     }
1015
1016     correction = 0.0;
1017     for (k = 0; k < nexp; k++) {
1018       correction += reminder[k];
1019     }
1020   }
1021 }
1022
1023 //_____________________________________________________________________________
1024 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
1025                                      , Int_t nexp)
1026 {
1027   //
1028   // Exponential filter "digital"
1029   // source will not be changed
1030   //
1031
1032   Int_t i        = 0;
1033   Int_t fAlphaL  = 0;
1034   Int_t fAlphaS  = 0;
1035   Int_t fTailPed = 0;
1036   Int_t iAlphaL  = 0;
1037   Int_t iAlphaS  = 0;
1038
1039   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
1040   // initialize (coefficient = alpha, rates = lambda)
1041
1042   Double_t dt = 0.1;
1043   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1044   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1045   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1046   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1047
1048   Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
1049   Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
1050   Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; //  9 bit paramter + fixed bits
1051   Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; //  9 bit paramter + fixed bits
1052
1053   if (nexp == 1) {
1054     fAlphaL = (Int_t) (c1 * 2048.0);
1055     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
1056   }
1057   if (nexp == 2) {
1058     fAlphaL = (Int_t) (c1 * 2048.0);
1059     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
1060     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
1061     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
1062   }
1063
1064   Double_t iAl = iAlphaL  / 2048.0;            // alpha  L: correspondence to floating point numbers
1065   Double_t iAs = iAlphaS  / 2048.0;            // alpha  S: correspondence to floating point numbers
1066   Double_t iLl = iLambdaL / 2048.0;            // lambda L: correspondence to floating point numbers
1067   Double_t iLs = iLambdaS / 2048.0;            // lambda S: correspondence to floating point numbers
1068
1069   Int_t h1;
1070   Int_t h2;
1071   Int_t rem1;
1072   Int_t rem2;
1073   Int_t correction;
1074   Int_t result;
1075   Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
1076
1077   Double_t xi = 1 - (iLl*iAs + iLs*iAl);             // Calculation of equilibrium values of the
1078   rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
1079   rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
1080
1081   // further initialization
1082   if ((rem1 + rem2) > 0x0FFF) {
1083     correction = 0x0FFF;
1084   }
1085   else {
1086     correction = (rem1 + rem2) & 0x0FFF;
1087   }
1088
1089   fTailPed = iFactor - correction;
1090
1091   for (i = 0; i < n; i++) {
1092
1093     result = (source[i]  - correction);
1094     if (result < 0) { // Too much undershoot
1095       result = 0;
1096     }
1097
1098     target[i] = result;
1099
1100     h1 = (rem1 + ((iAlphaL * result) >> 11));
1101     if (h1 > 0x0FFF) {
1102       h1 = 0x0FFF;
1103     }
1104     else {
1105       h1 &= 0x0FFF;
1106     }
1107
1108     h2 = (rem2 + ((iAlphaS * result) >> 11));
1109     if (h2 > 0x0FFF) {
1110       h2 = 0x0FFF;
1111     }
1112     else {
1113       h2 &= 0x0FFF;
1114     }
1115
1116     rem1 = (iLambdaL * h1 ) >> 11;
1117     rem2 = (iLambdaS * h2 ) >> 11;
1118
1119     if ((rem1 + rem2) > 0x0FFF) {
1120       correction = 0x0FFF;
1121     }
1122     else {
1123       correction = (rem1 + rem2) & 0x0FFF;
1124     }
1125
1126   }
1127
1128 }
1129
1130 //_____________________________________________________________________________
1131 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
1132                                       , Int_t n)
1133 {
1134   //
1135   // Exponential filter (M. Ivanov)
1136   // source will not be changed
1137   //
1138
1139   Int_t i = 0;
1140   Double_t sig1[100];
1141   Double_t sig2[100];
1142   Double_t sig3[100];
1143
1144   for (i = 0; i < n; i++) {
1145     sig1[i] = (Double_t)source[i];
1146   }
1147
1148   Float_t dt      = 0.1;
1149   Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
1150   Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
1151
1152   FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
1153   FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
1154
1155   for (i = 0; i < n; i++) {
1156     target[i] = sig3[i];
1157   }
1158
1159 }
1160
1161 //______________________________________________________________________________
1162 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
1163                                           , Double_t lambda, Int_t n)
1164 {
1165   //
1166   // Special filter (M. Ivanov)
1167   //
1168
1169   Int_t    i = 0;
1170   Double_t l = TMath::Exp(-lambda*0.5);
1171   Double_t in[1000];
1172   Double_t out[1000];
1173
1174   // Initialize in[] and out[] goes 0 ... 2*n+19
1175   for (i = 0; i < n*2+20; i++) {
1176     in[i] = out[i] = 0;
1177   }
1178
1179   // in[] goes 0, 1
1180   in[0] = ampin[0];
1181   in[1] = (ampin[0] + ampin[1]) * 0.5;
1182
1183   // Add charge to the end
1184   for (i = 0; i < 22; i++) {
1185     in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
1186   }
1187
1188   // Use arithmetic mean
1189   for (i = 1; i < n-1; i++) {
1190     in[2*i]   = ampin[i];    // in[] goes 2, 3, ... , 2*n-4, 2*n-3
1191     in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
1192   }
1193
1194   Double_t temp;
1195   out[2*n]    = in[2*n];
1196   temp        = 0;
1197   for (i = 2*n; i >= 0; i--) {
1198     out[i]    = in[i] + temp;
1199     temp      = l*(temp+in[i]);
1200   }
1201
1202   for (i = 0; i < n; i++){
1203     //ampout[i] = out[2*i+1];  // org
1204     ampout[i] = out[2*i];
1205   }
1206
1207 }
1208
1209 //______________________________________________________________________________
1210 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
1211                                             , Double_t norm, Double_t lambda
1212                                             , Int_t n)
1213 {
1214   //
1215   // Special filter (M. Ivanov)
1216   //
1217
1218   Int_t    i = 0;
1219
1220   Double_t l = TMath::Exp(-lambda*0.5);
1221   Double_t k = l*(1.0 - norm*lambda*0.5);
1222   Double_t in[1000];
1223   Double_t out[1000];
1224
1225   // Initialize in[] and out[] goes 0 ... 2*n+19
1226   for (i = 0; i < n*2+20; i++) {
1227     in[i] = out[i] = 0;
1228   }
1229
1230   // in[] goes 0, 1
1231   in[0] = ampin[0];
1232   in[1] = (ampin[0]+ampin[1])*0.5;
1233
1234   // Add charge to the end
1235   for (i =-2; i < 22; i++) {
1236     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
1237     in[2*(n-1)+i] = ampin[n-1];
1238   }
1239
1240   for (i = 1; i < n-2; i++) {
1241     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
1242     in[2*i]    = ampin[i];
1243     in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
1244     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
1245   }
1246
1247   Double_t temp;
1248   out[0] = in[0];
1249   temp   = in[0];
1250   for (i = 1; i <= 2*n; i++) {
1251     out[i] = in[i] + (k-l)*temp;
1252     temp   = in[i] +  k   *temp;
1253   }
1254
1255   for (i = 0; i < n; i++) {
1256     //ampout[i] = out[2*i+1];  // org
1257     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
1258     ampout[i] = TMath::Max(out[2*i],0.0);
1259   }
1260 }
1261
1262
1263 //_____________________________________________________________________________________
1264 //the following filter uses AliTRDtrapAlu-class
1265
1266 void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) {
1267   //static Int_t count = 0;
1268
1269   Double_t dt = 0.1;
1270   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1271   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1272   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1273   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1274
1275   nexp = 1;
1276
1277   //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual
1278   //parameters need to be adjusted
1279   AliTRDtrapAlu lambdaL;
1280   AliTRDtrapAlu lambdaS;
1281   AliTRDtrapAlu alphaL;
1282   AliTRDtrapAlu alphaS;
1283
1284   AliTRDtrapAlu correction;
1285   AliTRDtrapAlu result;
1286   AliTRDtrapAlu bufL;
1287   AliTRDtrapAlu bufS;
1288
1289   AliTRDtrapAlu bSource;
1290
1291   lambdaL.Init(1,11);
1292   lambdaS.Init(1,11);
1293   alphaL.Init(1,11);
1294   alphaS.Init(1,11);
1295
1296   //count=count+1;
1297
1298   lambdaL.AssignDouble(TMath::Exp(-dt/r1));
1299   lambdaS.AssignDouble(TMath::Exp(-dt/r2));
1300   alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants
1301   alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL
1302
1303   //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly
1304   correction.Init(10,2);
1305   result.Init(10,2);
1306   bufL.Init(10,2);
1307   bufS.Init(10,2);
1308   bSource.Init(10,2);
1309
1310   for(Int_t i = 0; i < n; i++) {
1311     bSource.AssignInt(source[i]);
1312     result = bSource - correction; // subtraction can produce an underflow
1313     if(result.GetSign() == kTRUE) {
1314       result.AssignInt(0);
1315     }
1316
1317     //target[i] = result.GetValuePre();  // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly
1318
1319     target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2
1320
1321     //printf("target-Wert zur Zeit %d : %d",i,target[i]);
1322     //printf("\n");
1323
1324     bufL  =  bufL + (result * alphaL);
1325     bufL  =  bufL * lambdaL;
1326
1327     bufS  =  bufS + (result * alphaS);
1328     bufS  =  bufS * lambdaS;  // eventually this should look like:
1329     // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more
1330
1331     correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF
1332   }
1333
1334
1335 }
1336
1337
1338
1339
1340
1341
1342
1343 //__________________________________________________________________________________
1344
1345
1346 // in order to use the Tracklets, please first
1347 // -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
1348 // -- set AliTRDfeeParam::fgkTFtype   to 3, in order to use the new tail cancellation filter
1349 //    currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
1350 // -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with           information about the Tracklet position (MCM, channel number)
1351
1352 // 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
1353
1354 void AliTRDmcmSim::Tracklet(){
1355     // tracklet calculation
1356     // if you use this code after a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
1357
1358   if(!CheckInitialized()){ return; }
1359
1360   Bool_t filtered = kTRUE;
1361
1362
1363
1364   AliTRDtrapAlu data;
1365   data.Init(10,2);
1366   if(fADCT[0][0]==-1){                      // check if filter was applied
1367     filtered = kFALSE;
1368     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1369       for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) {
1370         data.AssignInt(fADCR[iadc][iT]);
1371         fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called, this is done correctly
1372       }
1373     }
1374
1375   }
1376
1377   // the online ordering of mcm's is reverse to the TRAP-manual-ordering! reverse fADCT (to be consistent to TRAP), then do all calculations
1378   // reverse fADCT:
1379   Int_t** rev0 = new Int_t *[fNADC];
1380   Int_t** rev1 = new Int_t *[fNADC];
1381
1382   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1383     rev0[iadc] = new Int_t[fNTimeBin];
1384     rev1[iadc] = new Int_t[fNTimeBin];
1385     for( Int_t iT = 0; iT < fNTimeBin; iT++) {
1386       if( iadc <= fNADC-iadc-1 ) {
1387         rev0[iadc][iT]  = fADCT[fNADC-iadc-1][iT];
1388         rev1[iadc][iT]  = fADCT[iadc][iT];
1389         fADCT[iadc][iT] = rev0[iadc][iT];
1390       }
1391       else {
1392         rev0[iadc][iT]  = rev1[fNADC-iadc-1][iT];
1393         fADCT[iadc][iT] = rev0[iadc][iT];
1394       }
1395     }
1396   }
1397   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1398     delete[] rev0[iadc];
1399     delete[] rev1[iadc];
1400   }
1401
1402   delete[] rev0;
1403   delete[] rev1;
1404
1405   rev0 = NULL;
1406   rev1 = NULL;
1407
1408   // get the filtered pedestal; supports only electronic tail-cancellation filter
1409   AliTRDtrapAlu filPed;
1410   Int_t ep = 0;
1411   Int_t *ieffped = new Int_t[fNTimeBin];
1412   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1413     ieffped[iT] = ep;
1414   }
1415
1416   if( filtered == kTRUE ) {
1417     if( fFeeParam->IsPFon() ){
1418       ep = fFeeParam->GetPFeffectPedestal();
1419     }
1420     Int_t      nexp  = fFeeParam->GetTFnExp();
1421     Int_t  *isource  = new Int_t[fNTimeBin];
1422     filPed.Init(10,2);
1423     filPed.AssignInt(ep);
1424     Int_t epf = filPed.GetValue();
1425     for(Int_t iT = 0; iT < fNTimeBin; iT++){
1426       isource[iT] = ep;
1427       ieffped[iT] = epf;
1428     }
1429
1430     if( fFeeParam->IsTFon() ) {
1431       FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp);
1432     }
1433
1434     delete[] isource;
1435   }
1436
1437   //the following values should go to AliTRDfeeParam once they are defined; then they have to be read in properly
1438   //naming follows conventions in TRAP-manual
1439
1440
1441   Bool_t bVBY = kTRUE;                         // cluster-verification bypass
1442
1443   Double_t cQTParam = 0;                      // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10
1444   AliTRDtrapAlu cQTAlu;
1445   cQTAlu.Init(1,10,0,63);
1446   cQTAlu.AssignDouble(cQTParam);
1447   Int_t cQT = cQTAlu.GetValue();
1448
1449   // linear fit
1450   Int_t tFS = fFeeParam->GetLinearFitStart();  // linear fit start
1451   Int_t tFE = fFeeParam->GetLinearFitEnd();    // linear fit stop
1452
1453   // charge accumulators
1454   Int_t tQS0 = fFeeParam->GetQacc0Start();     // start-time for charge-accumulator 0
1455   Int_t tQE0 = fFeeParam->GetQacc0End();       // stop-time for charge-accumulator 0
1456   Int_t tQS1 = fFeeParam->GetQacc1Start();     // start-time for charge-accumulator 1
1457   Int_t tQE1 = fFeeParam->GetQacc1End();       // stop-time for charge-accumulator 1
1458   // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N
1459
1460   Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold
1461   AliTRDtrapAlu cTHAlu;
1462   cTHAlu.Init(12,2);
1463   cTHAlu.AssignDouble(cTHParam);
1464   Int_t cTH = cTHAlu.GetValue();                                 // cTH used for comparison
1465
1466   struct List_t {
1467     List_t *next;
1468     Int_t iadc;
1469     Int_t value;
1470   };
1471
1472   List_t selection[7];            // list with 7 elements
1473   List_t *list = NULL;
1474   List_t *listLeft = NULL;
1475
1476   Int_t* qsum = new Int_t[fNADC];
1477
1478   // fit sums
1479   AliTRDtrapAlu qsumAlu;
1480   qsumAlu.Init(12,2);           // charge sum will be 12+2 bits
1481   AliTRDtrapAlu dCOGAlu;
1482   dCOGAlu.Init(1,7,0,127);      // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT
1483   AliTRDtrapAlu yrawAlu;
1484   yrawAlu.Init(1,8,-1,255);
1485   AliTRDtrapAlu yAlu;
1486   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 -
1487   AliTRDtrapAlu xAlu;
1488   xAlu.Init(5,8);               // 8 past-comma bits because value will be added/multiplied to another value with this accuracy
1489   AliTRDtrapAlu xxAlu;
1490   xxAlu.Init(10,0);
1491   AliTRDtrapAlu yyAlu;
1492   yyAlu.Init(1,16,0,0xFFFF);    // maximum is 2^16-1; 16Bit for past-commas
1493   AliTRDtrapAlu xyAlu;
1494   xyAlu.Init(6,8);
1495   AliTRDtrapAlu XAlu;
1496   XAlu.Init(9,0);
1497   AliTRDtrapAlu XXAlu;
1498   XXAlu.Init(14,0);
1499   AliTRDtrapAlu YAlu;
1500   YAlu.Init(5,8);               // 14 bit, 1 is sign-bit; therefore only 13 bit
1501   AliTRDtrapAlu YYAlu;
1502   YYAlu.Init(5,16);
1503   AliTRDtrapAlu XYAlu;
1504   XYAlu.Init(8,8);              // 17 bit, 1 is sign-bit; therefore only 16 bit
1505   AliTRDtrapAlu qtruncAlu;
1506   qtruncAlu.Init(12,0);
1507   AliTRDtrapAlu QT0Alu;
1508   QT0Alu.Init(15,0);
1509   AliTRDtrapAlu QT1Alu;
1510   QT1Alu.Init(16,0);
1511
1512   AliTRDtrapAlu oneAlu;
1513   oneAlu.Init(1,8);
1514
1515
1516   AliTRDtrapAlu inverseNAlu;
1517   inverseNAlu.Init(1,8);        // simulates the LUT for 1/N
1518   AliTRDtrapAlu MeanChargeAlu;  // mean charge in ADC counts
1519   MeanChargeAlu.Init(8,0);
1520   AliTRDtrapAlu TotalChargeAlu;
1521   TotalChargeAlu.Init(17,8);
1522   //nr of post comma bits should be the same for inverseN and TotalCharge
1523
1524
1525   SetPosLUT();                    // initialize the position correction LUT for this MCM;
1526
1527
1528   // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
1529   Int_t *X   = new Int_t[fNADC-2];
1530   Int_t *XX  = new Int_t[fNADC-2];
1531   Int_t *Y   = new Int_t[fNADC-2];
1532   Int_t *YY  = new Int_t[fNADC-2];
1533   Int_t *XY  = new Int_t[fNADC-2];
1534   Int_t *N   = new Int_t[fNADC-2];
1535   Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
1536   Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
1537
1538   for (Int_t iCol = 0; iCol < fNADC-2; iCol++) {
1539
1540       // initialize fit-sums
1541       X[iCol]   = 0;
1542       XX[iCol]  = 0;
1543       Y[iCol]   = 0;
1544       YY[iCol]  = 0;
1545       XY[iCol]  = 0;
1546       N[iCol]   = 0;
1547       QT0[iCol] = 0;
1548       QT1[iCol] = 0;
1549   }
1550
1551
1552   filPed.Init(7,2);                         // convert filtered pedestal into 7+2Bits
1553
1554   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1555
1556     if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no?
1557
1558     // reset
1559     Int_t portChannel[4]   = {-1,-1,-1,-1};
1560     Int_t clusterCharge[4] = {0,0,0,0};
1561     Int_t leftCharge[4]    = {0,0,0,0};
1562     Int_t centerCharge[4]  = {0,0,0,0};
1563     Int_t rightCharge[4]   = {0,0,0,0};
1564
1565     Int_t mark = 0;
1566
1567     filPed.AssignFormatted(ieffped[iT]);   // no size-checking when using AssignFormatted; ieffped>=0
1568     filPed = filPed;                       // this checks the size
1569
1570     ieffped[iT] = filPed.GetValue();
1571
1572     for(Int_t i = 0; i<7; i++){
1573       selection[i].next       = NULL;
1574       selection[i].iadc       =   -1;     // value of -1: invalid adc
1575       selection[i].value      =    0;
1576
1577     }
1578     // selection[0] is starting list-element; just for pointing
1579
1580     // loop over inner adc's
1581     for (Int_t iCol = 1; iCol < fNADC-1; iCol++) {
1582
1583       Int_t left   = fADCT[iCol-1][iT];
1584       Int_t center = fADCT[iCol][iT];
1585       Int_t right  = fADCT[iCol+1][iT];
1586
1587       Int_t sum = left + center + right;            // cluster charge sum
1588       qsumAlu.AssignFormatted(sum);
1589       qsumAlu = qsumAlu;                        // size-checking; redundant
1590
1591       qsum[iCol] = qsumAlu.GetValue();
1592
1593       //hit detection and masking
1594       if(center>=left){
1595         if(center>right){
1596           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
1597             mark |= 1;                              // marker
1598           }
1599         }
1600       }
1601       mark = mark<<1;
1602     }
1603     mark = mark>>1;
1604
1605
1606     // get selection of 6 adc's and sort,starting with greatest values
1607
1608     //read three from right side and sort (primitive sorting algorithm)
1609     Int_t i = 0; // adc number
1610     Int_t j = 1; // selection number
1611     while(i<fNADC-2 && j<=3){
1612       i = i + 1;
1613       if( ((mark>>(i-1)) & 1) == 1) {
1614         selection[j].iadc  = fNADC-1-i;
1615         selection[j].value = qsum[fNADC-1-i]>>6;   // for hit-selection only the first 8 out of the 14 Bits are used for comparison
1616
1617         // insert into sorted list
1618         listLeft = &selection[0];
1619         list = listLeft->next;
1620
1621         if(list!=NULL) {
1622           while((list->next != NULL) && (selection[j].value <= list->value)){
1623             listLeft = list;
1624             list = list->next;
1625           }
1626
1627           if(selection[j].value<=list->value){
1628             selection[j].next = list->next;
1629             list->next = &selection[j];
1630           }
1631           else {
1632             listLeft->next = &selection[j];
1633             selection[j].next = list;
1634           }
1635         }
1636         else{
1637           listLeft->next = &selection[j];
1638           selection[j].next = list;
1639         }
1640
1641         j = j + 1;
1642       }
1643     }
1644
1645
1646     // read three from left side
1647     Int_t k = fNADC-2;
1648     while(k>i && j<=6) {
1649       if( ((mark>>(k-1)) & 1) == 1) {
1650         selection[j].iadc  = fNADC-1-k;
1651         selection[j].value = qsum[fNADC-1-k]>>6;
1652
1653         listLeft = &selection[0];
1654         list = listLeft->next;
1655
1656         if(list!=NULL){
1657           while((list->next != NULL) && (selection[j].value <= list->value)){
1658             listLeft = list;
1659             list = list->next;
1660           }
1661
1662           if(selection[j].value<=list->value){
1663             selection[j].next = list->next;
1664             list->next = &selection[j];
1665           }
1666           else {
1667             listLeft->next = &selection[j];
1668             selection[j].next = list;
1669           }
1670         }
1671         else{
1672           listLeft->next = &selection[j];
1673           selection[j].next = list;
1674         }
1675
1676         j = j + 1;
1677       }
1678       k = k - 1;
1679     }
1680
1681     // get the four with greatest charge-sum
1682     list = &selection[0];
1683     for(i = 0; i<4; i++){
1684       if(list->next == NULL) continue;
1685       list = list->next;
1686       if(list->iadc == -1) continue;
1687       Int_t adc = list->iadc;                              // channel number with selected hit
1688
1689       // the following arrays contain the four chosen channels in 1 time-bin
1690       portChannel[i]   = adc;
1691       clusterCharge[i] = qsum[adc];
1692       leftCharge[i]    = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
1693       centerCharge[i]  = fADCT[adc][iT] - ieffped[iT];
1694       rightCharge[i]   = fADCT[adc+1][iT] - ieffped[iT];
1695     }
1696
1697     // arithmetic unit
1698
1699     // cluster verification
1700     if(!bVBY){
1701       for(i = 0; i<4; i++){
1702         Int_t lr = leftCharge[i]*rightCharge[i]*1024;
1703         Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
1704         if (lr>=cc){
1705           portChannel[i]   = -1;                                 // set to invalid address
1706           clusterCharge[i] = 0;
1707         }
1708       }
1709     }
1710
1711     // fit-sums of valid channels
1712     // local hit position
1713     for(i = 0; i<4; i++){
1714       if (centerCharge[i] ==  0) {
1715         portChannel[i] = -1;
1716       }// prevent division by 0
1717
1718       if (portChannel[i]  == -1) continue;
1719
1720       Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
1721
1722       Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
1723       dCOG = (sign == kFALSE) ? dCOG : -dCOG;     // AssignDouble doesn't allow for signed doubles
1724       dCOGAlu.AssignDouble(dCOG);
1725       Int_t iLUTpos = dCOGAlu.GetValue();       // steers position in LUT
1726
1727       dCOG = dCOG/2;
1728       yrawAlu.AssignDouble(dCOG);
1729       Int_t iCOG = yrawAlu.GetValue();
1730       Int_t y = iCOG + fPosLUT[iLUTpos % 128];    // local position in pad-units
1731       yrawAlu.AssignFormatted(y);               // 0<y<1
1732       yAlu  = yrawAlu;                        // convert to 16 past-comma bits
1733
1734       if(sign == kTRUE) yAlu.SetSign(-1);       // buffer width of 9 bits; sign on real (not estimated) position
1735       xAlu.AssignInt(iT);                       // buffer width of 5 bits
1736
1737
1738       xxAlu = xAlu * xAlu;                  // buffer width of 10 bits -> fulfilled by x*x
1739
1740       yyAlu = yAlu * yAlu;                  // buffer width of 16 bits
1741
1742       xyAlu = xAlu * yAlu;                  // buffer width of 14 bits
1743
1744       Int_t adc = portChannel[i]-1;              // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
1745
1746       // calculate fit-sums recursively
1747       // interpretation of their bit-length is given as comment
1748
1749       // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
1750
1751       XAlu.AssignFormatted(X[adc]);
1752       XAlu = XAlu + xAlu;                   // buffer width of 9 bits
1753       X[adc] = XAlu.GetValue();
1754
1755       XXAlu.AssignFormatted(XX[adc]);
1756       XXAlu = XXAlu + xxAlu;                // buffer width of 14 bits
1757       XX[adc] = XXAlu.GetValue();
1758
1759       if (Y[adc] < 0) {
1760         YAlu.AssignFormatted(-Y[adc]);          // make sure that only positive values are assigned; sign-setting must be done by hand
1761         YAlu.SetSign(-1);
1762       }
1763       else {
1764         YAlu.AssignFormatted(Y[adc]);
1765         YAlu.SetSign(1);
1766       }
1767
1768       YAlu = YAlu + yAlu;                   // buffer width of 14 bits (8 past-comma);
1769       Y[adc] = YAlu.GetSignedValue();
1770
1771       YYAlu.AssignFormatted(YY[adc]);
1772       YYAlu = YYAlu + yyAlu;                // buffer width of 21 bits (16 past-comma)
1773       YY[adc] = YYAlu.GetValue();
1774
1775       if (XY[adc] < 0) {
1776         XYAlu.AssignFormatted(-XY[adc]);
1777         XYAlu.SetSign(-1);
1778       }
1779       else {
1780         XYAlu.AssignFormatted(XY[adc]);
1781         XYAlu.SetSign(1);
1782       }
1783
1784       XYAlu = XYAlu + xyAlu;                // buffer allows 17 bits (8 past-comma)
1785       XY[adc] = XYAlu.GetSignedValue();
1786
1787       N[adc]  = N[adc] + 1;
1788
1789
1790       // accumulated charge
1791       qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
1792       qtruncAlu = qsumAlu;
1793
1794       if(iT>=tQS0 && iT<=tQE0){
1795         QT0Alu.AssignFormatted(QT0[adc]);
1796         QT0Alu = QT0Alu + qtruncAlu;
1797         QT0[adc] = QT0Alu.GetValue();
1798         //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
1799       }
1800
1801       if(iT>=tQS1 && iT<=tQE1){
1802         QT1Alu.AssignFormatted(QT1[adc]);
1803         QT1Alu = QT1Alu + qtruncAlu;
1804         QT1[adc] = QT1Alu.GetValue();
1805         //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
1806       }
1807     }// i
1808
1809     // remapping is done!!
1810
1811   }//iT
1812
1813
1814
1815   // tracklet-assembly
1816
1817   // put into AliTRDfeeParam and take care that values are in proper range
1818   const Int_t cTCL = 1;      // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??)
1819   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%)
1820
1821   Int_t mPair   = 0;         // marker for possible tracklet pairs
1822   Int_t* hitSum = new Int_t[fNADC-3];
1823   // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18;
1824
1825   // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
1826   for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
1827     hitSum[iCol] = N[iCol] + N[iCol+1];
1828     if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
1829         mPair |= 1;         // mark as possible channel-pair
1830
1831     }
1832     mPair = mPair<<1;
1833   }
1834   mPair = mPair>>1;
1835
1836   List_t* selectPair = new List_t[fNADC-2];      // list with 18 elements (0..18) containing the left channel-nr and hit sums
1837                                                  // selectPair[18] is starting list-element just for pointing
1838   for(Int_t k = 0; k<fNADC-2; k++){
1839       selectPair[k].next       = NULL;
1840       selectPair[k].iadc       =   -1;           // invalid adc
1841       selectPair[k].value      =    0;
1842
1843     }
1844
1845  list = NULL;
1846  listLeft = NULL;
1847
1848   // read marker and sort according to hit-sum
1849
1850   Int_t adcL  = 0;            // left adc-channel-number (remapped)
1851   Int_t selNr = 0;            // current number in list
1852
1853   // insert marked channels into list and sort according to hit-sum
1854   while(adcL < fNADC-3 && selNr < fNADC-3){
1855
1856     if( ((mPair>>((fNADC-4)-(adcL))) & 1) == 1) {
1857       selectPair[selNr].iadc  = adcL;
1858       selectPair[selNr].value = hitSum[adcL];
1859
1860       listLeft = &selectPair[fNADC-3];
1861       list = listLeft->next;
1862
1863       if(list!=NULL) {
1864         while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
1865           listLeft = list;
1866           list = list->next;
1867         }
1868
1869         if(selectPair[selNr].value <= list->value){
1870           selectPair[selNr].next = list->next;
1871           list->next = &selectPair[selNr];
1872         }
1873         else {
1874           listLeft->next = &selectPair[selNr];
1875           selectPair[selNr].next = list;
1876         }
1877
1878       }
1879       else{
1880         listLeft->next = &selectPair[selNr];
1881         selectPair[selNr].next = list;
1882       }
1883
1884       selNr = selNr + 1;
1885     }
1886     adcL = adcL + 1;
1887   }
1888
1889   //select up to 4 channels with maximum number of hits
1890   Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
1891   Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
1892   list = &selectPair[fNADC-3];
1893
1894   for (Int_t i = 0; i<4; i++) {
1895     if(list->next == NULL) continue;
1896     list = list->next;
1897     if(list->iadc == -1) continue;
1898     lpairChannel[i] = list->iadc;        // channel number with selected hit
1899     rpairChannel[i] = lpairChannel[i]+1;
1900   }
1901
1902   // avoid submission of double tracklets
1903   for (Int_t i = 3; i>0; i--) {
1904     for (Int_t j = i-1; j>-1; j--) {
1905       if(lpairChannel[i] == rpairChannel[j]) {
1906         lpairChannel[i] = -1;
1907         rpairChannel[i] = -1;
1908         break;
1909       }
1910       /* if(rpairChannel[i] == lpairChannel[j]) {
1911         lpairChannel[i] = -1;
1912         rpairChannel[i] = -1;
1913         break;
1914         }*/
1915     }
1916   }
1917
1918   // merging of the fit-sums of the remainig channels
1919   // assume same data-word-width as for fit-sums for 1 channel
1920   // relative scales!
1921   Int_t mADC[4];
1922   Int_t mN[4];
1923   Int_t mQT0[4];
1924   Int_t mQT1[4];
1925   Int_t mX[4];
1926   Int_t mXX[4];
1927   Int_t mY[4];
1928   Int_t mYY[4];
1929   Int_t mXY[4];
1930   Int_t mOffset[4];
1931   Int_t mSlope[4];
1932   Int_t mMeanCharge[4];
1933   Int_t inverseN = 0;
1934   Double_t invN = 0;
1935   Int_t one = 0;
1936
1937   for (Int_t i = 0; i<4; i++){
1938     mADC[i] = -1;                        // set to invalid number
1939     mN[i]   =  0;
1940     mQT0[i] =  0;
1941     mQT1[i] =  0;
1942     mX[i]   =  0;
1943     mXX[i]  =  0;
1944     mY[i]   =  0;
1945     mYY[i]  =  0;
1946     mXY[i]  =  0;
1947     mOffset[i] = 0;
1948     mSlope[i]  = 0;
1949     mMeanCharge[i] = 0;
1950   }
1951
1952   oneAlu.AssignInt(1);
1953   one = oneAlu.GetValue();              // one with 8 past comma bits
1954
1955   for (Int_t i = 0; i<4; i++){
1956
1957
1958     mADC[i] = lpairChannel[i];          // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
1959                                         // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
1960     Int_t madc = mADC[i];
1961     if (madc == -1) continue;
1962
1963     YAlu.AssignInt(N[rpairChannel[i]]);
1964     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)
1965
1966     mN[i]    = hitSum[madc];
1967
1968     // 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
1969     if (N[madc+1] == 0) {
1970         mQT0[i] = QT0[madc];
1971         mQT1[i] = QT1[madc];
1972
1973     }
1974     else {
1975
1976         // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
1977
1978         mQT0[i]   = QT0[madc] + QT0[madc+1];
1979         QT0Alu.AssignFormatted(mQT0[i]);
1980         QT0Alu  = QT0Alu;                // size-check
1981         mQT0[i]   = QT0Alu.GetValue();     // write back
1982
1983         mQT1[i]   = QT1[madc] + QT1[madc+1];
1984         QT1Alu.AssignFormatted(mQT1[i]);
1985         QT1Alu  = QT1Alu;
1986         mQT1[i]   = QT1Alu.GetValue();
1987     }
1988
1989     // calculate the mean charge in adc values; later to be replaced by electron likelihood
1990     mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
1991     mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
1992     // 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
1993     invN = 1.0/(mN[i]);
1994     inverseNAlu.AssignDouble(invN);
1995     inverseN = inverseNAlu.GetValue();
1996     mMeanCharge[i] = mMeanCharge[i] * inverseN;  // now to be interpreted with 8 past-comma bits
1997     TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
1998     TotalChargeAlu = TotalChargeAlu;
1999     MeanChargeAlu = TotalChargeAlu;
2000     mMeanCharge[i] = MeanChargeAlu.GetValue();
2001
2002     // this check is not necessary; it is just for efficiency reasons
2003     if (N[madc+1] == 0) {
2004         mX[i]     =   X[madc];
2005         mXX[i]    =  XX[madc];
2006         mY[i]     =   Y[madc];
2007         mXY[i]    =  XY[madc];
2008         mYY[i]    =  YY[madc];
2009     }
2010     else {
2011
2012         mX[i]     =   X[madc] +  X[madc+1];
2013         XAlu.AssignFormatted(mX[i]);
2014         XAlu      = XAlu;
2015         mX[i]     = XAlu.GetValue();
2016
2017         mXX[i]    =  XX[madc] + XX[madc+1];
2018         XXAlu.AssignFormatted(mXX[i]);
2019         XXAlu     = XXAlu;
2020         mXX[i]    = XXAlu.GetValue();
2021
2022
2023         mY[i]     =   Y[madc] + Y[madc+1] + wpad;
2024         if (mY[i] < 0) {
2025             YAlu.AssignFormatted(-mY[i]);
2026             YAlu.SetSign(-1);
2027         }
2028         else {
2029             YAlu.AssignFormatted(mY[i]);
2030             YAlu.SetSign(1);
2031         }
2032         YAlu    = YAlu;
2033         mY[i]     = YAlu.GetSignedValue();
2034
2035         mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*one;    // multiplication by one to maintain the data format
2036
2037         if (mXY[i] < 0) {
2038             XYAlu.AssignFormatted(-mXY[i]);
2039             XYAlu.SetSign(-1);
2040         }
2041         else {
2042             XYAlu.AssignFormatted(mXY[i]);
2043             XYAlu.SetSign(1);
2044         }
2045         XYAlu   = XYAlu;
2046         mXY[i]    = XYAlu.GetSignedValue();
2047
2048         mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
2049         if (mYY[i] < 0) {
2050             YYAlu.AssignFormatted(-mYY[i]);
2051             YYAlu.SetSign(-1);
2052         }
2053         else {
2054             YYAlu.AssignFormatted(mYY[i]);
2055             YYAlu.SetSign(1);
2056         }
2057
2058         YYAlu   = YYAlu;
2059         mYY[i]    = YYAlu.GetSignedValue();
2060     }
2061
2062   }
2063
2064   // calculation of offset and slope from the merged fit-sums;
2065   // YY is needed for some error measure only; still to be done
2066   // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
2067
2068   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2069   // !!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   !!
2070   // !!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.    !!
2071   // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below).                                                                                             !!
2072   // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used.                                                                                          !!
2073   // !!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    !!
2074   // !!valid (in approximation) if tFS is close to the beginning of the drift region.                                                                                                                                                 !!
2075   // !!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                                                                    !!
2076   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2077
2078   // which formats should be chosen?
2079   AliTRDtrapAlu denomAlu;
2080   denomAlu.Init(20,8);
2081   AliTRDtrapAlu numAlu;
2082   numAlu.Init(20,8);
2083   // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
2084   // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
2085
2086   for (Int_t i = 0; i<4; i++) {
2087     if (mADC[i] == -1) continue;
2088
2089     Int_t num0  = (mN[i]*mXX[i]-mX[i]*mX[i]);
2090     if (num0 < 0) {
2091       denomAlu.AssignInt(-num0);    // num0 does not have to be interpreted as having past-comma bits -> AssignInt
2092       denomAlu.SetSign(-1);
2093     }
2094     else {
2095       denomAlu.AssignInt(num0);
2096       denomAlu.SetSign(1);
2097     }
2098
2099     Int_t num1  = mN[i]*mXY[i] - mX[i]*mY[i];
2100     if (num1 < 0) {
2101       numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
2102       numAlu.SetSign(-1);
2103     }
2104     else {
2105       numAlu.AssignFormatted(num1);
2106       numAlu.SetSign(1);
2107     }
2108     numAlu    = numAlu/denomAlu;
2109     mSlope[i]   = numAlu.GetSignedValue();
2110
2111     Int_t num2  = mXX[i]*mY[i] - mX[i]*mXY[i];
2112
2113     if (num2 < 0) {
2114       numAlu.AssignFormatted(-num2);
2115       numAlu.SetSign(-1);
2116     }
2117     else {
2118       numAlu.AssignFormatted(num2);
2119       numAlu.SetSign(1);
2120     }
2121
2122     numAlu    = numAlu/denomAlu;
2123
2124
2125     mOffset[i]  = numAlu.GetSignedValue();
2126     numAlu.SetSign(1);
2127     denomAlu.SetSign(1);
2128
2129
2130     //numAlu.AssignInt(mADC[i]+1);   // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)
2131     numAlu.AssignDouble((Double_t)mADC[i] + 1.5);      // numAlu has enough pre-comma place for that; correct trafo, best values
2132     mOffset[i]  = mOffset[i] + numAlu.GetValue();      // transform offset to a coord.system relative to chip; +1 to avoid neg. values
2133
2134     // up to here: adc-mapping according to TRAP manual and in line with pad-col mapping
2135     // reverse adc-counting to be again in line with the online mapping
2136     mADC[i]     = fNADC - 4 - mADC[i];                 // fNADC-4-mADC[i]: 0..17; remapping necessary;
2137     mADC[i]     = mADC[i] + 2;
2138     // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
2139   }
2140
2141   // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected;
2142   // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
2143
2144   // transform parameters to the local coordinate-system of a stack (used by GTU)
2145   AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
2146
2147   Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
2148   //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
2149
2150   // difference between width of inner and outer pads of a row is not accounted for;
2151
2152   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
2153   Double_t eCharge  = 0.3;                           // unit charge in (GeV/c)/m*T
2154   Double_t ptMin   = 2.3;                            // minimum transverse momentum (GeV/c); to be adjusted(?)
2155
2156   Double_t granularityOffset = 0.160;                // granularity for offset in mm
2157   Double_t granularitySlope  = 0.140;                // granularity for slope  in mm
2158
2159   // get the coordinates in SM-system; parameters:
2160
2161   Double_t zPos       =  (padPlane->GetRowPos(fRow))*10.0;  // z-position of the MCM; fRow is counted on a chamber; SM consists of 5
2162   // zPos is position of pad-borders;
2163   Double_t zOffset = 0.0;
2164   if ( fRow == 0 || fRow == 15 ) {
2165       zOffset = padPlane->GetLengthOPad();
2166   }
2167   else {
2168       zOffset = padPlane->GetLengthIPad();
2169   }
2170   zOffset = (-1.0) * zOffset/2.0;
2171   // turn zPos to be z-coordinate at middle of pad-row
2172   zPos = zPos + zOffset;
2173
2174
2175   Double_t xPos       =  0.0;                               // x-position of the upper border of the drift-chamber of actual layer
2176   Int_t    icol       =  0;                                 // column-number of adc-channel
2177   Double_t yPos[4];                                         // y-position of the pad to which ADC is connected
2178   Double_t dx         = 30.0;                               // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
2179   Double_t freqSample = fFeeParam->GetSamplingFrequency();  // retrieve the sampling frequency (10.019750 MHz)
2180   Double_t vdrift     = fCal->GetVdriftAverage(fChaId);     // averaged drift velocity for this detector (1.500000 cm/us)
2181   Int_t    nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
2182   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))
2183   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
2184   if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed
2185   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
2186   //Double_t lorAngle   =  7.0;                             // Lorentz-angle in degrees
2187   Double_t tiltAngle  = padPlane->GetTiltingAngle();        // sign-respecting tilting angle of pads in actual layer
2188   Double_t tiltTan    = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
2189   //Double_t lorTan     = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
2190
2191   Double_t alphaMax[4];                            // maximum deflection from the direction to the primary vertex; granularity of hit pads
2192   Double_t slopeMin[4];                            // local limits for the deflection
2193   Double_t slopeMax[4];
2194   Int_t   mslopeMin[4];                            // in granularity units; to be compared to mSlope[i]
2195   Int_t   mslopeMax[4];
2196
2197
2198   // x coord. of upper side of drift chambers in local SM-system (in mm)
2199   // obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
2200   // the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
2201   switch(fLayer)
2202     {
2203     case 0:
2204       xPos = 3003.0;
2205       break;
2206     case 1:
2207       xPos = 3129.0;
2208       break;
2209     case 2:
2210       xPos = 3255.0;
2211       break;
2212     case 3:
2213       xPos = 3381.0;
2214       break;
2215     case 4:
2216       xPos = 3507.0;
2217       break;
2218     case 5:
2219       xPos = 3633.0;
2220       break;
2221     }
2222
2223   // calculation of offset-correction n:
2224
2225   Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));
2226
2227   nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
2228   if (nCorrectOffset < 0) {
2229     numAlu.AssignInt(-nCorrectOffset);
2230     numAlu.SetSign(-1);
2231   }
2232   else {
2233     numAlu.AssignInt(nCorrectOffset);
2234     numAlu.SetSign(1);
2235   }
2236   nCorrectOffset = numAlu.GetSignedValue();
2237
2238   // the Lorentz correction to the offset
2239   Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
2240
2241
2242   lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
2243
2244   if(lorCorrectOffset < 0) {
2245       numAlu.AssignDouble(-lorCorrectOffset);
2246       numAlu.SetSign(-1);
2247   }
2248   else{
2249       numAlu.AssignDouble(lorCorrectOffset);
2250       numAlu.SetSign(1);
2251   }
2252
2253   Int_t mlorCorrectOffset = numAlu.GetSignedValue();
2254
2255
2256   Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
2257
2258   // calculation of slope-correction
2259
2260   // this is only true for tracks coming (approx.) from primary vertex
2261   // everything is evaluated for a tracklet covering the whole drift chamber
2262   Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
2263   // Double_t cCorrectSlope =  zPos/xPos*dx*tiltTan/granularitySlope;
2264   // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
2265
2266   if (cCorrectSlope < 0) {
2267       numAlu.AssignDouble(-cCorrectSlope);
2268       numAlu.SetSign(-1);
2269   }
2270   else {
2271       numAlu.AssignDouble(cCorrectSlope);
2272       numAlu.SetSign(1);
2273   }
2274   cCorrectSlope = numAlu.GetSignedValue();
2275
2276   // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
2277   // different pad-width of outer pads of a pad-plane not taken into account
2278   // 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)
2279   // 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
2280
2281
2282   Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope;  // >= 0.0
2283
2284   AliTRDtrapAlu correctAlu;
2285   correctAlu.Init(20,8);
2286
2287   AliTRDtrapAlu offsetAlu;
2288   offsetAlu.Init(13,0,-0x1000,0x0FFF);          // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
2289
2290   AliTRDtrapAlu slopeAlu;
2291   slopeAlu.Init(7,0,-0x40,0x3F);                // 7 bit-word;  2-complement (1 sign-bit);
2292
2293   for (Int_t i = 0; i<4; i++) {
2294
2295     if (mADC[i] == -1) continue;
2296
2297     icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
2298     yPos[i]   = (padPlane->GetColPos(icol))*10.0;
2299
2300
2301     // offset:
2302
2303     correctAlu.AssignDouble(mCorrectOffset);     // done because max. accuracy is 8 bit
2304     mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
2305     mOffset[i]  = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset));
2306     //mOffset[i]  = mOffset[i]*(-1);                   // adjust to direction of y-axes in online simulation
2307
2308     if (mOffset[i] < 0) {
2309       numAlu.AssignFormatted(-mOffset[i]);
2310       numAlu.SetSign(-1);
2311     }
2312     else {
2313       numAlu.AssignFormatted(mOffset[i]);
2314       numAlu.SetSign(1);
2315     }
2316
2317     offsetAlu = numAlu;
2318     mOffset[i]  = offsetAlu.GetSignedValue();
2319
2320
2321     // slope:
2322
2323     correctAlu.AssignDouble(mCorrectSlope);
2324     mCorrectSlope = correctAlu.GetValueWhole();
2325
2326     mSlope[i]   = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
2327
2328     if (mSlope[i] < 0) {
2329       numAlu.AssignFormatted(-mSlope[i]);
2330       numAlu.SetSign(-1);
2331     }
2332     else {
2333       numAlu.AssignFormatted(mSlope[i]);
2334       numAlu.SetSign(1);
2335     }
2336
2337     slopeAlu  = numAlu;     // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
2338     mSlope[i]   = slopeAlu.GetSignedValue();
2339
2340     // local (LTU) limits for the deflection
2341     // ATan returns angles in radian
2342     alphaMax[i]  = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
2343     slopeMin[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
2344     slopeMax[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
2345
2346     if (slopeMin[i] < 0) {
2347       slopeAlu.AssignDouble(-slopeMin[i]);
2348       slopeAlu.SetSign(-1);
2349     }
2350     else {
2351       slopeAlu.AssignDouble(slopeMin[i]);
2352       slopeAlu.SetSign(1);
2353     }
2354     mslopeMin[i] = slopeAlu.GetSignedValue();  // the borders should lie inside the range of mSlope -> usage of slopeAlu again
2355
2356     if (slopeMax[i] < 0) {
2357       slopeAlu.AssignDouble(-slopeMax[i]);
2358       slopeAlu.SetSign(-1);
2359     }
2360     else {
2361       slopeAlu.AssignDouble(slopeMax[i]);
2362       slopeAlu.SetSign(1);
2363     }
2364     mslopeMax[i] = slopeAlu.GetSignedValue();
2365   }
2366
2367   // suppress submission of tracks with low stiffness
2368   // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber)
2369
2370   // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
2371   // up to now they are sorted according to maximum hit sum
2372   // is the sorting really done in the TRAP-chip?
2373
2374   Int_t order[4] = {-1,-1,-1,-1};
2375   Int_t wordnr = 0;   // number of tracklet-words
2376
2377   for(Int_t j = 0; j < fMaxTracklets; j++) {
2378       //if( mADC[j] == -1) continue;
2379       if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
2380       wordnr++;
2381       if( wordnr-1 == 0) {
2382           order[0] = j;
2383           continue;
2384       }
2385       // wordnr-1>0, wordnr-1<4
2386       order[wordnr-1] = j;
2387       for( Int_t k = 0; k < wordnr-1; k++) {
2388           if( mOffset[j] < mOffset[order[k]] ) {
2389               for( Int_t l = wordnr-1; l > k; l-- ) {
2390                   order[l] = order[l-1];
2391               }
2392               order[k] = j;
2393               break;
2394           }
2395
2396       }
2397   }
2398
2399   // fill the bit-words in ascending order and without gaps
2400   UInt_t bitWord[4] = {0,0,0,0};                 // attention: unsigned int to have real 32 bits (no 2-complement)
2401   for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
2402       //Bool_t rem1 = kTRUE;
2403
2404     Int_t i = order[j];
2405     //bit-word is 2-complement and therefore without sign
2406     bitWord[j] =   1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
2407     //printf("\n");
2408     UInt_t shift  = 0;
2409     UInt_t shift2 = 0;
2410
2411
2412
2413
2414     /*printf("mean charge: %d\n",mMeanCharge[i]);
2415     printf("row: %d\n",fRow);
2416     printf("slope: %d\n",mSlope[i]);
2417     printf("pad position: %d\n",mOffset[i]);
2418     printf("channel: %d\n",mADC[i]);*/
2419
2420     // electron probability (currently not implemented; the mean charge is just scaled)
2421     shift = (UInt_t)mMeanCharge[i];
2422     for(Int_t iBit = 0; iBit < 8; iBit++) {
2423       bitWord[j]  = bitWord[j]<<1;
2424       bitWord[j] |= (shift>>(7-iBit))&1;
2425       //printf("0");
2426     }
2427
2428     // pad row
2429     shift = (UInt_t)fRow;
2430     for(Int_t iBit = 0; iBit < 4; iBit++) {
2431       bitWord[j]  = bitWord[j]<<1;
2432       bitWord[j] |= (shift>>(3-iBit))&1;
2433       //printf("%d", (fRow>>(3-iBit))&1);
2434     }
2435
2436     // deflection length
2437     if(mSlope[i] < 0) {
2438         shift = (UInt_t)(-mSlope[i]);
2439         // shift2 is 2-complement of shift
2440         shift2 = 1;
2441         for(Int_t iBit = 1; iBit < 7; iBit++) {
2442             shift2  = shift2<<1;
2443             shift2 |= (1- (((shift)>>(6-iBit))&1) );
2444             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2445         }
2446         shift2 = shift2 + 1;
2447         //printf("1");
2448         for(Int_t iBit = 0; iBit < 7; iBit++) {
2449             bitWord[j]  = bitWord[j]<<1;
2450             bitWord[j] |= (shift2>>(6-iBit))&1;
2451             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2452         }
2453     }
2454     else {
2455         shift = (UInt_t)(mSlope[i]);
2456         bitWord[j]  = bitWord[j]<<1;
2457         bitWord[j]   |= 0;
2458         //printf("0");
2459         for(Int_t iBit = 1; iBit < 7; iBit++) {
2460             bitWord[j]  = bitWord[j]<<1;
2461             bitWord[j] |= (shift>>(6-iBit))&1;
2462             //printf("%d",(mSlope[i]>>(6-iBit))&1);
2463         }
2464     }
2465
2466     // pad position
2467     if(mOffset[i] < 0) {
2468         shift = (UInt_t)(-mOffset[i]);
2469         shift2 = 1;
2470         for(Int_t iBit = 1; iBit < 13; iBit++) {
2471             shift2  = shift2<<1;
2472             shift2 |= (1-(((shift)>>(12-iBit))&1));
2473             //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
2474         }
2475         shift2 = shift2 + 1;
2476         //printf("1");
2477         for(Int_t iBit = 0; iBit < 13; iBit++) {
2478             bitWord[j]  = bitWord[j]<<1;
2479             bitWord[j] |= (shift2>>(12-iBit))&1;
2480             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2481         }
2482     }
2483     else {
2484         shift = (UInt_t)mOffset[i];
2485         bitWord[j] = bitWord[j]<<1;
2486         bitWord[j]   |= 0;
2487         //printf("0");
2488         for(Int_t iBit = 1; iBit < 13; iBit++) {
2489             bitWord[j]  = bitWord[j]<<1;
2490             bitWord[j] |= (shift>>(12-iBit))&1;
2491             //printf("%d",(mOffset[i]>>(12-iBit))&1);
2492         }
2493     }
2494
2495
2496
2497     //printf("bitWord: %u\n",bitWord[j]);
2498     //printf("adc: %d\n",mADC[i]);
2499     fMCMT[j] = bitWord[j];
2500   }
2501
2502   //printf("\n");
2503
2504
2505   delete [] qsum;
2506   delete [] ieffped;
2507
2508   delete [] X;
2509   delete [] XX;
2510   delete [] Y;
2511   delete [] YY;
2512   delete [] XY;
2513   delete [] N;
2514   delete [] QT0;
2515   delete [] QT1;
2516
2517   delete [] hitSum;
2518   delete [] selectPair;
2519
2520   delete padPlane;
2521
2522 //if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
2523
2524   if (!fFeeParam->GetMCTrackletOutput())
2525       return;
2526
2527   AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
2528   AliLog::SetFileOutput("../log/tracklet.log");
2529
2530   // testing for wordnr in order to speed up the simulation
2531   if (wordnr == 0)
2532     return;
2533
2534   UInt_t        *trackletWord = new UInt_t[fMaxTracklets];
2535   Int_t         *adcChannel   = new Int_t[fMaxTracklets];
2536   Int_t         *trackRef     = new Int_t[fMaxTracklets];
2537
2538   Int_t u = 0;
2539
2540   AliTRDdigitsManager *digman = new AliTRDdigitsManager();
2541   digman->ReadDigits(AliRunLoader::GetRunLoader()->GetLoader("TRDLoader")->TreeD());
2542   digman->SetUseDictionaries(kTRUE);
2543   AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
2544
2545   for (Int_t j = 0; j < fMaxTracklets; j++) {
2546       Int_t i = order[j];
2547       trackletWord[j] = 0;
2548       adcChannel[j] = -1;
2549       if (bitWord[j]!=0) {
2550           trackletWord[u] = bitWord[j];
2551           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
2552
2553 // Finding label of MC track
2554           TH1F *hTrkRef = new TH1F("trackref", "trackref", 100000, 0, 100000);
2555           Int_t track[3];
2556           Int_t padcol = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u]);
2557           Int_t padcol_ngb = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u] - 1);
2558           Int_t padrow = 4 * (fRobPos / 2) + fMcmPos / 4;
2559           Int_t det = 30 * fSector + 6 * fStack + fLayer;
2560           for(Int_t iTimebin = feeParam->GetLinearFitStart(); iTimebin < feeParam->GetLinearFitEnd(); iTimebin++) {
2561               track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2562               track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2563               track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2564               hTrkRef->Fill(track[0]);
2565               if (track[1] != track[0] && track[1] != -1)
2566                   hTrkRef->Fill(track[1]);
2567               if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2568                   hTrkRef->Fill(track[2]);
2569               if (padcol_ngb >= 0) {
2570                   track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2571                   track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2572                   track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2573                   hTrkRef->Fill(track[0]);
2574                   if (track[1] != track[0] && track[1] != -1)
2575                       hTrkRef->Fill(track[1]);
2576                   if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2577                       hTrkRef->Fill(track[2]);
2578               }
2579           }
2580           trackRef[u] = hTrkRef->GetMaximumBin() - 1;
2581           delete hTrkRef;
2582           u = u + 1;
2583       }
2584   }
2585
2586   AliDataLoader *dl = AliRunLoader::GetRunLoader()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
2587   if (!dl) {
2588     AliError("Could not get the tracklets data loader!");
2589   }
2590   else {
2591     TTree *trackletTree = dl->Tree();
2592     if (!trackletTree)
2593       dl->MakeTree();
2594     trackletTree = dl->Tree();
2595
2596    AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM();
2597    TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
2598    if (!trkbranch)
2599        trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
2600     trkbranch->SetAddress(&trkl);
2601
2602     for (Int_t iTracklet = 0; iTracklet < fMaxTracklets; iTracklet++) {
2603         if (trackletWord[iTracklet] == 0)
2604             continue;
2605         trkl->SetTrackletWord(trackletWord[iTracklet]);
2606         trkl->SetDetector(30*fSector + 6*fStack + fLayer);
2607         trkl->SetROB(fRobPos);
2608         trkl->SetMCM(fMcmPos);
2609         trkl->SetLabel(trackRef[iTracklet]);
2610         trackletTree->Fill();
2611     }
2612     delete trkl;
2613     dl->WriteData("OVERWRITE");
2614   }
2615
2616   delete [] trackletWord;
2617   delete [] adcChannel;
2618   delete [] trackRef;
2619   delete digman;
2620
2621   // to be done:
2622   // error measure for quality of fit (not necessarily needed for the trigger)
2623   // cluster quality threshold (not yet set)
2624   // electron probability
2625 }
2626 //_____________________________________________________________________________________
2627 void AliTRDmcmSim::GeneratefZSM1Dim()
2628 {
2629   //
2630   // Generate the array fZSM1Dim necessary
2631   // for the method ProduceRawStream
2632   //
2633
2634   // Fill the mapping
2635   // Supressed zeros indicated by -1 in digits array
2636   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ )
2637     {
2638       for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2639         {
2640
2641           if(fADCF[iadc][it]==-1)  // If is a supressed value
2642             {
2643               fZSM[iadc][it]=1;
2644             }
2645           else                    // Not suppressed
2646             {
2647               fZSM[iadc][it]=0;
2648             }
2649         }
2650     }
2651
2652   // Make the 1 dim projection
2653   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2654     {
2655       for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2656         {
2657           fZSM1Dim[iadc] &= fZSM[iadc][it];
2658         }
2659     }
2660 }
2661 //_______________________________________________________________________________________
2662 void AliTRDmcmSim::CopyArrays()
2663 {
2664   //
2665   // Initialize filtered data array with raw data
2666   // Method added for internal consistency
2667   //
2668
2669   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2670     {
2671       for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2672         {
2673           fADCF[iadc][it] = fADCR[iadc][it];
2674         }
2675     }
2676 }
2677 //_______________________________________________________________________________________
2678 void AliTRDmcmSim::StartfastZS(Int_t pads, Int_t timebins)
2679 {
2680   //
2681   // Initialize just the necessary elements to perform
2682   // the zero suppression in the digitizer
2683   //
2684
2685   fFeeParam  = AliTRDfeeParam::Instance();
2686   fSimParam  = AliTRDSimParam::Instance();
2687   fNADC      = pads;
2688   fNTimeBin  = timebins;
2689
2690   if( fADCR == NULL )
2691     {
2692       fADCR    = new Int_t *[fNADC];
2693       fADCF    = new Int_t *[fNADC];
2694       fADCT    = new Int_t *[fNADC];
2695       fZSM     = new Int_t *[fNADC];
2696       fZSM1Dim = new Int_t  [fNADC];
2697     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2698       {
2699         fADCR[iadc] = new Int_t[fNTimeBin];
2700         fADCF[iadc] = new Int_t[fNTimeBin];
2701         fADCT[iadc] = new Int_t[fNTimeBin];
2702         fZSM [iadc] = new Int_t[fNTimeBin];
2703       }
2704     }
2705
2706   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2707     {
2708       for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2709         {
2710           fADCR[iadc][it] =  0;
2711           fADCF[iadc][it] =  0;
2712           fADCT[iadc][it] = -1;
2713           fZSM [iadc][it] =  1;
2714         }
2715       fZSM1Dim[iadc] = 1;
2716     }
2717
2718   fInitialized = kTRUE;
2719 }
2720 //_______________________________________________________________________________________
2721 void AliTRDmcmSim::FlagDigitsArray(AliTRDarrayADC *tempdigs, Int_t valrow)
2722 {
2723   //
2724   // Modify the digits array to flag suppressed values
2725   //
2726
2727   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ )
2728     {
2729       for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2730         {
2731           if(fZSM[iadc][it]==1)
2732             {
2733               tempdigs->SetData(valrow,iadc,it,-1);
2734             }
2735         }
2736     }
2737 }
2738 //_______________________________________________________________________________________
2739 void AliTRDmcmSim::RestoreZeros()
2740 {
2741   //
2742   // Restore the zero-suppressed values (set as -1) to the value 0
2743   //
2744
2745   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ )
2746     {
2747       for( Int_t it = 0 ; it < fNTimeBin ; it++ )
2748         {
2749
2750           if(fADCF[iadc][it]==-1)  //if is a supressed zero, reset to zero
2751             {
2752               fADCF[iadc][it]=0;
2753               fADCR[iadc][it]=0;
2754             }
2755         }
2756     }
2757
2758 }
2759