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