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 = kFALSE ) // 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 // in order to use the Tracklets, please first
1259 // -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
1260 // -- set AliTRDfeeParam::fgkTFtype   to 3, in order to use the new tail cancellation filter
1261 //    currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
1262 // -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with           information about the Tracklet position (MCM, channel number)
1263
1264 // The code is designed such that the less possible calculations with AliTRDtrapAlu class-objects are performed; whenever possible calculations are done with doubles or integers and the results are transformed into the right format
1265
1266 void AliTRDmcmSim::Tracklet(){
1267     // tracklet calculation
1268     // if you use this code after 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-manual-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 go to AliTRDfeeParam once they are defined; then they 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   AliTRDtrapAlu oneAlu;
1425   oneAlu.Init(1,8);
1426
1427
1428   AliTRDtrapAlu inverseNAlu;
1429   inverseNAlu.Init(1,8);        // simulates the LUT for 1/N
1430   AliTRDtrapAlu MeanChargeAlu;  // mean charge in ADC counts
1431   MeanChargeAlu.Init(8,0);
1432   AliTRDtrapAlu TotalChargeAlu;
1433   TotalChargeAlu.Init(17,8);
1434   //nr of post comma bits should be the same for inverseN and TotalCharge
1435
1436
1437   SetPosLUT();                    // initialize the position correction LUT for this MCM;
1438
1439
1440   // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
1441   Int_t *X   = new Int_t[fNADC-2];
1442   Int_t *XX  = new Int_t[fNADC-2];
1443   Int_t *Y   = new Int_t[fNADC-2];
1444   Int_t *YY  = new Int_t[fNADC-2];
1445   Int_t *XY  = new Int_t[fNADC-2];
1446   Int_t *N   = new Int_t[fNADC-2];
1447   Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
1448   Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
1449
1450   for (Int_t iCol = 0; iCol < fNADC-2; iCol++) {
1451
1452       // initialize fit-sums
1453       X[iCol]   = 0;
1454       XX[iCol]  = 0;
1455       Y[iCol]   = 0;
1456       YY[iCol]  = 0;
1457       XY[iCol]  = 0;
1458       N[iCol]   = 0;
1459       QT0[iCol] = 0;
1460       QT1[iCol] = 0;
1461   }
1462
1463
1464   filPed.Init(7,2);                         // convert filtered pedestal into 7+2Bits
1465
1466   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1467
1468     if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no?
1469
1470     // reset
1471     Int_t portChannel[4]   = {-1,-1,-1,-1};
1472     Int_t clusterCharge[4] = {0,0,0,0};
1473     Int_t leftCharge[4]    = {0,0,0,0};
1474     Int_t centerCharge[4]  = {0,0,0,0};
1475     Int_t rightCharge[4]   = {0,0,0,0};
1476
1477     Int_t mark = 0;
1478
1479     filPed.AssignFormatted(ieffped[iT]);   // no size-checking when using AssignFormatted; ieffped>=0
1480     filPed = filPed;                       // this checks the size
1481
1482     ieffped[iT] = filPed.GetValue();
1483
1484     for(Int_t i = 0; i<7; i++){
1485       selection[i].next       = NULL;
1486       selection[i].iadc       =   -1;     // value of -1: invalid adc
1487       selection[i].value      =    0;
1488
1489     }
1490     // selection[0] is starting list-element; just for pointing
1491
1492     // loop over inner adc's
1493     for (Int_t iCol = 1; iCol < fNADC-1; iCol++) {
1494
1495       Int_t left   = fADCT[iCol-1][iT];
1496       Int_t center = fADCT[iCol][iT];
1497       Int_t right  = fADCT[iCol+1][iT];
1498
1499       Int_t sum = left + center + right;            // cluster charge sum
1500       qsumAlu.AssignFormatted(sum);
1501       qsumAlu = qsumAlu;                        // size-checking; redundant
1502
1503       qsum[iCol] = qsumAlu.GetValue();
1504
1505       //hit detection and masking
1506       if(center>=left){
1507         if(center>right){
1508           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
1509             mark |= 1;                              // marker
1510           }
1511         }
1512       }
1513       mark = mark<<1;
1514     }
1515     mark = mark>>1;
1516
1517
1518     // get selection of 6 adc's and sort,starting with greatest values
1519
1520     //read three from right side and sort (primitive sorting algorithm)
1521     Int_t i = 0; // adc number
1522     Int_t j = 1; // selection number
1523     while(i<fNADC-2 && j<=3){
1524       i = i + 1;
1525       if((mark>>(i-1)) & 1 == 1) {
1526         selection[j].iadc  = fNADC-1-i;
1527         selection[j].value = qsum[fNADC-1-i]>>6;   // for hit-selection only the first 8 out of the 14 Bits are used for comparison
1528
1529         // insert into sorted list
1530         listLeft = &selection[0];
1531         list = listLeft->next;
1532
1533         if(list!=NULL) {
1534           while((list->next != NULL) && (selection[j].value <= list->value)){
1535             listLeft = list;
1536             list = list->next;
1537           }
1538
1539           if(selection[j].value<=list->value){
1540             selection[j].next = list->next;
1541             list->next = &selection[j];
1542           }
1543           else {
1544             listLeft->next = &selection[j];
1545             selection[j].next = list;
1546           }
1547         }
1548         else{
1549           listLeft->next = &selection[j];
1550           selection[j].next = list;
1551         }
1552
1553         j = j + 1;
1554       }
1555     }
1556
1557
1558     // read three from left side
1559     Int_t k = fNADC-2;
1560     while(k>i && j<=6) {
1561       if((mark>>(k-1)) & 1 == 1) {
1562         selection[j].iadc  = fNADC-1-k;
1563         selection[j].value = qsum[fNADC-1-k]>>6;
1564
1565         listLeft = &selection[0];
1566         list = listLeft->next;
1567
1568         if(list!=NULL){
1569           while((list->next != NULL) && (selection[j].value <= list->value)){
1570             listLeft = list;
1571             list = list->next;
1572           }
1573
1574           if(selection[j].value<=list->value){
1575             selection[j].next = list->next;
1576             list->next = &selection[j];
1577           }
1578           else {
1579             listLeft->next = &selection[j];
1580             selection[j].next = list;
1581           }
1582         }
1583         else{
1584           listLeft->next = &selection[j];
1585           selection[j].next = list;
1586         }
1587
1588         j = j + 1;
1589       }
1590       k = k - 1;
1591     }
1592
1593     // get the four with greatest charge-sum
1594     list = &selection[0];
1595     for(i = 0; i<4; i++){
1596       if(list->next == NULL) continue;
1597       list = list->next;
1598       if(list->iadc == -1) continue;
1599       Int_t adc = list->iadc;                              // channel number with selected hit
1600
1601       // the following arrays contain the four chosen channels in 1 time-bin
1602       portChannel[i]   = adc;
1603       clusterCharge[i] = qsum[adc];
1604       leftCharge[i]    = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
1605       centerCharge[i]  = fADCT[adc][iT] - ieffped[iT];
1606       rightCharge[i]   = fADCT[adc+1][iT] - ieffped[iT];
1607     }
1608
1609     // arithmetic unit
1610
1611     // cluster verification
1612     if(!bVBY){
1613       for(i = 0; i<4; i++){
1614         Int_t lr = leftCharge[i]*rightCharge[i]*1024;
1615         Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
1616         if (lr>=cc){
1617           portChannel[i]   = -1;                                 // set to invalid address
1618           clusterCharge[i] = 0;
1619         }
1620       }
1621     }
1622
1623     // fit-sums of valid channels
1624     // local hit position
1625     for(i = 0; i<4; i++){
1626       if (centerCharge[i] ==  0) {
1627         portChannel[i] = -1;
1628       }// prevent division by 0
1629
1630       if (portChannel[i]  == -1) continue;
1631
1632       Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
1633
1634       Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
1635       dCOG = (sign == kFALSE) ? dCOG : -dCOG;     // AssignDouble doesn't allow for signed doubles
1636       dCOGAlu.AssignDouble(dCOG);
1637       Int_t iLUTpos = dCOGAlu.GetValue();       // steers position in LUT
1638
1639       dCOG = dCOG/2;
1640       yrawAlu.AssignDouble(dCOG);
1641       Int_t iCOG = yrawAlu.GetValue();
1642       Int_t y = iCOG + fPosLUT[iLUTpos % 128];    // local position in pad-units
1643       yrawAlu.AssignFormatted(y);               // 0<y<1
1644       yAlu  = yrawAlu;                        // convert to 16 past-comma bits
1645
1646       if(sign == kTRUE) yAlu.SetSign(-1);       // buffer width of 9 bits; sign on real (not estimated) position
1647       xAlu.AssignInt(iT);                       // buffer width of 5 bits
1648
1649
1650       xxAlu = xAlu * xAlu;                  // buffer width of 10 bits -> fulfilled by x*x
1651
1652       yyAlu = yAlu * yAlu;                  // buffer width of 16 bits
1653
1654       xyAlu = xAlu * yAlu;                  // buffer width of 14 bits
1655
1656       Int_t adc = portChannel[i]-1;              // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
1657
1658       // calculate fit-sums recursively
1659       // interpretation of their bit-length is given as comment
1660
1661       // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
1662
1663       XAlu.AssignFormatted(X[adc]);
1664       XAlu = XAlu + xAlu;                   // buffer width of 9 bits
1665       X[adc] = XAlu.GetValue();
1666
1667       XXAlu.AssignFormatted(XX[adc]);
1668       XXAlu = XXAlu + xxAlu;                // buffer width of 14 bits
1669       XX[adc] = XXAlu.GetValue();
1670
1671       if (Y[adc] < 0) {
1672         YAlu.AssignFormatted(-Y[adc]);          // make sure that only positive values are assigned; sign-setting must be done by hand
1673         YAlu.SetSign(-1);
1674       }
1675       else {
1676         YAlu.AssignFormatted(Y[adc]);
1677         YAlu.SetSign(1);
1678       }
1679
1680       YAlu = YAlu + yAlu;                   // buffer width of 14 bits (8 past-comma);
1681       Y[adc] = YAlu.GetSignedValue();
1682
1683       YYAlu.AssignFormatted(YY[adc]);
1684       YYAlu = YYAlu + yyAlu;                // buffer width of 21 bits (16 past-comma)
1685       YY[adc] = YYAlu.GetValue();
1686
1687       if (XY[adc] < 0) {
1688         XYAlu.AssignFormatted(-XY[adc]);
1689         XYAlu.SetSign(-1);
1690       }
1691       else {
1692         XYAlu.AssignFormatted(XY[adc]);
1693         XYAlu.SetSign(1);
1694       }
1695
1696       XYAlu = XYAlu + xyAlu;                // buffer allows 17 bits (8 past-comma)
1697       XY[adc] = XYAlu.GetSignedValue();
1698
1699       N[adc]  = N[adc] + 1;
1700
1701
1702       // accumulated charge
1703       qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
1704       qtruncAlu = qsumAlu;
1705
1706       if(iT>=tQS0 && iT<=tQE0){
1707         QT0Alu.AssignFormatted(QT0[adc]);
1708         QT0Alu = QT0Alu + qtruncAlu;
1709         QT0[adc] = QT0Alu.GetValue();
1710         //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
1711       }
1712
1713       if(iT>=tQS1 && iT<=tQE1){
1714         QT1Alu.AssignFormatted(QT1[adc]);
1715         QT1Alu = QT1Alu + qtruncAlu;
1716         QT1[adc] = QT1Alu.GetValue();
1717         //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
1718       }
1719     }// i
1720
1721     // remapping is done!!
1722
1723   }//iT
1724
1725
1726
1727   // tracklet-assembly
1728
1729   // put into AliTRDfeeParam and take care that values are in proper range
1730   const Int_t cTCL = 1;      // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??)
1731   const Int_t cTCT = 8;      // joint number of hits;     8<=TCT<=31; note that according to TRAP manual this number cannot be lower than 8; however it should be adjustable to the number of hits in the fit time range (40%)
1732
1733   Int_t mPair   = 0;         // marker for possible tracklet pairs
1734   Int_t* hitSum = new Int_t[fNADC-3];
1735   // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18;
1736
1737   // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
1738   for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
1739     hitSum[iCol] = N[iCol] + N[iCol+1];
1740     if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
1741         mPair |= 1;         // mark as possible channel-pair
1742
1743     }
1744     mPair = mPair<<1;
1745   }
1746   mPair = mPair>>1;
1747
1748   List_t* selectPair = new List_t[fNADC-2];      // list with 18 elements (0..18) containing the left channel-nr and hit sums
1749                                                  // selectPair[18] is starting list-element just for pointing
1750   for(Int_t k = 0; k<fNADC-2; k++){
1751       selectPair[k].next       = NULL;
1752       selectPair[k].iadc       =   -1;           // invalid adc
1753       selectPair[k].value      =    0;
1754
1755     }
1756
1757  list = NULL;
1758  listLeft = NULL;
1759
1760   // read marker and sort according to hit-sum
1761
1762   Int_t adcL  = 0;            // left adc-channel-number (remapped)
1763   Int_t selNr = 0;            // current number in list
1764
1765   // insert marked channels into list and sort according to hit-sum
1766   while(adcL < fNADC-3 && selNr < fNADC-3){
1767
1768     if((mPair>>((fNADC-4)-(adcL))) & 1 == 1) {
1769       selectPair[selNr].iadc  = adcL;
1770       selectPair[selNr].value = hitSum[adcL];
1771
1772       listLeft = &selectPair[fNADC-3];
1773       list = listLeft->next;
1774
1775       if(list!=NULL) {
1776         while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
1777           listLeft = list;
1778           list = list->next;
1779         }
1780
1781         if(selectPair[selNr].value <= list->value){
1782           selectPair[selNr].next = list->next;
1783           list->next = &selectPair[selNr];
1784         }
1785         else {
1786           listLeft->next = &selectPair[selNr];
1787           selectPair[selNr].next = list;
1788         }
1789
1790       }
1791       else{
1792         listLeft->next = &selectPair[selNr];
1793         selectPair[selNr].next = list;
1794       }
1795
1796       selNr = selNr + 1;
1797     }
1798     adcL = adcL + 1;
1799   }
1800
1801   //select up to 4 channels with maximum number of hits
1802   Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
1803   Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
1804   list = &selectPair[fNADC-3];
1805
1806   for (Int_t i = 0; i<4; i++) {
1807     if(list->next == NULL) continue;
1808     list = list->next;
1809     if(list->iadc == -1) continue;
1810     lpairChannel[i] = list->iadc;        // channel number with selected hit
1811     rpairChannel[i] = lpairChannel[i]+1;
1812   }
1813
1814   // avoid submission of double tracklets
1815   for (Int_t i = 3; i>0; i--) {
1816     for (Int_t j = i-1; j>-1; j--) {
1817       if(lpairChannel[i] == rpairChannel[j]) {
1818         lpairChannel[i] = -1;
1819         rpairChannel[i] = -1;
1820         break;
1821       }
1822       /* if(rpairChannel[i] == lpairChannel[j]) {
1823         lpairChannel[i] = -1;
1824         rpairChannel[i] = -1;
1825         break;
1826         }*/
1827     }
1828   }
1829
1830   // merging of the fit-sums of the remainig channels
1831   // assume same data-word-width as for fit-sums for 1 channel
1832   // relative scales!
1833   Int_t mADC[4];
1834   Int_t mN[4];
1835   Int_t mQT0[4];
1836   Int_t mQT1[4];
1837   Int_t mX[4];
1838   Int_t mXX[4];
1839   Int_t mY[4];
1840   Int_t mYY[4];
1841   Int_t mXY[4];
1842   Int_t mOffset[4];
1843   Int_t mSlope[4];
1844   Int_t mMeanCharge[4];
1845   Int_t inverseN = 0;
1846   Double_t invN = 0;
1847   Int_t one = 0;
1848
1849   for (Int_t i = 0; i<4; i++){
1850     mADC[i] = -1;                        // set to invalid number
1851     mN[i]   =  0;
1852     mQT0[i] =  0;
1853     mQT1[i] =  0;
1854     mX[i]   =  0;
1855     mXX[i]  =  0;
1856     mY[i]   =  0;
1857     mYY[i]  =  0;
1858     mXY[i]  =  0;
1859     mOffset[i] = 0;
1860     mSlope[i]  = 0;
1861     mMeanCharge[i] = 0;
1862   }
1863
1864   oneAlu.AssignInt(1);
1865   one = oneAlu.GetValue();              // one with 8 past comma bits
1866
1867   for (Int_t i = 0; i<4; i++){
1868
1869
1870     mADC[i] = lpairChannel[i];          // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
1871                                         // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
1872     Int_t madc = mADC[i];
1873     if (madc == -1) continue;
1874
1875     YAlu.AssignInt(N[rpairChannel[i]]);
1876     Int_t wpad  = YAlu.GetValue();       // enlarge hit counter of right channel by 8 past-comma bits; YAlu can have 5 pre-comma bits (values up to 63); hit counter<=nr of time bins (24)
1877
1878     mN[i]    = hitSum[madc];
1879
1880     // 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
1881     if (N[madc+1] == 0) {
1882         mQT0[i] = QT0[madc];
1883         mQT1[i] = QT1[madc];
1884
1885     }
1886     else {
1887
1888         // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
1889
1890         mQT0[i]   = QT0[madc] + QT0[madc+1];
1891         QT0Alu.AssignFormatted(mQT0[i]);
1892         QT0Alu  = QT0Alu;                // size-check
1893         mQT0[i]   = QT0Alu.GetValue();     // write back
1894
1895         mQT1[i]   = QT1[madc] + QT1[madc+1];
1896         QT1Alu.AssignFormatted(mQT1[i]);
1897         QT1Alu  = QT1Alu;
1898         mQT1[i]   = QT1Alu.GetValue();
1899     }
1900
1901     // calculate the mean charge in adc values; later to be replaced by electron likelihood
1902     mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
1903     mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
1904     // 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
1905     invN = 1.0/(mN[i]);
1906     inverseNAlu.AssignDouble(invN);
1907     inverseN = inverseNAlu.GetValue();
1908     mMeanCharge[i] = mMeanCharge[i] * inverseN;  // now to be interpreted with 8 past-comma bits
1909     TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
1910     TotalChargeAlu = TotalChargeAlu;
1911     MeanChargeAlu = TotalChargeAlu;
1912     mMeanCharge[i] = MeanChargeAlu.GetValue();
1913
1914     // this check is not necessary; it is just for efficiency reasons
1915     if (N[madc+1] == 0) {
1916         mX[i]     =   X[madc];
1917         mXX[i]    =  XX[madc];
1918         mY[i]     =   Y[madc];
1919         mXY[i]    =  XY[madc];
1920         mYY[i]    =  YY[madc];
1921     }
1922     else {
1923
1924         mX[i]     =   X[madc] +  X[madc+1];
1925         XAlu.AssignFormatted(mX[i]);
1926         XAlu      = XAlu;
1927         mX[i]     = XAlu.GetValue();
1928
1929         mXX[i]    =  XX[madc] + XX[madc+1];
1930         XXAlu.AssignFormatted(mXX[i]);
1931         XXAlu     = XXAlu;
1932         mXX[i]    = XXAlu.GetValue();
1933
1934
1935         mY[i]     =   Y[madc] + Y[madc+1] + wpad;
1936         if (mY[i] < 0) {
1937             YAlu.AssignFormatted(-mY[i]);
1938             YAlu.SetSign(-1);
1939         }
1940         else {
1941             YAlu.AssignFormatted(mY[i]);
1942             YAlu.SetSign(1);
1943         }
1944         YAlu    = YAlu;
1945         mY[i]     = YAlu.GetSignedValue();
1946
1947         mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*one;    // multiplication by one to maintain the data format
1948
1949         if (mXY[i] < 0) {
1950             XYAlu.AssignFormatted(-mXY[i]);
1951             XYAlu.SetSign(-1);
1952         }
1953         else {
1954             XYAlu.AssignFormatted(mXY[i]);
1955             XYAlu.SetSign(1);
1956         }
1957         XYAlu   = XYAlu;
1958         mXY[i]    = XYAlu.GetSignedValue();
1959
1960         mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
1961         if (mYY[i] < 0) {
1962             YYAlu.AssignFormatted(-mYY[i]);
1963             YYAlu.SetSign(-1);
1964         }
1965         else {
1966             YYAlu.AssignFormatted(mYY[i]);
1967             YYAlu.SetSign(1);
1968         }
1969
1970         YYAlu   = YYAlu;
1971         mYY[i]    = YYAlu.GetSignedValue();
1972     }
1973
1974   }
1975
1976   // calculation of offset and slope from the merged fit-sums;
1977   // YY is needed for some error measure only; still to be done
1978   // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
1979
1980   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1981   // !!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   !!
1982   // !!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.    !!
1983   // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below).                                                                                             !!
1984   // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used.                                                                                          !!
1985   // !!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    !!
1986   // !!valid (in approximation) if tFS is close to the beginning of the drift region.                                                                                                                                                 !!
1987   // !!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                                                                    !!
1988   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1989
1990   // which formats should be chosen?
1991   AliTRDtrapAlu denomAlu;
1992   denomAlu.Init(20,8);
1993   AliTRDtrapAlu numAlu;
1994   numAlu.Init(20,8);
1995   // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
1996   // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
1997
1998   for (Int_t i = 0; i<4; i++) {
1999     if (mADC[i] == -1) continue;
2000
2001     Int_t num0  = (mN[i]*mXX[i]-mX[i]*mX[i]);
2002     if (num0 < 0) {
2003       denomAlu.AssignInt(-num0);    // num0 does not have to be interpreted as having past-comma bits -> AssignInt
2004       denomAlu.SetSign(-1);
2005     }
2006     else {
2007       denomAlu.AssignInt(num0);
2008       denomAlu.SetSign(1);
2009     }
2010
2011     Int_t num1  = mN[i]*mXY[i] - mX[i]*mY[i];
2012     if (num1 < 0) {
2013       numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
2014       numAlu.SetSign(-1);
2015     }
2016     else {
2017       numAlu.AssignFormatted(num1);
2018       numAlu.SetSign(1);
2019     }
2020     numAlu    = numAlu/denomAlu;
2021     mSlope[i]   = numAlu.GetSignedValue();
2022
2023     Int_t num2  = mXX[i]*mY[i] - mX[i]*mXY[i];
2024
2025     if (num2 < 0) {
2026       numAlu.AssignFormatted(-num2);
2027       numAlu.SetSign(-1);
2028     }
2029     else {
2030       numAlu.AssignFormatted(num2);
2031       numAlu.SetSign(1);
2032     }
2033
2034     numAlu    = numAlu/denomAlu;
2035
2036
2037     mOffset[i]  = numAlu.GetSignedValue();
2038     numAlu.SetSign(1);
2039     denomAlu.SetSign(1);
2040
2041
2042     //numAlu.AssignInt(mADC[i]+1);   // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)
2043     numAlu.AssignDouble((Double_t)mADC[i] + 1.5);      // numAlu has enough pre-comma place for that; correct trafo, best values
2044     mOffset[i]  = mOffset[i] + numAlu.GetValue();      // transform offset to a coord.system relative to chip; +1 to avoid neg. values
2045
2046     // up to here: adc-mapping according to TRAP manual and in line with pad-col mapping
2047     // reverse adc-counting to be again in line with the online mapping
2048     mADC[i]     = fNADC - 4 - mADC[i];                 // fNADC-4-mADC[i]: 0..17; remapping necessary;
2049     mADC[i]     = mADC[i] + 2;
2050     // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
2051   }
2052
2053   // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected;
2054   // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
2055
2056   // transform parameters to the local coordinate-system of a stack (used by GTU)
2057   AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
2058
2059   Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
2060   //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
2061
2062   // difference between width of inner and outer pads of a row is not accounted for;
2063
2064   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
2065   Double_t eCharge  = 0.3;                           // unit charge in (GeV/c)/m*T
2066   Double_t ptMin   = 2.3;                            // minimum transverse momentum (GeV/c); to be adjusted(?)
2067
2068   Double_t granularityOffset = 0.160;                // granularity for offset in mm
2069   Double_t granularitySlope  = 0.140;                // granularity for slope  in mm
2070
2071   // get the coordinates in SM-system; parameters:
2072
2073   Double_t zPos       =  (padPlane->GetRowPos(fRow))*10.0;  // z-position of the MCM; fRow is counted on a chamber; SM consists of 5
2074   // zPos is position of pad-borders;
2075   Double_t zOffset = 0.0;
2076   if ( fRow == 0 || fRow == 15 ) {
2077       zOffset = padPlane->GetLengthOPad();
2078   }
2079   else {
2080       zOffset = padPlane->GetLengthIPad();
2081   }
2082   zOffset = (-1.0) * zOffset/2.0;
2083   // turn zPos to be z-coordinate at middle of pad-row
2084   zPos = zPos + zOffset;
2085
2086
2087   Double_t xPos       =  0.0;                               // x-position of the upper border of the drift-chamber of actual layer
2088   Int_t    icol       =  0;                                 // column-number of adc-channel
2089   Double_t yPos[4];                                         // y-position of the pad to which ADC is connected
2090   Double_t dx         = 30.0;                               // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
2091   Double_t freqSample = fFeeParam->GetSamplingFrequency();  // retrieve the sampling frequency (10.019750 MHz)
2092   Double_t vdrift     = fCal->GetVdriftAverage(fChaId);     // averaged drift velocity for this detector (1.500000 cm/us)
2093   Int_t    nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
2094   Int_t    nrOfAmplTimeBins  = 2;                           // the number of time bins between anode wire and cathode wires in ampl.region (3.5mm)(guess)(suppose v_drift+3.5cm/us there=>all clusters arrive at anode wire within one time bin (100ns))
2095   Int_t    nrOfOffsetCorrTimeBins = tFS - nrOfAmplTimeBins - 1; // -1 is  to be conservative; offset correction will not remove the shift but is supposed to improve it; if tFS = 5, 2 drift time bins before tFS are assumed
2096   if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed
2097   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
2098   //Double_t lorAngle   =  7.0;                             // Lorentz-angle in degrees
2099   Double_t tiltAngle  = padPlane->GetTiltingAngle();        // sign-respecting tilting angle of pads in actual layer
2100   Double_t tiltTan    = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
2101   //Double_t lorTan     = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
2102
2103   Double_t alphaMax[4];                            // maximum deflection from the direction to the primary vertex; granularity of hit pads
2104   Double_t slopeMin[4];                            // local limits for the deflection
2105   Double_t slopeMax[4];
2106   Int_t   mslopeMin[4];                            // in granularity units; to be compared to mSlope[i]
2107   Int_t   mslopeMax[4];
2108
2109
2110   // x coord. of upper side of drift chambers in local SM-system (in mm)
2111   // obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
2112   // the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
2113   switch(fLayer)
2114     {
2115     case 0:
2116       xPos = 3003.0;
2117       break;
2118     case 1:
2119       xPos = 3129.0;
2120       break;
2121     case 2:
2122       xPos = 3255.0;
2123       break;
2124     case 3:
2125       xPos = 3381.0;
2126       break;
2127     case 4:
2128       xPos = 3507.0;
2129       break;
2130     case 5:
2131       xPos = 3633.0;
2132       break;
2133     }
2134
2135   // calculation of offset-correction n:
2136
2137   Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));
2138
2139   nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
2140   if (nCorrectOffset < 0) {
2141     numAlu.AssignInt(-nCorrectOffset);
2142     numAlu.SetSign(-1);
2143   }
2144   else {
2145     numAlu.AssignInt(nCorrectOffset);
2146     numAlu.SetSign(1);
2147   }
2148   nCorrectOffset = numAlu.GetSignedValue();
2149
2150   // the Lorentz correction to the offset
2151   Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
2152
2153
2154   lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
2155
2156   if(lorCorrectOffset < 0) {
2157       numAlu.AssignDouble(-lorCorrectOffset);
2158       numAlu.SetSign(-1);
2159   }
2160   else{
2161       numAlu.AssignDouble(lorCorrectOffset);
2162       numAlu.SetSign(1);
2163   }
2164
2165   Int_t mlorCorrectOffset = numAlu.GetSignedValue();
2166
2167
2168   Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
2169
2170   // calculation of slope-correction
2171
2172   // this is only true for tracks coming (approx.) from primary vertex
2173   // everything is evaluated for a tracklet covering the whole drift chamber
2174   Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
2175   // Double_t cCorrectSlope =  zPos/xPos*dx*tiltTan/granularitySlope;
2176   // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
2177
2178   if (cCorrectSlope < 0) {
2179       numAlu.AssignDouble(-cCorrectSlope);
2180       numAlu.SetSign(-1);
2181   }
2182   else {
2183       numAlu.AssignDouble(cCorrectSlope);
2184       numAlu.SetSign(1);
2185   }
2186   cCorrectSlope = numAlu.GetSignedValue();
2187
2188   // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
2189   // different pad-width of outer pads of a pad-plane not taken into account
2190   // 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)
2191   // 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
2192
2193
2194   Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope;  // >= 0.0
2195
2196   AliTRDtrapAlu correctAlu;
2197   correctAlu.Init(20,8);
2198
2199   AliTRDtrapAlu offsetAlu;
2200   offsetAlu.Init(13,0,-0x1000,0x0FFF);          // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
2201
2202   AliTRDtrapAlu slopeAlu;
2203   slopeAlu.Init(7,0,-0x40,0x3F);                // 7 bit-word;  2-complement (1 sign-bit);
2204
2205   for (Int_t i = 0; i<4; i++) {
2206
2207     if (mADC[i] == -1) continue;
2208
2209     icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
2210     yPos[i]   = (padPlane->GetColPos(icol))*10.0;
2211
2212
2213     // offset:
2214
2215     correctAlu.AssignDouble(mCorrectOffset);     // done because max. accuracy is 8 bit
2216     mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
2217     mOffset[i]  = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset));
2218     //mOffset[i]  = mOffset[i]*(-1);                   // adjust to direction of y-axes in online simulation
2219
2220     if (mOffset[i] < 0) {
2221       numAlu.AssignFormatted(-mOffset[i]);
2222       numAlu.SetSign(-1);
2223     }
2224     else {
2225       numAlu.AssignFormatted(mOffset[i]);
2226       numAlu.SetSign(1);
2227     }
2228
2229     offsetAlu = numAlu;
2230     mOffset[i]  = offsetAlu.GetSignedValue();
2231
2232
2233     // slope:
2234
2235     correctAlu.AssignDouble(mCorrectSlope);
2236     mCorrectSlope = correctAlu.GetValueWhole();
2237
2238     mSlope[i]   = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
2239
2240     if (mSlope[i] < 0) {
2241       numAlu.AssignFormatted(-mSlope[i]);
2242       numAlu.SetSign(-1);
2243     }
2244     else {
2245       numAlu.AssignFormatted(mSlope[i]);
2246       numAlu.SetSign(1);
2247     }
2248
2249     slopeAlu  = numAlu;     // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
2250     mSlope[i]   = slopeAlu.GetSignedValue();
2251
2252     // local (LTU) limits for the deflection
2253     // ATan returns angles in radian
2254     alphaMax[i]  = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
2255     slopeMin[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
2256     slopeMax[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
2257
2258     if (slopeMin[i] < 0) {
2259       slopeAlu.AssignDouble(-slopeMin[i]);
2260       slopeAlu.SetSign(-1);
2261     }
2262     else {
2263       slopeAlu.AssignDouble(slopeMin[i]);
2264       slopeAlu.SetSign(1);
2265     }
2266     mslopeMin[i] = slopeAlu.GetSignedValue();  // the borders should lie inside the range of mSlope -> usage of slopeAlu again
2267
2268     if (slopeMax[i] < 0) {
2269       slopeAlu.AssignDouble(-slopeMax[i]);
2270       slopeAlu.SetSign(-1);
2271     }
2272     else {
2273       slopeAlu.AssignDouble(slopeMax[i]);
2274       slopeAlu.SetSign(1);
2275     }
2276     mslopeMax[i] = slopeAlu.GetSignedValue();
2277   }
2278
2279   // suppress submission of tracks with low stiffness
2280   // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber)
2281
2282   // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
2283   // up to now they are sorted according to maximum hit sum
2284   // is the sorting really done in the TRAP-chip?
2285
2286   Int_t order[4] = {-1,-1,-1,-1};
2287   Int_t wordnr = 0;   // number of tracklet-words
2288
2289   for(Int_t j = 0; j < fMaxTracklets; j++) {
2290       //if( mADC[j] == -1) continue;
2291       if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
2292       wordnr++;
2293       if( wordnr-1 == 0) {
2294           order[0] = j;
2295           continue;
2296       }
2297       // wordnr-1>0, wordnr-1<4
2298       order[wordnr-1] = j;
2299       for( Int_t k = 0; k < wordnr-1; k++) {
2300           if( mOffset[j] < mOffset[order[k]] ) {
2301               for( Int_t l = wordnr-1; l > k; l-- ) {
2302                   order[l] = order[l-1];
2303               }
2304               order[k] = j;
2305               break;
2306           }
2307
2308       }
2309   }
2310
2311   // fill the bit-words in ascending order and without gaps
2312   UInt_t bitWord[4] = {0,0,0,0};                 // attention: unsigned int to have real 32 bits (no 2-complement)
2313   for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
2314       //Bool_t rem1 = kTRUE;
2315
2316     Int_t i = order[j];
2317     //bit-word is 2-complement and therefore without sign
2318     bitWord[j] =   1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
2319     //printf("\n");
2320     UInt_t shift  = 0;
2321     UInt_t shift2 = 0;
2322
2323
2324
2325
2326     /*printf("mean charge: %d\n",mMeanCharge[i]);
2327     printf("row: %d\n",fRow);
2328     printf("slope: %d\n",mSlope[i]);
2329     printf("pad position: %d\n",mOffset[i]);
2330     printf("channel: %d\n",mADC[i]);*/
2331
2332     // electron probability (currently not implemented; the mean charge is just scaled)
2333     shift = (UInt_t)mMeanCharge[i];
2334     for(Int_t iBit = 0; iBit < 8; iBit++) {
2335       bitWord[j]  = bitWord[j]<<1;
2336       bitWord[j] |= (shift>>(7-iBit))&1;
2337       //printf("0");
2338     }
2339
2340     // pad row
2341     shift = (UInt_t)fRow;
2342     for(Int_t iBit = 0; iBit < 4; iBit++) {
2343       bitWord[j]  = bitWord[j]<<1;
2344       bitWord[j] |= (shift>>(3-iBit))&1;
2345       //printf("%d", (fRow>>(3-iBit))&1);
2346     }
2347
2348     // deflection length
2349     if(mSlope[i] < 0) {
2350         shift = (UInt_t)(-mSlope[i]);
2351         // shift2 is 2-complement of shift
2352         shift2 = 1;
2353         for(Int_t iBit = 1; iBit < 7; iBit++) {
2354             shift2  = shift2<<1;
2355             shift2 |= (1-((shift)>>(6-iBit))&1);
2356             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2357         }
2358         shift2 = shift2 + 1;
2359         //printf("1");
2360         for(Int_t iBit = 0; iBit < 7; iBit++) {
2361             bitWord[j]  = bitWord[j]<<1;
2362             bitWord[j] |= (shift2>>(6-iBit))&1;
2363             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2364         }
2365     }
2366     else {
2367         shift = (UInt_t)(mSlope[i]);
2368         bitWord[j]  = bitWord[j]<<1;
2369         bitWord[j]   |= 0;
2370         //printf("0");
2371         for(Int_t iBit = 1; iBit < 7; iBit++) {
2372             bitWord[j]  = bitWord[j]<<1;
2373             bitWord[j] |= (shift>>(6-iBit))&1;
2374             //printf("%d",(mSlope[i]>>(6-iBit))&1);
2375         }
2376     }
2377
2378     // pad position
2379     if(mOffset[i] < 0) {
2380         shift = (UInt_t)(-mOffset[i]);
2381         shift2 = 1;
2382         for(Int_t iBit = 1; iBit < 13; iBit++) {
2383             shift2  = shift2<<1;
2384             shift2 |= (1-((shift)>>(12-iBit))&1);
2385             //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
2386         }
2387         shift2 = shift2 + 1;
2388         //printf("1");
2389         for(Int_t iBit = 0; iBit < 13; iBit++) {
2390             bitWord[j]  = bitWord[j]<<1;
2391             bitWord[j] |= (shift2>>(12-iBit))&1;
2392             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2393         }
2394     }
2395     else {
2396         shift = (UInt_t)mOffset[i];
2397         bitWord[j] = bitWord[j]<<1;
2398         bitWord[j]   |= 0;
2399         //printf("0");
2400         for(Int_t iBit = 1; iBit < 13; iBit++) {
2401             bitWord[j]  = bitWord[j]<<1;
2402             bitWord[j] |= (shift>>(12-iBit))&1;
2403             //printf("%d",(mOffset[i]>>(12-iBit))&1);
2404         }
2405     }
2406
2407
2408
2409     //printf("bitWord: %u\n",bitWord[j]);
2410     //printf("adc: %d\n",mADC[i]);
2411     fMCMT[j] = bitWord[j];
2412   }
2413
2414   //printf("\n");
2415
2416
2417   delete [] qsum;
2418   delete [] ieffped;
2419
2420   delete [] X;
2421   delete [] XX;
2422   delete [] Y;
2423   delete [] YY;
2424   delete [] XY;
2425   delete [] N;
2426   delete [] QT0;
2427   delete [] QT1;
2428
2429   delete [] hitSum;
2430   delete [] selectPair;
2431
2432   delete padPlane;
2433
2434 //if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
2435
2436   if (!fFeeParam->GetMCTrackletOutput()) return;
2437
2438
2439   // structure: in the current directory a root-file called "TRD_readout_tree.root" is stored with subdirectories SMxx/sx (supermodule, stack);
2440   // in each of these subdirectories 6 trees according to layers are saved, called lx;
2441   // 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;
2442   // branch-name: mcmxxxwd;
2443   // another branch contains the channel-number (mcmxxxch)
2444
2445
2446   AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
2447   AliLog::SetFileOutput("../log/tracklet.log");
2448
2449
2450   UInt_t* trackletWord;
2451   Int_t*  adcChannel;
2452
2453   Int_t u = 0;
2454
2455   // testing for wordnr in order to speed up the simulation
2456
2457   if (wordnr == 0 && fNextEvent == 0) {
2458     return;
2459   }
2460
2461
2462   Int_t mcmNr = fRobPos * (fGeo->MCMmax()) + fMcmPos;
2463
2464   Char_t* SMName    = new Char_t[4];
2465   Char_t* stackName = new Char_t[2];
2466   Char_t* layerName = new Char_t[2];
2467
2468   Char_t* treeName  = new Char_t[2];
2469   Char_t* treeTitle = new Char_t[8];
2470   Char_t* branchNameWd = new Char_t[8]; // mcmxxxwd bit word
2471   Char_t* branchNameCh = new Char_t[8]; // mcmxxxch channel
2472
2473   Char_t* dirName = NULL;
2474   Char_t* treeFile  = NULL;
2475   Char_t* evFile = NULL;
2476   Char_t* curDir = new Char_t[255];
2477
2478   for (Int_t i = 0; i<255; i++) {
2479     curDir[i]='n';
2480   }
2481   sprintf(curDir,"%s",gSystem->BaseName(gSystem->WorkingDirectory()));
2482   Int_t rawEvent = 0;
2483   Int_t nrPos = 3;
2484
2485
2486   while(curDir[nrPos]!='n'){
2487
2488
2489     switch(curDir[nrPos]) {
2490     case '0':
2491       rawEvent = rawEvent*10 + 0;
2492       break;
2493     case '1':
2494       rawEvent = rawEvent*10 + 1;
2495       break;
2496     case '2':
2497       rawEvent = rawEvent*10 + 2;
2498       break;
2499     case '3':
2500       rawEvent = rawEvent*10 + 3;
2501       break;
2502     case '4':
2503       rawEvent = rawEvent*10 + 4;
2504       break;
2505     case '5':
2506       rawEvent = rawEvent*10 + 5;
2507       break;
2508     case '6':
2509       rawEvent = rawEvent*10 + 6;
2510       break;
2511     case '7':
2512       rawEvent = rawEvent*10 + 7;
2513       break;
2514     case '8':
2515       rawEvent = rawEvent*10 + 8;
2516       break;
2517     case '9':
2518       rawEvent = rawEvent*10 + 9;
2519       break;
2520
2521     }
2522     nrPos = nrPos + 1;
2523   }
2524   delete [] curDir;
2525
2526   //if (!gSystem->ChangeDirectory("../TRD_Tracklet")) {
2527    // gSystem->MakeDirectory("../TRD_Tracklet");
2528     //gSystem->ChangeDirectory("../TRD_Tracklet");
2529     //}
2530
2531   gSystem->ChangeDirectory("..");
2532
2533
2534   TFile *f = new TFile("TRD_readout_tree.root","update");
2535   TTree *tree          = NULL;
2536   TBranch *branch      = NULL;
2537   TBranch *branchCh   = NULL;
2538
2539   Int_t iEventNr = 0;
2540   Int_t dignr = 10; // nr of digits of a integer
2541   Int_t space = 1;  // additional char-space
2542
2543
2544   evFile = new Char_t[2+space];
2545   sprintf(evFile,"ev%d",iEventNr);
2546
2547
2548   while(f->cd(evFile)){
2549     iEventNr = iEventNr + 1;
2550     if (iEventNr/dignr > 0) {
2551       dignr = dignr * 10;
2552       space = space + 1;
2553     }
2554     delete [] evFile;
2555     evFile = NULL;
2556     evFile = new Char_t[2+space];
2557     sprintf(evFile,"ev%d",iEventNr);
2558   }
2559
2560   if(iEventNr == rawEvent) { fNextEvent = 1; } // new event
2561
2562   if (fNextEvent == 1) {
2563     fNextEvent = 0;
2564     // turn to head directory
2565     f->mkdir(evFile);
2566     f->cd(evFile);
2567     // create all subdirectories and trees in case of new event
2568
2569
2570     for (Int_t iSector = 0; iSector < 18; iSector++) {
2571
2572       if (iSector < 10) {
2573         sprintf(SMName,"SM0%d",iSector);
2574       }
2575       else {
2576         sprintf(SMName,"SM%d",iSector);
2577       }
2578
2579       for (Int_t iStack = 0; iStack < 5; iStack++) {
2580         sprintf(stackName,"s%d",iStack);
2581
2582         f->cd(evFile);
2583         if (iStack == 0) {
2584           gDirectory->mkdir(SMName);
2585         }
2586         gDirectory->cd(SMName);
2587         gDirectory->mkdir(stackName);
2588         gDirectory->cd(stackName);
2589
2590         for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2591           sprintf(layerName,"l%d",iLayer);
2592           sprintf(treeName,"%s",layerName);
2593           sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName);
2594           tree = new TTree(treeName,treeTitle);
2595           tree->Write("",TObject::kOverwrite);
2596           delete tree;
2597           tree = NULL;
2598         }
2599       }
2600     }
2601
2602
2603   }
2604
2605
2606   else {
2607     iEventNr = iEventNr - 1;
2608     dignr = dignr/10;
2609     if (iEventNr/dignr == 0) space = space - 1;
2610     delete [] evFile;
2611     evFile = NULL;
2612     evFile = new Char_t[2+space];
2613     sprintf(evFile,"ev%d",iEventNr);
2614   }
2615
2616   if (wordnr == 0) {
2617     delete [] SMName;
2618     delete [] stackName;
2619     delete [] layerName;
2620     delete [] treeName;
2621     delete [] treeTitle;
2622     delete [] branchNameWd;
2623     delete [] branchNameCh;
2624     delete [] evFile;
2625     f->Close();
2626     dirName   = new Char_t[6+space];
2627     //sprintf(dirName,"../raw%d",iEventNr);
2628     sprintf(dirName,"raw%d",iEventNr);
2629     gSystem->ChangeDirectory(dirName);
2630     delete [] dirName;
2631     return;
2632   }
2633
2634   dirName   = new Char_t[6+space];
2635   //sprintf(dirName,"../raw%d",iEventNr);
2636   sprintf(dirName,"raw%d",iEventNr);
2637
2638   f->cd(evFile);
2639
2640   if (fSector < 10) {
2641     sprintf(SMName,"SM0%d",fSector);
2642   }
2643   else {
2644     sprintf(SMName,"SM%d",fSector);
2645   }
2646   sprintf(stackName,"s%d",fStack);
2647   sprintf(layerName,"l%d",fLayer);
2648   sprintf(treeName,"%s",layerName);
2649   sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName);
2650
2651   treeFile = new Char_t[13+space];
2652   sprintf(treeFile,"%s/%s/%s/%s",evFile,SMName,stackName,treeName);
2653   delete [] evFile;
2654   evFile = NULL;
2655   gDirectory->cd(SMName);
2656   gDirectory->cd(stackName);
2657   tree = (TTree*)f->Get(treeFile);
2658   delete [] treeFile;
2659   treeFile = NULL;
2660
2661
2662   //make branch with number of words and fill
2663
2664   if (mcmNr < 10) {
2665     sprintf(branchNameWd,"mcm00%dwd",mcmNr);
2666     sprintf(branchNameCh,"mcm00%dch",mcmNr);
2667   }
2668   else {
2669     if (mcmNr < 100) {
2670       sprintf(branchNameWd,"mcm0%dwd",mcmNr);
2671       sprintf(branchNameCh,"mcm0%dch",mcmNr);
2672     }
2673     else {
2674       sprintf(branchNameWd,"mcm%dwd",mcmNr);
2675       sprintf(branchNameCh,"mcm%dch",mcmNr);
2676     }
2677   }
2678
2679
2680
2681   // fill the tracklet word; here wordnr > 0
2682
2683   trackletWord = new UInt_t[fMaxTracklets];
2684   adcChannel   = new Int_t[fMaxTracklets];
2685
2686   for (Int_t j = 0; j < fMaxTracklets; j++) {
2687       Int_t i = order[j];
2688       trackletWord[j] = 0;
2689       adcChannel[j] = -1;
2690       if (bitWord[j]!=0) {
2691           trackletWord[u] = bitWord[j];
2692           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
2693
2694           //fMCMT[u] = bitWord[j];
2695           u = u + 1;
2696       }
2697   }
2698
2699   branch = tree->GetBranch(branchNameWd);
2700   if(!branch) {
2701     //make branch and fill
2702     branch = tree->Branch(branchNameWd,trackletWord,"trackletWord[4]/i"); // 32 bit unsigned integer
2703     branch->Fill();
2704   }
2705
2706   branchCh = tree->GetBranch(branchNameCh);
2707   if(!branchCh) {
2708     //make branch and fill
2709     branchCh = tree->Branch(branchNameCh,adcChannel,"adcChannel[4]/i"); // 32 bit unsigned integer
2710     branchCh->Fill();
2711   }
2712
2713
2714   tree->Write("",TObject::kOverwrite);
2715
2716
2717   delete [] SMName;
2718   delete [] stackName;
2719   delete [] layerName;
2720   delete [] treeName;
2721   delete [] treeTitle;
2722   delete [] branchNameWd;
2723   delete [] branchNameCh;
2724   delete [] trackletWord;
2725   delete [] adcChannel;
2726
2727   f->Close();
2728   gSystem->ChangeDirectory(dirName);
2729   delete [] dirName;
2730
2731
2732
2733   // to be done:
2734   // error measure for quality of fit (not necessarily needed for the trigger)
2735   // cluster quality threshold (not yet set)
2736   // electron probability
2737
2738
2739
2740 }
2741
2742
2743
2744
2745
2746