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