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