TRD/AliTRDmcmSim.cxx

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