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