TRD/AliTRDmcmSim.cxx

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